TABLE OF CONTENTS
- ABINIT/dfpt_init_mag1
- ABINIT/dfpt_looppert
- ABINIT/dfpt_prtene
- ABINIT/eigen_meandege
- ABINIT/getcgqphase
- ABINIT/m_dfpt_loopert
ABINIT/dfpt_init_mag1 [ Functions ]
NAME
dfpt_init_mag1
FUNCTION
Initial guess of the first order magnetization/density for magnetic field perturbation. The first order magnetization is set so as to zero out the first order XC magnetic field, which should minimize the second order XC energy (without taking self-consistency into account).
INPUTS
ipert = perturbation type (works only for ipert==natom+5) idir = direction of the applied magnetic field cplex = complex or real first order density and magnetization nfft = dimension of the fft grid nspden= number of density matrix components nkxc = number of kxc components vxc0(nfft,nspden) = GS XC potential kxc0(nfft,nspden) = GS XC derivatives rhor0(nfft,nspden) = GS density matrix
OUTPUT
rhor1(cplex*nfft) = first order density magnetization guess
SOURCE
3233 subroutine dfpt_init_mag1(ipert,idir,rhor1,rhor0,cplex,nfft,nspden,vxc0,kxc0,nkxc) 3234 3235 !Arguments ------------------------------------ 3236 integer , intent(in) :: ipert,idir,cplex,nfft,nspden,nkxc 3237 real(dp), intent(in) :: vxc0(nfft,nspden),rhor0(nfft,nspden) 3238 real(dp), intent(in) :: kxc0(nfft,nkxc) 3239 real(dp), intent(out) :: rhor1(cplex*nfft,nspden) 3240 3241 !Local variables------------------------------- 3242 integer :: ipt 3243 real(dp) :: bxc0,bxc1 3244 real(dp) :: m1_norm,m0_norm 3245 real(dp) :: f_dot_m 3246 real(dp) :: mdir(3),fdir(3) 3247 3248 ! ************************************************************************* 3249 ABI_UNUSED(ipert) 3250 3251 if (nspden==2) then 3252 3253 if(cplex==1) then 3254 do ipt=1,nfft 3255 bxc1=half*(half*(kxc0(ipt,1)+kxc0(ipt,3))-kxc0(ipt,2)) ! d/dm Bxc 3256 !this overestimates the first order magnetization because of n1 not taken into account 3257 m1_norm=-half*(1/bxc1) 3258 rhor1(ipt,1)=zero ! rho_up+rho_dwn => charge density 3259 rhor1(ipt,2)=half*m1_norm ! rho_up=1/2(rho+m) => half*m 3260 end do 3261 else 3262 do ipt=1,cplex*nfft 3263 rhor1(ipt,:)=zero 3264 end do 3265 end if 3266 3267 else if(nspden==4) then 3268 3269 fdir=zero 3270 fdir(idir)= 1.0d0 3271 do ipt=1,nfft 3272 m0_norm=sqrt(rhor0(ipt,2)**2+rhor0(ipt,3)**2+rhor0(ipt,4)**2) 3273 mdir(1)=rhor0(ipt,2)/m0_norm 3274 mdir(2)=rhor0(ipt,3)/m0_norm 3275 mdir(3)=rhor0(ipt,4)/m0_norm 3276 f_dot_m=fdir(1)*mdir(1)+fdir(2)*mdir(2)+fdir(3)*mdir(3) ! projection of the field direction on m0 3277 3278 bxc1=half*(half*(kxc0(ipt,1)+kxc0(ipt,3))-kxc0(ipt,2)) ! d/dm Bxc 3279 m1_norm=(-half/bxc1)*f_dot_m ! get an estimate of the norm of m1 3280 3281 bxc0=-sqrt((half*(vxc0(ipt,1)-vxc0(ipt,2)))**2+vxc0(ipt,3)**2+vxc0(ipt,4)**2) 3282 if(cplex==1) then 3283 rhor1(ipt,1)=zero ! rho_up+rho_dwn => charge density 3284 rhor1(ipt,2)=m1_norm*mdir(1)-half*m0_norm/bxc0*(fdir(1)-f_dot_m*mdir(1)) ! m1x 3285 rhor1(ipt,3)=m1_norm*mdir(2)-half*m0_norm/bxc0*(fdir(2)-f_dot_m*mdir(2)) ! m1y 3286 rhor1(ipt,4)=m1_norm*mdir(3)-half*m0_norm/bxc0*(fdir(3)-f_dot_m*mdir(3)) ! m1z 3287 rhor1(ipt,:)=zero 3288 else 3289 rhor1(2*ipt-1,1)=zero ! Re rho_up+rho_dwn 3290 rhor1(2*ipt-1,2)=m1_norm*mdir(1)-half*m0_norm/bxc0*(fdir(1)-f_dot_m*mdir(1)) ! m1x 3291 rhor1(2*ipt-1,3)=m1_norm*mdir(2)-half*m0_norm/bxc0*(fdir(2)-f_dot_m*mdir(2)) ! m1x 3292 rhor1(2*ipt-1,4)=m1_norm*mdir(3)-half*m0_norm/bxc0*(fdir(3)-f_dot_m*mdir(3)) ! m1x 3293 rhor1(2*ipt ,1)=zero ! Im rho_up+rho_dwn 3294 rhor1(2*ipt ,2)=zero 3295 rhor1(2*ipt ,3)=zero 3296 rhor1(2*ipt ,4)=zero 3297 3298 rhor1(2*ipt-1,1)=zero; rhor1(2*ipt,1)=zero 3299 rhor1(2*ipt-1,2)=zero; rhor1(2*ipt,2)=zero 3300 rhor1(2*ipt-1,3)=zero; rhor1(2*ipt,3)=zero 3301 rhor1(2*ipt-1,4)=zero; rhor1(2*ipt,4)=zero 3302 3303 end if 3304 end do 3305 end if 3306 3307 end subroutine dfpt_init_mag1
ABINIT/dfpt_looppert [ Functions ]
NAME
dfpt_looppert
FUNCTION
Loop over perturbations
INPUTS
atindx(natom)=index table for atoms (see gstate.f) codvsn=code version cpus=cpu time limit in seconds dim_eigbrd=1 if eigbrd is to be computed dim_eig2nkq=1 if eig2nkq is to be computed doccde(mband*nkpt*nsppol)=derivative of occupancies wrt the energy dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset dyew(2,3,natom,3,natom)=Ewald part of the dynamical matrix dyfrlo(3,3,natom)=frozen wavefunctions local part of the dynamical matrix dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=frozen wavefunctions non-local part of the dynamical matrix dyfrx1(2,3,natom,3,natom)=frozen wf nonlin. xc core corr.(2) part of the dynamical matrix dyfrx2(3,3,natom)=frozen wf nonlin. xc core corr.(2) part of the dynamical matrix dyvdw(2,3,natom,3,natom*usevdw)=vdw DFT-D part of the dynamical matrix dyfr_cplex=1 if dyfrnl is real, 2 if it is complex dyfr_nondiag=1 if dyfrnl is non diagonal with respect to atoms; 0 otherwise eigbrd(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eigbrd)=broadening factors for the electronic eigenvalues eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq)=second derivatives of the electronic eigenvalues eltcore(6,6)=core contribution to the elastic tensor elteew(6+3*natom,6)=Ewald contribution to the elastic tensor eltfrhar(6,6)=Hartree contribution to the elastic tensor eltfrkin(6,6)=kinetic contribution to the elastic tensor eltfrloc(6+3*natom,6)=local psp contribution to the elastic tensor eltfrnl(6+3*natom,6)=non-local psp contribution to the elastic tensor eltfrxc(6+3*natom,6)=exchange-correlation contribution to the elastic tensor eltvdw(6+3*natom,6*usevdw)=vdw DFT-D part of the elastic tensor fermie=fermi energy (Hartree) iexit=index of "exit" on first line of file (0 if not found) indsym(4,nsym,natom)=indirect indexing array for atom labels kxc(nfftf,nkxc)=exchange and correlation kernel (see rhotoxc.f) mkmem =Number of k points treated by this node (GS data) mkqmem=Number of k+q points treated by this node (GS data) mk1mem=Number of k points treated by this node (RF data) mpert=maximum number of ipert mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor 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) nkpt=number of k points nkxc=second dimension of the kxc array nspden=number of spin-density components nsym=number of symmetry elements in space group occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh for the GS paw_ij(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 pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(my_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(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data pertsy(3,mpert)=set of perturbations that form a basis for all other perturbations prtbbb=if 1, bbb decomposition, also dimension d2bbb psps <type(pseudopotential_type)>=variables related to pseudopotentials rfpert(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. symq(4,2,nsym)=1 if symmetry preserves present qpoint. From littlegroup_q symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space) timrev=1 if time-reversal preserves the q wavevector; 0 otherwise. usecprj= 1 if cprj, cprjq, cprj1 arrays are stored in memory usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use vtrial(nfftf,nspden)=GS potential (Hartree) vxc(nfftf,nspden)=Exchange-Correlation GS potential (Hartree) vxcavg=average of vxc potential vxctau(nfftf,nspden,4*usekden)=derivative of e_xc with respect to kinetic energy density, for mGGA xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed) ddkfil(3)=unit numbers for the three possible ddk files d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some second order derivatives d2lo(2,mpert,3,mpert)=local contributions to the 2DTEs d2nl(2,mpert,3,mpert)=non-local contributions to the 2DTEs d2ovl(2,mpert,3,mpert*usepaw)=1st-order change of WF overlap contributions to the 2DTEs etotal=total energy (sum of 8 contributions) (hartree)
SOURCE
194 subroutine dfpt_looppert(atindx,blkflg,codvsn,cpus,dim_eigbrd,dim_eig2nkq,doccde,& 195 & ddkfil,dtfil,dtset,dyew,dyfrlo,dyfrnl,dyfrx1,dyfrx2,dyvdw,& 196 & dyfr_cplex,dyfr_nondiag,d2bbb,d2lo,d2nl,d2ovl,efmasdeg,efmasval,eigbrd,eig2nkq,& 197 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,& 198 & etotal,fermie,iexit,indsym,kxc,& 199 & mkmem,mkqmem,mk1mem,mpert,mpi_enreg,my_natom,nattyp,& 200 & nfftf,nhat,nkpt,nkxc,nspden,nsym,occ,& 201 & paw_an,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,& 202 & pertsy,prtbbb,psps,rfpert,rf2_dirs_from_rfpert_nl,rhog,rhor,symq,symrec,timrev,& 203 & usecprj,usevdw,vtrial,vxc,vxcavg,vxctau,xred,clflg,occ_rbz_pert,eigen0_pert,eigenq_pert,& 204 & eigen1_pert,nkpt_rbz,eigenq_fine,hdr_fine,hdr0) 205 206 !Arguments ------------------------------------ 207 integer, intent(in) :: dim_eigbrd,dim_eig2nkq,dyfr_cplex,dyfr_nondiag,mk1mem,mkmem,mkqmem,mpert 208 integer, intent(in) :: nfftf,nkpt,nkxc,nspden,nsym,prtbbb,timrev,usecprj,usevdw 209 integer, intent(out) :: iexit 210 integer, intent(inout) :: my_natom 211 real(dp), intent(in) :: cpus,vxcavg 212 real(dp), intent(inout) :: fermie 213 real(dp), intent(inout) :: etotal 214 character(len=8), intent(in) :: codvsn 215 type(MPI_type), intent(inout) :: mpi_enreg 216 type(datafiles_type), intent(in) :: dtfil 217 type(dataset_type), intent(in), target :: dtset 218 type(pawang_type),intent(in) :: pawang 219 type(pawfgr_type),intent(in) :: pawfgr 220 type(pseudopotential_type), intent(inout) :: psps 221 integer, intent(in) :: atindx(dtset%natom),indsym(4,nsym,dtset%natom) 222 integer, intent(in) :: nattyp(dtset%ntypat),pertsy(3,dtset%natom+6) 223 integer, intent(in) :: rfpert(mpert),rf2_dirs_from_rfpert_nl(3,3),symq(4,2,nsym),symrec(3,3,nsym) 224 integer, intent(out) :: ddkfil(3) 225 integer, intent(inout) :: blkflg(3,mpert,3,mpert) 226 integer, intent(out) :: clflg(3,mpert) 227 real(dp), intent(in) :: doccde(dtset%mband*nkpt*dtset%nsppol) 228 real(dp), intent(in) :: dyew(2,3,dtset%natom,3,dtset%natom) 229 real(dp), intent(in) :: dyfrlo(3,3,dtset%natom) 230 real(dp), intent(in) :: dyfrnl(dyfr_cplex,3,3,dtset%natom,1+(dtset%natom-1)*dyfr_nondiag) 231 real(dp), intent(in) :: dyfrx1(2,3,dtset%natom,3,dtset%natom) 232 real(dp), intent(in) :: dyfrx2(3,3,dtset%natom),dyvdw(2,3,dtset%natom,3,dtset%natom*usevdw) 233 real(dp), intent(in) :: eltcore(6,6),elteew(6+3*dtset%natom,6),eltfrhar(6,6) 234 real(dp), intent(in) :: eltfrkin(6,6),eltfrloc(6+3*dtset%natom,6) 235 real(dp), intent(in) :: eltfrnl(6+3*dtset%natom,6) 236 real(dp), intent(in) :: eltfrxc(6+3*dtset%natom,6),eltvdw(6+3*dtset%natom,6*usevdw) 237 real(dp), intent(in) :: kxc(nfftf,nkxc),nhat(nfftf,nspden) 238 real(dp), intent(in) :: occ(dtset%mband*nkpt*dtset%nsppol) 239 real(dp), intent(in) :: rhog(2,nfftf),rhor(nfftf,nspden),vxc(nfftf,nspden) 240 real(dp), intent(in) :: vtrial(nfftf,nspden) 241 real(dp), intent(inout) :: vxctau(nfftf,dtset%nspden,4*dtset%usekden) 242 real(dp), intent(inout) :: xred(3,dtset%natom) 243 real(dp), intent(inout) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)!vz_i 244 real(dp), intent(inout) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert) !vz_i 245 real(dp), intent(inout) :: d2ovl(2,3,mpert,3,mpert*psps%usepaw) !vz_i 246 real(dp), intent(out) :: eigbrd(2,dtset%mband*dtset%nsppol,nkpt,3,dtset%natom,3,dtset%natom*dim_eigbrd) 247 real(dp), intent(out) :: eig2nkq(2,dtset%mband*dtset%nsppol,nkpt,3,dtset%natom,3,dtset%natom*dim_eig2nkq) 248 type(efmasdeg_type),allocatable,intent(inout) :: efmasdeg(:) 249 type(efmasval_type),allocatable,intent(inout) :: efmasval(:,:) 250 type(paw_an_type),allocatable,target,intent(inout) :: paw_an(:) 251 type(paw_ij_type),allocatable,target,intent(inout) :: paw_ij(:) 252 type(pawfgrtab_type),allocatable,target,intent(inout) :: pawfgrtab(:) 253 type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw) 254 type(pawrhoij_type),allocatable,target,intent(inout) :: pawrhoij(:) 255 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 256 real(dp),pointer :: eigen1_pert(:,:,:) 257 real(dp),intent(out) :: occ_rbz_pert(:),eigen0_pert(:),eigenq_pert(:) 258 real(dp),pointer :: eigenq_fine(:,:,:) 259 integer, intent(out) :: nkpt_rbz 260 type(hdr_type),intent(out) :: hdr0,hdr_fine 261 262 !Local variables------------------------------- 263 !scalars 264 integer,parameter :: level=11,response=1,formeig1=1,master=0,fake_unit=-666 265 integer :: ask_accurate,band_index,bantot,bantot_rbz,bdeigrf,bdtot1_index,nsppol,nspinor,band2tot_index 266 integer :: bdtot_index,choice,cplex,cplex_rhoij,dim_eig2rf,formeig 267 integer :: gscase,g0term,iband,iblok,icase,icase_eq,idir,idir0,idir1,idir2,idir_eq,idir_dkdk,ierr 268 integer :: ii,ikpt,ikpt1,jband,initialized,iorder_cprj,ipert,ipert_cnt,ipert_eq,ipert_me,ireadwf0 269 integer :: iscf_mod,iscf_mod_save,isppol,istr,isym,mcg,mcgq,mcg1,mcprj,mcprjq,mband 270 integer :: mband_mem_rbz, mpert_ 271 integer :: mcgmq,mcg1mq,mpw1_mq !+/-q duplicates 272 integer :: maxidir,me,mgfftf,mkmem_rbz,mk1mem_rbz,mkqmem_rbz,mpw,mpw1,my_nkpt_rbz 273 integer :: n3xccc,nband_k,ncpgr,ndir,nkpt_eff,nkpt_max,nline_save,nmatel,npert_io,npert_me,nspden_rhoij 274 integer :: nstep_save,nsym1,ntypat,nwffile,nylmgr,nylmgr1,old_comm_atom,openexit,option,optorth,optthm,pertcase 275 integer :: qphase_rhoij,rdwr,rdwrpaw,spaceComm,smdelta,timrev_pert,timrev_kpt,to_compute_this_pert 276 integer :: useylmgr,useylmgr1,dfpt_scfcv_retcode,optn2 277 #ifdef HAVE_NETCDF 278 integer :: ncerr,ncid 279 #endif 280 real(dp) :: boxcut,dosdeltae,eberry,ecore,ecut_eff,ecutf,edocc,eei,eeig0,eew,efrhar,efrkin,efrloc 281 real(dp) :: efrnl,efrx1,efrx2,ehart,ehart01,ehart1,eii,ek,ek0,ek1,ek2,eloc0 282 real(dp) :: elpsp1,enl,enl0,enl1,end0,end1,entropy,enxc,eovl1,epaw1,evxctau0,evxctau1,evdw,exc1 283 real(dp) :: fsum,gsqcut,maxocc,nelectkq 284 real(dp) :: residm,tolwfr,tolwfr_save,toldfe_save,toldff_save,tolrff_save,tolvrs_save 285 real(dp) :: ucvol, eig1_r, eig1_i 286 real(dp) :: residm_mq !+/-q duplicates 287 logical,parameter :: paral_pert_inplace=.true.,remove_inv=.false. 288 logical :: first_entry,found_eq_gkk,has_nd,has_vxctau,t_exist,paral_atom,write_1wfk,init_rhor1 289 logical :: kramers_deg 290 character(len=fnlen) :: dscrpt,fiden1i,fiwf1i,fiwf1o,fiwfddk,fnamewff(4),gkkfilnam,fname,filnam 291 character(len=500) :: msg 292 type(crystal_t) :: crystal,ddb_crystal 293 type(dataset_type), pointer :: dtset_tmp 294 type(ebands_t) :: ebands_k,ebands_kq,gkk_ebands 295 type(ebands_t) :: ebands_kmq !+/-q duplicates 296 type(gkk_t) :: gkk2d 297 type(hdr_type) :: hdr,hdr_den,hdr_tmp 298 type(ddb_hdr_type) :: ddb_hdr, tmp_ddb_hdr 299 type(pawang_type) :: pawang1 300 type(wfk_t) :: ddk_f(4) 301 type(wvl_data) :: wvl 302 !arrays 303 integer :: eq_symop(3,3),ngfftf(18),file_index(4),rfdir(9),rf2dir(9),rf2_dir1(3),rf2_dir2(3) 304 integer,allocatable :: blkflg_save(:,:,:,:),dimcprj_srt(:),dyn(:),indkpt1(:),indkpt1_tmp(:) 305 integer,allocatable :: indsy1(:,:,:),irrzon1(:,:,:),istwfk_rbz(:),istwfk_pert(:,:,:) 306 integer,allocatable :: kg(:,:),kg1(:,:),nband_rbz(:),npwar1(:),npwarr(:),npwtot(:) 307 integer,allocatable :: kg1_mq(:,:),npwar1_mq(:),npwtot1_mq(:) !+q/-q duplicates 308 integer,allocatable :: npwtot1(:),npwar1_pert(:,:),npwarr_pert(:,:),npwtot_pert(:,:) 309 integer,allocatable :: pert_calc(:,:),pert_tmp(:,:),bz2ibz_smap(:,:) 310 integer,allocatable :: symaf1(:),symaf1_tmp(:),symrc1(:,:,:),symrl1(:,:,:),symrl1_tmp(:,:,:) 311 integer, pointer :: old_atmtab(:) 312 logical, allocatable :: distrb_flags(:,:,:) 313 real(dp) :: dielt(3,3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),tsec(2) 314 real(dp),allocatable :: buffer1(:,:,:,:,:),cg(:,:),cg1(:,:),cg1_active(:,:),cg1_3(:,:,:),cg0_pert(:,:) 315 real(dp),allocatable :: cg1_pert(:,:,:,:),cgq(:,:),gh0c1_pert(:,:,:,:) 316 real(dp),allocatable :: doccde_rbz(:),docckqde(:) 317 real(dp),allocatable :: gh1c_pert(:,:,:,:),eigen0(:),eigen0_copy(:),eigen1(:),eigen1_mean(:) 318 real(dp),allocatable :: eigenq(:),gh1c_set(:,:),gh0c1_set(:,:),kpq(:,:) 319 real(dp),allocatable :: kpq_rbz(:,:),kpt_rbz(:,:),occ_pert(:),occ_rbz(:),occkq(:),kpt_rbz_pert(:,:) 320 real(dp),allocatable :: occ_disk(:) 321 real(dp),allocatable :: vtrial_local(:,:) 322 real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phnons1(:,:,:),resid(:),rhog1(:,:) 323 real(dp),allocatable :: rhor1_save(:,:,:) 324 real(dp),allocatable :: rhor1(:,:),rho1wfg(:,:),rho1wfr(:,:),tnons1(:,:),tnons1_tmp(:,:) 325 real(dp),allocatable :: rhor1_pq(:,:),rhor1_mq(:,:),rhog1_pq(:,:),rhog1_mq(:,:) !+q/-q duplicates 326 real(dp),allocatable :: cg_mq(:,:),cg1_mq(:,:),resid_mq(:) ! 327 real(dp),allocatable :: cg1_active_mq(:,:),occk_mq(:) ! 328 real(dp),allocatable :: kmq(:,:),kmq_rbz(:,:),gh0c1_set_mq(:,:) ! 329 real(dp),allocatable :: eigen_mq(:),gh1c_set_mq(:,:),docckde_mq(:),eigen1_mq(:) ! 330 real(dp),allocatable :: vpsp1(:),work(:),wtk_folded(:),wtk_rbz(:),xccc3d1(:) 331 real(dp),allocatable :: ylm(:,:),ylm1(:,:),ylmgr(:,:,:),ylmgr1(:,:,:),zeff(:,:,:) 332 real(dp),allocatable :: phasecg(:,:),gauss(:,:) 333 real(dp),allocatable :: gkk(:,:,:,:,:) 334 logical :: has_cg1_3(3) 335 type(pawcprj_type),allocatable :: cprj(:,:),cprjq(:,:) 336 type(paw_ij_type),pointer :: paw_ij_pert(:) 337 type(paw_an_type),pointer :: paw_an_pert(:) 338 type(pawfgrtab_type),pointer :: pawfgrtab_pert(:) 339 type(pawrhoij_type),allocatable :: pawrhoij1(:) 340 type(pawrhoij_type),pointer :: pawrhoij_pert(:) 341 type(ddb_type) :: ddb 342 343 real(dp),allocatable :: doccde_tmp(:) 344 345 ! *********************************************************************** 346 347 DBG_ENTER("COLL") 348 349 call timab(141,1,tsec) 350 351 !Structured debugging if prtvol==-level 352 if(dtset%prtvol==-level)then 353 write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' dfpt_looppert : enter , debug mode ' 354 call wrtout(std_out,msg) 355 end if 356 357 dfpt_scfcv_retcode = -1 358 nsppol = dtset%nsppol; nspinor = dtset%nspinor 359 360 kramers_deg=.true. 361 if (dtset%tim1rev==0) then 362 kramers_deg=.false. 363 end if 364 365 has_vxctau = (size(vxctau) > 0) 366 evxctau0=zero; evxctau1=0 367 has_nd = ANY(ABS(dtset%nucdipmom(:,:))>tol8) 368 end0=zero; end1=zero 369 370 !Obtain dimensional translations in reciprocal space gprimd, 371 !metrics and unit cell volume, from rprimd. Also output rprimd, gprimd and ucvol 372 call mkrdim(dtset%acell_orig(1:3,1),dtset%rprim_orig(1:3,1:3,1),rprimd) 373 call metric(gmet,gprimd,std_out,rmet,rprimd,ucvol) 374 375 call crystal_init(dtset%amu_orig(:,1),crystal,dtset%spgroup,dtset%natom,dtset%npsp,& 376 & psps%ntypat,dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,& 377 & dtset%nspden==2.and.dtset%nsppol==1,remove_inv,psps%title,& 378 & symrel=dtset%symrel,tnons=dtset%tnons,symafm=dtset%symafm) 379 380 !Get FFT grid(s) sizes (be careful !) See NOTES in the comments at the beginning of respfn.F90 381 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then 382 mgfftf=pawfgr%mgfft;ngfftf(:)=pawfgr%ngfft(:) 383 ecutf=dtset%pawecutdg 384 else 385 mgfftf=dtset%mgfft;ngfftf(:)=dtset%ngfft(:) 386 ecutf=dtset%ecut 387 end if 388 ecut_eff=dtset%ecut*(dtset%dilatmx)**2 389 390 !Compute large sphere cut-off gsqcut 391 if (psps%usepaw==1) then 392 call wrtout(std_out,ch10//' FFT (fine) grid used for densities/potentials:') 393 end if 394 call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,dtset%qptn,ngfftf) 395 396 !Various initializations/allocations 397 iscf_mod=dtset%iscf 398 ntypat=psps%ntypat 399 nkpt_max=50;if (xmpi_paral==1) nkpt_max=-1 400 !TODO: this flag for paral_atom is ignored below 401 paral_atom=(dtset%natom/=my_natom) 402 cplex=2-timrev !cplex=2 ! DEBUG: impose cplex=2 403 first_entry=.true. 404 initialized=0 405 ecore=zero ; ek=zero ; ehart=zero ; enxc=zero ; eei=zero ; enl=zero ; eii=zero 406 d2bbb = zero 407 d2nl = zero 408 d2lo = zero 409 d2ovl = zero 410 clflg(:,:)=0 ! Array on calculated perturbations for eig2rf 411 if (psps%usepaw==1) then 412 ABI_MALLOC(dimcprj_srt,(dtset%natom)) 413 call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 414 end if 415 mband_mem_rbz = dtset%mband_mem 416 417 !Save values of SCF cycle parameters 418 iscf_mod_save = iscf_mod 419 nstep_save = dtset%nstep 420 nline_save = dtset%nline 421 tolwfr_save = dtset%tolwfr 422 toldfe_save = dtset%toldfe 423 toldff_save = dtset%toldff 424 tolrff_save = dtset%tolrff 425 tolvrs_save = dtset%tolvrs 426 427 !This dtset will be used in dfpt_scfcv to force non scf calculations for equivalent perturbations 428 nullify(dtset_tmp) 429 if (dtset%prepgkk/=0) then ! .and. dtset%use_nonscf_gkk==1) then !Later uncomment this - in scf case rhor1_save is used below only for testing 430 ABI_MALLOC(dtset_tmp,) 431 dtset_tmp = dtset%copy() 432 else 433 dtset_tmp => dtset 434 end if 435 436 !Generate the 1-dimensional phases 437 ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom)) 438 ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*dtset%natom)) 439 call getph(atindx,dtset%natom,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),ph1d,xred) 440 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then 441 call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred) 442 else 443 ph1df(:,:)=ph1d(:,:) 444 end if 445 446 !!Determine existence of pertubations and of pertubation symmetries 447 !!Create array with pertubations which have to be calculated 448 ! ABI_MALLOC(pert_tmp,(3*mpert)) 449 ! ipert_cnt=0 450 ! do ipert=1,mpert 451 ! do idir=1,3 452 ! if( rfpert(ipert)==1 .and. dtset%rfdir(idir) == 1 )then 453 ! if ((pertsy(idir,ipert)==1).or.& 454 !& ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.& 455 !& ((dtset%prepgkk == 1).and.(ipert <= dtset%natom)) ) then 456 ! ipert_cnt = ipert_cnt+1; 457 ! pert_tmp(ipert_cnt) = idir+(ipert-1)*3 458 ! else 459 ! write(msg, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,& 460 !& ' The perturbation idir=',idir,' ipert=',ipert,' is',ch10,& 461 !& ' symmetric of a previously calculated perturbation.',ch10,& 462 !& ' So, its SCF calculation is not needed.',ch10 463 ! call wrtout(std_out,msg,'COLL') 464 ! call wrtout(ab_out,msg,'COLL') 465 ! end if ! Test of existence of symmetry of perturbation 466 ! end if ! Test of existence of perturbation 467 ! end do 468 ! end do 469 ! ABI_MALLOC(pert_calc,(ipert_cnt)) 470 ! do icase=1,ipert_cnt 471 ! pert_calc(icase)=pert_tmp(icase) 472 ! end do 473 ! ABI_FREE(pert_tmp) 474 475 !Initialize rf2dir : 476 rf2_dir1(1:3)=dtset%rf2_pert1_dir(1:3) 477 rf2_dir2(1:3)=dtset%rf2_pert2_dir(1:3) 478 if (sum(rf2_dir1)==3.and.sum(rf2_dir2)==3.and.dtset%prepanl==1) then 479 ! Diagonal terms : 480 rf2dir(1) = rf2_dirs_from_rfpert_nl(1,1) 481 rf2dir(2) = rf2_dirs_from_rfpert_nl(2,2) 482 rf2dir(3) = rf2_dirs_from_rfpert_nl(3,3) 483 ! Upper triangular terms : 484 rf2dir(4) = rf2_dirs_from_rfpert_nl(2,3) 485 rf2dir(5) = rf2_dirs_from_rfpert_nl(1,3) 486 rf2dir(6) = rf2_dirs_from_rfpert_nl(1,2) 487 ! Lower triangular terms : 488 rf2dir(7) = rf2_dirs_from_rfpert_nl(3,2) 489 rf2dir(8) = rf2_dirs_from_rfpert_nl(3,1) 490 rf2dir(9) = rf2_dirs_from_rfpert_nl(2,1) 491 else 492 ! Diagonal terms : 493 rf2dir(1) = rf2_dir1(1)*rf2_dir2(1) 494 rf2dir(2) = rf2_dir1(2)*rf2_dir2(2) 495 rf2dir(3) = rf2_dir1(3)*rf2_dir2(3) 496 ! Upper triangular terms : 497 rf2dir(4) = rf2_dir1(2)*rf2_dir2(3) 498 rf2dir(5) = rf2_dir1(1)*rf2_dir2(3) 499 rf2dir(6) = rf2_dir1(1)*rf2_dir2(2) 500 ! Lower triangular terms : 501 rf2dir(7) = rf2_dir1(3)*rf2_dir2(2) 502 rf2dir(8) = rf2_dir1(3)*rf2_dir2(1) 503 rf2dir(9) = rf2_dir1(2)*rf2_dir2(1) 504 end if 505 506 !Determine existence of pertubations and of pertubation symmetries 507 !Create array with pertubations which have to be calculated 508 ABI_MALLOC(pert_tmp,(3,3*(dtset%natom+6)+18)) 509 ipert_cnt=0 510 do ipert=1,mpert 511 if (ipert<dtset%natom+10) then 512 maxidir = 3 513 rfdir(1:3) = dtset%rfdir(:) 514 rfdir(4:9) = 0 515 else 516 maxidir = 9 517 rfdir(1:9) = rf2dir(:) 518 end if 519 do idir=1,maxidir 520 to_compute_this_pert = 0 521 if(ipert<dtset%natom+10 .and. rfpert(ipert)==1 .and. rfdir(idir) == 1 ) then 522 if ((pertsy(idir,ipert)==1).or.& 523 & ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.& 524 & ((dtset%prepgkk == 1).and.(ipert <= dtset%natom)) ) then 525 to_compute_this_pert = 1 526 else 527 write(msg, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,& 528 ' The perturbation idir=',idir,' ipert=',ipert,' is',ch10,& 529 ' symmetric of a previously calculated perturbation.',ch10,& 530 ' So, its SCF calculation is not needed.',ch10 531 call wrtout([std_out, ab_out], msg) 532 end if ! Test of existence of symmetry of perturbation 533 else if (ipert==dtset%natom+11 .and. rfpert(ipert)==1 .and. rfdir(idir) == 1 ) then 534 to_compute_this_pert = 1 535 else if (ipert==dtset%natom+10 .and. rfpert(ipert)==1) then 536 if (dtset%rf2_dkdk==2 .or. dtset%rf2_dkdk==3) then 537 if (idir <= 3 .and. rfdir(idir) == 1) then 538 to_compute_this_pert = 1 539 else if (idir>=4.and.idir<=6) then 540 if (rfdir(idir) == 1 .or. rfdir(idir+3) == 1) to_compute_this_pert = 1 541 else if (idir>=7.and.idir<=9) then 542 if (rfdir(idir) == 1 .or. rfdir(idir-3) == 1) to_compute_this_pert = 1 543 end if 544 else 545 if (idir<=6) then 546 if (idir<=3) then 547 if (rfdir(idir) == 1) to_compute_this_pert = 1 548 else 549 if (rfdir(idir) == 1 .or. rfdir(idir+3) == 1) to_compute_this_pert = 1 550 end if 551 end if 552 end if 553 else if (ipert==dtset%natom+11 .and. rfpert(ipert)==1) then 554 if (idir <= 3 .and. rfdir(idir) == 1) then 555 to_compute_this_pert = 1 556 else if (idir>=4.and.idir<=6) then 557 if (rfdir(idir) == 1 .or. rfdir(idir+3) == 1) to_compute_this_pert = 1 558 else if (idir>=7.and.idir<=9) then 559 if (rfdir(idir) == 1 .or. rfdir(idir-3) == 1) to_compute_this_pert = 1 560 end if 561 end if 562 if (to_compute_this_pert /= 0) then 563 ipert_cnt = ipert_cnt+1; 564 pert_tmp(1,ipert_cnt) = ipert 565 pert_tmp(2,ipert_cnt) = idir 566 ! Store "pertcase" in pert_tmp(3,ipert_cnt) 567 if (ipert<dtset%natom+10) then 568 pert_tmp(3,ipert_cnt) = idir + (ipert-1)*3 569 else 570 pert_tmp(3,ipert_cnt) = idir + (ipert-dtset%natom-10)*9 + (dtset%natom+6)*3 571 end if 572 end if 573 end do ! idir 574 end do !ipert 575 ! do ipert=1,mpert 576 ! if (ipert<dtset%natom+10) then 577 ! maxidir = 3 578 ! rfdir(1:3) = dtset%rfdir(:) 579 ! rfdir(4:9) = 0 580 ! else 581 ! maxidir = 9 582 ! rfdir(1:9) = rf2dir(:) 583 ! end if 584 ! do idir=1,maxidir 585 ! if( rfpert(ipert)==1 .and. rfdir(idir) == 1 )then 586 ! to_compute_this_pert = 0 587 ! if (ipert>=dtset%natom+10) then 588 ! to_compute_this_pert = 1 589 ! else if ((pertsy(idir,ipert)==1).or.& 590 !& ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.& 591 !& ((dtset%prepgkk == 1).and.(ipert <= dtset%natom)) ) then 592 ! to_compute_this_pert = 1 593 ! end if 594 ! if (to_compute_this_pert /= 0) then 595 ! ipert_cnt = ipert_cnt+1; 596 ! pert_tmp(1,ipert_cnt) = ipert 597 ! pert_tmp(2,ipert_cnt) = idir 598 !! Store "pertcase" in pert_tmp(3,ipert_cnt) 599 ! if (ipert<dtset%natom+10) then 600 ! pert_tmp(3,ipert_cnt) = idir + (ipert-1)*3 601 ! else 602 ! pert_tmp(3,ipert_cnt) = idir + (ipert-dtset%natom-10)*9 + (dtset%natom+6)*3 603 ! end if 604 ! else 605 ! write(msg, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,& 606 !& ' The perturbation idir=',idir,' ipert=',ipert,' is',ch10,& 607 !& ' symmetric of a previously calculated perturbation.',ch10,& 608 !& ' So, its SCF calculation is not needed.',ch10 609 ! call wrtout(std_out,msg,'COLL') 610 ! call wrtout(ab_out,msg,'COLL') 611 ! end if ! Test of existence of symmetry of perturbation 612 ! end if ! Test of existence of perturbation 613 ! end do 614 ! end do 615 ABI_MALLOC(pert_calc,(3,ipert_cnt)) 616 do icase=1,ipert_cnt 617 pert_calc(:,icase)=pert_tmp(:,icase) 618 end do 619 ABI_FREE(pert_tmp) 620 621 if (dtset%prepgkk/=0) then ! .and. dtset%use_nonscf_gkk==1) then !Later uncomment this - in scf case rhor1_save is used below only for testing 622 ABI_MALLOC(rhor1_save,(cplex*nfftf,nspden,ipert_cnt)) 623 rhor1_save=zero 624 ABI_MALLOC(blkflg_save,(3,mpert,3,mpert)) 625 end if 626 627 ! Initialize quantities for netcdf print 628 ABI_MALLOC(eigen0_copy,(dtset%mband*nkpt*dtset%nsppol)) 629 eigen0_copy(:)=zero 630 631 ! SP : Retreval of the DDB information and computing of effective charge and 632 ! dielectric tensor 633 ABI_MALLOC(zeff,(3,3,dtset%natom)) 634 if (dtset%getddb .ne. 0 .or. dtset%irdddb .ne. 0 ) then 635 filnam = dtfil%filddbsin 636 call ddb%from_file(filnam, tmp_ddb_hdr, ddb_crystal, mpi_enreg%comm_world) 637 call tmp_ddb_hdr%free() 638 ! Get Dielectric Tensor and Effective Charges 639 ! (initialized to one_3D and zero if the derivatives are not available in the DDB file) 640 iblok = ddb%get_dielt_zeff(ddb_crystal,1,0,0,dielt,zeff) 641 call ddb_crystal%free() 642 call ddb%free() 643 end if 644 645 !%%%% Parallelization over perturbations %%%%% 646 !*Define file output/log file names 647 npert_io=ipert_cnt;if (dtset%nppert<=1) npert_io=0 648 call localfilnam(mpi_enreg%comm_pert,mpi_enreg%comm_cell_pert,mpi_enreg%comm_world,dtfil%filnam_ds,'_PRT',npert_io) 649 !Compute the number of perturbation done by the current cpu 650 if(mpi_enreg%paral_pert==1) then 651 npert_me = 0 ; ipert_me = 0 652 do icase=1,ipert_cnt 653 if (mpi_enreg%distrb_pert(icase)==mpi_enreg%me_pert) npert_me=npert_me +1 654 end do 655 end if 656 657 !*Redefine communicators 658 call set_pert_comm(mpi_enreg,dtset%nppert) 659 660 !*Redistribute PAW on-site data 661 nullify(old_atmtab,pawfgrtab_pert,pawrhoij_pert,paw_an_pert,paw_ij_pert) 662 if (paral_pert_inplace) then 663 call set_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,& 664 & paw_an,paw_ij,pawfgrtab,pawrhoij) 665 pawfgrtab_pert=>pawfgrtab ; pawrhoij_pert=>pawrhoij 666 paw_an_pert =>paw_an ; paw_ij_pert =>paw_ij 667 668 else 669 call set_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,& 670 & paw_an,paw_ij,pawfgrtab,pawrhoij,& 671 & paw_an_out=paw_an_pert,paw_ij_out=paw_ij_pert,& 672 & pawfgrtab_out=pawfgrtab_pert,pawrhoij_out=pawrhoij_pert) 673 674 end if 675 676 ! We can handle the time limit in dfpt_scfcv in a robust manner only if we have one perturbation. 677 if (ipert_cnt > 1 .or. mpi_enreg%paral_pert == 1) call disable_timelimit() 678 679 !Loop on perturbations 680 !========================================================================== 681 do icase=1,ipert_cnt 682 683 ! %%%% Parallelization over perturbations %%%%% 684 ! Select the perturbations treated by curent processor 685 if(mpi_enreg%paral_pert==1) then 686 if (mpi_enreg%distrb_pert(icase)/=mpi_enreg%me_pert) cycle 687 end if 688 689 ! Redefine output/log files 690 call localwrfile(mpi_enreg%comm_cell,icase,npert_io,mpi_enreg%paral_pert,0) 691 692 ! Set precision for FFT libs. 693 ii = fftcore_set_mixprec(dtset%mixprec) 694 695 !! Retrieve type and direction of the perturbation 696 ! if (pert_calc(icase) <= dtset%natom*3) then 697 ! idir = mod(pert_calc(icase),3) 698 ! if (idir==0) idir=3 699 ! ipert=( (pert_calc(icase)-idir) / 3 + 1) 700 ! else if (pert_calc(icase) <= dtset%natom*3+4) then 701 ! ipert = dtset%natom + ((pert_calc(icase) - 3*dtset%natom - 1) / 3) + 1 702 ! idir = mod(pert_calc(icase),3) 703 ! if (idir==0) idir=3 704 ! else 705 ! ipert = dtset%natom + ((pert_calc(icase) - 3*(dtset%natom+4) - 1) / 9) + 1 706 ! end if 707 ! pertcase=idir+(ipert-1)*3 708 709 ! Retrieve type and direction of the perturbation 710 ipert = pert_calc(1,icase) 711 idir = pert_calc(2,icase) 712 istr=idir 713 pertcase = pert_calc(3,icase) 714 715 ! Init MPI communicator 716 spaceComm=mpi_enreg%comm_cell 717 me=mpi_enreg%me_cell 718 719 ! ===== Describe the perturbation in output/log file 720 721 write(msg, '(a,80a,a,a,3f10.6)' ) ch10,('-',ii=1,80),ch10,& 722 ' Perturbation wavevector (in red.coord.) ',dtset%qptn(:) 723 call wrtout([std_out, ab_out],msg) 724 if(ipert>=1 .and. ipert<=dtset%natom)then 725 write(msg, '(a,i4,a,i4)' )' Perturbation : displacement of atom',ipert,' along direction',idir 726 call wrtout([std_out, ab_out], msg) 727 if(iscf_mod == -3)then 728 write(msg, '(a,a,a,a,a,a,a,a)' )ch10,& 729 ' dfpt_looppert : COMMENT -',ch10,& 730 ' The first-order density is imposed to be zero (iscf=-3).',ch10,& 731 ' Although this is strange in the case of phonons,',ch10,& 732 ' you are allowed to do so.' 733 call wrtout([std_out, ab_out], msg) 734 end if 735 else if(ipert==dtset%natom+1)then 736 write(msg,'(a,i4)')' Perturbation : derivative vs k along direction',idir 737 call wrtout([std_out, ab_out], msg) 738 if( iscf_mod /= -3 )then 739 write(msg, '(4a)' )ch10,& 740 ' dfpt_looppert : COMMENT -',ch10,& 741 ' In a d/dk calculation, iscf is set to -3 automatically.' 742 call wrtout([std_out, ab_out], msg) 743 iscf_mod=-3 744 end if 745 if( abs(dtset%dfpt_sciss) > 1.0d-8 )then 746 write(msg, '(a,a,a,a,f14.8,a,a)' )ch10,& 747 ' dfpt_looppert : WARNING -',ch10,& 748 ' Value of dfpt_sciss=',dtset%dfpt_sciss,ch10,& 749 ' Scissor with d/dk calculation : you are using a "naive" approach !' 750 call wrtout([std_out, ab_out], msg) 751 end if 752 else if(ipert==dtset%natom+2)then 753 write(msg, '(a,i4)' )' Perturbation : homogeneous electric field along direction',idir 754 call wrtout([std_out, ab_out], msg) 755 if( iscf_mod == -3 )then 756 write(msg, '(a,a,a,a,a,a)' )ch10,& 757 ' dfpt_looppert : COMMENT -',ch10,& 758 ' The first-order density is imposed to be zero (iscf=-3).',ch10,& 759 ' This corresponds to a calculation without local fields.' 760 call wrtout([std_out, ab_out], msg) 761 end if 762 else if(ipert==dtset%natom+10.or.ipert==dtset%natom+11)then 763 call rf2_getidirs(idir,idir1,idir2) 764 if(ipert==dtset%natom+10)then 765 write(msg,'(2(a,i1))') ' Perturbation : 2nd derivative wrt k, idir1 = ',idir1,& 766 ' idir2 = ',idir2 767 else 768 write(msg,'(2(a,i1),a)') ' Perturbation : 2nd derivative wrt k (idir1 =',idir1,& 769 ') and Efield (idir2 =',idir2,')' 770 end if 771 call wrtout([std_out, ab_out], msg) 772 if( iscf_mod /= -3 )then 773 write(msg, '(4a)' )ch10,& 774 ' dfpt_looppert : COMMENT -',ch10,& 775 ' In this case, iscf is set to -3 automatically.' 776 call wrtout([std_out, ab_out], msg) 777 iscf_mod=-3 778 end if 779 if( abs(dtset%dfpt_sciss) > 1.0d-8 )then 780 write(msg, '(a,a,a,a,f14.8,a,a)' )ch10,& 781 ' dfpt_looppert : WARNING -',ch10,& 782 ' Value of dfpt_sciss=',dtset%dfpt_sciss,ch10,& 783 ' Scissor with d/dk calculation : you are using a "naive" approach !' 784 call wrtout([std_out, ab_out], msg) 785 end if 786 ABI_MALLOC(occ_pert,(dtset%mband*nkpt*dtset%nsppol)) 787 occ_pert(:) = occ(:) - occ(1) 788 maxocc = maxval(abs(occ_pert)) 789 if (maxocc>1.0d-6.and.abs(maxocc-occ(1))>1.0d-6) then ! True if non-zero occupation numbers are not equal 790 write(msg, '(3a)' ) ' ipert=natom+10 or 11 does not work for a metallic system.',ch10,& 791 ' This perturbation will not be computed.' 792 ABI_WARNING(msg) 793 ABI_FREE(occ_pert) 794 cycle 795 end if 796 ABI_FREE(occ_pert) 797 else if(ipert>dtset%natom+11 .or. ipert<=0 )then 798 write(msg, '(a,i0,3a)' ) & 799 'ipert= ',ipert,' is outside the [1,dtset%natom+11] interval.',ch10,& 800 'This perturbation is not (yet) allowed.' 801 ABI_BUG(msg) 802 end if 803 ! Initialize the diverse parts of energy : 804 eew=zero ; evdw=zero ; efrloc=zero ; efrnl=zero ; efrx1=zero ; efrx2=zero 805 efrhar=zero ; efrkin=zero 806 if(ipert<=dtset%natom)then 807 eew=dyew(1,idir,ipert,idir,ipert) 808 if (usevdw==1) evdw=dyvdw(1,idir,ipert,idir,ipert) 809 efrloc=dyfrlo(idir,idir,ipert) 810 if (dyfr_nondiag==0) efrnl=dyfrnl(1,idir,idir,ipert,1) 811 if (dyfr_nondiag/=0) efrnl=dyfrnl(1,idir,idir,ipert,ipert) 812 efrx1=dyfrx1(1,idir,ipert,idir,ipert) 813 efrx2=dyfrx2(idir,idir,ipert) 814 else if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then 815 ! istr = 1,2,...,6 and indicates the cartesian strain component 816 if(ipert==dtset%natom+4) istr=idir+3 817 eii=eltcore(istr,istr) 818 eew=elteew(istr,istr) 819 if (usevdw==1) evdw=eltvdw(istr,istr) 820 efrhar=eltfrhar(istr,istr) 821 efrkin=eltfrkin(istr,istr) 822 efrloc=eltfrloc(istr,istr) 823 efrnl=eltfrnl(istr,istr) 824 efrx1=eltfrxc(istr,istr) 825 end if 826 827 ! Determine the subset of symmetry operations (nsym1 operations) 828 ! that leaves the perturbation invariant, and initialize corresponding arrays 829 ! symaf1, symrl1, tnons1 (and pawang1%zarot, if PAW).. 830 ABI_MALLOC(symaf1_tmp,(nsym)) 831 ABI_MALLOC(symrl1_tmp,(3,3,nsym)) 832 ABI_MALLOC(tnons1_tmp,(3,nsym)) 833 834 if (dtset%prepanl/=1.and.& 835 & dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.& 836 & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17) then 837 call littlegroup_pert(gprimd,idir,indsym,ab_out,ipert,dtset%natom,nsym,nsym1,2,& 838 & dtset%symafm,symaf1_tmp,symq,symrec,dtset%symrel,symrl1_tmp,0,dtset%tnons,tnons1_tmp) 839 else 840 nsym1 = 1 841 symaf1_tmp(1) = 1 842 symrl1_tmp(:,:,1) = dtset%symrel(:,:,1) 843 tnons1_tmp(:,1) = 0_dp 844 end if 845 ABI_MALLOC(indsy1,(4,nsym1,dtset%natom)) 846 ABI_MALLOC(symrc1,(3,3,nsym1)) 847 ABI_MALLOC(symaf1,(nsym1)) 848 ABI_MALLOC(symrl1,(3,3,nsym1)) 849 ABI_MALLOC(tnons1,(3,nsym1)) 850 symaf1(1:nsym1)=symaf1_tmp(1:nsym1) 851 symrl1(:,:,1:nsym1)=symrl1_tmp(:,:,1:nsym1) 852 tnons1(:,1:nsym1)=tnons1_tmp(:,1:nsym1) 853 ABI_FREE(symaf1_tmp) 854 ABI_FREE(symrl1_tmp) 855 ABI_FREE(tnons1_tmp) 856 857 ! Set up corresponding symmetry data 858 ABI_MALLOC(irrzon1,(dtset%nfft**(1-1/nsym1),2,(nspden/dtset%nsppol)-3*(nspden/4))) 859 ABI_MALLOC(phnons1,(2,dtset%nfft**(1-1/nsym1),(nspden/dtset%nsppol)-3*(nspden/4))) 860 call setsym(indsy1,irrzon1,1,dtset%natom,dtset%nfft,dtset%ngfft,nspden,dtset%nsppol,& 861 & nsym1,phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred) 862 if (psps%usepaw==1) then 863 ! Allocate/initialize only zarot in pawang1 datastructure 864 call pawang_init(pawang1,0,0,pawang%l_max-1,0,0,nsym1,0,0,0,0) 865 call setsym_ylm(gprimd,pawang1%l_max-1,pawang1%nsym,0,rprimd,symrc1,pawang1%zarot) 866 end if 867 868 ! Initialize k+q array 869 ABI_MALLOC(kpq,(3,nkpt)) 870 if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then 871 kpq(:,1:nkpt)=dtset%kptns(:,1:nkpt) ! Do not modify, needed for gfortran 872 else 873 do ikpt=1,nkpt 874 kpq(:,ikpt)=dtset%qptn(:)+dtset%kptns(:,ikpt) 875 end do 876 end if 877 ! In case wf1 at +q and -q are not related by time inversion symmetry, compute k-q as well for initializations 878 if (.not.kramers_deg) then 879 ABI_MALLOC(kmq,(3,nkpt)) 880 if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then 881 kmq(:,1:nkpt)=dtset%kptns(:,1:nkpt) ! Do not modify, needed for gfortran 882 else 883 do ikpt=1,nkpt 884 kmq(:,ikpt)=-dtset%qptn(:)+dtset%kptns(:,ikpt) ! kmq <= k-q 885 end do 886 end if 887 end if 888 889 ! Determine the subset of k-points needed in the "reduced Brillouin zone" and initialize other quantities 890 ABI_MALLOC(indkpt1_tmp,(nkpt)) 891 ABI_MALLOC(wtk_folded,(nkpt)) 892 ABI_MALLOC(bz2ibz_smap, (6, nkpt)) 893 indkpt1_tmp(:)=0 ; optthm=0 894 timrev_pert=timrev 895 if(dtset%ieig2rf>0) then 896 timrev_pert=0 897 call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,& 898 & 1,symrc1,timrev_pert,dtset%wtk,wtk_folded, bz2ibz_smap, xmpi_comm_self) 899 else 900 ! For the time being, the time reversal symmetry is not used 901 ! for ddk, elfd, mgfd perturbations. 902 timrev_pert=timrev 903 if(ipert==dtset%natom+1.or.ipert==dtset%natom+2.or.& 904 & ipert==dtset%natom+10.or.ipert==dtset%natom+11.or. & 905 & dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or. & 906 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17.or. & 907 & ipert==dtset%natom+5.or.dtset%prtfull1wf==1) timrev_pert=0 908 timrev_kpt = timrev_pert 909 910 !MR: Modified to agree with longwave driver 911 if(dtset%prepalw/=0) then 912 if (dtset%kptopt==2) timrev_pert=1 913 if (dtset%kptopt==3) timrev_pert=0 914 timrev_kpt = timrev_pert 915 !MR tmp: this has to be removed if perturbation-dependent spatial symmetries are 916 !implemented in the quadrupole and flexoelectrics routines 917 nsym1=1 918 919 if (dtset%rfstrs/=0.and.dtset%rfstrs_ref==0) then 920 write(msg,'(9a)')& 921 'If the outputs of this strain response function calculation are to be',ch10,& 922 'subsequently used as inputs of a longwave calculation the same energy',ch10,& 923 'reference has to be used in both cases. Otherwise wrong spatial-dispersion',ch10,& 924 'coefficients will be obtained',ch10,& 925 'Action: Put rfstrs_ref=1 ' 926 ABI_WARNING(msg) 927 end if 928 end if 929 930 ! The time reversal symmetry is not used for the BZ sampling when kptopt=3 or 4 931 if (dtset%kptopt==3.or.dtset%kptopt==4) timrev_kpt = 0 932 call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,& 933 nsym1,symrc1,timrev_kpt,dtset%wtk,wtk_folded, bz2ibz_smap, xmpi_comm_self) 934 end if 935 936 ABI_MALLOC(doccde_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol)) 937 ABI_MALLOC(indkpt1,(nkpt_rbz)) 938 ABI_MALLOC(istwfk_rbz,(nkpt_rbz)) 939 ABI_MALLOC(kpq_rbz,(3,nkpt_rbz)) 940 if (.not.kramers_deg) then 941 ABI_MALLOC(kmq_rbz,(3,nkpt_rbz)) 942 end if 943 ABI_MALLOC(kpt_rbz,(3,nkpt_rbz)) 944 ABI_MALLOC(nband_rbz,(nkpt_rbz*dtset%nsppol)) 945 ABI_MALLOC(occ_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol)) 946 ABI_MALLOC(occ_disk,(dtset%mband*nkpt_rbz*dtset%nsppol)) 947 ABI_MALLOC(wtk_rbz,(nkpt_rbz)) 948 indkpt1(:)=indkpt1_tmp(1:nkpt_rbz) 949 do ikpt=1,nkpt_rbz 950 istwfk_rbz(ikpt)=dtset%istwfk(indkpt1(ikpt)) 951 kpq_rbz(:,ikpt)=kpq(:,indkpt1(ikpt)) 952 kpt_rbz(:,ikpt)=dtset%kptns(:,indkpt1(ikpt)) 953 wtk_rbz(ikpt)=wtk_folded(indkpt1(ikpt)) 954 end do 955 if (.not.kramers_deg) then 956 do ikpt=1,nkpt_rbz 957 kmq_rbz(:,ikpt)=kmq(:,indkpt1(ikpt)) 958 end do 959 end if 960 ABI_FREE(indkpt1_tmp) 961 ABI_FREE(wtk_folded) 962 963 ! Transfer occ to occ_rbz and doccde to doccde_rbz : 964 ! this is a more delicate issue 965 ! NOTE : this takes into account that indkpt1 is ordered 966 ! MG: What about using occ(band,kpt,spin) ??? 967 ! GA: Would be better indeed, but I think indices neededed to be consistent 968 ! with second derivative of eigenvalues, which has spin index wrapped up 969 ! to avoid arrays with rank larger than 7. 970 bdtot_index=0;bdtot1_index=0 971 do isppol=1,dtset%nsppol 972 ikpt1=1 973 do ikpt=1,nkpt 974 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 975 ! Must test against ikpt1/=nkpt_rbz+1, before evaluate indkpt1(ikpt1) 976 if(ikpt1/=nkpt_rbz+1)then 977 if(ikpt==indkpt1(ikpt1))then 978 nband_rbz(ikpt1+(isppol-1)*nkpt_rbz)=nband_k 979 occ_rbz(1+bdtot1_index:nband_k+bdtot1_index) = occ(1+bdtot_index:nband_k+bdtot_index) 980 doccde_rbz(1+bdtot1_index:nband_k+bdtot1_index) = doccde(1+bdtot_index:nband_k+bdtot_index) 981 ikpt1=ikpt1+1 982 bdtot1_index=bdtot1_index+nband_k 983 end if 984 end if 985 bdtot_index=bdtot_index+nband_k 986 end do 987 end do 988 989 ! Compute maximum number of planewaves at k 990 call timab(142,1,tsec) 991 call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpt_rbz,mpi_enreg,mpw,nkpt_rbz) 992 call timab(142,2,tsec) 993 994 ! Allocate some k-dependent arrays at k 995 ABI_MALLOC(npwarr,(nkpt_rbz)) 996 ABI_MALLOC(npwtot,(nkpt_rbz)) 997 998 ! Determine distribution of k-points/bands over MPI processes 999 if (allocated(mpi_enreg%my_kpttab)) then 1000 ABI_FREE(mpi_enreg%my_kpttab) 1001 end if 1002 ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt_rbz)) 1003 if(xmpi_paral==1) then 1004 ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt_rbz,dtset%mband,dtset%nsppol)) 1005 call distrb2(dtset%mband,mband_mem_rbz,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,dtset%nsppol,mpi_enreg) 1006 else 1007 mpi_enreg%my_kpttab(:)=(/(ii,ii=1,nkpt_rbz)/) 1008 end if 1009 my_nkpt_rbz=maxval(mpi_enreg%my_kpttab) 1010 mkmem_rbz =my_nkpt_rbz ; mkqmem_rbz=my_nkpt_rbz ; mk1mem_rbz=my_nkpt_rbz 1011 ABI_UNUSED((/mkmem,mk1mem,mkqmem/)) 1012 1013 call initmpi_band(mkmem_rbz,mpi_enreg,nband_rbz,nkpt_rbz,dtset%nsppol) 1014 1015 ! given number of reduced kpt, store distribution of bands across procs 1016 ABI_MALLOC(distrb_flags,(nkpt_rbz,dtset%mband,dtset%nsppol)) 1017 distrb_flags = (mpi_enreg%proc_distrb == mpi_enreg%me_kpt) 1018 1019 ! Set up the basis sphere of planewaves at k 1020 ABI_MALLOC(kg,(3,mpw*mkmem_rbz)) 1021 call timab(143,1,tsec) 1022 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg,& 1023 & kpt_rbz,mkmem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,& 1024 & mpw,npwarr,npwtot,dtset%nsppol) 1025 call timab(143,2,tsec) 1026 1027 ! Set up the spherical harmonics (Ylm) at k 1028 useylmgr=0; option=0 ; nylmgr=0 1029 if (psps%useylm==1.and. & 1030 & (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or. & 1031 & (psps%usepaw==1.and.ipert==dtset%natom+2))) then 1032 useylmgr=1; option=1 ; nylmgr=3 1033 else if (psps%useylm==1.and.(ipert==dtset%natom+10.or.ipert==dtset%natom+11)) then 1034 useylmgr=1; option=2 ; nylmgr=9 1035 end if 1036 ABI_MALLOC(ylm,(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm)) 1037 ABI_MALLOC(ylmgr,(mpw*mkmem_rbz,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 1038 if (psps%useylm==1) then 1039 call initylmg(gprimd,kg,kpt_rbz,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,nband_rbz,nkpt_rbz,& 1040 & npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr) 1041 end if 1042 1043 ! Set up occupations for this perturbation 1044 if (dtset%ieig2rf>0) then 1045 if (.not.allocated(istwfk_pert)) then 1046 ABI_MALLOC(istwfk_pert,(nkpt,3,mpert)) 1047 ABI_MALLOC(occ_pert,(dtset%mband*nkpt*dtset%nsppol)) 1048 istwfk_pert(:,:,:)=0 ; occ_pert(:)=zero 1049 end if 1050 istwfk_pert(:,idir,ipert)=istwfk_rbz(:) 1051 occ_pert(:)=occ_rbz(:) 1052 end if 1053 if (dtset%efmas>0) then 1054 if (.not.allocated(istwfk_pert)) then 1055 ABI_MALLOC(istwfk_pert,(nkpt,3,mpert)) 1056 istwfk_pert(:,:,:)=0 1057 end if 1058 istwfk_pert(:,idir,ipert)=istwfk_rbz(:) 1059 end if 1060 1061 ! Print a separator in output file 1062 write(msg,'(3a)')ch10,'--------------------------------------------------------------------------------',ch10 1063 call wrtout(ab_out,msg) 1064 1065 ! Initialize band structure datatype at k 1066 ! TODO: occ_rbz and doccde_rbz are overwritten here, though they could be set from the input file, above 1067 ! need to check if this is correct (or indifferent) in all cases 1068 bantot_rbz=sum(nband_rbz(1:nkpt_rbz*dtset%nsppol)) 1069 ABI_MALLOC(eigen0,(bantot_rbz)) 1070 eigen0(:)=zero 1071 call ebands_init(bantot_rbz,ebands_k,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,& 1072 doccde_rbz,eigen0,istwfk_rbz,kpt_rbz,& 1073 nband_rbz,nkpt_rbz,npwarr,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,& 1074 dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 1075 dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 1076 ABI_FREE(eigen0) 1077 1078 ! Initialize header, update it with evolving variables 1079 gscase=0 ! A GS WF file is read 1080 call hdr_init(ebands_k,codvsn,dtset,hdr0,pawtab,gscase,psps,wvl%descr,& 1081 comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1082 1083 call hdr0%update(bantot_rbz,etotal,fermie,fermie,& ! CP: duplicated fermie to fit new def of hdr_update 1084 residm,rprimd,occ_rbz,pawrhoij_pert,xred,dtset%amu_orig(:,1),& 1085 comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1086 1087 ! Initialize GS wavefunctions at k 1088 ireadwf0=1; formeig=0 ; ask_accurate=1 ; optorth=0 1089 mcg=mpw*dtset%nspinor*mband_mem_rbz*mkmem_rbz*dtset%nsppol 1090 1091 call wrtout(std_out, sjoin(" Memory required for psi0_k, psi0_kq, psi1_kq: ", & 1092 ftoa(3 * two * mcg * dp * b2Mb, fmt="f8.1"), "[Mb] <<< MEM")) 1093 1094 if (one*mpw*dtset%nspinor*mband_mem_rbz*mkmem_rbz*dtset%nsppol > huge(1)) then 1095 write (msg,'(4a, 5(a,i0), 2a)')& 1096 "Default integer is not wide enough to store the size of the GS wavefunction array (WF0, mcg).",ch10,& 1097 "Action: increase the number of processors. Consider also OpenMP threads.",ch10,& 1098 "nspinor: ",dtset%nspinor, "mpw: ",mpw, "mband_mem_rbz: ",mband_mem_rbz, "mkmem_rbz: ",& 1099 mkmem_rbz, "nsppol: ",dtset%nsppol,ch10,& 1100 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode' 1101 ABI_ERROR(msg) 1102 end if 1103 ABI_MALLOC_OR_DIE(cg,(2,mcg), ierr) 1104 1105 ABI_MALLOC(eigen0,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1106 call timab(144,1,tsec) 1107 1108 ! Initialize the wave function type and read GS WFK 1109 call wfk_read_my_kptbands(dtfil%fnamewffk, distrb_flags, spacecomm, dtset%ecut*(dtset%dilatmx)**2,& 1110 & formeig, istwfk_rbz, kpt_rbz, mcg, dtset%mband, mband_mem_rbz, mkmem_rbz, mpw,& 1111 & dtset%natom, nkpt_rbz, npwarr, dtset%nspinor, dtset%nsppol, dtset%usepaw,& 1112 & cg, eigen=eigen0, occ=occ_disk) 1113 1114 call timab(144,2,tsec) 1115 1116 1117 ! Update energies GS energies at k 1118 call put_eneocc_vect(ebands_k, "eig", eigen0) 1119 1120 ! PAW: compute on-site projections of GS wavefunctions (cprj) (and derivatives) at k 1121 ncpgr=0 1122 ABI_MALLOC(cprj,(0,0)) 1123 if (psps%usepaw==1) then 1124 ncpgr=3 ! Valid for ipert<=natom (phonons), ipert=natom+2 (elec. field) 1125 ! or for ipert==natom+10,11 1126 if (ipert==dtset%natom+1) then 1127 if (dtset%orbmag.NE.0) then 1128 ncpgr=3 1129 else 1130 ncpgr=1 1131 end if 1132 end if 1133 if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) ncpgr=1 1134 if (usecprj==1) then 1135 ! distribute cprj by band as well 1136 ! NB: currently nsppol=2 is distributed in data (0s saved for other spin) 1137 ! but not in memory: all procs have nsppol 2 below 1138 mcprj=dtset%nspinor*mband_mem_rbz*mkmem_rbz*dtset%nsppol 1139 !mcprj=dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol 1140 ABI_FREE(cprj) 1141 ABI_MALLOC(cprj,(dtset%natom,mcprj)) 1142 call pawcprj_alloc(cprj,ncpgr,dimcprj_srt) 1143 if (ipert<=dtset%natom) then 1144 choice=2; iorder_cprj=0; idir0=0 1145 else if (ipert==dtset%natom+1) then 1146 if (dtset%orbmag.NE.0) then 1147 choice=5; iorder_cprj=0; idir0=0 1148 else 1149 choice=5; iorder_cprj=0; idir0=idir 1150 end if 1151 else if (ipert==dtset%natom+2) then 1152 choice=5; iorder_cprj=0; idir0=0 1153 else if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then 1154 choice=3; iorder_cprj=0; idir0=istr 1155 else if (ipert==dtset%natom+10.or.ipert==dtset%natom+11) then 1156 choice=5; iorder_cprj=0; idir0=0 ! Compute all first derivatives 1157 else 1158 choice=1; iorder_cprj=0; idir0=0 1159 end if 1160 call ctocprj(atindx,cg,choice,cprj,gmet,gprimd,-1,idir0,iorder_cprj,istwfk_rbz,& 1161 & kg,kpt_rbz,mcg,mcprj,dtset%mgfft,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,& 1162 & dtset%natom,nattyp,nband_rbz,dtset%natom,dtset%ngfft,nkpt_rbz,dtset%nloalg,& 1163 & npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,psps,& 1164 & rmet,dtset%typat,ucvol,dtfil%unpaw,xred,ylm,ylmgr) 1165 end if 1166 end if 1167 1168 ! Compute maximum number of planewaves at k+q 1169 ! Will be useful for both GS wfs at k+q and RF wavefunctions 1170 call timab(143,1,tsec) 1171 call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpq_rbz,mpi_enreg,mpw1,nkpt_rbz) 1172 if (.not.kramers_deg) then 1173 call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kmq_rbz,mpi_enreg,mpw1_mq,nkpt_rbz) 1174 !number of plane waves at k+q and k-q should be in principle the same to reconstruct rhor1_pq (?) 1175 !mpw1=max(mpw1,mpw1_tmp) 1176 else 1177 mpw1_mq=0 1178 end if 1179 call timab(143,2,tsec) 1180 1181 ! Allocate some arrays at k+q 1182 ABI_MALLOC(kg1,(3,mpw1*mk1mem_rbz)) 1183 ABI_MALLOC(npwar1,(nkpt_rbz)) 1184 ABI_MALLOC(npwtot1,(nkpt_rbz)) 1185 ! In case Kramers degeneracy is broken, do the same for k-q 1186 if (.not.kramers_deg) then 1187 ABI_MALLOC(kg1_mq,(3,mpw1_mq*mk1mem_rbz)) 1188 ABI_MALLOC(npwar1_mq,(nkpt_rbz)) 1189 ABI_MALLOC(npwtot1_mq,(nkpt_rbz)) 1190 end if 1191 1192 ! Set up the basis sphere of planewaves at k+q 1193 ! Will be useful for both GS wfs at k+q and RF wavefunctions 1194 call timab(142,1,tsec) 1195 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg1,& 1196 & kpq_rbz,mk1mem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw1,& 1197 & npwar1,npwtot1,dtset%nsppol) 1198 if (.not.kramers_deg) then 1199 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg1_mq,& 1200 & kmq_rbz,mk1mem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw1_mq,& 1201 & npwar1_mq,npwtot1_mq,dtset%nsppol) 1202 end if 1203 call timab(142,2,tsec) 1204 1205 ! Set up the spherical harmonics (Ylm) at k+q 1206 useylmgr1=0; option=0 ; nylmgr1=0 1207 if (psps%useylm==1.and. & 1208 & (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or. & 1209 & (psps%usepaw==1.and.ipert==dtset%natom+2))) then 1210 useylmgr1=1; option=1; nylmgr1=3 1211 else if (psps%useylm==1.and.(ipert==dtset%natom+10.or.ipert==dtset%natom+11)) then 1212 useylmgr1=1; option=2; nylmgr1=9 1213 end if 1214 ABI_MALLOC(ylm1,(mpw1*mk1mem_rbz,psps%mpsang*psps%mpsang*psps%useylm)) 1215 ABI_MALLOC(ylmgr1,(mpw1*mk1mem_rbz,nylmgr1,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 1216 if (psps%useylm==1) then 1217 call initylmg(gprimd,kg1,kpq_rbz,mk1mem_rbz,mpi_enreg,psps%mpsang,mpw1,nband_rbz,nkpt_rbz,& 1218 & npwar1,dtset%nsppol,option,rprimd,ylm1,ylmgr1) 1219 end if 1220 ! SPr: do the same for k-q if kramers_deg=.false. 1221 1222 ! Print a separator in output file 1223 write(msg, '(a,a)' )'--------------------------------------------------------------------------------',ch10 1224 call wrtout(ab_out,msg) 1225 1226 ! Initialize band structure datatype at k+q 1227 ABI_MALLOC(eigenq,(bantot_rbz)) 1228 ABI_MALLOC(doccde_tmp,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1229 eigenq(:)=zero 1230 call ebands_init(bantot_rbz,ebands_kq,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,& 1231 & doccde_tmp,eigenq,istwfk_rbz,kpq_rbz,& 1232 & nband_rbz,nkpt_rbz,npwar1,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,& 1233 & dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 1234 & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 1235 if (.not.kramers_deg) then 1236 eigenq(:)=zero 1237 call ebands_init(bantot_rbz,ebands_kmq,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,& 1238 & doccde_tmp,eigenq,istwfk_rbz,kmq_rbz,& 1239 & nband_rbz,nkpt_rbz,npwar1_mq,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,& 1240 & dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 1241 & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 1242 end if 1243 ABI_FREE(eigenq) 1244 ABI_FREE(doccde_tmp) 1245 1246 ! Initialize header 1247 call hdr_init(ebands_kq,codvsn,dtset,hdr,pawtab,pertcase,psps,wvl%descr, & 1248 comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab ) 1249 if (.not.kramers_deg) then 1250 call hdr_init(ebands_kmq,codvsn,dtset,hdr,pawtab,pertcase,psps,wvl%descr, & 1251 comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab ) 1252 end if 1253 1254 ! Initialize wavefunctions at k+q 1255 ! MG: Here it is possible to avoid the extra reading if the same k mesh can be used. 1256 ireadwf0=1 ; formeig=0 ; ask_accurate=1 ; optorth=0 1257 mcgq=mpw1*dtset%nspinor*mband_mem_rbz*mkqmem_rbz*dtset%nsppol 1258 !SPr: verified until here, add mcgq for -q 1259 if (one*mpw1*dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol > huge(1)) then 1260 write (msg,'(4a, 5(a,i0), 2a)')& 1261 "Default integer is not wide enough to store the size of the GS wavefunction array (WFKQ, mcgq).",ch10,& 1262 "Action: increase the number of processors. Consider also OpenMP threads.",ch10,& 1263 "nspinor: ",dtset%nspinor, "mpw1: ",mpw1, "mband: ",dtset%mband, "mkqmem_rbz: ",& 1264 mkqmem_rbz, "nsppol: ",dtset%nsppol,ch10,& 1265 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode' 1266 ABI_ERROR(msg) 1267 end if 1268 ABI_MALLOC_OR_DIE(cgq,(2,mcgq), ierr) 1269 1270 ABI_MALLOC(eigenq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1271 if (.not.kramers_deg) then 1272 !ABI_MALLOC_OR_DIE(cg_pq,(2,mcgq), ierr) 1273 !ABI_MALLOC(eigen_pq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1274 mcgmq=mpw1_mq*dtset%nspinor*mband_mem_rbz*mkqmem_rbz*dtset%nsppol 1275 ABI_MALLOC_OR_DIE(cg_mq,(2,mcgmq), ierr) 1276 1277 ABI_MALLOC(eigen_mq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1278 end if 1279 1280 !if (sum(dtset%qptn(1:3)**2)>=1.d-14) then ! non-zero q 1281 !TODO: for many other q this should be avoidable, in principle all if qptrlatt is a subgrid of kprtlatt 1282 ! Or at very least make a pointer instead of a full copy!!! 1283 if (dtfil%fnamewffq == dtfil%fnamewffk .and. sum(dtset%qptn(1:3)**2) < 1.d-14) then 1284 call wrtout(std_out, " qpt is Gamma, psi_k+q initialized from psi_k in memory") 1285 cgq = cg 1286 eigenq = eigen0 1287 else 1288 call timab(144,1,tsec) 1289 call wfk_read_my_kptbands(dtfil%fnamewffq, distrb_flags, spacecomm, dtset%ecut*(dtset%dilatmx)**2,& 1290 & formeig, istwfk_rbz, kpq_rbz, mcgq, dtset%mband, mband_mem_rbz, mkqmem_rbz, mpw1,& 1291 & dtset%natom, nkpt_rbz, npwar1, dtset%nspinor, dtset%nsppol, dtset%usepaw,& 1292 & cgq, eigen=eigenq, occ=occ_disk) 1293 call timab(144,2,tsec) 1294 1295 if (.not.kramers_deg) then 1296 !SPr: later "make" a separate WFQ file for "-q" 1297 call timab(144,1,tsec) 1298 call wfk_read_my_kptbands(dtfil%fnamewffq, distrb_flags, spacecomm,dtset%ecut*(dtset%dilatmx)**2, & 1299 & formeig, istwfk_rbz, kmq_rbz, mcgmq, dtset%mband, mband_mem_rbz, mkqmem_rbz, mpw1_mq,& 1300 & dtset%natom, nkpt_rbz, npwar1_mq, dtset%nspinor, dtset%nsppol, dtset%usepaw,& 1301 & cg_mq, eigen=eigen_mq, occ=occ_disk) 1302 call timab(144,2,tsec) 1303 1304 end if 1305 end if 1306 ABI_FREE(occ_disk) 1307 1308 ! Update energies GS energies at k + q 1309 call put_eneocc_vect(ebands_kq, "eig", eigenq) 1310 if (.not.kramers_deg) then 1311 call put_eneocc_vect(ebands_kmq, "eig", eigen_mq) 1312 end if 1313 1314 1315 ! PAW: compute on-site projections of GS wavefunctions (cprjq) (and derivatives) at k+q 1316 ABI_MALLOC(cprjq,(0,0)) 1317 if (psps%usepaw==1) then 1318 if (usecprj==1) then 1319 mcprjq=dtset%nspinor*mband_mem_rbz*mkqmem_rbz*dtset%nsppol 1320 !TODO : distribute cprj by band as well? 1321 !mcprjq=dtset%nspinor*mband_mem_rbz*mkqmem_rbz*dtset%nsppol 1322 ABI_FREE(cprjq) 1323 ABI_MALLOC(cprjq,(dtset%natom,mcprjq)) 1324 call pawcprj_alloc(cprjq,0,dimcprj_srt) 1325 if (ipert<=dtset%natom.and.(sum(dtset%qptn(1:3)**2)>=1.d-14)) then ! phonons at non-zero q 1326 choice=1 ; iorder_cprj=0 ; idir0=0 1327 call ctocprj(atindx,cgq,choice,cprjq,gmet,gprimd,-1,idir0,0,istwfk_rbz,& 1328 & kg1,kpq_rbz,mcgq,mcprjq,dtset%mgfft,mkqmem_rbz,mpi_enreg,psps%mpsang,mpw1,& 1329 & dtset%natom,nattyp,nband_rbz,dtset%natom,dtset%ngfft,nkpt_rbz,dtset%nloalg,& 1330 & npwar1,dtset%nspinor,dtset%nsppol,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,& 1331 & psps,rmet,dtset%typat,ucvol,dtfil%unpawq,xred,ylm1,ylmgr1) 1332 else if (mcprjq>0) then 1333 call pawcprj_copy(cprj,cprjq) 1334 end if 1335 end if 1336 end if 1337 1338 ! ===== Report on eigenq values 1339 if (dtset%ieig2rf>0.and.icase==ipert_cnt) then 1340 eigen0_pert(:) = eigen0(:) 1341 eigenq_pert(:) = eigenq(:) 1342 occ_rbz_pert(:) = occ_rbz(:) 1343 end if 1344 if (dtset%efmas>0.and.icase==ipert_cnt) then 1345 eigen0_pert(:) = eigen0(:) 1346 end if 1347 !call wrtout(std_out,ch10//' dfpt_looppert: eigenq array',"COLL") 1348 nkpt_eff=nkpt 1349 if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. nkpt>nkpt_max ) nkpt_eff=nkpt_max 1350 band_index=0 1351 do isppol=1,dtset%nsppol 1352 do ikpt=1,nkpt_rbz 1353 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 1354 if(ikpt<=nkpt_eff)then 1355 write(msg, '(a,i2,a,i5)' )' isppol=',isppol,', k point number',ikpt 1356 call wrtout(std_out,msg) 1357 do iband=1,nband_k,4 1358 write(msg, '(a,4es16.6)')' ',eigenq(iband+band_index:min(iband+3,nband_k)+band_index) 1359 call wrtout(std_out,msg) 1360 end do 1361 else if(ikpt==nkpt_eff+1)then 1362 write(msg,'(a,a)' )' respfn : prtvol=0, 1 or 2, stop printing eigenq.',ch10 1363 call wrtout(std_out,msg) 1364 end if 1365 band_index=band_index+nband_k 1366 end do 1367 end do 1368 1369 ! Generate occupation numbers for the reduced BZ at k+q 1370 ABI_MALLOC(docckqde,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1371 ABI_MALLOC(occkq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1372 if (.not.kramers_deg) then 1373 ABI_MALLOC(docckde_mq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1374 ABI_MALLOC(occk_mq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1375 end if 1376 1377 if(0<=dtset%occopt .and. dtset%occopt<=2)then 1378 ! Same occupation numbers at k and k+q (usually, insulating) 1379 occkq(:)=occ_rbz(:) 1380 docckqde(:)=zero ! docckqde is irrelevant in this case 1381 if(.not.kramers_deg) then 1382 occk_mq(:)=occ_rbz(:) 1383 docckde_mq(:)=zero 1384 end if 1385 else 1386 ! Metallic occupation numbers 1387 option=1 1388 dosdeltae=zero ! the DOS is not computed with option=1 1389 maxocc=two/(dtset%nspinor*dtset%nsppol) 1390 call getnel(docckqde,dosdeltae,eigenq,entropy,fermie,fermie,maxocc,dtset%mband,& 1391 nband_rbz,nelectkq,nkpt_rbz,dtset%nsppol,occkq,dtset%occopt,option,& 1392 dtset%tphysel,dtset%tsmear,fake_unit,wtk_rbz,1,dtset%nband(1)) ! CP: added 1, nband(0) to fit new def of getnel 1393 ! Compare nelect at k and nelelect at k+q 1394 write(msg, '(a,a,a,es16.6,a,es16.6,a)')& 1395 ' dfpt_looppert : total number of electrons, from k and k+q',ch10,& 1396 ' fully or partially occupied states are',dtset%nelect,' and',nelectkq,'.' 1397 call wrtout([std_out, ab_out], msg) 1398 if (.not.kramers_deg) then 1399 call getnel(docckde_mq,dosdeltae,eigen_mq,entropy,fermie,fermie,maxocc,dtset%mband,& 1400 nband_rbz,nelectkq,nkpt_rbz,dtset%nsppol,occk_mq,dtset%occopt,option,& 1401 dtset%tphysel,dtset%tsmear,fake_unit,wtk_rbz,1,dtset%nband(1)) 1402 ! Compare nelect at k and nelelect at k-q 1403 write(msg, '(a,a,a,es16.6,a,es16.6,a)')& 1404 ' dfpt_looppert : total number of electrons, from k and k-q',ch10,& 1405 ' fully or partially occupied states are',dtset%nelect,' and',nelectkq,'.' 1406 call wrtout([std_out, ab_out], msg) 1407 end if 1408 end if 1409 1410 ! Debug message 1411 if(dtset%prtvol==-level) call wrtout(std_out,'dfpt_looppert: initialisation of q part done.') 1412 1413 ! Initialisation of first-order wavefunctions 1414 write(msg,'(3a,i4)')' Initialisation of the first-order wave-functions :',ch10,' ireadwf=',dtfil%ireadwf 1415 call wrtout([std_out, ab_out], msg) 1416 call appdig(pertcase,dtfil%fnamewff1,fiwf1i) 1417 call appdig(pertcase,dtfil%fnameabo_1wf,fiwf1o) 1418 1419 ! Allocate 1st-order PAW occupancies (rhoij1) 1420 if (psps%usepaw==1) then 1421 ABI_MALLOC(pawrhoij1,(my_natom)) 1422 call pawrhoij_nullify(pawrhoij1) 1423 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,& 1424 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc) 1425 call pawrhoij_alloc(pawrhoij1,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,& 1426 & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,& 1427 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1428 if (cplex_rhoij/=hdr%pawrhoij(1)%cplex_rhoij.or.qphase_rhoij/=hdr%pawrhoij(1)%qphase) then 1429 ! Eventually reallocate hdr%pawrhoij 1430 call pawrhoij_free(hdr%pawrhoij) 1431 call pawrhoij_alloc(hdr%pawrhoij,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,& 1432 & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,& 1433 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1434 end if 1435 else 1436 ABI_MALLOC(pawrhoij1,(0)) 1437 end if 1438 1439 ! Initialize 1st-order wavefunctions 1440 formeig=1; ask_accurate=0; optorth=0 1441 ! NB: 4 Sept 2013: this was introducing a bug - for ieig2rf ==0 the dim was being set to 1 and passed to dfpt_vtowfk 1442 dim_eig2rf=0 1443 if ((dtset%ieig2rf > 0 .and. dtset%ieig2rf/=2) .or. dtset%efmas > 0) then 1444 dim_eig2rf=1 1445 end if 1446 mcg1=mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol 1447 if (one*mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol > huge(1)) then 1448 write (msg,'(4a, 5(a,i0), 2a)')& 1449 "Default integer is not wide enough to store the size of the GS wavefunction array (WFK1, mcg1).",ch10,& 1450 "Action: increase the number of processors. Consider also OpenMP threads.",ch10,& 1451 "nspinor: ",dtset%nspinor, "mpw1: ",mpw1, "mband: ",dtset%mband, "mk1mem_rbz: ",& 1452 mk1mem_rbz, "nsppol: ",dtset%nsppol,ch10,& 1453 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode' 1454 ABI_ERROR(msg) 1455 end if 1456 ABI_MALLOC_OR_DIE(cg1,(2,mcg1), ierr) 1457 ! space for all 3 ddk wavefunctions if call to orbmag will be needed 1458 if ( (dtset%orbmag .NE. 0) .AND. (dtset%rfddk .EQ. 1) .AND. (.NOT. ALLOCATED(cg1_3)) ) then 1459 ABI_MALLOC(cg1_3,(2,mcg1,3)) 1460 has_cg1_3(:) = .FALSE. 1461 end if 1462 if (.not.kramers_deg) then 1463 mcg1mq=mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol 1464 ABI_MALLOC_OR_DIE(cg1_mq,(2,mcg1mq), ierr) 1465 end if 1466 1467 ABI_MALLOC(cg1_active,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1468 ABI_MALLOC(gh1c_set,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1469 ABI_MALLOC(gh0c1_set,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1470 if (.not.kramers_deg) then 1471 ABI_MALLOC(cg1_active_mq,(2,mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1472 ABI_MALLOC(gh1c_set_mq,(2,mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1473 ABI_MALLOC(gh0c1_set_mq,(2,mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1474 end if 1475 ! XG090606 This is needed in the present 5.8.2 , for portability for the pathscale machine. 1476 ! However, it is due to a bug to be corrected by Paul Boulanger. When the bug will be corrected, 1477 ! this line should be removed. 1478 if(mk1mem_rbz/=0 .and. dtset%ieig2rf/=0)then 1479 cg1_active=zero 1480 gh1c_set=zero 1481 gh0c1_set=zero 1482 if (.not.kramers_deg) then 1483 cg1_active_mq=zero 1484 gh1c_set_mq=zero 1485 gh0c1_set_mq=zero 1486 end if 1487 end if 1488 ! MG TODO: 1489 ! eigen1 is not MPI-distributed, it scales as mband^2 and requires a lot of memory for dense k-meshes e.g. metals: 1490 ! this is what you get with nband 26 and ngkpt 42**3: 1491 ! 1492 ! <var=eigen1, A@m_dfpt_looppert.F90:1499, addr=0x150d756c0010, size_mb=391.204> 1493 ! 1494 ! we should try to MPI-distribute it over nkpt_rbs/nsppol or, at the very least, tell users to activate OpenMP. 1495 1496 ABI_MALLOC(eigen1,(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol)) 1497 ABI_MALLOC(resid,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1498 call timab(144,1,tsec) 1499 if ((file_exists(nctk_ncify(fiwf1i)) .or. file_exists(fiwf1i)) .and. (dtset%get1wf /= 0 .or. dtset%ird1wf /= 0)) then 1500 call wfk_read_my_kptbands(fiwf1i, distrb_flags, spacecomm, dtset%ecut*(dtset%dilatmx)**2,& 1501 formeig, istwfk_rbz, kpq_rbz, mcg1, dtset%mband, mband_mem_rbz, mk1mem_rbz, mpw1,& 1502 dtset%natom, nkpt_rbz, npwar1, dtset%nspinor, dtset%nsppol, dtset%usepaw,& 1503 cg1, eigen=eigen1, ask_accurate_=0) 1504 else 1505 cg1 = zero 1506 eigen1 = zero 1507 end if 1508 1509 call timab(144,2,tsec) 1510 if(.not.kramers_deg) then 1511 ABI_MALLOC(eigen1_mq,(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol)) 1512 ABI_MALLOC(resid_mq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1513 !initialize cg1_mq: 1514 call timab(144,1,tsec) 1515 if ((file_exists(nctk_ncify(fiwf1i)) .or. file_exists(fiwf1i)) .and. (dtset%get1wf > 0 .or. dtset%ird1wf > 0)) then 1516 call wfk_read_my_kptbands(fiwf1i, distrb_flags, spacecomm, dtset%ecut*(dtset%dilatmx)**2, & 1517 formeig, istwfk_rbz, kmq_rbz, mcg1mq, dtset%mband, mband_mem_rbz, mk1mem_rbz, mpw1_mq,& 1518 dtset%natom, nkpt_rbz, npwar1_mq, dtset%nspinor, dtset%nsppol, dtset%usepaw,& 1519 cg1_mq, eigen=eigen1_mq, ask_accurate_=0) 1520 else 1521 cg1_mq = zero 1522 eigen1_mq = zero 1523 end if 1524 call timab(144,2,tsec) 1525 end if 1526 1527 ! Eventually reytrieve 1st-order PAW occupancies from file header 1528 if (psps%usepaw==1.and.dtfil%ireadwf/=0) then 1529 call pawrhoij_copy(hdr%pawrhoij,pawrhoij1,comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 1530 end if 1531 1532 ! In case of electric field, or 2nd order perturbation : open the ddk (or dE) wf file(s) 1533 if ((ipert==dtset%natom+2.and.sum((dtset%qptn(1:3))**2)<1.0d-7.and. & 1534 & (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and. & 1535 & dtset%berryopt/= 7.and.dtset%berryopt/=14.and. & 1536 & dtset%berryopt/=16.and.dtset%berryopt/=17)) .or. & 1537 & ipert==dtset%natom+10.or.ipert==dtset%natom+11) then 1538 if (ipert<dtset%natom+10) then 1539 ! 1st order or berry phase, one direction 1540 nwffile=1 ! one direction needed 1541 file_index(1)=idir+dtset%natom*3 1542 fnamewff(1)=dtfil%fnamewffddk 1543 else if (ipert==dtset%natom+10) then 1544 ! 2nd order (k,k) 1545 if (idir<=3) then ! one direction needed 1546 nwffile=1 1547 file_index(1)=idir+dtset%natom*3 1548 fnamewff(1)=dtfil%fnamewffddk 1549 else ! two directions needed 1550 nwffile=2 1551 file_index(1)=idir1+dtset%natom*3 1552 file_index(2)=idir2+dtset%natom*3 1553 fnamewff(1)=dtfil%fnamewffddk 1554 fnamewff(2)=dtfil%fnamewffddk 1555 end if 1556 else if (ipert==dtset%natom+11) then 1557 ! 2nd order (k,E) 1558 nwffile=3 ! dk, dE and dkdk 1559 idir_dkdk = idir 1560 if(idir_dkdk>6) idir_dkdk = idir_dkdk - 3 1561 file_index(1)=idir_dkdk +(dtset%natom+6)*3 ! dkdk 1562 file_index(2)=idir2+(dtset%natom+1)*3 ! defld file (dir2) 1563 file_index(3)=idir1+dtset%natom*3 ! ddk file (dir1) 1564 fnamewff(1)=dtfil%fnamewffdkdk 1565 fnamewff(2)=dtfil%fnamewffdelfd 1566 fnamewff(3)=dtfil%fnamewffddk 1567 if (idir>3) then 1568 nwffile=4 1569 file_index(4)=idir2+dtset%natom*3 ! ddk file (dir2) 1570 fnamewff(4)=dtfil%fnamewffddk 1571 end if 1572 end if 1573 do ii=1,nwffile 1574 call appdig(file_index(ii),fnamewff(ii),fiwfddk) 1575 ! Checking the existence of data file 1576 if (.not. file_exists(fiwfddk)) then 1577 ! Trick needed to run Abinit test suite in netcdf mode. 1578 if (file_exists(nctk_ncify(fiwfddk))) then 1579 write(msg,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name." 1580 call wrtout(std_out,msg) 1581 fiwfddk = nctk_ncify(fiwfddk) 1582 end if 1583 if (.not. file_exists(fiwfddk)) then 1584 ABI_ERROR('Missing file: '//TRIM(fiwfddk)) 1585 end if 1586 end if 1587 write(msg,'(2a)')'- dfpt_looppert: read the DDK wavefunctions from file: ',trim(fiwfddk) 1588 call wrtout([std_out, ab_out],msg) 1589 ! Note that the unit number for these files is 50,51,52 or 53 (dtfil%unddk=50) 1590 call wfk_open_read(ddk_f(ii),fiwfddk,formeig1,dtset%iomode,dtfil%unddk+(ii-1), spacecomm) !xmpi_comm_self) 1591 end do 1592 end if 1593 1594 ! Get first-order local potentials and 1st-order core correction density change 1595 ! (do NOT include xccc3d1 in vpsp1 : this will be done in dfpt_scfcv because vpsp1 1596 ! might become spin-polarized) 1597 1598 n3xccc=0;if(psps%n1xccc/=0)n3xccc=nfftf 1599 ABI_MALLOC(xccc3d1,(cplex*n3xccc)) 1600 ABI_MALLOC(vpsp1,(cplex*nfftf)) 1601 1602 ! PAW: compute Vloc(1) and core(1) together in reciprocal space 1603 ! -------------------------------------------------------------- 1604 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 1605 ndir=1 1606 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,istr,ipert,& 1607 & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,& 1608 & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 1609 & atmrhor1=xccc3d1,atmvlocr1=vpsp1,optn_in=n3xccc/nfftf,optn2_in=1,vspl=psps%vlspl) 1610 else 1611 1612 ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp. 1613 ! ------------------------------------------------------------------------------ 1614 1615 if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then 1616 ! Section for strain perturbation 1617 1618 !To compute Absolute Deformation Potentials toghether with FxE tensor 1619 !the reference has to be the same as in the FxE routines 1620 g0term=0; if (dtset%rfstrs_ref==1) g0term=1 1621 1622 call vlocalstr(gmet,gprimd,gsqcut,istr,mgfftf,mpi_enreg,& 1623 & psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,ntypat,ph1df,psps%qgrid_vl,& 1624 & ucvol,psps%vlspl,vpsp1,g0term=g0term) 1625 else 1626 call dfpt_vlocal(atindx,cplex,gmet,gsqcut,idir,ipert,mpi_enreg,psps%mqgrid_vl,dtset%natom,& 1627 & nattyp,nfftf,ngfftf,ntypat,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,psps%qgrid_vl,& 1628 & dtset%qptn,ucvol,psps%vlspl,vpsp1,xred) 1629 !SPr: need vpsp1 for -q as well, but for magnetic field it's zero, to be done later 1630 end if 1631 1632 if(psps%n1xccc/=0)then 1633 call dfpt_mkcore(cplex,idir,ipert,dtset%natom,ntypat,ngfftf(1),psps%n1xccc,& 1634 & ngfftf(2),ngfftf(3),dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred) 1635 !SPr: same here, need xccc3d1 for -q as well for phonon pert.. to be done later 1636 end if ! psps%n1xccc/=0 1637 end if ! usepaw 1638 1639 eigen1(:)=zero; resid(:)=zero 1640 if(.not.kramers_deg) then 1641 eigen1_mq(:)=zero 1642 resid_mq(:)=zero 1643 end if 1644 ! Get starting charge density and Hartree + xc potential 1645 ABI_MALLOC(rhor1,(cplex*nfftf,nspden)) 1646 ABI_MALLOC(rhog1,(2,nfftf)) 1647 1648 if(.not.kramers_deg) then 1649 !Case when first order spinors at both +q and -q are not related by symmetry (time and/or space inversion) 1650 ABI_MALLOC(rhor1_pq,(cplex*nfftf,nspden)) 1651 ABI_MALLOC(rhog1_pq,(2,nfftf)) 1652 1653 ABI_MALLOC(rhor1_mq,(cplex*nfftf,nspden)) 1654 ABI_MALLOC(rhog1_mq,(2,nfftf)) 1655 end if 1656 1657 ! can we get this set of gkk matrices from previously calculated rhog1 through a non-scf calculation? 1658 found_eq_gkk=.false. 1659 if (dtset%prepgkk == 1 .and. ipert <= dtset%natom) then 1660 ! NOTE: this does not take into account combinations e.g. x+y -> z 1661 ! if rhor1 add linearly this could be done... 1662 do icase_eq = 1, icase-1 1663 idir_eq = mod(icase_eq,3) 1664 if (idir_eq==0) idir_eq=3 1665 ipert_eq = ( (icase_eq-idir_eq) / 3 + 1) 1666 1667 ! find sym which links old perturbation to present one 1668 do isym=1, nsym 1669 ! check that isym preserves qpt to begin with 1670 if (symq(4,1,isym) /= 1 .or. & 1671 & symq(1,1,isym) /= 0 .or. & 1672 & symq(2,1,isym) /= 0 .or. & 1673 & symq(3,1,isym) /= 0 ) cycle 1674 1675 eq_symop = dtset%symrel(:,:,isym) 1676 if (indsym(4,isym,ipert) == ipert_eq .and. & 1677 & abs(eq_symop(idir,idir_eq)) == 1 .and. & 1678 & sum(abs(eq_symop(:,idir_eq))) == 1) then 1679 found_eq_gkk = .true. 1680 exit 1681 end if 1682 end do ! isym 1683 if (found_eq_gkk) exit 1684 end do ! icase_eq 1685 end if ! check for prepgkk with symmetric pert 1686 1687 if (found_eq_gkk) then 1688 write (msg, '(a,l6,i6,2a,3i6,2a,3i6,a)') & 1689 ' found_eq_gkk,isym = ', found_eq_gkk, isym, ch10, & 1690 ' idir, ipert, icase = ', idir, ipert, icase, ch10, & 1691 ' idireq iperteq icaseeq = ', idir_eq, ipert_eq, icase_eq, ch10 1692 call wrtout(std_out,msg) 1693 ! 1694 ! Make density for present perturbation, which is symmetric of 1 or more previous perturbations: 1695 ! rotate 1DEN arrays with symrel(isym) to produce rhog1_eq rhor1_eq 1696 ! 1697 if (dtset%use_nonscf_gkk == 1) then 1698 call rotate_rho(cplex, timrev_pert, mpi_enreg, nfftf, ngfftf, nspden, & 1699 & rhor1_save(:,:,icase_eq), rhog1, rhor1, eq_symop, dtset%tnons(:,isym)) 1700 1701 rhor1 = rhor1 * eq_symop(idir,idir_eq) 1702 1703 ! TODO: rotate rhoij in PAW case 1704 1705 blkflg_save = blkflg 1706 dtset_tmp%iscf = -2 1707 iscf_mod = -2 1708 dtset_tmp%nstep = 1 1709 dtset_tmp%nline = 1 1710 if (abs(dtset_tmp%tolwfr) < 1.e-24) dtset_tmp%tolwfr = 1.e-24 1711 dtset_tmp%toldfe = zero 1712 dtset_tmp%toldff = zero 1713 dtset_tmp%tolrff = zero 1714 dtset_tmp%tolvrs = zero 1715 write (msg, '(a,i6,a)') ' NOTE: doing GKK calculation for icase ', icase, ' with non-SCF calculation' 1716 call wrtout(std_out,msg) 1717 !call wrtout(ab_out,msg,'COLL') ! decomment and update output files 1718 1719 else ! do not use non-scf shortcut, but save rotated 1DEN for comparison 1720 ! saves the rotated rho, for later comparison with the full SCF rhor1: comment lines below for iscf = -2 1721 call rotate_rho(cplex, timrev_pert, mpi_enreg, nfftf, ngfftf, nspden, & 1722 & rhor1_save(:,:,icase_eq), rhog1, rhor1_save(:,:,icase), eq_symop, & 1723 & dtset%tnons(:,isym)) 1724 rhor1_save(:,:,icase) = rhor1_save(:,:,icase) * eq_symop(idir,idir_eq) 1725 1726 end if ! force non scf calculation of other gkk, or not 1727 end if ! found an equiv perturbation for the gkk 1728 1729 if ( (dtfil%ireadwf==0 .and. iscf_mod/=-4 .and. dtset%get1den==0 .and. dtset%ird1den==0) & 1730 .or. (iscf_mod== -3 .and. ipert/=dtset%natom+11) ) then 1731 ! NOTE : For ipert==natom+11, we want to read the 1st order density from a previous calculation 1732 rhor1(:,:)=zero ; rhog1(:,:)=zero 1733 ! PAW: rhoij have been set to zero in call to pawrhoij_alloc above 1734 1735 init_rhor1 = ((ipert>=1 .and. ipert<=dtset%natom).or.ipert==dtset%natom+5) 1736 ! This section is needed in order to maintain the old behavior and pass the automatic tests 1737 if (psps%usepaw == 0) then 1738 init_rhor1 = init_rhor1 .and. all(psps%nctab(:ntypat)%has_tvale) 1739 else 1740 init_rhor1 = .False. 1741 end if 1742 1743 if (init_rhor1) then 1744 1745 if(ipert/=dtset%natom+5) then 1746 1747 ! Initialize rhor1 and rhog1 from the derivative of atomic densities/gaussians. 1748 ndir = 1; optn2 = 3 1749 if (psps%usepaw==1) then 1750 ! FIXME: Here there's a bug because has_tvale == 0 or 1 instead of 2 1751 if (all(pawtab(:ntypat)%has_tvale/=0)) optn2=2 1752 else if (psps%usepaw==0) then 1753 if (all(psps%nctab(:ntypat)%has_tvale)) optn2=2 1754 end if 1755 1756 if (optn2 == 3) then 1757 call wrtout(std_out," Initializing rhor1 from atom-centered gaussians") 1758 ABI_MALLOC(gauss,(2,ntypat)) 1759 call atom_gauss(ntypat, dtset%densty, psps%ziontypat, psps%znucltypat, gauss) 1760 1761 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,& 1762 mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,& 1763 ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 1764 atmrhor1=rhor1,optn_in=1,optn2_in=3,gauss=gauss) 1765 1766 ABI_FREE(gauss) 1767 else 1768 call wrtout(std_out," Initializing rhor1 from valence densities taken from pseudopotential files") 1769 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,& 1770 mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,& 1771 ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 1772 atmrhor1=rhor1,optn_in=1,optn2_in=2) 1773 end if 1774 1775 else 1776 ! Magnetic field perturbation 1777 call wrtout(std_out," Initializing rhor1 guess based on the ground state XC magnetic field") 1778 1779 call dfpt_init_mag1(ipert,idir,rhor1,rhor,cplex,nfftf,nspden,vxc,kxc,nkxc) 1780 1781 if(.not.kramers_deg) then 1782 rhor1_pq=rhor1 1783 call dfpt_init_mag1(ipert,idir,rhor1_mq,rhor,cplex,nfftf,nspden,vxc,kxc,nkxc) 1784 end if 1785 end if 1786 1787 call fourdp(cplex,rhog1,rhor1,-1,mpi_enreg,nfftf,1,ngfftf,0) 1788 if (.not.kramers_deg) then 1789 !call fourdp(cplex,rhog1_pq,rhor1_pq,-1,mpi_enreg,nfftf,1,ngfftf,0) 1790 rhog1_pq=rhog1 1791 call fourdp(cplex,rhog1_mq,rhor1_mq,-1,mpi_enreg,nfftf,1,ngfftf,0) 1792 end if 1793 end if 1794 1795 else 1796 ! rhor1 not being forced to 0.0 1797 if(iscf_mod>0) then 1798 ! cplex=2 gets the complex density, =1 only real part 1799 if (psps%usepaw==1) then 1800 ! Be careful: in PAW, rho does not include the 1st-order compensation density (to be added in dfpt_scfcv.F90) ! 1801 ABI_MALLOC(rho1wfg,(2,dtset%nfft)) 1802 ABI_MALLOC(rho1wfr,(dtset%nfft,nspden)) 1803 call dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon1,istwfk_rbz,& 1804 kg,kg1,dtset%mband,mband_mem_rbz,dtset%mgfft,mkmem_rbz,mk1mem_rbz,mpi_enreg,mpw,mpw1,nband_rbz,& 1805 dtset%nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,nsym1,& 1806 occ_rbz,phnons1,rho1wfg,rho1wfr,rprimd,symaf1,symrl1,tnons1,ucvol,wtk_rbz) 1807 call transgrid(cplex,mpi_enreg,nspden,+1,1,1,dtset%paral_kgb,pawfgr,rho1wfg,rhog1,rho1wfr,rhor1) 1808 ABI_FREE(rho1wfg) 1809 ABI_FREE(rho1wfr) 1810 else 1811 !SPr: need to modify dfpt_mkrho to taken into account q,-q and set proper formulas when +q and -q spinors are related 1812 call dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon1,istwfk_rbz,& 1813 kg,kg1,dtset%mband,mband_mem_rbz,dtset%mgfft,mkmem_rbz,mk1mem_rbz,mpi_enreg,mpw,mpw1,nband_rbz,& 1814 dtset%nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,nsym1,& 1815 occ_rbz,phnons1,rhog1,rhor1,rprimd,symaf1,symrl1,tnons1,ucvol,wtk_rbz) 1816 end if 1817 1818 else if (.not. found_eq_gkk) then 1819 ! negative iscf_mod and no symmetric rotation of rhor1 1820 ! Read rho1(r) from a disk file and broadcast data. 1821 rdwr=1;rdwrpaw=psps%usepaw;if(dtfil%ireadwf/=0) rdwrpaw=0 1822 if (ipert/=dtset%natom+11) then 1823 call appdig(pertcase,dtfil%fildens1in,fiden1i) 1824 else 1825 ! For ipert==natom+11, we want to read the 1st order density from a previous calculation 1826 call appdig(idir2+(dtset%natom+1)*3,dtfil%fildens1in,fiden1i) 1827 end if 1828 ! call appdig(pertcase,dtfil%fildens1in,fiden1i) 1829 call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor1, & 1830 hdr_den, pawrhoij1, spaceComm, check_hdr=hdr) 1831 etotal = hdr_den%etot; call hdr_den%free() 1832 1833 ! Compute up+down rho1(G) by fft 1834 ABI_MALLOC(work,(cplex*nfftf)) 1835 work(:)=rhor1(:,1) 1836 call fourdp(cplex,rhog1,work,-1,mpi_enreg,nfftf,1,ngfftf,0) 1837 ABI_FREE(work) 1838 end if ! rhor1 generated or read in from file 1839 1840 end if ! rhor1 set to 0 or read in from file 1841 1842 ! Check whether exiting was required by the user. 1843 ! If found then do not start minimization steps 1844 openexit=1 ; if(dtset%chkexit==0) openexit=0 1845 call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit) 1846 ! If immediate exit, and wavefunctions were not read, must zero eigenvalues 1847 if (iexit/=0) eigen1(:)=zero 1848 if (iexit/=0.and.(.not.kramers_deg)) eigen1_mq(:)=zero 1849 1850 if (iexit==0) then 1851 1852 ! Main calculation: get 1st-order wavefunctions from Sternheimer equation (SCF cycle) 1853 ! if ipert==natom+10 or natom+11 : get 2nd-order wavefunctions 1854 if (kramers_deg) then 1855 call dfpt_scfcv(atindx,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,& 1856 & dielt,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset_tmp,& 1857 & d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,& 1858 & ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,eloc0,elpsp1,& 1859 & end0,end1,enl0,enl1,eovl1,epaw1,etotal,evxctau0,evxctau1,evdw,exc1,& 1860 & fermie,gh0c1_set,gh1c_set,hdr,idir,& 1861 & indkpt1,indsy1,initialized,ipert,irrzon1,istwfk_rbz,& 1862 & kg,kg1,kpt_rbz,kxc,mband_mem_rbz,mgfftf,mkmem_rbz,mkqmem_rbz,mk1mem_rbz,& 1863 & mpert,mpi_enreg,mpw,mpw1,mpw1_mq,my_natom,& 1864 & nattyp,nband_rbz,ncpgr,nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,& 1865 & npwarr,npwar1,nspden,& 1866 & nsym1,n3xccc,occkq,occ_rbz,& 1867 & paw_an_pert,paw_ij_pert,pawang,pawang1,pawfgr,pawfgrtab_pert,pawrad,pawrhoij_pert,pawrhoij1,pawtab,& 1868 & pertcase,phnons1,ph1d,ph1df,prtbbb,psps,& 1869 & dtset%qptn,resid,residm,rhog,rhog1,& 1870 & rhor,rhor1,rprimd,symaf1,symrc1,symrl1,tnons1,& 1871 & usecprj,useylmgr,useylmgr1,ddk_f,vpsp1,vtrial,vxc,vxctau,& 1872 & wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,zeff,dfpt_scfcv_retcode,& 1873 & kramers_deg) 1874 else 1875 call dfpt_scfcv(atindx,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,& 1876 & dielt,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset_tmp,& 1877 & d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,& 1878 & ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,eloc0,elpsp1,& 1879 & end0,end1,enl0,enl1,eovl1,epaw1,etotal,evxctau0,evxctau1,evdw,exc1,& 1880 & fermie,gh0c1_set,gh1c_set,hdr,idir,& 1881 & indkpt1,indsy1,initialized,ipert,irrzon1,istwfk_rbz,& 1882 & kg,kg1,kpt_rbz,kxc,mband_mem_rbz,mgfftf,mkmem_rbz,mkqmem_rbz,mk1mem_rbz,& 1883 & mpert,mpi_enreg,mpw,mpw1,mpw1_mq,my_natom,& 1884 & nattyp,nband_rbz,ncpgr,nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,& 1885 & npwarr,npwar1,nspden,& 1886 & nsym1,n3xccc,occkq,occ_rbz,& 1887 & paw_an_pert,paw_ij_pert,pawang,pawang1,pawfgr,pawfgrtab_pert,pawrad,pawrhoij_pert,pawrhoij1,pawtab,& 1888 & pertcase,phnons1,ph1d,ph1df,prtbbb,psps,& 1889 & dtset%qptn,resid,residm,rhog,rhog1,& 1890 & rhor,rhor1,rprimd,symaf1,symrc1,symrl1,tnons1,& 1891 & usecprj,useylmgr,useylmgr1,ddk_f,vpsp1,vtrial,vxc,vxctau,& 1892 & wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,zeff,dfpt_scfcv_retcode,& 1893 & kramers_deg,& 1894 & cg_mq=cg_mq,cg1_mq=cg1_mq,cg1_active_mq=cg1_active_mq,docckde_mq=docckde_mq,eigen_mq=eigen_mq,& 1895 & eigen1_mq=eigen1_mq,gh0c1_set_mq=gh0c1_set_mq,gh1c_set_mq=gh1c_set_mq,& 1896 & kg1_mq=kg1_mq,npwar1_mq=npwar1_mq,occk_mq=occk_mq,resid_mq=resid_mq,residm_mq=residm_mq,& 1897 & rhog1_pq=rhog1_pq,rhog1_mq=rhog1_mq,rhor1_pq=rhor1_pq,rhor1_mq=rhor1_mq) 1898 end if 1899 1900 ! 2nd-order eigenvalues stuff 1901 if (dtset%ieig2rf>0) then 1902 if (first_entry) then 1903 nullify(eigen1_pert) 1904 first_entry = .false. 1905 end if 1906 if (.not.associated(eigen1_pert)) then 1907 ABI_MALLOC(eigen1_pert,(2*dtset%mband**2*nkpt*dtset%nsppol,3,mpert)) 1908 ABI_MALLOC_OR_DIE(cg1_pert,(2,mpw1*nspinor*mband_mem_rbz*mk1mem_rbz*nsppol*dim_eig2rf,3,mpert),ierr) 1909 ABI_MALLOC(gh0c1_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert)) 1910 ABI_MALLOC(gh1c_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert)) 1911 ABI_MALLOC(kpt_rbz_pert,(3,nkpt_rbz)) 1912 ABI_MALLOC(npwarr_pert,(nkpt_rbz,mpert)) 1913 ABI_MALLOC(npwar1_pert,(nkpt_rbz,mpert)) 1914 ABI_MALLOC(npwtot_pert,(nkpt_rbz,mpert)) 1915 eigen1_pert(:,:,:) = zero 1916 cg1_pert(:,:,:,:) = zero 1917 gh0c1_pert(:,:,:,:) = zero 1918 gh1c_pert(:,:,:,:) = zero 1919 npwar1_pert (:,:) = 0 1920 npwarr_pert (:,:) = 0 1921 kpt_rbz_pert = kpt_rbz 1922 end if 1923 clflg(idir,ipert)=1 1924 eigen1_pert(1:2*dtset%mband**2*nkpt_rbz*dtset%nsppol,idir,ipert) = eigen1(:) 1925 if(dtset%ieig2rf==1.or.dtset%ieig2rf==3.or.dtset%ieig2rf==4.or.dtset%ieig2rf==5) then 1926 cg1_pert(:,:,idir,ipert)=cg1_active(:,:) 1927 gh0c1_pert(:,:,idir,ipert)=gh0c1_set(:,:) 1928 gh1c_pert(:,:,idir,ipert)=gh1c_set(:,:) 1929 end if 1930 npwarr_pert(:,ipert)=npwarr(:) 1931 npwar1_pert(:,ipert)=npwar1(:) 1932 npwtot_pert(:,ipert)=npwtot(:) 1933 end if ! eig2rf 1934 1935 ! 2nd-order eigenvalues stuff for EFMAS 1936 if (dtset%efmas>0) then 1937 if (first_entry) then 1938 nullify(eigen1_pert) 1939 first_entry = .false. 1940 end if 1941 if (.not.associated(eigen1_pert)) then 1942 ABI_MALLOC(eigen1_pert,(2*dtset%mband**2*nkpt*dtset%nsppol,3,mpert)) 1943 ABI_MALLOC(cg1_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert)) 1944 ABI_MALLOC(gh0c1_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert)) 1945 ABI_MALLOC(gh1c_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert)) 1946 ABI_MALLOC(kpt_rbz_pert,(3,nkpt_rbz)) 1947 ABI_MALLOC(npwarr_pert,(nkpt_rbz,mpert)) 1948 eigen1_pert(:,:,:) = zero 1949 cg1_pert(:,:,:,:) = zero 1950 gh0c1_pert(:,:,:,:) = zero 1951 gh1c_pert(:,:,:,:) = zero 1952 npwarr_pert (:,:) = 0 1953 kpt_rbz_pert = kpt_rbz 1954 ABI_MALLOC(cg0_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1955 cg0_pert = cg 1956 end if 1957 eigen1_pert(1:2*dtset%mband**2*nkpt_rbz*dtset%nsppol,idir,ipert) = eigen1(:) 1958 cg1_pert(:,:,idir,ipert)=cg1_active(:,:) 1959 gh0c1_pert(:,:,idir,ipert)=gh0c1_set(:,:) 1960 gh1c_pert(:,:,idir,ipert)=gh1c_set(:,:) 1961 npwarr_pert(:,ipert)=npwarr(:) 1962 end if ! efmas 1963 ABI_FREE(gh1c_set) 1964 ABI_FREE(gh0c1_set) 1965 ABI_FREE(cg1_active) 1966 1967 !deallocate bit arrays for case without Kramers' degeneracy 1968 if(.not.kramers_deg) then 1969 ABI_FREE(gh1c_set_mq) 1970 ABI_FREE(gh0c1_set_mq) 1971 ABI_FREE(cg1_active_mq) 1972 end if 1973 1974 end if ! End of the check of hasty exit 1975 1976 call timab(146,1,tsec) 1977 1978 ! Print out message at the end of the iterations 1979 write(msg, '(80a,a,a,a,a)' ) ('=',ii=1,80),ch10,ch10,& 1980 & ' ----iterations are completed or convergence reached----',ch10 1981 call wrtout([std_out, ab_out], msg) 1982 1983 ! Print _gkk file for this perturbation 1984 if (dtset%prtgkk == 1) then 1985 call appdig(3*(ipert-1)+idir,dtfil%fnameabo_gkk,gkkfilnam) 1986 nmatel = dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol 1987 ABI_MALLOC(phasecg, (2, nmatel)) 1988 ! call getcgqphase(dtset, timrev, cg, mcg, cgq, mcgq, mpi_enreg, nkpt_rbz, npwarr, npwar1, phasecg) 1989 phasecg(1,:) = one 1990 phasecg(2,:) = zero 1991 ! NB: phasecg not actually used in outgkk for the moment (2013/08/15) 1992 call outgkk(bantot_rbz, nmatel,gkkfilnam,eigen0,eigen1,hdr0,hdr,mpi_enreg,phasecg) 1993 ABI_FREE(phasecg) 1994 1995 #ifdef HAVE_NETCDF 1996 ! Reshape eigen1 into gkk for netCDF output 1997 ABI_MALLOC_OR_DIE(gkk,(2*dtset%mband*dtset%nsppol,dtset%nkpt,1,1,dtset%mband), ierr) 1998 gkk(:,:,:,:,:) = zero 1999 mband = dtset%mband 2000 band_index = 0 2001 band2tot_index = 0 2002 do isppol=1,dtset%nsppol 2003 do ikpt =1,nkpt_rbz 2004 do iband=1,dtset%mband 2005 do jband=1,dtset%mband 2006 eig1_r = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index) 2007 eig1_i = eigen1(2*jband+(iband-1)*2*mband+band2tot_index) 2008 gkk(2*iband-1+2*band_index,ikpt,1,1,jband) = & 2009 & gkk(2*iband-1+2*band_index,ikpt,1,1,jband) + eig1_r 2010 gkk(2*iband+2*band_index,ikpt,1,1,jband) = & 2011 & gkk(2*iband+2*band_index,ikpt,1,1,jband) + eig1_i 2012 end do !jband 2013 end do !iband 2014 band2tot_index = band2tot_index + 2*mband**2 2015 end do !ikpt 2016 band_index = band_index + mband 2017 end do !isppol 2018 2019 ! Initialize ggk_ebands to write in the GKK.nc file 2020 ! MG FIXME: Here there's a bug because eigen0 is dimensioned with nkpt_rbz i.e. IBZ(q) 2021 ! but the ebands_t object is constructed with dimensions taken from hdr0 i.e. the IBZ(q=0). 2022 bantot= dtset%mband*dtset%nkpt*dtset%nsppol 2023 call ebands_init(bantot,gkk_ebands,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,& 2024 & doccde,eigen0,hdr0%istwfk,hdr0%kptns,& 2025 & hdr0%nband, hdr0%nkpt,hdr0%npwarr,hdr0%nsppol,hdr0%nspinor,& 2026 & hdr0%tphysel,hdr0%tsmear,hdr0%occopt,hdr0%occ,hdr0%wtk,& 2027 & hdr0%cellcharge, hdr0%kptopt, hdr0%kptrlatt_orig, hdr0%nshiftk_orig, hdr0%shiftk_orig, & 2028 & hdr0%kptrlatt, hdr0%nshiftk, hdr0%shiftk) 2029 2030 ! Init a gkk_t object 2031 call gkk_init(gkk,gkk2d,dtset%mband,dtset%nsppol,nkpt_rbz,1,1) 2032 2033 ! Write the netCDF file. 2034 if (me == master) then 2035 fname = strcat(gkkfilnam,".nc") 2036 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKK file") 2037 NCF_CHECK(crystal%ncwrite(ncid)) 2038 NCF_CHECK(ebands_ncwrite(gkk_ebands, ncid)) 2039 call gkk_ncwrite(gkk2d,dtset%qptn(:),dtset%wtq, ncid) 2040 NCF_CHECK(nf90_close(ncid)) 2041 end if 2042 2043 ! Free memory 2044 ABI_FREE(gkk) 2045 call gkk_free(gkk2d) 2046 call ebands_free(gkk_ebands) 2047 #endif 2048 end if 2049 2050 if (dtset%prepgkk == 1 .and. found_eq_gkk) then 2051 if (dtset%use_nonscf_gkk == 1) then 2052 ! Restore old values of SCF cycle parameters 2053 iscf_mod = iscf_mod_save 2054 dtset_tmp%iscf = iscf_mod_save 2055 dtset_tmp%nstep = nstep_save 2056 dtset_tmp%nline = nline_save 2057 dtset_tmp%tolwfr = tolwfr_save 2058 dtset_tmp%toldfe = toldfe_save 2059 dtset_tmp%toldff = toldff_save 2060 dtset_tmp%tolrff = tolrff_save 2061 dtset_tmp%tolvrs = tolvrs_save 2062 blkflg = blkflg_save ! this ensures we do not use the (unconverged) 2DTE from this non scf run 2063 ! Save density for present perturbation, for future use in symmetric perturbations 2064 rhor1_save(:,:,icase) = rhor1 2065 else 2066 write (msg, '(a,3E20.10)') 'norm diff = ', sum(abs(rhor1_save(:,:,icase) - rhor1)), & 2067 & sum(abs(rhor1)), sum(abs(rhor1_save(:,:,icase) - rhor1))/sum(abs(rhor1)) 2068 call wrtout(std_out, msg) 2069 end if 2070 end if 2071 2072 ! Write wavefunctions file only if convergence was not achieved. 2073 write_1wfk = .True. 2074 if (dtset%prtwf == 0) write_1wfk = .False. 2075 if (dtset%prtwf == -1 .and. dfpt_scfcv_retcode == 0) then 2076 write_1wfk = .False. 2077 call wrtout(ab_out," dfpt_looppert: DFPT cycle converged with prtwf=-1. Will skip output of the 1st-order WFK file.") 2078 end if 2079 2080 ! store DDK wavefunctions in memory for later call to orbmag 2081 ! only relevant for DDK pert with orbmag calculation 2082 if( (dtset%orbmag .NE. 0) .AND. (ipert .EQ. dtset%natom+1) ) then 2083 cg1_3(:,:,idir) = cg1(:,:) 2084 has_cg1_3(idir) = .TRUE. 2085 end if 2086 2087 2088 call outresid(dtset,kpt_rbz,dtset%mband, & 2089 & nband_rbz,nkpt_rbz,& 2090 & dtset%nsppol,resid) 2091 2092 if (write_1wfk) then 2093 ! Output 1st-order wavefunctions in file 2094 call wfk_write_my_kptbands(fiwf1o, distrb_flags, spacecomm, formeig, hdr, dtset%iomode, & 2095 & dtset%mband, mband_mem_rbz, mk1mem_rbz, dtset%mpw, nkpt_rbz, dtset%nspinor, dtset%nsppol, & 2096 & cg1, kg1, eigen1) 2097 2098 ! call outwf(cg1,dtset,psps,eigen1,fiwf1o,hdr,kg1,kpt_rbz,& 2099 !& dtset%mband,mcg1,mk1mem_rbz,mpi_enreg,mpw1,dtset%natom,nband_rbz,& 2100 !& nkpt_rbz,npwar1,dtset%nsppol,& 2101 !& occ_rbz,resid,response,dtfil%unwff2,wvl%wfs,wvl%descr) 2102 end if 2103 2104 2105 #ifdef HAVE_NETCDF 2106 ! Output DDK file in netcdf format. 2107 ! Can be used by optic instead of the 1WF file that is really huge. 2108 if (me == master .and. ipert == dtset%natom + 1 .and. dtset%prtevk == 1) then 2109 fname = strcat(dtfil%filnam_ds(4), "_EVK.nc" ) 2110 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EVK.nc file") 2111 ! Have to build hdr on k-grid with info about perturbation. 2112 call hdr_copy(hdr0, hdr_tmp) 2113 hdr_tmp%kptopt = dtset%kptopt 2114 hdr_tmp%pertcase = pertcase 2115 hdr_tmp%qptn = dtset%qptn(1:3) 2116 NCF_CHECK(hdr_tmp%ncwrite(ncid, 43, nc_define=.True.)) 2117 call hdr_tmp%free() 2118 NCF_CHECK(crystal%ncwrite(ncid)) 2119 NCF_CHECK(ebands_ncwrite(ebands_k, ncid)) 2120 ncerr = nctk_def_arrays(ncid, [ & 2121 nctkarr_t('h1_matrix_elements', "dp", & 2122 "two, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.) 2123 NCF_CHECK(ncerr) 2124 NCF_CHECK(nctk_set_datamode(ncid)) 2125 ncerr = nf90_put_var(ncid, nctk_idname(ncid, "h1_matrix_elements"), eigen1, & 2126 count=[2, dtset%mband, dtset%mband, nkpt_rbz, dtset%nsppol]) 2127 NCF_CHECK(ncerr) 2128 NCF_CHECK(nf90_close(ncid)) 2129 end if 2130 #endif 2131 2132 2133 ! If the perturbation is d/dk, evaluate the f-sum rule. 2134 if (ipert==dtset%natom+1 )then 2135 ! Note : the factor of two is related to the difference 2136 ! between Taylor expansion and perturbation expansion 2137 ! Note : this expression should be modified for ecutsm. 2138 ! Indeed, the present one will NOT tend to 1.0_dp. 2139 ek2=gmet(idir,idir)*(two_pi**2)*2.0_dp*dtset%nelect 2140 fsum=-ek1/ek2 2141 if(dtset%ecutsm<tol6)then 2142 write(msg, '(a,es20.10,a,a,es20.10)' ) & 2143 ' dfpt_looppert : ek2=',ek2,ch10,& 2144 ' f-sum rule ratio=',fsum 2145 else 2146 write(msg, '(a,es20.10,a,a,es20.10,a)' ) & 2147 ' dfpt_looppert : ek2=',ek2,ch10,& 2148 ' f-sum rule ratio=',fsum,' (note : ecutsm/=0)' 2149 end if 2150 call wrtout([std_out, ab_out] , msg) 2151 ! Write the diagonal elements of the dH/dk operator, after averaging over degenerate states 2152 ABI_MALLOC(eigen1_mean,(dtset%mband*nkpt_rbz*dtset%nsppol)) 2153 call eigen_meandege(eigen0,eigen1,eigen1_mean,dtset%mband,nband_rbz,nkpt_rbz,dtset%nsppol,1) 2154 option=4 2155 if (me == master) then 2156 call prteigrs(eigen1_mean,dtset%enunit,fermie,fermie,dtfil%fnametmp_1wf1_eig,ab_out,iscf_mod,kpt_rbz,dtset%kptopt,& 2157 & dtset%mband,nband_rbz,dtset%nbdbuf,nkpt_rbz,dtset%nnsclo,dtset%nsppol,occ_rbz,dtset%occopt,& 2158 & option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk_rbz) 2159 end if 2160 ABI_FREE(eigen1_mean) 2161 end if 2162 2163 ! Print the energies 2164 if (dtset%nline/=0 .or. dtset%nstep/=0)then 2165 call dfpt_prtene(dtset%berryopt,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,& 2166 & ehart01,ehart1,eii,ek0,ek1,eloc0,elpsp1,end0,end1,enl0,enl1,eovl1,epaw1,evxctau0,evxctau1,& 2167 & evdw,exc1,has_nd,has_vxctau,ab_out,ipert,dtset%natom,psps%usepaw,usevdw) 2168 end if 2169 2170 ! call orbmag if needed 2171 if ( (dtset%orbmag .NE. 0) .AND. (dtset%rfddk .EQ. 1) .AND. & 2172 & (COUNT(has_cg1_3) .EQ. 3) ) then 2173 2174 if ( .NOT. ALLOCATED(vtrial_local)) then 2175 ABI_MALLOC(vtrial_local,(nfftf,dtset%nspden)) 2176 end if 2177 vtrial_local = vtrial 2178 call orbmag(cg,cg1_3,cprj,dtset,eigen0,gsqcut,kg,mcg,mcg1,mcprj,mkmem_rbz,mpi_enreg,mpw,nfftf,ngfftf,& 2179 & npwarr,occ,paw_an,paw_ij,pawang,pawfgr,pawrad,pawtab,psps,rprimd,vtrial_local,vxctau,xred,ylm,ylmgr) 2180 if( ALLOCATED(vtrial_local) ) then 2181 ABI_FREE(vtrial_local) 2182 end if 2183 if( ALLOCATED(cg1_3) ) then 2184 ABI_FREE(cg1_3) 2185 has_cg1_3(:) = .FALSE. 2186 end if 2187 end if ! end call orbmag 2188 2189 if(mpi_enreg%paral_pert==1) then 2190 if (ipert_me < npert_me -1) then 2191 call hdr0%free() 2192 else 2193 eigen0_copy(1:dtset%mband*nkpt_rbz*dtset%nsppol) = eigen0 2194 end if 2195 ipert_me = ipert_me +1 2196 else 2197 if (icase == ipert_cnt) then 2198 eigen0_copy(1:dtset%mband*nkpt_rbz*dtset%nsppol) = eigen0 2199 else 2200 call hdr0%free() 2201 end if 2202 end if 2203 2204 !Release the temporary arrays (for k, k+q and 1st-order) 2205 ABI_FREE(cg) 2206 ABI_FREE(cgq) 2207 ABI_FREE(cg1) 2208 ABI_FREE(docckqde) 2209 if(.not.kramers_deg) then 2210 ABI_FREE(cg_mq) 2211 ABI_FREE(cg1_mq) 2212 ABI_FREE(docckde_mq) 2213 end if 2214 ABI_FREE(doccde_rbz) 2215 ABI_FREE(eigen0) 2216 ABI_FREE(eigenq) 2217 ABI_FREE(eigen1) 2218 ABI_FREE(kpq) 2219 if(.not.kramers_deg) then 2220 ABI_FREE(eigen_mq) 2221 ABI_FREE(eigen1_mq) 2222 ABI_FREE(kmq) 2223 end if 2224 ABI_FREE(indkpt1) 2225 ABI_FREE(indsy1) 2226 ABI_FREE(istwfk_rbz) 2227 ABI_FREE(irrzon1) 2228 ABI_FREE(kg) 2229 ABI_FREE(kg1) 2230 ABI_FREE(kpq_rbz) 2231 if(.not.kramers_deg) then 2232 ABI_FREE(kg1_mq) 2233 ABI_FREE(kmq_rbz) 2234 end if 2235 ABI_FREE(kpt_rbz) 2236 ABI_FREE(nband_rbz) 2237 ABI_FREE(npwarr) 2238 ABI_FREE(npwar1) 2239 ABI_FREE(npwtot) 2240 ABI_FREE(npwtot1) 2241 ABI_FREE(occkq) 2242 ABI_FREE(occ_rbz) 2243 ABI_FREE(phnons1) 2244 ABI_FREE(resid) 2245 ABI_FREE(rhog1) 2246 ABI_FREE(rhor1) 2247 ABI_FREE(symaf1) 2248 ABI_FREE(symrc1) 2249 ABI_FREE(symrl1) 2250 ABI_FREE(tnons1) 2251 ABI_FREE(wtk_rbz) 2252 ABI_FREE(xccc3d1) 2253 ABI_FREE(vpsp1) 2254 ABI_FREE(ylm) 2255 ABI_FREE(ylm1) 2256 ABI_FREE(ylmgr) 2257 ABI_FREE(ylmgr1) 2258 if(.not.kramers_deg) then 2259 ABI_FREE(npwar1_mq) 2260 ABI_FREE(npwtot1_mq) 2261 ABI_FREE(occk_mq) 2262 ABI_FREE(resid_mq) 2263 ABI_FREE(rhor1_pq) 2264 ABI_FREE(rhor1_mq) 2265 ABI_FREE(rhog1_pq) 2266 ABI_FREE(rhog1_mq) 2267 end if 2268 if (psps%usepaw==1) then 2269 call pawang_free(pawang1) 2270 call pawrhoij_free(pawrhoij1) 2271 if (usecprj==1) then 2272 call pawcprj_free(cprj) 2273 call pawcprj_free(cprjq) 2274 end if 2275 end if 2276 ABI_FREE(pawrhoij1) 2277 ABI_FREE(cprjq) 2278 ABI_FREE(cprj) 2279 if(xmpi_paral==1) then 2280 ABI_FREE(mpi_enreg%proc_distrb) 2281 ABI_FREE(mpi_enreg%my_kpttab) 2282 end if 2283 call hdr%free() 2284 2285 ! Clean band structure datatypes (should use it more in the future !) 2286 call ebands_free(ebands_k) 2287 call ebands_free(ebands_kq) 2288 if(.not.kramers_deg) call ebands_free(ebands_kmq) 2289 2290 ! %%%% Parallelization over perturbations %%%%% 2291 ! *Redefine output/log files 2292 call localredirect(mpi_enreg%comm_cell,mpi_enreg%comm_world,npert_io,mpi_enreg%paral_pert,0) 2293 2294 ABI_FREE(bz2ibz_smap) 2295 ABI_FREE(distrb_flags) 2296 2297 call timab(146,2,tsec) 2298 if(iexit/=0) exit 2299 end do ! End loop on perturbations 2300 ABI_FREE(zeff) 2301 2302 !%%%% Parallelization over perturbations %%%%% 2303 !*Restore default communicators 2304 call unset_pert_comm(mpi_enreg) 2305 !*Gather output/log files 2306 ABI_MALLOC(dyn,(npert_io)) 2307 if (npert_io>0) dyn=1 2308 call localrdfile(mpi_enreg%comm_pert,mpi_enreg%comm_world,.true.,npert_io,mpi_enreg%paral_pert,0,dyn) 2309 ABI_FREE(dyn) 2310 !*Restore PAW on-site data 2311 if (paral_pert_inplace) then 2312 call unset_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,paw_an,paw_ij,pawfgrtab,pawrhoij) 2313 else 2314 call unset_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,& 2315 & paw_an,paw_ij,pawfgrtab,pawrhoij,& 2316 & paw_an_out=paw_an_pert,paw_ij_out=paw_ij_pert,& 2317 & pawfgrtab_out=pawfgrtab_pert,pawrhoij_out=pawrhoij_pert) 2318 end if 2319 2320 !################################################################################# 2321 !Calculate the second-order eigenvalues for a wavevector Q 2322 2323 call timab(147,1,tsec) 2324 smdelta = dtset%smdelta 2325 bdeigrf = dtset%bdeigrf 2326 if(dtset%bdeigrf == -1) bdeigrf = dtset%mband 2327 2328 if(dtset%ieig2rf > 0) then 2329 2330 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 2331 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 2332 call wrtout(std_out,'Reading the dense grid WF file') 2333 ! We get the Abinit header of the file hdr_fine as ouput 2334 ! We get eigenq_fine(mband,hdr_fine%nkpt,hdr_fine%nsppol) as ouput 2335 fname = dtfil%fnameabi_wfkfine 2336 if (dtset%iomode == IO_MODE_ETSF) fname = nctk_ncify(dtfil%fnameabi_wfkfine) 2337 2338 call wfk_read_eigenvalues(fname,eigenq_fine,hdr_fine,mpi_enreg%comm_world) 2339 ABI_CHECK(SIZE(eigenq_fine,DIM=1)==Dtset%mband,"Size eigenq_fine != mband") 2340 end if 2341 ! DBSP ==> Has been changed to be able to make Bandstructure calculation 2342 ! if(dtset%kptopt==3 .or. dtset%kptopt==0)then 2343 if(dtset%kptopt==3 .or. dtset%kptopt==0 .or. dtset%kptopt < -4 .or. dtset%nsym==1) then 2344 !END 2345 if (dtset%nsym > 1) then ! .and. dtset%efmas==0) then 2346 ABI_ERROR("Symmetries are not implemented for temperature dependence calculations") 2347 end if 2348 write(std_out,*) 'Entering: eig2stern' 2349 if(smdelta>0)then 2350 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 2351 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 2352 call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,& 2353 & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,& 2354 & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,& 2355 & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,eigbrd,eigenq_fine,hdr_fine,hdr0) 2356 else 2357 call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,& 2358 & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,& 2359 & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,& 2360 & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,eigbrd) 2361 end if 2362 else 2363 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 2364 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 2365 call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,& 2366 & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,& 2367 & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,& 2368 & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,eigbrd,eigenq_fine,hdr_fine,hdr0) 2369 else 2370 call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,& 2371 & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,& 2372 & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,& 2373 & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset) 2374 end if 2375 end if 2376 call wrtout(std_out, 'Leaving: eig2stern') 2377 2378 2379 ! ------------------------- 2380 ! Output d2eig data to file 2381 ! ------------------------- 2382 if (dtset%ieig2rf==1.or.dtset%ieig2rf==2) then 2383 2384 ! GA: Here, mpert needs to be replaced by natom 2385 ! but why is mpert larger than natom in the first place? 2386 mpert_ = dtset%natom 2387 2388 ! Initialize perturbation flags 2389 ! GA: At the moment, they are all set to one 2390 ! Instead, they should be used to save individual perturbations 2391 ! to separate files and merge them after the loop. 2392 ABI_MALLOC(blkflg_save,(3,mpert_,3,mpert_)) 2393 blkflg_save = one 2394 2395 ! Initialize ddb object 2396 call ddb%init(dtset, 1, mpert_, & 2397 mband=bdeigrf,& 2398 nkpt=nkpt_rbz,& 2399 kpt=dtset%kptns(1:3,1:nkpt_rbz),& 2400 with_d2eig=.true.) 2401 2402 ! Create the ddb header 2403 dscrpt=' Note : temporary (transfer) database ' 2404 call ddb_hdr%init(dtset,psps,pawtab,dscrpt,1,& 2405 mpert=mpert_,& 2406 xred=xred,occ=occ_pert,& 2407 mband=bdeigrf / dtset%nsppol,& 2408 nkpt=nkpt_rbz,& 2409 kpt=dtset%kptns(:,1:nkpt_rbz)) 2410 2411 ! Set d2eig data 2412 call ddb%set_qpt(1, dtset%qptn) 2413 call ddb%set_d2eig_reshape(1, eig2nkq, blkflg_save) 2414 2415 call ddb_hdr%set_typ(ddb%nblok, ddb%typ) 2416 2417 ! Open the file and write header 2418 call ddb_hdr%open_write(dtfil%fnameabo_eigr2d, with_psps=1, comm=mpi_enreg%comm_world) 2419 2420 ! Write d2eig data block 2421 call ddb%write_d2eig(ddb_hdr, 1, comm=mpi_enreg%comm_world) 2422 2423 ! close and free memory 2424 call ddb_hdr%close() 2425 call ddb_hdr%free() 2426 call ddb%free() 2427 2428 if(smdelta>0) then 2429 ! write out _EIGI2D file 2430 2431 call ddb%init(dtset, 1, mpert_, & 2432 mband=bdeigrf,& 2433 nkpt=nkpt_rbz,& 2434 kpt=dtset%kptns(:,1:nkpt_rbz),& 2435 with_d2eig=.true.) 2436 2437 ! Create the ddb header 2438 dscrpt=' Note : temporary (transfer) database ' 2439 call ddb_hdr%init(dtset,psps,pawtab,dscrpt,1,& 2440 mpert=mpert_,& 2441 xred=xred,occ=occ_pert,& 2442 mband=bdeigrf / dtset%nsppol,& 2443 nkpt=nkpt_rbz,& 2444 kpt=dtset%kptns(:,1:nkpt_rbz)) 2445 2446 call ddb%set_qpt(1, dtset%qptn) 2447 call ddb%set_d2eig_reshape(1, eigbrd, blkflg_save, blktyp=BLKTYP_d2eig_im) 2448 2449 call ddb_hdr%set_typ(ddb%nblok, ddb%typ) 2450 2451 call ddb_hdr%open_write(dtfil%fnameabo_eigi2d, with_psps=1,comm=mpi_enreg%comm_world) 2452 call ddb%write_d2eig(ddb_hdr, 1, comm=mpi_enreg%comm_world) 2453 2454 call ddb_hdr%close() 2455 call ddb_hdr%free() 2456 call ddb%free() 2457 2458 end if !smdelta 2459 2460 ABI_FREE(blkflg_save) 2461 2462 end if !ieig2rf==1.or.ieig2rf==2 2463 else 2464 write(msg,'(3a)')& 2465 'K point grids must be the same for every perturbation: eig2stern not called',ch10,& 2466 'Action: Put kptopt=3 ' 2467 ABI_WARNING(msg) 2468 end if !kptopt 2469 ABI_FREE(gh1c_pert) 2470 ABI_FREE(gh0c1_pert) 2471 ABI_FREE(cg1_pert) 2472 ABI_FREE(kpt_rbz_pert) 2473 ABI_FREE(istwfk_pert) 2474 ABI_FREE(npwarr_pert) 2475 ABI_FREE(npwar1_pert) 2476 ABI_FREE(npwtot_pert) 2477 ABI_FREE(occ_pert) 2478 end if !if dtset%ieig2rf 2479 2480 ! Calculation of effective masses. 2481 if(dtset%efmas == 1) then 2482 call efmas_main(cg0_pert,cg1_pert,dim_eig2rf,dtset,efmasdeg,efmasval,eigen0_pert,& 2483 & eigen1_pert,gh0c1_pert,gh1c_pert,istwfk_pert,mpert,mpi_enreg,nkpt_rbz,npwarr_pert,rprimd) 2484 2485 ABI_FREE(gh1c_pert) 2486 ABI_FREE(gh0c1_pert) 2487 ABI_FREE(cg1_pert) 2488 ABI_FREE(istwfk_pert) 2489 ABI_FREE(npwarr_pert) 2490 ABI_FREE(cg0_pert) 2491 2492 if (dtset%prtefmas == 1 .and. me == master) then 2493 fname = strcat(dtfil%filnam_ds(4), "_EFMAS.nc") 2494 #ifdef HAVE_NETCDF 2495 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EFMAS file") 2496 NCF_CHECK(crystal%ncwrite(ncid)) 2497 !NCF_CHECK(ebands_ncwrite(ebands_k, ncid)) ! At this stage, ebands_k is not available 2498 call print_efmas(efmasdeg, efmasval, kpt_rbz_pert, ncid) 2499 NCF_CHECK(nf90_close(ncid)) 2500 #endif 2501 endif 2502 2503 call efmas_analysis(dtset,efmasdeg,efmasval,kpt_rbz_pert,mpi_enreg,nkpt_rbz,rprimd) 2504 ABI_FREE(kpt_rbz_pert) 2505 end if 2506 2507 !Free memory. 2508 if(dtset%ieig2rf /= 3 .and. dtset%ieig2rf /= 4 .and. dtset%ieig2rf /= 5) call hdr0%free() 2509 ABI_FREE(eigen0_copy) 2510 call crystal%free() 2511 2512 call timab(147,2,tsec) 2513 !###################################################################################### 2514 2515 !Get ddk file information, for later use in dfpt_dyout 2516 ddkfil(:)=0 2517 do idir=1,3 2518 file_index(1)=idir+dtset%natom*3 2519 call appdig(file_index(1),dtfil%fnamewffddk,fiwfddk) 2520 ! Check that ddk file exists 2521 t_exist = file_exists(fiwfddk) 2522 if (.not. t_exist) then 2523 ! Trick needed to run Abinit test suite in netcdf mode. 2524 t_exist = file_exists(nctk_ncify(fiwfddk)) 2525 if (t_exist) then 2526 write(msg,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name." 2527 call wrtout(std_out, msg) 2528 fiwfddk = nctk_ncify(fiwfddk) 2529 end if 2530 end if 2531 2532 ! If the file exists set ddkfil to a non-zero value 2533 if (t_exist) ddkfil(idir)=20+idir 2534 end do 2535 2536 ABI_FREE(ph1d) 2537 ABI_FREE(ph1df) 2538 ABI_FREE(pert_calc) 2539 if (psps%usepaw==1) then 2540 ABI_FREE(dimcprj_srt) 2541 end if 2542 2543 !destroy dtset_tmp 2544 if (dtset%prepgkk /= 0) then ! .and. dtset%use_nonscf_gkk == 1) then !Later uncomment this - in scf case rhor1_save is used below only for testing 2545 ABI_FREE(rhor1_save) 2546 ABI_FREE(blkflg_save) 2547 call dtset_tmp%free() 2548 ABI_FREE(dtset_tmp) 2549 end if 2550 2551 !In paral_pert-case some array's have to be reconstructed 2552 if(mpi_enreg%paral_pert==1) then 2553 ABI_MALLOC(buffer1,(2,3,mpert,3,mpert*(2+psps%usepaw))) 2554 buffer1(:,:,:,:,1:mpert)=d2lo(:,:,:,:,:) 2555 buffer1(:,:,:,:,1+mpert:2*mpert)=d2nl(:,:,:,:,:) 2556 if (psps%usepaw==1) then 2557 buffer1(:,:,:,:,1+2*mpert:3*mpert)=d2ovl(:,:,:,:,:) 2558 end if 2559 call xmpi_sum(buffer1,mpi_enreg%comm_pert,ierr) 2560 call xmpi_sum(blkflg,mpi_enreg%comm_pert,ierr) 2561 if(dtset%prtbbb==1) then 2562 call xmpi_sum(d2bbb,mpi_enreg%comm_pert,ierr) 2563 end if 2564 d2lo(:,:,:,:,:)=buffer1(:,:,:,:,1:mpert) 2565 d2nl(:,:,:,:,:)=buffer1(:,:,:,:,1+mpert:2*mpert) 2566 if (psps%usepaw==1) then 2567 d2ovl(:,:,:,:,:)=buffer1(:,:,:,:,1+2*mpert:3*mpert) 2568 end if 2569 ABI_FREE(buffer1) 2570 end if 2571 2572 if ( associated(old_atmtab)) then 2573 ABI_FREE(old_atmtab) 2574 nullify(old_atmtab) 2575 end if 2576 2577 call timab(141,2,tsec) 2578 2579 DBG_EXIT("COLL") 2580 2581 end subroutine dfpt_looppert
ABINIT/dfpt_prtene [ Functions ]
NAME
dfpt_prtene
FUNCTION
Print components of second derivative of total energy in nice format
INPUTS
eberry=energy associated with Berry phase edocc=correction to 2nd-order total energy coming from changes of occupation eeig0=0th-order eigenenergies part of 2nd-order total energy eew=Ewald part of 2nd-order total energy efrhar=hartree frozen-wavefunction part of 2nd-order tot. en. efrkin=kinetic frozen-wavefunction part of 2nd-order tot. en. efrloc=local psp. frozen-wavefunction part of 2nd-order tot. en. efrnl=nonlocal psp. frozen-wavefunction part of 2nd-order tot. en efrx1=xc core corr.(1) frozen-wavefunction part of 2nd-order tot. en efrx2=xc core corr.(2) frozen-wavefunction part of 2nd-order tot. en ehart01=inhomogeneous 1st-order Hartree part of 2nd-order total energy for strain perturbation only (zero otherwise, and not used) ehart1=1st-order Hartree part of 2nd-order total energy eii=pseudopotential core part of 2nd-order total energy ek0=0th-order kinetic energy part of 2nd-order total energy. ek1=1st-order kinetic energy part of 2nd-order total energy. eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy elpsp1=1st-order local pseudopot. part of 2nd-order total energy. end0=0th-order nuclear dipole part of 2nd-order total energy. end1=1st-order nuclear dipole part of 2nd-order total energy. enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy. enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy. eovl1=1st-order change of wave-functions overlap, part of 2nd-order energy PAW only - Eq(79) and Eq(80) of PRB 78, 035105 (2008) [[cite:Audouze2008]] epaw1=1st-order PAW on-site part of 2nd-order total energy. evxctau0=0th-order energy due to vxctau evxctau1=1st order energy due to vxctau evdw=DFT-D semi-empirical part of 2nd-order total energy exc1=1st-order exchange-correlation part of 2nd-order total energy has_nd=logical, whether nuclear dipole energies are present has_vxctau=logical, whether mgga vxctau energies are present iout=unit number to which output is written ipert=type of the perturbation natom=number of atoms in unit cell usepaw= 0 for non paw calculation; =1 for paw calculation usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use
OUTPUT
(only writing)
NOTES
all energies in Hartree
SOURCE
2816 subroutine dfpt_prtene(berryopt,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,& 2817 & ehart01,ehart1,eii,ek0,ek1,eloc0,elpsp1,end0,end1,enl0,enl1,eovl1,epaw1,evxctau0,evxctau1,& 2818 & evdw,exc1,has_nd,has_vxctau,iout,ipert,natom,usepaw,usevdw) 2819 2820 !Arguments ------------------------------- 2821 !scalars 2822 integer,intent(in) :: berryopt,iout,ipert,natom,usepaw,usevdw 2823 real(dp),intent(in) :: eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1 2824 real(dp),intent(in) :: efrx2,ehart01,ehart1,eii,ek0,ek1,eloc0,elpsp1,end0,end1,enl0,enl1 2825 real(dp),intent(in) :: eovl1,epaw1,evxctau0,evxctau1,evdw,exc1 2826 logical,intent(in) :: has_nd,has_vxctau 2827 2828 !Local variables ------------------------- 2829 !scalars 2830 integer :: nn 2831 logical :: berry_activated 2832 real(dp) :: enl1_effective,erelax,etotal 2833 character(len=10) :: numb 2834 character(len=10),parameter :: numbstr(20) = & 2835 & (/'One ','Two ','Three ','Four ','Five ', & 2836 & 'Six ','Seven ','Eight ','Nine ','Ten ', & 2837 & 'Eleven ','Twelve ','Thirteen ','Fourteen ','Fifteen ', & 2838 & 'Sixteen ','Seventeen ','Eighteen ','Nineteen ','Twenty '/) 2839 character(len=500) :: msg 2840 2841 ! ********************************************************************* 2842 2843 !Count and print the number of components of 2nd-order energy 2844 !MT feb 2015: this number is wrong! Should change it but 2845 ! need to change a lot of ref. files 2846 berry_activated=(berryopt== 4.or.berryopt== 6.or.berryopt== 7.or. & 2847 & berryopt==14.or.berryopt==16.or.berryopt==17) 2848 if (ipert==natom+1) nn=8 2849 if (ipert==natom+5) nn=8 2850 if (ipert==natom+2) nn=7 2851 if (ipert>=1.and.ipert<=natom) nn=13 2852 if (ipert==natom+3.or.ipert==natom+4) nn=17 2853 if (ipert==natom+2.and.berry_activated) nn=nn+1 2854 if (ipert==natom+10.or.ipert==natom+11) nn=1 ! means nothing, 2855 ! because we do not compute derivatives of the energy in this case 2856 if (usepaw==1) nn=nn+1 2857 if (usevdw==1) nn=nn+1 2858 write(msg, '(4a)' ) ch10,& 2859 & ' ',trim(numbstr(nn)),' components of 2nd-order total energy (hartree) are ' 2860 call wrtout(iout,msg) 2861 call wrtout(std_out,msg) 2862 2863 numb='1,2,3' 2864 write(msg, '(3a)' )& 2865 & ' ',trim(numb),': 0th-order hamiltonian combined with 1st-order wavefunctions' 2866 call wrtout(iout,msg) 2867 call wrtout(std_out,msg) 2868 write(msg, '(a,es17.8,a,es17.8,a,es17.8)' )& 2869 & ' kin0=',ek0, ' eigvalue=',eeig0,' local=',eloc0 2870 if (has_nd) then 2871 write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,es17.8)' )& 2872 & ' kin0=',ek0, ' eigvalue=',eeig0,' local=',eloc0,' nclr dpl0=',end0 2873 end if 2874 if (has_vxctau) then 2875 write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,es17.8)' )& 2876 & ' kin0=',ek0, ' eigvalue=',eeig0,' local=',eloc0,' evxctau0=',evxctau0 2877 end if 2878 call wrtout(iout,msg) 2879 call wrtout(std_out,msg) 2880 2881 numb='4,5,6';if( ipert==natom+3.or.ipert==natom+4) numb='4,5,6,7' 2882 write(msg, '(3a)' )& 2883 & ' ',trim(numb),': 1st-order hamiltonian combined with 1st and 0th-order wfs' 2884 call wrtout(iout,msg) 2885 call wrtout(std_out,msg) 2886 if(ipert/=natom+1.and.ipert/=natom+2)then 2887 write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,a)' ) & 2888 & ' loc psp =',elpsp1,' Hartree=',ehart1,' xc=',exc1,ch10,& 2889 & ' note that "loc psp" includes a xc core correction that could be resolved' 2890 else if(ipert==natom+1) then 2891 write(msg, '(a,es17.8,a,es17.8,a,es17.8)' ) & 2892 & ' kin1=',ek1, ' Hartree=',ehart1,' xc=',exc1 2893 if (has_nd) then 2894 write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,es17.8)' ) & 2895 & ' kin1=',ek1, ' Hartree=',ehart1,' xc=',exc1,' nclr dpl1=',end1 2896 end if 2897 if (has_vxctau) then 2898 write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,es17.8)' ) & 2899 & ' kin1=',ek1, ' Hartree=',ehart1,' xc=',exc1,' evxctau1=',evxctau1 2900 end if 2901 else if(ipert==natom+2) then 2902 write(msg, '(a,es17.8,a,es17.8,a,es17.8)' ) & 2903 & ' dotwf=',enl1, ' Hartree=',ehart1,' xc=',exc1 2904 end if 2905 if(ipert==natom+3 .or. ipert==natom+4) then 2906 write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,a,es17.8)' ) & 2907 & ' loc psp =',elpsp1,' Hartree=',ehart1,' xc=',exc1,ch10,& 2908 & ' kin1=',ek1 2909 end if 2910 call wrtout(iout,msg) 2911 call wrtout(std_out,msg) 2912 2913 enl1_effective=enl1;if (ipert==natom+2) enl1_effective=zero 2914 numb='7,8,9';if( ipert==natom+3.or.ipert==natom+4) numb='8,9,10' 2915 write(msg, '(5a,es17.8,a,es17.8,a,es17.8)' )& 2916 & ' ',trim(numb),': eventually, occupation + non-local contributions',ch10,& 2917 & ' edocc=',edocc,' enl0=',enl0,' enl1=',enl1_effective 2918 call wrtout(iout,msg) 2919 call wrtout(std_out,msg) 2920 2921 if (usepaw==1) then 2922 numb='10';if( ipert==natom+3.or.ipert==natom+4) numb='11' 2923 write(msg, '(3a,es17.8)' )& 2924 & ' ',trim(numb),': eventually, PAW "on-site" Hxc contribution: epaw1=',epaw1 2925 call wrtout(iout,msg) 2926 call wrtout(std_out,msg) 2927 end if 2928 2929 if(ipert/=natom+10 .and.ipert/=natom+11) then 2930 erelax=0.0_dp 2931 if(ipert>=1.and.ipert<=natom)then 2932 erelax=ek0+edocc+eeig0+eloc0+elpsp1+ehart1+exc1+enl0+enl1+epaw1 2933 else if(ipert==natom+1.or.ipert==natom+2)then 2934 ! JWZ: end0 and evxctau0 are included as "local" in getghc, while 2935 ! JWZ: end1 and evxctau1 are included in gvnlx1 (non local) in getgh1c 2936 ! JWZ: therefore, in erelax, end0 and evxctau0 are added in explicitly while 2937 ! JWZ: end1 and evxctau1 are already present in enl1 2938 erelax=ek0+edocc+eeig0+eloc0+ek1+ehart1+exc1+enl0+enl1+epaw1+end0+evxctau0 2939 else if(ipert==natom+3.or.ipert==natom+4)then 2940 erelax=ek0+edocc+eeig0+eloc0+ek1+elpsp1+ehart1+exc1+enl0+enl1+epaw1 2941 else if(ipert==natom+5)then 2942 erelax=ek0+edocc+eeig0+eloc0+ek1+elpsp1+ehart1+exc1+enl0+enl1+epaw1 2943 end if 2944 enl1_effective=enl1 2945 if (ipert==natom+1.or.ipert==natom+2) then 2946 if (1.0_dp+enl1/10.0_dp==1.0_dp) enl1_effective=zero 2947 end if 2948 2949 numb='1-9';if (usepaw==1) numb='1-10' 2950 if( ipert==natom+3.or.ipert==natom+4) then 2951 numb='1-10';if (usepaw==1) numb='1-11' 2952 end if 2953 write(msg, '(5a,es17.8)' )& 2954 & ' ',trim(numb),' gives the relaxation energy (to be shifted if some occ is /=2.0)',& 2955 & ch10,' erelax=',erelax 2956 call wrtout(iout,msg) 2957 call wrtout(std_out,msg) 2958 end if 2959 2960 if(ipert>=1.and.ipert<=natom)then 2961 2962 numb='10,11,12';if (usepaw==1) numb='11,12,13' 2963 write(msg, '(4a)' )& 2964 & ' ',trim(numb),' Non-relaxation contributions : ',& 2965 & 'frozen-wavefunctions and Ewald' 2966 call wrtout(iout,msg) 2967 call wrtout(std_out,msg) 2968 write(msg, '(a,es17.8,a,es17.8,a,es17.8)' ) & 2969 & ' fr.local=',efrloc,' fr.nonlo=',efrnl,' Ewald=',eew 2970 call wrtout(iout,msg) 2971 call wrtout(std_out,msg) 2972 2973 write(msg, '(a,es16.6)' )' dfpt_prtene : non-relax=',efrloc+efrnl+eew 2974 call wrtout(std_out,msg) 2975 2976 numb='13,14';if (usepaw==1) numb='14,15' 2977 write(msg, '(3a)' )& 2978 & ' ',trim(numb),' Frozen wf xc core corrections (1) and (2)' 2979 call wrtout(iout,msg) 2980 call wrtout(std_out,msg) 2981 write(msg, '(a,es17.8,a,es17.8)' ) & 2982 & ' frxc 1 =',efrx1,' frxc 2 =',efrx2 2983 call wrtout(iout,msg) 2984 call wrtout(std_out,msg) 2985 if (usepaw==1) then 2986 numb='16' 2987 write(msg, '(5a,es17.8)' )& 2988 & ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',& 2989 & ch10,' eovl1 =',eovl1 2990 call wrtout(iout,msg) 2991 call wrtout(std_out,msg) 2992 end if 2993 if (usevdw==1) then 2994 numb='15';if (usepaw==1) numb='17' 2995 write(msg, '(3a,es17.8)' )& 2996 & ' ',trim(numb),' Contribution from van der Waals DFT-D: evdw =',evdw 2997 call wrtout(iout,msg) 2998 call wrtout(std_out,msg) 2999 end if 3000 3001 write(msg, '(a)' )' Resulting in : ' 3002 call wrtout(iout,msg) 3003 call wrtout(std_out,msg) 3004 etotal=erelax+eew+efrloc+efrnl+efrx1+efrx2+evdw 3005 write(msg, '(a,e20.10,a,e22.12,a)' ) & 3006 & ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV' 3007 call wrtout(iout,msg) 3008 call wrtout(std_out,msg) 3009 write(msg, '(a,es20.10,a,es20.10,a)' ) & 3010 & ' (2DErelax=',erelax,' Ha. 2DEnonrelax=',etotal-erelax,' Ha)' 3011 call wrtout(iout,msg) 3012 call wrtout(std_out,msg) 3013 write(msg, '(a,es20.10,a,a)' ) & 3014 & ' ( non-var. 2DEtotal :',& 3015 & 0.5_dp*(elpsp1+enl1)+eovl1+eew+efrloc+efrnl+efrx1+efrx2+evdw,' Ha)',ch10 3016 call wrtout(iout,msg) 3017 call wrtout(std_out,msg) 3018 3019 else if(ipert==natom+1.or.ipert==natom+2)then 3020 if (ipert==natom+1.and.usepaw==1) then 3021 numb='11' 3022 write(msg, '(5a,es17.8)' )& 3023 & ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',& 3024 & ch10,' eovl1 =',eovl1 3025 call wrtout(iout,msg) 3026 call wrtout(std_out,msg) 3027 end if 3028 write(msg,*)' No Ewald or frozen-wf contrib.:',' the relaxation energy is the total one' 3029 if(berry_activated) then 3030 numb='10'; 3031 write(msg,'(3a,es20.10)')' ',trim(numb),' Berry phase energy :',eberry 3032 end if 3033 call wrtout(iout,msg) 3034 call wrtout(std_out,msg) 3035 etotal=erelax 3036 write(msg, '(a,e20.10,a,e22.12,a)' ) & 3037 & ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV' 3038 call wrtout(iout,msg) 3039 call wrtout(std_out,msg) 3040 write(msg, '(a,es20.10,a)' ) & 3041 & ' ( non-var. 2DEtotal :',0.5_dp*(ek1+enl1_effective)+eovl1,' Ha)' 3042 call wrtout(iout,msg) 3043 call wrtout(std_out,msg) 3044 3045 else if(ipert==natom+3 .or. ipert==natom+4) then 3046 numb='11,12,13';if (usepaw==1) numb='12,13,14' 3047 write(msg, '(4a)' )& 3048 & ' ',trim(numb),' Non-relaxation contributions : ','frozen-wavefunctions and Ewald' 3049 call wrtout(iout,msg) 3050 call wrtout(std_out,msg) 3051 write(msg, '(a,es17.8,a,es17.8,a,es17.8)' ) & 3052 & ' fr.hart=',efrhar,' fr.kin=',efrkin,' fr.loc=',efrloc 3053 call wrtout(iout,msg) 3054 call wrtout(std_out,msg) 3055 3056 numb='14,15,16';if (usepaw==1) numb='15,16,17' 3057 write(msg, '(4a)' )& 3058 & ' ',trim(numb),' Non-relaxation contributions : ','frozen-wavefunctions and Ewald' 3059 call wrtout(iout,msg) 3060 call wrtout(std_out,msg) 3061 write(msg, '(a,es17.8,a,es17.8,a,es17.8)' ) & 3062 & ' fr.nonl=',efrnl,' fr.xc=',efrx1,' Ewald=',eew 3063 call wrtout(iout,msg) 3064 call wrtout(std_out,msg) 3065 3066 numb='17';if (usepaw==1) numb='18' 3067 write(msg, '(4a)' )& 3068 & ' ',trim(numb),' Non-relaxation contributions : ','pseudopotential core energy' 3069 call wrtout(iout,msg) 3070 call wrtout(std_out,msg) 3071 write(msg, '(a,es17.8)' ) ' pspcore=',eii 3072 call wrtout(iout,msg) 3073 call wrtout(std_out,msg) 3074 if (usepaw==1) then 3075 numb='19' 3076 write(msg, '(5a,es17.8)' )& 3077 & ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',& 3078 & ch10,' eovl1 =',eovl1 3079 call wrtout(iout,msg) 3080 call wrtout(std_out,msg) 3081 end if 3082 if (usevdw==1) then 3083 numb='18';if (usepaw==1) numb='20' 3084 write(msg, '(3a,es17.8)' )& 3085 & ' ',trim(numb),' Contribution from van der Waals DFT-D: evdw =',evdw 3086 call wrtout(iout,msg) 3087 call wrtout(std_out,msg) 3088 end if 3089 3090 write(msg, '(a,es16.6)' )' dfpt_prtene : non-relax=',& 3091 & efrhar+efrkin+efrloc+efrnl+efrx1+eew+evdw 3092 call wrtout(std_out,msg) 3093 write(msg, '(a)' )' Resulting in : ' 3094 call wrtout(iout,msg) 3095 call wrtout(std_out,msg) 3096 etotal=erelax+efrhar+efrkin+efrloc+efrnl+efrx1+eew+eii+evdw 3097 write(msg, '(a,e20.10,a,e22.12,a)' ) & 3098 & ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV' 3099 call wrtout(iout,msg) 3100 call wrtout(std_out,msg) 3101 write(msg, '(a,es20.10,a,es20.10,a)' ) & 3102 & ' (2DErelax=',erelax,' Ha. 2DEnonrelax=',etotal-erelax,' Ha)' 3103 call wrtout(iout,msg) 3104 call wrtout(std_out,msg) 3105 write(msg, '(a,es20.10,a,a)' ) & 3106 & ' ( non-var. 2DEtotal :',& 3107 & 0.5_dp*(elpsp1+enl1+ek1+ehart01)+eovl1+& 3108 & efrhar+efrkin+efrloc+efrnl+efrx1+eew+eii+evdw,' Ha)',ch10 3109 call wrtout(iout,msg) 3110 call wrtout(std_out,msg) 3111 end if 3112 3113 end subroutine dfpt_prtene
ABINIT/eigen_meandege [ Functions ]
NAME
eigen_meandege
FUNCTION
This routine takes the mean values of the responses for the eigenstates that are degenerate in energy.
INPUTS
eigenresp((3-option)*mband**(3-option)*nkpt*nsppol)= input eigenresp
eigenrep(2*mband**2*nkpt*nsppol) for first-order derivatives of the eigenvalues
eigenrep(mband*nkpt*nsppol) for Fan or Debye-Waller second-order derivatives of the eigenvalues
mband= maximum number of bands
natom= number of atoms in the unit cell
nkpt= number of k-points
nsppol= 1 for unpolarized, 2 for spin-polarized
option= 1 for eigen(1), 2 for eigen(2) - Fan or Debye-Waller
OUTPUT
eigenresp_mean(mband*nkpt*nsppol)= eigenresp, averaged over degenerate states
SOURCE
3139 subroutine eigen_meandege(eigen0,eigenresp,eigenresp_mean,mband,nband,nkpt,nsppol,option) 3140 3141 !Arguments ------------------------------------ 3142 !scalars 3143 integer,intent(in) :: mband,nkpt,nsppol,option 3144 integer,intent(in) :: nband(nkpt*nsppol) 3145 3146 !arrays 3147 real(dp),intent(in) :: eigen0(mband*nkpt*nsppol) 3148 real(dp),intent(in) :: eigenresp((3-option)*mband**(3-option)*nkpt*nsppol) 3149 real(dp),intent(out) :: eigenresp_mean(mband*nkpt*nsppol) 3150 3151 !Local variables------------------------------- 3152 !scalars 3153 integer :: bdtot_index,bd2tot_index,iband,ii,ikpt,isppol,nband_k 3154 real(dp) :: eig0,mean 3155 character(len=500) :: msg 3156 !arrays 3157 3158 ! ********************************************************************* 3159 3160 if(option/=1 .and. option/=2)then 3161 write(msg, '(a,i0)' )' The argument option should be 1 or 2, while it is found that option=',option 3162 ABI_BUG(msg) 3163 end if 3164 3165 bdtot_index=0 ; bd2tot_index=0 3166 do isppol=1,nsppol 3167 do ikpt=1,nkpt 3168 nband_k=nband(ikpt+(isppol-1)*nkpt) 3169 if(option==1)then 3170 do iband=1,nband_k 3171 eigenresp_mean(iband+bdtot_index)=& 3172 & eigenresp(2*iband-1 + (iband-1)*2*nband_k + bd2tot_index) 3173 end do 3174 else if(option==2)then 3175 do iband=1,nband_k 3176 eigenresp_mean(iband+bdtot_index)=eigenresp(iband+bdtot_index) 3177 end do 3178 end if 3179 3180 ! Treat the case of degeneracies : take the mean of degenerate states 3181 if(nband_k>1)then 3182 eig0=eigen0(1+bdtot_index) 3183 ii=1 3184 do iband=2,nband_k 3185 if(eigen0(iband+bdtot_index)-eig0<tol8)then 3186 ii=ii+1 3187 else 3188 mean=sum(eigenresp_mean(iband-ii+bdtot_index:iband-1+bdtot_index))/ii 3189 eigenresp_mean(iband-ii+bdtot_index:iband-1+bdtot_index)=mean 3190 ii=1 3191 end if 3192 eig0=eigen0(iband+bdtot_index) 3193 if(iband==nband_k)then 3194 mean=sum(eigenresp_mean(iband-ii+1+bdtot_index:iband+bdtot_index))/ii 3195 eigenresp_mean(iband-ii+1+bdtot_index:iband+bdtot_index)=mean 3196 end if 3197 end do 3198 end if 3199 3200 bdtot_index=bdtot_index+nband_k 3201 bd2tot_index=bd2tot_index+2*nband_k**2 3202 end do 3203 end do 3204 3205 end subroutine eigen_meandege
ABINIT/getcgqphase [ Functions ]
NAME
getcgqphase
FUNCTION
extract phases from wave functions, to cancel contributions to gkk matrix elements
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset timrev = flag for use of time reversal symmetry cg = input wavefunctions mcg = dimension of cg = nspinor*mband*mpw*mkmem cgq = input wavefunctions at k+q mcgq = dimension of cgq = nspinor*mband*mpw*mkmem mpi_enreg = datastructure for mpi communication nkpt_rbz = number of k-points in reduced zone for present q point npwarr = array of numbers of plane waves for each k-point npwar1 = array of numbers of plane waves for each k+q point
OUTPUT
phasecg = phase of different wavefunction products <k,n | k+q,n'>
SOURCE
2608 subroutine getcgqphase(dtset, timrev, cg, mcg, cgq, mcgq, mpi_enreg, nkpt_rbz, npwarr, npwar1, phasecg) 2609 2610 !Arguments ------------------------------- 2611 ! scalars 2612 integer, intent(in) :: mcg, mcgq, timrev 2613 integer, intent(in) :: nkpt_rbz 2614 type(dataset_type), intent(in) :: dtset 2615 2616 ! arrays 2617 integer, intent(in) :: npwarr(nkpt_rbz) 2618 integer, intent(in) :: npwar1(nkpt_rbz) 2619 real(dp), intent(in) :: cg(2,mcg) 2620 real(dp), intent(in) :: cgq(2,mcgq) 2621 type(MPI_type), intent(in) :: mpi_enreg 2622 real(dp),intent(out) :: phasecg(2, dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol) 2623 2624 !Local variables ------------------------- 2625 ! local vars 2626 integer :: icg, icgq, isppol, ikpt, ipw 2627 integer :: istate, iband1, iband2 2628 integer :: npw_k, npw_q 2629 integer :: me, ierr, master, spaceComm, nprocs 2630 integer :: usepaw 2631 integer :: ddkflag, itrs, job, maxbd, mcg1_k, minbd, shiftbd 2632 real(dp) :: normsmat 2633 2634 integer, allocatable :: sflag_k(:) 2635 integer, allocatable :: pwind_k(:) 2636 2637 real(dp) :: cg1_dummy(1,1) 2638 real(dp) :: smat_inv_dummy(1,1,1) 2639 real(dp) :: smat_k_paw_dummy(1,1,1) 2640 real(dp) :: dtm_k_dummy(2) 2641 real(dp), allocatable :: smat_k(:,:,:) 2642 real(dp), allocatable :: pwnsfac_k(:,:) 2643 logical, allocatable :: my_kpt(:,:) 2644 !character(len=500) :: msg 2645 2646 ! ********************************************************************* 2647 2648 ABI_MALLOC(smat_k,(2,dtset%mband,dtset%mband)) 2649 ABI_MALLOC(sflag_k,(dtset%mband)) 2650 2651 !dummy use of timrev so abirules stops complaining. 2652 icg = timrev 2653 2654 !!MPI data for future use 2655 spaceComm=mpi_enreg%comm_cell 2656 nprocs=xmpi_comm_size(spaceComm) 2657 master=0 2658 me=mpi_enreg%me_kpt 2659 2660 ABI_MALLOC(my_kpt, (nkpt_rbz, dtset%nsppol)) 2661 my_kpt = .true. 2662 if (mpi_enreg%nproc_spkpt > 1) then 2663 do isppol = 1, dtset%nsppol 2664 do ikpt = 1, nkpt_rbz 2665 my_kpt(ikpt, isppol) = .not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,& 2666 & dtset%nband(ikpt),isppol,me)) 2667 end do 2668 end do 2669 end if 2670 2671 2672 !make trivial association of G vectors: we just want <psi_k| psi_k+q> 2673 !TODO: check this is correct wrt arrangement of kg vectors for k+q 2674 !looks ok : usually made in initberry, from the translations associated 2675 !to the symops, scalar product with the G vectors. The symop is the one 2676 !used to go from the irreducible k to the full zone k. In present context 2677 !we should be using only the reduced zone, and anyhow have the same k-grid 2678 !for the gkk matrix elements and for the cg here... 2679 ABI_MALLOC(pwind_k,(dtset%mpw)) 2680 ABI_MALLOC(pwnsfac_k,(4,dtset%mpw)) 2681 do ipw = 1, dtset%mpw 2682 pwind_k(ipw) = ipw 2683 pwnsfac_k(1,ipw) = one 2684 pwnsfac_k(2,ipw) = zero 2685 pwnsfac_k(3,ipw) = one 2686 pwnsfac_k(4,ipw) = zero 2687 end do 2688 2689 !flags for call to smatrix 2690 usepaw = 0 ! for now 2691 ddkflag = 0 2692 itrs = 0 2693 job = 0 2694 maxbd = 1 2695 mcg1_k = 1 2696 minbd = 1 2697 shiftbd = 1 2698 2699 !from overlap matrix for each wavefunction, extract phase 2700 icg = 0 2701 icgq = 0 2702 istate = 0 2703 2704 phasecg = zero 2705 do isppol = 1, dtset%nsppol 2706 do ikpt = 1, nkpt_rbz 2707 ! each proc only has certain k 2708 if (.not. my_kpt(ikpt, isppol)) then 2709 istate = istate + dtset%nband(ikpt)*dtset%nband(ikpt) 2710 cycle 2711 end if 2712 2713 npw_k = npwarr(ikpt) 2714 npw_q= npwar1(ikpt) 2715 2716 ! TODO: question: are the k-points in the ibz correctly ordered in cg and cgq? if not the icg below have to be adapted. 2717 sflag_k = 0 ! make sure all elements are calculated 2718 smat_k = zero 2719 2720 call smatrix(cg, cgq, cg1_dummy, ddkflag, dtm_k_dummy, icg, icgq,& 2721 & itrs, job, maxbd, mcg, mcgq, mcg1_k, minbd,dtset%mpw, dtset%mband, dtset%mband,& 2722 & npw_k, npw_q, dtset%nspinor, pwind_k, pwnsfac_k, sflag_k, shiftbd,& 2723 & smat_inv_dummy, smat_k, smat_k_paw_dummy, usepaw) 2724 2725 icg = icg + npw_k*dtset%nspinor*dtset%nband(ikpt) 2726 icgq = icgq + npw_q*dtset%nspinor*dtset%nband(ikpt) 2727 2728 do iband1 = 1, dtset%nband(ikpt) 2729 do iband2 = 1, dtset%nband(ikpt) 2730 istate = istate + 1 2731 ! normalise the overlap matrix element to get just the phase difference phi_k - phi_k+q 2732 normsmat = sqrt(smat_k(1,iband2, iband1)**2 & 2733 & + smat_k(2,iband2, iband1)**2) 2734 if (normsmat > tol12) then 2735 phasecg(:,istate) = smat_k(:,iband2, iband1) / normsmat 2736 ! NOTE: 21/9/2011 these appear to be always 1, i, or -i, to within 1.e-5 at worst! 2737 end if 2738 end do 2739 end do 2740 end do 2741 end do 2742 2743 !eventually do an mpi allreduce over the k-points for phasecg 2744 if (nprocs>1) then 2745 call xmpi_barrier(spaceComm) 2746 call xmpi_sum_master(phasecg,master,spaceComm,ierr) 2747 call xmpi_barrier(spaceComm) 2748 if (1==1) then 2749 call wrtout(std_out, 'In getcgqphase - contributions to phasecg collected') 2750 end if 2751 end if 2752 2753 ABI_FREE(sflag_k) 2754 ABI_FREE(smat_k) 2755 ABI_FREE(pwind_k) 2756 ABI_FREE(pwnsfac_k) 2757 ABI_FREE(my_kpt) 2758 2759 end subroutine getcgqphase
ABINIT/m_dfpt_loopert [ Modules ]
NAME
m_dfpt_loopert
FUNCTION
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (XG, DRH, MB, XW, MT, SPr, MJV) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_dfpt_loopert 23 24 use defs_basis 25 use m_dtset 26 use m_dtfil 27 use defs_wvltypes 28 use m_efmas_defs 29 use m_abicore 30 use m_xmpi 31 use m_errors 32 use m_wfk 33 use m_wffile 34 use m_io_redirect 35 use m_paral_pert 36 use m_nctk 37 use m_ddb 38 use m_wfd 39 use m_ddb_hdr 40 #ifdef HAVE_NETCDF 41 use netcdf 42 #endif 43 use m_hdr 44 use m_ebands 45 46 use defs_datatypes, only : pseudopotential_type, ebands_t 47 use defs_abitypes, only : MPI_type 48 use m_occ, only : getnel 49 use m_io_tools, only : file_exists 50 use m_time, only : timab 51 use m_fstrings, only : strcat, sjoin, ftoa 52 use m_geometry, only : mkrdim, metric, littlegroup_pert 53 use m_exit, only : exit_check, disable_timelimit 54 use m_atomdata, only : atom_gauss 55 use m_eig2d, only : eigr2d_init,eigr2d_t, eigr2d_ncwrite,eigr2d_free, & 56 gkk_t, gkk_init, gkk_ncwrite,gkk_free, outbsd, eig2stern 57 use m_crystal, only : crystal_init, crystal_t 58 use m_efmas, only : efmas_main, efmas_analysis, print_efmas 59 use m_fft, only : fourdp 60 use m_fftcore, only : fftcore_set_mixprec 61 use m_kg, only : getcut, getmpw, kpgio, getph 62 use m_iowf, only : outwf, outresid 63 use m_ioarr, only : read_rhor 64 use m_orbmag, only : orbmag 65 use m_pawang, only : pawang_type, pawang_init, pawang_free 66 use m_pawrad, only : pawrad_type 67 use m_pawtab, only : pawtab_type 68 use m_paw_an, only : paw_an_type 69 use m_paw_ij, only : paw_ij_type 70 use m_pawfgrtab, only : pawfgrtab_type 71 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_bcast, pawrhoij_copy, & 72 pawrhoij_nullify, pawrhoij_redistribute, pawrhoij_inquire_dim 73 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, pawcprj_getdim , pawcprj_output 74 use m_pawfgr, only : pawfgr_type 75 use m_paw_sphharm,only : setsym_ylm 76 use m_rf2, only : rf2_getidirs 77 use m_iogkk, only : outgkk 78 use m_inwffil, only : inwffil 79 use m_spacepar, only : rotate_rho, setsym 80 use m_initylmg, only : initylmg 81 use m_dfpt_scfcv, only : dfpt_scfcv 82 use m_dfpt_mkrho, only : dfpt_mkrho 83 use m_mpinfo, only : initmpi_band, distrb2, proc_distrb_cycle 84 use m_atm2fft, only : dfpt_atm2fft 85 use m_berrytk, only : smatrix 86 use m_common, only : prteigrs 87 use m_fourier_interpol, only : transgrid 88 use m_mkcore, only : dfpt_mkcore 89 use m_mklocl, only : dfpt_vlocal, vlocalstr 90 use m_cgprj, only : ctocprj 91 use m_symkpt, only : symkpt 92 93 implicit none 94 95 private