TABLE OF CONTENTS
- ABINIT/dfpt_mkvxc
- ABINIT/dfpt_mkvxc_noncoll
- ABINIT/dfpt_mkvxcgga
- ABINIT/dfpt_mkvxcgga_n0met
- ABINIT/dfpt_mkvxcggadq
- ABINIT/m_dfpt_mkvxc
ABINIT/dfpt_mkvxc [ Functions ]
NAME
dfpt_mkvxc
FUNCTION
Compute the first-order change of exchange-correlation potential due to atomic displacement: assemble the first-order density change with the frozen-core density change, then use the exchange-correlation kernel.
INPUTS
cplex= if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX ixc= choice of exchange-correlation scheme kxc(nfft,nkxc)=exchange and correlation kernel (see below) mpi_enreg=information about MPI parallelization 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 nhat1(cplex*nfft,2nspden*nhat1dim)= -PAW only- 1st-order compensation density nhat1dim= -PAW only- 1 if nhat1 array is used ; 0 otherwise nhat1gr(cplex*nfft,nspden,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 non_magnetic_xc= if true, handle density/potential as non-magnetic (even if it is) nspden=number of spin-density components n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used, otherwise, nfft option=if 0, work only with the XC core-correction, if 1, treat both density change and XC core correction if 2, treat only density change qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2). rhor1(cplex*nfft,nspden)=array for electron density in electrons/bohr**3. rprimd(3,3)=dimensional primitive translations in real space (bohr) usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc
OUTPUT
vxc1(cplex*nfft,nspden)=change in exchange-correlation potential (including core-correction, if applicable)
NOTES
Content of Kxc array:
===== if LDA
if nspden==1: kxc(:,1)= d2Exc/drho2
(kxc(:,2)= d2Exc/drho_up drho_dn)
if nspden>=2: kxc(:,1)= d2Exc/drho_up drho_up
kxc(:,2)= d2Exc/drho_up drho_dn
kxc(:,3)= d2Exc/drho_dn drho_dn
===== if GGA (or mGGA)
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)
Note about mGGA: 2nd derivatives involving Tau or Laplacian are not taken into account (yet)
SOURCE
129 subroutine dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,nhat1dim,nhat1gr,nhat1grdim,& 130 & nkxc,non_magnetic_xc,nspden,n3xccc,option,qphon,rhor1,rprimd,usexcnhat,vxc1,xccc3d1) 131 132 !Arguments ------------------------------------ 133 !scalars 134 integer,intent(in) :: cplex,ixc,n3xccc,nfft,nhat1dim,nhat1grdim 135 integer,intent(in) :: nkxc,nspden,option,usexcnhat 136 logical,intent(in) :: non_magnetic_xc 137 type(MPI_type),intent(in) :: mpi_enreg 138 !arrays 139 integer,intent(in) :: ngfft(18) 140 real(dp),intent(in),target :: nhat1(cplex*nfft,nspden*nhat1dim) 141 real(dp),intent(in),target :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim) 142 real(dp),intent(in) :: kxc(nfft,nkxc),qphon(3) 143 real(dp),intent(in),target :: rhor1(cplex*nfft,nspden) 144 real(dp),intent(in) :: rprimd(3,3),xccc3d1(cplex*n3xccc) 145 real(dp),intent(out) :: vxc1(cplex*nfft,nspden) 146 147 !Local variables------------------------------- 148 !scalars 149 integer :: ii,ir,ispden,nhat1dim_,nhat1rgdim_ 150 real(dp) :: rho1_dn,rho1_up,rho1im_dn,rho1im_up,rho1re_dn,rho1re_up 151 real(dp) :: spin_scale 152 !arrays 153 real(dp) :: gprimd(3,3),tsec(2) 154 real(dp), ABI_CONTIGUOUS pointer :: nhat1_(:,:),nhat1gr_(:,:,:),rhor1_(:,:) 155 156 ! ************************************************************************* 157 158 DBG_ENTER("COLL") 159 160 call timab(181,1,tsec) 161 162 if(nspden/=1 .and. nspden/=2) then 163 ABI_BUG('For nspden==4 please use dfpt_mkvxc_noncoll!') 164 end if 165 166 !Special case: no XC applied 167 if (ixc==0.or.nkxc==0) then 168 ABI_WARNING('Note that no xc is applied (ixc=0)') 169 vxc1=zero 170 return 171 end if 172 173 !Treat first LDA 174 if(nkxc==1.or.nkxc==3)then 175 176 ! PAW: eventually substract compensation density 177 if (option/=0) then 178 if ((usexcnhat==0.and.nhat1dim==1).or.(non_magnetic_xc)) then 179 ABI_MALLOC(rhor1_,(cplex*nfft,nspden)) 180 if (usexcnhat==0.and.nhat1dim==1) then 181 rhor1_(:,:)=rhor1(:,:)-nhat1(:,:) 182 else 183 rhor1_(:,:)=rhor1(:,:) 184 end if 185 if (non_magnetic_xc) then 186 if(nspden==2) rhor1_(:,2)=rhor1_(:,1)*half 187 if(nspden==4) rhor1_(:,2:4)=zero 188 end if 189 else 190 rhor1_ => rhor1 191 end if 192 end if 193 194 ! Case without non-linear core correction 195 if(n3xccc==0 .or. option==2)then 196 197 if(option==0)then ! No straight XC to compute 198 199 vxc1(:,:)=zero 200 201 else ! XC, without non-linear XC correction 202 203 ! Non-spin-polarized 204 if(nspden==1)then 205 if(cplex==1)then 206 do ir=1,nfft 207 vxc1(ir,1)=kxc(ir,1)*rhor1_(ir,1) 208 end do 209 else 210 do ir=1,nfft 211 vxc1(2*ir-1,1)=kxc(ir,1)*rhor1_(2*ir-1,1) 212 vxc1(2*ir ,1)=kxc(ir,1)*rhor1_(2*ir ,1) 213 end do 214 end if ! cplex==1 215 216 ! Spin-polarized 217 else 218 if(cplex==1)then 219 do ir=1,nfft 220 rho1_dn=rhor1_(ir,1)-rhor1_(ir,2) 221 vxc1(ir,1)=kxc(ir,1)*rhor1_(ir,2)+kxc(ir,2)*rho1_dn 222 vxc1(ir,2)=kxc(ir,2)*rhor1_(ir,2)+kxc(ir,3)*rho1_dn 223 end do 224 else 225 do ir=1,nfft 226 rho1re_dn=rhor1_(2*ir-1,1)-rhor1_(2*ir-1,2) 227 rho1im_dn=rhor1_(2*ir ,1)-rhor1_(2*ir ,2) 228 vxc1(2*ir-1,1)=kxc(ir,1)*rhor1_(2*ir-1,2)+kxc(ir,2)*rho1re_dn 229 vxc1(2*ir ,1)=kxc(ir,1)*rhor1_(2*ir ,2)+kxc(ir,2)*rho1im_dn 230 vxc1(2*ir-1,2)=kxc(ir,2)*rhor1_(2*ir-1,2)+kxc(ir,3)*rho1re_dn 231 vxc1(2*ir ,2)=kxc(ir,2)*rhor1_(2*ir ,2)+kxc(ir,3)*rho1im_dn 232 end do 233 end if ! cplex==1 234 end if ! nspden==1 235 236 end if ! option==0 237 238 ! Treat case with non-linear core correction 239 else 240 241 if(option==0)then 242 243 if(nspden==1)then 244 if(cplex==1)then 245 do ir=1,nfft 246 vxc1(ir,1)=kxc(ir,1)*xccc3d1(ir) 247 end do 248 else 249 do ir=1,nfft 250 vxc1(2*ir-1,1)=kxc(ir,1)*xccc3d1(2*ir-1) 251 vxc1(2*ir ,1)=kxc(ir,1)*xccc3d1(2*ir ) 252 end do 253 end if ! cplex==1 254 else 255 if(cplex==1)then 256 do ir=1,nfft 257 vxc1(ir,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(ir)*half 258 vxc1(ir,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(ir)*half 259 end do 260 else 261 do ir=1,nfft 262 vxc1(2*ir-1,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(2*ir-1)*half 263 vxc1(2*ir ,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(2*ir )*half 264 vxc1(2*ir-1,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(2*ir-1)*half 265 vxc1(2*ir ,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(2*ir )*half 266 end do 267 end if ! cplex==1 268 end if ! nspden==1 269 270 else ! option/=0 271 272 if(nspden==1)then 273 if(cplex==1)then 274 do ir=1,nfft 275 vxc1(ir,1)=kxc(ir,1)*(rhor1_(ir,1)+xccc3d1(ir)) 276 end do 277 else 278 do ir=1,nfft 279 vxc1(2*ir-1,1)=kxc(ir,1)*(rhor1_(2*ir-1,1)+xccc3d1(2*ir-1)) 280 vxc1(2*ir ,1)=kxc(ir,1)*(rhor1_(2*ir ,1)+xccc3d1(2*ir )) 281 end do 282 end if ! cplex==1 283 else 284 if(cplex==1)then 285 do ir=1,nfft 286 rho1_dn=rhor1_(ir,1)-rhor1_(ir,2) + xccc3d1(ir)*half 287 rho1_up=rhor1_(ir,2) + xccc3d1(ir)*half 288 vxc1(ir,1)=kxc(ir,1)*rho1_up+kxc(ir,2)*rho1_dn 289 vxc1(ir,2)=kxc(ir,2)*rho1_up+kxc(ir,3)*rho1_dn 290 end do 291 else 292 do ir=1,nfft 293 rho1re_dn=rhor1_(2*ir-1,1)-rhor1_(2*ir-1,2) + xccc3d1(2*ir-1)*half 294 rho1im_dn=rhor1_(2*ir ,1)-rhor1_(2*ir ,2) + xccc3d1(2*ir )*half 295 rho1re_up=rhor1_(2*ir-1,2) + xccc3d1(2*ir-1)*half 296 rho1im_up=rhor1_(2*ir ,2) + xccc3d1(2*ir )*half 297 vxc1(2*ir-1,1)=kxc(ir,1)*rho1re_up+kxc(ir,2)*rho1re_dn 298 vxc1(2*ir ,1)=kxc(ir,1)*rho1im_up+kxc(ir,2)*rho1im_dn 299 vxc1(2*ir-1,2)=kxc(ir,2)*rho1re_up+kxc(ir,3)*rho1re_dn 300 vxc1(2*ir ,2)=kxc(ir,2)*rho1im_up+kxc(ir,3)*rho1im_dn 301 end do 302 end if ! cplex==1 303 end if ! nspden==1 304 305 end if ! option==0 306 307 end if ! n3xccc==0 308 309 if (option/=0.and.((usexcnhat==0.and.nhat1dim==1).or.(non_magnetic_xc))) then 310 ABI_FREE(rhor1_) 311 end if 312 313 ! Treat GGA 314 else if (nkxc==7.or.nkxc==19) then 315 316 ! Transfer the data to spin-polarized storage 317 318 ! Treat the density change 319 ABI_MALLOC(rhor1_,(cplex*nfft,nspden)) 320 if (option==1 .or. option==2) then 321 if (nspden==1) then 322 do ir=1,cplex*nfft 323 rhor1_(ir,1)=rhor1(ir,1) 324 end do 325 else 326 if(non_magnetic_xc) then 327 do ir=1,cplex*nfft 328 rho1_dn=rhor1(ir,1)*half 329 rhor1_(ir,1)=rho1_dn 330 rhor1_(ir,2)=rho1_dn 331 end do 332 else 333 do ir=1,cplex*nfft 334 rho1_dn=rhor1(ir,1)-rhor1(ir,2) 335 rhor1_(ir,1)=rhor1(ir,2) 336 rhor1_(ir,2)=rho1_dn 337 end do 338 end if 339 end if 340 else 341 do ispden=1,nspden 342 do ir=1,cplex*nfft 343 rhor1_(ir,ispden)=zero 344 end do 345 end do 346 end if 347 348 if( (option==0 .or. option==1) .and. n3xccc/=0)then 349 spin_scale=one;if (nspden==2) spin_scale=half 350 do ispden=1,nspden 351 do ir=1,cplex*nfft 352 rhor1_(ir,ispden)=rhor1_(ir,ispden)+xccc3d1(ir)*spin_scale 353 end do 354 end do 355 end if 356 357 ! PAW: treat also compensation density (and gradients) 358 nhat1dim_=nhat1dim ; nhat1rgdim_=nhat1grdim 359 if (option/=0.and.nhat1dim==1.and.nspden==2) then 360 ABI_MALLOC(nhat1_,(cplex*nfft,nspden)) 361 if (non_magnetic_xc) then 362 do ir=1,cplex*nfft 363 rho1_dn=nhat1(ir,1)*half 364 nhat1_(ir,1:2)=rho1_dn 365 end do 366 else 367 do ir=1,cplex*nfft 368 rho1_dn=nhat1(ir,1)-nhat1(ir,2) 369 nhat1_(ir,1)=nhat1(ir,2) 370 nhat1_(ir,2)=rho1_dn 371 end do 372 end if 373 else if (option==0) then 374 ABI_MALLOC(nhat1_,(0,0)) 375 nhat1dim_=0 376 else 377 nhat1_ => nhat1 378 end if 379 if (option/=0.and.nhat1grdim==1.and.nspden==2) then 380 ABI_MALLOC(nhat1gr_,(cplex*nfft,nspden,3)) 381 if (non_magnetic_xc) then 382 do ii=1,3 383 do ir=1,cplex*nfft 384 rho1_dn=nhat1(ir,1)*half 385 nhat1gr_(ir,1:2,ii)=rho1_dn 386 end do 387 end do 388 else 389 do ii=1,3 390 do ir=1,cplex*nfft 391 rho1_dn=nhat1gr(ir,1,ii)-nhat1gr(ir,2,ii) 392 nhat1gr_(ir,1,ii)=nhat1gr(ir,2,ii) 393 nhat1gr_(ir,2,ii)=rho1_dn 394 end do 395 end do 396 end if 397 else if (option==0) then 398 ABI_MALLOC(nhat1gr_,(0,0,0)) 399 nhat1rgdim_=0 400 else 401 nhat1gr_ => nhat1gr 402 end if 403 404 call matr3inv(rprimd,gprimd) 405 406 call dfpt_mkvxcgga(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat1_,nhat1dim_,& 407 & nhat1gr_,nhat1rgdim_,nkxc,nspden,qphon,rhor1_,usexcnhat,vxc1) 408 409 ABI_FREE(rhor1_) 410 if ((option==0).or.(nhat1dim==1.and.nspden==2)) then 411 ABI_FREE(nhat1_) 412 end if 413 if ((option==0).or.(nhat1grdim==1.and.nspden==2)) then 414 ABI_FREE(nhat1gr_) 415 end if 416 417 else 418 ABI_BUG('Invalid nkxc!') 419 420 end if ! LDA or GGA 421 422 call timab(181,2,tsec) 423 424 DBG_EXIT("COLL") 425 426 end subroutine dfpt_mkvxc
ABINIT/dfpt_mkvxc_noncoll [ Functions ]
NAME
dfpt_mkvxc_noncoll
FUNCTION
Compute the first-order change of exchange-correlation potential due to atomic displacement for non-collinear spins: assemble the first-order density change with the frozen-core density change, then use the exchange-correlation kernel.
INPUTS
cplex= if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX ixc= choice of exchange-correlation scheme ixcrot= option for rotation of collinear spin potential to non collinear full matrix kxc(nfft,nkxc)=exchange and correlation kernel (see rhotoxc.F90) mpi_enreg=information about MPI parallelization 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 nhat(nfft,nspden*nhatdim)= -PAW only- GS compensation density nhatdim= -PAW only- 1 if nhat array is used ; 0 otherwise nhat1(cplex*nfft,nspden*nhat1dim)= -PAW only- 1st-order compensation density nhat1dim= -PAW only- 1 if nhat1 array is used ; 0 otherwise nhat1gr(cplex*nfft,nspden,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 non_magnetic_xc= if true, handle density/potential as non-magnetic (even if it is) nspden=number of spin-density components n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used, otherwise, nfft optnc=option for non-collinear magnetism (nspden=4): 1: the whole 2x2 Vres matrix is computed 2: only Vres^{11} and Vres^{22} are computed option=if 0, work only with the XC core-correction, if 1, treat both density change and XC core correction if 2, treat only density change qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2). rhor(nfft,nspden)=GS electron density in real space rhor1(cplex*nfft,nspden)=1st-order electron density in real space rprimd(3,3)=dimensional primitive translations in real space (bohr) usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc vxc(nfft,nspden)=GS XC potential
OUTPUT
vxc1(cplex*nfft,nspden)=change in exchange-correlation potential (including core-correction, if applicable)
SOURCE
770 subroutine dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat,nhatdim,nhat1,nhat1dim,& 771 & nhat1gr,nhat1grdim,nkxc,non_magnetic_xc,nspden,n3xccc,optnc,option,qphon,& 772 & rhor,rhor1,rprimd,usexcnhat,vxc,vxc1,xccc3d1,ixcrot) 773 774 !Arguments ------------------------------------ 775 !scalars 776 integer,intent(in) :: cplex,ixc,n3xccc,nfft,nhatdim,nhat1dim,nhat1grdim,optnc 777 integer,intent(in) :: nkxc,nspden,option,usexcnhat 778 logical,intent(in) :: non_magnetic_xc 779 type(MPI_type),intent(in) :: mpi_enreg 780 !arrays 781 integer,intent(in) :: ngfft(18) 782 real(dp),intent(in) :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim) 783 real(dp),intent(in) :: kxc(nfft,nkxc) 784 real(dp),intent(in) :: vxc(nfft,nspden) 785 real(dp),intent(in) :: nhat(nfft,nspden*nhatdim),nhat1(cplex*nfft,nspden*nhat1dim) 786 real(dp),intent(in),target :: rhor(nfft,nspden),rhor1(cplex*nfft,nspden) 787 real(dp),intent(in) :: qphon(3),rprimd(3,3),xccc3d1(cplex*n3xccc) 788 real(dp),intent(out) :: vxc1(cplex*nfft,nspden) 789 integer,optional,intent(in) :: ixcrot 790 !Local variables------------------------------- 791 !scalars 792 !arrays 793 real(dp) :: nhat1_zero(0,0),nhat1gr_zero(0,0,0),tsec(2) 794 real(dp),allocatable :: m_norm(:),rhor1_diag(:,:),vxc1_diag(:,:) 795 real(dp), ABI_CONTIGUOUS pointer :: mag(:,:),rhor_(:,:),rhor1_(:,:) 796 ! ************************************************************************* 797 798 ! Non-collinear magnetism 799 ! Has to locally "rotate" rho(r)^(1) (according to magnetization), 800 ! Compute Vxc(r)^(1) in the spin frame aligned with \vec{m} and rotate it back 801 802 DBG_ENTER("COLL") 803 ABI_UNUSED(nhat1gr) 804 805 call timab(181,1,tsec) 806 807 if(nspden/=4) then 808 ABI_BUG('only for nspden=4!') 809 end if 810 811 if(nkxc/=2*min(nspden,2)-1) then 812 ABI_BUG('nspden=4 works only with LSDA.') 813 end if 814 815 !Special case: no XC applied 816 if (ixc==0.or.nkxc==0) then 817 ABI_WARNING('Note that no xc is applied (ixc=0)') 818 vxc1(:,:)=zero 819 return 820 end if 821 822 823 824 !Treat first LDA 825 if(nkxc==1.or.nkxc==3)then 826 827 vxc1(:,:)=zero 828 829 ! PAW: possibly substract compensation density 830 if ((usexcnhat==0.and.nhatdim==1).or.(non_magnetic_xc)) then 831 ABI_MALLOC(rhor_,(nfft,nspden)) 832 if (usexcnhat==0.and.nhatdim==1) then 833 rhor_(:,:) =rhor(:,:)-nhat(:,:) 834 else 835 rhor_(:,:) =rhor(:,:) 836 end if 837 if (non_magnetic_xc) then 838 if(nspden==2) rhor_(:,2)=rhor_(:,1)*half 839 if(nspden==4) rhor_(:,2:4)=zero 840 end if 841 else 842 rhor_ => rhor 843 end if 844 if ((usexcnhat==0.and.nhat1dim==1).or.(non_magnetic_xc)) then 845 ABI_MALLOC(rhor1_,(cplex*nfft,nspden)) 846 if (usexcnhat==0.and.nhatdim==1) then 847 rhor1_(:,:)=rhor1(:,:)-nhat1(:,:) 848 else 849 rhor1_(:,:)=rhor1(:,:) 850 end if 851 if (non_magnetic_xc) then 852 if(nspden==2) rhor1_(:,2)=rhor1_(:,1)*half 853 if(nspden==4) rhor1_(:,2:4)=zero 854 end if 855 else 856 rhor1_ => rhor1 857 end if 858 859 ! Magnetization 860 mag => rhor_(:,2:4) 861 ABI_MALLOC(rhor1_diag,(cplex*nfft,2)) 862 ABI_MALLOC(vxc1_diag,(cplex*nfft,2)) 863 ABI_MALLOC(m_norm,(nfft)) 864 865 ! -- Rotate rho(r)^(1) 866 ! SPr: for option=0 the rhor is not used, only core density xccc3d1 867 ! rotate_mag is only to compute the m_norm 868 call rotate_mag(rhor1_,rhor1_diag,mag,nfft,cplex,mag_norm_out=m_norm,& 869 & rho_out_format=2) 870 871 ! -- Compute Vxc(r)^(1)=Kxc(r).rho(r)^(1)_rotated 872 ! Note for PAW: nhat has already been substracted; don't use it in dfpt_mkvxc 873 ! (put all nhat options to zero). 874 ! The collinear routine dfpt_mkvxc wants a general density built as (tr[rho],rho_upup) 875 call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1_zero,0,nhat1gr_zero,0,& 876 & nkxc,non_magnetic_xc,2,n3xccc,option,qphon,rhor1_diag,rprimd,0,vxc1_diag,xccc3d1) 877 878 !call test_rotations(0,1) 879 880 ! -- Rotate back Vxc(r)^(1) 881 if (optnc==1) then 882 if(present(ixcrot)) then 883 call rotate_back_mag_dfpt(option,vxc1_diag,vxc1,vxc,kxc,rhor1_,mag,nfft,cplex,& 884 & mag_norm_in=m_norm,rot_method=ixcrot) 885 else 886 call rotate_back_mag_dfpt(option,vxc1_diag,vxc1,vxc,kxc,rhor1_,mag,nfft,cplex,& 887 & mag_norm_in=m_norm) 888 end if 889 else 890 call rotate_back_mag(vxc1_diag,vxc1,mag,nfft,mag_norm_in=m_norm) 891 vxc1(:,3:4)=zero 892 end if 893 894 ABI_FREE(rhor1_diag) 895 ABI_FREE(vxc1_diag) 896 ABI_FREE(m_norm) 897 if ((usexcnhat==0.and.nhatdim==1).or.(non_magnetic_xc)) then 898 ABI_FREE(rhor_) 899 end if 900 if ((usexcnhat==0.and.nhat1dim==1).or.(non_magnetic_xc)) then 901 ABI_FREE(rhor1_) 902 end if 903 904 end if ! nkxc=1 or nkxc=3 905 906 call timab(181,2,tsec) 907 908 DBG_EXIT("COLL") 909 910 end subroutine dfpt_mkvxc_noncoll
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.
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=information 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
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)
SOURCE
497 subroutine dfpt_mkvxcgga(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,& 498 & nhat1,nhat1dim,nhat1gr,nhat1grdim,nkxc,& 499 & nspden,qphon,rhor1,usexcnhat,vxc1) 500 501 !Arguments ------------------------------------ 502 !scalars 503 integer,intent(in) :: cplex,nfft,nhat1dim,nhat1grdim,nkxc,nspden,usexcnhat 504 type(MPI_type),intent(in) :: mpi_enreg 505 !arrays 506 integer,intent(in) :: ngfft(18) 507 real(dp),intent(in) :: gprimd(3,3) 508 real(dp),intent(in) :: kxc(nfft,nkxc) 509 real(dp),intent(in) :: nhat1(cplex*nfft,nspden*nhat1dim) 510 real(dp),intent(in) :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim) 511 real(dp),intent(in) :: qphon(3) 512 real(dp),intent(in),target :: rhor1(cplex*nfft,nspden) 513 real(dp),intent(out) :: vxc1(cplex*nfft,nspden) 514 515 !Local variables------------------------------- 516 !scalars 517 integer :: ii,ir,ishift,ngrad,nspgrad,use_laplacian 518 logical :: test_nhat 519 real(dp) :: coeff_grho,coeff_grho_corr,coeff_grho_dn,coeff_grho_up 520 real(dp) :: coeffim_grho,coeffim_grho_corr,coeffim_grho_dn,coeffim_grho_up 521 real(dp) :: gradrho_gradrho1,gradrho_gradrho1_dn,gradrho_gradrho1_up 522 real(dp) :: gradrho_gradrho1im,gradrho_gradrho1im_dn,gradrho_gradrho1im_up 523 character(len=500) :: msg 524 !arrays 525 real(dp) :: r0(3),r0_dn(3),r0_up(3),r1(3),r1_dn(3),r1_up(3) 526 real(dp) :: r1im(3),r1im_dn(3),r1im_up(3) 527 real(dp),allocatable :: dnexcdn(:,:),rho1now(:,:,:) 528 real(dp),ABI_CONTIGUOUS pointer :: rhor1_ptr(:,:) 529 530 ! ************************************************************************* 531 532 DBG_ENTER("COLL") 533 534 if (nkxc/=12*min(nspden,2)-5) then 535 msg='Wrong nkxc value for GGA!' 536 ABI_BUG(msg) 537 end if 538 539 !metaGGA contributions are not taken into account here 540 use_laplacian=0 541 542 !PAW: substract 1st-order compensation density from 1st-order density 543 test_nhat=((nhat1dim==1).and.(usexcnhat==0.or.nhat1grdim==1)) 544 if (test_nhat) then 545 ABI_MALLOC(rhor1_ptr,(cplex*nfft,nspden)) 546 rhor1_ptr(:,:)=rhor1(:,:)-nhat1(:,:) 547 else 548 rhor1_ptr => rhor1 549 end if 550 551 !call filterpot(paral_kgb,cplex,gmet,gsqcut,nfft,ngfft,2,qphon,rhor1_ptr) 552 553 !Compute the gradients of the first-order density 554 !rho1now(:,:,1) contains the first-order density, and 555 !rho1now(:,:,2:4) contains the gradients of the first-order density 556 ishift=0 ; ngrad=2 557 ABI_MALLOC(rho1now,(cplex*nfft,nspden,ngrad*ngrad)) 558 call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,qphon,rhor1_ptr,rho1now) 559 560 !PAW: add "exact" gradients of compensation density 561 if (test_nhat.and.usexcnhat==1) then 562 rho1now(:,1:nspden,1)=rho1now(:,1:nspden,1)+nhat1(:,1:nspden) 563 end if 564 if (nhat1grdim==1) then 565 do ii=1,3 566 rho1now(:,1:nspden,ii+1)=rho1now(:,1:nspden,ii+1)+nhat1gr(:,1:nspden,ii) 567 end do 568 end if 569 if (test_nhat) then 570 ABI_FREE(rhor1_ptr) 571 end if 572 573 !Apply the XC kernel 574 nspgrad=2; if (nspden==2) nspgrad=5 575 ABI_MALLOC(dnexcdn,(cplex*nfft,nspgrad)) 576 577 if (cplex==1) then ! Treat real case first 578 if (nspden==1) then 579 do ir=1,nfft 580 r0(:)=kxc(ir,5:7) ; r1(:)=rho1now(ir,1,2:4) 581 gradrho_gradrho1=dot_product(r0,r1) 582 dnexcdn(ir,1)=kxc(ir,1)*rho1now(ir,1,1) + kxc(ir,3)*gradrho_gradrho1 583 coeff_grho=kxc(ir,3)*rho1now(ir,1,1) + kxc(ir,4)*gradrho_gradrho1 584 ! Reuse the storage in rho1now 585 rho1now(ir,1,2:4)=r1(:)*kxc(ir,2)+r0(:)*coeff_grho 586 end do 587 else 588 do ir=1,nfft 589 do ii=1,3 ! grad of spin-up ans spin_dn GS rho 590 r0_up(ii)=kxc(ir,13+2*ii);r0_dn(ii)=kxc(ir,12+2*ii)-kxc(ir,13+2*ii) 591 end do 592 r0(:)=r0_up(:)+r0_dn(:) ! grad of GS rho 593 r1_up(:)=rho1now(ir,1,2:4) ! grad of spin-up rho1 594 r1_dn(:)=rho1now(ir,2,2:4) ! grad of spin-down rho1 595 r1(:)=r1_up(:)+r1_dn(:) ! grad of GS rho1 596 gradrho_gradrho1_up=dot_product(r0_up,r1_up) 597 gradrho_gradrho1_dn=dot_product(r0_dn,r1_dn) 598 gradrho_gradrho1 =dot_product(r0,r1) 599 dnexcdn(ir,1)=kxc(ir, 1)*rho1now(ir,1,1) & 600 & +kxc(ir, 2)*rho1now(ir,2,1) & 601 & +kxc(ir, 6)*gradrho_gradrho1_up & 602 & +kxc(ir,11)*gradrho_gradrho1 603 dnexcdn(ir,2)=kxc(ir, 3)*rho1now(ir,2,1) & 604 & +kxc(ir, 2)*rho1now(ir,1,1) & 605 & +kxc(ir, 7)*gradrho_gradrho1_dn & 606 & +kxc(ir,12)*gradrho_gradrho1 607 coeff_grho_corr=kxc(ir,11)*rho1now(ir,1,1) & 608 & +kxc(ir,12)*rho1now(ir,2,1) & 609 & +kxc(ir,13)*gradrho_gradrho1 610 coeff_grho_up=kxc(ir,6)*rho1now(ir,1,1)+kxc(ir,8)*gradrho_gradrho1_up 611 coeff_grho_dn=kxc(ir,7)*rho1now(ir,2,1)+kxc(ir,9)*gradrho_gradrho1_dn 612 ! Reuse the storage in rho1now 613 rho1now(ir,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1_up(:) & 614 & +kxc(ir,10) *r1_dn(:) & 615 & +coeff_grho_up *r0_up(:) & 616 & +coeff_grho_corr *r0(:) 617 rho1now(ir,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1_dn(:) & 618 & +kxc(ir,10) *r1_up(:) & 619 & +coeff_grho_dn *r0_dn(:) & 620 & +coeff_grho_corr *r0(:) 621 end do 622 end if ! nspden 623 624 else ! if cplex==2 625 if (nspden==1) then 626 do ir=1,nfft 627 r0(:)=kxc(ir,5:7) 628 r1(:) =rho1now(2*ir-1,1,2:4) 629 r1im(:)=rho1now(2*ir ,1,2:4) 630 gradrho_gradrho1 =dot_product(r0,r1) 631 gradrho_gradrho1im=dot_product(r0,r1im) 632 dnexcdn(2*ir-1,1)=kxc(ir,1)*rho1now(2*ir-1,1,1) + kxc(ir,3)*gradrho_gradrho1 633 dnexcdn(2*ir ,1)=kxc(ir,1)*rho1now(2*ir ,1,1) + kxc(ir,3)*gradrho_gradrho1im 634 coeff_grho =kxc(ir,3)*rho1now(2*ir-1,1,1) + kxc(ir,4)*gradrho_gradrho1 635 coeffim_grho=kxc(ir,3)*rho1now(2*ir ,1,1) + kxc(ir,4)*gradrho_gradrho1im 636 ! Reuse the storage in rho1now 637 rho1now(2*ir-1,1,2:4)=r1(:) *kxc(ir,2)+r0(:)*coeff_grho 638 rho1now(2*ir ,1,2:4)=r1im(:)*kxc(ir,2)+r0(:)*coeffim_grho 639 end do 640 else 641 do ir=1,nfft 642 do ii=1,3 ! grad of spin-up ans spin_dn GS rho 643 r0_up(ii)=kxc(ir,13+2*ii);r0_dn(ii)=kxc(ir,12+2*ii)-kxc(ir,13+2*ii) 644 end do 645 r0(:)=r0_up(:)+r0_dn(:) ! grad of GS rho 646 r1_up(:)=rho1now(2*ir-1,1,2:4) ! grad of spin-up rho1 647 r1im_up(:)=rho1now(2*ir,1,2:4) ! grad of spin-up rho1 , im part 648 r1_dn(:)=rho1now(2*ir-1,2,2:4) ! grad of spin-down rho1 649 r1im_dn(:)=rho1now(2*ir,2,2:4) ! grad of spin-down rho1 , im part 650 r1(:)=r1_up(:)+r1_dn(:) ! grad of GS rho1 651 r1im(:)=r1im_up(:)+r1im_dn(:) ! grad of GS rho1, im part 652 gradrho_gradrho1_up =dot_product(r0_up,r1_up) 653 gradrho_gradrho1_dn =dot_product(r0_dn,r1_dn) 654 gradrho_gradrho1 =dot_product(r0,r1) 655 gradrho_gradrho1im_up=dot_product(r0_up,r1im_up) 656 gradrho_gradrho1im_dn=dot_product(r0_dn,r1im_dn) 657 gradrho_gradrho1im =dot_product(r0,r1im) 658 dnexcdn(2*ir-1,1)=kxc(ir, 1)*rho1now(2*ir-1,1,1) & 659 & +kxc(ir, 2)*rho1now(2*ir-1,2,1) & 660 & +kxc(ir, 6)*gradrho_gradrho1_up & 661 & +kxc(ir,11)*gradrho_gradrho1 662 dnexcdn(2*ir-1,2)=kxc(ir, 3)*rho1now(2*ir-1,2,1) & 663 & +kxc(ir, 2)*rho1now(2*ir-1,1,1) & 664 & +kxc(ir, 7)*gradrho_gradrho1_dn & 665 & +kxc(ir,12)*gradrho_gradrho1 666 dnexcdn(2*ir ,1)=kxc(ir, 1)*rho1now(2*ir ,1,1) & 667 & +kxc(ir, 2)*rho1now(2*ir ,2,1) & 668 & +kxc(ir, 6)*gradrho_gradrho1im_up & 669 & +kxc(ir,11)*gradrho_gradrho1im 670 dnexcdn(2*ir ,2)=kxc(ir, 3)*rho1now(2*ir ,2,1) & 671 & +kxc(ir, 2)*rho1now(2*ir ,1,1) & 672 & +kxc(ir, 7)*gradrho_gradrho1im_dn & 673 & +kxc(ir,12)*gradrho_gradrho1im 674 coeff_grho_corr =kxc(ir,11)*rho1now(2*ir-1,1,1) & 675 & +kxc(ir,12)*rho1now(2*ir-1,2,1) & 676 & +kxc(ir,13)*gradrho_gradrho1 677 coeffim_grho_corr=kxc(ir,11)*rho1now(2*ir ,1,1) & 678 & +kxc(ir,12)*rho1now(2*ir ,2,1) & 679 & +kxc(ir,13)*gradrho_gradrho1im 680 coeff_grho_up =kxc(ir,6)*rho1now(2*ir-1,1,1)+kxc(ir,8)*gradrho_gradrho1_up 681 coeff_grho_dn =kxc(ir,7)*rho1now(2*ir-1,2,1)+kxc(ir,9)*gradrho_gradrho1_dn 682 coeffim_grho_up=kxc(ir,6)*rho1now(2*ir ,1,1)+kxc(ir,8)*gradrho_gradrho1im_up 683 coeffim_grho_dn=kxc(ir,7)*rho1now(2*ir ,2,1)+kxc(ir,9)*gradrho_gradrho1im_dn 684 ! Reuse the storage in rho1now 685 rho1now(2*ir-1,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1_up(:) & 686 & +kxc(ir,10) *r1_dn(:) & 687 & +coeff_grho_up *r0_up(:) & 688 & +coeff_grho_corr*r0(:) 689 rho1now(2*ir-1,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1_dn(:) & 690 & +kxc(ir,10) *r1_up(:) & 691 & +coeff_grho_dn *r0_dn(:) & 692 & +coeff_grho_corr*r0(:) 693 rho1now(2*ir ,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1im_up(:) & 694 & +kxc(ir,10) *r1im_dn(:) & 695 & +coeffim_grho_up *r0_up(:) & 696 & +coeffim_grho_corr *r0(:) 697 rho1now(2*ir ,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1im_dn(:) & 698 & +kxc(ir,10) *r1im_up(:) & 699 & +coeffim_grho_dn *r0_dn(:) & 700 & +coeffim_grho_corr *r0(:) 701 end do 702 end if ! nspden 703 704 end if 705 706 vxc1(:,:)=zero 707 call xcpot(cplex,gprimd,ishift,use_laplacian,mpi_enreg,nfft,ngfft,ngrad,nspden,& 708 & nspgrad,qphon,depsxc=dnexcdn,rhonow=rho1now,vxc=vxc1) 709 710 !call filterpot(paral_kgb,cplex,gmet,gsqcut,nfft,ngfft,nspden,qphon,vxc1) 711 712 ABI_FREE(dnexcdn) 713 ABI_FREE(rho1now) 714 715 DBG_EXIT("COLL") 716 717 end subroutine dfpt_mkvxcgga
ABINIT/dfpt_mkvxcgga_n0met [ Functions ]
NAME
dfpt_mkvxcgga_n0met
FUNCTION
Compute the contribution to the second q-gradient of the metric perturbation that comes from gga XC potentials and depends only on ground state rho
INPUTS
beta= indicates the Cartesian direction of the metric perturbation cplex= if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX delta= indicates the Cartesian direction of the first q-gradient gamma= indicates the Cartesian direction of the second q-gradient 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=information about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT nkxc=second dimension of the kxc array nspden=number of spin-density components rho(cplex*nfft,2)=array for ground-state electron spin-density in electrons/bohr**3 (second index corresponds to spin-up and spin-down)
OUTPUT
vxc1(2*nfft,nspden)=change in exchange-correlation potential
NOTES
For the time being, a rather crude coding, to be optimized ...
Content of Kxc array:
Only works with nspden=1
===== 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)
SOURCE
1088 subroutine dfpt_mkvxcgga_n0met(beta,cplex,delta,gamma,gprimd,kxc,mpi_enreg,nfft,ngfft,& 1089 & nkxc,nspden,rhor,vxc1) 1090 1091 !Arguments ------------------------------------ 1092 !scalars 1093 integer,intent(in) :: beta,cplex,delta,gamma,nfft,nkxc,nspden 1094 type(MPI_type),intent(in) :: mpi_enreg 1095 !arrays 1096 integer,intent(in) :: ngfft(18) 1097 real(dp),intent(in) :: gprimd(3,3) 1098 real(dp),intent(in) :: kxc(nfft,nkxc) 1099 real(dp),intent(in),target :: rhor(cplex*nfft,nspden) 1100 real(dp),intent(out) :: vxc1(2*nfft,nspden) 1101 1102 !Local variables------------------------------- 1103 !scalars 1104 integer :: alpha,ii,ir,ishift,ngrad,nspgrad 1105 real(dp) :: delag,delad,delbd,delbg,deldg 1106 real(dp) :: gmodsq 1107 character(len=500) :: msg 1108 !arrays 1109 real(dp) :: r0(3) 1110 real(dp),allocatable :: dadgg(:,:),dadgtgn(:,:),gna(:,:),dadgngn_1(:,:),dadgngn_2(:,:) 1111 real(dp),allocatable :: dadgngn(:,:,:),kro_an(:,:,:),sumgrad(:,:,:) 1112 1113 ! ************************************************************************* 1114 1115 DBG_EXIT("COLL") 1116 1117 if (nkxc/=7) then 1118 msg='Wrong nkxc value for GGA in the longwave driver (optdriver=10)!' 1119 ABI_BUG(msg) 1120 end if 1121 1122 !Kronecker deltas 1123 delbd=0.0_dp; delbg=0.0_dp; deldg=0.0_dp 1124 if (beta==delta) delbd=1.0_dp 1125 if (beta==gamma) delbg=1.0_dp 1126 if (delta==gamma) deldg=1.0_dp 1127 1128 !Apply the XC kernel 1129 nspgrad=1 1130 ABI_MALLOC(dadgg,(cplex*nfft,nspgrad)) 1131 ABI_MALLOC(dadgtgn,(cplex*nfft,nspgrad)) 1132 ABI_MALLOC(gna,(cplex*nfft,nspgrad)) 1133 ABI_MALLOC(dadgngn_1,(cplex*nfft,nspgrad)) 1134 ABI_MALLOC(dadgngn_2,(cplex*nfft,nspgrad)) 1135 do ir=1,nfft 1136 r0(:)=kxc(ir,5:7) 1137 gmodsq=r0(1)**2+r0(2)**2+r0(3)**2 1138 dadgg(ir,1)=kxc(ir,4)*gmodsq*(delbd*r0(gamma)+delbg*r0(delta)) 1139 dadgtgn(ir,1)=two*kxc(ir,4)*r0(beta)*r0(delta)*r0(gamma) 1140 gna(ir,1)=(delbg*r0(delta)+delbd*r0(gamma)+two*deldg*r0(beta))*kxc(ir,2) 1141 dadgngn_1(ir,1)=delbd*kxc(ir,4)*rhor(ir,1)*r0(gamma) 1142 dadgngn_2(ir,1)=delbg*kxc(ir,4)*rhor(ir,1)*r0(delta) 1143 end do 1144 1145 !Incorporate the terms that do not need further treatment 1146 do ir=1,nfft 1147 ii=2*ir 1148 vxc1(ii-1,1)= -dadgg(ir,1)-dadgtgn(ir,1)-gna(ir,1) 1149 vxc1(ii,1)= zero 1150 end do 1151 ABI_FREE(dadgg) 1152 ABI_FREE(dadgtgn) 1153 ABI_FREE(gna) 1154 1155 !Build the last term whose gradient needs to be computed 1156 ABI_MALLOC(dadgngn,(cplex*nfft,nspgrad,3)) 1157 ABI_MALLOC(kro_an,(cplex*nfft,nspgrad,3)) 1158 ABI_MALLOC(sumgrad,(cplex*nfft,nspgrad,3)) 1159 do alpha=1,3 1160 delad=0.0_dp; delag=0.0_dp 1161 if (alpha==delta) delad=1.0_dp 1162 if (alpha==gamma) delag=1.0_dp 1163 do ir=1,nfft 1164 r0(:)=kxc(ir,5:7) 1165 dadgngn(ir,1,alpha)=(dadgngn_1(ir,1)+dadgngn_2(ir,1))*r0(alpha) 1166 kro_an(ir,1,alpha)=(delbd*delag+delbg*delad)*rhor(ir,1)*kxc(ir,2) 1167 sumgrad(ir,1,alpha)=dadgngn(ir,1,alpha)+kro_an(ir,1,alpha) 1168 end do 1169 end do 1170 1171 ABI_FREE(dadgngn_1) 1172 ABI_FREE(dadgngn_2) 1173 ABI_FREE(dadgngn) 1174 ABI_FREE(kro_an) 1175 1176 !Now the term whose sum over real-space derivatives has to be computed. 1177 !(Use the same routine as in the q-gradient of the XC kernel. It saves 1178 ! the gradient sum in the imaginary part of vxc1 and includes an additional 1179 ! two_pi factor. Need to fix this after the call.) 1180 ishift=0 ; ngrad=2 1181 call xcpotdq(sumgrad,cplex,gprimd,ishift,mpi_enreg,nfft, & 1182 & ngfft,ngrad,nspden,nspgrad,vxc1) 1183 1184 do ir=1,nfft 1185 ii=2*ir 1186 vxc1(ii-1,1)=vxc1(ii-1,1)+vxc1(ii,1)/two_pi 1187 vxc1(ii,1)=zero 1188 end do 1189 1190 1191 ABI_FREE(sumgrad) 1192 1193 end subroutine dfpt_mkvxcgga_n0met
ABINIT/dfpt_mkvxcggadq [ Functions ]
NAME
dfpt_mkvxcggadq
FUNCTION
Compute the first-order change of exchange-correlation potential in case of GGA functionals Use the q-gradient (Cartesian) of the exchange-correlation kernel.
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=information about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT nkxc=second dimension of the kxc array nspden=number of spin-density components qdirc= indicates the Cartesian direction of the q-gradient (1,2 or 3) rhor1tmp(cplex*nfft,2)=array for first-order electron spin-density in electrons/bohr**3 (second index corresponds to spin-up and spin-down)
OUTPUT
vxc1(2*nfft,nspden)=change in exchange-correlation potential
NOTES
For the time being, a rather crude coding, to be optimized ...
Content of Kxc array:
Only works with nspden=1
===== 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)
SOURCE
957 subroutine dfpt_mkvxcggadq(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,& 958 & nkxc,nspden,qdirc,rhor1,vxc1) 959 960 !Arguments ------------------------------------ 961 !scalars 962 integer,intent(in) :: cplex,nfft,nkxc,nspden,qdirc 963 type(MPI_type),intent(in) :: mpi_enreg 964 !arrays 965 integer,intent(in) :: ngfft(18) 966 real(dp),intent(in) :: gprimd(3,3) 967 real(dp),intent(in) :: kxc(nfft,nkxc) 968 real(dp),intent(in),target :: rhor1(cplex*nfft,nspden) 969 real(dp),intent(out) :: vxc1(2*nfft,nspden) 970 971 !Local variables------------------------------- 972 !scalars 973 integer :: ii,ir,ishift,ngrad,nspgrad 974 real(dp) :: gradrho_gradrho1 975 character(len=500) :: msg 976 !arrays 977 real(dp) :: qphon(3) 978 real(dp) :: r0(3),r1(3) 979 real(dp),allocatable :: ar1(:,:) 980 real(dp),allocatable :: a_gradi_r1(:,:) 981 real(dp),allocatable :: dadgradn_t1(:,:,:),dadgradn_t2(:,:) 982 real(dp),allocatable :: rho1now(:,:,:) 983 real(dp),ABI_CONTIGUOUS pointer :: rhor1_ptr(:,:) 984 985 ! ************************************************************************* 986 987 DBG_EXIT("COLL") 988 989 if (nkxc/=7) then 990 msg='Wrong nkxc value for GGA in the longwave driver (optdriver=10)!' 991 ABI_BUG(msg) 992 end if 993 994 !Compute the gradients of the first-order density 995 !rho1now(:,:,1) contains the first-order density, and 996 !rho1now(:,:,2:4) contains the gradients of the first-order density 997 ishift=0 ; ngrad=2 998 qphon(:)=zero 999 rhor1_ptr => rhor1 1000 ABI_MALLOC(rho1now,(cplex*nfft,nspden,ngrad*ngrad)) 1001 call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,qphon,rhor1_ptr,rho1now) 1002 1003 !Apply the XC kernel 1004 nspgrad=1 1005 ABI_MALLOC(ar1,(cplex*nfft,nspgrad)) 1006 ABI_MALLOC(a_gradi_r1,(cplex*nfft,nspgrad)) 1007 ABI_MALLOC(dadgradn_t1,(cplex*nfft,nspgrad,3)) 1008 ABI_MALLOC(dadgradn_t2,(cplex*nfft,nspgrad)) 1009 do ir=1,nfft 1010 r0(:)=kxc(ir,5:7); r1(:)=rho1now(ir,1,2:4) 1011 gradrho_gradrho1=dot_product(r0,r1) 1012 ar1(ir,1)=kxc(ir,2)*rho1now(ir,1,1) 1013 a_gradi_r1(ir,1)=kxc(ir,2)*r1(qdirc) 1014 dadgradn_t2(ir,1)=kxc(ir,4)*gradrho_gradrho1*r0(qdirc) 1015 dadgradn_t1(ir,1,:)=kxc(ir,4)*r0(:)*r0(qdirc)*rho1now(ir,1,1) 1016 end do 1017 do ii=1,3 1018 if (ii==qdirc) dadgradn_t1(:,1,ii)=dadgradn_t1(:,1,ii)+ar1(:,1) 1019 end do 1020 1021 !Incorporate the terms that do not need further treatment 1022 !(a -i factor is applied here) 1023 do ir=1,nfft 1024 ii=2*ir 1025 vxc1(ii-1,1)=zero 1026 vxc1(ii,1)= -a_gradi_r1(ir,1) -dadgradn_t2(ir,1) 1027 end do 1028 ABI_FREE(rho1now) 1029 ABI_FREE(a_gradi_r1) 1030 ABI_FREE(dadgradn_t2) 1031 ABI_FREE(ar1) 1032 1033 !Now the term whose sum over real-space derivatives has to be computed 1034 call xcpotdq(dadgradn_t1,cplex,gprimd,ishift,mpi_enreg,nfft, & 1035 & ngfft,ngrad,nspden,nspgrad,vxc1) 1036 1037 ABI_FREE(dadgradn_t1) 1038 1039 end subroutine dfpt_mkvxcggadq
ABINIT/m_dfpt_mkvxc [ Modules ]
NAME
m_dfpt_mkvxc
FUNCTION
COPYRIGHT
Copyright (C) 2001-2024 ABINIT group (XG, DRH, FR, EB, SPr) 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_dfpt_mkvxc 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 use m_xc_noncoll 28 29 use defs_abitypes, only : MPI_type 30 use m_time, only : timab 31 use m_symtk, only : matr3inv 32 use m_xctk, only : xcden, xcpot, xcpotdq 33 34 implicit none 35 36 private