TABLE OF CONTENTS
ABINIT/xcden [ Functions ]
NAME
xcden
FUNCTION
Prepare density data before calling local or semi-local xc evaluators.
NOTES
Can take into account a shift of the grid, for purpose of better accuracy. Can also compute the gradient of the density, for use with GGAs. Also eliminate eventual negative densities.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, DRH) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cplex=if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1) ishift : if ==0, do not shift the xc grid (usual case); if ==1, shift the xc grid nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ngrad : =1, only compute the density ; =2 also compute the gradient of the density. Note : ngrad**2 is also used to dimension rhonow nspden=number of spin-density components qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2). rhor(cplex*nfft,nspden)=real space electron density in electrons/bohr**3, on the unshifted grid (total in first half and spin-up in second half if nspden=2)
OUTPUT
rhonow(cplex*nfft,nspden,ngrad*ngrad)=electron (spin)-density in real space and eventually its gradient, either on the unshifted grid (if ishift==0, then equal to rhor),or on the shifted grid rhonow(:,:,1)=electron density in electrons/bohr**3 if ngrad==2 : rhonow(:,:,2:4)=gradient of electron density in electrons/bohr**4 OPTIONAL OUTPUT lrhonow(cplex*nfft,nspden)=Laplacian of the electron (spin)-density in real space in electrons/bohr**5 (in case of meta GGA)
PARENTS
afterscfloop,dfpt_mkvxcgga,dfpt_mkvxcstrgga,gammapositron_fft rhohxcpositron,rhotoxc
CHILDREN
fourdp,phase,ptabs_fourdp,timab
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 subroutine xcden (cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,paral_kgb,qphon,rhor,rhonow, & !Mandatory arguments 60 & lrhonow) !Optional arguments 61 62 use defs_basis 63 use defs_abitypes 64 use m_profiling_abi 65 use m_errors 66 67 use m_mpinfo, only : ptabs_fourdp 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'xcden' 73 use interfaces_18_timing 74 use interfaces_53_ffts 75 use interfaces_56_xc, except_this_one => xcden 76 !End of the abilint section 77 78 implicit none 79 80 !Arguments ------------------------------------ 81 !scalars 82 integer,intent(in) :: cplex,ishift,nfft,ngrad,nspden,paral_kgb 83 type(MPI_type),intent(in) :: mpi_enreg 84 !arrays 85 integer,intent(in) :: ngfft(18) 86 real(dp),intent(in) :: gprimd(3,3),qphon(3),rhor(cplex*nfft,nspden) 87 real(dp),intent(out) :: rhonow(cplex*nfft,nspden,ngrad*ngrad) 88 real(dp),intent(out),optional :: lrhonow(cplex*nfft,nspden) 89 90 !Local variables------------------------------- 91 !scalars 92 integer :: i1,i2,i3,id1,id2,id3,idir,ifft,ig1,ig2,ig3 93 integer :: ispden,n1,n2,n3,qeq0 94 real(dp) :: gc23_idir,gcart_idir,ph123i,ph123r,ph1i,ph1r,ph23i,ph23r,ph2i,ph2r 95 real(dp) :: ph3i,ph3r,work_im,work_re 96 character(len=500) :: message 97 !arrays 98 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 99 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 100 real(dp) :: tsec(2) 101 real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:),ph1(:),ph2(:),ph3(:) 102 real(dp),allocatable :: wkcmpx(:,:),work(:),workgr(:,:),worklp(:,:) 103 104 ! ************************************************************************* 105 106 !DEBUG 107 !write(std_out,*)' xcden : enter ' 108 !ENDDEBUG 109 110 if (ishift/=0 .and. ishift/=1) then 111 write(message, '(a,i0)' )'ishift must be 0 or 1 ; input was',ishift 112 MSG_BUG(message) 113 end if 114 115 if (ngrad/=1 .and. ngrad/=2) then 116 write(message, '(a,i0)' )'ngrad must be 1 or 2 ; input was',ngrad 117 MSG_BUG(message) 118 end if 119 120 !Keep local copy of fft dimensions 121 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 122 123 !Initialize computation of G in cartesian coordinates 124 id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2 125 126 !Get the distrib associated with this fft_grid 127 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 128 129 !Check whether q=0 130 qeq0=0 131 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1 132 133 if(ishift==0)then 134 135 ! Copy the input rhor in the new location. Will check on negative values later 136 137 do ispden=1,nspden 138 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,nfft,rhonow,rhor) 139 do ifft=1,cplex*nfft 140 rhonow(ifft,ispden,1)=rhor(ifft,ispden) 141 end do 142 end do 143 144 end if 145 146 if(ishift==1 .or. ngrad==2)then 147 148 ABI_ALLOCATE(wkcmpx,(2,nfft)) 149 ABI_ALLOCATE(work,(cplex*nfft)) 150 151 if(ishift==1)then 152 ! Precompute phases (The phases correspond to a shift of density on real space 153 ! grid from center at 0 0 0 to (1/2)*(1/n1,1/n2,1/n3).) 154 ABI_ALLOCATE(ph1,(2*n1)) 155 ABI_ALLOCATE(ph2,(2*n2)) 156 ABI_ALLOCATE(ph3,(2*n3)) 157 call phase(n1,ph1) 158 call phase(n2,ph2) 159 call phase(n3,ph3) 160 end if 161 162 do ispden=1,nspden 163 164 ! Obtain rho(G) in wkcmpx from input rho(r) 165 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,nfft,rhor,work) 166 do ifft=1,cplex*nfft 167 work(ifft)=rhor(ifft,ispden) 168 end do 169 170 call timab(82,1,tsec) 171 call fourdp(cplex,wkcmpx,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 172 call timab(82,2,tsec) 173 174 ! If shift is required, multiply now rho(G) by phase, then generate rho(r+delta) 175 if(ishift==1)then 176 !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,ph1i,ph1r,ph123i,ph123r,ph2i,ph2r,ph23i,ph23r,ph3i,ph3r,work_im,work_re) & 177 !$OMP&SHARED(n1,n2,n3,ph1,ph2,ph3,wkcmpx,mpi_enreg,fftn2_distrib) 178 do i3=1,n3 179 ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft) 180 ph3r=ph3(2*i3-1) 181 ph3i=ph3(2*i3 ) 182 do i2=1,n2 183 ph2r=ph2(2*i2-1) 184 ph2i=ph2(2*i2 ) 185 ph23r=ph2r*ph3r-ph2i*ph3i 186 ph23i=ph2i*ph3r+ph2r*ph3i 187 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 188 do i1=1,n1 189 ifft=ifft+1 190 ph1r=ph1(2*i1-1) 191 ph1i=ph1(2*i1 ) 192 ph123r=ph1r*ph23r-ph1i*ph23i 193 ph123i=ph1i*ph23r+ph1r*ph23i 194 ! Must use intermediate variables ! 195 work_re=ph123r*wkcmpx(1,ifft)-ph123i*wkcmpx(2,ifft) 196 work_im=ph123i*wkcmpx(1,ifft)+ph123r*wkcmpx(2,ifft) 197 wkcmpx(1,ifft)=work_re 198 wkcmpx(2,ifft)=work_im 199 end do 200 end if 201 end do 202 end do 203 call timab(82,1,tsec) 204 call fourdp(cplex,wkcmpx,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 205 call timab(82,2,tsec) 206 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,rhonow,work) 207 do ifft=1,cplex*nfft 208 rhonow(ifft,ispden,1)=work(ifft) 209 end do 210 end if 211 212 ! If gradient of the density is required, take care of the three components now 213 ! Note : this operation is applied on the eventually shifted rho(G) 214 if(ngrad==2)then 215 ABI_ALLOCATE(gcart1,(n1)) 216 ABI_ALLOCATE(gcart2,(n2)) 217 ABI_ALLOCATE(gcart3,(n3)) 218 ABI_ALLOCATE(workgr,(2,nfft)) 219 if (present(lrhonow)) then 220 ABI_ALLOCATE(worklp,(2,nfft)) 221 lrhonow(:,ispden)=zero 222 end if 223 do idir=1,3 224 225 do i1=1,n1 226 ig1=i1-(i1/id1)*n1-1 227 gcart1(i1)=gprimd(idir,1)*two_pi*(dble(ig1)+qphon(1)) 228 end do 229 ! Note that the G <-> -G symmetry must be maintained 230 if(mod(n1,2)==0 .and. qeq0==1)gcart1(n1/2+1)=zero 231 do i2=1,n2 232 ig2=i2-(i2/id2)*n2-1 233 gcart2(i2)=gprimd(idir,2)*two_pi*(dble(ig2)+qphon(2)) 234 end do 235 if(mod(n2,2)==0 .and. qeq0==1)gcart2(n2/2+1)=zero 236 do i3=1,n3 237 ig3=i3-(i3/id3)*n3-1 238 gcart3(i3)=gprimd(idir,3)*two_pi*(dble(ig3)+qphon(3)) 239 end do 240 if(mod(n3,2)==0 .and. qeq0==1)gcart3(n3/2+1)=zero 241 242 ! MG: Be careful here with OMP due to ifft. Disabled for the time being. 243 ! !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,gcart_idir,gc23_idir) & 244 ! !$OMP&SHARED(gcart1,gcart2,gcart3,n1,n2,n3,wkcmpx,workgr) 245 ifft = 0 246 do i3=1,n3 247 do i2=1,n2 248 gc23_idir=gcart2(i2)+gcart3(i3) 249 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 250 do i1=1,n1 251 ifft=ifft+1 252 gcart_idir=gc23_idir+gcart1(i1) 253 ! Multiply by i 2pi G(idir) 254 workgr(2,ifft)= gcart_idir*wkcmpx(1,ifft) 255 workgr(1,ifft)=-gcart_idir*wkcmpx(2,ifft) 256 ! Do the same to the gradient in order to get the laplacian 257 if (present(lrhonow)) then 258 worklp(2,ifft)= gcart_idir*workgr(1,ifft) 259 worklp(1,ifft)=-gcart_idir*workgr(2,ifft) 260 end if 261 end do 262 end if 263 end do 264 end do 265 call timab(82,1,tsec) 266 call fourdp(cplex,workgr,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 267 call timab(82,2,tsec) 268 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(idir,ispden,cplex,nfft,rhonow,work) 269 do ifft=1,cplex*nfft 270 rhonow(ifft,ispden,1+idir)=work(ifft) 271 end do 272 273 if (present(lrhonow)) then 274 call fourdp(cplex,worklp,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 275 do ifft=1,cplex*nfft 276 lrhonow(ifft,ispden)=lrhonow(ifft,ispden)+work(ifft) 277 end do 278 end if 279 280 end do 281 ABI_DEALLOCATE(gcart1) 282 ABI_DEALLOCATE(gcart2) 283 ABI_DEALLOCATE(gcart3) 284 ABI_DEALLOCATE(workgr) 285 if (allocated(worklp)) then 286 ABI_DEALLOCATE(worklp) 287 end if 288 end if 289 290 end do ! End loop on spins 291 292 ABI_DEALLOCATE(wkcmpx) 293 ABI_DEALLOCATE(work) 294 if(ishift==1) then 295 ABI_DEALLOCATE(ph1) 296 ABI_DEALLOCATE(ph2) 297 ABI_DEALLOCATE(ph3) 298 end if 299 300 end if ! End condition on ishift and ngrad 301 302 end subroutine xcden