TABLE OF CONTENTS
ABINIT/atm2fft [ Functions ]
NAME
atm2fft
FUNCTION
This routine sums atomic functions (density, kinetic 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 [kinetic] 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
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] 4: n^AT is the atomic PAW PS core kinetic density stored in array pawtab%ttaucorespl() 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...
SOURCE
181 subroutine atm2fft(atindx1,atmrho,atmvloc,dyfrn,dyfrv,eltfrn,gauss,gmet,gprimd,& 182 & grn,grv,gsqcut,mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,& 183 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,& 184 & psps,pawtab,ph1d,qgrid,qprtrb,rcut,rhog,rprimd,strn,strv,ucvol,usepaw,vg,vg1,vg1_core,vprtrb,vspl,& 185 & is2_in,comm_fft,me_g0,paral_kgb,distribfft) ! optional arguments 186 187 !Arguments ------------------------------------ 188 !scalars 189 integer,intent(in) :: mgfft,mqgrid,natom,nfft,ntypat,optatm,optdyfr,opteltfr 190 integer,intent(in) :: optgr,optn,optn2,optstr,optv,usepaw 191 integer,optional,intent(in) :: is2_in,me_g0,comm_fft,paral_kgb 192 real(dp),intent(in) :: rcut,gsqcut,ucvol 193 type(pseudopotential_type),target,intent(in) :: psps 194 type(distribfft_type),optional,intent(in),target :: distribfft 195 !arrays 196 integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),qprtrb(3) 197 real(dp),intent(in) :: gauss(2,ntypat*(optn2/3)),gmet(3,3),gprimd(3,3) 198 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid) 199 real(dp),intent(in) :: rhog(2,nfft*optv*max(optgr,optstr,optdyfr,opteltfr)) 200 real(dp),intent(inout) :: rprimd(3,3) 201 real(dp),intent(in) :: vg(2,nfft*optn*max(optgr,optstr,optdyfr,opteltfr)) 202 real(dp),intent(in) :: vg1(2,nfft*optn*opteltfr),vg1_core(2,nfft*optn*opteltfr) 203 real(dp),intent(in) :: vprtrb(2),vspl(mqgrid,2,ntypat*optv) 204 real(dp),intent(out) :: atmrho(nfft*optn) 205 real(dp),intent(inout) :: atmvloc(nfft*optv) 206 real(dp),intent(out) :: dyfrn(3,3,natom*optn*optdyfr),dyfrv(3,3,natom*optv*optdyfr) 207 real(dp),intent(out) :: eltfrn(6+3*natom,6) 208 real(dp),intent(inout) :: grn(3,natom*optn*optgr) 209 real(dp),intent(out) :: grv(3,natom*optv*optgr),strn(6*optn*optstr) 210 real(dp),intent(out) :: strv(6*optv*optstr) 211 type(pawtab_type),target,intent(in) :: pawtab(ntypat*usepaw) 212 213 !Local variables ------------------------------ 214 !scalars 215 integer,parameter :: im=2,re=1,icutcoul=3 216 integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ierr,ig1,ig1_,ig2,ig2_,ig3,ig3_,ii,is1,is2 217 integer :: itypat,jj,js,ka,kb,kd,kg,me_fft,my_comm_fft,ndir,n1,n2,n3,nproc_fft,paral_kgb_fft 218 integer :: shift1,shift2,shift3 219 logical :: have_g0 220 #ifdef FC_NVHPC 221 !Silly trick to prevent NVHPC optimization issue 222 logical :: nothing=.false. 223 #endif 224 real(dp),parameter :: tolfix=1.0000001_dp, vcutgeo(3)=zero 225 real(dp) :: aa,alf2pi2,bb,cc,cutoff,dbl_ig1,dbl_ig2,dbl_ig3,dd,dg1,dg2,d2g,diff 226 real(dp) :: dn_at,d2n_at,d2n_at2,dq,dq2div6,dqdiv6,dqm1,dv_at,ee,ff,gauss1,gauss2,gg,gmag,gsquar,n_at 227 real(dp) :: ph12i,ph12r,ph1i,ph1r,ph2i,ph2r,ph3i,ph3r,sfi,sfr,term,term1,term2,tmpni,tmpnr 228 real(dp) :: tmpvi,tmpvr,v_at,xnorm 229 character(len=500) :: message 230 type(distribfft_type),pointer :: my_distribfft 231 type(mpi_type) :: mpi_enreg_fft 232 !arrays 233 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 234 real(dp), ABI_CONTIGUOUS pointer :: tvalespl(:,:),tcorespl(:,:),ttaucorespl(:,:) 235 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 236 integer :: delta(6)=(/1,1,1,0,0,0/) 237 real(dp) :: dgm(3,3,6),d2gm(3,3,6,6),gcart(3),tsec(2) 238 real(dp),allocatable :: dyfrn_indx(:,:,:),dyfrv_indx(:,:,:),grn_indx(:,:) 239 real(dp),allocatable :: grv_indx(:,:),phim_igia(:),phre_igia(:),workn(:,:) 240 real(dp),allocatable :: gcutoff(:) 241 real(dp),allocatable :: workv(:,:) 242 243 ! ************************************************************************* 244 245 DBG_ENTER("COLL") 246 247 !Check optional arguments 248 if (present(comm_fft)) then 249 if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then 250 ABI_BUG(' Need paral_kgb and me_g0 with comm_fft !') 251 end if 252 end if 253 254 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 255 me_fft=ngfft(11) 256 nproc_fft=ngfft(10) 257 258 !Get the distrib associated with this fft_grid 259 if (present(distribfft)) then 260 my_distribfft => distribfft 261 else 262 ABI_MALLOC(my_distribfft,) 263 call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp') 264 end if 265 if (n2==my_distribfft%n2_coarse) then 266 fftn2_distrib => my_distribfft%tab_fftdp2_distrib 267 ffti2_local => my_distribfft%tab_fftdp2_local 268 else if (n2 == my_distribfft%n2_fine) then 269 fftn2_distrib => my_distribfft%tab_fftdp2dg_distrib 270 ffti2_local => my_distribfft%tab_fftdp2dg_local 271 else 272 ABI_BUG("Unable to find an allocated distrib for this fft grid") 273 end if 274 275 if (present(is2_in)) then 276 if(is2_in<1.or.is2_in>6) then 277 ABI_BUG("is2_in must be between 1 and 6") 278 else 279 ndir = 1 280 end if 281 ndir = -1 282 end if 283 284 !Zero out arrays to permit accumulation over atom types 285 if (optv==1.and.optatm==1) then 286 ABI_MALLOC(workv,(2,nfft)) 287 workv(:,:)=zero 288 end if 289 if (optn==1.and.optatm==1) then 290 ABI_MALLOC(workn,(2,nfft)) 291 workn(:,:)=zero 292 end if 293 if (optv==1.and.optgr==1) then 294 ABI_MALLOC(grv_indx,(3,natom)) 295 grv_indx(:,:)=zero 296 end if 297 if (optn==1.and.optgr==1) then 298 ABI_MALLOC(grn_indx,(3,natom)) 299 grn_indx(:,:)=zero 300 end if 301 if (optv==1.and.optdyfr==1) then 302 ABI_MALLOC(dyfrv_indx,(3,3,natom)) 303 dyfrv_indx(:,:,:)=zero 304 end if 305 if (optn==1.and.optdyfr==1) then 306 ABI_MALLOC(dyfrn_indx,(3,3,natom)) 307 dyfrn_indx(:,:,:)=zero 308 end if 309 if (optv==1.and.optstr==1) strv(:)=zero 310 if (optn==1.and.optstr==1) strn(:)=zero 311 if (opteltfr==1) eltfrn(:,:) = zero 312 313 !Compute 1st and 2nd derivatives of metric tensor wrt all strain components 314 !and store for use in inner loop below for elastic tensor. 315 if (opteltfr==1) then 316 dgm(:,:,:)=zero 317 d2gm(:,:,:,:)=zero 318 ! Loop over 2nd strain index 319 do is2=1,6 320 kg=idx(2*is2-1);kd=idx(2*is2) 321 do jj = 1,3 322 dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj)) 323 end do 324 325 ! Loop over 1st strain index, upper triangle only 326 do is1=1,is2 327 ka=idx(2*is1-1);kb=idx(2*is1) 328 do jj = 1,3 329 if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 330 & +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj) 331 if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 332 & +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj) 333 if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 334 & +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj) 335 if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)& 336 & +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj) 337 end do 338 end do !is1 339 end do !is2 340 end if 341 342 343 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 344 dqm1=1.0_dp/dq 345 dqdiv6=dq/6.0_dp 346 dq2div6=dq**2/6.0_dp 347 cutoff=gsqcut*tolfix 348 id1=n1/2+2 349 id2=n2/2+2 350 id3=n3/2+2 351 352 ABI_MALLOC(phre_igia,(natom)) 353 ABI_MALLOC(phim_igia,(natom)) 354 355 !Initialize Gcut-off array from m_termcutoff 356 !ABI_MALLOC(gcutoff,(nfft)) 357 call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,1,rcut,rprimd,vcutgeo) 358 359 ia1=1 360 do itypat=1,ntypat 361 ! ia1,ia2 sets range of loop over atoms: 362 ia2=ia1+nattyp(itypat)-1 363 ii=0 364 365 if (optn2==3)then 366 gauss1=gauss(1,itypat) 367 gauss2=gauss(2,itypat) 368 alf2pi2=(two_pi*gauss2)**2 369 end if 370 371 if (usepaw == 1) then 372 tcorespl => pawtab(itypat)%tcorespl 373 tvalespl => pawtab(itypat)%tvalespl 374 ttaucorespl => pawtab(itypat)%tcoretauspl 375 else 376 tcorespl => psps%nctab(itypat)%tcorespl 377 tvalespl => psps%nctab(itypat)%tvalespl 378 ttaucorespl => null() 379 end if 380 381 do i3=1,n3 382 ig3=i3-(i3/id3)*n3-1 383 ig3_=ig3;if (ig3_==(n3/2+1)) ig3_=0 384 do i2=1,n2 385 ig2=i2-(i2/id2)*n2-1 386 ig2_=ig2;if (ig2_==(n2/2+1)) ig2_=0 387 if(fftn2_distrib(i2)==me_fft) then 388 do i1=1,n1 389 ig1=i1-(i1/id1)*n1-1 390 ig1_=ig1;if (ig1_==(n1/2+1)) ig1_=0 391 ii=ii+1 392 gsquar=gsq_atm(ig1,ig2,ig3) 393 394 ! Skip G**2 outside cutoff: 395 if (gsquar<=cutoff) then 396 397 gmag=sqrt(gsquar) 398 have_g0=(ig1==0.and.ig2==0.and.ig3==0) 399 400 jj=1+int(gmag*dqm1) 401 diff=gmag-qgrid(jj) 402 403 ! Compute structure factor for all atoms of given type: 404 do ia=ia1,ia2 405 #ifdef FC_NVHPC 406 !Silly trick to prevent NVHPC optimization issue 407 if(nothing) write(100,*) shift1,shift2,shift3 408 #endif 409 shift1=1+n1+(ia-1)*(2*n1+1) 410 shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1) 411 shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1) 412 ph1r=ph1d(1,ig1+shift1);ph1i=ph1d(2,ig1+shift1) 413 ph2r=ph1d(1,ig2+shift2);ph2i=ph1d(2,ig2+shift2) 414 ph3r=ph1d(1,ig3+shift3);ph3i=ph1d(2,ig3+shift3) 415 ph12r=ph1r*ph2r-ph1i*ph2i 416 ph12i=ph1r*ph2i+ph1i*ph2r 417 phre_igia(ia)=ph12r*ph3r-ph12i*ph3i 418 phim_igia(ia)=ph12r*ph3i+ph12i*ph3r 419 end do 420 421 ! Assemble structure factors for this type of atom= sum[exp(-i.piG.R)] 422 if (optatm==1.or.optstr==1.or.opteltfr==1) then 423 sfr=zero;sfi=zero 424 do ia=ia1,ia2 425 sfr=sfr+phre_igia(ia) 426 sfi=sfi-phim_igia(ia) 427 end do 428 end if 429 430 ! Compute V^AT(G) and/or n^AT(G) for given type of atom 431 ! Evaluate spline fit: p. 86 Numerical Recipes, Press et al; 432 ! NOTE: error in book for sign of "aa" term in derivative; 433 ! ! also see splfit routine. 434 if (optv==1.or.optn2/=3) then 435 bb = diff*dqm1 436 aa = 1.0_dp-bb 437 cc = aa*(aa**2-1.0_dp)*dq2div6 438 dd = bb*(bb**2-1.0_dp)*dq2div6 439 end if 440 if (optv==1) then 441 if (have_g0) then 442 v_at=zero 443 else 444 v_at=(aa*vspl(jj,1,itypat)+bb*vspl(jj+1,1,itypat)+& 445 & cc*vspl(jj,2,itypat)+dd*vspl(jj+1,2,itypat)) & 446 & /gsquar * gcutoff(ii) 447 end if 448 end if 449 if (optn==1) then 450 if (optn2==1) then 451 n_at=(aa*tcorespl(jj,1)+bb*tcorespl(jj+1,1)+cc*tcorespl(jj,2)+dd*tcorespl(jj+1,2)) 452 else if (optn2==2) then 453 n_at=(aa*tvalespl(jj,1)+bb*tvalespl(jj+1,1)+cc*tvalespl(jj,2)+dd*tvalespl(jj+1,2)) 454 else if (optn2==3) then 455 n_at=gauss1*exp(-gsquar*alf2pi2) 456 else if (optn2==4) then 457 n_at=(aa*ttaucorespl(jj,1)+bb*ttaucorespl(jj+1,1)+cc*ttaucorespl(jj,2)+dd*ttaucorespl(jj+1,2)) 458 else 459 n_at=zero 460 end if 461 end if 462 463 ! Compute sum of local atomic potentials or densities 464 ! --------------------------------------------------- 465 if(optatm==1) then 466 ! Accumulate V^AT(G)*SF(G) or n^AT(G)*SF(G) 467 if (optv==1) then 468 workv(re,ii)=workv(re,ii)+sfr*v_at 469 workv(im,ii)=workv(im,ii)+sfi*v_at 470 end if 471 if (optn==1) then 472 workn(re,ii)=workn(re,ii)+sfr*n_at 473 workn(im,ii)=workn(im,ii)+sfi*n_at 474 end if 475 476 ! Compute contrib. to forces and/or frozen part of dyn. matrix 477 ! ------------------------------------------------------------- 478 else if (optgr==1.or.optdyfr==1) then 479 dbl_ig1=dble(ig1_);dbl_ig2=dble(ig2_);dbl_ig3=dble(ig3_) 480 ! Compute (2Pi)*V^AT(G)*rho(G) or (2Pi)*n^AT(G)*V(G) 481 if (optv==1) then 482 tmpvr=(two_pi*v_at)*rhog(re,ii) 483 tmpvi=(two_pi*v_at)*rhog(im,ii) 484 end if 485 if (optn==1) then 486 tmpnr=(two_pi*n_at)*vg(re,ii) 487 tmpni=(two_pi*n_at)*vg(im,ii) 488 end if 489 ! === contrib. to forces 490 if (optgr==1) then 491 ! Accumulate -(2Pi.G)*V^AT(G)*rho(G)*SF(G) 492 ! or -(2Pi)*n^AT(G)*V(G)*SF(G) into forces 493 if (optv==1) then 494 do ia=ia1,ia2 495 term=tmpvi*phre_igia(ia)+tmpvr*phim_igia(ia) 496 grv_indx(1,ia)=grv_indx(1,ia)-dbl_ig1*term 497 grv_indx(2,ia)=grv_indx(2,ia)-dbl_ig2*term 498 grv_indx(3,ia)=grv_indx(3,ia)-dbl_ig3*term 499 end do 500 end if 501 if (optn==1) then 502 do ia=ia1,ia2 503 term=tmpni*phre_igia(ia)+tmpnr*phim_igia(ia) 504 grn_indx(1,ia)=grn_indx(1,ia)-dbl_ig1*term 505 grn_indx(2,ia)=grn_indx(2,ia)-dbl_ig2*term 506 grn_indx(3,ia)=grn_indx(3,ia)-dbl_ig3*term 507 end do 508 end if 509 end if 510 ! === contrib. to frozen part of dyn. matrix 511 if (optdyfr==1) then 512 ! Accumulate -(2Pi^2.Gi.Gj)*V^AT(G)*rho(G)*SF(G) 513 ! or -(2Pi^2.Gi.Gj)*n^AT(G)*V(G)*SF(G) into dyn. matrix 514 if (optv==1) then 515 do ia=ia1,ia2 516 term=two_pi*(tmpvr*phre_igia(ia)-tmpvi*phim_igia(ia)) 517 dyfrv_indx(1,1,ia)=dyfrv_indx(1,1,ia)-dbl_ig1*dbl_ig1*term 518 dyfrv_indx(1,2,ia)=dyfrv_indx(1,2,ia)-dbl_ig1*dbl_ig2*term 519 dyfrv_indx(1,3,ia)=dyfrv_indx(1,3,ia)-dbl_ig1*dbl_ig3*term 520 dyfrv_indx(2,2,ia)=dyfrv_indx(2,2,ia)-dbl_ig2*dbl_ig2*term 521 dyfrv_indx(2,3,ia)=dyfrv_indx(2,3,ia)-dbl_ig2*dbl_ig3*term 522 dyfrv_indx(3,3,ia)=dyfrv_indx(3,3,ia)-dbl_ig3*dbl_ig3*term 523 end do 524 end if 525 if (optn==1) then 526 do ia=ia1,ia2 527 term=two_pi*(tmpnr*phre_igia(ia)-tmpni*phim_igia(ia)) 528 dyfrn_indx(1,1,ia)=dyfrn_indx(1,1,ia)-dbl_ig1*dbl_ig1*term 529 dyfrn_indx(1,2,ia)=dyfrn_indx(1,2,ia)-dbl_ig1*dbl_ig2*term 530 dyfrn_indx(1,3,ia)=dyfrn_indx(1,3,ia)-dbl_ig1*dbl_ig3*term 531 dyfrn_indx(2,2,ia)=dyfrn_indx(2,2,ia)-dbl_ig2*dbl_ig2*term 532 dyfrn_indx(2,3,ia)=dyfrn_indx(2,3,ia)-dbl_ig2*dbl_ig3*term 533 dyfrn_indx(3,3,ia)=dyfrn_indx(3,3,ia)-dbl_ig3*dbl_ig3*term 534 end do 535 end if 536 end if 537 end if 538 539 ! Compute (dV^AT(q)/dq)/q and/or (dn^AT(q)/dq)/q 540 ! For stress tensor and/or elastic tensor 541 ! --------------------------------- 542 if (optstr==1.or.opteltfr==1) then 543 ! Note: correction of Numerical Recipes sign error before (3._dp*aa**2-1._dp) 544 ! ee*dqm1 + ff*dqdiv6 is the best estimate of dV(q)/dq from splines 545 if (optv==1) then 546 if (have_g0) then 547 dv_at=zero 548 else 549 ee=vspl(jj+1,1,itypat)-vspl(jj,1,itypat) 550 ff=(3._dp*bb**2-1._dp)*vspl(jj+1,2,itypat)& 551 & -(3._dp*aa**2-1._dp)*vspl(jj ,2,itypat) 552 dv_at=((ee*dqm1+ff*dqdiv6)/gmag-2.0_dp*v_at)/gsquar 553 end if 554 end if 555 if (optn==1) then 556 if (have_g0) then 557 if (optn2==1) then 558 if (usepaw ==1) then 559 dn_at=pawtab(itypat)%dncdq0 560 else 561 dn_at=psps%nctab(itypat)%dncdq0 562 end if 563 else if (optn2==2) then 564 if (usepaw == 1) then 565 dn_at=pawtab(itypat)%dnvdq0 566 else 567 dn_at=psps%nctab(itypat)%dnvdq0 568 end if 569 else if (optn2==3) then 570 dn_at=-two*gauss1*alf2pi2 571 else if (optn2==4.and.usepaw==1) then 572 dn_at=pawtab(itypat)%dtaucdq0 573 end if 574 if (opteltfr==1) then 575 d2n_at = 0 576 end if 577 else 578 if (optn2==1) then 579 ee=tcorespl(jj+1,1)-tcorespl(jj,1) 580 ff=(3._dp*bb**2-1._dp)*tcorespl(jj+1,2) & 581 & -(3._dp*aa**2-1._dp)*tcorespl(jj,2) 582 ! Also get nc''(q) 583 if (opteltfr==1) gg=aa*tcorespl(jj,2)+bb*tcorespl(jj+1,2) 584 else if (optn2==2) then 585 ee=tvalespl(jj+1,1)-tvalespl(jj,1) 586 ff=(3._dp*bb**2-1._dp)*tvalespl(jj+1,2) & 587 & -(3._dp*aa**2-1._dp)*tvalespl(jj,2) 588 ! Also get nc''(q) 589 if (opteltfr==1) then 590 gg=aa*tvalespl(jj,2)+bb*tvalespl(jj+1,2) 591 end if 592 else if (optn2==3) then 593 dn_at=-two*gauss1*alf2pi2*exp(-gsquar*alf2pi2) 594 else if (optn2==4.and.usepaw==1) then 595 ee=ttaucorespl(jj+1,1)-ttaucorespl(jj,1) 596 ff=(3._dp*bb**2-1._dp)*ttaucorespl(jj+1,2) & 597 & -(3._dp*aa**2-1._dp)*ttaucorespl(jj,2) 598 ! Also get nc''(q) 599 if (opteltfr==1) gg=aa*ttaucorespl(jj,2)+bb*ttaucorespl(jj+1,2) 600 end if 601 dn_at = (ee*dqm1+ff*dqdiv6)/gmag 602 if (opteltfr==1) then 603 d2n_at = (gg-dn_at)/gsquar 604 d2n_at2 = gg/gmag**3 605 end if 606 end if 607 end if 608 609 ! Compute G in cartesian coordinates 610 gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+& 611 & gprimd(1,3)*dble(ig3) 612 gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+& 613 & gprimd(2,3)*dble(ig3) 614 gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+& 615 & gprimd(3,3)*dble(ig3) 616 end if 617 618 ! Compute contrib. to stress tensor 619 ! --------------------------------- 620 if (optstr==1)then 621 ! Accumulate -dV^AT/dG*rho(G)*SF(G)*Gi.Gj/G 622 ! or -dn^AT/dG*V(G)*SF(G)*Gi.Gj/G 623 ! into stress tensor 624 if (optv==1) then 625 term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi) 626 strv(1)=strv(1)-term*(dv_at*gcart(1)*gcart(1)+v_at) 627 strv(2)=strv(2)-term*(dv_at*gcart(2)*gcart(2)+v_at) 628 strv(3)=strv(3)-term*(dv_at*gcart(3)*gcart(3)+v_at) 629 strv(4)=strv(4)-term*dv_at*gcart(3)*gcart(2) 630 strv(5)=strv(5)-term*dv_at*gcart(3)*gcart(1) 631 strv(6)=strv(6)-term*dv_at*gcart(2)*gcart(1) 632 end if 633 if (optn==1) then 634 term=(vg(re,ii)*sfr+vg(im,ii)*sfi)*dn_at 635 strn(1)=strn(1)-term*gcart(1)*gcart(1) 636 strn(2)=strn(2)-term*gcart(2)*gcart(2) 637 strn(3)=strn(3)-term*gcart(3)*gcart(3) 638 strn(4)=strn(4)-term*gcart(3)*gcart(2) 639 strn(5)=strn(5)-term*gcart(3)*gcart(1) 640 strn(6)=strn(6)-term*gcart(2)*gcart(1) 641 end if 642 end if 643 644 ! Compute contrib. to elastic tensor 645 ! --------------------------------- 646 if (opteltfr==1) then 647 dbl_ig1=dble(ig1_);dbl_ig2=dble(ig2_);dbl_ig3=dble(ig3_) 648 ! if (ig1==0 .and. ig2==0 .and. ig3==0) cycle 649 ! Compute G*dG/Deps_{\gamme\delta} 650 dg2=0.5_dp*dgsqds_atm(ig1,ig2,ig3,is2_in) 651 652 term = vg (re,ii)*sfr + vg (im,ii)*sfi 653 term1 = vg1(re,ii)*sfr + vg1(im,ii)*sfi 654 term2 = vg1_core(re,ii)*sfr + vg1_core(im,ii)*sfi 655 656 do is1=1,6 657 658 ! Compute G*dG/Deps_{\alpha\beta} 659 dg1=0.5_dp*dgsqds_atm(ig1,ig2,ig3,is1) 660 ! Compute G^2*d2G/Deps_{alphabetagammadelta} 661 d2g=(0.25_dp*d2gsqds_atm(ig1,ig2,ig3,is1,is2_in)) 662 663 eltfrn(is1,is2_in) = eltfrn(is1,is2_in) + (term*(d2n_at*dg1*dg2 + d2g*dn_at)) 664 eltfrn(is1,is2_in) = eltfrn(is1,is2_in) + 0.5*(term1*dn_at*dg1) 665 eltfrn(is2_in,is1) = eltfrn(is2_in,is1) + 0.5*(term1*dn_at*dg1) 666 667 if(is2_in<=3)then 668 eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - term*dn_at*dg1 669 end if 670 if(is1<=3)then 671 eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - term*dn_at*dg2 672 end if 673 if(is2_in<=3.and.is1<=3)then 674 eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - (term1-term)*n_at 675 end if 676 end do 677 678 ! internal strain 679 do ia=ia1,ia2 680 js=7+3*(ia-1) 681 ! Compute -2pi*G*i*vxcis_core(G)*nat(G)*(exp(-iGr)) 682 term=(vg1_core(im,ii)*phre_igia(ia)+vg1_core(re,ii)*phim_igia(ia))*n_at*two_pi 683 eltfrn(js ,is2_in) = eltfrn(js ,is2_in) - dbl_ig1*term 684 eltfrn(js+1,is2_in) = eltfrn(js+1,is2_in) - dbl_ig2*term 685 eltfrn(js+2,is2_in) = eltfrn(js+2,is2_in) - dbl_ig3*term 686 687 ! Compute -2pi*G*i*vxc(G)*(dnat*dG/deps-delta*nat(G))*(exp(-iGr)) 688 term=(vg(im,ii)*phre_igia(ia)+vg(re,ii)*phim_igia(ia))*& 689 & (dn_at*dg2-delta(is2_in)*n_at)*two_pi 690 eltfrn(js ,is2_in) = eltfrn(js ,is2_in) - dbl_ig1*term 691 eltfrn(js+1,is2_in) = eltfrn(js+1,is2_in) - dbl_ig2*term 692 eltfrn(js+2,is2_in) = eltfrn(js+2,is2_in) - dbl_ig3*term 693 end do 694 end if 695 ! End skip G**2 outside cutoff: 696 end if 697 ! End loop on n1, n2, n3 698 end do 699 end if ! this plane is for me_fft 700 end do 701 end do 702 703 ! Symmetrize the dynamical matrix with respect to indices 704 if (optdyfr==1) then 705 if (optv==1) then 706 do ia=ia1,ia2 707 dyfrv_indx(2,1,ia)=dyfrv_indx(1,2,ia) 708 dyfrv_indx(3,1,ia)=dyfrv_indx(1,3,ia) 709 dyfrv_indx(3,2,ia)=dyfrv_indx(2,3,ia) 710 end do 711 end if 712 if (optn==1) then 713 do ia=ia1,ia2 714 dyfrn_indx(2,1,ia)=dyfrn_indx(1,2,ia) 715 dyfrn_indx(3,1,ia)=dyfrn_indx(1,3,ia) 716 dyfrn_indx(3,2,ia)=dyfrn_indx(2,3,ia) 717 end do 718 end if 719 end if 720 721 ia1=ia2+1 722 723 ! End loop on type of atoms 724 end do 725 726 ABI_FREE(phre_igia) 727 ABI_FREE(phim_igia) 728 ABI_FREE(gcutoff) 729 730 !Get local potential or density back to real space 731 if(optatm==1)then 732 ! Allow for the addition of a perturbing potential 733 if (optv==1) then 734 if ((vprtrb(1)**2+vprtrb(2)**2) > 1.d-30) then 735 ! Find the linear indices which correspond with the input wavevector qprtrb 736 ! The double modulus handles both i>=n and i<0, mapping into [0,n-1]; 737 ! then add 1 to get range [1,n] for each 738 i3=1+mod(n3+mod(qprtrb(3),n3),n3) 739 i2=1+mod(n2+mod(qprtrb(2),n2),n2) 740 i1=1+mod(n1+mod(qprtrb(1),n1),n1) 741 ! Compute the linear index in the 3 dimensional array 742 ii=i1+n1*((ffti2_local(i2)-1)+(n2/nproc_fft)*(i3-1)) 743 ! Add in the perturbation at G=qprtrb 744 workv(re,ii)=workv(re,ii)+0.5_dp*vprtrb(1) 745 workv(im,ii)=workv(im,ii)+0.5_dp*vprtrb(2) 746 ! Same thing for G=-qprtrb 747 i3=1+mod(n3+mod(-qprtrb(3),n3),n3) 748 i2=1+mod(n2+mod(-qprtrb(2),n2),n2) 749 i1=1+mod(n1+mod(-qprtrb(1),n1),n1) 750 ! ii=i1+n1*((i2-1)+n2*(i3-1)) 751 workv(re,ii)=workv(re,ii)+0.5_dp*vprtrb(1) 752 workv(im,ii)=workv(im,ii)-0.5_dp*vprtrb(2) 753 write(message, '(a,1p,2e12.4,a,0p,3i4,a)' )& 754 & ' atm2fft: perturbation of vprtrb=', vprtrb,& 755 & ' and q=',qprtrb,' has been added' 756 call wrtout(std_out,message,'COLL') 757 end if 758 end if 759 760 if (optv==1.or.optn==1) then 761 ! Create fake mpi_enreg to wrap fourdp 762 call initmpi_seq(mpi_enreg_fft) 763 ABI_FREE(mpi_enreg_fft%distribfft) 764 if (present(comm_fft)) then 765 call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb) 766 my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb 767 else 768 my_comm_fft=xmpi_comm_self;paral_kgb_fft=0; 769 mpi_enreg_fft%distribfft => my_distribfft 770 end if 771 ! Non-symetrized non-zero elements have to be nullified 772 ! Transform back to real space; divide by unit cell volume 773 xnorm=one/ucvol 774 if (optv==1) then 775 call zerosym(workv,2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft) 776 call fourdp(1,workv,atmvloc,1,mpi_enreg_fft,nfft,1,ngfft,0) 777 atmvloc(:)=atmvloc(:)*xnorm 778 ABI_FREE(workv) 779 end if 780 if (optn==1) then 781 call zerosym(workn,2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft) 782 call fourdp(1,workn,atmrho,1,mpi_enreg_fft,nfft,1,ngfft,0) 783 atmrho(:)=atmrho(:)*xnorm 784 ABI_FREE(workn) 785 end if 786 ! Destroy fake mpi_enreg 787 call unset_mpi_enreg_fft(mpi_enreg_fft) 788 end if 789 790 end if 791 792 !Additional treatment in case of parallelization 793 if (present(comm_fft)) then 794 if ((xmpi_comm_size(comm_fft)>1).and.& 795 & (optgr==1.or.optstr==1.or.optdyfr==1.or.opteltfr==1)) then 796 call timab(48,1,tsec) 797 if (optv==1) then 798 if (optgr==1)then 799 call xmpi_sum(grv_indx,comm_fft,ierr) 800 end if 801 if (optstr==1)then 802 call xmpi_sum(strv,comm_fft,ierr) 803 end if 804 if (optdyfr==1)then 805 call xmpi_sum(dyfrv_indx,comm_fft,ierr) 806 end if 807 end if 808 if (optn==1) then 809 if (optgr==1)then 810 call xmpi_sum(grn_indx,comm_fft,ierr) 811 end if 812 if (optstr==1)then 813 call xmpi_sum(strn,comm_fft,ierr) 814 end if 815 if (optdyfr==1)then 816 call xmpi_sum(dyfrn_indx,comm_fft,ierr) 817 end if 818 end if 819 call timab(48,2,tsec) 820 end if 821 end if 822 823 !Forces: re-order atoms 824 if (optgr==1) then 825 if (optv==1) then 826 do ia=1,natom 827 grv(1:3,atindx1(ia))=grv_indx(1:3,ia) 828 end do 829 ABI_FREE(grv_indx) 830 end if 831 if (optn==1) then 832 do ia=1,natom 833 grn(1:3,atindx1(ia))=grn_indx(1:3,ia) 834 end do 835 ABI_FREE(grn_indx) 836 end if 837 end if 838 839 !Elastic tensor: Fill in lower triangle 840 ! if (opteltfr==1) then 841 ! do is2=2,6 842 ! do is1=1,is2-1 843 ! eltfrn(is2,is1)=eltfrn(is1,is2) 844 ! end do 845 ! end do 846 ! end if 847 848 !Normalize stress tensor: 849 if (optstr==1) then 850 if (optv==1) then 851 strv(:)=strv(:)/ucvol 852 end if 853 if (optn==1) strn(:)=strn(:)/ucvol 854 end if 855 856 !Dynamical matrix: re-order atoms 857 if (optdyfr==1) then 858 if (optv==1) then 859 do ia=1,natom 860 dyfrv(1:3,1:3,atindx1(ia))=dyfrv_indx(1:3,1:3,ia) 861 end do 862 ABI_FREE(dyfrv_indx) 863 end if 864 if (optn==1) then 865 do ia=1,natom 866 dyfrn(1:3,1:3,atindx1(ia))=dyfrn_indx(1:3,1:3,ia) 867 end do 868 ABI_FREE(dyfrn_indx) 869 end if 870 end if 871 872 if (.not.present(distribfft)) then 873 call destroy_distribfft(my_distribfft) 874 ABI_FREE(my_distribfft) 875 end if 876 877 DBG_EXIT("COLL") 878 879 contains 880 881 function gsq_atm(i1,i2,i3) 882 883 real(dp) :: gsq_atm 884 integer,intent(in) :: i1,i2,i3 885 gsq_atm=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+dble(i3*i3)*gmet(3,3) & 886 & +two*(dble(i1*i2)*gmet(1,2)+dble(i2*i3)*gmet(2,3)+dble(i3*i1)*gmet(3,1)) 887 end function gsq_atm 888 889 function dgsqds_atm(i1,i2,i3,is) 890 !Define dG^2/ds based on G space metric derivative 891 real(dp) :: dgsqds_atm 892 integer,intent(in) :: i1,i2,i3,is 893 dgsqds_atm=dble(i1*i1)*dgm(1,1,is)+dble(i2*i2)*dgm(2,2,is)+& 894 & dble(i3*i3)*dgm(3,3,is)+& 895 & dble(i1*i2)*(dgm(1,2,is)+dgm(2,1,is))+& 896 & dble(i1*i3)*(dgm(1,3,is)+dgm(3,1,is))+& 897 & dble(i2*i3)*(dgm(2,3,is)+dgm(3,2,is)) 898 end function dgsqds_atm 899 900 function d2gsqds_atm(i1,i2,i3,is1,is2) 901 ! Define 2dG^2/ds1ds2 based on G space metric derivative 902 real(dp) :: d2gsqds_atm 903 integer,intent(in) :: i1,i2,i3,is1,is2 904 d2gsqds_atm=dble(i1*i1)*d2gm(1,1,is1,is2)+& 905 & dble(i2*i2)*d2gm(2,2,is1,is2)+dble(i3*i3)*d2gm(3,3,is1,is2)+& 906 & dble(i1*i2)*(d2gm(1,2,is1,is2)+d2gm(2,1,is1,is2))+& 907 & dble(i1*i3)*(d2gm(1,3,is1,is2)+d2gm(3,1,is1,is2))+& 908 & dble(i2*i3)*(d2gm(2,3,is1,is2)+d2gm(3,2,is1,is2)) 909 end function d2gsqds_atm 910 911 end subroutine atm2fft
ABINIT/dfpt_atm2fft [ Functions ]
NAME
dfpt_atm2fft
FUNCTION
This routine sums 1st-order atomic functions (density or potential) defined (in rec. space) on a radial grid to get global 1st-order quantities on the fine FFT grid. Possible options: optn=1: compute a sum of local 1st-order atomic densities optv=1: compute a sum of local 1st-order atomic potentials
INPUTS
atindx(natom)=index table for atoms ordered by type cplex: if 1, real space 1-order functions on FFT grid distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism eei=local pseudopotential part of total energy 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) idir=direction of atomic displacement (in case of phonons perturb.) used only if ndir=1 (see below) ipert=nindex of perturbation 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. ndir=number of directions of atomic displacement (in case of phonon): can be 1 (idir direction in then used) or 3 (all directions) 6 cartesian strain component 11,22,33,32,31,21 (in case strain perturbation) nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT nattyp(ntypat)=array describing how many atoms of each type in cell ntypat=number of types of atoms. optn,optn2,optv= (see NOTES below) paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present) 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. qphon(3)=wavevector of the phonon typat(natom)=type of each atom 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) xred(3,natom)=reduced atomic coordinates psps <type(pseudopotential_type)>=variables related to pseudopotentials
OUTPUT
======= if optv==1 ======= atmvlocr1(cplex*nfft,ndir)=sum of local 1st-order atomic potentials in real space atmvlocg1(2,nfft,ndir)=sum of local 1st-order atomic potentials in G-space ======= if optn==1 ======= --- if optatm==1 atmrhor1(cplex*nfft,ndir)=sum of 1st-order atomic densities in real space atmrhog1(2,nfft,ndir)=sum of 1st-order atomic densities in G-space
NOTES
Details on possible options: ============================ optv: controls the computation of a local 1st-order potential as sum of atomic potentials Vloc(r)=Sum_R[V1^AT(r-R)] optn: controls the computation of a 1st-order density as sum of atomic densities n(r)=Sum_R[n1^AT(r-R)] 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 Typical uses: ============= Computation of: - 1st-order local potential: optv=1 - 1st-order PS core density: optn=1, optn2=1
SOURCE
995 subroutine dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,& 996 & mgfft,mqgrid,natom,ndir,nfft,ngfft,ntypat,& 997 & ph1d,qgrid,qphon,typat,ucvol,usepaw,xred,psps,pawtab,& 998 & atmrhor1,atmrhog1,atmvlocr1,atmvlocg1,distribfft,gauss,comm_fft,me_g0,optn_in,& 999 & optn2_in,optv_in,paral_kgb,vspl) ! optional arguments 1000 1001 !Arguments ------------------------------------ 1002 !scalars 1003 integer,intent(in) :: cplex,idir,ipert,mgfft,mqgrid,natom,ndir,nfft,ntypat,usepaw 1004 integer,optional,intent(in) :: optn_in,optn2_in,optv_in 1005 integer,optional,intent(in) :: me_g0,comm_fft,paral_kgb 1006 real(dp),intent(in) :: gsqcut,ucvol 1007 type(pseudopotential_type),target,intent(in) :: psps 1008 type(distribfft_type),optional,intent(in),target :: distribfft 1009 !arrays 1010 integer,intent(in) :: atindx(natom),ngfft(18),typat(natom) 1011 real(dp),intent(in) :: gmet(3,3),gprimd(3,3) 1012 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid),qphon(3) 1013 real(dp),intent(in) :: xred(3,natom) 1014 real(dp),optional,intent(in) :: gauss(2,ntypat),vspl(mqgrid,2,ntypat) 1015 real(dp),optional,intent(out) :: atmrhor1(cplex*nfft,ndir) 1016 real(dp),optional,intent(out) :: atmrhog1(2,nfft,ndir) 1017 real(dp),optional,intent(out) :: atmvlocr1(cplex*nfft,ndir) 1018 real(dp),optional,intent(out) :: atmvlocg1(2,nfft,ndir) 1019 type(pawtab_type),target,intent(in) :: pawtab(ntypat*usepaw) 1020 1021 !Local variables ------------------------------ 1022 !scalars 1023 integer,parameter :: im=2,re=1 1024 integer :: i1,i2,i3,ia,ia1,ia2,iatm,iatom,id,id1,id2,id3 1025 integer :: ig1,ig1max,ig1min,ig2,ig2max,ig2min,ig3,ig3max,ig3min 1026 integer :: ii,itypat,jj,me_fft,my_comm_fft,n1,n2,n3,nattyp,nproc_fft,ntype,paral_kgb_fft 1027 integer :: optn,optv,optn2,shift1,shift2,shift3,type1,type2 1028 logical :: have_g0,qeq0,qeq05 1029 #ifdef FC_NVHPC 1030 !Silly trick to prevent NVHPC optimization issue 1031 logical :: nothing=.false. 1032 #endif 1033 real(dp),parameter :: tolfix=1.0000001_dp 1034 real(dp) :: aa,alf2pi2,bb,cc,cutoff,dd,diff,dq,dq2div6,dqdiv6,dqm1,ee,ff 1035 real(dp) :: gauss1,gauss2,gmag,gq1,gq2,gq3,gsquar,n_at,dn_at,ph12i,ph12r,ph1i 1036 real(dp) :: ph1r,ph2i,ph2r,ph3i,ph3r,phqim,phqre,qxred2pi 1037 real(dp) :: sfi,sfqi,sfqr,sfr,term_n,term_v,v_at,dv_at,xnorm 1038 type(distribfft_type),pointer :: my_distribfft 1039 type(mpi_type) :: mpi_enreg_fft 1040 !arrays 1041 integer :: eps1(6)=(/1,2,3,2,3,1/),eps2(6)=(/1,2,3,3,1,2/),jdir(ndir) 1042 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:) 1043 real(dp), ABI_CONTIGUOUS pointer :: tvalespl(:,:),tcorespl(:,:) 1044 real(dp) :: gq(6),gcart(3) 1045 real(dp),allocatable :: phim_igia(:),phre_igia(:),workn(:,:,:),workv(:,:,:) 1046 1047 !no_abirules 1048 !Define G^2 based on G space metric gmet. 1049 ! gsq(g1,g2,g3)=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+g3*g3*gmet(3,3) & 1050 ! & +two*(g1*g2*gmet(1,2)+g2*g3*gmet(2,3)+g3*g1*gmet(3,1)) 1051 ! ************************************************************************* 1052 1053 DBG_ENTER("COLL") 1054 1055 ! Check optional arguments 1056 if (present(comm_fft)) then 1057 if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then 1058 ABI_BUG('Need paral_kgb and me_g0 with comm_fft !') 1059 end if 1060 end if 1061 1062 if (present(gauss))then 1063 if (.not.present(optn2_in)) then 1064 optn2 = 3 1065 else 1066 if(optn2_in/=3)then 1067 ABI_BUG('optn2 must be set to 3!') 1068 else 1069 optn2 = optn2_in 1070 end if 1071 end if 1072 end if 1073 1074 if (present(atmrhor1))then 1075 if(.not.present(optn_in))then 1076 optn = 1 1077 else 1078 optn = optn_in 1079 end if 1080 if (.not.present(optn2_in)) then 1081 ABI_BUG('rho1 calculation need optn2 !') 1082 else 1083 optn2 = optn2_in 1084 end if 1085 else 1086 optn = 0 1087 optn2 = 0 1088 end if 1089 1090 if (present(atmvlocr1))then 1091 if(.not.present(optv_in))then 1092 optv = 1 1093 else 1094 if(optv_in/=1)then 1095 ABI_BUG('optv_in must be set to 1!') 1096 else 1097 optv = optv_in 1098 end if 1099 end if 1100 if(.not.present(vspl))then 1101 ABI_BUG('vloc1 calculation need vspl!') 1102 end if 1103 else 1104 optv = 0 1105 end if 1106 1107 if(ipert==natom+1.or.ipert==natom+2.or.ipert==natom+10.or.ipert==natom+11) then 1108 1109 ! (In case of d/dk or an electric/magnetic field) 1110 if (optn==1) then 1111 atmrhor1(1:cplex*nfft,1:ndir)=zero 1112 if (present(atmrhog1)) atmrhog1 = zero 1113 end if 1114 if (optv==1) then 1115 atmvlocr1(1:cplex*nfft,1:ndir)=zero 1116 if (present(atmvlocg1)) atmvlocg1 = zero 1117 end if 1118 1119 else 1120 1121 ! Useful quantities 1122 if (ipert/=natom+3.and.ipert/=natom+4) then 1123 iatom=ipert;iatm=atindx(iatom) 1124 itypat=typat(iatom) 1125 else 1126 !sum of all (strain pertubation) 1127 iatom = 1 1128 iatm = 1 1129 itypat = 1 1130 end if 1131 1132 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1133 me_fft=ngfft(11) 1134 nproc_fft=ngfft(10) 1135 if (ndir==1) then 1136 jdir(1)=idir 1137 else 1138 do id=1,ndir 1139 jdir(id)=id 1140 end do 1141 end if 1142 1143 ! Get the distrib associated with this fft_grid 1144 if (present(distribfft)) then 1145 my_distribfft => distribfft 1146 else 1147 ABI_MALLOC(my_distribfft,) 1148 call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp') 1149 end if 1150 if (n2==my_distribfft%n2_coarse) then 1151 fftn2_distrib => my_distribfft%tab_fftdp2_distrib 1152 else if (n2 == my_distribfft%n2_fine) then 1153 fftn2_distrib => my_distribfft%tab_fftdp2dg_distrib 1154 else 1155 ABI_BUG("Unable to find an allocated distrib for this fft grid") 1156 end if 1157 1158 qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) 1159 qeq05=(abs(abs(qphon(1))-half)<tol12.or. & 1160 & abs(abs(qphon(2))-half)<tol12.or. & 1161 & abs(abs(qphon(3))-half)<tol12) 1162 1163 if (nproc_fft>1.and.qeq05) then 1164 ABI_ERROR('not compatible with FFT parallelism') 1165 end if 1166 1167 if (optn2==3)then 1168 gauss1=gauss(1,itypat) 1169 gauss2=gauss(2,itypat) 1170 alf2pi2=(two_pi*gauss2)**2 1171 end if 1172 1173 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 1174 dqm1=one/dq 1175 dqdiv6=dq/six 1176 dq2div6=dq**2/six 1177 cutoff=gsqcut*tolfix 1178 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 1179 ig1max=-1;ig2max=-1;ig3max=-1 1180 ig1min=n1;ig2min=n2;ig3min=n3 1181 1182 ! Determination of phase qxred* 1183 qxred2pi=two_pi*(qphon(1)*xred(1,iatom)+ & 1184 & qphon(2)*xred(2,iatom)+ & 1185 & qphon(3)*xred(3,iatom) ) 1186 phqre=cos(qxred2pi) 1187 phqim=sin(qxred2pi) 1188 1189 ! Zero out temporary arrays 1190 if (optv==1) then 1191 ABI_MALLOC(workv,(2,nfft,ndir)) 1192 workv(:,:,:)=zero 1193 end if 1194 if (optn==1) then 1195 ABI_MALLOC(workn,(2,nfft,ndir)) 1196 workn(:,:,:)=zero 1197 end if 1198 1199 if (ipert==natom+3.or.ipert==natom+4) then 1200 ntype = ntypat 1201 ia1=1 1202 type1 = 1 1203 type2 = ntype 1204 ABI_MALLOC(phre_igia,(natom)) 1205 ABI_MALLOC(phim_igia,(natom)) 1206 else 1207 type1 = itypat 1208 type2 = itypat 1209 ntype = 1 1210 ABI_MALLOC(phre_igia,(iatm:iatm)) 1211 ABI_MALLOC(phim_igia,(iatm:iatm)) 1212 end if 1213 1214 1215 do itypat=type1,type2 1216 ! ia1,ia2 sets range of loop over atoms: 1217 if (ipert==natom+3.or.ipert==natom+4) then 1218 nattyp = count(typat(:)==itypat) 1219 ia2=ia1+nattyp-1 1220 else 1221 ia1 = iatm 1222 ia2 = iatm 1223 end if 1224 1225 if (usepaw == 1) then 1226 tcorespl => pawtab(itypat)%tcorespl 1227 tvalespl => pawtab(itypat)%tvalespl 1228 else 1229 tcorespl => psps%nctab(itypat)%tcorespl 1230 tvalespl => psps%nctab(itypat)%tvalespl 1231 end if 1232 1233 ii=0 1234 do i3=1,n3 1235 ig3=i3-(i3/id3)*n3-1 1236 gq3=dble(ig3)+qphon(3) 1237 gq(3)=gq3 1238 1239 do i2=1,n2 1240 if (fftn2_distrib(i2)==me_fft) then 1241 ig2=i2-(i2/id2)*n2-1 1242 gq2=dble(ig2)+qphon(2) 1243 gq(2)=gq2 1244 1245 do i1=1,n1 1246 ig1=i1-(i1/id1)*n1-1 1247 gq1=dble(ig1)+qphon(1) 1248 gq(1)=gq1 1249 1250 ii=ii+1 1251 ! gsquar=gsq(gq1,gq2,gq3) 1252 gsquar=gq1*gq1*gmet(1,1)+gq2*gq2*gmet(2,2)+gq3*gq3*gmet(3,3) & 1253 & +two*(gq1*gq2*gmet(1,2)+gq2*gq3*gmet(2,3)+gq3*gq1*gmet(3,1)) 1254 1255 1256 ! Skip G**2 outside cutoff: 1257 if (gsquar<=cutoff) then 1258 1259 ! Identify min/max indexes (to cancel unbalanced contributions later) 1260 if (qeq05) then 1261 ig1max=max(ig1max,ig1);ig1min=min(ig1min,ig1) 1262 ig2max=max(ig2max,ig2);ig2min=min(ig2min,ig2) 1263 ig3max=max(ig3max,ig3);ig3min=min(ig3min,ig3) 1264 end if 1265 1266 gmag=sqrt(gsquar) 1267 have_g0=(ig1==0.and.ig2==0.and.ig3==0.and.qeq0) 1268 1269 jj=1+int(gmag*dqm1) 1270 diff=gmag-qgrid(jj) 1271 1272 ! Compute structure factor 1273 phre_igia(:) = zero 1274 phim_igia(:) = zero 1275 1276 do ia=ia1,ia2 1277 #ifdef FC_NVHPC 1278 !Silly trick to prevent NVHPC optimization issue 1279 if(nothing) write(100,*) shift1,shift2,shift3 1280 #endif 1281 shift1=1+n1+(ia-1)*(2*n1+1) 1282 shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1) 1283 shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1) 1284 ph1r=ph1d(re,ig1+shift1);ph1i=ph1d(im,ig1+shift1) 1285 ph2r=ph1d(re,ig2+shift2);ph2i=ph1d(im,ig2+shift2) 1286 ph3r=ph1d(re,ig3+shift3);ph3i=ph1d(im,ig3+shift3) 1287 ph12r=ph1r*ph2r-ph1i*ph2i 1288 ph12i=ph1r*ph2i+ph1i*ph2r 1289 phre_igia(ia)=ph12r*ph3r-ph12i*ph3i 1290 phim_igia(ia)=ph12r*ph3i+ph12i*ph3r 1291 end do 1292 ! Compute V^AT(g+q) and/or n^AT(g+q) for given type of atom 1293 ! Evaluate spline fit: p. 86 Numerical Recipes, Press et al; 1294 ! Note the error in book for sign of "aa" term in derivative. 1295 if (optv==1.or.optn2/=3) then 1296 bb = diff*dqm1 1297 aa = one-bb 1298 cc = aa*(aa**2-one)*dq2div6 1299 dd = bb*(bb**2-one)*dq2div6 1300 end if 1301 1302 if (optv==1) then 1303 if (have_g0) then 1304 v_at=zero 1305 dv_at=zero 1306 else 1307 v_at=(aa*vspl(jj,1,itypat)+bb*vspl(jj+1,1,itypat)+& 1308 & cc*vspl(jj,2,itypat)+dd*vspl(jj+1,2,itypat))/gsquar 1309 1310 ! Also get (dV(q)/dq)/q: 1311 ! (note correction of Numerical Recipes sign error 1312 ! before (3._dp*aa**2-1._dp) 1313 ee= vspl(jj+1,1,itypat)-vspl(jj,1,itypat) 1314 ff= (3._dp*bb**2-1._dp)*vspl(jj+1,2,itypat) & 1315 & - (3._dp*aa**2-1._dp)*vspl(jj,2,itypat) 1316 dv_at = ( ( ee*dqm1 + ff*dqdiv6 )/gmag& 1317 & - 2.0_dp*v_at) / gsquar 1318 end if 1319 end if ! end optv 1320 1321 if (optn==1) then 1322 if (optn2==1) then 1323 n_at=(aa*tcorespl(jj,1)+bb*tcorespl(jj+1,1)+cc*tcorespl(jj,2)+dd*tcorespl(jj+1,2)) 1324 else if (optn2==2) then 1325 n_at=(aa*tvalespl(jj,1)+bb*tvalespl(jj+1,1)+cc*tvalespl(jj,2)+dd*tvalespl(jj+1,2)) 1326 else if (optn2==3) then 1327 n_at=gauss1*exp(-gsquar*alf2pi2) 1328 else 1329 n_at=zero 1330 end if 1331 1332 ! Also get (dn^AT(q)/dq)/q: 1333 if (have_g0) then 1334 if (optn2==1) then 1335 if (usepaw == 1) then 1336 dn_at=pawtab(itypat)%dncdq0 1337 else 1338 dn_at=psps%nctab(itypat)%dncdq0 1339 end if 1340 else if (optn2==2) then 1341 if (usepaw == 1) then 1342 dn_at=pawtab(itypat)%dnvdq0 1343 else 1344 dn_at=psps%nctab(itypat)%dnvdq0 1345 end if 1346 else if (optn2==3) then 1347 dn_at=-two*gauss1*alf2pi2 1348 end if 1349 else 1350 if (optn2==1) then 1351 ee=tcorespl(jj+1,1)-tcorespl(jj,1) 1352 ff=(3._dp*bb**2-1._dp)*tcorespl(jj+1,2) & 1353 & -(3._dp*aa**2-1._dp)*tcorespl(jj,2) 1354 else if (optn2==2) then 1355 ee=tvalespl(jj+1,1)-tvalespl(jj,1) 1356 ff=(3._dp*bb**2-1._dp)*tvalespl(jj+1,2) & 1357 & -(3._dp*aa**2-1._dp)*tvalespl(jj,2) 1358 else if (optn2==3) then 1359 dn_at=-two*gauss1*alf2pi2*exp(-gsquar*alf2pi2) 1360 else 1361 end if 1362 dn_at=(ee*dqm1+ff*dqdiv6)/gmag 1363 end if 1364 end if ! end optn 1365 1366 do id=1,ndir 1367 sfr=zero;sfi=zero 1368 do ia=ia1,ia2 1369 if (ipert==natom+3.or.ipert==natom+4) then 1370 ! sum[Exp(-i.2pi.g.xred)] 1371 sfr=sfr+phre_igia(ia) 1372 sfi=sfi-phim_igia(ia) 1373 else 1374 ! Exp(-i.2pi.g.xred) * -i.2pi.(g+q) 1375 sfr=-(two_pi*gq(jdir(id))*phim_igia(ia)) 1376 sfi=-(two_pi*gq(jdir(id))*phre_igia(ia)) 1377 end if 1378 end do 1379 1380 if (ipert/=natom+3.and.ipert/=natom+4) then 1381 ! Phonons case 1382 1383 ! Exp(-i.2pi.q.xred) => -i.2pi.(g+q).Exp(-i.2pi.(g+q).xred) 1384 sfqr= sfr*phqre+sfi*phqim 1385 sfqi=-sfr*phqim+sfi*phqre 1386 1387 if (optv == 1) then 1388 workv(re,ii,id) = workv(re,ii,id) + sfqr*v_at 1389 workv(im,ii,id) = workv(im,ii,id) + sfqi*v_at 1390 end if 1391 if (optn == 1) then 1392 workn(re,ii,id) = workn(re,ii,id) + sfqr*n_at 1393 workn(im,ii,id) = workn(im,ii,id) + sfqi*n_at 1394 end if 1395 1396 else 1397 ! Strain case 1398 1399 ! Compute G in cartesian coordinates 1400 gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+& 1401 & gprimd(1,3)*dble(ig3) 1402 gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+& 1403 & gprimd(2,3)*dble(ig3) 1404 gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+& 1405 & gprimd(3,3)*dble(ig3) 1406 1407 ! Accumulate -dV^AT/dG*rho(G)*SF(G)*Gi.Gj/G 1408 ! or -dn^AT/dG*V(G)*SF(G)*Gi.Gj/G 1409 if (optv==1) then 1410 if(jdir(id)<=3) then 1411 term_v = dv_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) + v_at 1412 else 1413 term_v = dv_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) 1414 end if 1415 workv(re,ii,id) = workv(re,ii,id) - (sfr*term_v) 1416 workv(im,ii,id) = workv(im,ii,id) - (sfi*term_v) 1417 end if 1418 1419 if (optn==1) then 1420 if(jdir(id)<=3) then 1421 term_n = dn_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) + n_at 1422 else 1423 term_n = dn_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) 1424 end if 1425 1426 workn(re,ii,id) = workn(re,ii,id) - (sfr*term_n) 1427 workn(im,ii,id) = workn(im,ii,id) - (sfi*term_n) 1428 end if 1429 1430 end if 1431 ! End loop on ndir 1432 1433 end do 1434 ! End skip G**2 outside cutoff 1435 end if 1436 ! End loop on n1, n2, n3 1437 end do 1438 end if ! this plane is selected 1439 end do 1440 end do 1441 ia1=ia2+1 1442 end do ! end loop itype 1443 1444 ABI_FREE(phre_igia) 1445 ABI_FREE(phim_igia) 1446 1447 if(ipert==natom+3.or.ipert==natom+4) then 1448 ! Set Vloc(G=0)=0: 1449 if (optv==1) then 1450 workv(re,1,:)=zero 1451 workv(im,1,:)=zero 1452 end if 1453 end if 1454 1455 ! Identify unbalanced g-vectors 1456 if (qeq05) then !This doesn't work in parallel 1457 ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2 1458 ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2 1459 ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2 1460 if (abs(abs(qphon(1))-half)<tol12) then 1461 if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max) 1462 if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min) 1463 end if 1464 if (abs(abs(qphon(2))-half)<tol12) then 1465 if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max) 1466 if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min) 1467 end if 1468 if (abs(abs(qphon(3))-half)<tol12) then 1469 if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max) 1470 if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min) 1471 end if 1472 end if 1473 1474 1475 ! Get 1st-order potential/density back to real space 1476 ! Non-symetrized non-zero elements have to be nullified 1477 ! Divide by unit cell volume 1478 1479 if (optv==1.or.optn==1) then 1480 1481 xnorm=one/ucvol 1482 ! Create fake mpi_enreg to wrap fourdp 1483 call initmpi_seq(mpi_enreg_fft) 1484 ABI_FREE(mpi_enreg_fft%distribfft) 1485 if (present(comm_fft)) then 1486 call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb) 1487 my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb 1488 else 1489 my_comm_fft=xmpi_comm_self;paral_kgb_fft=0; 1490 mpi_enreg_fft%distribfft => my_distribfft 1491 end if 1492 1493 if (optv==1) then 1494 do id=1,ndir 1495 ! Eliminate unbalanced g-vectors 1496 if (qeq0) then !q=0 1497 call zerosym(workv(:,:,id),2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft) 1498 else if (qeq05) then !q=1/2; this doesn't work in parallel 1499 call zerosym(workv(:,:,id),2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3) 1500 end if 1501 call fourdp(cplex,workv(:,:,id),atmvlocr1(:,id),1,mpi_enreg_fft,nfft,1,ngfft,0) 1502 atmvlocr1(:,id)=atmvlocr1(:,id)*xnorm 1503 end do 1504 1505 !if (present(atmvlocg1)) atmvlocg1 = workv 1506 ABI_FREE(workv) 1507 end if 1508 1509 if (optn==1) then 1510 do id=1,ndir 1511 ! Eliminate unbalanced g-vectors 1512 if (qeq0) then !q=0 1513 call zerosym(workn(:,:,id),2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft) 1514 else if (qeq05) then !q=1/2; this doesn't work in parallel 1515 call zerosym(workn(:,:,id),2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3) 1516 end if 1517 call fourdp(cplex,workn(:,:,id),atmrhor1(:,id),1,mpi_enreg_fft,nfft,1,ngfft,0) 1518 atmrhor1(:,id)=atmrhor1(:,id)*xnorm 1519 end do 1520 !if (present(atmrhog1)) atmrhog1 = workn 1521 ABI_FREE(workn) 1522 end if 1523 1524 ! Destroy fake mpi_enreg 1525 call unset_mpi_enreg_fft(mpi_enreg_fft) 1526 end if 1527 1528 if (.not.present(distribfft)) then 1529 call destroy_distribfft(my_distribfft) 1530 ABI_FREE(my_distribfft) 1531 end if 1532 1533 ! End the condition of non-electric-field 1534 end if 1535 1536 DBG_EXIT("COLL") 1537 1538 end subroutine dfpt_atm2fft
ABINIT/m_atm2fft [ Modules ]
NAME
m_atm2fft
FUNCTION
COPYRIGHT
Copyright (C) 1998-2022 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_atm2fft 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_xmpi 28 use m_distribfft 29 use m_dtset 30 31 use defs_abitypes, only : mpi_type 32 use m_time, only : timab 33 use defs_datatypes,only : pseudopotential_type 34 use m_gtermcutoff, only : termcutoff 35 use m_pawtab, only : pawtab_type 36 use m_fft, only : zerosym, fourdp 37 use m_mpinfo, only : set_mpi_enreg_fft, unset_mpi_enreg_fft, initmpi_seq 38 39 implicit none 40 41 private