TABLE OF CONTENTS
ABINIT/dfpt_mkvxcstr [ Functions ]
NAME
dfpt_mkvxcstr
FUNCTION
Compute the first-order change of exchange-correlation potential due to strain: 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 (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 idir=direction of the current perturbation ipert=type of the perturbation kxc(nfft,nkxc)=exchange and correlation kernel (see rhotoxc.f) mpi_enreg=information about MPI parallelization natom=number of atoms in cell. 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- compensation density nhat1(cplex*nfft,2nspden*usepaw)= -PAW only- 1st-order compensation density 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 strain-derivative frozen-wavefunction charge and the XC core-correction, if 1, treat both density change and XC core correction if 2, like 0 but multiply gradient strain derivative term by 2.0 for GGA. qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2). rhor(nfft,nspden)=array for GS electron density in electrons/bohr**3. rhor1(cplex*nfft,nspden)=array for electron density in electrons/bohr**3. rprimd(3,3)=dimensional primitive translations in real space (bohr) usepaw= 0 for non paw calculation; =1 for paw calculation usexcnhat= -PAW only- flag controling use of compensation density 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)
PARENTS
dfpt_eltfrxc,dfpt_nselt,dfpt_nstpaw,dfpt_rhotov
CHILDREN
dfpt_mkvxcstrgga,matr3inv,timab
SOURCE
58 #if defined HAVE_CONFIG_H 59 #include "config.h" 60 #endif 61 62 #include "abi_common.h" 63 64 65 subroutine dfpt_mkvxcstr(cplex,idir,ipert,kxc,mpi_enreg,natom,nfft,ngfft,nhat,nhat1,& 66 & nkxc,nspden,n3xccc,option,paral_kgb,qphon,rhor,rhor1,rprimd,usepaw,usexcnhat,vxc1,xccc3d1) 67 68 use defs_basis 69 use defs_abitypes 70 use m_errors 71 use m_profiling_abi 72 73 !This section has been created automatically by the script Abilint (TD). 74 !Do not modify the following lines by hand. 75 #undef ABI_FUNC 76 #define ABI_FUNC 'dfpt_mkvxcstr' 77 use interfaces_18_timing 78 use interfaces_32_util 79 use interfaces_72_response, except_this_one => dfpt_mkvxcstr 80 !End of the abilint section 81 82 implicit none 83 84 !Arguments ------------------------------------ 85 !scalars 86 integer,intent(in) :: cplex,idir,ipert,n3xccc,natom,nfft,nkxc,nspden,option 87 integer,intent(in) :: paral_kgb,usepaw,usexcnhat 88 type(MPI_type),intent(in) :: mpi_enreg 89 !arrays 90 integer,intent(in) :: ngfft(18) 91 real(dp),target,intent(in) :: nhat(nfft,nspden) 92 real(dp),target,intent(in) :: nhat1(cplex*nfft,nspden) 93 real(dp),intent(in) :: kxc(nfft,nkxc),qphon(3) 94 real(dp),target,intent(in) :: rhor(nfft,nspden),rhor1(cplex*nfft,nspden) 95 real(dp),intent(in) :: rprimd(3,3) 96 real(dp),intent(in) :: xccc3d1(cplex*n3xccc) 97 real(dp),intent(out) :: vxc1(cplex*nfft,nspden) 98 99 !Local variables------------------------------- 100 !scalars 101 integer :: ii,ir,istr 102 real(dp) :: rho1_dn,rho1_up,spin_scale,str_scale 103 character(len=500) :: message 104 !arrays 105 real(dp) :: gprimd(3,3),tsec(2) 106 real(dp),allocatable :: rhor1tmp(:,:),rhowk1(:,:) 107 real(dp),pointer :: rhor_(:,:),rhor1_(:,:) 108 109 ! ************************************************************************* 110 111 call timab(181,1,tsec) 112 113 if(nspden/=1 .and. nspden/=2) then 114 message = ' dfpt_mkvxc, Only for nspden==1 and 2.' 115 MSG_BUG(message) 116 end if 117 118 if (usepaw==1.and.usexcnhat==0) then 119 ABI_ALLOCATE(rhor_,(nfft,nspden)) 120 rhor_(:,:)=rhor(:,:)-nhat(:,:) 121 else 122 rhor_ => rhor 123 end if 124 125 if (usepaw==1.and.usexcnhat==0.and.option==1) then 126 ABI_ALLOCATE(rhor1_,(nfft,nspden)) 127 rhor1_(:,:)=rhor1(:,:)-nhat1(:,:) 128 else 129 rhor1_ => rhor1 130 end if 131 132 133 !Inhomogeneous term for diagonal strain 134 ABI_ALLOCATE(rhowk1,(nfft,nspden)) 135 if(option==0 .or. option==2) then 136 if(ipert==natom+3) then 137 rhowk1(:,:)=-rhor_(:,:) 138 else 139 rhowk1(:,:)=zero 140 end if 141 else if(option==1) then 142 if(ipert==natom+3) then 143 rhowk1(:,:)=rhor1_(:,:)-rhor_(:,:) 144 else 145 rhowk1(:,:)=rhor1_(:,:) 146 end if 147 end if 148 149 !Treat first LDA 150 if(nkxc==1.or.nkxc==3)then 151 152 ! Case without non-linear core correction 153 if(n3xccc==0)then 154 155 ! Non-spin-polarized 156 if(nspden==1)then 157 do ir=1,nfft 158 vxc1(ir,1)=kxc(ir,1)*rhowk1(ir,1) 159 end do 160 161 ! Spin-polarized 162 else 163 do ir=1,nfft 164 rho1_dn=rhowk1(ir,1)-rhowk1(ir,2) 165 vxc1(ir,1)=kxc(ir,1)*rhowk1(ir,2)+kxc(ir,2)*rho1_dn 166 vxc1(ir,2)=kxc(ir,2)*rhowk1(ir,2)+kxc(ir,3)*rho1_dn 167 end do 168 end if ! nspden==1 169 170 ! Treat case with non-linear core correction 171 else 172 if(nspden==1)then 173 do ir=1,nfft 174 vxc1(ir,1)=kxc(ir,1)*(rhowk1(ir,1)+xccc3d1(ir)) 175 end do 176 else 177 do ir=1,nfft 178 rho1_dn=rhowk1(ir,1)-rhowk1(ir,2) + xccc3d1(ir)*half 179 rho1_up=rhowk1(ir,2) + xccc3d1(ir)*half 180 vxc1(ir,1)=kxc(ir,1)*rho1_up+kxc(ir,2)*rho1_dn 181 vxc1(ir,2)=kxc(ir,2)*rho1_up+kxc(ir,3)*rho1_dn 182 end do 183 end if ! nspden==1 184 185 end if ! n3xccc==0 186 187 ! Treat GGA 188 else if (nkxc==7.or.nkxc==19) then 189 190 ! Generates gprimd and its strain derivative 191 ! Note that unlike the implicitly symmetric metric tensor strain 192 ! derivatives, we must explicltly symmetrize the strain derivative 193 ! here. 194 call matr3inv(rprimd,gprimd) 195 istr=idir + 3*(ipert-natom-3) 196 if(istr<1 .or. istr>6)then 197 write(message, '(a,i10,a,a,a)' )& 198 & 'Input dir gives istr=',istr,' not allowed.',ch10,& 199 & 'Possible values are 1,2,3,4,5,6 only.' 200 MSG_BUG(message) 201 end if 202 203 ! Rescalling needed for use in dfpt_eltfrxc for elastic tensor (not internal strain tensor). 204 str_scale=one;if(option==2) str_scale=two 205 206 ! FTransfer the data to spin-polarized storage 207 ABI_ALLOCATE(rhor1tmp,(cplex*nfft,nspden)) 208 if(nspden==1)then 209 do ir=1,cplex*nfft 210 rhor1tmp(ir,1)=rhowk1(ir,1) 211 end do 212 else 213 do ir=1,cplex*nfft 214 rho1_dn=rhowk1(ir,1)-rhowk1(ir,2) 215 rhor1tmp(ir,1)=rhowk1(ir,2) 216 rhor1tmp(ir,2)=rho1_dn 217 end do 218 end if ! nspden==1 219 if(n3xccc/=0)then 220 spin_scale=one;if (nspden==2) spin_scale=half 221 do ii=1,nspden 222 do ir=1,cplex*nfft 223 rhor1tmp(ir,ii)=rhor1tmp(ir,ii)+xccc3d1(ir)*spin_scale 224 end do 225 end do 226 end if 227 228 call dfpt_mkvxcstrgga(cplex,gprimd,istr,kxc,mpi_enreg,nfft,ngfft,nkxc,& 229 & nspden,paral_kgb,qphon,rhor1tmp,str_scale,vxc1) 230 ABI_DEALLOCATE(rhor1tmp) 231 232 else 233 MSG_BUG('Invalid nkxc!') 234 235 end if ! LDA or GGA 236 237 if (usepaw==1.and.usexcnhat==0) then 238 ABI_DEALLOCATE(rhor_) 239 end if 240 if (usepaw==1.and.usexcnhat==0.and.option==1) then 241 ABI_DEALLOCATE(rhor1_) 242 end if 243 244 ABI_DEALLOCATE(rhowk1) 245 246 call timab(181,2,tsec) 247 248 end subroutine dfpt_mkvxcstr