mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-21 06:25:15 -04:00
ccECP (#3940)
This commit is contained in:
Juerg Hutter
GitHub
parent
413afc8300
commit
71d5d6c00b
11 changed files with 1742 additions and 34 deletions
1001
data/ECP_POTENTIALS
1001
data/ECP_POTENTIALS
File diff suppressed because it is too large
Load diff
|
|
@ -934,6 +934,7 @@ list(
|
|||
common/eigenvalueproblems.F
|
||||
common/fparser.F
|
||||
common/gamma.F
|
||||
common/gfun.F
|
||||
common/kahan_sum.F
|
||||
common/lebedev.F
|
||||
common/linear_systems.F
|
||||
|
|
|
|||
|
|
@ -29,9 +29,11 @@ MODULE ai_overlap_ppl
|
|||
os_3center
|
||||
USE ai_overlap_debug, ONLY: init_os_overlap2,&
|
||||
os_overlap2
|
||||
USE gamma, ONLY: fgamma => fgamma_0
|
||||
USE gfun, ONLY: gfun_values
|
||||
USE kinds, ONLY: dp
|
||||
USE mathconstants, ONLY: fac,&
|
||||
pi
|
||||
USE mathconstants, ONLY: pi
|
||||
USE mathlib, ONLY: binomial
|
||||
USE orbital_pointers, ONLY: indco,&
|
||||
ncoset
|
||||
#include "../base/base_uses.f90"
|
||||
|
|
@ -44,7 +46,7 @@ MODULE ai_overlap_ppl
|
|||
|
||||
! *** Public subroutines ***
|
||||
|
||||
PUBLIC :: ppl_integral, ppl_integral_ri
|
||||
PUBLIC :: ecploc_integral, ppl_integral, ppl_integral_ri
|
||||
|
||||
CONTAINS
|
||||
|
||||
|
|
@ -182,6 +184,140 @@ CONTAINS
|
|||
DEALLOCATE (auxint)
|
||||
|
||||
END SUBROUTINE ppl_integral
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Calculation of three-center potential integrals <a|V(r)|b> over
|
||||
!> Cartesian Gaussian functions for the local part of ECP
|
||||
!> pseudopotential. Multiple terms C1-4 are possible.
|
||||
!>
|
||||
!> <a|V(ecploc)|b> = <a| C1/r*exp(-a1*r**2) + C2*exp(-a2*r**2) + C3*r*exp(-a3*r**2) +
|
||||
!> C4*r**2*exp(-a4*r**2)|b>
|
||||
!>
|
||||
!> \param la_max_set ...
|
||||
!> \param la_min_set ...
|
||||
!> \param npgfa ...
|
||||
!> \param rpgfa ...
|
||||
!> \param zeta ...
|
||||
!> \param lb_max_set ...
|
||||
!> \param lb_min_set ...
|
||||
!> \param npgfb ...
|
||||
!> \param rpgfb ...
|
||||
!> \param zetb ...
|
||||
!> \param nexp_ppl ...
|
||||
!> \param alpha_ppl ...
|
||||
!> \param nct_ppl ...
|
||||
!> \param cexp_ppl ...
|
||||
!> \param rpgfc ...
|
||||
!> \param rab ...
|
||||
!> \param dab ...
|
||||
!> \param rac ...
|
||||
!> \param dac ...
|
||||
!> \param rbc ...
|
||||
!> \param dbc ...
|
||||
!> \param vab ...
|
||||
!> \param s ...
|
||||
!> \param pab ...
|
||||
!> \param force_a ...
|
||||
!> \param force_b ...
|
||||
!> \param fs ...
|
||||
!> \param hab2 The derivative of the ppl integrals according to the weighting factors deltaR
|
||||
!> \param hab2_work ...
|
||||
!> \param deltaR Weighting factors for the derivatives wrt. nuclear positions
|
||||
!> \param iatom ...
|
||||
!> \param jatom ...
|
||||
!> \param katom ...
|
||||
!> \date 2025
|
||||
!> \author Juerg Hutter
|
||||
!> \version 1.0
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE ecploc_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
|
||||
lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
|
||||
nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, &
|
||||
rab, dab, rac, dac, rbc, dbc, vab, s, pab, &
|
||||
force_a, force_b, fs, hab2, hab2_work, &
|
||||
deltaR, iatom, jatom, katom)
|
||||
INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
|
||||
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
|
||||
INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
|
||||
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
|
||||
INTEGER, INTENT(IN) :: nexp_ppl
|
||||
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha_ppl
|
||||
INTEGER, DIMENSION(:), INTENT(IN) :: nct_ppl
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: cexp_ppl
|
||||
REAL(KIND=dp), INTENT(IN) :: rpgfc
|
||||
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
|
||||
REAL(KIND=dp), INTENT(IN) :: dab
|
||||
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
|
||||
REAL(KIND=dp), INTENT(IN) :: dac
|
||||
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rbc
|
||||
REAL(KIND=dp), INTENT(IN) :: dbc
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
|
||||
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
|
||||
OPTIONAL :: pab
|
||||
REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
|
||||
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
|
||||
OPTIONAL :: fs, hab2, hab2_work
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
|
||||
OPTIONAL :: deltaR
|
||||
INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
|
||||
|
||||
INTEGER :: iexp, ij, ipgf, jpgf, mmax, nexp
|
||||
REAL(KIND=dp) :: rho, sab, t, zetc
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: auxint
|
||||
REAL(KIND=dp), DIMENSION(3) :: pci
|
||||
|
||||
IF (PRESENT(pab)) THEN
|
||||
CPASSERT(PRESENT(force_a))
|
||||
CPASSERT(PRESENT(force_b))
|
||||
CPASSERT(PRESENT(fs))
|
||||
mmax = la_max_set + lb_max_set + 2
|
||||
force_a(:) = 0.0_dp
|
||||
force_b(:) = 0.0_dp
|
||||
ELSE IF (PRESENT(hab2)) THEN
|
||||
mmax = la_max_set + lb_max_set + 2
|
||||
ELSE
|
||||
mmax = la_max_set + lb_max_set
|
||||
END IF
|
||||
|
||||
ALLOCATE (auxint(0:mmax, npgfa*npgfb))
|
||||
auxint = 0._dp
|
||||
|
||||
! *** Calculate auxiliary integrals ***
|
||||
|
||||
DO ipgf = 1, npgfa
|
||||
! *** Screening ***
|
||||
IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
|
||||
DO jpgf = 1, npgfb
|
||||
! *** Screening ***
|
||||
IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
|
||||
(rpgfa(ipgf) + rpgfb(jpgf) < dab)) CYCLE
|
||||
ij = (ipgf - 1)*npgfb + jpgf
|
||||
rho = zeta(ipgf) + zetb(jpgf)
|
||||
pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
|
||||
sab = EXP(-(zeta(ipgf)*zetb(jpgf)/rho*dab*dab))
|
||||
t = rho*SUM(pci(:)*pci(:))
|
||||
|
||||
DO iexp = 1, nexp_ppl
|
||||
nexp = nct_ppl(iexp)
|
||||
zetc = alpha_ppl(iexp)
|
||||
CALL ecploc_aux(auxint(0:mmax, ij), mmax, t, rho, nexp, cexp_ppl(1, iexp), zetc)
|
||||
END DO
|
||||
|
||||
auxint(0:mmax, ij) = sab*auxint(0:mmax, ij)
|
||||
|
||||
END DO
|
||||
END DO
|
||||
|
||||
CALL os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
|
||||
lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
|
||||
rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
|
||||
vab2=hab2, vab2_work=hab2_work, &
|
||||
deltaR=deltaR, iatom=iatom, jatom=jatom, katom=katom)
|
||||
|
||||
DEALLOCATE (auxint)
|
||||
|
||||
END SUBROUTINE ecploc_integral
|
||||
! **************************************************************************************************
|
||||
!> \brief Calculation of two-center overlap integrals <a|c> over
|
||||
!> Cartesian Gaussian functions for the local part of the Goedecker
|
||||
|
|
@ -453,29 +589,135 @@ CONTAINS
|
|||
DO j = 0, MIN(i, pmax)
|
||||
kp = j
|
||||
ke = i - j
|
||||
auxint(i) = auxint(i) + expder(ke)*polder(kp)*choose(i, j)
|
||||
auxint(i) = auxint(i) + expder(ke)*polder(kp)*binomial(i, j)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
END SUBROUTINE ppl_aux
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param n ...
|
||||
!> \param k ...
|
||||
!> \return ...
|
||||
!> \param auxint ...
|
||||
!> \param mmax ...
|
||||
!> \param t ...
|
||||
!> \param rho ...
|
||||
!> \param nexp ...
|
||||
!> \param cexp ...
|
||||
!> \param zetc ...
|
||||
! **************************************************************************************************
|
||||
FUNCTION choose(n, k)
|
||||
SUBROUTINE ecploc_aux(auxint, mmax, t, rho, nexp, cexp, zetc)
|
||||
INTEGER, INTENT(IN) :: mmax
|
||||
REAL(KIND=dp), DIMENSION(0:mmax) :: auxint
|
||||
REAL(KIND=dp), INTENT(IN) :: t, rho
|
||||
INTEGER, INTENT(IN) :: nexp
|
||||
REAL(KIND=dp), INTENT(IN) :: cexp, zetc
|
||||
|
||||
INTEGER, INTENT(IN) :: n, k
|
||||
REAL(KIND=dp) :: choose
|
||||
INTEGER :: i, j, ke, kf
|
||||
REAL(KIND=dp) :: c0, c1, cc, cval, fa, fr, q, ts
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: expder, fdiff, funder, gfund
|
||||
|
||||
IF (n >= k) THEN
|
||||
choose = REAL(NINT(fac(n)/(fac(k)*fac(n - k))), KIND=dp)
|
||||
ELSE
|
||||
choose = 0.0_dp
|
||||
END IF
|
||||
q = rho + zetc
|
||||
fa = zetc/q
|
||||
fr = rho/q
|
||||
!
|
||||
ALLOCATE (expder(0:mmax), funder(0:mmax + 1))
|
||||
!
|
||||
SELECT CASE (nexp)
|
||||
CASE (0)
|
||||
cval = 2.0_dp*cexp/SQRT(q)*pi**1.5_dp*EXP(-t*fa)
|
||||
expder(0) = cval
|
||||
DO i = 1, mmax
|
||||
expder(i) = fa*expder(i - 1)
|
||||
END DO
|
||||
ts = fr*t
|
||||
ALLOCATE (gfund(0:mmax))
|
||||
CALL gfun_values(mmax, ts, gfund)
|
||||
|
||||
END FUNCTION choose
|
||||
funder(0) = gfund(0)
|
||||
DO i = 1, mmax
|
||||
funder(i) = 0.0_dp
|
||||
DO j = 0, i
|
||||
funder(i) = funder(i) + (-1)**j*binomial(i, j)*gfund(j)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
DEALLOCATE (gfund)
|
||||
DO i = 1, mmax
|
||||
funder(i) = fr**i*funder(i)
|
||||
END DO
|
||||
DO i = 0, mmax
|
||||
DO j = 0, i
|
||||
kf = j
|
||||
ke = i - j
|
||||
auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
|
||||
END DO
|
||||
END DO
|
||||
CASE (1)
|
||||
cval = cexp*2._dp*pi/q*EXP(-t*fa)
|
||||
expder(0) = cval
|
||||
DO i = 1, mmax
|
||||
expder(i) = fa*expder(i - 1)
|
||||
END DO
|
||||
ts = fr*t
|
||||
CALL fgamma(mmax, ts, funder)
|
||||
DO i = 1, mmax
|
||||
funder(i) = fr**i*funder(i)
|
||||
END DO
|
||||
DO i = 0, mmax
|
||||
DO j = 0, i
|
||||
kf = j
|
||||
ke = i - j
|
||||
auxint(i) = auxint(i) + expder(ke)*funder(kf)*binomial(i, j)
|
||||
END DO
|
||||
END DO
|
||||
CASE (2)
|
||||
cval = cexp*(pi/q)**1.5_dp*EXP(-t*fa)
|
||||
expder(0) = cval
|
||||
DO i = 1, mmax
|
||||
expder(i) = fa*expder(i - 1)
|
||||
END DO
|
||||
auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
|
||||
CASE (3)
|
||||
cval = 2.*pi*cexp/q**2*EXP(-t*fa)
|
||||
expder(0) = cval
|
||||
DO i = 1, mmax
|
||||
expder(i) = fa*expder(i - 1)
|
||||
END DO
|
||||
ts = fr*t
|
||||
CALL fgamma(mmax + 1, ts, funder)
|
||||
ALLOCATE (fdiff(0:mmax))
|
||||
fdiff(0) = (1.0_dp + ts)*funder(0) - ts*funder(1)
|
||||
DO i = 1, mmax
|
||||
fdiff(i) = fr**i*(-i*funder(i - 1) + (1.0_dp + ts)*funder(i) &
|
||||
+ i*funder(i) - ts*funder(i + 1))
|
||||
END DO
|
||||
DO i = 0, mmax
|
||||
DO j = 0, i
|
||||
kf = j
|
||||
ke = i - j
|
||||
auxint(i) = auxint(i) + expder(ke)*fdiff(kf)*binomial(i, j)
|
||||
END DO
|
||||
END DO
|
||||
DEALLOCATE (fdiff)
|
||||
CASE (4)
|
||||
cval = cexp/(4._dp*q**2)*(pi/q)**1.5_dp*EXP(-t*fa)
|
||||
expder(0) = cval
|
||||
DO i = 1, mmax
|
||||
expder(i) = fa*expder(i - 1)
|
||||
END DO
|
||||
c0 = 4._dp*rho/fa
|
||||
c1 = 6._dp*q + 4._dp*rho*t
|
||||
DO i = 0, mmax
|
||||
cc = -i*c0 + c1
|
||||
expder(i) = cc*expder(i)
|
||||
END DO
|
||||
auxint(0:mmax) = auxint(0:mmax) + expder(0:mmax)
|
||||
CASE DEFAULT
|
||||
CPABORT("nexp out of range [1..4]")
|
||||
END SELECT
|
||||
!
|
||||
DEALLOCATE (expder, funder)
|
||||
|
||||
END SUBROUTINE ecploc_aux
|
||||
! **************************************************************************************************
|
||||
|
||||
END MODULE ai_overlap_ppl
|
||||
|
|
|
|||
|
|
@ -2980,6 +2980,7 @@ CONTAINS
|
|||
SELECT CASE (nel)
|
||||
CASE DEFAULT
|
||||
CPABORT("Unknown Core State")
|
||||
CASE (0)
|
||||
CASE (2)
|
||||
potential%econf(0:3) = potential%econf(0:3) - ptable(2)%e_conv(0:3)
|
||||
CASE (10)
|
||||
|
|
|
|||
305
src/common/gfun.F
Normal file
305
src/common/gfun.F
Normal file
|
|
@ -0,0 +1,305 @@
|
|||
!--------------------------------------------------------------------------------------------------!
|
||||
! CP2K: A general program to perform molecular dynamics simulations !
|
||||
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
|
||||
! !
|
||||
! SPDX-License-Identifier: GPL-2.0-or-later !
|
||||
!--------------------------------------------------------------------------------------------------!
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Calculation of the G function G_n(t) for 1/R^2 operators
|
||||
!
|
||||
! (1) J. 0. Jensen, A. H. Cameri, C. P. Vlahacos, D. Zeroka, H. F. Hameka, C. N. Merrow,
|
||||
! Evaluation of one-electron integrals for arbitrary operators V(r) over cartesian Gaussians:
|
||||
! Application to inverse-square distance and Yukawa operators.
|
||||
! J. Comput. Chem. 14(8), 986 (1993).
|
||||
! doi: 10.1002/jcc.540140814
|
||||
! (2) B. Gao, A. J. Thorvaldsen, K. Ruud,
|
||||
! GEN1INT: A unified procedure for the evaluation of one-electron integrals over Gaussian
|
||||
! basis functions and their geometric derivatives.
|
||||
! Int. J. Quantum Chem. 111(4), 858 (2011).
|
||||
! doi: 10.1002/qua.22886
|
||||
! (3) libgrpp : specfun_gfun.c
|
||||
! (4) William Cody, Kathleen Paciorek, Henry Thacher,
|
||||
! Chebyshev Approximations for Dawson's Integral,
|
||||
! Mathematics of Computation,
|
||||
! Volume 24, Number 109, January 1970, pages 171-178.
|
||||
!
|
||||
!> \author JHU
|
||||
! **************************************************************************************************
|
||||
MODULE gfun
|
||||
|
||||
USE kinds, ONLY: dp
|
||||
#include "../base/base_uses.f90"
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
PRIVATE
|
||||
|
||||
! *** Global parameters ***
|
||||
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gfun'
|
||||
|
||||
PUBLIC :: gfun_values
|
||||
|
||||
CONTAINS
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param nmax ...
|
||||
!> \param t ...
|
||||
!> \param g ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE gfun_values(nmax, t, g)
|
||||
|
||||
INTEGER, INTENT(IN) :: nmax
|
||||
REAL(KIND=dp), INTENT(IN) :: t
|
||||
REAL(KIND=dp), DIMENSION(0:nmax), INTENT(OUT) :: g
|
||||
|
||||
INTEGER :: i
|
||||
REAL(KIND=dp) :: st
|
||||
|
||||
g = 0.0_dp
|
||||
|
||||
IF (t <= 12.0_dp) THEN
|
||||
! downward recursion
|
||||
g(nmax) = gfun_taylor(nmax, t);
|
||||
DO i = nmax, 1, -1
|
||||
g(i - 1) = (1.0_dp - 2.0_dp*t*g(i))/(2.0_dp*i - 1.0_dp)
|
||||
END DO
|
||||
ELSE
|
||||
! upward recursion
|
||||
st = SQRT(t)
|
||||
g(0) = daw(st)/st
|
||||
DO i = 0, nmax - 1
|
||||
g(i + 1) = (1.0_dp - (2.0_dp*i + 1.0_dp)*g(i))/(2.0_dp*t)
|
||||
END DO
|
||||
END IF
|
||||
|
||||
END SUBROUTINE gfun_values
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param n ...
|
||||
!> \param x ...
|
||||
!> \return ...
|
||||
! **************************************************************************************************
|
||||
FUNCTION gfun_taylor(n, x) RESULT(g)
|
||||
INTEGER, INTENT(IN) :: n
|
||||
REAL(KIND=dp), INTENT(IN) :: x
|
||||
REAL(KIND=dp) :: g
|
||||
|
||||
REAL(KIND=dp), PARAMETER :: eps = 1.E-15_dp
|
||||
|
||||
INTEGER :: k
|
||||
REAL(KIND=dp) :: ex, gk, tkk
|
||||
|
||||
ex = EXP(-x)
|
||||
tkk = 1.0_dp
|
||||
g = ex/REAL(2*n + 1, KIND=dp)
|
||||
DO k = 1, 100
|
||||
tkk = tkk*x/REAL(k, KIND=dp)
|
||||
gk = ex*tkk/REAL(2*n + 2*k + 1, KIND=dp)
|
||||
g = g + gk
|
||||
IF (gk < eps) EXIT
|
||||
END DO
|
||||
IF (gk > eps) THEN
|
||||
CPWARN("gfun_taylor did not converge")
|
||||
END IF
|
||||
|
||||
END FUNCTION gfun_taylor
|
||||
|
||||
!*****************************************************************************80
|
||||
!
|
||||
! DAW evaluates Dawson's integral function.
|
||||
!
|
||||
! Discussion:
|
||||
!
|
||||
! This routine evaluates Dawson's integral,
|
||||
!
|
||||
! F(x) = exp ( - x * x ) * Integral ( 0 <= t <= x ) exp ( t * t ) dt
|
||||
!
|
||||
! for a real argument x.
|
||||
!
|
||||
! Licensing:
|
||||
!
|
||||
! This code is distributed under the GNU LGPL license.
|
||||
!
|
||||
! Modified:
|
||||
!
|
||||
! 03 April 2007
|
||||
!
|
||||
! Author:
|
||||
!
|
||||
! Original FORTRAN77 version by William Cody.
|
||||
! FORTRAN90 version by John Burkardt.
|
||||
!
|
||||
! Reference:
|
||||
!
|
||||
! William Cody, Kathleen Paciorek, Henry Thacher,
|
||||
! Chebyshev Approximations for Dawson's Integral,
|
||||
! Mathematics of Computation,
|
||||
! Volume 24, Number 109, January 1970, pages 171-178.
|
||||
!
|
||||
! Parameters:
|
||||
!
|
||||
! Input, real ( kind = dp ) XX, the argument of the function.
|
||||
!
|
||||
! Output, real ( kind = dp ) DAW, the value of the function.
|
||||
!
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param xx ...
|
||||
!> \return ...
|
||||
! **************************************************************************************************
|
||||
FUNCTION daw(xx)
|
||||
|
||||
REAL(kind=dp) :: xx, daw
|
||||
|
||||
INTEGER :: i
|
||||
REAL(kind=dp) :: frac, one225, p1(10), p2(10), p3(10), &
|
||||
p4(10), q1(10), q2(9), q3(9), q4(9), &
|
||||
six25, sump, sumq, two5, w2, x, &
|
||||
xlarge, xmax, xsmall, y
|
||||
|
||||
!
|
||||
! Mathematical constants.
|
||||
!
|
||||
DATA six25/6.25D+00/
|
||||
DATA one225/12.25d0/
|
||||
DATA two5/25.0d0/
|
||||
!
|
||||
! Machine-dependent constants
|
||||
!
|
||||
DATA xsmall/1.05d-08/
|
||||
DATA xlarge/9.49d+07/
|
||||
DATA xmax/2.24d+307/
|
||||
!
|
||||
! Coefficients for R(9,9) approximation for |x| < 2.5
|
||||
!
|
||||
DATA p1/-2.69020398788704782410d-12, 4.18572065374337710778d-10, &
|
||||
-1.34848304455939419963d-08, 9.28264872583444852976d-07, &
|
||||
-1.23877783329049120592d-05, 4.07205792429155826266d-04, &
|
||||
-2.84388121441008500446d-03, 4.70139022887204722217d-02, &
|
||||
-1.38868086253931995101d-01, 1.00000000000000000004d+00/
|
||||
DATA q1/1.71257170854690554214d-10, 1.19266846372297253797d-08, &
|
||||
4.32287827678631772231d-07, 1.03867633767414421898d-05, &
|
||||
1.78910965284246249340d-04, 2.26061077235076703171d-03, &
|
||||
2.07422774641447644725d-02, 1.32212955897210128811d-01, &
|
||||
5.27798580412734677256d-01, 1.00000000000000000000d+00/
|
||||
!
|
||||
! Coefficients for R(9,9) approximation in J-fraction form
|
||||
! for x in [2.5, 3.5)
|
||||
!
|
||||
DATA p2/-1.70953804700855494930d+00, -3.79258977271042880786d+01, &
|
||||
2.61935631268825992835d+01, 1.25808703738951251885d+01, &
|
||||
-2.27571829525075891337d+01, 4.56604250725163310122d+00, &
|
||||
-7.33080089896402870750d+00, 4.65842087940015295573d+01, &
|
||||
-1.73717177843672791149d+01, 5.00260183622027967838d-01/
|
||||
DATA q2/1.82180093313514478378d+00, 1.10067081034515532891d+03, &
|
||||
-7.08465686676573000364d+00, 4.53642111102577727153d+02, &
|
||||
4.06209742218935689922d+01, 3.02890110610122663923d+02, &
|
||||
1.70641269745236227356d+02, 9.51190923960381458747d+02, &
|
||||
2.06522691539642105009d-01/
|
||||
!
|
||||
! Coefficients for R(9,9) approximation in J-fraction form
|
||||
! for x in [3.5, 5.0]
|
||||
!
|
||||
DATA p3/-4.55169503255094815112d+00, -1.86647123338493852582d+01, &
|
||||
-7.36315669126830526754d+00, -6.68407240337696756838d+01, &
|
||||
4.84507265081491452130d+01, 2.69790586735467649969d+01, &
|
||||
-3.35044149820592449072d+01, 7.50964459838919612289d+00, &
|
||||
-1.48432341823343965307d+00, 4.99999810924858824981d-01/
|
||||
DATA q3/4.47820908025971749852d+01, 9.98607198039452081913d+01, &
|
||||
1.40238373126149385228d+01, 3.48817758822286353588d+03, &
|
||||
-9.18871385293215873406d+00, 1.24018500009917163023d+03, &
|
||||
-6.88024952504512254535d+01, -2.31251575385145143070d+00, &
|
||||
2.50041492369922381761d-01/
|
||||
!
|
||||
! Coefficients for R(9,9) approximation in J-fraction form
|
||||
! for 5.0 < |x|.
|
||||
!
|
||||
DATA p4/-8.11753647558432685797d+00, -3.84043882477454453430d+01, &
|
||||
-2.23787669028751886675d+01, -2.88301992467056105854d+01, &
|
||||
-5.99085540418222002197d+00, -1.13867365736066102577d+01, &
|
||||
-6.52828727526980741590d+00, -4.50002293000355585708d+00, &
|
||||
-2.50000000088955834952d+00, 5.00000000000000488400d-01/
|
||||
DATA q4/2.69382300417238816428d+02, 5.04198958742465752861d+01, &
|
||||
6.11539671480115846173d+01, 2.08210246935564547889d+02, &
|
||||
1.97325365692316183531d+01, -1.22097010558934838708d+01, &
|
||||
-6.99732735041547247161d+00, -2.49999970104184464568d+00, &
|
||||
7.49999999999027092188d-01/
|
||||
|
||||
x = xx
|
||||
|
||||
IF (xlarge < ABS(x)) THEN
|
||||
|
||||
IF (ABS(x) <= xmax) THEN
|
||||
daw = 0.5D+00/x
|
||||
ELSE
|
||||
daw = 0.0D+00
|
||||
END IF
|
||||
|
||||
ELSE IF (ABS(x) < xsmall) THEN
|
||||
|
||||
daw = x
|
||||
|
||||
ELSE
|
||||
|
||||
y = x*x
|
||||
!
|
||||
! ABS(X) < 2.5.
|
||||
!
|
||||
IF (y < six25) THEN
|
||||
|
||||
sump = p1(1)
|
||||
sumq = q1(1)
|
||||
DO i = 2, 10
|
||||
sump = sump*y + p1(i)
|
||||
sumq = sumq*y + q1(i)
|
||||
END DO
|
||||
|
||||
daw = x*sump/sumq
|
||||
!
|
||||
! 2.5 <= ABS(X) < 3.5.
|
||||
!
|
||||
ELSE IF (y < one225) THEN
|
||||
|
||||
frac = 0.0D+00
|
||||
DO i = 1, 9
|
||||
frac = q2(i)/(p2(i) + y + frac)
|
||||
END DO
|
||||
|
||||
daw = (p2(10) + frac)/x
|
||||
!
|
||||
! 3.5 <= ABS(X) < 5.0.
|
||||
!
|
||||
ELSE IF (y < two5) THEN
|
||||
|
||||
frac = 0.0D+00
|
||||
DO i = 1, 9
|
||||
frac = q3(i)/(p3(i) + y + frac)
|
||||
END DO
|
||||
|
||||
daw = (p3(10) + frac)/x
|
||||
|
||||
ELSE
|
||||
!
|
||||
! 5.0 <= ABS(X) <= XLARGE.
|
||||
!
|
||||
w2 = 1.0D+00/x/x
|
||||
|
||||
frac = 0.0D+00
|
||||
DO i = 1, 9
|
||||
frac = q4(i)/(p4(i) + y + frac)
|
||||
END DO
|
||||
frac = p4(10) + frac
|
||||
|
||||
daw = (0.5D+00 + 0.5D+00*w2*frac)/x
|
||||
|
||||
END IF
|
||||
|
||||
END IF
|
||||
|
||||
END FUNCTION daw
|
||||
|
||||
END MODULE gfun
|
||||
|
|
@ -17,7 +17,8 @@
|
|||
! **************************************************************************************************
|
||||
MODULE core_ppl
|
||||
|
||||
USE ai_overlap_ppl, ONLY: ppl_integral,&
|
||||
USE ai_overlap_ppl, ONLY: ecploc_integral,&
|
||||
ppl_integral,&
|
||||
ppl_integral_ri
|
||||
USE atomic_kind_types, ONLY: atomic_kind_type,&
|
||||
get_atomic_kind_set
|
||||
|
|
@ -130,11 +131,12 @@ CONTAINS
|
|||
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
|
||||
INTEGER, DIMENSION(nexp_max) :: nct_ppl
|
||||
LOGICAL :: do_dR, doat, dokp, ecp_local, &
|
||||
ecp_semi_local, found, lpotextended, &
|
||||
ecp_semi_local, found, &
|
||||
libgrpp_local = .FALSE., lpotextended, &
|
||||
only_gaussians
|
||||
REAL(KIND=dp) :: alpha, atk0, atk1, dab, dac, dbc, f0, &
|
||||
ppl_radius
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qab, work
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hab2_w, ppl_fwork, ppl_work
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: hab, pab
|
||||
REAL(KIND=dp), ALLOCATABLE, &
|
||||
|
|
@ -232,13 +234,13 @@ CONTAINS
|
|||
!$OMP SHARED (ap_iterator, basis_set_list, calculate_forces, use_virial, &
|
||||
!$OMP matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
|
||||
!$OMP sab_orb, sac_ppl, nthread, ncoset, nkind, cell_to_index, &
|
||||
!$OMP ldsab, maxnset, maxder, do_dR, deltaR, doat, &
|
||||
!$OMP ldsab, maxnset, maxder, do_dR, deltaR, doat, libgrpp_local, &
|
||||
!$OMP maxlgto, nder, maxco, dokp, locks, natom) &
|
||||
!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, basis_set_a, basis_set_b, &
|
||||
!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, &
|
||||
!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
|
||||
!$OMP zetb, dab, irow, icol, h_block, found, iset, ncoa, &
|
||||
!$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, hab2, hab2_w, qab, &
|
||||
!$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, hab2, hab2_w, &
|
||||
!$OMP atk0, atk1, h1_1block, h1_2block, h1_3block, kkind, nseta, &
|
||||
!$OMP gth_potential, sgp_potential, alpha, cexp_ppl, lpotextended, &
|
||||
!$OMP ppl_radius, nexp_lpot, nexp_ppl, alpha_ppl, alpha_lpot, nct_ppl, &
|
||||
|
|
@ -278,7 +280,6 @@ CONTAINS
|
|||
IF (do_dR) THEN
|
||||
ALLOCATE (hab2(ldsab, ldsab, 4, maxnset, maxnset))
|
||||
ALLOCATE (hab2_w(ldsab, ldsab, 6))
|
||||
ALLOCATE (qab(ldsab, ldsab))
|
||||
ALLOCATE (ppl_fwork(ldai, ldai, maxder))
|
||||
END IF
|
||||
|
||||
|
|
@ -425,8 +426,7 @@ CONTAINS
|
|||
nct_ppl(1:nloc) = nrloc(1:nloc)
|
||||
alpha_ppl(1:nloc) = bloc(1:nloc)
|
||||
cval_ppl(1, 1:nloc) = aloc(1:nloc)
|
||||
only_gaussians = ALL(nct_ppl(1:nloc) == 2)
|
||||
IF (only_gaussians) nct_ppl(1:nloc) = nct_ppl(1:nloc) - 1
|
||||
only_gaussians = .FALSE.
|
||||
ELSE
|
||||
CALL get_potential(potential=sgp_potential, n_local=n_local, a_local=a_local, c_local=c_local)
|
||||
nexp_ppl = n_local
|
||||
|
|
@ -488,8 +488,7 @@ CONTAINS
|
|||
rab, dab, rac, dac, rbc, dbc, &
|
||||
hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
|
||||
force_a, force_b, ppl_fwork)
|
||||
ELSE
|
||||
|
||||
ELSEIF (libgrpp_local) THEN
|
||||
!$OMP CRITICAL(type1)
|
||||
CALL libgrpp_local_forces_ref(la_max(iset), la_min(iset), npgfa(iset), &
|
||||
rpgfa(:, iset), zeta(:, iset), &
|
||||
|
|
@ -500,6 +499,16 @@ CONTAINS
|
|||
hab(:, :, iset, jset), pab(:, :, iset, jset), &
|
||||
force_a, force_b)
|
||||
!$OMP END CRITICAL(type1)
|
||||
ELSE
|
||||
CALL ecploc_integral( &
|
||||
la_max(iset), la_min(iset), npgfa(iset), &
|
||||
rpgfa(:, iset), zeta(:, iset), &
|
||||
lb_max(jset), lb_min(jset), npgfb(jset), &
|
||||
rpgfb(:, jset), zetb(:, jset), &
|
||||
nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
|
||||
rab, dab, rac, dac, rbc, dbc, &
|
||||
hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
|
||||
force_a, force_b, ppl_fwork)
|
||||
END IF
|
||||
|
||||
IF (ecp_semi_local) THEN
|
||||
|
|
@ -565,7 +574,7 @@ CONTAINS
|
|||
nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
|
||||
rab, dab, rac, dac, rbc, dbc, hab(:, :, iset, jset), ppl_work)
|
||||
|
||||
ELSE
|
||||
ELSEIF (libgrpp_local) THEN
|
||||
!If the local part of the potential is more complex, we need libgrpp
|
||||
!$OMP CRITICAL(type1)
|
||||
CALL libgrpp_local_integrals(la_max(iset), la_min(iset), npgfa(iset), &
|
||||
|
|
@ -576,6 +585,14 @@ CONTAINS
|
|||
ppl_radius, rab, dab, rac, dac, dbc, &
|
||||
hab(:, :, iset, jset))
|
||||
!$OMP END CRITICAL(type1)
|
||||
ELSE
|
||||
CALL ecploc_integral( &
|
||||
la_max(iset), la_min(iset), npgfa(iset), &
|
||||
rpgfa(:, iset), zeta(:, iset), &
|
||||
lb_max(jset), lb_min(jset), npgfb(jset), &
|
||||
rpgfb(:, jset), zetb(:, jset), &
|
||||
nexp_ppl, alpha_ppl, nct_ppl, cval_ppl, ppl_radius, &
|
||||
rab, dab, rac, dac, rbc, dbc, hab(:, :, iset, jset), ppl_work)
|
||||
END IF
|
||||
|
||||
IF (ecp_semi_local) THEN
|
||||
|
|
@ -684,7 +701,7 @@ CONTAINS
|
|||
END DO
|
||||
END DO
|
||||
END IF
|
||||
IF (do_dR) DEALLOCATE (qab, hab2, ppl_fwork, hab2_w)
|
||||
IF (do_dR) DEALLOCATE (hab2, ppl_fwork, hab2_w)
|
||||
END DO ! slot
|
||||
|
||||
DEALLOCATE (hab, work, ppl_work)
|
||||
|
|
|
|||
|
|
@ -201,3 +201,23 @@ I D
|
|||
2 4.1442856 71.9688361
|
||||
2 0.9377235 9.3630657
|
||||
END LANL2DZ_ECP
|
||||
ccECP
|
||||
H nelec 0
|
||||
H ul
|
||||
1 21.24359508259891 1.00000000000000
|
||||
3 21.24359508259891 21.24359508259891
|
||||
2 21.77696655044365 -10.85192405303825
|
||||
H S
|
||||
2 1.000000000000000 0.00000000000000
|
||||
Cl nelec 10
|
||||
Cl ul
|
||||
1 7.944352 7.000000
|
||||
3 12.801261 55.610463
|
||||
2 6.296744 -22.860784
|
||||
Cl S
|
||||
2 17.908432 15.839234
|
||||
2 4.159880 44.469504
|
||||
Cl P
|
||||
2 7.931763 8.321946
|
||||
2 3.610412 24.044745
|
||||
END ccECP
|
||||
|
|
|
|||
55
tests/QS/regtest-ecp/HCl_ccECP.inp
Normal file
55
tests/QS/regtest-ecp/HCl_ccECP.inp
Normal file
|
|
@ -0,0 +1,55 @@
|
|||
&GLOBAL
|
||||
PRINT_LEVEL LOW
|
||||
PROJECT HCl
|
||||
RUN_TYPE ENERGY
|
||||
&END GLOBAL
|
||||
|
||||
&FORCE_EVAL
|
||||
&DFT
|
||||
POTENTIAL_FILE_NAME ./ECP_BASIS_POT
|
||||
&MGRID
|
||||
CUTOFF 300
|
||||
NGRIDS 5
|
||||
REL_CUTOFF 40
|
||||
&END MGRID
|
||||
&POISSON
|
||||
PERIODIC NONE
|
||||
PSOLVER MT
|
||||
&END POISSON
|
||||
&QS
|
||||
EPS_DEFAULT 1.0E-10
|
||||
METHOD GPW
|
||||
&END QS
|
||||
&SCF
|
||||
EPS_SCF 1.E-5
|
||||
MAX_SCF 25
|
||||
SCF_GUESS ATOMIC
|
||||
&END SCF
|
||||
&XC
|
||||
&XC_FUNCTIONAL LDA
|
||||
&END XC_FUNCTIONAL
|
||||
&END XC
|
||||
&END DFT
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC 6.0 6.0 6.0
|
||||
PERIODIC NONE
|
||||
&END CELL
|
||||
&COORD
|
||||
Cl 0.00000 0.00000 0.00000
|
||||
H 0.00000 0.00000 1.30000
|
||||
&END COORD
|
||||
&KIND H
|
||||
BASIS_SET DZV-GTH-PADE
|
||||
POTENTIAL ECP ccECP
|
||||
&END KIND
|
||||
&KIND Cl
|
||||
BASIS_SET DZVP-GTH-PADE
|
||||
POTENTIAL ECP ccECP
|
||||
&END KIND
|
||||
&TOPOLOGY
|
||||
&CENTER_COORDINATES
|
||||
&END CENTER_COORDINATES
|
||||
&END TOPOLOGY
|
||||
&END SUBSYS
|
||||
&END FORCE_EVAL
|
||||
65
tests/QS/regtest-ecp/HCl_force.inp
Normal file
65
tests/QS/regtest-ecp/HCl_force.inp
Normal file
|
|
@ -0,0 +1,65 @@
|
|||
&GLOBAL
|
||||
PRINT_LEVEL LOW
|
||||
PROJECT HCl
|
||||
RUN_TYPE DEBUG
|
||||
&END GLOBAL
|
||||
|
||||
&DEBUG
|
||||
CHECK_ATOM_FORCE 1 z
|
||||
DE 0.005
|
||||
DEBUG_DIPOLE .FALSE.
|
||||
DEBUG_FORCES .TRUE.
|
||||
DEBUG_POLARIZABILITY .FALSE.
|
||||
DEBUG_STRESS_TENSOR .FALSE.
|
||||
MAX_RELATIVE_ERROR 0.03
|
||||
&END DEBUG
|
||||
|
||||
&FORCE_EVAL
|
||||
&DFT
|
||||
POTENTIAL_FILE_NAME ./ECP_BASIS_POT
|
||||
&MGRID
|
||||
CUTOFF 300
|
||||
NGRIDS 5
|
||||
REL_CUTOFF 40
|
||||
&END MGRID
|
||||
&POISSON
|
||||
PERIODIC NONE
|
||||
PSOLVER MT
|
||||
&END POISSON
|
||||
&QS
|
||||
EPS_DEFAULT 1.0E-14
|
||||
METHOD GPW
|
||||
&END QS
|
||||
&SCF
|
||||
EPS_SCF 1.E-7
|
||||
MAX_SCF 50
|
||||
SCF_GUESS ATOMIC
|
||||
&END SCF
|
||||
&XC
|
||||
&XC_FUNCTIONAL LDA
|
||||
&END XC_FUNCTIONAL
|
||||
&END XC
|
||||
&END DFT
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC 6.0 6.0 6.0
|
||||
PERIODIC NONE
|
||||
&END CELL
|
||||
&COORD
|
||||
Cl 0.00000 0.00000 0.00000
|
||||
H 0.00000 0.00000 1.30000
|
||||
&END COORD
|
||||
&KIND H
|
||||
BASIS_SET DZV-GTH-PADE
|
||||
POTENTIAL ECP ccECP
|
||||
&END KIND
|
||||
&KIND Cl
|
||||
BASIS_SET DZVP-GTH-PADE
|
||||
POTENTIAL ECP ccECP
|
||||
&END KIND
|
||||
&TOPOLOGY
|
||||
&CENTER_COORDINATES
|
||||
&END CENTER_COORDINATES
|
||||
&END TOPOLOGY
|
||||
&END SUBSYS
|
||||
&END FORCE_EVAL
|
||||
|
|
@ -8,21 +8,20 @@
|
|||
BASIS_SET_FILE_NAME ./ECP_BASIS_POT
|
||||
POTENTIAL_FILE_NAME ./ECP_BASIS_POT
|
||||
&MGRID
|
||||
CUTOFF 400
|
||||
CUTOFF 300
|
||||
NGRIDS 5
|
||||
REL_CUTOFF 40
|
||||
&END MGRID
|
||||
&POISSON
|
||||
PERIODIC NONE
|
||||
PSOLVER WAVELET
|
||||
PSOLVER MT
|
||||
&END POISSON
|
||||
&QS
|
||||
EPS_DEFAULT 1.0E-12
|
||||
METHOD GAPW
|
||||
&END QS
|
||||
&SCF
|
||||
IGNORE_CONVERGENCE_FAILURE
|
||||
MAX_SCF 5
|
||||
MAX_SCF 50
|
||||
SCF_GUESS ATOMIC
|
||||
&END SCF
|
||||
&XC
|
||||
|
|
@ -32,7 +31,7 @@
|
|||
&END DFT
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC 10.0 10.0 10.0
|
||||
ABC 7.0 7.0 7.0
|
||||
PERIODIC NONE
|
||||
&END CELL
|
||||
&COORD
|
||||
|
|
|
|||
|
|
@ -1,3 +1,5 @@
|
|||
ICl_lanl2dz_gpw.inp 11 1.0E-11 -23.515911094537550
|
||||
Rn_stuttgart_gapw.inp 11 1.0E-11 -25.632304810373267
|
||||
#SbH3_def2_gapw.inp 11 1.0E-11 -212.357670903187426
|
||||
#SbH3_def2_gapw.inp 11 1.0E-11 -241.414769303399083
|
||||
HCl_ccECP.inp 11 1.0E-12 -15.457709431781252
|
||||
HCl_force.inp 0
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue