TABLE OF CONTENTS
- ABINIT/m_mklocl_realspace
- ABINIT/mklocl_realspace
- ABINIT/mklocl_wavelets
- mklocl_realspace/calcdVloc_mklocl
- mklocl_realspace/calcVloc_mklocl
- mklocl_realspace/createIonicPotential_new
- mklocl_realspace/ind_positions_mklocl
- mklocl_realspace/local_forces_new
- mklocl_realspace/vloc_zero_mklocl
- mklocl_wavelets/calcdVloc_wvl
- mklocl_wavelets/dvloc_zero_mklocl
- mklocl_wavelets/dvloc_zero_wvl
- mklocl_wavelets/local_forces_wvl
ABINIT/m_mklocl_realspace [ Modules ]
NAME
m_mklocl_realspace
FUNCTION
Routines related to the local part of the pseudopotentials. Computation is done in real space (useful for isolated systems).
COPYRIGHT
Copyright (C) 2013-2022 ABINIT group (TRangel, MT, DC) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
TODO
This module could be merged with m_mklocl
SOURCE
20 #if defined HAVE_CONFIG_H 21 #include "config.h" 22 #endif 23 24 #include "abi_common.h" 25 26 module m_mklocl_realspace 27 28 use defs_basis 29 use defs_wvltypes 30 use m_xmpi 31 use m_abicore 32 use m_errors 33 34 use defs_datatypes, only : pseudopotential_type 35 use defs_abitypes, only : MPI_type 36 use m_time, only : timab 37 use m_geometry, only : xred2xcart 38 use m_fft_mesh, only : mkgrid_fft 39 use m_mpinfo, only : ptabs_fourdp 40 use m_pawtab, only : pawtab_type 41 use m_paw_numeric, only : paw_splint, paw_splint_der 42 use m_psolver, only : psolver_hartree, psolver_kernel 43 use m_abi2big, only : wvl_rhov_abi2big 44 use m_wvl_wfs, only : derf_ab 45 use m_fft, only : fourdp 46 47 implicit none 48 49 private
ABINIT/mklocl_realspace [ Functions ]
NAME
mklocl_realspace
FUNCTION
This method is equivalent to mklocl_recipspace except that it uses real space pseudo-potentials. It is useful for isolated systems. Then the option 3 and 4 are not available for this implementation. Optionally compute: option=1 : local ionic potential throughout unit cell option=2 : contribution of local ionic potential to E gradient wrt xred
INPUTS
dtset <type(dataset_type)>=all input variables in this dataset mpi_enreg=information about MPI parallelization 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, see ~abinit/doc/variables/vargs.htm#ngfft nspden=number of spin-density components ntypat=number of types of atoms. option= (see above) pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials rhog(2,nfft)=electron density rho(G) (electrons/$\textrm{Bohr}^3$) rhor(nfft,nspden)=electron density in electrons/bohr**3. (needed if option==2 or if option==4) rprimd(3,3)=dimensional primitive translations in real space (bohr) ucvol=unit cell volume ($\textrm{Bohr}^3$). xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
(if option==1) vpsp(nfft)=local crystal pseudopotential in real space. (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.
SOURCE
99 subroutine mklocl_realspace(grtn,icoulomb,mpi_enreg,natom,nattyp,nfft,ngfft,nscforder, & 100 & nspden,ntypat,option,pawtab,psps,rhog,rhor,rprimd,typat,& 101 & ucvol,usewvl,vpsp,xred) 102 103 #if defined HAVE_BIGDFT 104 use BigDFT_API, only : coulomb_operator,deallocate_coulomb_operator 105 use defs_PSolver 106 #else 107 use defs_wvltypes, only : coulomb_operator 108 #endif 109 110 !Arguments ------------------------------------ 111 !scalars 112 integer,intent(in) :: natom,nfft,nspden,ntypat,option 113 real(dp),intent(in) :: ucvol 114 type(MPI_type),intent(in) :: mpi_enreg 115 type(pseudopotential_type),intent(in) :: psps 116 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 117 !arrays 118 integer,intent(in) :: icoulomb,nscforder,usewvl 119 integer,intent(in) :: nattyp(ntypat),ngfft(18),typat(natom) 120 real(dp),intent(in) :: rhog(2,nfft) 121 real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3) 122 real(dp),intent(in) :: xred(3,natom) 123 real(dp),intent(out) :: grtn(3,natom),vpsp(nfft) 124 125 !Local variables------------------------------- 126 character(len=1) :: geocode 127 !testing variables 128 !scalars 129 integer,parameter :: nStep=2 130 integer :: comm_fft,countParSeconde,i1,i2,i3 131 integer :: ia,ia1,ia2,igeo,ii,ind,itypat,ix,iy,iz,jj 132 integer :: kk,me_fft,n1,n2,n3,n3d,n_interpol 133 integer :: nproc_fft,tpsStart,tpsStop 134 real(dp),parameter :: min_rho_value=1.0d-12 135 real(dp) :: aa,bb,cc,dd,delta,deltaV,dr,dr2div6,invdr,r,vol_interpol,x,y,z,hgx,hgy,hgz,entmp 136 logical,parameter :: customRho=.false.,finiteDiff=.false.,testing=.false. 137 logical :: doIt 138 character(len=500) :: message 139 !arrays 140 integer :: ngfft_interpol(18) 141 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 142 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 143 real(dp) :: coord(3),coordXYZ(3),refValue(3),tsec(2) 144 real(dp),allocatable :: coordCart_interpol(:,:),coordRed_interpol(:,:) 145 real(dp),allocatable :: gridcart(:,:) 146 real(dp),allocatable :: grtn_cart_interpol(:,:),grtn_diff(:,:) 147 real(dp),allocatable :: rhog_interpol(:,:),rhog_testing(:,:),rhor_interpol(:) 148 real(dp),allocatable :: rhor_testing(:),rhor_work(:),xcart(:,:),vhartr(:),gxyz(:,:) 149 type(coulomb_operator):: kernel 150 151 ! ************************************************************************* 152 153 !Keep track of total time spent here 154 if (option==2) then 155 call timab(72,1,tsec) 156 end if 157 158 !Several constants (FFT sizes and parallelism) 159 n1 = ngfft(1) ; n2 = ngfft(2) ; n3 = ngfft(3) 160 nproc_fft = ngfft(10) ; me_fft = ngfft(11) 161 n3d = ngfft(13) !for parallel runs 162 if (nproc_fft==1) n3d=n3 !for serial runs 163 comm_fft=mpi_enreg%comm_fft 164 if(me_fft /= mpi_enreg%me_fft .or. nproc_fft /= mpi_enreg%nproc_fft) then 165 ABI_BUG("mpi_enreg%x_fft not equal to the corresponding values in ngfft") 166 end if 167 168 !Conditions for periodicity in the three directions 169 geocode='P' 170 if (icoulomb==1) geocode='F' 171 if (icoulomb==2) geocode='S' 172 173 !Get the distrib associated with this fft_grid 174 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 175 176 !Store xcart for each atom 177 ABI_MALLOC(xcart,(3, natom)) 178 call xred2xcart(natom, rprimd, xcart, xred) 179 !Store cartesian coordinates for each grid points 180 ABI_MALLOC(gridcart,(3, nfft)) 181 call mkgrid_fft(ffti3_local,fftn3_distrib,gridcart,nfft,ngfft,rprimd) 182 183 !Check whether all the PSP considered are of type GTH-HGH or PAW 184 doIt=.true. 185 do ii=1,psps%npsp 186 doIt=doIt .and.& 187 (psps%pspcod(ii)==2.or.psps%pspcod(ii)==3.or.psps%pspcod(ii)==10.or.psps%pspcod(ii)==7) 188 end do 189 190 !HGH-GTH/PAW treatment presumably starts here 191 if (doIt) then 192 193 ! Definition of the grid spacings as in the kernel routine 194 hgx = rprimd(1,1)/(n1) 195 hgy = rprimd(2,2)/(n2) 196 hgz = rprimd(3,3)/(n3) 197 198 199 !---------------------------------------------------------------------- 200 ! ----- Option 1: compute local ionic potential ----- 201 !---------------------------------------------------------------------- 202 if (option==1) then 203 204 call psolver_kernel( (/ hgx, hgy, hgz /), 2, icoulomb, me_fft, kernel, comm_fft, & 205 & (/n1,n2,n3/), nproc_fft, nscforder) 206 207 call createIonicPotential_new(fftn3_distrib,ffti3_local,& 208 & geocode,me_fft, nproc_fft, natom, & 209 & ntypat, typat, psps%gth_params%psppar, & 210 & int(psps%ziontypat), xcart,gridcart, hgx,hgy,hgz, & 211 & n1,n2,n3d,n3, kernel, vpsp, comm_fft,pawtab,psps%usepaw) 212 213 !---------------------------------------------------------------------- 214 ! ----- Option 2: compute forces induced by local ionic potential ----- 215 !---------------------------------------------------------------------- 216 else if (option == 2) then 217 218 ! Compute Hartree potential from rhor 219 ABI_MALLOC(vhartr,(nfft)) 220 call psolver_hartree(entmp, (/ hgx, hgy, hgz /), icoulomb, me_fft, comm_fft, nfft, & 221 & (/n1,n2,n3/), nproc_fft, nscforder, nspden, rhor, vhartr, usewvl) 222 223 ! Allocate temporary array for forces 224 ABI_MALLOC(gxyz,(3, natom)) 225 226 ! Calculate local part of the forces grtn (inspired from BigDFT routine) 227 call local_forces_new(fftn3_distrib,ffti3_local,geocode,me_fft, ntypat, natom, & 228 & typat, xcart, gridcart, psps%gth_params%psppar, int(psps%ziontypat), & 229 & hgx,hgy,hgz, n1,n2,n3,n3d, rhor,vhartr, gxyz, pawtab,psps%usepaw) 230 231 ! Forces should be in reduced coordinates. 232 do ia = 1, natom, 1 233 do igeo = 1, 3, 1 234 grtn(igeo, ia) = - rprimd(1, igeo) * gxyz(1, ia) & 235 & - rprimd(2, igeo) * gxyz(2, ia) & 236 & - rprimd(3, igeo) * gxyz(3, ia) 237 end do 238 end do 239 240 ! Deallocate local variables 241 ABI_FREE(vhartr) 242 ABI_FREE(gxyz) 243 end if 244 245 !---------------------------------------------------------------------- 246 ! ----- Section for the non-HGH/GTH/PAW pseudopotentials (testing) ---- 247 !---------------------------------------------------------------------- 248 else 249 250 if (testing) then 251 call system_clock(count_rate = countParSeconde) 252 call system_clock(tpsStart, count_rate = countParSeconde) 253 end if 254 255 ! dr is the r step in the sampling psps%vlspl 256 dr = psps%qgrid_vl(2) 257 invdr = 1._dp / dr 258 dr2div6 = dr * dr / 6._dp 259 260 if (option == 1) then 261 ! Set 0 in vpsp before summing 262 vpsp(:) = 0._dp 263 else if (option == 2) then 264 ! Allocate array to store cartesian gradient computed with 265 ! an interpolation of rhor 266 ABI_MALLOC(grtn_cart_interpol,(3, natom)) 267 grtn_cart_interpol(:, :) = 0._dp 268 269 n_interpol = nStep ** 3 270 ABI_MALLOC(coordRed_interpol,(3, nStep ** 3)) 271 ABI_MALLOC(coordCart_interpol,(3, nStep ** 3)) 272 273 if (testing .and. customRho) then 274 ! Use a custom rho instead of the self-consistent one. 275 ABI_MALLOC(rhor_testing,(nfft)) 276 ABI_MALLOC(rhog_testing,(2, nfft)) 277 end if 278 279 ABI_MALLOC(rhor_interpol,(nfft * n_interpol)) 280 ABI_MALLOC(rhor_work,(nfft * n_interpol)) 281 ABI_MALLOC(rhog_interpol,(2, nfft * n_interpol)) 282 283 if (testing .and. customRho) then 284 ! Testing only, changing rho with a centered gaussian 285 do ii = 1, nfft, 1 286 ! using the position of the first atom as center. 287 r = (gridcart(1, ii) - xcart(1, 1)) ** 2 + & 288 & (gridcart(2, ii) - xcart(2, 1)) ** 2 + & 289 & (gridcart(3, ii) - xcart(3, 1)) ** 2 290 rhor_testing(ii) = exp(-r/4._dp) 291 end do 292 ! Testing only, compute rhog_testing from rhor_testing 293 call fourdp(1,rhog_testing,rhor_testing,-1,mpi_enreg,nfft,1,ngfft,0) 294 end if 295 296 ! Compute the interpolation of rho, using a fourier transform 297 rhog_interpol(:, :) = 0._dp 298 ii = 0 299 do i3 = 1, n3, 1 300 if (i3 <= n3 / 2) then 301 iz = i3 302 else 303 iz = n3 * nStep - n3 + i3 304 end if 305 do i2 = 1, n2, 1 306 if (i2 <= n2 / 2) then 307 iy = i2 308 else 309 iy = n2 * nStep - n2 + i2 310 end if 311 do i1 = 1, n1, 1 312 ii = ii + 1 313 if (i1 <= n1 / 2) then 314 ix = i1 315 else 316 ix = n1 * nStep - n1 + i1 317 end if 318 jj = (iz - 1) * n2 * n1 * nStep ** 2 + (iy - 1) * n3 * nStep + ix 319 if (testing .and. customRho) then 320 rhog_interpol(:, jj) = rhog_testing(:, ii) 321 else 322 rhog_interpol(:, jj) = rhog(:, ii) 323 end if 324 end do 325 end do 326 end do 327 328 ! Compute the interpolation of rho from the Fourier transformation 329 ngfft_interpol(:) = ngfft(:) 330 ngfft_interpol(1:3) = (/ n1 * nStep, n2 * nStep, n3 * nStep /) 331 ngfft_interpol(4:6) = (/ n1 * nStep + 1, n2 * nStep + 1, n3 * nStep /) 332 call fourdp(1,rhog_interpol,rhor_work,1,mpi_enreg,nfft*n_interpol,1,ngfft_interpol,0) 333 334 ! Reorder rhor_interpol to be able to read it linearly 335 jj = 0 336 do i3 = 1, n3, 1 337 do i2 = 1, n2, 1 338 do i1 = 1, n1, 1 339 do iz = 1, nStep, 1 340 do iy = 1, nStep, 1 341 do ix = 1, nStep, 1 342 jj = jj + 1 343 kk = ((i3 - 1) * nStep + iz - 1) ! z coordinate in the interpolated grid 344 kk = kk * n1 * n2 * nStep ** 2 345 kk = kk + ((i2 - 1) * nStep + iy - 1) * n1 * nStep ! adding y coordinate 346 kk = kk + (i1 - 1) * nStep + ix ! adding x coordinate 347 rhor_interpol(jj) = rhor_work(kk) 348 end do 349 end do 350 end do 351 end do 352 end do 353 end do 354 ABI_FREE(rhor_work) 355 356 ! Compute grid access in the interpolated volume 357 ii = 0 358 do iz = 1, nStep, 1 359 z = real(iz - 1, dp) / real(nStep, dp) 360 do iy = 1, nStep, 1 361 y = real(iy - 1, dp) / real(nStep, dp) 362 do ix = 1, nStep, 1 363 x = real(ix - 1, dp) / real(nStep, dp) 364 ii = ii + 1 365 coordRed_interpol(:, ii) = (/ x, y, z /) 366 ! Assuming orthogonal box (should be change later) 367 coordCart_interpol(:, ii) = (/ x * rprimd(1, 1) / real(n1, dp), & 368 & y * rprimd(2, 2) / real(n2, dp), & 369 & z * rprimd(3, 3) / real(n3, dp) /) 370 end do 371 end do 372 end do 373 374 vol_interpol = 1._dp / real(nStep, dp) ** 3 375 ! Compute the coordinates (integer) of each atom and deduce 376 ! the max extens of the integral summation. 377 ! !! do ia = 1, natom, 1 378 ! !! coordAtom(1, ia) = int(xred(1, ia) * n1) + 1 379 ! !! coordAtom(2, ia) = int(xred(2, ia) * n2) + 1 380 ! !! coordAtom(3, ia) = int(xred(3, ia) * n3) + 1 381 ! !! end do 382 end if 383 384 if (testing .and. option == 2) then 385 call system_clock(tpsStop, count_rate = countParSeconde) 386 write(std_out,*) "Tps : ", real(tpsStop - tpsStart) / real(countParSeconde) 387 end if 388 389 ia1=1 390 do itypat = 1, ntypat, 1 391 ! ia1,ia2 sets range of loop over atoms: 392 ia2 = ia1 + nattyp(itypat) - 1 393 394 do ii = 1, nfft, 1 395 do ia = ia1, ia2, 1 396 if (option == 1) then 397 ! Compute the potential 398 ! r is the distance between grid point and atom 399 r = sqrt((gridcart(1, ii) - xcart(1, ia)) ** 2 + & 400 & (gridcart(2, ii) - xcart(2, ia)) ** 2 + & 401 & (gridcart(3, ii) - xcart(3, ia)) ** 2) 402 403 ! Coefficients needed to compute the spline. 404 jj = int(r * invdr) + 1 405 if (jj > psps%mqgrid_vl - 2) then 406 write(message, '(3a,i0,a,i0,a,a)' )& 407 & ' pseudo-potential local part sampling is not wide enough', ch10, & 408 & ' want to access position ', jj, ' whereas mqgrid_vl = ', psps%mqgrid_vl, ch10, & 409 & ' Action : no idea, contact developpers...' 410 ABI_ERROR(message) 411 end if 412 delta = r - psps%qgrid_vl(jj) 413 bb = delta * invdr 414 aa = 1._dp - bb 415 cc = aa * (aa ** 2 - 1._dp) * dr2div6 416 dd = bb * (bb ** 2 - 1._dp) * dr2div6 417 418 ! compute V(r) from the spline, jj and jj + 1 is braketting r in 419 ! the sampling 420 deltaV = aa * psps%vlspl(jj, 1, itypat) + bb * psps%vlspl(jj + 1, 1, itypat) + & 421 & cc * psps%vlspl(jj, 2, itypat) + dd * psps%vlspl(jj + 1, 2, itypat) 422 ! Add on grid point ii the contribution of atom ia 423 vpsp(ii) = vpsp(ii) + deltaV 424 else if (option == 2) then 425 ! Compute the forces, as gradient of energy (V(r).rho(r)) 426 427 ! Testing only - reference points 428 if (.false.) then 429 ! r is the distance between grid point and atom 430 r = sqrt((gridcart(1, ii) - xcart(1, ia)) ** 2 + & 431 & (gridcart(2, ii) - xcart(2, ia)) ** 2 + & 432 & (gridcart(3, ii) - xcart(3, ia)) ** 2) 433 434 ! Coefficients needed to compute the spline. 435 jj = int(r * invdr) + 1 436 delta = r - psps%qgrid_vl(jj) 437 bb = delta * invdr 438 aa = 1._dp - bb 439 cc = aa * (aa ** 2 - 1._dp) * dr2div6 440 dd = bb * (bb ** 2 - 1._dp) * dr2div6 441 442 ! When mesh position is on a node, forces are null. 443 if (r /= 0._dp) then 444 ! This value deltaV is the first derivative of V(r) taken at r. 445 deltaV = aa * psps%dvlspl(jj, 1, itypat) + bb * psps%dvlspl(jj + 1, 1, itypat) + & 446 & cc * psps%dvlspl(jj, 2, itypat) + dd * psps%dvlspl(jj + 1, 2, itypat) 447 ! We multiply by rho(r) to have an energy. 448 deltaV = deltaV * rhor(ii, 1) / r 449 refValue(:) = - deltaV * (gridcart(:, ii) - xcart(:, ia)) 450 grtn_cart_interpol(:, ia) = grtn_cart_interpol(:, ia) + refValue(:) 451 end if 452 end if 453 454 ! Compute the interpolation for the point ii 455 ind = (ii - 1) * n_interpol 456 do kk = 1, n_interpol, 1 457 ind = ind + 1 458 459 if (rhor_interpol(ind) > min_rho_value) then 460 ! Assume orthogonal box... 461 coordXYZ(1) = gridcart(1, ii) - xcart(1, ia) + coordCart_interpol(1, kk) 462 coordXYZ(2) = gridcart(2, ii) - xcart(2, ia) + coordCart_interpol(2, kk) 463 coordXYZ(3) = gridcart(3, ii) - xcart(3, ia) + coordCart_interpol(3, kk) 464 r = coordXYZ(1) ** 2 + coordXYZ(2) ** 2 + coordXYZ(3) ** 2 465 466 if (r /= 0._dp) then 467 r = sqrt(r) 468 ! Coefficients needed to compute the spline. 469 jj = int(r * invdr) + 1 470 delta = r - psps%qgrid_vl(jj) 471 bb = delta * invdr 472 aa = 1._dp - bb 473 cc = aa * (aa ** 2 - 1._dp) * dr2div6 474 dd = bb * (bb ** 2 - 1._dp) * dr2div6 475 deltaV = aa * psps%dvlspl(jj, 1, itypat) + & 476 & bb * psps%dvlspl(jj + 1, 1, itypat) + & 477 & cc * psps%dvlspl(jj, 2, itypat) + & 478 & dd * psps%dvlspl(jj + 1, 2, itypat) 479 deltaV = deltaV * rhor_interpol(ind) / r 480 grtn_cart_interpol(1, ia) = grtn_cart_interpol(1, ia) - deltaV * coordXYZ(1) 481 grtn_cart_interpol(2, ia) = grtn_cart_interpol(2, ia) - deltaV * coordXYZ(2) 482 grtn_cart_interpol(3, ia) = grtn_cart_interpol(3, ia) - deltaV * coordXYZ(3) 483 ! do igeo = 1, 3, 1 484 ! grtn_cart_interpol(igeo, ia) = grtn_cart_interpol(igeo, ia) - deltaV * coordXYZ(igeo) 485 ! end do 486 end if 487 end if 488 end do 489 490 ! ============= 491 ! Testing only 492 ! ============= 493 ! use of finite differences 494 if (finiteDiff) then 495 do igeo = 1, 3, 1 496 coord(:) = 0._dp 497 coord(igeo) = dr / 2.0_dp 498 r = sqrt((gridcart(1, ii) - xcart(1, ia) + coord(1)) ** 2 + & 499 & (gridcart(2, ii) - xcart(2, ia) + coord(2)) ** 2 + & 500 & (gridcart(3, ii) - xcart(3, ia) + coord(3)) ** 2) 501 502 ! Coefficients needed to compute the spline. 503 jj = int(r * invdr) + 1 504 delta = r - psps%qgrid_vl(jj) 505 bb = delta * invdr 506 aa = 1._dp - bb 507 cc = aa * (aa ** 2 - 1._dp) * dr2div6 508 dd = bb * (bb ** 2 - 1._dp) * dr2div6 509 510 deltaV = aa * psps%vlspl(jj, 1, itypat) + bb * psps%vlspl(jj + 1, 1, itypat) + & 511 & cc * psps%vlspl(jj, 2, itypat) + dd * psps%vlspl(jj + 1, 2, itypat) 512 513 514 coord(:) = 0._dp 515 coord(igeo) = -dr / 2.0_dp 516 r = sqrt((gridcart(1, ii) - xcart(1, ia) + coord(1)) ** 2 + & 517 & (gridcart(2, ii) - xcart(2, ia) + coord(2)) ** 2 + & 518 & (gridcart(3, ii) - xcart(3, ia) + coord(3)) ** 2) 519 520 ! Coefficients needed to compute the spline. 521 jj = int(r * invdr) + 1 522 delta = r - psps%qgrid_vl(jj) 523 bb = delta * invdr 524 aa = 1._dp - bb 525 cc = aa * (aa ** 2 - 1._dp) * dr2div6 526 dd = bb * (bb ** 2 - 1._dp) * dr2div6 527 528 deltaV = deltaV - (aa * psps%vlspl(jj, 1, itypat) + & 529 & bb * psps%vlspl(jj + 1, 1, itypat) + & 530 & cc * psps%vlspl(jj, 2, itypat) + & 531 & dd * psps%vlspl(jj + 1, 2, itypat)) 532 grtn_diff(igeo, ia) = grtn_diff(igeo, ia) - deltaV * rhor(ii, 1) / dr 533 end do 534 end if 535 ! ============= 536 ! Testing only 537 ! ============= 538 539 end if 540 end do 541 ! End loop over atoms of type itypat 542 end do 543 ! End loop over real space grid points 544 545 ia1 = ia2 + 1 546 end do 547 ! End loop over type of atoms 548 549 if(option==2)then 550 ! multiply the forces by the volume of a single box mesh. 551 grtn_cart_interpol(:, :) = grtn_cart_interpol(:, :) * & 552 & ucvol / real(n1 * n2 * n3, dp) * vol_interpol 553 ! Transform cartesian forces to reduce coordinates 554 do ia = 1, natom, 1 555 do igeo = 1, 3, 1 556 grtn(igeo, ia) = rprimd(1, igeo) * grtn_cart_interpol(1, ia) + & 557 & rprimd(2, igeo) * grtn_cart_interpol(2, ia) + & 558 & rprimd(3, igeo) * grtn_cart_interpol(3, ia) 559 end do 560 end do 561 ABI_FREE(rhor_interpol) 562 ABI_FREE(rhog_interpol) 563 ABI_FREE(coordRed_interpol) 564 ABI_FREE(coordCart_interpol) 565 if (testing .and. customRho) then 566 ABI_FREE(rhor_testing) 567 ABI_FREE(rhog_testing) 568 end if 569 570 if (testing) then 571 call system_clock(tpsStop, count_rate = countParSeconde) 572 write(std_out,*) "Tps : ", real(tpsStop - tpsStart) / real(countParSeconde) 573 write(std_out,*) grtn_cart_interpol 574 ABI_ERROR("Testing section!") 575 end if 576 577 end if 578 579 !------------------------- 580 end if ! GTH/HGH/PAW psps 581 582 !Release temporary memory 583 ABI_FREE(xcart) 584 ABI_FREE(gridcart) 585 586 !Close timing counters 587 if (option==2)then 588 call timab(72,2,tsec) 589 end if 590 591 end subroutine mklocl_realspace
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
INPUTS
efield (3)=external electric field mpi_enreg=information 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.
SOURCE
1337 subroutine mklocl_wavelets(efield, grtn, mpi_enreg, natom, nfft, & 1338 & nspden, option, rprimd, vpsp, wvl_den, wvl_descr, xcart) 1339 1340 use defs_wvltypes 1341 #if defined HAVE_BIGDFT 1342 use BigDFT_API, only : ELECTRONIC_DENSITY,createIonicPotential,local_forces 1343 use poisson_solver, only : H_potential 1344 #endif 1345 1346 !Arguments ------------------------------------ 1347 !scalars 1348 integer,intent(in) :: option, natom, nfft, nspden 1349 type(MPI_type),intent(in) :: mpi_enreg 1350 type(wvl_denspot_type), intent(inout) :: wvl_den 1351 type(wvl_internal_type), intent(in) :: wvl_descr 1352 !arrays 1353 real(dp),intent(in) :: rprimd(3,3),efield(3) 1354 real(dp),intent(inout) :: grtn(3,natom) 1355 real(dp), intent(inout) :: vpsp(nfft) 1356 real(dp),intent(inout) :: xcart(3,natom) 1357 1358 !Local variables------------------------------- 1359 #if defined HAVE_BIGDFT 1360 !scalars 1361 integer :: i,i1,i2,i3,ia,ierr,igeo,me,nproc,shift,comm 1362 real(dp) :: energ 1363 character(len=500) :: message 1364 !arrays 1365 real(dp) :: epot(3) 1366 real(dp) :: elecfield(3)=(/zero,zero,zero/) ! Not used here 1367 real(dp),allocatable :: gxyz(:,:),vhartr(:),rhov(:,:) 1368 #endif 1369 1370 ! ********************************************************************* 1371 1372 #if defined HAVE_BIGDFT 1373 1374 elecfield=zero !not used here 1375 1376 !Manage parallelism 1377 comm=mpi_enreg%comm_wvl 1378 nproc=xmpi_comm_size(comm) 1379 me=xmpi_comm_rank(comm) 1380 1381 !---------------------------------------------------------------------- 1382 ! ----- Option 1: compute local ionic potential ----- 1383 !---------------------------------------------------------------------- 1384 if (option == 1) then 1385 1386 ! We get the kernel for the Poisson solver (used to go from the ionic 1387 ! charge to the potential).If the kernel is uncomputed, it does it now. 1388 ! call psolver_kernel(wvl_den%denspot%dpbox%hgrids, 2, icoulomb, me, kernel, & 1389 !& comm, wvl_den%denspot%dpbox%ndims, nproc, nscforder) 1390 ! if (.not.associated(kernel%co%kernel)) then 1391 ! call psolver_kernel(wvl_den%denspot%dpbox%hgrids, 1, icoulomb, me, kernel, & 1392 !& comm, wvl_den%denspot%dpbox%ndims, nproc, nscforder) 1393 ! end if 1394 1395 message=ch10//' mklocl_wavelets: Create local potential from ions.' 1396 call wrtout(std_out,message,'COLL') 1397 1398 shift = 1 + wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) & 1399 & * wvl_den%denspot%dpbox%i3xcsh 1400 1401 ! Call the BigDFT routine 1402 call createIonicPotential(wvl_descr%atoms%astruct%geocode, me, nproc, (me == 0), wvl_descr%atoms, & 1403 & xcart, wvl_den%denspot%dpbox%hgrids(1), wvl_den%denspot%dpbox%hgrids(2), & 1404 & wvl_den%denspot%dpbox%hgrids(3), & 1405 & elecfield, wvl_descr%Glr%d%n1, wvl_descr%Glr%d%n2, wvl_descr%Glr%d%n3, & 1406 & wvl_den%denspot%dpbox%n3pi, wvl_den%denspot%dpbox%i3s + wvl_den%denspot%dpbox%i3xcsh, & 1407 & wvl_den%denspot%dpbox%ndims(1), wvl_den%denspot%dpbox%ndims(2), & 1408 & wvl_den%denspot%dpbox%ndims(3), wvl_den%denspot%pkernel, vpsp(shift), 0.d0,wvl_descr%rholoc) 1409 1410 ! Copy vpsp into bigdft object: 1411 call wvl_rhov_abi2big(1,vpsp,wvl_den%denspot%v_ext,shift=(shift-1)) 1412 1413 ! Eventually add the electric field 1414 if (maxval(efield) > tol12) then 1415 message=ch10//'mklocl_wavelets: Add the electric field.' 1416 call wrtout(std_out,message,'COLL') 1417 ! We add here the electric field since in BigDFT, the field must be on x... 1418 epot(:) = real(0.5, dp) * efield(:) * wvl_den%denspot%dpbox%hgrids(:) 1419 do i3 = 1, wvl_den%denspot%dpbox%n3pi, 1 1420 ia = (i3 - 1) * wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) 1421 do i2 = -14, 2 * wvl_descr%Glr%d%n2 + 16, 1 1422 i = ia + (i2 + 14) * wvl_den%denspot%dpbox%ndims(1) 1423 do i1 = -14, 2 * wvl_descr%Glr%d%n1 + 16, 1 1424 i = i + 1 1425 vpsp(shift + i) = vpsp(shift + i) + & 1426 & epot(1) * real(i1 - wvl_descr%Glr%d%n1, dp) + & 1427 & epot(2) * real(i2 - wvl_descr%Glr%d%n2, dp) + & 1428 & epot(3) * real(i3 - wvl_descr%Glr%d%n3, dp) 1429 end do 1430 end do 1431 end do 1432 end if 1433 1434 !---------------------------------------------------------------------- 1435 ! ----- Option 2: compute forces induced by local ionic potential ----- 1436 !---------------------------------------------------------------------- 1437 else if (option == 2) then 1438 1439 message=ch10//' mklocl_wavelets: compute local forces.' 1440 call wrtout(std_out,message,'COLL') 1441 1442 if (wvl_den%denspot%rhov_is/=ELECTRONIC_DENSITY) then 1443 message='denspot bigdft datstructure should contain rhor!' 1444 ABI_BUG(message) 1445 end if 1446 1447 ! Extract density rhor from bigDFT datastructure 1448 ABI_MALLOC(rhov,(nfft, nspden)) 1449 ABI_MALLOC(vhartr,(nfft)) 1450 rhov=zero ; vhartr=zero 1451 shift = wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) & 1452 & * wvl_den%denspot%dpbox%i3xcsh 1453 do i = 1, min(nfft,size(wvl_den%denspot%rhov)-shift) 1454 rhov(i, 1) = wvl_den%denspot%rhov(i + shift) 1455 vhartr(i) = wvl_den%denspot%rhov(i + shift) 1456 end do 1457 if (nspden == 2) then 1458 shift = shift + wvl_den%denspot%dpbox%ndims(1) * wvl_den%denspot%dpbox%ndims(2) & 1459 & * wvl_den%denspot%dpbox%n3d 1460 do i = 1, min(nfft,size(wvl_den%denspot%rhov)-shift) 1461 rhov(i, 2) = wvl_den%denspot%rhov(i + shift) 1462 vhartr(i) = vhartr(i) + wvl_den%denspot%rhov(i + shift) 1463 end do 1464 end if 1465 1466 ! Compute Hartree potential from rhor 1467 call H_potential('D',wvl_den%denspot%pkernel,vhartr,vhartr,energ,zero,.false.) 1468 1469 ! Allocate temporary array for forces 1470 ABI_MALLOC(gxyz,(3, natom)) 1471 1472 ! Calculate local part of the forces grtn (modified BigDFT routine) 1473 call local_forces_wvl(me,natom,xcart,& 1474 & wvl_den%denspot%dpbox%hgrids(1),& 1475 & wvl_den%denspot%dpbox%hgrids(2),& 1476 & wvl_den%denspot%dpbox%hgrids(3),& 1477 & wvl_descr%Glr%d%n1,wvl_descr%Glr%d%n2,wvl_descr%Glr%d%n3,& 1478 & wvl_den%denspot%dpbox%n3p,& 1479 & wvl_den%denspot%dpbox%i3s+wvl_den%denspot%dpbox%i3xcsh,& 1480 & wvl_den%denspot%dpbox%ndims(1),wvl_den%denspot%dpbox%ndims(2),& 1481 & rhov,vhartr,gxyz,wvl_descr) 1482 ! call local_forces(me, wvl_descr%atoms, xcart, & 1483 ! & wvl_den%denspot%dpbox%hgrids(1), wvl_den%denspot%dpbox%hgrids(2), & 1484 ! & wvl_den%denspot%dpbox%hgrids(3), & 1485 ! & wvl_descr%Glr%d%n1, wvl_descr%Glr%d%n2, wvl_descr%Glr%d%n3, & 1486 ! & wvl_den%denspot%dpbox%n3p, wvl_den%denspot%dpbox%i3s + wvl_den%denspot%dpbox%i3xcsh, & 1487 ! & wvl_den%denspot%dpbox%ndims(1), wvl_den%denspot%dpbox%ndims(2), & 1488 ! & rhov, vhartr,gxyz,locstrten,charge) 1489 1490 ! Pending: floc,locstrten and charge are not used here. 1491 ! Pending: put mpi_enreg%nscatterarr... in object denspot, initialize object, etc. 1492 1493 if (nproc > 1) then 1494 call xmpi_sum(gxyz, comm, ierr) 1495 end if 1496 1497 ! Forces should be in reduced coordinates. 1498 do ia = 1, natom, 1 1499 do igeo = 1, 3, 1 1500 grtn(igeo, ia) = - rprimd(1, igeo) * gxyz(1, ia) & 1501 & - rprimd(2, igeo) * gxyz(2, ia) & 1502 & - rprimd(3, igeo) * gxyz(3, ia) 1503 end do 1504 end do 1505 1506 ! Deallocate local variables 1507 ABI_FREE(vhartr) 1508 ABI_FREE(rhov) 1509 ABI_FREE(gxyz) 1510 1511 !---------------------------------------------------------------------- 1512 1513 else ! option switch 1514 message = 'Internal error, option should be 1 or 2!' 1515 ABI_ERROR(message) 1516 end if 1517 1518 #else 1519 BIGDFT_NOTENABLED_ERROR() 1520 if (.false.) write(std_out,*) option,natom,nfft,nspden,mpi_enreg%me,& 1521 & wvl_den%symObj,wvl_descr%h(1),rprimd(1,1),efield(1),grtn(1,1),vpsp(1),xcart(1,1) 1522 #endif 1523 1524 end subroutine mklocl_wavelets
mklocl_realspace/calcdVloc_mklocl [ Functions ]
[ Top ] [ mklocl_realspace ] [ Functions ]
NAME
calcdVloc_mklocl
FUNCTION
Compute 1st-derivative of long-range HGH local ionic potential (derf)
INPUTS
OUTPUT
SOURCE
1175 subroutine calcdVloc_mklocl(yy,xx,rloc,Z) 1176 1177 !Arguments ------------------------------------ 1178 !scalars 1179 real(dp),intent(in) :: xx,rloc,Z 1180 real(dp),intent(out) :: yy 1181 1182 !Local variables------------------------------- 1183 !scalars 1184 real(dp):: arg,tt 1185 1186 ! ************************************************************************* 1187 1188 arg=xx/(sqrt(2._dp)*rloc) 1189 call derf_ab(tt,arg) 1190 yy=(Z/(xx**2))* ( tt - 2._dp/sqrt(pi)*arg*exp(-arg**2) ) 1191 1192 end subroutine calcdVloc_mklocl
mklocl_realspace/calcVloc_mklocl [ Functions ]
[ Top ] [ mklocl_realspace ] [ Functions ]
NAME
calcVloc_mklocl
FUNCTION
INPUTS
OUTPUT
SOURCE
877 subroutine calcVloc_mklocl(yy,xx,rloc,Z) 878 879 !Arguments ------------------------------------ 880 !scalars 881 real(dp),intent(in) :: xx,rloc,Z 882 real(dp),intent(out) :: yy 883 884 !Local variables------------------------------- 885 !scalars 886 real(dp):: arg,tt 887 888 ! ************************************************************************* 889 890 arg=xx/(sqrt(2.0)*rloc) 891 call derf_ab(tt,arg) 892 yy=-Z/xx*tt 893 894 end subroutine calcVloc_mklocl
mklocl_realspace/createIonicPotential_new [ Functions ]
[ Top ] [ mklocl_realspace ] [ Functions ]
NAME
createIonicPotential_new
FUNCTION
INPUTS
OUTPUT
SOURCE
608 subroutine createIonicPotential_new(fftn3_distrib,ffti3_local,geocode,iproc,& 609 & nproc,nat,ntypes,iatype,psppar,nelpsp,rxyz,gridcart,& 610 & hxh,hyh,hzh,n1i,n2i,n3d,n3i,kernel,pot_ion,spaceworld,pawtab,usepaw) 611 612 use defs_wvltypes, only : coulomb_operator 613 614 !Arguments ------------------------------- 615 !scalars 616 integer, intent(in) :: iproc,nproc,ntypes,nat,n1i,n2i,n3i,n3d,spaceworld,usepaw 617 real(dp), intent(in) :: hxh,hyh,hzh 618 character(len=1), intent(in) :: geocode 619 type(coulomb_operator), intent(in) :: kernel 620 !arrays 621 integer, dimension(nat), intent(in) :: iatype 622 integer, dimension(ntypes), intent(in) :: nelpsp 623 integer, dimension(*), intent(in) ::fftn3_distrib,ffti3_local 624 real(dp), dimension(3,n1i*n2i*n3d), intent(in) :: gridcart 625 real(dp), dimension(0:4,0:6,ntypes), intent(in) :: psppar 626 real(dp), dimension(3,nat), intent(in) :: rxyz 627 real(dp), dimension(*), intent(inout) :: pot_ion 628 type(pawtab_type),intent(in) :: pawtab(ntypes*usepaw) 629 630 !Local variables ------------------------- 631 #if defined HAVE_BIGDFT 632 !scalars 633 integer :: iat,i1,i2,i3,j1,j2,j3,isx,isy,isz,iex,iey,iez,ierr,ityp 634 integer :: ind,nloc,iloc,i3loc,msz 635 logical :: gox,goy,goz,perx,pery,perz 636 real(dp) :: arg,charge,cutoff,ehart,eexcu,rholeaked,rholeaked_tot,rloc 637 real(dp) :: rx,ry,rz,rzero,r2,tt,tt_tot,vexcu,vhgh,x,xp,y,z 638 !arrays 639 real(dp) :: rr(1),vpaw(1) 640 real(dp),pointer :: rad(:),vloc(:),d2vloc(:) 641 #endif 642 643 ! ********************************************************************* 644 645 #if defined HAVE_BIGDFT 646 647 if(nproc<0)then 648 ABI_ERROR('nproc should not be negative') 649 end if 650 651 !Ionic charge (must be calculated for the PS active processes) 652 rholeaked=0._dp 653 !Ionic energy (can be calculated for all the processors) 654 655 !here we should insert the calculation of the ewald energy for the periodic BC case 656 !!! eion=0._dp 657 !!! do iat=1,nat 658 !!! ityp=iatype(iat) 659 !!! rx=rxyz(1,iat) 660 !!! ry=rxyz(2,iat) 661 !!! rz=rxyz(3,iat) 662 !!! ! ion-ion interaction 663 !!! do jat=1,iat-1 664 !!! dist=sqrt( (rx-rxyz(1,jat))**2+(ry-rxyz(2,jat))**2+(rz-rxyz(3,jat))**2 ) 665 !!! jtyp=iatype(jat) 666 !!! eion=eion+real(nelpsp(jtyp)*nelpsp(ityp),kind=dp)/dist 667 !!! enddo 668 !!! end do 669 !!! if (iproc.eq.0) write(std_out,'(1x,a,1pe22.14)') 'ion-ion interaction energy',eion 670 671 !Creates charge density arising from the ionic PSP cores 672 !the n3pi dimension indicates the number of planes trated by each processor in the FFT parallelisation 673 !for a plane wave treatment this value depends on whether the direct space is divided in planes or not 674 !I don't know this variable, which in the future must be inserted at the place of n3pi (LG) 675 !if n3pi=0 this means that the processors doesn't calculate anything 676 !if (n3pi >0 ) then 677 678 !conditions for periodicity in the three directions 679 perx=(geocode /= 'F') 680 pery=(geocode == 'P') 681 perz=(geocode /= 'F') 682 683 !this initialise the array to zero, it will work only if bigdft library is enabled 684 pot_ion(1:n1i*n2i*n3d)=zero 685 686 do iat=1,nat 687 ityp=iatype(iat) 688 rx=rxyz(1,iat) 689 ry=rxyz(2,iat) 690 rz=rxyz(3,iat) 691 692 rloc=psppar(0,0,ityp) 693 charge=real(nelpsp(ityp),kind=dp)/(2._dp*pi*sqrt(2._dp*pi)*rloc**3) 694 cutoff=10._dp*rloc 695 696 isx=floor((rx-cutoff)/hxh) 697 isy=floor((ry-cutoff)/hyh) 698 isz=floor((rz-cutoff)/hzh) 699 700 iex=ceiling((rx+cutoff)/hxh) 701 iey=ceiling((ry+cutoff)/hyh) 702 iez=ceiling((rz+cutoff)/hzh) 703 704 ! Calculate Ionic Density 705 ! using HGH parameters. 706 ! Eq. 1.104, T. Deutsch and L. Genovese, JDN. 12, 2011 707 do i3=isz,iez 708 z=real(i3,kind=dp)*hzh-rz 709 call ind_positions_mklocl(perz,i3,n3i,j3,goz) 710 if(fftn3_distrib(j3)==iproc) then 711 i3loc=ffti3_local(j3) 712 do i2=isy,iey 713 y=real(i2,kind=dp)*hyh-ry 714 call ind_positions_mklocl(pery,i2,n2i,j2,goy) 715 do i1=isx,iex 716 x=real(i1,kind=dp)*hxh-rx 717 call ind_positions_mklocl(perx,i1,n1i,j1,gox) 718 r2=x**2+y**2+z**2 719 if (goz .and. goy .and. gox ) then 720 ind=j1+(j2-1)*n1i+(i3loc-1)*n1i*n2i 721 r2=(gridcart(1,ind)-rx)**2+(gridcart(2,ind)-ry)**2+(gridcart(3,ind)-rz)**2 722 end if 723 arg=r2/rloc**2 724 xp=exp(-.5d0*arg) 725 if (goz .and. goy .and. gox ) then 726 pot_ion(ind)=pot_ion(ind)-xp*charge 727 else 728 rholeaked=rholeaked+xp*charge 729 end if 730 end do 731 end do 732 end if 733 end do 734 735 end do 736 737 !Check 738 tt=0._dp 739 do j3= 1,n3d 740 do i2= 1,n2i 741 do i1= 1,n1i 742 ind=i1+(i2-1)*n1i+(j3-1)*n1i*n2i 743 tt=tt+pot_ion(ind) 744 end do 745 end do 746 end do 747 748 tt=tt*hxh*hyh*hzh 749 rholeaked=rholeaked*hxh*hyh*hzh 750 751 call xmpi_sum(tt,tt_tot,spaceworld,ierr) 752 call xmpi_sum(rholeaked,rholeaked_tot,spaceworld,ierr) 753 754 if (iproc.eq.0) then 755 write(std_out,'(1x,a,f26.12,2x,1pe10.3)') & 756 & 'total ionic charge, leaked charge ',tt_tot,rholeaked_tot 757 end if 758 759 !Here the value of the datacode must be kept fixed 760 !there can be some problems when running this stuff in parallel, 761 ! if the ionic potential distribution does not agree with the 762 ! plane distribution which is supposed to hold for the Poisson Solver 763 call psolver(geocode,'D',iproc,nproc,n1i,n2i,n3i,0,hxh,hyh,hzh,& 764 & pot_ion,kernel%kernel,pot_ion,ehart,eexcu,vexcu,0._dp,.false.,1) 765 766 !Add the remaining short-range local terms 767 do iat=1,nat 768 ityp=iatype(iat) 769 770 rx=rxyz(1,iat) 771 ry=rxyz(2,iat) 772 rz=rxyz(3,iat) 773 774 ! determine number of local terms 775 rloc=psppar(0,0,ityp) 776 cutoff=10._dp*rloc 777 charge=real(nelpsp(ityp),kind=dp) 778 779 ! determine number of local terms (HGH pot) 780 nloc=0 781 do iloc=1,4 782 if (psppar(0,iloc,ityp).ne.0._dp) nloc=iloc 783 end do 784 785 ! PAW specifics 786 if (usepaw==1) then 787 msz=pawtab(ityp)%wvl%rholoc%msz 788 rad => pawtab(ityp)%wvl%rholoc%rad(1:msz) 789 vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,3) 790 d2vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,4) 791 rzero=rad(1);if (rzero<=1.d-10) rzero=rad(2) 792 end if 793 794 isx=floor((rx-cutoff)/hxh) 795 isy=floor((ry-cutoff)/hyh) 796 isz=floor((rz-cutoff)/hzh) 797 798 iex=ceiling((rx+cutoff)/hxh) 799 iey=ceiling((ry+cutoff)/hyh) 800 iez=ceiling((rz+cutoff)/hzh) 801 802 do i3=isz,iez 803 z=real(i3,kind=dp)*hzh-rz 804 call ind_positions_mklocl(perz,i3,n3i,j3,goz) 805 if(fftn3_distrib(j3) == iproc .and. goz) then !MPI 806 i3loc=ffti3_local(j3) 807 if (goz) then 808 do i2=isy,iey 809 y=real(i2,kind=dp)*hyh-ry 810 call ind_positions_mklocl(pery,i2,n2i,j2,goy) 811 if (goy) then 812 do i1=isx,iex 813 x=real(i1,kind=dp)*hxh-rx 814 call ind_positions_mklocl(perx,i1,n1i,j1,gox) 815 if (gox) then 816 ind=j1+(j2-1)*n1i+(i3loc-1)*n1i*n2i 817 r2=(gridcart(1,ind)-rx)**2+(gridcart(2,ind)-ry)**2+(gridcart(3,ind)-rz)**2 818 ! r2=x**2+y**2+z**2 819 820 ! HGH: V_S=gaussian potential of Eq. (9) in JCP 129, 014109(2008) 821 if (usepaw==0) then 822 if (nloc /= 0) then 823 arg=r2/rloc**2 824 xp=exp(-.5d0*arg) 825 tt=psppar(0,nloc,ityp) 826 do iloc=nloc-1,1,-1 827 tt=arg*tt+psppar(0,iloc,ityp) 828 end do 829 pot_ion(ind)=pot_ion(ind)+xp*tt 830 end if 831 832 ! PAW: V_PAW-V_L^HGH 833 else 834 rr(1)=sqrt(r2) 835 if (rr(1)>=rzero) then 836 call paw_splint(msz,rad,vloc,d2vloc,1,rr,vpaw) 837 call calcVloc_mklocl(vhgh,rr(1),rloc,charge) 838 pot_ion(ind)=pot_ion(ind)+vpaw(1)-vhgh 839 else 840 pot_ion(ind)=pot_ion(ind)+vloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc) 841 end if 842 end if 843 844 end if 845 end do 846 end if 847 end do 848 end if 849 end if 850 end do 851 852 end do !iat 853 #else 854 BIGDFT_NOTENABLED_ERROR() 855 if (.false.) write(std_out,*) geocode,iproc,nproc,ntypes,nat,n1i,n2i,n3i,n3d,spaceworld,usepaw,& 856 & hxh,hyh,hzh,iatype(1),nelpsp(1),fftn3_distrib(1),ffti3_local(1),gridcart(1,1),psppar(1,1,1),& 857 & rxyz(1,1),pot_ion(1),pawtab(1)%mesh_size,kernel%co 858 #endif 859 860 CONTAINS
mklocl_realspace/ind_positions_mklocl [ Functions ]
[ Top ] [ mklocl_realspace ] [ Functions ]
NAME
ind_positions_mklocl
FUNCTION
determine the index in which the potential must be inserted, following the BC determine also whether the index is inside or outside the box for free BC
INPUTS
OUTPUT
SOURCE
1279 subroutine ind_positions_mklocl(periodic,i,n,j,go) 1280 1281 1282 !Arguments ------------------------------- 1283 logical, intent(in) :: periodic 1284 integer, intent(in) :: i,n 1285 logical, intent(out) :: go 1286 integer, intent(out) :: j 1287 1288 ! ********************************************************************* 1289 1290 if (periodic) then 1291 go=.true. 1292 j=modulo(i-1,n)+1 1293 else 1294 j=i 1295 go=(i >= 1 .and. i <= n) 1296 end if 1297 1298 end subroutine ind_positions_mklocl
mklocl_realspace/local_forces_new [ Functions ]
[ Top ] [ mklocl_realspace ] [ Functions ]
NAME
local_forces_new
FUNCTION
INPUTS
OUTPUT
SOURCE
979 subroutine local_forces_new(fftn3_distrib,ffti3_local,& 980 geocode,iproc,ntypes,nat,iatype,rxyz,gridcart,psppar,nelpsp,hxh,hyh,hzh,& 981 n1,n2,n3,n3d,rho,pot,floc,pawtab,usepaw) 982 983 984 !Arguments ------------------------------- 985 !scalars 986 integer, intent(in) :: iproc,ntypes,nat,n1,n2,n3,n3d,usepaw 987 character(len=1), intent(in) :: geocode 988 real(dp), intent(in) :: hxh,hyh,hzh 989 !arrays 990 integer, dimension(*), intent(in) ::fftn3_distrib,ffti3_local 991 integer, dimension(nat), intent(in) :: iatype 992 integer, dimension(ntypes), intent(in) :: nelpsp 993 real(dp), dimension(3,n1*n2*n3d), intent(in) :: gridcart 994 real(dp), dimension(0:4,0:6,ntypes), intent(in) :: psppar 995 real(dp), dimension(3,nat), intent(in) :: rxyz 996 real(dp), dimension(*), intent(in) :: rho,pot 997 real(dp), dimension(3,nat), intent(out) :: floc 998 type(pawtab_type),intent(in) :: pawtab(ntypes*usepaw) 999 1000 !Local variables ------------------------- 1001 !scalars 1002 integer :: isx,isy,isz,iex,iey,iez,i1,i2,i3,j1,j2,j3,ind,iat,ityp,iloc,i3loc,msz,nloc 1003 logical :: perx,pery,perz,gox,goy,goz 1004 real(dp) :: arg,charge,cutoff,dvhgh,fxerf,fyerf,fzerf,fxgau,fygau,fzgau,forceleaked 1005 real(dp) :: forceloc,prefactor,rloc,rloc2,rhoel,rx,ry,rz,rzero,r2,tt,x,xp,y,z,Vel 1006 real(dp), dimension(4) :: cprime 1007 !arrays 1008 real(dp) :: dvpawdr(1),rr(1) 1009 real(dp),pointer :: rad(:),vloc(:),d2vloc(:) 1010 1011 ! ********************************************************************* 1012 1013 if (iproc == 0) write(std_out,'(1x,a)',advance='no')'Calculate local forces...' 1014 1015 !Conditions for periodicity in the three directions 1016 perx=(geocode /= 'F') 1017 pery=(geocode == 'P') 1018 perz=(geocode /= 'F') 1019 1020 forceleaked=zero 1021 1022 do iat=1,nat 1023 ityp=iatype(iat) 1024 ! Coordinates of the center 1025 rx=rxyz(1,iat) 1026 ry=rxyz(2,iat) 1027 rz=rxyz(3,iat) 1028 1029 ! Initialization of the forces 1030 ! ion-electron term, error function part 1031 fxerf=zero 1032 fyerf=zero 1033 fzerf=zero 1034 ! ion-electron term, gaussian part 1035 fxgau=zero 1036 fygau=zero 1037 fzgau=zero 1038 1039 ! Building array of coefficients of the derivative of the gaussian part 1040 cprime(1)=2._dp*psppar(0,2,ityp)-psppar(0,1,ityp) 1041 cprime(2)=4._dp*psppar(0,3,ityp)-psppar(0,2,ityp) 1042 cprime(3)=6._dp*psppar(0,4,ityp)-psppar(0,3,ityp) 1043 cprime(4)=-psppar(0,4,ityp) 1044 1045 ! Determine number of local terms (HGH pot) 1046 nloc=0 1047 do iloc=1,4 1048 if (psppar(0,iloc,ityp).ne.zero) nloc=iloc 1049 end do 1050 1051 ! Some constants depending on the atom type 1052 rloc=psppar(0,0,ityp) ; rloc2=rloc**2 1053 charge=real(nelpsp(ityp),kind=dp) 1054 prefactor=charge/(2._dp*pi*sqrt(2._dp*pi)*rloc**5) 1055 1056 ! PAW specifics 1057 if (usepaw==1) then 1058 msz=pawtab(ityp)%wvl%rholoc%msz 1059 rad => pawtab(ityp)%wvl%rholoc%rad(1:msz) 1060 vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,3) 1061 d2vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,4) 1062 rzero=rad(1);if (rad(1)<=1.d-10) rzero=rad(2) 1063 end if 1064 1065 ! Maximum extension of the gaussian 1066 cutoff=10._dp*rloc 1067 isx=floor((rx-cutoff)/hxh) 1068 isy=floor((ry-cutoff)/hyh) 1069 isz=floor((rz-cutoff)/hzh) 1070 iex=ceiling((rx+cutoff)/hxh) 1071 iey=ceiling((ry+cutoff)/hyh) 1072 iez=ceiling((rz+cutoff)/hzh) 1073 1074 ! Calculate the forces near the atom due to the gaussian 1075 ! and error function parts of the potential 1076 do i3=isz,iez 1077 z=real(i3,kind=dp)*hzh-rz 1078 call ind_positions_mklocl(perz,i3,n3,j3,goz) 1079 if(fftn3_distrib(j3)==iproc) then 1080 i3loc=ffti3_local(j3) 1081 do i2=isy,iey 1082 y=real(i2,kind=dp)*hyh-ry 1083 call ind_positions_mklocl(pery,i2,n2,j2,goy) 1084 do i1=isx,iex 1085 x=real(i1,kind=dp)*hxh-rx 1086 call ind_positions_mklocl(perx,i1,n1,j1,gox) 1087 1088 if (goz.and.goy.and.gox) then 1089 ind=j1+(j2-1)*n1+(i3loc-1)*n1*n2 1090 x=(gridcart(1,ind)-rx) 1091 y=(gridcart(2,ind)-ry) 1092 z=(gridcart(3,ind)-rz) 1093 r2=x**2+y**2+z**2 1094 xp=exp(-0.5_dp*r2/rloc**2) 1095 1096 ! Short range part 1097 rhoel=rho(ind) 1098 ! HGH: V_S^prime=gaussian 1099 if (usepaw==0) then 1100 if (nloc/=0) then 1101 arg=r2/rloc**2 1102 tt=cprime(nloc) 1103 do iloc=nloc-1,1,-1 1104 tt=arg*tt+cprime(iloc) 1105 end do 1106 forceloc=xp*tt*rhoel 1107 else 1108 forceloc=zero 1109 end if 1110 ! PAW: V_PAW^prime-V_L^prime 1111 else 1112 rr(1)=sqrt(r2) 1113 if (rr(1)>=rzero) then 1114 call paw_splint_der(msz,rad,vloc,d2vloc,1,rr,dvpawdr) 1115 call calcdVloc_mklocl(dvhgh,rr(1),rloc,charge) 1116 forceloc=rhoel*rloc2*(dvpawdr(1)-dvhgh)/rr(1) 1117 else 1118 forceloc=rhoel*rloc2*dvloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc) 1119 end if 1120 end if 1121 1122 fxgau=fxgau+forceloc*x 1123 fygau=fygau+forceloc*y 1124 fzgau=fzgau+forceloc*z 1125 1126 ! Long range part: error function 1127 Vel=pot(ind) 1128 fxerf=fxerf+xp*Vel*x 1129 fyerf=fyerf+xp*Vel*y 1130 fzerf=fzerf+xp*Vel*z 1131 1132 else if (nloc>0) then 1133 r2=x**2+y**2+z**2 1134 arg=r2/rloc**2 1135 xp=exp(-0.5_dp*arg) 1136 tt=cprime(nloc) 1137 do iloc=nloc-1,1,-1 1138 tt=arg*tt+cprime(iloc) 1139 end do 1140 forceleaked=forceleaked+xp*(1._dp+tt) 1141 end if 1142 end do 1143 end do 1144 end if 1145 end do 1146 1147 ! Final result of the forces 1148 floc(1,iat)=(hxh*hyh*hzh*prefactor)*fxerf+(hxh*hyh*hzh/rloc**2)*fxgau 1149 floc(2,iat)=(hxh*hyh*hzh*prefactor)*fyerf+(hxh*hyh*hzh/rloc**2)*fygau 1150 floc(3,iat)=(hxh*hyh*hzh*prefactor)*fzerf+(hxh*hyh*hzh/rloc**2)*fzgau 1151 1152 end do 1153 1154 forceleaked=forceleaked*prefactor*hxh*hyh*hzh 1155 if (iproc.eq.0) write(std_out,'(a,1pe12.5)') 'done. Leaked force: ',forceleaked 1156 1157 CONTAINS
mklocl_realspace/vloc_zero_mklocl [ Functions ]
[ Top ] [ mklocl_realspace ] [ Functions ]
NAME
vloc_zero_mklocl
FUNCTION
Use a quadratic interpolation to get limit of Vloc(x) at x->0
INPUTS
OUTPUT
SOURCE
912 function vloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc) 913 914 915 !Arguments ------------------------------------ 916 !scalars 917 integer,intent(in) :: msz 918 real(dp) :: vloc_zero_mklocl 919 real(dp),intent(in) :: charge,rloc 920 !arrays 921 real(dp) :: rad(msz),vloc(msz),d2vloc(msz) 922 923 !Local variables------------------------------- 924 !scalars 925 real(dp) :: y1,y2,y3,zz=0._dp 926 !arrays 927 real(dp) :: ll(3),xx(3),yy(3) 928 929 ! ************************************************************************* 930 931 !Select 3 points x1,x2,x3 near 0 932 if (rad(1)>1.d-10) then 933 xx(1:3)=rad(1:3) 934 else 935 xx(1:3)=rad(2:4) 936 end if 937 938 !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x 939 call paw_splint(msz,rad,vloc,d2vloc,3,xx,yy) 940 call calcVloc_mklocl(y1,xx(1),rloc,charge) 941 call calcVloc_mklocl(y2,xx(2),rloc,charge) 942 call calcVloc_mklocl(y3,xx(3),rloc,charge) 943 yy(1)= yy(1)-y1 944 yy(2)= yy(2)-y2 945 yy(3)= yy(3)-y3 946 947 !Find a polynomial of the form (z=0): 948 !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z) 949 950 !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3)) 951 ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3))) 952 !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3)) 953 ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3))) 954 !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2)) 955 ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2))) 956 957 vloc_zero_mklocl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3) 958 959 end function vloc_zero_mklocl
mklocl_wavelets/calcdVloc_wvl [ Functions ]
[ Top ] [ mklocl_wavelets ] [ Functions ]
NAME
calcdVloc_wvl
FUNCTION
Compute 1st-derivative of long-range HGH local ionic potential (derf)
INPUTS
OUTPUT
SOURCE
1763 subroutine calcdVloc_wvl(yy,xx,rloc,Z) 1764 1765 !Arguments ------------------------------------ 1766 !scalars 1767 real(dp),intent(in) :: xx,rloc,Z 1768 real(dp),intent(out) :: yy 1769 1770 !Local variables------------------------------- 1771 !scalars 1772 real(dp):: arg,tt 1773 1774 ! ************************************************************************* 1775 1776 arg=xx/(sqrt(2._dp)*rloc) 1777 call derf_ab(tt,arg) 1778 yy=(Z/(xx**2))* ( tt - 2._dp/sqrt(pi)*arg*exp(-arg**2) ) 1779 1780 end subroutine calcdVloc_wvl
mklocl_wavelets/dvloc_zero_mklocl [ Functions ]
[ Top ] [ mklocl_wavelets ] [ Functions ]
NAME
dvloc_zero_mklocl
FUNCTION
Use a quadratic interpolation to get limit of (1/x).dVloc(x)/dx at x->0
INPUTS
OUTPUT
SOURCE
1210 function dvloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc) 1211 1212 1213 !Arguments ------------------------------------ 1214 !scalars 1215 integer,intent(in) :: msz 1216 real(dp) :: dvloc_zero_mklocl 1217 real(dp),intent(in) :: charge,rloc 1218 !arrays 1219 real(dp) :: rad(msz),vloc(msz),d2vloc(msz) 1220 1221 !Local variables------------------------------- 1222 !scalars 1223 real(dp) :: y1,y2,y3,zz=0._dp 1224 !arrays 1225 real(dp) :: ll(3),xx(3),yy(3) 1226 1227 ! ************************************************************************* 1228 1229 !Select 3 points x1,x2,x3 near 0 1230 if (rad(1)>1.d-10) then 1231 xx(1:3)=rad(1:3) 1232 else 1233 xx(1:3)=rad(2:4) 1234 end if 1235 1236 !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x 1237 call paw_splint_der(msz,rad,vloc,d2vloc,3,xx,yy) 1238 call calcdVloc_mklocl(y1,xx(1),rloc,charge) 1239 call calcdVloc_mklocl(y2,xx(2),rloc,charge) 1240 call calcdVloc_mklocl(y3,xx(3),rloc,charge) 1241 yy(1)=(yy(1)-y1)/xx(1) 1242 yy(2)=(yy(2)-y2)/xx(2) 1243 yy(3)=(yy(3)-y3)/xx(3) 1244 1245 !Find a polynomial of the form (z=0): 1246 !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z) 1247 1248 !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3)) 1249 ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3))) 1250 !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3)) 1251 ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3))) 1252 !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2)) 1253 ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2))) 1254 1255 dvloc_zero_mklocl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3) 1256 1257 end function dvloc_zero_mklocl
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
INPUTS
OUTPUT
SOURCE
1798 function dvloc_zero_wvl(charge,rloc,msz,rad,vloc,d2vloc) 1799 1800 1801 !Arguments ------------------------------------ 1802 !scalars 1803 integer,intent(in) :: msz 1804 real(dp) :: dvloc_zero_wvl 1805 real(dp),intent(in) :: charge,rloc 1806 !arrays 1807 real(dp) :: rad(msz),vloc(msz),d2vloc(msz) 1808 1809 !Local variables------------------------------- 1810 !scalars 1811 real(dp) :: y1,y2,y3,zz=0._dp 1812 !arrays 1813 real(dp) :: ll(3),xx(3),yy(3) 1814 1815 ! ************************************************************************* 1816 1817 !Select 3 points x1,x2,x3 near 0 1818 if (rad(1)>1.d-10) then 1819 xx(1:3)=rad(1:3) 1820 else 1821 xx(1:3)=rad(2:4) 1822 end if 1823 1824 !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x 1825 call paw_splint_der(msz,rad,vloc,d2vloc,3,xx,yy) 1826 call calcdVloc_wvl(y1,xx(1),rloc,charge) 1827 call calcdVloc_wvl(y2,xx(2),rloc,charge) 1828 call calcdVloc_wvl(y3,xx(3),rloc,charge) 1829 yy(1)=(yy(1)-y1)/xx(1) 1830 yy(2)=(yy(2)-y2)/xx(2) 1831 yy(3)=(yy(3)-y3)/xx(3) 1832 1833 !Find a polynomial of the form (z=0): 1834 !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z) 1835 1836 !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3)) 1837 ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3))) 1838 !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3)) 1839 ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3))) 1840 !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2)) 1841 ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2))) 1842 1843 dvloc_zero_wvl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3) 1844 1845 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
SOURCE
1553 subroutine local_forces_wvl(iproc,natom,rxyz,hxh,hyh,hzh,n1,n2,n3,n3pi,i3s,n1i,n2i,& 1554 & rho,pot,floc,wvl) 1555 1556 use defs_wvltypes 1557 #if defined HAVE_BIGDFT 1558 use BigDFT_API, only : PSPCODE_PAW,ind_positions 1559 #endif 1560 1561 !Arguments ------------------------------- 1562 !scalars 1563 integer,intent(in) :: i3s,iproc,n1,n1i,n2,n2i,n3,n3pi,natom 1564 real(dp),intent(in) :: hxh,hyh,hzh 1565 type(wvl_internal_type),intent(in) :: wvl 1566 !arrays 1567 real(dp),intent(in) :: rxyz(3,natom) 1568 real(dp),dimension(*),intent(in) :: rho,pot 1569 real(dp),intent(out) :: floc(3,natom) 1570 1571 !Local variables ------------------------- 1572 #if defined HAVE_BIGDFT 1573 !scalars 1574 integer :: i1,i2,i3,iat,iex,iey,iez,iloc,ind,isx,isy,isz,ityp 1575 integer :: j1,j2,j3,msz=0 1576 integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3,nloc 1577 logical :: perx,pery,perz,gox,goy,goz,USE_PAW 1578 real(dp) :: arg,charge,cutoff,dvhgh,forceloc,fxerf,fyerf,fzerf,fxgau,fygau,fzgau 1579 real(dp) :: forceleaked,prefactor,r2,rhoel,rloc,rloc2,rx,ry,rz,rzero,tt,vel,x,xp,y,z 1580 !arrays 1581 real(dp) :: cprime(4),dvpawdr(1),rr(1) 1582 real(dp),pointer :: psppar(:,:),rad(:),vloc(:),d2vloc(:) 1583 #endif 1584 1585 ! ********************************************************************* 1586 1587 #if defined HAVE_BIGDFT 1588 1589 if (iproc==0) write(std_out,'(a)',advance='no') ' Calculate local forces...' 1590 1591 !PAW or NCPP ? 1592 USE_PAW=any(wvl%atoms%npspcode==PSPCODE_PAW) 1593 1594 !Conditions for periodicity in the three directions 1595 perx=(wvl%atoms%astruct%geocode /= 'F') 1596 pery=(wvl%atoms%astruct%geocode == 'P') 1597 perz=(wvl%atoms%astruct%geocode /= 'F') 1598 call ext_buffers(perx,nbl1,nbr1) 1599 call ext_buffers(pery,nbl2,nbr2) 1600 call ext_buffers(perz,nbl3,nbr3) 1601 1602 forceleaked=zero 1603 1604 do iat=1,natom 1605 ityp=wvl%atoms%astruct%iatype(iat) 1606 1607 ! Coordinates of the center 1608 rx=rxyz(1,iat) 1609 ry=rxyz(2,iat) 1610 rz=rxyz(3,iat) 1611 1612 ! Initialization of the forces 1613 ! ion-electron term, error function part 1614 fxerf=zero 1615 fyerf=zero 1616 fzerf=zero 1617 ! ion-electron term, gaussian part 1618 fxgau=zero 1619 fygau=zero 1620 fzgau=zero 1621 1622 ! Building array of coefficients of the derivative of the gaussian part 1623 psppar => wvl%atoms%psppar(:,:,ityp) 1624 cprime(1)=2._dp*wvl%atoms%psppar(0,2,ityp)-wvl%atoms%psppar(0,1,ityp) 1625 cprime(2)=4._dp*wvl%atoms%psppar(0,3,ityp)-wvl%atoms%psppar(0,2,ityp) 1626 cprime(3)=6._dp*wvl%atoms%psppar(0,4,ityp)-wvl%atoms%psppar(0,3,ityp) 1627 cprime(4)=-wvl%atoms%psppar(0,4,ityp) 1628 1629 ! Determine number of local terms (HGH pot) 1630 nloc=0 1631 do iloc=1,4 1632 if (wvl%atoms%psppar(0,iloc,ityp).ne.zero) nloc=iloc 1633 end do 1634 1635 ! Some constants depending on the atom type 1636 rloc=wvl%atoms%psppar(0,0,ityp) ; rloc2=rloc**2 1637 charge=real(wvl%atoms%nelpsp(ityp),kind=dp) 1638 prefactor=charge/(2._dp*pi*sqrt(2._dp*pi)*rloc**5) 1639 1640 ! PAW specifics 1641 if (USE_PAW) then 1642 msz=wvl%rholoc%msz(ityp) 1643 rad => wvl%rholoc%rad(1:msz,ityp) 1644 vloc => wvl%rholoc%d(1:msz,3,ityp) 1645 d2vloc => wvl%rholoc%d(1:msz,4,ityp) 1646 rzero=rad(1);if (rzero<=1.d-10) rzero=rad(2) 1647 end if 1648 1649 ! Maximum extension of the gaussian 1650 cutoff=10._dp*rloc 1651 isx=floor((rx-cutoff)/hxh) 1652 isy=floor((ry-cutoff)/hyh) 1653 isz=floor((rz-cutoff)/hzh) 1654 iex=ceiling((rx+cutoff)/hxh) 1655 iey=ceiling((ry+cutoff)/hyh) 1656 iez=ceiling((rz+cutoff)/hzh) 1657 1658 ! Calculate the forces near the atom due to the gaussian 1659 ! and error function parts of the potential 1660 if (n3pi>0) then 1661 do i3=isz,iez 1662 z=real(i3,kind=dp)*hzh-rz 1663 call ind_positions(perz,i3,n3,j3,goz) 1664 j3=j3+nbl3+1 1665 do i2=isy,iey 1666 y=real(i2,kind=dp)*hyh-ry 1667 call ind_positions(pery,i2,n2,j2,goy) 1668 do i1=isx,iex 1669 x=real(i1,kind=dp)*hxh-rx 1670 call ind_positions(perx,i1,n1,j1,gox) 1671 1672 r2=x**2+y**2+z**2 1673 xp=exp(-0.5_dp*r2/rloc2) 1674 1675 if ((j3>=i3s.and.j3<=i3s+n3pi-1) .and. (gox.and.goy)) then 1676 ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s+1-1)*n1i*n2i 1677 1678 ! Short range part 1679 rhoel=rho(ind) 1680 ! HGH: V_S^prime=gaussian 1681 if (.not.USE_PAW) then 1682 if (nloc/=0) then 1683 arg=r2/rloc2 1684 tt=cprime(nloc) 1685 do iloc=nloc-1,1,-1 1686 tt=arg*tt+cprime(iloc) 1687 end do 1688 forceloc=xp*tt*rhoel 1689 else 1690 forceloc=zero 1691 end if 1692 ! PAW: V_PAW^prime-V_L^prime 1693 else 1694 rr(1)=sqrt(r2) 1695 1696 if (rr(1)>=rzero) then 1697 call paw_splint_der(msz,rad,vloc,d2vloc,1,rr,dvpawdr) 1698 call calcdVloc_wvl(dvhgh,rr(1),rloc,charge) 1699 forceloc=rhoel*rloc2*(dvpawdr(1)-dvhgh)/rr(1) 1700 else 1701 forceloc=rhoel*rloc2*dvloc_zero_wvl(charge,rloc,msz,rad,vloc,d2vloc) 1702 end if 1703 end if 1704 1705 fxgau=fxgau+forceloc*x 1706 fygau=fygau+forceloc*y 1707 fzgau=fzgau+forceloc*z 1708 1709 ! Long range part: error function 1710 vel=pot(ind) 1711 fxerf=fxerf+xp*vel*x 1712 fyerf=fyerf+xp*vel*y 1713 fzerf=fzerf+xp*vel*z 1714 1715 else if ((.not.goz).and.(nloc>0)) then 1716 arg=r2/rloc2 1717 tt=cprime(nloc) 1718 do iloc=nloc-1,1,-1 1719 tt=arg*tt+cprime(iloc) 1720 end do 1721 forceleaked=forceleaked+prefactor*xp*tt*rho(1) 1722 end if 1723 1724 end do ! i1 1725 end do ! i2 1726 end do ! i3 1727 end if ! n3pi>0 1728 1729 ! Final result of the forces 1730 floc(1,iat)=(hxh*hyh*hzh*prefactor)*fxerf+(hxh*hyh*hzh/rloc2)*fxgau 1731 floc(2,iat)=(hxh*hyh*hzh*prefactor)*fyerf+(hxh*hyh*hzh/rloc2)*fygau 1732 floc(3,iat)=(hxh*hyh*hzh*prefactor)*fzerf+(hxh*hyh*hzh/rloc2)*fzgau 1733 1734 end do ! iat 1735 1736 forceleaked=forceleaked*hxh*hyh*hzh 1737 if (iproc.eq.0) write(std_out,'(a,1pe12.5)') 'done. Leaked force: ',forceleaked 1738 1739 #else 1740 BIGDFT_NOTENABLED_ERROR() 1741 if (.false.) write(std_out,*) i3s,iproc,n1,n1i,n2,n2i,n3,n3pi,natom,hxh,hyh,hzh,& 1742 & rxyz(1,1),floc(1,1),rho(1),pot(1),wvl%h(1) 1743 #endif 1744 1745 CONTAINS