TABLE OF CONTENTS
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.
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 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 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)
SIDE EFFECTS
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
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_dyxc1,dfpt_mkvxc_noncoll,dfpt_nstdy,dfpt_nstpaw,dfpt_rhotov dfptnl_loop,m_kxc,nres2vres
CHILDREN
dfpt_mkvxcgga,matr3inv,timab
SOURCE
96 #if defined HAVE_CONFIG_H 97 #include "config.h" 98 #endif 99 100 #include "abi_common.h" 101 102 subroutine dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,nhat1dim,nhat1gr,nhat1grdim,& 103 & nkxc,nspden,n3xccc,option,paral_kgb,qphon,rhor1,rprimd,usexcnhat,vxc1,xccc3d1) 104 105 use defs_basis 106 use defs_abitypes 107 use m_errors 108 use m_profiling_abi 109 110 !This section has been created automatically by the script Abilint (TD). 111 !Do not modify the following lines by hand. 112 #undef ABI_FUNC 113 #define ABI_FUNC 'dfpt_mkvxc' 114 use interfaces_18_timing 115 use interfaces_32_util 116 use interfaces_56_xc, except_this_one => dfpt_mkvxc 117 !End of the abilint section 118 119 implicit none 120 121 !Arguments ------------------------------------ 122 !scalars 123 integer,intent(in) :: cplex,ixc,n3xccc,nfft,nhat1dim,nhat1grdim 124 integer,intent(in) :: nkxc,nspden,option,paral_kgb,usexcnhat 125 type(MPI_type),intent(in) :: mpi_enreg 126 !arrays 127 integer,intent(in) :: ngfft(18) 128 real(dp),intent(in),target :: nhat1(cplex*nfft,nspden*nhat1dim) 129 real(dp),intent(in),target :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim) 130 real(dp),intent(in) :: kxc(nfft,nkxc),qphon(3) 131 real(dp),intent(in),target :: rhor1(cplex*nfft,nspden) 132 real(dp),intent(in) :: rprimd(3,3),xccc3d1(cplex*n3xccc) 133 real(dp),intent(out) :: vxc1(cplex*nfft,nspden) 134 135 !Local variables------------------------------- 136 !scalars 137 integer :: ii,ir,ispden,nhat1dim_,nhat1rgdim_ 138 real(dp) :: rho1_dn,rho1_up,rho1im_dn,rho1im_up,rho1re_dn,rho1re_up 139 real(dp) :: spin_scale 140 !arrays 141 real(dp) :: gprimd(3,3),tsec(2) 142 real(dp), ABI_CONTIGUOUS pointer :: nhat1_(:,:),nhat1gr_(:,:,:),rhor1_(:,:) 143 144 ! ************************************************************************* 145 146 DBG_ENTER("COLL") 147 148 call timab(181,1,tsec) 149 150 if(nspden/=1 .and. nspden/=2) then 151 MSG_BUG('For nspden==4 please use dfpt_mkvxc_noncoll!') 152 end if 153 154 !Special case: no XC applied 155 if (ixc==0.or.nkxc==0) then 156 MSG_WARNING('Note that no xc is applied (ixc=0)') 157 vxc1=zero 158 return 159 end if 160 161 !Treat first LDA 162 if(nkxc==1.or.nkxc==3)then 163 164 ! PAW: eventually substract compensation density 165 if (option/=0) then 166 if (usexcnhat==0.and.nhat1dim==1) then 167 ABI_ALLOCATE(rhor1_,(cplex*nfft,nspden)) 168 rhor1_(:,:)=rhor1(:,:)-nhat1(:,:) 169 else 170 rhor1_ => rhor1 171 end if 172 end if 173 174 ! Case without non-linear core correction 175 if(n3xccc==0 .or. option==2)then 176 177 if(option==0)then ! No straight XC to compute 178 179 vxc1(:,:)=zero 180 181 else ! XC, without non-linear XC correction 182 183 ! Non-spin-polarized 184 if(nspden==1)then 185 if(cplex==1)then 186 do ir=1,nfft 187 vxc1(ir,1)=kxc(ir,1)*rhor1_(ir,1) 188 end do 189 else 190 do ir=1,nfft 191 vxc1(2*ir-1,1)=kxc(ir,1)*rhor1_(2*ir-1,1) 192 vxc1(2*ir ,1)=kxc(ir,1)*rhor1_(2*ir ,1) 193 end do 194 end if ! cplex==1 195 196 ! Spin-polarized 197 else 198 if(cplex==1)then 199 do ir=1,nfft 200 rho1_dn=rhor1_(ir,1)-rhor1_(ir,2) 201 vxc1(ir,1)=kxc(ir,1)*rhor1_(ir,2)+kxc(ir,2)*rho1_dn 202 vxc1(ir,2)=kxc(ir,2)*rhor1_(ir,2)+kxc(ir,3)*rho1_dn 203 end do 204 else 205 do ir=1,nfft 206 rho1re_dn=rhor1_(2*ir-1,1)-rhor1_(2*ir-1,2) 207 rho1im_dn=rhor1_(2*ir ,1)-rhor1_(2*ir ,2) 208 vxc1(2*ir-1,1)=kxc(ir,1)*rhor1_(2*ir-1,2)+kxc(ir,2)*rho1re_dn 209 vxc1(2*ir ,1)=kxc(ir,1)*rhor1_(2*ir ,2)+kxc(ir,2)*rho1im_dn 210 vxc1(2*ir-1,2)=kxc(ir,2)*rhor1_(2*ir-1,2)+kxc(ir,3)*rho1re_dn 211 vxc1(2*ir ,2)=kxc(ir,2)*rhor1_(2*ir ,2)+kxc(ir,3)*rho1im_dn 212 end do 213 end if ! cplex==1 214 end if ! nspden==1 215 216 end if ! option==0 217 218 ! Treat case with non-linear core correction 219 else 220 221 if(option==0)then 222 223 if(nspden==1)then 224 if(cplex==1)then 225 do ir=1,nfft 226 vxc1(ir,1)=kxc(ir,1)*xccc3d1(ir) 227 end do 228 else 229 do ir=1,nfft 230 vxc1(2*ir-1,1)=kxc(ir,1)*xccc3d1(2*ir-1) 231 vxc1(2*ir ,1)=kxc(ir,1)*xccc3d1(2*ir ) 232 end do 233 end if ! cplex==1 234 else 235 if(cplex==1)then 236 do ir=1,nfft 237 vxc1(ir,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(ir)*half 238 vxc1(ir,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(ir)*half 239 end do 240 else 241 do ir=1,nfft 242 vxc1(2*ir-1,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(2*ir-1)*half 243 vxc1(2*ir ,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(2*ir )*half 244 vxc1(2*ir-1,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(2*ir-1)*half 245 vxc1(2*ir ,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(2*ir )*half 246 end do 247 end if ! cplex==1 248 end if ! nspden==1 249 250 else ! option/=0 251 252 if(nspden==1)then 253 if(cplex==1)then 254 do ir=1,nfft 255 vxc1(ir,1)=kxc(ir,1)*(rhor1_(ir,1)+xccc3d1(ir)) 256 end do 257 else 258 do ir=1,nfft 259 vxc1(2*ir-1,1)=kxc(ir,1)*(rhor1_(2*ir-1,1)+xccc3d1(2*ir-1)) 260 vxc1(2*ir ,1)=kxc(ir,1)*(rhor1_(2*ir ,1)+xccc3d1(2*ir )) 261 end do 262 end if ! cplex==1 263 else 264 if(cplex==1)then 265 do ir=1,nfft 266 rho1_dn=rhor1_(ir,1)-rhor1_(ir,2) + xccc3d1(ir)*half 267 rho1_up=rhor1_(ir,2) + xccc3d1(ir)*half 268 vxc1(ir,1)=kxc(ir,1)*rho1_up+kxc(ir,2)*rho1_dn 269 vxc1(ir,2)=kxc(ir,2)*rho1_up+kxc(ir,3)*rho1_dn 270 end do 271 else 272 do ir=1,nfft 273 rho1re_dn=rhor1_(2*ir-1,1)-rhor1_(2*ir-1,2) + xccc3d1(2*ir-1)*half 274 rho1im_dn=rhor1_(2*ir ,1)-rhor1_(2*ir ,2) + xccc3d1(2*ir )*half 275 rho1re_up=rhor1_(2*ir-1,2) + xccc3d1(2*ir-1)*half 276 rho1im_up=rhor1_(2*ir ,2) + xccc3d1(2*ir )*half 277 vxc1(2*ir-1,1)=kxc(ir,1)*rho1re_up+kxc(ir,2)*rho1re_dn 278 vxc1(2*ir ,1)=kxc(ir,1)*rho1im_up+kxc(ir,2)*rho1im_dn 279 vxc1(2*ir-1,2)=kxc(ir,2)*rho1re_up+kxc(ir,3)*rho1re_dn 280 vxc1(2*ir ,2)=kxc(ir,2)*rho1im_up+kxc(ir,3)*rho1im_dn 281 end do 282 end if ! cplex==1 283 end if ! nspden==1 284 285 end if ! option==0 286 287 end if ! n3xccc==0 288 289 if (option/=0.and.usexcnhat==0.and.nhat1dim==1) then 290 ABI_DEALLOCATE(rhor1_) 291 end if 292 293 ! Treat GGA 294 else if (nkxc==7.or.nkxc==19) then 295 296 ! Transfer the data to spin-polarized storage 297 298 ! Treat the density change 299 ABI_ALLOCATE(rhor1_,(cplex*nfft,nspden)) 300 if (option==1 .or. option==2) then 301 if (nspden==1) then 302 do ir=1,cplex*nfft 303 rhor1_(ir,1)=rhor1(ir,1) 304 end do 305 else 306 do ir=1,cplex*nfft 307 rho1_dn=rhor1(ir,1)-rhor1(ir,2) 308 rhor1_(ir,1)=rhor1(ir,2) 309 rhor1_(ir,2)=rho1_dn 310 end do 311 end if 312 else 313 do ispden=1,nspden 314 do ir=1,cplex*nfft 315 rhor1_(ir,ispden)=zero 316 end do 317 end do 318 end if 319 if( (option==0 .or. option==1) .and. n3xccc/=0)then 320 spin_scale=one;if (nspden==2) spin_scale=half 321 do ispden=1,nspden 322 do ir=1,cplex*nfft 323 rhor1_(ir,ispden)=rhor1_(ir,ispden)+xccc3d1(ir)*spin_scale 324 end do 325 end do 326 end if 327 328 ! PAW: treat also compensation density (and gradients) 329 nhat1dim_=nhat1dim ; nhat1rgdim_=nhat1grdim 330 if (option/=0.and.nhat1dim==1.and.nspden==2) then 331 ABI_ALLOCATE(nhat1_,(cplex*nfft,nspden)) 332 do ir=1,cplex*nfft 333 rho1_dn=nhat1(ir,1)-nhat1(ir,2) 334 nhat1_(ir,1)=nhat1(ir,2) 335 nhat1_(ir,2)=rho1_dn 336 end do 337 else if (option==0) then 338 ABI_ALLOCATE(nhat1_,(0,0)) 339 nhat1dim_=0 340 else 341 nhat1_ => nhat1 342 end if 343 if (option/=0.and.nhat1grdim==1.and.nspden==2) then 344 ABI_ALLOCATE(nhat1gr_,(cplex*nfft,nspden,3)) 345 do ii=1,3 346 do ir=1,cplex*nfft 347 rho1_dn=nhat1gr(ir,1,ii)-nhat1gr(ir,2,ii) 348 nhat1gr_(ir,1,ii)=nhat1gr(ir,2,ii) 349 nhat1gr_(ir,2,ii)=rho1_dn 350 end do 351 end do 352 else if (option==0) then 353 ABI_ALLOCATE(nhat1gr_,(0,0,0)) 354 nhat1rgdim_=0 355 else 356 nhat1gr_ => nhat1gr 357 end if 358 359 call matr3inv(rprimd,gprimd) 360 361 call dfpt_mkvxcgga(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat1_,nhat1dim_,& 362 & nhat1gr_,nhat1rgdim_,nkxc,nspden,paral_kgb,qphon,rhor1_,usexcnhat,vxc1) 363 364 ABI_DEALLOCATE(rhor1_) 365 if ((option==0).or.(nhat1dim==1.and.nspden==2)) then 366 ABI_DEALLOCATE(nhat1_) 367 end if 368 if ((option==0).or.(nhat1grdim==1.and.nspden==2)) then 369 ABI_DEALLOCATE(nhat1gr_) 370 end if 371 372 else 373 MSG_BUG('Invalid nkxc!') 374 375 end if ! LDA or GGA 376 377 call timab(181,2,tsec) 378 379 DBG_EXIT("COLL") 380 381 end subroutine dfpt_mkvxc