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