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