TABLE OF CONTENTS
ABINIT/dfpt_mkvxcstrgga [ Functions ]
NAME
dfpt_mkvxcstrgga
FUNCTION
Compute the first-order change of exchange-correlation potential for the strain perturbation in case of GGA functionals Use the exchange-correlation kernel.
COPYRIGHT
Copyright (C) 2001-2018 ABINIT group (DRH, XG) 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 gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1) istr=index of the strain perturbation (1..6) kxc(nfft,nkxc)=exchange and correlation kernel (see rhotoxc.f) 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 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) str_scale=scaling factor for gradient operator strain-derivative (1. or 2.)
OUTPUT
vxc1(cplex*nfft,nspden)=change in exchange-correlation potential
NOTES
Closely related to dfpt_mkvxcgga. 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_mkvxcstr
CHILDREN
xcden,xcpot
SOURCE
78 #if defined HAVE_CONFIG_H 79 #include "config.h" 80 #endif 81 82 #include "abi_common.h" 83 84 85 subroutine dfpt_mkvxcstrgga(cplex,gprimd,istr,kxc,mpi_enreg,nfft,ngfft,& 86 & nkxc,nspden,paral_kgb,qphon,rhor1tmp,str_scale,vxc1) 87 88 use defs_basis 89 use defs_abitypes 90 use m_profiling_abi 91 use m_errors 92 93 !This section has been created automatically by the script Abilint (TD). 94 !Do not modify the following lines by hand. 95 #undef ABI_FUNC 96 #define ABI_FUNC 'dfpt_mkvxcstrgga' 97 use interfaces_56_xc 98 !End of the abilint section 99 100 implicit none 101 102 !Arguments ------------------------------------ 103 !scalars 104 integer,intent(in) :: cplex,istr,nfft,nkxc,nspden,paral_kgb 105 real(dp) :: str_scale 106 type(MPI_type),intent(in) :: mpi_enreg 107 !arrays 108 integer,intent(in) :: ngfft(18) 109 real(dp),intent(in) :: gprimd(3,3),kxc(nfft,nkxc) 110 real(dp),intent(in) :: qphon(3),rhor1tmp(cplex*nfft,2) 111 real(dp),intent(out) :: vxc1(cplex*nfft,nspden) 112 !Local variables------------------------------- 113 !scalars 114 integer :: ii,ir,ishift,ispden,mgga,ngrad,nspgrad 115 real(dp) :: coeff_grho,coeff_grho_corr,coeff_grho_dn,coeff_grho_up 116 real(dp) :: gradrho_gradrho1,gradrho_gradrho1_dn,gradrho_gradrho1_up 117 character(len=500) :: msg 118 !arrays 119 real(dp) :: r0(3),r0_dn(3),r0_up(3),r1(3),r1_dn(3),r1_up(3) 120 real(dp),allocatable :: dnexcdn(:,:),rho1now(:,:,:),rhodgnow(:,:,:) 121 122 ! ************************************************************************* 123 124 DBG_ENTER("COLL") 125 126 if (nkxc/=12*min(nspden,2)-5) then 127 msg='Wrong nkxc value for GGA!' 128 MSG_BUG(msg) 129 end if 130 if (nspden>2) then 131 msg='Not compatible with non-collinear magnetism!' 132 MSG_ERROR(msg) 133 end if 134 135 !metaGGA contributions are not taken into account here 136 mgga=0 137 138 !if you uncomment the following line, you will have to modify 139 !the original function call to pass in gmet and gsqcut 140 !call filterpot(cplex,gmet,gsqcut,nfft,ngfft,2,qphon,rhor1tmp) 141 142 !Compute the gradients of the first-order density 143 !rho1now(:,:,1) contains the first-order density, and 144 !rho1now(:,:,2:4) contains the gradients of the first-order density 145 ishift=0 ; ngrad=2 146 ABI_ALLOCATE(rho1now,(cplex*nfft,nspden,ngrad*ngrad)) 147 call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,paral_kgb,qphon,rhor1tmp,rho1now) 148 149 !Calculate the 1st-order contribution to grad(n) from the strain derivative 150 ! acting on the gradient operator acting on the GS charge density, 151 !Simply use the following formula: 152 ! (dGprim/ds_alpha_beta)(i,j) = -half.( delta_alpha,i Gprim(beta,j) + delta_beta,i Gprim(alpha,j) ) 153 !To finally get: 154 ! (nabla)^(alpha,beta)_i[n] = -half ( delta_alpha,i nabla_beta[n] + delta_beta,i nabla_alpha[n] ) 155 ABI_ALLOCATE(rhodgnow,(cplex*nfft,nspden,3)) 156 rhodgnow(1:nfft,1:nspden,1:3)=zero 157 if (nspden==1) then 158 if (istr==1) rhodgnow(1:nfft,1,1)=- kxc(1:nfft,5) 159 if (istr==2) rhodgnow(1:nfft,1,2)=- kxc(1:nfft,6) 160 if (istr==3) rhodgnow(1:nfft,1,3)=- kxc(1:nfft,7) 161 if (istr==4) rhodgnow(1:nfft,1,2)=-half*kxc(1:nfft,7) 162 if (istr==4) rhodgnow(1:nfft,1,3)=-half*kxc(1:nfft,6) 163 if (istr==5) rhodgnow(1:nfft,1,1)=-half*kxc(1:nfft,7) 164 if (istr==5) rhodgnow(1:nfft,1,3)=-half*kxc(1:nfft,5) 165 if (istr==6) rhodgnow(1:nfft,1,1)=-half*kxc(1:nfft,6) 166 if (istr==6) rhodgnow(1:nfft,1,2)=-half*kxc(1:nfft,5) 167 else 168 if (istr==1) rhodgnow(1:nfft,1,1)=- kxc(1:nfft,15) 169 if (istr==2) rhodgnow(1:nfft,1,2)=- kxc(1:nfft,17) 170 if (istr==3) rhodgnow(1:nfft,1,3)=- kxc(1:nfft,19) 171 if (istr==4) rhodgnow(1:nfft,1,2)=-half*kxc(1:nfft,19) 172 if (istr==4) rhodgnow(1:nfft,1,3)=-half*kxc(1:nfft,17) 173 if (istr==5) rhodgnow(1:nfft,1,1)=-half*kxc(1:nfft,19) 174 if (istr==5) rhodgnow(1:nfft,1,3)=-half*kxc(1:nfft,15) 175 if (istr==6) rhodgnow(1:nfft,1,1)=-half*kxc(1:nfft,17) 176 if (istr==6) rhodgnow(1:nfft,1,2)=-half*kxc(1:nfft,15) 177 if (istr==1) rhodgnow(1:nfft,2,1)=- (kxc(1:nfft,14)-kxc(1:nfft,15)) 178 if (istr==2) rhodgnow(1:nfft,2,2)=- (kxc(1:nfft,16)-kxc(1:nfft,17)) 179 if (istr==3) rhodgnow(1:nfft,2,3)=- (kxc(1:nfft,18)-kxc(1:nfft,19)) 180 if (istr==4) rhodgnow(1:nfft,2,2)=-half*(kxc(1:nfft,18)-kxc(1:nfft,19)) 181 if (istr==4) rhodgnow(1:nfft,2,3)=-half*(kxc(1:nfft,16)-kxc(1:nfft,17)) 182 if (istr==5) rhodgnow(1:nfft,2,1)=-half*(kxc(1:nfft,18)-kxc(1:nfft,19)) 183 if (istr==5) rhodgnow(1:nfft,2,3)=-half*(kxc(1:nfft,14)-kxc(1:nfft,15)) 184 if (istr==6) rhodgnow(1:nfft,2,1)=-half*(kxc(1:nfft,16)-kxc(1:nfft,17)) 185 if (istr==6) rhodgnow(1:nfft,2,2)=-half*(kxc(1:nfft,14)-kxc(1:nfft,15)) 186 end if 187 188 !Add to the gradients of the first-order density 189 do ii=1,3 190 do ispden=1,nspden 191 rhodgnow(1:nfft,ispden,ii)=str_scale*rhodgnow(1:nfft,ispden,ii) 192 rho1now(1:nfft,ispden,1+ii)=rho1now(1:nfft,ispden,1+ii)+rhodgnow(1:nfft,ispden,ii) 193 end do 194 end do 195 196 !rho1now(:,:,1) contains the 1st-order density, and rho1now(:,:,2:4) contains the grads of the 1st-order density 197 198 !Apply the XC kernel 199 nspgrad=2; if (nspden==2) nspgrad=5 200 ABI_ALLOCATE(dnexcdn,(cplex*nfft,nspgrad)) 201 202 !== Non polarized 203 if (nspden==1) then 204 do ir=1,nfft 205 r0(:)=kxc(ir,5:7) ; r1(:)=rho1now(ir,1,2:4) 206 gradrho_gradrho1=dot_product(r0,r1) 207 dnexcdn(ir,1)=kxc(ir,1)*rho1now(ir,1,1) + kxc(ir,3)*gradrho_gradrho1 208 coeff_grho=kxc(ir,3)*rho1now(ir,1,1) + kxc(ir,4)*gradrho_gradrho1 209 ! Grad strain derivative contribution enters the following term with a 210 ! factor of two compared to above terms, so add it again. 211 r1(:)=r1(:)+rhodgnow(ir,1,1:3) 212 ! Reuse the storage in rho1now 213 rho1now(ir,1,2:4)=r1(:)*kxc(ir,2)+r0(:)*coeff_grho 214 end do 215 216 !== Spin-polarized 217 else ! nspden==2 218 do ir=1,nfft 219 do ii=1,3 ! grad of spin-up ans spin_dn GS rho 220 r0_up(ii)=kxc(ir,13+2*ii);r0_dn(ii)=kxc(ir,12+2*ii)-kxc(ir,13+2*ii) 221 end do 222 r0(:)=r0_up(:)+r0_dn(:) ! grad of GS rho 223 r1_up(:)=rho1now(ir,1,2:4) ! grad of spin-up rho1 224 r1_dn(:)=rho1now(ir,2,2:4) ! grad of spin-down rho1 225 r1(:)=r1_up(:)+r1_dn(:) ! grad of GS rho1 226 gradrho_gradrho1_up=dot_product(r0_up,r1_up) 227 gradrho_gradrho1_dn=dot_product(r0_dn,r1_dn) 228 gradrho_gradrho1 =dot_product(r0,r1) 229 230 dnexcdn(ir,1)=kxc(ir, 1)*rho1now(ir,1,1) & 231 & +kxc(ir, 2)*rho1now(ir,2,1) & 232 & +kxc(ir, 6)*gradrho_gradrho1_up & 233 & +kxc(ir,11)*gradrho_gradrho1 234 dnexcdn(ir,2)=kxc(ir, 3)*rho1now(ir,2,1) & 235 & +kxc(ir, 2)*rho1now(ir,1,1) & 236 & +kxc(ir, 7)*gradrho_gradrho1_dn & 237 & +kxc(ir,12)*gradrho_gradrho1 238 coeff_grho_corr=kxc(ir,11)*rho1now(ir,1,1) & 239 & +kxc(ir,12)*rho1now(ir,2,1) & 240 & +kxc(ir,13)*gradrho_gradrho1 241 coeff_grho_up=kxc(ir,6)*rho1now(ir,1,1)+kxc(ir,8)*gradrho_gradrho1_up 242 coeff_grho_dn=kxc(ir,7)*rho1now(ir,2,1)+kxc(ir,9)*gradrho_gradrho1_dn 243 244 ! Grad strain derivative contribution enters the following term with a 245 ! factor of two compared to above terms, so add it again. 246 r1_up(:)=r1_up(:)+rhodgnow(ir,1,1:3) 247 r1_dn(:)=r1_dn(:)+rhodgnow(ir,2,1:3) 248 249 ! Reuse the storage in rho1now 250 rho1now(ir,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1_up(:) & 251 & +kxc(ir,10) *r1_dn(:) & 252 & +coeff_grho_up *r0_up(:) & 253 & +coeff_grho_corr *r0(:) 254 rho1now(ir,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1_dn(:) & 255 & +kxc(ir,10) *r1_up(:) & 256 & +coeff_grho_dn *r0_dn(:) & 257 & +coeff_grho_corr *r0(:) 258 end do 259 260 end if ! nspden 261 ABI_DEALLOCATE(rhodgnow) 262 263 vxc1(:,:)=zero 264 call xcpot(cplex,dnexcdn,gprimd,ishift,mgga,mpi_enreg,nfft,ngfft,ngrad,nspden,& 265 & nspgrad,paral_kgb,qphon,rho1now,vxc1) 266 267 !if you uncomment the following line, you will have to modify 268 !the original function call to pass in gmet and gsqcut 269 !call filterpot(cplex,gmet,gsqcut,nfft,ngfft,nspden,qphon,vxc1) 270 271 ABI_DEALLOCATE(dnexcdn) 272 ABI_DEALLOCATE(rho1now) 273 274 end subroutine dfpt_mkvxcstrgga