mirror of https://gitlab.com/QEF/q-e.git
2862 lines
96 KiB
Fortran
2862 lines
96 KiB
Fortran
MODULE tsvdw_module
|
|
!
|
|
!! TS-vdW Code Version 14.0 (RAD/BS, Princeton University, February 2013).
|
|
!
|
|
!! SYNOPSIS:
|
|
!
|
|
!! * radial form of rhoA & drhoA mapped onto linear grid;
|
|
!! * atrho & rhosad on real-space mesh via linear interpolation;
|
|
!! * integration on spherical atomic domains (subsets of real-space mesh);
|
|
!! * quadratic veff derivatives computed linearly using sparse domain intersection algorithm.
|
|
!
|
|
!! All quantities necessary for the evaluation of the TS-vdW energy and forces are computed on the real-space
|
|
!! mesh using linear interpolation of the atomic pseudo-densities and their first derivatives which have been
|
|
!! mapped onto linear equispaced atomic grids from their original form computed on radial atomic grids via the
|
|
!! ATOMIC code.
|
|
!
|
|
USE cell_base, ONLY: h !h matrix for converting between r and s coordinates via r = h s
|
|
USE cell_base, ONLY: ainv !h^-1 matrix for converting between r and s coordinates via s = h^-1 r)
|
|
USE cell_base, ONLY: omega !cell volume (in au^3)
|
|
USE constants, ONLY: pi !pi in double-precision
|
|
USE fft_base, ONLY: dfftp !FFT derived data type
|
|
USE xc_lib, ONLY: xclib_get_id, &
|
|
xclib_dft_is_libxc
|
|
|
|
!correction to correlation utilized in functional
|
|
USE io_global, ONLY: stdout !print/write argument for standard output (to output file)
|
|
USE ions_base, ONLY: nat !number of total atoms (all atomic species)
|
|
USE ions_base, ONLY: nsp !number of unique atomic species
|
|
USE ions_base, ONLY: na !number of atoms within each atomic species
|
|
USE ions_base, ONLY: ityp !ityp(i):=type/species of ith atom
|
|
USE ions_base, ONLY: atm !atm(j):=name of jth atomic species (3 characters)
|
|
USE kinds, ONLY: DP !double-precision kind (selected_real_kind(14,200))
|
|
! the charge density is parallelized over the "band group" or processors
|
|
USE mp_bands, ONLY: nproc_bgrp !number of processors
|
|
USE mp_bands, ONLY: me_bgrp !processor number (0,1,...,nproc_bgrp-1)
|
|
USE mp_bands, ONLY: intra_bgrp_comm !standard MPI communicator
|
|
! atoms are parallelized over the "image group"
|
|
USE mp_images, ONLY: nproc_image !number of processors
|
|
USE mp_images, ONLY: me_image !processor number (0,1,...,nproc_image-1)
|
|
USE mp_images, ONLY: intra_image_comm !standard MPI communicator
|
|
USE mp, ONLY: mp_sum !MPI collection with sum
|
|
USE parallel_include !MPI header
|
|
USE uspp_param, ONLY: upf !atomic pseudo-potential data
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
SAVE
|
|
!
|
|
! PUBLIC variables
|
|
!
|
|
LOGICAL, PUBLIC :: vdw_isolated
|
|
!! isolated system control
|
|
REAL(DP), PUBLIC:: vdw_econv_thr
|
|
!! energy convergence threshold for periodic systems
|
|
REAL(DP), PUBLIC :: EtsvdW
|
|
!! the TS-vdW energy
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PUBLIC :: UtsvdW
|
|
!! the TS-vdW wavefunction forces (dispersion potential)
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: FtsvdW
|
|
!! the TS-vdW ionic forces (-dE/dR)
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: HtsvdW
|
|
!! the TS-vdW cell forces (dE/dh)
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PUBLIC :: VefftsvdW
|
|
!! the TS-vdW effective Hirshfeld volume
|
|
!
|
|
! PRIVATE variables
|
|
!
|
|
INTEGER, PARAMETER, PRIVATE :: NgpA=1000 !number of grid points for linear equispaced atomic grids (current value=1000pts)
|
|
INTEGER, PARAMETER, PRIVATE :: bsint=BIT_SIZE(NgpA) !integer bit size (for use in bit array manipulation)
|
|
INTEGER, PRIVATE :: me !processor number (1,2,...,nproc_image)
|
|
INTEGER, PRIVATE :: iproc !processor dummy index
|
|
INTEGER, PRIVATE :: nr1,nr2,nr3 !real space grid dimensions (global first, second, and third dimensions of the 3D grid)
|
|
INTEGER, PRIVATE :: nr1r,nr2r,nr3r !reduced real space grid dimensions (global first, second, and third dimensions of the 3D grid)
|
|
REAL(DP), PRIVATE :: ddamp !damping function parameter #1
|
|
REAL(DP), PRIVATE :: sR !damping function parameter #2
|
|
REAL(DP), PRIVATE :: spcutAmax !maximum radial cutoff for all atomic species
|
|
INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: nstates !number of atoms per processor
|
|
INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: sdispls !send displacement (offset) array
|
|
INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: rdispls !receive displacement (offset) array
|
|
INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: sendcount !send count array
|
|
INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: recvcount !receive count array
|
|
INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: istatus !MPI status array
|
|
INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: NsomegaA !number of points in the spherical atomic integration domain
|
|
INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: NsomegaAr !number of points in the reduced spherical atomic integration domain
|
|
INTEGER, DIMENSION(:), ALLOCATABLE, PRIVATE :: npair !number of unique atom pairs
|
|
INTEGER, DIMENSION(:,:), ALLOCATABLE, PRIVATE :: pair !unique atom pair overlap matrix
|
|
INTEGER, DIMENSION(:,:), ALLOCATABLE, PRIVATE :: gomegar !precursor to spherical atomic integration domain (intersection bit array)
|
|
INTEGER, DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: somegaA !spherical atomic integration domain
|
|
INTEGER, DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: somegaAr !reduced spherical atomic integration domain
|
|
INTEGER, DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: gomegaAr !reduced spherical atomic integration domain (intersection bit array)
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE:: predveffAdn !atomic dispersion potential prefactor
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: vfree !free atomic volumes for each atomic species
|
|
!GSz
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PUBLIC :: vfree_pub !free atomic volumes for each atomic species
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: dpfree !free atomic static dipole polarizability for each atomic species
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: R0free !free atomic vdW radius for each atomic species
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: C6AAfree !free atomic homonuclear C6 coefficient for each atomic species
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: veff !effective atomic volumes for each atom in the simulation cell
|
|
!GSz
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PUBLIC :: veff_pub !effective atomic volumes for each atom in the simulation cell
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: dpeff !effective atomic static dipole polarizability for each atom in the simulation cell
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: R0eff !effective atomic vdW radius for each atom in the simulation cell
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: C6AAeff !effective atomic homonuclear C6 coefficient for each atom in the simulation cell
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: rhosad !molecular pro-density (superposition of atomic densities) on real-space mesh
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE, PRIVATE :: rhotot !molecular charge density on real-space mesh
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE:: dveffAdn !the local copy of the TS-vdW wavefunction forces (dispersion potential)
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: spgrd !linear equispaced grid for each atomic species
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: sprho !atomic pseudo-density for each atomic species
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: spdrho !first derivative of atomic pseudo-density for each atomic species
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: spdata !linear grid cutoff (is,1) and linear grid spacing (is,2) for each atomic species
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: LIA !A coefficient for linear interpolation of rhoA
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: LIB !B coefficient for linear interpolation of rhoA
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: dLIA !A coefficient for linear interpolation of drhoA
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: dLIB !B coefficient for linear interpolation of drhoA
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: atxyz !Cartesian coordinates of ions adjusted according to PBC
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: C6ABfree !free atomic heteronuclear C6 coefficient for each atom pair
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE, PRIVATE :: C6ABeff !effective atomic heteronuclear C6 coefficient for each atom pair
|
|
REAL(DP), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: dveffdR !first derivative of effective volume wrt nuclear displacement
|
|
REAL(DP), DIMENSION(:,:,:), ALLOCATABLE, PRIVATE :: dveffdh !first derivative of effective volume wrt cell displacement
|
|
!
|
|
! PUBLIC subroutines
|
|
!
|
|
PUBLIC :: tsvdw_initialize
|
|
PUBLIC :: tsvdw_calculate
|
|
PUBLIC :: tsvdw_finalize
|
|
!
|
|
! PRIVATE subroutines
|
|
!
|
|
PRIVATE :: tsvdw_para_init
|
|
PRIVATE :: tsvdw_pbc
|
|
PRIVATE :: tsvdw_unique_pair
|
|
PRIVATE :: tsvdw_rhotot
|
|
PRIVATE :: tsvdw_screen
|
|
PRIVATE :: tsvdw_veff
|
|
PRIVATE :: tsvdw_dveff
|
|
PRIVATE :: tsvdw_effqnts
|
|
PRIVATE :: tsvdw_energy
|
|
PRIVATE :: tsvdw_wfforce
|
|
PRIVATE :: tsvdw_cleanup
|
|
PRIVATE :: Num1stDer
|
|
PRIVATE :: CubSplCoeff
|
|
PRIVATE :: GetVdWParam
|
|
!
|
|
!
|
|
CONTAINS
|
|
!
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_initialize()
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! TS-vdW initialization.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! Local Variables
|
|
!
|
|
LOGICAL :: uniform_grid=.FALSE.
|
|
INTEGER :: ip,iq,ir,is,it,NrgpA,NrgpintA,icutrA,Ndim
|
|
INTEGER :: iexch, icorr, igcx, igcc
|
|
REAL(DP) :: dxA,gfctrA,vref,eref,verr,d,dk1,dk2,dk3,num,den,drab,f1,f2,f3,L1,L2,L3
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE :: atgrdr,atgrdrab,atrhor,datrhor,d2atrhor,CSA,CSB,CSC,CSD
|
|
!
|
|
! Start of calculation banner...
|
|
!
|
|
WRITE(stdout,*)
|
|
WRITE(stdout,'(3X,"TS-vdW initialization")')
|
|
WRITE(stdout,'(3X,"---------------------")')
|
|
WRITE(stdout,*)
|
|
!
|
|
! Error messages for inconsistencies with current version of code...
|
|
!
|
|
!RAD: Have we missed any inconsistencies?
|
|
!
|
|
! Setup variables for use in TS-vdW module...
|
|
!
|
|
nr1=dfftp%nr1; nr2=dfftp%nr2; nr3=dfftp%nr3
|
|
nr1r=nr1/2; nr2r=nr2/2; nr3r=nr3/2
|
|
IF(MOD(nr1,2).EQ.1) nr1r=(nr1+1)/2
|
|
IF(MOD(nr2,2).EQ.1) nr2r=(nr2+1)/2
|
|
IF(MOD(nr3,2).EQ.1) nr3r=(nr3+1)/2
|
|
!
|
|
! Initialize the TS-vdW ionic forces, cell forces, and dispersion potential (wavefunction forces)...
|
|
!
|
|
ALLOCATE(FtsvdW(3,nat)); FtsvdW=0.0_DP
|
|
ALLOCATE(HtsvdW(3,3)); HtsvdW=0.0_DP
|
|
!
|
|
! Initialization of TS-vdW Hirshfeld effective volume public variable ... used in CP print_out.f90
|
|
!
|
|
ALLOCATE(VefftsvdW(nat)); VefftsvdW=0.0_DP
|
|
!
|
|
ALLOCATE(UtsvdW(dfftp%nnr)); UtsvdW=0.0_DP
|
|
!
|
|
! Set ddamp damping function parameter (set to 20 and functional independent)...
|
|
!
|
|
WRITE(stdout,'(3X,"Determining TS-vdW damping function parameters...")')
|
|
ddamp=20.0_DP
|
|
WRITE(stdout,'(5X,"ddamp = ",F9.6)') ddamp
|
|
!
|
|
! Set sR damping function parameter (functional dependent and currently only available for PBE & PBE0)...
|
|
!
|
|
iexch = xclib_get_id('LDA','EXCH')
|
|
icorr = xclib_get_id('LDA','CORR')
|
|
igcx = xclib_get_id('GGA','EXCH')
|
|
igcc = xclib_get_id('GGA','CORR')
|
|
!
|
|
IF ( xclib_dft_is_libxc('LDA','EXCH') .OR. xclib_dft_is_libxc('LDA','CORR') .OR.&
|
|
xclib_dft_is_libxc('GGA','EXCH') .OR. xclib_dft_is_libxc('GGA','CORR') ) &
|
|
CALL errore( 'tsvdw','TS-vdW sR parameter not available for libxc functionals', 1 )
|
|
!
|
|
IF( iexch==1 .AND. icorr==4 .AND. igcx==3 .AND. igcc==4) THEN
|
|
!
|
|
sR=0.94_DP !PBE=sla+pw+pbx+pbc
|
|
!
|
|
ELSEIF( iexch==6 .AND. icorr==4 .AND. igcx==8 .AND. igcc==4) THEN
|
|
!
|
|
sR=0.96_DP !PBE0=pb0x+pw+pb0x+pbc !RAD/BS: This line will not work in CP unless PBE0 code update funct.f90...
|
|
!
|
|
ELSE
|
|
!
|
|
CALL errore( 'tsvdw','TS-vdW sR parameter only available for PBE and PBE0 functionals...', 2 )
|
|
!
|
|
END IF
|
|
!
|
|
WRITE(stdout,'(5X,"sR = ",F9.6)') sR
|
|
!
|
|
! Allocate and initialize species-specific quantities...
|
|
!
|
|
ALLOCATE(vfree(nsp)); vfree=0.0_DP
|
|
IF(.NOT. ALLOCATED(vfree_pub)) ALLOCATE(vfree_pub(nsp)); vfree=0.0_DP
|
|
ALLOCATE(dpfree(nsp)); dpfree=0.0_DP
|
|
ALLOCATE(R0free(nsp)); R0free=0.0_DP
|
|
ALLOCATE(C6AAfree(nsp)); C6AAfree=0.0_DP
|
|
ALLOCATE(C6ABfree(nsp,nsp)); C6ABfree=0.0_DP
|
|
ALLOCATE(spdata(nsp,2)); spdata=0.0_DP
|
|
ALLOCATE(spgrd(nsp,0:NgpA)); spgrd=0.0_DP
|
|
ALLOCATE(sprho(nsp,0:NgpA)); sprho=0.0_DP
|
|
ALLOCATE(spdrho(nsp,0:NgpA)); spdrho=0.0_DP
|
|
ALLOCATE(LIA(nsp,0:NgpA)); LIA=0.0_DP
|
|
ALLOCATE(LIB(nsp,0:NgpA)); LIB=0.0_DP
|
|
ALLOCATE(dLIA(nsp,0:NgpA)); dLIA=0.0_DP
|
|
ALLOCATE(dLIB(nsp,0:NgpA)); dLIB=0.0_DP
|
|
!
|
|
spcutAmax=0.0_DP
|
|
!
|
|
! Loop over atomic species and extract species-dependent quantities to modular arrays...
|
|
!
|
|
DO is=1,nsp
|
|
!
|
|
! Obtain the radial grid and radial atomic pseudo-density from pseudo-potential file (via upf module) for
|
|
! the given atomic species. Convert the radial atomic pseudo-density to the real atomic pseudo-density using
|
|
! rho_real(r) = rho_radial(r) / (4*pi*r^2)...
|
|
!
|
|
WRITE(stdout,'(3X,"Initializing species # ",I3," with atomic symbol ",A3)') is,atm(is)
|
|
!
|
|
! Read in the number of grid points in radial mesh from upf...
|
|
!
|
|
NrgpA=upf(is)%mesh
|
|
!
|
|
! Transfer radial atomic grid (in upf) to local atgrdr array...
|
|
!
|
|
ALLOCATE(atgrdr(NrgpA)); atgrdr=0.0_DP
|
|
!
|
|
DO ir=1,NrgpA
|
|
!
|
|
atgrdr(ir)=upf(is)%r(ir)
|
|
!
|
|
END DO
|
|
!
|
|
! Transfer radial atomic grid spacing (in upf) to local atgrdrab array...
|
|
!
|
|
ALLOCATE(atgrdrab(NrgpA)); atgrdrab=0.0_DP
|
|
!
|
|
DO ir=1,NrgpA
|
|
!
|
|
atgrdrab(ir)=upf(is)%rab(ir)
|
|
!
|
|
END DO
|
|
!
|
|
! Determine whether radial grid is logarithmic/exponential or equispaced/uniform...
|
|
!
|
|
drab=atgrdrab(NrgpA)-atgrdrab(1)
|
|
uniform_grid = (DABS(drab).LT.(1.0E-6_DP))
|
|
IF (uniform_grid) WRITE(stdout,'(5X,"Equispaced/Uniform radial atomic grid detected...")')
|
|
!
|
|
! ----------------------------------------------------------------
|
|
! Logarithmic/Exponential grid (3 parameters: zmesh, xmin, dxA)
|
|
! ----------------------------------------------------------------
|
|
!
|
|
! For i = 1,2,...,NrgpA:
|
|
! r(i) = exp[xmin+(i-1)*dxA]/zmesh
|
|
! = exp[xmin]/zmesh * exp[(i-1)*dxA]
|
|
! = gfctrA * exp[(i-1)*dxA]
|
|
! rab(i) = r(i) * dxA
|
|
!
|
|
! Assumptions: grid does NOT start from zero (use simpson_cp90()).
|
|
!
|
|
! ---------------------------------------------
|
|
! Equispaced/Uniform grid (1 parameter: dxA)
|
|
! ---------------------------------------------
|
|
!
|
|
! For i = 1,2,...,NrgpA:
|
|
! r(i) = (i-1) * dxA
|
|
! rab(i) = dxA
|
|
!
|
|
! Assumptions: grid starts from zero (use simpson() for integration).
|
|
!
|
|
! Determine atomic radial grid parameters...
|
|
!
|
|
IF (uniform_grid.EQV..TRUE.) THEN
|
|
!
|
|
gfctrA=1.0_DP
|
|
dxA=atgrdrab(1)
|
|
!
|
|
ELSE
|
|
!
|
|
gfctrA=upf(is)%r(1)
|
|
dxA=DLOG(upf(is)%r(2)/upf(is)%r(1))
|
|
!
|
|
END IF
|
|
!
|
|
WRITE(stdout,'(5X,"Radial grid parameter: NrgpA is ",I5,".")') NrgpA
|
|
WRITE(stdout,'(5X,"Radial grid parameter: gfctrA is ",F9.6,".")') gfctrA
|
|
WRITE(stdout,'(5X,"Radial grid parameter: dxA is ",F9.6,".")') dxA
|
|
!
|
|
! Transfer radial atomic pseudo-density to atrhor array...
|
|
! Convert radial atomic pseudo-density to real atomic pseudo-density [n(r) = nrad(r)/(4*pi*r^2)]...
|
|
!
|
|
ALLOCATE(atrhor(NrgpA)); atrhor=0.0_DP
|
|
!
|
|
IF (uniform_grid.EQV..TRUE.) THEN
|
|
!
|
|
DO ir=2,NrgpA
|
|
!
|
|
atrhor(ir)=(upf(is)%rho_at(ir))/(4.0_DP*pi*atgrdr(ir)**(2.0_DP)) ! skip point at r=0...
|
|
!
|
|
END DO
|
|
!
|
|
! Quadratic extrapolation of the atomic density to r=0...
|
|
!
|
|
L1=((0.0_DP-atgrdr(3))*(0.0_DP-atgrdr(4)))/((atgrdr(2)-atgrdr(3))*(atgrdr(2)-atgrdr(4)))
|
|
L2=((0.0_DP-atgrdr(2))*(0.0_DP-atgrdr(4)))/((atgrdr(3)-atgrdr(2))*(atgrdr(3)-atgrdr(4)))
|
|
L3=((0.0_DP-atgrdr(2))*(0.0_DP-atgrdr(3)))/((atgrdr(4)-atgrdr(2))*(atgrdr(4)-atgrdr(3)))
|
|
atrhor(1)=L1*atrhor(2)+L2*atrhor(3)+L3*atrhor(4)
|
|
!
|
|
ELSE
|
|
!
|
|
DO ir=1,NrgpA
|
|
!
|
|
atrhor(ir)=(upf(is)%rho_at(ir))/(4.0_DP*pi*atgrdr(ir)**(2.0_DP))
|
|
!
|
|
END DO
|
|
!
|
|
END IF
|
|
!
|
|
! Set NrgpintA as the number of grid points (which must be odd) used during numerical integration using Simpson's rule...
|
|
!
|
|
IF (IAND(NrgpA,1).EQ.1) THEN
|
|
!
|
|
NrgpintA=NrgpA
|
|
!
|
|
ELSE
|
|
!
|
|
NrgpintA=NrgpA-1
|
|
!
|
|
END IF
|
|
!
|
|
! Compute the number of electrons (eref) for each atomic species via numerical integration
|
|
! of the atomic pseudo-density on the radial atomic grid using Simpson's rule...
|
|
!
|
|
eref=0.0_DP
|
|
!
|
|
DO ir=1,NrgpintA-2,2
|
|
!
|
|
f1=atrhor(ir )*atgrdrab(ir )*atgrdr(ir )**(2.0_DP) ! integrated quantity is rho
|
|
f2=atrhor(ir+1)*atgrdrab(ir+1)*atgrdr(ir+1)**(2.0_DP)
|
|
f3=atrhor(ir+2)*atgrdrab(ir+2)*atgrdr(ir+2)**(2.0_DP)
|
|
!
|
|
eref=eref+(f1+4.0_DP*f2+f3)
|
|
!
|
|
END DO
|
|
!
|
|
eref=(4.0_DP*pi/3.0_DP)*eref
|
|
WRITE(stdout,'(5X,"The number of valence electrons, eref, is ",F25.15,".")') eref
|
|
!
|
|
! Compute the reference free atom volume (vref) for each atomic species via numerical integration
|
|
! of the atomic pseudo-density on the radial atomic grid using Simpson's rule...
|
|
!
|
|
vref=0.0_DP
|
|
!
|
|
DO ir=1,NrgpintA-2,2
|
|
!
|
|
f1=atrhor(ir )*atgrdrab(ir )*atgrdr(ir )**(5.0_DP) ! integrated quantity is rho * r^3
|
|
f2=atrhor(ir+1)*atgrdrab(ir+1)*atgrdr(ir+1)**(5.0_DP)
|
|
f3=atrhor(ir+2)*atgrdrab(ir+2)*atgrdr(ir+2)**(5.0_DP)
|
|
!
|
|
vref=vref+(f1+4.0_DP*f2+f3)
|
|
!
|
|
END DO
|
|
!
|
|
vref=(4.0_DP*pi/3.0_DP)*vref
|
|
WRITE(stdout,'(5X,"The reference free atom volume, vref, is ",F25.15," bohr^3.")') vref
|
|
!
|
|
! Using the reference free atom volume, determine an acceptable radial grid cutoff value such that the
|
|
! free atom volume obtained using this cutoff does not deviate from the reference value by more than 1.0%.
|
|
!
|
|
WRITE(stdout,'(5X,"Determining intial radial grid cutoff...")')
|
|
!
|
|
DO iq=5,NrgpintA,2
|
|
!
|
|
vfree(is)=0.0_DP
|
|
verr=0.0_DP
|
|
!
|
|
DO ir=1,iq-2,2
|
|
!
|
|
f1=atrhor(ir )*atgrdrab(ir )*atgrdr(ir )**(5.0_DP) ! integrated quantity is rho * r^3
|
|
f2=atrhor(ir+1)*atgrdrab(ir+1)*atgrdr(ir+1)**(5.0_DP)
|
|
f3=atrhor(ir+2)*atgrdrab(ir+2)*atgrdr(ir+2)**(5.0_DP)
|
|
!
|
|
vfree(is)=vfree(is)+(f1+4.0_DP*f2+f3)
|
|
!
|
|
END DO
|
|
!
|
|
vfree(is)=(4.0_DP*pi/3.0_DP)*vfree(is)
|
|
verr=(vref-vfree(is))/vref*100.0_DP
|
|
!
|
|
IF (verr.LE.1.0_DP) THEN
|
|
!
|
|
icutrA=iq
|
|
!
|
|
WRITE(stdout,'(5X,"An acceptable radial grid cutoff was determined by retaining ",I4," of ",I4," radial grid points.")') &
|
|
icutrA,NrgpA
|
|
!
|
|
EXIT
|
|
!
|
|
END IF
|
|
!
|
|
END DO
|
|
!
|
|
WRITE(stdout,'(5X,"The magnitude of the atomic pseudo-density at the radial grid cutoff is ",ES13.6,".")') atrhor(icutrA)
|
|
WRITE(stdout,'(5X,"Using this radial grid cutoff value of ",F25.15," au:")') atgrdr(icutrA)
|
|
WRITE(stdout,'(5X,"The free atom volume computed with this cutoff is ",F25.15," bohr^3 with an error of ",F6.3,"%.")') &
|
|
vfree(is),verr
|
|
!
|
|
! Form 1st derivative of atrhor for input into cubic spline coefficient subroutine...
|
|
!
|
|
ALLOCATE(datrhor(NrgpA)); datrhor=0.0_DP
|
|
CALL Num1stDer(atgrdr,atrhor,NrgpA,dxA,datrhor)
|
|
!
|
|
! For logarithmic/exponential grid, transform linear derivative back to radial grid...
|
|
!
|
|
IF (.NOT.uniform_grid) THEN
|
|
!
|
|
DO ir=1,NrgpA
|
|
!
|
|
datrhor(ir)=datrhor(ir)/atgrdr(ir)
|
|
!
|
|
END DO
|
|
!
|
|
END IF
|
|
!
|
|
! Form the coefficients of the cubic spline interpolant (2nd derivatives) for the real atomic pseudo-density
|
|
! for use during cubic spline interpolation of the pseudo-density onto the linear equispaced atomic grid...
|
|
!
|
|
ALLOCATE(d2atrhor(NrgpA)); d2atrhor=0.0_DP
|
|
CALL CubSplCoeff(atgrdr,atrhor,NrgpA,datrhor,d2atrhor)
|
|
!
|
|
! Precompute cubic spline interpolation vectors (utilizing Taylor series form) via:
|
|
!
|
|
! y(x) = CSA + CSB*(x-x(k)) + CSC*(x-x(k))**2 + CSD*(x-x(k))**3
|
|
!
|
|
ALLOCATE(CSA(NrgpA)); CSA=0.0_DP
|
|
ALLOCATE(CSB(NrgpA)); CSB=0.0_DP
|
|
ALLOCATE(CSC(NrgpA)); CSC=0.0_DP
|
|
ALLOCATE(CSD(NrgpA)); CSD=0.0_DP
|
|
!
|
|
DO ir=1,NrgpA-1
|
|
!
|
|
! CSA(k) := y(k)
|
|
!
|
|
CSA(ir)=atrhor(ir)
|
|
!
|
|
! CSB(k) := delta(y)/delta(x) - 1/3*delta(x)*y''(k) - 1/6*delta(x)*y''(k+1)
|
|
!
|
|
CSB(ir)=(atrhor(ir+1)-atrhor(ir))/(atgrdr(ir+1)-atgrdr(ir))
|
|
CSB(ir)=CSB(ir)-((1.0_DP/3.0_DP)*(atgrdr(ir+1)-atgrdr(ir))*d2atrhor(ir))
|
|
CSB(ir)=CSB(ir)-((1.0_DP/6.0_DP)*(atgrdr(ir+1)-atgrdr(ir))*d2atrhor(ir+1))
|
|
!
|
|
! CSC(k) := 1/2*y''(k)
|
|
!
|
|
CSC(ir)=(1.0_DP/2.0_DP)*d2atrhor(ir)
|
|
!
|
|
! CSD(k) := 1/6*delta(y'')/delta(x)
|
|
!
|
|
CSD(ir)=((1.0_DP/6.0_DP)*(d2atrhor(ir+1)-d2atrhor(ir))/(atgrdr(ir+1)-atgrdr(ir)))
|
|
!
|
|
END DO
|
|
!
|
|
! Pack species-specific radial cutoff into (is,1) of spdata array...
|
|
!
|
|
spdata(is,1)=atgrdr(icutrA)
|
|
IF (spdata(is,1).GT.spcutAmax) spcutAmax=spdata(is,1)
|
|
!
|
|
! Compute and pack grid spacing of species-specific linear equispaced grid into (is,2) of spdata array...
|
|
!
|
|
spdata(is,2)=(atgrdr(icutrA)+1.0_DP)/DBLE(NgpA) !include additional buffer of 1 bohr...
|
|
WRITE(stdout,'(5X,"Linear grid spacing was computed as: ",F25.15," bohr.")') spdata(is,2)
|
|
!
|
|
! Form linear equispaced atomic grid (NOT including point at r=0) and pack into argument (is,:) of spgrd array...
|
|
!
|
|
DO ip=1,NgpA
|
|
!
|
|
spgrd(is,ip)=DBLE(ip)*spdata(is,2)
|
|
!
|
|
END DO
|
|
!
|
|
! Map atomic pseudo-density (currently on the radial atomic grid) onto linear equispaced atomic grid using
|
|
! cubic spline interpolation...Form first derivative of the atomic pseudo-density on the linear equispaced
|
|
! atomic grid via differentiation of the cubic spline interpolant...
|
|
!
|
|
DO ip=1,NgpA
|
|
!
|
|
d=spgrd(is,ip)
|
|
!
|
|
IF (uniform_grid.EQV..TRUE.) THEN
|
|
!
|
|
ir=INT(d/dxA)+1 !since the equispaced/uniform grid first point is at r=0...
|
|
!
|
|
ELSE
|
|
!
|
|
ir=FLOOR(DLOG(d*EXP(dxA)/gfctrA)/dxA)
|
|
!
|
|
END IF
|
|
!
|
|
dk1=d-atgrdr(ir); dk2=dk1*dk1; dk3=dk2*dk1
|
|
sprho(is,ip)=CSA(ir)+CSB(ir)*dk1+CSC(ir)*dk2+CSD(ir)*dk3 !Pack density into argument (is,:) of sprho array
|
|
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
|
|
!
|
|
END DO
|
|
!
|
|
! For computational efficiency during the remainder of the calculation, extrapolate sprho and spdrho to
|
|
! include the point at r=0 (this eliminates an if statement in crucial inner loops)...
|
|
! Use quadratic extrapolation to obtain these points...Hence the 0:NgpA dimension above...
|
|
!
|
|
spgrd(is,0)=0.0_DP !Extend linear grid to include point at r=0...
|
|
!
|
|
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)))
|
|
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)))
|
|
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)))
|
|
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...
|
|
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...
|
|
!
|
|
! Throughout the remainder of the code, to map the atomic quantities onto the real-space mesh, we will be
|
|
! utilizing the Taylor series form of linear interpolation, given by:
|
|
!
|
|
! y(x) = LIA + LIB*(x-x(k)) y'(x) = dLIA + dLIB*(x-x(k))
|
|
!
|
|
! for x(k) <= x <= x(k+1)...
|
|
!
|
|
DO ip=0,NgpA-1
|
|
!
|
|
! LIA(k) := y(k)
|
|
!
|
|
LIA(is,ip)=sprho(is,ip)
|
|
dLIA(is,ip)=spdrho(is,ip)
|
|
!
|
|
! LIB(k) := delta(y)/delta(x)
|
|
!
|
|
LIB(is,ip)=(sprho(is,ip+1)-sprho(is,ip))/(spgrd(is,ip+1)-spgrd(is,ip))
|
|
dLIB(is,ip)=(spdrho(is,ip+1)-spdrho(is,ip))/(spgrd(is,ip+1)-spgrd(is,ip))
|
|
!
|
|
END DO
|
|
!
|
|
! Populate reference free atom quantities...
|
|
!
|
|
CALL GetVdWParam(upf(is)%psd,C6AAfree(is),dpfree(is),R0free(is))
|
|
!
|
|
WRITE(stdout,'(5X,"The free atom static dipole polarizability is ",F13.6," bohr^3.")') dpfree(is)
|
|
WRITE(stdout,'(5X,"The free atom homonuclear C6 coefficient is ",F13.6," Hartree bohr^6.")') C6AAfree(is)
|
|
WRITE(stdout,'(5X,"The free atom vdW radius is ",F13.6," bohr.")') R0free(is)
|
|
!
|
|
! Clean-up all species-specific temporary arrays
|
|
!
|
|
IF (ALLOCATED(atgrdr)) DEALLOCATE(atgrdr)
|
|
IF (ALLOCATED(atgrdrab)) DEALLOCATE(atgrdrab)
|
|
IF (ALLOCATED(atrhor)) DEALLOCATE(atrhor)
|
|
IF (ALLOCATED(datrhor)) DEALLOCATE(datrhor)
|
|
IF (ALLOCATED(d2atrhor)) DEALLOCATE(d2atrhor)
|
|
IF (ALLOCATED(CSA)) DEALLOCATE(CSA)
|
|
IF (ALLOCATED(CSB)) DEALLOCATE(CSB)
|
|
IF (ALLOCATED(CSC)) DEALLOCATE(CSC)
|
|
IF (ALLOCATED(CSD)) DEALLOCATE(CSD)
|
|
!
|
|
END DO !is
|
|
!
|
|
! Compute free heteronuclear C6 coefficient matrix...
|
|
! C6ABfree(A,B)=[2*C6AAfree(A)*C6AAfree(B)]/[(dpfree(B)/dpfree(A))*C6AAfree(A)+(dpfree(A)/dpfree(B))*C6AAfree(B)]
|
|
!
|
|
DO is=1,nsp
|
|
!
|
|
DO it=1,nsp
|
|
!
|
|
num=2.0_DP*C6AAfree(is)*C6AAfree(it)
|
|
den=(dpfree(it)/dpfree(is))*C6AAfree(is)+(dpfree(is)/dpfree(it))*C6AAfree(it)
|
|
C6ABfree(is,it)=num/den
|
|
!
|
|
END DO
|
|
!
|
|
END DO
|
|
vfree_pub=vfree
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_initialize
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_calculate(tauin, rhor)
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Manages entire calculation of TS-vdW energy, wavefunction forces, and ion forces.
|
|
!
|
|
! TS-vdW Management Code: Manages entire calculation of TS-vdW energy, wavefunction forces, and ion forces via
|
|
! calls to PRIVATE subroutines below (called in each MD step). The calls to tsvdw_initialize and tsvdw_finalize
|
|
! are done once at the beginning (init_run) and the end (terminate_run).
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! I/O variables
|
|
!
|
|
REAL(DP), INTENT(IN) :: rhor(:)
|
|
REAL(DP) :: tauin(3,nat)
|
|
!
|
|
! Parallel initialization...
|
|
!
|
|
CALL tsvdw_para_init()
|
|
!
|
|
! Move all atoms into simulation cell by adjusting Cartesian coordinates according to PBCs...
|
|
!
|
|
CALL tsvdw_pbc(tauin)
|
|
!
|
|
! Compute unique atom pair list...
|
|
!
|
|
CALL tsvdw_unique_pair()
|
|
!
|
|
! Obtain molecular charge density given on the real-space mesh...
|
|
!
|
|
CALL tsvdw_rhotot( rhor )
|
|
!
|
|
! Determine spherical atomic integration domains and atom overlap (bit array)...
|
|
! Compute molecular pro-density (superposition of atomic densities) on the real-space mesh...
|
|
! Compute functional derivative of vdW energy wrt charge density (numerator only)...
|
|
!
|
|
CALL tsvdw_screen()
|
|
!
|
|
! Compute effective volume for each atom in the simulation cell...
|
|
! Complete functional derivative of vdW energy wrt charge density...
|
|
!
|
|
CALL tsvdw_veff()
|
|
!
|
|
! Calculate first derivative of veff wrt nuclear and cell displacements...
|
|
!
|
|
CALL tsvdw_dveff()
|
|
!
|
|
! Calculate effective quantities for each atom in the simulation cell...
|
|
!
|
|
CALL tsvdw_effqnts()
|
|
!
|
|
! Calculate total TS-vdW energy, dispersion potential prefactor, ionic forces, and cell forces...
|
|
!
|
|
CALL tsvdw_energy()
|
|
!
|
|
! Calculate total TS-vdW wavefunction forces (dispersion potential)...
|
|
!
|
|
CALL tsvdw_wfforce()
|
|
!
|
|
! Deallocate all arrays specific to tsvdw_calculate...
|
|
!
|
|
CALL tsvdw_cleanup()
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_calculate
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_para_init()
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! TS-vdW, parallel initialization
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
INTEGER :: i,j,k
|
|
!
|
|
me=me_image+1
|
|
!
|
|
ALLOCATE(nstates(nproc_image)); nstates=0
|
|
ALLOCATE(sdispls(nproc_image)); sdispls=0
|
|
ALLOCATE(sendcount(nproc_image)); sendcount=0
|
|
ALLOCATE(rdispls(nproc_image)); rdispls=0
|
|
ALLOCATE(recvcount(nproc_image)); recvcount=0
|
|
ALLOCATE(istatus(nproc_image)); istatus=0
|
|
!
|
|
! Assign workload of atoms over nproc_image processors
|
|
!
|
|
IF (nat.LE.nproc_image) THEN
|
|
!
|
|
DO i=1,nat
|
|
!
|
|
nstates(i)=1
|
|
!
|
|
END DO
|
|
!
|
|
ELSE
|
|
!
|
|
k=0
|
|
!
|
|
10 DO j=1,nproc_image
|
|
!
|
|
nstates(j)=nstates(j)+1
|
|
!
|
|
k=k+1
|
|
!
|
|
IF (k.GE.nat) GO TO 20
|
|
!
|
|
END DO
|
|
!
|
|
IF (k.LT.nat) GO TO 10
|
|
!
|
|
END IF
|
|
!
|
|
20 CONTINUE
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_para_init
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_pbc(tauin)
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Move all atoms into simulation cell by adjusting Cartesian coordinates according to PBCs.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! I/O variables
|
|
!
|
|
REAL(DP) :: tauin(3,nat)
|
|
!
|
|
! Local variables
|
|
!
|
|
INTEGER :: ia
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE :: atxyzs
|
|
!
|
|
! Initialization of PBC-adjusted Cartesian coordinates...
|
|
!
|
|
ALLOCATE(atxyz(3,nat)); atxyz=0.0_DP
|
|
ALLOCATE(atxyzs(3,nat)); atxyzs=0.0_DP
|
|
!
|
|
! Adjust Cartesian coordinates of ions according to periodic boundary conditions...
|
|
! N.B.: PBC are imposed here in the range [0,1)...
|
|
!
|
|
DO ia = 1, nat
|
|
!
|
|
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
|
|
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
|
|
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
|
|
!
|
|
atxyzs(1,ia)=atxyzs(1,ia)-FLOOR(atxyzs(1,ia)) ! impose PBC on s in range: [0,1)
|
|
atxyzs(2,ia)=atxyzs(2,ia)-FLOOR(atxyzs(2,ia)) ! impose PBC on s in range: [0,1)
|
|
atxyzs(3,ia)=atxyzs(3,ia)-FLOOR(atxyzs(3,ia)) ! impose PBC on s in range: [0,1)
|
|
!
|
|
atxyz(1,ia)=h(1,1)*atxyzs(1,ia)+h(1,2)*atxyzs(2,ia)+h(1,3)*atxyzs(3,ia) ! r = h s
|
|
atxyz(2,ia)=h(2,1)*atxyzs(1,ia)+h(2,2)*atxyzs(2,ia)+h(2,3)*atxyzs(3,ia) ! r = h s
|
|
atxyz(3,ia)=h(3,1)*atxyzs(1,ia)+h(3,2)*atxyzs(2,ia)+h(3,3)*atxyzs(3,ia) ! r = h s
|
|
!
|
|
END DO
|
|
!
|
|
IF (ALLOCATED(atxyzs)) DEALLOCATE(atxyzs)
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_pbc
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_unique_pair()
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Compute unique atom pair list.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! Local variables
|
|
!
|
|
INTEGER :: ia,ib,ias,ibs,ip,ir,i,j,k,jj,nj_max,nbmax,num,num1,jj_neib_of_i
|
|
REAL(DP) :: spcutA,spcutB,dAB(3),dAB2(3)
|
|
INTEGER, DIMENSION(:), ALLOCATABLE :: nj,overlap2
|
|
INTEGER, DIMENSION(:,:), ALLOCATABLE :: overlap
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE :: dABmic
|
|
!
|
|
CALL start_clock('tsvdw_pair')
|
|
!
|
|
! Allocate and initialize temporary arrays...
|
|
!
|
|
ALLOCATE(dABmic(nat)); dABmic=0.0_DP
|
|
ALLOCATE(overlap(nat,nat)); overlap=0
|
|
ALLOCATE(overlap2(nat)); overlap2=0
|
|
ALLOCATE(nj(nat)); nj=0
|
|
!
|
|
! Outer loop over atoms A to form non-unique atom pair overlap matrix...
|
|
!
|
|
DO ia=1,nat
|
|
!
|
|
nj(ia)=0; dABmic=0.0_DP
|
|
!
|
|
! Connect atom type with species-dependent quantities...
|
|
!
|
|
ias=ityp(ia)
|
|
!
|
|
! Transfer species-specific cutoff to spcutA...
|
|
!
|
|
spcutA=spdata(ias,1)
|
|
!
|
|
! Inner loop over atoms B...
|
|
!
|
|
DO ib=1,nat
|
|
!
|
|
IF(ib.NE.ia) THEN
|
|
!
|
|
! Connect atom type with species-dependent quantities...
|
|
!
|
|
ibs=ityp(ib)
|
|
!
|
|
! Transfer species-specific cutoff to spcutB...
|
|
!
|
|
spcutB=spdata(ibs,1)
|
|
!
|
|
! Compute distance between atom A and atom B (according to the minimum image convention)...
|
|
!
|
|
dAB(1)=atxyz(1,ia)-atxyz(1,ib) ! r_AB = r_A - r_B
|
|
dAB(2)=atxyz(2,ia)-atxyz(2,ib) ! r_AB = r_A - r_B
|
|
dAB(3)=atxyz(3,ia)-atxyz(3,ib) ! r_AB = r_A - r_B
|
|
!
|
|
dAB2(1)=ainv(1,1)*dAB(1)+ainv(1,2)*dAB(2)+ainv(1,3)*dAB(3) ! s_AB = h^-1 r_AB
|
|
dAB2(2)=ainv(2,1)*dAB(1)+ainv(2,2)*dAB(2)+ainv(2,3)*dAB(3) ! s_AB = h^-1 r_AB
|
|
dAB2(3)=ainv(3,1)*dAB(1)+ainv(3,2)*dAB(2)+ainv(3,3)*dAB(3) ! s_AB = h^-1 r_AB
|
|
!
|
|
dAB2(1)=dAB2(1)-IDNINT(dAB2(1)) ! impose MIC on s_AB in range: [-0.5,+0.5]
|
|
dAB2(2)=dAB2(2)-IDNINT(dAB2(2)) ! impose MIC on s_AB in range: [-0.5,+0.5]
|
|
dAB2(3)=dAB2(3)-IDNINT(dAB2(3)) ! impose MIC on s_AB in range: [-0.5,+0.5]
|
|
!
|
|
dAB(1)=h(1,1)*dAB2(1)+h(1,2)*dAB2(2)+h(1,3)*dAB2(3) ! r_AB = h s_AB (MIC)
|
|
dAB(2)=h(2,1)*dAB2(1)+h(2,2)*dAB2(2)+h(2,3)*dAB2(3) ! r_AB = h s_AB (MIC)
|
|
dAB(3)=h(3,1)*dAB2(1)+h(3,2)*dAB2(2)+h(3,3)*dAB2(3) ! r_AB = h s_AB (MIC)
|
|
!
|
|
dABmic(ib)=DSQRT(dAB(1)*dAB(1)+dAB(2)*dAB(2)+dAB(3)*dAB(3)) ! |r_A - r_B| (MIC)
|
|
!
|
|
IF(dABmic(ib).LT.(spcutA+spcutB)) THEN
|
|
!
|
|
nj(ia)=nj(ia)+1
|
|
overlap(nj(ia),ia)=ib
|
|
!
|
|
IF(nj(ia).EQ.1) THEN
|
|
!
|
|
overlap(nj(ia),ia)=ib
|
|
!
|
|
ELSE IF(dABmic(overlap(nj(ia)-1,ia)).LE.dABmic(ib)) THEN
|
|
!
|
|
overlap(nj(ia),ia)=ib
|
|
!
|
|
ELSE
|
|
!
|
|
overlap2(:)=0
|
|
!
|
|
DO ir=1,nj(ia)-1
|
|
!
|
|
IF(dABmic(overlap(ir,ia)).LT.dABmic(ib)) THEN
|
|
!
|
|
overlap2(ir)=overlap(ir,ia)
|
|
!
|
|
ELSE
|
|
!
|
|
overlap2(ir)=ib
|
|
!
|
|
DO ip=ir+1,nj(ia)
|
|
!
|
|
overlap2(ip)=overlap(ip-1,ia)
|
|
!
|
|
END DO
|
|
!
|
|
GO TO 30
|
|
!
|
|
END IF
|
|
!
|
|
END DO !ir
|
|
!
|
|
30 CONTINUE
|
|
!
|
|
DO ir=1,nj(ia)
|
|
!
|
|
overlap(ir,ia)=overlap2(ir)
|
|
!
|
|
END DO
|
|
!
|
|
END IF !nj(ia)
|
|
!
|
|
END IF !dABmic(j)
|
|
!
|
|
END IF !ia/=ib
|
|
!
|
|
END DO !ib
|
|
!
|
|
END DO !ia
|
|
!
|
|
IF (ALLOCATED(dABmic)) DEALLOCATE(dABmic)
|
|
IF (ALLOCATED(overlap2)) DEALLOCATE(overlap2)
|
|
!
|
|
! Now form unique atom pair overlap matrix...
|
|
!
|
|
nbmax=nat
|
|
!
|
|
ALLOCATE(pair(nbmax,nat)); pair=0
|
|
ALLOCATE(npair(nat)); npair=0
|
|
!
|
|
num=0; num1=0
|
|
!
|
|
DO j=1,nbmax
|
|
!
|
|
DO ia=1,nat
|
|
!
|
|
DO jj=1,nj(ia)
|
|
!
|
|
jj_neib_of_i=overlap(jj,ia)
|
|
!
|
|
IF(jj_neib_of_i.GT.0) THEN
|
|
!
|
|
pair(j,ia)=jj_neib_of_i
|
|
overlap(jj,ia)=0
|
|
num=num+1
|
|
!
|
|
DO k=1,nj(jj_neib_of_i)
|
|
!
|
|
IF(overlap(k,jj_neib_of_i).EQ.ia) THEN
|
|
!
|
|
overlap(k,jj_neib_of_i)=0
|
|
num1=num1+1
|
|
!
|
|
GO TO 40
|
|
!
|
|
END IF
|
|
!
|
|
END DO !k
|
|
!
|
|
END IF
|
|
!
|
|
END DO !jj
|
|
!
|
|
40 CONTINUE
|
|
!
|
|
END DO !ia
|
|
!
|
|
END DO !j
|
|
!
|
|
IF(num.NE.num1) THEN
|
|
!
|
|
CALL errore('tsvdw','ERROR: num .NE. num1...',1)
|
|
!
|
|
END IF
|
|
!
|
|
! Count number of unique atom pairs for each atom...
|
|
!
|
|
DO ia=1,nat
|
|
!
|
|
num=0
|
|
!
|
|
DO j=1,nbmax
|
|
!
|
|
IF(pair(j,ia).NE.0) num=num+1
|
|
!
|
|
END DO
|
|
!
|
|
npair(ia)=num
|
|
!
|
|
END DO
|
|
!
|
|
IF (ALLOCATED(overlap)) DEALLOCATE(overlap)
|
|
IF (ALLOCATED(nj)) DEALLOCATE(nj)
|
|
!
|
|
CALL stop_clock('tsvdw_pair')
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_unique_pair
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_rhotot( rhor )
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Obtain molecular charge density given on the real-space mesh.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
REAL(DP), INTENT(IN) :: rhor(:)
|
|
!
|
|
! Local variables
|
|
!
|
|
INTEGER :: ir,ierr
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE :: rhor_tmp
|
|
!
|
|
CALL start_clock('tsvdw_rhotot')
|
|
!
|
|
! Initialization of rhotot array (local copy of the real-space charge density)...
|
|
!
|
|
ALLOCATE(rhotot(nr1*nr2*nr3)); rhotot=0.0_DP
|
|
#if defined(__MPI)
|
|
!
|
|
! Initialization of rhor_tmp temporary buffers...
|
|
!
|
|
ALLOCATE(rhor_tmp(nr1*nr2*nr3)); rhor_tmp=0.0_DP
|
|
!
|
|
! Collect distributed rhor and broadcast to all processors...
|
|
!
|
|
DO iproc=1,nproc_bgrp
|
|
!
|
|
recvcount(iproc)=dfftp%nr3p(iproc)*nr1*nr2
|
|
!
|
|
END DO
|
|
!
|
|
rdispls(1) = 0
|
|
!
|
|
DO iproc=2,nproc_bgrp
|
|
!
|
|
rdispls(iproc)=rdispls(iproc-1)+recvcount(iproc-1)
|
|
!
|
|
END DO
|
|
!
|
|
CALL MPI_ALLGATHERV(rhor(1),dfftp%nr3p(me_bgrp+1)*nr1*nr2,&
|
|
MPI_DOUBLE_PRECISION,rhor_tmp(1),recvcount,rdispls,&
|
|
MPI_DOUBLE_PRECISION,intra_bgrp_comm,ierr)
|
|
!
|
|
! Transfer rhor temporary arrays to rhotot array...
|
|
!
|
|
rhotot=rhor_tmp
|
|
!
|
|
! Clean-up temporary arrays...
|
|
!
|
|
IF (ALLOCATED(rhor_tmp)) DEALLOCATE(rhor_tmp)
|
|
!
|
|
#else
|
|
rhotot(:) = rhor(:)
|
|
#endif
|
|
CALL stop_clock('tsvdw_rhotot')
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_rhotot
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_screen()
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Determine spherical atomic integration domains and atom overlap (bit array)...
|
|
!! Compute molecular pro-density (superposition of atomic densities) on the real-space mesh...
|
|
!! Compute functional derivative of vdW energy wrt charge density (numerator only).
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! Local variables
|
|
!
|
|
INTEGER :: ia,ias,Ntmp,Ntmpr,Npts,Nptsr,ir1,ir2,ir3,Ndim,Nints,off1,off1r,ioff,boff,iq,ir
|
|
REAL(DP) :: spcutA,spdA,dq(3),dqA(3),dqAs(3),dk1,rhoA
|
|
REAL(DP), ALLOCATABLE :: dqAmic(:,:,:),dveffAdntmp(:,:,:)
|
|
!
|
|
CALL start_clock('tsvdw_screen')
|
|
!
|
|
! Allocate and initialize gomegar array which contains (in a bit array) which atoms contribute to a given point
|
|
! on the reduced real-space grid for all atoms...
|
|
!
|
|
Nints=nat/bsint+1
|
|
ALLOCATE(gomegar(nr1r*nr2r*nr3r,Nints)); gomegar=0
|
|
!
|
|
! Allocate and initialize NsomegaA and NsomegaAr arrays which contains the number of points in the
|
|
! full and reduced spherical atomic integration domains for all atoms...
|
|
!
|
|
ALLOCATE(NsomegaA(nat)); NsomegaA=0
|
|
ALLOCATE(NsomegaAr(nat)); NsomegaAr=0
|
|
!
|
|
! Allocate and initialize somegaA and somegaAr arrays which contains the grid indices (along nr1,nr2,nr3) for each point in the
|
|
! full and reduced spherical atomic integration domains for each atom assigned to a given processor...
|
|
!
|
|
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)
|
|
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)
|
|
Ndim=MAX(1,nstates(me))
|
|
ALLOCATE(somegaA(Ntmp,3,Ndim)); somegaA=0
|
|
ALLOCATE(somegaAr(Ntmpr,3,Ndim)); somegaAr=0
|
|
!
|
|
! Allocate and initialize gomegaAr array which contains in a bit array all of the atoms that intersect with each point in the
|
|
! reduced spherical atomic integration domain for each atom assigned to a given processor...
|
|
!
|
|
ALLOCATE(gomegaAr(Ntmpr,Nints,Ndim)); gomegaAr=0
|
|
!
|
|
! Initialization of rhosad(r)...
|
|
!
|
|
ALLOCATE(rhosad(nr1*nr2*nr3)); rhosad=0.0_DP
|
|
!
|
|
! Initialization of dVA/dn(r)...
|
|
!
|
|
ALLOCATE(dveffAdn(Ntmp,Ndim)); dveffAdn=0.0_DP
|
|
!
|
|
DO iproc=1,nstates(me)
|
|
!
|
|
! Connect processor number with atom...
|
|
!
|
|
ia=me+nproc_image*(iproc-1)
|
|
!
|
|
! Connect atom type with species-dependent quantities...
|
|
!
|
|
ias=ityp(ia)
|
|
!
|
|
! Transfer species-specific cutoff to spcutA...
|
|
!
|
|
spcutA=spdata(ias,1)
|
|
!
|
|
! Precompute inverse of species-specific linear grid spacing (replaces / with * inside inner loop)...
|
|
!
|
|
spdA=1.0_DP/spdata(ias,2)
|
|
!
|
|
! Loop over grid points and determine if they belong to spherical atomic integration domain (if r < RcutA)...
|
|
!
|
|
Npts=0; Nptsr=0
|
|
!
|
|
ALLOCATE(dqAmic(nr1,nr2,nr3)); dqAmic=0.0_DP
|
|
ALLOCATE(dveffAdntmp(nr1,nr2,nr3)); dveffAdntmp=0.0_DP
|
|
!
|
|
!$omp parallel do private(dq,dqA,dqAs,ir,dk1,rhoA,off1,ioff,boff,off1r)
|
|
DO ir1=1,nr1
|
|
!
|
|
dq(1)=DBLE(ir1-1)/DBLE(nr1) ! s_i(1)
|
|
!
|
|
DO ir2=1,nr2
|
|
!
|
|
dq(2)=DBLE(ir2-1)/DBLE(nr2) ! s_i(2)
|
|
!
|
|
DO ir3=1,nr3
|
|
!
|
|
dq(3)=DBLE(ir3-1)/DBLE(nr3) ! s_i(3)
|
|
!
|
|
! Compute distance between grid point and atom according to minimum image convention (MIC)...
|
|
!
|
|
dqA(1)=h(1,1)*dq(1)+h(1,2)*dq(2)+h(1,3)*dq(3) ! r_i = h s_i
|
|
dqA(2)=h(2,1)*dq(1)+h(2,2)*dq(2)+h(2,3)*dq(3) ! r_i = h s_i
|
|
dqA(3)=h(3,1)*dq(1)+h(3,2)*dq(2)+h(3,3)*dq(3) ! r_i = h s_i
|
|
!
|
|
dqA(1)=dqA(1)-atxyz(1,ia) ! r_iA = r_i - r_A
|
|
dqA(2)=dqA(2)-atxyz(2,ia) ! r_iA = r_i - r_A
|
|
dqA(3)=dqA(3)-atxyz(3,ia) ! r_iA = r_i - r_A
|
|
!
|
|
dqAs(1)=ainv(1,1)*dqA(1)+ainv(1,2)*dqA(2)+ainv(1,3)*dqA(3) ! s_iA = h^-1 r_iA
|
|
dqAs(2)=ainv(2,1)*dqA(1)+ainv(2,2)*dqA(2)+ainv(2,3)*dqA(3) ! s_iA = h^-1 r_iA
|
|
dqAs(3)=ainv(3,1)*dqA(1)+ainv(3,2)*dqA(2)+ainv(3,3)*dqA(3) ! s_iA = h^-1 r_iA
|
|
!
|
|
dqAs(1)=dqAs(1)-IDNINT(dqAs(1)) ! impose MIC on s_iA in range: [-0.5,+0.5]
|
|
dqAs(2)=dqAs(2)-IDNINT(dqAs(2)) ! impose MIC on s_iA in range: [-0.5,+0.5]
|
|
dqAs(3)=dqAs(3)-IDNINT(dqAs(3)) ! impose MIC on s_iA in range: [-0.5,+0.5]
|
|
!
|
|
dqA(1)=h(1,1)*dqAs(1)+h(1,2)*dqAs(2)+h(1,3)*dqAs(3) ! r_iA = h s_iA (MIC)
|
|
dqA(2)=h(2,1)*dqAs(1)+h(2,2)*dqAs(2)+h(2,3)*dqAs(3) ! r_iA = h s_iA (MIC)
|
|
dqA(3)=h(3,1)*dqAs(1)+h(3,2)*dqAs(2)+h(3,3)*dqAs(3) ! r_iA = h s_iA (MIC)
|
|
!
|
|
dqAmic(ir1,ir2,ir3)=DSQRT(dqA(1)*dqA(1)+dqA(2)*dqA(2)+dqA(3)*dqA(3)) ! |r_i - r_A| (MIC)
|
|
!
|
|
! Screen grid point according to atomic radial cutoff...
|
|
!
|
|
IF (dqAmic(ir1,ir2,ir3).LE.spcutA) THEN
|
|
!
|
|
! Form rhosad(r) on the real-space mesh...
|
|
! 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...
|
|
!
|
|
! Determine the index in the atomic linear equispaced grid such that grd(ir) <= dqA <= grd(ir+1) and distance between dqA and grd(ir)...
|
|
!
|
|
ir=INT(dqAmic(ir1,ir2,ir3)*spdA)
|
|
dk1=dqAmic(ir1,ir2,ir3)-spgrd(ias,ir)
|
|
!
|
|
! Perform linear interpolation to obtain the value of the atomic pseudo-density at the given grid point...
|
|
!
|
|
rhoA=LIA(ias,ir)+LIB(ias,ir)*dk1
|
|
!
|
|
! Increment contribution to rhosad(r)...
|
|
!
|
|
off1=ir1+(ir2-1)*nr1+(ir3-1)*nr1*nr2 !global offset [nr1,nr2,nr3]
|
|
rhosad(off1)=rhosad(off1)+rhoA
|
|
!
|
|
! Form numerator of dVA/dn(r) only...
|
|
!
|
|
dveffAdntmp(ir1,ir2,ir3)=dqAmic(ir1,ir2,ir3)**(3.0_DP)*rhoA
|
|
!
|
|
! On reduced grid only, form screened somegaAr and gomegar...
|
|
!
|
|
IF ((MOD(ir1,2).EQ.1).AND.(MOD(ir2,2).EQ.1).AND.(MOD(ir3,2).EQ.1)) THEN
|
|
!
|
|
ioff=((ia-1)/bsint)+1 ! integer offset for gomegar bit array
|
|
boff=(ia-((ioff-1)*bsint))-1 ! bit offset for gomegar bit array
|
|
off1r=(ir1+1)/2+((ir2-1)/2)*nr1r+((ir3-1)/2)*nr1r*nr2r ! reduced global offset [nr1r,nr2r,nr3r]
|
|
!
|
|
gomegar(off1r,ioff)=IBSET(gomegar(off1r,ioff),boff)
|
|
!
|
|
END IF
|
|
!
|
|
END IF
|
|
!
|
|
END DO !ir3
|
|
!
|
|
END DO !ir2
|
|
!
|
|
END DO !ir1
|
|
!$omp end parallel do
|
|
!
|
|
DO ir1=1,nr1
|
|
!
|
|
DO ir2=1,nr2
|
|
!
|
|
DO ir3=1,nr3
|
|
!
|
|
! Screen grid point according to atomic radial cutoff...
|
|
!
|
|
IF (dqAmic(ir1,ir2,ir3).LE.spcutA) THEN
|
|
!
|
|
Npts=Npts+1
|
|
!
|
|
! Form screened somegaA...
|
|
!
|
|
somegaA(Npts,1,iproc)=ir1
|
|
somegaA(Npts,2,iproc)=ir2
|
|
somegaA(Npts,3,iproc)=ir3
|
|
!
|
|
dveffAdn(Npts,iproc)=dveffAdntmp(ir1,ir2,ir3)
|
|
!
|
|
! On reduced grid only, form screened somegaAr ...
|
|
!
|
|
IF ((MOD(ir1,2).EQ.1).AND.(MOD(ir2,2).EQ.1).AND.(MOD(ir3,2).EQ.1)) THEN
|
|
!
|
|
Nptsr=Nptsr+1
|
|
!
|
|
! Form reduced screened somegaAr...
|
|
!
|
|
somegaAr(Nptsr,1,iproc)=ir1
|
|
somegaAr(Nptsr,2,iproc)=ir2
|
|
somegaAr(Nptsr,3,iproc)=ir3
|
|
!
|
|
END IF
|
|
!
|
|
END IF
|
|
|
|
END DO !ir3
|
|
!
|
|
END DO !ir2
|
|
!
|
|
END DO !ir1
|
|
!
|
|
NsomegaA(ia)=Npts
|
|
NsomegaAr(ia)=Nptsr
|
|
!
|
|
IF (ALLOCATED(dqAmic)) DEALLOCATE(dqAmic)
|
|
IF (ALLOCATED(dveffAdntmp)) DEALLOCATE(dveffAdntmp)
|
|
!
|
|
END DO ! iproc
|
|
!
|
|
! Collect NsomegaA, NsomegaAr, gomegar, and rhosad over all processors and broadcast...
|
|
!
|
|
CALL mp_sum(NsomegaA,intra_image_comm)
|
|
CALL mp_sum(NsomegaAr,intra_image_comm)
|
|
CALL mp_sum(gomegar,intra_image_comm)
|
|
CALL mp_sum(rhosad,intra_image_comm)
|
|
!
|
|
! Decompose gomegar to gomegaAr to save on memory storage...
|
|
!
|
|
DO iproc=1,nstates(me)
|
|
!
|
|
! Connect processor number with atom...
|
|
!
|
|
ia=me+nproc_image*(iproc-1)
|
|
!
|
|
! Loop over points in the (pre-screened) reduced spherical atomic integration domain...
|
|
!
|
|
!$omp parallel do private(off1r,ir)
|
|
DO iq=1,NsomegaAr(ia)
|
|
!
|
|
DO ir=1,Nints
|
|
!
|
|
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]
|
|
gomegaAr(iq,ir,iproc)=gomegar(off1r,ir)
|
|
!
|
|
END DO
|
|
!
|
|
END DO !iq
|
|
!$omp end parallel do
|
|
!
|
|
END DO ! iproc
|
|
!
|
|
! Clean-up temporary arrays...
|
|
!
|
|
IF (ALLOCATED(gomegar)) DEALLOCATE(gomegar)
|
|
!
|
|
CALL stop_clock('tsvdw_screen')
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_screen
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_veff()
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Compute effective volume for each atom in the simulation cell.
|
|
!! Complete functional derivative of vdW energy wrt charge density.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! Local variables
|
|
!
|
|
INTEGER :: ia,iq,off1
|
|
REAL(DP) :: normr, veff_ia
|
|
!
|
|
CALL start_clock('tsvdw_veff')
|
|
!
|
|
! Initialization of effective volume...
|
|
!
|
|
ALLOCATE(veff(nat)); veff=0.0_DP
|
|
IF(.NOT. ALLOCATED(veff_pub)) ALLOCATE(veff_pub(nat)); veff_pub=0.0_DP
|
|
!
|
|
! Normalization factor for veff integral...
|
|
!
|
|
normr=omega/DBLE(nr1r*nr2r*nr3r)
|
|
!
|
|
! Loop over atoms in the simulation cell...
|
|
!
|
|
DO iproc=1,nstates(me)
|
|
!
|
|
! Connect processor number with atom...
|
|
!
|
|
ia=me+nproc_image*(iproc-1)
|
|
veff_ia=0.0_DP
|
|
!
|
|
! Loop over points in the (pre-screened) spherical atomic integration domain...
|
|
!
|
|
!$omp parallel do private(off1),reduction(+:veff_ia)
|
|
DO iq=1,NsomegaA(ia)
|
|
!
|
|
! Compute veff integrand and complete dispersion potential (functional derivative of veff(A) wrt charge density)...
|
|
!
|
|
! veff(A) = INT [|r-rA|^3*rhoA(|r-rA|)*rhotot(r)/rhosad(r)]
|
|
!
|
|
! dveff(A)/dn(r) = |r-rA|^3*rhoA(|r-rA|)/rhosad(r)
|
|
!
|
|
off1=somegaA(iq,1,iproc)+(somegaA(iq,2,iproc)-1)*nr1+(somegaA(iq,3,iproc)-1)*nr1*nr2 !global offset [nr1,nr2,nr3]
|
|
dveffAdn(iq,iproc)=dveffAdn(iq,iproc)/rhosad(off1)
|
|
!
|
|
! Increment veff...
|
|
!
|
|
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
|
|
!
|
|
veff_ia=veff_ia+(dveffAdn(iq,iproc)*rhotot(off1))
|
|
!
|
|
END IF
|
|
!
|
|
END DO !iq
|
|
!$omp end parallel do
|
|
!
|
|
! Apply final normalization to veff integral...
|
|
!
|
|
veff(ia)=normr*veff_ia
|
|
!
|
|
END DO !iproc
|
|
!
|
|
! Collect veff over all processors and broadcast...
|
|
!
|
|
CALL mp_sum(veff,intra_image_comm)
|
|
!
|
|
VefftsvdW = veff
|
|
veff_pub=veff
|
|
!
|
|
CALL stop_clock('tsvdw_veff')
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_veff
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_dveff()
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Calculate first derivative of veff wrt nuclear and cell displacements.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! Local variables
|
|
!
|
|
INTEGER :: ia,ib,ias,ibs,iq,ir,i,j,ipair,off1,ioff,boff
|
|
REAL(DP) :: spcutA,spcutB,spdA,spdB,dq(3),dqA(3),dqAs(3),dqB(3),dqBs(3),dqAmic,dqBmic,dABmic,normr
|
|
REAL(DP) :: dk1,dptmp1,dptmp2,dptmp3,dptmp4,dptmp5,rhoA,rhoB,drhoA,drhoB,dVAdRA(3),dVAdRB(3),dVBdRA(3)
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE :: predVAdRB
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE :: dqxyzr,dqAxyzs,predVBdRA
|
|
!
|
|
CALL start_clock('tsvdw_dveff')
|
|
!
|
|
! Initialization of the dveff/dR and dveff/dh arrays...
|
|
!
|
|
ALLOCATE(dveffdR(nat,nat,3)); dveffdR=0.0_DP
|
|
ALLOCATE(dveffdh(nat,3,3)); dveffdh=0.0_DP
|
|
!
|
|
! Normalization factor for dveff integrals...
|
|
!
|
|
normr=omega/DBLE(nr1r*nr2r*nr3r)
|
|
!
|
|
! Loop over atoms A in the simulation cell and compute dveffdR and dveffdh...
|
|
!
|
|
DO iproc=1,nstates(me)
|
|
!
|
|
! Connect processor number with atom...
|
|
!
|
|
ia=me+nproc_image*(iproc-1)
|
|
!
|
|
! Connect atom type with species-dependent quantities...
|
|
!
|
|
ias=ityp(ia)
|
|
!
|
|
! Transfer species-specific cutoff to spcutA...
|
|
!
|
|
spcutA=spdata(ias,1)
|
|
!
|
|
! Precompute inverse of species-specific linear grid spacing (replaces / with * during interpolation)...
|
|
!
|
|
spdA=1.0_DP/spdata(ias,2)
|
|
!
|
|
! Allocate and initialize atom-specific arrays...
|
|
!
|
|
ALLOCATE(dqxyzr(NsomegaAr(ia),3)); dqxyzr=0.0_DP
|
|
ALLOCATE(dqAxyzs(NsomegaAr(ia),3)); dqAxyzs=0.0_DP
|
|
ALLOCATE(predVAdRB(NsomegaAr(ia))); predVAdRB=0.0_DP
|
|
ALLOCATE(predVBdRA(NsomegaAr(ia),3)); predVBdRA=0.0_DP
|
|
!
|
|
! Initial loop over points in the (pre-screened) reduced spherical atomic integration domain for atom A to compute
|
|
! 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...
|
|
!
|
|
!$omp parallel do private(dq,dqA,dqAs,dqAmic,ir,dk1,rhoA,drhoA, &
|
|
!$omp off1,dptmp1,dptmp2,dptmp3,dptmp4,dVAdRA,dptmp5,i,j), &
|
|
!$omp reduction(+:dveffdh),reduction(+:dveffdR)
|
|
DO iq=1,NsomegaAr(ia)
|
|
!
|
|
! Compute global/cell reference frame Cartesian coordinates of given real-space grid point...
|
|
!
|
|
dq(1)=DBLE(somegaAr(iq,1,iproc)-1)/DBLE(nr1) ! s_i(1)
|
|
dq(2)=DBLE(somegaAr(iq,2,iproc)-1)/DBLE(nr2) ! s_i(2)
|
|
dq(3)=DBLE(somegaAr(iq,3,iproc)-1)/DBLE(nr3) ! s_i(3)
|
|
!
|
|
dqA(1)=h(1,1)*dq(1)+h(1,2)*dq(2)+h(1,3)*dq(3) ! r_i = h s_i
|
|
dqA(2)=h(2,1)*dq(1)+h(2,2)*dq(2)+h(2,3)*dq(3) ! r_i = h s_i
|
|
dqA(3)=h(3,1)*dq(1)+h(3,2)*dq(2)+h(3,3)*dq(3) ! r_i = h s_i
|
|
!
|
|
! Accumulate the Cartesian coordinates of the real-space grid point in the dqxyzr array:
|
|
!
|
|
! dqxyzr(:,1) := x-coordinate of grid point (global/cell reference frame)
|
|
! dqxyzr(:,2) := y-coordinate of grid point (global/cell reference frame)
|
|
! dqxyzr(:,3) := z-coordinate of grid point (global/cell reference frame)
|
|
!
|
|
dqxyzr(iq,1)=dqA(1)
|
|
dqxyzr(iq,2)=dqA(2)
|
|
dqxyzr(iq,3)=dqA(3)
|
|
!
|
|
! Compute distance between grid point and atom in scaled coordinates (s_iA) according to minimum image convention (MIC)...
|
|
!
|
|
dqA(1)=dqA(1)-atxyz(1,ia) ! r_iA = r_i - r_A
|
|
dqA(2)=dqA(2)-atxyz(2,ia) ! r_iA = r_i - r_A
|
|
dqA(3)=dqA(3)-atxyz(3,ia) ! r_iA = r_i - r_A
|
|
!
|
|
dqAs(1)=ainv(1,1)*dqA(1)+ainv(1,2)*dqA(2)+ainv(1,3)*dqA(3) ! s_iA = h^-1 r_iA
|
|
dqAs(2)=ainv(2,1)*dqA(1)+ainv(2,2)*dqA(2)+ainv(2,3)*dqA(3) ! s_iA = h^-1 r_iA
|
|
dqAs(3)=ainv(3,1)*dqA(1)+ainv(3,2)*dqA(2)+ainv(3,3)*dqA(3) ! s_iA = h^-1 r_iA
|
|
!
|
|
dqAs(1)=dqAs(1)-IDNINT(dqAs(1)) ! impose MIC on s_iA in range: [-0.5,+0.5]
|
|
dqAs(2)=dqAs(2)-IDNINT(dqAs(2)) ! impose MIC on s_iA in range: [-0.5,+0.5]
|
|
dqAs(3)=dqAs(3)-IDNINT(dqAs(3)) ! impose MIC on s_iA in range: [-0.5,+0.5]
|
|
!
|
|
! Accumulate the components of the s_i - s_A vector in the dqAxyzs array:
|
|
!
|
|
! dqAxyzs(:,1) := 1-coordinate of s_i - s_A vector (local/atom reference frame)
|
|
! dqAxyzs(:,2) := 2-coordinate of s_i - s_A vector (local/atom reference frame)
|
|
! dqAxyzs(:,3) := 3-coordinate of s_i - s_A vector (local/atom reference frame)
|
|
!
|
|
dqAxyzs(iq,1)=dqAs(1)
|
|
dqAxyzs(iq,2)=dqAs(2)
|
|
dqAxyzs(iq,3)=dqAs(3)
|
|
!
|
|
! Convert MIC distance components from scaled coordinates to Cartesian coordinates (s_iA -> r_iA)...
|
|
!
|
|
dqA(1)=h(1,1)*dqAs(1)+h(1,2)*dqAs(2)+h(1,3)*dqAs(3) ! r_iA = h s_iA (MIC)
|
|
dqA(2)=h(2,1)*dqAs(1)+h(2,2)*dqAs(2)+h(2,3)*dqAs(3) ! r_iA = h s_iA (MIC)
|
|
dqA(3)=h(3,1)*dqAs(1)+h(3,2)*dqAs(2)+h(3,3)*dqAs(3) ! r_iA = h s_iA (MIC)
|
|
!
|
|
dqAmic=DSQRT(dqA(1)*dqA(1)+dqA(2)*dqA(2)+dqA(3)*dqA(3)) ! |r_i - r_A| (MIC)
|
|
!
|
|
! Determine the index in the atomic linear equispaced grid such that grd(ir) <= dqA <= grd(ir+1) and distance between dqA and grd(ir)...
|
|
!
|
|
ir=INT(dqAmic*spdA)
|
|
dk1=dqAmic-spgrd(ias,ir)
|
|
!
|
|
! Perform linear interpolation to obtain the value of the atomic pseudo-density and its derivative at the given grid point...
|
|
!
|
|
rhoA=LIA(ias,ir)+LIB(ias,ir)*dk1 !rhoA at grid point via linear interpolation
|
|
drhoA=dLIA(ias,ir)+dLIB(ias,ir)*dk1 !drhoA at grid point via linear interpolation
|
|
!
|
|
! Compute global offset for rhosad(r) and rhotot(r), both computed on the real-space mesh...
|
|
!
|
|
off1=somegaAr(iq,1,iproc)+(somegaAr(iq,2,iproc)-1)*nr1+(somegaAr(iq,3,iproc)-1)*nr1*nr2 !global offset [nr1,nr2,nr3]
|
|
!
|
|
! Compute self-derivative dVA/dpA integrand for p={x,y,z}...
|
|
!
|
|
! 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|)]}
|
|
!
|
|
dptmp1=1.0_DP/rhosad(off1)
|
|
dptmp2=dqAmic*drhoA
|
|
dptmp3=dqAmic*rhotot(off1)*dptmp1
|
|
dptmp4=((rhoA*dptmp1-1.0_DP)*dptmp2-3.0_DP*rhoA)*dptmp3
|
|
!
|
|
dVAdRA=dqA*dptmp4 !dVA/dpA integrand/contribution for the given grid point...
|
|
!
|
|
! Increment self-derivative dVA/dpA for p={x,y,z}...
|
|
!
|
|
DO i=1,3
|
|
!
|
|
dveffdR(ia,ia,i)=dveffdR(ia,ia,i)+dVAdRA(i)
|
|
!
|
|
END DO !i
|
|
!
|
|
! Increment self-contribution to dVA/dhpq for p,q={x,y,z}...
|
|
!
|
|
! 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|)]}
|
|
!
|
|
DO i=1,3
|
|
!
|
|
DO j=1,3
|
|
!
|
|
dveffdh(ia,i,j)=dveffdh(ia,i,j)-dVAdRA(i)*dqAxyzs(iq,j)
|
|
!
|
|
END DO !j
|
|
!
|
|
END DO !i
|
|
!
|
|
! Precompute quantities necessary for dV(A)/dR(B) and dV(B)/dR(A)...
|
|
!
|
|
predVAdRB(iq)=dptmp1*rhoA*dqAmic*dqAmic*dptmp3
|
|
!
|
|
dptmp5=dptmp1*dptmp1*drhoA*rhotot(off1)
|
|
!
|
|
IF (dqAmic.LT.(1.0E-12_DP)) THEN
|
|
!
|
|
predVBdRA(iq,:)=dptmp5
|
|
!
|
|
ELSE
|
|
!
|
|
predVBdRA(iq,:)=dptmp5*dqA(:)/dqAmic
|
|
!
|
|
END IF
|
|
!
|
|
END DO !iq
|
|
!$omp end parallel do
|
|
!
|
|
! Inner loop over unique atom pairs B in the simulation cell to compute pair contributions to dveffdR and dveffdh...
|
|
!
|
|
!$omp parallel do private(dqB,dqBs,dqBmic,ir,dk1,rhoB,drhoB,dVAdRB,dVBdRA, &
|
|
!$omp i,j,ib,ibs,spcutB,spdB,ioff,boff), &
|
|
!$omp reduction(+:dveffdR),reduction(+:dveffdh)
|
|
DO ipair=1,npair(ia)
|
|
!
|
|
! Connect pair number with atom...
|
|
!
|
|
ib=pair(ipair,ia)
|
|
!
|
|
! Connect atom type with species-dependent quantities...
|
|
!
|
|
ibs=ityp(ib)
|
|
!
|
|
! Transfer species-specific cutoff to spcutB...
|
|
!
|
|
spcutB=spdata(ibs,1)
|
|
!
|
|
! Precompute inverse of species-specific linear grid spacing (replaces / with * during interpolation)...
|
|
!
|
|
spdB=1.0_DP/spdata(ibs,2)
|
|
!
|
|
! Determine atom B offsets for using gomegaAr bit array screening below...
|
|
!
|
|
ioff=((ib-1)/bsint)+1 ! integer offset for gomegaAr bit array
|
|
boff=(ib-((ioff-1)*bsint))-1 ! bit offset for gomegaAr bit array
|
|
!
|
|
! Inner loop over points in the (pre-screened) reduced spherical atomic integration domain for atom A to compute
|
|
! 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...
|
|
!
|
|
DO iq=1,NsomegaAr(ia)
|
|
!
|
|
! Determine if atom B contributes to the given point on the reduced spherical atomic integration domain on atom A (using gomegaAr bit array)...
|
|
!
|
|
IF (BTEST(gomegaAr(iq,ioff,iproc),boff)) THEN
|
|
!
|
|
! Compute distance between grid point and atom B according to minimum image convention (MIC)...
|
|
!
|
|
dqB(1)=dqxyzr(iq,1)-atxyz(1,ib) ! r_iB = r_i - r_B
|
|
dqB(2)=dqxyzr(iq,2)-atxyz(2,ib) ! r_iB = r_i - r_B
|
|
dqB(3)=dqxyzr(iq,3)-atxyz(3,ib) ! r_iB = r_i - r_B
|
|
!
|
|
dqBs(1)=ainv(1,1)*dqB(1)+ainv(1,2)*dqB(2)+ainv(1,3)*dqB(3) ! s_iB = h^-1 r_iB
|
|
dqBs(2)=ainv(2,1)*dqB(1)+ainv(2,2)*dqB(2)+ainv(2,3)*dqB(3) ! s_iB = h^-1 r_iB
|
|
dqBs(3)=ainv(3,1)*dqB(1)+ainv(3,2)*dqB(2)+ainv(3,3)*dqB(3) ! s_iB = h^-1 r_iB
|
|
!
|
|
dqBs(1)=dqBs(1)-IDNINT(dqBs(1)) ! impose MIC on s_iB in range: [-0.5,+0.5]
|
|
dqBs(2)=dqBs(2)-IDNINT(dqBs(2)) ! impose MIC on s_iB in range: [-0.5,+0.5]
|
|
dqBs(3)=dqBs(3)-IDNINT(dqBs(3)) ! impose MIC on s_iB in range: [-0.5,+0.5]
|
|
!
|
|
dqB(1)=h(1,1)*dqBs(1)+h(1,2)*dqBs(2)+h(1,3)*dqBs(3) ! r_iB = h s_iB (MIC)
|
|
dqB(2)=h(2,1)*dqBs(1)+h(2,2)*dqBs(2)+h(2,3)*dqBs(3) ! r_iB = h s_iB (MIC)
|
|
dqB(3)=h(3,1)*dqBs(1)+h(3,2)*dqBs(2)+h(3,3)*dqBs(3) ! r_iB = h s_iB (MIC)
|
|
!
|
|
dqBmic=DSQRT(dqB(1)*dqB(1)+dqB(2)*dqB(2)+dqB(3)*dqB(3)) ! |r_i - r_B| (MIC)
|
|
!
|
|
! Final screening based on the (pre-screened) spherical atomic integration domain on atom B...
|
|
!
|
|
IF (dqBmic.LE.spcutB) THEN
|
|
!
|
|
! Determine the index in the atomic linear equispaced grid such that grd(ir) <= dqB <= grd(ir+1) and distance between dqB and grd(ir)...
|
|
!
|
|
ir=INT(dqBmic*spdB)
|
|
dk1=dqBmic-spgrd(ibs,ir)
|
|
!
|
|
! Perform linear interpolation to obtain the value of the atomic pseudo-density and its derivative at the given grid point...
|
|
!
|
|
rhoB=LIA(ibs,ir)+LIB(ibs,ir)*dk1 !rhoB at grid point via linear interpolation
|
|
drhoB=dLIA(ibs,ir)+dLIB(ibs,ir)*dk1 !drhoB at grid point via linear interpolation
|
|
!
|
|
! Compute dVA/dpB integrand for p={x,y,z}...
|
|
!
|
|
! dVA/dpB = INT {(p-pB)/|r-rB|*[drho(|r-rB|)*|r-rA|^3*rho(|r-rA|)*rhotot(r)/rhosad(r)^2]}
|
|
!
|
|
IF (dqBmic.LT.(1.0E-12_DP)) THEN
|
|
!
|
|
dVAdRB(:)=predVAdRB(iq)*drhoB
|
|
!
|
|
ELSE
|
|
!
|
|
dVAdRB(:)=predVAdRB(iq)*drhoB*dqB(:)/dqBmic
|
|
!
|
|
END IF
|
|
!
|
|
! Increment non-self-derivative dVA/dpB for p={x,y,z}...
|
|
!
|
|
DO i=1,3
|
|
!
|
|
dveffdR(ia,ib,i)=dveffdR(ia,ib,i)+dVAdRB(i)
|
|
!
|
|
END DO !i
|
|
!
|
|
! Increment non-self-contribution to dVA/dhpq for p,q={x,y,z} from atom B...
|
|
!
|
|
! dVA/dhpq <-- INT {-(p-pB)*(qs-qsB)/|r-rB|*[drho(|r-rB|)*|r-rA|^3*rho(|r-rA|)*rhotot(r)/rhosad(r)^2]}
|
|
!
|
|
DO i=1,3
|
|
!
|
|
DO j=1,3
|
|
!
|
|
dveffdh(ia,i,j)=dveffdh(ia,i,j)-dVAdRB(i)*dqBs(j)
|
|
!
|
|
END DO !j
|
|
!
|
|
END DO !i
|
|
!
|
|
! Compute dVB/dpA integrand for p={x,y,z}...
|
|
!
|
|
! dVB/dpA = INT {(p-pA)/|r-rA|*[drho(|r-rA|)*|r-rB|^3*rho(|r-rB|)*rhotot(r)/rhosad(r)^2]}
|
|
!
|
|
dVBdRA(:)=predVBdRA(iq,:)*rhoB*dqBmic*dqBmic*dqBmic
|
|
!
|
|
! Increment non-self-derivative dVB/dpA for p={x,y,z} from atom A...
|
|
!
|
|
DO i=1,3
|
|
!
|
|
dveffdR(ib,ia,i)=dveffdR(ib,ia,i)+dVBdRA(i)
|
|
!
|
|
END DO !i
|
|
!
|
|
! Increment non-self-contribution to dVB/dhpq for p,q={x,y,z} from atom A...
|
|
!
|
|
! dVB/dhpq <-- INT {-(p-pA)*(qs-qsA)/|r-rA|*[drho(|r-rA|)*|r-rB|^3*rho(|r-rB|)*rhotot(r)/rhosad(r)^2]}
|
|
!
|
|
DO i=1,3
|
|
!
|
|
DO j=1,3
|
|
!
|
|
dveffdh(ib,i,j)=dveffdh(ib,i,j)-dVBdRA(i)*dqAxyzs(iq,j)
|
|
!
|
|
END DO !j
|
|
!
|
|
END DO !i
|
|
!
|
|
END IF
|
|
!
|
|
END IF !BTEST
|
|
!
|
|
END DO !iq
|
|
!
|
|
END DO !ipair
|
|
!$omp end parallel do
|
|
!
|
|
! Deallocate temporary arrays...
|
|
!
|
|
IF (ALLOCATED(dqxyzr)) DEALLOCATE(dqxyzr)
|
|
IF (ALLOCATED(dqAxyzs)) DEALLOCATE(dqAxyzs)
|
|
IF (ALLOCATED(predVAdRB)) DEALLOCATE(predVAdRB)
|
|
IF (ALLOCATED(predVBdRA)) DEALLOCATE(predVBdRA)
|
|
!
|
|
END DO !iproc
|
|
!
|
|
! Apply final normalization of dVA/dR integrals...
|
|
!
|
|
dveffdR=normr*dveffdR
|
|
!
|
|
! Apply final normalization of dVA/dhab integrals...
|
|
!
|
|
dveffdh=normr*dveffdh
|
|
!
|
|
! Collect dveffdR and dveffdh over all processors and broadcast...
|
|
!
|
|
CALL mp_sum(dveffdR,intra_image_comm)
|
|
CALL mp_sum(dveffdh,intra_image_comm)
|
|
!
|
|
CALL stop_clock('tsvdw_dveff')
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_dveff
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_effqnts()
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Calculate effective quantities for each atom in the simulation cell.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! Local variables
|
|
!
|
|
INTEGER :: ia,ib,ias,ibs
|
|
REAL(DP) :: vA,vB,num,den
|
|
!
|
|
! Initialization of base effective atomic quantities...
|
|
!
|
|
ALLOCATE(dpeff(nat)); dpeff=0.0_DP
|
|
ALLOCATE(R0eff(nat)); R0eff=0.0_DP
|
|
ALLOCATE(C6AAeff(nat)); C6AAeff=0.0_DP
|
|
ALLOCATE(C6ABeff(nat,nat)); C6ABeff=0.0_DP
|
|
!
|
|
! Population of base effective atomic quantities...
|
|
!
|
|
DO ia=1,nat
|
|
!
|
|
! Connect atom type with species-dependent quantities...
|
|
!
|
|
ias=ityp(ia)
|
|
!
|
|
! Precompute veff(A)/vfree(A) ratio...
|
|
!
|
|
vA=(veff(ia)/vfree(ias))
|
|
!
|
|
! Effective atomic static dipole polarizability array...
|
|
! dpeff(A)=[veff(A)/vfree(A)]*dpfree(A)
|
|
!
|
|
dpeff(ia)=vA*dpfree(ias)
|
|
!
|
|
! Effective atomic vdW radius array...
|
|
! R0eff(A)=[veff(A)/vfree(A)]^1/3*R0free(A)
|
|
!
|
|
R0eff(ia)=(vA**(1.0_DP/3.0_DP))*R0free(ias)
|
|
!
|
|
! Effective homonuclear C6 coefficient array...
|
|
! C6AAeff(A)=[veff(A)/vfree(A)]^2*C6AAfree(A)
|
|
!
|
|
C6AAeff(ia)=(vA**(2.0_DP))*C6AAfree(ias)
|
|
!
|
|
DO ib=1,nat
|
|
!
|
|
! Connect atom type with species-dependent quantities...
|
|
!
|
|
ibs=ityp(ib)
|
|
!
|
|
! Precompute veff(B)/vfree(B) ratio...
|
|
!
|
|
vB=(veff(ib)/vfree(ibs))
|
|
!
|
|
! Effective heteronuclear C6 coefficient matrix...
|
|
! C6ABeff(A,B)=(veff(A)/vfree(A))*(veff(B)/vfree(B))*C6ABfree(A,B)
|
|
!
|
|
C6ABeff(ia,ib)=(vA*vB)*C6ABfree(ias,ibs)
|
|
!
|
|
END DO !ib
|
|
!
|
|
END DO !ia
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_effqnts
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_energy()
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Calculate total TS-vdW energy, dispersion potential prefactor, ionic forces, and cell forces.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! Local variables
|
|
!
|
|
LOGICAL :: periodic_converged
|
|
INTEGER :: ia,ib,ic,ias,ibs,n_period,n1,n2,n3,i,j
|
|
|
|
REAL(DP) :: dAB(3),dAB2(3),dsAB(3),dABimg,dABimg2,dABimgn1,dABimgn2,dABimgn5,dABimgn6
|
|
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)
|
|
REAL(DP) :: EtsvdW_period,RAB0,edamp,fdamp,fdamp2,D1A,D1B,D2A,D2B,D12A,D12B,dptmp1,dptmp2,vtmp1(3),vtmp2(3)
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE :: predveffAdn_period
|
|
REAL(DP), DIMENSION(:,:), ALLOCATABLE :: FtsvdW_period
|
|
REAL(DP) :: HtsvdW_period(3,3)
|
|
!
|
|
CALL start_clock('tsvdw_energy')
|
|
!
|
|
! Initialize total TS-vdW energy, ion force, and cell force ...
|
|
!
|
|
EtsvdW=0.0_DP; FtsvdW=0.0_DP; HtsvdW=0.0_DP
|
|
!
|
|
! Allocate and initialize TS-vdW dispersion potential prefactor...
|
|
!
|
|
ALLOCATE(predveffAdn(nat)); predveffAdn=0.0_DP
|
|
!
|
|
! Allocate and initialize periodic contributions to TS-vdW ionic forces, cell forces, and dispersion potential prefactor...
|
|
!
|
|
ALLOCATE(FtsvdW_period(3,nat)); FtsvdW_period=0.0_DP
|
|
HtsvdW_period=0.0_DP
|
|
ALLOCATE(predveffAdn_period(nat)); predveffAdn_period=0.0_DP
|
|
!
|
|
! Precompute quantities outside all loops...
|
|
!
|
|
FDR0=ddamp/(2.0_DP*sR)
|
|
FDV0=-FDR0/3.0_DP
|
|
FCV0=0.5_DP
|
|
FRR0=-3.0_DP
|
|
!
|
|
! For periodic systems, converge the energy with respect to neighboring images...
|
|
!
|
|
n_period=0
|
|
periodic_converged=.FALSE.
|
|
!
|
|
DO WHILE (.NOT.periodic_converged)
|
|
!
|
|
EtsvdW_period=0.0_DP
|
|
FtsvdW_period=0.0_DP
|
|
HtsvdW_period=0.0_DP
|
|
predveffAdn_period=0.0_DP
|
|
!
|
|
! Outer loop over atoms A...
|
|
!
|
|
DO iproc=1,nstates(me)
|
|
!
|
|
! Connect processor number with atom...
|
|
!
|
|
ia=me+nproc_image*(iproc-1)
|
|
!
|
|
! Connect atom type with species-dependent quantities...
|
|
!
|
|
ias=ityp(ia)
|
|
!
|
|
! Precompute quantities outside loop over B...
|
|
!
|
|
FDV1=R0free(ias)/(vfree(ias)**(1.0_DP/3.0_DP)*veff(ia)**(2.0_DP/3.0_DP))
|
|
FCVA1=1.0_DP/vfree(ias)
|
|
FCVB1=veff(ia)*FCVA1
|
|
!
|
|
! Inner loop over atoms B...
|
|
!
|
|
|
|
!$omp parallel do private(ibs,RAB0,FRR2,FDR2,FDVA2,FDVB2,FCVB2,FCVA2, &
|
|
!$omp dAB,dAB2,FDVi,FDRi,FDRii,FCVi,FRRi,FRRii,n1,n2,n3,dsAB,dABimg2, &
|
|
!$omp dABimg,dABimgn1,dABimgn2,dABimgn5,dABimgn6,edamp,fdamp,fdamp2,dptmp1, &
|
|
!$omp dptmp2,i,j,vtmp1,vtmp2,D1A,D2A,D1B,D2B,D12A,D12B,ic), &
|
|
!$omp reduction(+:EtsvdW_period, FtsvdW_period, HtsvdW_period), &
|
|
!$omp reduction(+:predveffAdn_period)
|
|
DO ib=1,nat
|
|
!
|
|
! Connect atom type with species-dependent quantities...
|
|
!
|
|
ibs=ityp(ib)
|
|
!
|
|
! Compute RAB0 as the sum of the effective vdW radii of atoms A and B...
|
|
!
|
|
RAB0=R0eff(ia)+R0eff(ib)
|
|
!
|
|
! Precompute quantities outside loop over image cells...
|
|
!
|
|
FRR2=C6ABeff(ia,ib)
|
|
FDR2=FRR2/RAB0
|
|
FDVA2=FDR2/RAB0
|
|
FDVB2=FDVA2*R0free(ibs)/(vfree(ibs)**(1.0_DP/3.0_DP)*veff(ib)**(2.0_DP/3.0_DP))
|
|
FCVB2=C6ABfree(ias,ibs)/vfree(ibs)
|
|
FCVA2=FCVB2*veff(ib)
|
|
!
|
|
! Compute distance between atom A and atom B (according to the minimum image convention)...
|
|
!
|
|
dAB(1)=atxyz(1,ia)-atxyz(1,ib) ! r_AB = r_A - r_B
|
|
dAB(2)=atxyz(2,ia)-atxyz(2,ib) ! r_AB = r_A - r_B
|
|
dAB(3)=atxyz(3,ia)-atxyz(3,ib) ! r_AB = r_A - r_B
|
|
!
|
|
dAB2(1)=ainv(1,1)*dAB(1)+ainv(1,2)*dAB(2)+ainv(1,3)*dAB(3) ! s_AB = h^-1 r_AB
|
|
dAB2(2)=ainv(2,1)*dAB(1)+ainv(2,2)*dAB(2)+ainv(2,3)*dAB(3) ! s_AB = h^-1 r_AB
|
|
dAB2(3)=ainv(3,1)*dAB(1)+ainv(3,2)*dAB(2)+ainv(3,3)*dAB(3) ! s_AB = h^-1 r_AB
|
|
!
|
|
dAB2(1)=dAB2(1)-IDNINT(dAB2(1)) ! impose MIC on s_AB in range: [-0.5,+0.5]
|
|
dAB2(2)=dAB2(2)-IDNINT(dAB2(2)) ! impose MIC on s_AB in range: [-0.5,+0.5]
|
|
dAB2(3)=dAB2(3)-IDNINT(dAB2(3)) ! impose MIC on s_AB in range: [-0.5,+0.5]
|
|
!
|
|
! Initialize image-summed matrix elements...
|
|
!
|
|
FDVi=0.0_DP; FDRi=0.0_DP; FDRii=0.0_DP; FCVi=0.0_DP; FRRi=0.0_DP; FRRii=0.0_DP
|
|
!
|
|
! Loop over image cells...
|
|
!
|
|
DO n1=-n_period,n_period
|
|
!
|
|
DO n2=-n_period,n_period
|
|
!
|
|
DO n3=-n_period,n_period
|
|
!
|
|
IF ((ABS(n1).EQ.n_period).OR.(ABS(n2).EQ.n_period).OR.(ABS(n3).EQ.n_period)) THEN
|
|
!
|
|
! Recover MIC distance between atom A and atom B in crystal coordinates...
|
|
!
|
|
dsAB(1)=dAB2(1) ! s_AB (MIC)
|
|
dsAB(2)=dAB2(2) ! s_AB (MIC)
|
|
dsAB(3)=dAB2(3) ! s_AB (MIC)
|
|
!
|
|
! Increment MIC distance in crystal coordinates...
|
|
!
|
|
dsAB(1)=dsAB(1)+DBLE(n1) ! s_AB (incremented, MIC only if n_period == 0)
|
|
dsAB(2)=dsAB(2)+DBLE(n2) ! s_AB (incremented, MIC only if n_period == 0)
|
|
dsAB(3)=dsAB(3)+DBLE(n3) ! s_AB (incremented, MIC only if n_period == 0)
|
|
!
|
|
! Convert incremented distance back into cartesian coordinates...
|
|
!
|
|
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)
|
|
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)
|
|
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)
|
|
!
|
|
! Compute incremented distance between atom A and atom B...
|
|
!
|
|
dABimg2=dAB(1)*dAB(1)+dAB(2)*dAB(2)+dAB(3)*dAB(3)
|
|
dABimg=DSQRT(dABimg2)
|
|
!
|
|
! Precompute inverse powers of incremented distance between atom A and atom B...
|
|
!
|
|
IF ( dABimg > 0.0_dp ) THEN
|
|
dABimgn1=1.0_DP/dABimg
|
|
dABimgn2=dABimgn1*dABimgn1
|
|
dABimgn5=dABimgn2*dABimgn2*dABimgn1
|
|
dABimgn6=dABimgn5*dABimgn1
|
|
ELSE
|
|
dABimgn1=0.0_DP
|
|
dABimgn2=0.0_DP
|
|
dABimgn5=0.0_DP
|
|
dABimgn6=0.0_DP
|
|
END IF
|
|
!
|
|
! Precompute damping function (fdamp) and damping function exponential (edamp)...
|
|
!
|
|
edamp=EXP(-ddamp*(dABimg/(sR*RAB0)-1.0_DP))
|
|
fdamp=1.0_DP/(1.0_DP+edamp)
|
|
fdamp2=fdamp*fdamp
|
|
!
|
|
! Apply delta[ia;ib] x delta[n1,n2,n3;0,0,0] conditional...
|
|
!
|
|
IF (n_period.EQ.0.AND.ia.EQ.ib) THEN
|
|
!
|
|
! Do not include self-interaction in the simulation cell...
|
|
!
|
|
FDVi=FDVi+0.0_DP
|
|
FDRi=FDRi+0.0_DP
|
|
FDRii=FDRii+0.0_DP
|
|
FCVi=FCVi+0.0_DP
|
|
FRRi=FRRi+0.0_DP
|
|
FRRii=FRRii+0.0_DP
|
|
!
|
|
ELSE
|
|
!
|
|
! Increment image-summed matrix elements...
|
|
!
|
|
dptmp1=edamp*fdamp2*dABimgn5
|
|
FDVi=FDVi+dptmp1
|
|
!
|
|
dptmp2=fdamp*dABimgn6
|
|
FCVi=FCVi+dptmp2
|
|
!
|
|
dptmp1=dptmp1*dABimgn2
|
|
dptmp2=dptmp2*dABimgn2
|
|
!
|
|
DO i=1,3
|
|
!
|
|
vtmp1(i)=dptmp1*dAB(i)
|
|
FDRi(i)=FDRi(i)+vtmp1(i)
|
|
!
|
|
vtmp2(i)=dptmp2*dAB(i)
|
|
FRRi(i)=FRRi(i)+vtmp2(i)
|
|
!
|
|
DO j=1,3
|
|
!
|
|
FDRii(i,j)=FDRii(i,j)+vtmp1(i)*dsAB(j)
|
|
FRRii(i,j)=FRRii(i,j)+vtmp2(i)*dsAB(j)
|
|
!
|
|
END DO
|
|
!
|
|
END DO
|
|
!
|
|
END IF
|
|
!
|
|
END IF !n_period conditional
|
|
!
|
|
END DO !n3
|
|
!
|
|
END DO !n2
|
|
!
|
|
END DO !n1
|
|
!
|
|
! Increment period energy via EtsvdWAB = - 1/2 * C6ABeff * FAB...
|
|
!
|
|
EtsvdW_period=EtsvdW_period-(FCV0*FRR2*FCVi)
|
|
!
|
|
! Increment dispersion potential (predveffAdn) prefactor...
|
|
!
|
|
! predveffAdn(A) = (d * R0freeA * C6ABeff * edamp * fdamp^2) / (6 * sR * vfreeA^1/3 * veffA^2/3 * RAB0^2 * RAB^5)
|
|
! - (C6ABfree * veffB * fdamp) / (2 * vfreeA * vfreeB * RAB^6)
|
|
!
|
|
! predveffAdn(B) = (d * R0freeB * C6ABeff * edamp * fdamp^2) / (6 * sR * vfreeB^1/3 * veffB^2/3 * RAB0^2 * RAB^5)
|
|
! - (C6ABfree * veffA * fdamp) / (2 * vfreeA * vfreeB * RAB^6)
|
|
!
|
|
predveffAdn_period(ia)=predveffAdn_period(ia)-(FDV0*FDV1*FDVA2*FDVi+FCV0*FCVA1*FCVA2*FCVi)
|
|
predveffAdn_period(ib)=predveffAdn_period(ib)-(FDV0*FDVB2*FDVi+FCV0*FCVB1*FCVB2*FCVi)
|
|
!
|
|
! Increment effective volume derivative contributions to ionic and cell forces...
|
|
!
|
|
! (dfdamp/dVA) --> D1A = - (d * R0freeA * C6ABeff * edamp * fdamp^2) / (6 * sR * vfreeA^1/3 * veffA^2/3 * RAB0^2 * RAB^5)
|
|
!
|
|
! (dfdamp/dVB) --> D1B = - (d * R0freeB * C6ABeff * edamp * fdamp^2) / (6 * sR * vfreeB^1/3 * veffB^2/3 * RAB0^2 * RAB^5)
|
|
!
|
|
! (dC6AB/dVA) --> D2A = (C6ABfree * veffB * fdamp) / (2 * vfreeA * vfreeB * RAB^6)
|
|
!
|
|
! (dC6AB/dVB) --> D2B = (C6ABfree * veffA * fdamp) / (2 * vfreeA * vfreeB * RAB^6)
|
|
!
|
|
D1A=FDV0*FDV1*FDVA2*FDVi ! (dfdamp/dVA)
|
|
D2A=FCV0*FCVA1*FCVA2*FCVi ! (dC6AB/dVA)
|
|
!
|
|
D1B=FDV0*FDVB2*FDVi ! (dfdamp/dVB)
|
|
D2B=FCV0*FCVB1*FCVB2*FCVi ! (dC6AB/dVB)
|
|
!
|
|
D12A=D1A+D2A; D12B=D1B+D2B
|
|
!
|
|
DO i=1,3
|
|
!
|
|
DO ic=1,nat
|
|
!
|
|
FtsvdW_period(i,ic)=FtsvdW_period(i,ic)+(dveffdR(ia,ic,i)*D12A+dveffdR(ib,ic,i)*D12B)
|
|
!
|
|
END DO
|
|
!
|
|
DO j=1,3
|
|
!
|
|
HtsvdW_period(i,j)=HtsvdW_period(i,j)+(dveffdh(ia,i,j)*D12A+dveffdh(ib,i,j)*D12B)
|
|
!
|
|
END DO
|
|
!
|
|
END DO
|
|
!
|
|
! Increment RAB derivative contributions to ionic and cell forces...
|
|
!
|
|
! (dfdamp/dRA) --> D1A = (d * C6ABeff * edamp * fdamp^2) / (2 * sR * RAB0 * RAB^7)
|
|
!
|
|
! (dfdamp/dRB) --> D1B = - D1A
|
|
!
|
|
! (dRAB^-6/dRA) --> D2A = - (3 * C6ABeff * fdamp) / (RAB^8)
|
|
!
|
|
! (dRAB^-6/dRB) --> D2B = - D2A
|
|
!
|
|
D1A=FDR0*FDR2 ! (dfdamp/dRA)
|
|
D2A=FRR0*FRR2 ! (dRAB^-6/dRA)
|
|
!
|
|
! N.B.: Manually zero out the force contribution from an atom in the simulation cell and
|
|
! any of its images (this applies to distance derivatives NOT volume derivatives)...
|
|
!
|
|
IF (ia.NE.ib) THEN
|
|
!
|
|
DO i=1,3
|
|
!
|
|
FtsvdW_period(i,ia)=FtsvdW_period(i,ia)+(D1A*FDRi(i)+D2A*FRRi(i))
|
|
FtsvdW_period(i,ib)=FtsvdW_period(i,ib)+(-D1A*FDRi(i)-D2A*FRRi(i))
|
|
!
|
|
DO j=1,3
|
|
!
|
|
HtsvdW_period(i,j)=HtsvdW_period(i,j)+(D1A*FDRii(i,j)+D2A*FRRii(i,j))
|
|
!
|
|
END DO
|
|
!
|
|
END DO
|
|
!
|
|
END IF
|
|
!
|
|
END DO !ib
|
|
!$omp end parallel do
|
|
!
|
|
END DO !iproc
|
|
!
|
|
! Synchronize n_period contribution from all processors...
|
|
!
|
|
CALL mp_sum(EtsvdW_period,intra_image_comm)
|
|
CALL mp_sum(FtsvdW_period,intra_image_comm)
|
|
CALL mp_sum(HtsvdW_period,intra_image_comm)
|
|
CALL mp_sum(predveffAdn_period,intra_image_comm)
|
|
!
|
|
! Increment total quantities...
|
|
!
|
|
EtsvdW=EtsvdW+EtsvdW_period !EvdW
|
|
FtsvdW=FtsvdW+FtsvdW_period !(-dE/dR)
|
|
HtsvdW=HtsvdW-HtsvdW_period !(dE/dh)
|
|
predveffAdn=predveffAdn+predveffAdn_period !(dE/dVA) & (dE/dVB)
|
|
!
|
|
! DEBUGGING
|
|
!WRITE(stdout,'(I10,2F25.12)') n_period,EtsvdW_period,EtsvdW
|
|
! DEBUGGING
|
|
!
|
|
! Periodic convergence loop conditionals...
|
|
!
|
|
IF ( vdw_isolated .OR. (ABS(EtsvdW_period) <= vdw_econv_thr) ) periodic_converged=.TRUE.
|
|
!
|
|
n_period=n_period+1
|
|
!
|
|
END DO !convergence loop
|
|
!
|
|
! Deallocate temporary arrays...
|
|
!
|
|
IF (ALLOCATED(FtsvdW_period)) DEALLOCATE(FtsvdW_period)
|
|
IF (ALLOCATED(predveffAdn_period)) DEALLOCATE(predveffAdn_period)
|
|
!
|
|
CALL stop_clock('tsvdw_energy')
|
|
!
|
|
!! DEBUGGING
|
|
!WRITE(stdout,'(3X,"<START veff ARRAY>")')
|
|
!DO ia=1,nat
|
|
! WRITE(stdout,'(5X,I5,F25.12)') ia,veff(ia)
|
|
!END DO
|
|
!WRITE(stdout,'(3X,"<END veff ARRAY>")')
|
|
!!
|
|
!WRITE(stdout,'(3X,"<START FtsvdW ARRAY>")')
|
|
!DO ia=1,nat
|
|
! WRITE(stdout,'(5X,I5,3F25.12)') ia,FtsvdW(:,ia)
|
|
!END DO
|
|
!WRITE(stdout,'(3X,"<END FtsvdW ARRAY>")')
|
|
!!
|
|
!WRITE(stdout,'(3X,"<START HtsvdW ARRAY>")')
|
|
!DO i=1,3
|
|
! WRITE(stdout,'(5X,3F25.12)') HtsvdW(i,:)
|
|
!END DO
|
|
!WRITE(stdout,'(3X,"<END HtsvdW ARRAY>")')
|
|
!! DEBUGGING
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_energy
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_wfforce()
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Calculate total TS-vdW wavefunction forces (dispersion potential).
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! Local variables
|
|
!
|
|
INTEGER :: ia,ip,iq,off1
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE :: UtsvdWA
|
|
!
|
|
CALL start_clock('tsvdw_wfforce')
|
|
!
|
|
! Initialization of UtsvdwA array...
|
|
!
|
|
ALLOCATE(UtsvdWA(nr1*nr2*nr3)); UtsvdWA=0.0_DP
|
|
!
|
|
! Loop over atoms and populate UtsvdWA from predveffAdn and dveffAdn...
|
|
!
|
|
DO iproc=1,nstates(me)
|
|
!
|
|
! Connect processor number with atom...
|
|
!
|
|
ia=me+nproc_image*(iproc-1)
|
|
!
|
|
! Loop over points in the (pre-screened) spherical atomic integration domain...
|
|
!
|
|
!$omp parallel do private(off1)
|
|
DO iq=1,NsomegaA(ia)
|
|
!
|
|
off1=somegaA(iq,1,iproc)+(somegaA(iq,2,iproc)-1)*nr1+(somegaA(iq,3,iproc)-1)*nr1*nr2 !global offset [nr1,nr2,nr3]
|
|
UtsvdWA(off1)=UtsvdWA(off1)+predveffAdn(ia)*dveffAdn(iq,iproc)
|
|
!
|
|
END DO !iq
|
|
!$omp end parallel do
|
|
!
|
|
END DO !iproc
|
|
!
|
|
! Collect UtsvdWA over all processors and broadcast...
|
|
!
|
|
CALL mp_sum(UtsvdWA,intra_image_comm)
|
|
!
|
|
! Partition out dispersion potential consistent with slabs of the charge density...
|
|
!
|
|
IF (dfftp%nr3p(me_bgrp+1).NE.0) THEN
|
|
!
|
|
!$omp parallel do
|
|
DO ip=1,dfftp%nr3p(me_bgrp+1)*nr1*nr2
|
|
!
|
|
UtsvdW(ip)=UtsvdWA(ip+rdispls(me_bgrp+1))
|
|
!
|
|
END DO
|
|
!$omp end parallel do
|
|
!
|
|
END IF
|
|
!
|
|
! Deallocate temporary arrays...
|
|
!
|
|
IF (ALLOCATED(UtsvdWA)) DEALLOCATE(UtsvdWA)
|
|
!
|
|
CALL stop_clock('tsvdw_wfforce')
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_wfforce
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_cleanup()
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Deallocate all arrays specific to \texttt{tsvdw\_calculate}.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! Deallocate tsvdw_calculate specific arrays...
|
|
!
|
|
IF (ALLOCATED(atxyz)) DEALLOCATE(atxyz)
|
|
IF (ALLOCATED(rhosad)) DEALLOCATE(rhosad)
|
|
IF (ALLOCATED(rhotot)) DEALLOCATE(rhotot)
|
|
IF (ALLOCATED(veff)) DEALLOCATE(veff)
|
|
IF (ALLOCATED(dpeff)) DEALLOCATE(dpeff)
|
|
IF (ALLOCATED(R0eff)) DEALLOCATE(R0eff)
|
|
IF (ALLOCATED(C6AAeff)) DEALLOCATE(C6AAeff)
|
|
IF (ALLOCATED(C6ABeff)) DEALLOCATE(C6ABeff)
|
|
IF (ALLOCATED(dveffdR)) DEALLOCATE(dveffdR)
|
|
IF (ALLOCATED(dveffdh)) DEALLOCATE(dveffdh)
|
|
IF (ALLOCATED(somegaA)) DEALLOCATE(somegaA)
|
|
IF (ALLOCATED(somegaAr)) DEALLOCATE(somegaAr)
|
|
IF (ALLOCATED(gomegaAr)) DEALLOCATE(gomegaAr)
|
|
IF (ALLOCATED(NsomegaA)) DEALLOCATE(NsomegaA)
|
|
IF (ALLOCATED(NsomegaAr)) DEALLOCATE(NsomegaAr)
|
|
IF (ALLOCATED(nstates)) DEALLOCATE(nstates)
|
|
IF (ALLOCATED(sdispls)) DEALLOCATE(sdispls)
|
|
IF (ALLOCATED(sendcount)) DEALLOCATE(sendcount)
|
|
IF (ALLOCATED(rdispls)) DEALLOCATE(rdispls)
|
|
IF (ALLOCATED(recvcount)) DEALLOCATE(recvcount)
|
|
IF (ALLOCATED(istatus)) DEALLOCATE(istatus)
|
|
IF (ALLOCATED(npair)) DEALLOCATE(npair)
|
|
IF (ALLOCATED(pair)) DEALLOCATE(pair)
|
|
IF (ALLOCATED(dveffAdn)) DEALLOCATE(dveffAdn)
|
|
IF (ALLOCATED(predveffAdn)) DEALLOCATE(predveffAdn)
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_cleanup
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE tsvdw_finalize()
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Deallocate module-specific arrays.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
IF (ALLOCATED(UtsvdW)) DEALLOCATE(UtsvdW)
|
|
IF (ALLOCATED(FtsvdW)) DEALLOCATE(FtsvdW)
|
|
IF (ALLOCATED(HtsvdW)) DEALLOCATE(HtsvdW)
|
|
IF (ALLOCATED(VefftsvdW))DEALLOCATE(VefftsvdW)
|
|
IF (ALLOCATED(vfree)) DEALLOCATE(vfree)
|
|
IF (ALLOCATED(dpfree)) DEALLOCATE(dpfree)
|
|
IF (ALLOCATED(R0free)) DEALLOCATE(R0free)
|
|
IF (ALLOCATED(C6AAfree)) DEALLOCATE(C6AAfree)
|
|
IF (ALLOCATED(C6ABfree)) DEALLOCATE(C6ABfree)
|
|
IF (ALLOCATED(spgrd)) DEALLOCATE(spgrd)
|
|
IF (ALLOCATED(sprho)) DEALLOCATE(sprho)
|
|
IF (ALLOCATED(spdrho)) DEALLOCATE(spdrho)
|
|
IF (ALLOCATED(spdata)) DEALLOCATE(spdata)
|
|
IF (ALLOCATED(LIA)) DEALLOCATE(LIA)
|
|
IF (ALLOCATED(LIB)) DEALLOCATE(LIB)
|
|
IF (ALLOCATED(dLIA)) DEALLOCATE(dLIA)
|
|
IF (ALLOCATED(dLIB)) DEALLOCATE(dLIB)
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE tsvdw_finalize
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE Num1stDer(r,f,N,h,df)
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Compute first derivative on linear mesh and then transform back to radial/exponential grid.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! I/O variables
|
|
!
|
|
INTEGER :: N
|
|
REAL(DP) :: h,r(N),f(N),df(N)
|
|
!
|
|
! Local variables
|
|
!
|
|
INTEGER, PARAMETER :: Ndp=7 !using 7-point formulae...
|
|
INTEGER :: ip,ir
|
|
INTEGER :: A1(Ndp),A2(Ndp),A3(Ndp),A4(Ndp),A5(Ndp),A6(Ndp),A7(Ndp)
|
|
REAL(DP) :: dsum
|
|
!
|
|
! Populate Bickley coefficient vectors (Math. Gaz.,v25,p19-27,1941) according to cases...
|
|
!
|
|
DATA A1/ -1764_DP, 4320_DP, -5400_DP, 4800_DP, -2700_DP, 864_DP, -120_DP /
|
|
DATA A2/ -120_DP, -924_DP, 1800_DP, -1200_DP, 600_DP, -180_DP, 24_DP /
|
|
DATA A3/ 24_DP, -288_DP, -420_DP, 960_DP, -360_DP, 96_DP, -12_DP /
|
|
DATA A4/ -12_DP, 108_DP, -540_DP, 0_DP, 540_DP, -108_DP, 12_DP /
|
|
DATA A5/ 12_DP, -96_DP, 360_DP, -960_DP, 420_DP, 288_DP, -24_DP /
|
|
DATA A6/ -24_DP, 180_DP, -600_DP, 1200_DP, -1800_DP, 924_DP, 120_DP /
|
|
DATA A7/ 120_DP, -864_DP, 2700_DP, -4800_DP, 5400_DP, -4320_DP, 1764_DP /
|
|
!
|
|
! Compute first derivative on linear mesh and then transform back to radial/exponential grid...
|
|
!
|
|
DO ir=1,N
|
|
!
|
|
dsum=0.0_DP
|
|
!
|
|
! Deal with different cases one-by-one...
|
|
!
|
|
IF (ir.EQ.1) THEN
|
|
DO ip=1,Ndp
|
|
dsum=dsum+A1(ip)*f(ir-1+ip)
|
|
END DO
|
|
ELSE IF (ir.EQ.2) THEN
|
|
DO ip=1,Ndp
|
|
dsum=dsum+A2(ip)*f(ir-2+ip)
|
|
END DO
|
|
ELSE IF (ir.EQ.3) THEN
|
|
DO ip=1,Ndp
|
|
dsum=dsum+A3(ip)*f(ir-3+ip)
|
|
END DO
|
|
ELSE IF (ir.GE.4.AND.ir.LE.N-3) THEN
|
|
DO ip=1,Ndp
|
|
dsum=dsum+A4(ip)*f(ir-4+ip)
|
|
END DO
|
|
ELSE IF (ir.EQ.N-2) THEN
|
|
DO ip=1,Ndp
|
|
dsum=dsum+A5(ip)*f(ir-5+ip)
|
|
END DO
|
|
ELSE IF (ir.EQ.N-1) THEN
|
|
DO ip=1,Ndp
|
|
dsum=dsum+A6(ip)*f(ir-6+ip)
|
|
END DO
|
|
ELSE IF (ir.EQ.N) THEN
|
|
DO ip=1,Ndp
|
|
dsum=dsum+A7(ip)*f(ir-7+ip)
|
|
END DO
|
|
ELSE
|
|
WRITE(stdout,'("Error in Num1stDer subroutine...")')
|
|
END IF
|
|
!
|
|
! Final Normalization...
|
|
!
|
|
df(ir)=dsum/(720.0_DP*h)
|
|
!
|
|
END DO !ir
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE Num1stDer
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE CubSplCoeff(r,f,N,df,d2f)
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Compute second derivatives at each of the atomic radial grid points using the cubic spline methodology (i.e., smooth &
|
|
!! continuous piecewise first and second derivatives). These second derivatives will be utilized during cubic spline interpolation
|
|
!! as a higher accuracy alternative to linear interpolation during the construction of the linear atomic grids. The two-parameter
|
|
!! boundary conditions that will be utilized below are known as a clamped cubic spline in that the first derivative at both the
|
|
!! first and last grid point were computed numerically and provided as input.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! I/O variables
|
|
!
|
|
INTEGER :: N
|
|
REAL(DP) :: r(N),f(N),df(N),d2f(N)
|
|
!
|
|
! Local variables
|
|
!
|
|
INTEGER :: i,j
|
|
REAL(DP) :: dy1,dyn,p,q,un,qn
|
|
REAL(DP), DIMENSION(:), ALLOCATABLE :: work
|
|
!
|
|
!
|
|
ALLOCATE(work(N)); work=0.0_DP
|
|
!
|
|
d2f=0.0_DP
|
|
!
|
|
! Enforce 'clamped' boundary condition at the first radial grid point...
|
|
!
|
|
dy1=df(1)
|
|
d2f(1)=-0.5_DP
|
|
work(1)=(3.0_DP/(r(2)-r(1)))*((f(2)-f(1))/(r(2)-r(1))-dy1)
|
|
!
|
|
! Decomposition loop of the tridiagonal algorithm for the second derivatives...
|
|
!
|
|
DO i=2,N-1
|
|
p=(r(i)-r(i-1))/(r(i+1)-r(i-1))
|
|
q=p*d2f(i-1)+2.0_DP
|
|
d2f(i)=(p-1.0_DP)/q
|
|
work(i)=(f(i+1)-f(i))/(r(i+1)-r(i))-(f(i)-f(i-1))/(r(i)-r(i-1))
|
|
work(i)=(6.0_DP*work(i)/(r(i+1)-r(i-1))-p*work(i-1))/q
|
|
END DO
|
|
!
|
|
! Enforce 'clamped' boundary condition at the last radial grid point...
|
|
!
|
|
dyn=df(N)
|
|
qn=0.5_DP
|
|
un=(3.0_DP/(r(N)-r(N-1)))*(dyn-(f(N)-f(N-1))/(r(N)-r(N-1)))
|
|
d2f(N)=(un-qn*work(N-1))/(qn*d2f(N-1)+1.0_DP)
|
|
!
|
|
! Back substitution loop of the tridiagonal algorithm for the second derivatives...
|
|
!
|
|
DO j=N-1,1,-1
|
|
d2f(j)=d2f(j)*d2f(j+1)+work(j)
|
|
END DO
|
|
!
|
|
! Clean-up and return home...
|
|
!
|
|
DEALLOCATE(work)
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE CubSplCoeff
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
SUBROUTINE GetVdWParam(atom,C6,alpha,R0)
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!! Get VdW parameters for each atom.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
! I/O variables
|
|
!
|
|
CHARACTER(LEN=2) :: atom
|
|
REAL(DP) :: C6,alpha,R0
|
|
!
|
|
SELECT CASE (atom)
|
|
!
|
|
CASE ('H')
|
|
alpha=4.500000_DP
|
|
C6=6.500000_DP
|
|
R0=3.100000_DP
|
|
!
|
|
CASE ('He')
|
|
alpha=1.380000_DP
|
|
C6=1.460000_DP
|
|
R0=2.650000_DP
|
|
!
|
|
CASE ('Li')
|
|
alpha=164.200000_DP
|
|
C6=1387.000000_DP
|
|
R0=4.160000_DP
|
|
!
|
|
CASE ('Be')
|
|
alpha=38.000000_DP
|
|
C6=214.000000_DP
|
|
R0=4.170000_DP
|
|
!
|
|
CASE ('B')
|
|
alpha=21.000000_DP
|
|
C6=99.500000_DP
|
|
R0=3.890000_DP
|
|
!
|
|
CASE ('C')
|
|
alpha=12.000000_DP
|
|
C6=46.600000_DP
|
|
R0=3.590000_DP
|
|
!
|
|
CASE ('N')
|
|
alpha=7.400000_DP
|
|
C6=24.200000_DP
|
|
R0=3.340000_DP
|
|
!
|
|
CASE ('O')
|
|
alpha=5.400000_DP
|
|
C6=15.600000_DP
|
|
R0=3.190000_DP
|
|
!
|
|
CASE ('F')
|
|
alpha=3.800000_DP
|
|
C6=9.520000_DP
|
|
R0=3.040000_DP
|
|
!
|
|
CASE ('Ne')
|
|
alpha=2.670000_DP
|
|
C6=6.380000_DP
|
|
R0=2.910000_DP
|
|
!
|
|
CASE ('Na')
|
|
alpha=162.700000_DP
|
|
C6=1556.000000_DP
|
|
R0=3.730000_DP
|
|
!
|
|
CASE ('Mg')
|
|
alpha=71.000000_DP
|
|
C6=627.000000_DP
|
|
R0=4.270000_DP
|
|
!
|
|
CASE ('Al')
|
|
alpha=60.000000_DP
|
|
C6=528.000000_DP
|
|
R0=4.330000_DP
|
|
!
|
|
CASE ('Si')
|
|
alpha=37.000000_DP
|
|
C6=305.000000_DP
|
|
R0=4.200000_DP
|
|
!
|
|
CASE ('P')
|
|
alpha=25.000000_DP
|
|
C6=185.000000_DP
|
|
R0=4.010000_DP
|
|
!
|
|
CASE ('S')
|
|
alpha=19.600000_DP
|
|
C6=134.000000_DP
|
|
R0=3.860000_DP
|
|
!
|
|
CASE ('Cl')
|
|
alpha=15.000000_DP
|
|
C6=94.600000_DP
|
|
R0=3.710000_DP
|
|
!
|
|
CASE ('Ar')
|
|
alpha=11.100000_DP
|
|
C6=64.300000_DP
|
|
R0=3.550000_DP
|
|
!
|
|
CASE ('K')
|
|
alpha=292.900000_DP
|
|
C6=3897.000000_DP
|
|
R0=3.710000_DP
|
|
!
|
|
CASE ('Ca')
|
|
alpha=160.000000_DP
|
|
C6=2221.000000_DP
|
|
R0=4.650000_DP
|
|
!
|
|
CASE ('Sc')
|
|
alpha=120.000000_DP
|
|
C6=1383.000000_DP
|
|
R0=4.590000_DP
|
|
!
|
|
CASE ('Ti')
|
|
alpha=98.000000_DP
|
|
C6=1044.000000_DP
|
|
R0=4.510000_DP
|
|
!
|
|
CASE ('V')
|
|
alpha=84.000000_DP
|
|
C6=832.000000_DP
|
|
R0=4.440000_DP
|
|
!
|
|
CASE ('Cr')
|
|
alpha=78.000000_DP
|
|
C6=602.000000_DP
|
|
R0=3.990000_DP
|
|
!
|
|
CASE ('Mn')
|
|
alpha=63.000000_DP
|
|
C6=552.000000_DP
|
|
R0=3.970000_DP
|
|
!
|
|
CASE ('Fe')
|
|
alpha=56.000000_DP
|
|
C6=482.000000_DP
|
|
R0=4.230000_DP
|
|
!
|
|
CASE ('Co')
|
|
alpha=50.000000_DP
|
|
C6=408.000000_DP
|
|
R0=4.180000_DP
|
|
!
|
|
CASE ('Ni')
|
|
alpha=48.000000_DP
|
|
C6=373.000000_DP
|
|
R0=3.820000_DP
|
|
!
|
|
CASE ('Cu')
|
|
alpha=42.000000_DP
|
|
C6=253.000000_DP
|
|
R0=3.760000_DP
|
|
!
|
|
CASE ('Zn')
|
|
alpha=40.000000_DP
|
|
C6=284.000000_DP
|
|
R0=4.020000_DP
|
|
!
|
|
CASE ('Ga')
|
|
alpha=60.000000_DP
|
|
C6=498.000000_DP
|
|
R0=4.190000_DP
|
|
!
|
|
CASE ('Ge')
|
|
alpha=41.000000_DP
|
|
C6=354.000000_DP
|
|
R0=4.200000_DP
|
|
!
|
|
CASE ('As')
|
|
alpha=29.000000_DP
|
|
C6=246.000000_DP
|
|
R0=4.110000_DP
|
|
!
|
|
CASE ('Se')
|
|
alpha=25.000000_DP
|
|
C6=210.000000_DP
|
|
R0=4.040000_DP
|
|
!
|
|
CASE ('Br')
|
|
alpha=20.000000_DP
|
|
C6=162.000000_DP
|
|
R0=3.930000_DP
|
|
!
|
|
CASE ('Kr')
|
|
alpha=16.800000_DP
|
|
C6=129.600000_DP
|
|
R0=3.820000_DP
|
|
!
|
|
CASE ('Rb')
|
|
alpha=319.200000_DP
|
|
C6=4691.000000_DP
|
|
R0=3.720000_DP
|
|
!
|
|
CASE ('Sr')
|
|
alpha=199.000000_DP
|
|
C6=3170.000000_DP
|
|
R0=4.540000_DP
|
|
!
|
|
CASE ('Y')
|
|
alpha=126.7370_DP
|
|
C6=1968.580_DP
|
|
R0=4.81510_DP
|
|
!
|
|
CASE ('Zr')
|
|
alpha=119.97_DP
|
|
C6=1677.91_DP
|
|
R0=4.53_DP
|
|
!
|
|
CASE ('Nb')
|
|
alpha=101.603_DP
|
|
C6=1263.61_DP
|
|
R0=4.2365_DP
|
|
!
|
|
CASE ('Mo')
|
|
alpha=88.4225785_DP
|
|
C6=1028.73_DP
|
|
R0=4.099_DP
|
|
!
|
|
CASE ('Tc')
|
|
alpha=80.083_DP
|
|
C6=1390.87_DP
|
|
R0=4.076_DP
|
|
!
|
|
CASE ('Ru')
|
|
alpha=65.8950_DP
|
|
C6=609.754_DP
|
|
R0=3.99530_DP
|
|
!
|
|
CASE ('Rh')
|
|
alpha=56.1_DP
|
|
C6=469.0_DP
|
|
R0=3.95_DP
|
|
!
|
|
CASE ('Pd')
|
|
alpha=23.680000_DP
|
|
C6=157.500000_DP
|
|
R0=3.66000_DP
|
|
!
|
|
CASE ('Ag')
|
|
alpha=50.600000_DP
|
|
C6=339.000000_DP
|
|
R0=3.820000_DP
|
|
!
|
|
CASE ('Cd')
|
|
alpha=39.7_DP
|
|
C6=452.0_DP
|
|
R0=3.99_DP
|
|
!
|
|
CASE ('In')
|
|
alpha=70.22000_DP
|
|
C6=707.046000_DP
|
|
R0=4.23198000_DP
|
|
!
|
|
CASE ('Sn')
|
|
alpha=55.9500_DP
|
|
C6=587.41700_DP
|
|
R0=4.303000_DP
|
|
!
|
|
CASE ('Sb')
|
|
alpha=43.671970_DP
|
|
C6=459.322_DP
|
|
R0=4.2760_DP
|
|
!
|
|
CASE ('Te')
|
|
alpha=37.65_DP
|
|
C6=396.0_DP
|
|
R0=4.22_DP
|
|
!
|
|
CASE ('I')
|
|
alpha=35.000000_DP
|
|
C6=385.000000_DP
|
|
R0=4.170000_DP
|
|
!
|
|
CASE ('Xe')
|
|
alpha=27.300000_DP
|
|
C6=285.900000_DP
|
|
R0=4.080000_DP
|
|
!
|
|
CASE ('Cs')
|
|
alpha=427.12_DP
|
|
C6=6582.08_DP
|
|
R0=3.78_DP
|
|
!
|
|
CASE ('Ba')
|
|
alpha=275.0_DP
|
|
C6=5727.0_DP
|
|
R0=4.77_DP
|
|
!
|
|
CASE ('Hf')
|
|
alpha=99.52_DP
|
|
C6=1274.8_DP
|
|
R0=4.21_DP
|
|
!
|
|
CASE ('Ta')
|
|
alpha=82.53_DP
|
|
C6=1019.92_DP
|
|
R0=4.15_DP
|
|
!
|
|
CASE ('W')
|
|
alpha=71.041_DP
|
|
C6=847.93_DP
|
|
R0=4.08_DP
|
|
!
|
|
CASE ('Re')
|
|
alpha=63.04_DP
|
|
C6=710.2_DP
|
|
R0=4.02_DP
|
|
!
|
|
CASE ('Os')
|
|
alpha=55.055_DP
|
|
C6=596.67_DP
|
|
R0=3.84_DP
|
|
!
|
|
CASE ('Ir')
|
|
alpha=42.51_DP
|
|
C6=359.1_DP
|
|
R0=4.00_DP
|
|
!
|
|
CASE ('Pt')
|
|
alpha=39.68_DP
|
|
C6=347.1_DP
|
|
R0=3.92_DP
|
|
!
|
|
CASE ('Au')
|
|
alpha=36.5_DP
|
|
C6=298.0_DP
|
|
R0=3.86_DP
|
|
!
|
|
CASE ('Hg')
|
|
alpha=33.9_DP
|
|
C6=392.0_DP
|
|
R0=3.98_DP
|
|
!
|
|
CASE ('Tl')
|
|
alpha=69.92_DP
|
|
C6=717.44_DP
|
|
R0=3.91_DP
|
|
!
|
|
CASE ('Pb')
|
|
alpha=61.8_DP
|
|
C6=697.0_DP
|
|
R0=4.31_DP
|
|
!
|
|
CASE ('Bi')
|
|
alpha=49.02_DP
|
|
C6=571.0_DP
|
|
R0=4.32_DP
|
|
!
|
|
CASE ('Po')
|
|
alpha=45.013_DP
|
|
C6=530.92_DP
|
|
R0=4.097_DP
|
|
!
|
|
CASE ('At')
|
|
alpha=38.93_DP
|
|
C6=457.53_DP
|
|
R0=4.07_DP
|
|
!
|
|
CASE ('Rn')
|
|
alpha=33.54_DP
|
|
C6=390.63_DP
|
|
R0=4.23_DP
|
|
!
|
|
CASE DEFAULT
|
|
!
|
|
CALL errore('tsvdw','Reference free atom parameters not available for requested atom type...',1)
|
|
!
|
|
END SELECT
|
|
!
|
|
RETURN
|
|
!
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
END SUBROUTINE GetVdWParam
|
|
!--------------------------------------------------------------------------------------------------------------
|
|
!
|
|
!
|
|
END MODULE tsvdw_module
|
|
|