TABLE OF CONTENTS
ABINIT/dfptnl_loop [ Functions ]
NAME
dfptnl_loop
FUNCTION
Loop over the perturbations j1, j2 and j3
COPYRIGHT
Copyright (C) 2018-2022 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,use_gpu_cuda=dtset%use_gpu_cuda,paw_ij=paw_ij0) 299 300 ABI_MALLOC(vpsp1,(cplex*nfftf)) 301 ABI_MALLOC(xccc3d1,(cplex*nfftf)) 302 ABI_MALLOC(xccc3d2,(cplex*nfftf)) 303 ABI_MALLOC(xccc3d3,(cplex*nfftf)) 304 ABI_MALLOC(vhartr1_i2pert,(cplex*nfftf)) 305 ABI_MALLOC(vxc1_i2pert,(cplex*nfftf,dtset%nspden)) 306 ABI_MALLOC(vtrial1_i2pert,(cplex*nfftf,dtset%nspden)) 307 308 ABI_MALLOC(vresid_dum,(0,0)) 309 ! PAW stuff 310 usexcnhat = 0 311 nhat1grdim=0 312 ABI_MALLOC(nhat1gr,(0,0,0)) 313 nhat1gr(:,:,:) = zero 314 rdwrpaw=psps%usepaw 315 !Allocate 1st-order PAW occupancies (rhoij1) 316 if (psps%usepaw==1) then 317 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,& 318 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc) 319 ABI_MALLOC(pawrhoij1_i1pert,(natom)) 320 ABI_MALLOC(pawrhoij1_i2pert,(natom)) 321 ABI_MALLOC(pawrhoij1_i3pert,(natom)) 322 call pawrhoij_nullify(pawrhoij1_i1pert) 323 call pawrhoij_nullify(pawrhoij1_i2pert) 324 call pawrhoij_nullify(pawrhoij1_i3pert) 325 call pawrhoij_alloc(pawrhoij1_i1pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,& 326 & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 327 call pawrhoij_alloc(pawrhoij1_i2pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,& 328 & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 329 call pawrhoij_alloc(pawrhoij1_i3pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,& 330 & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 331 else 332 ABI_MALLOC(pawrhoij1_i1pert,(0)) 333 ABI_MALLOC(pawrhoij1_i2pert,(0)) 334 ABI_MALLOC(pawrhoij1_i3pert,(0)) 335 end if 336 337 mcg=mpw*nspinor*mband*mkmem*nsppol 338 339 !Allocations/initializations for PAW only 340 if(psps%usepaw==1) then 341 usexcnhat=maxval(pawtab(:)%usexcnhat) 342 ! 1st-order compensation density 343 ABI_MALLOC(nhat1_i1pert,(cplex*nfftf,dtset%nspden)) 344 nhat1_i1pert=zero 345 ABI_MALLOC(nhat1_i2pert,(cplex*nfftf,dtset%nspden)) 346 nhat1_i2pert=zero 347 ABI_MALLOC(nhat1_i3pert,(cplex*nfftf,dtset%nspden)) 348 nhat1_i3pert=zero 349 350 ! 1st-order arrays/variables related to the PAW spheres 351 ABI_MALLOC(paw_an1_i2pert,(natom)) 352 ABI_MALLOC(paw_ij1_i2pert,(natom)) 353 call paw_an_nullify(paw_an1_i2pert) 354 call paw_ij_nullify(paw_ij1_i2pert) 355 356 has_dijfr=1 357 has_diju=merge(0,1,dtset%usepawu==0) 358 359 call paw_an_init(paw_an1_i2pert,dtset%natom,dtset%ntypat,0,0,dtset%nspden,cplex,dtset%pawxcdev,& 360 & dtset%typat,pawang,pawtab,has_vxc=1,& 361 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 362 363 call paw_ij_init(paw_ij1_i2pert,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,0,dtset%natom,& 364 & dtset%ntypat,dtset%typat,pawtab,& 365 & has_dij=1,has_dijhartree=1,has_dijfr=has_dijfr,has_dijU=has_diju,& 366 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 367 else 368 ABI_MALLOC(nhat1_i1pert,(0,0)) 369 ABI_MALLOC(nhat1_i2pert,(0,0)) 370 ABI_MALLOC(nhat1_i3pert,(0,0)) 371 ABI_MALLOC(paw_an1_i2pert,(0)) 372 ABI_MALLOC(paw_ij1_i2pert,(0)) 373 end if ! PAW 374 375 n3xccc=0;if(psps%n1xccc/=0)n3xccc=nfftf 376 non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 377 378 !Loop over the perturbations j1, j2, j3 379 380 pert1case = 0 ; pert2case = 0 ; pert3case = 0 381 382 do i1pert = 1, mpert 383 do i1dir = 1, 3 384 385 if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then 386 387 pert1case = i1dir + (i1pert-1)*3 388 counter = pert1case 389 call appdig(pert1case,dtfil%fnamewff1,fiwf1i) 390 391 call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,& 392 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 393 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 394 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 395 & dtset%nsppol,dtset%nsym,& 396 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 397 & dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i,wvl) 398 399 if (ireadwf==1) then 400 call WffClose (wff1,ierr) 401 end if 402 403 flag1 = 0 404 rho1r1(:,:) = zero 405 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 406 call appdig(pert1case,dtfil%fildens1in,fiden1i) 407 408 call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho1r1, & 409 hdr_den, pawrhoij1_i1pert, comm_cell, check_hdr=hdr) 410 call hdr_den%free() 411 end if 412 413 xccc3d1(:) = zero 414 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 415 ndir=1 416 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i1dir,i1pert,& 417 & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,& 418 & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 419 & atmrhor1=xccc3d1,optn_in=n3xccc/nfftf,optn2_in=1,optv_in=0,vspl=psps%vlspl) 420 else 421 ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp. 422 ! ------------------------------------------------------------------------------ 423 if(psps%n1xccc/=0)then 424 call dfpt_mkcore(cplex,i1dir,i1pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,& 425 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred) 426 end if ! psps%n1xccc/=0 427 end if ! usepaw 428 429 do i3pert = 1, mpert 430 do i3dir = 1, 3 431 432 if ((maxval(rfpert(i1dir,i1pert,:,:,i3dir,i3pert))==1)) then 433 434 pert3case = i3dir + (i3pert-1)*3 435 counter = 100*pert3case + pert1case 436 call appdig(pert3case,dtfil%fnamewff1,fiwf3i) 437 438 call inwffil(ask_accurate,cg3,dtset,dtset%ecut,ecut_eff,eigen3,dtset%exchn2n3d,& 439 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 440 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 441 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 442 & dtset%nsppol,dtset%nsym,& 443 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 444 & dtfil%unkg1,wff3,wfft3,dtfil%unwff3,& 445 & fiwf3i,wvl) 446 if (ireadwf==1) then 447 call WffClose (wff3,ierr) 448 end if 449 450 flag3 = 0 451 rho3r1(:,:) = zero 452 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 453 454 call appdig(pert3case,dtfil%fildens1in,fiden1i) 455 456 call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho3r1, & 457 hdr_den, pawrhoij1_i3pert, comm_cell, check_hdr=hdr) 458 call hdr_den%free() 459 end if 460 461 xccc3d3(:) = zero 462 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 463 ndir=1 464 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i3dir,i3pert,& 465 & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,& 466 & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 467 & atmrhor1=xccc3d3,optn_in=n3xccc/nfftf,optn2_in=1,optv_in=0,vspl=psps%vlspl) 468 else 469 ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp. 470 ! ------------------------------------------------------------------------------ 471 if(psps%n1xccc/=0)then 472 call dfpt_mkcore(cplex,i3dir,i3pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,& 473 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d3,xred) 474 end if ! psps%n1xccc/=0 475 end if ! usepaw 476 477 do i2pert = 1, mpert 478 do i2dir = 1, 3 479 480 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 481 482 blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 483 484 npert_phon = 0 485 if(i1pert<=dtset%natom) npert_phon = npert_phon + 1 486 if(i2pert<=dtset%natom) npert_phon = npert_phon + 1 487 if(i3pert<=dtset%natom) npert_phon = npert_phon + 1 488 if (npert_phon>1) then 489 ABI_ERROR("dfptnl_loop is available with at most one phonon perturbation. Change your input!") 490 end if 491 492 pert2case = i2dir + (i2pert-1)*3 493 counter = 100*pert2case + pert2case 494 call appdig(pert2case,dtfil%fnamewff1,fiwf2i) 495 496 call inwffil(ask_accurate,cg2,dtset,dtset%ecut,ecut_eff,eigen2,dtset%exchn2n3d,& 497 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 498 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 499 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 500 & dtset%nsppol,dtset%nsym,& 501 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 502 & dtfil%unkg1,wff2,wfft2,dtfil%unwff2,& 503 & fiwf2i,wvl) 504 if (ireadwf==1) then 505 call WffClose (wff2,ierr) 506 end if 507 508 ! Read the first-order densities from disk-files 509 rho2r1(:,:) = zero ; rho2g1(:,:) = zero 510 511 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 512 513 call appdig(pert2case,dtfil%fildens1in,fiden1i) 514 515 call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho2r1, & 516 hdr_den, pawrhoij1_i2pert , comm_cell, check_hdr=hdr) 517 call hdr_den%free() 518 519 ! Compute up+down rho1(G) by fft 520 ABI_MALLOC(work,(cplex*nfftf)) 521 work(:)=rho2r1(:,1) 522 call fourdp(cplex,rho2g1,work,-1,mpi_enreg,nfftf,1,ngfftf,0) 523 ABI_FREE(work) 524 525 end if 526 527 xccc3d2(:)=zero ; vpsp1(:)=zero 528 ! PAW: compute Vloc(1) and core(1) together in reciprocal space 529 ! -------------------------------------------------------------- 530 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 531 ndir=1 532 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i2dir,i2pert,& 533 & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,& 534 & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 535 & atmrhor1=xccc3d2,atmvlocr1=vpsp1,optn_in=n3xccc/nfftf,optn2_in=1,vspl=psps%vlspl) 536 ! PAW only: we sometimes have to compute 1st-order compensation density 537 ! and eventually add it to density from 1st-order WFs 538 ! ---------------------------------------------------------------------- 539 if (psps%usepaw==1) then 540 541 !Force the computation of nhatfr 542 do iatom=1,dtset%natom 543 pawfgrtab(iatom)%nhatfr_allocated = 0 544 pawfgrtab(iatom)%nhatfr = zero 545 end do 546 547 ! This portion of code works only when npert_phon<=1 548 if (i1pert<=natom.and.usexcnhat==0) then 549 call pawnhatfr(0,i1dir,i1pert,1,dtset%natom,nspden,psps%ntypat,& 550 & pawang,pawfgrtab(i1pert),pawrhoij(i1pert),pawtab,rprimd,& 551 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 552 end if 553 if (i2pert<=natom) then 554 call pawnhatfr(0,i2dir,i2pert,1,dtset%natom,nspden,psps%ntypat,& 555 & pawang,pawfgrtab(i2pert),pawrhoij(i2pert),pawtab,rprimd,& 556 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 557 end if 558 if (i3pert<=natom.and.usexcnhat==0) then 559 call pawnhatfr(0,i3dir,i3pert,1,dtset%natom,nspden,psps%ntypat,& 560 & pawang,pawfgrtab(i3pert),pawrhoij(i3pert),pawtab,rprimd,& 561 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 562 end if 563 564 if (usexcnhat==0) then 565 566 call pawmknhat(dummy_real,cplex,0,i1dir,i1pert,0,gprimd,natom,dtset%natom,& 567 & nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i1pert,& 568 & pawrhoij1_i1pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,& 569 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 570 if (flag1==0) then 571 rho1r1(:,:) = rho1r1(:,:) - nhat1_i1pert(:,:) 572 flag1 = 1 573 end if 574 575 call pawmknhat(dummy_real,cplex,0,i3dir,i3pert,0,gprimd,natom,dtset%natom,& 576 & nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i3pert,& 577 & pawrhoij1_i3pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,& 578 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 579 if (flag3==0) then 580 rho3r1(:,:) = rho3r1(:,:) - nhat1_i3pert(:,:) 581 flag3 = 1 582 end if 583 584 end if 585 586 call pawmknhat(dummy_real,cplex,0,i2dir,i2pert,0,gprimd,natom,dtset%natom,& 587 & nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i2pert,& 588 & pawrhoij1_i2pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,& 589 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 590 591 end if 592 593 else 594 595 ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp. 596 ! ------------------------------------------------------------------------------ 597 if(psps%n1xccc/=0)then 598 call dfpt_mkcore(cplex,i2dir,i2pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,& 599 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d2,xred) 600 end if ! psps%n1xccc/=0 601 602 call dfpt_vlocal(atindx,cplex,gmet,gsqcut,i2dir,i2pert,mpi_enreg,psps%mqgrid_vl,dtset%natom,& 603 & nattyp,nfftf,ngfftf,psps%ntypat,n1,n2,n3,ph1df,psps%qgrid_vl,& 604 & dtset%qptn,ucvol,psps%vlspl,vpsp1,xred) 605 606 end if ! usepaw 607 608 option=1;optene=0 609 call dfpt_rhotov(cplex,dummy_real,dummy_real,dummy_real,dummy_real,dummy_real,& 610 & gsqcut,i2dir,i2pert,dtset%ixc,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,nhat,& 611 & nhat1_i2pert,nhat1gr,nhat1grdim,nkxc,nspden,n3xccc,non_magnetic_xc,optene,option,& 612 & dtset%qptn,rhog,rho2g1,rhor,rho2r1,rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1_i2pert,& 613 & vpsp1,vresid_dum,dummy_real,vtrial1_i2pert,vxc,vxc1_i2pert,xccc3d2,dtset%ixcrot) 614 615 if (psps%usepaw==1.and.usexcnhat==0) then 616 rho2r1(:,:) = rho2r1(:,:) - nhat1_i2pert(:,:) 617 end if 618 619 if (psps%usepaw==1)then 620 call paw_an_reset_flags(paw_an1_i2pert) ! Force the recomputation of on-site potentials 621 call paw_ij_reset_flags(paw_ij1_i2pert,all=.true.) ! Force the recomputation of Dij 622 optfr=0 623 call pawdijfr(gprimd,i2dir,i2pert,natom,natom,nfftf,ngfftf,nspden,nsppol,& 624 & psps%ntypat,optfr,paw_ij1_i2pert,pawang,pawfgrtab,pawrad,pawtab,cplex,qphon,& 625 & rprimd,ucvol,vpsp1,vtrial,vxc,xred,& 626 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 627 628 ! Computation of "on-site" first-order potentials, first-order densities 629 option=1 630 call pawdenpot(dummy_real,dummy_real2,dummy_real3,i2pert,dtset%ixc,natom,dtset%natom,& 631 & nspden,psps%ntypat,dtset%nucdipmom,& 632 & 0,option,paw_an1_i2pert,paw_an0,paw_ij1_i2pert,pawang,& 633 & dtset%pawprtvol,pawrad,pawrhoij1_i2pert,dtset%pawspnorb,pawtab,dtset%pawxcdev,& 634 & dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp, & 635 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 636 ! First-order Dij computation 637 ! call timab(561,1,tsec) 638 if (has_dijfr>0) then 639 !vpsp1 contribution to Dij already stored in frozen part of Dij 640 ABI_MALLOC(vtrial1_tmp,(cplex*nfftf,nspden)) 641 vtrial1_tmp=vtrial1_i2pert 642 do ii=1,min(nspden,2) 643 vtrial1_tmp(:,ii)=vtrial1_tmp(:,ii)-vpsp1(:) 644 end do 645 else 646 vtrial1_tmp => vtrial1_i2pert 647 end if 648 call pawdij(cplex,dtset%enunit,gprimd,i2pert,natom,dtset%natom,& 649 & nfftf,nfftotf,dtset%nspden,psps%ntypat,paw_an1_i2pert,paw_ij1_i2pert,pawang,& 650 & pawfgrtab,dtset%pawprtvol,pawrad,pawrhoij1_i2pert,dtset%pawspnorb,pawtab,& 651 & dtset%pawxcdev,qphon,dtset%spnorbscl,ucvol,dtset%cellcharge(1),vtrial1_tmp,vxc1_i2pert,xred,& 652 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 653 if (has_dijfr>0) then 654 ABI_FREE(vtrial1_tmp) 655 end if 656 call symdij(gprimd,indsy1,i2pert,natom,dtset%natom,nsym1,psps%ntypat,0,& 657 & paw_ij1_i2pert,pawang1,dtset%pawprtvol,pawtab,rprimd,symaf1,symrc1, & 658 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,& 659 & qphon=qphon) 660 ! call timab(561,2,tsec) 661 662 end if ! end usepaw section 663 664 nwffile = 1 665 file_index(1) = i2dir + 3*(i2pert-1) 666 fnamewff(1) = dtfil%fnamewff1 667 668 if (i2pert==natom+2) then 669 670 nwffile = 3 671 file_index(2) = i2dir+natom*3 672 fnamewff(2) = dtfil%fnamewffddk 673 ! As npert_phon<=1 and i2pert==natom+2, i1pert or i3pert is necessarly equal to natom+2 674 if (i3pert==natom+2) then 675 second_idir = i3dir 676 else if (i1pert==natom+2) then 677 second_idir = i1dir 678 else 679 ABI_BUG(" i1pert or i3pert is supposed to be equal to natom+2, which is not the case here.") 680 end if 681 call rf2_getidir(i2dir,second_idir,idir_dkde) 682 file_index(3) = idir_dkde+9+(dtset%natom+6)*3 683 fnamewff(3) = dtfil%fnamewffdkde 684 685 if (npert_phon==1.and.psps%usepaw==1.and.second_idir/=i2dir) then 686 nwffile = 5 687 file_index(4) = second_idir+natom*3 688 fnamewff(4) = dtfil%fnamewffddk 689 call rf2_getidir(second_idir,i2dir,idir_dkde) ! i2dir and second_idir are reversed 690 file_index(5) = idir_dkde+9+(dtset%natom+6)*3 691 fnamewff(5) = dtfil%fnamewffdkde 692 end if 693 694 end if 695 696 do ii=1,nwffile 697 call appdig(file_index(ii),fnamewff(ii),fiwfddk) 698 ! Checking the existence of data file 699 if (.not. file_exists(fiwfddk)) then 700 ! Trick needed to run Abinit test suite in netcdf mode. 701 if (file_exists(nctk_ncify(fiwfddk))) then 702 write(message,"(3a)")"- File: ",trim(fiwfddk),& 703 " does not exist but found netcdf file with similar name." 704 call wrtout(std_out,message,'COLL') 705 fiwfddk = nctk_ncify(fiwfddk) 706 end if 707 if (.not. file_exists(fiwfddk)) then 708 ABI_ERROR('Missing file: '//TRIM(fiwfddk)) 709 end if 710 end if 711 write(message,'(2a)')'-dfptnl_loop : read the wavefunctions from file: ',trim(fiwfddk) 712 call wrtout(std_out,message,'COLL') 713 call wrtout(ab_out,message,'COLL') 714 ! Note that the unit number for these files is 50,51,52 or 53 (dtfil%unddk=50) 715 call wfk_open_read(ddk_f(ii),fiwfddk,1,dtset%iomode,dtfil%unddk+(ii-1),mpi_enreg%comm_cell) 716 end do 717 718 ! Perform DFPT part of the 3dte calculation 719 call timab(513,1,tsec) 720 ! NOTE : eigen2 equals zero here 721 722 call dfptnl_pert(atindx,cg,cg1,cg2,cg3,cplex,dtfil,dtset,d3etot,eigen0,gs_hamkq,k3xc,indsy1,i1dir,& 723 & i2dir,i3dir,i1pert,i2pert,i3pert,kg,mband,mgfft,mkmem,mk1mem,mpert,mpi_enreg,& 724 & mpsang,mpw,natom,nattyp,nfftf,nfftotf,ngfftf,nkpt,nk3xc,nspden,nspinor,nsppol,nsym1,npwarr,occ,& 725 & pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawtab,pawrhoij,pawrhoij1_i1pert,pawrhoij1_i2pert,pawrhoij1_i3pert,& 726 & paw_an0,paw_an1_i2pert,paw_ij1_i2pert,ph1d,psps,rho1r1,rho2r1,rho3r1,& 727 & rprimd,symaf1,symrc1,ucvol,vtrial,vhartr1_i2pert,vtrial1_i2pert,vxc1_i2pert,& 728 & ddk_f,xccc3d1,xccc3d2,xccc3d3,xred,& 729 & d3etot_1,d3etot_2,d3etot_3,d3etot_4,d3etot_5,d3etot_6,d3etot_7,d3etot_8,d3etot_9) 730 call timab(513,2,tsec) 731 732 733 ! Eventually close the dot file 734 do ii=1,nwffile 735 call ddk_f(ii)%close() 736 end do 737 738 ! if (psps%usepaw==1) then 739 ! do ii=1,natom 740 ! pawfgrtab(ii)%nhatfr = zero 741 ! end do 742 ! end if 743 744 end if ! rfpert 745 end do ! i2dir 746 end do ! i2pert 747 748 end if ! rfpert 749 end do ! i3dir 750 end do ! i3pert 751 752 end if ! rfpert 753 end do ! i1dir 754 end do ! i1pert 755 756 !More memory cleaning 757 call gs_hamkq%free() 758 759 ABI_FREE(cg1) 760 ABI_FREE(cg2) 761 ABI_FREE(cg3) 762 ABI_FREE(eigen1) 763 ABI_FREE(eigen2) 764 ABI_FREE(eigen3) 765 ABI_FREE(rho1r1) 766 ABI_FREE(rho2r1) 767 ABI_FREE(rho2g1) 768 ABI_FREE(rho3r1) 769 ABI_FREE(nhat1gr) 770 ABI_FREE(vresid_dum) 771 ABI_FREE(vtrial1_i2pert) 772 ABI_FREE(vxc1_i2pert) 773 ABI_FREE(vhartr1_i2pert) 774 ABI_FREE(vpsp1) 775 ABI_FREE(xccc3d1) 776 ABI_FREE(xccc3d2) 777 ABI_FREE(xccc3d3) 778 if (psps%usepaw==1) then 779 call pawrhoij_free(pawrhoij1_i1pert) 780 call pawrhoij_free(pawrhoij1_i2pert) 781 call pawrhoij_free(pawrhoij1_i3pert) 782 ABI_FREE(nhat1_i1pert) 783 ABI_FREE(nhat1_i2pert) 784 ABI_FREE(nhat1_i3pert) 785 call paw_an_free(paw_an1_i2pert) 786 call paw_ij_free(paw_ij1_i2pert) 787 ABI_FREE(paw_an1_i2pert) 788 ABI_FREE(paw_ij1_i2pert) 789 end if 790 ABI_FREE(pawrhoij1_i1pert) 791 ABI_FREE(pawrhoij1_i2pert) 792 ABI_FREE(pawrhoij1_i3pert) 793 794 call timab(503,2,tsec) 795 796 DBG_EXIT("COLL") 797 798 end subroutine dfptnl_loop
ABINIT/m_dfptnl_loop [ Modules ]
NAME
m_dfptnl_loop
FUNCTION
COPYRIGHT
Copyright (C) 2018-2022 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