TABLE OF CONTENTS
ABINIT/atm2fft [ Functions ]
NAME
atm2fft
FUNCTION
This routine sums atomic functions (density or potential) defined (in rec. space) on a radial grid to get global quantities on the fine FFT grid. It can also compute contribution to energy derivatives of these atomic functions. Possible options: optn=1: compute a sum of local atomic densities or contrib. to energy derivatives optv=1: compute a sum of local atomic potentials or contrib. to energy derivatives optatm =1: computes sum of atomic potentials/densities optgr =1: computes contribution of atomic pot./dens. to forces optstr =1: computes contribution of atomic pot./dens. to stress tensor optdyfr=1: computes contribution of atomic pot./dens. to frozen part of dyn. matrix
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (FJ, MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism gauss(2,ntypat)= params for gaussian atm density (optn2=3) for each atom type gmet(3,3)=reciprocal space metric gprimd(3,3)=reciprocal space dimensional primitive translations gsqcut=cutoff on |G|^2: see setup1 for definition (doubled sphere) kxc(2,nfft)=exchange and correlation kernel me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present) mgfft=maximum size of 1D FFTs comm_fft=--optional-- MPI communicator over FFT components mqgrid=number of grid pts in q array for f(q) spline. natom=number of atoms in unit cell. nattyp(ntypat)=number of 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 ntypat=number of types of atoms. optatm,optdyfr,optgr,optn,optn2,optstr,optv= (see NOTES below) paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present) psps <type(pseudopotential_type)>=variables related to pseudopotentials pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information qgrid(mqgrid)=q grid for spline from 0 to qmax. qprtrb(3)= integer wavevector of possible perturbing potential in basis of reciprocal lattice translations rhog(2,nfft)=electron density rho(G) in reciprocal space (used only if optv=1 and (optgr=1 or optstr=1 or optdyfr=1)) ucvol=unit cell volume usepaw= 0 for non paw calculation; =1 for paw calculation vspl(mqgrid,2,ntypat)=q^2 v(q) spline of an atomic potential (used only if optv=1) vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero, perturbing potential is added of the form V(G)=(vprtrb(1)+I*vprtrb(2))/2 at the values G=qprtrb and (vprtrb(1)-I*vprtrb(2))/2 at G=-qprtrb vg(2,nfft)= potential V(G) in reciprocal space (used only if optn=1 and (optgr=1 or optstr=1 or optdyfr=1)) vg1(2,nfft)= 1st-order potential V(G) in reciprocal space (used only if opteltfr==1) vg1_core(2,nfft)= 1-order potential V(G) in reciprocal space with only core contribution (used only if opteltfr==1)
OUTPUT
======= if optv==1 ======= ============================ --- if optatm==1 atmvloc(nfft)=sum of local atomic potentials in real space --- if optgr==1 grv(3,natom)=contribution of atomic potentials to forces --- if optstr==1 strv(6)=contribution of atomic potentials to stress tensor cart. coordinates, symmetric tensor, 6 comp. in order 11,22,33,32,31,21 --- if optdyfr==1 dyfrv(3,3,natom)=contribution of atomic potentials to frozen part of dyn. matrix ======= if optn==1 ======= ============================ --- if optatm==1 atmrho(nfft)=sum of atomic densities in real space --- if optgr==1 grn(3,natom)=contribution of atomic densities to forces --- if optstr==1 strn(6)=contribution of atomic densities to stress tensor cart. coordinates, symmetric tensor, 6 comp. in order 11,22,33,32,31,21 --- if optdyfr==1 dyfrn(3,3,natom)=contribution of atomic densities to frozen part of dyn. matrix --- if opteltfr==1 eltfrn(6+3*natom,6)=contribution of atomic density to frozen part of stress tensor
NOTES
Details on possible options: ============================ optv: controls the computation of a local potential as sum of atomic potentials Vloc(r)=Sum_R[V^AT(r-R)] or its contributions to energy derivatives, i.e. derivatives of Int[Vloc(r).rho(r).dr] rho(r) is stored in reciprocal space in array rhog() V^AT is stored in reciprocal space in array vspl (in practice vspl(q)=q^2.V^AT(q)) optn: controls the computation of a density as sum of atomic densities n(r)=Sum_R[n^AT(r-R)] or its contributions to energy derivatives, i.e. derivatives of Int[n(r).V(r).dr] V(r) is stored in reciprocal space in array vg() n^AT is stored in reciprocal space: if optn2=1: n^AT is the atomic PAW PS core density stored in array pawtab%tcorespl() 2: n^AT is the atomic PAW PS valence density stored in array pawtab%tvalespl() 3: n^AT is a gaussian density: n(g)=gauss(1,ityp)*exp[-(gauss(2,ityp)*G)^2] Note: optv and optn can be activated together Options controlling which contrib. to Etot derivatives are computed: optatm =1: computes Vloc(r) or n(r) as sum of atomic potentials/densities optgr =1: computes contribution of atomic Vloc(r) or n(r) to forces optstr =1: computes contribution of atomic Vloc(r) or n(r) to stress tensor optdyfr =1: computes contribution of atomic Vloc(r) or n(r) to fr part of dyn. matrix opteltfr=1: computes contribution of atomic Vloc(r) or n(r) to elastic tensor Note: optatm, optgr, optstr, optelfr and optdyfr can be activated together Typical uses: ============= Computation of: - local potential: optv=1, optatm=1 - contrib. of local potential to Etot derivatives: optv=1, rhog=total valence density optgr=1 or optstr=1 or optdyfr=1 - PS core density: optn=1, optn2=1, optatm=1 - contrib. of NLCC to Etot derivatives: optn=1, optn2=1, vg=XC potential optgr=1 or optstr=1 or optdyfr=1 - sum of atomic valence densities: optn=1, optn2=2 or 3, optatm=1 - correction of forces due to potential residual: optn=1, optn2=2 or 3, optgr=1 vg=potential residual etc...
PARENTS
dfpt_dyfro,dfpt_eltfrxc,extraprho,forces,fresidrsp,prcref,prcref_PMA respfn,setvtr,stress
CHILDREN
destroy_distribfft,fourdp,init_distribfft_seq,initmpi_seq set_mpi_enreg_fft,timab,unset_mpi_enreg_fft,wrtout,xmpi_sum,zerosym
SOURCE
145 #if defined HAVE_CONFIG_H 146 #include "config.h" 147 #endif 148 149 #include "abi_common.h" 150 151 subroutine atm2fft(atindx1,atmrho,atmvloc,dyfrn,dyfrv,eltfrn,gauss,gmet,gprimd,& 152 & grn,grv,gsqcut,mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,& 153 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,& 154 & psps,pawtab,ph1d,qgrid,qprtrb,rhog,strn,strv,ucvol,usepaw,vg,vg1,vg1_core,vprtrb,vspl,& 155 & is2_in,comm_fft,me_g0,paral_kgb,distribfft) ! optional arguments 156 157 use m_profiling_abi 158 use m_errors 159 use defs_basis 160 use defs_abitypes 161 use m_errors 162 use m_xmpi 163 164 use defs_datatypes,only : pseudopotential_type 165 use m_pawtab, only : pawtab_type 166 use m_distribfft, only : distribfft_type 167 use m_fft, only : zerosym 168 use m_mpinfo, only : set_mpi_enreg_fft, unset_mpi_enreg_fft 169 170 !This section has been created automatically by the script Abilint (TD). 171 !Do not modify the following lines by hand. 172 #undef ABI_FUNC 173 #define ABI_FUNC 'atm2fft' 174 use interfaces_14_hidewrite 175 use interfaces_18_timing 176 use interfaces_51_manage_mpi 177 use interfaces_53_ffts 178 !End of the abilint section 179 180 implicit none 181 182 !Arguments ------------------------------------ 183 !scalars 184 integer,intent(in) :: mgfft,mqgrid,natom,nfft,ntypat,optatm,optdyfr,opteltfr 185 integer,intent(in) :: optgr,optn,optn2,optstr,optv,usepaw 186 integer,optional,intent(in) :: is2_in,me_g0,comm_fft,paral_kgb 187 real(dp),intent(in) :: gsqcut,ucvol 188 type(pseudopotential_type),target,intent(in) :: psps 189 type(distribfft_type),optional,intent(in),target :: distribfft 190 !arrays 191 integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),qprtrb(3) 192 real(dp),intent(in) :: gauss(2,ntypat*(optn2/3)),gmet(3,3),gprimd(3,3) 193 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid) 194 real(dp),intent(in) :: rhog(2,nfft*optv*max(optgr,optstr,optdyfr,opteltfr)) 195 real(dp),intent(in) :: vg(2,nfft*optn*max(optgr,optstr,optdyfr,opteltfr)) 196 real(dp),intent(in) :: vg1(2,nfft*optn*opteltfr),vg1_core(2,nfft*optn*opteltfr) 197 real(dp),intent(in) :: vprtrb(2),vspl(mqgrid,2,ntypat*optv) 198 real(dp),intent(out) :: atmrho(nfft*optn) 199 real(dp),intent(inout) :: atmvloc(nfft*optv) 200 real(dp),intent(out) :: dyfrn(3,3,natom*optn*optdyfr),dyfrv(3,3,natom*optv*optdyfr) 201 real(dp),intent(out) :: eltfrn(6+3*natom,6) 202 real(dp),intent(inout) :: grn(3,natom*optn*optgr) 203 real(dp),intent(out) :: grv(3,natom*optv*optgr),strn(6*optn*optstr) 204 real(dp),intent(out) :: strv(6*optv*optstr) 205 type(pawtab_type),target,intent(in) :: pawtab(ntypat*usepaw) 206 207 !Local variables ------------------------------ 208 !scalars 209 integer,parameter :: im=2,re=1 210 integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ierr,ig1,ig1_,ig2,ig2_,ig3,ig3_,ii,is1,is2 211 integer :: itypat,jj,js,ka,kb,kd,kg,me_fft,my_comm_fft,ndir,n1,n2,n3,nproc_fft,paral_kgb_fft 212 integer :: shift1,shift2,shift3 213 logical :: have_g0 214 real(dp),parameter :: tolfix=1.0000001_dp 215 real(dp) :: aa,alf2pi2,bb,cc,cutoff,dbl_ig1,dbl_ig2,dbl_ig3,dd,dg1,dg2,d2g,diff 216 real(dp) :: dn_at,d2n_at,d2n_at2,dq,dq2div6,dqdiv6,dqm1,dv_at,ee,ff,gauss1,gauss2,gg,gmag,gsquar,n_at 217 real(dp) :: ph12i,ph12r,ph1i,ph1r,ph2i,ph2r,ph3i,ph3r,sfi,sfr,term,term1,term2,tmpni,tmpnr 218 real(dp) :: tmpvi,tmpvr,v_at,xnorm 219 character(len=500) :: message 220 type(distribfft_type),pointer :: my_distribfft 221 type(mpi_type) :: mpi_enreg_fft 222 !arrays 223 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 224 real(dp), ABI_CONTIGUOUS pointer :: tvalespl(:,:),tcorespl(:,:) 225 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 226 integer :: delta(6)=(/1,1,1,0,0,0/) 227 real(dp) :: dgm(3,3,6),d2gm(3,3,6,6),gcart(3),tsec(2) 228 real(dp),allocatable :: dyfrn_indx(:,:,:),dyfrv_indx(:,:,:),grn_indx(:,:) 229 real(dp),allocatable :: grv_indx(:,:),phim_igia(:),phre_igia(:),workn(:,:) 230 real(dp),allocatable :: workv(:,:) 231 232 ! ************************************************************************* 233 234 DBG_ENTER("COLL") 235 236 !Check optional arguments 237 if (present(comm_fft)) then 238 if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then 239 MSG_BUG(' Need paral_kgb and me_g0 with comm_fft !') 240 end if 241 end if 242 243 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 244 245 me_fft=ngfft(11) 246 nproc_fft=ngfft(10) 247 248 249 !Get the distrib associated with this fft_grid 250 if (present(distribfft)) then 251 my_distribfft => distribfft 252 else 253 ABI_DATATYPE_ALLOCATE(my_distribfft,) 254 call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp') 255 end if 256 if (n2==my_distribfft%n2_coarse) then 257 fftn2_distrib => my_distribfft%tab_fftdp2_distrib 258 ffti2_local => my_distribfft%tab_fftdp2_local 259 else if (n2 == my_distribfft%n2_fine) then 260 fftn2_distrib => my_distribfft%tab_fftdp2dg_distrib 261 ffti2_local => my_distribfft%tab_fftdp2dg_local 262 else 263 MSG_BUG("Unable to find an allocated distrib for this fft grid") 264 end if 265 266 if (present(is2_in)) then 267 if(is2_in<1.or.is2_in>6) then 268 MSG_BUG("is2_in must be between 1 and 6") 269 else 270 ndir = 1 271 end if 272 ndir = -1 273 end if 274 275 !Zero out arrays to permit accumulation over atom types 276 if (optv==1.and.optatm==1) then 277 ABI_ALLOCATE(workv,(2,nfft)) 278 workv(:,:)=zero 279 end if 280 if (optn==1.and.optatm==1) then 281 ABI_ALLOCATE(workn,(2,nfft)) 282 workn(:,:)=zero 283 end if 284 if (optv==1.and.optgr==1) then 285 ABI_ALLOCATE(grv_indx,(3,natom)) 286 grv_indx(:,:)=zero 287 end if 288 if (optn==1.and.optgr==1) then 289 ABI_ALLOCATE(grn_indx,(3,natom)) 290 grn_indx(:,:)=zero 291 end if 292 if (optv==1.and.optdyfr==1) then 293 ABI_ALLOCATE(dyfrv_indx,(3,3,natom)) 294 dyfrv_indx(:,:,:)=zero 295 end if 296 if (optn==1.and.optdyfr==1) then 297 ABI_ALLOCATE(dyfrn_indx,(3,3,natom)) 298 dyfrn_indx(:,:,:)=zero 299 end if 300 if (optv==1.and.optstr==1) strv(:)=zero 301 if (optn==1.and.optstr==1) strn(:)=zero 302 if (opteltfr==1) eltfrn(:,:) = zero 303 304 !Compute 1st and 2nd derivatives of metric tensor wrt all strain components 305 !and store for use in inner loop below for elastic tensor. 306 if (opteltfr==1) then 307 dgm(:,:,:)=zero 308 d2gm(:,:,:,:)=zero 309 ! Loop over 2nd strain index 310 do is2=1,6 311 kg=idx(2*is2-1);kd=idx(2*is2) 312 do jj = 1,3 313 dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj)) 314 end do 315 316 ! Loop over 1st strain index, upper triangle only 317 do is1=1,is2 318 ka=idx(2*is1-1);kb=idx(2*is1) 319 do jj = 1,3 320 if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 321 & +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj) 322 if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 323 & +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj) 324 if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 325 & +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj) 326 if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 327 & +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj) 328 end do 329 end do !is1 330 end do !is2 331 end if 332 333 334 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 335 dqm1=1.0_dp/dq 336 dqdiv6=dq/6.0_dp 337 dq2div6=dq**2/6.0_dp 338 cutoff=gsqcut*tolfix 339 id1=n1/2+2 340 id2=n2/2+2 341 id3=n3/2+2 342 343 ABI_ALLOCATE(phre_igia,(natom)) 344 ABI_ALLOCATE(phim_igia,(natom)) 345 346 ia1=1 347 do itypat=1,ntypat 348 ! ia1,ia2 sets range of loop over atoms: 349 ia2=ia1+nattyp(itypat)-1 350 ii=0 351 352 if (optn2==3)then 353 gauss1=gauss(1,itypat) 354 gauss2=gauss(2,itypat) 355 alf2pi2=(two_pi*gauss2)**2 356 end if 357 358 if (usepaw == 1) then 359 tcorespl => pawtab(itypat)%tcorespl 360 tvalespl => pawtab(itypat)%tvalespl 361 else 362 tcorespl => psps%nctab(itypat)%tcorespl 363 tvalespl => psps%nctab(itypat)%tvalespl 364 end if 365 366 do i3=1,n3 367 ig3=i3-(i3/id3)*n3-1 368 ig3_=ig3;if (ig3_==(n3/2+1)) ig3_=0 369 do i2=1,n2 370 ig2=i2-(i2/id2)*n2-1 371 ig2_=ig2;if (ig2_==(n2/2+1)) ig2_=0 372 if(fftn2_distrib(i2)==me_fft) then 373 do i1=1,n1 374 ig1=i1-(i1/id1)*n1-1 375 ig1_=ig1;if (ig1_==(n1/2+1)) ig1_=0 376 ii=ii+1 377 gsquar=gsq_atm(ig1,ig2,ig3) 378 379 ! Skip G**2 outside cutoff: 380 if (gsquar<=cutoff) then 381 382 gmag=sqrt(gsquar) 383 have_g0=(ig1==0.and.ig2==0.and.ig3==0) 384 385 jj=1+int(gmag*dqm1) 386 diff=gmag-qgrid(jj) 387 388 ! Compute structure factor for all atoms of given type: 389 do ia=ia1,ia2 390 shift1=1+n1+(ia-1)*(2*n1+1) 391 shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1) 392 shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1) 393 ph1r=ph1d(1,ig1+shift1);ph1i=ph1d(2,ig1+shift1) 394 ph2r=ph1d(1,ig2+shift2);ph2i=ph1d(2,ig2+shift2) 395 ph3r=ph1d(1,ig3+shift3);ph3i=ph1d(2,ig3+shift3) 396 ph12r=ph1r*ph2r-ph1i*ph2i 397 ph12i=ph1r*ph2i+ph1i*ph2r 398 phre_igia(ia)=ph12r*ph3r-ph12i*ph3i 399 phim_igia(ia)=ph12r*ph3i+ph12i*ph3r 400 end do 401 402 ! Assemble structure factors for this type of atom= sum[exp(-i.piG.R)] 403 if (optatm==1.or.optstr==1.or.opteltfr==1) then 404 sfr=zero;sfi=zero 405 do ia=ia1,ia2 406 sfr=sfr+phre_igia(ia) 407 sfi=sfi-phim_igia(ia) 408 end do 409 end if 410 411 ! Compute V^AT(G) and/or n^AT(G) for given type of atom 412 ! Evaluate spline fit: p. 86 Numerical Recipes, Press et al; 413 ! NOTE: error in book for sign of "aa" term in derivative; 414 ! ! also see splfit routine. 415 if (optv==1.or.optn2/=3) then 416 bb = diff*dqm1 417 aa = 1.0_dp-bb 418 cc = aa*(aa**2-1.0_dp)*dq2div6 419 dd = bb*(bb**2-1.0_dp)*dq2div6 420 end if 421 if (optv==1) then 422 if (have_g0) then 423 v_at=zero 424 else 425 v_at=(aa*vspl(jj,1,itypat)+bb*vspl(jj+1,1,itypat)+& 426 & cc*vspl(jj,2,itypat)+dd*vspl(jj+1,2,itypat)) & 427 & /gsquar 428 end if 429 end if 430 if (optn==1) then 431 if (optn2==1) then 432 n_at=(aa*tcorespl(jj,1)+bb*tcorespl(jj+1,1)+cc*tcorespl(jj,2)+dd*tcorespl(jj+1,2)) 433 else if (optn2==2) then 434 n_at=(aa*tvalespl(jj,1)+bb*tvalespl(jj+1,1)+cc*tvalespl(jj,2)+dd*tvalespl(jj+1,2)) 435 else if (optn2==3) then 436 n_at=gauss1*exp(-gsquar*alf2pi2) 437 else 438 n_at=zero 439 end if 440 end if 441 442 ! Compute sum of local atomic potentials or densities 443 ! --------------------------------------------------- 444 if(optatm==1) then 445 ! Accumulate V^AT(G)*SF(G) or n^AT(G)*SF(G) 446 if (optv==1) then 447 workv(re,ii)=workv(re,ii)+sfr*v_at 448 workv(im,ii)=workv(im,ii)+sfi*v_at 449 end if 450 if (optn==1) then 451 workn(re,ii)=workn(re,ii)+sfr*n_at 452 workn(im,ii)=workn(im,ii)+sfi*n_at 453 end if 454 455 ! Compute contrib. to forces and/or frozen part of dyn. matrix 456 ! ------------------------------------------------------------- 457 else if (optgr==1.or.optdyfr==1) then 458 dbl_ig1=dble(ig1_);dbl_ig2=dble(ig2_);dbl_ig3=dble(ig3_) 459 ! Compute (2Pi)*V^AT(G)*rho(G) or (2Pi)*n^AT(G)*V(G) 460 if (optv==1) then 461 tmpvr=(two_pi*v_at)*rhog(re,ii) 462 tmpvi=(two_pi*v_at)*rhog(im,ii) 463 end if 464 if (optn==1) then 465 tmpnr=(two_pi*n_at)*vg(re,ii) 466 tmpni=(two_pi*n_at)*vg(im,ii) 467 end if 468 ! === contrib. to forces 469 if (optgr==1) then 470 ! Accumulate -(2Pi.G)*V^AT(G)*rho(G)*SF(G) 471 ! or -(2Pi)*n^AT(G)*V(G)*SF(G) into forces 472 if (optv==1) then 473 do ia=ia1,ia2 474 term=tmpvi*phre_igia(ia)+tmpvr*phim_igia(ia) 475 grv_indx(1,ia)=grv_indx(1,ia)-dbl_ig1*term 476 grv_indx(2,ia)=grv_indx(2,ia)-dbl_ig2*term 477 grv_indx(3,ia)=grv_indx(3,ia)-dbl_ig3*term 478 end do 479 end if 480 if (optn==1) then 481 do ia=ia1,ia2 482 term=tmpni*phre_igia(ia)+tmpnr*phim_igia(ia) 483 grn_indx(1,ia)=grn_indx(1,ia)-dbl_ig1*term 484 grn_indx(2,ia)=grn_indx(2,ia)-dbl_ig2*term 485 grn_indx(3,ia)=grn_indx(3,ia)-dbl_ig3*term 486 end do 487 end if 488 end if 489 ! === contrib. to frozen part of dyn. matrix 490 if (optdyfr==1) then 491 ! Accumulate -(2Pi^2.Gi.Gj)*V^AT(G)*rho(G)*SF(G) 492 ! or -(2Pi^2.Gi.Gj)*n^AT(G)*V(G)*SF(G) into dyn. matrix 493 if (optv==1) then 494 do ia=ia1,ia2 495 term=two_pi*(tmpvr*phre_igia(ia)-tmpvi*phim_igia(ia)) 496 dyfrv_indx(1,1,ia)=dyfrv_indx(1,1,ia)-dbl_ig1*dbl_ig1*term 497 dyfrv_indx(1,2,ia)=dyfrv_indx(1,2,ia)-dbl_ig1*dbl_ig2*term 498 dyfrv_indx(1,3,ia)=dyfrv_indx(1,3,ia)-dbl_ig1*dbl_ig3*term 499 dyfrv_indx(2,2,ia)=dyfrv_indx(2,2,ia)-dbl_ig2*dbl_ig2*term 500 dyfrv_indx(2,3,ia)=dyfrv_indx(2,3,ia)-dbl_ig2*dbl_ig3*term 501 dyfrv_indx(3,3,ia)=dyfrv_indx(3,3,ia)-dbl_ig3*dbl_ig3*term 502 end do 503 end if 504 if (optn==1) then 505 do ia=ia1,ia2 506 term=two_pi*(tmpnr*phre_igia(ia)-tmpni*phim_igia(ia)) 507 dyfrn_indx(1,1,ia)=dyfrn_indx(1,1,ia)-dbl_ig1*dbl_ig1*term 508 dyfrn_indx(1,2,ia)=dyfrn_indx(1,2,ia)-dbl_ig1*dbl_ig2*term 509 dyfrn_indx(1,3,ia)=dyfrn_indx(1,3,ia)-dbl_ig1*dbl_ig3*term 510 dyfrn_indx(2,2,ia)=dyfrn_indx(2,2,ia)-dbl_ig2*dbl_ig2*term 511 dyfrn_indx(2,3,ia)=dyfrn_indx(2,3,ia)-dbl_ig2*dbl_ig3*term 512 dyfrn_indx(3,3,ia)=dyfrn_indx(3,3,ia)-dbl_ig3*dbl_ig3*term 513 end do 514 end if 515 end if 516 end if 517 518 ! Compute (dV^AT(q)/dq)/q and/or (dn^AT(q)/dq)/q 519 ! For stress tensor and/or elastic tensor 520 ! --------------------------------- 521 if (optstr==1.or.opteltfr==1) then 522 ! Note: correction of Numerical Recipes sign error before (3._dp*aa**2-1._dp) 523 ! ee*dqm1 + ff*dqdiv6 is the best estimate of dV(q)/dq from splines 524 if (optv==1) then 525 if (have_g0) then 526 dv_at=zero 527 else 528 ee=vspl(jj+1,1,itypat)-vspl(jj,1,itypat) 529 ff=(3._dp*bb**2-1._dp)*vspl(jj+1,2,itypat)& 530 & -(3._dp*aa**2-1._dp)*vspl(jj ,2,itypat) 531 dv_at=((ee*dqm1+ff*dqdiv6)/gmag-2.0_dp*v_at)/gsquar 532 end if 533 end if 534 if (optn==1) then 535 if (have_g0) then 536 if (optn2==1) then 537 if (usepaw ==1) then 538 dn_at=pawtab(itypat)%dncdq0 539 else 540 dn_at=psps%nctab(itypat)%dncdq0 541 end if 542 else if (optn2==2) then 543 if (usepaw == 1) then 544 dn_at=pawtab(itypat)%dnvdq0 545 else 546 dn_at=psps%nctab(itypat)%dnvdq0 547 end if 548 else if (optn2==3) then 549 dn_at=-two*gauss1*alf2pi2 550 end if 551 if (opteltfr==1) then 552 d2n_at = 0 553 end if 554 else 555 if (optn2==1) then 556 ee=tcorespl(jj+1,1)-tcorespl(jj,1) 557 ff=(3._dp*bb**2-1._dp)*tcorespl(jj+1,2) & 558 & -(3._dp*aa**2-1._dp)*tcorespl(jj,2) 559 ! Also get nc''(q) 560 if (opteltfr==1) then 561 gg=aa*tcorespl(jj,2)+bb*tcorespl(jj+1,2) 562 end if 563 else if (optn2==2) then 564 ee=tvalespl(jj+1,1)-tvalespl(jj,1) 565 ff=(3._dp*bb**2-1._dp)*tvalespl(jj+1,2) & 566 & -(3._dp*aa**2-1._dp)*tvalespl(jj,2) 567 ! Also get nc''(q) 568 if (opteltfr==1) then 569 gg=aa*tvalespl(jj,2)+bb*tvalespl(jj+1,2) 570 end if 571 else if (optn2==3) then 572 dn_at=-two*gauss1*alf2pi2*exp(-gsquar*alf2pi2) 573 else 574 end if 575 dn_at = (ee*dqm1+ff*dqdiv6)/gmag 576 if (opteltfr==1) then 577 d2n_at = (gg-dn_at)/gsquar 578 d2n_at2 = gg/gmag**3 579 end if 580 end if 581 end if 582 583 ! Compute G in cartesian coordinates 584 gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+& 585 & gprimd(1,3)*dble(ig3) 586 gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+& 587 & gprimd(2,3)*dble(ig3) 588 gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+& 589 & gprimd(3,3)*dble(ig3) 590 end if 591 592 ! Compute contrib. to stress tensor 593 ! --------------------------------- 594 if (optstr==1)then 595 ! Accumulate -dV^AT/dG*rho(G)*SF(G)*Gi.Gj/G 596 ! or -dn^AT/dG*V(G)*SF(G)*Gi.Gj/G 597 ! into stress tensor 598 if (optv==1) then 599 term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi) 600 strv(1)=strv(1)-term*(dv_at*gcart(1)*gcart(1)+v_at) 601 strv(2)=strv(2)-term*(dv_at*gcart(2)*gcart(2)+v_at) 602 strv(3)=strv(3)-term*(dv_at*gcart(3)*gcart(3)+v_at) 603 strv(4)=strv(4)-term*dv_at*gcart(3)*gcart(2) 604 strv(5)=strv(5)-term*dv_at*gcart(3)*gcart(1) 605 strv(6)=strv(6)-term*dv_at*gcart(2)*gcart(1) 606 end if 607 if (optn==1) then 608 term=(vg(re,ii)*sfr+vg(im,ii)*sfi)*dn_at 609 strn(1)=strn(1)-term*gcart(1)*gcart(1) 610 strn(2)=strn(2)-term*gcart(2)*gcart(2) 611 strn(3)=strn(3)-term*gcart(3)*gcart(3) 612 strn(4)=strn(4)-term*gcart(3)*gcart(2) 613 strn(5)=strn(5)-term*gcart(3)*gcart(1) 614 strn(6)=strn(6)-term*gcart(2)*gcart(1) 615 end if 616 end if 617 618 ! Compute contrib. to elastic tensor 619 ! --------------------------------- 620 if (opteltfr==1) then 621 dbl_ig1=dble(ig1_);dbl_ig2=dble(ig2_);dbl_ig3=dble(ig3_) 622 ! if (ig1==0 .and. ig2==0 .and. ig3==0) cycle 623 ! Compute G*dG/Deps_{\gamme\delta} 624 dg2=0.5_dp*dgsqds_atm(ig1,ig2,ig3,is2_in) 625 626 term = vg (re,ii)*sfr + vg (im,ii)*sfi 627 term1 = vg1(re,ii)*sfr + vg1(im,ii)*sfi 628 term2 = vg1_core(re,ii)*sfr + vg1_core(im,ii)*sfi 629 630 do is1=1,6 631 632 ! Compute G*dG/Deps_{\alpha\beta} 633 dg1=0.5_dp*dgsqds_atm(ig1,ig2,ig3,is1) 634 ! Compute G^2*d2G/Deps_{alphabetagammadelta} 635 d2g=(0.25_dp*d2gsqds_atm(ig1,ig2,ig3,is1,is2_in)) 636 637 eltfrn(is1,is2_in) = eltfrn(is1,is2_in) + (term*(d2n_at*dg1*dg2 + d2g*dn_at)) 638 eltfrn(is1,is2_in) = eltfrn(is1,is2_in) + 0.5*(term1*dn_at*dg1) 639 eltfrn(is2_in,is1) = eltfrn(is2_in,is1) + 0.5*(term1*dn_at*dg1) 640 641 if(is2_in<=3)then 642 eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - term*dn_at*dg1 643 end if 644 if(is1<=3)then 645 eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - term*dn_at*dg2 646 end if 647 if(is2_in<=3.and.is1<=3)then 648 eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - (term1-term)*n_at 649 end if 650 end do 651 652 ! internal strain 653 do ia=ia1,ia2 654 js=7+3*(ia-1) 655 ! Compute -2pi*G*i*vxcis_core(G)*nat(G)*(exp(-iGr)) 656 term=(vg1_core(im,ii)*phre_igia(ia)+vg1_core(re,ii)*phim_igia(ia))*n_at*two_pi 657 eltfrn(js ,is2_in) = eltfrn(js ,is2_in) - dbl_ig1*term 658 eltfrn(js+1,is2_in) = eltfrn(js+1,is2_in) - dbl_ig2*term 659 eltfrn(js+2,is2_in) = eltfrn(js+2,is2_in) - dbl_ig3*term 660 661 ! Compute -2pi*G*i*vxc(G)*(dnat*dG/deps-delta*nat(G))*(exp(-iGr)) 662 term=(vg(im,ii)*phre_igia(ia)+vg(re,ii)*phim_igia(ia))*& 663 & (dn_at*dg2-delta(is2_in)*n_at)*two_pi 664 eltfrn(js ,is2_in) = eltfrn(js ,is2_in) - dbl_ig1*term 665 eltfrn(js+1,is2_in) = eltfrn(js+1,is2_in) - dbl_ig2*term 666 eltfrn(js+2,is2_in) = eltfrn(js+2,is2_in) - dbl_ig3*term 667 end do 668 end if 669 ! End skip G**2 outside cutoff: 670 end if 671 ! End loop on n1, n2, n3 672 end do 673 end if ! this plane is for me_fft 674 end do 675 end do 676 677 ! Symmetrize the dynamical matrix with respect to indices 678 if (optdyfr==1) then 679 if (optv==1) then 680 do ia=ia1,ia2 681 dyfrv_indx(2,1,ia)=dyfrv_indx(1,2,ia) 682 dyfrv_indx(3,1,ia)=dyfrv_indx(1,3,ia) 683 dyfrv_indx(3,2,ia)=dyfrv_indx(2,3,ia) 684 end do 685 end if 686 if (optn==1) then 687 do ia=ia1,ia2 688 dyfrn_indx(2,1,ia)=dyfrn_indx(1,2,ia) 689 dyfrn_indx(3,1,ia)=dyfrn_indx(1,3,ia) 690 dyfrn_indx(3,2,ia)=dyfrn_indx(2,3,ia) 691 end do 692 end if 693 end if 694 695 ia1=ia2+1 696 697 ! End loop on type of atoms 698 end do 699 700 ABI_DEALLOCATE(phre_igia) 701 ABI_DEALLOCATE(phim_igia) 702 703 !Get local potential or density back to real space 704 if(optatm==1)then 705 ! Allow for the addition of a perturbing potential 706 if ((optv==1).and.(vprtrb(1)**2+vprtrb(2)**2) > 1.d-30) then 707 ! Find the linear indices which correspond with the input wavevector qprtrb 708 ! The double modulus handles both i>=n and i<0, mapping into [0,n-1]; 709 ! then add 1 to get range [1,n] for each 710 i3=1+mod(n3+mod(qprtrb(3),n3),n3) 711 i2=1+mod(n2+mod(qprtrb(2),n2),n2) 712 i1=1+mod(n1+mod(qprtrb(1),n1),n1) 713 ! Compute the linear index in the 3 dimensional array 714 ii=i1+n1*((ffti2_local(i2)-1)+(n2/nproc_fft)*(i3-1)) 715 ! Add in the perturbation at G=qprtrb 716 workv(re,ii)=workv(re,ii)+0.5_dp*vprtrb(1) 717 workv(im,ii)=workv(im,ii)+0.5_dp*vprtrb(2) 718 ! Same thing for G=-qprtrb 719 i3=1+mod(n3+mod(-qprtrb(3),n3),n3) 720 i2=1+mod(n2+mod(-qprtrb(2),n2),n2) 721 i1=1+mod(n1+mod(-qprtrb(1),n1),n1) 722 ! ii=i1+n1*((i2-1)+n2*(i3-1)) 723 workv(re,ii)=workv(re,ii)+0.5_dp*vprtrb(1) 724 workv(im,ii)=workv(im,ii)-0.5_dp*vprtrb(2) 725 write(message, '(a,1p,2e12.4,a,0p,3i4,a)' )& 726 & ' atm2fft: perturbation of vprtrb=', vprtrb,& 727 & ' and q=',qprtrb,' has been added' 728 call wrtout(std_out,message,'COLL') 729 end if 730 731 if (optv==1.or.optn==1) then 732 ! Create fake mpi_enreg to wrap fourdp 733 call initmpi_seq(mpi_enreg_fft) 734 ABI_DATATYPE_DEALLOCATE(mpi_enreg_fft%distribfft) 735 if (present(comm_fft)) then 736 call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb) 737 my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb 738 else 739 my_comm_fft=xmpi_comm_self;paral_kgb_fft=0; 740 mpi_enreg_fft%distribfft => my_distribfft 741 end if 742 ! Non-symetrized non-zero elements have to be nullified 743 ! Transform back to real space; divide by unit cell volume 744 xnorm=one/ucvol 745 if (optv==1) then 746 call zerosym(workv,2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft) 747 call fourdp(1,workv,atmvloc,1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0) 748 atmvloc(:)=atmvloc(:)*xnorm 749 ABI_DEALLOCATE(workv) 750 end if 751 if (optn==1) then 752 call zerosym(workn,2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft) 753 call fourdp(1,workn,atmrho,1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0) 754 atmrho(:)=atmrho(:)*xnorm 755 ABI_DEALLOCATE(workn) 756 end if 757 ! Destroy fake mpi_enreg 758 call unset_mpi_enreg_fft(mpi_enreg_fft) 759 end if 760 761 end if 762 763 !Additional treatment in case of parallelization 764 if (present(comm_fft)) then 765 if ((xmpi_comm_size(comm_fft)>1).and.& 766 & (optgr==1.or.optstr==1.or.optdyfr==1.or.opteltfr==1)) then 767 call timab(48,1,tsec) 768 if (optv==1) then 769 if (optgr==1)then 770 call xmpi_sum(grv_indx,comm_fft,ierr) 771 end if 772 if (optstr==1)then 773 call xmpi_sum(strv,comm_fft,ierr) 774 end if 775 if (optdyfr==1)then 776 call xmpi_sum(dyfrv_indx,comm_fft,ierr) 777 end if 778 end if 779 if (optn==1) then 780 if (optgr==1)then 781 call xmpi_sum(grn_indx,comm_fft,ierr) 782 end if 783 if (optstr==1)then 784 call xmpi_sum(strn,comm_fft,ierr) 785 end if 786 if (optdyfr==1)then 787 call xmpi_sum(dyfrn_indx,comm_fft,ierr) 788 end if 789 end if 790 call timab(48,2,tsec) 791 end if 792 end if 793 794 !Forces: re-order atoms 795 if (optgr==1) then 796 if (optv==1) then 797 do ia=1,natom 798 grv(1:3,atindx1(ia))=grv_indx(1:3,ia) 799 end do 800 ABI_DEALLOCATE(grv_indx) 801 end if 802 if (optn==1) then 803 do ia=1,natom 804 grn(1:3,atindx1(ia))=grn_indx(1:3,ia) 805 end do 806 ABI_DEALLOCATE(grn_indx) 807 end if 808 end if 809 810 !Elastic tensor: Fill in lower triangle 811 ! if (opteltfr==1) then 812 ! do is2=2,6 813 ! do is1=1,is2-1 814 ! eltfrn(is2,is1)=eltfrn(is1,is2) 815 ! end do 816 ! end do 817 ! end if 818 819 !Normalize stress tensor: 820 if (optstr==1) then 821 if (optv==1) then 822 strv(:)=strv(:)/ucvol 823 end if 824 if (optn==1) strn(:)=strn(:)/ucvol 825 end if 826 827 !Dynamical matrix: re-order atoms 828 if (optdyfr==1) then 829 if (optv==1) then 830 do ia=1,natom 831 dyfrv(1:3,1:3,atindx1(ia))=dyfrv_indx(1:3,1:3,ia) 832 end do 833 ABI_DEALLOCATE(dyfrv_indx) 834 end if 835 if (optn==1) then 836 do ia=1,natom 837 dyfrn(1:3,1:3,atindx1(ia))=dyfrn_indx(1:3,1:3,ia) 838 end do 839 ABI_DEALLOCATE(dyfrn_indx) 840 end if 841 end if 842 843 if (.not.present(distribfft)) then 844 call destroy_distribfft(my_distribfft) 845 ABI_DATATYPE_DEALLOCATE(my_distribfft) 846 end if 847 848 DBG_EXIT("COLL") 849 850 contains 851 852 function gsq_atm(i1,i2,i3) 853 854 855 !This section has been created automatically by the script Abilint (TD). 856 !Do not modify the following lines by hand. 857 #undef ABI_FUNC 858 #define ABI_FUNC 'gsq_atm' 859 !End of the abilint section 860 861 real(dp) :: gsq_atm 862 integer,intent(in) :: i1,i2,i3 863 gsq_atm=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+dble(i3*i3)*gmet(3,3) & 864 & +two*(dble(i1*i2)*gmet(1,2)+dble(i2*i3)*gmet(2,3)+dble(i3*i1)*gmet(3,1)) 865 end function gsq_atm 866 867 function dgsqds_atm(i1,i2,i3,is) 868 !Define dG^2/ds based on G space metric derivative 869 870 !This section has been created automatically by the script Abilint (TD). 871 !Do not modify the following lines by hand. 872 #undef ABI_FUNC 873 #define ABI_FUNC 'dgsqds_atm' 874 !End of the abilint section 875 876 real(dp) :: dgsqds_atm 877 integer,intent(in) :: i1,i2,i3,is 878 dgsqds_atm=dble(i1*i1)*dgm(1,1,is)+dble(i2*i2)*dgm(2,2,is)+& 879 & dble(i3*i3)*dgm(3,3,is)+& 880 & dble(i1*i2)*(dgm(1,2,is)+dgm(2,1,is))+& 881 & dble(i1*i3)*(dgm(1,3,is)+dgm(3,1,is))+& 882 & dble(i2*i3)*(dgm(2,3,is)+dgm(3,2,is)) 883 end function dgsqds_atm 884 885 function d2gsqds_atm(i1,i2,i3,is1,is2) 886 ! Define 2dG^2/ds1ds2 based on G space metric derivative 887 888 !This section has been created automatically by the script Abilint (TD). 889 !Do not modify the following lines by hand. 890 #undef ABI_FUNC 891 #define ABI_FUNC 'd2gsqds_atm' 892 !End of the abilint section 893 894 real(dp) :: d2gsqds_atm 895 integer,intent(in) :: i1,i2,i3,is1,is2 896 d2gsqds_atm=dble(i1*i1)*d2gm(1,1,is1,is2)+& 897 & dble(i2*i2)*d2gm(2,2,is1,is2)+dble(i3*i3)*d2gm(3,3,is1,is2)+& 898 & dble(i1*i2)*(d2gm(1,2,is1,is2)+d2gm(2,1,is1,is2))+& 899 & dble(i1*i3)*(d2gm(1,3,is1,is2)+d2gm(3,1,is1,is2))+& 900 & dble(i2*i3)*(d2gm(2,3,is1,is2)+d2gm(3,2,is1,is2)) 901 end function d2gsqds_atm 902 903 end subroutine atm2fft