TABLE OF CONTENTS
ABINIT/respfn [ Functions ]
NAME
respfn
FUNCTION
Primary routine for conducting DFT calculations of Response functions.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (XG, DRH, MT) 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
codvsn=code version cpui=initial cpu time dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset | mband=maximum number of bands | mgfft=maximum single fft dimension | mkmem=Number of k points treated by this node | mpw=maximum number of planewaves in basis sphere (large number) | natom=number of atoms in unit cell | nfft=(effective) number of FFT grid points (for this processor) | nkpt=number of k points | nspden=number of spin-density components | nsppol=number of channels for spin-polarization (1 or 2) | nsym=number of symmetry elements in space group mkmems(3)=array containing the tree values of mkmem (see above) (k-GS, k+q-GS and RF) mpi_enreg=information about MPI parallelization npwtot(nkpt)=number of planewaves in basis and boundary at each k point xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
etotal=total energy (sum of 7 or 8 contributions) (hartree)
SIDE EFFECTS
iexit=index of "exit" on first line of file (0 if not found)
occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point
Occupations number may have been read from a previous dataset...
pawang <type(pawang_type)>=paw angular mesh and related data
pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
Some dimensions in pawrad have been set in driver.f
pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
Some dimensions in pawtab have been set in driver.f
psps <type(pseudopotential_type)>=variables related to pseudopotentials
Before entering the first time in respfn, a significant part of psps
has been initialized: the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,
mpsso,mgrid,ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp
All the remaining components of psps are to be initialized in the call
to pspini. The next time the code enters respfn, psps might be identical
to the one of the previous dtset, in which case, no reinitialisation
is scheduled in pspini.f .
results_respfn <type(results_respfn_type)>=stores some results of respfn calls
NOTES
USE OF FFT GRIDS:
=================
In case of PAW:
---------------
Two FFT grids are used:
- A "coarse" FFT grid (defined by ecut)
for the application of the Hamiltonian on the plane waves basis.
It is defined by nfft, ngfft, mgfft, ...
Hamiltonian, wave-functions, density related to WFs (rhor here), ...
are expressed on this grid.
- A "fine" FFT grid (defined) by ecutdg)
for the computation of the density inside PAW spheres.
It is defined by nfftf, ngfftf, mgfftf, ...
Total density, potentials, ...
are expressed on this grid.
In case of norm-conserving:
---------------------------
- Only the usual FFT grid (defined by ecut) is used.
It is defined by nfft, ngfft, mgfft, ...
For compatibility reasons, (nfftf,ngfftf,mgfftf)
are set equal to (nfft,ngfft,mgfft) in that case.
PARENTS
driver
CHILDREN
alloc_hamilt_gpu,atm2fft,check_kxc,chkpawovlp,chkph3,crystal_free crystal_init,d2frnl,d2sym3,ddb_hdr_free,ddb_hdr_init,ddb_hdr_open_write dealloc_hamilt_gpu,dfpt_dyfro,dfpt_dyout,dfpt_dyxc1,dfpt_eltfrhar dfpt_eltfrkin,dfpt_eltfrloc,dfpt_eltfrxc,dfpt_ewald,dfpt_gatherdy dfpt_looppert,dfpt_phfrq,dfpt_prtph,ebands_free,efmasdeg_free_array efmasfr_free_array,eig2tot,eigen_meandege,elph2_fanddw,elt_ewald exit_check,fourdp,getcut,getph,hartre,hdr_free,hdr_init,hdr_update initrhoij,initylmg,inwffil,irreducible_set_pert,kpgio,littlegroup_q matr3inv,mkcore,mklocl,mkrho,newocc,nhatgrid,outddbnc,paw_an_free paw_an_init,paw_an_nullify,paw_gencond,paw_ij_free,paw_ij_init paw_ij_nullify,pawdenpot,pawdij,pawexpiqr,pawfgr_destroy,pawfgr_init pawfgrtab_free,pawfgrtab_init,pawinit,pawmknhat,pawpuxinit pawrhoij_alloc,pawrhoij_bcast,pawrhoij_copy,pawrhoij_free pawrhoij_nullify,pawtab_get_lsize,prteigrs,pspini,q0dy3_apply q0dy3_calc,read_rhor,rhotoxc,setsym,setsymrhoij,setup1,status,symdij symmetrize_xred,sytens,timab,transgrid,vdw_dftd2,vdw_dftd3,wffclose wings3,wrtloctens,wrtout,xcdata_init,xmpi_bcast
SOURCE
105 #if defined HAVE_CONFIG_H 106 #include "config.h" 107 #endif 108 109 #include "abi_common.h" 110 111 subroutine respfn(codvsn,cpui,dtfil,dtset,etotal,iexit,& 112 & mkmems,mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,psps,results_respfn,xred) 113 114 use defs_basis 115 use defs_datatypes 116 use defs_abitypes 117 use defs_wvltypes 118 use m_efmas_defs 119 use m_profiling_abi 120 use m_xmpi 121 use m_exit 122 use m_wffile 123 use m_errors 124 use m_ebands 125 use m_results_respfn 126 use m_hdr 127 use m_crystal 128 use m_xcdata 129 130 use m_fstrings, only : strcat 131 use m_kpts, only : symkchk 132 use m_dynmat, only : chkph3, d2sym3, q0dy3_apply, q0dy3_calc, wings3, dfpt_phfrq, sytens 133 use m_ddb, only : DDB_VERSION 134 use m_ddb_hdr, only : ddb_hdr_type, ddb_hdr_init, ddb_hdr_free, ddb_hdr_open_write 135 use m_occ, only : newocc 136 use m_efmas, only : efmasdeg_free_array, efmasfr_free_array 137 use m_wfk, only : wfk_read_eigenvalues 138 use m_ioarr, only : read_rhor 139 use m_pawang, only : pawang_type 140 use m_pawrad, only : pawrad_type 141 use m_pawtab, only : pawtab_type,pawtab_get_lsize 142 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 143 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 144 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free 145 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy, & 146 & pawrhoij_bcast, pawrhoij_nullify, pawrhoij_get_nspden 147 use m_pawdij, only : pawdij, symdij 148 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 149 use m_paw_finegrid,only : pawexpiqr 150 use m_paw_dmft, only : paw_dmft_type 151 use m_kg, only : getcut, getph, kpgio 152 153 !This section has been created automatically by the script Abilint (TD). 154 !Do not modify the following lines by hand. 155 #undef ABI_FUNC 156 #define ABI_FUNC 'respfn' 157 use interfaces_14_hidewrite 158 use interfaces_18_timing 159 use interfaces_32_util 160 use interfaces_41_geometry 161 use interfaces_41_xc_lowlevel 162 #if defined HAVE_GPU_CUDA 163 use interfaces_52_manage_cuda 164 #endif 165 use interfaces_53_ffts 166 use interfaces_56_recipspace 167 use interfaces_56_xc 168 use interfaces_64_psp 169 use interfaces_65_paw 170 use interfaces_67_common 171 use interfaces_72_response 172 use interfaces_79_seqpar_mpi 173 use interfaces_95_drive, except_this_one => respfn 174 !End of the abilint section 175 176 implicit none 177 178 !Arguments ------------------------------------ 179 integer,intent(inout) :: iexit 180 real(dp),intent(in) :: cpui 181 real(dp),intent(inout) :: etotal !vz_i 182 character(len=6),intent(in) :: codvsn 183 type(MPI_type),intent(inout) :: mpi_enreg 184 type(datafiles_type),intent(in) :: dtfil 185 type(dataset_type),intent(in) :: dtset 186 type(pawang_type),intent(inout) :: pawang 187 type(pseudopotential_type),intent(inout) :: psps 188 integer,intent(in) :: mkmems(3) 189 integer,intent(inout) :: npwtot(dtset%nkpt) 190 real(dp),intent(inout) :: xred(3,dtset%natom) 191 real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 192 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 193 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 194 type(results_respfn_type),intent(inout) :: results_respfn 195 196 !Local variables------------------------------- 197 integer,parameter :: formeig=0,level=10 198 integer,parameter :: response=1,syuse=0,master=0,cplex1=1 199 integer :: nk3xc 200 integer :: analyt,ask_accurate,band_index,bantot,bdeigrf,coredens_method,cplex 201 integer :: dim_eig2nkq,dim_eigbrd,dyfr_cplex,dyfr_nondiag,gnt_option 202 integer :: gscase,has_dijnd,has_kxc,iatom,iatom_tot,iband,idir,ider,ierr,ifft,ii,ikpt,indx 203 integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert 204 integer :: initialized,ipert,ipert2,ireadwf0,iscf,iscf_eff,ispden,isppol 205 integer :: itypat,izero,mcg,me,mgfftf,mk1mem,mkqmem,mpert,mu 206 integer :: my_natom,n1,natom,n3xccc,nband_k,nfftf,nfftot,nfftotf,nhatdim,nhatgrdim 207 integer :: nkpt_eff,nkpt_max,nkpt_rbz,nkxc,nkxc1,nspden_rhoij,ntypat,nzlmopt,openexit 208 integer :: optcut,option,optgr0,optgr1,optgr2,optorth,optrad 209 integer :: optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv 210 integer :: outd2,pawbec,pawpiezo,prtbbb,psp_gencond,qzero,rdwr,rdwrpaw 211 integer :: req_cplex_dij,rfasr,rfddk,rfelfd,rfphon,rfstrs,rfuser,rf2_dkdk,rf2_dkde,rfmagn 212 integer :: spaceworld,sumg0,sz1,sz2,tim_mkrho,timrev,usecprj,usevdw 213 integer :: usexcnhat,use_sym,vloc_method 214 logical :: has_full_piezo,has_allddk,paral_atom,qeq0,use_nhat_gga,call_pawinit,non_magnetic_xc 215 real(dp) :: boxcut,compch_fft,compch_sph,cpus,ecore,ecut_eff,ecutdg_eff,ecutf 216 real(dp) :: eei,eew,ehart,eii,ek,enl,entropy,enxc 217 real(dp) :: epaw,epawdc,etot,evdw,fermie,gsqcut,gsqcut_eff,gsqcutc_eff,qphnrm,residm 218 real(dp) :: ucvol,vxcavg 219 character(len=fnlen) :: dscrpt 220 character(len=fnlen) :: filename 221 character(len=500) :: message 222 type(ebands_t) :: bstruct 223 type(hdr_type) :: hdr,hdr_fine,hdr0,hdr_den 224 type(ddb_hdr_type) :: ddb_hdr 225 type(paw_dmft_type) :: paw_dmft 226 type(pawfgr_type) :: pawfgr 227 type(wffile_type) :: wffgs,wfftgs 228 type(wvl_data) :: wvl 229 type(crystal_t) :: Crystal 230 type(xcdata_type) :: xcdata 231 integer :: ddkfil(3),ngfft(18),ngfftf(18),rfdir(3),rf2_dirs_from_rfpert_nl(3,3) 232 integer,allocatable :: atindx(:),atindx1(:),blkflg(:,:,:,:),blkflgfrx1(:,:,:,:),blkflg1(:,:,:,:) 233 integer,allocatable :: blkflg2(:,:,:,:),carflg(:,:,:,:),clflg(:,:),indsym(:,:,:) 234 integer,allocatable :: irrzon(:,:,:),kg(:,:),l_size_atm(:),nattyp(:),npwarr(:) 235 integer,allocatable :: pertsy(:,:),rfpert(:),rfpert_nl(:,:,:,:,:,:),symq(:,:,:),symrec(:,:,:) 236 real(dp) :: dum_gauss(0),dum_dyfrn(0),dum_dyfrv(0),dum_eltfrxc(0) 237 real(dp) :: dum_grn(0),dum_grv(0),dum_rhog(0),dum_vg(0) 238 real(dp) :: dummy6(6),gmet(3,3),gmet_for_kg(3,3),gprimd(3,3),gprimd_for_kg(3,3),qphon(3) 239 real(dp) :: rmet(3,3),rprimd(3,3),rprimd_for_kg(3,3),strn_dummy6(6),strv_dummy6(6),strsxc(6),tsec(2) 240 real(dp),parameter :: k0(3)=(/zero,zero,zero/) 241 real(dp),allocatable :: amass(:),becfrnl(:,:,:),cg(:,:),d2bbb(:,:,:,:,:,:),d2cart(:,:,:,:,:) 242 real(dp),allocatable :: d2cart_bbb(:,:,:,:,:,:),d2eig0(:,:,:,:,:) 243 real(dp),allocatable :: d2k0(:,:,:,:,:),d2lo(:,:,:,:,:),d2loc0(:,:,:,:,:) 244 real(dp),allocatable :: d2matr(:,:,:,:,:),d2nfr(:,:,:,:,:),d2nl(:,:,:,:,:),d2ovl(:,:,:,:,:) 245 real(dp),allocatable :: d2nl0(:,:,:,:,:),d2nl1(:,:,:,:,:),d2tmp(:,:,:,:,:),d2vn(:,:,:,:,:) 246 real(dp),allocatable :: displ(:),doccde(:) 247 real(dp),allocatable :: dyew(:,:,:,:,:),dyewq0(:,:,:),dyfrlo(:,:,:),dyfrlo_indx(:,:,:) 248 real(dp),allocatable :: dyfrnl(:,:,:,:,:),dyfrwf(:,:,:,:,:),dyfrx1(:,:,:,:,:),dyvdw(:,:,:,:,:) 249 real(dp),allocatable :: dyfrx2(:,:,:),eigen0(:),eigval(:),eigvec(:) 250 real(dp),allocatable :: eig2nkq(:,:,:,:,:,:,:),eigbrd(:,:,:,:,:,:,:) 251 real(dp),allocatable :: eigen_fan(:),eigen_ddw(:),eigen_fanddw(:) 252 real(dp),allocatable :: eigen_fan_mean(:),eigen_ddw_mean(:) 253 real(dp),allocatable :: eltcore(:,:),elteew(:,:),eltfrhar(:,:),eltfrkin(:,:) 254 real(dp),allocatable :: eltfrloc(:,:),eltfrnl(:,:),eltfrxc(:,:),eltvdw(:,:),grtn_indx(:,:) 255 real(dp),allocatable :: grxc(:,:),kxc(:,:),nhat(:,:),nhatgr(:,:,:) 256 real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phfrq(:),phnons(:,:,:),piezofrnl(:,:) 257 real(dp),allocatable :: rhog(:,:),rhor(:,:),rhowfg(:,:),rhowfr(:,:) 258 real(dp),allocatable :: vhartr(:),vpsp(:),vtrial(:,:) 259 real(dp),allocatable :: vxc(:,:),work(:),xccc3d(:),ylm(:,:),ylmgr(:,:,:) 260 real(dp),pointer :: eigenq_fine(:,:,:),eigen1_pert(:,:,:) 261 real(dp),allocatable :: eigen0_pert(:),eigenq_pert(:),occ_rbz_pert(:) 262 type(efmasdeg_type),allocatable :: efmasdeg(:) 263 type(efmasfr_type),allocatable :: efmasfr(:,:) 264 type(paw_an_type),allocatable :: paw_an(:) 265 type(paw_ij_type),allocatable :: paw_ij(:) 266 type(pawfgrtab_type),allocatable,save :: pawfgrtab(:) 267 type(pawrhoij_type),allocatable :: pawrhoij(:),pawrhoij_read(:) 268 ! *********************************************************************** 269 270 DBG_ENTER("COLL") 271 272 call timab(132,1,tsec) 273 call timab(133,1,tsec) 274 275 call status(0,dtfil%filstat,iexit,level,'init ') 276 277 ! Initialise non_magnetic_xc for rhohxc 278 non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14) 279 280 !Some data for parallelism 281 nkpt_max=50;if(xmpi_paral==1)nkpt_max=-1 282 my_natom=mpi_enreg%my_natom 283 paral_atom=(my_natom/=dtset%natom) 284 !Define FFT grid(s) sizes (be careful !) 285 !See NOTES in the comments at the beginning of this file. 286 call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf) 287 288 !Structured debugging if dtset%prtvol==-level 289 if(dtset%prtvol==-level)then 290 write(message,'(80a,a,a)') ('=',ii=1,80),ch10,' respfn : enter , debug mode ' 291 call wrtout(std_out,message,'COLL') 292 end if 293 294 !Option input variables 295 iscf=dtset%iscf 296 297 !Respfn input variables 298 rfasr=dtset%rfasr ; rfdir(1:3)=dtset%rfdir(1:3) 299 rfddk=dtset%rfddk ; rfelfd=dtset%rfelfd ; rfmagn=dtset%rfmagn 300 rfphon=dtset%rfphon ; rfstrs=dtset%rfstrs 301 rfuser=dtset%rfuser ; rf2_dkdk=dtset%rf2_dkdk ; rf2_dkde=dtset%rf2_dkde 302 303 pawbec=0 ; if(psps%usepaw==1.and.(rfphon==1.or.(rfelfd==1.or.rfelfd==3))) pawbec=1 304 pawpiezo=0; if(psps%usepaw==1.and.(rfstrs/=0.or.(rfelfd==1.or.rfelfd==3))) pawpiezo=1 305 !AM 10152015 -- WARNING --- the full calculation of the piezoelectric tensor 306 !from electric field perturbation is only available 307 !if nsym/=1 (strain perturbation is not symmetrized): 308 has_full_piezo=.False. ; if(pawpiezo==1.and.dtset%nsym==1) has_full_piezo=.True. 309 usevdw=0;if (dtset%vdw_xc>=5.and.dtset%vdw_xc<=7) usevdw=1 310 !mkmem variables (mkmem is already argument) 311 mkqmem=mkmems(2) ; mk1mem=mkmems(3) 312 313 ntypat=psps%ntypat 314 natom=dtset%natom 315 nfftot=product(ngfft(1:3)) 316 nfftotf=product(ngfftf(1:3)) 317 318 call status(0,dtfil%filstat,iexit,level,'call setup1 ') 319 320 !LIKELY TO BE TAKEN AWAY 321 initialized=0 322 ek=zero ; ehart=zero ; enxc=zero ; eei=zero ; enl=zero 323 eii=zero ; eew=zero ; ecore=zero 324 325 !Set up for iterations 326 ABI_ALLOCATE(amass,(natom)) 327 call setup1(dtset%acell_orig(1:3,1),amass,dtset%amu_orig(:,1),bantot,dtset,& 328 & ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,& 329 & natom,ngfftf,ngfft,dtset%nkpt,dtset%nsppol,& 330 & response,rmet,dtset%rprim_orig(1:3,1:3,1),rprimd,ucvol,psps%usepaw) 331 332 !In some cases (e.g. getcell/=0), the plane wave vectors have 333 ! to be generated from the original simulation cell 334 rprimd_for_kg=rprimd 335 if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=dtset%rprimd_orig(:,:,1) 336 call matr3inv(rprimd_for_kg,gprimd_for_kg) 337 gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg) 338 339 !Define the set of admitted perturbations 340 mpert=natom+7 341 if (rf2_dkdk>0.or.rf2_dkde>0) mpert=natom+11 342 343 !Initialize the list of perturbations rfpert 344 ABI_ALLOCATE(rfpert,(mpert)) 345 rfpert(:)=0 346 if(rfphon==1)rfpert(dtset%rfatpol(1):dtset%rfatpol(2))=1 347 348 if(rfddk==1)rfpert(natom+1)=1 349 if(rfddk==2)rfpert(natom+6)=1 350 351 if(rf2_dkdk/=0)rfpert(natom+10)=1 352 if(rf2_dkde/=0)rfpert(natom+11)=1 353 354 if(rfelfd==1.or.rfelfd==2)rfpert(natom+1)=1 355 if(rfelfd==1.or.rfelfd==3)rfpert(natom+2)=1 356 357 if(rfstrs==1.or.rfstrs==3)rfpert(natom+3)=1 358 if(rfstrs==2.or.rfstrs==3)rfpert(natom+4)=1 359 360 if(rfuser==1.or.rfuser==3)rfpert(natom+6)=1 361 if(rfuser==2.or.rfuser==3)rfpert(natom+7)=1 362 363 if(rfmagn==1) rfpert(natom+5)=1 364 365 qeq0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2<1.d-14) 366 367 !Init spaceworld 368 spaceworld=mpi_enreg%comm_cell 369 me = xmpi_comm_rank(spaceworld) 370 371 !Set up the basis sphere of planewaves 372 ABI_ALLOCATE(kg,(3,dtset%mpw*dtset%mkmem)) 373 ABI_ALLOCATE(npwarr,(dtset%nkpt)) 374 call status(0,dtfil%filstat,iexit,level,'call kpgio(1) ') 375 call kpgio(ecut_eff,dtset%exchn2n3d,gmet_for_kg,dtset%istwfk,kg,& 376 & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw,npwarr,npwtot,& 377 & dtset%nsppol) 378 379 !Set up the Ylm for each k point 380 ABI_ALLOCATE(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)) 381 if (rfstrs/=0.or.pawbec==1.or.pawpiezo==1.or.dtset%efmas>0) then 382 ABI_ALLOCATE(ylmgr,(dtset%mpw*dtset%mkmem,9,psps%mpsang*psps%mpsang*psps%useylm)) 383 else 384 ABI_ALLOCATE(ylmgr,(0,0,psps%useylm)) 385 end if 386 if (psps%useylm==1) then 387 call status(0,dtfil%filstat,iexit,level,'call initylmg ') 388 option=0 389 if (rfstrs/=0.or.pawbec==1.or.pawpiezo==1.or.dtset%efmas>0) option=2 390 call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,& 391 & npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr) 392 end if 393 394 call timab(133,2,tsec) 395 call timab(134,1,tsec) 396 397 !Open and read pseudopotential files 398 call status(0,dtfil%filstat,iexit,level,'call pspini(1)') 399 call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,& 400 & psps,rprimd,comm_mpi=mpi_enreg%comm_cell) 401 402 call timab(134,2,tsec) 403 call timab(135,1,tsec) 404 405 !Initialize band structure datatype 406 bstruct = ebands_from_dtset(dtset, npwarr) 407 408 !Initialize PAW atomic occupancies 409 if (psps%usepaw==1) then 410 ABI_DATATYPE_ALLOCATE(pawrhoij,(my_natom)) 411 call pawrhoij_nullify(pawrhoij) 412 call initrhoij(dtset%pawcpxocc,dtset%lexexch,dtset%lpawu, & 413 & my_natom,natom,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat,& 414 & pawrhoij,dtset%pawspnorb,pawtab,dtset%spinat,dtset%typat,& 415 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 416 else 417 ABI_DATATYPE_ALLOCATE(pawrhoij,(0)) 418 end if 419 420 !Initialize header 421 gscase=0 422 call hdr_init(bstruct,codvsn,dtset,hdr,pawtab,gscase,psps,wvl%descr, & 423 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 424 425 !Update header, with evolving variables, when available 426 !Here, rprimd, xred and occ are available 427 etot=hdr%etot ; fermie=hdr%fermie ; residm=hdr%residm 428 !If parallelism over atom, hdr is distributed 429 call hdr_update(hdr,bantot,etot,fermie,& 430 & residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1), & 431 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 432 433 !Clean band structure datatype (should use it more in the future !) 434 call ebands_free(bstruct) 435 436 call status(0,dtfil%filstat,iexit,level,'call inwffil(1') 437 438 !Initialize wavefunction files and wavefunctions. 439 ireadwf0=1 440 441 mcg=dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol 442 ABI_STAT_ALLOCATE(cg,(2,mcg), ierr) 443 ABI_CHECK(ierr==0, "out-of-memory in cg") 444 445 ABI_ALLOCATE(eigen0,(dtset%mband*dtset%nkpt*dtset%nsppol)) 446 eigen0(:)=zero ; ask_accurate=1 447 optorth=0 448 449 hdr%rprimd=rprimd_for_kg ! We need the rprimd that was used to generate de G vectors 450 call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,& 451 & formeig,hdr,ireadwf0,dtset%istwfk,kg,dtset%kptns,& 452 & dtset%localrdwf,dtset%mband,mcg,dtset%mkmem,mpi_enreg,dtset%mpw,& 453 & dtset%nband,ngfft,dtset%nkpt,npwarr,dtset%nsppol,dtset%nsym,& 454 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 455 & dtfil%unkg,wffgs,wfftgs,dtfil%unwffgs,dtfil%fnamewffk,wvl) 456 hdr%rprimd=rprimd 457 458 !Close wffgs, if it was ever opened (in inwffil) 459 if (ireadwf0==1) then 460 call WffClose(wffgs,ierr) 461 end if 462 463 if (psps%usepaw==1.and.ireadwf0==1) then 464 ! if parallelism, pawrhoij is distributed, hdr%pawrhoij is not 465 call pawrhoij_copy(hdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,& 466 & mpi_atmtab=mpi_enreg%my_atmtab) 467 end if 468 469 call timab(135,2,tsec) 470 call timab(136,1,tsec) 471 472 !Report on eigen0 values ! Should use prteigrs.F90 473 write(message, '(a,a)' ) 474 call wrtout(std_out,ch10//' respfn : eigen0 array','COLL') 475 nkpt_eff=dtset%nkpt 476 if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. dtset%nkpt>nkpt_max ) nkpt_eff=nkpt_max 477 band_index=0 478 do isppol=1,dtset%nsppol 479 do ikpt=1,dtset%nkpt 480 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 481 if(ikpt<=nkpt_eff)then 482 write(message, '(a,i2,a,i5)' )' isppol=',isppol,', k point number',ikpt 483 call wrtout(std_out,message,'COLL') 484 do iband=1,nband_k,4 485 write(message, '(a,4es16.6)')' ',eigen0(iband+band_index:min(iband+3,nband_k)+band_index) 486 call wrtout(std_out,message,'COLL') 487 end do 488 else if(ikpt==nkpt_eff+1)then 489 write(message,'(a,a)' )' respfn : prtvol=0, 1 or 2, stop printing eigen0.',ch10 490 call wrtout(std_out,message,'COLL') 491 end if 492 band_index=band_index+nband_k 493 end do 494 end do 495 496 !Allocation for forces and atomic positions (should be taken away, also argument ... ) 497 ABI_ALLOCATE(grxc,(3,natom)) 498 499 !Do symmetry stuff 500 ABI_ALLOCATE(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 501 ABI_ALLOCATE(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 502 ABI_ALLOCATE(indsym,(4,dtset%nsym,natom)) 503 ABI_ALLOCATE(symrec,(3,3,dtset%nsym)) 504 irrzon=0;indsym=0;symrec=0;phnons=zero 505 !If the density is to be computed by mkrho, need irrzon and phnons 506 iscf_eff=0;if(dtset%getden==0)iscf_eff=1 507 call setsym(indsym,irrzon,iscf_eff,natom,& 508 & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,& 509 & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred) 510 511 !Symmetrize atomic coordinates over space group elements: 512 call symmetrize_xred(indsym,natom,dtset%nsym,dtset%symrel,dtset%tnons,xred) 513 514 !Examine the symmetries of the q wavevector 515 ABI_ALLOCATE(symq,(4,2,dtset%nsym)) 516 timrev=1 517 518 ! By default use symmetries. 519 use_sym = 1 520 if (dtset%prtgkk == 1)then 521 use_sym = 0 522 call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol,use_sym=use_sym) 523 else 524 call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol) 525 end if 526 527 !Generate an index table of atoms, in order for them to be used 528 !type after type. 529 ABI_ALLOCATE(atindx,(natom)) 530 ABI_ALLOCATE(atindx1,(natom)) 531 ABI_ALLOCATE(nattyp,(ntypat)) 532 indx=1 533 do itypat=1,ntypat 534 nattyp(itypat)=0 535 do iatom=1,natom 536 if(dtset%typat(iatom)==itypat)then 537 atindx(iatom)=indx 538 atindx1(indx)=iatom 539 indx=indx+1 540 nattyp(itypat)=nattyp(itypat)+1 541 end if 542 end do 543 end do 544 545 !Here allocation of GPU for fft calculations 546 #if defined HAVE_GPU_CUDA 547 if (dtset%use_gpu_cuda==1) then 548 call alloc_hamilt_gpu(atindx1,dtset,gprimd,mpi_enreg,nattyp,npwarr,0,psps,dtset%use_gpu_cuda) 549 end if 550 #endif 551 552 !Compute structure factor phases for current atomic pos: 553 ABI_ALLOCATE(ph1d,(2,3*(2*dtset%mgfft+1)*natom)) 554 ABI_ALLOCATE(ph1df,(2,3*(2*mgfftf+1)*natom)) 555 call getph(atindx,natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred) 556 557 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then 558 call getph(atindx,natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred) 559 else 560 ph1df(:,:)=ph1d(:,:) 561 end if 562 563 !Compute occupation numbers and fermi energy, in case occupation scheme is metallic. 564 ABI_ALLOCATE(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol)) 565 if( dtset%occopt>=3.and.dtset%occopt<=8 ) then 566 567 call newocc(doccde,eigen0,entropy,fermie,dtset%spinmagntarget,dtset%mband,dtset%nband,& 568 & dtset%nelect,dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,dtset%occopt,dtset%prtvol,dtset%stmbias,& 569 & dtset%tphysel,dtset%tsmear,dtset%wtk) 570 571 ! Update fermie and occ 572 etot=hdr%etot ; residm=hdr%residm 573 call hdr_update(hdr,bantot,etot,fermie,residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1)) 574 575 else 576 ! doccde is irrelevant in this case 577 doccde(:)=zero 578 end if 579 580 !Recompute first large sphere cut-off gsqcut, without taking into account dilatmx 581 ecutf=dtset%ecut 582 if (psps%usepaw==1) then 583 ecutf=dtset%pawecutdg 584 call wrtout(std_out,ch10//' FFT (fine) grid used in SCF cycle:','COLL') 585 end if 586 587 call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,k0,ngfftf) 588 589 !PAW: 1- Initialize values for several arrays depending only on atomic data 590 !2- Check overlap 591 !3- Identify FFT points in spheres and compute g_l(r).Y_lm(r) (and exp(-i.q.r) if needed) 592 !4- Allocate PAW specific arrays 593 !5- Compute perturbed local potential inside spheres 594 !6- Eventually open temporary storage files 595 if(psps%usepaw==1) then 596 ! 1-Initialize values for several arrays depending only on atomic data 597 gnt_option=1 598 if (dtset%pawxcdev==2.or.dtset%rfphon/=0.or.dtset%rfstrs/=0.or.dtset%rfelfd==1.or.& 599 dtset%rfelfd==3.or.dtset%rf2_dkde==1) gnt_option=2 600 601 ! Test if we have to call pawinit 602 call paw_gencond(Dtset,gnt_option,"test",call_pawinit) 603 604 if (psp_gencond==1.or.call_pawinit) then 605 ! Some gen-cond have to be added... 606 call timab(553,1,tsec) 607 call pawinit(gnt_option,zero,zero,dtset%pawlcutd,dtset%pawlmix,& 608 & psps%mpsang,dtset%pawnphi,dtset%nsym,dtset%pawntheta,& 609 & pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%xclevel,dtset%usepotzero) 610 call setsymrhoij(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,& 611 & rprimd,symrec,pawang%zarot) 612 613 ! Update internal values 614 call paw_gencond(Dtset,gnt_option,"save",call_pawinit) 615 616 call timab(553,2,tsec) 617 else 618 if (pawtab(1)%has_kij ==1) pawtab(1:psps%ntypat)%has_kij =2 619 if (pawtab(1)%has_nabla==1) pawtab(1:psps%ntypat)%has_nabla=2 620 end if 621 psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore) 622 call setsymrhoij(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,rprimd,symrec,pawang%zarot) 623 pawtab(:)%usepawu=0 624 pawtab(:)%useexexch=0 625 pawtab(:)%exchmix=zero 626 ! if (dtset%usepawu>0.or.dtset%useexexch>0) then 627 call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,& 628 & dtset%jpawu,dtset%lexexch,dtset%lpawu,ntypat,pawang,dtset%pawprtvol,pawrad,& 629 & pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu) 630 ! end if 631 compch_fft=-1.d5;compch_sph=-1.d5 632 usexcnhat=maxval(pawtab(:)%usexcnhat) 633 usecprj=dtset%pawusecp 634 ! 2-Check overlap 635 call chkpawovlp(natom,psps%ntypat,dtset%pawovlp,pawtab,rmet,dtset%typat,xred) 636 ! 3-Identify FFT points in spheres and compute g_l(r).Y_lm(r) and exp(-i.q.r) 637 ABI_DATATYPE_ALLOCATE(pawfgrtab,(my_natom)) 638 if (my_natom>0) then 639 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,& 640 & mpi_atmtab=mpi_enreg%my_atmtab) 641 call pawfgrtab_init(pawfgrtab,1,l_size_atm,pawrhoij(1)%nspden,dtset%typat,& 642 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 643 ABI_DEALLOCATE(l_size_atm) 644 end if 645 use_nhat_gga=(dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0) 646 optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm 647 if (use_nhat_gga) then 648 optgr1=dtset%pawstgylm 649 if (rfphon==1) optgr2=1 650 end if 651 if (rfphon==1.or.rfstrs/=0.or.rfelfd==3.or.rf2_dkde==1) then ! LB2016-11-28 : Why not rfelfd==1? 652 if (optgr1==0) optgr1=dtset%pawstgylm 653 if (optgr2==0) optgr2=dtset%pawstgylm 654 if (optrad==0.and.(.not.qeq0.or.rfstrs/=0.or.rfelfd==3.or.rf2_dkde==1)) optrad=1 ! LB2016-11-28 : Why not rfelfd==1? 655 end if 656 if (rfelfd==1.or.rfelfd==3.or.rf2_dkde==1) then 657 if (optgr1==0) optgr1=dtset%pawstgylm 658 end if 659 call status(0,dtfil%filstat,iexit,level,'call nhatgrid ') 660 call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfftf,psps%ntypat,& 661 & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,& 662 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab ) 663 ! Compute exp(iq.r) factors around the atoms 664 if (.not.qeq0) then 665 do iatom=1,my_natom 666 iatom_tot=iatom; if (paral_atom) iatom_tot=mpi_enreg%my_atmtab(iatom) 667 if (allocated(pawfgrtab(iatom)%expiqr)) then 668 ABI_DEALLOCATE(pawfgrtab(iatom)%expiqr) 669 end if 670 ABI_ALLOCATE(pawfgrtab(iatom)%expiqr,(2,pawfgrtab(iatom)%nfgd)) 671 call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,pawfgrtab(iatom)%nfgd,dtset%qptn,& 672 & pawfgrtab(iatom)%rfgd,xred(:,iatom_tot)) 673 pawfgrtab(iatom)%expiqr_allocated=1 674 end do 675 end if 676 ! 4-Allocate PAW specific arrays 677 ABI_DATATYPE_ALLOCATE(paw_an,(my_natom)) 678 ABI_DATATYPE_ALLOCATE(paw_ij,(my_natom)) 679 call paw_an_nullify(paw_an) 680 call paw_ij_nullify(paw_ij) 681 has_kxc=0;nkxc1=0;cplex=1 682 has_dijnd=0; req_cplex_dij=1 683 if(any(abs(dtset%nucdipmom)>tol8)) then 684 has_dijnd=1; req_cplex_dij=2 685 end if 686 if (rfphon/=0.or.rfelfd==1.or.rfelfd==3.or.rfstrs/=0.or.rf2_dkde/=0) then 687 has_kxc=1;nkxc1=2*min(dtset%nspden,2)-1 ! LDA 688 if(dtset%xclevel==2.and.dtset%nspden==1) nkxc1=7 ! GGA non-polarized 689 if(dtset%xclevel==2.and.dtset%nspden==2) nkxc1=19 ! GGA polarized 690 if (dtset%nspden==4) nkxc1=nkxc1+3 ! Non-coll.: need to store 3 additional arrays in kxc 691 end if 692 call paw_an_init(paw_an,dtset%natom,dtset%ntypat,nkxc1,dtset%nspden,& 693 & cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1,has_vxc_ex=1,has_kxc=has_kxc,& 694 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 695 call paw_ij_init(paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%pawspnorb,& 696 & natom,dtset%ntypat,dtset%typat,pawtab,has_dij=1,has_dijhartree=1,has_dijnd=has_dijnd,& 697 & has_dijso=1,has_pawu_occ=1,has_exexch_pot=1,req_cplex_dij=req_cplex_dij,& 698 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 699 700 else ! PAW vs NCPP 701 usexcnhat=0;usecprj=0 702 use_nhat_gga=.false. 703 ABI_DATATYPE_ALLOCATE(paw_an,(0)) 704 ABI_DATATYPE_ALLOCATE(paw_ij,(0)) 705 ABI_DATATYPE_ALLOCATE(pawfgrtab,(0)) 706 end if 707 708 ABI_ALLOCATE(rhog,(2,nfftf)) 709 ABI_ALLOCATE(rhor,(nfftf,dtset%nspden)) 710 711 !Read ground-state charge density from diskfile in case getden /= 0 712 !or compute it from wfs that were read previously : rhor as well as rhog 713 714 if (dtset%getden /= 0 .or. dtset%irdden /= 0) then 715 ! Read rho1(r) from a disk file and broadcast data. 716 ! This part is not compatible with MPI-FFT (note single_proc=.True. below) 717 718 rdwr=1;rdwrpaw=psps%usepaw;if(ireadwf0/=0) rdwrpaw=0 719 if (rdwrpaw/=0) then 720 ABI_DATATYPE_ALLOCATE(pawrhoij_read,(natom)) 721 call pawrhoij_nullify(pawrhoij_read) 722 nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb) 723 call pawrhoij_alloc(pawrhoij_read,dtset%pawcpxocc,nspden_rhoij,dtset%nspinor,& 724 & dtset%nsppol,dtset%typat,pawtab=pawtab) 725 else 726 ABI_DATATYPE_ALLOCATE(pawrhoij_read,(0)) 727 end if 728 729 ! MT july 2013: Should we read rhoij from the density file ? 730 call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor, & 731 hdr_den, pawrhoij_read, spaceworld, check_hdr=hdr) 732 etotal = hdr_den%etot; call hdr_free(hdr_den) 733 734 if (rdwrpaw/=0) then 735 call pawrhoij_bcast(pawrhoij_read,hdr%pawrhoij,0,spaceworld) 736 call pawrhoij_free(pawrhoij_read) 737 end if 738 ABI_DATATYPE_DEALLOCATE(pawrhoij_read) 739 740 ! Compute up+down rho(G) by fft 741 ABI_ALLOCATE(work,(nfftf)) 742 work(:)=rhor(:,1) 743 call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 744 ABI_DEALLOCATE(work) 745 746 else 747 izero=0 748 ! Obtain the charge density from read wfs 749 ! Be careful: in PAW, compensation density has to be added ! 750 tim_mkrho=4 751 paw_dmft%use_sc_dmft=0 ! respfn with dmft not implemented 752 paw_dmft%use_dmft=0 ! respfn with dmft not implemented 753 if (psps%usepaw==1) then 754 ABI_ALLOCATE(rhowfg,(2,dtset%nfft)) 755 ABI_ALLOCATE(rhowfr,(dtset%nfft,dtset%nspden)) 756 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,& 757 & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs) 758 call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor) 759 ABI_DEALLOCATE(rhowfg) 760 ABI_DEALLOCATE(rhowfr) 761 else 762 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,& 763 & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs) 764 end if 765 end if ! getden 766 767 !In PAW, compensation density has eventually to be added 768 nhatgrdim=0;nhatdim=0 769 ABI_ALLOCATE(nhatgr,(0,0,0)) 770 if (psps%usepaw==1.and. ((usexcnhat==0).or.(dtset%getden==0).or.dtset%xclevel==2)) then 771 nhatdim=1 772 ABI_ALLOCATE(nhat,(nfftf,dtset%nspden)) 773 call timab(558,1,tsec) 774 nhatgrdim=0;if (dtset%xclevel==2.and.dtset%pawnhatxc>0) nhatgrdim=usexcnhat 775 ider=2*nhatgrdim 776 if (nhatgrdim>0) then 777 ABI_DEALLOCATE(nhatgr) 778 ABI_ALLOCATE(nhatgr,(nfftf,dtset%nspden,3)) 779 end if 780 izero=0;cplex=1;ipert=0;idir=0;qphon(:)=zero 781 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,natom,& 782 & nfftf,ngfftf,nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,& 783 & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred, & 784 & mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom) 785 if (dtset%getden==0) then 786 rhor(:,:)=rhor(:,:)+nhat(:,:) 787 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 788 end if 789 call timab(558,2,tsec) 790 else 791 ABI_ALLOCATE(nhat,(0,0)) 792 end if 793 794 !The GS irrzon and phnons were only needed to symmetrize the GS density 795 ABI_DEALLOCATE(irrzon) 796 ABI_DEALLOCATE(phnons) 797 798 !jmb 2012 write(std_out,'(a)')' ' ! needed to make ibm6_xlf12 pass tests. No idea why this works. JWZ 5 Sept 2011 799 !Will compute now the total potential 800 801 !Compute local ionic pseudopotential vpsp and core electron density xccc3d: 802 n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf 803 ABI_ALLOCATE(xccc3d,(n3xccc)) 804 ABI_ALLOCATE(vpsp,(nfftf)) 805 806 !Determine by which method the local ionic potential and/or 807 ! the pseudo core charge density have to be computed 808 !Local ionic potential: 809 ! Method 1: PAW ; Method 2: Norm-conserving PP 810 vloc_method=1;if (psps%usepaw==0) vloc_method=2 811 !Pseudo core charge density: 812 ! Method 1: PAW, nc_xccc_gspace ; Method 2: Norm-conserving PP 813 coredens_method=1;if (psps%usepaw==0) coredens_method=2 814 if (psps%nc_xccc_gspace==1) coredens_method=1 815 if (psps%nc_xccc_gspace==0) coredens_method=2 816 817 !Local ionic potential and/or pseudo core charge by method 1 818 if (vloc_method==1.or.coredens_method==1) then 819 call timab(562,1,tsec) 820 optv=0;if (vloc_method==1) optv=1 821 optn=0;if (coredens_method==1) optn=n3xccc/nfftf 822 optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1 823 call atm2fft(atindx1,xccc3d,vpsp,dum_dyfrn,dum_dyfrv,dum_eltfrxc,dum_gauss,gmet,gprimd,& 824 & dum_grn,dum_grv,gsqcut,mgfftf,psps%mqgrid_vl,natom,nattyp,nfftf,ngfftf,& 825 & ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1df,psps%qgrid_vl,& 826 & dtset%qprtrb,dum_rhog,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dum_vg,dum_vg,dum_vg,dtset%vprtrb,psps%vlspl) 827 call timab(562,2,tsec) 828 end if 829 830 !Local ionic potential by method 2 831 if (vloc_method==2) then 832 option=1 833 ABI_ALLOCATE(dyfrlo_indx,(3,3,natom)) 834 ABI_ALLOCATE(grtn_indx,(3,natom)) 835 call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,& 836 & grtn_indx,gsqcut,dummy6,mgfftf,mpi_enreg,natom,nattyp,& 837 & nfftf,ngfftf,dtset%nspden,ntypat,option,pawtab,ph1df,psps,& 838 & dtset%qprtrb,rhog,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred) 839 ABI_DEALLOCATE(dyfrlo_indx) 840 ABI_DEALLOCATE(grtn_indx) 841 end if 842 843 !Pseudo core electron density by method 2 844 if (coredens_method==2.and.psps%n1xccc/=0) then 845 option=1 846 ABI_ALLOCATE(dyfrx2,(3,3,natom)) 847 ABI_ALLOCATE(vxc,(0,0)) ! dummy 848 call mkcore(dummy6,dyfrx2,grxc,mpi_enreg,natom,nfftf,dtset%nspden,ntypat,& 849 & ngfftf(1),psps%n1xccc,ngfftf(2),ngfftf(3),option,rprimd,dtset%typat,ucvol,vxc,& 850 & psps%xcccrc,psps%xccc1d,xccc3d,xred) 851 ABI_DEALLOCATE(dyfrx2) 852 ABI_DEALLOCATE(vxc) ! dummy 853 end if 854 855 !Set up hartree and xc potential. Compute kxc here. 856 ABI_ALLOCATE(vhartr,(nfftf)) 857 858 call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,rhog,rprimd,vhartr) 859 860 option=2 ; nk3xc=1 861 nkxc=2*min(dtset%nspden,2)-1;if(dtset%xclevel==2)nkxc=12*min(dtset%nspden,2)-5 862 call check_kxc(dtset%ixc,dtset%optdriver) 863 ABI_ALLOCATE(kxc,(nfftf,nkxc)) 864 ABI_ALLOCATE(vxc,(nfftf,dtset%nspden)) 865 866 _IBM6("Before rhotoxc") 867 868 call xcdata_init(xcdata,dtset=dtset) 869 call rhotoxc(enxc,kxc,mpi_enreg,nfftf,ngfftf,& 870 & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor,& 871 & rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,vhartr=vhartr) 872 873 !Compute local + Hxc potential, and subtract mean potential. 874 ABI_ALLOCATE(vtrial,(nfftf,dtset%nspden)) 875 do ispden=1,min(dtset%nspden,2) 876 do ifft=1,nfftf 877 vtrial(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft) 878 end do 879 end do 880 if (dtset%nspden==4) then 881 do ispden=3,4 882 do ifft=1,nfftf 883 vtrial(ifft,ispden)=vxc(ifft,ispden) 884 end do 885 end do 886 end if 887 ABI_DEALLOCATE(vhartr) 888 889 890 if(dtset%prtvol==-level)then 891 call wrtout(std_out,' respfn: ground-state density and potential set up.','COLL') 892 end if 893 894 !PAW: compute Dij quantities (psp strengths) 895 if (psps%usepaw==1)then 896 cplex=1;ipert=0;option=1 897 nzlmopt=0;if (dtset%pawnzlm>0) nzlmopt=-1 898 899 call pawdenpot(compch_sph,epaw,epawdc,ipert,dtset%ixc,my_natom,natom,dtset%nspden,& 900 & ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang,dtset%pawprtvol,& 901 & pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,& 902 & dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp, & 903 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 904 905 call timab(561,1,tsec) 906 call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,natom,nfftf,nfftotf,& 907 & dtset%nspden,ntypat,paw_an,paw_ij,pawang,pawfgrtab,dtset%pawprtvol,& 908 & pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,k0,& 909 & dtset%spnorbscl,ucvol,dtset%charge,vtrial,vxc,xred,nucdipmom=dtset%nucdipmom,& 910 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 911 call symdij(gprimd,indsym,ipert,my_natom,natom,dtset%nsym,ntypat,0,& 912 & paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,& 913 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 914 call timab(561,2,tsec) 915 916 end if 917 918 !-----2. Frozen-wavefunctions and Ewald(q=0) parts of 2DTE 919 920 dyfr_nondiag=0;if (psps%usepaw==1.and.rfphon==1) dyfr_nondiag=1 921 dyfr_cplex=1;if (psps%usepaw==1.and.rfphon==1.and.(.not.qeq0)) dyfr_cplex=2 922 ABI_ALLOCATE(dyew,(2,3,natom,3,natom)) 923 ABI_ALLOCATE(dyewq0,(3,3,natom)) 924 ABI_ALLOCATE(dyfrlo,(3,3,natom)) 925 ABI_ALLOCATE(dyfrx2,(3,3,natom)) 926 ABI_ALLOCATE(dyfrnl,(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)) 927 ABI_ALLOCATE(dyfrwf,(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)) 928 ABI_ALLOCATE(becfrnl,(3,natom,3*pawbec)) 929 ABI_ALLOCATE(piezofrnl,(6,3*pawpiezo)) 930 ABI_ALLOCATE(dyvdw,(2,3,natom,3,natom*usevdw)) 931 dyew(:,:,:,:,:)=zero 932 dyewq0(:,:,:)=zero 933 dyfrnl(:,:,:,:,:)=zero 934 dyfrwf(:,:,:,:,:)=zero 935 dyfrlo(:,:,:)=zero 936 dyfrx2(:,:,:)=zero 937 if (usevdw==1) dyvdw(:,:,:,:,:)=zero 938 if (pawbec==1) becfrnl(:,:,:)=zero 939 if (pawpiezo==1) piezofrnl(:,:)=zero 940 941 ABI_ALLOCATE(eltcore,(6,6)) 942 ABI_ALLOCATE(elteew,(6+3*natom,6)) 943 ABI_ALLOCATE(eltfrhar,(6,6)) 944 ABI_ALLOCATE(eltfrnl,(6+3*natom,6)) 945 ABI_ALLOCATE(eltfrloc,(6+3*natom,6)) 946 ABI_ALLOCATE(eltfrkin,(6,6)) 947 ABI_ALLOCATE(eltfrxc,(6+3*natom,6)) 948 ABI_ALLOCATE(eltvdw,(6+3*natom,6*usevdw)) 949 eltcore(:,:)=zero 950 elteew(:,:)=zero 951 eltfrnl(:,:)=zero 952 eltfrloc(:,:)=zero 953 eltfrkin(:,:)=zero 954 eltfrhar(:,:)=zero 955 eltfrxc(:,:)=zero 956 if (usevdw==1) eltvdw(:,:)=zero 957 958 !Section common to all perturbations 959 !Compute the nonlocal part of the elastic tensor and/or dynamical matrix 960 if (rfstrs/=0.or.rfphon==1.or.dtset%efmas>0.or.pawbec==1.or.pawpiezo==1)then 961 call d2frnl(becfrnl,cg,dtfil,dtset,dyfrnl,dyfr_cplex,& 962 & dyfr_nondiag,efmasdeg,efmasfr,eigen0,eltfrnl,gsqcut,has_allddk,indsym,kg,mgfftf,& 963 & mpi_enreg,psps%mpsang,my_natom,natom,nfftf,ngfft,ngfftf,& 964 & npwarr,occ,paw_ij,pawang,pawbec,pawfgrtab,pawpiezo,pawrad,& 965 & pawrhoij,pawtab,ph1d,ph1df,piezofrnl,psps,rprimd,rfphon,& 966 & rfstrs,symrec,vtrial,vxc,xred,ylm,ylmgr) 967 end if 968 969 !No more need of these local derivatives 970 if (rfphon==1.and.psps%usepaw==1.and.(.not.use_nhat_gga)) then 971 do iatom=1,my_natom 972 if (allocated(pawfgrtab(iatom)%gylmgr2)) then 973 ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2) 974 end if 975 pawfgrtab(iatom)%gylmgr2_allocated=0 976 end do 977 end if 978 979 !Section for the atomic displacement/electric field perturbations 980 if (rfphon==1) then 981 982 ! Compute the local of the dynamical matrix 983 ! dyfrnl has not yet been symmetrized, but will be in the next routine 984 call dfpt_dyfro(atindx1,dyfrnl,dyfrlo,dyfrwf,dyfrx2,dyfr_cplex,dyfr_nondiag,& 985 & gmet,gprimd,gsqcut,indsym,mgfftf,mpi_enreg,psps%mqgrid_vl,& 986 & natom,nattyp, nfftf,ngfftf,dtset%nspden,dtset%nsym,ntypat,& 987 & psps%n1xccc,n3xccc,dtset%paral_kgb,psps,pawtab,ph1df,psps%qgrid_vl,& 988 & dtset%qptn,rhog,rprimd,symq,symrec,dtset%typat,ucvol,& 989 & psps%usepaw,psps%vlspl,vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred) 990 991 _IBM6("Before dfpt_ewald") 992 993 ! Compute Ewald (q=0) contribution 994 sumg0=0;qphon(:)=zero 995 call status(0,dtfil%filstat,iexit,level,'call dfpt_ewald(1)') 996 call dfpt_ewald(dyew,gmet,my_natom,natom,qphon,rmet,sumg0,dtset%typat,ucvol,xred,psps%ziontypat,& 997 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 998 option=1 999 call q0dy3_calc(natom,dyewq0,dyew,option) 1000 ! Calculate the DFT-D2 vdw part of the dynamical matrix 1001 if (usevdw==1.and.dtset%vdw_xc==5) then 1002 call vdw_dftd2(evdw,dtset%ixc,natom,ntypat,1,dtset%typat,rprimd,dtset%vdw_tol,& 1003 & xred,psps%znuclpsp,dyn_vdw_dftd2=dyvdw,qphon=dtset%qptn) 1004 end if 1005 ! Calculate the DFT-D3/D3(BJ) vdw part of the dynamical matrix 1006 if (usevdw==1.and.(dtset%vdw_xc==6.or.dtset%vdw_xc==7)) then 1007 call vdw_dftd3(evdw,dtset%ixc,natom,ntypat,1,dtset%typat,rprimd,& 1008 & dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,xred,psps%znuclpsp,& 1009 & dyn_vdw_dftd3=dyvdw,qphon=dtset%qptn) 1010 end if 1011 !The frozen-wavefunction part of the dynamical matrix is now: 1012 ! d2frnl: non-local contribution 1013 ! dyfrlo: local contribution 1014 ! dyfrx2: 2nd order xc core correction contribution 1015 ! dyew : Ewald contribution 1016 ! dyvdw : vdw DFT-D contribution 1017 ! dyfrwf: all contributions 1018 ! In case of PAW, it misses a term coming from the perturbed overlap operator 1019 end if 1020 1021 !Section for the strain perturbation 1022 if(rfstrs/=0) then 1023 1024 ! Verify that k-point set has full space-group symmetry; otherwise exit 1025 call status(0,dtfil%filstat,iexit,level,'call symkchk ') 1026 timrev=1 1027 if (symkchk(dtset%kptns,dtset%nkpt,dtset%nsym,symrec,timrev,message) /= 0) then 1028 MSG_ERROR(message) 1029 end if 1030 1031 ! Calculate the kinetic part of the elastic tensor 1032 call dfpt_eltfrkin(cg,eltfrkin,dtset%ecut,dtset%ecutsm,dtset%effmass_free,& 1033 & dtset%istwfk,kg,dtset%kptns,dtset%mband,dtset%mgfft,dtset%mkmem,mpi_enreg,& 1034 & dtset%mpw,dtset%nband,dtset%nkpt,ngfft,npwarr,& 1035 & dtset%nspinor,dtset%nsppol,occ,rprimd,dtset%wtk) 1036 1037 ! Calculate the hartree part of the elastic tensor 1038 call dfpt_eltfrhar(eltfrhar,rprimd,gsqcut,mpi_enreg,nfftf,ngfftf,rhog) 1039 1040 ! Calculate the xc part of the elastic tensor 1041 call dfpt_eltfrxc(atindx,dtset,eltfrxc,enxc,gsqcut,kxc,mpi_enreg,mgfftf,& 1042 & nattyp,nfftf,ngfftf,ngfftf,nhat,nkxc,n3xccc,pawtab,ph1df,psps,rhor,rprimd,& 1043 & usexcnhat,vxc,xccc3d,xred) 1044 1045 ! Calculate the local potential part of the elastic tensor 1046 call dfpt_eltfrloc(atindx,eltfrloc,gmet,gprimd,gsqcut,mgfftf,mpi_enreg,psps%mqgrid_vl,& 1047 & natom,nattyp,nfftf,ngfftf,ntypat,ph1df,psps%qgrid_vl,rhog,psps%vlspl) 1048 1049 ! Calculate the Ewald part of the elastic tensor 1050 call elt_ewald(elteew,gmet,gprimd,my_natom,natom,ntypat,rmet,rprimd,& 1051 & dtset%typat,ucvol,xred,psps%ziontypat,& 1052 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1053 1054 ! Calculate the DFT-D2 vdw part of the elastic tensor 1055 if (usevdw==1.and.dtset%vdw_xc==5) then 1056 option=1; if (rfphon==1) option=0 1057 call vdw_dftd2(evdw,dtset%ixc,natom,ntypat,option,dtset%typat,rprimd,dtset%vdw_tol,& 1058 & xred,psps%znuclpsp,elt_vdw_dftd2=eltvdw) 1059 end if 1060 ! Calculate the DFT-D3/D3(BJ) vdw part of the elastic tensor 1061 if (usevdw==1.and.(dtset%vdw_xc==6.or.dtset%vdw_xc==7)) then 1062 option=1; if (rfphon==1) option=0 1063 call vdw_dftd3(evdw,dtset%ixc,natom,ntypat,option,dtset%typat,rprimd,& 1064 & dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,xred,psps%znuclpsp,& 1065 & elt_vdw_dftd3=eltvdw) 1066 end if 1067 ! Calculate the psp core energy part of elastic tensor (trivial) 1068 eltcore(1:3,1:3)=ecore/ucvol 1069 1070 !The frozen-wavefunction part of the elastic tensor is now: 1071 ! eltfrnl: non-local contribution 1072 ! eltfrloc: local contribution 1073 ! eltfrkin: kinetic contribution 1074 ! eltfrhar: Hartree contribution 1075 ! eltfrx: XC contribution 1076 ! eltcore: psps core contribution 1077 ! elteew: Ewald contribution 1078 ! eltvdw: vdw DFT-D contribution 1079 ! In case of PAW, it misses a term coming from the perturbed overlap operator 1080 end if 1081 1082 ABI_DEALLOCATE(vpsp) 1083 ABI_DEALLOCATE(xccc3d) 1084 1085 if(dtset%prtvol==-level)then 1086 call wrtout(std_out,' respfn: frozen wavef. and Ewald(q=0) part of 2DTE done.','COLL') 1087 end if 1088 1089 call timab(136,2,tsec) 1090 1091 !-----3. Initialisation of 1st response, taking into account the q vector. 1092 1093 call timab(137,1,tsec) 1094 1095 write(message,'(3a)')ch10,' ==> initialize data related to q vector <== ',ch10 1096 call wrtout(std_out,message,'COLL') 1097 call wrtout(ab_out,message,'COLL') 1098 1099 qphon(:)=dtset%qptn(:) 1100 sumg0=1 1101 1102 !Treat the case of q=0 or q too close to 0 1103 qzero=0 1104 if(qeq0)then 1105 qphon(:)=zero 1106 write(message,'(3a)')& 1107 & ' respfn : the norm of the phonon wavelength (as input) was small (<1.d-7).',ch10,& 1108 & ' q has been set exactly to (0 0 0)' 1109 call wrtout(std_out,message,'COLL') 1110 sumg0=0 1111 qzero=1 1112 else 1113 if(rfelfd/=0 .or. rfstrs/=0 .or. rfddk /= 0 .or. rf2_dkdk /= 0 .or. rf2_dkde /= 0) then 1114 ! Temporarily, ... 1115 write(message, '(a,a,a,3es16.6,a,a,6(a,i2),a,a,a)' )ch10,& 1116 & 'The treatment of non-zero wavevector q is restricted to phonons.',& 1117 & 'However, the input normalized qpt is',qphon(:),',',ch10,& 1118 & 'while rfelfd=',rfelfd,', rfddk=',rfddk,', rf2_dkdk=',rf2_dkdk,', rf2_dkde=',rf2_dkde,& 1119 & ' and rfstrs=',rfstrs,'.',ch10,& 1120 & 'Action: change qpt, or rfelfd, or rfstrs in the input file.' 1121 MSG_ERROR(message) 1122 else if(rfasr.eq.2)then 1123 write(message,'(2a)')ch10,' rfasr=2 not allowed with q/=0 => rfasr was reset to 0.' 1124 MSG_WARNING(message) 1125 rfasr=0 1126 end if 1127 end if 1128 1129 _IBM6("Before irreducible_set_pert") 1130 1131 !Determine the symmetrical perturbations 1132 ABI_ALLOCATE(pertsy,(3,natom+6)) 1133 call irreducible_set_pert(indsym,natom+6,natom,dtset%nsym,pertsy,rfdir,rfpert,symq,symrec,dtset%symrel) 1134 write(message,'(a)') ' The list of irreducible perturbations for this q vector is:' 1135 call wrtout(ab_out,message,'COLL') 1136 call wrtout(std_out,message,'COLL') 1137 ii=1 1138 do ipert=1,natom+6 1139 do idir=1,3 1140 if(rfpert(ipert)==1.and.rfdir(idir)==1)then 1141 if( pertsy(idir,ipert)==1 )then 1142 write(message, '(i5,a,i2,a,i4)' )ii,') idir=',idir,' ipert=',ipert 1143 call wrtout(ab_out,message,'COLL') 1144 call wrtout(std_out,message,'COLL') 1145 ii=ii+1 1146 end if 1147 end if 1148 end do 1149 end do 1150 1151 !test if the user left default rfdir 0 0 0 1152 if (ii==1 .and. rf2_dkdk==0 .and. rf2_dkde==0) then 1153 write(message,'(5a)')ch10,& 1154 & ' WARNING: no perturbations to be done at this q-point.',ch10,& 1155 & ' You may have forgotten to set the rfdir or rfatpol variables. Continuing normally.',ch10 1156 call wrtout(ab_out,message,'COLL') 1157 MSG_WARNING(message) 1158 end if 1159 1160 if (dtset%prepanl==1.and.(rf2_dkdk/=0 .or. rf2_dkde/=0)) then 1161 ABI_ALLOCATE(rfpert_nl,(3,natom+2,3,natom+2,3,natom+2)) 1162 rfpert_nl = 0 1163 rfpert_nl(:,natom+2,:,natom+2,:,natom+2) = 1 1164 rfpert_nl(:,1:natom,:,natom+2,:,natom+2) = 1 1165 rfpert_nl(:,natom+2,:,1:natom,:,natom+2) = 1 1166 rfpert_nl(:,natom+2,:,natom+2,:,1:natom) = 1 1167 call sytens(indsym,natom+2,natom,dtset%nsym,rfpert_nl,symrec,dtset%symrel) 1168 write(message, '(a,a,a,a,a)' ) ch10, & 1169 & ' The list of irreducible elements of the Raman and non-linear',& 1170 & ch10,' optical susceptibility tensors is:',ch10 1171 call wrtout(std_out,message,'COLL') 1172 1173 write(message,'(12x,a)')& 1174 & 'i1pert i1dir i2pert i2dir i3pert i3dir' 1175 call wrtout(std_out,message,'COLL') 1176 n1 = 0 1177 rf2_dirs_from_rfpert_nl(:,:) = 0 1178 do i1pert = 1, natom + 2 1179 do i1dir = 1, 3 1180 do i2pert = 1, natom + 2 1181 do i2dir = 1, 3 1182 do i3pert = 1, natom + 2 1183 do i3dir = 1,3 1184 if (rfpert_nl(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 1185 n1 = n1 + 1 1186 write(message,'(2x,i4,a,6(5x,i3))') n1,')', & 1187 & i1pert,i1dir,i2pert,i2dir,i3pert,i3dir 1188 call wrtout(std_out,message,'COLL') 1189 if (i2pert==natom+2) then 1190 if (i3pert==natom+2) then 1191 rf2_dirs_from_rfpert_nl(i3dir,i2dir) = 1 1192 else if (i1pert==natom+2) then 1193 rf2_dirs_from_rfpert_nl(i1dir,i2dir) = 1 1194 end if 1195 end if 1196 end if 1197 end do 1198 end do 1199 end do 1200 end do 1201 end do 1202 end do 1203 write(message,'(a,a)') ch10,ch10 1204 call wrtout(std_out,message,'COLL') 1205 1206 write(message,'(a)') 'rf2_dirs_from_rfpert_nl :' 1207 call wrtout(std_out,message,'COLL') 1208 do i1dir = 1, 3 1209 do i2dir = 1, 3 1210 write(message,'(3(a,i1))') ' ',i1dir,' ',i2dir,' : ',rf2_dirs_from_rfpert_nl(i1dir,i2dir) 1211 call wrtout(std_out,message,'COLL') 1212 end do 1213 end do 1214 end if 1215 1216 !Contribution to the dynamical matrix from ion-ion energy 1217 if(rfphon==1)then 1218 call status(0,dtfil%filstat,iexit,level,'call dfpt_ewald(2)') 1219 call dfpt_ewald(dyew,gmet,my_natom,natom,qphon,rmet,sumg0,dtset%typat,ucvol,xred,psps%ziontypat, & 1220 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1221 call q0dy3_apply(natom,dyewq0,dyew) 1222 end if 1223 1224 !1-order contribution of the xc core correction to the dynamical matrix 1225 ABI_ALLOCATE(dyfrx1,(2,3,natom,3,natom)) 1226 dyfrx1(:,:,:,:,:)=zero 1227 if(rfphon==1.and.psps%n1xccc/=0)then 1228 ABI_ALLOCATE(blkflgfrx1,(3,natom,3,natom)) 1229 !FR non-collinear magnetism 1230 if (dtset%nspden==4) then 1231 call dfpt_dyxc1(atindx,blkflgfrx1,dyfrx1,gmet,gsqcut,dtset%ixc,kxc,mgfftf,mpert,mpi_enreg,& 1232 & psps%mqgrid_vl,natom,nfftf,ngfftf,nkxc,dtset%nspden,& 1233 & ntypat,psps%n1xccc,dtset%paral_kgb,psps,pawtab,ph1df,psps%qgrid_vl,qphon,& 1234 & rfdir,rfpert,rprimd,timrev,dtset%typat,ucvol,psps%usepaw,psps%xcccrc,psps%xccc1d,xred,rhor=rhor,vxc=vxc) 1235 else 1236 call dfpt_dyxc1(atindx,blkflgfrx1,dyfrx1,gmet,gsqcut,dtset%ixc,kxc,mgfftf,mpert,mpi_enreg,& 1237 & psps%mqgrid_vl,natom,nfftf,ngfftf,nkxc,dtset%nspden,& 1238 & ntypat,psps%n1xccc,dtset%paral_kgb,psps,pawtab,ph1df,psps%qgrid_vl,qphon,& 1239 & rfdir,rfpert,rprimd,timrev,dtset%typat,ucvol,psps%usepaw,psps%xcccrc,psps%xccc1d,xred) 1240 end if 1241 end if 1242 1243 !Deallocate the arrays that were needed only for the frozen wavefunction part 1244 ABI_DEALLOCATE(ph1d) 1245 ABI_DEALLOCATE(ph1df) 1246 ABI_DEALLOCATE(cg) 1247 ABI_DEALLOCATE(kg) 1248 ABI_DEALLOCATE(npwarr) 1249 if(xmpi_paral==1) then 1250 ABI_DEALLOCATE(mpi_enreg%proc_distrb) 1251 end if 1252 1253 ABI_ALLOCATE(blkflg,(3,mpert,3,mpert)) 1254 ABI_ALLOCATE(d2eig0,(2,3,mpert,3,mpert)) 1255 ABI_ALLOCATE(d2k0,(2,3,mpert,3,mpert)) 1256 ABI_ALLOCATE(d2lo,(2,3,mpert,3,mpert)) 1257 ABI_ALLOCATE(d2loc0,(2,3,mpert,3,mpert)) 1258 ABI_ALLOCATE(d2nfr,(2,3,mpert,3,mpert)) 1259 ABI_ALLOCATE(d2nl,(2,3,mpert,3,mpert)) 1260 ABI_ALLOCATE(d2nl0,(2,3,mpert,3,mpert)) 1261 ABI_ALLOCATE(d2nl1,(2,3,mpert,3,mpert)) 1262 ABI_ALLOCATE(d2vn,(2,3,mpert,3,mpert)) 1263 ABI_ALLOCATE(d2ovl,(2,3,mpert,3,mpert*psps%usepaw)) 1264 blkflg(:,:,:,:)=0 1265 d2eig0(:,:,:,:,:)=zero ; d2k0(:,:,:,:,:)=zero 1266 d2lo(:,:,:,:,:)=zero ; d2loc0(:,:,:,:,:)=zero 1267 d2nfr(:,:,:,:,:)=zero ; d2nl(:,:,:,:,:)=zero 1268 d2nl0(:,:,:,:,:)=zero ; d2nl1(:,:,:,:,:)=zero 1269 d2vn(:,:,:,:,:)=zero 1270 if (psps%usepaw==1) d2ovl(:,:,:,:,:)=zero 1271 1272 prtbbb=dtset%prtbbb 1273 ABI_ALLOCATE(d2bbb,(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)) 1274 ABI_ALLOCATE(d2cart_bbb,(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)) 1275 if(prtbbb==1)then 1276 d2cart_bbb(:,:,:,:,:,:)=zero 1277 d2bbb(:,:,:,:,:,:)=zero 1278 end if 1279 1280 dim_eig2nkq = 0 1281 if(dtset%ieig2rf /= 0) dim_eig2nkq = 1 1282 ABI_ALLOCATE(eig2nkq,(2,dtset%mband*dtset%nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq)) 1283 dim_eigbrd=0 1284 if(dtset%ieig2rf /= 0 .and. dtset%smdelta>0 ) dim_eigbrd = 1 1285 ABI_ALLOCATE(eigbrd,(2,dtset%mband*dtset%nsppol,dtset%nkpt,3,natom,3,natom*dim_eigbrd)) 1286 1287 call timab(137,2,tsec) 1288 1289 1290 !Check whether exiting was required by the user. 1291 !If found then do not start minimization steps 1292 !At this first call to exit_check, initialize cpus 1293 cpus=dtset%cpus 1294 if(abs(cpus)>1.0d-5)cpus=cpus+cpui 1295 openexit=1; if(dtset%chkexit==0) openexit=0 1296 call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit) 1297 1298 !TEMPORARY: for testing purpose only 1299 ! if (rfstrs/=0.and.dtset%usepaw==1) iexit=1 1300 1301 _IBM6("Before dfpt_looppert") 1302 1303 if (iexit==0) then 1304 ! ####################################################################### 1305 write(message,'(a,80a)')ch10,('=',mu=1,80) 1306 call wrtout(ab_out,message,'COLL') 1307 call wrtout(std_out,message,'COLL') 1308 1309 ddkfil(:)=0 1310 1311 ! MGNAG: WHY THIS? all these variables are declared as optional pointers in dfpt_looppert! 1312 ! but they are allocated here so why using pointers! Moreover OPTIONAL arguments MUST 1313 ! be passed by keyword for better clarity and robustness! 1314 ! People should learn how to program Fort90 before being allowed to change the code! 1315 ! v5[26] crashes in dfpt_looppert 1316 ! The best solution would be using a datatype to gather the results! 1317 ABI_ALLOCATE(clflg,(3,mpert)) 1318 if(dtset%ieig2rf > 0) then 1319 ABI_ALLOCATE(eigen0_pert,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1320 ABI_ALLOCATE(eigenq_pert,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1321 ABI_ALLOCATE(occ_rbz_pert,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1322 end if 1323 if(dtset%efmas > 0) then 1324 ABI_ALLOCATE(eigen0_pert,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1325 end if 1326 ! Note that kg, cg, eigen0, mpw and npwarr are NOT passed to dfpt_looppert : 1327 ! they will be reinitialized for each perturbation, with an eventually 1328 ! reduced set of k point, thanks to the use of symmetry operations. 1329 call dfpt_looppert(atindx,blkflg,codvsn,cpus,dim_eigbrd,dim_eig2nkq,doccde,& 1330 & ddkfil,dtfil,dtset,dyew,dyfrlo,dyfrnl,dyfrx1,dyfrx2,dyvdw,dyfr_cplex,dyfr_nondiag,& 1331 & d2bbb,d2lo,d2nl,d2ovl,efmasdeg,efmasfr,eigbrd,eig2nkq,& 1332 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,& 1333 & etotal,fermie,iexit,indsym,kxc,& 1334 & dtset%mkmem,mkqmem,mk1mem,mpert,mpi_enreg,my_natom,nattyp,& 1335 & nfftf,nhat,dtset%nkpt,nkxc,dtset%nspden,dtset%nsym,occ,& 1336 & paw_an,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,& 1337 & pertsy,prtbbb,psps,rfpert,rf2_dirs_from_rfpert_nl,rhog,rhor,symq,symrec,timrev,& 1338 & usecprj,usevdw,vtrial,vxc,vxcavg,xred,clflg,occ_rbz_pert,eigen0_pert,eigenq_pert,& 1339 & eigen1_pert,nkpt_rbz,eigenq_fine,hdr_fine,hdr0) 1340 1341 ! ##################################################################### 1342 end if 1343 1344 call timab(138,1,tsec) 1345 1346 write(message, '(80a,a,a,a,a)' ) ('=',mu=1,80),ch10,ch10,& 1347 & ' ---- first-order wavefunction calculations are completed ----',ch10 1348 call wrtout(ab_out,message,'COLL') 1349 call wrtout(std_out,message,'COLL') 1350 1351 ABI_DEALLOCATE(vxc) 1352 1353 if (dtset%prepanl==1.and.(rf2_dkdk/=0 .or. rf2_dkde/=0)) then 1354 ABI_DEALLOCATE(rfpert_nl) 1355 end if 1356 1357 !Output of the localization tensor 1358 if ( rfpert(natom+1) /= 0 .and. (me == 0) .and. dtset%occopt<=2) then 1359 call wrtloctens(blkflg,d2bbb,d2nl,dtset%mband,mpert,natom,dtset%prtbbb,rprimd,psps%usepaw) 1360 end if 1361 1362 !The perturbation natom+1 was only an auxiliary perturbation, 1363 !needed to construct the electric field response, so its flag is now set to 0. 1364 !rfpert(natom+1)=0 1365 1366 !Were 2DTE computed ? 1367 if(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0 .or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 .and. rfuser==0 .and. rfmagn==0)then 1368 1369 write(message,'(a,a)' )ch10,' respfn : d/dk was computed, but no 2DTE, so no DDB output.' 1370 call wrtout(std_out,message,'COLL') 1371 call wrtout(ab_out,message,'COLL') 1372 1373 ! If 2DTE were computed, only one processor must output them and compute 1374 ! frequencies. 1375 else if(me==0)then 1376 1377 write(message,'(a,a)' )ch10,' ==> Compute Derivative Database <== ' 1378 call wrtout(std_out,message,'COLL') 1379 call wrtout(ab_out,message,'COLL') 1380 1381 ! In the RESPFN code, dfpt_nstdy and stady3 were called here 1382 d2nfr(:,:,:,:,:)=d2lo(:,:,:,:,:)+d2nl(:,:,:,:,:) 1383 if (psps%usepaw==1) d2nfr(:,:,:,:,:)=d2nfr(:,:,:,:,:)+d2ovl(:,:,:,:,:) 1384 1385 ! In case of bbb decomposition 1386 if(prtbbb==1)then 1387 ABI_ALLOCATE(blkflg1,(3,mpert,3,mpert)) 1388 ABI_ALLOCATE(blkflg2,(3,mpert,3,mpert)) 1389 blkflg2(:,:,:,:) = blkflg(:,:,:,:) 1390 do ipert = 1, mpert 1391 do ipert2 = 1, mpert 1392 if ((ipert /= natom + 2).and.(ipert>natom).and.(ipert2/=natom+2)) then 1393 blkflg2(:,ipert2,:,ipert) = 0 1394 end if 1395 end do 1396 end do 1397 ABI_ALLOCATE(d2tmp,(2,3,mpert,3,mpert)) 1398 do iband = 1,dtset%mband 1399 d2tmp(:,:,:,:,:)=zero 1400 blkflg1(:,:,:,:) = blkflg2(:,:,:,:) 1401 d2tmp(:,:,natom+2,:,:) = d2bbb(:,:,:,:,iband,iband) 1402 call d2sym3(blkflg1,d2tmp,indsym,mpert,natom,dtset%nsym,qphon,symq,& 1403 & symrec,dtset%symrel,timrev) 1404 d2bbb(:,:,:,:,iband,iband) = d2tmp(:,:,natom+2,:,:) 1405 end do 1406 ABI_DEALLOCATE(blkflg1) 1407 ABI_DEALLOCATE(blkflg2) 1408 ABI_DEALLOCATE(d2tmp) 1409 end if 1410 1411 ! Complete the d2nfr matrix by symmetrization of the existing elements 1412 call d2sym3(blkflg,d2nfr,indsym,mpert,natom,dtset%nsym,qphon,symq,symrec,dtset%symrel,timrev) 1413 1414 if(rfphon==1.and.psps%n1xccc/=0)then 1415 ! Complete the dyfrx1 matrix by symmetrization of the existing elements 1416 call d2sym3(blkflgfrx1,dyfrx1,indsym,natom,natom,dtset%nsym,qphon,symq,symrec,dtset%symrel,timrev) 1417 end if 1418 1419 ! Note that there is a bug in d2sym3 which will set some elements of 1420 ! blkflg to 1 even when no corresponding symmetry-related element 1421 ! has been computed. This has the effect of producing spurious extra 1422 ! output lines in the 2nd-order matrix listing in the .out file 1423 ! and in the DDB file. The suprious matrix elements are all zero, 1424 ! so this is primarily an annoyance.(DRH) 1425 1426 1427 ! Add the frozen-wf (dyfrwf) part to the ewald part (dyew), 1428 ! the part 1 of the frozen wf part of the xc core correction 1429 ! (dyfrx1) and the non-frozen part (dynfr) to get the second-order 1430 ! derivative matrix (d2matr), then 1431 ! take account of the non-cartesian coordinates (d2cart). 1432 ABI_ALLOCATE(d2cart,(2,3,mpert,3,mpert)) 1433 ABI_ALLOCATE(carflg,(3,mpert,3,mpert)) 1434 ABI_ALLOCATE(d2matr,(2,3,mpert,3,mpert)) 1435 outd2=1 1436 call status(0,dtfil%filstat,iexit,level,'call dfpt_gatherdy ') 1437 call dfpt_gatherdy(becfrnl,dtset%berryopt,blkflg,carflg,& 1438 & dyew,dyfrwf,dyfrx1,dyfr_cplex,dyfr_nondiag,dyvdw,d2bbb,d2cart,d2cart_bbb,d2matr,d2nfr,& 1439 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,& 1440 & gprimd,dtset%mband,mpert,natom,ntypat,outd2,pawbec,pawpiezo,piezofrnl,dtset%prtbbb,& 1441 & rfasr,rfpert,rprimd,dtset%typat,ucvol,usevdw,psps%ziontypat) 1442 1443 dscrpt=' Note : temporary (transfer) database ' 1444 1445 ! Initialize the header of the DDB file 1446 call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,& 1447 & 1,xred=xred,occ=occ,ngfft=ngfft) 1448 1449 ! Open the formatted derivative database file, and write the header 1450 call status(0,dtfil%filstat,iexit,level,'call ddb_hdr_open_write') 1451 call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_ddb, dtfil%unddb) 1452 1453 call ddb_hdr_free(ddb_hdr) 1454 1455 ! Output of the dynamical matrix (master only) 1456 call status(0,dtfil%filstat,iexit,level,'call dfpt_dyout ') 1457 call dfpt_dyout(becfrnl,dtset%berryopt,blkflg,carflg,dtfil%unddb,ddkfil,dyew,dyfrlo,& 1458 & dyfrnl,dyfrx1,dyfrx2,dyfr_cplex,dyfr_nondiag,dyvdw,d2cart,d2cart_bbb,d2eig0,& 1459 & d2k0,d2lo,d2loc0,d2matr,d2nl,d2nl0,d2nl1,d2ovl,d2vn,& 1460 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,& 1461 & has_full_piezo,has_allddk,ab_out,dtset%mband,mpert,natom,ntypat,& 1462 & outd2,pawbec,pawpiezo,piezofrnl,dtset%prtbbb,dtset%prtvol,qphon,qzero,& 1463 & dtset%typat,rfdir,rfpert,rfphon,rfstrs,psps%usepaw,usevdw,psps%ziontypat) 1464 1465 close(dtfil%unddb) 1466 1467 #ifdef HAVE_NETCDF 1468 ! Output dynamical matrix in NetCDF format. 1469 call crystal_init(dtset%amu_orig(:,1), Crystal, & 1470 & dtset%spgroup, dtset%natom, dtset%npsp, psps%ntypat, & 1471 & dtset%nsym, rprimd, dtset%typat, xred, dtset%ziontypat, dtset%znucl, 1, & 1472 & dtset%nspden==2.and.dtset%nsppol==1, .false., hdr%title, & 1473 & dtset%symrel, dtset%tnons, dtset%symafm) 1474 1475 filename = strcat(dtfil%filnam_ds(4),"_DDB.nc") 1476 1477 call outddbnc(filename, mpert, d2matr, blkflg, dtset%qptn, Crystal) 1478 1479 call crystal_free(Crystal) 1480 1481 #endif 1482 1483 1484 ! In case of phonons, diagonalize the dynamical matrix 1485 if(rfphon==1)then 1486 1487 ! First, suppress the 'wings' elements, 1488 ! for which the diagonal element is not known 1489 call wings3(carflg,d2cart,mpert) 1490 1491 ! Check the analyticity of the dynamical matrix 1492 analyt=0 1493 if (rfpert(natom+2)==0 .or. rfpert(natom+2)==2 .or. sumg0==1 ) analyt=1 1494 1495 ! Diagonalize the analytic part 1496 ABI_ALLOCATE(displ,(2*3*natom*3*natom)) 1497 ABI_ALLOCATE(eigval,(3*natom)) 1498 ABI_ALLOCATE(eigvec,(2*3*natom*3*natom)) 1499 ABI_ALLOCATE(phfrq,(3*natom)) 1500 qphnrm=one 1501 call dfpt_phfrq(dtset%amu_orig(:,1),displ,d2cart,eigval,eigvec,indsym,mpert,& 1502 & dtset%nsym,natom,dtset%nsym,ntypat,phfrq,qphnrm,qphon,& 1503 & dtset%rprimd_orig(1:3,1:3,1),0,dtset%symrel,dtset%symafm,dtset%typat,ucvol) 1504 1505 ! Print the phonon frequencies 1506 call dfpt_prtph(displ,0,dtset%enunit,ab_out,natom,phfrq,qphnrm,qphon) 1507 1508 ! Check the completeness of the dynamical matrix and eventually send a warning 1509 call chkph3(carflg,0,mpert,natom) 1510 end if ! end case of phonons 1511 end if !end me == 0 1512 1513 !Compute the other terms for AHC dynamic and AHC full 1514 if (.not.(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0.or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 .and. rfuser==0 & 1515 & .and. rfmagn==0)) then 1516 if(rfphon==1) then ! AHC can only be computed in case of phonons 1517 1518 ! Stuff for parallelism 1519 if(master /= me) then 1520 ABI_ALLOCATE(phfrq,(3*natom)) 1521 ABI_ALLOCATE(displ,(2*3*natom*3*natom)) 1522 end if 1523 call xmpi_bcast (phfrq,master,mpi_enreg%comm_cell,ierr) !Broadcast phfrq and displ 1524 call xmpi_bcast (displ,master,mpi_enreg%comm_cell,ierr) !to all processes 1525 1526 if(dtset%ieig2rf == 3 .or. dtset%ieig2rf == 4 .or. dtset%ieig2rf == 5 ) then 1527 bdeigrf = dtset%bdeigrf 1528 if(dtset%bdeigrf == -1) bdeigrf = dtset%mband 1529 ! if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 1530 ! & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 1531 ! write(std_out,*)'Reading the dense grid WF file' 1532 ! call wfk_read_eigenvalues(dtfil%fnameabi_wfkfine,eigenq_fine,hdr_fine,mpi_enreg%comm_world) 1533 ! ABI_CHECK(SIZE(eigenq_fine,DIM=1)==Dtset%mband,"Size eigenq_fine != mband") 1534 ! endif 1535 if(dtset%kptopt==3 .or. dtset%kptopt==0 .or. dtset%nsym==1)then 1536 write(std_out,*) 'Entering: eig2tot' 1537 if(dtset%smdelta>0)then 1538 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 1539 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 1540 call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,& 1541 & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,& 1542 & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,& 1543 & occ_rbz_pert,hdr0,eigbrd,eigenq_fine,hdr_fine) 1544 else 1545 call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,& 1546 & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,& 1547 & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,& 1548 & occ_rbz_pert,hdr0,eigbrd) 1549 end if 1550 else 1551 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 1552 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 1553 call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,& 1554 & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,& 1555 & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,& 1556 & occ_rbz_pert,hdr0,eigbrd,eigenq_fine,hdr_fine) 1557 else 1558 call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,& 1559 & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,& 1560 & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,& 1561 & occ_rbz_pert,hdr0) 1562 end if 1563 end if 1564 write(std_out,*) 'Leaving: eig2tot' 1565 end if 1566 end if 1567 if (dtset%ieig2rf > 0) then 1568 ABI_DEALLOCATE(eigen0_pert) 1569 ABI_DEALLOCATE(eigenq_pert) 1570 ABI_DEALLOCATE(occ_rbz_pert) 1571 ABI_DEALLOCATE(eigen1_pert) 1572 call hdr_free(hdr0) 1573 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 1574 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 1575 ! call hdr_free(hdr0) 1576 call hdr_free(hdr_fine) 1577 ABI_DEALLOCATE(eigenq_fine) 1578 end if 1579 end if ! ieig2rf == 3 or %ieig2rf == 4 or %ieig2rf == 5 1580 end if ! rfphon==1 1581 end if 1582 if(dtset%efmas>0) then 1583 ABI_DEALLOCATE(eigen0_pert) 1584 ABI_DEALLOCATE(eigen1_pert) 1585 end if 1586 ABI_DEALLOCATE(doccde) 1587 1588 1589 if(me==0)then 1590 if (.not.(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0 .or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 .and.rfuser==0 & 1591 & .and. rfmagn==0) )then 1592 if(rfphon==1)then 1593 ! Compute and print the T=0 Fan, and possibly DDW contributions to the eigenenergies. 1594 if(dtset%ieig2rf > 0) then 1595 write(message, '(80a,9a)' ) ('=',mu=1,80),ch10,ch10,& 1596 & ' ---- T=0 shift of eigenenergies due to electron-phonon interation at q ---- ',ch10,& 1597 & ' Warning : the total shift must be computed through anaddb, ',ch10,& 1598 & ' here, only the contribution of one q point is printed. ',ch10,& 1599 & ' Print first the electronic eigenvalues, then the q-dependent Fan shift of eigenvalues.' 1600 call wrtout(ab_out,message,'COLL') 1601 call wrtout(std_out, message,'COLL') 1602 1603 if(qeq0)then 1604 write(message, '(a)' )& 1605 & ' Phonons at gamma, also compute the Diagonal Debye-Waller shift of eigenvalues.' 1606 call wrtout(ab_out,message,'COLL') 1607 call wrtout(std_out,message,'COLL') 1608 end if 1609 1610 write(message, '(a)' ) ' ' 1611 call wrtout(ab_out,message,'COLL') 1612 call wrtout(std_out,message,'COLL') 1613 1614 call prteigrs(eigen0,dtset%enunit,fermie,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,& 1615 & dtset%mband,dtset%nband,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,3,0,dtset%prtvol,& 1616 & eigen0,zero,zero,dtset%wtk) 1617 1618 write(message, '(a)' ) ch10 1619 call wrtout(ab_out,message,'COLL') 1620 call wrtout(std_out,message,'COLL') 1621 1622 ! Compute and print Fan contribution 1623 ABI_ALLOCATE(eigen_fan,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1624 ABI_ALLOCATE(eigen_fan_mean,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1625 call elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_fan,gprimd,& 1626 & dtset%mband,natom,dtset%nkpt,dtset%nsppol,1,phfrq,dtset%prtvol) 1627 call eigen_meandege(eigen0,eigen_fan,eigen_fan_mean,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,2) 1628 call prteigrs(eigen_fan_mean,dtset%enunit,fermie,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,& 1629 & dtset%mband,dtset%nband,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,5,0,dtset%prtvol,& 1630 & eigen0,zero,zero,dtset%wtk) 1631 1632 if(qeq0 .or. dtset%getgam_eig2nkq>0)then 1633 1634 write(message, '(a)' ) ch10 1635 call wrtout(ab_out,message,'COLL') 1636 call wrtout(std_out,message,'COLL') 1637 1638 ! Compute and print Diagonal Debye-Waller contribution 1639 ABI_ALLOCATE(eigen_ddw,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1640 ABI_ALLOCATE(eigen_ddw_mean,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1641 if(qeq0)then 1642 call elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_ddw,gprimd,& 1643 & dtset%mband,natom,dtset%nkpt,dtset%nsppol,2,phfrq,dtset%prtvol) 1644 if(results_respfn%gam_jdtset == -dtset%jdtset)then 1645 sz1=dtset%mband*dtset%nsppol 1646 sz2=natom*dim_eig2nkq 1647 ABI_ALLOCATE(results_respfn%gam_eig2nkq,(2,sz1,dtset%nkpt,3,natom,3,sz2)) 1648 results_respfn%gam_eig2nkq(:,:,:,:,:,:,:)=eig2nkq(:,:,:,:,:,:,:) 1649 results_respfn%gam_jdtset=dtset%jdtset 1650 end if 1651 else if(dtset%getgam_eig2nkq>0)then 1652 if(results_respfn%gam_jdtset==dtset%getgam_eig2nkq)then 1653 call elph2_fanddw(dim_eig2nkq,displ,results_respfn%gam_eig2nkq,eigen_ddw,& 1654 & gprimd,dtset%mband,natom,dtset%nkpt,dtset%nsppol,2,phfrq,dtset%prtvol) 1655 else 1656 write(message,'(a,i0,2a,i0,2a)')& 1657 & 'results_respfn%gam_jdtset=',results_respfn%gam_jdtset,ch10,& 1658 & 'dtset%getgam_eig2nkq=',dtset%getgam_eig2nkq,ch10,& 1659 & 'So, it seems eig2nkq at gamma has not yet been computed, while it is needed now.' 1660 MSG_BUG(message) 1661 end if 1662 end if 1663 call eigen_meandege(eigen0,eigen_ddw,eigen_ddw_mean,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,2) 1664 call prteigrs(eigen_ddw_mean,dtset%enunit,fermie,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,& 1665 & dtset%mband,dtset%nband,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,6,0,dtset%prtvol,& 1666 & eigen0,zero,zero,dtset%wtk) 1667 1668 write(message, '(a)' ) ch10 1669 call wrtout(ab_out,message,'COLL') 1670 call wrtout(std_out,message,'COLL') 1671 1672 ! Print sum of mean Fan and DDW 1673 ABI_ALLOCATE(eigen_fanddw,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1674 eigen_fanddw=eigen_fan_mean+eigen_ddw_mean 1675 call prteigrs(eigen_fanddw,dtset%enunit,fermie,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,& 1676 & dtset%mband,dtset%nband,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,7,0,dtset%prtvol,& 1677 & eigen0,zero,zero,dtset%wtk) 1678 1679 ABI_DEALLOCATE(eigen_ddw) 1680 ABI_DEALLOCATE(eigen_ddw_mean) 1681 ABI_DEALLOCATE(eigen_fanddw) 1682 1683 end if 1684 1685 ABI_DEALLOCATE(eigen_fan) 1686 ABI_DEALLOCATE(eigen_fan_mean) 1687 end if 1688 1689 ! In case of a non-analytical part, 1690 ! get the phonon frequencies for three different directions (in cartesian coordinates) 1691 if(analyt==0)then 1692 qphnrm=zero 1693 do idir=1,3 1694 ! Need to know the corresponding dielectric constant 1695 if(carflg(idir,natom+2,idir,natom+2)==1)then 1696 qphon(:)=zero ; qphon(idir)=one 1697 ! Get the phonon frequencies 1698 call dfpt_phfrq(dtset%amu_orig(:,1),displ,d2cart,eigval,eigvec,indsym,mpert,& 1699 & dtset%nsym,natom,dtset%nsym,ntypat,phfrq,qphnrm,qphon,& 1700 & dtset%rprimd_orig(1:3,1:3,1),0,dtset%symrel,dtset%symafm,dtset%typat,ucvol) 1701 ! Print the phonon frequencies 1702 call dfpt_prtph(displ,0,dtset%enunit,ab_out,natom,phfrq,qphnrm,qphon) 1703 ! Check the completeness of the dynamical matrix 1704 ! and eventually send a warning 1705 call chkph3(carflg,idir,mpert,natom) 1706 end if 1707 end do 1708 if (idir < 4) then 1709 qphon(idir)=zero 1710 end if 1711 end if 1712 1713 ABI_DEALLOCATE(displ) 1714 ABI_DEALLOCATE(eigval) 1715 ABI_DEALLOCATE(eigvec) 1716 ABI_DEALLOCATE(phfrq) 1717 end if ! rfphon == 1 1718 ABI_DEALLOCATE(carflg) 1719 ABI_DEALLOCATE(d2cart) 1720 ABI_DEALLOCATE(d2matr) 1721 end if ! End condition on if.not. 1722 end if ! master node 1723 1724 !Deallocate arrays 1725 if (allocated(displ)) then 1726 ABI_DEALLOCATE(displ) 1727 end if 1728 if (allocated(eigval)) then 1729 ABI_DEALLOCATE(eigval) 1730 end if 1731 if (allocated(eigvec)) then 1732 ABI_DEALLOCATE(eigvec) 1733 end if 1734 if (allocated(phfrq)) then 1735 ABI_DEALLOCATE(phfrq) 1736 end if 1737 1738 ABI_DEALLOCATE(clflg) 1739 ABI_DEALLOCATE(amass) 1740 ABI_DEALLOCATE(atindx) 1741 ABI_DEALLOCATE(atindx1) 1742 ABI_DEALLOCATE(blkflg) 1743 ABI_DEALLOCATE(dyew) 1744 ABI_DEALLOCATE(dyewq0) 1745 ABI_DEALLOCATE(dyfrlo) 1746 ABI_DEALLOCATE(dyfrnl) 1747 ABI_DEALLOCATE(dyfrwf) 1748 ABI_DEALLOCATE(dyfrx1) 1749 ABI_DEALLOCATE(dyfrx2) 1750 ABI_DEALLOCATE(dyvdw) 1751 ABI_DEALLOCATE(d2bbb) 1752 ABI_DEALLOCATE(d2cart_bbb) 1753 ABI_DEALLOCATE(d2eig0) 1754 ABI_DEALLOCATE(d2k0) 1755 ABI_DEALLOCATE(d2lo) 1756 ABI_DEALLOCATE(d2loc0) 1757 ABI_DEALLOCATE(d2nfr) 1758 ABI_DEALLOCATE(d2nl) 1759 ABI_DEALLOCATE(d2nl0) 1760 ABI_DEALLOCATE(d2nl1) 1761 ABI_DEALLOCATE(d2ovl) 1762 ABI_DEALLOCATE(d2vn) 1763 ABI_DEALLOCATE(eigen0) 1764 ABI_DEALLOCATE(eig2nkq) 1765 ABI_DEALLOCATE(eigbrd) 1766 ABI_DEALLOCATE(eltcore) 1767 ABI_DEALLOCATE(elteew) 1768 ABI_DEALLOCATE(eltfrhar) 1769 ABI_DEALLOCATE(eltfrnl) 1770 ABI_DEALLOCATE(eltfrloc) 1771 ABI_DEALLOCATE(eltfrkin) 1772 ABI_DEALLOCATE(eltfrxc) 1773 ABI_DEALLOCATE(eltvdw) 1774 ABI_DEALLOCATE(becfrnl) 1775 ABI_DEALLOCATE(piezofrnl) 1776 call efmasdeg_free_array(efmasdeg) 1777 call efmasfr_free_array(efmasfr) 1778 ABI_DEALLOCATE(grxc) 1779 ABI_DEALLOCATE(indsym) 1780 ABI_DEALLOCATE(kxc) 1781 ABI_DEALLOCATE(nattyp) 1782 ABI_DEALLOCATE(pertsy) 1783 ABI_DEALLOCATE(rfpert) 1784 ABI_DEALLOCATE(rhog) 1785 ABI_DEALLOCATE(rhor) 1786 ABI_DEALLOCATE(symq) 1787 ABI_DEALLOCATE(symrec) 1788 ABI_DEALLOCATE(vtrial) 1789 ABI_DEALLOCATE(ylm) 1790 ABI_DEALLOCATE(ylmgr) 1791 call pawfgr_destroy(pawfgr) 1792 if (psps%usepaw==1) then 1793 call pawrhoij_free(pawrhoij) 1794 call paw_an_free(paw_an) 1795 call paw_ij_free(paw_ij) 1796 call pawfgrtab_free(pawfgrtab) 1797 end if 1798 ABI_DEALLOCATE(nhat) 1799 ABI_DEALLOCATE(nhatgr) 1800 ABI_DATATYPE_DEALLOCATE(pawrhoij) 1801 ABI_DATATYPE_DEALLOCATE(paw_an) 1802 ABI_DATATYPE_DEALLOCATE(paw_ij) 1803 ABI_DATATYPE_DEALLOCATE(pawfgrtab) 1804 if(rfphon==1.and.psps%n1xccc/=0)then 1805 ABI_DEALLOCATE(blkflgfrx1) 1806 end if 1807 1808 !Clean the header 1809 call hdr_free(hdr) 1810 1811 !Clean GPU data 1812 #if defined HAVE_GPU_CUDA 1813 if (dtset%use_gpu_cuda==1) then 1814 call dealloc_hamilt_gpu(0,dtset%use_gpu_cuda) 1815 end if 1816 #endif 1817 1818 call status(0,dtfil%filstat,iexit,level,'exit') 1819 1820 call timab(138,2,tsec) 1821 call timab(132,2,tsec) 1822 1823 DBG_EXIT("COLL") 1824 1825 end subroutine respfn