TABLE OF CONTENTS
ABINIT/xcpot [ Functions ]
NAME
xcpot
FUNCTION
Process data after calling local or semi-local xc evaluators The derivative of Exc with respect to the density, gradient of density in case of GGAs, and Laplacian of density in case of meta-GGA are contained in depsxc(:,:). In case of GGAs (and meta-GGAs) the gradient of the density contained in rhonow(:,:,2:4) is already multiplied by the local partial derivative of the XC functional. Can take into account a shift of the grid, for purpose of better accuracy
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 depsxc(cplex*nfft,nspgrad)=derivative of Exc with respect to the (spin-)density, or to the norm of the gradient of the (spin-)density, further divided by the norm of the gradient of the (spin-)density The different components of depsxc will be for nspden=1, depsxc(:,1)=d(rho.exc)/d(rho) and if ngrad=2, depsxc(:,2)=1/2*1/|grad rho_up|*d(n.exc)/d(|grad rho_up|) +1/|grad rho|*d(rho.exc)/d(|grad rho|) and if mgga=1, depsxc(:,3)=d(rho.exc)/d(lapl rho) for nspden>=2, depsxc(:,1)=d(rho.exc)/d(rho_up) depsxc(:,2)=d(rho.exc)/d(rho_down) and if ngrad=2, depsxc(:,3)=1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|) depsxc(:,4)=1/|grad rho_down|*d(rho.exc)/d(|grad rho_down|) depsxc(:,5)=1/|grad rho|*d(rho.exc)/d(|grad rho|) and if mgga=1, depsxc(:,6)=d(rho.exc)/d(lapl rho_up) depsxc(:,7)=d(rho.exc)/d(lapl rho_down) 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 mgga : 1 if we use a meta-GGA functional. 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 take into account derivative wrt the density ; =2, also take into account derivative wrt the gradient of the density. nspden=number of spin-density components nspgrad=number of spin-density and spin-density-gradient components qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2). rhonow(cplex*nfft,nspden,ngrad*ngrad)=electron (spin)-density in real space and eventually its gradient already multiplied by the local partial derivative of the XC functional, 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 el./bohr**4, times local partial derivative of the functional, as required by the GGA In this routine, rhonow is used only in the GGA case (ngrad=2).
OUTPUT
(see side effects)
SIDE EFFECTS
Input/Output: vxc(cplex*nfft,nspden)=xc potential (spin up in first half and spin down in second half if nspden>=2). Contribution from the present shifted or unshifted grid is ADDED to the input vxc data. vxctau(cplex*nfft,nspden,4)=derivative of XC energy density with respect to kinetic energy density (depsxcdtau). The arrays vxctau(nfft,nspden,4) contains also the gradient of vxctau (gvxctau) which will be computed here in vxctau(:,:,2:4).
PARENTS
dfpt_mkvxcgga,dfpt_mkvxcstrgga,rhotoxc
CHILDREN
fourdp,phase,ptabs_fourdp,timab
SOURCE
81 #if defined HAVE_CONFIG_H 82 #include "config.h" 83 #endif 84 85 #include "abi_common.h" 86 87 subroutine xcpot (cplex,depsxc,gprimd,ishift,mgga,mpi_enreg,nfft,ngfft,ngrad,nspden,& 88 & nspgrad,paral_kgb,qphon,rhonow,vxc,& 89 & vxctau) ! optional argument 90 91 use defs_basis 92 use defs_abitypes 93 use m_profiling_abi 94 use m_errors 95 96 use m_mpinfo, only : ptabs_fourdp 97 98 !This section has been created automatically by the script Abilint (TD). 99 !Do not modify the following lines by hand. 100 #undef ABI_FUNC 101 #define ABI_FUNC 'xcpot' 102 use interfaces_18_timing 103 use interfaces_53_ffts 104 use interfaces_56_xc, except_this_one => xcpot 105 !End of the abilint section 106 107 implicit none 108 109 !Arguments ------------------------------------ 110 !scalars 111 integer,intent(in) :: cplex,ishift,mgga,nfft,ngrad,nspden,nspgrad,paral_kgb 112 type(MPI_type),intent(in) :: mpi_enreg 113 !arrays 114 integer,intent(in) :: ngfft(18) 115 real(dp),intent(in) :: depsxc(cplex*nfft,nspgrad),gprimd(3,3),qphon(3) 116 real(dp),intent(in) :: rhonow(cplex*nfft,nspden,ngrad*ngrad) 117 real(dp),intent(inout) :: vxc(cplex*nfft,nspden) !vz_i 118 real(dp),intent(inout),optional :: vxctau(cplex*nfft,nspden,4) 119 120 !Local variables------------------------------- 121 !scalars 122 integer :: i1,i2,i3,id1,id2,id3,idir,ifft,ig1,ig2,ig3,ispden,n1,n2,n3,qeq0 123 real(dp),parameter :: lowden=1.d-14,precis=1.d-15 124 real(dp) :: gc23_idir,gcart_idir,ph123i,ph123r,ph1i,ph1r,ph23i,ph23r,ph2i,ph2r 125 real(dp) :: ph3i,ph3r,work_im,work_re 126 character(len=500) :: message 127 !arrays 128 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 129 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 130 real(dp) :: tsec(2) 131 real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:),ph1(:),ph2(:),ph3(:) 132 real(dp),allocatable :: wkcmpx(:,:),wkcmpxtau(:,:) 133 real(dp),allocatable :: work(:),workgr(:,:),worklp(:,:),worktau(:,:) 134 135 ! ************************************************************************* 136 137 if (ishift/=0 .and. ishift/=1) then 138 write(message, '(a,i0)' )' ishift must be 0 or 1 ; input was',ishift 139 MSG_BUG(message) 140 end if 141 142 if (ngrad/=1 .and. ngrad/=2 ) then 143 write(message, '(a,i0)' )' ngrad must be 1 or 2 ; input was',ngrad 144 MSG_BUG(message) 145 end if 146 147 !Keep local copy of fft dimensions 148 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 149 150 !Initialize computation of G in cartesian coordinates 151 id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2 152 153 !Get the distrib associated with this fft_grid 154 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 155 156 !Check whether q=0 157 qeq0=0 158 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1 159 160 if(ishift==0)then ! Add the value of depsxc to vxc 161 do ispden=1,min(nspden,2) 162 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,depsxc,nfft,vxc,ispden) 163 do ifft=1,cplex*nfft 164 vxc(ifft,ispden)=vxc(ifft,ispden)+depsxc(ifft,ispden) 165 end do 166 end do 167 168 end if 169 170 !If the grid is shifted, or if gradient corrections are present, there must be FFTs. 171 if(ishift==1 .or. ngrad==2)then 172 173 ABI_ALLOCATE( wkcmpx,(2,nfft)) 174 ABI_ALLOCATE(work,(cplex*nfft)) 175 176 if(ishift==1)then 177 ABI_ALLOCATE( ph1,(2*n1)) 178 ABI_ALLOCATE(ph2,(2*n2)) 179 ABI_ALLOCATE(ph3,(2*n3)) 180 ! Precompute phases (The phases correspond to a shift of density on real space 181 ! grid from center at 0 0 0 to (1/2)*(1/n1,1/n2,1/n3).) 182 call phase(n1,ph1) 183 call phase(n2,ph2) 184 call phase(n3,ph3) 185 end if 186 187 do ispden=1,min(nspden,2) 188 189 ! Initialize wkcmpx either to 0 or to the shifted vxc value 190 if(ishift==0)then 191 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(nfft,wkcmpx) 192 do ifft=1,nfft 193 wkcmpx(:,ifft)=zero 194 end do 195 196 else 197 ! Obtain depsxc(G)*phase in wkcmpx from input depsxc(r+delta) 198 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,depsxc,ispden,nfft,work) 199 do ifft=1,cplex*nfft 200 work(ifft)=depsxc(ifft,ispden) 201 end do 202 call timab(82,1,tsec) 203 call fourdp(cplex,wkcmpx,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 204 call timab(82,2,tsec) 205 end if 206 207 ! If gradient correction is present, take care of the three components now 208 ! Note : this operation is done on the eventually shifted grid 209 if(ngrad==2)then 210 ABI_ALLOCATE(gcart1,(n1)) 211 ABI_ALLOCATE(gcart2,(n2)) 212 ABI_ALLOCATE(gcart3,(n3)) 213 ABI_ALLOCATE(workgr,(2,nfft)) 214 if(mgga==1) then 215 ABI_ALLOCATE(worklp,(2,nfft)) 216 end if 217 if(present(vxctau)) then 218 ABI_ALLOCATE(worktau,(2,nfft)) 219 ABI_ALLOCATE(wkcmpxtau,(2,nfft)) 220 end if 221 222 do idir=1,3 223 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,idir,ispden,nfft,rhonow,work) 224 do ifft=1,cplex*nfft 225 work(ifft)=rhonow(ifft,ispden,1+idir) 226 end do 227 228 call timab(82,1,tsec) 229 call fourdp(cplex,workgr,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 230 call timab(82,2,tsec) 231 232 ! IF Meta-GGA then take care of the laplacian term involved. 233 ! Note : this operation is done on the eventually shifted grid 234 if(mgga==1)then 235 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,idir,ispden,nspden,nfft,depsxc,work) 236 do ifft=1,cplex*nfft 237 if(nspden==1)then 238 work(ifft)=depsxc(ifft,2+ispden) 239 else if(nspden==2)then 240 work(ifft)=depsxc(ifft,5+ispden) 241 end if 242 end do 243 244 call timab(82,1,tsec) 245 call fourdp(cplex,worklp,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 246 call timab(82,2,tsec) 247 end if 248 249 if(present(vxctau))then 250 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxctau,work) 251 do ifft=1,cplex*nfft 252 work(ifft)=vxctau(ifft,ispden,1) 253 end do 254 255 call timab(82,1,tsec) 256 call fourdp(cplex,worktau,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0) 257 call timab(82,2,tsec) 258 end if ! present vxctau 259 260 do i1=1,n1 261 ig1=i1-(i1/id1)*n1-1 262 gcart1(i1)=gprimd(idir,1)*two_pi*(dble(ig1)+qphon(1)) 263 end do 264 ! Note that the G <-> -G symmetry must be maintained 265 if(mod(n1,2)==0 .and. qeq0==1)gcart1(n1/2+1)=zero 266 do i2=1,n2 267 ig2=i2-(i2/id2)*n2-1 268 gcart2(i2)=gprimd(idir,2)*two_pi*(dble(ig2)+qphon(2)) 269 end do 270 if(mod(n2,2)==0 .and. qeq0==1)gcart2(n2/2+1)=zero 271 do i3=1,n3 272 ig3=i3-(i3/id3)*n3-1 273 gcart3(i3)=gprimd(idir,3)*two_pi*(dble(ig3)+qphon(3)) 274 end do 275 if(mod(n3,2)==0 .and. qeq0==1)gcart3(n3/2+1)=zero 276 277 ! !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,gc23_idir,gcart_idir) & 278 ! !$OMP&SHARED(gcart1,gcart2,gcart3,n1,n2,n3,wkcmpx,workgr) 279 ifft = 0 280 do i3=1,n3 281 do i2=1,n2 282 gc23_idir=gcart2(i2)+gcart3(i3) 283 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 284 do i1=1,n1 285 ifft=ifft+1 286 gcart_idir=gc23_idir+gcart1(i1) 287 ! Multiply by - i 2pi G(idir) and accumulate in wkcmpx 288 wkcmpx(1,ifft)=wkcmpx(1,ifft)+gcart_idir*workgr(2,ifft) 289 wkcmpx(2,ifft)=wkcmpx(2,ifft)-gcart_idir*workgr(1,ifft) 290 if(mgga==1)then 291 ! Multiply by - i 2pi G(idir) and accumulate in wkcmpx 292 wkcmpx(1,ifft)=wkcmpx(1,ifft)-gcart_idir**2*worklp(1,ifft) 293 wkcmpx(2,ifft)=wkcmpx(2,ifft)-gcart_idir**2*worklp(2,ifft) 294 end if 295 if(present(vxctau))then 296 wkcmpxtau(1,ifft)= gcart_idir*worktau(2,ifft) 297 wkcmpxtau(2,ifft)=-gcart_idir*worktau(1,ifft) 298 end if 299 end do 300 end if 301 end do 302 end do 303 304 if(present(vxctau))then 305 call timab(82,1,tsec) 306 call fourdp(cplex,wkcmpxtau,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 307 call timab(82,2,tsec) 308 309 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxctau,work) 310 do ifft=1,cplex*nfft 311 vxctau(ifft,ispden,1+idir)=work(ifft) 312 end do 313 end if ! present vxctau 314 315 end do ! enddo idir 316 ABI_DEALLOCATE(gcart1) 317 ABI_DEALLOCATE(gcart2) 318 ABI_DEALLOCATE(gcart3) 319 ABI_DEALLOCATE(workgr) 320 if(mgga==1)then 321 ABI_DEALLOCATE(worklp) 322 end if 323 if(present(vxctau)) then 324 ABI_DEALLOCATE(worktau) 325 ABI_DEALLOCATE(wkcmpxtau) 326 end if 327 end if 328 329 ! wkcmpx(:,:) contains now the full exchange-correlation potential, but 330 ! eventually for the shifted grid 331 332 if(ishift==1)then 333 ! Take away the phase to get depsxc(G) 334 ifft=0 335 do i3=1,n3 336 ph3r=ph3(2*i3-1) 337 ph3i=ph3(2*i3 ) 338 do i2=1,n2 339 ph2r=ph2(2*i2-1) 340 ph2i=ph2(2*i2 ) 341 ph23r=ph2r*ph3r-ph2i*ph3i 342 ph23i=ph2i*ph3r+ph2r*ph3i 343 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 344 do i1=1,n1 345 ifft=ifft+1 346 ph1r=ph1(2*i1-1) 347 ph1i=ph1(2*i1 ) 348 ph123r=ph1r*ph23r-ph1i*ph23i 349 ph123i=ph1i*ph23r+ph1r*ph23i 350 ! Multiply by phase. Must use intermediate variables ! 351 work_re= ph123r*wkcmpx(1,ifft)+ph123i*wkcmpx(2,ifft) 352 work_im=-ph123i*wkcmpx(1,ifft)+ph123r*wkcmpx(2,ifft) 353 wkcmpx(1,ifft)=work_re 354 wkcmpx(2,ifft)=work_im 355 end do 356 end if 357 end do 358 end do 359 end if 360 361 call timab(82,1,tsec) 362 call fourdp(cplex,wkcmpx,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 363 call timab(82,2,tsec) 364 365 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxc,work) 366 do ifft=1,cplex*nfft 367 vxc(ifft,ispden)=vxc(ifft,ispden)+work(ifft) 368 end do 369 370 end do ! End loop on spins 371 372 if(ishift==1) then 373 ABI_DEALLOCATE(ph1) 374 ABI_DEALLOCATE(ph2) 375 ABI_DEALLOCATE(ph3) 376 end if 377 ABI_DEALLOCATE(wkcmpx) 378 ABI_DEALLOCATE(work) 379 380 end if ! End condition on ishift 381 382 end subroutine xcpot