TABLE OF CONTENTS
- ABINIT/mklocl_wavelets
- mklocl_wavelets/calcdVloc_wvl
- mklocl_wavelets/dvloc_zero_wvl
- mklocl_wavelets/local_forces_wvl
ABINIT/mklocl_wavelets [ Functions ]
NAME
mklocl_wavelets
FUNCTION
Compute the ionic local potential when the pseudo-potentials are GTH, using the special decomposition of these pseudo. The resulting potential is computed with free boundary conditions. It gives the same result than mklocl_realspace for the GTH pseudo only with a different way to compute the potential. Optionally compute : option=1 : local ionic potential throughout unit cell option=2 : contribution of local ionic potential to E gradient wrt xred
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DC,TRangel,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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
efield (3)=external electric field mpi_enreg=informations about MPI parallelization natom=number of atoms nfft=size of vpsp (local potential) nspden=number of spin-density components option=type of calculation (potential, forces, ...) rprimd(3,3)=dimensional primitive translations in real space (bohr) xcart(3,natom)=cartesian atomic coordinates. wvl_den=density-potential BigDFT object wvl_descr=wavelet BigDFT object
OUTPUT
(if option==1) vpsp(nfft)=the potential resulting from the ionic
density of charge.
(if option==2) grtn(3,natom)=grads of Etot wrt tn. These gradients are in
reduced coordinates. Multiply them by rprimd to get
gradients in cartesian coordinates.
PARENTS
mklocl,wvl_wfsinp_scratch
CHILDREN
calcdvloc_wvl,derf_ab,paw_splint_der
SOURCE
51 #if defined HAVE_CONFIG_H 52 #include "config.h" 53 #endif 54 55 #include "abi_common.h" 56 57 subroutine mklocl_wavelets(efield, grtn, mpi_enreg, natom, nfft, & 58 & nspden, option, rprimd, vpsp, wvl_den, wvl_descr, xcart) 59 60 use defs_basis 61 use defs_abitypes 62 use defs_wvltypes 63 use m_abi2big, only : wvl_rhov_abi2big 64 use m_profiling_abi 65 use m_xmpi 66 use m_errors 67 #if defined HAVE_BIGDFT 68 use BigDFT_API, only : ELECTRONIC_DENSITY,createIonicPotential,local_forces 69 use poisson_solver, only : H_potential 70 #endif 71 72 !This section has been created automatically by the script Abilint (TD). 73 !Do not modify the following lines by hand. 74 #undef ABI_FUNC 75 #define ABI_FUNC 'mklocl_wavelets' 76 use interfaces_14_hidewrite 77 use interfaces_67_common, except_this_one => mklocl_wavelets 78 !End of the abilint section 79 80 implicit none 81 82 !Arguments ------------------------------------ 83 !scalars 84 integer,intent(in) :: option, natom, nfft, nspden 85 type(MPI_type),intent(in) :: mpi_enreg 86 type(wvl_denspot_type), intent(inout) :: wvl_den 87 type(wvl_internal_type), intent(in) :: wvl_descr 88 !arrays 89 real(dp),intent(in) :: rprimd(3,3),efield(3) 90 real(dp),intent(inout) :: grtn(3,natom) 91 real(dp), intent(inout) :: vpsp(nfft) 92 real(dp),intent(inout) :: xcart(3,natom) 93 94 !Local variables------------------------------- 95 #if defined HAVE_BIGDFT 96 !scalars 97 integer :: i,i1,i2,i3,ia,ierr,igeo,me,nproc,shift,comm 98 real(dp) :: energ 99 character(len=500) :: message 100 !arrays 101 real(dp) :: epot(3) 102 real(dp) :: elecfield(3)=(/zero,zero,zero/) ! Not used here 103 real(dp),allocatable :: gxyz(:,:),vhartr(:),rhov(:,:) 104 #endif 105 106 ! ********************************************************************* 107 108 #if defined HAVE_BIGDFT 109 110 elecfield=zero !not used here 111 112 !Manage parallelism 113 comm=mpi_enreg%comm_wvl 114 nproc=xmpi_comm_size(comm) 115 me=xmpi_comm_rank(comm) 116 117 !---------------------------------------------------------------------- 118 ! ----- Option 1: compute local ionic potential ----- 119 !---------------------------------------------------------------------- 120 if (option == 1) then 121 122 ! We get the kernel for the Poisson solver (used to go from the ionic 123 ! charge to the potential).If the kernel is uncomputed, it does it now. 124 ! call psolver_kernel(wvl_den%denspot%dpbox%hgrids, 2, icoulomb, me, kernel, & 125 !& comm, wvl_den%denspot%dpbox%ndims, nproc, nscforder) 126 ! if (.not.associated(kernel%co%kernel)) then 127 ! call psolver_kernel(wvl_den%denspot%dpbox%hgrids, 1, icoulomb, me, kernel, & 128 !& comm, wvl_den%denspot%dpbox%ndims, nproc, nscforder) 129 ! end if 130 131 message=ch10//' mklocl_wavelets: Create local potential from ions.' 132 call wrtout(std_out,message,'COLL') 133 134 shift = 1 + wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) & 135 & * wvl_den%denspot%dpbox%i3xcsh 136 137 ! Call the BigDFT routine 138 call createIonicPotential(wvl_descr%atoms%astruct%geocode, me, nproc, (me == 0), wvl_descr%atoms, & 139 & xcart, wvl_den%denspot%dpbox%hgrids(1), wvl_den%denspot%dpbox%hgrids(2), & 140 & wvl_den%denspot%dpbox%hgrids(3), & 141 & elecfield, wvl_descr%Glr%d%n1, wvl_descr%Glr%d%n2, wvl_descr%Glr%d%n3, & 142 & wvl_den%denspot%dpbox%n3pi, wvl_den%denspot%dpbox%i3s + wvl_den%denspot%dpbox%i3xcsh, & 143 & wvl_den%denspot%dpbox%ndims(1), wvl_den%denspot%dpbox%ndims(2), & 144 & wvl_den%denspot%dpbox%ndims(3), wvl_den%denspot%pkernel, vpsp(shift), 0.d0,wvl_descr%rholoc) 145 146 ! Copy vpsp into bigdft object: 147 call wvl_rhov_abi2big(1,vpsp,wvl_den%denspot%v_ext,shift=(shift-1)) 148 149 ! Eventually add the electric field 150 if (maxval(efield) > tol12) then 151 message=ch10//'mklocl_wavelets: Add the electric field.' 152 call wrtout(std_out,message,'COLL') 153 ! We add here the electric field since in BigDFT, the field must be on x... 154 epot(:) = real(0.5, dp) * efield(:) * wvl_den%denspot%dpbox%hgrids(:) 155 do i3 = 1, wvl_den%denspot%dpbox%n3pi, 1 156 ia = (i3 - 1) * wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) 157 do i2 = -14, 2 * wvl_descr%Glr%d%n2 + 16, 1 158 i = ia + (i2 + 14) * wvl_den%denspot%dpbox%ndims(1) 159 do i1 = -14, 2 * wvl_descr%Glr%d%n1 + 16, 1 160 i = i + 1 161 vpsp(shift + i) = vpsp(shift + i) + & 162 & epot(1) * real(i1 - wvl_descr%Glr%d%n1, dp) + & 163 & epot(2) * real(i2 - wvl_descr%Glr%d%n2, dp) + & 164 & epot(3) * real(i3 - wvl_descr%Glr%d%n3, dp) 165 end do 166 end do 167 end do 168 end if 169 170 !---------------------------------------------------------------------- 171 ! ----- Option 2: compute forces induced by local ionic potential ----- 172 !---------------------------------------------------------------------- 173 else if (option == 2) then 174 175 message=ch10//' mklocl_wavelets: compute local forces.' 176 call wrtout(std_out,message,'COLL') 177 178 if (wvl_den%denspot%rhov_is/=ELECTRONIC_DENSITY) then 179 message='denspot bigdft datstructure should contain rhor!' 180 MSG_BUG(message) 181 end if 182 183 ! Extract density rhor from bigDFT datastructure 184 ABI_ALLOCATE(rhov,(nfft, nspden)) 185 ABI_ALLOCATE(vhartr,(nfft)) 186 shift = wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) & 187 & * wvl_den%denspot%dpbox%i3xcsh 188 do i = 1, nfft 189 rhov(i, 1) = wvl_den%denspot%rhov(i + shift) 190 vhartr(i) = wvl_den%denspot%rhov(i + shift) 191 end do 192 if (nspden == 2) then 193 shift = shift + wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) & 194 & * wvl_den%denspot%dpbox%n3d 195 do i = 1, nfft 196 rhov(i, 2) = wvl_den%denspot%rhov(i + shift) 197 vhartr(i) = vhartr(i) + wvl_den%denspot%rhov(i + shift) 198 end do 199 end if 200 201 ! Compute Hartree potential from rhor 202 call H_potential('D',wvl_den%denspot%pkernel,vhartr,vhartr,energ,zero,.false.) 203 204 ! Allocate temporary array for forces 205 ABI_ALLOCATE(gxyz,(3, natom)) 206 207 ! Calculate local part of the forces grtn (modified BigDFT routine) 208 call local_forces_wvl(me,natom,xcart,& 209 & wvl_den%denspot%dpbox%hgrids(1),& 210 & wvl_den%denspot%dpbox%hgrids(2),& 211 & wvl_den%denspot%dpbox%hgrids(3),& 212 & wvl_descr%Glr%d%n1,wvl_descr%Glr%d%n2,wvl_descr%Glr%d%n3,& 213 & wvl_den%denspot%dpbox%n3p,& 214 & wvl_den%denspot%dpbox%i3s+wvl_den%denspot%dpbox%i3xcsh,& 215 & wvl_den%denspot%dpbox%ndims(1),wvl_den%denspot%dpbox%ndims(2),& 216 & rhov,vhartr,gxyz,wvl_descr) 217 ! call local_forces(me, wvl_descr%atoms, xcart, & 218 ! & wvl_den%denspot%dpbox%hgrids(1), wvl_den%denspot%dpbox%hgrids(2), & 219 ! & wvl_den%denspot%dpbox%hgrids(3), & 220 ! & wvl_descr%Glr%d%n1, wvl_descr%Glr%d%n2, wvl_descr%Glr%d%n3, & 221 ! & wvl_den%denspot%dpbox%n3p, wvl_den%denspot%dpbox%i3s + wvl_den%denspot%dpbox%i3xcsh, & 222 ! & wvl_den%denspot%dpbox%ndims(1), wvl_den%denspot%dpbox%ndims(2), & 223 ! & rhov, vhartr,gxyz,locstrten,charge) 224 225 ! Pending: floc,locstrten and charge are not used here. 226 ! Pending: put mpi_enreg%nscatterarr... in object denspot, initialize object, etc. 227 228 if (nproc > 1) then 229 call xmpi_sum(gxyz, comm, ierr) 230 end if 231 232 ! Forces should be in reduced coordinates. 233 do ia = 1, natom, 1 234 do igeo = 1, 3, 1 235 grtn(igeo, ia) = - rprimd(1, igeo) * gxyz(1, ia) & 236 & - rprimd(2, igeo) * gxyz(2, ia) & 237 & - rprimd(3, igeo) * gxyz(3, ia) 238 end do 239 end do 240 241 ! Deallocate local variables 242 ABI_DEALLOCATE(vhartr) 243 ABI_DEALLOCATE(rhov) 244 ABI_DEALLOCATE(gxyz) 245 246 !---------------------------------------------------------------------- 247 248 else ! option switch 249 message = 'Internal error, option should be 1 or 2!' 250 MSG_ERROR(message) 251 end if 252 253 #else 254 BIGDFT_NOTENABLED_ERROR() 255 if (.false.) write(std_out,*) option,natom,nfft,nspden,mpi_enreg%me,& 256 & wvl_den%symObj,wvl_descr%h(1),rprimd(1,1),efield(1),grtn(1,1),vpsp(1),xcart(1,1) 257 #endif 258 259 end subroutine mklocl_wavelets
mklocl_wavelets/calcdVloc_wvl [ Functions ]
[ Top ] [ mklocl_wavelets ] [ Functions ]
NAME
calcdVloc_wvl
FUNCTION
Compute 1st-derivative of long-range HGH local ionic potential (derf)
COPYRIGHT
Copyright (C) 2016-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
OUTPUT
PARENTS
mklocl_wavelets
CHILDREN
calcdvloc_wvl,derf_ab,paw_splint_der
SOURCE
533 subroutine calcdVloc_wvl(yy,xx,rloc,Z) 534 535 use defs_basis 536 537 !This section has been created automatically by the script Abilint (TD). 538 !Do not modify the following lines by hand. 539 #undef ABI_FUNC 540 #define ABI_FUNC 'calcdVloc_wvl' 541 use interfaces_43_wvl_wrappers 542 !End of the abilint section 543 544 implicit none 545 !Arguments ------------------------------------ 546 !scalars 547 real(dp),intent(in) :: xx,rloc,Z 548 real(dp),intent(out) :: yy 549 550 !Local variables------------------------------- 551 !scalars 552 real(dp):: arg,tt 553 554 ! ************************************************************************* 555 556 arg=xx/(sqrt(2._dp)*rloc) 557 call derf_ab(tt,arg) 558 yy=(Z/(xx**2))* ( tt - 2._dp/sqrt(pi)*arg*exp(-arg**2) ) 559 560 end subroutine calcdVloc_wvl
mklocl_wavelets/dvloc_zero_wvl [ Functions ]
[ Top ] [ mklocl_wavelets ] [ Functions ]
NAME
dvloc_zero_wvl
FUNCTION
Use a quadratic interpolation to get limit of (1/x).dVloc(x)/dx at x->0
COPYRIGHT
Copyright (C) 2013-2018 ABINIT group (TRangel,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
OUTPUT
PARENTS
SOURCE
586 function dvloc_zero_wvl(charge,rloc,msz,rad,vloc,d2vloc) 587 588 use defs_basis 589 590 !This section has been created automatically by the script Abilint (TD). 591 !Do not modify the following lines by hand. 592 #undef ABI_FUNC 593 #define ABI_FUNC 'dvloc_zero_wvl' 594 !End of the abilint section 595 596 implicit none 597 598 !Arguments ------------------------------------ 599 !scalars 600 integer,intent(in) :: msz 601 real(dp) :: dvloc_zero_wvl 602 real(dp),intent(in) :: charge,rloc 603 !arrays 604 real(dp) :: rad(msz),vloc(msz),d2vloc(msz) 605 606 !Local variables------------------------------- 607 !scalars 608 real(dp) :: y1,y2,y3,zz=0._dp 609 !arrays 610 real(dp) :: ll(3),xx(3),yy(3) 611 612 ! ************************************************************************* 613 614 !Select 3 points x1,x2,x3 near 0 615 if (rad(1)>1.d-10) then 616 xx(1:3)=rad(1:3) 617 else 618 xx(1:3)=rad(2:4) 619 end if 620 621 !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x 622 call paw_splint_der(msz,rad,vloc,d2vloc,3,xx,yy) 623 call calcdVloc_wvl(y1,xx(1),rloc,charge) 624 call calcdVloc_wvl(y2,xx(2),rloc,charge) 625 call calcdVloc_wvl(y3,xx(3),rloc,charge) 626 yy(1)=(yy(1)-y1)/xx(1) 627 yy(2)=(yy(2)-y2)/xx(2) 628 yy(3)=(yy(3)-y3)/xx(3) 629 630 !Find a polynomial of the form (z=0): 631 !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z) 632 633 !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3)) 634 ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3))) 635 !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3)) 636 ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3))) 637 !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2)) 638 ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2))) 639 640 dvloc_zero_wvl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3) 641 642 end function dvloc_zero_wvl
mklocl_wavelets/local_forces_wvl [ Functions ]
[ Top ] [ mklocl_wavelets ] [ Functions ]
NAME
local_forces_wvl
FUNCTION
INPUTS
hxh,hyh,hzh=wavelet grid spacings iproc=current MPI process number n1,n2,n3=number of wavelet points in each direction n1i,n2i=size of intermediate xy wvl grid n3pi=number of xy wvl planes handled by current MPI process i3s=starting index of local potential for current MPI process natom=number of atoms pot(*)=Hartree ionic potential rho(*)=electronic density rxyz(3,natom)=cartesian coordinates of atoms wvl=wavelet BigDFT object
OUTPUT
floc(3,natom)=local ionic potential contribution to forces
PARENTS
mklocl_wavelets
CHILDREN
calcdvloc_wvl,derf_ab,paw_splint_der
SOURCE
294 #if defined HAVE_CONFIG_H 295 #include "config.h" 296 #endif 297 298 #include "abi_common.h" 299 300 subroutine local_forces_wvl(iproc,natom,rxyz,hxh,hyh,hzh,n1,n2,n3,n3pi,i3s,n1i,n2i,& 301 & rho,pot,floc,wvl) 302 303 use defs_basis 304 use defs_wvltypes 305 use m_errors 306 use m_paw_numeric, only : paw_splint_der 307 #if defined HAVE_BIGDFT 308 use BigDFT_API, only : PSPCODE_PAW,ind_positions 309 #endif 310 311 !This section has been created automatically by the script Abilint (TD). 312 !Do not modify the following lines by hand. 313 #undef ABI_FUNC 314 #define ABI_FUNC 'local_forces_wvl' 315 !End of the abilint section 316 317 implicit none 318 319 !Arguments ------------------------------- 320 !scalars 321 integer,intent(in) :: i3s,iproc,n1,n1i,n2,n2i,n3,n3pi,natom 322 real(dp),intent(in) :: hxh,hyh,hzh 323 type(wvl_internal_type),intent(in) :: wvl 324 !arrays 325 real(dp),intent(in) :: rxyz(3,natom) 326 real(dp),dimension(*),intent(in) :: rho,pot 327 real(dp),intent(out) :: floc(3,natom) 328 329 !Local variables ------------------------- 330 #if defined HAVE_BIGDFT 331 !scalars 332 integer :: i1,i2,i3,iat,iex,iey,iez,iloc,ind,isx,isy,isz,ityp 333 integer :: j1,j2,j3,msz=0 334 integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3,nloc 335 logical :: perx,pery,perz,gox,goy,goz,USE_PAW 336 real(dp) :: arg,charge,cutoff,dvhgh,forceloc,fxerf,fyerf,fzerf,fxgau,fygau,fzgau 337 real(dp) :: forceleaked,prefactor,r2,rhoel,rloc,rloc2,rx,ry,rz,rzero,tt,vel,x,xp,y,z 338 !arrays 339 real(dp) :: cprime(4),dvpawdr(1),rr(1) 340 real(dp),pointer :: psppar(:,:),rad(:),vloc(:),d2vloc(:) 341 #endif 342 343 ! ********************************************************************* 344 345 #if defined HAVE_BIGDFT 346 347 if (iproc==0) write(std_out,'(a)',advance='no') ' Calculate local forces...' 348 349 !PAW or NCPP ? 350 USE_PAW=any(wvl%atoms%npspcode==PSPCODE_PAW) 351 352 !Conditions for periodicity in the three directions 353 perx=(wvl%atoms%astruct%geocode /= 'F') 354 pery=(wvl%atoms%astruct%geocode == 'P') 355 perz=(wvl%atoms%astruct%geocode /= 'F') 356 call ext_buffers(perx,nbl1,nbr1) 357 call ext_buffers(pery,nbl2,nbr2) 358 call ext_buffers(perz,nbl3,nbr3) 359 360 forceleaked=zero 361 362 do iat=1,natom 363 ityp=wvl%atoms%astruct%iatype(iat) 364 365 ! Coordinates of the center 366 rx=rxyz(1,iat) 367 ry=rxyz(2,iat) 368 rz=rxyz(3,iat) 369 370 ! Initialization of the forces 371 ! ion-electron term, error function part 372 fxerf=zero 373 fyerf=zero 374 fzerf=zero 375 ! ion-electron term, gaussian part 376 fxgau=zero 377 fygau=zero 378 fzgau=zero 379 380 ! Building array of coefficients of the derivative of the gaussian part 381 psppar => wvl%atoms%psppar(:,:,ityp) 382 cprime(1)=2._dp*wvl%atoms%psppar(0,2,ityp)-wvl%atoms%psppar(0,1,ityp) 383 cprime(2)=4._dp*wvl%atoms%psppar(0,3,ityp)-wvl%atoms%psppar(0,2,ityp) 384 cprime(3)=6._dp*wvl%atoms%psppar(0,4,ityp)-wvl%atoms%psppar(0,3,ityp) 385 cprime(4)=-wvl%atoms%psppar(0,4,ityp) 386 387 ! Determine number of local terms (HGH pot) 388 nloc=0 389 do iloc=1,4 390 if (wvl%atoms%psppar(0,iloc,ityp).ne.zero) nloc=iloc 391 end do 392 393 ! Some constants depending on the atom type 394 rloc=wvl%atoms%psppar(0,0,ityp) ; rloc2=rloc**2 395 charge=real(wvl%atoms%nelpsp(ityp),kind=dp) 396 prefactor=charge/(2._dp*pi*sqrt(2._dp*pi)*rloc**5) 397 398 ! PAW specifics 399 if (USE_PAW) then 400 msz=wvl%rholoc%msz(ityp) 401 rad => wvl%rholoc%rad(1:msz,ityp) 402 vloc => wvl%rholoc%d(1:msz,3,ityp) 403 d2vloc => wvl%rholoc%d(1:msz,4,ityp) 404 rzero=rad(1);if (rzero<=1.d-10) rzero=rad(2) 405 end if 406 407 ! Maximum extension of the gaussian 408 cutoff=10._dp*rloc 409 isx=floor((rx-cutoff)/hxh) 410 isy=floor((ry-cutoff)/hyh) 411 isz=floor((rz-cutoff)/hzh) 412 iex=ceiling((rx+cutoff)/hxh) 413 iey=ceiling((ry+cutoff)/hyh) 414 iez=ceiling((rz+cutoff)/hzh) 415 416 ! Calculate the forces near the atom due to the gaussian 417 ! and error function parts of the potential 418 if (n3pi>0) then 419 do i3=isz,iez 420 z=real(i3,kind=dp)*hzh-rz 421 call ind_positions(perz,i3,n3,j3,goz) 422 j3=j3+nbl3+1 423 do i2=isy,iey 424 y=real(i2,kind=dp)*hyh-ry 425 call ind_positions(pery,i2,n2,j2,goy) 426 do i1=isx,iex 427 x=real(i1,kind=dp)*hxh-rx 428 call ind_positions(perx,i1,n1,j1,gox) 429 430 r2=x**2+y**2+z**2 431 xp=exp(-0.5_dp*r2/rloc2) 432 433 if ((j3>=i3s.and.j3<=i3s+n3pi-1) .and. (gox.and.goy)) then 434 ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s+1-1)*n1i*n2i 435 436 ! Short range part 437 rhoel=rho(ind) 438 ! HGH: V_S^prime=gaussian 439 if (.not.USE_PAW) then 440 if (nloc/=0) then 441 arg=r2/rloc2 442 tt=cprime(nloc) 443 do iloc=nloc-1,1,-1 444 tt=arg*tt+cprime(iloc) 445 end do 446 forceloc=xp*tt*rhoel 447 else 448 forceloc=zero 449 end if 450 ! PAW: V_PAW^prime-V_L^prime 451 else 452 rr(1)=sqrt(r2) 453 454 if (rr(1)>=rzero) then 455 call paw_splint_der(msz,rad,vloc,d2vloc,1,rr,dvpawdr) 456 call calcdVloc_wvl(dvhgh,rr(1),rloc,charge) 457 forceloc=rhoel*rloc2*(dvpawdr(1)-dvhgh)/rr(1) 458 else 459 forceloc=rhoel*rloc2*dvloc_zero_wvl(charge,rloc,msz,rad,vloc,d2vloc) 460 end if 461 end if 462 463 fxgau=fxgau+forceloc*x 464 fygau=fygau+forceloc*y 465 fzgau=fzgau+forceloc*z 466 467 ! Long range part: error function 468 vel=pot(ind) 469 fxerf=fxerf+xp*vel*x 470 fyerf=fyerf+xp*vel*y 471 fzerf=fzerf+xp*vel*z 472 473 else if ((.not.goz).and.(nloc>0)) then 474 arg=r2/rloc2 475 tt=cprime(nloc) 476 do iloc=nloc-1,1,-1 477 tt=arg*tt+cprime(iloc) 478 end do 479 forceleaked=forceleaked+prefactor*xp*tt*rho(1) 480 end if 481 482 end do ! i1 483 end do ! i2 484 end do ! i3 485 end if ! n3pi>0 486 487 ! Final result of the forces 488 floc(1,iat)=(hxh*hyh*hzh*prefactor)*fxerf+(hxh*hyh*hzh/rloc2)*fxgau 489 floc(2,iat)=(hxh*hyh*hzh*prefactor)*fyerf+(hxh*hyh*hzh/rloc2)*fygau 490 floc(3,iat)=(hxh*hyh*hzh*prefactor)*fzerf+(hxh*hyh*hzh/rloc2)*fzgau 491 492 end do ! iat 493 494 forceleaked=forceleaked*hxh*hyh*hzh 495 if (iproc.eq.0) write(std_out,'(a,1pe12.5)') 'done. Leaked force: ',forceleaked 496 497 #else 498 BIGDFT_NOTENABLED_ERROR() 499 if (.false.) write(std_out,*) i3s,iproc,n1,n1i,n2,n2i,n3,n3pi,natom,hxh,hyh,hzh,& 500 & rxyz(1,1),floc(1,1),rho(1),pot(1),wvl%h(1) 501 #endif 502 503 CONTAINS