TABLE OF CONTENTS
- ABINIT/m_dfptlw_pert
- ABINIT/m_dfptlw_pert/dfptlw_pert
- ABINIT/m_dfptlw_pert/lw_elecstic
- ABINIT/m_dfptlw_pert/preca_ffnl
ABINIT/m_dfptlw_pert [ Modules ]
NAME
m_dfptlw_pert
FUNCTION
COPYRIGHT
Copyright (C) 2022-2024 ABINIT group (MR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
PARENTS
CHILDREN
SOURCE
20 #if defined HAVE_CONFIG_H 21 #include "config.h" 22 #endif 23 24 #include "abi_common.h" 25 26 module m_dfptlw_pert 27 28 use defs_basis 29 use defs_abitypes 30 use defs_datatypes 31 use m_dtset 32 use m_dtfil 33 use m_errors 34 use m_profiling_abi 35 use m_hamiltonian 36 use m_cgtools 37 use m_pawcprj 38 use m_pawfgr 39 use m_wfk 40 use m_xmpi 41 use m_getgh1c 42 use m_mklocl 43 use m_initylmg, only : initylmg 44 use m_fstrings, only : itoa, sjoin 45 use m_io_tools, only : file_exists 46 use m_time, only : cwtime 47 use m_kg, only : mkkpg 48 use m_mpinfo, only : proc_distrb_cycle 49 use m_dfptlw_wf 50 use m_dfpt_mkvxc, only : dfpt_mkvxcggadq 51 use m_dfptlw_nv, only : dfptlw_geom 52 use m_spacepar, only : hartredq 53 use m_cgtools, only : dotprod_vn 54 use m_mkffnl, only : mkffnl 55 56 implicit none 57 58 public :: dfptlw_pert 59 public :: preca_ffnl 60 61 private 62 63 ! ************************************************************************* 64 65 contains
ABINIT/m_dfptlw_pert/dfptlw_pert [ Functions ]
NAME
dfptlw_pert
FUNCTION
Compute first-order response function contributions to the spatial-dispersion 3rd order energy derivatives of the longwave driver. The main inputs are : - GS WFs and Hamiltonian (cg,gs_hamkq) - 1st-order WFs for two perturbations i1pert/i1dir,i2pert/i2dir (cg1,cg2) - 1st-order Local+SCF potentials for i1pert and i2pert - 1st-order WFs DDK and 2nd-order WF D2_DKDK (d2_dkdk_f)
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (MR) 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
cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = array for planewave
coefficients of wavefunctions
cg1 = first derivative of cg with respect the perturbation i1pert
cg2 = first derivative of cg with respect the perturbation i2pert
cplex= if 1, real space 1-order functions on FFT grid are REAL,
if 2, COMPLEX
dimffnl= third dimension of ffnl
dtset <type(dataset_type)>=all input variables for this dataset
eigen1(2*mband*mband*nkpt*nsppol)=1st-order eigenvalues for i1pert,i1dir (hartree)
eigen2(2*mband*mband*nkpt*nsppol)=1st-order eigenvalues for i2pert,i2dir (hartree)
ffnl(dtset%mkmem,dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives
gmet(3,3)=reciprocal space metric tensor in bohr**-2
gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
gsqcut=large sphere cut-off
i1dir,i2dir,i3dir=directions of the corresponding perturbations
i1pert,i2pert,i3pert = type of perturbation that has to be computed
kg(3,mpw*mkmem_rbz)=reduced planewave coordinates
kxc(nfft,nkxc)=exchange and correlation kernel
mband = maximum number of bands
mkmem_rbz = maximum number of k points which can fit in core memory
mk1mem = maximum number of k points for first-order WF
which can fit in core memory
mpert =maximum number of ipert
mpi_enreg=MPI-parallelisation information
mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
mpw = maximum number of planewaves in basis sphere (large number)
natom = number of atoms in unit cell
n1dq= third dimension of vlocal1_i1pertdq
n2dq= third dimension of vlocal1_i2pertdq
nfft= number of FFT grid points (for this proc)
ngfft(1:18)=integer array with FFT box dimensions and other
nkpt = number of k points
nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed.
nspden = number of spin-density components
nspinor = number of spinorial components of the wavefunctions
nsppol = number of channels for spin-polarization (1 or 2)
npwarr(nkpt) = array holding npw for each k point
nylmgr=second dimension of ylmgr_k
occ(mband*nkpt*nsppol) = occupation number for each band and k
pawfgr <type(pawfgr_type)>=fine grid parameters and related data
psps <type(pseudopotential_type)> = variables related to pseudopotentials
rho1g1(2,nfft)=G-space RF electron density in electrons/bohr**3 (i1pert)
rho1r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i1pert)
rho2r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i2pert)
rmet(3,3)=real space metric tensor in bohr**2
rprimd(3,3) = dimensional primitive translations (bohr)
samepert= .true. if i1pert=i2pert and i1dir=i2dir
ucvol=volume of the unit cell
useylmgr= if 1 use the derivative of spherical harmonics
vpsp1_i1pertdq(cplex*nfft,nspden,n1dq)= local potential of first-order
gradient Hamiltonian for i1pert
vpsp1_i1pertdq(cplex*nfft,nspden,n1dq)= local potential of second-order
gradient Hamiltonian for i1pert
vpsp1_i1pertdq_geom(cplex*nfft,nspden,3)= local potential of first-order
gradient Hamiltonian for i1pert wrp to i3dir and i2dir
vpsp1_i2pertdq(cplex*nfft,nspden,n2dq)= local potential of first-order
gradient Hamiltonian for i2pert
ddk_f = wf files
d2_dkdk_f = wf files
ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics
ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics
OUTPUT
d3etot(2,3,mpert,3,mpert,3,mpert) = third derivatives of the energy tensor d3etot_t4(2,n2dq)= t4 term which might need to be converted to type-II d3etot_t5(2,n1dq)= t5 term which might need to be converted to type-II d3etot_tgeom(2,2)= Geometric term which needs to be converted to type-II
SIDE EFFECTS
TO DO!
PARENTS
m_dfptlw_loop
CHILDREN
dotprod_vn
SOURCE
168 subroutine dfptlw_pert(cg,cg1,cg2,cplex,d3etot,d3etot_t4,d3etot_t5,d3etot_tgeom,& 169 & dimffnl,dtset,eigen1,eigen2,ffnl,gmet,gs_hamkq,gsqcut,i1dir,i2dir,i3dir,& 170 & i1pert,i2pert,i3pert,kg,kxc,mband,mkmem_rbz,mk1mem,mpert,mpi_enreg,mpsang,mpw,natom,& 171 & n1dq,n2dq,nfft,ngfft,nkpt,nkxc,& 172 & nspden,nspinor,nsppol,npwarr,nylmgr,occ,pawfgr,psps,rho1g1,rho1r1,rho2r1,rmet,rprimd,samepert,& 173 & ucvol,useylmgr,vpsp1_i1pertdq,vpsp1_i1pertdqdq,vpsp1_i1pertdq_geom,vpsp1_i2pertdq,ddk_f,d2_dkdk_f,d2_dkdk_f2,ylm,ylmgr) 174 175 !Arguments ------------------------------------ 176 !scalars 177 integer,intent(in) :: cplex,dimffnl,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,mband 178 integer,intent(in) :: mk1mem,mkmem_rbz,mpert,mpsang,mpw,natom,n1dq,n2dq,nfft,nkpt,nkxc,nspden 179 integer,intent(in) :: nspinor,nsppol,nylmgr,useylmgr 180 real(dp),intent(in) :: gsqcut,ucvol 181 logical,intent(in) :: samepert 182 type(MPI_type),intent(inout) :: mpi_enreg 183 type(dataset_type),intent(in) :: dtset 184 type(pseudopotential_type),intent(in) :: psps 185 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 186 type(pawfgr_type),intent(in) :: pawfgr 187 type(wfk_t),intent(inout) :: ddk_f,d2_dkdk_f, d2_dkdk_f2 188 189 !arrays 190 integer,intent(in) :: kg(3,mpw*mkmem_rbz),ngfft(18),npwarr(nkpt) 191 real(dp),intent(in) :: eigen1(2*mband*mband*nkpt*nsppol) 192 real(dp),intent(in) :: eigen2(2*mband*mband*nkpt*nsppol) 193 real(dp),intent(in) :: ffnl(mkmem_rbz,mpw,dimffnl,psps%lmnmax,psps%ntypat) 194 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) 195 real(dp),intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol) 196 real(dp),intent(in) :: cg2(2,mpw*nspinor*mband*mk1mem*nsppol) 197 real(dp),intent(in) :: gmet(3,3),kxc(nfft,nkxc) 198 real(dp),intent(in) :: occ(mband*nkpt*nsppol) 199 real(dp),intent(in) :: rho1g1(2,nfft),rho1r1(cplex*nfft,dtset%nspden) 200 real(dp),intent(in) :: rho2r1(cplex*nfft,dtset%nspden) 201 real(dp),intent(in) :: rmet(3,3),rprimd(3,3) 202 real(dp),intent(in) :: vpsp1_i1pertdq(2*nfft,nspden,n1dq) 203 real(dp),intent(in) :: vpsp1_i1pertdqdq(2*nfft,nspden,n2dq) 204 real(dp),intent(in) :: vpsp1_i1pertdq_geom(2*nfft,nspden,3) 205 real(dp),intent(in) :: vpsp1_i2pertdq(2*nfft,nspden,n2dq) 206 real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert) 207 real(dp),intent(out) :: d3etot_t4(2,n2dq),d3etot_t5(2,n1dq) 208 real(dp),intent(out) :: d3etot_tgeom(2,n2dq) 209 real(dp),intent(in) :: ylm(mpw*mk1mem,psps%mpsang*psps%mpsang*psps%useylm) 210 real(dp),intent(in) :: ylmgr(mpw*mk1mem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 211 212 !Variables ------------------------------------ 213 !scalars 214 integer :: bandtot,bd2tot,icg,idq,ierr,ii,ikc,ikg,ikpt,ilm,isppol,istwf_k,me,n1,n2,n3,n4,n5,n6 215 integer :: nband_k,npw_k,spaceworld,tim_getgh1c 216 integer :: usepaw 217 real(dp) :: tmpim,tmpre,wtk_k 218 real(dp) :: cpu, wall, gflops 219 character(len=1000) :: msg 220 logical :: with_nonlocal_i1pert, with_nonlocal_i2pert 221 !arrays 222 integer,allocatable :: kg_k(:,:) 223 real(dp) :: d3etot_t1(2),d3etot_t1_k(2) 224 real(dp) :: d3etot_t2(2),d3etot_t2_k(2) 225 real(dp) :: d3etot_t3(2),d3etot_t3_k(2) 226 real(dp) :: d3etot_t4_k(2,n2dq) 227 real(dp) :: d3etot_t5_k(2,n1dq) 228 real(dp) :: d3etot_tgeom_k(2,n2dq) 229 real(dp) :: d3etot_telec(2) 230 real(dp) :: e3tot(2),kpt(3) 231 real(dp),allocatable :: eig1_k(:),eig2_k(:),occ_k(:) 232 real(dp),allocatable :: ylm_k(:,:),ylmgr_k(:,:,:) 233 real(dp),allocatable :: ffnl_k(:,:,:,:) 234 235 ! ************************************************************************* 236 237 DBG_ENTER("COLL") 238 239 write(msg,'(2a,3(a,i2,a,i1))') ch10,'LONGWAVE : ',& 240 ' perts : ',i1pert,'.',i1dir,' / ',i2pert,'.',i2dir,' / ',i3pert,'.',i3dir 241 call wrtout(std_out,msg,'COLL') 242 call wrtout(ab_out,msg,'COLL') 243 244 !Init parallelism 245 spaceworld=mpi_enreg%comm_cell 246 me=mpi_enreg%me_kpt 247 248 !Additional definitions 249 tim_getgh1c=0 250 usepaw=dtset%usepaw 251 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 252 n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6) 253 with_nonlocal_i1pert=.true. ; if (i1pert==natom+2) with_nonlocal_i1pert=.false. 254 with_nonlocal_i2pert=.true. ; if (i2pert==natom+2) with_nonlocal_i2pert=.false. 255 256 !Initialize d3etot parts 257 d3etot_t1=zero 258 d3etot_t2=zero 259 d3etot_t3=zero 260 d3etot_t4=zero 261 d3etot_t5=zero 262 d3etot_telec=zero 263 d3etot_tgeom=zero 264 265 !Calculate the electrostatic contribution 266 call lw_elecstic(cplex,d3etot_telec,gmet,gs_hamkq%gprimd,gsqcut,& 267 & i3dir,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,rho1g1,rho1r1,rho2r1,ucvol) 268 269 !Loop over spins 270 bandtot = 0 271 bd2tot = 0 272 icg=0 273 do isppol = 1, nsppol 274 275 ! Loop over k-points 276 ikg = 0 277 ikc = 0 278 do ikpt = 1, nkpt 279 280 nband_k = dtset%nband(ikpt+(isppol-1)*nkpt) 281 npw_k = npwarr(ikpt) 282 istwf_k = dtset%istwfk(ikpt) 283 284 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,mpi_enreg%me)) then 285 bandtot = bandtot + nband_k 286 bd2tot = bd2tot + 2*nband_k**2 287 cycle ! Skip the rest of the k-point loop 288 end if 289 ikc= ikc + 1 290 291 ABI_MALLOC(occ_k,(nband_k)) 292 occ_k(:) = occ(1+bandtot:nband_k+bandtot) 293 wtk_k = dtset%wtk(ikpt) 294 kpt(:) = dtset%kptns(:,ikpt) 295 296 ABI_MALLOC(eig1_k,(2*nband_k**2)) 297 ABI_MALLOC(eig2_k,(2*nband_k**2)) 298 ABI_MALLOC(kg_k,(3,npw_k)) 299 ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm)) 300 ABI_MALLOC(ylmgr_k,(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 301 ABI_MALLOC(ffnl_k,(npw_k,dimffnl,psps%lmnmax,psps%ntypat)) 302 303 !Get plane-wave vectors and related data at k 304 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 305 if (dtset%ffnl_lw==1) then 306 if (psps%useylm==1) then 307 do ilm=1,psps%mpsang*psps%mpsang 308 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 309 end do 310 if (useylmgr==1) then 311 do ilm=1,psps%mpsang*psps%mpsang 312 do ii=1,nylmgr 313 ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm) 314 end do 315 end do 316 end if 317 end if 318 else if (dtset%ffnl_lw==0) then 319 ffnl_k(1:npw_k,:,:,:)=ffnl(ikc,1:npw_k,:,:,:) 320 end if 321 322 !Get matrix elements for uniform perturbations 323 eig1_k(:)=eigen1(1+bd2tot:2*nband_k**2+bd2tot) 324 eig2_k(:)=eigen2(1+bd2tot:2*nband_k**2+bd2tot) 325 326 !Compute the stationary terms of d3etot depending on response functions 327 call dfpt_1wf(cg,cg1,cg2,cplex,ddk_f,d2_dkdk_f,d2_dkdk_f2,d3etot_t1_k,d3etot_t2_k,d3etot_t3_k,& 328 & d3etot_t4_k,d3etot_t5_k,dimffnl,dtset,eig1_k,eig2_k,ffnl_k,gs_hamkq,icg,& 329 & i1dir,i2dir,i3dir,i1pert,i2pert,ikpt,isppol,istwf_k,& 330 & kg_k,kpt,mkmem_rbz,mpi_enreg,mpw,natom,nband_k,& 331 & n1dq,n2dq,nfft,ngfft,npw_k,nspden,nsppol,nylmgr,occ_k,& 332 & pawfgr,psps,rmet,rprimd,samepert,useylmgr,& 333 & vpsp1_i1pertdq,vpsp1_i2pertdq,& 334 & wtk_k,ylm_k,ylmgr_k) 335 336 ! Add the contribution from each k-point. 337 d3etot_t1=d3etot_t1 + d3etot_t1_k 338 d3etot_t2=d3etot_t2 + d3etot_t2_k 339 d3etot_t3=d3etot_t3 + d3etot_t3_k 340 d3etot_t4=d3etot_t4 + d3etot_t4_k 341 d3etot_t5=d3etot_t5 + d3etot_t5_k 342 343 !Compute the nonvariational geometric term 344 call cwtime(cpu, wall, gflops, "start") 345 if (i1pert<=natom.and.(i2pert==natom+3.or.i2pert==natom+4)) then 346 call dfptlw_geom(cg,d3etot_tgeom_k,dimffnl,dtset, & 347 & ffnl_k,gs_hamkq,icg, & 348 & i1dir,i2dir,i3dir,i1pert,i2pert,ikpt, & 349 & isppol,istwf_k,kg_k,kpt,mkmem_rbz,mpi_enreg,natom,mpw,nband_k,n2dq,nfft, & 350 & ngfft,npw_k,nspden,nsppol,nylmgr,occ_k, & 351 & psps,rmet,rprimd,useylmgr,vpsp1_i1pertdqdq,vpsp1_i1pertdq_geom,wtk_k,ylm_k,ylmgr_k) 352 353 !Add the contribution from each k-point 354 d3etot_tgeom=d3etot_tgeom + d3etot_tgeom_k 355 end if 356 call cwtime(cpu, wall, gflops, "stop") 357 358 ! Keep track of total number of bands 359 bandtot = bandtot + nband_k 360 bd2tot = bd2tot + 2*nband_k**2 361 362 ! Shift arrays memory 363 icg=icg+npw_k*dtset%nspinor*nband_k 364 ikg=ikg+npw_k 365 366 ABI_FREE(eig1_k) 367 ABI_FREE(eig2_k) 368 ABI_FREE(occ_k) 369 ABI_FREE(kg_k) 370 ABI_FREE(ylm_k) 371 ABI_FREE(ylmgr_k) 372 ABI_FREE(ffnl_k) 373 374 end do !ikpt 375 376 end do !isppol 377 378 379 !=== MPI communications ================== 380 if (xmpi_paral==1) then 381 call xmpi_sum(d3etot_t1,spaceworld,ierr) 382 call xmpi_sum(d3etot_t2,spaceworld,ierr) 383 call xmpi_sum(d3etot_t3,spaceworld,ierr) 384 call xmpi_sum(d3etot_t4,spaceworld,ierr) 385 call xmpi_sum(d3etot_t5,spaceworld,ierr) 386 call xmpi_sum(d3etot_tgeom,spaceworld,ierr) 387 end if 388 389 !Apply +i or -i in case of strain perturbation. 390 if (i1pert==natom+3.or.i1pert==natom+4) then 391 tmpre=d3etot_telec(1);tmpim=d3etot_telec(2) ; d3etot_telec(1)=tmpim;d3etot_telec(2)=-tmpre 392 tmpre=d3etot_t1(1);tmpim=d3etot_t1(2) ; d3etot_t1(1)=tmpim;d3etot_t1(2)=-tmpre 393 tmpre=d3etot_t2(1);tmpim=d3etot_t2(2) ; d3etot_t2(1)=tmpim;d3etot_t2(2)=-tmpre 394 tmpre=d3etot_t3(1);tmpim=d3etot_t3(2) ; d3etot_t3(1)=tmpim;d3etot_t3(2)=-tmpre 395 do idq=1,n2dq 396 tmpre=d3etot_t4(1,idq);tmpim=d3etot_t4(2,idq) ; d3etot_t4(1,idq)=tmpim;d3etot_t4(2,idq)=-tmpre 397 end do 398 do idq=1,n1dq 399 tmpre=d3etot_t5(1,idq);tmpim=d3etot_t5(2,idq) ; d3etot_t5(1,idq)=tmpim;d3etot_t5(2,idq)=-tmpre 400 end do 401 end if 402 if (i2pert==natom+3.or.i2pert==natom+4) then 403 tmpre=d3etot_telec(1);tmpim=d3etot_telec(2) ; d3etot_telec(1)=-tmpim;d3etot_telec(2)=tmpre 404 tmpre=d3etot_t1(1);tmpim=d3etot_t1(2) ; d3etot_t1(1)=-tmpim;d3etot_t1(2)=tmpre 405 tmpre=d3etot_t2(1);tmpim=d3etot_t2(2) ; d3etot_t2(1)=-tmpim;d3etot_t2(2)=tmpre 406 tmpre=d3etot_t3(1);tmpim=d3etot_t3(2) ; d3etot_t3(1)=-tmpim;d3etot_t3(2)=tmpre 407 do idq=1,n2dq 408 tmpre=d3etot_t4(1,idq);tmpim=d3etot_t4(2,idq) ; d3etot_t4(1,idq)=-tmpim;d3etot_t4(2,idq)=tmpre 409 if (i1pert<=natom) then 410 tmpre=d3etot_tgeom(1,idq);tmpim=d3etot_tgeom(2,idq) ; d3etot_tgeom(1,idq)=-tmpim;d3etot_tgeom(2,idq)=tmpre 411 end if 412 end do 413 do idq=1,n1dq 414 tmpre=d3etot_t5(1,idq);tmpim=d3etot_t5(2,idq) ; d3etot_t5(1,idq)=-tmpim;d3etot_t5(2,idq)=tmpre 415 end do 416 end if 417 418 !Join all the contributions in e3tot except t4 and t5 which may need to be 419 !converted to type-II in case of strain perturbation. 420 !Apply here the two factor to the stationary wf1 contributions 421 !(see PRB 105, 064101 (2022)) 422 d3etot_t1(:)=two*d3etot_t1(:) 423 d3etot_t2(:)=two*d3etot_t2(:) 424 d3etot_t3(:)=two*d3etot_t3(:) 425 d3etot_t4(:,:)=two*d3etot_t4(:,:) 426 d3etot_t5(:,:)=two*d3etot_t5(:,:) 427 e3tot(:)=d3etot_t1(:)+d3etot_t2(:)+d3etot_t3(:)+d3etot_telec(:) 428 429 430 !Before printing, set small contributions to zero 431 !Real parts 432 if (abs(d3etot_t1(1))<tol8) d3etot_t1(1)= zero 433 if (abs(d3etot_t2(1))<tol8) d3etot_t2(1)= zero 434 if (abs(d3etot_t3(1))<tol8) d3etot_t3(1)= zero 435 do idq=1,n2dq 436 if (abs(d3etot_t4(1,idq))<tol8) d3etot_t4(1,idq)= zero 437 if (abs(d3etot_tgeom(1,idq))<tol8) d3etot_tgeom(1,idq)= zero 438 end do 439 do idq=1,n1dq 440 if (abs(d3etot_t5(1,idq))<tol8) d3etot_t5(1,idq)= zero 441 end do 442 if (abs(d3etot_telec(1))<tol8) d3etot_telec(1)= zero 443 if (abs(e3tot(1)) <tol8) e3tot(1)= zero 444 445 !Imaginary parts 446 if (abs(d3etot_t1(2))<tol8) d3etot_t1(2)= zero 447 if (abs(d3etot_t2(2))<tol8) d3etot_t2(2)= zero 448 if (abs(d3etot_t3(2))<tol8) d3etot_t3(2)= zero 449 do idq=1,n2dq 450 if (abs(d3etot_t4(2,idq))<tol8) d3etot_t4(2,idq)= zero 451 if (abs(d3etot_tgeom(2,idq))<tol8) d3etot_tgeom(2,idq)= zero 452 end do 453 do idq=1,n1dq 454 if (abs(d3etot_t5(2,idq))<tol8) d3etot_t5(2,idq)= zero 455 end do 456 if (abs(d3etot_telec(2))<tol8) d3etot_telec(2)= zero 457 if (abs(e3tot(2)) <tol8) e3tot(2)= zero 458 459 if (dtset%prtvol>=10) then 460 write(msg,'(4(a,2(a,f18.8)),a)') & 461 ch10,' d3etot_telec = ',d3etot_telec(1), ',',d3etot_telec(2),& 462 ch10,' d3etot_t1 = ',d3etot_t1(1), ',',d3etot_t1(2),& 463 ch10,' d3etot_t2 = ',d3etot_t2(1), ',',d3etot_t2(2),& 464 ch10,' d3etot_t3 = ',d3etot_t3(1), ',',d3etot_t3(2) 465 call wrtout(std_out,msg,'COLL') 466 call wrtout(ab_out,msg,'COLL') 467 if (n2dq==1) then 468 write(msg,'(2(a,f18.8))') & 469 ' d3etot_t4 = ',d3etot_t4(1,1), ',',d3etot_t4(2,1) 470 else if (n2dq==2) then 471 write(msg,'(2(2(a,f18.8),a))') & 472 ' d3etot_t4(dw shear) = ',d3etot_t4(1,1), ',',d3etot_t4(2,1),ch10,& 473 ' d3etot_t4(up shear) = ',d3etot_t4(1,2), ',',d3etot_t4(2,2) 474 end if 475 call wrtout(std_out,msg,'COLL') 476 call wrtout(ab_out,msg,'COLL') 477 if (n1dq==1) then 478 write(msg,'(2(a,f18.8))') & 479 ' d3etot_t5 = ',d3etot_t5(1,1), ',',d3etot_t5(2,1) 480 else if (n1dq==2) then 481 write(msg,'(2(2(a,f18.8),a))') & 482 ' d3etot_t5(dw shear) = ',d3etot_t5(1,1), ',',d3etot_t5(2,1),ch10,& 483 ' d3etot_t5(up shear) = ',d3etot_t5(1,2), ',',d3etot_t5(2,2) 484 end if 485 call wrtout(std_out,msg,'COLL') 486 call wrtout(ab_out,msg,'COLL') 487 if (i1pert<=natom.and.(i2pert==natom+3.or.i2pert==natom+4)) then 488 if (n2dq==1) then 489 write(msg,'(2(a,f18.8))') & 490 ' d3etot_tgeom = ',d3etot_tgeom(1,1), ',',d3etot_tgeom(2,1) 491 else if (n2dq==2) then 492 write(msg,'(2(2(a,f18.8),a))') & 493 'd3etot_tgeom(dw shear) = ',d3etot_tgeom(1,1), ',',d3etot_tgeom(2,1),ch10,& 494 'd3etot_tgeom(up shear) = ',d3etot_tgeom(1,2), ',',d3etot_tgeom(2,2) 495 end if 496 call wrtout(std_out,msg,'COLL') 497 call wrtout(ab_out,msg,'COLL') 498 end if 499 end if 500 501 d3etot(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=e3tot(:) 502 503 !Deallocations 504 505 DBG_EXIT("COLL") 506 507 end subroutine dfptlw_pert
ABINIT/m_dfptlw_pert/lw_elecstic [ Functions ]
NAME
lw_elecstic
FUNCTION
This routine calculates the electrostatic term of the spatial-dispersion third-order energy derivative for a couple of perturbations and a gradient direction.
COPYRIGHT
Copyright (C) 2022-2024 ABINIT group (MR) 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 gmet(3,3)=reciprocal space metric tensor in bohr**-2 gprimd(3,3)=reciprocal space dimensional primitive translations gsqcut=large sphere cut-off i3dir= directions of the 3th perturbations kxc(nfft,nkxc)=exchange and correlation kernel mpi_enreg=information about MPI parallelization nfft= number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed. nspden = number of spin-density components rho1g1(2,nfft)=G-space RF electron density in electrons/bohr**3 (i1pert) rho1r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i1pert) rho2r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i2pert) ucvol=volume of the unit cell
OUTPUT
d3etot_telec(2)= Electrostatic term of the third-order energy derivative
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
SOURCE
556 #if defined HAVE_CONFIG_H 557 #include "config.h" 558 #endif 559 560 #include "abi_common.h" 561 562 563 subroutine lw_elecstic(cplex,d3etot_telec,gmet,gprimd,gsqcut,& 564 & i3dir,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,rho1g1,rho1r1,rho2r1,ucvol) 565 566 use defs_basis 567 use m_errors 568 use m_profiling_abi 569 570 implicit none 571 572 !Arguments ------------------------------------ 573 integer,intent(in) :: cplex,i3dir 574 integer,intent(in) :: nfft,nkxc,nspden 575 real(dp),intent(in) :: gsqcut,ucvol 576 type(MPI_type),intent(inout) :: mpi_enreg 577 !arrays 578 integer,intent(in) :: ngfft(18) 579 real(dp),intent(in) :: gmet(3,3),gprimd(3,3) 580 real(dp),intent(in) :: rho1g1(2,nfft),rho1r1(cplex*nfft,nspden) 581 real(dp),intent(in) :: rho2r1(cplex*nfft,nspden),kxc(nfft,nkxc) 582 real(dp),intent(out) :: d3etot_telec(2) 583 584 !Local variables------------------------------- 585 !scalars 586 integer :: ii,jj,nfftot,qcar 587 real(dp) :: doti,dotr 588 589 !arrays 590 real(dp),allocatable :: rhor1_cplx(:,:) 591 real(dp),allocatable :: vxc1dq(:,:),vxc1dq_car(:,:,:),vqgradhart(:) 592 593 ! ************************************************************************* 594 595 DBG_ENTER("COLL") 596 597 !If GGA xc first calculate the Cartesian q gradient of the xc kernel 598 if (nkxc == 7) then 599 ABI_MALLOC(vxc1dq,(2*nfft,nspden)) 600 ABI_MALLOC(vxc1dq_car,(2*nfft,nspden,3)) 601 do qcar=1,3 602 call dfpt_mkvxcggadq(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,qcar,rho1r1,vxc1dq) 603 vxc1dq_car(:,:,qcar)=vxc1dq(:,:) 604 end do 605 end if 606 607 !Calculate the q gradient of the Hartree potential 608 ABI_MALLOC(vqgradhart,(2*nfft)) 609 call hartredq(2,gmet,gsqcut,mpi_enreg,nfft,ngfft,i3dir,rho1g1,vqgradhart) 610 611 !If GGA convert the gradient of xc kernel to reduced coordinates and incorporate it to the Hartree part 612 if (nkxc == 7) then 613 vxc1dq=zero 614 do qcar=1,3 615 vxc1dq(:,:)=vxc1dq(:,:) + gprimd(qcar,i3dir) * vxc1dq_car(:,:,qcar) 616 end do 617 vqgradhart(:)=vqgradhart(:)+vxc1dq(:,1) 618 ABI_FREE(vxc1dq_car) 619 end if 620 621 !Calculate the electrostatic energy term with the i2pert density response 622 !I need a complex density for the dotprod_vn 623 ABI_MALLOC(rhor1_cplx,(2*nfft,nspden)) 624 rhor1_cplx=zero 625 do ii=1,nfft 626 jj=ii*2 627 rhor1_cplx(jj-1,:)=rho2r1(ii,:) 628 end do 629 630 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 631 call dotprod_vn(2,rhor1_cplx,dotr,doti,nfft,nfftot,nspden,2,vqgradhart,ucvol) 632 633 d3etot_telec(1)=dotr 634 d3etot_telec(2)=doti 635 636 !Deallocations 637 ABI_SFREE(vxc1dq) 638 ABI_FREE(vqgradhart) 639 ABI_FREE(rhor1_cplx) 640 641 DBG_EXIT("COLL") 642 643 end subroutine lw_elecstic
ABINIT/m_dfptlw_pert/preca_ffnl [ Functions ]
NAME
preca_ffnl
FUNCTION
Calculates the nonlocal form factors and derivatives for all the atoms and k points.
COPYRIGHT
Copyright (C) 2022-2024 ABINIT group (MR) 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
dimffnl= second dimension of ffnl gmet(3,3)= reciprocal-space metric tensor gprimd(3,3)= dimensional reciprocal space primitive translations (b^-1) ider= if 1 first order derivatives of ffnl are calculated if 2 first and second order derivatives of ffnl are calculated idir0= variable that controls the way in which the derivatives of ffnl are calculated and saved kg(3,mpw)=integer coordinates of G vectors in basis sphere kptns(3,nkpt)=k points in terms of reciprocal translations mband= masimum number of bands mkmem= maximum number of k points which can fit in core memory mpi_enreg=information about MPI parallelization mpw = maximum number of planewaves in basis sphere (large number) nkpt = number of k point npwarr(nkpt)=array holding npw for each k point nylmgr=second dimension of ylmgr psps <type(pseudopotential_type)> = variables related to pseudopotentials rmet(3,3)= real-space metric tensor useylmgr= if 1 use the derivative of spherical harmonics ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics
OUTPUT
ffnl(mkmem,npw_k,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
SOURCE
696 #if defined HAVE_CONFIG_H 697 #include "config.h" 698 #endif 699 700 #include "abi_common.h" 701 702 703 subroutine preca_ffnl(dimffnl,ffnl,gmet,gprimd,ider,idir0,kg,kptns,mband,mkmem,mpi_enreg,mpw,nkpt, & 704 & npwarr,nylmgr,psps,rmet,useylmgr,ylm,ylmgr) 705 706 use defs_basis 707 use m_errors 708 use m_profiling_abi 709 710 implicit none 711 712 !Arguments ------------------------------------ 713 !scalars 714 integer , intent(in) :: dimffnl,ider,idir0,mband,mkmem,mpw,nkpt,nylmgr,useylmgr 715 type(pseudopotential_type),intent(in) :: psps 716 type(MPI_type),intent(in) :: mpi_enreg 717 !arrays 718 integer,intent(in) :: kg(3,mpw*mkmem) 719 integer,intent(in) :: npwarr(nkpt) 720 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kptns(3,nkpt),rmet(3,3) 721 real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm) 722 real(dp),intent(in) :: ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 723 real(dp),intent(out) :: ffnl(mkmem,mpw,dimffnl,psps%lmnmax,psps%ntypat) 724 725 !Local variables------------------------------- 726 !scalars 727 integer :: ii,ikc,ikg,ikpt,ilm,nkpg,npw_k 728 character(len=500) :: msg 729 !arrays 730 integer,allocatable :: kg_k(:,:) 731 real(dp) :: kpt(3) 732 real(dp),allocatable :: ffnl_k(:,:,:,:),kpg_k(:,:) 733 real(dp),allocatable :: ylm_k(:,:),ylmgr_k(:,:,:),ylmgr_k_part(:,:,:) 734 735 ! ************************************************************************* 736 737 DBG_ENTER("COLL") 738 739 !Loop over k-points 740 ikg=0 741 ikc=0 742 do ikpt = 1, nkpt 743 744 npw_k = npwarr(ikpt) 745 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,1,mpi_enreg%me)) then 746 cycle ! Skip the rest of the k-point loop 747 end if 748 ikc= ikc + 1 749 750 ABI_MALLOC(kg_k,(3,npw_k)) 751 ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 752 ABI_MALLOC(ylmgr_k,(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 753 754 kpt(:)= kptns(:,ikpt) 755 756 !Get plane-wave vectors and related data at k 757 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 758 if (psps%useylm==1) then 759 do ilm=1,psps%mpsang*psps%mpsang 760 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 761 end do 762 if (useylmgr==1) then 763 do ilm=1,psps%mpsang*psps%mpsang 764 do ii=1,nylmgr 765 ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm) 766 end do 767 end do 768 end if 769 end if 770 771 if (dimffnl==2.or.dimffnl==4) then 772 ABI_MALLOC(ylmgr_k_part,(npw_k,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 773 ylmgr_k_part(:,:,:)=ylmgr_k(:,1:3,:) 774 else if (dimffnl==10) then 775 ABI_MALLOC(ylmgr_k_part,(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 776 ylmgr_k_part(:,:,:)=ylmgr_k(:,:,:) 777 else 778 msg='wrong size for ffnl via dimffnl!' 779 ABI_BUG(msg) 780 end if 781 782 783 nkpg=0 784 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 785 ABI_MALLOC(ffnl_k,(npw_k,dimffnl,psps%lmnmax,psps%ntypat)) 786 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl_k,psps%ffspl,gmet,gprimd,ider,idir0,& 787 & psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,& 788 & npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k_part) 789 790 ffnl(ikc,1:npw_k,:,:,:)=ffnl_k(1:npw_k,:,:,:) 791 792 ABI_FREE(kg_k) 793 ABI_FREE(ylm_k) 794 ABI_FREE(ylmgr_k) 795 ABI_FREE(ylmgr_k_part) 796 ABI_FREE(ffnl_k) 797 ABI_FREE(kpg_k) 798 799 !Shift arrays memory 800 ikg=ikg+npw_k 801 802 end do 803 804 DBG_EXIT("COLL") 805 806 end subroutine preca_ffnl