TABLE OF CONTENTS
ABINIT/dfpt_rhotov [ Functions ]
NAME
dfpt_rhotov
FUNCTION
This routine is called to compute, from a given 1st-order total density - the trial (local) 1st-order potential and/or the residual potential, - some contributions to the 2nd-order energy
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (XG, DRH, MT, 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cplex: if 1, real space 1-order WF on FFT grid are REAL; if 2, COMPLEX gsqcut=cutoff on (k+G)^2 (bohr^-2) idir=direction of atomic displacement (=1,2 or 3 : displacement of atom ipert along the 1st, 2nd or 3rd axis). ipert=type of the perturbation ixc= choice of exchange-correlation scheme kxc(nfft,nkxc)=exchange-correlation kernel 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 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 array kxc, see rhotoxc.f for a description nspden=number of spin-density components n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used optene=0: the contributions to the 2nd order energy are not computed 1: the contributions to the 2nd order energy are computed optres=0: the trial potential residual is computed ; the input potential value is kept 1: the new value of the trial potential is computed in place of the input value paral_kgb=flag controlling (k,g,bands) parallelization qphon(3)=reduced coordinates for the phonon wavelength rhog(2,nfft)=array for Fourier transform of GS electron density rhog1(2,nfft)=RF electron density in reciprocal space rhor(nfft,nspden)=array for GS electron density in electrons/bohr**3. rhor1(cplex*nfft,nspden)=RF electron density in real space (electrons/bohr**3). rprimd(3,3)=dimensional primitive translations in real space (bohr) ucvol=unit cell volume in ($\textrm{bohr}^{3}$) usepaw= 0 for non paw calculation; =1 for paw calculation usexcnhat= -PAW only- flag controling use of compensation density in Vxc vpsp1(cplex*nfft)=first-order derivative of the ionic potential xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc
OUTPUT
vhartr1(cplex*nfft)=1-order Hartree potential (not output if size=0)
vxc1(cplex*nfft,nspden)= 1st-order XC potential (not output if size=0)
==== if optene==1
ehart01=inhomogeneous 1st-order Hartree part of 2nd-order total energy
ehart1=1st-order Hartree part of 2nd-order total energy
exc1=1st-order exchange-correlation part of 2nd-order total energy
elpsp1=1st-order local pseudopot. part of 2nd-order total energy.
==== if optres==0
vresid1(cplex*nfft,nspden)=potential residual
vres2=square of the norm of the residual
SIDE EFFECTS
==== if optres==1
vtrial1(cplex*nfft,nspden)= new value of 1st-order trial potential
PARENTS
dfpt_scfcv
CHILDREN
dfpt_mkvxc,dfpt_mkvxc_noncoll,dfpt_mkvxcstr,dfpt_v1zeeman,dotprod_vn hartre,hartrestr,sqnorm_v,timab
SOURCE
79 #if defined HAVE_CONFIG_H 80 #include "config.h" 81 #endif 82 83 #include "abi_common.h" 84 85 86 subroutine dfpt_rhotov(cplex,ehart01,ehart1,elpsp1,exc1,elmag1,gsqcut,idir,ipert,& 87 & ixc,kxc,mpi_enreg,natom,nfft,ngfft,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,nspden,n3xccc,& 88 & optene,optres,paral_kgb,qphon,rhog,rhog1,rhor,rhor1,rprimd,ucvol,& 89 & usepaw,usexcnhat,vhartr1,vpsp1,vresid1,vres2,vtrial1,vxc,vxc1,xccc3d1,ixcrot) 90 91 use defs_basis 92 use defs_abitypes 93 use m_profiling_abi 94 use m_errors 95 use m_cgtools 96 97 !This section has been created automatically by the script Abilint (TD). 98 !Do not modify the following lines by hand. 99 #undef ABI_FUNC 100 #define ABI_FUNC 'dfpt_rhotov' 101 use interfaces_18_timing 102 use interfaces_56_xc 103 use interfaces_72_response, except_this_one => dfpt_rhotov 104 !End of the abilint section 105 106 implicit none 107 108 !Arguments ------------------------------------ 109 !scalars 110 integer,intent(in) :: cplex,idir,ipert,ixc,n3xccc,natom,nfft,nhat1grdim,nkxc,nspden 111 integer,intent(in) :: optene,optres,paral_kgb,usepaw,usexcnhat,ixcrot 112 real(dp),intent(in) :: gsqcut,ucvol 113 real(dp),intent(inout) :: ehart01 !vz_i 114 real(dp),intent(out) :: vres2 115 type(MPI_type),intent(in) :: mpi_enreg 116 !arrays 117 real(dp),intent(in) :: kxc(nfft,nkxc) 118 real(dp),intent(in) :: vxc(nfft,nspden) 119 real(dp),intent(in) :: nhat(nfft,nspden) 120 real(dp),intent(in) :: nhat1(cplex*nfft,nspden) !vz_d 121 real(dp),intent(in) :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim) 122 real(dp),intent(in) :: qphon(3),rhog(2,nfft) 123 real(dp),intent(in) :: rhog1(2,nfft) 124 real(dp),target,intent(in) :: rhor(nfft,nspden),rhor1(cplex*nfft,nspden) 125 real(dp),intent(in) :: rprimd(3,3),vpsp1(cplex*nfft) 126 real(dp),intent(in) :: xccc3d1(cplex*n3xccc) 127 real(dp),intent(inout) :: vtrial1(cplex*nfft,nspden),elpsp1,ehart1,exc1,elmag1 128 real(dp),intent(out) :: vresid1(cplex*nfft,nspden) 129 real(dp),target,intent(out) :: vhartr1(:),vxc1(:,:) 130 131 !Local variables------------------------------- 132 !scalars 133 integer :: ifft,ispden,nfftot,option 134 integer :: optnc,nkxc_cur 135 logical :: vhartr1_allocated,vxc1_allocated 136 real(dp) :: doti,elpsp10 137 !arrays 138 integer,intent(in) :: ngfft(18) 139 real(dp) :: tsec(20) 140 real(dp),allocatable :: rhor1_nohat(:,:),vhartr01(:),vxc1val(:,:) 141 real(dp),pointer :: rhor1_(:,:),vhartr1_(:),vxc1_(:,:),v1zeeman(:,:) 142 143 ! ********************************************************************* 144 145 call timab(157,1,tsec) 146 147 !FR EB SPr 148 if (nspden==4) then 149 if(usepaw==1) then 150 MSG_ERROR('DFPT with nspden=4 works only for norm-conserving psp!') 151 end if 152 end if 153 154 !Get size of FFT grid 155 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 156 157 !Eventually allocate temporary memory space 158 vhartr1_allocated=(size(vhartr1)>0) 159 if (vhartr1_allocated) then 160 vhartr1_ => vhartr1 161 else 162 ABI_ALLOCATE(vhartr1_,(cplex*nfft)) 163 end if 164 vxc1_allocated=(size(vxc1)>0) 165 if (vxc1_allocated) then 166 vxc1_ => vxc1 167 else 168 ABI_ALLOCATE(vxc1_,(cplex*nfft,nspden)) 169 end if 170 171 !If needed, store pseudo density without charge compensation 172 if (usepaw==1.and.usexcnhat==0) then 173 ABI_ALLOCATE(rhor1_,(cplex*nfft,nspden)) 174 rhor1_(:,:)=rhor1(:,:)-nhat1(:,:) 175 else 176 rhor1_ => rhor1 177 end if 178 179 180 if(ipert==natom+5)then 181 ABI_ALLOCATE(v1zeeman,(cplex*nfft,nspden)) 182 call dfpt_v1zeeman(nspden,nfft,cplex,idir,v1zeeman) 183 end if 184 185 !------ Compute 1st-order Hartree potential (and energy) ---------------------- 186 187 call hartre(cplex,gsqcut,0,mpi_enreg,nfft,ngfft,paral_kgb,rhog1,rprimd,vhartr1_,qpt=qphon) 188 189 if (optene>0) then 190 call dotprod_vn(cplex,rhor1,ehart1,doti,nfft,nfftot,1,1,vhartr1_,ucvol) 191 end if 192 193 if (optene>0) ehart01=zero 194 if(ipert==natom+3 .or. ipert==natom+4) then 195 ABI_ALLOCATE(vhartr01,(cplex*nfft)) 196 call hartrestr(gsqcut,idir,ipert,mpi_enreg,natom,nfft,ngfft,paral_kgb,rhog,rprimd,vhartr01) 197 if (optene>0) then 198 call dotprod_vn(cplex,rhor1,ehart01,doti,nfft,nfftot,1,1,vhartr01,ucvol) 199 ehart01=two*ehart01 200 ehart1=ehart1+ehart01 201 end if 202 ! Note that there is a factor 2.0_dp difference with the similar GS formula 203 vhartr1_(:)=vhartr1_(:)+vhartr01(:) 204 205 ABI_DEALLOCATE(vhartr01) 206 end if 207 208 !------ Compute 1st-order XC potential (and energy) ---------------------- 209 !(including the XC core correction) 210 211 !Compute Vxc^(1) (with or without valence contribution according to options) 212 option=0;if (optene==0) option=1 213 if(ipert==natom+3.or.ipert==natom+4) then 214 call dfpt_mkvxcstr(cplex,idir,ipert,kxc,mpi_enreg,natom,nfft,ngfft,nhat,& 215 & nhat1,nkxc,nspden,n3xccc,option,paral_kgb,qphon,rhor,rhor1,rprimd,& 216 & usepaw,usexcnhat,vxc1_,xccc3d1) 217 else 218 ! FR EB non-collinear magnetism 219 ! the second nkxc should be nkxc_cur (see 67_common/nres2vres.F90) 220 if (nspden==4) then 221 optnc=1 222 nkxc_cur=nkxc ! TODO: remove nkxc_cur? 223 224 call dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,& 225 & nspden,n3xccc,optnc,option,paral_kgb,qphon,rhor,rhor1,rprimd,usexcnhat,vxc,vxc1_,xccc3d1,ixcrot=ixcrot) 226 227 else 228 call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,& 229 & nspden,n3xccc,option,paral_kgb,qphon,rhor1,rprimd,usexcnhat,vxc1_,xccc3d1) 230 end if !nspden==4 231 end if 232 233 !Compute local contribution to 2nd-order energy (includes Vxc and Vpsp and Vmag) 234 if (optene>0) then 235 if (usepaw==0) then 236 call dotprod_vn(cplex,rhor1,elpsp10,doti,nfft,nfftot,nspden,1,vxc1_,ucvol) 237 call dotprod_vn(cplex,rhor1,elpsp1 ,doti,nfft,nfftot,1 ,1,vpsp1,ucvol) 238 if (ipert==natom+5) then 239 call dotprod_vn(cplex,rhor1,elmag1 ,doti,nfft,nfftot,nspden,1,v1zeeman,ucvol) 240 end if 241 else 242 if (usexcnhat/=0) then 243 ABI_ALLOCATE(rhor1_nohat,(cplex*nfft,1)) 244 rhor1_nohat(:,1)=rhor1(:,1)-nhat1(:,1) 245 call dotprod_vn(cplex,rhor1 ,elpsp10,doti,nfft,nfftot,nspden,1,vxc1_,ucvol) 246 call dotprod_vn(cplex,rhor1_nohat,elpsp1 ,doti,nfft,nfftot,1 ,1,vpsp1,ucvol) 247 ABI_DEALLOCATE(rhor1_nohat) 248 else 249 call dotprod_vn(cplex,rhor1_,elpsp10,doti,nfft,nfftot,nspden,1,vxc1_,ucvol) 250 call dotprod_vn(cplex,rhor1_,elpsp1 ,doti,nfft,nfftot,1 ,1,vpsp1,ucvol) 251 end if 252 end if 253 254 ! Note that there is a factor 2 difference with the similar GS formula 255 elpsp1=two*(elpsp1+elpsp10) 256 end if 257 258 259 !Compute XC valence contribution exc1 and complete eventually Vxc^(1) 260 if (optene>0) then 261 ABI_ALLOCATE(vxc1val,(cplex*nfft,nspden)) 262 vxc1val=zero 263 option=2 264 !FR SPr EB non-collinear magnetism 265 if (nspden==4) then 266 optnc=1 267 nkxc_cur=nkxc 268 call dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,& 269 & nspden,n3xccc,optnc,option,paral_kgb,qphon,rhor,rhor1,rprimd,usexcnhat,vxc,vxc1val,xccc3d1,ixcrot=ixcrot) 270 else 271 call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,& 272 & nspden,n3xccc,option,paral_kgb,qphon,rhor1,rprimd,usexcnhat,vxc1val,xccc3d1) 273 end if !nspden==4 274 vxc1_(:,:)=vxc1_(:,:)+vxc1val(:,:) 275 call dotprod_vn(cplex,rhor1_,exc1,doti,nfft,nfftot,nspden,1,vxc1val,ucvol) 276 ABI_DEALLOCATE(vxc1val) 277 end if 278 279 if (usepaw==1.and.usexcnhat==0) then 280 ABI_DEALLOCATE(rhor1_) 281 end if 282 283 !DEBUG (do not take away) 284 !Compute NSC energy ensc1 associated with rhor1 in vtrial1, for debugging purposes 285 !call dotprod_vn(cplex,rhor1,ensc1,doti,nfft,nfftot,nspden,1,vtrial1,ucvol) 286 !write(std_out,*)' ek0+eeig0+eloc0=',ek0+eeig0+eloc0 287 !write(std_out,*)' ensc1=',ensc1 288 !Compute NSC energy associated with vtrial1, for debugging purposes 289 !call dotprod_vn(cplex,rhor1,ensc1,doti,mpi_enreg,nfft,nfftot,nspden,1,vtrial1,ucvol) 290 !ensc1=ensc1+half*enl1 291 !write(std_out,*)' dfpt_rhotov : check NSC energy, diff=',& 292 !& ek0+edocc+eeig0+eloc0+enl0+ensc1 293 !write(std_out,*)' evarNSC=',ek0+edocc+eeig0+eloc0+enl0 294 !write(std_out,*)' ensc1,exc1=',ensc1,exc1 295 !ENDDEBUG 296 297 !Here, vhartr1 contains Hartree potential, vpsp1 contains local psp, 298 !while vxc1 contain xc potential 299 300 !------ Produce residual vector and square of norm of it ------------- 301 !(only if requested ; if optres==0) 302 if (optres==0) then 303 !$OMP PARALLEL DO COLLAPSE(2) 304 do ispden=1,min(nspden,2) 305 do ifft=1,cplex*nfft 306 vresid1(ifft,ispden)=vhartr1_(ifft)+vxc1_(ifft,ispden)+vpsp1(ifft)-vtrial1(ifft,ispden) 307 end do 308 end do 309 if(nspden==4)then 310 !$OMP PARALLEL DO COLLAPSE(2) 311 do ispden=3,4 312 do ifft=1,cplex*nfft 313 vresid1(ifft,ispden)=vxc1_(ifft,ispden)-vtrial1(ifft,ispden) 314 end do 315 end do 316 end if 317 318 if (ipert==natom+5) then 319 vresid1 = vresid1 + v1zeeman 320 end if 321 ! Compute square norm vres2 of potential residual vresid 322 call sqnorm_v(cplex,nfft,vres2,nspden,optres,vresid1) 323 324 else 325 326 ! ------ Produce new value of trial potential------------- 327 ! (only if requested ; if optres==1) 328 329 !$OMP PARALLEL DO COLLAPSE(2) 330 do ispden=1,min(nspden,2) 331 do ifft=1,cplex*nfft 332 vtrial1(ifft,ispden)=vhartr1_(ifft)+vxc1_(ifft,ispden)+vpsp1(ifft) 333 end do 334 end do 335 if(nspden==4)then 336 !$OMP PARALLEL DO COLLAPSE(2) 337 do ispden=3,4 338 do ifft=1,cplex*nfft 339 vtrial1(ifft,ispden)=vxc1_(ifft,ispden) 340 end do 341 end do 342 end if 343 344 if (ipert==natom+5) then 345 vtrial1 = vtrial1 + v1zeeman 346 end if 347 348 end if 349 350 !Release temporary memory space 351 if (.not.vhartr1_allocated) then 352 ABI_DEALLOCATE(vhartr1_) 353 end if 354 if (.not.vxc1_allocated) then 355 ABI_DEALLOCATE(vxc1_) 356 end if 357 358 if (ipert==natom+5) then 359 ABI_DEALLOCATE(v1zeeman) 360 end if 361 362 call timab(157,2,tsec) 363 364 end subroutine dfpt_rhotov