1MODULE tsvdw_module 2! 3!---------------------------------------------------------------------------------------------------------------- 4! TS-vdW Code Version 14.0 (RAD/BS, Princeton University, February 2013) 5!---------------------------------------------------------------------------------------------------------------- 6! All quantities necessary for the evaluation of the TS-vdW energy and forces are computed on the real-space 7! mesh using linear interpolation of the atomic pseudo-densities and their first derivatives which have been 8! mapped onto linear equispaced atomic grids from their original form computed on radial atomic grids via the 9! ATOMIC code. 10!---------------------------------------------------------------------------------------------------------------- 11! SYNOPSIS: radial form of rhoA & drhoA mapped onto linear grid; 12! atrho & rhosad on real-space mesh via linear interpolation; 13! integration on spherical atomic domains (subsets of real-space mesh); 14! quadratic veff derivatives computed linearly using sparse domain intersection algorithm. 15!---------------------------------------------------------------------------------------------------------------- 16! 17USE cell_base, ONLY: h !h matrix for converting between r and s coordinates via r = h s 18USE cell_base, ONLY: ainv !h^-1 matrix for converting between r and s coordinates via s = h^-1 r) 19USE cell_base, ONLY: omega !cell volume (in au^3) 20USE constants, ONLY: pi !pi in double-precision 21USE fft_base, ONLY: dfftp !FFT derived data type 22USE funct, ONLY: get_iexch !retrieves type of exchange utilized in functional 23USE funct, ONLY: get_icorr !retrieves type of correlation utilized in functional 24USE funct, ONLY: get_igcx !retrieves type of gradient correction to exchange utilized in functional 25USE funct, ONLY: get_igcc !retrieves type of gradient correction to correlation utilized in functional 26USE io_global, ONLY: stdout !print/write argument for standard output (to output file) 27USE ions_base, ONLY: nat !number of total atoms (all atomic species) 28USE ions_base, ONLY: nsp !number of unique atomic species 29USE ions_base, ONLY: na !number of atoms within each atomic species 30USE ions_base, ONLY: ityp !ityp(i):=type/species of ith atom 31USE ions_base, ONLY: atm !atm(j):=name of jth atomic species (3 characters) 32USE kinds, ONLY: DP !double-precision kind (selected_real_kind(14,200)) 33! the charge density is parallelized over the "band group" or processors 34USE mp_bands, ONLY: nproc_bgrp !number of processors 35USE mp_bands, ONLY: me_bgrp !processor number (0,1,...,nproc_bgrp-1) 36USE mp_bands, ONLY: intra_bgrp_comm !standard MPI communicator 37! atoms are parallelized over the "image group" 38USE mp_images, ONLY: nproc_image !number of processors 39USE mp_images, ONLY: me_image !processor number (0,1,...,nproc_image-1) 40USE mp_images, ONLY: intra_image_comm !standard MPI communicator 41USE mp, ONLY: mp_sum !MPI collection with sum 42USE parallel_include !MPI header 43USE uspp_param, ONLY: upf !atomic pseudo-potential data 44! 45IMPLICIT NONE 46! 47SAVE 48! 49! PUBLIC variables 50! 51LOGICAL, PUBLIC :: vdw_isolated ! isolated system control 52REAL(DP), PUBLIC:: vdw_econv_thr ! energy convergence threshold for periodic systems 53REAL(DP), PUBLIC :: EtsvdW !the TS-vdW energy 54REAL(DP), DIMENSION(:), ALLOCATABLE, PUBLIC :: UtsvdW !the TS-vdW wavefunction forces (dispersion potential) 55REAL(DP), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: FtsvdW !the TS-vdW ionic forces (-dE/dR) 56REAL(DP), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: HtsvdW !the TS-vdW cell forces (dE/dh) 57REAL(DP), DIMENSION(:), ALLOCATABLE, PUBLIC :: VefftsvdW !the TS-vdW effective Hirshfeld volume 58! 59! PRIVATE variables 60! 61INTEGER, PARAMETER, PRIVATE :: NgpA=1000 !number of grid points for linear equispaced atomic grids (current value=1000pts) 62INTEGER, PARAMETER, PRIVATE :: bsint=BIT_SIZE(NgpA) !integer bit size (for use in bit array manipulation) 63INTEGER, PRIVATE :: me !processor number (1,2,...,nproc_image) 64INTEGER, PRIVATE :: iproc !processor dummy index 65INTEGER, PRIVATE :: nr1,nr2,nr3 !real space grid dimensions (global first, second, and third dimensions of the 3D grid) 66INTEGER, PRIVATE :: nr1r,nr2r,nr3r !reduced real space grid dimensions (global first, second, and third dimensions of the 3D grid) 67REAL(DP), PRIVATE :: ddamp !damping function parameter #1 68REAL(DP), PRIVATE :: sR !damping function parameter #2 69REAL(DP), PRIVATE :: spcutAmax !maximum radial cutoff for all atomic species 70INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: nstates !number of atoms per processor 71INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: sdispls !send displacement (offset) array 72INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: rdispls !receive displacement (offset) array 73INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: sendcount !send count array 74INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: recvcount !receive count array 75INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: istatus !MPI status array 76INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: NsomegaA !number of points in the spherical atomic integration domain 77INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: NsomegaAr !number of points in the reduced spherical atomic integration domain 78INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: npair !number of unique atom pairs 79INTEGER, DIMENSION(:,:), ALLOCATABLE, PRIVATE :: pair !unique atom pair overlap matrix 80INTEGER, DIMENSION(:,:), ALLOCATABLE, PRIVATE :: gomegar !precursor to spherical atomic integration domain (intersection bit array) 81INTEGER, DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: somegaA !spherical atomic integration domain 82INTEGER, DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: somegaAr !reduced spherical atomic integration domain 83INTEGER, DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: gomegaAr !reduced spherical atomic integration domain (intersection bit array) 84REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE:: predveffAdn !atomic dispersion potential prefactor 85REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: vfree !free atomic volumes for each atomic species 86REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: dpfree !free atomic static dipole polarizability for each atomic species 87REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: R0free !free atomic vdW radius for each atomic species 88REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: C6AAfree !free atomic homonuclear C6 coefficient for each atomic species 89REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: veff !effective atomic volumes for each atom in the simulation cell 90REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: dpeff !effective atomic static dipole polarizability for each atom in the simulation cell 91REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: R0eff !effective atomic vdW radius for each atom in the simulation cell 92REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: C6AAeff !effective atomic homonuclear C6 coefficient for each atom in the simulation cell 93REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: rhosad !molecular pro-density (superposition of atomic densities) on real-space mesh 94REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: rhotot !molecular charge density on real-space mesh 95REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE:: dveffAdn !the local copy of the TS-vdW wavefunction forces (dispersion potential) 96REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: spgrd !linear equispaced grid for each atomic species 97REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: sprho !atomic pseudo-density for each atomic species 98REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: spdrho !first derivative of atomic pseudo-density for each atomic species 99REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: spdata !linear grid cutoff (is,1) and linear grid spacing (is,2) for each atomic species 100REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: LIA !A coefficient for linear interpolation of rhoA 101REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: LIB !B coefficient for linear interpolation of rhoA 102REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: dLIA !A coefficient for linear interpolation of drhoA 103REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: dLIB !B coefficient for linear interpolation of drhoA 104REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: atxyz !Cartesian coordinates of ions adjusted according to PBC 105REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: C6ABfree !free atomic heteronuclear C6 coefficient for each atom pair 106REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: C6ABeff !effective atomic heteronuclear C6 coefficient for each atom pair 107REAL(DP), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: dveffdR !first derivative of effective volume wrt nuclear displacement 108REAL(DP), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: dveffdh !first derivative of effective volume wrt cell displacement 109! 110! PUBLIC subroutines 111! 112PUBLIC :: tsvdw_initialize 113PUBLIC :: tsvdw_calculate 114PUBLIC :: tsvdw_finalize 115! 116! PRIVATE subroutines 117! 118PRIVATE :: tsvdw_para_init 119PRIVATE :: tsvdw_pbc 120PRIVATE :: tsvdw_unique_pair 121PRIVATE :: tsvdw_rhotot 122PRIVATE :: tsvdw_screen 123PRIVATE :: tsvdw_veff 124PRIVATE :: tsvdw_dveff 125PRIVATE :: tsvdw_effqnts 126PRIVATE :: tsvdw_energy 127PRIVATE :: tsvdw_wfforce 128PRIVATE :: tsvdw_cleanup 129PRIVATE :: Num1stDer 130PRIVATE :: CubSplCoeff 131PRIVATE :: GetVdWParam 132 ! 133 ! 134 CONTAINS 135 ! 136 ! 137 !-------------------------------------------------------------------------------------------------------------- 138 SUBROUTINE tsvdw_initialize() 139 !-------------------------------------------------------------------------------------------------------------- 140 ! 141 IMPLICIT NONE 142 ! 143 ! Local Variables 144 ! 145 LOGICAL :: uniform_grid=.FALSE. 146 INTEGER :: ip,iq,ir,is,it,NrgpA,NrgpintA,icutrA,Ndim 147 REAL(DP) :: dxA,gfctrA,vref,eref,verr,d,dk1,dk2,dk3,num,den,drab,f1,f2,f3,L1,L2,L3 148 REAL(DP), DIMENSION(:), ALLOCATABLE :: atgrdr,atgrdrab,atrhor,datrhor,d2atrhor,CSA,CSB,CSC,CSD 149 ! 150 ! Start of calculation banner... 151 ! 152 WRITE(stdout,*) 153 WRITE(stdout,'(3X,"TS-vdW initialization")') 154 WRITE(stdout,'(3X,"---------------------")') 155 WRITE(stdout,*) 156 ! 157 ! Error messages for inconsistencies with current version of code... 158 ! 159 !RAD: Have we missed any inconsistencies? 160 ! 161 ! Setup variables for use in TS-vdW module... 162 ! 163 nr1=dfftp%nr1; nr2=dfftp%nr2; nr3=dfftp%nr3 164 nr1r=nr1/2; nr2r=nr2/2; nr3r=nr3/2 165 IF(MOD(nr1,2).EQ.1) nr1r=(nr1+1)/2 166 IF(MOD(nr2,2).EQ.1) nr2r=(nr2+1)/2 167 IF(MOD(nr3,2).EQ.1) nr3r=(nr3+1)/2 168 ! 169 ! Initialize the TS-vdW ionic forces, cell forces, and dispersion potential (wavefunction forces)... 170 ! 171 ALLOCATE(FtsvdW(3,nat)); FtsvdW=0.0_DP 172 ALLOCATE(HtsvdW(3,3)); HtsvdW=0.0_DP 173 ! 174 ! Initialization of TS-vdW Hirshfeld effective volume public variable ... used in CP print_out.f90 175 ! 176 ALLOCATE(VefftsvdW(nat)); VefftsvdW=0.0_DP 177 ! 178 ALLOCATE(UtsvdW(dfftp%nnr)); UtsvdW=0.0_DP 179 ! 180 ! Set ddamp damping function parameter (set to 20 and functional independent)... 181 ! 182 WRITE(stdout,'(3X,"Determining TS-vdW damping function parameters...")') 183 ddamp=20.0_DP 184 WRITE(stdout,'(5X,"ddamp = ",F9.6)') ddamp 185 ! 186 ! Set sR damping function parameter (functional dependent and currently only available for PBE & PBE0)... 187 ! 188 IF (get_iexch().EQ.1.AND.get_icorr().EQ.4.AND.get_igcx().EQ.3.AND.get_igcc().EQ.4) THEN 189 ! 190 sR=0.94_DP !PBE=sla+pw+pbx+pbc 191 ! 192 ELSE IF (get_iexch().EQ.6.AND.get_icorr().EQ.4.AND.get_igcx().EQ.8.AND.get_igcc().EQ.4) THEN 193 ! 194 sR=0.96_DP !PBE0=pb0x+pw+pb0x+pbc !RAD/BS: This line will not work in CP unless PBE0 code update funct.f90... 195 ! 196 ELSE 197 ! 198 CALL errore('tsvdw','TS-vdW sR parameter only available for PBE and PBE0 functionals...',1) 199 ! 200 END IF 201 ! 202 WRITE(stdout,'(5X,"sR = ",F9.6)') sR 203 ! 204 ! Allocate and initialize species-specific quantities... 205 ! 206 ALLOCATE(vfree(nsp)); vfree=0.0_DP 207 ALLOCATE(dpfree(nsp)); dpfree=0.0_DP 208 ALLOCATE(R0free(nsp)); R0free=0.0_DP 209 ALLOCATE(C6AAfree(nsp)); C6AAfree=0.0_DP 210 ALLOCATE(C6ABfree(nsp,nsp)); C6ABfree=0.0_DP 211 ALLOCATE(spdata(nsp,2)); spdata=0.0_DP 212 ALLOCATE(spgrd(nsp,0:NgpA)); spgrd=0.0_DP 213 ALLOCATE(sprho(nsp,0:NgpA)); sprho=0.0_DP 214 ALLOCATE(spdrho(nsp,0:NgpA)); spdrho=0.0_DP 215 ALLOCATE(LIA(nsp,0:NgpA)); LIA=0.0_DP 216 ALLOCATE(LIB(nsp,0:NgpA)); LIB=0.0_DP 217 ALLOCATE(dLIA(nsp,0:NgpA)); dLIA=0.0_DP 218 ALLOCATE(dLIB(nsp,0:NgpA)); dLIB=0.0_DP 219 ! 220 spcutAmax=0.0_DP 221 ! 222 ! Loop over atomic species and extract species-dependent quantities to modular arrays... 223 ! 224 DO is=1,nsp 225 ! 226 ! Obtain the radial grid and radial atomic pseudo-density from pseudo-potential file (via upf module) for 227 ! the given atomic species. Convert the radial atomic pseudo-density to the real atomic pseudo-density using 228 ! rho_real(r) = rho_radial(r) / (4*pi*r^2)... 229 ! 230 WRITE(stdout,'(3X,"Initializing species # ",I3," with atomic symbol ",A3)') is,atm(is) 231 ! 232 ! Read in the number of grid points in radial mesh from upf... 233 ! 234 NrgpA=upf(is)%mesh 235 ! 236 ! Transfer radial atomic grid (in upf) to local atgrdr array... 237 ! 238 ALLOCATE(atgrdr(NrgpA)); atgrdr=0.0_DP 239 ! 240 DO ir=1,NrgpA 241 ! 242 atgrdr(ir)=upf(is)%r(ir) 243 ! 244 END DO 245 ! 246 ! Transfer radial atomic grid spacing (in upf) to local atgrdrab array... 247 ! 248 ALLOCATE(atgrdrab(NrgpA)); atgrdrab=0.0_DP 249 ! 250 DO ir=1,NrgpA 251 ! 252 atgrdrab(ir)=upf(is)%rab(ir) 253 ! 254 END DO 255 ! 256 ! Determine whether radial grid is logarithmic/exponential or equispaced/uniform... 257 ! 258 drab=atgrdrab(NrgpA)-atgrdrab(1) 259 uniform_grid = (DABS(drab).LT.(1.0E-6_DP)) 260 IF (uniform_grid) WRITE(stdout,'(5X,"Equispaced/Uniform radial atomic grid detected...")') 261 ! 262 ! ---------------------------------------------------------------- 263 ! Logarithmic/Exponential grid (3 parameters: zmesh, xmin, dxA) 264 ! ---------------------------------------------------------------- 265 ! 266 ! For i = 1,2,...,NrgpA: 267 ! r(i) = exp[xmin+(i-1)*dxA]/zmesh 268 ! = exp[xmin]/zmesh * exp[(i-1)*dxA] 269 ! = gfctrA * exp[(i-1)*dxA] 270 ! rab(i) = r(i) * dxA 271 ! 272 ! Assumptions: grid does NOT start from zero (use simpson_cp90()). 273 ! 274 ! --------------------------------------------- 275 ! Equispaced/Uniform grid (1 parameter: dxA) 276 ! --------------------------------------------- 277 ! 278 ! For i = 1,2,...,NrgpA: 279 ! r(i) = (i-1) * dxA 280 ! rab(i) = dxA 281 ! 282 ! Assumptions: grid starts from zero (use simpson() for integration). 283 ! 284 ! Determine atomic radial grid parameters... 285 ! 286 IF (uniform_grid.EQV..TRUE.) THEN 287 ! 288 gfctrA=1.0_DP 289 dxA=atgrdrab(1) 290 ! 291 ELSE 292 ! 293 gfctrA=upf(is)%r(1) 294 dxA=DLOG(upf(is)%r(2)/upf(is)%r(1)) 295 ! 296 END IF 297 ! 298 WRITE(stdout,'(5X,"Radial grid parameter: NrgpA is ",I5,".")') NrgpA 299 WRITE(stdout,'(5X,"Radial grid parameter: gfctrA is ",F9.6,".")') gfctrA 300 WRITE(stdout,'(5X,"Radial grid parameter: dxA is ",F9.6,".")') dxA 301 ! 302 ! Transfer radial atomic pseudo-density to atrhor array... 303 ! Convert radial atomic pseudo-density to real atomic pseudo-density [n(r) = nrad(r)/(4*pi*r^2)]... 304 ! 305 ALLOCATE(atrhor(NrgpA)); atrhor=0.0_DP 306 ! 307 IF (uniform_grid.EQV..TRUE.) THEN 308 ! 309 DO ir=2,NrgpA 310 ! 311 atrhor(ir)=(upf(is)%rho_at(ir))/(4.0_DP*pi*atgrdr(ir)**(2.0_DP)) ! skip point at r=0... 312 ! 313 END DO 314 ! 315 ! Quadratic extrapolation of the atomic density to r=0... 316 ! 317 L1=((0.0_DP-atgrdr(3))*(0.0_DP-atgrdr(4)))/((atgrdr(2)-atgrdr(3))*(atgrdr(2)-atgrdr(4))) 318 L2=((0.0_DP-atgrdr(2))*(0.0_DP-atgrdr(4)))/((atgrdr(3)-atgrdr(2))*(atgrdr(3)-atgrdr(4))) 319 L3=((0.0_DP-atgrdr(2))*(0.0_DP-atgrdr(3)))/((atgrdr(4)-atgrdr(2))*(atgrdr(4)-atgrdr(3))) 320 atrhor(1)=L1*atrhor(2)+L2*atrhor(3)+L3*atrhor(4) 321 ! 322 ELSE 323 ! 324 DO ir=1,NrgpA 325 ! 326 atrhor(ir)=(upf(is)%rho_at(ir))/(4.0_DP*pi*atgrdr(ir)**(2.0_DP)) 327 ! 328 END DO 329 ! 330 END IF 331 ! 332 ! Set NrgpintA as the number of grid points (which must be odd) used during numerical integration using Simpson's rule... 333 ! 334 IF (IAND(NrgpA,1).EQ.1) THEN 335 ! 336 NrgpintA=NrgpA 337 ! 338 ELSE 339 ! 340 NrgpintA=NrgpA-1 341 ! 342 END IF 343 ! 344 ! Compute the number of electrons (eref) for each atomic species via numerical integration 345 ! of the atomic pseudo-density on the radial atomic grid using Simpson's rule... 346 ! 347 eref=0.0_DP 348 ! 349 DO ir=1,NrgpintA-2,2 350 ! 351 f1=atrhor(ir )*atgrdrab(ir )*atgrdr(ir )**(2.0_DP) ! integrated quantity is rho 352 f2=atrhor(ir+1)*atgrdrab(ir+1)*atgrdr(ir+1)**(2.0_DP) 353 f3=atrhor(ir+2)*atgrdrab(ir+2)*atgrdr(ir+2)**(2.0_DP) 354 ! 355 eref=eref+(f1+4.0_DP*f2+f3) 356 ! 357 END DO 358 ! 359 eref=(4.0_DP*pi/3.0_DP)*eref 360 WRITE(stdout,'(5X,"The number of valence electrons, eref, is ",F25.15,".")') eref 361 ! 362 ! Compute the reference free atom volume (vref) for each atomic species via numerical integration 363 ! of the atomic pseudo-density on the radial atomic grid using Simpson's rule... 364 ! 365 vref=0.0_DP 366 ! 367 DO ir=1,NrgpintA-2,2 368 ! 369 f1=atrhor(ir )*atgrdrab(ir )*atgrdr(ir )**(5.0_DP) ! integrated quantity is rho * r^3 370 f2=atrhor(ir+1)*atgrdrab(ir+1)*atgrdr(ir+1)**(5.0_DP) 371 f3=atrhor(ir+2)*atgrdrab(ir+2)*atgrdr(ir+2)**(5.0_DP) 372 ! 373 vref=vref+(f1+4.0_DP*f2+f3) 374 ! 375 END DO 376 ! 377 vref=(4.0_DP*pi/3.0_DP)*vref 378 WRITE(stdout,'(5X,"The reference free atom volume, vref, is ",F25.15," bohr^3.")') vref 379 ! 380 ! Using the reference free atom volume, determine an acceptable radial grid cutoff value such that the 381 ! free atom volume obtained using this cutoff does not deviate from the reference value by more than 1.0%. 382 ! 383 WRITE(stdout,'(5X,"Determining intial radial grid cutoff...")') 384 ! 385 DO iq=5,NrgpintA,2 386 ! 387 vfree(is)=0.0_DP 388 verr=0.0_DP 389 ! 390 DO ir=1,iq-2,2 391 ! 392 f1=atrhor(ir )*atgrdrab(ir )*atgrdr(ir )**(5.0_DP) ! integrated quantity is rho * r^3 393 f2=atrhor(ir+1)*atgrdrab(ir+1)*atgrdr(ir+1)**(5.0_DP) 394 f3=atrhor(ir+2)*atgrdrab(ir+2)*atgrdr(ir+2)**(5.0_DP) 395 ! 396 vfree(is)=vfree(is)+(f1+4.0_DP*f2+f3) 397 ! 398 END DO 399 ! 400 vfree(is)=(4.0_DP*pi/3.0_DP)*vfree(is) 401 verr=(vref-vfree(is))/vref*100.0_DP 402 ! 403 IF (verr.LE.1.0_DP) THEN 404 ! 405 icutrA=iq 406 ! 407 WRITE(stdout,'(5X,"An acceptable radial grid cutoff was determined by retaining ",I4," of ",I4," radial grid points.")') & 408 icutrA,NrgpA 409 ! 410 EXIT 411 ! 412 END IF 413 ! 414 END DO 415 ! 416 WRITE(stdout,'(5X,"The magnitude of the atomic pseudo-density at the radial grid cutoff is ",ES13.6,".")') atrhor(icutrA) 417 WRITE(stdout,'(5X,"Using this radial grid cutoff value of ",F25.15," au:")') atgrdr(icutrA) 418 WRITE(stdout,'(5X,"The free atom volume computed with this cutoff is ",F25.15," bohr^3 with an error of ",F6.3,"%.")') & 419 vfree(is),verr 420 ! 421 ! Form 1st derivative of atrhor for input into cubic spline coefficient subroutine... 422 ! 423 ALLOCATE(datrhor(NrgpA)); datrhor=0.0_DP 424 CALL Num1stDer(atgrdr,atrhor,NrgpA,dxA,datrhor) 425 ! 426 ! For logarithmic/exponential grid, transform linear derivative back to radial grid... 427 ! 428 IF (.NOT.uniform_grid) THEN 429 ! 430 DO ir=1,NrgpA 431 ! 432 datrhor(ir)=datrhor(ir)/atgrdr(ir) 433 ! 434 END DO 435 ! 436 END IF 437 ! 438 ! Form the coefficients of the cubic spline interpolant (2nd derivatives) for the real atomic pseudo-density 439 ! for use during cubic spline interpolation of the pseudo-density onto the linear equispaced atomic grid... 440 ! 441 ALLOCATE(d2atrhor(NrgpA)); d2atrhor=0.0_DP 442 CALL CubSplCoeff(atgrdr,atrhor,NrgpA,datrhor,d2atrhor) 443 ! 444 ! Precompute cubic spline interpolation vectors (utilizing Taylor series form) via: 445 ! 446 ! y(x) = CSA + CSB*(x-x(k)) + CSC*(x-x(k))**2 + CSD*(x-x(k))**3 447 ! 448 ALLOCATE(CSA(NrgpA)); CSA=0.0_DP 449 ALLOCATE(CSB(NrgpA)); CSB=0.0_DP 450 ALLOCATE(CSC(NrgpA)); CSC=0.0_DP 451 ALLOCATE(CSD(NrgpA)); CSD=0.0_DP 452 ! 453 DO ir=1,NrgpA-1 454 ! 455 ! CSA(k) := y(k) 456 ! 457 CSA(ir)=atrhor(ir) 458 ! 459 ! CSB(k) := delta(y)/delta(x) - 1/3*delta(x)*y''(k) - 1/6*delta(x)*y''(k+1) 460 ! 461 CSB(ir)=(atrhor(ir+1)-atrhor(ir))/(atgrdr(ir+1)-atgrdr(ir)) 462 CSB(ir)=CSB(ir)-((1.0_DP/3.0_DP)*(atgrdr(ir+1)-atgrdr(ir))*d2atrhor(ir)) 463 CSB(ir)=CSB(ir)-((1.0_DP/6.0_DP)*(atgrdr(ir+1)-atgrdr(ir))*d2atrhor(ir+1)) 464 ! 465 ! CSC(k) := 1/2*y''(k) 466 ! 467 CSC(ir)=(1.0_DP/2.0_DP)*d2atrhor(ir) 468 ! 469 ! CSD(k) := 1/6*delta(y'')/delta(x) 470 ! 471 CSD(ir)=((1.0_DP/6.0_DP)*(d2atrhor(ir+1)-d2atrhor(ir))/(atgrdr(ir+1)-atgrdr(ir))) 472 ! 473 END DO 474 ! 475 ! Pack species-specific radial cutoff into (is,1) of spdata array... 476 ! 477 spdata(is,1)=atgrdr(icutrA) 478 IF (spdata(is,1).GT.spcutAmax) spcutAmax=spdata(is,1) 479 ! 480 ! Compute and pack grid spacing of species-specific linear equispaced grid into (is,2) of spdata array... 481 ! 482 spdata(is,2)=(atgrdr(icutrA)+1.0_DP)/DBLE(NgpA) !include additional buffer of 1 bohr... 483 WRITE(stdout,'(5X,"Linear grid spacing was computed as: ",F25.15," bohr.")') spdata(is,2) 484 ! 485 ! Form linear equispaced atomic grid (NOT including point at r=0) and pack into argument (is,:) of spgrd array... 486 ! 487 DO ip=1,NgpA 488 ! 489 spgrd(is,ip)=DBLE(ip)*spdata(is,2) 490 ! 491 END DO 492 ! 493 ! Map atomic pseudo-density (currently on the radial atomic grid) onto linear equispaced atomic grid using 494 ! cubic spline interpolation...Form first derivative of the atomic pseudo-density on the linear equispaced 495 ! atomic grid via differentiation of the cubic spline interpolant... 496 ! 497 DO ip=1,NgpA 498 ! 499 d=spgrd(is,ip) 500 ! 501 IF (uniform_grid.EQV..TRUE.) THEN 502 ! 503 ir=INT(d/dxA)+1 !since the equispaced/uniform grid first point is at r=0... 504 ! 505 ELSE 506 ! 507 ir=FLOOR(DLOG(d*EXP(dxA)/gfctrA)/dxA) 508 ! 509 END IF 510 ! 511 dk1=d-atgrdr(ir); dk2=dk1*dk1; dk3=dk2*dk1 512 sprho(is,ip)=CSA(ir)+CSB(ir)*dk1+CSC(ir)*dk2+CSD(ir)*dk3 !Pack density into argument (is,:) of sprho array 513 spdrho(is,ip)=CSB(ir)+2.0_DP*CSC(ir)*dk1+3.0_DP*CSD(ir)*dk2 !Pack density derivative into argument (is,:) of spdrho array 514 ! 515 END DO 516 ! 517 ! For computational efficiency during the remainder of the calculation, extrapolate sprho and spdrho to 518 ! include the point at r=0 (this eliminates an if statement in crucial inner loops)... 519 ! Use quadratic extrapolation to obtain these points...Hence the 0:NgpA dimension above... 520 ! 521 spgrd(is,0)=0.0_DP !Extend linear grid to include point at r=0... 522 ! 523 L1=((0.0_DP-spgrd(is,2))*(0.0_DP-spgrd(is,3)))/((spgrd(is,1)-spgrd(is,2))*(spgrd(is,1)-spgrd(is,3))) 524 L2=((0.0_DP-spgrd(is,1))*(0.0_DP-spgrd(is,3)))/((spgrd(is,2)-spgrd(is,1))*(spgrd(is,2)-spgrd(is,3))) 525 L3=((0.0_DP-spgrd(is,1))*(0.0_DP-spgrd(is,2)))/((spgrd(is,3)-spgrd(is,1))*(spgrd(is,3)-spgrd(is,2))) 526 sprho(is,0)=L1*sprho(is,1)+L2*sprho(is,2)+L3*sprho(is,3) !Extend atomic pseudo-density to include point at r=0... 527 spdrho(is,0)=L1*spdrho(is,1)+L2*spdrho(is,2)+L3*spdrho(is,3) !Extend atomic pseudo-density derivative to include point at r=0... 528 ! 529 ! Throughout the remainder of the code, to map the atomic quantities onto the real-space mesh, we will be 530 ! utilizing the Taylor series form of linear interpolation, given by: 531 ! 532 ! y(x) = LIA + LIB*(x-x(k)) y'(x) = dLIA + dLIB*(x-x(k)) 533 ! 534 ! for x(k) <= x <= x(k+1)... 535 ! 536 DO ip=0,NgpA-1 537 ! 538 ! LIA(k) := y(k) 539 ! 540 LIA(is,ip)=sprho(is,ip) 541 dLIA(is,ip)=spdrho(is,ip) 542 ! 543 ! LIB(k) := delta(y)/delta(x) 544 ! 545 LIB(is,ip)=(sprho(is,ip+1)-sprho(is,ip))/(spgrd(is,ip+1)-spgrd(is,ip)) 546 dLIB(is,ip)=(spdrho(is,ip+1)-spdrho(is,ip))/(spgrd(is,ip+1)-spgrd(is,ip)) 547 ! 548 END DO 549 ! 550 ! Populate reference free atom quantities... 551 ! 552 CALL GetVdWParam(upf(is)%psd,C6AAfree(is),dpfree(is),R0free(is)) 553 ! 554 WRITE(stdout,'(5X,"The free atom static dipole polarizability is ",F13.6," bohr^3.")') dpfree(is) 555 WRITE(stdout,'(5X,"The free atom homonuclear C6 coefficient is ",F13.6," Hartree bohr^6.")') C6AAfree(is) 556 WRITE(stdout,'(5X,"The free atom vdW radius is ",F13.6," bohr.")') R0free(is) 557 ! 558 ! Clean-up all species-specific temporary arrays 559 ! 560 IF (ALLOCATED(atgrdr)) DEALLOCATE(atgrdr) 561 IF (ALLOCATED(atgrdrab)) DEALLOCATE(atgrdrab) 562 IF (ALLOCATED(atrhor)) DEALLOCATE(atrhor) 563 IF (ALLOCATED(datrhor)) DEALLOCATE(datrhor) 564 IF (ALLOCATED(d2atrhor)) DEALLOCATE(d2atrhor) 565 IF (ALLOCATED(CSA)) DEALLOCATE(CSA) 566 IF (ALLOCATED(CSB)) DEALLOCATE(CSB) 567 IF (ALLOCATED(CSC)) DEALLOCATE(CSC) 568 IF (ALLOCATED(CSD)) DEALLOCATE(CSD) 569 ! 570 END DO !is 571 ! 572 ! Compute free heteronuclear C6 coefficient matrix... 573 ! C6ABfree(A,B)=[2*C6AAfree(A)*C6AAfree(B)]/[(dpfree(B)/dpfree(A))*C6AAfree(A)+(dpfree(A)/dpfree(B))*C6AAfree(B)] 574 ! 575 DO is=1,nsp 576 ! 577 DO it=1,nsp 578 ! 579 num=2.0_DP*C6AAfree(is)*C6AAfree(it) 580 den=(dpfree(it)/dpfree(is))*C6AAfree(is)+(dpfree(is)/dpfree(it))*C6AAfree(it) 581 C6ABfree(is,it)=num/den 582 ! 583 END DO 584 ! 585 END DO 586 ! 587 RETURN 588 ! 589 !-------------------------------------------------------------------------------------------------------------- 590 END SUBROUTINE tsvdw_initialize 591 !-------------------------------------------------------------------------------------------------------------- 592 ! 593 !-------------------------------------------------------------------------------------------------------------- 594 SUBROUTINE tsvdw_calculate(tauin, rhor) 595 !-------------------------------------------------------------------------------------------------------------- 596 ! TS-vdW Management Code: Manages entire calculation of TS-vdW energy, wavefunction forces, and ion forces via 597 ! calls to PRIVATE subroutines below (called in each MD step). The calls to tsvdw_initialize and tsvdw_finalize 598 ! are done once at the beginning (init_run) and the end (terminate_run). 599 !-------------------------------------------------------------------------------------------------------------- 600 ! 601 IMPLICIT NONE 602 ! 603 ! I/O variables 604 ! 605 REAL(DP), INTENT(IN) :: rhor(:) 606 REAL(DP) :: tauin(3,nat) 607 ! 608 ! Parallel initialization... 609 ! 610 CALL tsvdw_para_init() 611 ! 612 ! Move all atoms into simulation cell by adjusting Cartesian coordinates according to PBCs... 613 ! 614 CALL tsvdw_pbc(tauin) 615 ! 616 ! Compute unique atom pair list... 617 ! 618 CALL tsvdw_unique_pair() 619 ! 620 ! Obtain molecular charge density given on the real-space mesh... 621 ! 622 CALL tsvdw_rhotot( rhor ) 623 ! 624 ! Determine spherical atomic integration domains and atom overlap (bit array)... 625 ! Compute molecular pro-density (superposition of atomic densities) on the real-space mesh... 626 ! Compute functional derivative of vdW energy wrt charge density (numerator only)... 627 ! 628 CALL tsvdw_screen() 629 ! 630 ! Compute effective volume for each atom in the simulation cell... 631 ! Complete functional derivative of vdW energy wrt charge density... 632 ! 633 CALL tsvdw_veff() 634 ! 635 ! Calculate first derivative of veff wrt nuclear and cell displacements... 636 ! 637 CALL tsvdw_dveff() 638 ! 639 ! Calculate effective quantities for each atom in the simulation cell... 640 ! 641 CALL tsvdw_effqnts() 642 ! 643 ! Calculate total TS-vdW energy, dispersion potential prefactor, ionic forces, and cell forces... 644 ! 645 CALL tsvdw_energy() 646 ! 647 ! Calculate total TS-vdW wavefunction forces (dispersion potential)... 648 ! 649 CALL tsvdw_wfforce() 650 ! 651 ! Deallocate all arrays specific to tsvdw_calculate... 652 ! 653 CALL tsvdw_cleanup() 654 ! 655 RETURN 656 ! 657 !-------------------------------------------------------------------------------------------------------------- 658 END SUBROUTINE tsvdw_calculate 659 !-------------------------------------------------------------------------------------------------------------- 660 ! 661 !-------------------------------------------------------------------------------------------------------------- 662 SUBROUTINE tsvdw_para_init() 663 !-------------------------------------------------------------------------------------------------------------- 664 ! 665 IMPLICIT NONE 666 ! 667 INTEGER :: i,j,k 668 ! 669 me=me_image+1 670 ! 671 ALLOCATE(nstates(nproc_image)); nstates=0 672 ALLOCATE(sdispls(nproc_image)); sdispls=0 673 ALLOCATE(sendcount(nproc_image)); sendcount=0 674 ALLOCATE(rdispls(nproc_image)); rdispls=0 675 ALLOCATE(recvcount(nproc_image)); recvcount=0 676 ALLOCATE(istatus(nproc_image)); istatus=0 677 ! 678 ! Assign workload of atoms over nproc_image processors 679 ! 680 IF (nat.LE.nproc_image) THEN 681 ! 682 DO i=1,nat 683 ! 684 nstates(i)=1 685 ! 686 END DO 687 ! 688 ELSE 689 ! 690 k=0 691 ! 69210 DO j=1,nproc_image 693 ! 694 nstates(j)=nstates(j)+1 695 ! 696 k=k+1 697 ! 698 IF (k.GE.nat) GO TO 20 699 ! 700 END DO 701 ! 702 IF (k.LT.nat) GO TO 10 703 ! 704 END IF 705 ! 706 20 CONTINUE 707 ! 708 RETURN 709 ! 710 !-------------------------------------------------------------------------------------------------------------- 711 END SUBROUTINE tsvdw_para_init 712 !-------------------------------------------------------------------------------------------------------------- 713 ! 714 !-------------------------------------------------------------------------------------------------------------- 715 SUBROUTINE tsvdw_pbc(tauin) 716 !-------------------------------------------------------------------------------------------------------------- 717 ! 718 IMPLICIT NONE 719 ! 720 ! I/O variables 721 ! 722 REAL(DP) :: tauin(3,nat) 723 ! 724 ! Local variables 725 ! 726 INTEGER :: ia 727 REAL(DP), DIMENSION(:,:), ALLOCATABLE :: atxyzs 728 ! 729 ! Initialization of PBC-adjusted Cartesian coordinates... 730 ! 731 ALLOCATE(atxyz(3,nat)); atxyz=0.0_DP 732 ALLOCATE(atxyzs(3,nat)); atxyzs=0.0_DP 733 ! 734 ! Adjust Cartesian coordinates of ions according to periodic boundary conditions... 735 ! N.B.: PBC are imposed here in the range [0,1)... 736 ! 737 DO ia = 1, nat 738 ! 739 atxyzs(1,ia)=ainv(1,1)*tauin(1,ia)+ainv(1,2)*tauin(2,ia)+ainv(1,3)*tauin(3,ia) ! s = h^-1 r 740 atxyzs(2,ia)=ainv(2,1)*tauin(1,ia)+ainv(2,2)*tauin(2,ia)+ainv(2,3)*tauin(3,ia) ! s = h^-1 r 741 atxyzs(3,ia)=ainv(3,1)*tauin(1,ia)+ainv(3,2)*tauin(2,ia)+ainv(3,3)*tauin(3,ia) ! s = h^-1 r 742 ! 743 atxyzs(1,ia)=atxyzs(1,ia)-FLOOR(atxyzs(1,ia)) ! impose PBC on s in range: [0,1) 744 atxyzs(2,ia)=atxyzs(2,ia)-FLOOR(atxyzs(2,ia)) ! impose PBC on s in range: [0,1) 745 atxyzs(3,ia)=atxyzs(3,ia)-FLOOR(atxyzs(3,ia)) ! impose PBC on s in range: [0,1) 746 ! 747 atxyz(1,ia)=h(1,1)*atxyzs(1,ia)+h(1,2)*atxyzs(2,ia)+h(1,3)*atxyzs(3,ia) ! r = h s 748 atxyz(2,ia)=h(2,1)*atxyzs(1,ia)+h(2,2)*atxyzs(2,ia)+h(2,3)*atxyzs(3,ia) ! r = h s 749 atxyz(3,ia)=h(3,1)*atxyzs(1,ia)+h(3,2)*atxyzs(2,ia)+h(3,3)*atxyzs(3,ia) ! r = h s 750 ! 751 END DO 752 ! 753 IF (ALLOCATED(atxyzs)) DEALLOCATE(atxyzs) 754 ! 755 RETURN 756 ! 757 !-------------------------------------------------------------------------------------------------------------- 758 END SUBROUTINE tsvdw_pbc 759 !-------------------------------------------------------------------------------------------------------------- 760 ! 761 !-------------------------------------------------------------------------------------------------------------- 762 SUBROUTINE tsvdw_unique_pair() 763 !-------------------------------------------------------------------------------------------------------------- 764 ! 765 IMPLICIT NONE 766 ! 767 ! Local variables 768 ! 769 INTEGER :: ia,ib,ias,ibs,ip,ir,i,j,k,jj,nj_max,nbmax,num,num1,jj_neib_of_i 770 REAL(DP) :: spcutA,spcutB,dAB(3),dAB2(3) 771 INTEGER, DIMENSION(:), ALLOCATABLE :: nj,overlap2 772 INTEGER, DIMENSION(:,:), ALLOCATABLE :: overlap 773 REAL(DP), DIMENSION(:), ALLOCATABLE :: dABmic 774 ! 775 CALL start_clock('tsvdw_pair') 776 ! 777 ! Allocate and initialize temporary arrays... 778 ! 779 ALLOCATE(dABmic(nat)); dABmic=0.0_DP 780 ALLOCATE(overlap(nat,nat)); overlap=0 781 ALLOCATE(overlap2(nat)); overlap2=0 782 ALLOCATE(nj(nat)); nj=0 783 ! 784 ! Outer loop over atoms A to form non-unique atom pair overlap matrix... 785 ! 786 DO ia=1,nat 787 ! 788 nj(ia)=0; dABmic=0.0_DP 789 ! 790 ! Connect atom type with species-dependent quantities... 791 ! 792 ias=ityp(ia) 793 ! 794 ! Transfer species-specific cutoff to spcutA... 795 ! 796 spcutA=spdata(ias,1) 797 ! 798 ! Inner loop over atoms B... 799 ! 800 DO ib=1,nat 801 ! 802 IF(ib.NE.ia) THEN 803 ! 804 ! Connect atom type with species-dependent quantities... 805 ! 806 ibs=ityp(ib) 807 ! 808 ! Transfer species-specific cutoff to spcutB... 809 ! 810 spcutB=spdata(ibs,1) 811 ! 812 ! Compute distance between atom A and atom B (according to the minimum image convention)... 813 ! 814 dAB(1)=atxyz(1,ia)-atxyz(1,ib) ! r_AB = r_A - r_B 815 dAB(2)=atxyz(2,ia)-atxyz(2,ib) ! r_AB = r_A - r_B 816 dAB(3)=atxyz(3,ia)-atxyz(3,ib) ! r_AB = r_A - r_B 817 ! 818 dAB2(1)=ainv(1,1)*dAB(1)+ainv(1,2)*dAB(2)+ainv(1,3)*dAB(3) ! s_AB = h^-1 r_AB 819 dAB2(2)=ainv(2,1)*dAB(1)+ainv(2,2)*dAB(2)+ainv(2,3)*dAB(3) ! s_AB = h^-1 r_AB 820 dAB2(3)=ainv(3,1)*dAB(1)+ainv(3,2)*dAB(2)+ainv(3,3)*dAB(3) ! s_AB = h^-1 r_AB 821 ! 822 dAB2(1)=dAB2(1)-IDNINT(dAB2(1)) ! impose MIC on s_AB in range: [-0.5,+0.5] 823 dAB2(2)=dAB2(2)-IDNINT(dAB2(2)) ! impose MIC on s_AB in range: [-0.5,+0.5] 824 dAB2(3)=dAB2(3)-IDNINT(dAB2(3)) ! impose MIC on s_AB in range: [-0.5,+0.5] 825 ! 826 dAB(1)=h(1,1)*dAB2(1)+h(1,2)*dAB2(2)+h(1,3)*dAB2(3) ! r_AB = h s_AB (MIC) 827 dAB(2)=h(2,1)*dAB2(1)+h(2,2)*dAB2(2)+h(2,3)*dAB2(3) ! r_AB = h s_AB (MIC) 828 dAB(3)=h(3,1)*dAB2(1)+h(3,2)*dAB2(2)+h(3,3)*dAB2(3) ! r_AB = h s_AB (MIC) 829 ! 830 dABmic(ib)=DSQRT(dAB(1)*dAB(1)+dAB(2)*dAB(2)+dAB(3)*dAB(3)) ! |r_A - r_B| (MIC) 831 ! 832 IF(dABmic(ib).LT.(spcutA+spcutB)) THEN 833 ! 834 nj(ia)=nj(ia)+1 835 overlap(nj(ia),ia)=ib 836 ! 837 IF(nj(ia).EQ.1) THEN 838 ! 839 overlap(nj(ia),ia)=ib 840 ! 841 ELSE IF(dABmic(overlap(nj(ia)-1,ia)).LE.dABmic(ib)) THEN 842 ! 843 overlap(nj(ia),ia)=ib 844 ! 845 ELSE 846 ! 847 overlap2(:)=0 848 ! 849 DO ir=1,nj(ia)-1 850 ! 851 IF(dABmic(overlap(ir,ia)).LT.dABmic(ib)) THEN 852 ! 853 overlap2(ir)=overlap(ir,ia) 854 ! 855 ELSE 856 ! 857 overlap2(ir)=ib 858 ! 859 DO ip=ir+1,nj(ia) 860 ! 861 overlap2(ip)=overlap(ip-1,ia) 862 ! 863 END DO 864 ! 865 GO TO 30 866 ! 867 END IF 868 ! 869 END DO !ir 870 ! 87130 CONTINUE 872 ! 873 DO ir=1,nj(ia) 874 ! 875 overlap(ir,ia)=overlap2(ir) 876 ! 877 END DO 878 ! 879 END IF !nj(ia) 880 ! 881 END IF !dABmic(j) 882 ! 883 END IF !ia/=ib 884 ! 885 END DO !ib 886 ! 887 END DO !ia 888 ! 889 IF (ALLOCATED(dABmic)) DEALLOCATE(dABmic) 890 IF (ALLOCATED(overlap2)) DEALLOCATE(overlap2) 891 ! 892 ! Now form unique atom pair overlap matrix... 893 ! 894 nbmax=nat 895 ! 896 ALLOCATE(pair(nbmax,nat)); pair=0 897 ALLOCATE(npair(nat)); npair=0 898 ! 899 num=0; num1=0 900 ! 901 DO j=1,nbmax 902 ! 903 DO ia=1,nat 904 ! 905 DO jj=1,nj(ia) 906 ! 907 jj_neib_of_i=overlap(jj,ia) 908 ! 909 IF(jj_neib_of_i.GT.0) THEN 910 ! 911 pair(j,ia)=jj_neib_of_i 912 overlap(jj,ia)=0 913 num=num+1 914 ! 915 DO k=1,nj(jj_neib_of_i) 916 ! 917 IF(overlap(k,jj_neib_of_i).EQ.ia) THEN 918 ! 919 overlap(k,jj_neib_of_i)=0 920 num1=num1+1 921 ! 922 GO TO 40 923 ! 924 END IF 925 ! 926 END DO !k 927 ! 928 END IF 929 ! 930 END DO !jj 931 ! 93240 CONTINUE 933 ! 934 END DO !ia 935 ! 936 END DO !j 937 ! 938 IF(num.NE.num1) THEN 939 ! 940 CALL errore('tsvdw','ERROR: num .NE. num1...',1) 941 ! 942 END IF 943 ! 944 ! Count number of unique atom pairs for each atom... 945 ! 946 DO ia=1,nat 947 ! 948 num=0 949 ! 950 DO j=1,nbmax 951 ! 952 IF(pair(j,ia).NE.0) num=num+1 953 ! 954 END DO 955 ! 956 npair(ia)=num 957 ! 958 END DO 959 ! 960 IF (ALLOCATED(overlap)) DEALLOCATE(overlap) 961 IF (ALLOCATED(nj)) DEALLOCATE(nj) 962 ! 963 CALL stop_clock('tsvdw_pair') 964 ! 965 RETURN 966 ! 967 !-------------------------------------------------------------------------------------------------------------- 968 END SUBROUTINE tsvdw_unique_pair 969 !-------------------------------------------------------------------------------------------------------------- 970 ! 971 !-------------------------------------------------------------------------------------------------------------- 972 SUBROUTINE tsvdw_rhotot( rhor ) 973 !-------------------------------------------------------------------------------------------------------------- 974 ! 975 IMPLICIT NONE 976 ! 977 REAL(DP), INTENT(IN) :: rhor(:) 978 ! 979 ! Local variables 980 ! 981 INTEGER :: ir,ierr 982 REAL(DP), DIMENSION(:), ALLOCATABLE :: rhor_tmp 983 ! 984 CALL start_clock('tsvdw_rhotot') 985 ! 986 ! Initialization of rhotot array (local copy of the real-space charge density)... 987 ! 988 ALLOCATE(rhotot(nr1*nr2*nr3)); rhotot=0.0_DP 989#if defined(__MPI) 990 ! 991 ! Initialization of rhor_tmp temporary buffers... 992 ! 993 ALLOCATE(rhor_tmp(nr1*nr2*nr3)); rhor_tmp=0.0_DP 994 ! 995 ! Collect distributed rhor and broadcast to all processors... 996 ! 997 DO iproc=1,nproc_bgrp 998 ! 999 recvcount(iproc)=dfftp%nr3p(iproc)*nr1*nr2 1000 ! 1001 END DO 1002 ! 1003 rdispls(1) = 0 1004 ! 1005 DO iproc=2,nproc_bgrp 1006 ! 1007 rdispls(iproc)=rdispls(iproc-1)+recvcount(iproc-1) 1008 ! 1009 END DO 1010 ! 1011 CALL MPI_ALLGATHERV(rhor(1),dfftp%nr3p(me_bgrp+1)*nr1*nr2,& 1012 MPI_DOUBLE_PRECISION,rhor_tmp(1),recvcount,rdispls,& 1013 MPI_DOUBLE_PRECISION,intra_bgrp_comm,ierr) 1014 ! 1015 ! Transfer rhor temporary arrays to rhotot array... 1016 ! 1017 rhotot=rhor_tmp 1018 ! 1019 ! Clean-up temporary arrays... 1020 ! 1021 IF (ALLOCATED(rhor_tmp)) DEALLOCATE(rhor_tmp) 1022 ! 1023#else 1024 rhotot(:) = rhor(:) 1025#endif 1026 CALL stop_clock('tsvdw_rhotot') 1027 ! 1028 RETURN 1029 ! 1030 !-------------------------------------------------------------------------------------------------------------- 1031 END SUBROUTINE tsvdw_rhotot 1032 !-------------------------------------------------------------------------------------------------------------- 1033 ! 1034 !-------------------------------------------------------------------------------------------------------------- 1035 SUBROUTINE tsvdw_screen() 1036 !-------------------------------------------------------------------------------------------------------------- 1037 ! 1038 IMPLICIT NONE 1039 ! 1040 ! Local variables 1041 ! 1042 INTEGER :: ia,ias,Ntmp,Ntmpr,Npts,Nptsr,ir1,ir2,ir3,Ndim,Nints,off1,off1r,ioff,boff,iq,ir 1043 REAL(DP) :: spcutA,spdA,dq(3),dqA(3),dqAs(3),dk1,rhoA 1044 REAL(DP), ALLOCATABLE :: dqAmic(:,:,:),dveffAdntmp(:,:,:) 1045 ! 1046 CALL start_clock('tsvdw_screen') 1047 ! 1048 ! Allocate and initialize gomegar array which contains (in a bit array) which atoms contribute to a given point 1049 ! on the reduced real-space grid for all atoms... 1050 ! 1051 Nints=nat/bsint+1 1052 ALLOCATE(gomegar(nr1r*nr2r*nr3r,Nints)); gomegar=0 1053 ! 1054 ! Allocate and initialize NsomegaA and NsomegaAr arrays which contains the number of points in the 1055 ! full and reduced spherical atomic integration domains for all atoms... 1056 ! 1057 ALLOCATE(NsomegaA(nat)); NsomegaA=0 1058 ALLOCATE(NsomegaAr(nat)); NsomegaAr=0 1059 ! 1060 ! Allocate and initialize somegaA and somegaAr arrays which contains the grid indices (along nr1,nr2,nr3) for each point in the 1061 ! full and reduced spherical atomic integration domains for each atom assigned to a given processor... 1062 ! 1063 Ntmp=INT(1.10_DP*((4.0_DP/3.0_DP)*pi*spcutAmax**(3.0_DP)/omega)*(nr1*nr2*nr3)) ! Number of points in the full sphere + 10% buffer (to be safe) 1064 Ntmpr=INT(1.10_DP*((4.0_DP/3.0_DP)*pi*spcutAmax**(3.0_DP)/omega)*(nr1r*nr2r*nr3r)) ! Number of points in the reduced sphere + 10% buffer (to be safe) 1065 Ndim=MAX(1,nstates(me)) 1066 ALLOCATE(somegaA(Ntmp,3,Ndim)); somegaA=0 1067 ALLOCATE(somegaAr(Ntmpr,3,Ndim)); somegaAr=0 1068 ! 1069 ! Allocate and initialize gomegaAr array which contains in a bit array all of the atoms that intersect with each point in the 1070 ! reduced spherical atomic integration domain for each atom assigned to a given processor... 1071 ! 1072 ALLOCATE(gomegaAr(Ntmpr,Nints,Ndim)); gomegaAr=0 1073 ! 1074 ! Initialization of rhosad(r)... 1075 ! 1076 ALLOCATE(rhosad(nr1*nr2*nr3)); rhosad=0.0_DP 1077 ! 1078 ! Initialization of dVA/dn(r)... 1079 ! 1080 ALLOCATE(dveffAdn(Ntmp,Ndim)); dveffAdn=0.0_DP 1081 ! 1082 DO iproc=1,nstates(me) 1083 ! 1084 ! Connect processor number with atom... 1085 ! 1086 ia=me+nproc_image*(iproc-1) 1087 ! 1088 ! Connect atom type with species-dependent quantities... 1089 ! 1090 ias=ityp(ia) 1091 ! 1092 ! Transfer species-specific cutoff to spcutA... 1093 ! 1094 spcutA=spdata(ias,1) 1095 ! 1096 ! Precompute inverse of species-specific linear grid spacing (replaces / with * inside inner loop)... 1097 ! 1098 spdA=1.0_DP/spdata(ias,2) 1099 ! 1100 ! Loop over grid points and determine if they belong to spherical atomic integration domain (if r < RcutA)... 1101 ! 1102 Npts=0; Nptsr=0 1103 ! 1104 ALLOCATE(dqAmic(nr1,nr2,nr3)); dqAmic=0.0_DP 1105 ALLOCATE(dveffAdntmp(nr1,nr2,nr3)); dveffAdntmp=0.0_DP 1106 ! 1107!$omp parallel do private(dq,dqA,dqAs,ir,dk1,rhoA,off1,ioff,boff,off1r) 1108 DO ir1=1,nr1 1109 ! 1110 dq(1)=DBLE(ir1-1)/DBLE(nr1) ! s_i(1) 1111 ! 1112 DO ir2=1,nr2 1113 ! 1114 dq(2)=DBLE(ir2-1)/DBLE(nr2) ! s_i(2) 1115 ! 1116 DO ir3=1,nr3 1117 ! 1118 dq(3)=DBLE(ir3-1)/DBLE(nr3) ! s_i(3) 1119 ! 1120 ! Compute distance between grid point and atom according to minimum image convention (MIC)... 1121 ! 1122 dqA(1)=h(1,1)*dq(1)+h(1,2)*dq(2)+h(1,3)*dq(3) ! r_i = h s_i 1123 dqA(2)=h(2,1)*dq(1)+h(2,2)*dq(2)+h(2,3)*dq(3) ! r_i = h s_i 1124 dqA(3)=h(3,1)*dq(1)+h(3,2)*dq(2)+h(3,3)*dq(3) ! r_i = h s_i 1125 ! 1126 dqA(1)=dqA(1)-atxyz(1,ia) ! r_iA = r_i - r_A 1127 dqA(2)=dqA(2)-atxyz(2,ia) ! r_iA = r_i - r_A 1128 dqA(3)=dqA(3)-atxyz(3,ia) ! r_iA = r_i - r_A 1129 ! 1130 dqAs(1)=ainv(1,1)*dqA(1)+ainv(1,2)*dqA(2)+ainv(1,3)*dqA(3) ! s_iA = h^-1 r_iA 1131 dqAs(2)=ainv(2,1)*dqA(1)+ainv(2,2)*dqA(2)+ainv(2,3)*dqA(3) ! s_iA = h^-1 r_iA 1132 dqAs(3)=ainv(3,1)*dqA(1)+ainv(3,2)*dqA(2)+ainv(3,3)*dqA(3) ! s_iA = h^-1 r_iA 1133 ! 1134 dqAs(1)=dqAs(1)-IDNINT(dqAs(1)) ! impose MIC on s_iA in range: [-0.5,+0.5] 1135 dqAs(2)=dqAs(2)-IDNINT(dqAs(2)) ! impose MIC on s_iA in range: [-0.5,+0.5] 1136 dqAs(3)=dqAs(3)-IDNINT(dqAs(3)) ! impose MIC on s_iA in range: [-0.5,+0.5] 1137 ! 1138 dqA(1)=h(1,1)*dqAs(1)+h(1,2)*dqAs(2)+h(1,3)*dqAs(3) ! r_iA = h s_iA (MIC) 1139 dqA(2)=h(2,1)*dqAs(1)+h(2,2)*dqAs(2)+h(2,3)*dqAs(3) ! r_iA = h s_iA (MIC) 1140 dqA(3)=h(3,1)*dqAs(1)+h(3,2)*dqAs(2)+h(3,3)*dqAs(3) ! r_iA = h s_iA (MIC) 1141 ! 1142 dqAmic(ir1,ir2,ir3)=DSQRT(dqA(1)*dqA(1)+dqA(2)*dqA(2)+dqA(3)*dqA(3)) ! |r_i - r_A| (MIC) 1143 ! 1144 ! Screen grid point according to atomic radial cutoff... 1145 ! 1146 IF (dqAmic(ir1,ir2,ir3).LE.spcutA) THEN 1147 ! 1148 ! Form rhosad(r) on the real-space mesh... 1149 ! N.B. This algorithm only works when the images of a given atom are greater than the radial grid cutoff values for ALL atomic species... 1150 ! 1151 ! Determine the index in the atomic linear equispaced grid such that grd(ir) <= dqA <= grd(ir+1) and distance between dqA and grd(ir)... 1152 ! 1153 ir=INT(dqAmic(ir1,ir2,ir3)*spdA) 1154 dk1=dqAmic(ir1,ir2,ir3)-spgrd(ias,ir) 1155 ! 1156 ! Perform linear interpolation to obtain the value of the atomic pseudo-density at the given grid point... 1157 ! 1158 rhoA=LIA(ias,ir)+LIB(ias,ir)*dk1 1159 ! 1160 ! Increment contribution to rhosad(r)... 1161 ! 1162 off1=ir1+(ir2-1)*nr1+(ir3-1)*nr1*nr2 !global offset [nr1,nr2,nr3] 1163 rhosad(off1)=rhosad(off1)+rhoA 1164 ! 1165 ! Form numerator of dVA/dn(r) only... 1166 ! 1167 dveffAdntmp(ir1,ir2,ir3)=dqAmic(ir1,ir2,ir3)**(3.0_DP)*rhoA 1168 ! 1169 ! On reduced grid only, form screened somegaAr and gomegar... 1170 ! 1171 IF ((MOD(ir1,2).EQ.1).AND.(MOD(ir2,2).EQ.1).AND.(MOD(ir3,2).EQ.1)) THEN 1172 ! 1173 ioff=((ia-1)/bsint)+1 ! integer offset for gomegar bit array 1174 boff=(ia-((ioff-1)*bsint))-1 ! bit offset for gomegar bit array 1175 off1r=(ir1+1)/2+((ir2-1)/2)*nr1r+((ir3-1)/2)*nr1r*nr2r ! reduced global offset [nr1r,nr2r,nr3r] 1176 ! 1177 gomegar(off1r,ioff)=IBSET(gomegar(off1r,ioff),boff) 1178 ! 1179 END IF 1180 ! 1181 END IF 1182 ! 1183 END DO !ir3 1184 ! 1185 END DO !ir2 1186 ! 1187 END DO !ir1 1188!$omp end parallel do 1189 ! 1190 DO ir1=1,nr1 1191 ! 1192 DO ir2=1,nr2 1193 ! 1194 DO ir3=1,nr3 1195 ! 1196 ! Screen grid point according to atomic radial cutoff... 1197 ! 1198 IF (dqAmic(ir1,ir2,ir3).LE.spcutA) THEN 1199 ! 1200 Npts=Npts+1 1201 ! 1202 ! Form screened somegaA... 1203 ! 1204 somegaA(Npts,1,iproc)=ir1 1205 somegaA(Npts,2,iproc)=ir2 1206 somegaA(Npts,3,iproc)=ir3 1207 ! 1208 dveffAdn(Npts,iproc)=dveffAdntmp(ir1,ir2,ir3) 1209 ! 1210 ! On reduced grid only, form screened somegaAr ... 1211 ! 1212 IF ((MOD(ir1,2).EQ.1).AND.(MOD(ir2,2).EQ.1).AND.(MOD(ir3,2).EQ.1)) THEN 1213 ! 1214 Nptsr=Nptsr+1 1215 ! 1216 ! Form reduced screened somegaAr... 1217 ! 1218 somegaAr(Nptsr,1,iproc)=ir1 1219 somegaAr(Nptsr,2,iproc)=ir2 1220 somegaAr(Nptsr,3,iproc)=ir3 1221 ! 1222 END IF 1223 ! 1224 END IF 1225 1226 END DO !ir3 1227 ! 1228 END DO !ir2 1229 ! 1230 END DO !ir1 1231 ! 1232 NsomegaA(ia)=Npts 1233 NsomegaAr(ia)=Nptsr 1234 ! 1235 IF (ALLOCATED(dqAmic)) DEALLOCATE(dqAmic) 1236 IF (ALLOCATED(dveffAdntmp)) DEALLOCATE(dveffAdntmp) 1237 ! 1238 END DO ! iproc 1239 ! 1240 ! Collect NsomegaA, NsomegaAr, gomegar, and rhosad over all processors and broadcast... 1241 ! 1242 CALL mp_sum(NsomegaA,intra_image_comm) 1243 CALL mp_sum(NsomegaAr,intra_image_comm) 1244 CALL mp_sum(gomegar,intra_image_comm) 1245 CALL mp_sum(rhosad,intra_image_comm) 1246 ! 1247 ! Decompose gomegar to gomegaAr to save on memory storage... 1248 ! 1249 DO iproc=1,nstates(me) 1250 ! 1251 ! Connect processor number with atom... 1252 ! 1253 ia=me+nproc_image*(iproc-1) 1254 ! 1255 ! Loop over points in the (pre-screened) reduced spherical atomic integration domain... 1256 ! 1257!$omp parallel do private(off1r,ir) 1258 DO iq=1,NsomegaAr(ia) 1259 ! 1260 DO ir=1,Nints 1261 ! 1262 off1r=(somegaAr(iq,1,iproc)+1)/2+((somegaAr(iq,2,iproc)-1)/2)*nr1r+((somegaAr(iq,3,iproc)-1)/2)*nr1r*nr2r ! reduced global offset [nr1r,nr2r,nr3r] 1263 gomegaAr(iq,ir,iproc)=gomegar(off1r,ir) 1264 ! 1265 END DO 1266 ! 1267 END DO !iq 1268!$omp end parallel do 1269 ! 1270 END DO ! iproc 1271 ! 1272 ! Clean-up temporary arrays... 1273 ! 1274 IF (ALLOCATED(gomegar)) DEALLOCATE(gomegar) 1275 ! 1276 CALL stop_clock('tsvdw_screen') 1277 ! 1278 RETURN 1279 ! 1280 !-------------------------------------------------------------------------------------------------------------- 1281 END SUBROUTINE tsvdw_screen 1282 !-------------------------------------------------------------------------------------------------------------- 1283 ! 1284 !-------------------------------------------------------------------------------------------------------------- 1285 SUBROUTINE tsvdw_veff() 1286 !-------------------------------------------------------------------------------------------------------------- 1287 ! 1288 IMPLICIT NONE 1289 ! 1290 ! Local variables 1291 ! 1292 INTEGER :: ia,iq,off1 1293 REAL(DP) :: normr 1294 ! 1295 CALL start_clock('tsvdw_veff') 1296 ! 1297 ! Initialization of effective volume... 1298 ! 1299 ALLOCATE(veff(nat)); veff=0.0_DP 1300 ! 1301 ! Normalization factor for veff integral... 1302 ! 1303 normr=omega/DBLE(nr1r*nr2r*nr3r) 1304 ! 1305 ! Loop over atoms in the simulation cell... 1306 ! 1307 DO iproc=1,nstates(me) 1308 ! 1309 ! Connect processor number with atom... 1310 ! 1311 ia=me+nproc_image*(iproc-1) 1312 ! 1313 ! Loop over points in the (pre-screened) spherical atomic integration domain... 1314 ! 1315!$omp parallel do private(off1),reduction(+:veff) 1316 DO iq=1,NsomegaA(ia) 1317 ! 1318 ! Compute veff integrand and complete dispersion potential (functional derivative of veff(A) wrt charge density)... 1319 ! 1320 ! veff(A) = INT [|r-rA|^3*rhoA(|r-rA|)*rhotot(r)/rhosad(r)] 1321 ! 1322 ! dveff(A)/dn(r) = |r-rA|^3*rhoA(|r-rA|)/rhosad(r) 1323 ! 1324 off1=somegaA(iq,1,iproc)+(somegaA(iq,2,iproc)-1)*nr1+(somegaA(iq,3,iproc)-1)*nr1*nr2 !global offset [nr1,nr2,nr3] 1325 dveffAdn(iq,iproc)=dveffAdn(iq,iproc)/rhosad(off1) 1326 ! 1327 ! Increment veff... 1328 ! 1329 IF ((MOD(somegaA(iq,1,iproc),2).EQ.1).AND.(MOD(somegaA(iq,2,iproc),2).EQ.1).AND.(MOD(somegaA(iq,3,iproc),2).EQ.1)) THEN 1330 ! 1331 veff(ia)=veff(ia)+(dveffAdn(iq,iproc)*rhotot(off1)) 1332 ! 1333 END IF 1334 ! 1335 END DO !iq 1336!$omp end parallel do 1337 ! 1338 ! Apply final normalization to veff integral... 1339 ! 1340 veff(ia)=normr*veff(ia) 1341 ! 1342 END DO !iproc 1343 ! 1344 ! Collect veff over all processors and broadcast... 1345 ! 1346 CALL mp_sum(veff,intra_image_comm) 1347 ! 1348 VefftsvdW = veff 1349 ! 1350 CALL stop_clock('tsvdw_veff') 1351 ! 1352 RETURN 1353 ! 1354 !-------------------------------------------------------------------------------------------------------------- 1355 END SUBROUTINE tsvdw_veff 1356 !-------------------------------------------------------------------------------------------------------------- 1357 ! 1358 !-------------------------------------------------------------------------------------------------------------- 1359 SUBROUTINE tsvdw_dveff() 1360 !-------------------------------------------------------------------------------------------------------------- 1361 ! 1362 IMPLICIT NONE 1363 ! 1364 ! Local variables 1365 ! 1366 INTEGER :: ia,ib,ias,ibs,iq,ir,i,j,ipair,off1,ioff,boff 1367 REAL(DP) :: spcutA,spcutB,spdA,spdB,dq(3),dqA(3),dqAs(3),dqB(3),dqBs(3),dqAmic,dqBmic,dABmic,normr 1368 REAL(DP) :: dk1,dptmp1,dptmp2,dptmp3,dptmp4,dptmp5,rhoA,rhoB,drhoA,drhoB,dVAdRA(3),dVAdRB(3),dVBdRA(3) 1369 REAL(DP), DIMENSION(:), ALLOCATABLE :: predVAdRB 1370 REAL(DP), DIMENSION(:,:), ALLOCATABLE :: dqxyzr,dqAxyzs,predVBdRA 1371 ! 1372 CALL start_clock('tsvdw_dveff') 1373 ! 1374 ! Initialization of the dveff/dR and dveff/dh arrays... 1375 ! 1376 ALLOCATE(dveffdR(nat,nat,3)); dveffdR=0.0_DP 1377 ALLOCATE(dveffdh(nat,3,3)); dveffdh=0.0_DP 1378 ! 1379 ! Normalization factor for dveff integrals... 1380 ! 1381 normr=omega/DBLE(nr1r*nr2r*nr3r) 1382 ! 1383 ! Loop over atoms A in the simulation cell and compute dveffdR and dveffdh... 1384 ! 1385 DO iproc=1,nstates(me) 1386 ! 1387 ! Connect processor number with atom... 1388 ! 1389 ia=me+nproc_image*(iproc-1) 1390 ! 1391 ! Connect atom type with species-dependent quantities... 1392 ! 1393 ias=ityp(ia) 1394 ! 1395 ! Transfer species-specific cutoff to spcutA... 1396 ! 1397 spcutA=spdata(ias,1) 1398 ! 1399 ! Precompute inverse of species-specific linear grid spacing (replaces / with * during interpolation)... 1400 ! 1401 spdA=1.0_DP/spdata(ias,2) 1402 ! 1403 ! Allocate and initialize atom-specific arrays... 1404 ! 1405 ALLOCATE(dqxyzr(NsomegaAr(ia),3)); dqxyzr=0.0_DP 1406 ALLOCATE(dqAxyzs(NsomegaAr(ia),3)); dqAxyzs=0.0_DP 1407 ALLOCATE(predVAdRB(NsomegaAr(ia))); predVAdRB=0.0_DP 1408 ALLOCATE(predVBdRA(NsomegaAr(ia),3)); predVBdRA=0.0_DP 1409 ! 1410 ! Initial loop over points in the (pre-screened) reduced spherical atomic integration domain for atom A to compute 1411 ! self-derivative (dV(A)/dR(A)), quantities necessary for dV(A)/dR(B) and dV(B)/dR(A), and self-contribution to dV(A)/dh... 1412 ! 1413!$omp parallel do private(dq,dqA,dqAs,dqAmic,ir,dk1,rhoA,drhoA, & 1414!$omp off1,dptmp1,dptmp2,dptmp3,dptmp4,dVAdRA,dptmp5,i,j), & 1415!$omp reduction(-:dveffdh),reduction(+:dveffdR) 1416 DO iq=1,NsomegaAr(ia) 1417 ! 1418 ! Compute global/cell reference frame Cartesian coordinates of given real-space grid point... 1419 ! 1420 dq(1)=DBLE(somegaAr(iq,1,iproc)-1)/DBLE(nr1) ! s_i(1) 1421 dq(2)=DBLE(somegaAr(iq,2,iproc)-1)/DBLE(nr2) ! s_i(2) 1422 dq(3)=DBLE(somegaAr(iq,3,iproc)-1)/DBLE(nr3) ! s_i(3) 1423 ! 1424 dqA(1)=h(1,1)*dq(1)+h(1,2)*dq(2)+h(1,3)*dq(3) ! r_i = h s_i 1425 dqA(2)=h(2,1)*dq(1)+h(2,2)*dq(2)+h(2,3)*dq(3) ! r_i = h s_i 1426 dqA(3)=h(3,1)*dq(1)+h(3,2)*dq(2)+h(3,3)*dq(3) ! r_i = h s_i 1427 ! 1428 ! Accumulate the Cartesian coordinates of the real-space grid point in the dqxyzr array: 1429 ! 1430 ! dqxyzr(:,1) := x-coordinate of grid point (global/cell reference frame) 1431 ! dqxyzr(:,2) := y-coordinate of grid point (global/cell reference frame) 1432 ! dqxyzr(:,3) := z-coordinate of grid point (global/cell reference frame) 1433 ! 1434 dqxyzr(iq,1)=dqA(1) 1435 dqxyzr(iq,2)=dqA(2) 1436 dqxyzr(iq,3)=dqA(3) 1437 ! 1438 ! Compute distance between grid point and atom in scaled coordinates (s_iA) according to minimum image convention (MIC)... 1439 ! 1440 dqA(1)=dqA(1)-atxyz(1,ia) ! r_iA = r_i - r_A 1441 dqA(2)=dqA(2)-atxyz(2,ia) ! r_iA = r_i - r_A 1442 dqA(3)=dqA(3)-atxyz(3,ia) ! r_iA = r_i - r_A 1443 ! 1444 dqAs(1)=ainv(1,1)*dqA(1)+ainv(1,2)*dqA(2)+ainv(1,3)*dqA(3) ! s_iA = h^-1 r_iA 1445 dqAs(2)=ainv(2,1)*dqA(1)+ainv(2,2)*dqA(2)+ainv(2,3)*dqA(3) ! s_iA = h^-1 r_iA 1446 dqAs(3)=ainv(3,1)*dqA(1)+ainv(3,2)*dqA(2)+ainv(3,3)*dqA(3) ! s_iA = h^-1 r_iA 1447 ! 1448 dqAs(1)=dqAs(1)-IDNINT(dqAs(1)) ! impose MIC on s_iA in range: [-0.5,+0.5] 1449 dqAs(2)=dqAs(2)-IDNINT(dqAs(2)) ! impose MIC on s_iA in range: [-0.5,+0.5] 1450 dqAs(3)=dqAs(3)-IDNINT(dqAs(3)) ! impose MIC on s_iA in range: [-0.5,+0.5] 1451 ! 1452 ! Accumulate the components of the s_i - s_A vector in the dqAxyzs array: 1453 ! 1454 ! dqAxyzs(:,1) := 1-coordinate of s_i - s_A vector (local/atom reference frame) 1455 ! dqAxyzs(:,2) := 2-coordinate of s_i - s_A vector (local/atom reference frame) 1456 ! dqAxyzs(:,3) := 3-coordinate of s_i - s_A vector (local/atom reference frame) 1457 ! 1458 dqAxyzs(iq,1)=dqAs(1) 1459 dqAxyzs(iq,2)=dqAs(2) 1460 dqAxyzs(iq,3)=dqAs(3) 1461 ! 1462 ! Convert MIC distance components from scaled coordinates to Cartesian coordinates (s_iA -> r_iA)... 1463 ! 1464 dqA(1)=h(1,1)*dqAs(1)+h(1,2)*dqAs(2)+h(1,3)*dqAs(3) ! r_iA = h s_iA (MIC) 1465 dqA(2)=h(2,1)*dqAs(1)+h(2,2)*dqAs(2)+h(2,3)*dqAs(3) ! r_iA = h s_iA (MIC) 1466 dqA(3)=h(3,1)*dqAs(1)+h(3,2)*dqAs(2)+h(3,3)*dqAs(3) ! r_iA = h s_iA (MIC) 1467 ! 1468 dqAmic=DSQRT(dqA(1)*dqA(1)+dqA(2)*dqA(2)+dqA(3)*dqA(3)) ! |r_i - r_A| (MIC) 1469 ! 1470 ! Determine the index in the atomic linear equispaced grid such that grd(ir) <= dqA <= grd(ir+1) and distance between dqA and grd(ir)... 1471 ! 1472 ir=INT(dqAmic*spdA) 1473 dk1=dqAmic-spgrd(ias,ir) 1474 ! 1475 ! Perform linear interpolation to obtain the value of the atomic pseudo-density and its derivative at the given grid point... 1476 ! 1477 rhoA=LIA(ias,ir)+LIB(ias,ir)*dk1 !rhoA at grid point via linear interpolation 1478 drhoA=dLIA(ias,ir)+dLIB(ias,ir)*dk1 !drhoA at grid point via linear interpolation 1479 ! 1480 ! Compute global offset for rhosad(r) and rhotot(r), both computed on the real-space mesh... 1481 ! 1482 off1=somegaAr(iq,1,iproc)+(somegaAr(iq,2,iproc)-1)*nr1+(somegaAr(iq,3,iproc)-1)*nr1*nr2 !global offset [nr1,nr2,nr3] 1483 ! 1484 ! Compute self-derivative dVA/dpA integrand for p={x,y,z}... 1485 ! 1486 ! dVA/dpA = INT {(p-pA)*|r-rA|*rhotot(r)/rhosad(r)*[|r-rA|*rho(|r-rA|)*drho(|r-rA|)/rhosad(r)-|r-rA|*drho(|r-rA|)-3*rho(|r-rA|)]} 1487 ! 1488 dptmp1=1.0_DP/rhosad(off1) 1489 dptmp2=dqAmic*drhoA 1490 dptmp3=dqAmic*rhotot(off1)*dptmp1 1491 dptmp4=((rhoA*dptmp1-1.0_DP)*dptmp2-3.0_DP*rhoA)*dptmp3 1492 ! 1493 dVAdRA=dqA*dptmp4 !dVA/dpA integrand/contribution for the given grid point... 1494 ! 1495 ! Increment self-derivative dVA/dpA for p={x,y,z}... 1496 ! 1497 DO i=1,3 1498 ! 1499 dveffdR(ia,ia,i)=dveffdR(ia,ia,i)+dVAdRA(i) 1500 ! 1501 END DO !i 1502 ! 1503 ! Increment self-contribution to dVA/dhpq for p,q={x,y,z}... 1504 ! 1505 ! dVA/dhpq <-- INT {-(p-pA)*(qs-qsA)*|r-rA|*rhotot(r)/rhosad(r)*[|r-rA|*rho(|r-rA|)*drho(|r-rA|)/rhosad(r)-|r-rA|*drho(|r-rA|)-3*rho(|r-rA|)]} 1506 ! 1507 DO i=1,3 1508 ! 1509 DO j=1,3 1510 ! 1511 dveffdh(ia,i,j)=dveffdh(ia,i,j)-dVAdRA(i)*dqAxyzs(iq,j) 1512 ! 1513 END DO !j 1514 ! 1515 END DO !i 1516 ! 1517 ! Precompute quantities necessary for dV(A)/dR(B) and dV(B)/dR(A)... 1518 ! 1519 predVAdRB(iq)=dptmp1*rhoA*dqAmic*dqAmic*dptmp3 1520 ! 1521 dptmp5=dptmp1*dptmp1*drhoA*rhotot(off1) 1522 ! 1523 IF (dqAmic.LT.(1.0E-12_DP)) THEN 1524 ! 1525 predVBdRA(iq,:)=dptmp5 1526 ! 1527 ELSE 1528 ! 1529 predVBdRA(iq,:)=dptmp5*dqA(:)/dqAmic 1530 ! 1531 END IF 1532 ! 1533 END DO !iq 1534!$omp end parallel do 1535 ! 1536 ! Inner loop over unique atom pairs B in the simulation cell to compute pair contributions to dveffdR and dveffdh... 1537 ! 1538!$omp parallel do private(dqB,dqBs,dqBmic,ir,dk1,rhoB,drhoB,dVAdRB,dVBdRA, & 1539!$omp i,j,ib,ibs,spcutB,spdB,ioff,boff), & 1540!$omp reduction(+:dveffdR),reduction(-:dveffdh) 1541 DO ipair=1,npair(ia) 1542 ! 1543 ! Connect pair number with atom... 1544 ! 1545 ib=pair(ipair,ia) 1546 ! 1547 ! Connect atom type with species-dependent quantities... 1548 ! 1549 ibs=ityp(ib) 1550 ! 1551 ! Transfer species-specific cutoff to spcutB... 1552 ! 1553 spcutB=spdata(ibs,1) 1554 ! 1555 ! Precompute inverse of species-specific linear grid spacing (replaces / with * during interpolation)... 1556 ! 1557 spdB=1.0_DP/spdata(ibs,2) 1558 ! 1559 ! Determine atom B offsets for using gomegaAr bit array screening below... 1560 ! 1561 ioff=((ib-1)/bsint)+1 ! integer offset for gomegaAr bit array 1562 boff=(ib-((ioff-1)*bsint))-1 ! bit offset for gomegaAr bit array 1563 ! 1564 ! Inner loop over points in the (pre-screened) reduced spherical atomic integration domain for atom A to compute 1565 ! non-self-derivatives (dV(A)/dR(B) and dV(B)/dR(A)) and non-self-contributions to dV(A)/dh and dV(B)/dh in the overlapping integration domain... 1566 ! 1567 DO iq=1,NsomegaAr(ia) 1568 ! 1569 ! Determine if atom B contributes to the given point on the reduced spherical atomic integration domain on atom A (using gomegaAr bit array)... 1570 ! 1571 IF (BTEST(gomegaAr(iq,ioff,iproc),boff)) THEN 1572 ! 1573 ! Compute distance between grid point and atom B according to minimum image convention (MIC)... 1574 ! 1575 dqB(1)=dqxyzr(iq,1)-atxyz(1,ib) ! r_iB = r_i - r_B 1576 dqB(2)=dqxyzr(iq,2)-atxyz(2,ib) ! r_iB = r_i - r_B 1577 dqB(3)=dqxyzr(iq,3)-atxyz(3,ib) ! r_iB = r_i - r_B 1578 ! 1579 dqBs(1)=ainv(1,1)*dqB(1)+ainv(1,2)*dqB(2)+ainv(1,3)*dqB(3) ! s_iB = h^-1 r_iB 1580 dqBs(2)=ainv(2,1)*dqB(1)+ainv(2,2)*dqB(2)+ainv(2,3)*dqB(3) ! s_iB = h^-1 r_iB 1581 dqBs(3)=ainv(3,1)*dqB(1)+ainv(3,2)*dqB(2)+ainv(3,3)*dqB(3) ! s_iB = h^-1 r_iB 1582 ! 1583 dqBs(1)=dqBs(1)-IDNINT(dqBs(1)) ! impose MIC on s_iB in range: [-0.5,+0.5] 1584 dqBs(2)=dqBs(2)-IDNINT(dqBs(2)) ! impose MIC on s_iB in range: [-0.5,+0.5] 1585 dqBs(3)=dqBs(3)-IDNINT(dqBs(3)) ! impose MIC on s_iB in range: [-0.5,+0.5] 1586 ! 1587 dqB(1)=h(1,1)*dqBs(1)+h(1,2)*dqBs(2)+h(1,3)*dqBs(3) ! r_iB = h s_iB (MIC) 1588 dqB(2)=h(2,1)*dqBs(1)+h(2,2)*dqBs(2)+h(2,3)*dqBs(3) ! r_iB = h s_iB (MIC) 1589 dqB(3)=h(3,1)*dqBs(1)+h(3,2)*dqBs(2)+h(3,3)*dqBs(3) ! r_iB = h s_iB (MIC) 1590 ! 1591 dqBmic=DSQRT(dqB(1)*dqB(1)+dqB(2)*dqB(2)+dqB(3)*dqB(3)) ! |r_i - r_B| (MIC) 1592 ! 1593 ! Final screening based on the (pre-screened) spherical atomic integration domain on atom B... 1594 ! 1595 IF (dqBmic.LE.spcutB) THEN 1596 ! 1597 ! Determine the index in the atomic linear equispaced grid such that grd(ir) <= dqB <= grd(ir+1) and distance between dqB and grd(ir)... 1598 ! 1599 ir=INT(dqBmic*spdB) 1600 dk1=dqBmic-spgrd(ibs,ir) 1601 ! 1602 ! Perform linear interpolation to obtain the value of the atomic pseudo-density and its derivative at the given grid point... 1603 ! 1604 rhoB=LIA(ibs,ir)+LIB(ibs,ir)*dk1 !rhoB at grid point via linear interpolation 1605 drhoB=dLIA(ibs,ir)+dLIB(ibs,ir)*dk1 !drhoB at grid point via linear interpolation 1606 ! 1607 ! Compute dVA/dpB integrand for p={x,y,z}... 1608 ! 1609 ! dVA/dpB = INT {(p-pB)/|r-rB|*[drho(|r-rB|)*|r-rA|^3*rho(|r-rA|)*rhotot(r)/rhosad(r)^2]} 1610 ! 1611 IF (dqBmic.LT.(1.0E-12_DP)) THEN 1612 ! 1613 dVAdRB(:)=predVAdRB(iq)*drhoB 1614 ! 1615 ELSE 1616 ! 1617 dVAdRB(:)=predVAdRB(iq)*drhoB*dqB(:)/dqBmic 1618 ! 1619 END IF 1620 ! 1621 ! Increment non-self-derivative dVA/dpB for p={x,y,z}... 1622 ! 1623 DO i=1,3 1624 ! 1625 dveffdR(ia,ib,i)=dveffdR(ia,ib,i)+dVAdRB(i) 1626 ! 1627 END DO !i 1628 ! 1629 ! Increment non-self-contribution to dVA/dhpq for p,q={x,y,z} from atom B... 1630 ! 1631 ! dVA/dhpq <-- INT {-(p-pB)*(qs-qsB)/|r-rB|*[drho(|r-rB|)*|r-rA|^3*rho(|r-rA|)*rhotot(r)/rhosad(r)^2]} 1632 ! 1633 DO i=1,3 1634 ! 1635 DO j=1,3 1636 ! 1637 dveffdh(ia,i,j)=dveffdh(ia,i,j)-dVAdRB(i)*dqBs(j) 1638 ! 1639 END DO !j 1640 ! 1641 END DO !i 1642 ! 1643 ! Compute dVB/dpA integrand for p={x,y,z}... 1644 ! 1645 ! dVB/dpA = INT {(p-pA)/|r-rA|*[drho(|r-rA|)*|r-rB|^3*rho(|r-rB|)*rhotot(r)/rhosad(r)^2]} 1646 ! 1647 dVBdRA(:)=predVBdRA(iq,:)*rhoB*dqBmic*dqBmic*dqBmic 1648 ! 1649 ! Increment non-self-derivative dVB/dpA for p={x,y,z} from atom A... 1650 ! 1651 DO i=1,3 1652 ! 1653 dveffdR(ib,ia,i)=dveffdR(ib,ia,i)+dVBdRA(i) 1654 ! 1655 END DO !i 1656 ! 1657 ! Increment non-self-contribution to dVB/dhpq for p,q={x,y,z} from atom A... 1658 ! 1659 ! dVB/dhpq <-- INT {-(p-pA)*(qs-qsA)/|r-rA|*[drho(|r-rA|)*|r-rB|^3*rho(|r-rB|)*rhotot(r)/rhosad(r)^2]} 1660 ! 1661 DO i=1,3 1662 ! 1663 DO j=1,3 1664 ! 1665 dveffdh(ib,i,j)=dveffdh(ib,i,j)-dVBdRA(i)*dqAxyzs(iq,j) 1666 ! 1667 END DO !j 1668 ! 1669 END DO !i 1670 ! 1671 END IF 1672 ! 1673 END IF !BTEST 1674 ! 1675 END DO !iq 1676 ! 1677 END DO !ipair 1678!$omp end parallel do 1679 ! 1680 ! Deallocate temporary arrays... 1681 ! 1682 IF (ALLOCATED(dqxyzr)) DEALLOCATE(dqxyzr) 1683 IF (ALLOCATED(dqAxyzs)) DEALLOCATE(dqAxyzs) 1684 IF (ALLOCATED(predVAdRB)) DEALLOCATE(predVAdRB) 1685 IF (ALLOCATED(predVBdRA)) DEALLOCATE(predVBdRA) 1686 ! 1687 END DO !iproc 1688 ! 1689 ! Apply final normalization of dVA/dR integrals... 1690 ! 1691 dveffdR=normr*dveffdR 1692 ! 1693 ! Apply final normalization of dVA/dhab integrals... 1694 ! 1695 dveffdh=normr*dveffdh 1696 ! 1697 ! Collect dveffdR and dveffdh over all processors and broadcast... 1698 ! 1699 CALL mp_sum(dveffdR,intra_image_comm) 1700 CALL mp_sum(dveffdh,intra_image_comm) 1701 ! 1702 CALL stop_clock('tsvdw_dveff') 1703 ! 1704 RETURN 1705 ! 1706 !-------------------------------------------------------------------------------------------------------------- 1707 END SUBROUTINE tsvdw_dveff 1708 !-------------------------------------------------------------------------------------------------------------- 1709 ! 1710 !-------------------------------------------------------------------------------------------------------------- 1711 SUBROUTINE tsvdw_effqnts() 1712 !-------------------------------------------------------------------------------------------------------------- 1713 ! 1714 IMPLICIT NONE 1715 ! 1716 ! Local variables 1717 ! 1718 INTEGER :: ia,ib,ias,ibs 1719 REAL(DP) :: vA,vB,num,den 1720 ! 1721 ! Initialization of base effective atomic quantities... 1722 ! 1723 ALLOCATE(dpeff(nat)); dpeff=0.0_DP 1724 ALLOCATE(R0eff(nat)); R0eff=0.0_DP 1725 ALLOCATE(C6AAeff(nat)); C6AAeff=0.0_DP 1726 ALLOCATE(C6ABeff(nat,nat)); C6ABeff=0.0_DP 1727 ! 1728 ! Population of base effective atomic quantities... 1729 ! 1730 DO ia=1,nat 1731 ! 1732 ! Connect atom type with species-dependent quantities... 1733 ! 1734 ias=ityp(ia) 1735 ! 1736 ! Precompute veff(A)/vfree(A) ratio... 1737 ! 1738 vA=(veff(ia)/vfree(ias)) 1739 ! 1740 ! Effective atomic static dipole polarizability array... 1741 ! dpeff(A)=[veff(A)/vfree(A)]*dpfree(A) 1742 ! 1743 dpeff(ia)=vA*dpfree(ias) 1744 ! 1745 ! Effective atomic vdW radius array... 1746 ! R0eff(A)=[veff(A)/vfree(A)]^1/3*R0free(A) 1747 ! 1748 R0eff(ia)=(vA**(1.0_DP/3.0_DP))*R0free(ias) 1749 ! 1750 ! Effective homonuclear C6 coefficient array... 1751 ! C6AAeff(A)=[veff(A)/vfree(A)]^2*C6AAfree(A) 1752 ! 1753 C6AAeff(ia)=(vA**(2.0_DP))*C6AAfree(ias) 1754 ! 1755 DO ib=1,nat 1756 ! 1757 ! Connect atom type with species-dependent quantities... 1758 ! 1759 ibs=ityp(ib) 1760 ! 1761 ! Precompute veff(B)/vfree(B) ratio... 1762 ! 1763 vB=(veff(ib)/vfree(ibs)) 1764 ! 1765 ! Effective heteronuclear C6 coefficient matrix... 1766 ! C6ABeff(A,B)=(veff(A)/vfree(A))*(veff(B)/vfree(B))*C6ABfree(A,B) 1767 ! 1768 C6ABeff(ia,ib)=(vA*vB)*C6ABfree(ias,ibs) 1769 ! 1770 END DO !ib 1771 ! 1772 END DO !ia 1773 ! 1774 RETURN 1775 ! 1776 !-------------------------------------------------------------------------------------------------------------- 1777 END SUBROUTINE tsvdw_effqnts 1778 !-------------------------------------------------------------------------------------------------------------- 1779 ! 1780 !-------------------------------------------------------------------------------------------------------------- 1781 SUBROUTINE tsvdw_energy() 1782 !-------------------------------------------------------------------------------------------------------------- 1783 ! 1784 IMPLICIT NONE 1785 ! 1786 ! Local variables 1787 ! 1788 LOGICAL :: periodic_converged 1789 INTEGER :: ia,ib,ic,ias,ibs,n_period,n1,n2,n3,i,j 1790 1791 REAL(DP) :: dAB(3),dAB2(3),dsAB(3),dABimg,dABimg2,dABimgn1,dABimgn2,dABimgn5,dABimgn6 1792 REAL(DP) :: FDV0,FDR0,FCV0,FRR0,FDV1,FCVA1,FCVB1,FDVA2,FDVB2,FDR2,FRR2,FCVA2,FCVB2,FDVi,FDRi(3),FDRii(3,3),FCVi,FRRi(3),FRRii(3,3) 1793 REAL(DP) :: EtsvdW_period,RAB0,edamp,fdamp,fdamp2,D1A,D1B,D2A,D2B,D12A,D12B,dptmp1,dptmp2,vtmp1(3),vtmp2(3) 1794 REAL(DP), DIMENSION(:), ALLOCATABLE :: predveffAdn_period 1795 REAL(DP), DIMENSION(:,:), ALLOCATABLE :: FtsvdW_period,HtsvdW_period 1796 ! 1797 CALL start_clock('tsvdw_energy') 1798 ! 1799 ! Initialize total TS-vdW energy, ion force, and cell force ... 1800 ! 1801 EtsvdW=0.0_DP; FtsvdW=0.0_DP; HtsvdW=0.0_DP 1802 ! 1803 ! Allocate and initialize TS-vdW dispersion potential prefactor... 1804 ! 1805 ALLOCATE(predveffAdn(nat)); predveffAdn=0.0_DP 1806 ! 1807 ! Allocate and initialize periodic contributions to TS-vdW ionic forces, cell forces, and dispersion potential prefactor... 1808 ! 1809 ALLOCATE(FtsvdW_period(3,nat)); FtsvdW_period=0.0_DP 1810 ALLOCATE(HtsvdW_period(3,3)); HtsvdW_period=0.0_DP 1811 ALLOCATE(predveffAdn_period(nat)); predveffAdn_period=0.0_DP 1812 ! 1813 ! Precompute quantities outside all loops... 1814 ! 1815 FDR0=ddamp/(2.0_DP*sR) 1816 FDV0=-FDR0/3.0_DP 1817 FCV0=0.5_DP 1818 FRR0=-3.0_DP 1819 ! 1820 ! For periodic systems, converge the energy with respect to neighboring images... 1821 ! 1822 n_period=0 1823 periodic_converged=.FALSE. 1824 ! 1825 DO WHILE (.NOT.periodic_converged) 1826 ! 1827 EtsvdW_period=0.0_DP 1828 FtsvdW_period=0.0_DP 1829 HtsvdW_period=0.0_DP 1830 predveffAdn_period=0.0_DP 1831 ! 1832 ! Outer loop over atoms A... 1833 ! 1834 DO iproc=1,nstates(me) 1835 ! 1836 ! Connect processor number with atom... 1837 ! 1838 ia=me+nproc_image*(iproc-1) 1839 ! 1840 ! Connect atom type with species-dependent quantities... 1841 ! 1842 ias=ityp(ia) 1843 ! 1844 ! Precompute quantities outside loop over B... 1845 ! 1846 FDV1=R0free(ias)/(vfree(ias)**(1.0_DP/3.0_DP)*veff(ia)**(2.0_DP/3.0_DP)) 1847 FCVA1=1.0_DP/vfree(ias) 1848 FCVB1=veff(ia)*FCVA1 1849 ! 1850 ! Inner loop over atoms B... 1851 ! 1852 1853!$omp parallel private(ibs,RAB0,FRR2,FDR2,FDVA2,FDVB2,FCVB2,FCVA2, & 1854!$omp dAB,dAB2,FDVi,FDRi,FDRii,FCVi,FRRi,FRRii,n1,n2,n3,dsAB,dABimg2, & 1855!$omp dABimg,dABimgn1,dABimgn2,dABimgn5,dABimgn6,edamp,fdamp,fdamp2,dptmp1, & 1856!$omp dptmp2,i,j,vtmp1,vtmp2,D1A,D2A,D1B,D2B,D12A,D12B,ic), & 1857!$omp reduction(-:EtsvdW_period),reduction(+:FtsvdW_period), & 1858!$omp reduction(+:HtsvdW_period),reduction(-:predveffAdn_period) 1859!$omp do 1860 DO ib=1,nat 1861 ! 1862 ! Connect atom type with species-dependent quantities... 1863 ! 1864 ibs=ityp(ib) 1865 ! 1866 ! Compute RAB0 as the sum of the effective vdW radii of atoms A and B... 1867 ! 1868 RAB0=R0eff(ia)+R0eff(ib) 1869 ! 1870 ! Precompute quantities outside loop over image cells... 1871 ! 1872 FRR2=C6ABeff(ia,ib) 1873 FDR2=FRR2/RAB0 1874 FDVA2=FDR2/RAB0 1875 FDVB2=FDVA2*R0free(ibs)/(vfree(ibs)**(1.0_DP/3.0_DP)*veff(ib)**(2.0_DP/3.0_DP)) 1876 FCVB2=C6ABfree(ias,ibs)/vfree(ibs) 1877 FCVA2=FCVB2*veff(ib) 1878 ! 1879 ! Compute distance between atom A and atom B (according to the minimum image convention)... 1880 ! 1881 dAB(1)=atxyz(1,ia)-atxyz(1,ib) ! r_AB = r_A - r_B 1882 dAB(2)=atxyz(2,ia)-atxyz(2,ib) ! r_AB = r_A - r_B 1883 dAB(3)=atxyz(3,ia)-atxyz(3,ib) ! r_AB = r_A - r_B 1884 ! 1885 dAB2(1)=ainv(1,1)*dAB(1)+ainv(1,2)*dAB(2)+ainv(1,3)*dAB(3) ! s_AB = h^-1 r_AB 1886 dAB2(2)=ainv(2,1)*dAB(1)+ainv(2,2)*dAB(2)+ainv(2,3)*dAB(3) ! s_AB = h^-1 r_AB 1887 dAB2(3)=ainv(3,1)*dAB(1)+ainv(3,2)*dAB(2)+ainv(3,3)*dAB(3) ! s_AB = h^-1 r_AB 1888 ! 1889 dAB2(1)=dAB2(1)-IDNINT(dAB2(1)) ! impose MIC on s_AB in range: [-0.5,+0.5] 1890 dAB2(2)=dAB2(2)-IDNINT(dAB2(2)) ! impose MIC on s_AB in range: [-0.5,+0.5] 1891 dAB2(3)=dAB2(3)-IDNINT(dAB2(3)) ! impose MIC on s_AB in range: [-0.5,+0.5] 1892 ! 1893 ! Initialize image-summed matrix elements... 1894 ! 1895 FDVi=0.0_DP; FDRi=0.0_DP; FDRii=0.0_DP; FCVi=0.0_DP; FRRi=0.0_DP; FRRii=0.0_DP 1896 ! 1897 ! Loop over image cells... 1898 ! 1899 DO n1=-n_period,n_period 1900 ! 1901 DO n2=-n_period,n_period 1902 ! 1903 DO n3=-n_period,n_period 1904 ! 1905 IF ((ABS(n1).EQ.n_period).OR.(ABS(n2).EQ.n_period).OR.(ABS(n3).EQ.n_period)) THEN 1906 ! 1907 ! Recover MIC distance between atom A and atom B in crystal coordinates... 1908 ! 1909 dsAB(1)=dAB2(1) ! s_AB (MIC) 1910 dsAB(2)=dAB2(2) ! s_AB (MIC) 1911 dsAB(3)=dAB2(3) ! s_AB (MIC) 1912 ! 1913 ! Increment MIC distance in crystal coordinates... 1914 ! 1915 dsAB(1)=dsAB(1)+DBLE(n1) ! s_AB (incremented, MIC only if n_period == 0) 1916 dsAB(2)=dsAB(2)+DBLE(n2) ! s_AB (incremented, MIC only if n_period == 0) 1917 dsAB(3)=dsAB(3)+DBLE(n3) ! s_AB (incremented, MIC only if n_period == 0) 1918 ! 1919 ! Convert incremented distance back into cartesian coordinates... 1920 ! 1921 dAB(1)=h(1,1)*dsAB(1)+h(1,2)*dsAB(2)+h(1,3)*dsAB(3) ! r_AB = h s_AB (MIC only if n_period == 0) 1922 dAB(2)=h(2,1)*dsAB(1)+h(2,2)*dsAB(2)+h(2,3)*dsAB(3) ! r_AB = h s_AB (MIC only if n_period == 0) 1923 dAB(3)=h(3,1)*dsAB(1)+h(3,2)*dsAB(2)+h(3,3)*dsAB(3) ! r_AB = h s_AB (MIC only if n_period == 0) 1924 ! 1925 ! Compute incremented distance between atom A and atom B... 1926 ! 1927 dABimg2=dAB(1)*dAB(1)+dAB(2)*dAB(2)+dAB(3)*dAB(3) 1928 dABimg=DSQRT(dABimg2) 1929 ! 1930 ! Precompute inverse powers of incremented distance between atom A and atom B... 1931 ! 1932 IF ( dABimg > 0.0_dp ) THEN 1933 dABimgn1=1.0_DP/dABimg 1934 dABimgn2=dABimgn1*dABimgn1 1935 dABimgn5=dABimgn2*dABimgn2*dABimgn1 1936 dABimgn6=dABimgn5*dABimgn1 1937 ELSE 1938 dABimgn1=0.0_DP 1939 dABimgn2=0.0_DP 1940 dABimgn5=0.0_DP 1941 dABimgn6=0.0_DP 1942 END IF 1943 ! 1944 ! Precompute damping function (fdamp) and damping function exponential (edamp)... 1945 ! 1946 edamp=EXP(-ddamp*(dABimg/(sR*RAB0)-1.0_DP)) 1947 fdamp=1.0_DP/(1.0_DP+edamp) 1948 fdamp2=fdamp*fdamp 1949 ! 1950 ! Apply delta[ia;ib] x delta[n1,n2,n3;0,0,0] conditional... 1951 ! 1952 IF (n_period.EQ.0.AND.ia.EQ.ib) THEN 1953 ! 1954 ! Do not include self-interaction in the simulation cell... 1955 ! 1956 FDVi=FDVi+0.0_DP 1957 FDRi=FDRi+0.0_DP 1958 FDRii=FDRii+0.0_DP 1959 FCVi=FCVi+0.0_DP 1960 FRRi=FRRi+0.0_DP 1961 FRRii=FRRii+0.0_DP 1962 ! 1963 ELSE 1964 ! 1965 ! Increment image-summed matrix elements... 1966 ! 1967 dptmp1=edamp*fdamp2*dABimgn5 1968 FDVi=FDVi+dptmp1 1969 ! 1970 dptmp2=fdamp*dABimgn6 1971 FCVi=FCVi+dptmp2 1972 ! 1973 dptmp1=dptmp1*dABimgn2 1974 dptmp2=dptmp2*dABimgn2 1975 ! 1976 DO i=1,3 1977 ! 1978 vtmp1(i)=dptmp1*dAB(i) 1979 FDRi(i)=FDRi(i)+vtmp1(i) 1980 ! 1981 vtmp2(i)=dptmp2*dAB(i) 1982 FRRi(i)=FRRi(i)+vtmp2(i) 1983 ! 1984 DO j=1,3 1985 ! 1986 FDRii(i,j)=FDRii(i,j)+vtmp1(i)*dsAB(j) 1987 FRRii(i,j)=FRRii(i,j)+vtmp2(i)*dsAB(j) 1988 ! 1989 END DO 1990 ! 1991 END DO 1992 ! 1993 END IF 1994 ! 1995 END IF !n_period conditional 1996 ! 1997 END DO !n3 1998 ! 1999 END DO !n2 2000 ! 2001 END DO !n1 2002 ! 2003 ! Increment period energy via EtsvdWAB = - 1/2 * C6ABeff * FAB... 2004 ! 2005 EtsvdW_period=EtsvdW_period-(FCV0*FRR2*FCVi) 2006 ! 2007 ! Increment dispersion potential (predveffAdn) prefactor... 2008 ! 2009 ! predveffAdn(A) = (d * R0freeA * C6ABeff * edamp * fdamp^2) / (6 * sR * vfreeA^1/3 * veffA^2/3 * RAB0^2 * RAB^5) 2010 ! - (C6ABfree * veffB * fdamp) / (2 * vfreeA * vfreeB * RAB^6) 2011 ! 2012 ! predveffAdn(B) = (d * R0freeB * C6ABeff * edamp * fdamp^2) / (6 * sR * vfreeB^1/3 * veffB^2/3 * RAB0^2 * RAB^5) 2013 ! - (C6ABfree * veffA * fdamp) / (2 * vfreeA * vfreeB * RAB^6) 2014 ! 2015 predveffAdn_period(ia)=predveffAdn_period(ia)-(FDV0*FDV1*FDVA2*FDVi+FCV0*FCVA1*FCVA2*FCVi) 2016 predveffAdn_period(ib)=predveffAdn_period(ib)-(FDV0*FDVB2*FDVi+FCV0*FCVB1*FCVB2*FCVi) 2017 ! 2018 ! Increment effective volume derivative contributions to ionic and cell forces... 2019 ! 2020 ! (dfdamp/dVA) --> D1A = - (d * R0freeA * C6ABeff * edamp * fdamp^2) / (6 * sR * vfreeA^1/3 * veffA^2/3 * RAB0^2 * RAB^5) 2021 ! 2022 ! (dfdamp/dVB) --> D1B = - (d * R0freeB * C6ABeff * edamp * fdamp^2) / (6 * sR * vfreeB^1/3 * veffB^2/3 * RAB0^2 * RAB^5) 2023 ! 2024 ! (dC6AB/dVA) --> D2A = (C6ABfree * veffB * fdamp) / (2 * vfreeA * vfreeB * RAB^6) 2025 ! 2026 ! (dC6AB/dVB) --> D2B = (C6ABfree * veffA * fdamp) / (2 * vfreeA * vfreeB * RAB^6) 2027 ! 2028 D1A=FDV0*FDV1*FDVA2*FDVi ! (dfdamp/dVA) 2029 D2A=FCV0*FCVA1*FCVA2*FCVi ! (dC6AB/dVA) 2030 ! 2031 D1B=FDV0*FDVB2*FDVi ! (dfdamp/dVB) 2032 D2B=FCV0*FCVB1*FCVB2*FCVi ! (dC6AB/dVB) 2033 ! 2034 D12A=D1A+D2A; D12B=D1B+D2B 2035 ! 2036 DO i=1,3 2037 ! 2038 DO ic=1,nat 2039 ! 2040 FtsvdW_period(i,ic)=FtsvdW_period(i,ic)+(dveffdR(ia,ic,i)*D12A+dveffdR(ib,ic,i)*D12B) 2041 ! 2042 END DO 2043 ! 2044 DO j=1,3 2045 ! 2046 HtsvdW_period(i,j)=HtsvdW_period(i,j)+(dveffdh(ia,i,j)*D12A+dveffdh(ib,i,j)*D12B) 2047 ! 2048 END DO 2049 ! 2050 END DO 2051 ! 2052 ! Increment RAB derivative contributions to ionic and cell forces... 2053 ! 2054 ! (dfdamp/dRA) --> D1A = (d * C6ABeff * edamp * fdamp^2) / (2 * sR * RAB0 * RAB^7) 2055 ! 2056 ! (dfdamp/dRB) --> D1B = - D1A 2057 ! 2058 ! (dRAB^-6/dRA) --> D2A = - (3 * C6ABeff * fdamp) / (RAB^8) 2059 ! 2060 ! (dRAB^-6/dRB) --> D2B = - D2A 2061 ! 2062 D1A=FDR0*FDR2 ! (dfdamp/dRA) 2063 D2A=FRR0*FRR2 ! (dRAB^-6/dRA) 2064 ! 2065 ! N.B.: Manually zero out the force contribution from an atom in the simulation cell and 2066 ! any of its images (this applies to distance derivatives NOT volume derivatives)... 2067 ! 2068 IF (ia.NE.ib) THEN 2069 ! 2070 DO i=1,3 2071 ! 2072 FtsvdW_period(i,ia)=FtsvdW_period(i,ia)+(D1A*FDRi(i)+D2A*FRRi(i)) 2073 FtsvdW_period(i,ib)=FtsvdW_period(i,ib)+(-D1A*FDRi(i)-D2A*FRRi(i)) 2074 ! 2075 DO j=1,3 2076 ! 2077 HtsvdW_period(i,j)=HtsvdW_period(i,j)+(D1A*FDRii(i,j)+D2A*FRRii(i,j)) 2078 ! 2079 END DO 2080 ! 2081 END DO 2082 ! 2083 END IF 2084 ! 2085 END DO !ib 2086!$omp end do 2087!$omp end parallel 2088 ! 2089 END DO !iproc 2090 ! 2091 ! Synchronize n_period contribution from all processors... 2092 ! 2093 CALL mp_sum(EtsvdW_period,intra_image_comm) 2094 CALL mp_sum(FtsvdW_period,intra_image_comm) 2095 CALL mp_sum(HtsvdW_period,intra_image_comm) 2096 CALL mp_sum(predveffAdn_period,intra_image_comm) 2097 ! 2098 ! Increment total quantities... 2099 ! 2100 EtsvdW=EtsvdW+EtsvdW_period !EvdW 2101 FtsvdW=FtsvdW+FtsvdW_period !(-dE/dR) 2102 HtsvdW=HtsvdW-HtsvdW_period !(dE/dh) 2103 predveffAdn=predveffAdn+predveffAdn_period !(dE/dVA) & (dE/dVB) 2104 ! 2105 ! DEBUGGING 2106 !WRITE(stdout,'(I10,2F25.12)') n_period,EtsvdW_period,EtsvdW 2107 ! DEBUGGING 2108 ! 2109 ! Periodic convergence loop conditionals... 2110 ! 2111 IF ( vdw_isolated .OR. (ABS(EtsvdW_period) <= vdw_econv_thr) ) periodic_converged=.TRUE. 2112 ! 2113 n_period=n_period+1 2114 ! 2115 END DO !convergence loop 2116 ! 2117 ! Deallocate temporary arrays... 2118 ! 2119 IF (ALLOCATED(FtsvdW_period)) DEALLOCATE(FtsvdW_period) 2120 IF (ALLOCATED(HtsvdW_period)) DEALLOCATE(HtsvdW_period) 2121 IF (ALLOCATED(predveffAdn_period)) DEALLOCATE(predveffAdn_period) 2122 ! 2123 CALL stop_clock('tsvdw_energy') 2124 ! 2125 !! DEBUGGING 2126 !WRITE(stdout,'(3X,"<START veff ARRAY>")') 2127 !DO ia=1,nat 2128 ! WRITE(stdout,'(5X,I5,F25.12)') ia,veff(ia) 2129 !END DO 2130 !WRITE(stdout,'(3X,"<END veff ARRAY>")') 2131 !! 2132 !WRITE(stdout,'(3X,"<START FtsvdW ARRAY>")') 2133 !DO ia=1,nat 2134 ! WRITE(stdout,'(5X,I5,3F25.12)') ia,FtsvdW(:,ia) 2135 !END DO 2136 !WRITE(stdout,'(3X,"<END FtsvdW ARRAY>")') 2137 !! 2138 !WRITE(stdout,'(3X,"<START HtsvdW ARRAY>")') 2139 !DO i=1,3 2140 ! WRITE(stdout,'(5X,3F25.12)') HtsvdW(i,:) 2141 !END DO 2142 !WRITE(stdout,'(3X,"<END HtsvdW ARRAY>")') 2143 !! DEBUGGING 2144 ! 2145 RETURN 2146 ! 2147 !-------------------------------------------------------------------------------------------------------------- 2148 END SUBROUTINE tsvdw_energy 2149 !-------------------------------------------------------------------------------------------------------------- 2150 ! 2151 !-------------------------------------------------------------------------------------------------------------- 2152 SUBROUTINE tsvdw_wfforce() 2153 !-------------------------------------------------------------------------------------------------------------- 2154 ! 2155 IMPLICIT NONE 2156 ! 2157 ! Local variables 2158 ! 2159 INTEGER :: ia,ip,iq,off1 2160 REAL(DP), DIMENSION(:), ALLOCATABLE :: UtsvdWA 2161 ! 2162 CALL start_clock('tsvdw_wfforce') 2163 ! 2164 ! Initialization of UtsvdwA array... 2165 ! 2166 ALLOCATE(UtsvdWA(nr1*nr2*nr3)); UtsvdWA=0.0_DP 2167 ! 2168 ! Loop over atoms and populate UtsvdWA from predveffAdn and dveffAdn... 2169 ! 2170 DO iproc=1,nstates(me) 2171 ! 2172 ! Connect processor number with atom... 2173 ! 2174 ia=me+nproc_image*(iproc-1) 2175 ! 2176 ! Loop over points in the (pre-screened) spherical atomic integration domain... 2177 ! 2178!$omp parallel do private(off1) 2179 DO iq=1,NsomegaA(ia) 2180 ! 2181 off1=somegaA(iq,1,iproc)+(somegaA(iq,2,iproc)-1)*nr1+(somegaA(iq,3,iproc)-1)*nr1*nr2 !global offset [nr1,nr2,nr3] 2182 UtsvdWA(off1)=UtsvdWA(off1)+predveffAdn(ia)*dveffAdn(iq,iproc) 2183 ! 2184 END DO !iq 2185!$omp end parallel do 2186 ! 2187 END DO !iproc 2188 ! 2189 ! Collect UtsvdWA over all processors and broadcast... 2190 ! 2191 CALL mp_sum(UtsvdWA,intra_image_comm) 2192 ! 2193 ! Partition out dispersion potential consistent with slabs of the charge density... 2194 ! 2195 IF (dfftp%nr3p(me_bgrp+1).NE.0) THEN 2196 ! 2197!$omp parallel do 2198 DO ip=1,dfftp%nr3p(me_bgrp+1)*nr1*nr2 2199 ! 2200 UtsvdW(ip)=UtsvdWA(ip+rdispls(me_bgrp+1)) 2201 ! 2202 END DO 2203!$omp end parallel do 2204 ! 2205 END IF 2206 ! 2207 ! Deallocate temporary arrays... 2208 ! 2209 IF (ALLOCATED(UtsvdWA)) DEALLOCATE(UtsvdWA) 2210 ! 2211 CALL stop_clock('tsvdw_wfforce') 2212 ! 2213 RETURN 2214 ! 2215 !-------------------------------------------------------------------------------------------------------------- 2216 END SUBROUTINE tsvdw_wfforce 2217 !-------------------------------------------------------------------------------------------------------------- 2218 ! 2219 !-------------------------------------------------------------------------------------------------------------- 2220 SUBROUTINE tsvdw_cleanup() 2221 !-------------------------------------------------------------------------------------------------------------- 2222 ! 2223 IMPLICIT NONE 2224 ! 2225 ! Deallocate tsvdw_calculate specific arrays... 2226 ! 2227 IF (ALLOCATED(atxyz)) DEALLOCATE(atxyz) 2228 IF (ALLOCATED(rhosad)) DEALLOCATE(rhosad) 2229 IF (ALLOCATED(rhotot)) DEALLOCATE(rhotot) 2230 IF (ALLOCATED(veff)) DEALLOCATE(veff) 2231 IF (ALLOCATED(dpeff)) DEALLOCATE(dpeff) 2232 IF (ALLOCATED(R0eff)) DEALLOCATE(R0eff) 2233 IF (ALLOCATED(C6AAeff)) DEALLOCATE(C6AAeff) 2234 IF (ALLOCATED(C6ABeff)) DEALLOCATE(C6ABeff) 2235 IF (ALLOCATED(dveffdR)) DEALLOCATE(dveffdR) 2236 IF (ALLOCATED(dveffdh)) DEALLOCATE(dveffdh) 2237 IF (ALLOCATED(somegaA)) DEALLOCATE(somegaA) 2238 IF (ALLOCATED(somegaAr)) DEALLOCATE(somegaAr) 2239 IF (ALLOCATED(gomegaAr)) DEALLOCATE(gomegaAr) 2240 IF (ALLOCATED(NsomegaA)) DEALLOCATE(NsomegaA) 2241 IF (ALLOCATED(NsomegaAr)) DEALLOCATE(NsomegaAr) 2242 IF (ALLOCATED(nstates)) DEALLOCATE(nstates) 2243 IF (ALLOCATED(sdispls)) DEALLOCATE(sdispls) 2244 IF (ALLOCATED(sendcount)) DEALLOCATE(sendcount) 2245 IF (ALLOCATED(rdispls)) DEALLOCATE(rdispls) 2246 IF (ALLOCATED(recvcount)) DEALLOCATE(recvcount) 2247 IF (ALLOCATED(istatus)) DEALLOCATE(istatus) 2248 IF (ALLOCATED(npair)) DEALLOCATE(npair) 2249 IF (ALLOCATED(pair)) DEALLOCATE(pair) 2250 IF (ALLOCATED(dveffAdn)) DEALLOCATE(dveffAdn) 2251 IF (ALLOCATED(predveffAdn)) DEALLOCATE(predveffAdn) 2252 ! 2253 RETURN 2254 ! 2255 !-------------------------------------------------------------------------------------------------------------- 2256 END SUBROUTINE tsvdw_cleanup 2257 !-------------------------------------------------------------------------------------------------------------- 2258 ! 2259 !-------------------------------------------------------------------------------------------------------------- 2260 SUBROUTINE tsvdw_finalize() 2261 !-------------------------------------------------------------------------------------------------------------- 2262 ! 2263 IMPLICIT NONE 2264 ! 2265 ! Deallocate module-specific arrays... 2266 ! 2267 IF (ALLOCATED(UtsvdW)) DEALLOCATE(UtsvdW) 2268 IF (ALLOCATED(FtsvdW)) DEALLOCATE(FtsvdW) 2269 IF (ALLOCATED(HtsvdW)) DEALLOCATE(HtsvdW) 2270 IF (ALLOCATED(VefftsvdW))DEALLOCATE(VefftsvdW) 2271 IF (ALLOCATED(vfree)) DEALLOCATE(vfree) 2272 IF (ALLOCATED(dpfree)) DEALLOCATE(dpfree) 2273 IF (ALLOCATED(R0free)) DEALLOCATE(R0free) 2274 IF (ALLOCATED(C6AAfree)) DEALLOCATE(C6AAfree) 2275 IF (ALLOCATED(C6ABfree)) DEALLOCATE(C6ABfree) 2276 IF (ALLOCATED(spgrd)) DEALLOCATE(spgrd) 2277 IF (ALLOCATED(sprho)) DEALLOCATE(sprho) 2278 IF (ALLOCATED(spdrho)) DEALLOCATE(spdrho) 2279 IF (ALLOCATED(spdata)) DEALLOCATE(spdata) 2280 IF (ALLOCATED(LIA)) DEALLOCATE(LIA) 2281 IF (ALLOCATED(LIB)) DEALLOCATE(LIB) 2282 IF (ALLOCATED(dLIA)) DEALLOCATE(dLIA) 2283 IF (ALLOCATED(dLIB)) DEALLOCATE(dLIB) 2284 ! 2285 RETURN 2286 ! 2287 !-------------------------------------------------------------------------------------------------------------- 2288 END SUBROUTINE tsvdw_finalize 2289 !-------------------------------------------------------------------------------------------------------------- 2290 ! 2291 !-------------------------------------------------------------------------------------------------------------- 2292 SUBROUTINE Num1stDer(r,f,N,h,df) 2293 !-------------------------------------------------------------------------------------------------------------- 2294 ! 2295 IMPLICIT NONE 2296 ! 2297 ! I/O variables 2298 ! 2299 INTEGER :: N 2300 REAL(DP) :: h,r(N),f(N),df(N) 2301 ! 2302 ! Local variables 2303 ! 2304 INTEGER, PARAMETER :: Ndp=7 !using 7-point formulae... 2305 INTEGER :: ip,ir 2306 INTEGER :: A1(Ndp),A2(Ndp),A3(Ndp),A4(Ndp),A5(Ndp),A6(Ndp),A7(Ndp) 2307 REAL(DP) :: dsum 2308 ! 2309 ! Populate Bickley coefficient vectors (Math. Gaz.,v25,p19-27,1941) according to cases... 2310 ! 2311 DATA A1/ -1764_DP, 4320_DP, -5400_DP, 4800_DP, -2700_DP, 864_DP, -120_DP / 2312 DATA A2/ -120_DP, -924_DP, 1800_DP, -1200_DP, 600_DP, -180_DP, 24_DP / 2313 DATA A3/ 24_DP, -288_DP, -420_DP, 960_DP, -360_DP, 96_DP, -12_DP / 2314 DATA A4/ -12_DP, 108_DP, -540_DP, 0_DP, 540_DP, -108_DP, 12_DP / 2315 DATA A5/ 12_DP, -96_DP, 360_DP, -960_DP, 420_DP, 288_DP, -24_DP / 2316 DATA A6/ -24_DP, 180_DP, -600_DP, 1200_DP, -1800_DP, 924_DP, 120_DP / 2317 DATA A7/ 120_DP, -864_DP, 2700_DP, -4800_DP, 5400_DP, -4320_DP, 1764_DP / 2318 ! 2319 ! Compute first derivative on linear mesh and then transform back to radial/exponential grid... 2320 ! 2321 DO ir=1,N 2322 ! 2323 dsum=0.0_DP 2324 ! 2325 ! Deal with different cases one-by-one... 2326 ! 2327 IF (ir.EQ.1) THEN 2328 DO ip=1,Ndp 2329 dsum=dsum+A1(ip)*f(ir-1+ip) 2330 END DO 2331 ELSE IF (ir.EQ.2) THEN 2332 DO ip=1,Ndp 2333 dsum=dsum+A2(ip)*f(ir-2+ip) 2334 END DO 2335 ELSE IF (ir.EQ.3) THEN 2336 DO ip=1,Ndp 2337 dsum=dsum+A3(ip)*f(ir-3+ip) 2338 END DO 2339 ELSE IF (ir.GE.4.AND.ir.LE.N-3) THEN 2340 DO ip=1,Ndp 2341 dsum=dsum+A4(ip)*f(ir-4+ip) 2342 END DO 2343 ELSE IF (ir.EQ.N-2) THEN 2344 DO ip=1,Ndp 2345 dsum=dsum+A5(ip)*f(ir-5+ip) 2346 END DO 2347 ELSE IF (ir.EQ.N-1) THEN 2348 DO ip=1,Ndp 2349 dsum=dsum+A6(ip)*f(ir-6+ip) 2350 END DO 2351 ELSE IF (ir.EQ.N) THEN 2352 DO ip=1,Ndp 2353 dsum=dsum+A7(ip)*f(ir-7+ip) 2354 END DO 2355 ELSE 2356 WRITE(stdout,'("Error in Num1stDer subroutine...")') 2357 END IF 2358 ! 2359 ! Final Normalization... 2360 ! 2361 df(ir)=dsum/(720.0_DP*h) 2362 ! 2363 END DO !ir 2364 ! 2365 RETURN 2366 ! 2367 !-------------------------------------------------------------------------------------------------------------- 2368 END SUBROUTINE Num1stDer 2369 !-------------------------------------------------------------------------------------------------------------- 2370 ! 2371 !-------------------------------------------------------------------------------------------------------------- 2372 SUBROUTINE CubSplCoeff(r,f,N,df,d2f) 2373 !-------------------------------------------------------------------------------------------------------------- 2374 ! 2375 IMPLICIT NONE 2376 ! 2377 ! I/O variables 2378 ! 2379 INTEGER :: N 2380 REAL(DP) :: r(N),f(N),df(N),d2f(N) 2381 ! 2382 ! Local variables 2383 ! 2384 INTEGER :: i,j 2385 REAL(DP) :: dy1,dyn,p,q,un,qn 2386 REAL(DP), DIMENSION(:), ALLOCATABLE :: work 2387 ! 2388 ! ---------------------------------------------------------------------------------------------------------------------------------- 2389 ! SYNOPSIS: Compute second derivatives at each of the atomic radial grid points using the cubic spline methodology (i.e., smooth & 2390 ! continuous piecewise first and second derivatives). These second derivatives will be utilized during cubic spline interpolation 2391 ! as a higher accuracy alternative to linear interpolation during the construction of the linear atomic grids. The two-parameter 2392 ! boundary conditions that will be utilized below are known as a clamped cubic spline in that the first derivative at both the 2393 ! first and last grid point were computed numerically and provided as input... 2394 ! ---------------------------------------------------------------------------------------------------------------------------------- 2395 ! 2396 ALLOCATE(work(N)); work=0.0_DP 2397 ! 2398 d2f=0.0_DP 2399 ! 2400 ! Enforce 'clamped' boundary condition at the first radial grid point... 2401 ! 2402 dy1=df(1) 2403 d2f(1)=-0.5_DP 2404 work(1)=(3.0_DP/(r(2)-r(1)))*((f(2)-f(1))/(r(2)-r(1))-dy1) 2405 ! 2406 ! Decomposition loop of the tridiagonal algorithm for the second derivatives... 2407 ! 2408 DO i=2,N-1 2409 p=(r(i)-r(i-1))/(r(i+1)-r(i-1)) 2410 q=p*d2f(i-1)+2.0_DP 2411 d2f(i)=(p-1.0_DP)/q 2412 work(i)=(f(i+1)-f(i))/(r(i+1)-r(i))-(f(i)-f(i-1))/(r(i)-r(i-1)) 2413 work(i)=(6.0_DP*work(i)/(r(i+1)-r(i-1))-p*work(i-1))/q 2414 END DO 2415 ! 2416 ! Enforce 'clamped' boundary condition at the last radial grid point... 2417 ! 2418 dyn=df(N) 2419 qn=0.5_DP 2420 un=(3.0_DP/(r(N)-r(N-1)))*(dyn-(f(N)-f(N-1))/(r(N)-r(N-1))) 2421 d2f(N)=(un-qn*work(N-1))/(qn*d2f(N-1)+1.0_DP) 2422 ! 2423 ! Back substitution loop of the tridiagonal algorithm for the second derivatives... 2424 ! 2425 DO j=N-1,1,-1 2426 d2f(j)=d2f(j)*d2f(j+1)+work(j) 2427 END DO 2428 ! 2429 ! Clean-up and return home... 2430 ! 2431 DEALLOCATE(work) 2432 ! 2433 RETURN 2434 ! 2435 !-------------------------------------------------------------------------------------------------------------- 2436 END SUBROUTINE CubSplCoeff 2437 !-------------------------------------------------------------------------------------------------------------- 2438 ! 2439 !-------------------------------------------------------------------------------------------------------------- 2440 SUBROUTINE GetVdWParam(atom,C6,alpha,R0) 2441 !-------------------------------------------------------------------------------------------------------------- 2442 ! 2443 IMPLICIT NONE 2444 ! 2445 ! I/O variables 2446 ! 2447 CHARACTER(LEN=2) :: atom 2448 REAL(DP) :: C6,alpha,R0 2449 ! 2450 SELECT CASE (atom) 2451 ! 2452 CASE ('H') 2453 alpha=4.500000_DP 2454 C6=6.500000_DP 2455 R0=3.100000_DP 2456 ! 2457 CASE ('He') 2458 alpha=1.380000_DP 2459 C6=1.460000_DP 2460 R0=2.650000_DP 2461 ! 2462 CASE ('Li') 2463 alpha=164.200000_DP 2464 C6=1387.000000_DP 2465 R0=4.160000_DP 2466 ! 2467 CASE ('Be') 2468 alpha=38.000000_DP 2469 C6=214.000000_DP 2470 R0=4.170000_DP 2471 ! 2472 CASE ('B') 2473 alpha=21.000000_DP 2474 C6=99.500000_DP 2475 R0=3.890000_DP 2476 ! 2477 CASE ('C') 2478 alpha=12.000000_DP 2479 C6=46.600000_DP 2480 R0=3.590000_DP 2481 ! 2482 CASE ('N') 2483 alpha=7.400000_DP 2484 C6=24.200000_DP 2485 R0=3.340000_DP 2486 ! 2487 CASE ('O') 2488 alpha=5.400000_DP 2489 C6=15.600000_DP 2490 R0=3.190000_DP 2491 ! 2492 CASE ('F') 2493 alpha=3.800000_DP 2494 C6=9.520000_DP 2495 R0=3.040000_DP 2496 ! 2497 CASE ('Ne') 2498 alpha=2.670000_DP 2499 C6=6.380000_DP 2500 R0=2.910000_DP 2501 ! 2502 CASE ('Na') 2503 alpha=162.700000_DP 2504 C6=1556.000000_DP 2505 R0=3.730000_DP 2506 ! 2507 CASE ('Mg') 2508 alpha=71.000000_DP 2509 C6=627.000000_DP 2510 R0=4.270000_DP 2511 ! 2512 CASE ('Al') 2513 alpha=60.000000_DP 2514 C6=528.000000_DP 2515 R0=4.330000_DP 2516 ! 2517 CASE ('Si') 2518 alpha=37.000000_DP 2519 C6=305.000000_DP 2520 R0=4.200000_DP 2521 ! 2522 CASE ('P') 2523 alpha=25.000000_DP 2524 C6=185.000000_DP 2525 R0=4.010000_DP 2526 ! 2527 CASE ('S') 2528 alpha=19.600000_DP 2529 C6=134.000000_DP 2530 R0=3.860000_DP 2531 ! 2532 CASE ('Cl') 2533 alpha=15.000000_DP 2534 C6=94.600000_DP 2535 R0=3.710000_DP 2536 ! 2537 CASE ('Ar') 2538 alpha=11.100000_DP 2539 C6=64.300000_DP 2540 R0=3.550000_DP 2541 ! 2542 CASE ('K') 2543 alpha=292.900000_DP 2544 C6=3897.000000_DP 2545 R0=3.710000_DP 2546 ! 2547 CASE ('Ca') 2548 alpha=160.000000_DP 2549 C6=2221.000000_DP 2550 R0=4.650000_DP 2551 ! 2552 CASE ('Sc') 2553 alpha=120.000000_DP 2554 C6=1383.000000_DP 2555 R0=4.590000_DP 2556 ! 2557 CASE ('Ti') 2558 alpha=98.000000_DP 2559 C6=1044.000000_DP 2560 R0=4.510000_DP 2561 ! 2562 CASE ('V') 2563 alpha=84.000000_DP 2564 C6=832.000000_DP 2565 R0=4.440000_DP 2566 ! 2567 CASE ('Cr') 2568 alpha=78.000000_DP 2569 C6=602.000000_DP 2570 R0=3.990000_DP 2571 ! 2572 CASE ('Mn') 2573 alpha=63.000000_DP 2574 C6=552.000000_DP 2575 R0=3.970000_DP 2576 ! 2577 CASE ('Fe') 2578 alpha=56.000000_DP 2579 C6=482.000000_DP 2580 R0=4.230000_DP 2581 ! 2582 CASE ('Co') 2583 alpha=50.000000_DP 2584 C6=408.000000_DP 2585 R0=4.180000_DP 2586 ! 2587 CASE ('Ni') 2588 alpha=48.000000_DP 2589 C6=373.000000_DP 2590 R0=3.820000_DP 2591 ! 2592 CASE ('Cu') 2593 alpha=42.000000_DP 2594 C6=253.000000_DP 2595 R0=3.760000_DP 2596 ! 2597 CASE ('Zn') 2598 alpha=40.000000_DP 2599 C6=284.000000_DP 2600 R0=4.020000_DP 2601 ! 2602 CASE ('Ga') 2603 alpha=60.000000_DP 2604 C6=498.000000_DP 2605 R0=4.190000_DP 2606 ! 2607 CASE ('Ge') 2608 alpha=41.000000_DP 2609 C6=354.000000_DP 2610 R0=4.200000_DP 2611 ! 2612 CASE ('As') 2613 alpha=29.000000_DP 2614 C6=246.000000_DP 2615 R0=4.110000_DP 2616 ! 2617 CASE ('Se') 2618 alpha=25.000000_DP 2619 C6=210.000000_DP 2620 R0=4.040000_DP 2621 ! 2622 CASE ('Br') 2623 alpha=20.000000_DP 2624 C6=162.000000_DP 2625 R0=3.930000_DP 2626 ! 2627 CASE ('Kr') 2628 alpha=16.800000_DP 2629 C6=129.600000_DP 2630 R0=3.820000_DP 2631 ! 2632 CASE ('Rb') 2633 alpha=319.200000_DP 2634 C6=4691.000000_DP 2635 R0=3.720000_DP 2636 ! 2637 CASE ('Sr') 2638 alpha=199.000000_DP 2639 C6=3170.000000_DP 2640 R0=4.540000_DP 2641 ! 2642 CASE ('Y') 2643 alpha=126.7370_DP 2644 C6=1968.580_DP 2645 R0=4.81510_DP 2646 ! 2647 CASE ('Zr') 2648 alpha=119.97_DP 2649 C6=1677.91_DP 2650 R0=4.53_DP 2651 ! 2652 CASE ('Nb') 2653 alpha=101.603_DP 2654 C6=1263.61_DP 2655 R0=4.2365_DP 2656 ! 2657 CASE ('Mo') 2658 alpha=88.4225785_DP 2659 C6=1028.73_DP 2660 R0=4.099_DP 2661 ! 2662 CASE ('Tc') 2663 alpha=80.083_DP 2664 C6=1390.87_DP 2665 R0=4.076_DP 2666 ! 2667 CASE ('Ru') 2668 alpha=65.8950_DP 2669 C6=609.754_DP 2670 R0=3.99530_DP 2671 ! 2672 CASE ('Rh') 2673 alpha=56.1_DP 2674 C6=469.0_DP 2675 R0=3.95_DP 2676 ! 2677 CASE ('Pd') 2678 alpha=23.680000_DP 2679 C6=157.500000_DP 2680 R0=3.66000_DP 2681 ! 2682 CASE ('Ag') 2683 alpha=50.600000_DP 2684 C6=339.000000_DP 2685 R0=3.820000_DP 2686 ! 2687 CASE ('Cd') 2688 alpha=39.7_DP 2689 C6=452.0_DP 2690 R0=3.99_DP 2691 ! 2692 CASE ('In') 2693 alpha=70.22000_DP 2694 C6=707.046000_DP 2695 R0=4.23198000_DP 2696 ! 2697 CASE ('Sn') 2698 alpha=55.9500_DP 2699 C6=587.41700_DP 2700 R0=4.303000_DP 2701 ! 2702 CASE ('Sb') 2703 alpha=43.671970_DP 2704 C6=459.322_DP 2705 R0=4.2760_DP 2706 ! 2707 CASE ('Te') 2708 alpha=37.65_DP 2709 C6=396.0_DP 2710 R0=4.22_DP 2711 ! 2712 CASE ('I') 2713 alpha=35.000000_DP 2714 C6=385.000000_DP 2715 R0=4.170000_DP 2716 ! 2717 CASE ('Xe') 2718 alpha=27.300000_DP 2719 C6=285.900000_DP 2720 R0=4.080000_DP 2721 ! 2722 CASE ('Cs') 2723 alpha=427.12_DP 2724 C6=6582.08_DP 2725 R0=3.78_DP 2726 ! 2727 CASE ('Ba') 2728 alpha=275.0_DP 2729 C6=5727.0_DP 2730 R0=4.77_DP 2731 ! 2732 CASE ('Hf') 2733 alpha=99.52_DP 2734 C6=1274.8_DP 2735 R0=4.21_DP 2736 ! 2737 CASE ('Ta') 2738 alpha=82.53_DP 2739 C6=1019.92_DP 2740 R0=4.15_DP 2741 ! 2742 CASE ('W') 2743 alpha=71.041_DP 2744 C6=847.93_DP 2745 R0=4.08_DP 2746 ! 2747 CASE ('Re') 2748 alpha=63.04_DP 2749 C6=710.2_DP 2750 R0=4.02_DP 2751 ! 2752 CASE ('Os') 2753 alpha=55.055_DP 2754 C6=596.67_DP 2755 R0=3.84_DP 2756 ! 2757 CASE ('Ir') 2758 alpha=42.51_DP 2759 C6=359.1_DP 2760 R0=4.00_DP 2761 ! 2762 CASE ('Pt') 2763 alpha=39.68_DP 2764 C6=347.1_DP 2765 R0=3.92_DP 2766 ! 2767 CASE ('Au') 2768 alpha=36.5_DP 2769 C6=298.0_DP 2770 R0=3.86_DP 2771 ! 2772 CASE ('Hg') 2773 alpha=33.9_DP 2774 C6=392.0_DP 2775 R0=3.98_DP 2776 ! 2777 CASE ('Tl') 2778 alpha=69.92_DP 2779 C6=717.44_DP 2780 R0=3.91_DP 2781 ! 2782 CASE ('Pb') 2783 alpha=61.8_DP 2784 C6=697.0_DP 2785 R0=4.31_DP 2786 ! 2787 CASE ('Bi') 2788 alpha=49.02_DP 2789 C6=571.0_DP 2790 R0=4.32_DP 2791 ! 2792 CASE ('Po') 2793 alpha=45.013_DP 2794 C6=530.92_DP 2795 R0=4.097_DP 2796 ! 2797 CASE ('At') 2798 alpha=38.93_DP 2799 C6=457.53_DP 2800 R0=4.07_DP 2801 ! 2802 CASE ('Rn') 2803 alpha=33.54_DP 2804 C6=390.63_DP 2805 R0=4.23_DP 2806 ! 2807 CASE DEFAULT 2808 ! 2809 CALL errore('tsvdw','Reference free atom parameters not available for requested atom type...',1) 2810 ! 2811 END SELECT 2812 ! 2813 RETURN 2814 ! 2815 !-------------------------------------------------------------------------------------------------------------- 2816 END SUBROUTINE GetVdWParam 2817 !-------------------------------------------------------------------------------------------------------------- 2818 ! 2819 ! 2820END MODULE tsvdw_module 2821 2822