1
2! Copyright (C) 2017 T. Mueller, J. K. Dewhurst, S. Sharma and E. K. U. Gross.
3! This file is distributed under the terms of the GNU General Public License.
4! See the file COPYING for license details.
5
6subroutine rhomaguk(ik0,lock,evecu)
7use modmain
8use modulr
9use modomp
10implicit none
11! arguments
12integer, intent(in) :: ik0
13integer(8), intent(in) :: lock(nqpt)
14complex(8), intent(in) :: evecu(nstulr,nstulr)
15! local variables
16integer ik,ikpa,ist,i,j
17integer ngk0,ispn,is,ias
18integer npc,ir,nthd
19real(8) wm,wi
20! automatic arrays
21integer idx(nstsv)
22complex(8) zfft(nqpt),wfir(ngtc,nspinor)
23! allocatable arrays
24integer(8), allocatable :: lockl(:)
25complex(8), allocatable :: apwalm(:,:,:,:),evecfv(:,:),evecsv(:,:)
26complex(8), allocatable :: wfmt(:,:,:,:),wfgk(:,:,:)
27complex(8), allocatable :: wfrmt(:,:,:,:),wfrgk(:,:,:)
28! central k-point
29ik=(ik0-1)*nkpa+1
30! number of G+k-vectors for central k-point
31ngk0=ngk(1,ik)
32! initialise the local OpenMP locks
33allocate(lockl(nqpt))
34do ir=1,nqpt
35  call omp_init_lock(lockl(ir))
36end do
37! get the eigenvectors from file
38allocate(evecfv(nmatmax,nstfv),evecsv(nstsv,nstsv))
39call getevecfv(filext,ik,vkl(:,ik),vgkl(:,:,:,ik),evecfv)
40call getevecsv(filext,ik,vkl(:,ik),evecsv)
41! find the matching coefficients
42allocate(apwalm(ngkmax,apwordmax,lmmaxapw,natmtot))
43call match(ngk0,vgkc(:,:,1,ik),gkc(:,1,ik),sfacgk(:,:,1,ik),apwalm)
44! index to all states
45do ist=1,nstsv
46  idx(ist)=ist
47end do
48allocate(wfmt(npcmtmax,natmtot,nspinor,nstsv),wfgk(ngk0,nspinor,nstsv))
49call genwfsv(.false.,.true.,nstsv,idx,ngdgc,igfc,ngk0,igkig(:,1,ik),apwalm, &
50 evecfv,evecsv,wfmt,ngk0,wfgk)
51deallocate(apwalm,evecfv,evecsv)
52allocate(wfrmt(npcmtmax,natmtot,nspinor,nqpt),wfrgk(ngk0,nspinor,nqpt))
53! loop over ultra long-range states
54do j=1,nstulr
55  wm=occulr(j,ik0)
56  if (abs(wm).lt.epsocc) cycle
57  wm=wm*wkpt(ik)
58  wi=wm/omega
59! zero the ultra long-range wavefunctions
60  call holdthd(2,nthd)
61!$OMP PARALLEL SECTIONS DEFAULT(SHARED) &
62!$OMP PRIVATE(ir,ispn) &
63!$OMP NUM_THREADS(nthd)
64!$OMP SECTION
65  do ir=1,nqpt
66    do ispn=1,nspinor
67      call wfmt0(wfrmt(:,:,ispn,ir))
68    end do
69  end do
70!$OMP SECTION
71  wfrgk(:,:,:)=0.d0
72!$OMP END PARALLEL SECTIONS
73  call freethd(nthd)
74! parallel loop over second-variational states
75  call holdthd(nstsv,nthd)
76!$OMP PARALLEL DO DEFAULT(SHARED) &
77!$OMP PRIVATE(zfft,ikpa,i,ir) &
78!$OMP NUM_THREADS(nthd)
79  do ist=1,nstsv
80    zfft(:)=0.d0
81! loop over kappa-points
82    do ikpa=1,nkpa
83      i=(ikpa-1)*nstsv+ist
84! store the wavefunction in Q-space
85      zfft(iqfft(ikpa))=evecu(i,j)
86    end do
87! Fourier transform to R-space
88    call zfftifc(3,ngridq,1,zfft)
89! loop over R-points
90    do ir=1,nqpt
91      call omp_set_lock(lockl(ir))
92      call wfadd(zfft(ir),wfmt(:,:,:,ist),wfgk(:,:,ist),wfrmt(:,:,:,ir), &
93       wfrgk(:,:,ir))
94      call omp_unset_lock(lockl(ir))
95    end do
96  end do
97!$OMP END PARALLEL DO
98  call freethd(nthd)
99! parallel loop over R-points
100  call holdthd(nqpt,nthd)
101!$OMP PARALLEL DO DEFAULT(SHARED) &
102!$OMP PRIVATE(wfir,ispn,ias,is,npc) &
103!$OMP NUM_THREADS(nthd)
104  do ir=1,nqpt
105    do ispn=1,nspinor
106! Fourier transform the interstitial part to real-space
107      wfir(:,ispn)=0.d0
108      wfir(igfc(igkig(1:ngk0,1,ik)),ispn)=wfrgk(1:ngk0,ispn,ir)
109      call zfftifc(3,ngdgc,1,wfir(:,ispn))
110    end do
111! add to the density and magnetisation
112    call omp_set_lock(lock(ir))
113    do ias=1,natmtot
114      is=idxis(ias)
115      npc=npcmt(is)
116      if (spinpol) then
117        if (ncmag) then
118          call rmk1(npc,wm,wfrmt(:,ias,1,ir),wfrmt(:,ias,2,ir), &
119           rhormt(:,ias,ir),magrmt(:,ias,1,ir),magrmt(:,ias,2,ir), &
120           magrmt(:,ias,3,ir))
121        else
122          call rmk2(npc,wm,wfrmt(:,ias,1,ir),wfrmt(:,ias,2,ir), &
123           rhormt(:,ias,ir),magrmt(:,ias,1,ir))
124        end if
125      else
126        call rmk3(npc,wm,wfrmt(:,ias,1,ir),rhormt(:,ias,ir))
127      end if
128    end do
129    if (spinpol) then
130      if (ncmag) then
131        call rmk1(ngtc,wi,wfir,wfir(:,2),rhorir(:,ir),magrir(:,1,ir), &
132         magrir(:,2,ir),magrir(:,3,ir))
133      else
134        call rmk2(ngtc,wi,wfir,wfir(:,2),rhorir(:,ir),magrir(:,1,ir))
135      end if
136    else
137      call rmk3(ngtc,wi,wfir,rhorir(:,ir))
138    end if
139    call omp_unset_lock(lock(ir))
140! end loop over R-points
141  end do
142!$OMP END PARALLEL DO
143  call freethd(nthd)
144! end loop over long-range states
145end do
146! destroy the local OpenMP locks
147do ir=1,nqpt
148  call omp_destroy_lock(lockl(ir))
149end do
150deallocate(lockl)
151deallocate(wfmt,wfgk,wfrmt,wfrgk)
152return
153
154contains
155
156pure subroutine wfmt0(wfmt)
157implicit none
158! arguments
159complex(8), intent(out) :: wfmt(npcmtmax,natmtot)
160! local variables
161integer is,ias
162do ias=1,natmtot
163  is=idxis(ias)
164  wfmt(1:npcmt(is),ias)=0.d0
165end do
166end subroutine
167
168subroutine wfadd(za,wfmt1,wfgk1,wfmt2,wfgk2)
169implicit none
170! arguments
171complex(8), intent(in) :: za
172complex(8), intent(in) :: wfmt1(npcmtmax,natmtot,nspinor)
173complex(8), intent(in) :: wfgk1(ngk0,nspinor)
174complex(8), intent(inout) :: wfmt2(npcmtmax,natmtot,nspinor)
175complex(8), intent(inout) :: wfgk2(ngk0,nspinor)
176! local variables
177integer ispn,is,ias
178do ispn=1,nspinor
179  do ias=1,natmtot
180    is=idxis(ias)
181    call zaxpy(npcmt(is),za,wfmt1(:,ias,ispn),1,wfmt2(:,ias,ispn),1)
182  end do
183end do
184call zaxpy(ngk0*nspinor,za,wfgk1,1,wfgk2,1)
185end subroutine
186
187pure subroutine rmk1(n,wo,wf1,wf2,rho,mag1,mag2,mag3)
188implicit none
189! arguments
190integer, intent(in) :: n
191real(8), intent(in) :: wo
192complex(8), intent(in) :: wf1(n),wf2(n)
193real(8), intent(inout) :: rho(n),mag1(n),mag2(n),mag3(n)
194! local variables
195integer i
196real(8) wo2,t1,t2
197real(8) a1,b1,a2,b2
198wo2=2.d0*wo
199!$OMP SIMD PRIVATE(a1,b1,a2,b2,t1,t2) SIMDLEN(8)
200do i=1,n
201  a1=dble(wf1(i)); b1=aimag(wf1(i))
202  a2=dble(wf2(i)); b2=aimag(wf2(i))
203  t1=a1**2+b1**2; t2=a2**2+b2**2
204  mag1(i)=mag1(i)+wo2*(a1*a2+b1*b2)
205  mag2(i)=mag2(i)+wo2*(a1*b2-b1*a2)
206  mag3(i)=mag3(i)+wo*(t1-t2)
207  rho(i)=rho(i)+wo*(t1+t2)
208end do
209end subroutine
210
211pure subroutine rmk2(n,wo,wf1,wf2,rho,mag)
212implicit none
213! arguments
214integer, intent(in) :: n
215real(8), intent(in) :: wo
216complex(8), intent(in) :: wf1(n),wf2(n)
217real(8), intent(inout) :: rho(n),mag(n)
218! local variables
219integer i
220real(8) t1,t2
221!$OMP SIMD PRIVATE(t1,t2) SIMDLEN(8)
222do i=1,n
223  t1=dble(wf1(i))**2+aimag(wf1(i))**2
224  t2=dble(wf2(i))**2+aimag(wf2(i))**2
225  mag(i)=mag(i)+wo*(t1-t2)
226  rho(i)=rho(i)+wo*(t1+t2)
227end do
228end subroutine
229
230pure subroutine rmk3(n,wo,wf,rho)
231implicit none
232! arguments
233integer, intent(in) :: n
234real(8), intent(in) :: wo
235complex(8), intent(in) :: wf(n)
236real(8), intent(inout) :: rho(n)
237rho(:)=rho(:)+wo*(dble(wf(:))**2+aimag(wf(:))**2)
238end subroutine
239
240end subroutine
241
242