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
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 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 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 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
SOURCE
108 subroutine dfpt_rhotov(cplex,ehart01,ehart1,elpsp1,exc1,elmag1,gsqcut,idir,ipert,& 109 & ixc,kxc,mpi_enreg,natom,nfft,ngfft,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,nspden,n3xccc,& 110 & non_magnetic_xc,optene,optres,qphon,rhog,rhog1,rhor,rhor1,rprimd,ucvol,& 111 & usepaw,usexcnhat,vhartr1,vpsp1,vresid1,vres2,vtrial1,vxc,vxc1,xccc3d1,ixcrot) 112 113 !Arguments ------------------------------------ 114 !scalars 115 integer,intent(in) :: cplex,idir,ipert,ixc,n3xccc,natom,nfft,nhat1grdim,nkxc,nspden 116 integer,intent(in) :: optene,optres,usepaw,usexcnhat,ixcrot 117 logical,intent(in) :: non_magnetic_xc 118 real(dp),intent(in) :: gsqcut,ucvol 119 real(dp),intent(inout) :: ehart01 120 real(dp),intent(out) :: vres2 121 type(MPI_type),intent(in) :: mpi_enreg 122 !arrays 123 real(dp),intent(in) :: kxc(nfft,nkxc) 124 real(dp),intent(in) :: vxc(nfft,nspden) 125 real(dp),intent(in) :: nhat(nfft,nspden) 126 real(dp),intent(in) :: nhat1(cplex*nfft,nspden) !vz_d 127 real(dp),intent(in) :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim) 128 real(dp),intent(in) :: qphon(3),rhog(2,nfft) 129 real(dp),intent(in) :: rhog1(2,nfft) 130 real(dp),target,intent(in) :: rhor(nfft,nspden),rhor1(cplex*nfft,nspden) 131 real(dp),intent(in) :: rprimd(3,3),vpsp1(cplex*nfft) 132 real(dp),intent(in) :: xccc3d1(cplex*n3xccc) 133 real(dp),intent(inout) :: vtrial1(cplex*nfft,nspden),elpsp1,ehart1,exc1,elmag1 134 real(dp),intent(out) :: vresid1(cplex*nfft,nspden) 135 real(dp),target,intent(out) :: vhartr1(:),vxc1(:,:) 136 137 !Local variables------------------------------- 138 !scalars 139 integer :: ifft,ispden,nfftot,option 140 integer :: optnc,nkxc_cur 141 logical :: vhartr1_allocated,vxc1_allocated 142 real(dp) :: doti,elpsp10 143 !arrays 144 integer,intent(in) :: ngfft(18) 145 real(dp) :: tsec(20) 146 real(dp),parameter :: dummyvgeo(3)=zero 147 real(dp),allocatable :: rhor1_nohat(:,:),vhartr01(:),vxc1val(:,:) 148 real(dp),pointer :: rhor1_(:,:),vhartr1_(:),vxc1_(:,:),v1zeeman(:,:) 149 150 ! ********************************************************************* 151 152 call timab(157,1,tsec) 153 154 !FR EB SPr 155 if (nspden==4) then 156 if(usepaw==1) then 157 ABI_ERROR('DFPT with nspden=4 works only for norm-conserving psp!') 158 end if 159 end if 160 161 !Get size of FFT grid 162 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 163 164 !Eventually allocate temporary memory space 165 vhartr1_allocated=(size(vhartr1)>0) 166 if (vhartr1_allocated) then 167 vhartr1_ => vhartr1 168 else 169 ABI_MALLOC(vhartr1_,(cplex*nfft)) 170 end if 171 vxc1_allocated=(size(vxc1)>0) 172 if (vxc1_allocated) then 173 vxc1_ => vxc1 174 else 175 ABI_MALLOC(vxc1_,(cplex*nfft,nspden)) 176 end if 177 178 !If needed, store pseudo density without charge compensation 179 if (usepaw==1.and.usexcnhat==0) then 180 ABI_MALLOC(rhor1_,(cplex*nfft,nspden)) 181 rhor1_(:,:)=rhor1(:,:)-nhat1(:,:) 182 else 183 rhor1_ => rhor1 184 end if 185 186 187 if(ipert==natom+5)then 188 ABI_MALLOC(v1zeeman,(cplex*nfft,nspden)) 189 call dfpt_v1zeeman(nspden,nfft,cplex,idir,v1zeeman) 190 end if 191 192 !------ Compute 1st-order Hartree potential (and energy) ---------------------- 193 194 call hartre(cplex,gsqcut,3,0,mpi_enreg,nfft,ngfft,1,zero,rhog1,rprimd,dummyvgeo,vhartr1_,qpt=qphon) 195 196 if (optene>0) then 197 call dotprod_vn(cplex,rhor1,ehart1,doti,nfft,nfftot,1,1,vhartr1_,ucvol) 198 end if 199 200 if (optene>0) ehart01=zero 201 if(ipert==natom+3 .or. ipert==natom+4) then 202 ABI_MALLOC(vhartr01,(cplex*nfft)) 203 call hartrestr(gsqcut,idir,ipert,mpi_enreg,natom,nfft,ngfft,rhog,rprimd,vhartr01) 204 if (optene>0) then 205 call dotprod_vn(cplex,rhor1,ehart01,doti,nfft,nfftot,1,1,vhartr01,ucvol) 206 ehart01=two*ehart01 207 ehart1=ehart1+ehart01 208 end if 209 ! Note that there is a factor 2.0_dp difference with the similar GS formula 210 vhartr1_(:)=vhartr1_(:)+vhartr01(:) 211 212 ABI_FREE(vhartr01) 213 end if 214 215 !------ Compute 1st-order XC potential (and energy) ---------------------- 216 !(including the XC core correction) 217 218 !Compute Vxc^(1) (with or without valence contribution according to options) 219 option=0;if (optene==0) option=1 220 if(ipert==natom+3.or.ipert==natom+4) then 221 call dfpt_mkvxcstr(cplex,idir,ipert,kxc,mpi_enreg,natom,nfft,ngfft,nhat,& 222 & nhat1,nkxc,non_magnetic_xc,nspden,n3xccc,option,qphon,rhor,rhor1,rprimd,& 223 & usepaw,usexcnhat,vxc1_,xccc3d1) 224 else 225 ! FR EB non-collinear magnetism 226 ! the second nkxc should be nkxc_cur (see 67_common/nres2vres.F90) 227 if (nspden==4) then 228 optnc=1 229 nkxc_cur=nkxc ! TODO: remove nkxc_cur? 230 231 call dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,& 232 & non_magnetic_xc,nspden,n3xccc,optnc,option,qphon,rhor,rhor1,rprimd,usexcnhat,vxc,vxc1_,xccc3d1,ixcrot=ixcrot) 233 234 else 235 call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,& 236 & non_magnetic_xc,nspden,n3xccc,option,qphon,rhor1,rprimd,usexcnhat,vxc1_,xccc3d1) 237 end if !nspden==4 238 end if 239 240 !Compute local contribution to 2nd-order energy (includes Vxc and Vpsp and Vmag) 241 if (optene>0) then 242 if (usepaw==0) then 243 call dotprod_vn(cplex,rhor1,elpsp10,doti,nfft,nfftot,nspden,1,vxc1_,ucvol) 244 call dotprod_vn(cplex,rhor1,elpsp1 ,doti,nfft,nfftot,1 ,1,vpsp1,ucvol) 245 if (ipert==natom+5) then 246 call dotprod_vn(cplex,rhor1,elmag1 ,doti,nfft,nfftot,nspden,1,v1zeeman,ucvol) 247 end if 248 else 249 if (usexcnhat/=0) then 250 ABI_MALLOC(rhor1_nohat,(cplex*nfft,1)) 251 rhor1_nohat(:,1)=rhor1(:,1)-nhat1(:,1) 252 call dotprod_vn(cplex,rhor1 ,elpsp10,doti,nfft,nfftot,nspden,1,vxc1_,ucvol) 253 call dotprod_vn(cplex,rhor1_nohat,elpsp1 ,doti,nfft,nfftot,1 ,1,vpsp1,ucvol) 254 ABI_FREE(rhor1_nohat) 255 else 256 call dotprod_vn(cplex,rhor1_,elpsp10,doti,nfft,nfftot,nspden,1,vxc1_,ucvol) 257 call dotprod_vn(cplex,rhor1_,elpsp1 ,doti,nfft,nfftot,1 ,1,vpsp1,ucvol) 258 end if 259 end if 260 261 ! Note that there is a factor 2 difference with the similar GS formula 262 elpsp1=two*(elpsp1+elpsp10) 263 end if 264 265 266 !Compute XC valence contribution exc1 and complete eventually Vxc^(1) 267 if (optene>0) then 268 ABI_MALLOC(vxc1val,(cplex*nfft,nspden)) 269 vxc1val=zero 270 option=2 271 !FR SPr EB non-collinear magnetism 272 if (nspden==4) then 273 optnc=1 274 nkxc_cur=nkxc 275 call dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat,usepaw,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,& 276 & non_magnetic_xc,nspden,n3xccc,optnc,option,qphon,rhor,rhor1,rprimd,usexcnhat,vxc,vxc1val,xccc3d1,ixcrot=ixcrot) 277 else 278 call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,usepaw,nhat1gr,nhat1grdim,nkxc,& 279 & non_magnetic_xc,nspden,n3xccc,option,qphon,rhor1,rprimd,usexcnhat,vxc1val,xccc3d1) 280 end if !nspden==4 281 vxc1_(:,:)=vxc1_(:,:)+vxc1val(:,:) 282 call dotprod_vn(cplex,rhor1_,exc1,doti,nfft,nfftot,nspden,1,vxc1val,ucvol) 283 ABI_FREE(vxc1val) 284 end if 285 286 if (usepaw==1.and.usexcnhat==0) then 287 ABI_FREE(rhor1_) 288 end if 289 290 !DEBUG (do not take away) 291 !Compute NSC energy ensc1 associated with rhor1 in vtrial1, for debugging purposes 292 !call dotprod_vn(cplex,rhor1,ensc1,doti,nfft,nfftot,nspden,1,vtrial1,ucvol) 293 !write(std_out,*)' ek0+eeig0+eloc0=',ek0+eeig0+eloc0 294 !write(std_out,*)' ensc1=',ensc1 295 !Compute NSC energy associated with vtrial1, for debugging purposes 296 !call dotprod_vn(cplex,rhor1,ensc1,doti,mpi_enreg,nfft,nfftot,nspden,1,vtrial1,ucvol) 297 !ensc1=ensc1+half*enl1 298 !write(std_out,*)' dfpt_rhotov : check NSC energy, diff=',& 299 !& ek0+edocc+eeig0+eloc0+enl0+ensc1 300 !write(std_out,*)' evarNSC=',ek0+edocc+eeig0+eloc0+enl0 301 !write(std_out,*)' ensc1,exc1=',ensc1,exc1 302 !ENDDEBUG 303 304 !Here, vhartr1 contains Hartree potential, vpsp1 contains local psp, 305 !while vxc1 contain xc potential 306 307 !------ Produce residual vector and square of norm of it ------------- 308 !(only if requested ; if optres==0) 309 if (optres==0) then 310 !$OMP PARALLEL DO COLLAPSE(2) 311 do ispden=1,min(nspden,2) 312 do ifft=1,cplex*nfft 313 vresid1(ifft,ispden)=vhartr1_(ifft)+vxc1_(ifft,ispden)+vpsp1(ifft)-vtrial1(ifft,ispden) 314 end do 315 end do 316 if(nspden==4)then 317 !$OMP PARALLEL DO COLLAPSE(2) 318 do ispden=3,4 319 do ifft=1,cplex*nfft 320 vresid1(ifft,ispden)=vxc1_(ifft,ispden)-vtrial1(ifft,ispden) 321 end do 322 end do 323 end if 324 325 if (ipert==natom+5) then 326 vresid1 = vresid1 + v1zeeman 327 end if 328 ! Compute square norm vres2 of potential residual vresid 329 call sqnorm_v(cplex,nfft,vres2,nspden,optres,vresid1) 330 331 else 332 333 ! ------ Produce new value of trial potential------------- 334 ! (only if requested ; if optres==1) 335 336 !$OMP PARALLEL DO COLLAPSE(2) 337 do ispden=1,min(nspden,2) 338 do ifft=1,cplex*nfft 339 vtrial1(ifft,ispden)=vhartr1_(ifft)+vxc1_(ifft,ispden)+vpsp1(ifft) 340 end do 341 end do 342 if(nspden==4)then 343 !$OMP PARALLEL DO COLLAPSE(2) 344 do ispden=3,4 345 do ifft=1,cplex*nfft 346 vtrial1(ifft,ispden)=vxc1_(ifft,ispden) 347 end do 348 end do 349 end if 350 351 if (ipert==natom+5) then 352 vtrial1 = vtrial1 + v1zeeman 353 end if 354 355 end if 356 357 !Release temporary memory space 358 if (.not.vhartr1_allocated) then 359 ABI_FREE(vhartr1_) 360 end if 361 if (.not.vxc1_allocated) then 362 ABI_FREE(vxc1_) 363 end if 364 365 if (ipert==natom+5) then 366 ABI_FREE(v1zeeman) 367 end if 368 369 call timab(157,2,tsec) 370 371 end subroutine dfpt_rhotov
ABINIT/dfpt_v1zeeman [ Functions ]
NAME
dfpt_v1zeeman
FUNCTION
Calculate 1st order Zeeman potential = -vec{\sigma}.\vec{b}, where
sigma is the vector of Pauli matrices and \vec{b} is the unit
vector indicating the perturbing field direction.
INPUTS
nspden = number of density matrix components
nfft = numbder of fft grid points
cplex = complex or real density matrix
idir = direction of the perturbing field in Cartesian frame
1: along x
2: along y
3: along z
4: identity matrix at each fft point is returned (for density-density response)
OUTPUT
v1zeeman(nfft*cplex,nspden)= 1st order Zeeman potential, or Identity matrix (electrostatic potential) for idir=4
SIDE EFFECTS
NOTES
The definition of components of the potential matrix differ depending on cplex for nspden=4: For cplex=1, the potential is defined as (V_upup,V_dndn,Re[V_updn],Im[V_updn]) For cplex=2, the definition is (V_upup,V_dndn,V_updn,i.V_updn)
SOURCE
406 subroutine dfpt_v1zeeman(nspden,nfft,cplex,idir,v1zeeman) 407 408 !Arguments ------------------------------------ 409 integer , intent(in) :: idir,nfft,cplex,nspden 410 real(dp), intent(inout) :: v1zeeman(cplex*nfft,nspden) 411 412 !Local variables------------------------------- 413 integer :: ifft 414 !character(len=500) :: msg 415 416 ! ************************************************************************* 417 418 DBG_ENTER("COLL") 419 420 ! if (option/=1 .and. option/=2 ) then 421 ! write(msg,'(3a,i0)')& 422 !& 'The argument option should be 1 or 2,',ch10,& 423 !& 'however, option=',option 424 ! ABI_BUG(msg) 425 ! end if 426 ! 427 ! if (sizein<1) then 428 ! write(msg,'(3a,i0)')& 429 !& 'The argument sizein should be a positive number,',ch10,& 430 !& 'however, sizein=',sizein 431 ! ABI_ERROR(msg) 432 ! end if 433 434 DBG_EXIT("COLL") 435 436 select case(cplex) 437 case(1) 438 if (nspden==4) then 439 if(idir==3)then ! Zeeman field along the 3rd axis (z) 440 v1zeeman(:,1)=-0.5d0 441 v1zeeman(:,2)=+0.5d0 442 v1zeeman(:,3)= 0.0d0 443 v1zeeman(:,4)= 0.0d0 444 else if(idir==2)then ! Zeeman field along the 2nd axis (y) 445 v1zeeman(:,1)= 0.0d0 446 v1zeeman(:,2)= 0.0d0 447 v1zeeman(:,3)= 0.0d0 448 v1zeeman(:,4)=+0.5d0 449 else ! Zeeman field along the 1st axis (x) 450 v1zeeman(:,1)= 0.0d0 451 v1zeeman(:,2)= 0.0d0 452 v1zeeman(:,3)=-0.5d0 453 v1zeeman(:,4)= 0.0d0 454 end if 455 else if (nspden==2) then 456 v1zeeman(:,1)=-0.5e0 457 v1zeeman(:,2)= 0.5e0 458 else 459 v1zeeman(:,1)= 0.0e0 460 end if 461 case(2) 462 if (nspden==2) then 463 do ifft=1,nfft 464 v1zeeman(2*ifft-1,1) =-0.5e0 465 v1zeeman(2*ifft ,1) = 0.0e0 466 v1zeeman(2*ifft-1,2) = 0.5e0 467 v1zeeman(2*ifft ,2) = 0.0e0 468 end do 469 else if (nspden==4) then 470 select case(idir) 471 case(1) !along x, v1=-sigma_x 472 do ifft=1,nfft 473 v1zeeman(2*ifft-1,1)= 0.0e0 !Re[V^11] 474 v1zeeman(2*ifft ,1)= 0.0e0 !Im[V^11] 475 v1zeeman(2*ifft-1,2)= 0.0e0 !Re[V^22] 476 v1zeeman(2*ifft ,2)= 0.0e0 !Im[V^22] 477 v1zeeman(2*ifft-1,3)=-0.5e0 !Re[V^12] 478 v1zeeman(2*ifft ,3)= 0.0e0 !Im[V^12] 479 v1zeeman(2*ifft-1,4)= 0.0e0 !Re[i.V^21]=Im[V^12] 480 v1zeeman(2*ifft ,4)=-0.5e0 !Im[i.V^21]=Re[V^12] 481 end do 482 case(2) !along y, v1 = -sigma_y 483 do ifft=1,nfft 484 v1zeeman(2*ifft-1,1)= 0.0e0 !Re[V^11] 485 v1zeeman(2*ifft ,1)= 0.0e0 !Im[V^11] 486 v1zeeman(2*ifft-1,2)= 0.0e0 !Re[V^22] 487 v1zeeman(2*ifft ,2)= 0.0e0 !Im[V^22] 488 v1zeeman(2*ifft-1,3)= 0.0e0 !Re[V^12] 489 v1zeeman(2*ifft ,3)=+0.5e0 !Im[V^12] 490 v1zeeman(2*ifft-1,4)=+0.5e0 !Re[i.V^21]=Im[V^12] 491 v1zeeman(2*ifft ,4)= 0.0e0 !Im[i.V^21]=Re[V^12] 492 end do 493 case(3) 494 do ifft=1,nfft 495 v1zeeman(2*ifft-1,1)=-0.5e0 !Re[V^11] 496 v1zeeman(2*ifft ,1)= 0.0e0 !Im[V^11] 497 v1zeeman(2*ifft-1,2)= 0.5e0 !Re[V^22] 498 v1zeeman(2*ifft ,2)= 0.0e0 !Im[V^22] 499 v1zeeman(2*ifft-1,3)= 0.0e0 !Re[V^12] 500 v1zeeman(2*ifft ,3)= 0.0e0 !Im[V^12] 501 v1zeeman(2*ifft-1,4)= 0.0e0 !Re[i.V^21] 502 v1zeeman(2*ifft ,4)= 0.0e0 !Im[i.V^21] 503 end do 504 end select 505 end if 506 end select !cplex 507 508 end subroutine dfpt_v1zeeman
ABINIT/m_dfpt_rhotov [ Modules ]
NAME
m_dfpt_rhotov
FUNCTION
COPYRIGHT
Copyright (C) 1999-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_dfpt_rhotov 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_cgtools 28 29 use defs_abitypes, only : MPI_type 30 use m_time, only : timab 31 use m_spacepar, only : hartrestr, hartre 32 use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxc_noncoll 33 use m_dfpt_mkvxcstr, only : dfpt_mkvxcstr 34 35 implicit none 36 37 private