TABLE OF CONTENTS
- ABINIT/m_dfptlw_loop
- ABINIT/m_dfptlw_loop/dfptlw_loop
- ABINIT/m_dfptlw_loop/dfptlw_typeIproc
- ABINIT/m_dfptlw_loop/read_1eig
ABINIT/m_dfptlw_loop [ Modules ]
NAME
m_dfptlw_loop
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 .
NOTES
PARENTS
CHILDREN
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 module m_dfptlw_loop 28 29 use defs_basis 30 use defs_wvltypes 31 use m_errors 32 use m_abicore 33 use m_hdr 34 use m_nctk 35 use m_wffile 36 use m_wfk 37 use m_dtset 38 use m_dtfil 39 40 use defs_datatypes, only : pseudopotential_type 41 use defs_abitypes, only : MPI_type 42 use m_time, only : timab 43 use m_io_tools, only : file_exists,iomode_from_fname,get_unit 44 use m_kg, only : getcut,getph 45 use m_inwffil, only : inwffil 46 use m_fft, only : fourdp 47 use m_ioarr, only : read_rhor 48 use m_hamiltonian, only : gs_hamiltonian_type, init_hamiltonian 49 use m_pawdij, only : pawdij, pawdijfr, symdij 50 use m_pawfgr, only : pawfgr_type 51 use m_pawfgrtab, only : pawfgrtab_type 52 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags 53 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags, paw_ij_print 54 use m_pawang, only : pawang_type 55 use m_pawrad, only : pawrad_type 56 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_nullify, & 57 & pawrhoij_io, pawrhoij_inquire_dim 58 use m_paw_nhat, only : pawmknhat,pawnhatfr 59 use m_paw_denpot, only : pawdenpot 60 use m_pawtab, only : pawtab_type 61 use m_rf2, only : rf2_getidir 62 use m_initylmg, only : initylmg 63 use m_atm2fft, only : dfpt_atm2fft 64 use m_dfpt_mkvxc, only : dfpt_mkvxc 65 use m_dfpt_rhotov, only : dfpt_rhotov 66 use m_mkcore, only : dfpt_mkcore 67 use m_mklocl, only : dfpt_vlocal, vlocalstr,dfpt_vlocaldq,dfpt_vlocaldqdq,dfpt_vmetdqdq 68 use m_dfptlw_pert, only : dfptlw_pert 69 use m_dynmat, only : cart39 70 use m_xmpi 71 72 implicit none 73 74 private
ABINIT/m_dfptlw_loop/dfptlw_loop [ Functions ]
NAME
dfptlw_loop
FUNCTION
Loop over two perturbations j1, j2 and a q gradient
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave coefficients of wavefunctions d3e_pert1(mpert)=array with the i1pert cases to calculate d3e_pert2(mpert)=array with the i2pert cases to calculate dimffnl= third dimension of ffnl dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset 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 gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1) kg(3,mpw*mkmem)=reduced planewave coordinates kxc(nfftf,nkxc)=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) nkpt = number of k points nkxc=second dimension of the array kxc, see rhotoxc.f for a description 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 pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data 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. rmet(3,3)=real space metric tensor in bohr**2 rprimd(3,3)=dimensional primitive translations (bohr) ucvol = unit cell volume (bohr^3) useylmgr= if 1 use the derivative of spherical harmonics xred(3,natom) = reduced atomic coordinates 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
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
SIDE EFFECTS
hdr <type(hdr_type)>=the header of wf, den and pot files
PARENTS
m_longwave
CHILDREN
SOURCE
155 subroutine dfptlw_loop(atindx,blkflg,cg,d3e_pert1,d3e_pert2,d3etot,dimffnl,dtfil,dtset,& 156 & ffnl,gmet,gprimd,& 157 & hdr,kg,kxc,mband,mgfft,mkmem,mk1mem,& 158 & mpert,mpi_enreg,mpw,natom,nattyp,ngfftf,nfftf,nkpt,nkxc,nspinor,nsppol,& 159 & npwarr,nylmgr,occ,& 160 & pawfgr,pawtab,& 161 & psps,rfpert,rhog,rhor,rmet,rprimd,ucvol,useylmgr,xred,ylm,ylmgr) 162 163 implicit none 164 165 !Arguments ------------------------------------ 166 !scalars 167 integer,intent(in) :: dimffnl,mband,mgfft,mk1mem,mkmem,mpert,mpw,natom,nfftf 168 integer,intent(in) :: nkpt,nkxc,nspinor,nsppol,nylmgr,useylmgr 169 real(dp),intent(in) :: ucvol 170 type(MPI_type),intent(inout) :: mpi_enreg 171 type(datafiles_type),intent(in) :: dtfil 172 type(dataset_type),intent(in) :: dtset 173 type(hdr_type),intent(inout) :: hdr 174 type(pawfgr_type),intent(in) :: pawfgr 175 type(pseudopotential_type),intent(in) :: psps 176 177 !arrays 178 integer,intent(in) :: atindx(natom),d3e_pert1(mpert),d3e_pert2(mpert) 179 integer,intent(in) :: kg(3,mk1mem*mpw) 180 integer,intent(in) :: nattyp(psps%ntypat),ngfftf(18),npwarr(nkpt) 181 integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert) 182 integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) 183 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gmet(3,3) 184 real(dp),intent(in) :: ffnl(mkmem,mpw,dimffnl,psps%lmnmax,psps%ntypat) 185 real(dp),intent(in) :: gprimd(3,3),kxc(nfftf,nkxc) 186 real(dp),intent(in) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden),rmet(3,3),rprimd(3,3) 187 real(dp),intent(in) :: xred(3,natom) 188 real(dp),intent(inout) :: occ(mband*nkpt*nsppol) 189 real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert) 190 real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm) 191 real(dp),intent(in) :: ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 192 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 193 194 !Local variables------------------------------- 195 !scalars 196 integer :: alpha,ask_accurate,beta,comm_cell,cplex,delta,dkdk_index,formeig,gamma 197 integer :: ia1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,idir_dkdk 198 integer :: idq,ierr,ii,ireadwf,istr,itypat,mcg,me,mpsang 199 integer :: n1,n2,n3,n1dq,n2dq,nhat1grdim,nfftotf,nspden,n3xccc 200 integer :: optgeom,opthartdqdq,optorth,pawread 201 integer :: pert1case,pert2case,pert3case,timrev,usexcnhat 202 real(dp) :: boxcut,delad,delag,delbd,delbg,ecut,ecut_eff,gsqcut 203 logical :: samepert 204 character(len=500) :: message 205 character(len=fnlen) :: fiden1i,fiwf1i,fiwf2i,fiwfddk,fiwfdkdk 206 type(gs_hamiltonian_type) :: gs_hamkq 207 type(wffile_type) :: wff1,wff2,wfft1,wfft2 208 type(wfk_t) :: ddk_f,d2_dkdk_f,d2_dkdk_f2 209 type(wvl_data) :: wvl 210 type(hdr_type) :: hdr_den 211 !arrays 212 integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/) 213 real(dp),allocatable :: cg1(:,:),cg2(:,:) 214 real(dp),allocatable :: d3etot_t4(:,:),d3etot_t5(:,:),d3etot_tgeom(:,:),eigen1(:),eigen2(:) 215 real(dp),allocatable :: nhat1(:,:),ph1d(:,:) 216 real(dp),allocatable :: rho1g1(:,:),rho1r1(:,:) 217 real(dp),allocatable :: rho2g1(:,:),rho2r1(:,:) 218 real(dp),allocatable :: t4_typeI(:,:,:,:,:,:),t4_typeII(:,:,:,:,:,:,:) 219 real(dp),allocatable :: t5_typeI(:,:,:,:,:,:),t5_typeII(:,:,:,:,:,:,:) 220 real(dp),allocatable :: tgeom_typeI(:,:,:,:,:,:),tgeom_typeII(:,:,:,:,:,:,:) 221 real(dp),allocatable :: vhart1dqdq(:),vpsp1dqdq(:) 222 real(dp),allocatable :: vpsp1_i1pertdq(:,:,:),vpsp1_i2pertdq(:,:,:) 223 real(dp),allocatable :: vpsp1_i1pertdq_geom(:,:,:), vpsp1_i1pertdqdq(:,:,:) 224 real(dp),allocatable :: vxc1dqdq(:),work(:) 225 type(pawrhoij_type),allocatable :: pawrhoij_read(:) 226 227 228 ! ************************************************************************* 229 230 DBG_ENTER("COLL") 231 232 !Init parallelism 233 comm_cell=mpi_enreg%comm_cell 234 me=mpi_enreg%me_kpt 235 236 !Various initializations 237 timrev = 1 ! as q=0 238 cplex = 2 - timrev 239 nspden = dtset%nspden 240 ecut=dtset%ecut 241 ecut_eff = ecut*(dtset%dilatmx)**2 242 mpsang = psps%mpsang 243 optorth=1;if (psps%usepaw==1) optorth=0 244 opthartdqdq=1 245 246 247 ABI_MALLOC(cg1,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol)) 248 ABI_MALLOC(cg2,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol)) 249 ABI_MALLOC(eigen1,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol)) 250 ABI_MALLOC(eigen2,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol)) 251 ABI_MALLOC(rho1r1,(cplex*nfftf,dtset%nspden)) 252 ABI_MALLOC(rho2r1,(cplex*nfftf,dtset%nspden)) 253 ABI_MALLOC(rho1g1,(2,nfftf)) 254 ABI_MALLOC(rho2g1,(2,nfftf)) 255 256 ask_accurate=1 ; formeig = 1 ; ireadwf = 1 257 n1=ngfftf(1) ; n2=ngfftf(2) ; n3=ngfftf(3) 258 nfftotf=n1*n2*n3 259 260 !Allocations for type-I terms 261 ABI_MALLOC(t4_typeII,(2,3,mpert,3,mpert,3,mpert)) 262 t4_typeII(:,:,:,:,:,:,:)=zero 263 if (d3e_pert2(natom+3)==1.or.d3e_pert2(natom+4)==1) then 264 ABI_MALLOC(t4_typeI,(2,3,mpert,3,3,3)) 265 t4_typeI(:,:,:,:,:,:)=zero 266 end if 267 ABI_MALLOC(t5_typeII,(2,3,mpert,3,mpert,3,mpert)) 268 t5_typeII(:,:,:,:,:,:,:)=zero 269 if (d3e_pert1(natom+3)==1.or.d3e_pert1(natom+4)==1) then 270 ABI_MALLOC(t5_typeI,(2,3,mpert,3,3,3)) 271 t5_typeI(:,:,:,:,:,:)=zero 272 end if 273 ABI_MALLOC(tgeom_typeII,(2,3,mpert,3,mpert,3,mpert)) 274 tgeom_typeII(:,:,:,:,:,:,:)=zero 275 if (any(d3e_pert1(1:natom)==1).and.(d3e_pert2(natom+3)==1.or.d3e_pert2(natom+4)==1)) then 276 ABI_MALLOC(tgeom_typeI,(2,3,mpert,3,3,3)) 277 tgeom_typeI(:,:,:,:,:,:)=zero 278 end if 279 280 !Compute large sphere cut-off gsqcut 281 call getcut(boxcut,ecut,gmet,gsqcut,dtset%iboxcut,std_out,dtset%qptn,dtset%ngfft) 282 283 !Generate the 1-dimensional phases 284 ABI_MALLOC(ph1d,(2,3*(2*mgfft+1)*dtset%natom)) 285 call getph(atindx,dtset%natom,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),ph1d,xred) 286 287 !==== Initialize most of the Hamiltonian (and derivative) ==== 288 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 289 !2) Perform the setup needed for the non-local factors: 290 !3) Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk. 291 call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,dtset%nsppol,nspden,dtset%natom,& 292 & dtset%typat,xred,dtset%nfft,mgfft,dtset%ngfft,rprimd,dtset%nloalg,ph1d=ph1d,& 293 & gpu_option=dtset%gpu_option) 294 295 !Specific allocations for strain-gradient perturbation 296 if (dtset%lw_flexo==1.or.dtset%lw_flexo==2.or.dtset%lw_flexo==4) then 297 ABI_MALLOC(vhart1dqdq,(2*nfftf)) 298 ABI_MALLOC(vpsp1dqdq,(2*nfftf)) 299 ABI_MALLOC(vxc1dqdq,(2*nfftf)) 300 end if 301 302 !This is necessary to deactivate paw options in the dfpt_rhotov routine 303 ABI_MALLOC(pawrhoij_read,(0)) 304 usexcnhat=0 305 n3xccc=0 306 pawread=0 307 nhat1grdim=0 308 ABI_MALLOC(nhat1,(cplex*dtset%nfft,nspden)) 309 nhat1=zero 310 311 mcg=mpw*nspinor*mband*mkmem*nsppol 312 313 pert1case = 0 ; pert2case = 0 ; pert3case = 0 314 315 do i1pert = 1, mpert 316 do i1dir = 1, 3 317 318 if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then 319 320 pert1case = i1dir + (i1pert-1)*3 321 call appdig(pert1case,dtfil%fnamewff1,fiwf1i) 322 323 call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,& 324 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 325 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 326 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 327 & dtset%nsppol,dtset%nsym,& 328 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 329 & dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i,wvl) 330 331 332 if (ireadwf==1) then 333 call WffClose (wff1,ierr) 334 end if 335 336 call read_1eig(eigen1,formeig,mband,nkpt,nsppol,fiwf1i) 337 338 rho1r1(:,:) = zero; rho1g1(:,:) = zero 339 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 340 call appdig(pert1case,dtfil%fildens1in,fiden1i) 341 342 call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, psps%usepaw, mpi_enreg, rho1r1, & 343 hdr_den, pawrhoij_read, comm_cell, check_hdr=hdr) 344 call hdr_den%free() 345 end if 346 347 !Perform FFT rhor1 to rhog1 348 ABI_MALLOC(work,(cplex*nfftf)) 349 work(:)=rho1r1(:,1) 350 call fourdp(cplex,rho1g1,work,-1,mpi_enreg,dtset%nfft,1,dtset%ngfft,0) 351 ABI_FREE(work) 352 353 !Allocate the first-order gradient local potential 354 if (i1pert <= natom+3) then 355 n1dq=1 356 ABI_MALLOC(vpsp1_i1pertdq,(2*nfftf,dtset%nspden,n1dq)) 357 else if (i1pert == natom+4) then 358 n1dq=2 359 ABI_MALLOC(vpsp1_i1pertdq,(2*nfftf,dtset%nspden,n1dq)) 360 else 361 n1dq=1 362 end if 363 ABI_MALLOC(d3etot_t5,(2,n1dq)) 364 365 do i2pert = 1, mpert 366 do i2dir = 1, 3 367 368 if ((maxval(rfpert(i1dir,i1pert,i2dir,i2pert,:,:))==1)) then 369 370 pert2case = i2dir + (i2pert-1)*3 371 call appdig(pert2case,dtfil%fnamewff1,fiwf2i) 372 373 call inwffil(ask_accurate,cg2,dtset,dtset%ecut,ecut_eff,eigen2,dtset%exchn2n3d,& 374 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,& 375 & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,& 376 & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,& 377 & dtset%nsppol,dtset%nsym,& 378 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 379 & dtfil%unkg1,wff2,wfft2,dtfil%unwff2,fiwf2i,wvl) 380 381 if (ireadwf==1) then 382 call WffClose (wff2,ierr) 383 end if 384 385 call read_1eig(eigen2,formeig,mband,nkpt,nsppol,fiwf2i) 386 387 if (i1pert==i2pert.and.i1dir==i2dir) then 388 samepert=.true. 389 else 390 samepert=.false. 391 end if 392 393 rho2r1(:,:) = zero; rho2g1(:,:) = zero 394 if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then 395 call appdig(pert2case,dtfil%fildens1in,fiden1i) 396 397 call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, psps%usepaw, mpi_enreg, rho2r1, & 398 hdr_den, pawrhoij_read, comm_cell, check_hdr=hdr) 399 call hdr_den%free() 400 end if 401 402 if (.not.samepert) then 403 !Perform FFT rhor1 to rhog1 404 ABI_MALLOC(work,(cplex*nfftf)) 405 work(:)=rho2r1(:,1) 406 call fourdp(cplex,rho2g1,work,-1,mpi_enreg,dtset%nfft,1,dtset%ngfft,0) 407 ABI_FREE(work) 408 end if !samepert 409 410 !Allocate the first-order gradient local potential 411 if (i2pert <= natom+3) then 412 n2dq=1 413 ABI_MALLOC(vpsp1_i2pertdq,(2*nfftf,dtset%nspden,n2dq)) 414 else if (i2pert == natom+4) then 415 n2dq=2 416 ABI_MALLOC(vpsp1_i2pertdq,(2*nfftf,dtset%nspden,n2dq)) 417 else 418 n2dq=1 419 end if 420 ABI_MALLOC(d3etot_t4,(2,n2dq)) 421 ABI_MALLOC(d3etot_tgeom,(2,n2dq)) 422 423 !Calculate the first-order gradient local potential that enters the geometric term 424 ABI_MALLOC(vpsp1_i1pertdq_geom,(2*nfftf,dtset%nspden,3)) 425 if (i1pert <= natom .and. (i2pert == natom+3.or.i2pert == natom+4)) then 426 427 !calculate the second of the two first-gradient directions 428 do ii=1,3 429 call dfpt_vlocaldq(atindx,2,gmet,gsqcut,i1dir,i1pert,mpi_enreg, & 430 & psps%mqgrid_vl,dtset%natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3, & 431 & ph1d,ii,psps%qgrid_vl,dtset%qptn,ucvol,psps%vlspl,vpsp1_i1pertdq_geom(:,1,ii)) 432 end do 433 end if 434 435 !Allocate the second-gradient array 436 ABI_MALLOC(vpsp1_i1pertdqdq,(2*nfftf,dtset%nspden,n2dq)) 437 438 do i3pert = 1, mpert 439 do i3dir = 1, 3 440 441 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 442 443 blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 444 445 !Calculate local potentials for first-order gradient Hamiltonians 446 !gradient of i1pert: 447 if (i1pert<=natom) then 448 !Get q-gradient of first-order local part of the pseudopotential 449 call dfpt_vlocaldq(atindx,2,gmet,gsqcut,i1dir,i1pert,mpi_enreg, & 450 & psps%mqgrid_vl,dtset%natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3, & 451 & ph1d,i3dir,psps%qgrid_vl,dtset%qptn,ucvol,psps%vlspl,vpsp1_i1pertdq(:,1,1)) 452 453 if (i2pert == natom+3.or.i2pert == natom+4) then 454 gamma=i3dir 455 do idq= 1, n2dq 456 if (i2pert==natom+3) then 457 istr=i2dir 458 else 459 istr=idq*3+i2dir 460 endif 461 delta=idx(2*istr) 462 call dfpt_vlocaldqdq(atindx,2,gs_hamkq%gmet,gsqcut,i1dir,i1pert,mpi_enreg, & 463 & psps%mqgrid_vl,dtset%natom,& 464 & nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3, & 465 & ph1d,gamma,delta,psps%qgrid_vl,& 466 & dtset%qptn,ucvol,psps%vlspl,vpsp1_i1pertdqdq(:,1,idq)) 467 end do 468 end if 469 470 else if (i1pert==natom+3.or.i1pert==natom+4) then 471 istr=i1dir; if (i1pert==natom+4) istr=3+i1dir 472 !Get 2nd q-gradient of first-order local part of the pseudopotential and of the Hartree 473 !(and XC if GGA) contribution from ground state density 474 call dfpt_vmetdqdq(2,gmet,gprimd,gsqcut,istr,i1pert,kxc,mpi_enreg, & 475 & psps%mqgrid_vl,natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3,& 476 & nkxc,nspden,opthartdqdq,ph1d,i3dir,psps%qgrid_vl,& 477 & dtset%qptn,rhog,rhor,ucvol,psps%vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq) 478 vpsp1_i1pertdq(:,1,1)=vhart1dqdq(:)+vpsp1dqdq(:)+vxc1dqdq(:) 479 if (i1pert==natom+4) then 480 !Here we need to calculate both extradiagonal shear-strains 481 !because the second gradient of the metric perturbation is 482 !type-I, i.e., non symmetric with respect to the 483 !permutation of the strain indexes. 484 istr=6+i1dir 485 call dfpt_vmetdqdq(2,gmet,gprimd,gsqcut,istr,i1pert,kxc,mpi_enreg, & 486 & psps%mqgrid_vl,natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3,& 487 & nkxc,nspden,opthartdqdq,ph1d,i3dir,psps%qgrid_vl,& 488 & dtset%qptn,rhog,rhor,ucvol,psps%vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq) 489 vpsp1_i1pertdq(:,1,2)=vhart1dqdq(:)+vpsp1dqdq(:)+vxc1dqdq(:) 490 end if 491 end if 492 493 if (.not.samepert) then 494 !gradient of i2pert: 495 if (i2pert<=natom) then 496 !Get q-gradient of first-order local part of the pseudopotential 497 call dfpt_vlocaldq(atindx,2,gmet,gsqcut,i2dir,i2pert,mpi_enreg, & 498 & psps%mqgrid_vl,dtset%natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3, & 499 & ph1d,i3dir,psps%qgrid_vl,dtset%qptn,ucvol,psps%vlspl,vpsp1_i2pertdq(:,1,1)) 500 else if (i2pert==natom+3.or.i2pert==natom+4) then 501 istr=i2dir; if (i2pert==natom+4) istr=3+i2dir 502 !Get 2nd q-gradient of first-order local part of the pseudopotential and of the Hartree 503 !(and XC if GGA) contribution from ground state density 504 call dfpt_vmetdqdq(2,gmet,gprimd,gsqcut,istr,i2pert,kxc,mpi_enreg, & 505 & psps%mqgrid_vl,natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3,& 506 & nkxc,nspden,opthartdqdq,ph1d,i3dir,psps%qgrid_vl,& 507 & dtset%qptn,rhog,rhor,ucvol,psps%vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq) 508 vpsp1_i2pertdq(:,1,1)=vhart1dqdq(:)+vpsp1dqdq(:)+vxc1dqdq(:) 509 if (i2pert==natom+4) then 510 !Here we need to calculate both extradiagonal shear-strains 511 !because the second gradient of the metric perturbation is 512 !type-I, i.e., non symmetric with respect to the 513 !permutation of the strain indexes. 514 istr=6+i2dir 515 call dfpt_vmetdqdq(2,gmet,gprimd,gsqcut,istr,i2pert,kxc,mpi_enreg, & 516 & psps%mqgrid_vl,natom,nattyp,dtset%nfft,dtset%ngfft,dtset%ntypat,n1,n2,n3,& 517 & nkxc,nspden,opthartdqdq,ph1d,i3dir,psps%qgrid_vl,& 518 & dtset%qptn,rhog,rhor,ucvol,psps%vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq) 519 vpsp1_i2pertdq(:,1,2)=vhart1dqdq(:)+vpsp1dqdq(:)+vxc1dqdq(:) 520 end if 521 end if 522 end if !samepert 523 524 !Prepare ddk wf file 525 pert3case = i3dir + natom*3 526 call appdig(pert3case,dtfil%fnamewffddk,fiwfddk) 527 ! Checking the existence of data file 528 if (.not. file_exists(fiwfddk)) then 529 ! Trick needed to run Abinit test suite in netcdf mode. 530 if (file_exists(nctk_ncify(fiwfddk))) then 531 write(message,"(3a)")"- File: ",trim(fiwfddk),& 532 " does not exist but found netcdf file with similar name." 533 call wrtout(std_out,message,'COLL') 534 fiwfddk = nctk_ncify(fiwfddk) 535 end if 536 if (.not. file_exists(fiwfddk)) then 537 ABI_ERROR('Missing file: '//TRIM(fiwfddk)) 538 end if 539 end if 540 write(message,'(2a)')'-dfptlw_loop : read the ddk wavefunctions from file: ',trim(fiwfddk) 541 call wrtout(std_out,message,'COLL') 542 !call wrtout(ab_out,message,'COLL') 543 ! Note that the unit number for these files is 50,51,52 or 53 (dtfil%unddk=50) 544 call wfk_open_read(ddk_f,fiwfddk,1,dtset%iomode,dtfil%unddk,mpi_enreg%comm_cell) 545 546 !Prepare d2_dkdk wf file 547 !For i1pert 548 if (i1pert==natom+2) then 549 call rf2_getidir(i1dir,i3dir,idir_dkdk) 550 !if (idir_dkdk>6) idir_dkdk=idir_dkdk-3 551 dkdk_index=idir_dkdk+(dtset%natom+6)*3 552 call appdig(dkdk_index,dtfil%fnamewffdkdk,fiwfdkdk) 553 !Check that d2_dkdk file exists and open it 554 if (.not. file_exists(fiwfdkdk)) then 555 ! Trick needed to run Abinit test suite in netcdf mode. 556 if (file_exists(nctk_ncify(fiwfdkdk))) then 557 write(message,"(3a)")"- File: ",trim(fiwfdkdk),& 558 " does not exist but found netcdf file with similar name." 559 call wrtout(std_out,message,'COLL') 560 fiwfdkdk = nctk_ncify(fiwfdkdk) 561 end if 562 if (.not. file_exists(fiwfdkdk)) then 563 ABI_ERROR('Missing file: '//TRIM(fiwfdkdk)) 564 end if 565 end if 566 write(message,'(2a)')'-dfptlw_loop : read the d2_dkdk wavefunctions from file: ',trim(fiwfdkdk) 567 call wrtout(std_out,message,'COLL') 568 !call wrtout(ab_out,message,'COLL') 569 call wfk_open_read(d2_dkdk_f,fiwfdkdk,1,dtset%iomode,dtfil%unddk+1,mpi_enreg%comm_cell) 570 571 end if 572 573 !Prepare d2_dkdk wf file 574 !For i2pert 575 if (i2pert==natom+2.and..not.samepert) then 576 call rf2_getidir(i2dir,i3dir,idir_dkdk) 577 !if (idir_dkdk>6) idir_dkdk=idir_dkdk-3 578 dkdk_index=idir_dkdk+(dtset%natom+6)*3 579 call appdig(dkdk_index,dtfil%fnamewffdkdk,fiwfdkdk) 580 !Check that d2_dkdk file exists and open it 581 if (.not. file_exists(fiwfdkdk)) then 582 ! Trick needed to run Abinit test suite in netcdf mode. 583 if (file_exists(nctk_ncify(fiwfdkdk))) then 584 write(message,"(3a)")"- File: ",trim(fiwfdkdk),& 585 " does not exist but found netcdf file with similar name." 586 call wrtout(std_out,message,'COLL') 587 fiwfdkdk = nctk_ncify(fiwfdkdk) 588 end if 589 if (.not. file_exists(fiwfdkdk)) then 590 ABI_ERROR('Missing file: '//TRIM(fiwfdkdk)) 591 end if 592 end if 593 write(message,'(2a)')'-dfptlw_loop : read the d2_dkdk wavefunctions from file: ',trim(fiwfdkdk) 594 call wrtout(std_out,message,'COLL') 595 !call wrtout(ab_out,message,'COLL') 596 call wfk_open_read(d2_dkdk_f2,fiwfdkdk,1,dtset%iomode,dtfil%unddk+2,mpi_enreg%comm_cell) 597 end if 598 599 !Perform the longwave DFPT part of the 3dte calculation 600 call dfptlw_pert(cg,cg1,cg2,cplex,d3etot,d3etot_t4,d3etot_t5,d3etot_tgeom,dimffnl,dtset, & 601 & eigen1,eigen2,ffnl,gmet,gs_hamkq,gsqcut,i1dir,& 602 & i2dir,i3dir,i1pert,i2pert,i3pert,kg,kxc,mband,mkmem,mk1mem,mpert,mpi_enreg,& 603 & mpsang,mpw,natom,n1dq,n2dq,nfftf,ngfftf,nkpt,nkxc,nspden,nspinor,nsppol,npwarr,nylmgr,occ,& 604 & pawfgr,psps,rho1g1,rho1r1,rho2r1,rmet,rprimd,samepert,ucvol,useylmgr,& 605 & vpsp1_i1pertdq,vpsp1_i1pertdqdq,vpsp1_i1pertdq_geom,vpsp1_i2pertdq,& 606 & ddk_f,d2_dkdk_f,d2_dkdk_f2,ylm,ylmgr) 607 608 !close ddk file 609 call ddk_f%close() 610 611 !close d2_dkdk file (i1pert) 612 if (i1pert==natom+2) call d2_dkdk_f%close() 613 614 ! Close d2_dkdk file (i2pert) 615 if (i2pert==natom+2.and..not.samepert) call d2_dkdk_f2%close() 616 617 !Save the type-I terms 618 if (i2pert==natom+3.or.i2pert==natom+4) then 619 gamma=i3dir 620 do idq=1,n2dq 621 if (i2pert==natom+3) then 622 istr=i2dir 623 else 624 istr=idq*3+i2dir 625 endif 626 beta=idx(2*istr-1); delta=idx(2*istr) 627 t4_typeI(:,i1dir,i1pert,beta,delta,gamma)=d3etot_t4(:,idq) 628 end do 629 else 630 t4_typeII(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=d3etot_t4(:,1) 631 end if 632 633 if (i1pert==natom+3.or.i1pert==natom+4) then 634 gamma=i3dir 635 do idq=1,n1dq 636 if (i1pert==natom+3) then 637 istr=i1dir 638 else 639 istr=idq*3+i1dir 640 endif 641 beta=idx(2*istr-1); delta=idx(2*istr) 642 t5_typeI(:,i2dir,i2pert,beta,delta,gamma)=d3etot_t5(:,idq) 643 end do 644 else 645 t5_typeII(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=d3etot_t5(:,1) 646 end if 647 648 649 if (i1pert<=natom.and.(i2pert==natom+3.or.i2pert==natom+4)) then 650 alpha=i1dir 651 gamma=i3dir 652 do idq=1,n2dq 653 if (i2pert==natom+3) then 654 istr=i2dir 655 else 656 istr=idq*3+i2dir 657 endif 658 beta=idx(2*istr-1); delta=idx(2*istr) 659 tgeom_typeI(:,i1dir,i1pert,beta,delta,gamma)=d3etot_tgeom(:,idq) 660 661 !Incorporate here the G=0 contribution of the geometric term 662 ia1=0 663 itypat=0 664 do ii=1,dtset%ntypat 665 ia1=ia1+nattyp(ii) 666 if (atindx(i1pert)<=ia1.and.itypat==0) itypat=ii 667 end do 668 delad=zero ; if (alpha==delta) delad=one 669 delbd=zero ; if (beta==delta) delbd=one 670 delag=zero ; if (alpha==gamma) delag=one 671 delbg=zero ; if (beta==gamma) delbg=one 672 673 tgeom_typeI(1,i1dir,i1pert,beta,delta,gamma)= & 674 & tgeom_typeI(1,i1dir,i1pert,beta,delta,gamma) + & 675 & pi*pi*rhog(1,1)*psps%vlspl(1,2,itypat)*(delag*delbd+delad*delbg) 676 end do 677 else 678 tgeom_typeII(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=d3etot_tgeom(:,1) 679 end if 680 681 end if ! rfpert 682 end do ! ir3dir 683 end do ! ir3pert 684 685 ABI_SFREE(vpsp1_i2pertdq) 686 ABI_FREE(vpsp1_i1pertdq_geom) 687 ABI_FREE(vpsp1_i1pertdqdq) 688 ABI_FREE(d3etot_t4) 689 ABI_FREE(d3etot_tgeom) 690 691 end if ! rfpert 692 end do ! i2dir 693 end do ! i2pert 694 695 ABI_SFREE(vpsp1_i1pertdq) 696 ABI_FREE(d3etot_t5) 697 698 end if ! rfpert 699 end do ! i1dir 700 end do ! i1pert 701 702 !More memory cleaning 703 call gs_hamkq%free() 704 705 ABI_FREE(cg1) 706 ABI_FREE(cg2) 707 ABI_FREE(eigen1) 708 ABI_FREE(eigen2) 709 ABI_FREE(rho1r1) 710 ABI_FREE(rho2r1) 711 ABI_FREE(rho1g1) 712 ABI_FREE(rho2g1) 713 ABI_FREE(nhat1) 714 ABI_FREE(pawrhoij_read) 715 ABI_FREE(ph1d) 716 717 if (dtset%lw_flexo==1.or.dtset%lw_flexo==2.or.dtset%lw_flexo==4) then 718 ABI_FREE(vhart1dqdq) 719 ABI_FREE(vpsp1dqdq) 720 ABI_FREE(vxc1dqdq) 721 end if 722 723 !Treatment of T4 and T5 terms that have a q-gradient of a rf Hamiltonian 724 !they need to be converted to type-II for strain perturbation 725 if (d3e_pert2(natom+3)==1.or.d3e_pert2(natom+4)==1) then 726 optgeom=0 727 call dfptlw_typeIproc(blkflg,gprimd,optgeom,mpert,natom,rfpert,rprimd,t4_typeI,t4_typeII) 728 end if 729 730 if (d3e_pert1(natom+3)==1.or.d3e_pert1(natom+4)==1) then 731 optgeom=0 732 call dfptlw_typeIproc(blkflg,gprimd,optgeom,mpert,natom,rfpert,rprimd,t5_typeI,t5_typeII) 733 end if 734 735 !Tgeom has to be converted to type-II 736 !To do it we need to convert the three involve indexes to cartessian. Then, 737 !after type-II conversion the q-gradient index is back converted to reduced. 738 if (any(d3e_pert1(1:natom)==1).and.(d3e_pert2(natom+3)==1.or.d3e_pert2(natom+4)==1)) then 739 optgeom=1 740 call dfptlw_typeIproc(blkflg,gprimd,optgeom,mpert,natom,rfpert,rprimd,tgeom_typeI,tgeom_typeII) 741 end if 742 743 !Incorporate T4, T5 and Tgeom to d3etot 744 d3etot(:,:,:,:,:,:,:)= d3etot(:,:,:,:,:,:,:) + & 745 & t4_typeII(:,:,:,:,:,:,:) + & 746 & t5_typeII(:,:,:,:,:,:,:) + & 747 & tgeom_typeII(:,:,:,:,:,:,:) 748 749 !Anounce end of spatial-dispersion calculation 750 write(message, '(a,a,a,a)' ) ch10,ch10,& 751 & ' -- Spatial-dispersion 3rd-order derivatives completed -- ',ch10 752 call wrtout(std_out,message,'COLL') 753 call wrtout(ab_out,message,'COLL') 754 755 !Deallocations 756 ABI_FREE(t4_typeII) 757 ABI_FREE(t5_typeII) 758 ABI_FREE(tgeom_typeII) 759 ABI_SFREE(t4_typeI) 760 ABI_SFREE(t5_typeI) 761 ABI_SFREE(tgeom_typeI) 762 763 DBG_EXIT("COLL") 764 765 end subroutine dfptlw_loop
ABINIT/m_dfptlw_loop/dfptlw_typeIproc [ Functions ]
NAME
dfptlw_typeIproc
FUNCTION
Process type-I terms and convert them to type-II in the d3etot mixed (reduced/cartessian) coordinates.
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
blkflg(3,mpert,3,mpert) = flags for each element of the 3DTE gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1) optgeom= if 1 do special treatment for the geometric term mpert =maximum number of ipert natom = number of atoms in unit cell 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 rprimd(3,3)=dimensional primitive translations (bohr) t_typeI(2,3,mpert,3,3,3)= Input type-I tensor
OUTPUT
t_typeII(2,3,mpert,3,mpert,3,mpert)= type-II tensor converted to the mixed
coordinates.
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
SOURCE
809 #if defined HAVE_CONFIG_H 810 #include "config.h" 811 #endif 812 813 #include "abi_common.h" 814 815 816 subroutine dfptlw_typeIproc(blkflg,gprimd,optgeom,mpert,natom,rfpert,rprimd,t_typeI,& 817 & t_typeII) 818 819 use defs_basis 820 use m_errors 821 use m_profiling_abi 822 use m_dynmat, only : cart39 823 824 implicit none 825 826 !Arguments ------------------------------------ 827 !scalars 828 integer,intent(in) :: mpert,natom,optgeom 829 !arrays 830 integer,intent(in) :: blkflg(3,mpert,3,mpert,3,mpert) 831 integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert) 832 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3) 833 real(dp),intent(inout) :: t_typeI(2,3,mpert,3,3,3) 834 real(dp),intent(inout) :: t_typeII(2,3,mpert,3,mpert,3,mpert) 835 836 !Local variables------------------------------- 837 !scalar 838 integer :: beta,delta,gamma,ii 839 integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,istr 840 real(dp) :: fac 841 !arrays 842 integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/) 843 integer :: flg1(3),flg2(3) 844 real(dp) :: vec1(3),vec2(3) 845 846 ! ************************************************************************* 847 848 DBG_ENTER("COLL") 849 850 if (optgeom==1) then 851 !Transform the metric perturbation direction 852 !(treat it as an atomic displacement) 853 flg1(:)=1 854 do i1pert=1,natom 855 do i1dir=1,3 856 do gamma=1,3 857 do ii=1,2 858 do delta=1,3 859 do beta=1,3 860 vec1(beta)=t_typeI(ii,i1dir,i1pert,beta,delta,gamma) 861 end do 862 call cart39(flg1,flg2,gprimd,i1pert,natom,rprimd,vec1,vec2) 863 do beta=1,3 864 t_typeI(ii,i1dir,i1pert,beta,delta,gamma)=vec2(beta) 865 end do 866 end do 867 end do 868 end do 869 end do 870 end do 871 872 !Transform the second q-gradient direction 873 !(treat it as an electric field) 874 do i1pert=1,natom 875 do i1dir=1,3 876 do gamma=1,3 877 do ii=1,2 878 do beta=1,3 879 do delta=1,3 880 vec1(delta)=t_typeI(ii,i1dir,i1pert,beta,delta,gamma) 881 end do 882 call cart39(flg1,flg2,gprimd,natom+2,natom,rprimd,vec1,vec2) 883 do delta=1,3 884 t_typeI(ii,i1dir,i1pert,beta,delta,gamma)=vec2(delta) 885 end do 886 end do 887 end do 888 end do 889 end do 890 end do 891 892 !Transform the first q-gradient direction 893 !(treat it as an electric field) 894 do i1pert=1,natom 895 do i1dir=1,3 896 do ii=1,2 897 do beta=1,3 898 do delta=1,3 899 do gamma=1,3 900 vec1(gamma)=t_typeI(ii,i1dir,i1pert,beta,delta,gamma) 901 end do 902 call cart39(flg1,flg2,gprimd,natom+2,natom,rprimd,vec1,vec2) 903 do gamma=1,3 904 t_typeI(ii,i1dir,i1pert,beta,delta,gamma)=vec2(gamma) 905 end do 906 end do 907 end do 908 end do 909 end do 910 end do 911 912 end if 913 914 fac=two_pi ** 2 915 i3pert= natom+8 916 do i1pert = 1, mpert 917 do i1dir = 1, 3 918 if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then 919 do i2pert = natom+3, natom+4 920 do i2dir = 1, 3 921 istr=(i2pert-natom-3)*3+i2dir 922 beta=idx(2*istr-1); delta=idx(2*istr) 923 do ii=1,2 924 925 !Transform into type-II 926 do i3dir=1,3 927 gamma=i3dir 928 t_typeII(ii,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)= & 929 & t_typeI(ii,i1dir,i1pert,beta,delta,gamma) + & 930 & t_typeI(ii,i1dir,i1pert,delta,gamma,beta) - & 931 & t_typeI(ii,i1dir,i1pert,gamma,beta,delta) 932 end do ! i3dir 933 934 !Transform i3dir into reduced coordinates 935 do i3dir=1,3 936 vec1(i3dir)=t_typeII(ii,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 937 flg1(i3dir)=blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 938 end do 939 call cart39(flg1,flg2,transpose(rprimd),natom+2,natom,transpose(gprimd),vec1,vec2) 940 do i3dir=1,3 941 t_typeII(ii,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=vec2(i3dir)*fac 942 end do 943 944 end do ! ii 945 end do ! i2dir 946 end do ! i2pert 947 end if ! rfpert 948 end do ! i1dir 949 end do ! i1pert 950 951 DBG_EXIT("COLL") 952 953 end subroutine dfptlw_typeIproc
ABINIT/m_dfptlw_loop/read_1eig [ Functions ]
NAME
read_1eig
FUNCTION
Reads all the first-order energies from a given _1WF file. Data is read by the master and broadcasted.
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
formeig option (format of the eigenvalues and occupations) : 0 => ground-state format 1 => respfn format mband=maximum number of bands nkpt= number of k points nsppol=1 for unpolarized, 2 for spin-polarized wffnm=name (character data) of file for input wavefunctions.
OUTPUT
eigen(2*mband*mband*nkpt*nsppol)=matrix of eigenvalues
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
SOURCE
993 #if defined HAVE_CONFIG_H 994 #include "config.h" 995 #endif 996 997 #include "abi_common.h" 998 999 1000 subroutine read_1eig(eigen,formeig,mband,nkpt,nsppol,wffnm) 1001 1002 implicit none 1003 1004 !Arguments ------------------------------------ 1005 !scalars 1006 integer,intent(in) :: formeig,mband,nkpt,nsppol 1007 character(len=*),intent(inout) :: wffnm 1008 !arrays 1009 real(dp),intent(out) :: eigen((2*mband)**formeig*mband*nkpt*nsppol) 1010 1011 !Local variables------------------------------- 1012 !scalar 1013 integer :: bd2tot,comm,ierr,ik_bz,iomode,isppol,master,my_rank 1014 type(wfk_t) :: Wfk1 1015 type(hdr_type) :: hdr1 1016 character(len=500) :: msg 1017 !arrays 1018 real(dp),allocatable :: eig_buffer(:) 1019 1020 ! ************************************************************************* 1021 1022 DBG_ENTER("COLL") 1023 1024 comm = xmpi_world 1025 master = 0 1026 my_rank = xmpi_comm_rank(comm) 1027 1028 ! Master opens the 1WF file 1029 if (my_rank == master) then 1030 iomode = iomode_from_fname(wffnm) 1031 1032 !Check that atdis file exists and open it 1033 if (.not. file_exists(wffnm)) then 1034 ! Trick needed to run Abinit test suite in netcdf mode. 1035 if (file_exists(nctk_ncify(wffnm))) then 1036 write(msg,"(3a)")"- File: ",trim(wffnm),& 1037 " does not exist but found netcdf file with similar name." 1038 call wrtout(std_out,msg,'COLL') 1039 wffnm = nctk_ncify(wffnm) 1040 end if 1041 if (.not. file_exists(wffnm)) then 1042 ABI_ERROR('Missing file: '//TRIM(wffnm)) 1043 end if 1044 end if 1045 write(msg,'(a,a)')'-open 1wf file :',trim(wffnm) 1046 call wrtout(std_out,msg,'COLL') 1047 1048 call wfk_open_read(Wfk1, wffnm, formeig, iomode, get_unit(), xmpi_comm_self, Hdr_out=hdr1) 1049 end if 1050 1051 ! Master Broadcasts the header to all procs in comm 1052 call hdr1%bcast(master, my_rank, comm) 1053 1054 !Allocate buffer for MPI communicatio with max dimensions. 1055 ABI_MALLOC(eig_buffer,((2*mband)**formeig*mband*nsppol)) 1056 1057 bd2tot = 0 1058 do isppol=1,nsppol 1059 do ik_bz=1,hdr1%nkpt 1060 1061 !Master reads and broadcasts 1062 if (my_rank == master) then 1063 call Wfk1%read_band_block([1,mband], ik_bz, isppol, xmpio_single, & 1064 & eig_k=eig_buffer) 1065 end if 1066 1067 call xmpi_bcast(eig_buffer, master, comm, ierr) 1068 1069 eigen(1+bd2tot:2*mband**2+bd2tot)=eig_buffer(:) 1070 1071 !Keep track of total number of bands 1072 bd2tot = bd2tot + 2*mband**2 1073 1074 end do 1075 end do 1076 1077 if (my_rank == master) call Wfk1%close() 1078 call hdr1%free() 1079 1080 ABI_FREE(eig_buffer) 1081 1082 DBG_EXIT("COLL") 1083 1084 1085 end subroutine read_1eig