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,"")') !DO ia=1,nat ! WRITE(stdout,'(5X,I5,F25.12)') ia,veff(ia) !END DO !WRITE(stdout,'(3X,"")') !! !WRITE(stdout,'(3X,"")') !DO ia=1,nat ! WRITE(stdout,'(5X,I5,3F25.12)') ia,FtsvdW(:,ia) !END DO !WRITE(stdout,'(3X,"")') !! !WRITE(stdout,'(3X,"")') !DO i=1,3 ! WRITE(stdout,'(5X,3F25.12)') HtsvdW(i,:) !END DO !WRITE(stdout,'(3X,"")') !! 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