File: rhomaguk.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (262 lines) | stat: -rw-r--r-- 6,837 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262

! Copyright (C) 2017 T. Mueller, J. K. Dewhurst, S. Sharma and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine rhomaguk(ik0,lock,evecu)
use modmain
use modulr
use modomp
implicit none
! arguments
integer, intent(in) :: ik0
integer(8), intent(in) :: lock(nqpt)
complex(8), intent(in) :: evecu(nstulr,nstulr)
! local variables
integer ik,ist,ikpa
integer ngk0,igk,ifg
integer ispn,is,ias
integer npc,ir,i,j
integer nthd,ithd
real(8) wm,wi
! automatic arrays
integer idx(nstsv)
! allocatable arrays
integer(8), allocatable :: lockl(:)
complex(8), allocatable :: apwalm(:,:,:,:),evecfv(:,:),evecsv(:,:)
complex(8), allocatable :: wfmt(:,:,:,:),wfgk(:,:,:)
complex(8), allocatable :: wfmtu(:,:,:,:),wfgku(:,:,:)
complex(8), allocatable :: wfir(:,:),zfft(:,:)
! central k-point
ik=(ik0-1)*nkpa+1
! number of G+k-vectors for central k-point
ngk0=ngk(1,ik)
! initialise the local OpenMP locks
allocate(lockl(nqpt))
do ir=1,nqpt
  call omp_init_lock(lockl(ir))
end do
! get the eigenvectors from file
allocate(evecfv(nmatmax,nstfv),evecsv(nstsv,nstsv))
call getevecfv(filext,ik,vkl(:,ik),vgkl(:,:,:,ik),evecfv)
call getevecsv(filext,ik,vkl(:,ik),evecsv)
! find the matching coefficients
allocate(apwalm(ngkmax,apwordmax,lmmaxapw,natmtot))
call match(ngk0,gkc(:,1,ik),tpgkc(:,:,1,ik),sfacgk(:,:,1,ik),apwalm)
! index to all states
do ist=1,nstsv
  idx(ist)=ist
end do
allocate(wfmt(npcmtmax,natmtot,nspinor,nstsv),wfgk(ngk0,nspinor,nstsv))
call genwfsv(.false.,.true.,nstsv,idx,ngridg,igfft,ngk(1,ik),igkig(:,1,ik), &
 apwalm,evecfv,evecsv,wfmt,ngk0,wfgk)
deallocate(apwalm,evecfv,evecsv)
allocate(wfmtu(npcmtmax,natmtot,nspinor,nqpt),wfgku(ngk0,nspinor,nqpt))
! loop over ultra long-range states
do j=1,nstulr
  wm=wkpt(ik)*occulr(j,ik0)
  wi=wm/omega
  if (abs(wm).lt.epsocc) cycle
! zero the ultra long-range wavefunctions
  call omp_hold(2,nthd)
!$OMP PARALLEL SECTIONS DEFAULT(SHARED) &
!$OMP PRIVATE(ir,ispn) &
!$OMP NUM_THREADS(nthd)
!$OMP SECTION
  do ir=1,nqpt
    do ispn=1,nspinor
      call wfmt0(wfmtu(:,:,ispn,ir))
    end do
  end do
!$OMP SECTION
  wfgku(:,:,:)=0.d0
!$OMP END PARALLEL SECTIONS
  call omp_free(nthd)
! parallel loop over second-variational states
  call omp_hold(nstsv,nthd)
  allocate(zfft(nqpt,0:nthd-1))
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(ithd,ikpa,i,ir) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do ist=1,nstsv
    ithd=omp_get_thread_num()
    zfft(:,ithd)=0.d0
! loop over kappa-points
    do ikpa=1,nkpa
      i=(ikpa-1)*nstsv+ist
! store the wavefunction in Q-space
      zfft(iqfft(ikpa),ithd)=evecu(i,j)
    end do
! Fourier transform to R-space
    call zfftifc(3,ngridq,1,zfft(:,ithd))
! loop over R-points
    do ir=1,nqpt
      call omp_set_lock(lockl(ir))
      call wfadd(zfft(ir,ithd),wfmt(:,:,:,ist),wfgk(:,:,ist),wfmtu(:,:,:,ir), &
       wfgku(:,:,ir))
      call omp_unset_lock(lockl(ir))
    end do
  end do
!$OMP END DO
!$OMP END PARALLEL
  deallocate(zfft)
  call omp_free(nthd)
