TABLE OF CONTENTS
ABINIT/m_xctk [ Modules ]
NAME
m_xctk
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_xctk 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 28 use defs_abitypes, only : MPI_type 29 use m_time, only : timab 30 use m_mpinfo, only : ptabs_fourdp 31 use m_fft_mesh, only : phase 32 use m_fft, only : fourdp 33 34 implicit none 35 36 private
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.
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)
SOURCE
84 subroutine xcden (cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,qphon,rhor,rhonow, & !Mandatory arguments 85 & lrhonow) !Optional arguments 86 87 !Arguments ------------------------------------ 88 !scalars 89 integer,intent(in) :: cplex,ishift,nfft,ngrad,nspden 90 type(MPI_type),intent(in) :: mpi_enreg 91 !arrays 92 integer,intent(in) :: ngfft(18) 93 real(dp),intent(in) :: gprimd(3,3),qphon(3),rhor(cplex*nfft,nspden) 94 real(dp),intent(out) :: rhonow(cplex*nfft,nspden,ngrad*ngrad) 95 real(dp),intent(out),optional :: lrhonow(cplex*nfft,nspden) 96 97 !Local variables------------------------------- 98 !scalars 99 integer :: i1,i2,i3,id1,id2,id3,idir,ifft,ig1,ig2,ig3 100 integer :: ispden,n1,n2,n3,qeq0 101 real(dp) :: gc23_idir,gcart_idir,ph123i,ph123r,ph1i,ph1r,ph23i,ph23r,ph2i,ph2r 102 real(dp) :: ph3i,ph3r,work_im,work_re 103 character(len=500) :: message 104 !arrays 105 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 106 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 107 real(dp) :: tsec(2) 108 real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:),ph1(:),ph2(:),ph3(:) 109 real(dp),allocatable :: wkcmpx(:,:),work(:),workgr(:,:),worklp(:,:) 110 111 ! ************************************************************************* 112 113 !DEBUG 114 !write(std_out,*)' xcden : enter ' 115 !ENDDEBUG 116 117 118 if (ishift/=0 .and. ishift/=1) then 119 write(message, '(a,i0)' )'ishift must be 0 or 1 ; input was',ishift 120 ABI_BUG(message) 121 end if 122 123 if (ngrad/=1 .and. ngrad/=2) then 124 write(message, '(a,i0)' )'ngrad must be 1 or 2 ; input was',ngrad 125 ABI_BUG(message) 126 end if 127 128 !Keep local copy of fft dimensions 129 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 130 131 !Initialize computation of G in cartesian coordinates 132 id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2 133 134 !Get the distrib associated with this fft_grid 135 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 136 137 !Check whether q=0 138 qeq0=0 139 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1 140 141 if(ishift==0)then 142 143 ! Copy the input rhor in the new location. Will check on negative values later 144 145 do ispden=1,nspden 146 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,nfft,rhonow,rhor) 147 do ifft=1,cplex*nfft 148 rhonow(ifft,ispden,1)=rhor(ifft,ispden) 149 end do 150 end do 151 152 end if 153 154 if(ishift==1 .or. ngrad==2)then 155 156 ABI_MALLOC(wkcmpx,(2,nfft)) 157 ABI_MALLOC(work,(cplex*nfft)) 158 159 if(ishift==1)then 160 ! Precompute phases (The phases correspond to a shift of density on real space 161 ! grid from center at 0 0 0 to (1/2)*(1/n1,1/n2,1/n3).) 162 ABI_MALLOC(ph1,(2*n1)) 163 ABI_MALLOC(ph2,(2*n2)) 164 ABI_MALLOC(ph3,(2*n3)) 165 call phase(n1,ph1) 166 call phase(n2,ph2) 167 call phase(n3,ph3) 168 end if 169 170 do ispden=1,nspden 171 172 ! Obtain rho(G) in wkcmpx from input rho(r) 173 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,nfft,rhor,work) 174 do ifft=1,cplex*nfft 175 work(ifft)=rhor(ifft,ispden) 176 end do 177 178 call timab(82,1,tsec) 179 call fourdp(cplex,wkcmpx,work,-1,mpi_enreg,nfft,1,ngfft,0) 180 call timab(82,2,tsec) 181 182 ! If shift is required, multiply now rho(G) by phase, then generate rho(r+delta) 183 if(ishift==1)then 184 !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,ph1i,ph1r,ph123i,ph123r,ph2i,ph2r,ph23i,ph23r,ph3i,ph3r,work_im,work_re) & 185 !$OMP&SHARED(n1,n2,n3,ph1,ph2,ph3,wkcmpx,mpi_enreg,fftn2_distrib) 186 do i3=1,n3 187 ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft) 188 ph3r=ph3(2*i3-1) 189 ph3i=ph3(2*i3 ) 190 do i2=1,n2 191 ph2r=ph2(2*i2-1) 192 ph2i=ph2(2*i2 ) 193 ph23r=ph2r*ph3r-ph2i*ph3i 194 ph23i=ph2i*ph3r+ph2r*ph3i 195 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 196 do i1=1,n1 197 ifft=ifft+1 198 ph1r=ph1(2*i1-1) 199 ph1i=ph1(2*i1 ) 200 ph123r=ph1r*ph23r-ph1i*ph23i 201 ph123i=ph1i*ph23r+ph1r*ph23i 202 ! Must use intermediate variables ! 203 work_re=ph123r*wkcmpx(1,ifft)-ph123i*wkcmpx(2,ifft) 204 work_im=ph123i*wkcmpx(1,ifft)+ph123r*wkcmpx(2,ifft) 205 wkcmpx(1,ifft)=work_re 206 wkcmpx(2,ifft)=work_im 207 end do 208 end if 209 end do 210 end do 211 call timab(82,1,tsec) 212 call fourdp(cplex,wkcmpx,work,1,mpi_enreg,nfft,1,ngfft,0) 213 call timab(82,2,tsec) 214 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,rhonow,work) 215 do ifft=1,cplex*nfft 216 rhonow(ifft,ispden,1)=work(ifft) 217 end do 218 end if 219 220 ! If gradient of the density is required, take care of the three components now 221 ! Note : this operation is applied on the eventually shifted rho(G) 222 if(ngrad==2)then 223 ABI_MALLOC(gcart1,(n1)) 224 ABI_MALLOC(gcart2,(n2)) 225 ABI_MALLOC(gcart3,(n3)) 226 ABI_MALLOC(workgr,(2,nfft)) 227 if (present(lrhonow)) then 228 ABI_MALLOC(worklp,(2,nfft)) 229 lrhonow(:,ispden)=zero 230 end if 231 do idir=1,3 232 233 do i1=1,n1 234 ig1=i1-(i1/id1)*n1-1 235 gcart1(i1)=gprimd(idir,1)*two_pi*(dble(ig1)+qphon(1)) 236 end do 237 ! Note that the G <-> -G symmetry must be maintained 238 if(mod(n1,2)==0 .and. qeq0==1)gcart1(n1/2+1)=zero 239 do i2=1,n2 240 ig2=i2-(i2/id2)*n2-1 241 gcart2(i2)=gprimd(idir,2)*two_pi*(dble(ig2)+qphon(2)) 242 end do 243 if(mod(n2,2)==0 .and. qeq0==1)gcart2(n2/2+1)=zero 244 do i3=1,n3 245 ig3=i3-(i3/id3)*n3-1 246 gcart3(i3)=gprimd(idir,3)*two_pi*(dble(ig3)+qphon(3)) 247 end do 248 if(mod(n3,2)==0 .and. qeq0==1)gcart3(n3/2+1)=zero 249 250 ! MG: Be careful here with OMP due to ifft. Disabled for the time being. 251 ! !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,gcart_idir,gc23_idir) & 252 ! !$OMP&SHARED(gcart1,gcart2,gcart3,n1,n2,n3,wkcmpx,workgr) 253 ifft = 0 254 do i3=1,n3 255 do i2=1,n2 256 gc23_idir=gcart2(i2)+gcart3(i3) 257 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 258 do i1=1,n1 259 ifft=ifft+1 260 gcart_idir=gc23_idir+gcart1(i1) 261 ! Multiply by i 2pi G(idir) 262 workgr(2,ifft)= gcart_idir*wkcmpx(1,ifft) 263 workgr(1,ifft)=-gcart_idir*wkcmpx(2,ifft) 264 ! Do the same to the gradient in order to get the laplacian 265 if (present(lrhonow)) then 266 worklp(2,ifft)= gcart_idir*workgr(1,ifft) 267 worklp(1,ifft)=-gcart_idir*workgr(2,ifft) 268 end if 269 end do 270 end if 271 end do 272 end do 273 call timab(82,1,tsec) 274 call fourdp(cplex,workgr,work,1,mpi_enreg,nfft,1,ngfft,0) 275 call timab(82,2,tsec) 276 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(idir,ispden,cplex,nfft,rhonow,work) 277 do ifft=1,cplex*nfft 278 rhonow(ifft,ispden,1+idir)=work(ifft) 279 end do 280 281 if (present(lrhonow)) then 282 call fourdp(cplex,worklp,work,1,mpi_enreg,nfft,1,ngfft,0) 283 do ifft=1,cplex*nfft 284 lrhonow(ifft,ispden)=lrhonow(ifft,ispden)+work(ifft) 285 end do 286 end if 287 288 end do 289 ABI_FREE(gcart1) 290 ABI_FREE(gcart2) 291 ABI_FREE(gcart3) 292 ABI_FREE(workgr) 293 if (allocated(worklp)) then 294 ABI_FREE(worklp) 295 end if 296 end if 297 298 end do ! End loop on spins 299 if (allocated(wkcmpx)) then 300 ABI_FREE(wkcmpx) 301 end if 302 ABI_FREE(work) 303 if(ishift==1) then 304 ABI_FREE(ph1) 305 ABI_FREE(ph2) 306 ABI_FREE(ph3) 307 end if 308 309 end if ! End condition on ishift and ngrad 310 311 end subroutine xcden
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
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 use_laplacian=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 use_laplacian=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 use_laplacian : 1 if we use a functional depending on the laplacian of the density 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 (all optional: [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).
SOURCE
379 subroutine xcpot (cplex,gprimd,ishift,use_laplacian,mpi_enreg,nfft,ngfft,ngrad,nspden,& 380 & nspgrad,qphon,& 381 & depsxc,rhonow,vxc,vxctau) ! optional argument 382 383 !Arguments ------------------------------------ 384 !scalars 385 integer,intent(in) :: cplex,ishift,nfft,ngrad,nspden,nspgrad,use_laplacian 386 type(MPI_type),intent(in) :: mpi_enreg 387 !arrays 388 integer,intent(in) :: ngfft(18) 389 real(dp),intent(in),optional :: rhonow(cplex*nfft,nspden,ngrad*ngrad) 390 real(dp),intent(in),optional :: depsxc(cplex*nfft,nspgrad),gprimd(3,3),qphon(3) 391 real(dp),intent(inout),optional :: vxc(cplex*nfft,nspden) 392 real(dp),intent(inout),optional :: vxctau(cplex*nfft,nspden,4) 393 394 !Local variables------------------------------- 395 !scalars 396 integer :: i1,i2,i3,id1,id2,id3,idir,ifft,ig1,ig2,ig3,ispden,n1,n2,n3,qeq0 397 real(dp),parameter :: lowden=1.d-14,precis=1.d-15 398 real(dp) :: gc23_idir,gcart_idir,ph123i,ph123r,ph1i,ph1r,ph23i,ph23r,ph2i,ph2r 399 real(dp) :: ph3i,ph3r,work_im,work_re 400 character(len=500) :: message 401 !arrays 402 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 403 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 404 logical :: with_vxc,with_vxctau 405 real(dp) :: tsec(2) 406 real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:),ph1(:),ph2(:),ph3(:) 407 real(dp),allocatable :: wkcmpx(:,:),wkcmpxtau(:,:) 408 real(dp),allocatable :: work(:),workgr(:,:),worklp(:,:),worktau(:,:) 409 410 ! ************************************************************************* 411 412 if (ishift/=0 .and. ishift/=1) then 413 write(message, '(a,i0)' )' ishift must be 0 or 1 ; input was',ishift 414 ABI_BUG(message) 415 end if 416 417 if (ngrad/=1 .and. ngrad/=2 ) then 418 write(message, '(a,i0)' )' ngrad must be 1 or 2 ; input was',ngrad 419 ABI_BUG(message) 420 end if 421 422 with_vxc=present(vxc) ; with_vxctau=present(vxctau) 423 if (with_vxc) then 424 if ((.not.present(rhonow)).or.(.not.present(depsxc))) then 425 message='need rhonow or depsxc!' 426 ABI_BUG(message) 427 end if 428 end if 429 430 !Keep local copy of fft dimensions 431 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 432 433 !Initialize computation of G in cartesian coordinates 434 id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2 435 436 !Get the distrib associated with this fft_grid 437 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 438 439 !Check whether q=0 440 qeq0=0;if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1 441 442 if(with_vxc.and.ishift==0)then ! Add the value of depsxc to vxc 443 do ispden=1,min(nspden,2) 444 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,depsxc,nfft,vxc,ispden) 445 do ifft=1,cplex*nfft 446 vxc(ifft,ispden)=vxc(ifft,ispden)+depsxc(ifft,ispden) 447 end do 448 end do 449 end if 450 451 !If the grid is shifted, or if gradient corrections are present, there must be FFTs. 452 if(ishift==1 .or. ngrad==2)then 453 454 if(with_vxc.or.with_vxctau) then 455 ABI_MALLOC(work,(cplex*nfft)) 456 end if 457 if (with_vxc) then 458 ABI_MALLOC(wkcmpx,(2,nfft)) 459 end if 460 461 if(ishift==1)then 462 ABI_MALLOC(ph1,(2*n1)) 463 ABI_MALLOC(ph2,(2*n2)) 464 ABI_MALLOC(ph3,(2*n3)) 465 ! Precompute phases (The phases correspond to a shift of density on real space 466 ! grid from center at 0 0 0 to (1/2)*(1/n1,1/n2,1/n3).) 467 call phase(n1,ph1) 468 call phase(n2,ph2) 469 call phase(n3,ph3) 470 end if 471 472 do ispden=1,min(nspden,2) 473 474 ! Initialize wkcmpx either to 0 or to the shifted vxc value 475 if (with_vxc) then 476 if(ishift==0)then 477 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(nfft,wkcmpx) 478 do ifft=1,nfft 479 wkcmpx(:,ifft)=zero 480 end do 481 else 482 ! Obtain depsxc(G)*phase in wkcmpx from input depsxc(r+delta) 483 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,depsxc,ispden,nfft,work) 484 do ifft=1,cplex*nfft 485 work(ifft)=depsxc(ifft,ispden) 486 end do 487 call timab(82,1,tsec) 488 call fourdp(cplex,wkcmpx,work,-1,mpi_enreg,nfft,1,ngfft,0) 489 call timab(82,2,tsec) 490 end if 491 end if 492 493 ! If gradient correction is present, take care of the three components now 494 ! Note : this operation is done on the eventually shifted grid 495 if (ngrad==2) then 496 ABI_MALLOC(gcart1,(n1)) 497 ABI_MALLOC(gcart2,(n2)) 498 ABI_MALLOC(gcart3,(n3)) 499 if (with_vxc) then 500 ABI_MALLOC(workgr,(2,nfft)) 501 if (use_laplacian==1) then 502 ABI_MALLOC(worklp,(2,nfft)) 503 end if 504 end if 505 if (with_vxctau) then 506 ABI_MALLOC(worktau,(2,nfft)) 507 ABI_MALLOC(wkcmpxtau,(2,nfft)) 508 end if 509 510 do idir=1,3 511 512 if (with_vxc) then 513 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,idir,ispden,nfft,rhonow,work) 514 do ifft=1,cplex*nfft 515 work(ifft)=rhonow(ifft,ispden,1+idir) 516 end do 517 call timab(82,1,tsec) 518 call fourdp(cplex,workgr,work,-1,mpi_enreg,nfft,1,ngfft,0) 519 call timab(82,2,tsec) 520 521 ! IF Meta-GGA then take care of the laplacian term involved. 522 ! Note : this operation is done on the eventually shifted grid 523 if(use_laplacian==1)then 524 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,idir,ispden,nspden,nfft,depsxc,work) 525 do ifft=1,cplex*nfft 526 if(nspden==1)then 527 work(ifft)=depsxc(ifft,2+ispden) 528 else if(nspden==2)then 529 work(ifft)=depsxc(ifft,5+ispden) 530 end if 531 end do 532 call timab(82,1,tsec) 533 call fourdp(cplex,worklp,work,-1,mpi_enreg,nfft,1,ngfft,0) 534 call timab(82,2,tsec) 535 end if 536 end if 537 538 if(with_vxctau)then 539 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxctau,work) 540 do ifft=1,cplex*nfft 541 work(ifft)=vxctau(ifft,ispden,1) 542 end do 543 call timab(82,1,tsec) 544 call fourdp(cplex,worktau,work,-1,mpi_enreg,nfft,1,ngfft,0) 545 call timab(82,2,tsec) 546 end if ! present vxctau 547 548 do i1=1,n1 549 ig1=i1-(i1/id1)*n1-1 550 gcart1(i1)=gprimd(idir,1)*two_pi*(dble(ig1)+qphon(1)) 551 end do 552 ! Note that the G <-> -G symmetry must be maintained 553 if(mod(n1,2)==0 .and. qeq0==1)gcart1(n1/2+1)=zero 554 do i2=1,n2 555 ig2=i2-(i2/id2)*n2-1 556 gcart2(i2)=gprimd(idir,2)*two_pi*(dble(ig2)+qphon(2)) 557 end do 558 if(mod(n2,2)==0 .and. qeq0==1)gcart2(n2/2+1)=zero 559 do i3=1,n3 560 ig3=i3-(i3/id3)*n3-1 561 gcart3(i3)=gprimd(idir,3)*two_pi*(dble(ig3)+qphon(3)) 562 end do 563 if(mod(n3,2)==0 .and. qeq0==1)gcart3(n3/2+1)=zero 564 565 ! !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,gc23_idir,gcart_idir) & 566 ! !$OMP&SHARED(gcart1,gcart2,gcart3,n1,n2,n3,wkcmpx,workgr) 567 ifft = 0 568 do i3=1,n3 569 do i2=1,n2 570 gc23_idir=gcart2(i2)+gcart3(i3) 571 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 572 do i1=1,n1 573 ifft=ifft+1 574 gcart_idir=gc23_idir+gcart1(i1) 575 if(with_vxc)then 576 ! Multiply by - i 2pi G(idir) and accumulate in wkcmpx 577 wkcmpx(1,ifft)=wkcmpx(1,ifft)+gcart_idir*workgr(2,ifft) 578 wkcmpx(2,ifft)=wkcmpx(2,ifft)-gcart_idir*workgr(1,ifft) 579 if(use_laplacian==1)then 580 ! Multiply by - i 2pi G(idir) and accumulate in wkcmpx 581 wkcmpx(1,ifft)=wkcmpx(1,ifft)-gcart_idir**2*worklp(1,ifft) 582 wkcmpx(2,ifft)=wkcmpx(2,ifft)-gcart_idir**2*worklp(2,ifft) 583 end if 584 end if 585 if(with_vxctau)then 586 wkcmpxtau(1,ifft)= gcart_idir*worktau(2,ifft) 587 wkcmpxtau(2,ifft)=-gcart_idir*worktau(1,ifft) 588 end if 589 end do 590 end if 591 end do 592 end do 593 594 if (with_vxctau) then 595 call timab(82,1,tsec) 596 call fourdp(cplex,wkcmpxtau,work,1,mpi_enreg,nfft,1,ngfft,0) 597 call timab(82,2,tsec) 598 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxctau,work) 599 do ifft=1,cplex*nfft 600 vxctau(ifft,ispden,1+idir)=work(ifft) 601 end do 602 end if 603 604 end do ! enddo idir 605 606 ABI_FREE(gcart1) 607 ABI_FREE(gcart2) 608 ABI_FREE(gcart3) 609 if (with_vxc) then 610 ABI_FREE(workgr) 611 if (use_laplacian==1) then 612 ABI_FREE(worklp) 613 end if 614 end if 615 if (with_vxctau) then 616 ABI_FREE(worktau) 617 ABI_FREE(wkcmpxtau) 618 end if 619 620 end if 621 622 ! wkcmpx(:,:) contains now the full exchange-correlation potential, but 623 ! eventually for the shifted grid 624 625 if (with_vxc) then 626 if(ishift==1)then 627 ! Take away the phase to get depsxc(G) 628 ifft=0 629 do i3=1,n3 630 ph3r=ph3(2*i3-1) 631 ph3i=ph3(2*i3 ) 632 do i2=1,n2 633 ph2r=ph2(2*i2-1) 634 ph2i=ph2(2*i2 ) 635 ph23r=ph2r*ph3r-ph2i*ph3i 636 ph23i=ph2i*ph3r+ph2r*ph3i 637 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 638 do i1=1,n1 639 ifft=ifft+1 640 ph1r=ph1(2*i1-1) 641 ph1i=ph1(2*i1 ) 642 ph123r=ph1r*ph23r-ph1i*ph23i 643 ph123i=ph1i*ph23r+ph1r*ph23i 644 ! Multiply by phase. Must use intermediate variables ! 645 work_re= ph123r*wkcmpx(1,ifft)+ph123i*wkcmpx(2,ifft) 646 work_im=-ph123i*wkcmpx(1,ifft)+ph123r*wkcmpx(2,ifft) 647 wkcmpx(1,ifft)=work_re 648 wkcmpx(2,ifft)=work_im 649 end do 650 end if 651 end do 652 end do 653 end if 654 655 call timab(82,1,tsec) 656 call fourdp(cplex,wkcmpx,work,1,mpi_enreg,nfft,1,ngfft,0) 657 call timab(82,2,tsec) 658 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxc,work) 659 do ifft=1,cplex*nfft 660 vxc(ifft,ispden)=vxc(ifft,ispden)+work(ifft) 661 end do 662 end if 663 664 end do ! End loop on spins 665 666 if(ishift==1) then 667 ABI_FREE(ph1) 668 ABI_FREE(ph2) 669 ABI_FREE(ph3) 670 end if 671 if(with_vxc) then 672 ABI_FREE(wkcmpx) 673 end if 674 if(with_vxc.or.with_vxctau) then 675 ABI_FREE(work) 676 end if 677 678 end if ! End condition on ishift/ngrad 679 680 end subroutine xcpot
ABINIT/xcpotdq [ Functions ]
NAME
xcpotdq
FUNCTION
Equivalent to xcpot for the q-derivative of the GGA xc kernel.
INPUTS
agradn(cplex*nfft,nspgrad,3)=kxc(:,4)*gradrho(:,:)*gradrho(:,qdir)*rho1(*) 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 (not implemented) 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
OUTPUT
vxc(cplex*nfft,nspden)]=q-derivative of the GGA xc potential.
At input already includes three terms.
SOURCE
709 subroutine xcpotdq (agradn,cplex,gprimd,ishift,mpi_enreg, & 710 & nfft,ngfft,ngrad,nspden,nspgrad,vxc) 711 712 !Arguments ------------------------------------ 713 !scalars 714 integer,intent(in) :: cplex,ishift,nfft,ngrad,nspden,nspgrad 715 type(MPI_type),intent(in) :: mpi_enreg 716 !arrays 717 integer,intent(in) :: ngfft(18) 718 real(dp),intent(in) :: agradn(cplex*nfft,nspgrad,3) 719 real(dp),intent(in) :: gprimd(3,3) 720 real(dp),intent(inout) :: vxc(2*nfft,nspden) 721 722 !Local variables------------------------------- 723 !scalars 724 integer :: i1,i2,i3,id1,id2,id3,idir,ifft,ig1,ig2,ig3,ispden,n1,n2,n3 725 real(dp),parameter :: lowden=1.d-14,precis=1.d-15 726 real(dp) :: gc23_idir,gcart_idir 727 character(len=500) :: message 728 !arrays 729 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 730 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 731 real(dp) :: tsec(2) 732 real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:) 733 real(dp),allocatable :: wkcmpx(:,:) 734 real(dp),allocatable :: work(:),workgr(:,:) 735 736 ! ************************************************************************* 737 738 if (ishift/=0) then 739 write(message, '(a,i0)' )' ishift must be 0 ; input was',ishift 740 ABI_BUG(message) 741 end if 742 743 if (ngrad/=2) then 744 write(message, '(a,i0)' )' ngrad must be 2 ; input was',ngrad 745 ABI_BUG(message) 746 end if 747 748 !Keep local copy of fft dimensions 749 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 750 751 !Initialize computation of G in cartesian coordinates 752 id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2 753 754 !Get the distrib associated with this fft_grid 755 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 756 757 !Compute the real-space gradient of de second term 758 ABI_MALLOC(work,(cplex*nfft)) 759 ABI_MALLOC(wkcmpx,(2,nfft)) 760 ABI_MALLOC(workgr,(2,nfft)) 761 762 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(nfft,wkcmpx) 763 do ifft=1,nfft 764 wkcmpx(:,ifft)=zero 765 end do 766 767 ! Obtain agradn(G)*phase in wkcmpx from input agradn(r) 768 ispden=1 769 ABI_MALLOC(gcart1,(n1)) 770 ABI_MALLOC(gcart2,(n2)) 771 ABI_MALLOC(gcart3,(n3)) 772 do idir=1, 3 773 774 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,idir,agradn,ispden,nfft,work) 775 do ifft=1,cplex*nfft 776 work(ifft)=agradn(ifft,ispden,idir) 777 end do 778 call timab(82,1,tsec) 779 call fourdp(cplex,workgr,work,-1,mpi_enreg,nfft,1,ngfft,0) 780 call timab(82,2,tsec) 781 782 do i1=1,n1 783 ig1=i1-(i1/id1)*n1-1 784 gcart1(i1)=gprimd(idir,1)*two_pi*dble(ig1) 785 end do 786 !Note that the G <-> -G symmetry must be maintained 787 if(mod(n1,2)==0) gcart1(n1/2+1)=zero 788 do i2=1,n2 789 ig2=i2-(i2/id2)*n2-1 790 gcart2(i2)=gprimd(idir,2)*two_pi*dble(ig2) 791 end do 792 if(mod(n2,2)==0) gcart2(n2/2+1)=zero 793 do i3=1,n3 794 ig3=i3-(i3/id3)*n3-1 795 gcart3(i3)=gprimd(idir,3)*two_pi*dble(ig3) 796 end do 797 if(mod(n3,2)==0) gcart3(n3/2+1)=zero 798 799 ! !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,gc23_idir,gcart_idir) & 800 ! !$OMP&SHARED(gcart1,gcart2,gcart3,n1,n2,n3,wkcmpx,workgr) 801 ifft = 0 802 do i3=1,n3 803 do i2=1,n2 804 gc23_idir=gcart2(i2)+gcart3(i3) 805 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 806 do i1=1,n1 807 ifft=ifft+1 808 gcart_idir=gc23_idir+gcart1(i1) 809 ! Multiply by -i 2pi G(idir) and accumulate in wkcmpx 810 wkcmpx(1,ifft)=wkcmpx(1,ifft)+gcart_idir*workgr(2,ifft) 811 wkcmpx(2,ifft)=wkcmpx(2,ifft)-gcart_idir*workgr(1,ifft) 812 end do 813 end if 814 end do 815 end do 816 817 end do 818 819 ABI_FREE(gcart1) 820 ABI_FREE(gcart2) 821 ABI_FREE(gcart3) 822 ABI_FREE(workgr) 823 824 call timab(82,1,tsec) 825 call fourdp(cplex,wkcmpx,work,1,mpi_enreg,nfft,1,ngfft,0) 826 call timab(82,2,tsec) 827 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vxc,work) 828 do ifft=1,nfft 829 vxc(2*ifft,ispden)=vxc(2*ifft,ispden)+work(ifft) 830 !Apply here the two pi factor 831 vxc(2*ifft,ispden)=vxc(2*ifft,ispden)*two_pi 832 end do 833 ABI_FREE(wkcmpx) 834 ABI_FREE(work) 835 836 end subroutine xcpotdq