TABLE OF CONTENTS
ABINIT/dfptnl_loop [ Functions ]
NAME
dfptnl_loop
FUNCTION
Loop over the perturbations j1, j2 and j3
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (LB) 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
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave coefficients of wavefunctions dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree) gmet(3,3)=reciprocal space metric tensor in bohr**-2 gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1) gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere--appropriate for charge density rho(G), Hartree potential, and pseudopotentials kg(3,mpw*mkmem)=reduced planewave coordinates kxc(nfftf,nkxc)=exchange-correlation kernel k3xc(nfftf,nk3xc)=third-order exchange-correlation kernel mband = maximum number of bands mgfft = maximum single fft dimension mkmem = Number of k points treated by this node. mk1mem = Number of k points for first-order WF treated by this node. mpert =maximum number of ipert mpi_enreg=MPI-parallelisation information mpw = maximum number of planewaves in basis sphere (large number) natom = number of atoms in unit cell nattyp(ntypat)= # atoms of each type. nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90) ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid (see NOTES in respfn.F90) nhat=compensation charge density on fine rectangular grid nkpt = number of k points nkxc=second dimension of the array kxc, see rhotoxc.f for a description nk3xc=second dimension of the array k3xc 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 occ(mband*nkpt*nsppol) = occupation number for each band and k paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh paw_ij0(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS pawang <type(pawang_type)>=paw angular mesh and related data pawang1 <type(pawang_type)>=pawang datastructure containing only the symmetries preserving the perturbation pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information ph1df(2,3*(2*mgfftf+1)*natom)=one-dimensional structure factor information (fine grid) psps <type(pseudopotential_type)> = variables related to pseudopotentials rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations that have to be computed 1 -> element has to be computed explicitely -1 -> use symmetry operations to obtain the corresponding element rhog(2,nfftf)=array for Fourier transform of GS electron density rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3. rprimd(3,3)=dimensional primitive translations (bohr) ucvol = unit cell volume (bohr^3) usecprj= 1 if cprj, cprjq, cprj1 arrays are stored in memory vtrial(nfftf,nspden)=GS Vtrial(r). vxc(nfftf,nspden)=Exchange-Correlation GS potential (Hartree) xred(3,natom) = reduced atomic coordinates nsym1=number of symmetry elements in space group consistent with perturbation indsy1(4,nsym1,natom)=indirect indexing array for atom labels symaf1(nsym1)=anti(ferromagnetic) part of symmetry operations symrc1(3,3,nsym1)=symmetry operations in reciprocal space
OUTPUT
blkflg(3,mpert,3,mpert) = flags for each element of the 3DTE
(=1 if computed)
d3etot(2,3,mpert,3,mpert,3,mpert) = third derivatives of the energy tensor
= \sum_{i=1}^9 d3etot_i
d3etot_1(2,3,mpert,3,mpert,3,mpert) = 1st term of d3etot
d3etot_2(2,3,mpert,3,mpert,3,mpert) = 2nd term of d3etot
d3etot_3(2,3,mpert,3,mpert,3,mpert) = 3rd term of d3etot
d3etot_4(2,3,mpert,3,mpert,3,mpert) = 4th term of d3etot
d3etot_5(2,3,mpert,3,mpert,3,mpert) = 5th term of d3etot
d3etot_6(2,3,mpert,3,mpert,3,mpert) = 6th term of d3etot
d3etot_7(2,3,mpert,3,mpert,3,mpert) = 7th term of d3etot
d3etot_8(2,3,mpert,3,mpert,3,mpert) = 8th term of d3etot
d3etot_9(2,3,mpert,3,mpert,3,mpert) = 9th term of d3etot
SIDE EFFECTS
hdr <type(hdr_type)>=the header of wf, den and pot files
SOURCE
130 subroutine dfptnl_loop(atindx,blkflg,cg,dtfil,dtset,d3etot,eigen0,gmet,gprimd,gsqcut, & 131 & hdr,kg,kxc,k3xc,mband,mgfft,mgfftf,mkmem,mk1mem,& 132 & mpert,mpi_enreg,mpw,natom,nattyp,ngfftf,nfftf,nhat,nkpt,nkxc,nk3xc,nspinor,nsppol,& 133 & npwarr,occ,paw_an0,paw_ij0,& 134 & pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,& 135 & ph1d,ph1df,psps,rfpert,rhog,rhor,rprimd,ucvol,usecprj,vtrial,vxc,xred,& 136 & nsym1,indsy1,symaf1,symrc1,& 137 & d3etot_1,d3etot_2,d3etot_3,d3etot_4,d3etot_5,d3etot_6,d3etot_7,d3etot_8,d3etot_9) 138 139 use defs_basis 140 use defs_wvltypes 141 use m_errors 142 use m_abicore 143 use m_hdr 144 use m_nctk 145 use m_wffile 146 use m_wfk 147 use m_dtset 148 use m_dtfil 149 150 use defs_datatypes, only : pseudopotential_type 151 use defs_abitypes, only : MPI_type 152 use m_time, only : timab 153 use m_io_tools, only : file_exists 154 use m_kg, only : getph 155 use m_inwffil, only : inwffil 156 use m_fft, only : fourdp 157 use m_ioarr, only : read_rhor 158 use m_hamiltonian, only : gs_hamiltonian_type, init_hamiltonian 159 use m_pawdij, only : pawdij, pawdijfr, symdij 160 use m_pawfgr, only : pawfgr_type 161 use m_pawfgrtab, only : pawfgrtab_type 162 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags 163 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags, paw_ij_print 164 use m_pawang, only : pawang_type 165 use m_pawrad, only : pawrad_type 166 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_nullify, & 167 & pawrhoij_io, pawrhoij_inquire_dim 168 use m_paw_nhat, only : pawmknhat,pawnhatfr 169 use m_paw_denpot, only : pawdenpot 170 use m_pawtab, only : pawtab_type 171 use m_rf2, only : rf2_getidir 172 use m_initylmg, only : initylmg 173 use m_atm2fft, only : dfpt_atm2fft 174 use m_dfpt_mkvxc, only : dfpt_mkvxc 175 use m_dfpt_rhotov, only : dfpt_rhotov 176 use m_mkcore, only : dfpt_mkcore 177 use m_mklocl, only : dfpt_vlocal 178 use m_dfptnl_pert, only : dfptnl_pert 179 180 !Arguments ------------------------------------ 181 !scalars 182 integer,intent(in) :: mband,mgfft,mgfftf,mk1mem,mkmem,mpert,mpw,natom,nfftf 183 integer,intent(in) :: nk3xc,nkpt,nkxc,nspinor,nsppol,nsym1,usecprj 184 real(dp),intent(in) :: gsqcut,ucvol 185 type(MPI_type),intent(inout) :: mpi_enreg 186 type(datafiles_type),intent(in) :: dtfil 187 type(dataset_type),intent(in) :: dtset 188 type(hdr_type),intent(inout) :: hdr 189 type(pawang_type),intent(inout) :: pawang,pawang1 190 type(pawfgr_type),intent(in) :: pawfgr 191 type(pseudopotential_type),intent(in) :: psps 192 193 !arrays 194 integer,intent(in) :: atindx(natom),kg(3,mk1mem*mpw) 195 integer,intent(in) :: nattyp(psps%ntypat),ngfftf(18),npwarr(nkpt) 196 integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert) 197 integer,intent(in) :: indsy1(4,nsym1,dtset%natom),symaf1(nsym1),symrc1(3,3,nsym1) 198 integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) !vz_i 199 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gmet(3,3) 200 real(dp),intent(in) :: eigen0(dtset%mband*dtset%nkpt*dtset%nsppol) 201 real(dp),intent(in) :: gprimd(3,3),k3xc(nfftf,nk3xc),kxc(nfftf,nkxc) 202 real(dp),intent(in) :: nhat(nfftf,dtset%nspden) 203 real(dp),intent(in) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden),rprimd(3,3) 204 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),ph1df(2,3*(2*mgfftf+1)*natom) 205 real(dp),intent(in) :: vtrial(nfftf,dtset%nspden),xred(3,natom) 206 real(dp),intent(in) :: vxc(nfftf,dtset%nspden) 207 real(dp),intent(inout) :: occ(mband*nkpt*nsppol) 208 real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert) !vz_i 209 real(dp),intent(inout) :: d3etot_1(2,3,mpert,3,mpert,3,mpert) 210 real(dp),intent(inout) :: d3etot_2(2,3,mpert,3,mpert,3,mpert) 211 real(dp),intent(inout) :: d3etot_3(2,3,mpert,3,mpert,3,mpert) 212 real(dp),intent(inout) :: d3etot_4(2,3,mpert,3,mpert,3,mpert) 213 real(dp),intent(inout) :: d3etot_5(2,3,mpert,3,mpert,3,mpert) 214 real(dp),intent(inout) :: d3etot_6(2,3,mpert,3,mpert,3,mpert) 215 real(dp),intent(inout) :: d3etot_7(2,3,mpert,3,mpert,3,mpert) 216 real(dp),intent(inout) :: d3etot_8(2,3,mpert,3,mpert,3,mpert) 217 real(dp),intent(inout) :: d3etot_9(2,3,mpert,3,mpert,3,mpert) 218 type(pawfgrtab_type),intent(inout) :: pawfgrtab(natom*psps%usepaw) 219 type(pawrhoij_type),intent(in) :: pawrhoij(natom*psps%usepaw) 220 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 221 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 222 type(paw_an_type),intent(in) :: paw_an0(natom*psps%usepaw) 223 type(paw_ij_type),intent(in) :: paw_ij0(natom*psps%usepaw) 224 225 !Local variables------------------------------- 226 !scalars 227 integer,parameter :: level=51 228 integer :: ask_accurate,comm_cell,counter,cplex,cplex_rhoij,formeig,flag1,flag3 229 integer :: has_dijfr,has_diju 230 integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,iatom,idir_dkde,ierr,ii,ireadwf 231 integer :: mcg,mpsang,n1,n2,n3,n3xccc,ndir,nfftotf,nhat1grdim,npert_phon,nspden,nspden_rhoij,nwffile 232 integer :: option,optene,optfr,optorth,pert1case,pert2case,pert3case 233 integer :: qphase_rhoij,rdwrpaw,second_idir,timrev,usexcnhat 234 logical :: non_magnetic_xc 235 real(dp) :: dummy_real,dummy_real2, dummy_real3, ecut_eff 236 character(len=500) :: message 237 character(len=fnlen) :: fiden1i,fiwf1i,fiwf2i,fiwf3i,fiwfddk,fnamewff(5) 238 type(gs_hamiltonian_type) :: gs_hamkq 239 type(wffile_type) :: wff1,wff2,wff3,wfft1,wfft2,wfft3 240 type(wfk_t) :: ddk_f(5) 241 type(wvl_data) :: wvl 242 type(hdr_type) :: hdr_den 243 !arrays 244 integer :: file_index(5) 245 real(dp) :: qphon(3),tsec(2) 246 real(dp),allocatable :: cg1(:,:),cg2(:,:),cg3(:,:),eigen1(:),eigen2(:),eigen3(:) 247 real(dp),allocatable :: nhat1_i1pert(:,:),nhat1_i2pert(:,:),nhat1_i3pert(:,:) 248 real(dp),allocatable :: nhat1gr(:,:,:),vresid_dum(:,:) 249 real(dp),allocatable :: rho1r1(:,:) 250 real(dp),allocatable :: rho2g1(:,:),rho2r1(:,:),rho3r1(:,:),vhartr1_i2pert(:) 251 real(dp),allocatable :: vpsp1(:),vxc1_i2pert(:,:),work(:) 252 real(dp),allocatable,target :: vtrial1_i2pert(:,:) 253 real(dp),pointer :: vtrial1_tmp(:,:) 254 real(dp),allocatable :: xccc3d1(:),xccc3d2(:),xccc3d3(:) 255 type(pawrhoij_type),allocatable :: pawrhoij1_i1pert(:),pawrhoij1_i2pert(:),pawrhoij1_i3pert(:) 256 type(paw_an_type),allocatable :: paw_an1_i2pert(:) 257 type(paw_ij_type),allocatable :: paw_ij1_i2pert(:) 258 259 ! *********************************************************************** 260 261 DBG_ENTER("COLL") 262 263 call timab(503,1,tsec) 264 265 comm_cell = mpi_enreg%comm_cell 266 267 timrev = 1 ! as q=0 268 cplex = 2 - timrev 269 nspden = dtset%nspden 270 ecut_eff = (dtset%ecut)*(dtset%dilatmx)**2 271 mpsang = psps%mpsang 272 optorth=1;if (psps%usepaw==1) optorth=0 273 274 qphon(:)=zero 275 276 ABI_MALLOC(cg1,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol)) 277 ABI_MALLOC(cg2,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol)) 278 ABI_MALLOC(cg3,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol)) 279 ABI_MALLOC(eigen1,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol)) 280 ABI_MALLOC(eigen2,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol)) 281 ABI_MALLOC(eigen3,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol)) 282 ABI_MALLOC(rho1r1,(cplex*nfftf,dtset%nspden)) 283 ABI_MALLOC(rho2r1,(cplex*nfftf,dtset%nspden)) 284 ABI_MALLOC(rho2g1,(2,nfftf)) 285 ABI_MALLOC(rho3r1,(cplex*nfftf,dtset%nspden)) 286 287 ask_accurate=1 ; formeig = 1 ; ireadwf = 1 288 n1=ngfftf(1) ; n2=ngfftf(2) ; n3=ngfftf(3) 289 nfftotf=n1*n2*n3 290 291 !==== Initialize most of the Hamiltonian (and derivative) ==== 292 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 293 !2) Perform the setup needed for the non-local factors: 294 !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk. 295 !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients. 296 call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,natom,& 297 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,& 298 & usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,paw_ij=paw_ij0,& 299 & gpu_option=dtset%gpu_option) 300 301 ABI_MALLOC(vpsp1,(cplex*nfftf)) 302 ABI_MALLOC(xccc3d1,(cplex*nfftf)) 303 ABI_MALLOC(xccc3d2,(cplex*nfftf)) 304 ABI_MALLOC(xccc3d3,(cplex*nfftf)) 305 ABI_MALLOC(vhartr1_i2pert,(cplex*nfftf)) 306 ABI_MALLOC(vxc1_i2pert,(cplex*nfftf,dtset%nspden)) 307 ABI_MALLOC(vtrial1_i2pert,(cplex*nfftf,dtset%nspden)) 308 309 ABI_MALLOC(vresid_dum,(0,0)) 310 ! PAW stuff 311 usexcnhat = 0 312 nhat1grdim=0 313 ABI_MALLOC(nhat1gr,(0,0,0)) 314 nhat1gr(:,:,:) = zero 315 rdwrpaw=psps%usepaw 316 !Allocate 1st-order PAW occupancies (rhoij1) 317 if (psps%usepaw==1) then 318 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,& 319 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc) 320 ABI_MALLOC(pawrhoij1_i1pert,(natom)) 321 ABI_MALLOC(pawrhoij1_i2pert,(natom)) 322 ABI_MALLOC(pawrhoij1_i3pert,(natom)) 323 call pawrhoij_nullify(pawrhoij1_i1pert) 324 call pawrhoij_nullify(pawrhoij1_i2pert) 325 call pawrhoij_nullify(pawrhoij1_i3pert) 326 call pawrhoij_alloc(pawrhoij1_i1pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,& 327 & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 328 call pawrhoij_alloc(pawrhoij1_i2pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,& 329 & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 330 call pawrhoij_alloc(pawrhoij1_i3pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,& 331 & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 332 else 333 ABI_MALLOC(pawrhoij1_i1pert,(0)) 334 ABI_MALLOC(pawrhoij1_i2pert,(0)) 335 ABI_MALLOC(pawrhoij1_i3pert,(0)) 336 end if 337 338 mcg=mpw*nspinor*mband*mkmem*nsppol 339 340 !Allocations/initializations for PAW only 341 if(psps%usepaw==1) then 342 usexcnhat=maxval(pawtab(:)%usexcnhat) 343 ! 1st-order compensation density 344 ABI_MALLOC(nhat1_i1pert,(cplex*nfftf,dtset%nspden)) 345 nhat1_i1pert=zero 346 ABI_MALLOC(nhat1_i2pert,(cplex*nfftf,dtset%nspden)) 347 nhat1_i2pert=zero 348 ABI_MALLOC(nhat1_i3pert,(cplex*nfftf,dtset%nspden)) 349 nhat1_i3pert=zero 350 351 ! 1st-order arrays/variables related to the PAW spheres 352 ABI_MALLOC(paw_an1_i2pert,(natom)) 353 ABI_MALLOC(paw_ij1_i2pert,(natom)) 354 call paw_an_nullify(paw_an1_i2pert) 355 call paw_ij_nullify(paw_ij1_i2pert) 356 357 has_dijfr=1 358 has_diju=merge(0,1,dtset%usepawu==0) 359 360 call paw_an_init(paw_an1_i2pert,dtset%natom,dtset%ntypat,0,0,dtset%nspden,cplex,dtset%pawxcdev,& 361 & dtset%typat,pawang,pawtab,has_vxc=1,& 362 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 363 364 call paw_ij_init(paw_ij1_i2pert,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,0,dtset%natom,& 365 & dtset%ntypat,dtset%typat,pawtab,& 366 & has_dij=1,has_dijhartree=1,has_dijfr=has_dijfr,has_dijU=has_diju,& 367 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 368 else 369 ABI_MALLOC(nhat1_i1pert,(0,0)) 370 ABI_MALLOC(nhat1_i2pert,(0,0)) 371 ABI_MALLOC(nhat1_i3pert,(0,0)) 372 ABI_MALLOC(paw_an1_i2pert,(0)) 373 ABI_MALLOC(paw_ij1_i2pert,(0)) 374 end if ! PAW 375 376 n3xccc=0;if(psps%n1xccc/=0)n3xccc=nfftf 377 non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 378 379 !Loop over the perturbations j1, j2, j3 380 381 pert1case = 0 ; pert2case = 0 ; pert3case = 0 382 383 do i1pert = 1, mpert 384 do i1dir = 1, 3 385 386 if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then 387 388 pert1case = i1dir + (i1pert-1)*3 389 counter = pert1case 390 call appdig(pert1case,dtfil%fnamewff1,fiwf1i) 391 392 call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,& 393 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 394 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 395 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 396 & dtset%nsppol,dtset%nsym,& 397 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 398 & dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i,wvl) 399 400 if (ireadwf==1) then 401 call WffClose (wff1,ierr) 402 end if 403 404 flag1 = 0 405 rho1r1(:,:) = zero 406 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 407 call appdig(pert1case,dtfil%fildens1in,fiden1i) 408 409 call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho1r1, & 410 hdr_den, pawrhoij1_i1pert, comm_cell, check_hdr=hdr) 411 call hdr_den%free() 412 end if 413 414 xccc3d1(:) = zero 415 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 416 ndir=1 417 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i1dir,i1pert,& 418 & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,& 419 & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 420 & atmrhor1=xccc3d1,optn_in=n3xccc/nfftf,optn2_in=1,optv_in=0,vspl=psps%vlspl) 421 else 422 ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp. 423 ! ------------------------------------------------------------------------------ 424 if(psps%n1xccc/=0)then 425 call dfpt_mkcore(cplex,i1dir,i1pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,& 426 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred) 427 end if ! psps%n1xccc/=0 428 end if ! usepaw 429 430 do i3pert = 1, mpert 431 do i3dir = 1, 3 432 433 if ((maxval(rfpert(i1dir,i1pert,:,:,i3dir,i3pert))==1)) then 434 435 pert3case = i3dir + (i3pert-1)*3 436 counter = 100*pert3case + pert1case 437 call appdig(pert3case,dtfil%fnamewff1,fiwf3i) 438 439 call inwffil(ask_accurate,cg3,dtset,dtset%ecut,ecut_eff,eigen3,dtset%exchn2n3d,& 440 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 441 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 442 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 443 & dtset%nsppol,dtset%nsym,& 444 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 445 & dtfil%unkg1,wff3,wfft3,dtfil%unwff3,& 446 & fiwf3i,wvl) 447 if (ireadwf==1) then 448 call WffClose (wff3,ierr) 449 end if 450 451 flag3 = 0 452 rho3r1(:,:) = zero 453 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 454 455 call appdig(pert3case,dtfil%fildens1in,fiden1i) 456 457 call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho3r1, & 458 hdr_den, pawrhoij1_i3pert, comm_cell, check_hdr=hdr) 459 call hdr_den%free() 460 end if 461 462 xccc3d3(:) = zero 463 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 464 ndir=1 465 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i3dir,i3pert,& 466 & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,& 467 & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 468 & atmrhor1=xccc3d3,optn_in=n3xccc/nfftf,optn2_in=1,optv_in=0,vspl=psps%vlspl) 469 else 470 ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp. 471 ! ------------------------------------------------------------------------------ 472 if(psps%n1xccc/=0)then 473 call dfpt_mkcore(cplex,i3dir,i3pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,& 474 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d3,xred) 475 end if ! psps%n1xccc/=0 476 end if ! usepaw 477 478 do i2pert = 1, mpert 479 do i2dir = 1, 3 480 481 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 482 483 blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 484 485 npert_phon = 0 486 if(i1pert<=dtset%natom) npert_phon = npert_phon + 1 487 if(i2pert<=dtset%natom) npert_phon = npert_phon + 1 488 if(i3pert<=dtset%natom) npert_phon = npert_phon + 1 489 if (npert_phon>1) then 490 ABI_ERROR("dfptnl_loop is available with at most one phonon perturbation. Change your input!") 491 end if 492 493 pert2case = i2dir + (i2pert-1)*3 494 counter = 100*pert2case + pert2case 495 call appdig(pert2case,dtfil%fnamewff1,fiwf2i) 496 497 call inwffil(ask_accurate,cg2,dtset,dtset%ecut,ecut_eff,eigen2,dtset%exchn2n3d,& 498 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 499 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 500 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 501 & dtset%nsppol,dtset%nsym,& 502 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 503 & dtfil%unkg1,wff2,wfft2,dtfil%unwff2,& 504 & fiwf2i,wvl) 505 if (ireadwf==1) then 506 call WffClose (wff2,ierr) 507 end if 508 509 ! Read the first-order densities from disk-files 510 rho2r1(:,:) = zero ; rho2g1(:,:) = zero 511 512 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 513 514 call appdig(pert2case,dtfil%fildens1in,fiden1i) 515 516 call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho2r1, & 517 hdr_den, pawrhoij1_i2pert , comm_cell, check_hdr=hdr) 518 call hdr_den%free() 519 520 ! Compute up+down rho1(G) by fft 521 ABI_MALLOC(work,(cplex*nfftf)) 522 work(:)=rho2r1(:,1) 523 call fourdp(cplex,rho2g1,work,-1,mpi_enreg,nfftf,1,ngfftf,0) 524 ABI_FREE(work) 525 526 end if 527 528 xccc3d2(:)=zero ; vpsp1(:)=zero 529 ! PAW: compute Vloc(1) and core(1) together in reciprocal space 530 ! -------------------------------------------------------------- 531 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 532 ndir=1 533 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i2dir,i2pert,& 534 & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,& 535 & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 536 & atmrhor1=xccc3d2,atmvlocr1=vpsp1,optn_in=n3xccc/nfftf,optn2_in=1,vspl=psps%vlspl) 537 ! PAW only: we sometimes have to compute 1st-order compensation density 538 ! and eventually add it to density from 1st-order WFs 539 ! ---------------------------------------------------------------------- 540 if (psps%usepaw==1) then 541 542 !Force the computation of nhatfr 543 do iatom=1,dtset%natom 544 pawfgrtab(iatom)%nhatfr_allocated = 0 545 pawfgrtab(iatom)%nhatfr = zero 546 end do 547 548 ! This portion of code works only when npert_phon<=1 549 if (i1pert<=natom.and.usexcnhat==0) then 550 call pawnhatfr(0,i1dir,i1pert,1,dtset%natom,nspden,psps%ntypat,& 551 & pawang,pawfgrtab(i1pert),pawrhoij(i1pert),pawtab,rprimd,& 552 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 553 end if 554 if (i2pert<=natom) then 555 call pawnhatfr(0,i2dir,i2pert,1,dtset%natom,nspden,psps%ntypat,& 556 & pawang,pawfgrtab(i2pert),pawrhoij(i2pert),pawtab,rprimd,& 557 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 558 end if 559 if (i3pert<=natom.and.usexcnhat==0) then 560 call pawnhatfr(0,i3dir,i3pert,1,dtset%natom,nspden,psps%ntypat,& 561 & pawang,pawfgrtab(i3pert),pawrhoij(i3pert),pawtab,rprimd,& 562 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 563 end if 564 565 if (usexcnhat==0) then 566 567 call pawmknhat(dummy_real,cplex,0,i1dir,i1pert,0,gprimd,natom,dtset%natom,& 568 & nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i1pert,& 569 & pawrhoij1_i1pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,& 570 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 571 if (flag1==0) then 572 rho1r1(:,:) = rho1r1(:,:) - nhat1_i1pert(:,:) 573 flag1 = 1 574 end if 575 576 call pawmknhat(dummy_real,cplex,0,i3dir,i3pert,0,gprimd,natom,dtset%natom,& 577 & nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i3pert,& 578 & pawrhoij1_i3pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,& 579 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 580 if (flag3==0) then 581 rho3r1(:,:) = rho3r1(:,:) - nhat1_i3pert(:,:) 582 flag3 = 1 583 end if 584 585 end if 586 587 call pawmknhat(dummy_real,cplex,0,i2dir,i2pert,0,gprimd,natom,dtset%natom,& 588 & nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i2pert,& 589 & pawrhoij1_i2pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,& 590 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 591 592 end if 593 594 else 595 596 ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp. 597 ! ------------------------------------------------------------------------------ 598 if(psps%n1xccc/=0)then 599 call dfpt_mkcore(cplex,i2dir,i2pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,& 600 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d2,xred) 601 end if ! psps%n1xccc/=0 602 603 call dfpt_vlocal(atindx,cplex,gmet,gsqcut,i2dir,i2pert,mpi_enreg,psps%mqgrid_vl,dtset%natom,& 604 & nattyp,nfftf,ngfftf,psps%ntypat,n1,n2,n3,ph1df,psps%qgrid_vl,& 605 & dtset%qptn,ucvol,psps%vlspl,vpsp1,xred) 606 607 end if ! usepaw 608 609 option=1;optene=0 610 call dfpt_rhotov(cplex,dummy_real,dummy_real,dummy_real,dummy_real,dummy_real,& 611 & gsqcut,i2dir,i2pert,dtset%ixc,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,nhat,& 612 & nhat1_i2pert,nhat1gr,nhat1grdim,nkxc,nspden,n3xccc,non_magnetic_xc,optene,option,& 613 & dtset%qptn,rhog,rho2g1,rhor,rho2r1,rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1_i2pert,& 614 & vpsp1,vresid_dum,dummy_real,vtrial1_i2pert,vxc,vxc1_i2pert,xccc3d2,dtset%ixcrot) 615 616 if (psps%usepaw==1.and.usexcnhat==0) then 617 rho2r1(:,:) = rho2r1(:,:) - nhat1_i2pert(:,:) 618 end if 619 620 if (psps%usepaw==1)then 621 call paw_an_reset_flags(paw_an1_i2pert) ! Force the recomputation of on-site potentials 622 call paw_ij_reset_flags(paw_ij1_i2pert,all=.true.) ! Force the recomputation of Dij 623 optfr=0 624 call pawdijfr(gprimd,i2dir,i2pert,natom,natom,nfftf,ngfftf,nspden,nsppol,& 625 & psps%ntypat,optfr,paw_ij1_i2pert,pawang,pawfgrtab,pawrad,pawtab,cplex,qphon,& 626 & rprimd,ucvol,vpsp1,vtrial,vxc,xred,& 627 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 628 629 ! Computation of "on-site" first-order potentials, first-order densities 630 option=1 631 call pawdenpot(dummy_real,dummy_real2,dummy_real3,i2pert,dtset%ixc,natom,dtset%natom,& 632 & nspden,psps%ntypat,dtset%nucdipmom,& 633 & 0,option,paw_an1_i2pert,paw_an0,paw_ij1_i2pert,pawang,& 634 & dtset%pawprtvol,pawrad,pawrhoij1_i2pert,dtset%pawspnorb,pawtab,dtset%pawxcdev,& 635 & dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,ucvol,psps%znuclpsp, & 636 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 637 ! First-order Dij computation 638 ! call timab(561,1,tsec) 639 if (has_dijfr>0) then 640 !vpsp1 contribution to Dij already stored in frozen part of Dij 641 ABI_MALLOC(vtrial1_tmp,(cplex*nfftf,nspden)) 642 vtrial1_tmp=vtrial1_i2pert 643 do ii=1,min(nspden,2) 644 vtrial1_tmp(:,ii)=vtrial1_tmp(:,ii)-vpsp1(:) 645 end do 646 else 647 vtrial1_tmp => vtrial1_i2pert 648 end if 649 call pawdij(cplex,dtset%enunit,gprimd,i2pert,natom,dtset%natom,& 650 & nfftf,nfftotf,dtset%nspden,psps%ntypat,paw_an1_i2pert,paw_ij1_i2pert,pawang,& 651 & pawfgrtab,dtset%pawprtvol,pawrad,pawrhoij1_i2pert,dtset%pawspnorb,pawtab,& 652 & dtset%pawxcdev,qphon,dtset%spnorbscl,ucvol,dtset%cellcharge(1),vtrial1_tmp,vxc1_i2pert,xred,& 653 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 654 if (has_dijfr>0) then 655 ABI_FREE(vtrial1_tmp) 656 end if 657 call symdij(gprimd,indsy1,i2pert,natom,dtset%natom,nsym1,psps%ntypat,0,& 658 & paw_ij1_i2pert,pawang1,dtset%pawprtvol,pawtab,rprimd,symaf1,symrc1, & 659 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,& 660 & qphon=qphon) 661 ! call timab(561,2,tsec) 662 663 end if ! end usepaw section 664 665 nwffile = 1 666 file_index(1) = i2dir + 3*(i2pert-1) 667 fnamewff(1) = dtfil%fnamewff1 668 669 if (i2pert==natom+2) then 670 671 nwffile = 3 672 file_index(2) = i2dir+natom*3 673 fnamewff(2) = dtfil%fnamewffddk 674 ! As npert_phon<=1 and i2pert==natom+2, i1pert or i3pert is necessarly equal to natom+2 675 if (i3pert==natom+2) then 676 second_idir = i3dir 677 else if (i1pert==natom+2) then 678 second_idir = i1dir 679 else 680 ABI_BUG(" i1pert or i3pert is supposed to be equal to natom+2, which is not the case here.") 681 end if 682 call rf2_getidir(i2dir,second_idir,idir_dkde) 683 file_index(3) = idir_dkde+9+(dtset%natom+6)*3 684 fnamewff(3) = dtfil%fnamewffdkde 685 686 if (npert_phon==1.and.psps%usepaw==1.and.second_idir/=i2dir) then 687 nwffile = 5 688 file_index(4) = second_idir+natom*3 689 fnamewff(4) = dtfil%fnamewffddk 690 call rf2_getidir(second_idir,i2dir,idir_dkde) ! i2dir and second_idir are reversed 691 file_index(5) = idir_dkde+9+(dtset%natom+6)*3 692 fnamewff(5) = dtfil%fnamewffdkde 693 end if 694 695 end if 696 697 do ii=1,nwffile 698 call appdig(file_index(ii),fnamewff(ii),fiwfddk) 699 ! Checking the existence of data file 700 if (.not. file_exists(fiwfddk)) then 701 ! Trick needed to run Abinit test suite in netcdf mode. 702 if (file_exists(nctk_ncify(fiwfddk))) then 703 write(message,"(3a)")"- File: ",trim(fiwfddk),& 704 " does not exist but found netcdf file with similar name." 705 call wrtout(std_out,message,'COLL') 706 fiwfddk = nctk_ncify(fiwfddk) 707 end if 708 if (.not. file_exists(fiwfddk)) then 709 ABI_ERROR('Missing file: '//TRIM(fiwfddk)) 710 end if 711 end if 712 write(message,'(2a)')'-dfptnl_loop : read the wavefunctions from file: ',trim(fiwfddk) 713 call wrtout(std_out,message,'COLL') 714 call wrtout(ab_out,message,'COLL') 715 ! Note that the unit number for these files is 50,51,52 or 53 (dtfil%unddk=50) 716 call wfk_open_read(ddk_f(ii),fiwfddk,1,dtset%iomode,dtfil%unddk+(ii-1),mpi_enreg%comm_cell) 717 end do 718 719 ! Perform DFPT part of the 3dte calculation 720 call timab(513,1,tsec) 721 ! NOTE : eigen2 equals zero here 722 723 call dfptnl_pert(atindx,cg,cg1,cg2,cg3,cplex,dtfil,dtset,d3etot,eigen0,gs_hamkq,k3xc,indsy1,i1dir,& 724 & i2dir,i3dir,i1pert,i2pert,i3pert,kg,mband,mgfft,mkmem,mk1mem,mpert,mpi_enreg,& 725 & mpsang,mpw,natom,nattyp,nfftf,nfftotf,ngfftf,nkpt,nk3xc,nspden,nspinor,nsppol,nsym1,npwarr,occ,& 726 & pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawtab,pawrhoij,pawrhoij1_i1pert,pawrhoij1_i2pert,pawrhoij1_i3pert,& 727 & paw_an0,paw_an1_i2pert,paw_ij1_i2pert,ph1d,psps,rho1r1,rho2r1,rho3r1,& 728 & rprimd,symaf1,symrc1,ucvol,vtrial,vhartr1_i2pert,vtrial1_i2pert,vxc1_i2pert,& 729 & ddk_f,xccc3d1,xccc3d2,xccc3d3,xred,& 730 & d3etot_1,d3etot_2,d3etot_3,d3etot_4,d3etot_5,d3etot_6,d3etot_7,d3etot_8,d3etot_9) 731 call timab(513,2,tsec) 732 733 734 ! Eventually close the dot file 735 do ii=1,nwffile 736 call ddk_f(ii)%close() 737 end do 738 739 ! if (psps%usepaw==1) then 740 ! do ii=1,natom 741 ! pawfgrtab(ii)%nhatfr = zero 742 ! end do 743 ! end if 744 745 end if ! rfpert 746 end do ! i2dir 747 end do ! i2pert 748 749 end if ! rfpert 750 end do ! i3dir 751 end do ! i3pert 752 753 end if ! rfpert 754 end do ! i1dir 755 end do ! i1pert 756 757 !More memory cleaning 758 call gs_hamkq%free() 759 760 ABI_FREE(cg1) 761 ABI_FREE(cg2) 762 ABI_FREE(cg3) 763 ABI_FREE(eigen1) 764 ABI_FREE(eigen2) 765 ABI_FREE(eigen3) 766 ABI_FREE(rho1r1) 767 ABI_FREE(rho2r1) 768 ABI_FREE(rho2g1) 769 ABI_FREE(rho3r1) 770 ABI_FREE(nhat1gr) 771 ABI_FREE(vresid_dum) 772 ABI_FREE(vtrial1_i2pert) 773 ABI_FREE(vxc1_i2pert) 774 ABI_FREE(vhartr1_i2pert) 775 ABI_FREE(vpsp1) 776 ABI_FREE(xccc3d1) 777 ABI_FREE(xccc3d2) 778 ABI_FREE(xccc3d3) 779 if (psps%usepaw==1) then 780 call pawrhoij_free(pawrhoij1_i1pert) 781 call pawrhoij_free(pawrhoij1_i2pert) 782 call pawrhoij_free(pawrhoij1_i3pert) 783 ABI_FREE(nhat1_i1pert) 784 ABI_FREE(nhat1_i2pert) 785 ABI_FREE(nhat1_i3pert) 786 call paw_an_free(paw_an1_i2pert) 787 call paw_ij_free(paw_ij1_i2pert) 788 ABI_FREE(paw_an1_i2pert) 789 ABI_FREE(paw_ij1_i2pert) 790 end if 791 ABI_FREE(pawrhoij1_i1pert) 792 ABI_FREE(pawrhoij1_i2pert) 793 ABI_FREE(pawrhoij1_i3pert) 794 795 call timab(503,2,tsec) 796 797 DBG_EXIT("COLL") 798 799 end subroutine dfptnl_loop
ABINIT/m_dfptnl_loop [ Modules ]
NAME
m_dfptnl_loop
FUNCTION
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (LB) 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
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_dfptnl_loop 22 23 implicit none 24 25 private