! parallel loop over R-points
  call omp_hold(nqpt,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(wfir,ispn,igk,ifg) &
!$OMP PRIVATE(ias,is,npc) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
  do ir=1,nqpt
    allocate(wfir(ngtot,nspinor))
    do ispn=1,nspinor
! Fourier transform the interstitial part to real-space
      wfir(:,ispn)=0.d0
      do igk=1,ngk0
        ifg=igfft(igkig(igk,1,ik))
        wfir(ifg,ispn)=wfgku(igk,ispn,ir)
      end do
      call zfftifc(3,ngridg,1,wfir(:,ispn))
    end do
! add to the density and magnetisation
    call omp_set_lock(lock(ir))
    do ias=1,natmtot
      is=idxis(ias)
      npc=npcmt(is)
      if (spinpol) then
        if (ncmag) then
          call rmk1(npc,wm,wfmtu(:,ias,1,ir),wfmtu(:,ias,2,ir), &
           rhormt(:,ias,ir),magrmt(:,ias,1,ir),magrmt(:,ias,2,ir), &
           magrmt(:,ias,3,ir))
        else
          call rmk2(npc,wm,wfmtu(:,ias,1,ir),wfmtu(:,ias,2,ir), &
           rhormt(:,ias,ir),magrmt(:,ias,1,ir))
        end if
      else
        call rmk3(npc,wm,wfmtu(:,ias,1,ir),rhormt(:,ias,ir))
      end if
    end do
    if (spinpol) then
      if (ncmag) then
        call rmk1(ngtot,wi,wfir(:,1),wfir(:,2),rhorir(:,ir),magrir(:,1,ir), &
         magrir(:,2,ir),magrir(:,3,ir))
      else
        call rmk2(ngtot,wi,wfir(:,1),wfir(:,2),rhorir(:,ir),magrir(:,1,ir))
      end if
    else
      call rmk3(ngtot,wi,wfir(:,1),rhorir(:,ir))
    end if
    call omp_unset_lock(lock(ir))
    deallocate(wfir)
! end loop over R-points
  end do
!$OMP END DO
!$OMP END PARALLEL
  call omp_free(nthd)
! end loop over long-range states
end do
! destroy the local OpenMP locks
do ir=1,nqpt
  call omp_destroy_lock(lockl(ir))
end do
deallocate(lockl)
deallocate(wfmt,wfgk,wfmtu,wfgku)
return

contains

subroutine wfmt0(wfmt)
implicit none
! arguments
complex(8), intent(out) :: wfmt(npcmtmax,natmtot)
! local variables
integer is,ias
do ias=1,natmtot
  is=idxis(ias)
  wfmt(1:npcmt(is),ias)=0.d0
end do
return
end subroutine

subroutine wfadd(za,wfmt1,wfgk1,wfmt2,wfgk2)
implicit none
! arguments
complex(8), intent(in) :: za
complex(8), intent(in) :: wfmt1(npcmtmax,natmtot,nspinor)
complex(8), intent(in) :: wfgk1(ngk0,nspinor)
complex(8), intent(inout) :: wfmt2(npcmtmax,natmtot,nspinor)
complex(8), intent(inout) :: wfgk2(ngk0,nspinor)
! local variables
integer ispn,is,ias
do ispn=1,nspinor
  do ias=1,natmtot
    is=idxis(ias)
    call zaxpy(npcmt(is),za,wfmt1(:,ias,ispn),1,wfmt2(:,ias,ispn),1)
  end do
end do
do ispn=1,nspinor
  call zaxpy(ngk0,za,wfgk1(:,ispn),1,wfgk2(:,ispn),1)
end do
return
end subroutine

subroutine rmk1(n,wo,wf1,wf2,rho,mag1,mag2,mag3)
implicit none
! arguments
integer, intent(in) :: n
real(8), intent(in) :: wo
complex(8), intent(in) :: wf1(n),wf2(n)
real(8), intent(inout) :: rho(n),mag1(n),mag2(n),mag3(n)
! local variables
integer i
real(8) wo2,t1,t2
complex(8) z1,z2
wo2=2.d0*wo
do i=1,n
  z1=wf1(i)
  z2=wf2(i)
  t1=dble(z1)**2+aimag(z1)**2
  t2=dble(z2)**2+aimag(z2)**2
  z1=conjg(z1)*z2
  rho(i)=rho(i)+wo*(t1+t2)
  mag1(i)=mag1(i)+wo2*dble(z1)
  mag2(i)=mag2(i)+wo2*aimag(z1)
  mag3(i)=mag3(i)+wo*(t1-t2)
end do
return
end subroutine

subroutine rmk2(n,wo,wf1,wf2,rho,mag)
implicit none
! arguments
integer, intent(in) :: n
real(8), intent(in) :: wo
complex(8), intent(in) :: wf1(n),wf2(n)
real(8), intent(inout) :: rho(n),mag(n)
! local variables
integer i
real(8) t1,t2
do i=1,n
  t1=dble(wf1(i))**2+aimag(wf1(i))**2
  t2=dble(wf2(i))**2+aimag(wf2(i))**2
  rho(i)=rho(i)+wo*(t1+t2)
  mag(i)=mag(i)+wo*(t1-t2)
end do
return
end subroutine

subroutine rmk3(n,wo,wf,rho)
implicit none
! arguments
integer, intent(in) :: n
real(8), intent(in) :: wo
complex(8), intent(in) :: wf(n)
real(8), intent(inout) :: rho(n)
rho(:)=rho(:)+wo*(dble(wf(:))**2+aimag(wf(:))**2)
return
end subroutine

end subroutine