TABLE OF CONTENTS
ABINIT/pawgrnl [ Functions ]
NAME
pawgrnl
FUNCTION
PAW: Add to GRadients of total energy due to non-local term of Hamiltonian the contribution due to Dij derivatives In particular, compute contribution to forces, stresses, dyn. matrix Remember: Vnl=Sum_ij[|p_i>Dij<p_j|]
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (FJ, 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx dimnhat=second dimension of array nhat (0 or # of spin components) distribfft<type(distribfft_type)>=--optional-- contains all the informations related to the FFT parallelism and plane sharing dyfr_cplex=1 if dyfrnl is real, 2 if it is complex gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere mgfft=maximum size of 1D FFTs me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present) mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms comm_fft=--optional-- MPI communicator over FFT components (=mpi_comm_grid is not present) mpi_comm_grid=--optional-- MPI communicator over real space grid components (=comm_fft is not present) my_natom=number of atoms treated by current processor natom=total number of atoms in cell nattyp(ntypat)=array describing how many atoms of each type in cell nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nhat(nfft,dimnhat)=compensation charge density on rectangular grid in real space nspden=number of spin-density components nsym=number of symmetries in space group ntypat=number of types of atoms optgr= 1 if gradients with respect to atomic position(s) have to be computed optgr2= 1 if 2nd gradients with respect to atomic position(s) have to be computed optstr= 1 if gradients with respect to strain(s) have to be computed optstr2= 1 if 2nd gradients with respect to strain(s) have to be computed paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present) pawang <type(pawang_type)>=paw angular mesh and related data pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) information psps <type(pseudopotential_type)>=variables related to pseudopotentials qphon(3)=wavevector of the phonon rprimd(3,3)=dimensional primitive translations in real space (bohr) symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates typat(natom)=types of atoms ucvol=unit cell volume vtrial(nfft,nspden)= total local potential vxc(nfft,nspden)=XC potential xred(3,natom)=reduced dimensionless atomic coordinates
SIDE EFFECTS
At input, this terms contain contribution from non-local projectors derivatives At output, they are updated with the contribution of Dij derivatives ==== if optgr=1 ==== grnl(3*natom) =gradients of NL energy wrt atomic coordinates ==== if optstr=1 ==== nlstr(6) =gradients of NL energy wrt strains ==== if optgr2=1 ==== dyfrnl(dyfr_cplex,3,3,natom,natom) =2nd gradients of NL energy wrt atomic coordinates ==== if optstr=2 ==== eltfrnl(6+3*natom,6)=non-symmetrized non-local contribution to the elastic tensor
NOTES
In the case of parallelisation over atoms and calculation of dynamical matrix (optgr2=1) several data are gathered and no more distributed inside this routine.
PARENTS
d2frnl,etotfor,forstr
CHILDREN
convert_notation,destroy_distribfft,dfpt_atm2fft,free_my_atmtab get_my_atmtab,init_distribfft_seq,metric,pawexpiqr,pawfgrtab_free pawfgrtab_gather,pawfgrtab_nullify,pawgylm,pawrfgd_fft,pawrhoij_free pawrhoij_gather,pawrhoij_nullify,stresssym,xmpi_sum
SOURCE
88 #if defined HAVE_CONFIG_H 89 #include "config.h" 90 #endif 91 92 #include "abi_common.h" 93 94 subroutine pawgrnl(atindx1,dimnhat,dyfrnl,dyfr_cplex,eltfrnl,grnl,gsqcut,mgfft,my_natom,natom,& 95 & nattyp,nfft,ngfft,nhat,nlstr,nspden,nsym,ntypat,optgr,optgr2,optstr,optstr2,& 96 & pawang,pawfgrtab,pawrhoij,pawtab,ph1d,psps,qphon,rprimd,symrec,typat,ucvol,vtrial,vxc,xred,& 97 & mpi_atmtab,comm_atom,comm_fft,mpi_comm_grid,me_g0,paral_kgb,distribfft) ! optional arguments (parallelism) 98 99 use defs_basis 100 use defs_datatypes 101 use defs_abitypes 102 use m_profiling_abi 103 use m_xmpi 104 use m_errors 105 106 use m_geometry, only : metric 107 use m_distribfft, only : distribfft_type,init_distribfft_seq,destroy_distribfft 108 use m_pawang, only : pawang_type 109 use m_pawtab, only : pawtab_type 110 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_nullify, pawfgrtab_gather 111 use m_pawrhoij, only : pawrhoij_type, pawrhoij_free, pawrhoij_gather, pawrhoij_nullify 112 use m_paw_finegrid, only : pawgylm, pawrfgd_fft, pawexpiqr 113 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 114 115 !This section has been created automatically by the script Abilint (TD). 116 !Do not modify the following lines by hand. 117 #undef ABI_FUNC 118 #define ABI_FUNC 'pawgrnl' 119 use interfaces_41_geometry 120 use interfaces_64_psp 121 use interfaces_65_paw, except_this_one => pawgrnl 122 !End of the abilint section 123 124 implicit none 125 126 !Arguments ------------------------------------ 127 !scalars 128 integer,intent(in) :: dimnhat,dyfr_cplex,mgfft,my_natom,natom,nfft,nspden,nsym,ntypat 129 integer,intent(in) :: optgr,optgr2,optstr,optstr2 130 integer,optional,intent(in) :: me_g0,comm_atom,comm_fft,mpi_comm_grid,paral_kgb 131 real(dp),intent(in) :: gsqcut,ucvol 132 type(distribfft_type),optional,target,intent(in) :: distribfft 133 type(pawang_type),intent(in) :: pawang 134 type(pseudopotential_type),intent(in) :: psps 135 !arrays 136 integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18) 137 integer,intent(in) :: symrec(3,3,nsym),typat(natom) 138 integer,optional,target,intent(in) :: mpi_atmtab(:) 139 real(dp),intent(in) :: nhat(nfft,dimnhat),ph1d(2,3*(2*mgfft+1)*natom),qphon(3) 140 real(dp),intent(in) :: rprimd(3,3),vxc(nfft,nspden),xred(3,natom) 141 real(dp),intent(in),target :: vtrial(nfft,nspden) 142 real(dp),intent(inout) :: dyfrnl(dyfr_cplex,3,3,natom,natom*optgr2) 143 real(dp),intent(inout) :: eltfrnl(6+3*natom,6),grnl(3*natom*optgr) 144 real(dp),intent(inout) :: nlstr(6*optstr) 145 type(pawfgrtab_type),target,intent(inout) :: pawfgrtab(:) 146 type(pawrhoij_type),target,intent(inout) :: pawrhoij(:) 147 type(pawtab_type),intent(in) :: pawtab(ntypat) 148 149 !Local variables------------------------------- 150 !scalars 151 integer :: bufind,bufsiz,cplex,dimvtrial,eps_alpha,eps_beta,eps_gamma,eps_delta,iatm,iatom 152 integer :: iatom_pawfgrtab,iatom_pawrhoij,iatom_tot,iatshft,ic,idiag,idir,ier,ilm,indx,irhoij 153 integer :: isel,ishift_grhoij,ishift_gr,ishift2_gr,ishift_gr2,ishift_str,ishift_str2,ishift_str2is,ispden 154 integer :: ispvtr,itypat,jatom,jatom_tot,jatm,jc,jrhoij,jtypat,klm,klmn,klmn1,ll,lm_size 155 integer :: lm_sizej,lmax,lmin,lmn2_size,me_fft,mu,mua,mub,mushift,my_me_g0,my_comm_atom,my_comm_fft 156 integer :: my_comm_grid,my_paral_kgb,n1,n2,n3,nfftot,nfgd,nfgd_jatom 157 integer :: ngrad,ngrad_nondiag,ngradp,ngradp_nondiag,ngrhat,nsploop 158 integer :: opt1,opt2,opt3,qne0,usexcnhat 159 logical,parameter :: save_memory=.true. 160 logical :: has_phase,my_atmtab_allocated 161 logical :: paral_atom,paral_atom_pawfgrtab,paral_atom_pawrhoij,paral_grid 162 real(dp) :: dlt_tmp,fact_ucvol,grhat_x,hatstr_diag,rcut_jatom,ro,ro_d,ucvol_ 163 character(len=500) :: msg 164 type(distribfft_type),pointer :: my_distribfft 165 type(pawfgrtab_type),pointer :: pawfgrtab_iatom,pawfgrtab_jatom 166 type(pawrhoij_type),pointer :: pawrhoij_iatom,pawrhoij_jatom 167 !arrays 168 integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/) 169 integer,parameter :: eps1(6)=(/1,2,3,2,3,1/),eps2(6)=(/1,2,3,3,1,2/) 170 integer,parameter :: mu9(9)=(/1,2,3,4,5,6,4,5,6/) 171 integer,allocatable :: atindx(:),atm_indx(:),mu4(:) 172 integer,allocatable,target :: ifftsph_tmp(:) 173 integer,ABI_CONTIGUOUS pointer :: ffti3_local(:),fftn3_distrib(:),ifft_jatom(:) 174 integer, pointer :: my_atmtab(:) 175 real(dp) :: gmet(3,3),gprimd(3,3),hatstr(6),rdum(1),rdum2(1),rmet(3,3),tmp(12) 176 real(dp) :: work1(dyfr_cplex,3,3),work2(dyfr_cplex,3,3) 177 real(dp),allocatable :: buf(:,:),buf1(:),dyfr(:,:,:,:,:),eltfr(:,:) 178 real(dp),allocatable :: grhat_tmp(:,:),grhat_tmp2(:,:),hatgr(:) 179 real(dp),allocatable :: prod(:,:),prodp(:,:),vloc(:),vpsp1_gr(:,:),vpsp1_str(:,:) 180 real(dp),allocatable,target :: rfgd_tmp(:,:) 181 real(dp),ABI_CONTIGUOUS pointer :: gylm_jatom(:,:),gylmgr_jatom(:,:,:),gylmgr2_jatom(:,:,:),expiqr_jatom(:,:) 182 real(dp),ABI_CONTIGUOUS pointer :: rfgd_jatom(:,:),vtrial_(:,:) 183 type(coeff2_type),allocatable :: prod_nondiag(:),prodp_nondiag(:) 184 type(pawfgrtab_type),pointer :: pawfgrtab_(:),pawfgrtab_tot(:) 185 type(pawrhoij_type),pointer :: pawrhoij_(:),pawrhoij_tot(:) 186 187 ! ************************************************************************* 188 189 DBG_ENTER("COLL") 190 191 !Compatibility tests 192 qne0=0;if (qphon(1)**2+qphon(2)**2+qphon(3)**2>=1.d-15) qne0=1 193 if (my_natom>0) then 194 if ((optgr2==1.or.optstr2==1).and.pawrhoij(1)%ngrhoij==0) then 195 msg='Inconsistency between variables optgr2/optstr2 and ngrhoij!' 196 MSG_BUG(msg) 197 end if 198 if (pawfgrtab(1)%rfgd_allocated==0) then 199 if ((optgr2==1.and.qne0==1).or.optstr2==1) then 200 MSG_BUG('pawfgrtab()%rfgd array must be allocated!') 201 end if 202 end if 203 end if 204 205 !---------------------------------------------------------------------- 206 !Parallelism setup 207 208 !Set up parallelism over atoms 209 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 210 paral_atom_pawfgrtab=(size(pawfgrtab)/=natom) 211 paral_atom_pawrhoij=(size(pawrhoij)/=natom) 212 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 213 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 214 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 215 if (paral_atom) then 216 ABI_ALLOCATE(atm_indx,(natom)) 217 atm_indx=-1 218 do iatom=1,my_natom 219 atm_indx(my_atmtab(iatom))=iatom 220 end do 221 end if 222 223 !Set up parallelism over real space grid and/or FFT 224 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nfftot=n1*n2*n3 225 my_comm_grid=xmpi_comm_self;my_comm_fft=xmpi_comm_self;me_fft=0 226 my_me_g0=1;my_paral_kgb=0;paral_grid=.false.;nullify(my_distribfft) 227 if (present(mpi_comm_grid).or.present(comm_fft)) then 228 if (present(mpi_comm_grid)) my_comm_grid=mpi_comm_grid 229 if (present(comm_fft)) my_comm_fft=comm_fft 230 if (.not.present(mpi_comm_grid)) my_comm_grid=comm_fft 231 if (.not.present(comm_fft)) my_comm_fft=mpi_comm_grid 232 paral_grid=(xmpi_comm_size(my_comm_grid)>1) 233 me_fft=xmpi_comm_rank(my_comm_fft) 234 end if 235 if (optgr2==1.or.optstr2==1) then 236 if (present(comm_fft)) then 237 if ((.not.present(paral_kgb)).or.(.not.present(me_g0)).or.(.not.present(distribfft))) then 238 MSG_BUG(' Need paral_kgb, me_g0 and distribfft with comm_fft !') 239 end if 240 my_me_g0=me_g0;my_paral_kgb=paral_kgb 241 my_distribfft => distribfft 242 else 243 ABI_DATATYPE_ALLOCATE(my_distribfft,) 244 call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp') 245 end if 246 if (n2 == my_distribfft%n2_coarse) then 247 fftn3_distrib => my_distribfft%tab_fftdp3_distrib 248 ffti3_local => my_distribfft%tab_fftdp3_local 249 else 250 fftn3_distrib => my_distribfft%tab_fftdp3dg_distrib 251 ffti3_local => my_distribfft%tab_fftdp3dg_local 252 end if 253 else 254 nullify(my_distribfft,fftn3_distrib,ffti3_local) 255 end if 256 257 !---------------------------------------------------------------------- 258 !Initializations 259 260 !Compute different geometric tensors 261 !ucvol is not computed here but provided as input arg 262 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol_) 263 fact_ucvol=ucvol/dble(nfftot) 264 265 !Retrieve local potential according to the use of nhat in XC 266 usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat) 267 if (usexcnhat==0) then 268 ABI_ALLOCATE(vtrial_,(nfft,1)) 269 dimvtrial=1 270 !$OMP PARALLEL DO PRIVATE(ic) SHARED(nfft,vtrial,vtrial_,vxc) 271 do ic=1,nfft 272 vtrial_(ic,1)=vtrial(ic,1)-vxc(ic,1) 273 end do 274 else 275 dimvtrial=nspden 276 vtrial_ => vtrial 277 end if 278 279 !Initializations and allocations 280 ngrhat=0;ngrad=0;ngradp=0;ngrad_nondiag=0;ngradp_nondiag=0 281 ishift_grhoij=0;ishift_gr=0;ishift_gr2=0;ishift_str=0;ishift_str2=0;ishift_str2is=0;ishift2_gr=0 282 cplex=1;if (qne0==1) cplex=2 283 if (optgr==1) then 284 ABI_ALLOCATE(hatgr,(3*natom)) 285 hatgr=zero 286 ngrad=ngrad+3 287 ngrhat=ngrhat+3 288 ishift_gr2=ishift_gr2+3 289 end if 290 if (optgr2==1) then 291 mu=min(dyfr_cplex,cplex) 292 ngrad =ngrad +9 293 ngradp=ngradp+3 294 ngrad_nondiag = ngrad_nondiag +9*mu 295 ngradp_nondiag= ngradp_nondiag+3*mu 296 ngrhat= ngrhat+9*mu 297 end if 298 if (optstr==1) then 299 hatstr=zero 300 ngrad=ngrad+6 301 ngrhat=ngrhat+6 302 ishift_gr=ishift_gr+6 303 ishift_gr2=ishift_gr2+6 304 ishift_str2=ishift_str2+6 305 ishift_str2is = ishift_str2is+6 306 end if 307 if (optstr2==1) then 308 ngrad =ngrad+6*(6+3) 309 ngradp=ngradp+(6+3) 310 ngrad_nondiag =ngrad_nondiag+6*(6+3) 311 ngradp_nondiag=ngradp_nondiag+3 312 ishift2_gr=ishift2_gr+3 313 ngrhat=ngrhat+6*(6+3) 314 ishift_gr=ishift_gr+(6+3) 315 ishift_gr2=ishift_gr2+6*(6+3) 316 ishift_str2is=ishift_str2is+36 317 ishift_grhoij = 6 318 end if 319 320 nsploop=nspden;if (dimvtrial<nspden) nsploop=2 321 if (optgr2/=1.and.optstr2/=1) then 322 ABI_ALLOCATE(grhat_tmp,(ngrhat,1)) 323 else 324 ABI_ALLOCATE(grhat_tmp,(ngrhat,natom)) 325 grhat_tmp=zero 326 ABI_DATATYPE_ALLOCATE(prod_nondiag,(natom)) 327 ABI_DATATYPE_ALLOCATE(prodp_nondiag,(natom)) 328 ABI_ALLOCATE(atindx,(natom)) 329 if(optgr2==1.or.optstr2==1)then 330 ABI_ALLOCATE(vpsp1_gr,(cplex*nfft,3)) 331 vpsp1_gr(:,:)= zero 332 end if 333 if (optgr2==1) then 334 ABI_ALLOCATE(dyfr,(dyfr_cplex,3,3,natom,natom)) 335 dyfr=zero 336 end if 337 if (optstr2==1) then 338 ABI_ALLOCATE(vpsp1_str,(cplex*nfft,6)) 339 ABI_ALLOCATE(grhat_tmp2,(18,natom)) 340 ABI_ALLOCATE(eltfr,(6+3*natom,6)) 341 eltfr=zero 342 end if 343 ABI_ALLOCATE(mu4,(4)) 344 atindx(:)=0 345 do iatom=1,natom 346 iatm=0 347 do while (atindx(iatom)==0.and.iatm<natom) 348 iatm=iatm+1;if (atindx1(iatm)==iatom) atindx(iatom)=iatm 349 end do 350 end do 351 end if 352 353 !The computation of dynamical matrix and elastic tensor requires the knowledge of 354 !g_l(r-R).Y_lm(r-R) and derivatives for all atoms 355 !Compute them here, except memory saving is activated 356 if ((.not.save_memory).and.(optgr2==1.or.optstr2==1)) then 357 do jatom=1,size(pawfgrtab) 358 jatom_tot=jatom;if (paral_atom_pawfgrtab) jatom_tot=my_atmtab(jatom) 359 pawfgrtab_jatom => pawfgrtab(jatom) 360 lm_sizej=pawfgrtab_jatom%l_size**2 361 opt1=0;opt2=0;opt3=0 362 if (pawfgrtab_jatom%gylm_allocated==0) then 363 if (allocated(pawfgrtab_jatom%gylm)) then 364 ABI_DEALLOCATE(pawfgrtab_jatom%gylm) 365 end if 366 ABI_ALLOCATE(pawfgrtab_jatom%gylm,(pawfgrtab_jatom%nfgd,lm_sizej)) 367 pawfgrtab_jatom%gylm_allocated=2;opt1=1 368 end if 369 if (pawfgrtab_jatom%gylmgr_allocated==0) then 370 if (allocated(pawfgrtab_jatom%gylmgr)) then 371 ABI_DEALLOCATE(pawfgrtab_jatom%gylmgr) 372 end if 373 ABI_ALLOCATE(pawfgrtab_jatom%gylmgr,(3,pawfgrtab_jatom%nfgd,lm_sizej)) 374 pawfgrtab_jatom%gylmgr_allocated=2;opt2=1 375 end if 376 if (opt1+opt2+opt3>0) then 377 call pawgylm(pawfgrtab_jatom%gylm,pawfgrtab_jatom%gylmgr,& 378 & pawfgrtab_jatom%gylmgr2,lm_sizej,pawfgrtab_jatom%nfgd,& 379 & opt1,opt2,opt3,pawtab(typat(jatom_tot)),pawfgrtab_jatom%rfgd) 380 end if 381 if (optgr2==1.and.qne0==1) then 382 if (pawfgrtab_jatom%expiqr_allocated==0) then 383 if (allocated(pawfgrtab_jatom%expiqr)) then 384 ABI_DEALLOCATE(pawfgrtab_jatom%expiqr) 385 end if 386 pawfgrtab_jatom%expiqr_allocated=2 387 ABI_ALLOCATE(pawfgrtab_jatom%expiqr,(2,nfgd)) 388 call pawexpiqr(pawfgrtab_jatom%expiqr,gprimd,pawfgrtab_jatom%nfgd,& 389 & qphon,pawfgrtab_jatom%rfgd,xred(:,jatom_tot)) 390 end if 391 end if 392 end do 393 end if 394 395 !The computation of dynamical matrix and elastic tensor might require some communications 396 if ((optgr2==1.or.optstr2==1).and.paral_atom.and.paral_atom_pawfgrtab.and.(.not.save_memory)) then 397 ABI_DATATYPE_ALLOCATE(pawfgrtab_tot,(natom)) 398 call pawfgrtab_nullify(pawfgrtab_tot) 399 call pawfgrtab_gather(pawfgrtab,pawfgrtab_tot,my_comm_atom,ier,mpi_atmtab=my_atmtab) 400 else 401 pawfgrtab_tot => pawfgrtab 402 end if 403 if ((optgr2==1.or.optstr2==1).and.paral_atom.and.paral_atom_pawrhoij) then 404 ABI_DATATYPE_ALLOCATE(pawrhoij_tot,(natom)) 405 call pawrhoij_nullify(pawrhoij_tot) 406 call pawrhoij_gather(pawrhoij,pawrhoij_tot,-1,my_comm_atom, & 407 & with_rhoijres=.false.,with_rhoij_=.false.,with_lmnmix=.false.) 408 else 409 pawrhoij_tot => pawrhoij 410 end if 411 412 if (save_memory) then 413 pawfgrtab_ => pawfgrtab 414 pawrhoij_ => pawrhoij 415 else 416 pawfgrtab_ => pawfgrtab_tot 417 pawrhoij_ => pawrhoij_tot 418 end if 419 420 !---------------------------------------------------------------------- 421 !Loops over types and atoms 422 423 iatshft=0 424 do itypat=1,ntypat 425 426 lmn2_size=pawtab(itypat)%lmn2_size 427 lm_size=pawtab(itypat)%lcut_size**2 428 429 do iatm=iatshft+1,iatshft+nattyp(itypat) 430 431 iatom_tot=atindx1(iatm) 432 iatom=iatom_tot 433 if (paral_atom) then 434 if (save_memory.or.(optgr2/=1.and.optstr2/=1)) iatom=atm_indx(iatom_tot) 435 end if 436 437 if (iatom==-1) cycle 438 iatom_pawfgrtab=iatom_tot;if (paral_atom_pawfgrtab) iatom_pawfgrtab=iatom 439 iatom_pawrhoij =iatom_tot;if (paral_atom_pawrhoij) iatom_pawrhoij =iatom 440 pawfgrtab_iatom => pawfgrtab_(iatom_pawfgrtab) 441 pawrhoij_iatom => pawrhoij_(iatom_pawrhoij) 442 443 idiag=1;if (optgr2==1.or.optstr2==1) idiag=iatm 444 nfgd=pawfgrtab_iatom%nfgd 445 446 ABI_ALLOCATE(vloc,(nfgd)) 447 if (ngrad>0) then 448 ABI_ALLOCATE(prod,(ngrad,lm_size)) 449 end if 450 if (ngradp>0) then 451 ABI_ALLOCATE(prodp,(ngradp,lm_size)) 452 end if 453 if (ngrad_nondiag>0.and.ngradp_nondiag>0) then 454 do jatm=1,natom 455 jtypat=typat(atindx1(jatm)) 456 lm_sizej=pawtab(jtypat)%lcut_size**2 457 ABI_ALLOCATE(prod_nondiag(jatm)%value,(ngrad_nondiag,lm_sizej)) 458 ABI_ALLOCATE(prodp_nondiag(jatm)%value,(ngradp_nondiag,lm_sizej)) 459 end do 460 end if 461 462 grhat_tmp=zero 463 if(optstr2==1) grhat_tmp2=zero 464 465 ! ------------------------------------------------------------------ 466 ! Compute some useful data 467 468 ! Eventually compute g_l(r).Y_lm(r) derivatives for the current atom (if not already done) 469 if ((optgr==1.or.optstr==1).and.(optgr2/=1).and.(optstr2/=1)) then 470 if (pawfgrtab_iatom%gylmgr_allocated==0) then 471 if (allocated(pawfgrtab_iatom%gylmgr)) then 472 ABI_DEALLOCATE(pawfgrtab_iatom%gylmgr) 473 end if 474 ABI_ALLOCATE(pawfgrtab_iatom%gylmgr,(3,pawfgrtab_iatom%nfgd,lm_size)) 475 pawfgrtab_iatom%gylmgr_allocated=2 476 call pawgylm(rdum,pawfgrtab_iatom%gylmgr,rdum2,lm_size,pawfgrtab_iatom%nfgd,& 477 & 0,1,0,pawtab(itypat),pawfgrtab_iatom%rfgd) 478 end if 479 480 end if 481 if (optgr2==1.or.optstr2==1) then 482 opt1=0;opt2=0;opt3=0 483 if (pawfgrtab_iatom%gylm_allocated==0) then 484 if (allocated(pawfgrtab_iatom%gylm)) then 485 ABI_DEALLOCATE(pawfgrtab_iatom%gylm) 486 end if 487 ABI_ALLOCATE(pawfgrtab_iatom%gylm,(pawfgrtab_iatom%nfgd,lm_size)) 488 pawfgrtab_iatom%gylm_allocated=2;opt1=1 489 end if 490 if (pawfgrtab_iatom%gylmgr_allocated==0) then 491 if (allocated(pawfgrtab_iatom%gylmgr)) then 492 ABI_DEALLOCATE(pawfgrtab_iatom%gylmgr) 493 end if 494 ABI_ALLOCATE(pawfgrtab_iatom%gylmgr,(3,pawfgrtab_iatom%nfgd,lm_size)) 495 pawfgrtab_iatom%gylmgr_allocated=2;opt2=1 496 end if 497 if (pawfgrtab_iatom%gylmgr2_allocated==0) then 498 if (allocated(pawfgrtab_iatom%gylmgr2)) then 499 ABI_DEALLOCATE(pawfgrtab_iatom%gylmgr2) 500 end if 501 ABI_ALLOCATE(pawfgrtab_iatom%gylmgr2,(6,pawfgrtab_iatom%nfgd,lm_size)) 502 pawfgrtab_iatom%gylmgr2_allocated=2;opt3=1 503 end if 504 if (opt1+opt2+opt3>0) then 505 call pawgylm(pawfgrtab_iatom%gylm,pawfgrtab_iatom%gylmgr,& 506 & pawfgrtab_iatom%gylmgr2,lm_size,pawfgrtab_iatom%nfgd,& 507 & opt1,opt2,opt3,pawtab(itypat),pawfgrtab_iatom%rfgd) 508 end if 509 end if 510 511 ! Eventually compute exp(-i.q.r) factors for the current atom (if not already done) 512 if (optgr2==1.and.qne0==1.and.(pawfgrtab_iatom%expiqr_allocated==0)) then 513 if (allocated(pawfgrtab_iatom%expiqr)) then 514 ABI_DEALLOCATE(pawfgrtab_iatom%expiqr) 515 end if 516 ABI_ALLOCATE(pawfgrtab_iatom%expiqr,(2,nfgd)) 517 call pawexpiqr(pawfgrtab_iatom%expiqr,gprimd,nfgd,qphon,& 518 & pawfgrtab_iatom%rfgd,xred(:,iatom)) 519 pawfgrtab_iatom%expiqr_allocated=2 520 end if 521 has_phase=(optgr2==1.and.pawfgrtab_iatom%expiqr_allocated/=0) 522 523 ! Eventually compute 1st-order potential 524 if (optgr2==1.or.optstr2==1) then 525 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,iatom_tot,& 526 & mgfft,psps%mqgrid_vl,natom,3,nfft,ngfft,ntypat,ph1d,& 527 & psps%qgrid_vl,qphon,typat,ucvol,psps%usepaw,xred,psps,pawtab,atmvlocr1=vpsp1_gr,& 528 & vspl=psps%vlspl,comm_fft=my_comm_fft,me_g0=my_me_g0,& 529 & paral_kgb=my_paral_kgb,distribfft=my_distribfft) 530 if (cplex==1) then 531 do ic=1,nfft 532 tmp(1:3)=vpsp1_gr(ic,1:3) 533 do mu=1,3 534 vpsp1_gr(ic,mu)=-(gprimd(mu,1)*tmp(1)+gprimd(mu,2)*tmp(2)+gprimd(mu,3)*tmp(3)) 535 end do 536 end do 537 else ! cplex=2 538 do ic=1,nfft 539 jc=2*ic;tmp(1:3)=vpsp1_gr(jc-1,1:3);tmp(4:6)=vpsp1_gr(jc,1:3) 540 do mu=1,3 541 vpsp1_gr(jc-1,mu)=-(gprimd(mu,1)*tmp(1)+gprimd(mu,2)*tmp(2)+gprimd(mu,3)*tmp(3)) 542 vpsp1_gr(jc ,mu)=-(gprimd(mu,1)*tmp(4)+gprimd(mu,2)*tmp(5)+gprimd(mu,3)*tmp(6)) 543 end do 544 end do 545 end if 546 end if 547 if (optstr2==1) then 548 vpsp1_str(:,:) = zero 549 call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,natom+3,& 550 & mgfft,psps%mqgrid_vl,natom,6,nfft,ngfft,ntypat,& 551 & ph1d,psps%qgrid_vl,qphon,typat,ucvol,psps%usepaw,xred,psps,pawtab,atmvlocr1=vpsp1_str,& 552 & vspl=psps%vlspl,comm_fft=my_comm_fft,me_g0=my_me_g0,& 553 & paral_kgb=my_paral_kgb,distribfft=my_distribfft) 554 end if 555 556 ! ------------------------------------------------------------------ 557 ! Loop over spin components 558 559 do ispden=1,nsploop 560 561 ! ----- Retrieve potential (subtle if nspden=4 ;-) 562 if (nspden/=4) then 563 ispvtr=min(dimvtrial,ispden) 564 do ic=1,nfgd 565 jc = pawfgrtab_iatom%ifftsph(ic) 566 vloc(ic)=vtrial_(jc,ispvtr) 567 end do 568 else 569 if (ispden==1) then 570 ispvtr=min(dimvtrial,2) 571 do ic=1,nfgd 572 jc=pawfgrtab_iatom%ifftsph(ic) 573 vloc(ic)=half*(vtrial_(jc,1)+vtrial_(jc,ispvtr)) 574 end do 575 else if (ispden==4) then 576 ispvtr=min(dimvtrial,2) 577 do ic=1,nfgd 578 jc=pawfgrtab_iatom%ifftsph(ic) 579 vloc(ic)=half*(vtrial_(jc,1)-vtrial_(jc,ispvtr)) 580 end do 581 else if (ispden==2) then 582 ispvtr=min(dimvtrial,3) 583 do ic=1,nfgd 584 jc=pawfgrtab_iatom%ifftsph(ic) 585 vloc(ic)=vtrial_(jc,ispvtr) 586 end do 587 else ! ispden=3 588 ispvtr=min(dimvtrial,4) 589 do ic=1,nfgd 590 jc=pawfgrtab_iatom%ifftsph(ic) 591 vloc(ic)=-vtrial_(jc,ispvtr) 592 end do 593 end if 594 end if 595 596 ! ----------------------------------------------------------------------- 597 ! ----- Compute projected scalars (integrals of vloc and Q_ij^hat) ------ 598 ! ----- and/or their derivatives ---------------------------------------- 599 600 if (ngrad>0) prod=zero 601 if (ngradp>0) prodp=zero 602 603 ! ==== Contribution to forces ==== 604 if (optgr==1) then 605 do ilm=1,lm_size 606 do ic=1,pawfgrtab_iatom%nfgd 607 do mu=1,3 608 prod(mu+ishift_gr,ilm)=prod(mu+ishift_gr,ilm)-& 609 & vloc(ic)*pawfgrtab_iatom%gylmgr(mu,ic,ilm) 610 end do 611 end do 612 end do 613 end if ! optgr 614 615 ! ==== Contribution to stresses ==== 616 if (optstr==1) then 617 do ilm=1,lm_size 618 do ic=1,pawfgrtab_iatom%nfgd 619 jc=pawfgrtab_iatom%ifftsph(ic) 620 do mu=1,6 621 mua=alpha(mu);mub=beta(mu) 622 prod(mu+ishift_str,ilm)=prod(mu+ishift_str,ilm) & 623 & +half*vloc(ic)& 624 & *(pawfgrtab_iatom%gylmgr(mua,ic,ilm)*pawfgrtab_iatom%rfgd(mub,ic)& 625 & +pawfgrtab_iatom%gylmgr(mub,ic,ilm)*pawfgrtab_iatom%rfgd(mua,ic)) 626 end do 627 end do 628 end do 629 end if ! optstr 630 631 ! ==== Diagonal contribution to frozen wf part of dyn. matrix ==== 632 if (optgr2==1) then 633 ! Diagonal contribution 634 do ilm=1,lm_size 635 do ic=1,pawfgrtab_iatom%nfgd 636 do mu=1,9 637 prod(ishift_gr2+mu,ilm)=prod(ishift_gr2+mu,ilm) & 638 & +half*vloc(ic)*pawfgrtab_iatom%gylmgr2(mu9(mu),ic,ilm) 639 end do 640 do mu=1,3 641 prodp(ishift_gr+mu,ilm)=prodp(ishift_gr+mu,ilm) & 642 & -vloc(ic)*pawfgrtab_iatom%gylmgr(mu,ic,ilm) 643 end do 644 end do 645 end do 646 end if ! optgr2 647 648 ! ==== Diagonal contribution to elastic tensor ==== 649 if (optstr2==1) then 650 do ilm=1,lm_size 651 do ic=1,pawfgrtab_iatom%nfgd 652 mu=1 653 jc=pawfgrtab_iatom%ifftsph(ic) 654 do mua=1,6 655 eps_alpha=eps1(mua);eps_beta=eps2(mua); 656 do mub=1,6 657 eps_gamma=eps1(mub);eps_delta=eps2(mub); 658 mu4 = zero 659 call convert_notation(mu4,eps_alpha,eps_beta,eps_gamma,eps_delta) 660 ! v_loc*d2glylm 661 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) + half*half*vloc(ic)*( & 662 & pawfgrtab_iatom%rfgd(eps_beta,ic)*pawfgrtab_iatom%rfgd(eps_gamma,ic)*& 663 & pawfgrtab_iatom%gylmgr2(mu4(2),ic,ilm)& 664 & +pawfgrtab_iatom%rfgd(eps_alpha,ic)*pawfgrtab_iatom%rfgd(eps_gamma,ic)*& 665 pawfgrtab_iatom%gylmgr2(mu4(4),ic,ilm)& 666 & +pawfgrtab_iatom%rfgd(eps_beta,ic) *pawfgrtab_iatom%rfgd(eps_delta,ic)*& 667 pawfgrtab_iatom%gylmgr2(mu4(1),ic,ilm)& 668 & +pawfgrtab_iatom%rfgd(eps_alpha,ic)*pawfgrtab_iatom%rfgd(eps_delta,ic)*& 669 pawfgrtab_iatom%gylmgr2(mu4(3),ic,ilm)) 670 if(eps_gamma==eps_beta)then 671 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) & 672 & +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic)) 673 end if 674 if(eps_gamma==eps_alpha)then 675 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) & 676 & +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic)) 677 end if 678 if(eps_delta==eps_beta)then 679 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) & 680 & +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic)) 681 end if 682 if(eps_delta==eps_alpha)then 683 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) & 684 & +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic)) 685 end if 686 ! d(vloc)/d(eps_gammadelta) * d(gylm)/d(eps_alphabeta) 687 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm)& 688 & +vpsp1_str(jc,mub)*half*(& 689 & (pawfgrtab_iatom%gylmgr(eps_alpha,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic)& 690 & +pawfgrtab_iatom%gylmgr(eps_beta,ic,ilm) *pawfgrtab_iatom%rfgd(eps_alpha,ic))) 691 ! d(vloc)/d(eps_alphabeta) * d(gylm)/d(eps_gammadelta) 692 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm)& 693 & +vpsp1_str(jc,mua)*half*(& 694 & (pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_delta,ic)& 695 & +pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_gamma,ic))) 696 ! delta_alphabeta * dv_loc/depsgammadelta * (gylm) 697 if (mua<=3) then 698 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) & 699 & +vpsp1_str(jc,mub)*pawfgrtab_iatom%gylm(ic,ilm) 700 end if 701 ! delta_gammadelta * dv_loc/depsalphabeta * (gylm) 702 if (mub<=3) then 703 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) & 704 & +vpsp1_str(jc,mua) * pawfgrtab_iatom%gylm(ic,ilm) 705 end if 706 ! delta_gammadelta * v_loc * d(gylm)/d(eps_alphabeta) 707 if (mub<=3) then 708 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) & 709 & +half*vloc(ic)& 710 & *(pawfgrtab_iatom%gylmgr(eps_beta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic)& 711 & + pawfgrtab_iatom%gylmgr(eps_alpha,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic)) 712 end if 713 ! delta_alphabeta * v_loc * d(gylm)/d(eps_gammadelta) 714 if (mua<=3) then 715 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) & 716 & +half*vloc(ic)& 717 & *(pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_delta,ic)& 718 & + pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_gamma,ic)) 719 end if 720 ! delta_gammadelta delta_alphabeta * v_loc * (gylm) 721 if (mua<=3.and.mub<=3) then 722 prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) & 723 & +vloc(ic)*pawfgrtab_iatom%gylm(ic,ilm) 724 end if 725 mu=mu+1 726 end do !end loop mub 727 end do !end loop mua 728 ! vloc * d(gylm)/d(eps_alphabeta) 729 do mu=1,6 730 mua=alpha(mu);mub=beta(mu) 731 prodp(ishift_str2+mu,ilm)=prodp(ishift_str2+mu,ilm)& 732 & +half*vloc(ic)*& 733 & (pawfgrtab_iatom%gylmgr(mua,ic,ilm)*pawfgrtab_iatom%rfgd(mub,ic)& 734 & +pawfgrtab_iatom%gylmgr(mub,ic,ilm)*pawfgrtab_iatom%rfgd(mua,ic)) 735 ! d(vloc)/d(eps_alphabeta or gammadelta) * gylm 736 prodp(ishift_str2+mu,ilm)=prodp(ishift_str2+mu,ilm)& 737 & +vpsp1_str(jc,mu)*pawfgrtab_iatom%gylm(ic,ilm) 738 ! delta_alphabeta * vloc * gylm 739 if (mu<=3) then 740 prodp(ishift_str2+mu,ilm)=prodp(ishift_str2+mu,ilm)& 741 & +vloc(ic)*pawfgrtab_iatom%gylm(ic,ilm) 742 end if 743 744 ! INTERNAL STRAIN CONTRIBUTION: 745 do idir=1,3 746 ! v_loc*d2glylm/dR contribution: 747 eps_alpha=alpha(mu);eps_beta=beta(mu); 748 call convert_notation(mu4,eps_alpha,eps_beta,idir,idir) 749 prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)& 750 & -half*vloc(ic)& 751 & *(pawfgrtab_iatom%gylmgr2(mu4(3),ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic)& 752 & +pawfgrtab_iatom%gylmgr2(mu4(1),ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic)) 753 if (idir==eps_beta)then 754 prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)& 755 & -half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_alpha,ic,ilm)) 756 end if 757 if (idir==eps_alpha)then 758 prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)& 759 & -half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_beta,ic,ilm)) 760 end if 761 ! delta_gammadelta * v_loc * d(gylm)/dR 762 if (mu<=3) then 763 prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)-& 764 vloc(ic)*pawfgrtab_iatom%gylmgr(idir,ic,ilm) 765 end if 766 ! dv_loc/deps_alph_beta * d(gylm)/dR 767 prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)-& 768 vpsp1_str(jc,mu)*pawfgrtab_iatom%gylmgr(idir,ic,ilm) 769 end do 770 end do 771 do idir=1,3 772 ! v_loc * d(gylm)/dR 773 prodp(6+idir,ilm) = prodp(6+idir,ilm)-vloc(ic)*pawfgrtab_iatom%gylmgr(idir,ic,ilm) 774 end do !end loop idir 775 ! END INTERNAL STRAIN CONTRIBUTION 776 777 end do 778 end do 779 end if !optstr2 780 781 ! Off-diagonal contributions 782 if (optgr2==1.or.optstr2==1) then 783 do jatm=1,natom 784 jatom_tot=atindx1(jatm);jtypat=typat(jatom_tot) 785 jatom=jatom_tot;if (paral_atom.and.save_memory) jatom=atm_indx(jatom_tot) 786 lm_sizej=pawtab(jtypat)%lcut_size**2 787 788 ! Retrieve data for the atom j 789 if (save_memory.and.jatom/=iatom) then 790 rcut_jatom=pawtab(jtypat)%rshp 791 call pawrfgd_fft(ifftsph_tmp,gmet,n1,n2,n3,nfgd_jatom,rcut_jatom,rfgd_tmp,rprimd,& 792 & ucvol,xred(:,jatom_tot),fft_distrib=fftn3_distrib,fft_index=ffti3_local,me_fft=me_fft) 793 ifft_jatom => ifftsph_tmp ; rfgd_jatom => rfgd_tmp 794 ABI_ALLOCATE(gylm_jatom,(nfgd_jatom,lm_sizej)) 795 ABI_ALLOCATE(gylmgr_jatom,(3,nfgd_jatom,lm_sizej)) 796 opt1=1;opt2=1;opt3=0;gylmgr2_jatom=>gylmgr_jatom 797 call pawgylm(gylm_jatom,gylmgr_jatom,gylmgr2_jatom,lm_sizej,nfgd_jatom,& 798 & opt1,opt2,opt3,pawtab(typat(jatom_tot)),rfgd_jatom) 799 if (optgr2==1.and.qne0==1) then 800 ABI_ALLOCATE(expiqr_jatom,(2,nfgd_jatom)) 801 call pawexpiqr(expiqr_jatom,gprimd,nfgd_jatom,qphon,rfgd_jatom,xred(:,jatom_tot)) 802 end if 803 else 804 pawfgrtab_jatom => pawfgrtab_tot(jatom) 805 nfgd_jatom = pawfgrtab_jatom%nfgd 806 ifft_jatom => pawfgrtab_jatom%ifftsph 807 rfgd_jatom => pawfgrtab_jatom%rfgd 808 gylm_jatom => pawfgrtab_jatom%gylm 809 gylmgr_jatom => pawfgrtab_jatom%gylmgr 810 gylmgr2_jatom => pawfgrtab_jatom%gylmgr2 811 expiqr_jatom => pawfgrtab_jatom%expiqr 812 end if 813 814 ! ==== Off-diagonal contribution to frozen wf part of dyn. matrix ==== 815 if (optgr2==1) then 816 mu = min(dyfr_cplex,cplex) 817 prod_nondiag(jatm)%value(ishift_gr2+1:ishift_gr2+(9*mu),:) = zero 818 prodp_nondiag(jatm)%value(ishift2_gr+1:ishift2_gr +(3*mu),:) = zero 819 if (has_phase.or.cplex==2) then 820 if (dyfr_cplex==1.or.cplex==1) then 821 do ilm=1,lm_sizej 822 do ic=1,nfgd_jatom 823 jc=2*ifft_jatom(ic) 824 tmp(1:3)=vpsp1_gr(jc-1,1:3)*expiqr_jatom(1,ic) & 825 & -vpsp1_gr(jc ,1:3)*expiqr_jatom(2,ic) 826 do mu=1,9 827 mua=alpha(mu);mub=beta(mu) 828 prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+mu,ilm) & 829 & +tmp(mua)*gylmgr_jatom(mub,ic,ilm) 830 end do 831 do mu=1,3 832 prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm) & 833 & -tmp(mu)*gylm_jatom(ic,ilm) 834 end do 835 end do 836 end do 837 else 838 do ilm=1,lm_sizej 839 do ic=1,nfgd_jatom 840 jc=2*ifft_jatom(ic) 841 tmp(1:3)=vpsp1_gr(jc-1,1:3)*expiqr_jatom(1,ic) & 842 & -vpsp1_gr(jc ,1:3)*expiqr_jatom(2,ic) 843 tmp(4:6)=vpsp1_gr(jc-1,1:3)*expiqr_jatom(2,ic) & 844 & +vpsp1_gr(jc ,1:3)*expiqr_jatom(1,ic) 845 do mu=1,9 846 mua=alpha(mu);mub=beta(mu) 847 prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+mu,ilm) & 848 & +tmp(mua )*gylmgr_jatom(mub,ic,ilm) 849 prod_nondiag(jatm)%value(ishift_gr2+9+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+9+mu,ilm) & 850 & +tmp(3+mua)*gylmgr_jatom(mub,ic,ilm) 851 end do 852 do mu=1,3 853 prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm) & 854 & -tmp( mu)*gylm_jatom(ic,ilm) 855 prodp_nondiag(jatm)%value(ishift2_gr+3+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+3+mu,ilm) & 856 & -tmp(3+mu)*gylm_jatom(ic,ilm) 857 end do 858 end do 859 end do 860 end if 861 else ! no phase 862 do ilm=1,lm_sizej 863 do ic=1,nfgd_jatom 864 jc=ifft_jatom(ic) 865 do mu=1,9 866 mua=alpha(mu);mub=beta(mu) 867 prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+mu,ilm) & 868 & +vpsp1_gr(jc,mua)*gylmgr_jatom(mub,ic,ilm) 869 end do 870 do mu=1,3 871 prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm) & 872 & -vpsp1_gr(jc,mu)*gylm_jatom(ic,ilm) 873 end do 874 end do 875 end do 876 end if 877 end if ! optgr2 878 879 ! ==== Off-diagonal contribution to elastic tensor ==== 880 if (optstr2==1) then 881 prod_nondiag(jatm)%value(ishift_str2is+1:ishift_str2is+18,:)=zero 882 prodp_nondiag(jatm)%value(1:3,:)=zero 883 do ilm=1,lm_sizej 884 do ic=1,nfgd_jatom 885 mu=1;jc=ifft_jatom(ic) 886 ! INTERNAL STRAIN CONTRIBUTION: 887 do mua=1,6 888 eps_alpha=eps1(mua);eps_beta=eps2(mua); 889 ! d(-vloc)/dR * d(gylm)/d(eps_gamma_delta) 890 do idir=1,3 891 prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)=& 892 & prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)& 893 & -vpsp1_gr(jc,idir)*half& 894 & *(gylmgr_jatom(eps_alpha,ic,ilm) * rfgd_jatom(eps_beta ,ic)& 895 & +gylmgr_jatom(eps_beta ,ic,ilm) * rfgd_jatom(eps_alpha,ic)) 896 ! delta_alphabeta * d(-v_loc/dr) * gylm 897 if (mua<=3) then 898 prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)=& 899 & prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)& 900 & -vpsp1_gr(jc,idir)*gylm_jatom(ic,ilm) 901 end if 902 end do ! dir 903 end do ! mua 904 do idir=1,3 905 ! d(-v_loc/dr) * gylm 906 prodp_nondiag(jatm)%value(idir,ilm) = prodp_nondiag(jatm)%value(idir,ilm)& 907 & -vpsp1_gr(jc,idir)*gylm_jatom(ic,ilm) 908 end do !end loop idir 909 ! END INTERNAL STRAIN CONTRIBUTION 910 end do 911 end do 912 end if ! optstr2 913 914 ! Release temp memory allocated for atom j 915 if (save_memory.and.jatom/=iatom) then 916 ABI_DEALLOCATE(ifftsph_tmp) 917 ABI_DEALLOCATE(rfgd_tmp) 918 ABI_DEALLOCATE(gylm_jatom) 919 ABI_DEALLOCATE(gylmgr_jatom) 920 if (optgr2==1.and.qne0==1) then 921 ABI_DEALLOCATE(expiqr_jatom) 922 end if 923 end if 924 925 end do ! loop on atoms j 926 end if ! optgr2 or optstr2 927 928 ! --- Apply scaling factor on integrals --- 929 if (ngrad >0) prod (:,:)=prod (:,:)*fact_ucvol 930 if (ngradp>0) prodp(:,:)=prodp(:,:)*fact_ucvol 931 if (ngrad_nondiag>0) then 932 do jatm=1,natom 933 prod_nondiag(jatm)%value(:,:)=prod_nondiag(jatm)%value(:,:)*fact_ucvol 934 end do 935 end if 936 if (ngradp_nondiag>0) then 937 do jatm=1,natom 938 prodp_nondiag(jatm)%value(:,:)=prodp_nondiag(jatm)%value(:,:)*fact_ucvol 939 end do 940 end if 941 942 ! --- Reduction in case of parallelization --- 943 if (paral_grid) then 944 if (ngrad>0) then 945 call xmpi_sum(prod,my_comm_grid,ier) 946 end if 947 if (ngradp>0) then 948 call xmpi_sum(prodp,my_comm_grid,ier) 949 end if 950 if (ngrad_nondiag>0.or.ngradp_nondiag>0) then 951 bufsiz=0;bufind=0 952 do jatm=1,natom 953 jtypat=typat(atindx1(jatm)) 954 bufsiz=bufsiz+pawtab(jtypat)%lcut_size**2 955 end do 956 ABI_ALLOCATE(buf,(ngrad_nondiag+ngradp_nondiag,bufsiz)) 957 do jatm=1,natom 958 jtypat=typat(atindx1(jatm)) 959 lm_sizej=pawtab(jtypat)%lcut_size**2 960 if (ngrad_nondiag> 0) buf(1:ngrad_nondiag,bufind+1:bufind+lm_sizej)= & 961 & prod_nondiag(jatm)%value(:,:) 962 if (ngradp_nondiag>0) buf(ngrad_nondiag+1:ngrad_nondiag+ngradp_nondiag, & 963 & bufind+1:bufind+lm_sizej)=prodp_nondiag(jatm)%value(:,:) 964 bufind=bufind+lm_sizej*(ngrad_nondiag+ngradp_nondiag) 965 end do 966 call xmpi_sum(buf,my_comm_grid,ier) 967 bufind=0 968 do jatm=1,natom 969 jtypat=typat(atindx1(jatm)) 970 lm_sizej=pawtab(jtypat)%lcut_size**2 971 if (ngrad> 0) prod_nondiag(jatm)%value(:,:)= & 972 & buf(1:ngrad_nondiag,bufind+1:bufind+lm_sizej) 973 if (ngradp>0) prodp_nondiag(jatm)%value(:,:)= & 974 & buf(ngrad_nondiag+1:ngrad_nondiag+ngradp_nondiag,bufind+1:bufind+lm_sizej) 975 bufind=bufind+lm_sizej*(ngrad_nondiag+ngradp_nondiag) 976 end do 977 ABI_DEALLOCATE(buf) 978 end if 979 end if 980 981 ! ---------------------------------------------------------------- 982 ! Compute final sums (i.e. derivatives of Sum_ij[rho_ij.Intg{Qij.Vloc}] 983 984 ! ---- Compute terms common to all gradients 985 jrhoij=1 986 do irhoij=1,pawrhoij_iatom%nrhoijsel 987 klmn=pawrhoij_iatom%rhoijselect(irhoij) 988 klm =pawtab(itypat)%indklmn(1,klmn) 989 lmin=pawtab(itypat)%indklmn(3,klmn) 990 lmax=pawtab(itypat)%indklmn(4,klmn) 991 ro =pawrhoij_iatom%rhoijp(jrhoij,ispden) 992 ro_d=ro*pawtab(itypat)%dltij(klmn) 993 do ll=lmin,lmax,2 994 do ilm=ll**2+1,(ll+1)**2 995 isel=pawang%gntselect(ilm,klm) 996 if (isel>0) then 997 grhat_x=ro_d*pawtab(itypat)%qijl(ilm,klmn) 998 do mu=1,ngrad 999 grhat_tmp(mu,idiag)=grhat_tmp(mu,idiag)+grhat_x*prod(mu,ilm) 1000 end do 1001 end if 1002 end do 1003 end do 1004 jrhoij=jrhoij+pawrhoij_iatom%cplex 1005 end do 1006 1007 ! ---- Add additional (diagonal) terms for dynamical matrix 1008 ! ---- Terms including rhoij derivatives 1009 if (optgr2==1) then 1010 klmn1=1 1011 do klmn=1,lmn2_size 1012 klm =pawtab(itypat)%indklmn(1,klmn) 1013 lmin=pawtab(itypat)%indklmn(3,klmn) 1014 lmax=pawtab(itypat)%indklmn(4,klmn) 1015 dlt_tmp=pawtab(itypat)%dltij(klmn) 1016 do ll=lmin,lmax,2 1017 do ilm=ll**2+1,(ll+1)**2 1018 isel=pawang%gntselect(ilm,klm) 1019 if (isel>0) then 1020 ro_d= dlt_tmp*pawtab(itypat)%qijl(ilm,klmn) 1021 do mu=1,9 1022 mua=alpha(mu);mub=beta(mu) 1023 grhat_tmp(ishift_gr2+mu,idiag)=grhat_tmp(ishift_gr2+mu,idiag)& 1024 & +ro_d*pawrhoij_iatom%grhoij(ishift_grhoij+mua,klmn1,ispden)*prodp(mub+ishift_gr,ilm) 1025 end do 1026 end if 1027 end do 1028 end do 1029 klmn1=klmn1+pawrhoij_iatom%cplex 1030 end do ! klmn 1031 end if ! optgr2 1032 1033 ! ---- Add additional (diagonal) terms for elastic tensor 1034 ! ---- Terms including rhoij derivatives 1035 if (optstr2==1)then 1036 klmn1=1 1037 do klmn=1,lmn2_size 1038 klm =pawtab(itypat)%indklmn(1,klmn) 1039 lmin=pawtab(itypat)%indklmn(3,klmn) 1040 lmax=pawtab(itypat)%indklmn(4,klmn) 1041 dlt_tmp=pawtab(itypat)%dltij(klmn) 1042 do ll=lmin,lmax,2 1043 do ilm=ll**2+1,(ll+1)**2 1044 isel=pawang%gntselect(ilm,klm) 1045 if (isel>0) then 1046 ro_d=dlt_tmp*pawtab(itypat)%qijl(ilm,klmn) 1047 mu=1 1048 do mua=1,6 1049 do mub=1,6 1050 grhat_tmp(ishift_str2+mu,iatm)= grhat_tmp(ishift_str2+mu,iatm)& 1051 & +ro_d*pawrhoij_iatom%grhoij(mub,klmn1,ispden)*prodp(mua,ilm) 1052 grhat_tmp(ishift_str2+mu,iatm)= grhat_tmp(ishift_str2+mu,iatm)& 1053 & +ro_d*pawrhoij_iatom%grhoij(mua,klmn1,ispden)*prodp(mub,ilm) 1054 mu=mu+1 1055 end do 1056 ! INTERNAL STRAIN CONTRIBUTION 1057 do idir=1,3 1058 grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm) = grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm)& 1059 & +ro_d*pawrhoij_iatom%grhoij(ishift_grhoij+idir,klmn1,ispden)*prodp(mua,ilm) 1060 grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm) = grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm)& 1061 & +ro_d*pawrhoij_iatom%grhoij(mua,klmn1,ispden)*prodp(6+idir,ilm) 1062 end do 1063 end do 1064 end if 1065 end do 1066 end do 1067 klmn1=klmn1+pawrhoij_iatom%cplex 1068 end do 1069 end if ! optstr2 1070 1071 ! ---- Add off-diagonal additional contributions for second gradients 1072 if (optgr2==1.or.optstr2==1) then 1073 do jatm=1,natom 1074 jatom_tot=atindx1(jatm);jtypat=typat(jatom_tot) 1075 pawrhoij_jatom => pawrhoij_tot(jatom_tot) 1076 1077 ! ---- Dynamical matrix 1078 if (optgr2==1) then 1079 1080 ! Off-diagonal term including rhoij 1081 if (dyfr_cplex==1.or.cplex==1) then 1082 jrhoij=1 1083 do irhoij=1,pawrhoij_jatom%nrhoijsel 1084 klmn=pawrhoij_jatom%rhoijselect(irhoij) 1085 klm =pawtab(jtypat)%indklmn(1,klmn) 1086 lmin=pawtab(jtypat)%indklmn(3,klmn) 1087 lmax=pawtab(jtypat)%indklmn(4,klmn) 1088 ro =pawrhoij_jatom%rhoijp(jrhoij,ispden) 1089 ro_d=ro*pawtab(jtypat)%dltij(klmn) 1090 do ll=lmin,lmax,2 1091 do ilm=ll**2+1,(ll+1)**2 1092 isel=pawang%gntselect(ilm,klm) 1093 if (isel>0) then 1094 grhat_x=ro_d*pawtab(jtypat)%qijl(ilm,klmn) 1095 do mu=1,9 1096 grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm) & 1097 & +grhat_x*prod_nondiag(jatm)%value(ishift_gr2+mu,ilm) 1098 end do 1099 end if 1100 end do 1101 end do 1102 jrhoij=jrhoij+pawrhoij_jatom%cplex 1103 end do 1104 else 1105 jrhoij=1;mushift=ishift_gr2+9 1106 do irhoij=1,pawrhoij_jatom%nrhoijsel 1107 klmn=pawrhoij_jatom%rhoijselect(irhoij) 1108 klm =pawtab(jtypat)%indklmn(1,klmn) 1109 lmin=pawtab(jtypat)%indklmn(3,klmn) 1110 lmax=pawtab(jtypat)%indklmn(4,klmn) 1111 ro =pawrhoij_jatom%rhoijp(jrhoij,ispden) 1112 ro_d=ro*pawtab(jtypat)%dltij(klmn) 1113 do ll=lmin,lmax,2 1114 do ilm=ll**2+1,(ll+1)**2 1115 isel=pawang%gntselect(ilm,klm) 1116 if (isel>0) then 1117 grhat_x=ro_d*pawtab(jtypat)%qijl(ilm,klmn) 1118 do mu=1,9 1119 grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm)& 1120 & +grhat_x*prod_nondiag(jatm)%value(ishift_gr2+mu,ilm) 1121 grhat_tmp(mushift+mu,jatm)=grhat_tmp(mushift+mu,jatm)& 1122 & +grhat_x*prod_nondiag(jatm)%value(ishift_gr2+9+mu,ilm) 1123 end do 1124 end if 1125 end do 1126 end do 1127 jrhoij=jrhoij+pawrhoij_jatom%cplex 1128 end do 1129 end if 1130 1131 ! Off-diagonal term including rhoij derivative 1132 if (dyfr_cplex==1.or.cplex==1) then 1133 klmn1=1 1134 do klmn=1,pawrhoij_jatom%lmn2_size 1135 klm =pawtab(jtypat)%indklmn(1,klmn) 1136 lmin=pawtab(jtypat)%indklmn(3,klmn) 1137 lmax=pawtab(jtypat)%indklmn(4,klmn) 1138 dlt_tmp=pawtab(jtypat)%dltij(klmn) 1139 do ll=lmin,lmax,2 1140 do ilm=ll**2+1,(ll+1)**2 1141 isel=pawang%gntselect(ilm,klm) 1142 if (isel>0) then 1143 ro_d=dlt_tmp*pawtab(jtypat)%qijl(ilm,klmn) 1144 do mu=1,9 1145 mua=alpha(mu);mub=beta(mu) 1146 grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm) & 1147 & +ro_d*pawrhoij_jatom%grhoij(ishift_grhoij+mua,klmn1,ispden) & 1148 & *prodp_nondiag(jatm)%value(ishift2_gr+mub,ilm) 1149 end do 1150 end if 1151 end do 1152 end do 1153 klmn1=klmn1+pawrhoij_jatom%cplex 1154 end do ! klmn 1155 else ! ngradp_nondiag>=6 1156 klmn1=1;mushift=ishift_gr2+9 1157 do klmn=1,pawrhoij_jatom%lmn2_size 1158 klm =pawtab(jtypat)%indklmn(1,klmn) 1159 lmin=pawtab(jtypat)%indklmn(3,klmn) 1160 lmax=pawtab(jtypat)%indklmn(4,klmn) 1161 dlt_tmp=pawtab(jtypat)%dltij(klmn) 1162 do ll=lmin,lmax,2 1163 do ilm=ll**2+1,(ll+1)**2 1164 isel=pawang%gntselect(ilm,klm) 1165 if (isel>0) then 1166 ro_d=dlt_tmp*pawtab(jtypat)%qijl(ilm,klmn) 1167 do mu=1,9 1168 mua=alpha(mu);mub=beta(mu) 1169 grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm) & 1170 & +ro_d*pawrhoij_jatom%grhoij(ishift_grhoij+mua,klmn1,ispden) & 1171 & *prodp_nondiag(jatm)%value(ishift2_gr+mub,ilm) 1172 grhat_tmp(mushift+mu,jatm)=grhat_tmp(mushift+mu,jatm) & 1173 & +ro_d*pawrhoij_jatom%grhoij(ishift_grhoij+mua,klmn1,ispden) & 1174 & *prodp_nondiag(jatm)%value(ishift2_gr+3+mub,ilm) 1175 end do 1176 end if 1177 end do 1178 end do 1179 klmn1=klmn1+pawrhoij_jatom%cplex 1180 end do 1181 end if 1182 end if ! optgr2 1183 1184 ! ---- Elastic tensor 1185 if (optstr2==1)then 1186 1187 ! Off-diagonal term including rhoij 1188 jrhoij=1; 1189 do irhoij=1,pawrhoij_jatom%nrhoijsel 1190 klmn=pawrhoij_jatom%rhoijselect(irhoij) 1191 klm =pawtab(jtypat)%indklmn(1,klmn) 1192 lmin=pawtab(jtypat)%indklmn(3,klmn) 1193 lmax=pawtab(jtypat)%indklmn(4,klmn) 1194 ro =pawrhoij_jatom%rhoijp(jrhoij,ispden) 1195 ro_d=ro*pawtab(jtypat)%dltij(klmn) 1196 do ll=lmin,lmax,2 1197 do ilm=ll**2+1,(ll+1)**2 1198 isel=pawang%gntselect(ilm,klm) 1199 if (isel>0) then 1200 grhat_x=ro_d*pawtab(jtypat)%qijl(ilm,klmn) 1201 do mu=1,18 1202 grhat_tmp2(mu,jatm)=grhat_tmp2(mu,jatm) & 1203 & +grhat_x*prod_nondiag(jatm)%value(ishift_str2is+mu,ilm) 1204 end do 1205 end if 1206 end do 1207 end do 1208 jrhoij=jrhoij+pawrhoij_jatom%cplex 1209 end do 1210 ! Off-diagonal term including rhoij derivative 1211 klmn1=1 1212 do klmn=1,pawrhoij_jatom%lmn2_size 1213 klm =pawtab(jtypat)%indklmn(1,klmn) 1214 lmin=pawtab(jtypat)%indklmn(3,klmn) 1215 lmax=pawtab(jtypat)%indklmn(4,klmn) 1216 dlt_tmp=pawtab(jtypat)%dltij(klmn) 1217 do ll=lmin,lmax,2 1218 do ilm=ll**2+1,(ll+1)**2 1219 isel=pawang%gntselect(ilm,klm) 1220 if (isel>0) then 1221 ro_d=dlt_tmp*pawtab(jtypat)%qijl(ilm,klmn) 1222 mu=1 1223 do mua=1,6 1224 do idir=1,3 1225 grhat_tmp2((mua-1)*3+idir,jatm) = grhat_tmp2((mua-1)*3+idir,jatm) & 1226 & +ro_d*pawrhoij_jatom%grhoij(mua,klmn1,ispden) & 1227 & *prodp_nondiag(jatm)%value(idir,ilm) 1228 end do 1229 end do 1230 end if 1231 end do 1232 end do 1233 klmn1=klmn1+pawrhoij_jatom%cplex 1234 end do 1235 end if ! optstr2 1236 1237 end do ! jatm 1238 end if ! optgr2 or optstr2 1239 1240 ! ---------------------------------------------------------------- 1241 ! End of loop over spin components 1242 1243 end do ! ispden 1244 1245 ! Eventually free temporary space for g_l(r).Y_lm(r) factors 1246 if (pawfgrtab_iatom%gylm_allocated==2) then 1247 ABI_DEALLOCATE(pawfgrtab_iatom%gylm) 1248 ABI_ALLOCATE(pawfgrtab_iatom%gylm,(0,0)) 1249 pawfgrtab_iatom%gylm_allocated=0 1250 end if 1251 if (pawfgrtab_iatom%gylmgr_allocated==2) then 1252 ABI_DEALLOCATE(pawfgrtab_iatom%gylmgr) 1253 ABI_ALLOCATE(pawfgrtab_iatom%gylmgr,(0,0,0)) 1254 pawfgrtab_iatom%gylmgr_allocated=0 1255 end if 1256 if (pawfgrtab_iatom%gylmgr2_allocated==2) then 1257 ABI_DEALLOCATE(pawfgrtab_iatom%gylmgr2) 1258 ABI_ALLOCATE(pawfgrtab_iatom%gylmgr2,(0,0,0)) 1259 pawfgrtab_iatom%gylmgr2_allocated=0 1260 end if 1261 if (pawfgrtab_iatom%expiqr_allocated==2) then 1262 ABI_DEALLOCATE(pawfgrtab_iatom%expiqr) 1263 ABI_ALLOCATE(pawfgrtab_iatom%expiqr,(0,0)) 1264 pawfgrtab_iatom%expiqr_allocated=0 1265 end if 1266 1267 ! ---------------------------------------------------------------- 1268 ! Copy results in corresponding arrays 1269 1270 ! ==== Forces ==== 1271 ! Convert from cartesian to reduced coordinates 1272 if (optgr==1) then 1273 mushift=3*(iatm-1) 1274 tmp(1:3)=grhat_tmp(ishift_gr+1:ishift_gr+3,idiag) 1275 do mu=1,3 1276 hatgr(mu+mushift)=rprimd(1,mu)*tmp(1)+rprimd(2,mu)*tmp(2)+rprimd(3,mu)*tmp(3) 1277 end do 1278 end if 1279 1280 ! ==== Stresses ==== 1281 if (optstr==1) then 1282 hatstr(1:6)=hatstr(1:6)+grhat_tmp(ishift_str+1:ishift_str+6,idiag) 1283 end if 1284 1285 ! ==== Frozen wf part of dyn. matrix ==== 1286 if (optgr2==1) then 1287 do jatm=1,natom 1288 do mu=1,9 1289 mua=alpha(mu);mub=beta(mu) 1290 dyfr(1,mub,mua,jatm,iatm)=grhat_tmp(ishift_gr2+mu,jatm) 1291 end do 1292 if (dyfr_cplex==2.and.cplex==2) then 1293 mushift=ishift_gr2+9 1294 do mu=1,9 1295 mua=alpha(mu);mub=beta(mu) 1296 dyfr(2,mub,mua,jatm,iatm)=grhat_tmp(mushift+mu,jatm) 1297 end do 1298 end if 1299 end do 1300 end if 1301 1302 ! ==== Elastic tensor ==== 1303 if (optstr2==1) then 1304 eltfr(1:6,1:6)=eltfr(1:6,1:6)+ & 1305 & reshape(grhat_tmp(ishift_str2+1:ishift_str2+36,iatm),(/6,6/)) 1306 ! Convert internal Strain in reduced coordinates 1307 do mua = 1,6 1308 tmp(1:3)=grhat_tmp(ishift_str2is+(mua-1)*3+1:ishift_str2is+(mua-1)*3+3,iatm) 1309 do idir=1,3 1310 eltfr(6+(iatm-1)*3+idir,mua)=eltfr(6+(iatm-1)*3+idir,mua)+ & 1311 & (rprimd(1,idir)*tmp(1)+rprimd(2,idir)*tmp(2)+rprimd(3,idir)*tmp(3)) 1312 end do 1313 do jatm=1,natom 1314 tmp(1:3)=grhat_tmp2((mua-1)*3+1:(mua-1)*3+3,jatm) 1315 do idir=1,3 1316 eltfr(6+(iatm-1)*3+idir,mua)=eltfr(6+(iatm-1)*3+idir,mua)+ & 1317 & (rprimd(1,idir)*tmp(1)+rprimd(2,idir)*tmp(2)+rprimd(3,idir)*tmp(3)) 1318 end do 1319 end do 1320 end do 1321 end if 1322 1323 ! ---------------------------------------------------------------- 1324 ! End loops on types and atoms 1325 1326 ABI_DEALLOCATE(vloc) 1327 if (ngrad>0) then 1328 ABI_DEALLOCATE(prod) 1329 end if 1330 if (ngradp>0) then 1331 ABI_DEALLOCATE(prodp) 1332 end if 1333 if (optgr2==1.or.optstr2==1) then 1334 do jatm=1,natom 1335 ABI_DEALLOCATE(prod_nondiag(jatm)%value) 1336 ABI_DEALLOCATE(prodp_nondiag(jatm)%value) 1337 end do 1338 end if 1339 end do 1340 iatshft=iatshft+nattyp(itypat) 1341 end do 1342 1343 !Reduction in case of parallelisation over atoms 1344 if (paral_atom) then 1345 bufsiz=3*natom*optgr+6*optstr 1346 if (save_memory) bufsiz=bufsiz+9*dyfr_cplex*natom**2*optgr2+6*(6+3*natom)*optstr2 1347 if (bufsiz>0) then 1348 ABI_ALLOCATE(buf1,(bufsiz)) 1349 if (optgr==1) buf1(1:3*natom)=hatgr(1:3*natom) 1350 indx=optgr*3*natom 1351 if (optstr==1) buf1(indx+1:indx+6)=hatstr(1:6) 1352 indx=indx+optstr*6 1353 if (save_memory) then 1354 if (optgr2==1) then 1355 buf1(indx+1:indx+9*dyfr_cplex*natom**2)= & 1356 & reshape(dyfr,(/9*dyfr_cplex*natom**2/)) 1357 indx=indx+9*dyfr_cplex*natom**2 1358 end if 1359 if (optstr2==1) then 1360 buf1(indx+1:indx+6*(6+3*natom))= & 1361 & reshape(eltfr,(/6*(6+3*natom)/)) 1362 indx=indx+6*(6+3*natom) 1363 end if 1364 end if 1365 call xmpi_sum(buf1,my_comm_atom,ier) 1366 if (optgr==1) hatgr(1:3*natom)=buf1(1:3*natom) 1367 indx=optgr*3*natom 1368 if (optstr==1) hatstr(1:6)=buf1(indx+1:indx+6) 1369 indx=indx+optstr*6 1370 if (save_memory) then 1371 if (optgr2==1) then 1372 dyfr(1:dyfr_cplex,1:3,1:3,1:natom,1:natom)= & 1373 & reshape(buf1(indx+1:indx+9*dyfr_cplex*natom**2),(/dyfr_cplex,3,3,natom,natom/)) 1374 indx=indx+9*dyfr_cplex*natom**2 1375 end if 1376 if (optstr2==1) then 1377 eltfr(1:6+3*natom,1:6)= & 1378 & reshape(buf1(indx+1:indx+6*(6+3*natom)),(/6+3*natom,6/)) 1379 indx=indx+6*(6+3*natom) 1380 end if 1381 end if 1382 ABI_DEALLOCATE(buf1) 1383 end if 1384 end if 1385 1386 !Deallocate additional memory 1387 ABI_DEALLOCATE(grhat_tmp) 1388 if (optgr2==1.or.optstr2==1) then 1389 ABI_DEALLOCATE(mu4) 1390 ABI_DEALLOCATE(atindx) 1391 if (optgr2==1.or.optstr2==1) then 1392 ABI_DEALLOCATE(vpsp1_gr) 1393 end if 1394 if (optstr2==1) then 1395 ABI_DEALLOCATE(grhat_tmp2) 1396 ABI_DEALLOCATE(vpsp1_str) 1397 end if 1398 ABI_DATATYPE_DEALLOCATE(prod_nondiag) 1399 ABI_DATATYPE_DEALLOCATE(prodp_nondiag) 1400 if (.not.save_memory) then 1401 do jatom=1,size(pawfgrtab) 1402 pawfgrtab_jatom => pawfgrtab(jatom) 1403 if (pawfgrtab(jatom)%gylm_allocated==2) then 1404 ABI_DEALLOCATE(pawfgrtab(jatom)%gylm) 1405 ABI_ALLOCATE(pawfgrtab(jatom)%gylm,(0,0)) 1406 pawfgrtab(jatom)%gylm_allocated=0 1407 end if 1408 if (pawfgrtab(jatom)%gylmgr_allocated==2) then 1409 ABI_DEALLOCATE(pawfgrtab(jatom)%gylmgr) 1410 ABI_ALLOCATE(pawfgrtab(jatom)%gylmgr,(0,0,0)) 1411 pawfgrtab(jatom)%gylmgr_allocated=0 1412 end if 1413 if (pawfgrtab(jatom)%gylmgr2_allocated==2) then 1414 ABI_DEALLOCATE(pawfgrtab(jatom)%gylmgr2) 1415 ABI_ALLOCATE(pawfgrtab(jatom)%gylmgr2,(0,0,0)) 1416 pawfgrtab(jatom)%gylmgr2_allocated=0 1417 end if 1418 if (pawfgrtab(jatom)%expiqr_allocated==2) then 1419 ABI_DEALLOCATE(pawfgrtab(jatom)%expiqr) 1420 ABI_ALLOCATE(pawfgrtab(jatom)%expiqr,(0,0)) 1421 pawfgrtab(jatom)%expiqr_allocated=0 1422 end if 1423 end do 1424 end if 1425 if (paral_atom) then 1426 if ((.not.save_memory).and.paral_atom_pawfgrtab) then 1427 call pawfgrtab_free(pawfgrtab_tot) 1428 ABI_DATATYPE_DEALLOCATE(pawfgrtab_tot) 1429 end if 1430 if (paral_atom_pawrhoij) then 1431 call pawrhoij_free(pawrhoij_tot) 1432 ABI_DATATYPE_DEALLOCATE(pawrhoij_tot) 1433 end if 1434 end if 1435 end if 1436 1437 !---------------------------------------------------------------------- 1438 !Update non-local gardients 1439 1440 !===== Update forces ===== 1441 if (optgr==1) then 1442 grnl(1:3*natom)=grnl(1:3*natom)+hatgr(1:3*natom) 1443 ABI_DEALLOCATE(hatgr) 1444 end if 1445 1446 !===== Convert stresses (add diag and off-diag contributions) ===== 1447 if (optstr==1) then 1448 1449 ! Has to compute int[nhat*vtrial] 1450 hatstr_diag=zero 1451 if (nspden==1.or.dimvtrial==1) then 1452 do ic=1,nfft 1453 hatstr_diag=hatstr_diag+vtrial_(ic,1)*nhat(ic,1) 1454 end do 1455 else if (nspden==2) then 1456 do ic=1,nfft 1457 hatstr_diag=hatstr_diag+vtrial_(ic,1)*nhat(ic,2)+vtrial_(ic,2)*(nhat(ic,1)-nhat(ic,2)) 1458 end do 1459 else if (nspden==4) then 1460 do ic=1,nfft 1461 hatstr_diag=hatstr_diag+half*(vtrial_(ic,1)*(nhat(ic,1)+nhat(ic,4)) & 1462 & +vtrial_(ic,2)*(nhat(ic,1)-nhat(ic,4))) & 1463 & +vtrial_(ic,3)*nhat(ic,2)+vtrial_(ic,4)*nhat(ic,3) 1464 end do 1465 end if 1466 hatstr_diag=hatstr_diag*fact_ucvol 1467 if (paral_grid) then 1468 call xmpi_sum(hatstr_diag,my_comm_grid,ier) 1469 end if 1470 1471 ! Convert hat contribution 1472 1473 hatstr(1:3)=(hatstr(1:3)+hatstr_diag)/ucvol 1474 hatstr(4:6)= hatstr(4:6)/ucvol 1475 1476 ! Add to already computed NL contrib 1477 nlstr(1:6)=nlstr(1:6)+hatstr(1:6) 1478 1479 ! Apply symmetries 1480 call stresssym(gprimd,nsym,nlstr,symrec) 1481 end if 1482 1483 !===== Convert dynamical matrix (from cartesian to reduced coordinates) ===== 1484 if (optgr2==1) then 1485 do iatm=1,natom 1486 do jatm=1,natom 1487 do mua=1,3 1488 do mub=1,3 1489 work1(1,mua,mub)=dyfr(1,mub,mua,jatm,iatm)+dyfr(1,mua,mub,iatm,jatm) 1490 end do 1491 end do 1492 if (dyfr_cplex==2) then 1493 do mua=1,3 1494 do mub=1,3 1495 work1(2,mua,mub)=dyfr(2,mub,mua,jatm,iatm)-dyfr(2,mua,mub,iatm,jatm) 1496 end do 1497 end do 1498 end if 1499 do mu=1,3 1500 work2(:,:,mu)=rprimd(1,mu)*work1(:,:,1)+rprimd(2,mu)*work1(:,:,2)+rprimd(3,mu)*work1(:,:,3) 1501 end do 1502 do mub=1,3 1503 do mua=1,3 1504 dyfrnl(:,mua,mub,jatm,iatm)=dyfrnl(:,mua,mub,jatm,iatm) & ! Already contains NL projectors contribution 1505 & +rprimd(1,mua)*work2(:,1,mub) & 1506 & +rprimd(2,mua)*work2(:,2,mub) & 1507 & +rprimd(3,mua)*work2(:,3,mub) 1508 end do 1509 end do 1510 end do 1511 end do 1512 ABI_DEALLOCATE(dyfr) 1513 end if 1514 1515 !===== Update elastic tensor ===== 1516 if (optstr2==1) then 1517 eltfrnl(1:6+3*natom,1:6)=eltfrnl(1:6+3*natom,1:6)+eltfr(1:6+3*natom,1:6) 1518 ABI_DEALLOCATE(eltfr) 1519 end if 1520 1521 !---------------------------------------------------------------------- 1522 !End 1523 1524 !Destroy temporary space 1525 if (usexcnhat==0) then 1526 ABI_DEALLOCATE(vtrial_) 1527 end if 1528 1529 !Destroy atom tables used for parallelism 1530 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1531 if (paral_atom) then 1532 ABI_DEALLOCATE(atm_indx) 1533 end if 1534 1535 !Destroy FFT tables used for parallelism 1536 if ((optgr2==1.or.optstr2==1).and.(.not.present(comm_fft))) then 1537 call destroy_distribfft(my_distribfft) 1538 ABI_DATATYPE_DEALLOCATE(my_distribfft) 1539 end if 1540 1541 DBG_ENTER("COLL") 1542 1543 1544 end subroutine pawgrnl