TABLE OF CONTENTS
- ABINIT/mklocl_realspace
- 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/dvloc_zero_mklocl
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 usefull 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
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR,TRangel) 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
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.
SIDE EFFECTS
PARENTS
mklocl
CHILDREN
SOURCE
57 #if defined HAVE_CONFIG_H 58 #include "config.h" 59 #endif 60 61 #include "abi_common.h" 62 63 subroutine mklocl_realspace(grtn,icoulomb,mpi_enreg,natom,nattyp,nfft,ngfft,nscforder, & 64 & nspden,ntypat,option,pawtab,psps,rhog,rhor,rprimd,typat,& 65 & ucvol,usewvl,vpsp,xred) 66 67 use defs_basis 68 use defs_datatypes 69 use defs_abitypes 70 use m_xmpi 71 use m_profiling_abi 72 use m_errors 73 74 use m_geometry, only : xred2xcart 75 use m_mpinfo, only : ptabs_fourdp 76 use m_pawtab, only : pawtab_type 77 use m_paw_numeric, only : paw_splint 78 #if defined HAVE_BIGDFT 79 use BigDFT_API, only : coulomb_operator,deallocate_coulomb_operator 80 use defs_PSolver 81 #else 82 use defs_wvltypes, only : coulomb_operator 83 #endif 84 85 !This section has been created automatically by the script Abilint (TD). 86 !Do not modify the following lines by hand. 87 #undef ABI_FUNC 88 #define ABI_FUNC 'mklocl_realspace' 89 use interfaces_18_timing 90 use interfaces_53_ffts 91 use interfaces_62_poisson 92 use interfaces_67_common, except_this_one => mklocl_realspace 93 !End of the abilint section 94 95 implicit none 96 97 !Arguments ------------------------------------ 98 !scalars 99 integer,intent(in) :: natom,nfft,nspden,ntypat,option 100 real(dp),intent(in) :: ucvol 101 type(MPI_type),intent(in) :: mpi_enreg 102 type(pseudopotential_type),intent(in) :: psps 103 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 104 !arrays 105 integer,intent(in) :: icoulomb,nscforder,usewvl 106 integer,intent(in) :: nattyp(ntypat),ngfft(18),typat(natom) 107 real(dp),intent(in) :: rhog(2,nfft) 108 real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3) 109 real(dp),intent(in) :: xred(3,natom) 110 real(dp),intent(out) :: grtn(3,natom),vpsp(nfft) 111 112 !Local variables------------------------------- 113 character(len=1) :: geocode 114 !testing variables 115 !scalars 116 integer,parameter :: nStep=2 117 integer :: comm_fft,countParSeconde,i1,i2,i3 118 integer :: ia,ia1,ia2,igeo,ii,ind,itypat,ix,iy,iz,jj 119 integer :: kk,me_fft,n1,n2,n3,n3d,n_interpol 120 integer :: nproc_fft,tpsStart,tpsStop 121 real(dp),parameter :: min_rho_value=1.0d-12 122 real(dp) :: aa,bb,cc,dd,delta,deltaV,dr,dr2div6,invdr,r,vol_interpol,x,y,z,hgx,hgy,hgz,entmp 123 logical,parameter :: customRho=.false.,finiteDiff=.false.,testing=.false. 124 logical :: doIt 125 character(len=500) :: message 126 !arrays 127 integer :: ngfft_interpol(18) 128 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 129 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 130 real(dp) :: coord(3),coordXYZ(3),refValue(3),tsec(2) 131 real(dp),allocatable :: coordCart_interpol(:,:),coordRed_interpol(:,:) 132 real(dp),allocatable :: gridcart(:,:) 133 real(dp),allocatable :: grtn_cart_interpol(:,:),grtn_diff(:,:) 134 real(dp),allocatable :: rhog_interpol(:,:),rhog_testing(:,:),rhor_interpol(:) 135 real(dp),allocatable :: rhor_testing(:),rhor_work(:),xcart(:,:),vhartr(:),gxyz(:,:) 136 type(coulomb_operator):: kernel 137 138 ! ************************************************************************* 139 140 !Keep track of total time spent here 141 if (option==2) then 142 call timab(72,1,tsec) 143 end if 144 145 !Several constants (FFT sizes and parallelism) 146 n1 = ngfft(1) ; n2 = ngfft(2) ; n3 = ngfft(3) 147 nproc_fft = ngfft(10) ; me_fft = ngfft(11) 148 n3d = ngfft(13) !for parallel runs 149 if (nproc_fft==1) n3d=n3 !for serial runs 150 comm_fft=mpi_enreg%comm_fft 151 if(me_fft /= mpi_enreg%me_fft .or. nproc_fft /= mpi_enreg%nproc_fft) then 152 MSG_BUG("mpi_enreg%x_fft not equal to the corresponding values in ngfft") 153 end if 154 155 !Conditions for periodicity in the three directions 156 geocode='P' 157 if (icoulomb==1) geocode='F' 158 if (icoulomb==2) geocode='S' 159 160 !Get the distrib associated with this fft_grid 161 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 162 163 !Store xcart for each atom 164 ABI_ALLOCATE(xcart,(3, natom)) 165 call xred2xcart(natom, rprimd, xcart, xred) 166 !Store cartesian coordinates for each grid points 167 ABI_ALLOCATE(gridcart,(3, nfft)) 168 call mkgrid_fft(ffti3_local,fftn3_distrib,gridcart,nfft,ngfft,rprimd) 169 170 !Check whether all the PSP considered are of type GTH-HGH or PAW 171 doIt=.true. 172 do ii=1,psps%npsp 173 doIt=doIt .and.& 174 (psps%pspcod(ii)==2.or.psps%pspcod(ii)==3.or.psps%pspcod(ii)==10.or.psps%pspcod(ii)==7) 175 end do 176 177 !HGH-GTH/PAW treatment presumably starts here 178 if (doIt) then 179 180 ! Definition of the grid spacings as in the kernel routine 181 hgx = rprimd(1,1)/(n1) 182 hgy = rprimd(2,2)/(n2) 183 hgz = rprimd(3,3)/(n3) 184 185 186 !---------------------------------------------------------------------- 187 ! ----- Option 1: compute local ionic potential ----- 188 !---------------------------------------------------------------------- 189 if (option==1) then 190 191 call psolver_kernel( (/ hgx, hgy, hgz /), 2, icoulomb, me_fft, kernel, comm_fft, & 192 & (/n1,n2,n3/), nproc_fft, nscforder) 193 194 call createIonicPotential_new(fftn3_distrib,ffti3_local,& 195 & geocode,me_fft, nproc_fft, natom, & 196 & ntypat, typat, psps%gth_params%psppar, & 197 & int(psps%ziontypat), xcart,gridcart, hgx,hgy,hgz, & 198 & n1,n2,n3d,n3, kernel, vpsp, comm_fft,pawtab,psps%usepaw) 199 200 !---------------------------------------------------------------------- 201 ! ----- Option 2: compute forces induced by local ionic potential ----- 202 !---------------------------------------------------------------------- 203 else if (option == 2) then 204 205 ! Compute Hartree potential from rhor 206 ABI_ALLOCATE(vhartr,(nfft)) 207 call psolver_hartree(entmp, (/ hgx, hgy, hgz /), icoulomb, me_fft, comm_fft, nfft, & 208 & (/n1,n2,n3/), nproc_fft, nscforder, nspden, rhor, vhartr, usewvl) 209 210 ! Allocate temporary array for forces 211 ABI_ALLOCATE(gxyz,(3, natom)) 212 213 ! Calculate local part of the forces grtn (inspired from BigDFT routine) 214 call local_forces_new(fftn3_distrib,ffti3_local,geocode,me_fft, ntypat, natom, & 215 & typat, xcart, gridcart, psps%gth_params%psppar, int(psps%ziontypat), & 216 & hgx,hgy,hgz, n1,n2,n3,n3d, rhor,vhartr, gxyz, pawtab,psps%usepaw) 217 218 ! Forces should be in reduced coordinates. 219 do ia = 1, natom, 1 220 do igeo = 1, 3, 1 221 grtn(igeo, ia) = - rprimd(1, igeo) * gxyz(1, ia) & 222 & - rprimd(2, igeo) * gxyz(2, ia) & 223 & - rprimd(3, igeo) * gxyz(3, ia) 224 end do 225 end do 226 227 ! Deallocate local variables 228 ABI_DEALLOCATE(vhartr) 229 ABI_DEALLOCATE(gxyz) 230 end if 231 232 !---------------------------------------------------------------------- 233 ! ----- Section for the non-HGH/GTH/PAW pseudopotentials (testing) ---- 234 !---------------------------------------------------------------------- 235 else 236 237 if (testing) then 238 call system_clock(count_rate = countParSeconde) 239 call system_clock(tpsStart, count_rate = countParSeconde) 240 end if 241 242 ! dr is the r step in the sampling psps%vlspl 243 dr = psps%qgrid_vl(2) 244 invdr = 1._dp / dr 245 dr2div6 = dr * dr / 6._dp 246 247 if (option == 1) then 248 ! Set 0 in vpsp before summing 249 vpsp(:) = 0._dp 250 else if (option == 2) then 251 ! Allocate array to store cartesian gradient computed with 252 ! an interpolation of rhor 253 ABI_ALLOCATE(grtn_cart_interpol,(3, natom)) 254 grtn_cart_interpol(:, :) = 0._dp 255 256 n_interpol = nStep ** 3 257 ABI_ALLOCATE(coordRed_interpol,(3, nStep ** 3)) 258 ABI_ALLOCATE(coordCart_interpol,(3, nStep ** 3)) 259 260 if (testing .and. customRho) then 261 ! Use a custom rho instead of the self-consistent one. 262 ABI_ALLOCATE(rhor_testing,(nfft)) 263 ABI_ALLOCATE(rhog_testing,(2, nfft)) 264 end if 265 266 ABI_ALLOCATE(rhor_interpol,(nfft * n_interpol)) 267 ABI_ALLOCATE(rhor_work,(nfft * n_interpol)) 268 ABI_ALLOCATE(rhog_interpol,(2, nfft * n_interpol)) 269 270 if (testing .and. customRho) then 271 ! Testing only, changing rho with a centered gaussian 272 do ii = 1, nfft, 1 273 ! using the position of the first atom as center. 274 r = (gridcart(1, ii) - xcart(1, 1)) ** 2 + & 275 & (gridcart(2, ii) - xcart(2, 1)) ** 2 + & 276 & (gridcart(3, ii) - xcart(3, 1)) ** 2 277 rhor_testing(ii) = exp(-r/4._dp) 278 end do 279 ! Testing only, compute rhog_testing from rhor_testing 280 call fourdp(1,rhog_testing,rhor_testing,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0) 281 end if 282 283 ! Compute the interpolation of rho, using a fourier transform 284 rhog_interpol(:, :) = 0._dp 285 ii = 0 286 do i3 = 1, n3, 1 287 if (i3 <= n3 / 2) then 288 iz = i3 289 else 290 iz = n3 * nStep - n3 + i3 291 end if 292 do i2 = 1, n2, 1 293 if (i2 <= n2 / 2) then 294 iy = i2 295 else 296 iy = n2 * nStep - n2 + i2 297 end if 298 do i1 = 1, n1, 1 299 ii = ii + 1 300 if (i1 <= n1 / 2) then 301 ix = i1 302 else 303 ix = n1 * nStep - n1 + i1 304 end if 305 jj = (iz - 1) * n2 * n1 * nStep ** 2 + (iy - 1) * n3 * nStep + ix 306 if (testing .and. customRho) then 307 rhog_interpol(:, jj) = rhog_testing(:, ii) 308 else 309 rhog_interpol(:, jj) = rhog(:, ii) 310 end if 311 end do 312 end do 313 end do 314 315 ! Compute the interpolation of rho from the Fourier transformation 316 ngfft_interpol(:) = ngfft(:) 317 ngfft_interpol(1:3) = (/ n1 * nStep, n2 * nStep, n3 * nStep /) 318 ngfft_interpol(4:6) = (/ n1 * nStep + 1, n2 * nStep + 1, n3 * nStep /) 319 call fourdp(1,rhog_interpol,rhor_work,1,mpi_enreg,nfft*n_interpol,ngfft_interpol,mpi_enreg%paral_kgb,0) 320 321 ! Reorder rhor_interpol to be able to read it linearly 322 jj = 0 323 do i3 = 1, n3, 1 324 do i2 = 1, n2, 1 325 do i1 = 1, n1, 1 326 do iz = 1, nStep, 1 327 do iy = 1, nStep, 1 328 do ix = 1, nStep, 1 329 jj = jj + 1 330 kk = ((i3 - 1) * nStep + iz - 1) ! z coordinate in the interpolated grid 331 kk = kk * n1 * n2 * nStep ** 2 332 kk = kk + ((i2 - 1) * nStep + iy - 1) * n1 * nStep ! adding y coordinate 333 kk = kk + (i1 - 1) * nStep + ix ! adding x coordinate 334 rhor_interpol(jj) = rhor_work(kk) 335 end do 336 end do 337 end do 338 end do 339 end do 340 end do 341 ABI_DEALLOCATE(rhor_work) 342 343 ! Compute grid access in the interpolated volume 344 ii = 0 345 do iz = 1, nStep, 1 346 z = real(iz - 1, dp) / real(nStep, dp) 347 do iy = 1, nStep, 1 348 y = real(iy - 1, dp) / real(nStep, dp) 349 do ix = 1, nStep, 1 350 x = real(ix - 1, dp) / real(nStep, dp) 351 ii = ii + 1 352 coordRed_interpol(:, ii) = (/ x, y, z /) 353 ! Assuming orthogonal box (should be change later) 354 coordCart_interpol(:, ii) = (/ x * rprimd(1, 1) / real(n1, dp), & 355 & y * rprimd(2, 2) / real(n2, dp), & 356 & z * rprimd(3, 3) / real(n3, dp) /) 357 end do 358 end do 359 end do 360 361 vol_interpol = 1._dp / real(nStep, dp) ** 3 362 ! Compute the coordinates (integer) of each atom and deduce 363 ! the max extens of the integral summation. 364 ! !! do ia = 1, natom, 1 365 ! !! coordAtom(1, ia) = int(xred(1, ia) * n1) + 1 366 ! !! coordAtom(2, ia) = int(xred(2, ia) * n2) + 1 367 ! !! coordAtom(3, ia) = int(xred(3, ia) * n3) + 1 368 ! !! end do 369 end if 370 371 if (testing .and. option == 2) then 372 call system_clock(tpsStop, count_rate = countParSeconde) 373 write(std_out,*) "Tps : ", real(tpsStop - tpsStart) / real(countParSeconde) 374 end if 375 376 ia1=1 377 do itypat = 1, ntypat, 1 378 ! ia1,ia2 sets range of loop over atoms: 379 ia2 = ia1 + nattyp(itypat) - 1 380 381 do ii = 1, nfft, 1 382 do ia = ia1, ia2, 1 383 if (option == 1) then 384 ! Compute the potential 385 ! r is the distance between grid point and atom 386 r = sqrt((gridcart(1, ii) - xcart(1, ia)) ** 2 + & 387 & (gridcart(2, ii) - xcart(2, ia)) ** 2 + & 388 & (gridcart(3, ii) - xcart(3, ia)) ** 2) 389 390 ! Coefficients needed to compute the spline. 391 jj = int(r * invdr) + 1 392 if (jj > psps%mqgrid_vl - 2) then 393 write(message, '(3a,i0,a,i0,a,a)' )& 394 & ' pseudo-potential local part sampling is not wide enough', ch10, & 395 & ' want to access position ', jj, ' whereas mqgrid_vl = ', psps%mqgrid_vl, ch10, & 396 & ' Action : no idea, contact developpers...' 397 MSG_ERROR(message) 398 end if 399 delta = r - psps%qgrid_vl(jj) 400 bb = delta * invdr 401 aa = 1._dp - bb 402 cc = aa * (aa ** 2 - 1._dp) * dr2div6 403 dd = bb * (bb ** 2 - 1._dp) * dr2div6 404 405 ! compute V(r) from the spline, jj and jj + 1 is braketting r in 406 ! the sampling 407 deltaV = aa * psps%vlspl(jj, 1, itypat) + bb * psps%vlspl(jj + 1, 1, itypat) + & 408 & cc * psps%vlspl(jj, 2, itypat) + dd * psps%vlspl(jj + 1, 2, itypat) 409 ! Add on grid point ii the contribution of atom ia 410 vpsp(ii) = vpsp(ii) + deltaV 411 else if (option == 2) then 412 ! Compute the forces, as gradient of energy (V(r).rho(r)) 413 414 ! Testing only - reference points 415 if (.false.) then 416 ! r is the distance between grid point and atom 417 r = sqrt((gridcart(1, ii) - xcart(1, ia)) ** 2 + & 418 & (gridcart(2, ii) - xcart(2, ia)) ** 2 + & 419 & (gridcart(3, ii) - xcart(3, ia)) ** 2) 420 421 ! Coefficients needed to compute the spline. 422 jj = int(r * invdr) + 1 423 delta = r - psps%qgrid_vl(jj) 424 bb = delta * invdr 425 aa = 1._dp - bb 426 cc = aa * (aa ** 2 - 1._dp) * dr2div6 427 dd = bb * (bb ** 2 - 1._dp) * dr2div6 428 429 ! When mesh position is on a node, forces are null. 430 if (r /= 0._dp) then 431 ! This value deltaV is the first derivative of V(r) taken at r. 432 deltaV = aa * psps%dvlspl(jj, 1, itypat) + bb * psps%dvlspl(jj + 1, 1, itypat) + & 433 & cc * psps%dvlspl(jj, 2, itypat) + dd * psps%dvlspl(jj + 1, 2, itypat) 434 ! We multiply by rho(r) to have an energy. 435 deltaV = deltaV * rhor(ii, 1) / r 436 refValue(:) = - deltaV * (gridcart(:, ii) - xcart(:, ia)) 437 grtn_cart_interpol(:, ia) = grtn_cart_interpol(:, ia) + refValue(:) 438 end if 439 end if 440 441 ! Compute the interpolation for the point ii 442 ind = (ii - 1) * n_interpol 443 do kk = 1, n_interpol, 1 444 ind = ind + 1 445 446 if (rhor_interpol(ind) > min_rho_value) then 447 ! Assume orthogonal box... 448 coordXYZ(1) = gridcart(1, ii) - xcart(1, ia) + coordCart_interpol(1, kk) 449 coordXYZ(2) = gridcart(2, ii) - xcart(2, ia) + coordCart_interpol(2, kk) 450 coordXYZ(3) = gridcart(3, ii) - xcart(3, ia) + coordCart_interpol(3, kk) 451 r = coordXYZ(1) ** 2 + coordXYZ(2) ** 2 + coordXYZ(3) ** 2 452 453 if (r /= 0._dp) then 454 r = sqrt(r) 455 ! Coefficients needed to compute the spline. 456 jj = int(r * invdr) + 1 457 delta = r - psps%qgrid_vl(jj) 458 bb = delta * invdr 459 aa = 1._dp - bb 460 cc = aa * (aa ** 2 - 1._dp) * dr2div6 461 dd = bb * (bb ** 2 - 1._dp) * dr2div6 462 deltaV = aa * psps%dvlspl(jj, 1, itypat) + & 463 & bb * psps%dvlspl(jj + 1, 1, itypat) + & 464 & cc * psps%dvlspl(jj, 2, itypat) + & 465 & dd * psps%dvlspl(jj + 1, 2, itypat) 466 deltaV = deltaV * rhor_interpol(ind) / r 467 grtn_cart_interpol(1, ia) = grtn_cart_interpol(1, ia) - deltaV * coordXYZ(1) 468 grtn_cart_interpol(2, ia) = grtn_cart_interpol(2, ia) - deltaV * coordXYZ(2) 469 grtn_cart_interpol(3, ia) = grtn_cart_interpol(3, ia) - deltaV * coordXYZ(3) 470 ! do igeo = 1, 3, 1 471 ! grtn_cart_interpol(igeo, ia) = grtn_cart_interpol(igeo, ia) - deltaV * coordXYZ(igeo) 472 ! end do 473 end if 474 end if 475 end do 476 477 ! ============= 478 ! Testing only 479 ! ============= 480 ! use of finite differences 481 if (finiteDiff) then 482 do igeo = 1, 3, 1 483 coord(:) = 0._dp 484 coord(igeo) = dr / 2.0_dp 485 r = sqrt((gridcart(1, ii) - xcart(1, ia) + coord(1)) ** 2 + & 486 & (gridcart(2, ii) - xcart(2, ia) + coord(2)) ** 2 + & 487 & (gridcart(3, ii) - xcart(3, ia) + coord(3)) ** 2) 488 489 ! Coefficients needed to compute the spline. 490 jj = int(r * invdr) + 1 491 delta = r - psps%qgrid_vl(jj) 492 bb = delta * invdr 493 aa = 1._dp - bb 494 cc = aa * (aa ** 2 - 1._dp) * dr2div6 495 dd = bb * (bb ** 2 - 1._dp) * dr2div6 496 497 deltaV = aa * psps%vlspl(jj, 1, itypat) + bb * psps%vlspl(jj + 1, 1, itypat) + & 498 & cc * psps%vlspl(jj, 2, itypat) + dd * psps%vlspl(jj + 1, 2, itypat) 499 500 501 coord(:) = 0._dp 502 coord(igeo) = -dr / 2.0_dp 503 r = sqrt((gridcart(1, ii) - xcart(1, ia) + coord(1)) ** 2 + & 504 & (gridcart(2, ii) - xcart(2, ia) + coord(2)) ** 2 + & 505 & (gridcart(3, ii) - xcart(3, ia) + coord(3)) ** 2) 506 507 ! Coefficients needed to compute the spline. 508 jj = int(r * invdr) + 1 509 delta = r - psps%qgrid_vl(jj) 510 bb = delta * invdr 511 aa = 1._dp - bb 512 cc = aa * (aa ** 2 - 1._dp) * dr2div6 513 dd = bb * (bb ** 2 - 1._dp) * dr2div6 514 515 deltaV = deltaV - (aa * psps%vlspl(jj, 1, itypat) + & 516 & bb * psps%vlspl(jj + 1, 1, itypat) + & 517 & cc * psps%vlspl(jj, 2, itypat) + & 518 & dd * psps%vlspl(jj + 1, 2, itypat)) 519 grtn_diff(igeo, ia) = grtn_diff(igeo, ia) - deltaV * rhor(ii, 1) / dr 520 end do 521 end if 522 ! ============= 523 ! Testing only 524 ! ============= 525 526 end if 527 end do 528 ! End loop over atoms of type itypat 529 end do 530 ! End loop over real space grid points 531 532 ia1 = ia2 + 1 533 end do 534 ! End loop over type of atoms 535 536 if(option==2)then 537 ! multiply the forces by the volume of a single box mesh. 538 grtn_cart_interpol(:, :) = grtn_cart_interpol(:, :) * & 539 & ucvol / real(n1 * n2 * n3, dp) * vol_interpol 540 ! Transform cartesian forces to reduce coordinates 541 do ia = 1, natom, 1 542 do igeo = 1, 3, 1 543 grtn(igeo, ia) = rprimd(1, igeo) * grtn_cart_interpol(1, ia) + & 544 & rprimd(2, igeo) * grtn_cart_interpol(2, ia) + & 545 & rprimd(3, igeo) * grtn_cart_interpol(3, ia) 546 end do 547 end do 548 ABI_DEALLOCATE(rhor_interpol) 549 ABI_DEALLOCATE(rhog_interpol) 550 ABI_DEALLOCATE(coordRed_interpol) 551 ABI_DEALLOCATE(coordCart_interpol) 552 if (testing .and. customRho) then 553 ABI_DEALLOCATE(rhor_testing) 554 ABI_DEALLOCATE(rhog_testing) 555 end if 556 557 if (testing) then 558 call system_clock(tpsStop, count_rate = countParSeconde) 559 write(std_out,*) "Tps : ", real(tpsStop - tpsStart) / real(countParSeconde) 560 write(std_out,*) grtn_cart_interpol 561 MSG_ERROR("Testing section!") 562 end if 563 564 end if 565 566 !------------------------- 567 end if ! GTH/HGH/PAW psps 568 569 !Release temporary memory 570 ABI_DEALLOCATE(xcart) 571 ABI_DEALLOCATE(gridcart) 572 573 !Close timing counters 574 if (option==2)then 575 call timab(72,2,tsec) 576 end if 577 578 end subroutine mklocl_realspace
mklocl_realspace/calcdVloc_mklocl [ Functions ]
[ Top ] [ mklocl_realspace ] [ Functions ]
NAME
calcdVloc_mklocl
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_realspace
CHILDREN
SOURCE
1246 subroutine calcdVloc_mklocl(yy,xx,rloc,Z) 1247 1248 use defs_basis 1249 1250 !This section has been created automatically by the script Abilint (TD). 1251 !Do not modify the following lines by hand. 1252 #undef ABI_FUNC 1253 #define ABI_FUNC 'calcdVloc_mklocl' 1254 use interfaces_43_wvl_wrappers 1255 !End of the abilint section 1256 1257 implicit none 1258 !Arguments ------------------------------------ 1259 !scalars 1260 real(dp),intent(in) :: xx,rloc,Z 1261 real(dp),intent(out) :: yy 1262 1263 !Local variables------------------------------- 1264 !scalars 1265 real(dp):: arg,tt 1266 1267 ! ************************************************************************* 1268 1269 arg=xx/(sqrt(2._dp)*rloc) 1270 call derf_ab(tt,arg) 1271 yy=(Z/(xx**2))* ( tt - 2._dp/sqrt(pi)*arg*exp(-arg**2) ) 1272 1273 end subroutine calcdVloc_mklocl
mklocl_realspace/calcVloc_mklocl [ Functions ]
[ Top ] [ mklocl_realspace ] [ Functions ]
NAME
calcVloc_mklocl
FUNCTION
COPYRIGHT
Copyright (C) 2013-2018 ABINIT group (TRangel) 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_realspace
CHILDREN
SOURCE
894 subroutine calcVloc_mklocl(yy,xx,rloc,Z) 895 896 use defs_basis 897 898 !This section has been created automatically by the script Abilint (TD). 899 !Do not modify the following lines by hand. 900 #undef ABI_FUNC 901 #define ABI_FUNC 'calcVloc_mklocl' 902 use interfaces_43_wvl_wrappers 903 !End of the abilint section 904 905 implicit none 906 !Arguments ------------------------------------ 907 !scalars 908 real(dp),intent(in) :: xx,rloc,Z 909 real(dp),intent(out) :: yy 910 911 !Local variables------------------------------- 912 !scalars 913 real(dp):: arg,tt 914 915 ! ************************************************************************* 916 917 arg=xx/(sqrt(2.0)*rloc) 918 call derf_ab(tt,arg) 919 yy=-Z/xx*tt 920 921 end subroutine calcVloc_mklocl
mklocl_realspace/createIonicPotential_new [ Functions ]
[ Top ] [ mklocl_realspace ] [ Functions ]
NAME
createIonicPotential_new
FUNCTION
INPUTS
OUTPUT
PARENTS
mklocl_realspace
CHILDREN
SOURCE
599 subroutine createIonicPotential_new(fftn3_distrib,ffti3_local,geocode,iproc,& 600 & nproc,nat,ntypes,iatype,psppar,nelpsp,rxyz,gridcart,& 601 & hxh,hyh,hzh,n1i,n2i,n3d,n3i,kernel,pot_ion,spaceworld,pawtab,usepaw) 602 603 use defs_basis 604 use defs_datatypes 605 use m_profiling_abi 606 use m_errors 607 use m_xmpi 608 use defs_wvltypes, only : coulomb_operator 609 use m_pawtab, only : pawtab_type 610 use m_paw_numeric, only : paw_splint 611 612 !This section has been created automatically by the script Abilint (TD). 613 !Do not modify the following lines by hand. 614 #undef ABI_FUNC 615 #define ABI_FUNC 'createIonicPotential_new' 616 use interfaces_67_common, except_this_one => createIonicPotential_new 617 !End of the abilint section 618 619 implicit none 620 621 !Arguments ------------------------------- 622 !scalars 623 integer, intent(in) :: iproc,nproc,ntypes,nat,n1i,n2i,n3i,n3d,spaceworld,usepaw 624 real(dp), intent(in) :: hxh,hyh,hzh 625 character(len=1), intent(in) :: geocode 626 type(coulomb_operator), intent(in) :: kernel 627 !arrays 628 integer, dimension(nat), intent(in) :: iatype 629 integer, dimension(ntypes), intent(in) :: nelpsp 630 integer, dimension(*), intent(in) ::fftn3_distrib,ffti3_local 631 real(dp), dimension(3,n1i*n2i*n3d), intent(in) :: gridcart 632 real(dp), dimension(0:4,0:6,ntypes), intent(in) :: psppar 633 real(dp), dimension(3,nat), intent(in) :: rxyz 634 real(dp), dimension(*), intent(inout) :: pot_ion 635 type(pawtab_type),intent(in) :: pawtab(ntypes*usepaw) 636 637 !Local variables ------------------------- 638 #if defined HAVE_BIGDFT 639 !scalars 640 integer :: iat,i1,i2,i3,j1,j2,j3,isx,isy,isz,iex,iey,iez,ierr,ityp 641 integer :: ind,nloc,iloc,i3loc,msz 642 logical :: gox,goy,goz,perx,pery,perz 643 real(dp) :: arg,charge,cutoff,ehart,eexcu,rholeaked,rholeaked_tot,rloc 644 real(dp) :: rx,ry,rz,rzero,r2,tt,tt_tot,vexcu,vhgh,x,xp,y,z 645 !arrays 646 real(dp) :: rr(1),vpaw(1) 647 real(dp),pointer :: rad(:),vloc(:),d2vloc(:) 648 #endif 649 650 ! ********************************************************************* 651 652 #if defined HAVE_BIGDFT 653 654 if(nproc<0)then 655 MSG_ERROR('nproc should not be negative') 656 end if 657 658 !Ionic charge (must be calculated for the PS active processes) 659 rholeaked=0._dp 660 !Ionic energy (can be calculated for all the processors) 661 662 !here we should insert the calculation of the ewald energy for the periodic BC case 663 !!! eion=0._dp 664 !!! do iat=1,nat 665 !!! ityp=iatype(iat) 666 !!! rx=rxyz(1,iat) 667 !!! ry=rxyz(2,iat) 668 !!! rz=rxyz(3,iat) 669 !!! ! ion-ion interaction 670 !!! do jat=1,iat-1 671 !!! dist=sqrt( (rx-rxyz(1,jat))**2+(ry-rxyz(2,jat))**2+(rz-rxyz(3,jat))**2 ) 672 !!! jtyp=iatype(jat) 673 !!! eion=eion+real(nelpsp(jtyp)*nelpsp(ityp),kind=dp)/dist 674 !!! enddo 675 !!! end do 676 !!! if (iproc.eq.0) write(std_out,'(1x,a,1pe22.14)') 'ion-ion interaction energy',eion 677 678 !Creates charge density arising from the ionic PSP cores 679 !the n3pi dimension indicates the number of planes trated by each processor in the FFT parallelisation 680 !for a plane wave treatment this value depends on whether the direct space is divided in planes or not 681 !I don't know this variable, which in the future must be inserted at the place of n3pi (LG) 682 !if n3pi=0 this means that the processors doesn't calculate anything 683 !if (n3pi >0 ) then 684 685 !conditions for periodicity in the three directions 686 perx=(geocode /= 'F') 687 pery=(geocode == 'P') 688 perz=(geocode /= 'F') 689 690 !this initialise the array to zero, it will work only if bigdft library is enabled 691 pot_ion(1:n1i*n2i*n3d)=zero 692 693 do iat=1,nat 694 ityp=iatype(iat) 695 rx=rxyz(1,iat) 696 ry=rxyz(2,iat) 697 rz=rxyz(3,iat) 698 699 rloc=psppar(0,0,ityp) 700 charge=real(nelpsp(ityp),kind=dp)/(2._dp*pi*sqrt(2._dp*pi)*rloc**3) 701 cutoff=10._dp*rloc 702 703 isx=floor((rx-cutoff)/hxh) 704 isy=floor((ry-cutoff)/hyh) 705 isz=floor((rz-cutoff)/hzh) 706 707 iex=ceiling((rx+cutoff)/hxh) 708 iey=ceiling((ry+cutoff)/hyh) 709 iez=ceiling((rz+cutoff)/hzh) 710 711 ! Calculate Ionic Density 712 ! using HGH parameters. 713 ! Eq. 1.104, T. Deutsch and L. Genovese, JDN. 12, 2011 714 do i3=isz,iez 715 z=real(i3,kind=dp)*hzh-rz 716 call ind_positions_mklocl(perz,i3,n3i,j3,goz) 717 if(fftn3_distrib(j3)==iproc) then 718 i3loc=ffti3_local(j3) 719 do i2=isy,iey 720 y=real(i2,kind=dp)*hyh-ry 721 call ind_positions_mklocl(pery,i2,n2i,j2,goy) 722 do i1=isx,iex 723 x=real(i1,kind=dp)*hxh-rx 724 call ind_positions_mklocl(perx,i1,n1i,j1,gox) 725 r2=x**2+y**2+z**2 726 if (goz .and. goy .and. gox ) then 727 ind=j1+(j2-1)*n1i+(i3loc-1)*n1i*n2i 728 r2=(gridcart(1,ind)-rx)**2+(gridcart(2,ind)-ry)**2+(gridcart(3,ind)-rz)**2 729 end if 730 arg=r2/rloc**2 731 xp=exp(-.5d0*arg) 732 if (goz .and. goy .and. gox ) then 733 pot_ion(ind)=pot_ion(ind)-xp*charge 734 else 735 rholeaked=rholeaked+xp*charge 736 end if 737 end do 738 end do 739 end if 740 end do 741 742 end do 743 744 !Check 745 tt=0._dp 746 do j3= 1,n3d 747 do i2= 1,n2i 748 do i1= 1,n1i 749 ind=i1+(i2-1)*n1i+(j3-1)*n1i*n2i 750 tt=tt+pot_ion(ind) 751 end do 752 end do 753 end do 754 755 tt=tt*hxh*hyh*hzh 756 rholeaked=rholeaked*hxh*hyh*hzh 757 758 call xmpi_sum(tt,tt_tot,spaceworld,ierr) 759 call xmpi_sum(rholeaked,rholeaked_tot,spaceworld,ierr) 760 761 if (iproc.eq.0) then 762 write(std_out,'(1x,a,f26.12,2x,1pe10.3)') & 763 & 'total ionic charge, leaked charge ',tt_tot,rholeaked_tot 764 end if 765 766 !Here the value of the datacode must be kept fixed 767 !there can be some problems when running this stuff in parallel, 768 ! if the ionic potential distribution does not agree with the 769 ! plane distribution which is supposed to hold for the Poisson Solver 770 call psolver(geocode,'D',iproc,nproc,n1i,n2i,n3i,0,hxh,hyh,hzh,& 771 & pot_ion,kernel%kernel,pot_ion,ehart,eexcu,vexcu,0._dp,.false.,1) 772 773 !Add the remaining short-range local terms 774 do iat=1,nat 775 ityp=iatype(iat) 776 777 rx=rxyz(1,iat) 778 ry=rxyz(2,iat) 779 rz=rxyz(3,iat) 780 781 ! determine number of local terms 782 rloc=psppar(0,0,ityp) 783 cutoff=10._dp*rloc 784 charge=real(nelpsp(ityp),kind=dp) 785 786 ! determine number of local terms (HGH pot) 787 nloc=0 788 do iloc=1,4 789 if (psppar(0,iloc,ityp).ne.0._dp) nloc=iloc 790 end do 791 792 ! PAW specifics 793 if (usepaw==1) then 794 msz=pawtab(ityp)%wvl%rholoc%msz 795 rad => pawtab(ityp)%wvl%rholoc%rad(1:msz) 796 vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,3) 797 d2vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,4) 798 rzero=rad(1);if (rzero<=1.d-10) rzero=rad(2) 799 end if 800 801 isx=floor((rx-cutoff)/hxh) 802 isy=floor((ry-cutoff)/hyh) 803 isz=floor((rz-cutoff)/hzh) 804 805 iex=ceiling((rx+cutoff)/hxh) 806 iey=ceiling((ry+cutoff)/hyh) 807 iez=ceiling((rz+cutoff)/hzh) 808 809 do i3=isz,iez 810 z=real(i3,kind=dp)*hzh-rz 811 call ind_positions_mklocl(perz,i3,n3i,j3,goz) 812 if(fftn3_distrib(j3) == iproc .and. goz) then !MPI 813 i3loc=ffti3_local(j3) 814 if (goz) then 815 do i2=isy,iey 816 y=real(i2,kind=dp)*hyh-ry 817 call ind_positions_mklocl(pery,i2,n2i,j2,goy) 818 if (goy) then 819 do i1=isx,iex 820 x=real(i1,kind=dp)*hxh-rx 821 call ind_positions_mklocl(perx,i1,n1i,j1,gox) 822 if (gox) then 823 ind=j1+(j2-1)*n1i+(i3loc-1)*n1i*n2i 824 r2=(gridcart(1,ind)-rx)**2+(gridcart(2,ind)-ry)**2+(gridcart(3,ind)-rz)**2 825 ! r2=x**2+y**2+z**2 826 827 ! HGH: V_S=gaussian potential of Eq. (9) in JCP 129, 014109(2008) 828 if (usepaw==0) then 829 if (nloc /= 0) then 830 arg=r2/rloc**2 831 xp=exp(-.5d0*arg) 832 tt=psppar(0,nloc,ityp) 833 do iloc=nloc-1,1,-1 834 tt=arg*tt+psppar(0,iloc,ityp) 835 end do 836 pot_ion(ind)=pot_ion(ind)+xp*tt 837 end if 838 839 ! PAW: V_PAW-V_L^HGH 840 else 841 rr(1)=sqrt(r2) 842 if (rr(1)>=rzero) then 843 call paw_splint(msz,rad,vloc,d2vloc,1,rr,vpaw) 844 call calcVloc_mklocl(vhgh,rr(1),rloc,charge) 845 pot_ion(ind)=pot_ion(ind)+vpaw(1)-vhgh 846 else 847 pot_ion(ind)=pot_ion(ind)+vloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc) 848 end if 849 end if 850 851 end if 852 end do 853 end if 854 end do 855 end if 856 end if 857 end do 858 859 end do !iat 860 #else 861 BIGDFT_NOTENABLED_ERROR() 862 if (.false.) write(std_out,*) geocode,iproc,nproc,ntypes,nat,n1i,n2i,n3i,n3d,spaceworld,usepaw,& 863 & hxh,hyh,hzh,iatype(1),nelpsp(1),fftn3_distrib(1),ffti3_local(1),gridcart(1,1),psppar(1,1,1),& 864 & rxyz(1,1),pot_ion(1),pawtab(1)%mesh_size,kernel%co 865 #endif 866 867 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
PARENTS
mklocl_realspace
CHILDREN
SOURCE
1382 subroutine ind_positions_mklocl(periodic,i,n,j,go) 1383 1384 use m_profiling_abi 1385 1386 !This section has been created automatically by the script Abilint (TD). 1387 !Do not modify the following lines by hand. 1388 #undef ABI_FUNC 1389 #define ABI_FUNC 'ind_positions_mklocl' 1390 !End of the abilint section 1391 1392 implicit none 1393 1394 !Arguments ------------------------------- 1395 logical, intent(in) :: periodic 1396 integer, intent(in) :: i,n 1397 logical, intent(out) :: go 1398 integer, intent(out) :: j 1399 1400 ! ********************************************************************* 1401 1402 if (periodic) then 1403 go=.true. 1404 j=modulo(i-1,n)+1 1405 else 1406 j=i 1407 go=(i >= 1 .and. i <= n) 1408 end if 1409 1410 end subroutine ind_positions_mklocl
mklocl_realspace/local_forces_new [ Functions ]
[ Top ] [ mklocl_realspace ] [ Functions ]
NAME
local_forces_new
FUNCTION
INPUTS
OUTPUT
PARENTS
mklocl_realspace
CHILDREN
SOURCE
1027 subroutine local_forces_new(fftn3_distrib,ffti3_local,& 1028 geocode,iproc,ntypes,nat,iatype,rxyz,gridcart,psppar,nelpsp,hxh,hyh,hzh,& 1029 n1,n2,n3,n3d,rho,pot,floc,pawtab,usepaw) 1030 1031 use defs_basis 1032 use m_profiling_abi 1033 use m_pawtab, only : pawtab_type 1034 use m_paw_numeric, only : paw_splint_der 1035 1036 !This section has been created automatically by the script Abilint (TD). 1037 !Do not modify the following lines by hand. 1038 #undef ABI_FUNC 1039 #define ABI_FUNC 'local_forces_new' 1040 use interfaces_67_common, except_this_one => local_forces_new 1041 !End of the abilint section 1042 1043 implicit none 1044 1045 !Arguments ------------------------------- 1046 !scalars 1047 integer, intent(in) :: iproc,ntypes,nat,n1,n2,n3,n3d,usepaw 1048 character(len=1), intent(in) :: geocode 1049 real(dp), intent(in) :: hxh,hyh,hzh 1050 !arrays 1051 integer, dimension(*), intent(in) ::fftn3_distrib,ffti3_local 1052 integer, dimension(nat), intent(in) :: iatype 1053 integer, dimension(ntypes), intent(in) :: nelpsp 1054 real(dp), dimension(3,n1*n2*n3d), intent(in) :: gridcart 1055 real(dp), dimension(0:4,0:6,ntypes), intent(in) :: psppar 1056 real(dp), dimension(3,nat), intent(in) :: rxyz 1057 real(dp), dimension(*), intent(in) :: rho,pot 1058 real(dp), dimension(3,nat), intent(out) :: floc 1059 type(pawtab_type),intent(in) :: pawtab(ntypes*usepaw) 1060 1061 !Local variables ------------------------- 1062 !scalars 1063 integer :: isx,isy,isz,iex,iey,iez,i1,i2,i3,j1,j2,j3,ind,iat,ityp,iloc,i3loc,msz,nloc 1064 logical :: perx,pery,perz,gox,goy,goz 1065 real(dp) :: arg,charge,cutoff,dvhgh,fxerf,fyerf,fzerf,fxgau,fygau,fzgau,forceleaked 1066 real(dp) :: forceloc,prefactor,rloc,rloc2,rhoel,rx,ry,rz,rzero,r2,tt,x,xp,y,z,Vel 1067 real(dp), dimension(4) :: cprime 1068 !arrays 1069 real(dp) :: dvpawdr(1),rr(1) 1070 real(dp),pointer :: rad(:),vloc(:),d2vloc(:) 1071 1072 ! ********************************************************************* 1073 1074 if (iproc == 0) write(std_out,'(1x,a)',advance='no')'Calculate local forces...' 1075 1076 !Conditions for periodicity in the three directions 1077 perx=(geocode /= 'F') 1078 pery=(geocode == 'P') 1079 perz=(geocode /= 'F') 1080 1081 forceleaked=zero 1082 1083 do iat=1,nat 1084 ityp=iatype(iat) 1085 ! Coordinates of the center 1086 rx=rxyz(1,iat) 1087 ry=rxyz(2,iat) 1088 rz=rxyz(3,iat) 1089 1090 ! Initialization of the forces 1091 ! ion-electron term, error function part 1092 fxerf=zero 1093 fyerf=zero 1094 fzerf=zero 1095 ! ion-electron term, gaussian part 1096 fxgau=zero 1097 fygau=zero 1098 fzgau=zero 1099 1100 ! Building array of coefficients of the derivative of the gaussian part 1101 cprime(1)=2._dp*psppar(0,2,ityp)-psppar(0,1,ityp) 1102 cprime(2)=4._dp*psppar(0,3,ityp)-psppar(0,2,ityp) 1103 cprime(3)=6._dp*psppar(0,4,ityp)-psppar(0,3,ityp) 1104 cprime(4)=-psppar(0,4,ityp) 1105 1106 ! Determine number of local terms (HGH pot) 1107 nloc=0 1108 do iloc=1,4 1109 if (psppar(0,iloc,ityp).ne.zero) nloc=iloc 1110 end do 1111 1112 ! Some constants depending on the atom type 1113 rloc=psppar(0,0,ityp) ; rloc2=rloc**2 1114 charge=real(nelpsp(ityp),kind=dp) 1115 prefactor=charge/(2._dp*pi*sqrt(2._dp*pi)*rloc**5) 1116 1117 ! PAW specifics 1118 if (usepaw==1) then 1119 msz=pawtab(ityp)%wvl%rholoc%msz 1120 rad => pawtab(ityp)%wvl%rholoc%rad(1:msz) 1121 vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,3) 1122 d2vloc => pawtab(ityp)%wvl%rholoc%d(1:msz,4) 1123 rzero=rad(1);if (rad(1)<=1.d-10) rzero=rad(2) 1124 end if 1125 1126 ! Maximum extension of the gaussian 1127 cutoff=10._dp*rloc 1128 isx=floor((rx-cutoff)/hxh) 1129 isy=floor((ry-cutoff)/hyh) 1130 isz=floor((rz-cutoff)/hzh) 1131 iex=ceiling((rx+cutoff)/hxh) 1132 iey=ceiling((ry+cutoff)/hyh) 1133 iez=ceiling((rz+cutoff)/hzh) 1134 1135 ! Calculate the forces near the atom due to the gaussian 1136 ! and error function parts of the potential 1137 do i3=isz,iez 1138 z=real(i3,kind=dp)*hzh-rz 1139 call ind_positions_mklocl(perz,i3,n3,j3,goz) 1140 if(fftn3_distrib(j3)==iproc) then 1141 i3loc=ffti3_local(j3) 1142 do i2=isy,iey 1143 y=real(i2,kind=dp)*hyh-ry 1144 call ind_positions_mklocl(pery,i2,n2,j2,goy) 1145 do i1=isx,iex 1146 x=real(i1,kind=dp)*hxh-rx 1147 call ind_positions_mklocl(perx,i1,n1,j1,gox) 1148 1149 if (goz.and.goy.and.gox) then 1150 ind=j1+(j2-1)*n1+(i3loc-1)*n1*n2 1151 x=(gridcart(1,ind)-rx) 1152 y=(gridcart(2,ind)-ry) 1153 z=(gridcart(3,ind)-rz) 1154 r2=x**2+y**2+z**2 1155 xp=exp(-0.5_dp*r2/rloc**2) 1156 1157 ! Short range part 1158 rhoel=rho(ind) 1159 ! HGH: V_S^prime=gaussian 1160 if (usepaw==0) then 1161 if (nloc/=0) then 1162 arg=r2/rloc**2 1163 tt=cprime(nloc) 1164 do iloc=nloc-1,1,-1 1165 tt=arg*tt+cprime(iloc) 1166 end do 1167 forceloc=xp*tt*rhoel 1168 else 1169 forceloc=zero 1170 end if 1171 ! PAW: V_PAW^prime-V_L^prime 1172 else 1173 rr(1)=sqrt(r2) 1174 if (rr(1)>=rzero) then 1175 call paw_splint_der(msz,rad,vloc,d2vloc,1,rr,dvpawdr) 1176 call calcdVloc_mklocl(dvhgh,rr(1),rloc,charge) 1177 forceloc=rhoel*rloc2*(dvpawdr(1)-dvhgh)/rr(1) 1178 else 1179 forceloc=rhoel*rloc2*dvloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc) 1180 end if 1181 end if 1182 1183 fxgau=fxgau+forceloc*x 1184 fygau=fygau+forceloc*y 1185 fzgau=fzgau+forceloc*z 1186 1187 ! Long range part: error function 1188 Vel=pot(ind) 1189 fxerf=fxerf+xp*Vel*x 1190 fyerf=fyerf+xp*Vel*y 1191 fzerf=fzerf+xp*Vel*z 1192 1193 else if (nloc>0) then 1194 r2=x**2+y**2+z**2 1195 arg=r2/rloc**2 1196 xp=exp(-0.5_dp*arg) 1197 tt=cprime(nloc) 1198 do iloc=nloc-1,1,-1 1199 tt=arg*tt+cprime(iloc) 1200 end do 1201 forceleaked=forceleaked+xp*(1._dp+tt) 1202 end if 1203 end do 1204 end do 1205 end if 1206 end do 1207 1208 ! Final result of the forces 1209 floc(1,iat)=(hxh*hyh*hzh*prefactor)*fxerf+(hxh*hyh*hzh/rloc**2)*fxgau 1210 floc(2,iat)=(hxh*hyh*hzh*prefactor)*fyerf+(hxh*hyh*hzh/rloc**2)*fygau 1211 floc(3,iat)=(hxh*hyh*hzh*prefactor)*fzerf+(hxh*hyh*hzh/rloc**2)*fzgau 1212 1213 end do 1214 1215 forceleaked=forceleaked*prefactor*hxh*hyh*hzh 1216 if (iproc.eq.0) write(std_out,'(a,1pe12.5)') 'done. Leaked force: ',forceleaked 1217 1218 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
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
947 function vloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc) 948 949 use defs_basis 950 951 !This section has been created automatically by the script Abilint (TD). 952 !Do not modify the following lines by hand. 953 #undef ABI_FUNC 954 #define ABI_FUNC 'vloc_zero_mklocl' 955 !End of the abilint section 956 957 implicit none 958 959 !Arguments ------------------------------------ 960 !scalars 961 integer,intent(in) :: msz 962 real(dp) :: vloc_zero_mklocl 963 real(dp),intent(in) :: charge,rloc 964 !arrays 965 real(dp) :: rad(msz),vloc(msz),d2vloc(msz) 966 967 !Local variables------------------------------- 968 !scalars 969 real(dp) :: y1,y2,y3,zz=0._dp 970 !arrays 971 real(dp) :: ll(3),xx(3),yy(3) 972 973 ! ************************************************************************* 974 975 !Select 3 points x1,x2,x3 near 0 976 if (rad(1)>1.d-10) then 977 xx(1:3)=rad(1:3) 978 else 979 xx(1:3)=rad(2:4) 980 end if 981 982 !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x 983 call paw_splint(msz,rad,vloc,d2vloc,3,xx,yy) 984 call calcVloc_mklocl(y1,xx(1),rloc,charge) 985 call calcVloc_mklocl(y2,xx(2),rloc,charge) 986 call calcVloc_mklocl(y3,xx(3),rloc,charge) 987 yy(1)= yy(1)-y1 988 yy(2)= yy(2)-y2 989 yy(3)= yy(3)-y3 990 991 !Find a polynomial of the form (z=0): 992 !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z) 993 994 !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3)) 995 ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3))) 996 !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3)) 997 ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3))) 998 !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2)) 999 ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2))) 1000 1001 vloc_zero_mklocl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3) 1002 1003 end function vloc_zero_mklocl
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
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
1299 function dvloc_zero_mklocl(charge,rloc,msz,rad,vloc,d2vloc) 1300 1301 use defs_basis 1302 1303 !This section has been created automatically by the script Abilint (TD). 1304 !Do not modify the following lines by hand. 1305 #undef ABI_FUNC 1306 #define ABI_FUNC 'dvloc_zero_mklocl' 1307 !End of the abilint section 1308 1309 implicit none 1310 1311 !Arguments ------------------------------------ 1312 !scalars 1313 integer,intent(in) :: msz 1314 real(dp) :: dvloc_zero_mklocl 1315 real(dp),intent(in) :: charge,rloc 1316 !arrays 1317 real(dp) :: rad(msz),vloc(msz),d2vloc(msz) 1318 1319 !Local variables------------------------------- 1320 !scalars 1321 real(dp) :: y1,y2,y3,zz=0._dp 1322 !arrays 1323 real(dp) :: ll(3),xx(3),yy(3) 1324 1325 ! ************************************************************************* 1326 1327 !Select 3 points x1,x2,x3 near 0 1328 if (rad(1)>1.d-10) then 1329 xx(1:3)=rad(1:3) 1330 else 1331 xx(1:3)=rad(2:4) 1332 end if 1333 1334 !Find the corresponding values of y=(V^PAW(x)-V^HGH(x))/x 1335 call paw_splint_der(msz,rad,vloc,d2vloc,3,xx,yy) 1336 call calcdVloc_mklocl(y1,xx(1),rloc,charge) 1337 call calcdVloc_mklocl(y2,xx(2),rloc,charge) 1338 call calcdVloc_mklocl(y3,xx(3),rloc,charge) 1339 yy(1)=(yy(1)-y1)/xx(1) 1340 yy(2)=(yy(2)-y2)/xx(2) 1341 yy(3)=(yy(3)-y3)/xx(3) 1342 1343 !Find a polynomial of the form (z=0): 1344 !P(z) = y1.L1(z) + y2.L2(z) + y3.L3(z) 1345 1346 !L1(z) = (z-x2)(z-x3)/((x1-x2)(x1-x3)) 1347 ll(1)=(zz-xx(2))*(zz-xx(3))/((xx(1)-xx(2))*(xx(1)-xx(3))) 1348 !L2(z) = (z-x1)(z-x3)/((x2-x1)(x2-x3)) 1349 ll(2)=(zz-xx(1))*(zz-xx(3))/((xx(2)-xx(1))*(xx(2)-xx(3))) 1350 !L3(z) = (z-x1)(z-x2)/((x3-x1)(x3-x2)) 1351 ll(3)=(zz-xx(1))*(zz-xx(2))/((xx(3)-xx(1))*(xx(3)-xx(2))) 1352 1353 dvloc_zero_mklocl=yy(1)*ll(1)+yy(2)*ll(2)+yy(3)*ll(3) 1354 1355 end function dvloc_zero_mklocl