TABLE OF CONTENTS
ABINIT/d2frnl [ Functions ]
NAME
d2frnl
FUNCTION
Compute the frozen-wavefunction non-local contribution for response functions (strain and/or phonon)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GM, AR, MB, MT, AM) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnuC.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol)=<G|Cnk>=Fourier coefficients of WF dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset 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 eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere has_allddk= True if all ddk file are present on disk indsym(4,nsym,natom)=indirect indexing array for atom labels kg(3,mpw*mkmem)=work array for coordinates of G vectors in basis primitive translations mgfftf=maximum size of 1D FFTs for the fine FFT grid (PAW) mpi_enreg=information about MPI parallelization mpsang= 1+maximum angular momentum for nonlocal pseudopotentials my_natom=number of atoms treated by current processor natom=number of atoms in unit cell nfftf= -PAW ONLY- number of FFT grid points for the fine grid (nfftf=nfft for norm-conserving potential runs) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid (ngs_rbzfftf=ngfft for norm-conserving potential runs) npwarr(nkpt)=number of planewaves at each k point, and boundary ntypat=integer specification of atom type (1, 2, ...) occ(mband*nkpt*nsppol)=occupation numbers of bands (usually 2) at each k point rfphon=1 if non local contribution of dynamical matrix have to be computed rfstrs!=0 if non local contribution of elastic tensor have to be computed paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawang <type(pawang_type)>=paw angular mesh and related data pawbec= flag for the computation of Born Effective Charge within PAW ; set to 1 if yes pawpiezo= flag for the computation of piezoelectric tensor within PAW ; set to 1 if yes pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=phase information related to structure factor ph1df(2,3*(2*mgfftf+1)*natom)=phase information related to structure factor on the fine FFT grid (PAW) psps <type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3)=dimensional real space primitive translations (bohr) symrec(3,3,nsym)=symmetries in reciprocal space (dimensionless) vtrial(nfftf,nspden)=total potential (Hartree+XC+loc) vxc(nfftf,nspden)=XC potential xred(3,natom)=reduced coordinates of atoms (dimensionless) ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylmgr(mpw*mkmem,9,mpsang*mpsang*useylm)= gradients of real spherical harmonics for each G and k point
OUTPUT
becfrnl(3,natom,3*pawbec)=NL frozen contribution to Born Effective Charges (PAW only) (3,natom) = derivative wr to the displ. of one atom in one direction (3) = derivative wr to electric field in one direction piezofrnl(3,6*pawpiezo)=NL frozen contribution to piezoelectric tensor (PAW only) dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)= non-symmetrized non-local contribution to the dynamical matrix If NCPP, it depends on one atom If PAW, it depends on two atoms eltfrnl(6+3*natom,6)=non-symmetrized non-local contribution to the elastic tensor
SIDE EFFECTS
===== if psps%usepaw==1 pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid pawfgrtab(:)%gylmgr2 are deallocated here pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data (gradients of rhoij for each atom with respect to atomic positions are computed here)
PARENTS
respfn
CHILDREN
appdig,check_degeneracies,destroy_hamiltonian,dotprod_g init_hamiltonian,load_k_hamiltonian,load_spin_hamiltonian,metric,mkffnl mkkin,mkkpg,nonlop,paw_ij_free,paw_ij_init,paw_ij_nullify paw_ij_reset_flags,pawaccrhoij,pawcprj_alloc,pawcprj_free,pawdij2e1kb pawdijfr,pawfgrtab_free,pawfgrtab_init,pawgrnl,pawrhoij_free pawrhoij_gather,pawrhoij_nullify,pawtab_get_lsize,strconv,symrhoij timab,wfk_close,wfk_open_read,wfk_read_bks,wrtout,xmpi_sum
SOURCE
95 #if defined HAVE_CONFIG_H 96 #include "config.h" 97 #endif 98 99 #include "abi_common.h" 100 101 subroutine d2frnl(becfrnl,cg,dtfil,dtset,dyfrnl,dyfr_cplex,dyfr_nondiag,efmasdeg,efmasfr,eigen,eltfrnl,& 102 & gsqcut,has_allddk,indsym,kg,mgfftf,mpi_enreg,mpsang,my_natom,natom,nfftf,ngfft,ngfftf,npwarr,& 103 & occ,paw_ij,pawang,pawbec,pawfgrtab,pawpiezo,pawrad,pawrhoij,pawtab,ph1d,ph1df,piezofrnl,psps,& 104 & rprimd,rfphon,rfstrs,symrec,vtrial,vxc,xred,ylm,ylmgr) 105 106 use defs_basis 107 use defs_datatypes 108 use defs_abitypes 109 use m_profiling_abi 110 use m_xmpi 111 use m_errors 112 use m_cgtools 113 use m_nctk 114 use m_hamiltonian 115 use m_efmas_defs 116 use m_wfk 117 118 use m_geometry, only : metric 119 use m_efmas, only : check_degeneracies 120 use m_io_tools, only : file_exists 121 use m_hdr, only : hdr_skip 122 use m_pawang, only : pawang_type 123 use m_pawrad, only : pawrad_type 124 use m_pawtab, only : pawtab_type,pawtab_get_lsize 125 use m_pawfgrtab,only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free 126 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify,& 127 & paw_ij_reset_flags 128 use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_free, pawrhoij_gather,& 129 & pawrhoij_nullify, symrhoij 130 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get,& 131 & pawcprj_copy, pawcprj_free 132 use m_pawdij, only : pawdijfr 133 use m_kg, only : mkkin 134 135 !This section has been created automatically by the script Abilint (TD). 136 !Do not modify the following lines by hand. 137 #undef ABI_FUNC 138 #define ABI_FUNC 'd2frnl' 139 use interfaces_14_hidewrite 140 use interfaces_18_timing 141 use interfaces_32_util 142 use interfaces_41_geometry 143 use interfaces_65_paw 144 use interfaces_66_nonlocal 145 !End of the abilint section 146 147 implicit none 148 149 !Arguments ------------------------------------ 150 !scalars 151 integer,intent(in) :: dyfr_cplex,dyfr_nondiag,mgfftf,mpsang,my_natom,natom 152 integer,intent(in) :: nfftf,pawbec,pawpiezo,rfphon,rfstrs 153 real(dp),intent(in) :: gsqcut 154 type(MPI_type),intent(in) :: mpi_enreg 155 type(datafiles_type),intent(in) :: dtfil 156 type(dataset_type),intent(in) :: dtset 157 type(pawang_type),intent(in) :: pawang 158 type(pseudopotential_type),intent(in) :: psps 159 !arrays 160 integer,intent(in) :: indsym(4,dtset%nsym,natom),kg(3,dtset%mpw*dtset%mkmem) 161 integer,intent(in) :: ngfft(18),ngfftf(18),npwarr(dtset%nkpt) 162 integer,intent(in) :: symrec(3,3,dtset%nsym) 163 real(dp),intent(in) :: cg(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol) 164 real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 165 real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 166 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*natom) 167 real(dp),intent(in) :: ph1df(2,3*(2*mgfftf+1)*natom),rprimd(3,3) 168 real(dp),intent(in) :: vxc(nfftf,dtset%nspden),xred(3,natom) 169 real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,mpsang*mpsang*psps%useylm) 170 real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,9,mpsang*mpsang*psps%useylm) 171 real(dp),intent(in),target :: vtrial(nfftf,dtset%nspden) 172 real(dp),intent(out) :: becfrnl(3,natom,3*pawbec),piezofrnl(6,3*pawpiezo) 173 real(dp),intent(out) :: dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag) 174 real(dp),intent(out) :: eltfrnl(6+3*natom,6) 175 logical,intent(inout):: has_allddk 176 type(efmasdeg_type),allocatable,intent(out):: efmasdeg(:) 177 type(efmasfr_type),allocatable,intent(out):: efmasfr(:,:) 178 type(paw_ij_type),intent(in) :: paw_ij(my_natom) 179 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 180 type(pawrad_type),intent(in) :: pawrad(psps%ntypat) 181 type(pawrhoij_type),intent(inout),target :: pawrhoij(my_natom*psps%usepaw) 182 type(pawtab_type),intent(in) :: pawtab(psps%ntypat) 183 184 !Local variables------------------------------- 185 !scalars 186 integer,parameter :: formeig1=1,usecprj=0 187 integer :: bdtot_index,bufdim,choice_bec2,choice_bec54,choice_efmas,choice_phon,choice_strs,choice_piez3,choice_piez55 188 integer :: cplex,cplx,cpopt,cpopt_bec,ddkcase 189 integer :: dimffnl,dimffnl_str,dimnhat,ia,iatom,iashift,iband,jband,ibg,icg,icplx,ideg,ider,idir 190 integer :: ider_str,idir_ffnl,idir_str,ielt,ieltx,ierr,ii,ikg,ikpt,ilm,ipw 191 integer :: ispinor,isppol,istwf_k,isub,itypat,jj,jsub,klmn,master,me,mu 192 integer :: my_comm_atom,n1,n2,n3,nband_k,ncpgr,nfftot,ngrhoij,nkpg,nnlout_bec1,nnlout_bec2,nnlout_efmas 193 integer :: nnlout_piez1,nnlout_piez2,nnlout_phon,nnlout_strs,npw_,npw_k,nsp,nsploop,nu 194 integer :: optgr,optgr2,option,option_rhoij,optstr,optstr2,paw_opt,paw_opt_1,paw_opt_3,paw_opt_efmas 195 integer :: shift_rhoij,signs,signs_field,spaceworld,sz2,sz3,tim_nonlop 196 real(dp) :: arg,eig_k,enl,enlk,occ_k,ucvol,wtk_k 197 logical :: has_ddk_file,need_becfr,need_efmas,need_piezofr,paral_atom,t_test,usetimerev 198 character(len=500) :: msg 199 type(gs_hamiltonian_type) :: gs_ham 200 !arrays 201 integer :: ik_ddk(3),ddkfil(3) 202 integer,allocatable :: dimlmn(:),kg_k(:,:),l_size_atm(:) 203 integer,pointer :: my_atmtab(:) 204 real(dp) :: dotprod(2),dummy(0),gmet(3,3),gprimd(3,3),grhoij(3),kpoint(3),nonlop_dum(1,1) 205 real(dp) :: rmet(3,3),tsec(2) 206 real(dp),allocatable :: becfrnl_tmp(:,:,:),becfrnlk(:,:,:),becij(:,:,:,:),cg_left(:,:) 207 real(dp),allocatable :: cwavef(:,:),ddk(:,:),ddkinpw(:,:,:),dyfrnlk(:,:) 208 real(dp),allocatable :: elt_work(:,:),eltfrnlk(:,:),enlout_bec1(:),enlout_bec2(:),enlout_efmas(:) 209 real(dp),allocatable :: enlout_piez1(:),enlout_piez2(:),enlout_phon(:),enlout_strs(:) 210 real(dp),allocatable :: gh2c(:,:),gs2c(:,:) 211 real(dp),allocatable :: kpg_k(:,:),mpibuf(:),nhat_dum(:,:),piezofrnlk(:,:),ph3d(:,:,:) 212 real(dp),allocatable :: svectout(:,:),ylm_k(:,:),ylmgr_k(:,:,:) 213 real(dp),allocatable,target :: ffnl(:,:,:,:),ffnl_str(:,:,:,:,:) 214 character(len=fnlen) :: fiwfddk(3) 215 type(paw_ij_type),allocatable :: paw_ij_tmp(:) 216 type(pawcprj_type),allocatable,target :: cwaveprj(:,:) 217 type(pawfgrtab_type),allocatable :: pawfgrtab_tmp(:) 218 type(pawrhoij_type),pointer :: pawrhoij_tot(:) 219 type(wfk_t) :: ddkfiles(3) 220 221 ! ************************************************************************* 222 223 DBG_ENTER("COLL") 224 225 call timab(159,1,tsec) 226 227 write(msg,'(3a)')ch10,' ==> Calculation of the frozen part of the second order derivatives, this can take some time...',ch10 228 call wrtout(std_out,msg,'COLL') 229 230 !Set up parallelism 231 spaceworld=mpi_enreg%comm_cell 232 me=mpi_enreg%me_kpt 233 master=0 234 paral_atom=(my_natom/=natom) 235 my_comm_atom=mpi_enreg%comm_atom 236 my_atmtab=>mpi_enreg%my_atmtab 237 238 !Compute gmet, gprimd and ucvol from rprimd 239 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 240 241 !If needed, check for ddk files (used for effective charges) 242 if (pawbec==1.or.pawpiezo==1) then 243 ddkfil(:)=0 244 do ii=1,3 245 ddkcase=ii+natom*3 246 call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk(ii)) 247 t_test = file_exists(fiwfddk(ii)) 248 ! Trick needed to run Abinit test suite in netcdf mode. 249 if (.not. t_test .and. file_exists(nctk_ncify(fiwfddk(ii)))) then 250 t_test = .True.; fiwfddk(ii) = nctk_ncify(fiwfddk(ii)) 251 write(msg,"(3a)")"- File: ",trim(fiwfddk(ii))," does not exist but found netcdf file with similar name." 252 call wrtout(std_out,msg,'COLL') 253 end if 254 if (t_test) ddkfil(ii)=20+ii ! Note the use of unit numbers 21, 22 and 23 255 end do 256 has_ddk_file=(any(ddkfil(:)>0)) 257 has_allddk =(all(ddkfil(:)>0)) 258 else 259 has_ddk_file=.FALSE. 260 has_allddk =.FALSE. 261 end if 262 263 if(pawbec==1.or.pawpiezo==1.and.has_ddk_file) then 264 if(.not.has_allddk) then 265 write(msg,'(5a)')ch10,& 266 & ' WARNING: All ddk perturbations are needed to compute',ch10,& 267 & ' the frozen part of Born effective charges and/or piezoelectric tensor.',ch10 268 call wrtout(std_out,msg,'COLL') 269 else 270 write(msg,'(5a)')ch10,& 271 & ' All ddk perturbations are available.',ch10,& 272 & ' The frozen part of Born effective charges and/or piezoelectric tensor will be computed',ch10 273 call wrtout(std_out,msg,'COLL') 274 end if 275 end if 276 277 need_becfr=(pawbec==1.and.has_ddk_file) 278 need_piezofr=(pawpiezo==1.and.has_ddk_file) 279 280 !Initialization of frozen non local array 281 if(rfphon==1) then 282 dyfrnl(:,:,:,:,:)=zero 283 ABI_ALLOCATE(dyfrnlk,(6,natom)) 284 end if 285 if(rfstrs/=0)then 286 eltfrnl(:,:)=zero;enl=zero 287 ABI_ALLOCATE(eltfrnlk,(6+3*natom,6)) 288 end if 289 if (need_becfr) then 290 becfrnl(:,:,:)=zero 291 ABI_ALLOCATE(becfrnlk,(3,natom,3)) 292 end if 293 if (need_piezofr) then 294 piezofrnl(:,:)=zero 295 ABI_ALLOCATE(piezofrnlk,(6,3)) 296 end if 297 need_efmas=dtset%efmas>0 298 if(need_efmas.and.(rfphon==1.or.rfstrs/=0.or.need_becfr.or.need_piezofr)) then 299 write(msg,'(5a)')ch10,& 300 & ' ERROR: Efmas calculation is incompatible with phonons, elastic tensor, Born effective charges,',ch10,& 301 & ' and piezoelectric tensor calculations. Please revise your input.',ch10 302 MSG_ERROR(msg) 303 end if 304 305 !Common initialization 306 bdtot_index=0;ibg=0;icg=0 307 nsploop=dtset%nsppol;if (dtset%nspden==4) nsploop=4 308 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 309 nfftot=ngfftf(1)*ngfftf(2)*ngfftf(3) 310 311 !Common data for "nonlop" routine 312 tim_nonlop=6 313 signs=1 ; signs_field = 2 ; eig_k=zero ; idir=0 314 ! ffnl are in cartesian coordinates for EFMAS (idir==4), 315 ! in contrast to reduced coordinates for the other responses (idir==0). 316 idir_ffnl=0 ; if(need_efmas) idir_ffnl=4 317 choice_phon=0;choice_strs=0 318 if(rfphon==1)then 319 shift_rhoij=0 320 choice_phon=4 321 nnlout_phon=max(1,6*natom) 322 ABI_ALLOCATE(enlout_phon,(nnlout_phon)) 323 end if 324 if(rfstrs/=0)then 325 shift_rhoij=6 326 choice_strs=6 327 nnlout_strs=6*(3*natom+6) 328 ABI_ALLOCATE(enlout_strs,(nnlout_strs)) 329 end if 330 if (psps%usepaw==0) then 331 paw_opt=0 ; cpopt=-1 332 else 333 paw_opt=2 ; cpopt=1+2*usecprj 334 end if 335 if(need_piezofr)then 336 choice_piez3 = 3 337 choice_piez55 = 55 338 nnlout_piez1 = 6 339 nnlout_piez2 = 36 340 paw_opt_1 = 1 341 paw_opt_3 = 3 342 ABI_ALLOCATE(enlout_piez1,(nnlout_piez1)) 343 ABI_ALLOCATE(enlout_piez2,(nnlout_piez2)) 344 end if 345 if (need_becfr) then 346 choice_bec2=2 ; choice_bec54=54 347 nnlout_bec1=max(1,3*natom) ; nnlout_bec2=max(1,18*natom); 348 paw_opt_1=1 ; paw_opt_3=3 ; cpopt_bec=-1 349 ABI_ALLOCATE(enlout_bec1,(nnlout_bec1)) 350 ABI_ALLOCATE(enlout_bec2,(nnlout_bec2)) 351 else 352 choice_bec2=0 ; choice_bec54=0 353 nnlout_bec1=0 354 end if 355 if(need_efmas) then 356 ABI_MALLOC(enlout_efmas,(0)) 357 ABI_DATATYPE_ALLOCATE(efmasdeg,(dtset%nkpt)) 358 ABI_DATATYPE_ALLOCATE(efmasfr,(dtset%nkpt,dtset%mband)) 359 end if 360 361 !Initialize Hamiltonian (k-independent terms) 362 call init_hamiltonian(gs_ham,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,natom,& 363 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,& 364 & paw_ij=paw_ij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,& 365 & usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom) 366 367 !===== PAW specific section 368 if (psps%usepaw==1) then 369 370 ! Define several sizes & flags 371 ncpgr=0;ngrhoij=0 372 if(rfphon==1)then 373 ncpgr=3;ngrhoij=3 374 end if 375 if(rfphon==1.or.need_becfr)then 376 ncpgr=6;ngrhoij=6 377 end if 378 if(rfstrs/=0.and.rfphon==1)then 379 ncpgr=9;ngrhoij=9 380 end if 381 if(rfstrs/=0.or.need_piezofr)then 382 ncpgr=9;ngrhoij=9 383 end if 384 385 ! If PAW and Born Eff. Charges, one has to compute some additional data: 386 ! For each atom and for electric field direction k: 387 ! becij(k)=<Phi_i|r_k-R_k|Phi_j>-<tPhi_i|r_k-R_k|tPhi_j> + sij.R_k 388 if (need_becfr.or.need_piezofr) then 389 ABI_ALLOCATE(becij,(gs_ham%dimekb1,gs_ham%dimekb2,dtset%nspinor**2,3)) 390 becij=zero 391 ABI_DATATYPE_ALLOCATE(paw_ij_tmp,(my_natom)) 392 ABI_DATATYPE_ALLOCATE(pawfgrtab_tmp,(my_natom)) 393 call paw_ij_nullify(paw_ij_tmp) 394 cplex=1;nsp=1 ! Force nsppol/nspden to 1 because Dij^(1) due to electric field is spin-independent 395 call paw_ij_init(paw_ij_tmp,cplex,dtset%nspinor,nsp,nsp,dtset%pawspnorb,natom,psps%ntypat,& 396 & dtset%typat,pawtab,has_dijfr=1,comm_atom=my_comm_atom,mpi_atmtab=my_atmtab ) 397 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,mpi_atmtab=my_atmtab) 398 call pawfgrtab_init(pawfgrtab_tmp,1,l_size_atm,dtset%nspden,dtset%typat,& 399 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 400 ABI_DEALLOCATE(l_size_atm) 401 do ii=1,3 ! Loop over direction of electric field 402 call paw_ij_reset_flags(paw_ij_tmp,all=.True.) 403 call pawdijfr(cplex,gprimd,ii,natom+2,my_natom,natom,nfftf,ngfftf,nsp,psps%ntypat,& 404 & 0,paw_ij_tmp,pawang,pawfgrtab_tmp,pawrad,pawtab,& 405 & (/zero,zero,zero/),rprimd,ucvol,vtrial,vtrial,vxc,xred,& 406 & comm_atom=my_comm_atom, mpi_atmtab=my_atmtab ) ! vtrial not used here 407 do isppol=1,dtset%nspinor**2 408 call pawdij2e1kb(paw_ij_tmp(:),nsp,my_comm_atom,e1kbfr=becij(:,:,:,ii),mpi_atmtab=my_atmtab) 409 end do 410 end do 411 call paw_ij_free(paw_ij_tmp) 412 call pawfgrtab_free(pawfgrtab_tmp) 413 ABI_DATATYPE_DEALLOCATE(paw_ij_tmp) 414 ABI_DATATYPE_DEALLOCATE(pawfgrtab_tmp) 415 end if 416 417 ! PAW occupancies: need to communicate when paral atom is activated 418 if (paral_atom) then 419 ABI_DATATYPE_ALLOCATE(pawrhoij_tot,(natom)) 420 call pawrhoij_nullify(pawrhoij_tot) 421 call pawrhoij_gather(pawrhoij,pawrhoij_tot,-1,my_comm_atom) 422 else 423 pawrhoij_tot => pawrhoij 424 end if 425 426 ! Projected WF (cprj) and PAW occupancies (& gradients) 427 ABI_DATATYPE_ALLOCATE(cwaveprj,(natom,dtset%nspinor)) 428 call pawcprj_alloc(cwaveprj,ncpgr,gs_ham%dimcprj) 429 do iatom=1,natom 430 sz2=pawrhoij_tot(iatom)%cplex*pawrhoij_tot(iatom)%lmn2_size 431 sz3=pawrhoij_tot(iatom)%nspden 432 ABI_ALLOCATE(pawrhoij_tot(iatom)%grhoij,(ngrhoij,sz2,sz3)) 433 pawrhoij_tot(iatom)%ngrhoij=ngrhoij 434 pawrhoij_tot(iatom)%grhoij=zero 435 end do 436 usetimerev=(dtset%kptopt>0.and.dtset%kptopt<3) 437 438 else 439 ABI_DATATYPE_ALLOCATE(cwaveprj,(0,0)) 440 end if !PAW 441 442 !If needed, manage ddk files 443 !Open ddk WF file(s) 444 if (need_becfr.or.need_piezofr) then 445 do ii=1,3 ! Loop over elect. field directions 446 if (ddkfil(ii)/=0) then 447 write(msg, '(a,a)') '-open ddk wf file :',trim(fiwfddk(ii)) 448 call wrtout(std_out,msg,'COLL') 449 call wrtout(ab_out,msg,'COLL') 450 call wfk_open_read(ddkfiles(ii),fiwfddk(ii),formeig1,dtset%iomode,ddkfil(ii),spaceworld) 451 end if 452 end do 453 end if 454 455 !LOOP OVER SPINS 456 do isppol=1,dtset%nsppol 457 458 ! Continue to initialize the Hamiltonian (PAW DIJ coefficients) 459 call load_spin_hamiltonian(gs_ham,isppol,with_nonlocal=.true.) 460 461 ! Rewind (k+G) data if needed 462 ikg=0 463 464 ! Loop over k points 465 do ikpt=1,dtset%nkpt 466 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 467 istwf_k=dtset%istwfk(ikpt) 468 npw_k=npwarr(ikpt) 469 wtk_k=dtset%wtk(ikpt) 470 kpoint(:)=dtset%kptns(:,ikpt) 471 472 ! Skip this k-point if not the proper processor 473 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 474 bdtot_index=bdtot_index+nband_k 475 cycle 476 end if 477 478 ! If needed, manage ddk files 479 if (need_becfr.or.need_piezofr) then 480 do ii=1,3 ! Loop over elect. field directions 481 if (ddkfil(ii)/=0)then 482 ! Number of k points to skip in the full set of k pointsp 483 ik_ddk(ii) = wfk_findk(ddkfiles(ii), kpoint) 484 ABI_CHECK(ik_ddk(ii) /= -1, "Cannot find k-point in DDK") 485 npw_ = ddkfiles(ii)%hdr%npwarr(ik_ddk(ii)) 486 if (npw_/=npw_k) then 487 write(unit=msg,fmt='(a,i3,a,i5,a,i3,a,a,i5,a,a,i5)')& 488 & 'For isppol = ',isppol,', ikpt = ',ikpt,' and idir = ',ii,ch10,& 489 & 'the number of plane waves in the ddk file is equal to', npw_,ch10,& 490 & 'while it should be ',npw_k 491 MSG_ERROR(msg) 492 end if 493 494 end if 495 end do 496 end if 497 498 ABI_ALLOCATE(cwavef,(2,npw_k*dtset%nspinor)) 499 if (need_becfr.or.need_piezofr) then 500 ABI_ALLOCATE(svectout,(2,npw_k*dtset%nspinor)) 501 end if 502 if (need_efmas) then 503 ABI_MALLOC(cg_left,(2,npw_k*dtset%nspinor)) 504 ABI_MALLOC(gh2c,(2,npw_k*dtset%nspinor)) 505 ABI_MALLOC(gs2c,(2,npw_k*dtset%nspinor)) 506 end if 507 508 ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang*psps%useylm)) 509 if(rfstrs/=0.or.need_becfr.or.need_piezofr.or.need_efmas)then 510 ABI_ALLOCATE(ylmgr_k,(npw_k,9,mpsang*mpsang*psps%useylm)) 511 else 512 ABI_ALLOCATE(ylmgr_k,(0,0,0)) 513 end if 514 515 ABI_ALLOCATE(kg_k,(3,npw_k)) 516 kg_k(:,:) = 0 517 !$OMP PARALLEL DO 518 do ipw=1,npw_k 519 kg_k(1,ipw)=kg(1,ipw+ikg) 520 kg_k(2,ipw)=kg(2,ipw+ikg) 521 kg_k(3,ipw)=kg(3,ipw+ikg) 522 end do 523 if (psps%useylm==1) then 524 !SOMP PARALLEL DO COLLAPSE(2) 525 do ilm=1,mpsang*mpsang 526 do ipw=1,npw_k 527 ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm) 528 end do 529 end do 530 if(rfstrs/=0.or.need_becfr.or.need_piezofr.or.need_efmas)then 531 !SOMP PARALLEL DO COLLAPSE(3) 532 do ilm=1,mpsang*mpsang 533 do ii=1,9 534 do ipw=1,npw_k 535 ylmgr_k(ipw,ii,ilm)=ylmgr(ipw+ikg,ii,ilm) 536 end do 537 end do 538 end do 539 end if 540 end if 541 542 cplex=2;if (istwf_k>1) cplex=1 543 544 ! Compute (k+G) vectors (only if useylm=1) 545 nkpg=0 546 if (rfstrs/=0.or.need_efmas.or.pawpiezo==1) nkpg=3*dtset%nloalg(3) 547 ABI_ALLOCATE(kpg_k,(npw_k,nkpg)) 548 if (nkpg>0) then 549 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 550 end if 551 552 !EFMAS: Compute second order derivatives w/r to k for all direction for this k-point. 553 if (need_efmas) then 554 ABI_ALLOCATE(ddkinpw,(npw_k,3,3)) 555 do mu=1,3 556 do nu=1,3 557 ! call d2kpg(ddkinpw(:,mu,nu),dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,mu,nu,kg_k,kpoint,npw_k) 558 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,ddkinpw(:,mu,nu),kpoint,npw_k,mu,nu) 559 end do 560 end do 561 end if 562 563 ! Compute nonlocal form factors ffnl at all (k+G): 564 ider=0;dimffnl=1; 565 if(need_becfr) then 566 ider=1;dimffnl=4 567 end if 568 if(rfstrs/=0.or.need_piezofr.or.need_efmas)then 569 ider=2;dimffnl=3+7*psps%useylm 570 end if 571 ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,psps%ntypat)) 572 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,& 573 & gmet,gprimd,ider,idir_ffnl,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,& 574 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,& 575 & psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k) 576 577 ! For piezoelectric tensor need additional ffnl derivatives 578 if(need_piezofr)then 579 ider_str=1 ; dimffnl_str=2 580 ABI_ALLOCATE(ffnl_str,(npw_k,dimffnl_str,psps%lmnmax,psps%ntypat,6)) 581 do mu=1,6 !loop over strain 582 idir_str=-mu 583 call mkffnl(psps%dimekb,dimffnl_str,psps%ekb,ffnl_str(:,:,:,:,mu),& 584 & psps%ffspl,gmet,gprimd,ider_str,idir_str,psps%indlmn,kg_k,kpg_k,& 585 & kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,& 586 & psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k) 587 end do 588 end if 589 590 ! Load k-dependent part in the Hamiltonian datastructure 591 ABI_ALLOCATE(ph3d,(2,npw_k,gs_ham%matblk)) 592 call load_k_hamiltonian(gs_ham,kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,& 593 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_ph3d=.true.) 594 595 ! Initialize contributions from current k point 596 if(rfphon==1) dyfrnlk(:,:)=zero 597 if(rfstrs/=0)then 598 enlk=zero;eltfrnlk(:,:)=zero 599 end if 600 if (need_becfr) becfrnlk(:,:,:)=zero 601 if (need_piezofr) piezofrnlk(:,:)=zero 602 if(need_efmas) then 603 call check_degeneracies(efmasdeg(ikpt),dtset%efmas_bands(:,ikpt),nband_k,eigen(bdtot_index+1:bdtot_index+nband_k), & 604 & dtset%efmas_deg_tol) 605 do ideg=1,efmasdeg(ikpt)%ndegs 606 if(efmasdeg(ikpt)%treated(ideg)) then 607 ABI_MALLOC(efmasfr(ikpt,ideg)%ch2c,(efmasdeg(ikpt)%deg_dim(ideg),efmasdeg(ikpt)%deg_dim(ideg),3,3)) 608 efmasfr(ikpt,ideg)%ch2c=zero 609 else 610 ABI_MALLOC(efmasfr(ikpt,ideg)%ch2c,(0,0,0,0)) 611 end if 612 end do 613 end if 614 615 ! Loop over bands 616 do iband=1,nband_k 617 618 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me)) then 619 cycle 620 end if 621 622 occ_k=occ(iband+bdtot_index) 623 cwavef(:,1:npw_k*dtset%nspinor) = cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 624 625 ! Compute non-local contributions from n,k 626 if (psps%usepaw==1) eig_k=eigen(iband+bdtot_index) 627 628 ! === Dynamical matrix 629 if(rfphon==1) then 630 call nonlop(choice_phon,cpopt,cwaveprj,enlout_phon,gs_ham,idir,(/eig_k/),mpi_enreg,1,& 631 & nnlout_phon,paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef) 632 ! Accumulate non-local contributions from n,k 633 dyfrnlk(:,:)=dyfrnlk(:,:)+occ_k*reshape(enlout_phon(:),(/6,natom/)) 634 end if 635 636 ! === Elastic tensor 637 if(rfstrs/=0) then 638 call nonlop(choice_strs,cpopt,cwaveprj,enlout_strs,gs_ham,idir,(/eig_k/),mpi_enreg,1,& 639 & nnlout_strs,paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef) 640 ! Accumulate non-local contribut ions from n,k 641 eltfrnlk(:,:)=eltfrnlk(:,:)+occ_k*reshape(enlout_strs(:),(/3*natom+6,6/)) 642 end if !endo if strs 643 644 ! PAW: accumulate gradients of rhoij 645 !EFMAS: Bug with efmas currently; to be looked into... 646 if (psps%usepaw==1.and.(.not.need_efmas)) then 647 call pawaccrhoij(gs_ham%atindx,cplex,cwaveprj,cwaveprj,0,isppol,natom,& 648 & natom,dtset%nspinor,occ_k,3,pawrhoij_tot,usetimerev,wtk_k) 649 end if 650 651 ! PAW: Compute frozen contribution to piezo electric tensor 652 if (need_piezofr) then 653 do ii=1,3 ! Loop over elect. field directions 654 call nonlop(choice_piez3,cpopt,cwaveprj,enlout_piez1,gs_ham,0,(/zero/),mpi_enreg,1,& 655 & nnlout_piez1,paw_opt_1,signs,nonlop_dum,tim_nonlop,cwavef,cwavef,enl=becij(:,:,:,ii)) 656 piezofrnlk(:,ii)=piezofrnlk(:,ii)+occ_k*enlout_piez1(:) 657 end do !end do ii 658 end if 659 660 ! PAW: Compute frozen contribution to Born Effective Charges 661 if (need_becfr) then 662 do ii=1,3 ! Loop over elect. field directions 663 call nonlop(choice_bec2,cpopt,cwaveprj,enlout_bec1,gs_ham,0,(/zero/),mpi_enreg,1,& 664 & nnlout_bec1,paw_opt_1,signs,nonlop_dum,tim_nonlop,cwavef,cwavef,enl=becij(:,:,:,ii)) 665 becfrnlk(:,:,ii)=becfrnlk(:,:,ii)+occ_k*reshape(enlout_bec1(:),(/3,natom/)) 666 end do !end do ii 667 end if 668 669 if (need_becfr.or.need_piezofr) then 670 do ii=1,3 ! Loop over elect. field directions 671 ! Not able to compute if ipert=(Elect. field) and no ddk WF file 672 if (ddkfil(ii)==0) cycle 673 ! Read ddk wave function 674 ABI_ALLOCATE(ddk,(2,npw_k*dtset%nspinor)) 675 if (ddkfil(ii)/=0) then 676 call wfk_read_bks(ddkfiles(ii), iband, ik_ddk(ii), isppol, xmpio_single, cg_bks=ddk) 677 ! Multiply ddk by +i 678 do jj=1,npw_k*dtset%nspinor 679 arg=ddk(1,jj) 680 ddk(1,jj)=-ddk(2,jj);ddk(2,jj)=arg 681 end do 682 else 683 ddk=zero 684 end if 685 686 if(need_becfr)then 687 do iatom=1,natom !Loop over atom 688 ia=gs_ham%atindx(iatom) 689 do mu=1,3 !loop over atom direction 690 call nonlop(choice_bec2,cpopt_bec,cwaveprj,enlout_bec1,gs_ham,mu,(/zero/),& 691 & mpi_enreg,1,nnlout_bec1,paw_opt_3,signs_field,svectout,tim_nonlop,& 692 & cwavef,cwavef,iatom_only=iatom) 693 call dotprod_g(dotprod(1),dotprod(2),istwf_k,npw_k*dtset%nspinor,2,svectout,ddk,& 694 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 695 becfrnlk(mu,ia,ii)=becfrnlk(mu,ia,ii)+occ_k*dotprod(1) 696 end do 697 end do 698 end if 699 700 if(need_piezofr)then 701 do mu=1,6 !loop over strain 702 call load_k_hamiltonian(gs_ham,ffnl_k=ffnl_str(:,:,:,:,mu)) 703 call nonlop(choice_piez3,cpopt,cwaveprj,enlout_piez1,gs_ham,mu,(/zero/),mpi_enreg,1,& 704 & nnlout_piez1,paw_opt_3,signs_field,svectout,tim_nonlop,cwavef,svectout) 705 call dotprod_g(dotprod(1),dotprod(2),istwf_k,npw_k*dtset%nspinor,2,svectout,ddk,& 706 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 707 piezofrnlk(mu,ii)=piezofrnlk(mu,ii)+occ_k*dotprod(1) 708 end do 709 call load_k_hamiltonian(gs_ham,ffnl_k=ffnl) 710 end if 711 712 ABI_DEALLOCATE(ddk) 713 end do ! End loop ddk file 714 end if 715 716 if(need_piezofr)then 717 enlout_piez2 = zero 718 call nonlop(choice_piez55,cpopt,cwaveprj,enlout_piez2,gs_ham,0,(/zero/),mpi_enreg,1,& 719 & nnlout_piez2,paw_opt_3,signs,nonlop_dum,tim_nonlop,cwavef,cwavef) 720 ! Multiply enlout by +i 721 iashift = 1 722 do mu=1,6 ! strain 723 do nu=1,3 ! k 724 piezofrnlk(mu,nu)=piezofrnlk(mu,nu)-occ_k*(enlout_piez2(iashift+1)) ! Real part 725 ! piezofrnlk(mu,nu)=piezofrnlk(mu,nu)+occ_k*(enlout_piez2(iashift ))! Imaginary part 726 iashift = iashift + 2 727 end do 728 end do 729 end if 730 731 if(need_becfr)then 732 call nonlop(choice_bec54,cpopt,cwaveprj,enlout_bec2,gs_ham,0,(/zero/),mpi_enreg,1,& 733 & nnlout_bec2,paw_opt_3,signs,nonlop_dum,tim_nonlop,cwavef,cwavef) 734 ! Multiply enlout by +i 735 iashift = 1 736 do iatom=1,natom ! atm 737 do mu=1,3 ! atm pos. 738 do nu=1,3 ! k 739 becfrnlk(mu,iatom,nu)=becfrnlk(mu,iatom,nu)-occ_k*(enlout_bec2(iashift+1)) ! Real part 740 ! becfrnlk(mu,iatom,nu)=becfrnlk(mu,iatom,nu)+occ_k*(enlout_bec2(iashift ))! Imaginary part 741 iashift = iashift + 2 742 end do 743 end do 744 end do 745 end if 746 747 if(need_efmas) then 748 if(iband>=efmasdeg(ikpt)%band_range(1) .and. iband<=(efmasdeg(ikpt)%band_range(2))) then 749 choice_efmas=8; signs=2 750 cpopt=-1 !To prevent re-use of stored dgxdt, which are not for all direction required for EFMAS. 751 paw_opt_efmas=0; if(psps%usepaw/=0) paw_opt_efmas=4 !To get both gh2c and gs2c 752 nnlout_efmas=0; tim_nonlop=0 ! No tim_nonlop for efmas, currently. 753 do mu=1,3 754 do nu=1,3 755 idir=3*(mu-1)+nu !xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9, (xyz,xyz)=(mu,nu) 756 gh2c=zero; gs2c=zero 757 call nonlop(choice_efmas,cpopt,cwaveprj,enlout_efmas,gs_ham,idir,(/eig_k/),mpi_enreg,& 758 1,nnlout_efmas,paw_opt_efmas,signs,gs2c,tim_nonlop,cwavef,gh2c) 759 do ispinor=1,dtset%nspinor 760 ii = 1+(ispinor-1)*npw_k 761 do icplx=1,2 762 gh2c(icplx,ii:ispinor*npw_k) = gh2c(icplx,ii:ispinor*npw_k) + & 763 & ddkinpw(1:npw_k,mu,nu)*cwavef(icplx,ii:ispinor*npw_k) 764 end do 765 end do 766 gh2c = gh2c - eig_k*gs2c 767 ideg = efmasdeg(ikpt)%ideg(iband) 768 do jband=efmasdeg(ikpt)%degs_bounds(1,ideg),efmasdeg(ikpt)%degs_bounds(2,ideg) 769 cg_left(:,1:npw_k*dtset%nspinor) = cg(:,1+(jband-1)*npw_k*dtset%nspinor+icg:jband*npw_k*dtset%nspinor+icg) 770 dotprod=0 771 call dotprod_g(dotprod(1),dotprod(2),istwf_k,npw_k*dtset%nspinor,2,cg_left,gh2c,mpi_enreg%me_g0,& 772 & mpi_enreg%comm_spinorfft) 773 isub = iband-efmasdeg(ikpt)%degl(ideg) 774 jsub = jband-efmasdeg(ikpt)%degl(ideg) 775 efmasfr(ikpt,ideg)%ch2c(jsub,isub,mu,nu)=cmplx(dotprod(1),dotprod(2),kind=dpc) 776 end do 777 end do 778 end do 779 end if 780 end if 781 782 end do ! End of loop on bands 783 784 if(rfphon==1) then 785 do iatom=1,natom 786 ia=iatom;if (dyfr_nondiag==0) ia=1 787 dyfrnl(1,1,1,iatom,ia)=dyfrnl(1,1,1,iatom,ia)+wtk_k*dyfrnlk(1,iatom) 788 dyfrnl(1,2,2,iatom,ia)=dyfrnl(1,2,2,iatom,ia)+wtk_k*dyfrnlk(2,iatom) 789 dyfrnl(1,3,3,iatom,ia)=dyfrnl(1,3,3,iatom,ia)+wtk_k*dyfrnlk(3,iatom) 790 dyfrnl(1,2,3,iatom,ia)=dyfrnl(1,2,3,iatom,ia)+wtk_k*dyfrnlk(4,iatom) 791 dyfrnl(1,1,3,iatom,ia)=dyfrnl(1,1,3,iatom,ia)+wtk_k*dyfrnlk(5,iatom) 792 dyfrnl(1,1,2,iatom,ia)=dyfrnl(1,1,2,iatom,ia)+wtk_k*dyfrnlk(6,iatom) 793 end do 794 end if ! end if rfphon 795 if(rfstrs/=0) then 796 eltfrnl(:,:)=eltfrnl(:,:)+dtset%wtk(ikpt)*eltfrnlk(:,:) 797 end if 798 if(need_becfr) then 799 becfrnl(:,:,:)=becfrnl(:,:,:)+dtset%wtk(ikpt)*becfrnlk(:,:,:) 800 end if 801 if(need_piezofr) then 802 piezofrnl(:,:)=piezofrnl(:,:)+dtset%wtk(ikpt)*piezofrnlk(:,:) 803 end if 804 ! Increment indexes 805 bdtot_index=bdtot_index+nband_k 806 if (dtset%mkmem/=0) then 807 ibg=ibg+nband_k*dtset%nspinor 808 icg=icg+npw_k*dtset%nspinor*nband_k 809 ikg=ikg+npw_k 810 end if 811 812 ABI_DEALLOCATE(ffnl) 813 ABI_DEALLOCATE(kpg_k) 814 ABI_DEALLOCATE(ph3d) 815 ABI_DEALLOCATE(ylm_k) 816 ABI_DEALLOCATE(ylmgr_k) 817 ABI_DEALLOCATE(cwavef) 818 ABI_DEALLOCATE(kg_k) 819 if (need_becfr.or.need_piezofr) then 820 ABI_DEALLOCATE(svectout) 821 end if 822 if (need_piezofr) then 823 ABI_DEALLOCATE(ffnl_str) 824 end if 825 if (need_efmas) then 826 ABI_DEALLOCATE(ddkinpw) 827 ABI_DEALLOCATE(cg_left) 828 ABI_DEALLOCATE(gh2c) 829 ABI_DEALLOCATE(gs2c) 830 end if 831 832 end do ! End loops on isppol and ikpt 833 end do 834 if(rfphon==1) then 835 ABI_DEALLOCATE(dyfrnlk) 836 ABI_DEALLOCATE(enlout_phon) 837 end if 838 if(rfstrs/=0) then 839 ABI_DEALLOCATE(eltfrnlk) 840 ABI_DEALLOCATE(enlout_strs) 841 end if 842 if (need_becfr) then 843 ABI_DEALLOCATE(becfrnlk) 844 ABI_DEALLOCATE(enlout_bec1) 845 ABI_DEALLOCATE(enlout_bec2) 846 end if 847 if(need_piezofr)then 848 ABI_DEALLOCATE(enlout_piez1) 849 ABI_DEALLOCATE(enlout_piez2) 850 ABI_DEALLOCATE(piezofrnlk) 851 end if 852 if(need_efmas) then 853 ABI_DEALLOCATE(enlout_efmas) 854 end if 855 if (psps%usepaw==1) then 856 if (need_becfr.or.need_piezofr) then 857 ABI_DEALLOCATE(becij) 858 end if 859 call pawcprj_free(cwaveprj) 860 end if 861 ABI_DATATYPE_DEALLOCATE(cwaveprj) 862 863 !Fill in lower triangle of matrixes 864 if (rfphon==1) then 865 do iatom=1,natom 866 ia=iatom;if (dyfr_nondiag==0) ia=1 867 dyfrnl(1,3,2,iatom,ia)=dyfrnl(1,2,3,iatom,ia) 868 dyfrnl(1,3,1,iatom,ia)=dyfrnl(1,1,3,iatom,ia) 869 dyfrnl(1,2,1,iatom,ia)=dyfrnl(1,1,2,iatom,ia) 870 end do 871 end if 872 if(rfstrs/=0)then 873 do jj=2,6 874 do ii=1,jj-1 875 eltfrnl(jj,ii)=eltfrnl(ii,jj) 876 end do 877 end do 878 end if 879 880 !Parallel case: accumulate (n,k) contributions 881 if (xmpi_paral==1) then 882 call timab(48,1,tsec) 883 ! Accumulate dyfrnl 884 if(rfphon==1)then 885 call xmpi_sum(dyfrnl,spaceworld,ierr) 886 end if 887 ! Accumulate eltfrnl. 888 if(rfstrs/=0)then 889 call xmpi_sum(eltfrnl,spaceworld,ierr) 890 end if 891 ! Accumulate becfrnl 892 if (need_becfr) then 893 call xmpi_sum(becfrnl,spaceworld,ierr) 894 end if 895 ! Accumulate piezofrnl 896 if (need_piezofr) then 897 call xmpi_sum(piezofrnl,spaceworld,ierr) 898 end if 899 900 ! PAW: accumulate gradients of rhoij 901 if (psps%usepaw==1) then 902 ABI_ALLOCATE(dimlmn,(natom)) 903 dimlmn(1:natom)=pawrhoij_tot(1:natom)%cplex*pawrhoij_tot(1:natom)%lmn2_size 904 bufdim=ncpgr*sum(dimlmn)*nsploop 905 ABI_ALLOCATE(mpibuf,(bufdim)) 906 ii=0;mpibuf=zero 907 do iatom=1,natom 908 do isppol=1,nsploop 909 do mu=1,ncpgr 910 mpibuf(ii+1:ii+dimlmn(iatom))=pawrhoij_tot(iatom)%grhoij(mu,1:dimlmn(iatom),isppol) 911 ii=ii+dimlmn(iatom) 912 end do 913 end do 914 end do 915 call xmpi_sum(mpibuf,spaceworld,ierr) 916 ii=0 917 do iatom=1,natom 918 do isppol=1,nsploop 919 do mu=1,ncpgr 920 pawrhoij_tot(iatom)%grhoij(mu,1:dimlmn(iatom),isppol)=mpibuf(ii+1:ii+dimlmn(iatom)) 921 ii=ii+dimlmn(iatom) 922 end do 923 end do 924 end do 925 ABI_DEALLOCATE(mpibuf) 926 ABI_DEALLOCATE(dimlmn) 927 end if 928 call timab(48,2,tsec) 929 end if 930 931 !====== PAW: Additional steps 932 if (psps%usepaw==1) then 933 934 ! Symmetrize rhoij gradients and transfer to cartesian (reciprocal space) coord. 935 ! This symetrization is necessary in the antiferromagnetic case... 936 if (rfphon==1.and.rfstrs==0) then 937 option_rhoij=2;option=0 938 call symrhoij(pawrhoij_tot,pawrhoij_tot,option_rhoij,gprimd,indsym,0,natom,dtset%nsym,& 939 & psps%ntypat,option,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,dtset%typat,& 940 & comm_atom=my_comm_atom, mpi_atmtab=my_atmtab) 941 else if (rfphon==1.and.rfstrs==1) then 942 option_rhoij=23;option=0 943 call symrhoij(pawrhoij_tot,pawrhoij_tot,option_rhoij,gprimd,indsym,0,natom,dtset%nsym,& 944 & psps%ntypat,option,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,dtset%typat,& 945 & comm_atom=my_comm_atom, mpi_atmtab=my_atmtab) 946 end if 947 948 ! Translate coordinates 949 do iatom=1,natom 950 cplx=pawrhoij_tot(iatom)%cplex 951 do isppol=1,nsploop 952 do klmn=1,pawrhoij_tot(iatom)%lmn2_size 953 do ii=1,cplx 954 if(rfphon==1.or.rfstrs/=0)then 955 grhoij(1:3)=pawrhoij_tot(iatom)%grhoij(shift_rhoij+1:shift_rhoij+3,cplx*(klmn-1)+ii,isppol) 956 do mu=1,3 957 pawrhoij_tot(iatom)%grhoij(shift_rhoij+mu,cplx*(klmn-1)+ii,isppol)=gprimd(mu,1)*grhoij(1)& 958 & +gprimd(mu,2)*grhoij(2)+gprimd(mu,3)*grhoij(3) 959 end do 960 end if 961 if(rfstrs/=0)then 962 call strconv(pawrhoij_tot(iatom)%grhoij(1:6,cplx*(klmn-1)+ii,isppol),gprimd,& 963 & pawrhoij_tot(iatom)%grhoij(1:6,cplx*(klmn-1)+ii,isppol)) 964 end if 965 end do 966 end do 967 end do 968 end do 969 970 ! In case of elastic tensor computation, add diagonal contribution: 971 ! -delta_{alphabeta} rhoi_{ij} to drhoij/d_eps 972 if(rfstrs/=0)then 973 do iatom=1,natom 974 cplx=pawrhoij_tot(iatom)%cplex 975 do isppol=1,nsploop 976 do nu=1,pawrhoij_tot(iatom)%nrhoijsel 977 klmn=pawrhoij_tot(iatom)%rhoijselect(nu) 978 do ii=1,cplx 979 pawrhoij_tot(iatom)%grhoij(1:3,cplx*(klmn-1)+ii,isppol)= & 980 & pawrhoij_tot(iatom)%grhoij(1:3,cplx*(klmn-1)+ii,isppol)& 981 & -pawrhoij_tot(iatom)%rhoijp(cplx*(nu-1)+ii,isppol) 982 end do 983 end do 984 end do 985 end do 986 end if 987 988 ! Add gradients due to Dij derivatives to dynamical matrix/stress tensor 989 dimnhat=0;optgr=0;optgr2=0;optstr=0;optstr2=0 990 if (rfphon==1) optgr2=1 991 if (rfstrs/=0) optstr2=1 992 ABI_ALLOCATE(nhat_dum,(1,0)) 993 call pawgrnl(gs_ham%atindx1,dimnhat,dyfrnl,dyfr_cplex,eltfrnl,dummy,gsqcut,mgfftf,my_natom,natom,& 994 & gs_ham%nattyp,nfftf,ngfftf,nhat_dum,dummy,dtset%nspden,dtset%nsym,psps%ntypat,optgr,optgr2,optstr,optstr2,& 995 & pawang,pawfgrtab,pawrhoij_tot,pawtab,ph1df,psps,dtset%qptn,rprimd,symrec,dtset%typat,ucvol,vtrial,vxc,xred,& 996 & mpi_atmtab=my_atmtab,comm_atom=my_comm_atom) 997 ABI_DEALLOCATE(nhat_dum) 998 end if !PAW 999 1000 !The indexing array atindx is used to reestablish the correct order of atoms 1001 if (rfstrs/=0)then 1002 ABI_ALLOCATE(elt_work,(6+3*natom,6)) 1003 elt_work(1:6,1:6)=eltfrnl(1:6,1:6) 1004 do ia=1,natom 1005 ielt=7+3*(ia-1) 1006 ieltx=7+3*(gs_ham%atindx(ia)-1) 1007 elt_work(ielt:ielt+2,1:6)=eltfrnl(ieltx:ieltx+2,1:6) 1008 end do 1009 eltfrnl(:,:)=elt_work(:,:) 1010 ABI_DEALLOCATE(elt_work) 1011 end if 1012 1013 !Born Effective Charges and PAW: 1014 !1-Re-order atoms -- 2-Add diagonal contribution from rhoij 1015 !3-Multiply by -1 because that the effective charges 1016 ! are minus the second derivatives of the energy 1017 if (need_becfr) then 1018 ABI_ALLOCATE(becfrnl_tmp,(3,natom,3)) 1019 becfrnl_tmp=-becfrnl 1020 do ia=1,natom ! Atom (sorted by type) 1021 iatom=gs_ham%atindx1(ia) ! Atom (not sorted) 1022 itypat=dtset%typat(iatom) 1023 do ii=1,3 ! Direction of electric field 1024 do jj=1,3 ! Direction of atom 1025 becfrnl(jj,iatom,ii)=becfrnl_tmp(jj,ia,ii) 1026 end do 1027 end do 1028 end do 1029 ABI_DEALLOCATE(becfrnl_tmp) 1030 end if 1031 1032 !Piezoelectric Tensor 1033 !-Multiply by -1 because that the piezoelectric tensor 1034 ! are minus the second derivatives of the energy 1035 if (need_piezofr) then 1036 piezofrnl=-piezofrnl 1037 end if 1038 1039 !Close the ddk files 1040 do ii=1,3 1041 call wfk_close(ddkfiles(ii)) 1042 end do 1043 1044 !Release now useless memory 1045 if (psps%usepaw==1) then 1046 do iatom=1,natom 1047 ABI_DEALLOCATE(pawrhoij_tot(iatom)%grhoij) 1048 pawrhoij_tot(iatom)%ngrhoij=0 1049 end do 1050 if (paral_atom) then 1051 call pawrhoij_free(pawrhoij_tot) 1052 ABI_DATATYPE_DEALLOCATE(pawrhoij_tot) 1053 end if 1054 end if 1055 call destroy_hamiltonian(gs_ham) 1056 call timab(159,2,tsec) 1057 1058 write(msg,'(3a)')ch10,' ==> Calculation of the frozen part of the second order derivative done',ch10 1059 call wrtout(std_out,msg,'COLL') 1060 1061 DBG_EXIT("COLL") 1062 1063 end subroutine d2frnl