TABLE OF CONTENTS
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
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (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
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
PARENTS
dfpt_dyxc1,dfpt_eltfrxc,dfpt_looppert,dfpt_nstpaw,pawgrnl
CHILDREN
destroy_distribfft,fourdp,init_distribfft_seq,initmpi_seq set_mpi_enreg_fft,unset_mpi_enreg_fft,zerosym
SOURCE
96 #if defined HAVE_CONFIG_H 97 #include "config.h" 98 #endif 99 100 #include "abi_common.h" 101 102 subroutine dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,& 103 & mgfft,mqgrid,natom,ndir,nfft,ngfft,ntypat,& 104 & ph1d,qgrid,qphon,typat,ucvol,usepaw,xred,psps,pawtab,& 105 & atmrhor1,atmrhog1,atmvlocr1,atmvlocg1,distribfft,gauss,comm_fft,me_g0,optn_in,& 106 & optn2_in,optv_in,paral_kgb,vspl) ! optional arguments 107 108 use defs_basis 109 use defs_abitypes 110 use m_profiling_abi 111 use m_errors 112 113 use m_xmpi, only : xmpi_comm_self, xmpi_comm_size 114 use defs_datatypes, only : pseudopotential_type 115 use m_pawtab, only : pawtab_type 116 use m_distribfft, only : distribfft_type 117 use m_fft, only : zerosym 118 use m_mpinfo, only : set_mpi_enreg_fft, unset_mpi_enreg_fft 119 120 !This section has been created automatically by the script Abilint (TD). 121 !Do not modify the following lines by hand. 122 #undef ABI_FUNC 123 #define ABI_FUNC 'dfpt_atm2fft' 124 use interfaces_51_manage_mpi 125 use interfaces_53_ffts 126 !End of the abilint section 127 128 implicit none 129 130 !Arguments ------------------------------------ 131 !scalars 132 integer,intent(in) :: cplex,idir,ipert,mgfft,mqgrid,natom,ndir,nfft,ntypat,usepaw 133 integer,optional,intent(in) :: optn_in,optn2_in,optv_in 134 integer,optional,intent(in) :: me_g0,comm_fft,paral_kgb 135 real(dp),intent(in) :: gsqcut,ucvol 136 type(pseudopotential_type),target,intent(in) :: psps 137 type(distribfft_type),optional,intent(in),target :: distribfft 138 !arrays 139 integer,intent(in) :: atindx(natom),ngfft(18),typat(natom) 140 real(dp),intent(in) :: gmet(3,3),gprimd(3,3) 141 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid),qphon(3) 142 real(dp),intent(in) :: xred(3,natom) 143 real(dp),optional,intent(in) :: gauss(2,ntypat),vspl(mqgrid,2,ntypat) 144 real(dp),optional,intent(out) :: atmrhor1(cplex*nfft,ndir) 145 real(dp),optional,intent(out) :: atmrhog1(2,nfft,ndir) 146 real(dp),optional,intent(out) :: atmvlocr1(cplex*nfft,ndir) 147 real(dp),optional,intent(out) :: atmvlocg1(2,nfft,ndir) 148 type(pawtab_type),target,intent(in) :: pawtab(ntypat*usepaw) 149 150 !Local variables ------------------------------ 151 !scalars 152 integer,parameter :: im=2,re=1 153 integer :: i1,i2,i3,ia,ia1,ia2,iatm,iatom,id,id1,id2,id3 154 integer :: ig1,ig1max,ig1min,ig2,ig2max,ig2min,ig3,ig3max,ig3min 155 integer :: ii,itypat,jj,me_fft,my_comm_fft,n1,n2,n3,nattyp,nproc_fft,ntype,paral_kgb_fft 156 integer :: optn,optv,optn2,shift1,shift2,shift3,type1,type2 157 logical :: have_g0,qeq0,qeq05 158 real(dp),parameter :: tolfix=1.0000001_dp 159 real(dp) :: aa,alf2pi2,bb,cc,cutoff,dd,diff,dq,dq2div6,dqdiv6,dqm1,ee,ff 160 real(dp) :: gauss1,gauss2,gmag,gq1,gq2,gq3,gsquar,n_at,dn_at,ph12i,ph12r,ph1i 161 real(dp) :: ph1r,ph2i,ph2r,ph3i,ph3r,phqim,phqre,qxred2pi 162 real(dp) :: sfi,sfqi,sfqr,sfr,term_n,term_v,v_at,dv_at,xnorm 163 type(distribfft_type),pointer :: my_distribfft 164 type(mpi_type) :: mpi_enreg_fft 165 !arrays 166 integer :: eps1(6)=(/1,2,3,2,3,1/),eps2(6)=(/1,2,3,3,1,2/),jdir(ndir) 167 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:) 168 real(dp), ABI_CONTIGUOUS pointer :: tvalespl(:,:),tcorespl(:,:) 169 real(dp) :: gq(6),gcart(3) 170 real(dp),allocatable :: phim_igia(:),phre_igia(:),workn(:,:,:),workv(:,:,:) 171 172 !no_abirules 173 !Define G^2 based on G space metric gmet. 174 ! gsq(g1,g2,g3)=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+g3*g3*gmet(3,3) & 175 ! & +two*(g1*g2*gmet(1,2)+g2*g3*gmet(2,3)+g3*g1*gmet(3,1)) 176 ! ************************************************************************* 177 178 DBG_ENTER("COLL") 179 180 ! Check optional arguments 181 if (present(comm_fft)) then 182 if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then 183 MSG_BUG('Need paral_kgb and me_g0 with comm_fft !') 184 end if 185 end if 186 187 if (present(gauss))then 188 if (.not.present(optn2_in)) then 189 optn2 = 3 190 else 191 if(optn2_in/=3)then 192 MSG_BUG('optn2 must be set to 3!') 193 else 194 optn2 = optn2_in 195 end if 196 end if 197 end if 198 199 if (present(atmrhor1))then 200 if(.not.present(optn_in))then 201 optn = 1 202 else 203 optn = optn_in 204 end if 205 if (.not.present(optn2_in)) then 206 MSG_BUG('rho1 calculation need optn2 !') 207 else 208 optn2 = optn2_in 209 end if 210 else 211 optn = 0 212 optn2 = 0 213 end if 214 215 if (present(atmvlocr1))then 216 if(.not.present(optv_in))then 217 optv = 1 218 else 219 if(optv_in/=1)then 220 MSG_BUG('optv_in must be set to 1!') 221 else 222 optv = optv_in 223 end if 224 end if 225 if(.not.present(vspl))then 226 MSG_BUG('vloc1 calculation need vspl!') 227 end if 228 else 229 optv = 0 230 end if 231 232 if(ipert==natom+1.or.ipert==natom+2.or.ipert==natom+10.or.ipert==natom+11) then 233 234 ! (In case of d/dk or an electric/magnetic field) 235 if (optn==1) then 236 atmrhor1(1:cplex*nfft,1:ndir)=zero 237 if (present(atmrhog1)) atmrhog1 = zero 238 end if 239 if (optv==1) then 240 atmvlocr1(1:cplex*nfft,1:ndir)=zero 241 if (present(atmvlocg1)) atmvlocg1 = zero 242 end if 243 244 else 245 246 ! Useful quantities 247 if (ipert/=natom+3.and.ipert/=natom+4) then 248 iatom=ipert;iatm=atindx(iatom) 249 itypat=typat(iatom) 250 else 251 !sum of all (strain pertubation) 252 iatom = 1 253 iatm = 1 254 itypat = 1 255 end if 256 257 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 258 me_fft=ngfft(11) 259 nproc_fft=ngfft(10) 260 if (ndir==1) then 261 jdir(1)=idir 262 else 263 do id=1,ndir 264 jdir(id)=id 265 end do 266 end if 267 268 ! Get the distrib associated with this fft_grid 269 if (present(distribfft)) then 270 my_distribfft => distribfft 271 else 272 ABI_DATATYPE_ALLOCATE(my_distribfft,) 273 call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp') 274 end if 275 if (n2==my_distribfft%n2_coarse) then 276 fftn2_distrib => my_distribfft%tab_fftdp2_distrib 277 else if (n2 == my_distribfft%n2_fine) then 278 fftn2_distrib => my_distribfft%tab_fftdp2dg_distrib 279 else 280 MSG_BUG("Unable to find an allocated distrib for this fft grid") 281 end if 282 283 qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) 284 qeq05=(abs(abs(qphon(1))-half)<tol12.or. & 285 & abs(abs(qphon(2))-half)<tol12.or. & 286 & abs(abs(qphon(3))-half)<tol12) 287 288 if (nproc_fft>1.and.qeq05) then 289 MSG_ERROR('not compatible with FFT parallelism') 290 end if 291 292 if (optn2==3)then 293 gauss1=gauss(1,itypat) 294 gauss2=gauss(2,itypat) 295 alf2pi2=(two_pi*gauss2)**2 296 end if 297 298 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 299 dqm1=one/dq 300 dqdiv6=dq/six 301 dq2div6=dq**2/six 302 cutoff=gsqcut*tolfix 303 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 304 ig1max=-1;ig2max=-1;ig3max=-1 305 ig1min=n1;ig2min=n2;ig3min=n3 306 307 ! Determination of phase qxred* 308 qxred2pi=two_pi*(qphon(1)*xred(1,iatom)+ & 309 & qphon(2)*xred(2,iatom)+ & 310 & qphon(3)*xred(3,iatom) ) 311 phqre=cos(qxred2pi) 312 phqim=sin(qxred2pi) 313 314 ! Zero out temporary arrays 315 if (optv==1) then 316 ABI_ALLOCATE(workv,(2,nfft,ndir)) 317 workv(:,:,:)=zero 318 end if 319 if (optn==1) then 320 ABI_ALLOCATE(workn,(2,nfft,ndir)) 321 workn(:,:,:)=zero 322 end if 323 324 if (ipert==natom+3.or.ipert==natom+4) then 325 ntype = ntypat 326 ia1=1 327 type1 = 1 328 type2 = ntype 329 ABI_ALLOCATE(phre_igia,(natom)) 330 ABI_ALLOCATE(phim_igia,(natom)) 331 else 332 type1 = itypat 333 type2 = itypat 334 ntype = 1 335 ABI_ALLOCATE(phre_igia,(iatm:iatm)) 336 ABI_ALLOCATE(phim_igia,(iatm:iatm)) 337 end if 338 339 340 do itypat=type1,type2 341 ! ia1,ia2 sets range of loop over atoms: 342 if (ipert==natom+3.or.ipert==natom+4) then 343 nattyp = count(typat(:)==itypat) 344 ia2=ia1+nattyp-1 345 else 346 ia1 = iatm 347 ia2 = iatm 348 end if 349 350 if (usepaw == 1) then 351 tcorespl => pawtab(itypat)%tcorespl 352 tvalespl => pawtab(itypat)%tvalespl 353 else 354 tcorespl => psps%nctab(itypat)%tcorespl 355 tvalespl => psps%nctab(itypat)%tvalespl 356 end if 357 358 ii=0 359 do i3=1,n3 360 ig3=i3-(i3/id3)*n3-1 361 gq3=dble(ig3)+qphon(3) 362 gq(3)=gq3 363 364 do i2=1,n2 365 if (fftn2_distrib(i2)==me_fft) then 366 ig2=i2-(i2/id2)*n2-1 367 gq2=dble(ig2)+qphon(2) 368 gq(2)=gq2 369 370 do i1=1,n1 371 ig1=i1-(i1/id1)*n1-1 372 gq1=dble(ig1)+qphon(1) 373 gq(1)=gq1 374 375 ii=ii+1 376 ! gsquar=gsq(gq1,gq2,gq3) 377 gsquar=gq1*gq1*gmet(1,1)+gq2*gq2*gmet(2,2)+gq3*gq3*gmet(3,3) & 378 & +two*(gq1*gq2*gmet(1,2)+gq2*gq3*gmet(2,3)+gq3*gq1*gmet(3,1)) 379 380 381 ! Skip G**2 outside cutoff: 382 if (gsquar<=cutoff) then 383 384 ! Identify min/max indexes (to cancel unbalanced contributions later) 385 if (qeq05) then 386 ig1max=max(ig1max,ig1);ig1min=min(ig1min,ig1) 387 ig2max=max(ig2max,ig2);ig2min=min(ig2min,ig2) 388 ig3max=max(ig3max,ig3);ig3min=min(ig3min,ig3) 389 end if 390 391 gmag=sqrt(gsquar) 392 have_g0=(ig1==0.and.ig2==0.and.ig3==0.and.qeq0) 393 394 jj=1+int(gmag*dqm1) 395 diff=gmag-qgrid(jj) 396 397 ! Compute structure factor 398 phre_igia(:) = zero 399 phim_igia(:) = zero 400 401 do ia=ia1,ia2 402 shift1=1+n1+(ia-1)*(2*n1+1) 403 shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1) 404 shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1) 405 ph1r=ph1d(re,ig1+shift1);ph1i=ph1d(im,ig1+shift1) 406 ph2r=ph1d(re,ig2+shift2);ph2i=ph1d(im,ig2+shift2) 407 ph3r=ph1d(re,ig3+shift3);ph3i=ph1d(im,ig3+shift3) 408 ph12r=ph1r*ph2r-ph1i*ph2i 409 ph12i=ph1r*ph2i+ph1i*ph2r 410 phre_igia(ia)=ph12r*ph3r-ph12i*ph3i 411 phim_igia(ia)=ph12r*ph3i+ph12i*ph3r 412 end do 413 ! Compute V^AT(g+q) and/or n^AT(g+q) for given type of atom 414 ! Evaluate spline fit: p. 86 Numerical Recipes, Press et al; 415 ! Note the error in book for sign of "aa" term in derivative. 416 if (optv==1.or.optn2/=3) then 417 bb = diff*dqm1 418 aa = one-bb 419 cc = aa*(aa**2-one)*dq2div6 420 dd = bb*(bb**2-one)*dq2div6 421 end if 422 423 if (optv==1) then 424 if (have_g0) then 425 v_at=zero 426 dv_at=zero 427 else 428 v_at=(aa*vspl(jj,1,itypat)+bb*vspl(jj+1,1,itypat)+& 429 & cc*vspl(jj,2,itypat)+dd*vspl(jj+1,2,itypat))/gsquar 430 431 ! Also get (dV(q)/dq)/q: 432 ! (note correction of Numerical Recipes sign error 433 ! before (3._dp*aa**2-1._dp) 434 ee= vspl(jj+1,1,itypat)-vspl(jj,1,itypat) 435 ff= (3._dp*bb**2-1._dp)*vspl(jj+1,2,itypat) & 436 & - (3._dp*aa**2-1._dp)*vspl(jj,2,itypat) 437 dv_at = ( ( ee*dqm1 + ff*dqdiv6 )/gmag& 438 & - 2.0_dp*v_at) / gsquar 439 end if 440 end if ! end optv 441 442 if (optn==1) then 443 if (optn2==1) then 444 n_at=(aa*tcorespl(jj,1)+bb*tcorespl(jj+1,1)+cc*tcorespl(jj,2)+dd*tcorespl(jj+1,2)) 445 else if (optn2==2) then 446 n_at=(aa*tvalespl(jj,1)+bb*tvalespl(jj+1,1)+cc*tvalespl(jj,2)+dd*tvalespl(jj+1,2)) 447 else if (optn2==3) then 448 n_at=gauss1*exp(-gsquar*alf2pi2) 449 else 450 n_at=zero 451 end if 452 453 ! Also get (dn^AT(q)/dq)/q: 454 if (have_g0) then 455 if (optn2==1) then 456 if (usepaw == 1) then 457 dn_at=pawtab(itypat)%dncdq0 458 else 459 dn_at=psps%nctab(itypat)%dncdq0 460 end if 461 else if (optn2==2) then 462 if (usepaw == 1) then 463 dn_at=pawtab(itypat)%dnvdq0 464 else 465 dn_at=psps%nctab(itypat)%dnvdq0 466 end if 467 else if (optn2==3) then 468 dn_at=-two*gauss1*alf2pi2 469 end if 470 else 471 if (optn2==1) then 472 ee=tcorespl(jj+1,1)-tcorespl(jj,1) 473 ff=(3._dp*bb**2-1._dp)*tcorespl(jj+1,2) & 474 & -(3._dp*aa**2-1._dp)*tcorespl(jj,2) 475 else if (optn2==2) then 476 ee=tvalespl(jj+1,1)-tvalespl(jj,1) 477 ff=(3._dp*bb**2-1._dp)*tvalespl(jj+1,2) & 478 & -(3._dp*aa**2-1._dp)*tvalespl(jj,2) 479 else if (optn2==3) then 480 dn_at=-two*gauss1*alf2pi2*exp(-gsquar*alf2pi2) 481 else 482 end if 483 dn_at=(ee*dqm1+ff*dqdiv6)/gmag 484 end if 485 end if ! end optn 486 487 do id=1,ndir 488 sfr=zero;sfi=zero 489 do ia=ia1,ia2 490 if (ipert==natom+3.or.ipert==natom+4) then 491 ! sum[Exp(-i.2pi.g.xred)] 492 sfr=sfr+phre_igia(ia) 493 sfi=sfi-phim_igia(ia) 494 else 495 ! Exp(-i.2pi.g.xred) * -i.2pi.(g+q) 496 sfr=-(two_pi*gq(jdir(id))*phim_igia(ia)) 497 sfi=-(two_pi*gq(jdir(id))*phre_igia(ia)) 498 end if 499 end do 500 501 if (ipert/=natom+3.and.ipert/=natom+4) then 502 ! Phonons case 503 504 ! Exp(-i.2pi.q.xred) => -i.2pi.(g+q).Exp(-i.2pi.(g+q).xred) 505 sfqr= sfr*phqre+sfi*phqim 506 sfqi=-sfr*phqim+sfi*phqre 507 508 if (optv == 1) then 509 workv(re,ii,id) = workv(re,ii,id) + sfqr*v_at 510 workv(im,ii,id) = workv(im,ii,id) + sfqi*v_at 511 end if 512 if (optn == 1) then 513 workn(re,ii,id) = workn(re,ii,id) + sfqr*n_at 514 workn(im,ii,id) = workn(im,ii,id) + sfqi*n_at 515 end if 516 517 else 518 ! Strain case 519 520 ! Compute G in cartesian coordinates 521 gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+& 522 & gprimd(1,3)*dble(ig3) 523 gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+& 524 & gprimd(2,3)*dble(ig3) 525 gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+& 526 & gprimd(3,3)*dble(ig3) 527 528 ! Accumulate -dV^AT/dG*rho(G)*SF(G)*Gi.Gj/G 529 ! or -dn^AT/dG*V(G)*SF(G)*Gi.Gj/G 530 if (optv==1) then 531 if(jdir(id)<=3) then 532 term_v = dv_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) + v_at 533 else 534 term_v = dv_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) 535 end if 536 workv(re,ii,id) = workv(re,ii,id) - (sfr*term_v) 537 workv(im,ii,id) = workv(im,ii,id) - (sfi*term_v) 538 end if 539 540 if (optn==1) then 541 if(jdir(id)<=3) then 542 term_n = dn_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) + n_at 543 else 544 term_n = dn_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) 545 end if 546 547 workn(re,ii,id) = workn(re,ii,id) - (sfr*term_n) 548 workn(im,ii,id) = workn(im,ii,id) - (sfi*term_n) 549 end if 550 551 end if 552 ! End loop on ndir 553 554 end do 555 ! End skip G**2 outside cutoff 556 end if 557 ! End loop on n1, n2, n3 558 end do 559 end if ! this plane is selected 560 end do 561 end do 562 ia1=ia2+1 563 end do ! end loop itype 564 565 ABI_DEALLOCATE(phre_igia) 566 ABI_DEALLOCATE(phim_igia) 567 568 if(ipert==natom+3.or.ipert==natom+4) then 569 ! Set Vloc(G=0)=0: 570 if (optv==1) then 571 workv(re,1,:)=zero 572 workv(im,1,:)=zero 573 end if 574 end if 575 576 ! Identify unbalanced g-vectors 577 if (qeq05) then !This doesn't work in parallel 578 ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2 579 ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2 580 ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2 581 if (abs(abs(qphon(1))-half)<tol12) then 582 if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max) 583 if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min) 584 end if 585 if (abs(abs(qphon(2))-half)<tol12) then 586 if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max) 587 if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min) 588 end if 589 if (abs(abs(qphon(3))-half)<tol12) then 590 if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max) 591 if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min) 592 end if 593 end if 594 595 596 ! Get 1st-order potential/density back to real space 597 ! Non-symetrized non-zero elements have to be nullified 598 ! Divide by unit cell volume 599 600 if (optv==1.or.optn==1) then 601 602 xnorm=one/ucvol 603 ! Create fake mpi_enreg to wrap fourdp 604 call initmpi_seq(mpi_enreg_fft) 605 ABI_DATATYPE_DEALLOCATE(mpi_enreg_fft%distribfft) 606 if (present(comm_fft)) then 607 call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb) 608 my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb 609 else 610 my_comm_fft=xmpi_comm_self;paral_kgb_fft=0; 611 mpi_enreg_fft%distribfft => my_distribfft 612 end if 613 614 if (optv==1) then 615 do id=1,ndir 616 ! Eliminate unbalanced g-vectors 617 if (qeq0) then !q=0 618 call zerosym(workv(:,:,id),2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft) 619 else if (qeq05) then !q=1/2; this doesn't work in parallel 620 call zerosym(workv(:,:,id),2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3) 621 end if 622 call fourdp(cplex,workv(:,:,id),atmvlocr1(:,id),1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0) 623 atmvlocr1(:,id)=atmvlocr1(:,id)*xnorm 624 end do 625 626 !if (present(atmvlocg1)) atmvlocg1 = workv 627 ABI_DEALLOCATE(workv) 628 end if 629 630 if (optn==1) then 631 do id=1,ndir 632 ! Eliminate unbalanced g-vectors 633 if (qeq0) then !q=0 634 call zerosym(workn(:,:,id),2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft) 635 else if (qeq05) then !q=1/2; this doesn't work in parallel 636 call zerosym(workn(:,:,id),2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3) 637 end if 638 call fourdp(cplex,workn(:,:,id),atmrhor1(:,id),1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0) 639 atmrhor1(:,id)=atmrhor1(:,id)*xnorm 640 end do 641 !if (present(atmrhog1)) atmrhog1 = workn 642 ABI_DEALLOCATE(workn) 643 end if 644 645 ! Destroy fake mpi_enreg 646 call unset_mpi_enreg_fft(mpi_enreg_fft) 647 end if 648 649 if (.not.present(distribfft)) then 650 call destroy_distribfft(my_distribfft) 651 ABI_DATATYPE_DEALLOCATE(my_distribfft) 652 end if 653 654 ! End the condition of non-electric-field 655 end if 656 657 DBG_EXIT("COLL") 658 659 end subroutine dfpt_atm2fft