TABLE OF CONTENTS
ABINIT/dfpt_looppert [ Functions ]
NAME
dfpt_looppert
FUNCTION
Loop over perturbations
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (XG, DRH, MB, XW, MT,SPr) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) 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)=boradening 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=informations 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 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)
PARENTS
respfn
CHILDREN
appdig,atom_gauss,crystal_free,crystal_init,ctocprj,ddb_free ddb_from_file,ddb_hdr_free,ddb_hdr_init,ddb_hdr_open_write,dfpt_atm2fft dfpt_init_mag1,dfpt_mkcore,dfpt_mkrho,dfpt_prtene,dfpt_scfcv dfpt_vlocal,disable_timelimit,distrb2,dtset_copy,dtset_free,ebands_free ebands_init,efmas_main,eig2stern,eigen_meandege,eigr2d_free,eigr2d_init eigr2d_ncwrite,exit_check,fourdp,getcgqphase,getcut,getmpw,getnel,getph gkk_free,gkk_init,gkk_ncwrite,hdr_copy,hdr_free,hdr_init,hdr_update initmpi_band,initylmg,inwffil,kpgio,littlegroup_pert,localfilnam localrdfile,localredirect,localwrfile,metric,mkrdim,outbsd,outgkk,outwf pawang_free,pawang_init,pawcprj_alloc,pawcprj_copy,pawcprj_free pawcprj_getdim,pawrhoij_alloc,pawrhoij_copy,pawrhoij_free pawrhoij_nullify,prteigrs,put_eneocc_vect,read_rhor,rf2_getidirs rotate_rho,set_pert_comm,set_pert_paw,setsym,setsymrhoij,status,symkpt timab,transgrid,unset_pert_comm,unset_pert_paw,vlocalstr,wffclose wfk_open_read,wfk_read_eigenvalues,wrtout,xmpi_sum
SOURCE
116 #if defined HAVE_CONFIG_H 117 #include "config.h" 118 #endif 119 120 #include "abi_common.h" 121 122 subroutine dfpt_looppert(atindx,blkflg,codvsn,cpus,dim_eigbrd,dim_eig2nkq,doccde,& 123 & ddkfil,dtfil,dtset,dyew,dyfrlo,dyfrnl,dyfrx1,dyfrx2,dyvdw,& 124 & dyfr_cplex,dyfr_nondiag,d2bbb,d2lo,d2nl,d2ovl,efmasdeg,efmasfr,eigbrd,eig2nkq,& 125 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,& 126 & etotal,fermie,iexit,indsym,kxc,& 127 & mkmem,mkqmem,mk1mem,mpert,mpi_enreg,my_natom,nattyp,& 128 & nfftf,nhat,nkpt,nkxc,nspden,nsym,occ,& 129 & paw_an,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,& 130 & pertsy,prtbbb,psps,rfpert,rf2_dirs_from_rfpert_nl,rhog,rhor,symq,symrec,timrev,& 131 & usecprj,usevdw,vtrial,vxc,vxcavg,xred,clflg,occ_rbz_pert,eigen0_pert,eigenq_pert,& 132 & eigen1_pert,nkpt_rbz,eigenq_fine,hdr_fine,hdr0) 133 134 use defs_basis 135 use defs_datatypes 136 use defs_abitypes 137 use defs_wvltypes 138 use m_efmas_defs 139 use m_profiling_abi 140 use m_xmpi 141 use m_errors 142 use m_wfk 143 use m_wffile 144 use m_io_redirect 145 use m_paral_pert 146 use m_abi_etsf 147 use m_nctk 148 use m_ddb 149 #ifdef HAVE_NETCDF 150 use netcdf 151 #endif 152 use m_hdr 153 use m_ebands 154 155 use m_occ, only : getnel 156 use m_ddb_hdr, only : ddb_hdr_type, ddb_hdr_init, ddb_hdr_free, ddb_hdr_open_write 157 use m_io_tools, only : file_exists 158 use m_fstrings, only : strcat 159 use m_geometry, only : mkrdim, metric 160 use m_exit, only : exit_check, disable_timelimit 161 use m_atomdata, only : atom_gauss 162 use m_eig2d, only : eigr2d_init,eigr2d_t, eigr2d_ncwrite,eigr2d_free, & 163 & gkk_t, gkk_init, gkk_ncwrite,gkk_free 164 use m_crystal, only : crystal_init, crystal_free, crystal_t 165 use m_crystal_io, only : crystal_ncwrite 166 use m_efmas, only : efmas_main 167 use m_kg, only : getcut, getmpw, kpgio, getph 168 use m_dtset, only : dtset_copy, dtset_free 169 use m_iowf, only : outwf 170 use m_ioarr, only : read_rhor 171 use m_pawang, only : pawang_type, pawang_init, pawang_free 172 use m_pawrad, only : pawrad_type 173 use m_pawtab, only : pawtab_type 174 use m_paw_an, only : paw_an_type 175 use m_paw_ij, only : paw_ij_type 176 use m_pawfgrtab, only : pawfgrtab_type 177 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_bcast, pawrhoij_copy, & 178 & pawrhoij_nullify, pawrhoij_redistribute, pawrhoij_get_nspden 179 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, pawcprj_getdim 180 use m_pawfgr, only : pawfgr_type 181 use m_rf2, only : rf2_getidirs 182 183 !This section has been created automatically by the script Abilint (TD). 184 !Do not modify the following lines by hand. 185 #undef ABI_FUNC 186 #define ABI_FUNC 'dfpt_looppert' 187 use interfaces_14_hidewrite 188 use interfaces_18_timing 189 use interfaces_32_util 190 use interfaces_41_geometry 191 use interfaces_51_manage_mpi 192 use interfaces_53_ffts 193 use interfaces_56_recipspace 194 use interfaces_64_psp 195 use interfaces_65_paw 196 use interfaces_66_nonlocal 197 use interfaces_67_common 198 use interfaces_72_response 199 use interfaces_79_seqpar_mpi 200 !End of the abilint section 201 202 implicit none 203 204 !Arguments ------------------------------------ 205 integer, intent(in) :: dim_eigbrd,dim_eig2nkq,dyfr_cplex,dyfr_nondiag,mk1mem,mkmem,mkqmem,mpert 206 integer, intent(in) :: nfftf,nkpt,nkxc,nspden,nsym,prtbbb,timrev,usecprj,usevdw 207 integer, intent(out) :: iexit 208 integer, intent(inout) :: my_natom 209 real(dp), intent(in) :: cpus,vxcavg 210 real(dp), intent(inout) :: fermie 211 real(dp), intent(inout) :: etotal 212 character(len=6), intent(in) :: codvsn 213 type(MPI_type), intent(inout) :: mpi_enreg 214 type(datafiles_type), intent(in) :: dtfil 215 type(dataset_type), intent(in), target :: dtset 216 type(pawang_type),intent(in) :: pawang 217 type(pawfgr_type),intent(in) :: pawfgr 218 type(pseudopotential_type), intent(inout) :: psps 219 integer, intent(in) :: atindx(dtset%natom),indsym(4,nsym,dtset%natom) 220 integer, intent(in) :: nattyp(dtset%ntypat),pertsy(3,dtset%natom+6) 221 integer, intent(in) :: rfpert(mpert),rf2_dirs_from_rfpert_nl(3,3),symq(4,2,nsym),symrec(3,3,nsym) 222 integer, intent(out) :: ddkfil(3) 223 integer, intent(inout) :: blkflg(3,mpert,3,mpert) 224 integer, intent(out) :: clflg(3,mpert) 225 real(dp), intent(in) :: doccde(dtset%mband*nkpt*dtset%nsppol) 226 real(dp), intent(in) :: dyew(2,3,dtset%natom,3,dtset%natom) 227 real(dp), intent(in) :: dyfrlo(3,3,dtset%natom) 228 real(dp), intent(in) :: dyfrnl(dyfr_cplex,3,3,dtset%natom,1+(dtset%natom-1)*dyfr_nondiag) 229 real(dp), intent(in) :: dyfrx1(2,3,dtset%natom,3,dtset%natom) 230 real(dp), intent(in) :: dyfrx2(3,3,dtset%natom),dyvdw(2,3,dtset%natom,3,dtset%natom*usevdw) 231 real(dp), intent(in) :: eltcore(6,6),elteew(6+3*dtset%natom,6),eltfrhar(6,6) 232 real(dp), intent(in) :: eltfrkin(6,6),eltfrloc(6+3*dtset%natom,6) 233 real(dp), intent(in) :: eltfrnl(6+3*dtset%natom,6) 234 real(dp), intent(in) :: eltfrxc(6+3*dtset%natom,6),eltvdw(6+3*dtset%natom,6*usevdw) 235 real(dp), intent(in) :: kxc(nfftf,nkxc),nhat(nfftf,nspden) 236 real(dp), intent(in) :: occ(dtset%mband*nkpt*dtset%nsppol) 237 real(dp), intent(in) :: rhog(2,nfftf),rhor(nfftf,nspden),vxc(nfftf,nspden) 238 real(dp), intent(in) :: vtrial(nfftf,nspden) 239 real(dp), intent(inout) :: xred(3,dtset%natom) 240 real(dp), intent(inout) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)!vz_i 241 real(dp), intent(inout) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert) !vz_i 242 real(dp), intent(inout) :: d2ovl(2,3,mpert,3,mpert*psps%usepaw) !vz_i 243 real(dp), intent(out) :: eigbrd(2,dtset%mband*dtset%nsppol,nkpt,3,dtset%natom,3,dtset%natom*dim_eigbrd) 244 real(dp), intent(out) :: eig2nkq(2,dtset%mband*dtset%nsppol,nkpt,3,dtset%natom,3,dtset%natom*dim_eig2nkq) 245 type(efmasdeg_type),allocatable,intent(in) :: efmasdeg(:) 246 type(efmasfr_type),allocatable,intent(in) :: efmasfr(:,:) 247 type(paw_an_type),allocatable,target,intent(inout) :: paw_an(:) 248 type(paw_ij_type),allocatable,target,intent(inout) :: paw_ij(:) 249 type(pawfgrtab_type),allocatable,target,intent(inout) :: pawfgrtab(:) 250 type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw) 251 type(pawrhoij_type),allocatable,target,intent(inout) :: pawrhoij(:) 252 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 253 real(dp),pointer :: eigen1_pert(:,:,:) 254 real(dp),intent(out) :: occ_rbz_pert(:),eigen0_pert(:),eigenq_pert(:) 255 real(dp),pointer :: eigenq_fine(:,:,:) 256 integer, intent(out) :: nkpt_rbz 257 type(hdr_type),intent(out) :: hdr0,hdr_fine 258 259 !Local variables------------------------------- 260 !scalars 261 integer,parameter :: level=11,response=1,formeig1=1,master=0,fake_unit=-666 262 integer :: ask_accurate,band_index,bantot,bantot_rbz,bdeigrf,bdtot1_index,nsppol,nspinor,band2tot_index 263 integer :: bdtot_index,choice,cplex,cplex_rhoij,dim_eig2rf,formeig 264 integer :: gscase,iband,iblok,icase,icase_eq,idir,idir0,idir1,idir2,idir_eq,idir_dkdk,ierr 265 integer :: ii,ikpt,ikpt1,jband,initialized,iorder_cprj,ipert,ipert_cnt,ipert_eq,ipert_me,ireadwf0 266 integer :: iscf_mod,iscf_mod_save,isppol,istr,isym,mcg,mcgq,mcg1,mcprj,mcprjq,mband 267 integer :: mcgmq,mcg1mq,mpw1_mq !+/-q duplicates 268 integer :: maxidir,me,mgfftf,mkmem_rbz,mk1mem_rbz,mkqmem_rbz,mpw,mpw1,my_nkpt_rbz 269 integer :: n3xccc,nband_k,ncpgr,ndir,nkpt_eff,nkpt_max,nline_save,nmatel,npert_io,npert_me,nspden_rhoij 270 integer :: nstep_save,nsym1,ntypat,nwffile,nylmgr,nylmgr1,old_comm_atom,openexit,option,optorth,optthm,pertcase 271 integer :: rdwr,rdwrpaw,spaceComm,smdelta,timrev_pert,timrev_kpt,to_compute_this_pert 272 integer :: unitout,useylmgr,useylmgr1,vrsddb,dfpt_scfcv_retcode,optn2 273 #ifdef HAVE_NETCDF 274 integer :: ncerr,ncid 275 #endif 276 real(dp) :: boxcut,dosdeltae,eberry,ecore,ecut_eff,ecutf,edocc,eei,eeig0,eew,efrhar,efrkin,efrloc 277 real(dp) :: efrnl,efrx1,efrx2,ehart,ehart01,ehart1,eii,ek,ek0,ek1,ek2,eloc0 278 real(dp) :: elpsp1,enl,enl0,enl1,entropy,enxc,eovl1,epaw1,evdw,exc1,fsum,gsqcut,maxocc,nelectkq 279 real(dp) :: residm,tolwfr,tolwfr_save,toldfe_save,toldff_save,tolrff_save,tolvrs_save 280 real(dp) :: ucvol, eig1_r, eig1_i 281 real(dp) :: residm_mq !+/-q duplicates 282 logical,parameter :: paral_pert_inplace=.true.,remove_inv=.false. 283 logical :: first_entry,found_eq_gkk,t_exist,paral_atom,write_1wfk,init_rhor1 284 logical :: kramers_deg 285 character(len=fnlen) :: dscrpt,fiden1i,fiwf1i,fiwf1o,fiwfddk,fnamewff(4),gkkfilnam,fname,filnam 286 character(len=500) :: message 287 type(crystal_t) :: crystal,ddb_crystal 288 type(dataset_type), pointer :: dtset_tmp 289 type(ebands_t) :: ebands_k,ebands_kq,gkk_ebands 290 type(ebands_t) :: ebands_kmq !+/-q duplicates 291 type(eigr2d_t) :: eigr2d,eigi2d 292 type(gkk_t) :: gkk2d 293 type(hdr_type) :: hdr,hdr_den,hdr_tmp 294 type(ddb_hdr_type) :: ddb_hdr 295 type(pawang_type) :: pawang1 296 type(wffile_type) :: wff1,wffgs,wffkq,wffnow,wfftgs,wfftkq 297 type(wfk_t) :: ddk_f(4) 298 type(wvl_data) :: wvl 299 !arrays 300 integer :: eq_symop(3,3),ngfftf(18),file_index(4),rfdir(9),rf2dir(9),rf2_dir1(3),rf2_dir2(3) 301 integer,allocatable :: blkflg_save(:,:,:,:),dimcprj_srt(:),dummy(:),dyn(:),indkpt1(:),indkpt1_tmp(:) 302 integer,allocatable :: indsy1(:,:,:),irrzon1(:,:,:),istwfk_rbz(:),istwfk_pert(:,:,:) 303 integer,allocatable :: kg(:,:),kg1(:,:),nband_rbz(:),npwar1(:),npwarr(:),npwtot(:) 304 integer,allocatable :: kg1_mq(:,:),npwar1_mq(:),npwtot1_mq(:) !+q/-q duplicates 305 integer,allocatable :: npwtot1(:),npwar1_pert(:,:),npwarr_pert(:,:),npwtot_pert(:,:) 306 integer,allocatable :: pert_calc(:,:),pert_tmp(:,:) 307 integer,allocatable :: symaf1(:),symaf1_tmp(:),symrc1(:,:,:),symrl1(:,:,:),symrl1_tmp(:,:,:) 308 integer, pointer :: old_atmtab(:) 309 real(dp) :: dielt(3,3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),tsec(2) 310 real(dp),allocatable :: buffer1(:,:,:,:,:),cg(:,:),cg1(:,:),cg1_active(:,:),cg0_pert(:,:) 311 real(dp),allocatable :: cg1_pert(:,:,:,:),cgq(:,:),gh0c1_pert(:,:,:,:) 312 real(dp),allocatable :: doccde_rbz(:),docckqde(:) 313 real(dp),allocatable :: gh1c_pert(:,:,:,:),eigen0(:),eigen0_copy(:),eigen1(:),eigen1_mean(:) 314 real(dp),allocatable :: eigenq(:),gh1c_set(:,:),gh0c1_set(:,:),kpq(:,:) 315 real(dp),allocatable :: kpq_rbz(:,:),kpt_rbz(:,:),occ_pert(:),occ_rbz(:),occkq(:),kpt_rbz_pert(:,:) 316 real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phnons1(:,:,:),resid(:),rhog1(:,:) 317 real(dp),allocatable :: rhor1_save(:,:,:) 318 real(dp),allocatable :: rhor1(:,:),rho1wfg(:,:),rho1wfr(:,:),tnons1(:,:),tnons1_tmp(:,:) 319 real(dp),allocatable :: rhor1_pq(:,:),rhor1_mq(:,:),rhog1_pq(:,:),rhog1_mq(:,:) !+q/-q duplicates 320 real(dp),allocatable :: cg_mq(:,:),cg1_mq(:,:),resid_mq(:) ! 321 real(dp),allocatable :: cg1_active_mq(:,:),occk_mq(:) ! 322 real(dp),allocatable :: kmq(:,:),kmq_rbz(:,:),gh0c1_set_mq(:,:) ! 323 real(dp),allocatable :: eigen_mq(:),gh1c_set_mq(:,:),docckde_mq(:),eigen1_mq(:) ! 324 real(dp),allocatable :: vpsp1(:),work(:),wtk_folded(:),wtk_rbz(:),xccc3d1(:) 325 real(dp),allocatable :: ylm(:,:),ylm1(:,:),ylmgr(:,:,:),ylmgr1(:,:,:),zeff(:,:,:) 326 real(dp),allocatable :: phasecg(:,:),gauss(:,:) 327 real(dp),allocatable :: gkk(:,:,:,:,:) 328 type(pawcprj_type),allocatable :: cprj(:,:),cprjq(:,:) 329 type(paw_ij_type),pointer :: paw_ij_pert(:) 330 type(paw_an_type),pointer :: paw_an_pert(:) 331 type(pawfgrtab_type),pointer :: pawfgrtab_pert(:) 332 type(pawrhoij_type),allocatable :: pawrhoij1(:) 333 type(pawrhoij_type),pointer :: pawrhoij_pert(:) 334 type(ddb_type) :: ddb 335 336 ! *********************************************************************** 337 338 DBG_ENTER("COLL") 339 340 _IBM6("In dfpt_looppert") 341 342 call timab(141,1,tsec) 343 344 call status(0,dtfil%filstat,iexit,level,'init ') 345 346 !Structured debugging if prtvol==-level 347 if(dtset%prtvol==-level)then 348 write(message,'(80a,a,a)') ('=',ii=1,80),ch10,' dfpt_looppert : enter , debug mode ' 349 call wrtout(std_out,message,'COLL') 350 end if 351 352 dfpt_scfcv_retcode = -1 353 nsppol = dtset%nsppol; nspinor = dtset%nspinor 354 355 kramers_deg=.true. 356 if (dtset%tim1rev==0) then 357 kramers_deg=.false. 358 end if 359 360 !Obtain dimensional translations in reciprocal space gprimd, 361 !metrics and unit cell volume, from rprimd. Also output rprimd, gprimd and ucvol 362 call mkrdim(dtset%acell_orig(1:3,1),dtset%rprim_orig(1:3,1:3,1),rprimd) 363 call metric(gmet,gprimd,std_out,rmet,rprimd,ucvol) 364 365 call crystal_init(dtset%amu_orig(:,1),crystal,dtset%spgroup,dtset%natom,dtset%npsp,& 366 & psps%ntypat,dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,& 367 & dtset%nspden==2.and.dtset%nsppol==1,remove_inv,psps%title,& 368 & symrel=dtset%symrel,tnons=dtset%tnons,symafm=dtset%symafm) 369 370 !Get FFT grid(s) sizes (be careful !) See NOTES in the comments at the beginning of respfn.F90 371 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then 372 mgfftf=pawfgr%mgfft;ngfftf(:)=pawfgr%ngfft(:) 373 ecutf=dtset%pawecutdg 374 else 375 mgfftf=dtset%mgfft;ngfftf(:)=dtset%ngfft(:) 376 ecutf=dtset%ecut 377 end if 378 ecut_eff=dtset%ecut*(dtset%dilatmx)**2 379 380 !Compute large sphere cut-off gsqcut 381 if (psps%usepaw==1) then 382 call wrtout(std_out,ch10//' FFT (fine) grid used for densities/potentials:','COLL') 383 end if 384 call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,dtset%qptn,ngfftf) 385 386 !Various initializations/allocations 387 iscf_mod=dtset%iscf 388 ntypat=psps%ntypat 389 nkpt_max=50;if (xmpi_paral==1) nkpt_max=-1 390 paral_atom=(dtset%natom/=my_natom) 391 cplex=2-timrev !cplex=2 ! DEBUG: impose cplex=2 392 first_entry=.true. 393 initialized=0 394 ecore=zero ; ek=zero ; ehart=zero ; enxc=zero ; eei=zero ; enl=zero ; eii=zero 395 clflg(:,:)=0 ! Array on calculated perturbations for eig2rf 396 if (psps%usepaw==1) then 397 ABI_ALLOCATE(dimcprj_srt,(dtset%natom)) 398 call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 399 end if 400 401 !Save values of SCF cycle parameters 402 iscf_mod_save = iscf_mod 403 nstep_save = dtset%nstep 404 nline_save = dtset%nline 405 tolwfr_save = dtset%tolwfr 406 toldfe_save = dtset%toldfe 407 toldff_save = dtset%toldff 408 tolrff_save = dtset%tolrff 409 tolvrs_save = dtset%tolvrs 410 411 !This dtset will be used in dfpt_scfcv to force non scf calculations for equivalent perturbations 412 nullify(dtset_tmp) 413 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 414 ABI_DATATYPE_ALLOCATE(dtset_tmp,) 415 call dtset_copy(dtset_tmp, dtset) 416 else 417 dtset_tmp => dtset 418 end if 419 420 !Generate the 1-dimensional phases 421 ABI_ALLOCATE(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom)) 422 ABI_ALLOCATE(ph1df,(2,3*(2*mgfftf+1)*dtset%natom)) 423 call getph(atindx,dtset%natom,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),ph1d,xred) 424 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then 425 call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred) 426 else 427 ph1df(:,:)=ph1d(:,:) 428 end if 429 430 !!Determine existence of pertubations and of pertubation symmetries 431 !!Create array with pertubations which have to be calculated 432 ! ABI_ALLOCATE(pert_tmp,(3*mpert)) 433 ! ipert_cnt=0 434 ! do ipert=1,mpert 435 ! do idir=1,3 436 ! if( rfpert(ipert)==1 .and. dtset%rfdir(idir) == 1 )then 437 ! if ((pertsy(idir,ipert)==1).or.& 438 !& ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.& 439 !& ((dtset%prepgkk == 1).and.(ipert <= dtset%natom)) ) then 440 ! ipert_cnt = ipert_cnt+1; 441 ! pert_tmp(ipert_cnt) = idir+(ipert-1)*3 442 ! else 443 ! write(message, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,& 444 !& ' The perturbation idir=',idir,' ipert=',ipert,' is',ch10,& 445 !& ' symmetric of a previously calculated perturbation.',ch10,& 446 !& ' So, its SCF calculation is not needed.',ch10 447 ! call wrtout(std_out,message,'COLL') 448 ! call wrtout(ab_out,message,'COLL') 449 ! end if ! Test of existence of symmetry of perturbation 450 ! end if ! Test of existence of perturbation 451 ! end do 452 ! end do 453 ! ABI_ALLOCATE(pert_calc,(ipert_cnt)) 454 ! do icase=1,ipert_cnt 455 ! pert_calc(icase)=pert_tmp(icase) 456 ! end do 457 ! ABI_DEALLOCATE(pert_tmp) 458 459 !Initialize rf2dir : 460 rf2_dir1(1:3)=dtset%rf2_pert1_dir(1:3) 461 rf2_dir2(1:3)=dtset%rf2_pert2_dir(1:3) 462 if (sum(rf2_dir1)==3.and.sum(rf2_dir2)==3.and.dtset%prepanl==1) then 463 ! Diagonal terms : 464 rf2dir(1) = rf2_dirs_from_rfpert_nl(1,1) 465 rf2dir(2) = rf2_dirs_from_rfpert_nl(2,2) 466 rf2dir(3) = rf2_dirs_from_rfpert_nl(3,3) 467 ! Upper triangular terms : 468 rf2dir(4) = rf2_dirs_from_rfpert_nl(2,3) 469 rf2dir(5) = rf2_dirs_from_rfpert_nl(1,3) 470 rf2dir(6) = rf2_dirs_from_rfpert_nl(1,2) 471 ! Lower triangular terms : 472 rf2dir(7) = rf2_dirs_from_rfpert_nl(3,2) 473 rf2dir(8) = rf2_dirs_from_rfpert_nl(3,1) 474 rf2dir(9) = rf2_dirs_from_rfpert_nl(2,1) 475 else 476 ! Diagonal terms : 477 rf2dir(1) = rf2_dir1(1)*rf2_dir2(1) 478 rf2dir(2) = rf2_dir1(2)*rf2_dir2(2) 479 rf2dir(3) = rf2_dir1(3)*rf2_dir2(3) 480 ! Upper triangular terms : 481 rf2dir(4) = rf2_dir1(2)*rf2_dir2(3) 482 rf2dir(5) = rf2_dir1(1)*rf2_dir2(3) 483 rf2dir(6) = rf2_dir1(1)*rf2_dir2(2) 484 ! Lower triangular terms : 485 rf2dir(7) = rf2_dir1(3)*rf2_dir2(2) 486 rf2dir(8) = rf2_dir1(3)*rf2_dir2(1) 487 rf2dir(9) = rf2_dir1(2)*rf2_dir2(1) 488 end if 489 490 !Determine existence of pertubations and of pertubation symmetries 491 !Create array with pertubations which have to be calculated 492 ABI_ALLOCATE(pert_tmp,(3,3*(dtset%natom+6)+18)) 493 ipert_cnt=0 494 do ipert=1,mpert 495 if (ipert<dtset%natom+10) then 496 maxidir = 3 497 rfdir(1:3) = dtset%rfdir(:) 498 rfdir(4:9) = 0 499 else 500 maxidir = 9 501 rfdir(1:9) = rf2dir(:) 502 end if 503 do idir=1,maxidir 504 to_compute_this_pert = 0 505 if(ipert<dtset%natom+10 .and. rfpert(ipert)==1 .and. rfdir(idir) == 1 ) then 506 if ((pertsy(idir,ipert)==1).or.& 507 & ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.& 508 & ((dtset%prepgkk == 1).and.(ipert <= dtset%natom)) ) then 509 to_compute_this_pert = 1 510 else 511 write(message, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,& 512 & ' The perturbation idir=',idir,' ipert=',ipert,' is',ch10,& 513 & ' symmetric of a previously calculated perturbation.',ch10,& 514 & ' So, its SCF calculation is not needed.',ch10 515 call wrtout(std_out,message,'COLL') 516 call wrtout(ab_out,message,'COLL') 517 end if ! Test of existence of symmetry of perturbation 518 else if (ipert==dtset%natom+11 .and. rfpert(ipert)==1 .and. rfdir(idir) == 1 ) then 519 to_compute_this_pert = 1 520 else if (ipert==dtset%natom+10 .and. rfpert(ipert)==1 .and. idir <= 6) then 521 if (idir<=3) then 522 if (rfdir(idir) == 1) to_compute_this_pert = 1 523 else 524 if (rfdir(idir) == 1 .or. rfdir(idir+3) == 1) to_compute_this_pert = 1 525 end if 526 end if 527 if (to_compute_this_pert /= 0) then 528 ipert_cnt = ipert_cnt+1; 529 pert_tmp(1,ipert_cnt) = ipert 530 pert_tmp(2,ipert_cnt) = idir 531 ! Store "pertcase" in pert_tmp(3,ipert_cnt) 532 if (ipert<dtset%natom+10) then 533 pert_tmp(3,ipert_cnt) = idir + (ipert-1)*3 534 else 535 pert_tmp(3,ipert_cnt) = idir + (ipert-dtset%natom-10)*9 + (dtset%natom+6)*3 536 end if 537 end if 538 end do ! idir 539 end do !ipert 540 ! do ipert=1,mpert 541 ! if (ipert<dtset%natom+10) then 542 ! maxidir = 3 543 ! rfdir(1:3) = dtset%rfdir(:) 544 ! rfdir(4:9) = 0 545 ! else 546 ! maxidir = 9 547 ! rfdir(1:9) = rf2dir(:) 548 ! end if 549 ! do idir=1,maxidir 550 ! if( rfpert(ipert)==1 .and. rfdir(idir) == 1 )then 551 ! to_compute_this_pert = 0 552 ! if (ipert>=dtset%natom+10) then 553 ! to_compute_this_pert = 1 554 ! else if ((pertsy(idir,ipert)==1).or.& 555 !& ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.& 556 !& ((dtset%prepgkk == 1).and.(ipert <= dtset%natom)) ) then 557 ! to_compute_this_pert = 1 558 ! end if 559 ! if (to_compute_this_pert /= 0) then 560 ! ipert_cnt = ipert_cnt+1; 561 ! pert_tmp(1,ipert_cnt) = ipert 562 ! pert_tmp(2,ipert_cnt) = idir 563 !! Store "pertcase" in pert_tmp(3,ipert_cnt) 564 ! if (ipert<dtset%natom+10) then 565 ! pert_tmp(3,ipert_cnt) = idir + (ipert-1)*3 566 ! else 567 ! pert_tmp(3,ipert_cnt) = idir + (ipert-dtset%natom-10)*9 + (dtset%natom+6)*3 568 ! end if 569 ! else 570 ! write(message, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,& 571 !& ' The perturbation idir=',idir,' ipert=',ipert,' is',ch10,& 572 !& ' symmetric of a previously calculated perturbation.',ch10,& 573 !& ' So, its SCF calculation is not needed.',ch10 574 ! call wrtout(std_out,message,'COLL') 575 ! call wrtout(ab_out,message,'COLL') 576 ! end if ! Test of existence of symmetry of perturbation 577 ! end if ! Test of existence of perturbation 578 ! end do 579 ! end do 580 ABI_ALLOCATE(pert_calc,(3,ipert_cnt)) 581 do icase=1,ipert_cnt 582 pert_calc(:,icase)=pert_tmp(:,icase) 583 end do 584 ABI_DEALLOCATE(pert_tmp) 585 586 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 587 ABI_ALLOCATE(rhor1_save,(cplex*nfftf,nspden,ipert_cnt)) 588 rhor1_save=zero 589 ABI_ALLOCATE(blkflg_save,(3,mpert,3,mpert)) 590 end if 591 592 ! Initialize quantities for netcdf print 593 ABI_ALLOCATE(eigen0_copy,(dtset%mband*nkpt*dtset%nsppol)) 594 eigen0_copy(:)=zero 595 596 ! SP : Retreval of the DDB information and computing of effective charge and 597 ! dielectric tensor 598 ABI_ALLOCATE(zeff,(3,3,dtset%natom)) 599 if (dtset%getddb .ne. 0 .or. dtset%irdddb .ne. 0 ) then 600 filnam = dtfil%filddbsin !'test_DDB' 601 ABI_ALLOCATE(dummy,(dtset%natom)) 602 call ddb_from_file(ddb,filnam,1,dtset%natom,0,dummy,ddb_crystal,mpi_enreg%comm_world) 603 ! Get Dielectric Tensor and Effective Charges 604 ! (initialized to one_3D and zero if the derivatives are not available in 605 ! the DDB file) 606 iblok = ddb_get_dielt_zeff(ddb,ddb_crystal,1,0,0,dielt,zeff) 607 call crystal_free(ddb_crystal) 608 call ddb_free(ddb) 609 ABI_DEALLOCATE(dummy) 610 end if 611 612 !%%%% Parallelization over perturbations %%%%% 613 !*Define file output/log file names 614 npert_io=ipert_cnt;if (dtset%nppert<=1) npert_io=0 615 call localfilnam(mpi_enreg%comm_pert,mpi_enreg%comm_cell_pert,mpi_enreg%comm_world,dtfil%filnam_ds,'_PRT',npert_io) 616 !Compute the number of perturbation done by the current cpu 617 if(mpi_enreg%paral_pert==1) then 618 npert_me = 0 ; ipert_me = 0 619 do icase=1,ipert_cnt 620 if (mpi_enreg%distrb_pert(icase)==mpi_enreg%me_pert) npert_me=npert_me +1 621 end do 622 end if 623 624 !*Redefine communicators 625 call set_pert_comm(mpi_enreg,dtset%nppert) 626 627 !*Redistribute PAW on-site data 628 nullify(old_atmtab,pawfgrtab_pert,pawrhoij_pert,paw_an_pert,paw_ij_pert) 629 if (paral_pert_inplace) then 630 call set_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,& 631 & paw_an,paw_ij,pawfgrtab,pawrhoij) 632 pawfgrtab_pert=>pawfgrtab ; pawrhoij_pert=>pawrhoij 633 paw_an_pert =>paw_an ; paw_ij_pert =>paw_ij 634 635 else 636 call set_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,& 637 & paw_an,paw_ij,pawfgrtab,pawrhoij,& 638 & paw_an_out=paw_an_pert,paw_ij_out=paw_ij_pert,& 639 & pawfgrtab_out=pawfgrtab_pert,pawrhoij_out=pawrhoij_pert) 640 641 end if 642 643 ! We can handle the time limit in dfpt_scfcv in a robust manner only if we have one perturbation. 644 if (ipert_cnt > 1 .or. mpi_enreg%paral_pert == 1) call disable_timelimit() 645 646 !Loop on perturbations 647 !========================================================================== 648 do icase=1,ipert_cnt 649 _IBM6("In loop on perts") 650 651 ! %%%% Parallelization over perturbations %%%%% 652 ! Select the perturbations treated by curent processor 653 if(mpi_enreg%paral_pert==1) then 654 if (mpi_enreg%distrb_pert(icase)/=mpi_enreg%me_pert) cycle 655 end if 656 657 ! *Redefine output/log files 658 call localwrfile(mpi_enreg%comm_cell,icase,npert_io,mpi_enreg%paral_pert,0) 659 660 !! Retrieve type and direction of the perturbation 661 ! if (pert_calc(icase) <= dtset%natom*3) then 662 ! idir = mod(pert_calc(icase),3) 663 ! if (idir==0) idir=3 664 ! ipert=( (pert_calc(icase)-idir) / 3 + 1) 665 ! else if (pert_calc(icase) <= dtset%natom*3+4) then 666 ! ipert = dtset%natom + ((pert_calc(icase) - 3*dtset%natom - 1) / 3) + 1 667 ! idir = mod(pert_calc(icase),3) 668 ! if (idir==0) idir=3 669 ! else 670 ! ipert = dtset%natom + ((pert_calc(icase) - 3*(dtset%natom+4) - 1) / 9) + 1 671 ! end if 672 ! pertcase=idir+(ipert-1)*3 673 674 ! Retrieve type and direction of the perturbation 675 ipert = pert_calc(1,icase) 676 idir = pert_calc(2,icase) 677 istr=idir 678 pertcase = pert_calc(3,icase) 679 680 ! Init MPI communicator 681 spaceComm=mpi_enreg%comm_cell 682 me=mpi_enreg%me_cell 683 684 ! ===== Describe the perturbation in output/log file 685 _IBM6("IBM6 before print perturbation") 686 687 write(message, '(a,80a,a,a,3f10.6)' ) ch10,('-',ii=1,80),ch10,& 688 & ' Perturbation wavevector (in red.coord.) ',dtset%qptn(:) 689 call wrtout(std_out,message,'COLL') 690 call wrtout(ab_out,message,'COLL') 691 if(ipert>=1 .and. ipert<=dtset%natom)then 692 write(message, '(a,i4,a,i4)' )' Perturbation : displacement of atom',ipert,' along direction',idir 693 call wrtout(std_out,message,'COLL') 694 call wrtout(ab_out,message,'COLL') 695 if(iscf_mod == -3)then 696 write(message, '(a,a,a,a,a,a,a,a)' )ch10,& 697 & ' dfpt_looppert : COMMENT -',ch10,& 698 & ' The first-order density is imposed to be zero (iscf=-3).',ch10,& 699 & ' Although this is strange in the case of phonons,',ch10,& 700 & ' you are allowed to do so.' 701 call wrtout(std_out,message,'COLL') 702 call wrtout(ab_out,message,'COLL') 703 end if 704 else if(ipert==dtset%natom+1)then 705 write(message,'(a,i4)')' Perturbation : derivative vs k along direction',idir 706 call wrtout(std_out,message,'COLL') 707 call wrtout(ab_out,message,'COLL') 708 if( iscf_mod /= -3 )then 709 write(message, '(4a)' )ch10,& 710 & ' dfpt_looppert : COMMENT -',ch10,& 711 & ' In a d/dk calculation, iscf is set to -3 automatically.' 712 call wrtout(std_out,message,'COLL') 713 call wrtout(ab_out,message,'COLL') 714 iscf_mod=-3 715 end if 716 if( abs(dtset%dfpt_sciss) > 1.0d-8 )then 717 write(message, '(a,a,a,a,f14.8,a,a)' )ch10,& 718 & ' dfpt_looppert : WARNING -',ch10,& 719 & ' Value of dfpt_sciss=',dtset%dfpt_sciss,ch10,& 720 & ' Scissor with d/dk calculation : you are using a "naive" approach !' 721 call wrtout(std_out,message,'COLL') 722 call wrtout(ab_out,message,'COLL') 723 end if 724 else if(ipert==dtset%natom+2)then 725 write(message, '(a,i4)' )' Perturbation : homogeneous electric field along direction',idir 726 call wrtout(std_out,message,'COLL') 727 call wrtout(ab_out,message,'COLL') 728 if( iscf_mod == -3 )then 729 write(message, '(a,a,a,a,a,a)' )ch10,& 730 & ' dfpt_looppert : COMMENT -',ch10,& 731 & ' The first-order density is imposed to be zero (iscf=-3).',ch10,& 732 & ' This corresponds to a calculation without local fields.' 733 call wrtout(std_out,message,'COLL') 734 call wrtout(ab_out,message,'COLL') 735 end if 736 else if(ipert==dtset%natom+10.or.ipert==dtset%natom+11)then 737 call rf2_getidirs(idir,idir1,idir2) 738 if(ipert==dtset%natom+10)then 739 write(message,'(2(a,i1))') ' Perturbation : 2nd derivative wrt k, idir1 = ',idir1,& 740 & ' idir2 = ',idir2 741 else 742 write(message,'(2(a,i1),a)') ' Perturbation : 2nd derivative wrt k (idir1 =',idir1,& 743 & ') and Efield (idir2 =',idir2,')' 744 end if 745 call wrtout(std_out,message,'COLL') 746 call wrtout(ab_out,message,'COLL') 747 if( iscf_mod /= -3 )then 748 write(message, '(4a)' )ch10,& 749 & ' dfpt_looppert : COMMENT -',ch10,& 750 & ' In this case, iscf is set to -3 automatically.' 751 call wrtout(std_out,message,'COLL') 752 call wrtout(ab_out,message,'COLL') 753 iscf_mod=-3 754 end if 755 if( abs(dtset%dfpt_sciss) > 1.0d-8 )then 756 write(message, '(a,a,a,a,f14.8,a,a)' )ch10,& 757 & ' dfpt_looppert : WARNING -',ch10,& 758 & ' Value of dfpt_sciss=',dtset%dfpt_sciss,ch10,& 759 & ' Scissor with d/dk calculation : you are using a "naive" approach !' 760 call wrtout(std_out,message,'COLL') 761 call wrtout(ab_out,message,'COLL') 762 end if 763 ABI_ALLOCATE(occ_pert,(dtset%mband*nkpt*dtset%nsppol)) 764 occ_pert(:) = occ(:) - occ(1) 765 maxocc = maxval(abs(occ_pert)) 766 if (maxocc>1.0d-6.and.abs(maxocc-occ(1))>1.0d-6) then ! True if non-zero occupation numbers are not equal 767 write(message, '(3a)' ) ' ipert=natom+10 or 11 does not work for a metallic system.',ch10,& 768 ' This perturbation will not be computed.' 769 MSG_WARNING(message) 770 ABI_DEALLOCATE(occ_pert) 771 cycle 772 end if 773 if (ipert==dtset%natom+11.and.minval(occ)<1.0d-6) then 774 write(message, '(3a)' ) ' ipert=natom+11 does not work with empty bands.',ch10,& 775 ' This perturbation will not be computed.' 776 MSG_WARNING(message) 777 ABI_DEALLOCATE(occ_pert) 778 cycle 779 end if 780 ABI_DEALLOCATE(occ_pert) 781 else if(ipert>dtset%natom+11 .or. ipert<=0 )then 782 write(message, '(a,i4,a,a,a)' ) & 783 & ' ipert=',ipert,' is outside the [1,dtset%natom+11] interval.',ch10,& 784 & ' This perturbation is not (yet) allowed.' 785 MSG_BUG(message) 786 end if 787 ! Initialize the diverse parts of energy : 788 eew=zero ; evdw=zero ; efrloc=zero ; efrnl=zero ; efrx1=zero ; efrx2=zero 789 efrhar=zero ; efrkin=zero 790 if(ipert<=dtset%natom)then 791 eew=dyew(1,idir,ipert,idir,ipert) 792 if (usevdw==1) evdw=dyvdw(1,idir,ipert,idir,ipert) 793 efrloc=dyfrlo(idir,idir,ipert) 794 if (dyfr_nondiag==0) efrnl=dyfrnl(1,idir,idir,ipert,1) 795 if (dyfr_nondiag/=0) efrnl=dyfrnl(1,idir,idir,ipert,ipert) 796 efrx1=dyfrx1(1,idir,ipert,idir,ipert) 797 efrx2=dyfrx2(idir,idir,ipert) 798 else if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then 799 ! istr = 1,2,...,6 and indicates the cartesian strain component 800 if(ipert==dtset%natom+4) istr=idir+3 801 eii=eltcore(istr,istr) 802 eew=elteew(istr,istr) 803 if (usevdw==1) evdw=eltvdw(istr,istr) 804 efrhar=eltfrhar(istr,istr) 805 efrkin=eltfrkin(istr,istr) 806 efrloc=eltfrloc(istr,istr) 807 efrnl=eltfrnl(istr,istr) 808 efrx1=eltfrxc(istr,istr) 809 end if 810 811 ! Determine the subset of symmetry operations (nsym1 operations) 812 ! that leaves the perturbation invariant, and initialize corresponding arrays 813 ! symaf1, symrl1, tnons1 (and pawang1%zarot, if PAW).. 814 ABI_ALLOCATE(symaf1_tmp,(nsym)) 815 ABI_ALLOCATE(symrl1_tmp,(3,3,nsym)) 816 ABI_ALLOCATE(tnons1_tmp,(3,nsym)) 817 if (dtset%prepanl/=1.and.& 818 & dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.& 819 & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17) then 820 call littlegroup_pert(gprimd,idir,indsym,ab_out,ipert,dtset%natom,nsym,nsym1,2,& 821 & dtset%symafm,symaf1_tmp,symq,symrec,dtset%symrel,symrl1_tmp,0,dtset%tnons,tnons1_tmp) 822 else 823 nsym1 = 1 824 symaf1_tmp(1) = 1 825 symrl1_tmp(:,:,1) = dtset%symrel(:,:,1) 826 tnons1_tmp(:,1) = 0_dp 827 end if 828 ABI_ALLOCATE(indsy1,(4,nsym1,dtset%natom)) 829 ABI_ALLOCATE(symrc1,(3,3,nsym1)) 830 ABI_ALLOCATE(symaf1,(nsym1)) 831 ABI_ALLOCATE(symrl1,(3,3,nsym1)) 832 ABI_ALLOCATE(tnons1,(3,nsym1)) 833 symaf1(1:nsym1)=symaf1_tmp(1:nsym1) 834 symrl1(:,:,1:nsym1)=symrl1_tmp(:,:,1:nsym1) 835 tnons1(:,1:nsym1)=tnons1_tmp(:,1:nsym1) 836 ABI_DEALLOCATE(symaf1_tmp) 837 ABI_DEALLOCATE(symrl1_tmp) 838 ABI_DEALLOCATE(tnons1_tmp) 839 840 ! Set up corresponding symmetry data 841 ABI_ALLOCATE(irrzon1,(dtset%nfft**(1-1/nsym1),2,(nspden/dtset%nsppol)-3*(nspden/4))) 842 ABI_ALLOCATE(phnons1,(2,dtset%nfft**(1-1/nsym1),(nspden/dtset%nsppol)-3*(nspden/4))) 843 call setsym(indsy1,irrzon1,1,dtset%natom,dtset%nfft,dtset%ngfft,nspden,dtset%nsppol,& 844 & nsym1,phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred) 845 if (psps%usepaw==1) then 846 ! Allocate/initialize only zarot in pawang1 datastructure 847 call pawang_init(pawang1,0,pawang%l_max-1,0,nsym1,0,1,0,0,0) 848 call setsymrhoij(gprimd,pawang1%l_max-1,pawang1%nsym,0,rprimd,symrc1,pawang1%zarot) 849 end if 850 851 ! Initialize k+q array 852 ABI_ALLOCATE(kpq,(3,nkpt)) 853 if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then 854 kpq(:,1:nkpt)=dtset%kptns(:,1:nkpt) ! Do not modify, needed for gfortran 855 else 856 do ikpt=1,nkpt 857 kpq(:,ikpt)=dtset%qptn(:)+dtset%kptns(:,ikpt) 858 end do 859 end if 860 ! In case wf1 at +q and -q are not related by time inversion symmetry, compute k-q as well 861 ! for initializations 862 if (.not.kramers_deg) then 863 ABI_ALLOCATE(kmq,(3,nkpt)) 864 if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then 865 kmq(:,1:nkpt)=dtset%kptns(:,1:nkpt) ! Do not modify, needed for gfortran 866 else 867 do ikpt=1,nkpt 868 kmq(:,ikpt)=-dtset%qptn(:)+dtset%kptns(:,ikpt) ! kmq <= k-q 869 end do 870 end if 871 end if 872 873 ! Determine the subset of k-points needed in the "reduced Brillouin zone", 874 ! and initialize other quantities 875 ABI_ALLOCATE(indkpt1_tmp,(nkpt)) 876 ABI_ALLOCATE(wtk_folded,(nkpt)) 877 indkpt1_tmp(:)=0 ; optthm=0 878 timrev_pert=timrev 879 if(dtset%ieig2rf>0) then 880 timrev_pert=0 881 call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,& 882 & 1,symrc1,timrev_pert,dtset%wtk,wtk_folded) 883 else 884 ! For the time being, the time reversal symmetry is not used 885 ! for ddk, elfd, mgfd perturbations. 886 timrev_pert=timrev 887 if(ipert==dtset%natom+1.or.ipert==dtset%natom+2.or.& 888 & ipert==dtset%natom+10.or.ipert==dtset%natom+11.or. & 889 & dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or. & 890 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17.or. & 891 & ipert==dtset%natom+5) timrev_pert=0 892 timrev_kpt = timrev_pert 893 ! The time reversal symmetry is not used for the BZ sampling when kptopt=3 or 4 894 if (dtset%kptopt==3.or.dtset%kptopt==4) timrev_kpt = 0 895 call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,& 896 nsym1,symrc1,timrev_kpt,dtset%wtk,wtk_folded) 897 end if 898 899 ABI_ALLOCATE(doccde_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol)) 900 ABI_ALLOCATE(indkpt1,(nkpt_rbz)) 901 ABI_ALLOCATE(istwfk_rbz,(nkpt_rbz)) 902 ABI_ALLOCATE(kpq_rbz,(3,nkpt_rbz)) 903 if (.not.kramers_deg) then 904 ABI_ALLOCATE(kmq_rbz,(3,nkpt_rbz)) 905 end if 906 ABI_ALLOCATE(kpt_rbz,(3,nkpt_rbz)) 907 ABI_ALLOCATE(nband_rbz,(nkpt_rbz*dtset%nsppol)) 908 ABI_ALLOCATE(occ_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol)) 909 ABI_ALLOCATE(wtk_rbz,(nkpt_rbz)) 910 indkpt1(:)=indkpt1_tmp(1:nkpt_rbz) 911 do ikpt=1,nkpt_rbz 912 istwfk_rbz(ikpt)=dtset%istwfk(indkpt1(ikpt)) 913 kpq_rbz(:,ikpt)=kpq(:,indkpt1(ikpt)) 914 kpt_rbz(:,ikpt)=dtset%kptns(:,indkpt1(ikpt)) 915 wtk_rbz(ikpt)=wtk_folded(indkpt1(ikpt)) 916 end do 917 if (.not.kramers_deg) then 918 do ikpt=1,nkpt_rbz 919 kmq_rbz(:,ikpt)=kmq(:,indkpt1(ikpt)) 920 end do 921 end if 922 ABI_DEALLOCATE(indkpt1_tmp) 923 ABI_DEALLOCATE(wtk_folded) 924 925 ! Transfer occ to occ_rbz and doccde to doccde_rbz : 926 ! this is a more delicate issue 927 ! NOTE : this takes into account that indkpt1 is ordered 928 ! MG: What about using occ(band,kpt,spin) ??? 929 bdtot_index=0;bdtot1_index=0 930 do isppol=1,dtset%nsppol 931 ikpt1=1 932 do ikpt=1,nkpt 933 nband_k=dtset%nband(ikpt+(isppol-1)*nkpt) 934 ! Must test against ikpt1/=nkpt_rbz+1, before evaluate indkpt1(ikpt1) 935 if(ikpt1/=nkpt_rbz+1)then 936 if(ikpt==indkpt1(ikpt1))then 937 nband_rbz(ikpt1+(isppol-1)*nkpt_rbz)=nband_k 938 occ_rbz(1+bdtot1_index:nband_k+bdtot1_index) = occ(1+bdtot_index:nband_k+bdtot_index) 939 doccde_rbz(1+bdtot1_index:nband_k+bdtot1_index) = doccde(1+bdtot_index:nband_k+bdtot_index) 940 ikpt1=ikpt1+1 941 bdtot1_index=bdtot1_index+nband_k 942 end if 943 end if 944 bdtot_index=bdtot_index+nband_k 945 end do 946 end do 947 948 _IBM6("IBM6 in dfpt_looppert before getmpw") 949 950 ! Compute maximum number of planewaves at k 951 call timab(142,1,tsec) 952 call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpt_rbz,mpi_enreg,mpw,nkpt_rbz) 953 call timab(142,2,tsec) 954 955 ! Allocate some k-dependent arrays at k 956 ABI_ALLOCATE(kg,(3,mpw*nkpt_rbz)) 957 ABI_ALLOCATE(npwarr,(nkpt_rbz)) 958 ABI_ALLOCATE(npwtot,(nkpt_rbz)) 959 960 ! Determine distribution of k-points/bands over MPI processes 961 if (allocated(mpi_enreg%my_kpttab)) then 962 ABI_DEALLOCATE(mpi_enreg%my_kpttab) 963 end if 964 ABI_ALLOCATE(mpi_enreg%my_kpttab,(nkpt_rbz)) 965 if(xmpi_paral==1) then 966 ABI_ALLOCATE(mpi_enreg%proc_distrb,(nkpt_rbz,dtset%mband,dtset%nsppol)) 967 call distrb2(dtset%mband,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,dtset%nsppol,mpi_enreg) 968 else 969 mpi_enreg%my_kpttab(:)=(/(ii,ii=1,nkpt_rbz)/) 970 end if 971 my_nkpt_rbz=maxval(mpi_enreg%my_kpttab) 972 call initmpi_band(mpi_enreg,nband_rbz,nkpt_rbz,dtset%nsppol) 973 mkmem_rbz =my_nkpt_rbz ; mkqmem_rbz=my_nkpt_rbz ; mk1mem_rbz=my_nkpt_rbz 974 ABI_UNUSED((/mkmem,mk1mem,mkqmem/)) 975 976 _IBM6("IBM6 before kpgio") 977 978 ! Set up the basis sphere of planewaves at k 979 call timab(143,1,tsec) 980 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg,& 981 & kpt_rbz,mkmem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw,npwarr,npwtot,dtset%nsppol) 982 call timab(143,2,tsec) 983 984 ! Set up the spherical harmonics (Ylm) at k 985 useylmgr=0; option=0 ; nylmgr=0 986 if (psps%useylm==1.and. & 987 & (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or. & 988 & (psps%usepaw==1.and.ipert==dtset%natom+2))) then 989 useylmgr=1; option=1 ; nylmgr=3 990 else if (psps%useylm==1.and.(ipert==dtset%natom+10.or.ipert==dtset%natom+11)) then 991 useylmgr=1; option=2 ; nylmgr=9 992 end if 993 ABI_ALLOCATE(ylm,(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm)) 994 ABI_ALLOCATE(ylmgr,(mpw*mkmem_rbz,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 995 if (psps%useylm==1) then 996 call initylmg(gprimd,kg,kpt_rbz,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,nband_rbz,nkpt_rbz,& 997 & npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr) 998 end if 999 1000 _IBM6("Before ieig2rf > 0") 1001 1002 ! Set up occupations for this perturbation 1003 if (dtset%ieig2rf>0) then 1004 if (.not.allocated(istwfk_pert)) then 1005 ABI_ALLOCATE(istwfk_pert,(nkpt,3,mpert)) 1006 ABI_ALLOCATE(occ_pert,(dtset%mband*nkpt*dtset%nsppol)) 1007 istwfk_pert(:,:,:)=0 ; occ_pert(:)=zero 1008 end if 1009 istwfk_pert(:,idir,ipert)=istwfk_rbz(:) 1010 occ_pert(:)=occ_rbz(:) 1011 end if 1012 if (dtset%efmas>0) then 1013 if (.not.allocated(istwfk_pert)) then 1014 ABI_ALLOCATE(istwfk_pert,(nkpt,3,mpert)) 1015 istwfk_pert(:,:,:)=0 1016 end if 1017 istwfk_pert(:,idir,ipert)=istwfk_rbz(:) 1018 end if 1019 1020 ! Print a separator in output file 1021 write(message,'(3a)')ch10,'--------------------------------------------------------------------------------',ch10 1022 call wrtout(ab_out,message,'COLL') 1023 1024 ! Initialize band structure datatype at k 1025 bantot_rbz=sum(nband_rbz(1:nkpt_rbz*dtset%nsppol)) 1026 ABI_ALLOCATE(eigen0,(bantot_rbz)) 1027 eigen0(:)=zero 1028 call ebands_init(bantot_rbz,ebands_k,dtset%nelect,doccde_rbz,eigen0,istwfk_rbz,kpt_rbz,& 1029 & nband_rbz,nkpt_rbz,npwarr,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,& 1030 & dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 1031 & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 1032 ABI_DEALLOCATE(eigen0) 1033 1034 ! Initialize header, update it with evolving variables 1035 gscase=0 ! A GS WF file is read 1036 call hdr_init(ebands_k,codvsn,dtset,hdr0,pawtab,gscase,psps,wvl%descr,& 1037 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1038 1039 call hdr_update(hdr0,bantot_rbz,etotal,fermie,& 1040 & residm,rprimd,occ_rbz,pawrhoij_pert,xred,dtset%amu_orig(:,1),& 1041 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1042 1043 _IBM6("before inwffil") 1044 1045 ! Initialize GS wavefunctions at k 1046 ireadwf0=1; formeig=0 ; ask_accurate=1 ; optorth=0 1047 mcg=mpw*dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol 1048 if (one*mpw*dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol > huge(1)) then 1049 write (message,'(4a, 5(a,i0), 2a)')& 1050 & "Default integer is not wide enough to store the size of the GS wavefunction array (WF0, mcg).",ch10,& 1051 & "Action: increase the number of processors. Consider also OpenMP threads.",ch10,& 1052 & "nspinor: ",dtset%nspinor, "mpw: ",mpw, "mband: ",dtset%mband, "mkmem_rbz: ",& 1053 & mkmem_rbz, "nsppol: ",dtset%nsppol,ch10,& 1054 & 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode' 1055 MSG_ERROR(message) 1056 end if 1057 ABI_STAT_ALLOCATE(cg,(2,mcg), ierr) 1058 ABI_CHECK(ierr==0, "out-of-memory in cg") 1059 1060 ABI_ALLOCATE(eigen0,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1061 call timab(144,1,tsec) 1062 call status(0,dtfil%filstat,iexit,level,'call inwffil-k') 1063 call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,& 1064 & formeig,hdr0,ireadwf0,istwfk_rbz,kg,& 1065 & kpt_rbz,dtset%localrdwf,dtset%mband,mcg,& 1066 & mkmem_rbz,mpi_enreg,mpw,nband_rbz,dtset%ngfft,nkpt_rbz,npwarr,& 1067 & dtset%nsppol,nsym,occ_rbz,optorth,dtset%symafm,& 1068 & dtset%symrel,dtset%tnons,dtfil%unkg,wffgs,wfftgs,& 1069 & dtfil%unwffgs,dtfil%fnamewffk,wvl) 1070 call timab(144,2,tsec) 1071 ! Close wffgs%unwff, if it was ever opened (in inwffil) 1072 if (ireadwf0==1) then 1073 call WffClose(wffgs,ierr) 1074 end if 1075 ! Update energies GS energies at k 1076 call put_eneocc_vect(ebands_k, "eig", eigen0) 1077 1078 ! PAW: compute on-site projections of GS wavefunctions (cprj) (and derivatives) at k 1079 ncpgr=0 1080 ABI_DATATYPE_ALLOCATE(cprj,(0,0)) 1081 if (psps%usepaw==1) then 1082 ncpgr=3 ! Valid for ipert<=natom (phonons), ipert=natom+2 (elec. field) 1083 ! or for ipert==natom+10,11 1084 if (ipert==dtset%natom+1) ncpgr=1 1085 if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) ncpgr=1 1086 if (usecprj==1) then 1087 mcprj=dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol 1088 ABI_DATATYPE_DEALLOCATE(cprj) 1089 ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,mcprj)) 1090 call pawcprj_alloc(cprj,ncpgr,dimcprj_srt) 1091 if (ipert<=dtset%natom) then 1092 choice=2; iorder_cprj=0; idir0=0 1093 else if (ipert==dtset%natom+1) then 1094 choice=5; iorder_cprj=0; idir0=idir 1095 else if (ipert==dtset%natom+2) then 1096 choice=5; iorder_cprj=0; idir0=0 1097 else if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then 1098 choice=3; iorder_cprj=0; idir0=istr 1099 else if (ipert==dtset%natom+10.or.ipert==dtset%natom+11) then 1100 choice=5; iorder_cprj=0; idir0=0 ! Compute all first derivatives 1101 else 1102 choice=1; iorder_cprj=0; idir0=0 1103 end if 1104 call ctocprj(atindx,cg,choice,cprj,gmet,gprimd,-1,idir0,iorder_cprj,istwfk_rbz,& 1105 & kg,kpt_rbz,mcg,mcprj,dtset%mgfft,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,& 1106 & dtset%natom,nattyp,nband_rbz,dtset%natom,dtset%ngfft,nkpt_rbz,dtset%nloalg,& 1107 & npwarr,dtset%nspinor,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,psps,& 1108 & rmet,dtset%typat,ucvol,dtfil%unpaw,xred,ylm,ylmgr) 1109 end if 1110 end if 1111 1112 ! Compute maximum number of planewaves at k+q 1113 ! Will be useful for both GS wfs at k+q and RF wavefunctions 1114 call timab(143,1,tsec) 1115 call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpq_rbz,mpi_enreg,mpw1,nkpt_rbz) 1116 if (.not.kramers_deg) then 1117 call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kmq_rbz,mpi_enreg,mpw1_mq,nkpt_rbz) 1118 !number of plane waves at k+q and k-q should be in principle the same to reconstruct rhor1_pq (?) 1119 !mpw1=max(mpw1,mpw1_tmp) 1120 else 1121 mpw1_mq=0 1122 end if 1123 call timab(143,2,tsec) 1124 1125 ! Allocate some arrays at k+q 1126 ABI_ALLOCATE(kg1,(3,mpw1*mk1mem_rbz)) 1127 ABI_ALLOCATE(npwar1,(nkpt_rbz)) 1128 ABI_ALLOCATE(npwtot1,(nkpt_rbz)) 1129 ! In case Kramers degeneracy is broken, do the same for k-q 1130 if (.not.kramers_deg) then 1131 ABI_ALLOCATE(kg1_mq,(3,mpw1_mq*mk1mem_rbz)) 1132 ABI_ALLOCATE(npwar1_mq,(nkpt_rbz)) 1133 ABI_ALLOCATE(npwtot1_mq,(nkpt_rbz)) 1134 end if 1135 1136 ! Set up the basis sphere of planewaves at k+q 1137 ! Will be useful for both GS wfs at k+q and RF wavefunctions 1138 call timab(142,1,tsec) 1139 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg1,& 1140 & kpq_rbz,mk1mem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw1,& 1141 & npwar1,npwtot1,dtset%nsppol) 1142 if (.not.kramers_deg) then 1143 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg1_mq,& 1144 & kmq_rbz,mk1mem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw1_mq,& 1145 & npwar1_mq,npwtot1_mq,dtset%nsppol) 1146 end if 1147 call timab(142,2,tsec) 1148 1149 ! Set up the spherical harmonics (Ylm) at k+q 1150 useylmgr1=0; option=0 ; ; nylmgr1=0 1151 if (psps%useylm==1.and. & 1152 & (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or. & 1153 & (psps%usepaw==1.and.ipert==dtset%natom+2))) then 1154 useylmgr1=1; option=1 ; ; nylmgr1=3 1155 else if (psps%useylm==1.and.(ipert==dtset%natom+10.or.ipert==dtset%natom+11)) then 1156 useylmgr1=1; option=2 ; ; nylmgr1=9 1157 end if 1158 ABI_ALLOCATE(ylm1,(mpw1*mk1mem_rbz,psps%mpsang*psps%mpsang*psps%useylm)) 1159 ABI_ALLOCATE(ylmgr1,(mpw1*mk1mem_rbz,nylmgr1,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 1160 if (psps%useylm==1) then 1161 call initylmg(gprimd,kg1,kpq_rbz,mk1mem_rbz,mpi_enreg,psps%mpsang,mpw1,nband_rbz,nkpt_rbz,& 1162 & npwar1,dtset%nsppol,option,rprimd,ylm1,ylmgr1) 1163 end if 1164 ! SPr: do the same for k-q if kramers_deg=.false. 1165 1166 ! Print a separator in output file 1167 write(message, '(a,a)' )'--------------------------------------------------------------------------------',ch10 1168 call wrtout(ab_out,message,'COLL') 1169 1170 ! Initialize band structure datatype at k+q 1171 ABI_ALLOCATE(eigenq,(bantot_rbz)) 1172 eigenq(:)=zero 1173 call ebands_init(bantot_rbz,ebands_kq,dtset%nelect,doccde_rbz,eigenq,istwfk_rbz,kpq_rbz,& 1174 & nband_rbz,nkpt_rbz,npwar1,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,& 1175 & dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 1176 & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 1177 if (.not.kramers_deg) then 1178 eigenq(:)=zero 1179 call ebands_init(bantot_rbz,ebands_kmq,dtset%nelect,doccde_rbz,eigenq,istwfk_rbz,kmq_rbz,& 1180 & nband_rbz,nkpt_rbz,npwar1_mq,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,& 1181 & dtset%charge, dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, & 1182 & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk) 1183 end if 1184 ABI_DEALLOCATE(eigenq) 1185 1186 ! Initialize header 1187 call hdr_init(ebands_kq,codvsn,dtset,hdr,pawtab,pertcase,psps,wvl%descr, & 1188 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab ) 1189 if (.not.kramers_deg) then 1190 call hdr_init(ebands_kmq,codvsn,dtset,hdr,pawtab,pertcase,psps,wvl%descr, & 1191 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab ) 1192 end if 1193 1194 ! Initialize wavefunctions at k+q 1195 ! MG: Here it is possible to avoid the extra reading if the same k mesh can be used. 1196 ireadwf0=1 ; formeig=0 ; ask_accurate=1 ; optorth=0 1197 mcgq=mpw1*dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol 1198 !SPr: verified until here, add mcgq for -q 1199 if (one*mpw1*dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol > huge(1)) then 1200 write (message,'(4a, 5(a,i0), 2a)')& 1201 & "Default integer is not wide enough to store the size of the GS wavefunction array (WFKQ, mcgq).",ch10,& 1202 & "Action: increase the number of processors. Consider also OpenMP threads.",ch10,& 1203 & "nspinor: ",dtset%nspinor, "mpw1: ",mpw1, "mband: ",dtset%mband, "mkqmem_rbz: ",& 1204 & mkqmem_rbz, "nsppol: ",dtset%nsppol,ch10,& 1205 & 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode' 1206 MSG_ERROR(message) 1207 end if 1208 1209 ABI_STAT_ALLOCATE(cgq,(2,mcgq), ierr) 1210 ABI_CHECK(ierr==0, "out-of-memory in cgq") 1211 ABI_ALLOCATE(eigenq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1212 if (.not.kramers_deg) then 1213 !ABI_STAT_ALLOCATE(cg_pq,(2,mcgq), ierr) 1214 !ABI_CHECK(ierr==0, "out-of-memory in cgmq") 1215 !ABI_ALLOCATE(eigen_pq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1216 mcgmq=mpw1_mq*dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol 1217 ABI_STAT_ALLOCATE(cg_mq,(2,mcgmq), ierr) 1218 ABI_CHECK(ierr==0, "out-of-memory in cgmq") 1219 ABI_ALLOCATE(eigen_mq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1220 end if 1221 1222 !if (sum(dtset%qptn(1:3)**2)>=1.d-14) then ! non-zero q 1223 if (dtfil%fnamewffq == dtfil%fnamewffk .and. sum(dtset%qptn(1:3)**2) < 1.d-14) then 1224 call wrtout(std_out, "qpt is Gamma, psi_k+q initialized from psi_k in memory") 1225 cgq = cg 1226 eigenq = eigen0 1227 else 1228 call timab(144,1,tsec) 1229 call status(0,dtfil%filstat,iexit,level,'call inwffilkq') 1230 call inwffil(ask_accurate,cgq,dtset,dtset%ecut,ecut_eff,eigenq,dtset%exchn2n3d,& 1231 & formeig,hdr,& 1232 & ireadwf0,istwfk_rbz,kg1,kpq_rbz,dtset%localrdwf,dtset%mband,mcgq,& 1233 & mkqmem_rbz,mpi_enreg,mpw1,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1,& 1234 & dtset%nsppol,nsym,occ_rbz,optorth,& 1235 & dtset%symafm,dtset%symrel,dtset%tnons,& 1236 & dtfil%unkg1,wffkq,wfftkq,dtfil%unwffkq,dtfil%fnamewffq,wvl) 1237 call timab(144,2,tsec) 1238 ! Close dtfil%unwffkq, if it was ever opened (in inwffil) 1239 if (ireadwf0==1) then 1240 call WffClose(wffkq,ierr) 1241 end if 1242 if (.not.kramers_deg) then 1243 !SPr: later "make" a separate WFQ file for "-q" 1244 call timab(144,1,tsec) 1245 call status(0,dtfil%filstat,iexit,level,'call inwffilkq') 1246 call inwffil(ask_accurate,cg_mq,dtset,dtset%ecut,ecut_eff,eigen_mq,dtset%exchn2n3d,& 1247 & formeig,hdr,& 1248 & ireadwf0,istwfk_rbz,kg1_mq,kmq_rbz,dtset%localrdwf,dtset%mband,mcgmq,& 1249 & mkqmem_rbz,mpi_enreg,mpw1_mq,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1_mq,& 1250 & dtset%nsppol,nsym,occ_rbz,optorth,& 1251 & dtset%symafm,dtset%symrel,dtset%tnons,& 1252 & dtfil%unkg1,wffkq,wfftkq,dtfil%unwffkq,dtfil%fnamewffq,wvl) 1253 call timab(144,2,tsec) 1254 ! Close dtfil%unwffkq, if it was ever opened (in inwffil) 1255 if (ireadwf0==1) then 1256 call WffClose(wffkq,ierr) 1257 end if 1258 end if 1259 end if 1260 ! Update energies GS energies at k + q 1261 call put_eneocc_vect(ebands_kq, "eig", eigenq) 1262 if (.not.kramers_deg) then 1263 call put_eneocc_vect(ebands_kmq, "eig", eigen_mq) 1264 end if 1265 1266 ! PAW: compute on-site projections of GS wavefunctions (cprjq) (and derivatives) at k+q 1267 ABI_DATATYPE_ALLOCATE(cprjq,(0,0)) 1268 if (psps%usepaw==1) then 1269 if (usecprj==1) then 1270 call status(0,dtfil%filstat,iexit,level,'call ctocprjq') 1271 mcprjq=dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol 1272 ABI_DATATYPE_DEALLOCATE(cprjq) 1273 ABI_DATATYPE_ALLOCATE(cprjq,(dtset%natom,mcprjq)) 1274 call pawcprj_alloc(cprjq,0,dimcprj_srt) 1275 if (ipert<=dtset%natom.and.(sum(dtset%qptn(1:3)**2)>=1.d-14)) then ! phonons at non-zero q 1276 choice=1 ; iorder_cprj=0 ; idir0=0 1277 call ctocprj(atindx,cgq,choice,cprjq,gmet,gprimd,-1,idir0,0,istwfk_rbz,& 1278 & kg1,kpq_rbz,mcgq,mcprjq,dtset%mgfft,mkqmem_rbz,mpi_enreg,psps%mpsang,mpw1,& 1279 & dtset%natom,nattyp,nband_rbz,dtset%natom,dtset%ngfft,nkpt_rbz,dtset%nloalg,& 1280 & npwar1,dtset%nspinor,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,& 1281 & psps,rmet,dtset%typat,ucvol,dtfil%unpawq,xred,ylm1,ylmgr1) 1282 else if (mcprjq>0) then 1283 call pawcprj_copy(cprj,cprjq) 1284 end if 1285 end if 1286 end if 1287 1288 ! ===== Report on eigenq values 1289 if (dtset%ieig2rf>0.and.icase==ipert_cnt) then 1290 eigen0_pert(:) = eigen0(:) 1291 eigenq_pert(:) = eigenq(:) 1292 occ_rbz_pert(:) = occ_rbz(:) 1293 end if 1294 if (dtset%efmas>0.and.icase==ipert_cnt) then 1295 eigen0_pert(:) = eigen0(:) 1296 end if 1297 call wrtout(std_out,ch10//' dfpt_looppert : eigenq array',"COLL") 1298 nkpt_eff=nkpt 1299 if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. nkpt>nkpt_max ) nkpt_eff=nkpt_max 1300 band_index=0 1301 do isppol=1,dtset%nsppol 1302 do ikpt=1,nkpt_rbz 1303 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 1304 if(ikpt<=nkpt_eff)then 1305 write(message, '(a,i2,a,i5)' )' isppol=',isppol,', k point number',ikpt 1306 call wrtout(std_out,message,'COLL') 1307 do iband=1,nband_k,4 1308 write(message, '(a,4es16.6)')' ',eigenq(iband+band_index:min(iband+3,nband_k)+band_index) 1309 call wrtout(std_out,message,'COLL') 1310 end do 1311 else if(ikpt==nkpt_eff+1)then 1312 write(message,'(a,a)' )' respfn : prtvol=0, 1 or 2, stop printing eigenq.',ch10 1313 call wrtout(std_out,message,'COLL') 1314 end if 1315 band_index=band_index+nband_k 1316 end do 1317 end do 1318 1319 ! Generate occupation numbers for the reduced BZ at k+q 1320 ABI_ALLOCATE(docckqde,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1321 ABI_ALLOCATE(occkq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1322 if (.not.kramers_deg) then 1323 ABI_ALLOCATE(docckde_mq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1324 ABI_ALLOCATE(occk_mq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1325 end if 1326 1327 if(0<=dtset%occopt .and. dtset%occopt<=2)then 1328 ! Same occupation numbers at k and k+q (usually, insulating) 1329 occkq(:)=occ_rbz(:) 1330 docckqde(:)=zero ! docckqde is irrelevant in this case 1331 if(.not.kramers_deg) then 1332 occk_mq(:)=occ_rbz(:) 1333 docckde_mq(:)=zero 1334 end if 1335 else 1336 ! Metallic occupation numbers 1337 option=1 1338 dosdeltae=zero ! the DOS is not computed with option=1 1339 maxocc=two/(dtset%nspinor*dtset%nsppol) 1340 call getnel(docckqde,dosdeltae,eigenq,entropy,fermie,maxocc,dtset%mband,& 1341 & nband_rbz,nelectkq,nkpt_rbz,dtset%nsppol,occkq,dtset%occopt,option,& 1342 & dtset%tphysel,dtset%tsmear,fake_unit,wtk_rbz) 1343 ! Compare nelect at k and nelelect at k+q 1344 write(message, '(a,a,a,es16.6,a,es16.6,a)')& 1345 & ' dfpt_looppert : total number of electrons, from k and k+q',ch10,& 1346 & ' fully or partially occupied states are',dtset%nelect,' and',nelectkq,'.' 1347 call wrtout(ab_out,message,'COLL') 1348 call wrtout(std_out,message,'COLL') 1349 if (.not.kramers_deg) then 1350 call getnel(docckde_mq,dosdeltae,eigen_mq,entropy,fermie,maxocc,dtset%mband,& 1351 & nband_rbz,nelectkq,nkpt_rbz,dtset%nsppol,occk_mq,dtset%occopt,option,& 1352 & dtset%tphysel,dtset%tsmear,fake_unit,wtk_rbz) 1353 ! Compare nelect at k and nelelect at k-q 1354 write(message, '(a,a,a,es16.6,a,es16.6,a)')& 1355 & ' dfpt_looppert : total number of electrons, from k and k-q',ch10,& 1356 & ' fully or partially occupied states are',dtset%nelect,' and',nelectkq,'.' 1357 call wrtout(ab_out,message,'COLL') 1358 call wrtout(std_out,message,'COLL') 1359 end if 1360 end if 1361 1362 ! Debug message 1363 if(dtset%prtvol==-level)then 1364 call wrtout(std_out,'dfpt_looppert : initialisation of q part done. ','COLL') 1365 end if 1366 1367 ! Initialisation of first-order wavefunctions 1368 write(message,'(3a,i4)')& 1369 & ' Initialisation of the first-order wave-functions :',ch10,& 1370 & ' ireadwf=',dtfil%ireadwf 1371 call wrtout(std_out,message,'COLL') 1372 call wrtout(ab_out,message,'COLL') 1373 call appdig(pertcase,dtfil%fnamewff1,fiwf1i) 1374 call appdig(pertcase,dtfil%fnameabo_1wf,fiwf1o) 1375 1376 ! Allocate 1st-order PAW occupancies (rhoij1) 1377 if (psps%usepaw==1) then 1378 ABI_DATATYPE_ALLOCATE(pawrhoij1,(my_natom)) 1379 call pawrhoij_nullify(pawrhoij1) 1380 cplex_rhoij=max(cplex,dtset%pawcpxocc) 1381 nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb) 1382 call pawrhoij_alloc(pawrhoij1,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,& 1383 & dtset%typat,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1384 if (cplex_rhoij/=hdr%pawrhoij(1)%cplex) then 1385 ! Eventually reallocate hdr%pawrhoij 1386 call pawrhoij_free(hdr%pawrhoij) 1387 call pawrhoij_alloc(hdr%pawrhoij,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,& 1388 & dtset%typat,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1389 end if 1390 else 1391 ABI_DATATYPE_ALLOCATE(pawrhoij1,(0)) 1392 end if 1393 1394 ! Initialize 1st-order wavefunctions 1395 formeig=1; ask_accurate=0; optorth=0 1396 ! NB: 4 Sept 2013: this was introducing a bug - for ieig2rf ==0 the dim was being set to 1 and passed to dfpt_vtowfk 1397 dim_eig2rf=0 1398 if ((dtset%ieig2rf > 0 .and. dtset%ieig2rf/=2) .or. dtset%efmas > 0) then 1399 dim_eig2rf=1 1400 end if 1401 mcg1=mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol 1402 if (one*mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol > huge(1)) then 1403 write (message,'(4a, 5(a,i0), 2a)')& 1404 & "Default integer is not wide enough to store the size of the GS wavefunction array (WFK1, mcg1).",ch10,& 1405 & "Action: increase the number of processors. Consider also OpenMP threads.",ch10,& 1406 & "nspinor: ",dtset%nspinor, "mpw1: ",mpw1, "mband: ",dtset%mband, "mk1mem_rbz: ",& 1407 & mk1mem_rbz, "nsppol: ",dtset%nsppol,ch10,& 1408 & 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode' 1409 MSG_ERROR(message) 1410 end if 1411 ABI_STAT_ALLOCATE(cg1,(2,mcg1), ierr) 1412 ABI_CHECK(ierr==0, "out of memory in cg1") 1413 if (.not.kramers_deg) then 1414 mcg1mq=mpw1_mq*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol 1415 ABI_STAT_ALLOCATE(cg1_mq,(2,mcg1mq), ierr) 1416 ABI_CHECK(ierr==0, "out of memory in cg1_mq") 1417 end if 1418 1419 ABI_ALLOCATE(cg1_active,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1420 ABI_ALLOCATE(gh1c_set,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1421 ABI_ALLOCATE(gh0c1_set,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1422 if (.not.kramers_deg) then 1423 ABI_ALLOCATE(cg1_active_mq,(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1424 ABI_ALLOCATE(gh1c_set_mq,(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1425 ABI_ALLOCATE(gh0c1_set_mq,(2,mpw1_mq*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1426 end if 1427 ! XG090606 This is needed in the present 5.8.2 , for portability for the pathscale machine. 1428 ! However, it is due to a bug to be corrected by Paul Boulanger. When the bug will be corrected, 1429 ! this line should be removed. 1430 if(mk1mem_rbz/=0 .and. dtset%ieig2rf/=0)then 1431 cg1_active=zero 1432 gh1c_set=zero 1433 gh0c1_set=zero 1434 if (.not.kramers_deg) then 1435 cg1_active_mq=zero 1436 gh1c_set_mq=zero 1437 gh0c1_set_mq=zero 1438 end if 1439 end if 1440 ABI_ALLOCATE(eigen1,(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol)) 1441 ABI_ALLOCATE(resid,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1442 call timab(144,1,tsec) 1443 call status(pertcase,dtfil%filstat,iexit,level,'call inwffil ') 1444 call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,& 1445 & formeig,hdr,& 1446 & dtfil%ireadwf,istwfk_rbz,kg1,kpq_rbz,dtset%localrdwf,& 1447 & dtset%mband,mcg1,mk1mem_rbz,mpi_enreg,mpw1,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1,& 1448 & dtset%nsppol,nsym1,occ_rbz,optorth,& 1449 & symaf1,symrl1,tnons1,dtfil%unkg1,wff1,wffnow,dtfil%unwff1,& 1450 & fiwf1i,wvl) 1451 call timab(144,2,tsec) 1452 ! Close wff1 (filename fiwf1i), if it was ever opened (in inwffil) 1453 if (dtfil%ireadwf==1) then 1454 call WffClose(wff1,ierr) 1455 end if 1456 if(.not.kramers_deg) then 1457 ABI_ALLOCATE(eigen1_mq,(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol)) 1458 ABI_ALLOCATE(resid_mq,(dtset%mband*nkpt_rbz*dtset%nsppol)) 1459 !initialize cg1_mq: 1460 call timab(144,1,tsec) 1461 call status(pertcase,dtfil%filstat,iexit,level,'call inwffil ') 1462 call inwffil(ask_accurate,cg1_mq,dtset,dtset%ecut,ecut_eff,eigen1_mq,dtset%exchn2n3d,& 1463 & formeig,hdr,& 1464 & dtfil%ireadwf,istwfk_rbz,kg1_mq,kmq_rbz,dtset%localrdwf,& 1465 & dtset%mband,mcg1mq,mk1mem_rbz,mpi_enreg,mpw1_mq,nband_rbz,dtset%ngfft,nkpt_rbz,npwar1_mq,& 1466 & dtset%nsppol,nsym1,occ_rbz,optorth,& 1467 & symaf1,symrl1,tnons1,dtfil%unkg1,wff1,wffnow,dtfil%unwff1,& 1468 & fiwf1i,wvl) 1469 call timab(144,2,tsec) 1470 ! Close wff1 (filename fiwf1i), if it was ever opened (in inwffil) 1471 if (dtfil%ireadwf==1) then 1472 call WffClose(wff1,ierr) 1473 end if 1474 end if 1475 1476 ! Eventually reytrieve 1st-order PAW occupancies from file header 1477 if (psps%usepaw==1.and.dtfil%ireadwf/=0) then 1478 call pawrhoij_copy(hdr%pawrhoij,pawrhoij1,comm_atom=mpi_enreg%comm_atom , & 1479 & mpi_atmtab=mpi_enreg%my_atmtab) 1480 end if 1481 1482 ! In case of electric field, or 2nd order perturbation : open the ddk (or dE) wf file(s) 1483 if ((ipert==dtset%natom+2.and.sum((dtset%qptn(1:3))**2)<1.0d-7.and. & 1484 & (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and. & 1485 & dtset%berryopt/= 7.and.dtset%berryopt/=14.and. & 1486 & dtset%berryopt/=16.and.dtset%berryopt/=17)) .or. & 1487 & ipert==dtset%natom+10.or.ipert==dtset%natom+11) then 1488 if (ipert<dtset%natom+10) then 1489 ! 1st order or berry phase, one direction 1490 nwffile=1 ! one direction needed 1491 file_index(1)=idir+dtset%natom*3 1492 fnamewff(1)=dtfil%fnamewffddk 1493 else if (ipert==dtset%natom+10) then 1494 ! 2nd order (k,k) 1495 if (idir<=3) then ! one direction needed 1496 nwffile=1 1497 file_index(1)=idir+dtset%natom*3 1498 fnamewff(1)=dtfil%fnamewffddk 1499 else ! two directions needed 1500 nwffile=2 1501 file_index(1)=idir1+dtset%natom*3 1502 file_index(2)=idir2+dtset%natom*3 1503 fnamewff(1)=dtfil%fnamewffddk 1504 fnamewff(2)=dtfil%fnamewffddk 1505 end if 1506 else if (ipert==dtset%natom+11) then 1507 ! 2nd order (k,E) 1508 nwffile=3 ! dk, dE and dkdk 1509 idir_dkdk = idir 1510 if(idir_dkdk>6) idir_dkdk = idir_dkdk - 3 1511 file_index(1)=idir_dkdk +(dtset%natom+6)*3 ! dkdk 1512 file_index(2)=idir2+(dtset%natom+1)*3 ! defld file (dir2) 1513 file_index(3)=idir1+dtset%natom*3 ! ddk file (dir1) 1514 fnamewff(1)=dtfil%fnamewffdkdk 1515 fnamewff(2)=dtfil%fnamewffdelfd 1516 fnamewff(3)=dtfil%fnamewffddk 1517 if (idir>3) then 1518 nwffile=4 1519 file_index(4)=idir2+dtset%natom*3 ! ddk file (dir2) 1520 fnamewff(4)=dtfil%fnamewffddk 1521 end if 1522 end if 1523 do ii=1,nwffile 1524 call appdig(file_index(ii),fnamewff(ii),fiwfddk) 1525 ! Checking the existence of data file 1526 if (.not. file_exists(fiwfddk)) then 1527 ! Trick needed to run Abinit test suite in netcdf mode. 1528 if (file_exists(nctk_ncify(fiwfddk))) then 1529 write(message,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name." 1530 call wrtout(std_out,message,'COLL') 1531 fiwfddk = nctk_ncify(fiwfddk) 1532 end if 1533 if (.not. file_exists(fiwfddk)) then 1534 MSG_ERROR('Missing file: '//TRIM(fiwfddk)) 1535 end if 1536 end if 1537 write(message,'(2a)')'-dfpt_looppert : read the wavefunctions from file: ',trim(fiwfddk) 1538 call wrtout(std_out,message,'COLL') 1539 call wrtout(ab_out,message,'COLL') 1540 ! Note that the unit number for these files is 50,51,52 or 53 (dtfil%unddk=50) 1541 call wfk_open_read(ddk_f(ii),fiwfddk,formeig1,dtset%iomode,dtfil%unddk+(ii-1),spaceComm) 1542 end do 1543 end if 1544 1545 ! Get first-order local potentials and 1st-order core correction density change 1546 ! (do NOT include xccc3d1 in vpsp1 : this will be done in dfpt_scfcv because vpsp1 1547 ! might become spin-polarized) 1548 1549 n3xccc=0;if(psps%n1xccc/=0)n3xccc=nfftf 1550 ABI_ALLOCATE(xccc3d1,(cplex*n3xccc)) 1551 ABI_ALLOCATE(vpsp1,(cplex*nfftf)) 1552 1553 ! PAW: compute Vloc(1) and core(1) together in reciprocal space 1554 ! -------------------------------------------------------------- 1555 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 1556 ndir=1 1557 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,istr,ipert,& 1558 & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,& 1559 & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 1560 & atmrhor1=xccc3d1,atmvlocr1=vpsp1,optn_in=n3xccc/nfftf,optn2_in=1,vspl=psps%vlspl) 1561 else 1562 1563 ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp. 1564 ! ------------------------------------------------------------------------------ 1565 1566 if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then 1567 ! Section for strain perturbation 1568 call vlocalstr(gmet,gprimd,gsqcut,istr,mgfftf,mpi_enreg,& 1569 & psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,ntypat,dtset%paral_kgb,ph1df,psps%qgrid_vl,& 1570 & ucvol,psps%vlspl,vpsp1) 1571 else 1572 call dfpt_vlocal(atindx,cplex,gmet,gsqcut,idir,ipert,mpi_enreg,psps%mqgrid_vl,dtset%natom,& 1573 & nattyp,nfftf,ngfftf,ntypat,ngfftf(1),ngfftf(2),ngfftf(3),dtset%paral_kgb,ph1df,psps%qgrid_vl,& 1574 & dtset%qptn,ucvol,psps%vlspl,vpsp1,xred) 1575 !SPr: need vpsp1 for -q as well, but for magnetic field it's zero, to be done later 1576 end if 1577 1578 if(psps%n1xccc/=0)then 1579 call dfpt_mkcore(cplex,idir,ipert,dtset%natom,ntypat,ngfftf(1),psps%n1xccc,& 1580 & ngfftf(2),ngfftf(3),dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred) 1581 !SPr: same here, need xccc3d1 for -q as well for phonon pert.. to be done later 1582 end if ! psps%n1xccc/=0 1583 end if ! usepaw 1584 1585 eigen1(:)=zero; resid(:)=zero 1586 if(.not.kramers_deg) then 1587 eigen1_mq(:)=zero 1588 resid_mq(:)=zero 1589 end if 1590 ! Get starting charge density and Hartree + xc potential 1591 ABI_ALLOCATE(rhor1,(cplex*nfftf,nspden)) 1592 ABI_ALLOCATE(rhog1,(2,nfftf)) 1593 1594 if(.not.kramers_deg) then 1595 !Case when first order spinors at both +q and -q are not related by symmetry (time and/or space inversion) 1596 ABI_ALLOCATE(rhor1_pq,(cplex*nfftf,nspden)) 1597 ABI_ALLOCATE(rhog1_pq,(2,nfftf)) 1598 1599 ABI_ALLOCATE(rhor1_mq,(cplex*nfftf,nspden)) 1600 ABI_ALLOCATE(rhog1_mq,(2,nfftf)) 1601 end if 1602 1603 ! can we get this set of gkk matrices from previously calculated rhog1 through a non-scf calculation? 1604 found_eq_gkk=.false. 1605 if (dtset%prepgkk == 1 .and. ipert <= dtset%natom) then 1606 ! NOTE: this does not take into account combinations e.g. x+y -> z 1607 ! if rhor1 add linearly this could be done... 1608 do icase_eq = 1, icase-1 1609 idir_eq = mod(icase_eq,3) 1610 if (idir_eq==0) idir_eq=3 1611 ipert_eq = ( (icase_eq-idir_eq) / 3 + 1) 1612 1613 ! find sym which links old perturbation to present one 1614 do isym=1, nsym 1615 ! check that isym preserves qpt to begin with 1616 if (symq(4,1,isym) /= 1 .or. & 1617 & symq(1,1,isym) /= 0 .or. & 1618 & symq(2,1,isym) /= 0 .or. & 1619 & symq(3,1,isym) /= 0 ) cycle 1620 1621 eq_symop = dtset%symrel(:,:,isym) 1622 if (indsym(4,isym,ipert) == ipert_eq .and. & 1623 & abs(eq_symop(idir,idir_eq)) == 1 .and. & 1624 & sum(abs(eq_symop(:,idir_eq))) == 1) then 1625 found_eq_gkk = .true. 1626 exit 1627 end if 1628 end do ! isym 1629 if (found_eq_gkk) exit 1630 end do ! icase_eq 1631 end if ! check for prepgkk with symmetric pert 1632 1633 if (found_eq_gkk) then 1634 write (message, '(a,l6,i6,2a,3i6,2a,3i6,a)') & 1635 & ' found_eq_gkk,isym = ', found_eq_gkk, isym, ch10, & 1636 & ' idir, ipert, icase = ', idir, ipert, icase, ch10, & 1637 & ' idireq iperteq icaseeq = ', idir_eq, ipert_eq, icase_eq, ch10 1638 call wrtout(std_out,message,'COLL') 1639 ! 1640 ! Make density for present perturbation, which is symmetric of 1 or more previous perturbations: 1641 ! rotate 1DEN arrays with symrel(isym) to produce rhog1_eq rhor1_eq 1642 ! 1643 if (dtset%use_nonscf_gkk == 1) then 1644 call rotate_rho(cplex, timrev_pert, mpi_enreg, nfftf, ngfftf, nspden, & 1645 & rhor1_save(:,:,icase_eq), rhog1, rhor1, eq_symop, dtset%tnons(:,isym)) 1646 1647 rhor1 = rhor1 * eq_symop(idir,idir_eq) 1648 1649 ! TODO: rotate rhoij in PAW case 1650 1651 blkflg_save = blkflg 1652 dtset_tmp%iscf = -2 1653 iscf_mod = -2 1654 dtset_tmp%nstep = 1 1655 dtset_tmp%nline = 1 1656 if (abs(dtset_tmp%tolwfr) < 1.e-24) dtset_tmp%tolwfr = 1.e-24 1657 dtset_tmp%toldfe = zero 1658 dtset_tmp%toldff = zero 1659 dtset_tmp%tolrff = zero 1660 dtset_tmp%tolvrs = zero 1661 write (message, '(a,i6,a)') ' NOTE: doing GKK calculation for icase ', icase, ' with non-SCF calculation' 1662 call wrtout(std_out,message,'COLL') 1663 !call wrtout(ab_out,message,'COLL') ! decomment and update output files 1664 1665 else ! do not use non-scf shortcut, but save rotated 1DEN for comparison 1666 ! saves the rotated rho, for later comparison with the full SCF rhor1: comment lines below for iscf = -2 1667 call rotate_rho(cplex, timrev_pert, mpi_enreg, nfftf, ngfftf, nspden, & 1668 & rhor1_save(:,:,icase_eq), rhog1, rhor1_save(:,:,icase), eq_symop, & 1669 & dtset%tnons(:,isym)) 1670 rhor1_save(:,:,icase) = rhor1_save(:,:,icase) * eq_symop(idir,idir_eq) 1671 1672 end if ! force non scf calculation of other gkk, or not 1673 end if ! found an equiv perturbation for the gkk 1674 1675 if ( (dtfil%ireadwf==0 .and. iscf_mod/=-4 .and. dtset%get1den==0 .and. dtset%ird1den==0) & 1676 .or. (iscf_mod== -3 .and. ipert/=dtset%natom+11) ) then 1677 ! NOTE : For ipert==natom+11, we want to read the 1st order density from a previous calculation 1678 rhor1(:,:)=zero ; rhog1(:,:)=zero 1679 ! PAW: rhoij have been set to zero in call to pawrhoij_alloc above 1680 1681 init_rhor1 = ((ipert>=1 .and. ipert<=dtset%natom).or.ipert==dtset%natom+5) 1682 ! This section is needed in order to maintain the old behavior and pass the automatic tests 1683 if (psps%usepaw == 0) then 1684 init_rhor1 = init_rhor1 .and. all(psps%nctab(:ntypat)%has_tvale) 1685 else 1686 init_rhor1 = .False. 1687 end if 1688 1689 if (init_rhor1) then 1690 1691 if(ipert/=dtset%natom+5) then 1692 1693 ! Initialize rhor1 and rhog1 from the derivative of atomic densities/gaussians. 1694 ndir = 1; optn2 = 3 1695 if (psps%usepaw==1) then 1696 ! FIXME: Here there's a bug because has_tvale == 0 or 1 instead of 2 1697 if (all(pawtab(:ntypat)%has_tvale/=0)) optn2=2 1698 else if (psps%usepaw==0) then 1699 if (all(psps%nctab(:ntypat)%has_tvale)) optn2=2 1700 end if 1701 1702 if (optn2 == 3) then 1703 call wrtout(std_out," Initializing rhor1 from atom-centered gaussians", "COLL") 1704 ABI_ALLOCATE(gauss,(2,ntypat)) 1705 call atom_gauss(ntypat, dtset%densty, psps%ziontypat, psps%znucltypat, gauss) 1706 1707 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,& 1708 mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,& 1709 ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 1710 atmrhor1=rhor1,optn_in=1,optn2_in=3,gauss=gauss) 1711 1712 ABI_FREE(gauss) 1713 else 1714 call wrtout(std_out," Initializing rhor1 from valence densities taken from pseudopotential files", "COLL") 1715 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,& 1716 mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,& 1717 ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,& 1718 atmrhor1=rhor1,optn_in=1,optn2_in=2) 1719 end if 1720 1721 else 1722 !Magnetic field perturbation 1723 call wrtout(std_out," Initializing rhor1 guess based on the ground state XC magnetic field", "COLL") 1724 1725 call dfpt_init_mag1(idir,rhor1,rhor,cplex,nfftf,nspden,vxc,kxc,nkxc) 1726 1727 if(.not.kramers_deg) then 1728 rhor1_pq=rhor1 1729 call dfpt_init_mag1(idir,rhor1_mq,rhor,cplex,nfftf,nspden,vxc,kxc,nkxc) 1730 end if 1731 end if 1732 1733 call fourdp(cplex,rhog1,rhor1,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 1734 if (.not.kramers_deg) then 1735 !call fourdp(cplex,rhog1_pq,rhor1_pq,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 1736 rhog1_pq=rhog1 1737 call fourdp(cplex,rhog1_mq,rhor1_mq,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 1738 end if 1739 end if 1740 1741 else ! rhor1 not being forced to 0.0 1742 if(iscf_mod>0) then 1743 ! cplex=2 gets the complex density, =1 only real part 1744 if (psps%usepaw==1) then 1745 ! Be careful: in PAW, rho does not include the 1st-order compensation 1746 ! density (to be added in dfpt_scfcv.F90) ! 1747 ABI_ALLOCATE(rho1wfg,(2,dtset%nfft)) 1748 ABI_ALLOCATE(rho1wfr,(dtset%nfft,nspden)) 1749 call dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon1,istwfk_rbz,& 1750 & kg,kg1,dtset%mband,dtset%mgfft,mkmem_rbz,mk1mem_rbz,mpi_enreg,mpw,mpw1,nband_rbz,& 1751 & dtset%nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,nsym1,& 1752 & occ_rbz,dtset%paral_kgb,phnons1,rho1wfg,rho1wfr,rprimd,symaf1,symrl1,ucvol,& 1753 & wtk_rbz) 1754 call transgrid(cplex,mpi_enreg,nspden,+1,1,1,dtset%paral_kgb,pawfgr,rho1wfg,rhog1,rho1wfr,rhor1) 1755 ABI_DEALLOCATE(rho1wfg) 1756 ABI_DEALLOCATE(rho1wfr) 1757 else 1758 !SPr: need to modify dfpt_mkrho to taken into account q,-q and set proper formulas when +q and -q spinors are 1759 ! related 1760 call dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon1,istwfk_rbz,& 1761 & kg,kg1,dtset%mband,dtset%mgfft,mkmem_rbz,mk1mem_rbz,mpi_enreg,mpw,mpw1,nband_rbz,& 1762 & dtset%nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,nsym1,& 1763 & occ_rbz,dtset%paral_kgb,phnons1,rhog1,rhor1,rprimd,symaf1,symrl1,ucvol,& 1764 & wtk_rbz) 1765 end if 1766 1767 else if (.not. found_eq_gkk) then 1768 ! negative iscf_mod and no symmetric rotation of rhor1 1769 ! Read rho1(r) from a disk file and broadcast data. 1770 rdwr=1;rdwrpaw=psps%usepaw;if(dtfil%ireadwf/=0) rdwrpaw=0 1771 if (ipert/=dtset%natom+11) then 1772 call appdig(pertcase,dtfil%fildens1in,fiden1i) 1773 else ! For ipert==natom+11, we want to read the 1st order density from a previous calculation 1774 call appdig(idir2+(dtset%natom+1)*3,dtfil%fildens1in,fiden1i) 1775 end if 1776 ! call appdig(pertcase,dtfil%fildens1in,fiden1i) 1777 call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor1, & 1778 hdr_den, pawrhoij1, spaceComm, check_hdr=hdr) 1779 etotal = hdr_den%etot; call hdr_free(hdr_den) 1780 1781 ! Compute up+down rho1(G) by fft 1782 ABI_ALLOCATE(work,(cplex*nfftf)) 1783 work(:)=rhor1(:,1) 1784 call fourdp(cplex,rhog1,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 1785 ABI_DEALLOCATE(work) 1786 end if ! rhor1 generated or read in from file 1787 1788 end if ! rhor1 set to 0 or read in from file 1789 1790 ! Check whether exiting was required by the user. 1791 ! If found then do not start minimization steps 1792 openexit=1 ; if(dtset%chkexit==0) openexit=0 1793 call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit) 1794 ! If immediate exit, and wavefunctions were not read, must zero eigenvalues 1795 if (iexit/=0) eigen1(:)=zero 1796 if (iexit/=0.and.(.not.kramers_deg)) eigen1_mq(:)=zero 1797 1798 if (iexit==0) then 1799 _IBM6("before dfpt_scfcv") 1800 1801 ! Main calculation: get 1st-order wavefunctions from Sternheimer equation (SCF cycle) 1802 ! if ipert==natom+10 or natom+11 : get 2nd-order wavefunctions 1803 if (kramers_deg) then 1804 call dfpt_scfcv(atindx,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,& 1805 & dielt,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset_tmp,& 1806 & d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,& 1807 & ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,eloc0,elpsp1,& 1808 & enl0,enl1,eovl1,epaw1,etotal,evdw,exc1,fermie,gh0c1_set,gh1c_set,hdr,idir,& 1809 & indkpt1,indsy1,initialized,ipert,irrzon1,istwfk_rbz,& 1810 & kg,kg1,kpt_rbz,kxc,mgfftf,mkmem_rbz,mkqmem_rbz,mk1mem_rbz,& 1811 & mpert,mpi_enreg,mpw,mpw1,mpw1_mq,my_natom,& 1812 & nattyp,nband_rbz,ncpgr,nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,& 1813 & npwarr,npwar1,nspden,& 1814 & nsym1,n3xccc,occkq,occ_rbz,& 1815 & paw_an_pert,paw_ij_pert,pawang,pawang1,pawfgr,pawfgrtab_pert,pawrad,pawrhoij_pert,pawrhoij1,pawtab,& 1816 & pertcase,phnons1,ph1d,ph1df,prtbbb,psps,& 1817 & dtset%qptn,resid,residm,rhog,rhog1,& 1818 & rhor,rhor1,rprimd,symaf1,symrc1,symrl1,& 1819 & usecprj,useylmgr,useylmgr1,ddk_f,vpsp1,vtrial,vxc,& 1820 & wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,zeff,dfpt_scfcv_retcode,& 1821 & kramers_deg) 1822 else 1823 call dfpt_scfcv(atindx,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,& 1824 & dielt,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset_tmp,& 1825 & d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,& 1826 & ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,eloc0,elpsp1,& 1827 & enl0,enl1,eovl1,epaw1,etotal,evdw,exc1,fermie,gh0c1_set,gh1c_set,hdr,idir,& 1828 & indkpt1,indsy1,initialized,ipert,irrzon1,istwfk_rbz,& 1829 & kg,kg1,kpt_rbz,kxc,mgfftf,mkmem_rbz,mkqmem_rbz,mk1mem_rbz,& 1830 & mpert,mpi_enreg,mpw,mpw1,mpw1_mq,my_natom,& 1831 & nattyp,nband_rbz,ncpgr,nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,& 1832 & npwarr,npwar1,nspden,& 1833 & nsym1,n3xccc,occkq,occ_rbz,& 1834 & paw_an_pert,paw_ij_pert,pawang,pawang1,pawfgr,pawfgrtab_pert,pawrad,pawrhoij_pert,pawrhoij1,pawtab,& 1835 & pertcase,phnons1,ph1d,ph1df,prtbbb,psps,& 1836 & dtset%qptn,resid,residm,rhog,rhog1,& 1837 & rhor,rhor1,rprimd,symaf1,symrc1,symrl1,& 1838 & usecprj,useylmgr,useylmgr1,ddk_f,vpsp1,vtrial,vxc,& 1839 & wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,zeff,dfpt_scfcv_retcode,& 1840 & kramers_deg,& 1841 & cg_mq=cg_mq,cg1_mq=cg1_mq,cg1_active_mq=cg1_active_mq,docckde_mq=docckde_mq,eigen_mq=eigen_mq,& 1842 & eigen1_mq=eigen1_mq,gh0c1_set_mq=gh0c1_set_mq,gh1c_set_mq=gh1c_set_mq,& 1843 & kg1_mq=kg1_mq,npwar1_mq=npwar1_mq,occk_mq=occk_mq,resid_mq=resid_mq,residm_mq=residm_mq,& 1844 & rhog1_pq=rhog1_pq,rhog1_mq=rhog1_mq,rhor1_pq=rhor1_pq,rhor1_mq=rhor1_mq) 1845 end if 1846 1847 _IBM6("after dfpt_scfcv") 1848 1849 ! 2nd-order eigenvalues stuff 1850 if (dtset%ieig2rf>0) then 1851 if (first_entry) then 1852 nullify(eigen1_pert) 1853 first_entry = .false. 1854 end if 1855 if (.not.associated(eigen1_pert)) then 1856 ABI_ALLOCATE(eigen1_pert,(2*dtset%mband**2*nkpt*dtset%nsppol,3,mpert)) 1857 ABI_STAT_ALLOCATE(cg1_pert,(2,mpw1*nspinor*dtset%mband*mk1mem_rbz*nsppol*dim_eig2rf,3,mpert),ierr) 1858 ABI_CHECK(ierr==0, "out-of-memory in cg1_pert") 1859 ABI_ALLOCATE(gh0c1_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert)) 1860 ABI_ALLOCATE(gh1c_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert)) 1861 ABI_ALLOCATE(kpt_rbz_pert,(3,nkpt_rbz)) 1862 ABI_ALLOCATE(npwarr_pert,(nkpt_rbz,mpert)) 1863 ABI_ALLOCATE(npwar1_pert,(nkpt_rbz,mpert)) 1864 ABI_ALLOCATE(npwtot_pert,(nkpt_rbz,mpert)) 1865 eigen1_pert(:,:,:) = zero 1866 cg1_pert(:,:,:,:) = zero 1867 gh0c1_pert(:,:,:,:) = zero 1868 gh1c_pert(:,:,:,:) = zero 1869 npwar1_pert (:,:) = 0 1870 npwarr_pert (:,:) = 0 1871 kpt_rbz_pert = kpt_rbz 1872 end if 1873 clflg(idir,ipert)=1 1874 eigen1_pert(1:2*dtset%mband**2*nkpt_rbz*dtset%nsppol,idir,ipert) = eigen1(:) 1875 if(dtset%ieig2rf==1.or.dtset%ieig2rf==3.or.dtset%ieig2rf==4.or.dtset%ieig2rf==5) then 1876 cg1_pert(:,:,idir,ipert)=cg1_active(:,:) 1877 gh0c1_pert(:,:,idir,ipert)=gh0c1_set(:,:) 1878 gh1c_pert(:,:,idir,ipert)=gh1c_set(:,:) 1879 end if 1880 npwarr_pert(:,ipert)=npwarr(:) 1881 npwar1_pert(:,ipert)=npwar1(:) 1882 npwtot_pert(:,ipert)=npwtot(:) 1883 end if 1884 ! 2nd-order eigenvalues stuff for EFMAS 1885 if (dtset%efmas>0) then 1886 if (first_entry) then 1887 nullify(eigen1_pert) 1888 first_entry = .false. 1889 end if 1890 if (.not.associated(eigen1_pert)) then 1891 ABI_ALLOCATE(eigen1_pert,(2*dtset%mband**2*nkpt*dtset%nsppol,3,mpert)) 1892 ABI_ALLOCATE(cg1_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert)) 1893 ABI_ALLOCATE(gh0c1_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert)) 1894 ABI_ALLOCATE(gh1c_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert)) 1895 ABI_ALLOCATE(kpt_rbz_pert,(3,nkpt_rbz)) 1896 ABI_ALLOCATE(npwarr_pert,(nkpt_rbz,mpert)) 1897 eigen1_pert(:,:,:) = zero 1898 cg1_pert(:,:,:,:) = zero 1899 gh0c1_pert(:,:,:,:) = zero 1900 gh1c_pert(:,:,:,:) = zero 1901 npwarr_pert (:,:) = 0 1902 kpt_rbz_pert = kpt_rbz 1903 ABI_ALLOCATE(cg0_pert,(2,mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol*dim_eig2rf)) 1904 cg0_pert = cg 1905 end if 1906 eigen1_pert(1:2*dtset%mband**2*nkpt_rbz*dtset%nsppol,idir,ipert) = eigen1(:) 1907 cg1_pert(:,:,idir,ipert)=cg1_active(:,:) 1908 gh0c1_pert(:,:,idir,ipert)=gh0c1_set(:,:) 1909 gh1c_pert(:,:,idir,ipert)=gh1c_set(:,:) 1910 npwarr_pert(:,ipert)=npwarr(:) 1911 end if 1912 ABI_DEALLOCATE(gh1c_set) 1913 ABI_DEALLOCATE(gh0c1_set) 1914 ABI_DEALLOCATE(cg1_active) 1915 if(.not.kramers_deg) then 1916 ABI_DEALLOCATE(gh1c_set_mq) 1917 ABI_DEALLOCATE(gh0c1_set_mq) 1918 ABI_DEALLOCATE(cg1_active_mq) 1919 end if 1920 1921 end if ! End of the check of hasty exit 1922 1923 call timab(146,1,tsec) 1924 1925 ! Print out message at the end of the iterations 1926 write(message, '(80a,a,a,a,a)' ) ('=',ii=1,80),ch10,ch10,& 1927 & ' ----iterations are completed or convergence reached----',ch10 1928 call wrtout(ab_out,message,'COLL') 1929 call wrtout(std_out, message,'COLL') 1930 1931 ! Print _gkk file for this perturbation 1932 if (dtset%prtgkk == 1) then 1933 call appdig(3*(ipert-1)+idir,dtfil%fnameabo_gkk,gkkfilnam) 1934 nmatel = dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol 1935 ABI_ALLOCATE(phasecg, (2, nmatel)) 1936 call getcgqphase(dtset, timrev, cg, mcg, cgq, mcgq, mpi_enreg, nkpt_rbz, npwarr, npwar1, phasecg) 1937 phasecg(1,:) = one 1938 phasecg(2,:) = zero 1939 ! NB: phasecg not actually used in outgkk for the moment (2013/08/15) 1940 call outgkk(bantot_rbz, nmatel,gkkfilnam,eigen0,eigen1,hdr0,hdr,mpi_enreg,phasecg) 1941 ABI_DEALLOCATE(phasecg) 1942 1943 #ifdef HAVE_NETCDF 1944 ! Reshape eigen1 into gkk for netCDF output 1945 ABI_STAT_ALLOCATE(gkk,(2*dtset%mband*dtset%nsppol,dtset%nkpt,1,1,dtset%mband), ierr) 1946 ABI_CHECK(ierr==0, "out-of-memory in gkk") 1947 gkk(:,:,:,:,:) = zero 1948 mband = dtset%mband 1949 band_index = 0 1950 band2tot_index = 0 1951 do isppol=1,dtset%nsppol 1952 do ikpt =1,nkpt_rbz 1953 do iband=1,dtset%mband 1954 do jband=1,dtset%mband 1955 eig1_r = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index) 1956 eig1_i = eigen1(2*jband+(iband-1)*2*mband+band2tot_index) 1957 gkk(2*iband-1+2*band_index,ikpt,1,1,jband) = & 1958 & gkk(2*iband-1+2*band_index,ikpt,1,1,jband) + eig1_r 1959 gkk(2*iband+2*band_index,ikpt,1,1,jband) = & 1960 & gkk(2*iband+2*band_index,ikpt,1,1,jband) + eig1_i 1961 end do !jband 1962 end do !iband 1963 band2tot_index = band2tot_index + 2*mband**2 1964 end do !ikpt 1965 band_index = band_index + mband 1966 end do !isppol 1967 1968 ! Initialize ggk_ebands to write in the GKK.nc file 1969 ! MG FIXME: Here there's a bug because eigen0 is dimensioned with nkpt_rbz i.e. IBZ(q) 1970 ! but the ebands_t object is constructed with dimensions taken from hdr0 i.e. the IBZ(q=0). 1971 bantot= dtset%mband*dtset%nkpt*dtset%nsppol 1972 call ebands_init(bantot,gkk_ebands,dtset%nelect,doccde,eigen0,hdr0%istwfk,hdr0%kptns,& 1973 & hdr0%nband, hdr0%nkpt,hdr0%npwarr,hdr0%nsppol,hdr0%nspinor,& 1974 & hdr0%tphysel,hdr0%tsmear,hdr0%occopt,hdr0%occ,hdr0%wtk,& 1975 & hdr0%charge, hdr0%kptopt, hdr0%kptrlatt_orig, hdr0%nshiftk_orig, hdr0%shiftk_orig, & 1976 & hdr0%kptrlatt, hdr0%nshiftk, hdr0%shiftk) 1977 1978 ! Init a gkk_t object 1979 call gkk_init(gkk,gkk2d,dtset%mband,dtset%nsppol,nkpt_rbz,1,1) 1980 1981 ! Write the netCDF file. 1982 fname = strcat(gkkfilnam,".nc") 1983 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKK file") 1984 NCF_CHECK(crystal_ncwrite(crystal, ncid)) 1985 NCF_CHECK(ebands_ncwrite(gkk_ebands, ncid)) 1986 call gkk_ncwrite(gkk2d,dtset%qptn(:),dtset%wtq, ncid) 1987 NCF_CHECK(nf90_close(ncid)) 1988 1989 ! Free memory 1990 ABI_DEALLOCATE(gkk) 1991 call gkk_free(gkk2d) 1992 call ebands_free(gkk_ebands) 1993 #endif 1994 end if 1995 1996 if (dtset%prepgkk == 1 .and. found_eq_gkk) then 1997 if (dtset%use_nonscf_gkk == 1) then 1998 ! Restore old values of SCF cycle parameters 1999 iscf_mod = iscf_mod_save 2000 dtset_tmp%iscf = iscf_mod_save 2001 dtset_tmp%nstep = nstep_save 2002 dtset_tmp%nline = nline_save 2003 dtset_tmp%tolwfr = tolwfr_save 2004 dtset_tmp%toldfe = toldfe_save 2005 dtset_tmp%toldff = toldff_save 2006 dtset_tmp%tolrff = tolrff_save 2007 dtset_tmp%tolvrs = tolvrs_save 2008 blkflg = blkflg_save ! this ensures we do not use the (unconverged) 2DTE from this non scf run 2009 ! Save density for present perturbation, for future use in symmetric perturbations 2010 rhor1_save(:,:,icase) = rhor1 2011 else 2012 write (message, '(a,3E20.10)') 'norm diff = ', sum(abs(rhor1_save(:,:,icase) - rhor1)), & 2013 & sum(abs(rhor1)), sum(abs(rhor1_save(:,:,icase) - rhor1))/sum(abs(rhor1)) 2014 call wrtout(std_out, message,'COLL') 2015 end if 2016 end if 2017 2018 ! Write wavefunctions file only if convergence was not achieved. 2019 write_1wfk = .True. 2020 if (dtset%prtwf==-1 .and. dfpt_scfcv_retcode == 0) then 2021 write_1wfk = .False. 2022 call wrtout(ab_out," dfpt_looppert: DFPT cycle converged with prtwf=-1. Will skip output of the 1st-order WFK file.","COLL") 2023 end if 2024 2025 if (write_1wfk) then 2026 ! Output 1st-order wavefunctions in file 2027 call status(pertcase,dtfil%filstat,iexit,level,'call outwf ') 2028 call outwf(cg1,dtset,psps,eigen1,fiwf1o,hdr,kg1,kpt_rbz,& 2029 & dtset%mband,mcg1,mk1mem_rbz,mpi_enreg,mpw1,dtset%natom,nband_rbz,& 2030 & nkpt_rbz,npwar1,dtset%nsppol,& 2031 & occ_rbz,resid,response,dtfil%unwff2,wvl%wfs,wvl%descr) 2032 end if 2033 2034 #ifdef HAVE_NETCDF 2035 ! Output DDK file in netcdf format. 2036 if (me == master .and. ipert == dtset%natom + 1) then 2037 fname = strcat(dtfil%filnam_ds(4), "_DDK.nc") 2038 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating DDK.nc file") 2039 ! Have to build hdr on k-grid with info about perturbation. 2040 call hdr_copy(hdr0, hdr_tmp) 2041 hdr_tmp%kptopt = dtset%kptopt 2042 hdr_tmp%pertcase = pertcase 2043 hdr_tmp%qptn = dtset%qptn(1:3) 2044 NCF_CHECK(hdr_ncwrite(hdr_tmp, ncid, 43, nc_define=.True.)) 2045 call hdr_free(hdr_tmp) 2046 NCF_CHECK(crystal_ncwrite(crystal, ncid)) 2047 NCF_CHECK(ebands_ncwrite(ebands_k, ncid)) 2048 ncerr = nctk_def_arrays(ncid, [ & 2049 nctkarr_t('h1_matrix_elements', "dp", & 2050 "two, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.) 2051 NCF_CHECK(ncerr) 2052 NCF_CHECK(nctk_set_datamode(ncid)) 2053 ncerr = nf90_put_var(ncid, nctk_idname(ncid, "h1_matrix_elements"), eigen1, & 2054 count=[2, dtset%mband, dtset%mband, nkpt_rbz, dtset%nsppol]) 2055 NCF_CHECK(ncerr) 2056 NCF_CHECK(nf90_close(ncid)) 2057 end if 2058 #endif 2059 2060 2061 ! If the perturbation is d/dk, evaluate the f-sum rule. 2062 if (ipert==dtset%natom+1 )then 2063 ! Note : the factor of two is related to the difference 2064 ! between Taylor expansion and perturbation expansion 2065 ! Note : this expression should be modified for ecutsm. 2066 ! Indeed, the present one will NOT tend to 1.0_dp. 2067 ek2=gmet(idir,idir)*(two_pi**2)*2.0_dp*dtset%nelect 2068 fsum=-ek1/ek2 2069 if(dtset%ecutsm<tol6)then 2070 write(message, '(a,es20.10,a,a,es20.10)' ) & 2071 & ' dfpt_looppert : ek2=',ek2,ch10,& 2072 & ' f-sum rule ratio=',fsum 2073 else 2074 write(message, '(a,es20.10,a,a,es20.10,a)' ) & 2075 & ' dfpt_looppert : ek2=',ek2,ch10,& 2076 & ' f-sum rule ratio=',fsum,' (note : ecutsm/=0)' 2077 end if 2078 call wrtout(std_out,message,'COLL') 2079 call wrtout(ab_out,message,'COLL') 2080 ! Write the diagonal elements of the dH/dk operator, after averaging over degenerate states 2081 ABI_ALLOCATE(eigen1_mean,(dtset%mband*nkpt_rbz*dtset%nsppol)) 2082 call eigen_meandege(eigen0,eigen1,eigen1_mean,dtset%mband,nband_rbz,nkpt_rbz,dtset%nsppol,1) 2083 option=4 2084 call prteigrs(eigen1_mean,dtset%enunit,fermie,dtfil%fnametmp_1wf1_eig,ab_out,iscf_mod,kpt_rbz,dtset%kptopt,& 2085 & dtset%mband,nband_rbz,nkpt_rbz,dtset%nnsclo,dtset%nsppol,occ_rbz,dtset%occopt,& 2086 & option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk_rbz) 2087 ABI_DEALLOCATE(eigen1_mean) 2088 end if 2089 2090 ! Print the energies 2091 if (dtset%nline/=0 .or. dtset%nstep/=0)then 2092 call dfpt_prtene(dtset%berryopt,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,& 2093 & ehart01,ehart1,eii,ek0,ek1,eloc0,elpsp1,enl0,enl1,eovl1,epaw1,evdw,exc1,ab_out,& 2094 & ipert,dtset%natom,psps%usepaw,usevdw) 2095 end if 2096 2097 if(mpi_enreg%paral_pert==1) then 2098 if (ipert_me < npert_me -1) then 2099 call hdr_free(hdr0) 2100 else 2101 eigen0_copy(1:dtset%mband*nkpt_rbz*dtset%nsppol) = eigen0 2102 end if 2103 ipert_me = ipert_me +1 2104 else 2105 if (icase == ipert_cnt) then 2106 eigen0_copy(1:dtset%mband*nkpt_rbz*dtset%nsppol) = eigen0 2107 else 2108 call hdr_free(hdr0) 2109 end if 2110 end if 2111 2112 !Release the temporary arrays (for k, k+q and 1st-order) 2113 ABI_DEALLOCATE(cg) 2114 ABI_DEALLOCATE(cgq) 2115 ABI_DEALLOCATE(cg1) 2116 ABI_DEALLOCATE(docckqde) 2117 if(.not.kramers_deg) then 2118 ABI_DEALLOCATE(cg_mq) 2119 ABI_DEALLOCATE(cg1_mq) 2120 ABI_DEALLOCATE(docckde_mq) 2121 end if 2122 ABI_DEALLOCATE(doccde_rbz) 2123 ABI_DEALLOCATE(eigen0) 2124 ABI_DEALLOCATE(eigenq) 2125 ABI_DEALLOCATE(eigen1) 2126 ABI_DEALLOCATE(kpq) 2127 if(.not.kramers_deg) then 2128 ABI_DEALLOCATE(eigen_mq) 2129 ABI_DEALLOCATE(eigen1_mq) 2130 ABI_DEALLOCATE(kmq) 2131 end if 2132 ABI_DEALLOCATE(indkpt1) 2133 ABI_DEALLOCATE(indsy1) 2134 ABI_DEALLOCATE(istwfk_rbz) 2135 ABI_DEALLOCATE(irrzon1) 2136 ABI_DEALLOCATE(kg) 2137 ABI_DEALLOCATE(kg1) 2138 ABI_DEALLOCATE(kpq_rbz) 2139 if(.not.kramers_deg) then 2140 ABI_DEALLOCATE(kg1_mq) 2141 ABI_DEALLOCATE(kmq_rbz) 2142 end if 2143 ABI_DEALLOCATE(kpt_rbz) 2144 ABI_DEALLOCATE(nband_rbz) 2145 ABI_DEALLOCATE(npwarr) 2146 ABI_DEALLOCATE(npwar1) 2147 ABI_DEALLOCATE(npwtot) 2148 ABI_DEALLOCATE(npwtot1) 2149 ABI_DEALLOCATE(occkq) 2150 ABI_DEALLOCATE(occ_rbz) 2151 ABI_DEALLOCATE(phnons1) 2152 ABI_DEALLOCATE(resid) 2153 ABI_DEALLOCATE(rhog1) 2154 ABI_DEALLOCATE(rhor1) 2155 ABI_DEALLOCATE(symaf1) 2156 ABI_DEALLOCATE(symrc1) 2157 ABI_DEALLOCATE(symrl1) 2158 ABI_DEALLOCATE(tnons1) 2159 ABI_DEALLOCATE(wtk_rbz) 2160 ABI_DEALLOCATE(xccc3d1) 2161 ABI_DEALLOCATE(vpsp1) 2162 ABI_DEALLOCATE(ylm) 2163 ABI_DEALLOCATE(ylm1) 2164 ABI_DEALLOCATE(ylmgr) 2165 ABI_DEALLOCATE(ylmgr1) 2166 if(.not.kramers_deg) then 2167 ABI_DEALLOCATE(npwar1_mq) 2168 ABI_DEALLOCATE(npwtot1_mq) 2169 ABI_DEALLOCATE(occk_mq) 2170 ABI_DEALLOCATE(resid_mq) 2171 ABI_DEALLOCATE(rhor1_pq) 2172 ABI_DEALLOCATE(rhor1_mq) 2173 ABI_DEALLOCATE(rhog1_pq) 2174 ABI_DEALLOCATE(rhog1_mq) 2175 end if 2176 if (psps%usepaw==1) then 2177 call pawang_free(pawang1) 2178 call pawrhoij_free(pawrhoij1) 2179 if (usecprj==1) then 2180 call pawcprj_free(cprj) 2181 call pawcprj_free(cprjq) 2182 end if 2183 end if 2184 ABI_DATATYPE_DEALLOCATE(pawrhoij1) 2185 ABI_DATATYPE_DEALLOCATE(cprjq) 2186 ABI_DATATYPE_DEALLOCATE(cprj) 2187 if(xmpi_paral==1) then 2188 ABI_DEALLOCATE(mpi_enreg%proc_distrb) 2189 ABI_DEALLOCATE(mpi_enreg%my_kpttab) 2190 end if 2191 call hdr_free(hdr) 2192 2193 ! Clean band structure datatypes (should use it more in the future !) 2194 call ebands_free(ebands_k) 2195 call ebands_free(ebands_kq) 2196 if(.not.kramers_deg) then 2197 call ebands_free(ebands_kmq) 2198 end if 2199 2200 ! %%%% Parallelization over perturbations %%%%% 2201 ! *Redefine output/log files 2202 call localredirect(mpi_enreg%comm_cell,mpi_enreg%comm_world,npert_io,& 2203 & mpi_enreg%paral_pert,0) 2204 2205 call timab(146,2,tsec) 2206 if(iexit/=0) exit 2207 end do ! End loop on perturbations 2208 ABI_DEALLOCATE(zeff) 2209 2210 !%%%% Parallelization over perturbations %%%%% 2211 !*Restore default communicators 2212 call unset_pert_comm(mpi_enreg) 2213 !*Gather output/log files 2214 ABI_ALLOCATE(dyn,(npert_io)) 2215 if (npert_io>0) dyn=1 2216 call localrdfile(mpi_enreg%comm_pert,mpi_enreg%comm_world,.true.,npert_io,& 2217 & mpi_enreg%paral_pert,0,dyn) 2218 ABI_DEALLOCATE(dyn) 2219 !*Restore PAW on-site data 2220 if (paral_pert_inplace) then 2221 call unset_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,& 2222 & paw_an,paw_ij,pawfgrtab,pawrhoij) 2223 else 2224 call unset_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,& 2225 & paw_an,paw_ij,pawfgrtab,pawrhoij,& 2226 & paw_an_out=paw_an_pert,paw_ij_out=paw_ij_pert,& 2227 & pawfgrtab_out=pawfgrtab_pert,pawrhoij_out=pawrhoij_pert) 2228 end if 2229 2230 !################################################################################# 2231 !Calculate the second-order eigenvalues for a wavevector Q 2232 2233 call timab(147,1,tsec) 2234 smdelta = dtset%smdelta 2235 bdeigrf = dtset%bdeigrf 2236 if(dtset%bdeigrf == -1) bdeigrf = dtset%mband 2237 2238 if(dtset%ieig2rf > 0) then 2239 2240 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 2241 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 2242 call wrtout(std_out,'Reading the dense grid WF file',"COLL") 2243 ! We get the Abinit header of the file hdr_fine as ouput 2244 ! We get eigenq_fine(mband,hdr_fine%nkpt,hdr_fine%nsppol) as ouput 2245 fname = dtfil%fnameabi_wfkfine 2246 if (dtset%iomode == IO_MODE_ETSF) fname = nctk_ncify(dtfil%fnameabi_wfkfine) 2247 2248 call wfk_read_eigenvalues(fname,eigenq_fine,hdr_fine,mpi_enreg%comm_world) 2249 ABI_CHECK(SIZE(eigenq_fine,DIM=1)==Dtset%mband,"Size eigenq_fine != mband") 2250 end if 2251 ! DBSP ==> Has been changed to be able to make Bandstructure calculation 2252 ! if(dtset%kptopt==3 .or. dtset%kptopt==0)then 2253 if(dtset%kptopt==3 .or. dtset%kptopt==0 .or. dtset%kptopt < -4 .or. dtset%nsym==1) then 2254 !END 2255 if (dtset%nsym > 1) then ! .and. dtset%efmas==0) then 2256 MSG_ERROR("Symmetries are not implemented for temperature dependence calculations") 2257 end if 2258 write(std_out,*) 'Entering: eig2stern' 2259 if(smdelta>0)then 2260 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 2261 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 2262 call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,& 2263 & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,& 2264 & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,& 2265 & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,eigbrd,eigenq_fine,hdr_fine,hdr0) 2266 else 2267 call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,& 2268 & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,& 2269 & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,& 2270 & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,eigbrd) 2271 end if 2272 else 2273 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 2274 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 2275 call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,& 2276 & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,& 2277 & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,& 2278 & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,eigbrd,eigenq_fine,hdr_fine,hdr0) 2279 else 2280 call eig2stern(occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,& 2281 & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,& 2282 & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,& 2283 & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset) 2284 end if 2285 end if 2286 call wrtout(std_out, 'Leaving: eig2stern', "COLL") 2287 2288 if (dtset%ieig2rf==1.or.dtset%ieig2rf==2) then 2289 if (me==master) then 2290 ! print _EIGR2D file for this perturbation 2291 dscrpt=' Note : temporary (transfer) database ' 2292 unitout = dtfil%unddb 2293 vrsddb=100401 2294 2295 call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,& 2296 & 1,xred=xred,occ=occ_pert) 2297 2298 call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_eigr2d, dtfil%unddb) 2299 2300 call outbsd(bdeigrf,dtset,eig2nkq,dtset%natom,nkpt_rbz,unitout) 2301 ! print _EIGI2D file for this perturbation 2302 if(smdelta>0) then 2303 2304 unitout = dtfil%unddb 2305 call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_eigi2d, unitout) 2306 2307 call outbsd(bdeigrf,dtset,eigbrd,dtset%natom,nkpt_rbz,unitout) 2308 end if !smdelta 2309 2310 call ddb_hdr_free(ddb_hdr) 2311 2312 ! Output of the EIGR2D.nc file. 2313 fname = strcat(dtfil%filnam_ds(4),"_EIGR2D.nc") 2314 2315 ! Electronic band energies. 2316 bantot= dtset%mband*dtset%nkpt*dtset%nsppol 2317 call ebands_init(bantot,gkk_ebands,dtset%nelect,doccde,eigen0_pert,hdr0%istwfk,hdr0%kptns,& 2318 & hdr0%nband, hdr0%nkpt,hdr0%npwarr,hdr0%nsppol,hdr0%nspinor,& 2319 & hdr0%tphysel,hdr0%tsmear,hdr0%occopt,hdr0%occ,hdr0%wtk,& 2320 & hdr0%charge, hdr0%kptopt, hdr0%kptrlatt_orig, hdr0%nshiftk_orig, hdr0%shiftk_orig, & 2321 & hdr0%kptrlatt, hdr0%nshiftk, hdr0%shiftk) 2322 2323 ! Second order derivative EIGR2D (real and Im) 2324 call eigr2d_init(eig2nkq,eigr2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom) 2325 #ifdef HAVE_NETCDF 2326 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EIGR2D file") 2327 NCF_CHECK(crystal_ncwrite(crystal, ncid)) 2328 NCF_CHECK(ebands_ncwrite(gkk_ebands, ncid)) 2329 call eigr2d_ncwrite(eigr2d,dtset%qptn(:),dtset%wtq,ncid) 2330 NCF_CHECK(nf90_close(ncid)) 2331 #endif 2332 if(smdelta>0) then 2333 ! Output of the EIGI2D.nc file. 2334 fname = strcat(dtfil%filnam_ds(4),"_EIGI2D.nc") 2335 ! Broadening EIGI2D (real and Im) 2336 call eigr2d_init(eigbrd,eigi2d,dtset%mband,hdr0%nsppol,nkpt_rbz,dtset%natom) 2337 #ifdef HAVE_NETCDF 2338 NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EIGI2D file") 2339 NCF_CHECK(crystal_ncwrite(crystal, ncid)) 2340 NCF_CHECK(ebands_ncwrite(gkk_ebands, ncid)) 2341 call eigr2d_ncwrite(eigi2d,dtset%qptn(:),dtset%wtq,ncid) 2342 NCF_CHECK(nf90_close(ncid)) 2343 #endif 2344 end if 2345 end if 2346 end if !ieig2rf==1.or.ieig2rf==2 2347 else 2348 write(message,'(3a)')& 2349 & 'K point grids must be the same for every perturbation: eig2stern not called',ch10,& 2350 & 'Action: Put kptopt=3 ' 2351 MSG_WARNING(message) 2352 end if !kptopt 2353 ABI_DEALLOCATE(gh1c_pert) 2354 ABI_DEALLOCATE(gh0c1_pert) 2355 ABI_DEALLOCATE(cg1_pert) 2356 ABI_DEALLOCATE(kpt_rbz_pert) 2357 ABI_DEALLOCATE(istwfk_pert) 2358 ABI_DEALLOCATE(npwarr_pert) 2359 ABI_DEALLOCATE(npwar1_pert) 2360 ABI_DEALLOCATE(npwtot_pert) 2361 ABI_DEALLOCATE(occ_pert) 2362 end if !if dtset%ieig2rf 2363 2364 !Calculation of effective masses. 2365 if(dtset%efmas == 1) then 2366 call efmas_main(cg0_pert,cg1_pert,dim_eig2rf,dtset,efmasdeg,efmasfr,eigen0_pert,& 2367 & eigen1_pert,gh0c1_pert,gh1c_pert,istwfk_pert,& 2368 & kpt_rbz_pert,mpert,mpi_enreg,& 2369 & nkpt_rbz,npwarr_pert,rprimd) 2370 ABI_DEALLOCATE(gh1c_pert) 2371 ABI_DEALLOCATE(gh0c1_pert) 2372 ABI_DEALLOCATE(cg1_pert) 2373 ABI_DEALLOCATE(kpt_rbz_pert) 2374 ABI_DEALLOCATE(istwfk_pert) 2375 ABI_DEALLOCATE(npwarr_pert) 2376 ABI_DEALLOCATE(cg0_pert) 2377 end if 2378 2379 2380 !Free memory. 2381 if(dtset%ieig2rf /= 3 .and. dtset%ieig2rf /= 4 .and. dtset%ieig2rf /= 5) then 2382 call hdr_free(hdr0) 2383 end if 2384 ABI_DEALLOCATE(eigen0_copy) 2385 call crystal_free(crystal) 2386 2387 ! GKK stuff (deprecated) 2388 call ebands_free(gkk_ebands) 2389 call eigr2d_free(eigr2d) 2390 call eigr2d_free(eigi2d) 2391 2392 call timab(147,2,tsec) 2393 !###################################################################################### 2394 2395 !Get ddk file information, for later use in dfpt_dyout 2396 ddkfil(:)=0 2397 do idir=1,3 2398 file_index(1)=idir+dtset%natom*3 2399 call appdig(file_index(1),dtfil%fnamewffddk,fiwfddk) 2400 ! Check that ddk file exists 2401 t_exist = file_exists(fiwfddk) 2402 if (.not. t_exist) then 2403 ! Trick needed to run Abinit test suite in netcdf mode. 2404 t_exist = file_exists(nctk_ncify(fiwfddk)) 2405 if (t_exist) then 2406 write(message,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name." 2407 call wrtout(std_out,message,'COLL') 2408 fiwfddk = nctk_ncify(fiwfddk) 2409 end if 2410 end if 2411 2412 ! If the file exists set ddkfil to a non-zero value 2413 if (t_exist) ddkfil(idir)=20+idir 2414 end do 2415 2416 ABI_DEALLOCATE(ph1d) 2417 ABI_DEALLOCATE(ph1df) 2418 ABI_DEALLOCATE(pert_calc) 2419 if (psps%usepaw==1) then 2420 ABI_DEALLOCATE(dimcprj_srt) 2421 end if 2422 2423 !destroy dtset_tmp 2424 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 2425 ABI_DEALLOCATE(rhor1_save) 2426 ABI_DEALLOCATE(blkflg_save) 2427 call dtset_free (dtset_tmp) 2428 ABI_DATATYPE_DEALLOCATE(dtset_tmp) 2429 end if 2430 2431 !In paral_pert-case some array's have to be reconstructed 2432 if(mpi_enreg%paral_pert==1) then 2433 ABI_ALLOCATE(buffer1,(2,3,mpert,3,mpert*(2+psps%usepaw))) 2434 buffer1(:,:,:,:,1:mpert)=d2lo(:,:,:,:,:) 2435 buffer1(:,:,:,:,1+mpert:2*mpert)=d2nl(:,:,:,:,:) 2436 if (psps%usepaw==1) then 2437 buffer1(:,:,:,:,1+2*mpert:3*mpert)=d2ovl(:,:,:,:,:) 2438 end if 2439 call xmpi_sum(buffer1,mpi_enreg%comm_pert,ierr) 2440 call xmpi_sum(blkflg,mpi_enreg%comm_pert,ierr) 2441 if(dtset%prtbbb==1) then 2442 call xmpi_sum(d2bbb,mpi_enreg%comm_pert,ierr) 2443 end if 2444 d2lo(:,:,:,:,:)=buffer1(:,:,:,:,1:mpert) 2445 d2nl(:,:,:,:,:)=buffer1(:,:,:,:,1+mpert:2*mpert) 2446 if (psps%usepaw==1) then 2447 d2ovl(:,:,:,:,:)=buffer1(:,:,:,:,1+2*mpert:3*mpert) 2448 end if 2449 ABI_DEALLOCATE(buffer1) 2450 end if 2451 2452 if ( associated(old_atmtab)) then 2453 ABI_DEALLOCATE(old_atmtab) 2454 nullify(old_atmtab) 2455 end if 2456 2457 call status(0,dtfil%filstat,iexit,level,'exit') 2458 2459 call timab(141,2,tsec) 2460 2461 DBG_EXIT("COLL") 2462 2463 end subroutine dfpt_looppert