TABLE OF CONTENTS
ABINIT/dfpt_mkvxcgga [ Functions ]
NAME
dfpt_mkvxcgga
FUNCTION
Compute the first-order change of exchange-correlation potential in case of GGA functionals Use the exchange-correlation kernel.
COPYRIGHT
Copyright (C) 2001-2018 ABINIT group (XG, 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 .
INPUTS
cplex= if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX gmet(3,3)=metrix tensor in G space in Bohr**-2. gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1) gsqcut=cutoff value on G**2 for sphere inside fft box. kxc(nfft,nkxc)=exchange and correlation kernel (see below) mpi_enreg=informations about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT nhat1(cplex*nfft,2*nhat1dim)= -PAW only- 1st-order compensation density nhat1dim= -PAW only- 1 if nhat1 array is used ; 0 otherwise nhat1gr(cplex*nfft,2,3*nhat1grdim)= -PAW only- gradients of 1st-order compensation density nhat1grdim= -PAW only- 1 if nhat1gr array is used ; 0 otherwise nkxc=second dimension of the kxc array nspden=number of spin-density components qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2). rhor1tmp(cplex*nfft,2)=array for first-order electron spin-density in electrons/bohr**3 (second index corresponds to spin-up and spin-down) usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
OUTPUT
vxc1(cplex*nfft,nspden)=change in exchange-correlation potential
SIDE EFFECTS
NOTES
For the time being, a rather crude coding, to be optimized ...
Content of Kxc array:
===== if GGA
if nspden==1:
kxc(:,1)= d2Exc/drho2
kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
kxc(:,5)= gradx(rho)
kxc(:,6)= grady(rho)
kxc(:,7)= gradz(rho)
if nspden>=2:
kxc(:,1)= d2Exc/drho_up drho_up
kxc(:,2)= d2Exc/drho_up drho_dn
kxc(:,3)= d2Exc/drho_dn drho_dn
kxc(:,4)= 1/|grad(rho_up)| dEx/d|grad(rho_up)|
kxc(:,5)= 1/|grad(rho_dn)| dEx/d|grad(rho_dn)|
kxc(:,6)= 1/|grad(rho_up)| d2Ex/d|grad(rho_up)| drho_up
kxc(:,7)= 1/|grad(rho_dn)| d2Ex/d|grad(rho_dn)| drho_dn
kxc(:,8)= 1/|grad(rho_up)| * d/d|grad(rho_up)| ( 1/|grad(rho_up)| dEx/d|grad(rho_up)| )
kxc(:,9)= 1/|grad(rho_dn)| * d/d|grad(rho_dn)| ( 1/|grad(rho_dn)| dEx/d|grad(rho_dn)| )
kxc(:,10)=1/|grad(rho)| dEc/d|grad(rho)|
kxc(:,11)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_up
kxc(:,12)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_dn
kxc(:,13)=1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dEc/d|grad(rho)| )
kxc(:,14)=gradx(rho_up)
kxc(:,15)=gradx(rho_dn)
kxc(:,16)=grady(rho_up)
kxc(:,17)=grady(rho_dn)
kxc(:,18)=gradz(rho_up)
kxc(:,19)=gradz(rho_dn)
PARENTS
dfpt_mkvxc
CHILDREN
xcden,xcpot
SOURCE
84 #if defined HAVE_CONFIG_H 85 #include "config.h" 86 #endif 87 88 #include "abi_common.h" 89 90 subroutine dfpt_mkvxcgga(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,& 91 & nhat1,nhat1dim,nhat1gr,nhat1grdim,nkxc,& 92 & nspden,paral_kgb,qphon,rhor1,usexcnhat,vxc1) 93 94 use defs_basis 95 use defs_abitypes 96 use m_profiling_abi 97 use m_errors 98 99 !This section has been created automatically by the script Abilint (TD). 100 !Do not modify the following lines by hand. 101 #undef ABI_FUNC 102 #define ABI_FUNC 'dfpt_mkvxcgga' 103 use interfaces_56_xc, except_this_one => dfpt_mkvxcgga 104 !End of the abilint section 105 106 implicit none 107 108 !Arguments ------------------------------------ 109 !scalars 110 integer,intent(in) :: cplex,nfft,nhat1dim,nhat1grdim,nkxc,nspden,paral_kgb,usexcnhat 111 type(MPI_type),intent(in) :: mpi_enreg 112 !arrays 113 integer,intent(in) :: ngfft(18) 114 real(dp),intent(in) :: gprimd(3,3) 115 real(dp),intent(in) :: kxc(nfft,nkxc) 116 real(dp),intent(in) :: nhat1(cplex*nfft,nspden*nhat1dim) 117 real(dp),intent(in) :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim) 118 real(dp),intent(in) :: qphon(3) 119 real(dp),intent(in),target :: rhor1(cplex*nfft,nspden) 120 real(dp),intent(out) :: vxc1(cplex*nfft,nspden) 121 122 !Local variables------------------------------- 123 !scalars 124 integer :: ii,ir,ishift,mgga,ngrad,nspgrad 125 logical :: test_nhat 126 real(dp) :: coeff_grho,coeff_grho_corr,coeff_grho_dn,coeff_grho_up 127 real(dp) :: coeffim_grho,coeffim_grho_corr,coeffim_grho_dn,coeffim_grho_up 128 real(dp) :: gradrho_gradrho1,gradrho_gradrho1_dn,gradrho_gradrho1_up 129 real(dp) :: gradrho_gradrho1im,gradrho_gradrho1im_dn,gradrho_gradrho1im_up 130 character(len=500) :: msg 131 !arrays 132 real(dp) :: r0(3),r0_dn(3),r0_up(3),r1(3),r1_dn(3),r1_up(3) 133 real(dp) :: r1im(3),r1im_dn(3),r1im_up(3) 134 real(dp),allocatable :: dnexcdn(:,:),rho1now(:,:,:) 135 real(dp),ABI_CONTIGUOUS pointer :: rhor1_ptr(:,:) 136 137 ! ************************************************************************* 138 139 DBG_ENTER("COLL") 140 141 if (nkxc/=12*min(nspden,2)-5) then 142 msg='Wrong nkxc value for GGA!' 143 MSG_BUG(msg) 144 end if 145 146 !metaGGA contributions are not taken into account here 147 mgga=0 148 149 !PAW: substract 1st-order compensation density from 1st-order density 150 test_nhat=((nhat1dim==1).and.(usexcnhat==0.or.nhat1grdim==1)) 151 if (test_nhat) then 152 ABI_ALLOCATE(rhor1_ptr,(cplex*nfft,nspden)) 153 rhor1_ptr(:,:)=rhor1(:,:)-nhat1(:,:) 154 else 155 rhor1_ptr => rhor1 156 end if 157 158 !call filterpot(paral_kgb,cplex,gmet,gsqcut,nfft,ngfft,2,qphon,rhor1_ptr) 159 160 !Compute the gradients of the first-order density 161 !rho1now(:,:,1) contains the first-order density, and 162 !rho1now(:,:,2:4) contains the gradients of the first-order density 163 ishift=0 ; ngrad=2 164 ABI_ALLOCATE(rho1now,(cplex*nfft,nspden,ngrad*ngrad)) 165 call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,paral_kgb,qphon,rhor1_ptr,rho1now) 166 167 !PAW: add "exact" gradients of compensation density 168 if (test_nhat.and.usexcnhat==1) then 169 rho1now(:,1:nspden,1)=rho1now(:,1:nspden,1)+nhat1(:,1:nspden) 170 end if 171 if (nhat1grdim==1) then 172 do ii=1,3 173 rho1now(:,1:nspden,ii+1)=rho1now(:,1:nspden,ii+1)+nhat1gr(:,1:nspden,ii) 174 end do 175 end if 176 if (test_nhat) then 177 ABI_DEALLOCATE(rhor1_ptr) 178 end if 179 180 !Apply the XC kernel 181 nspgrad=2; if (nspden==2) nspgrad=5 182 ABI_ALLOCATE(dnexcdn,(cplex*nfft,nspgrad)) 183 184 if (cplex==1) then ! Treat real case first 185 if (nspden==1) then 186 do ir=1,nfft 187 r0(:)=kxc(ir,5:7) ; r1(:)=rho1now(ir,1,2:4) 188 gradrho_gradrho1=dot_product(r0,r1) 189 dnexcdn(ir,1)=kxc(ir,1)*rho1now(ir,1,1) + kxc(ir,3)*gradrho_gradrho1 190 coeff_grho=kxc(ir,3)*rho1now(ir,1,1) + kxc(ir,4)*gradrho_gradrho1 191 ! Reuse the storage in rho1now 192 rho1now(ir,1,2:4)=r1(:)*kxc(ir,2)+r0(:)*coeff_grho 193 end do 194 else 195 do ir=1,nfft 196 do ii=1,3 ! grad of spin-up ans spin_dn GS rho 197 r0_up(ii)=kxc(ir,13+2*ii);r0_dn(ii)=kxc(ir,12+2*ii)-kxc(ir,13+2*ii) 198 end do 199 r0(:)=r0_up(:)+r0_dn(:) ! grad of GS rho 200 r1_up(:)=rho1now(ir,1,2:4) ! grad of spin-up rho1 201 r1_dn(:)=rho1now(ir,2,2:4) ! grad of spin-down rho1 202 r1(:)=r1_up(:)+r1_dn(:) ! grad of GS rho1 203 gradrho_gradrho1_up=dot_product(r0_up,r1_up) 204 gradrho_gradrho1_dn=dot_product(r0_dn,r1_dn) 205 gradrho_gradrho1 =dot_product(r0,r1) 206 dnexcdn(ir,1)=kxc(ir, 1)*rho1now(ir,1,1) & 207 & +kxc(ir, 2)*rho1now(ir,2,1) & 208 & +kxc(ir, 6)*gradrho_gradrho1_up & 209 & +kxc(ir,11)*gradrho_gradrho1 210 dnexcdn(ir,2)=kxc(ir, 3)*rho1now(ir,2,1) & 211 & +kxc(ir, 2)*rho1now(ir,1,1) & 212 & +kxc(ir, 7)*gradrho_gradrho1_dn & 213 & +kxc(ir,12)*gradrho_gradrho1 214 coeff_grho_corr=kxc(ir,11)*rho1now(ir,1,1) & 215 & +kxc(ir,12)*rho1now(ir,2,1) & 216 & +kxc(ir,13)*gradrho_gradrho1 217 coeff_grho_up=kxc(ir,6)*rho1now(ir,1,1)+kxc(ir,8)*gradrho_gradrho1_up 218 coeff_grho_dn=kxc(ir,7)*rho1now(ir,2,1)+kxc(ir,9)*gradrho_gradrho1_dn 219 ! Reuse the storage in rho1now 220 rho1now(ir,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1_up(:) & 221 & +kxc(ir,10) *r1_dn(:) & 222 & +coeff_grho_up *r0_up(:) & 223 & +coeff_grho_corr *r0(:) 224 rho1now(ir,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1_dn(:) & 225 & +kxc(ir,10) *r1_up(:) & 226 & +coeff_grho_dn *r0_dn(:) & 227 & +coeff_grho_corr *r0(:) 228 end do 229 end if ! nspden 230 231 else ! if cplex==2 232 if (nspden==1) then 233 do ir=1,nfft 234 r0(:)=kxc(ir,5:7) 235 r1(:) =rho1now(2*ir-1,1,2:4) 236 r1im(:)=rho1now(2*ir ,1,2:4) 237 gradrho_gradrho1 =dot_product(r0,r1) 238 gradrho_gradrho1im=dot_product(r0,r1im) 239 dnexcdn(2*ir-1,1)=kxc(ir,1)*rho1now(2*ir-1,1,1) + kxc(ir,3)*gradrho_gradrho1 240 dnexcdn(2*ir ,1)=kxc(ir,1)*rho1now(2*ir ,1,1) + kxc(ir,3)*gradrho_gradrho1im 241 coeff_grho =kxc(ir,3)*rho1now(2*ir-1,1,1) + kxc(ir,4)*gradrho_gradrho1 242 coeffim_grho=kxc(ir,3)*rho1now(2*ir ,1,1) + kxc(ir,4)*gradrho_gradrho1im 243 ! Reuse the storage in rho1now 244 rho1now(2*ir-1,1,2:4)=r1(:) *kxc(ir,2)+r0(:)*coeff_grho 245 rho1now(2*ir ,1,2:4)=r1im(:)*kxc(ir,2)+r0(:)*coeffim_grho 246 end do 247 else 248 do ir=1,nfft 249 do ii=1,3 ! grad of spin-up ans spin_dn GS rho 250 r0_up(ii)=kxc(ir,13+2*ii);r0_dn(ii)=kxc(ir,12+2*ii)-kxc(ir,13+2*ii) 251 end do 252 r0(:)=r0_up(:)+r0_dn(:) ! grad of GS rho 253 r1_up(:)=rho1now(2*ir-1,1,2:4) ! grad of spin-up rho1 254 r1im_up(:)=rho1now(2*ir,1,2:4) ! grad of spin-up rho1 , im part 255 r1_dn(:)=rho1now(2*ir-1,2,2:4) ! grad of spin-down rho1 256 r1im_dn(:)=rho1now(2*ir,2,2:4) ! grad of spin-down rho1 , im part 257 r1(:)=r1_up(:)+r1_dn(:) ! grad of GS rho1 258 r1im(:)=r1im_up(:)+r1im_dn(:) ! grad of GS rho1, im part 259 gradrho_gradrho1_up =dot_product(r0_up,r1_up) 260 gradrho_gradrho1_dn =dot_product(r0_dn,r1_dn) 261 gradrho_gradrho1 =dot_product(r0,r1) 262 gradrho_gradrho1im_up=dot_product(r0_up,r1im_up) 263 gradrho_gradrho1im_dn=dot_product(r0_dn,r1im_dn) 264 gradrho_gradrho1im =dot_product(r0,r1im) 265 dnexcdn(2*ir-1,1)=kxc(ir, 1)*rho1now(2*ir-1,1,1) & 266 & +kxc(ir, 2)*rho1now(2*ir-1,2,1) & 267 & +kxc(ir, 6)*gradrho_gradrho1_up & 268 & +kxc(ir,11)*gradrho_gradrho1 269 dnexcdn(2*ir-1,2)=kxc(ir, 3)*rho1now(2*ir-1,2,1) & 270 & +kxc(ir, 2)*rho1now(2*ir-1,1,1) & 271 & +kxc(ir, 7)*gradrho_gradrho1_dn & 272 & +kxc(ir,12)*gradrho_gradrho1 273 dnexcdn(2*ir ,1)=kxc(ir, 1)*rho1now(2*ir ,1,1) & 274 & +kxc(ir, 2)*rho1now(2*ir ,2,1) & 275 & +kxc(ir, 6)*gradrho_gradrho1im_up & 276 & +kxc(ir,11)*gradrho_gradrho1im 277 dnexcdn(2*ir ,2)=kxc(ir, 3)*rho1now(2*ir ,2,1) & 278 & +kxc(ir, 2)*rho1now(2*ir ,1,1) & 279 & +kxc(ir, 7)*gradrho_gradrho1im_dn & 280 & +kxc(ir,12)*gradrho_gradrho1im 281 coeff_grho_corr =kxc(ir,11)*rho1now(2*ir-1,1,1) & 282 & +kxc(ir,12)*rho1now(2*ir-1,2,1) & 283 & +kxc(ir,13)*gradrho_gradrho1 284 coeffim_grho_corr=kxc(ir,11)*rho1now(2*ir ,1,1) & 285 & +kxc(ir,12)*rho1now(2*ir ,2,1) & 286 & +kxc(ir,13)*gradrho_gradrho1im 287 coeff_grho_up =kxc(ir,6)*rho1now(2*ir-1,1,1)+kxc(ir,8)*gradrho_gradrho1_up 288 coeff_grho_dn =kxc(ir,7)*rho1now(2*ir-1,2,1)+kxc(ir,9)*gradrho_gradrho1_dn 289 coeffim_grho_up=kxc(ir,6)*rho1now(2*ir ,1,1)+kxc(ir,8)*gradrho_gradrho1im_up 290 coeffim_grho_dn=kxc(ir,7)*rho1now(2*ir ,2,1)+kxc(ir,9)*gradrho_gradrho1im_dn 291 ! Reuse the storage in rho1now 292 rho1now(2*ir-1,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1_up(:) & 293 & +kxc(ir,10) *r1_dn(:) & 294 & +coeff_grho_up *r0_up(:) & 295 & +coeff_grho_corr*r0(:) 296 rho1now(2*ir-1,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1_dn(:) & 297 & +kxc(ir,10) *r1_up(:) & 298 & +coeff_grho_dn *r0_dn(:) & 299 & +coeff_grho_corr*r0(:) 300 rho1now(2*ir ,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1im_up(:) & 301 & +kxc(ir,10) *r1im_dn(:) & 302 & +coeffim_grho_up *r0_up(:) & 303 & +coeffim_grho_corr *r0(:) 304 rho1now(2*ir ,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1im_dn(:) & 305 & +kxc(ir,10) *r1im_up(:) & 306 & +coeffim_grho_dn *r0_dn(:) & 307 & +coeffim_grho_corr *r0(:) 308 end do 309 end if ! nspden 310 311 end if 312 313 vxc1(:,:)=zero 314 call xcpot(cplex,dnexcdn,gprimd,ishift,mgga,mpi_enreg,nfft,ngfft,ngrad,nspden,& 315 & nspgrad,paral_kgb,qphon,rho1now,vxc1) 316 317 !call filterpot(paral_kgb,cplex,gmet,gsqcut,nfft,ngfft,nspden,qphon,vxc1) 318 319 ABI_DEALLOCATE(dnexcdn) 320 ABI_DEALLOCATE(rho1now) 321 322 DBG_EXIT("COLL") 323 324 end subroutine dfpt_mkvxcgga