mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-28 22:25:32 -04:00
ECP nuclear gradients
This commit is contained in:
abussy
Augustin Bussy
parent
644f5a76d6
commit
48d3695c3c
8 changed files with 1053 additions and 41 deletions
|
|
@ -31,8 +31,10 @@ MODULE core_ppl
|
|||
sgp_potential_type
|
||||
USE kinds, ONLY: dp,&
|
||||
int_8
|
||||
USE libgrpp_integrals, ONLY: libgrpp_local_integral,&
|
||||
libgrpp_semilocal_integral
|
||||
USE libgrpp_integrals, ONLY: libgrpp_local_forces_ref,&
|
||||
libgrpp_local_integrals,&
|
||||
libgrpp_semilocal_forces_ref,&
|
||||
libgrpp_semilocal_integrals
|
||||
USE lri_environment_types, ONLY: lri_kind_type
|
||||
USE orbital_pointers, ONLY: init_orbital_pointers,&
|
||||
ncoset
|
||||
|
|
@ -475,12 +477,31 @@ CONTAINS
|
|||
hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
|
||||
force_a, force_b, ppl_fwork)
|
||||
ELSE
|
||||
CPABORT("ECP gradients NYI")
|
||||
|
||||
!$OMP CRITICAL(type1)
|
||||
CALL libgrpp_local_forces_ref(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, cval_ppl(1, :), nct_ppl, &
|
||||
ppl_radius, rab, dab, rac, dac, dbc, &
|
||||
hab(:, :, iset, jset), pab(:, :, iset, jset), &
|
||||
force_a, force_b)
|
||||
!$OMP END CRITICAL(type1)
|
||||
END IF
|
||||
|
||||
IF (ecp_semi_local) THEN
|
||||
! semi local ECP part -- forces
|
||||
CPABORT("ECP gradients NYI")
|
||||
|
||||
!$OMP CRITICAL(type2)
|
||||
CALL libgrpp_semilocal_forces_ref(la_max(iset), la_min(iset), npgfa(iset), &
|
||||
rpgfa(:, iset), zeta(:, iset), &
|
||||
lb_max(jset), lb_min(jset), npgfb(jset), &
|
||||
rpgfb(:, jset), zetb(:, jset), &
|
||||
slmax, npot, bpot, apot, nrpot, &
|
||||
ppl_radius, rab, dab, rac, dac, dbc, &
|
||||
hab(:, :, iset, jset), pab(:, :, iset, jset), &
|
||||
force_a, force_b)
|
||||
!$OMP END CRITICAL(type2)
|
||||
END IF
|
||||
! *** The derivatives w.r.t. atomic center c are ***
|
||||
! *** calculated using the translational invariance ***
|
||||
|
|
@ -535,26 +556,26 @@ CONTAINS
|
|||
ELSE
|
||||
!If the local part of the potential is more complex, we need libgrpp
|
||||
!$OMP CRITICAL(type1)
|
||||
CALL libgrpp_local_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, cval_ppl(1, :), nct_ppl, &
|
||||
ppl_radius, rab, dab, rac, dac, dbc, &
|
||||
hab(:, :, iset, jset))
|
||||
CALL libgrpp_local_integrals(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, cval_ppl(1, :), nct_ppl, &
|
||||
ppl_radius, rab, dab, rac, dac, dbc, &
|
||||
hab(:, :, iset, jset))
|
||||
!$OMP END CRITICAL(type1)
|
||||
END IF
|
||||
|
||||
IF (ecp_semi_local) THEN
|
||||
! semi local ECP part
|
||||
!$OMP CRITICAL(type2)
|
||||
CALL libgrpp_semilocal_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), &
|
||||
slmax, npot, bpot, apot, nrpot, &
|
||||
ppl_radius, rab, dab, rac, dac, dbc, &
|
||||
hab(:, :, iset, jset))
|
||||
CALL libgrpp_semilocal_integrals(la_max(iset), la_min(iset), npgfa(iset), &
|
||||
rpgfa(:, iset), zeta(:, iset), &
|
||||
lb_max(jset), lb_min(jset), npgfb(jset), &
|
||||
rpgfb(:, jset), zetb(:, jset), &
|
||||
slmax, npot, bpot, apot, nrpot, &
|
||||
ppl_radius, rab, dab, rac, dac, dbc, &
|
||||
hab(:, :, iset, jset))
|
||||
!$OMP END CRITICAL(type2)
|
||||
END IF
|
||||
END IF
|
||||
|
|
|
|||
|
|
@ -12,10 +12,12 @@
|
|||
MODULE libgrpp_integrals
|
||||
USE kinds, ONLY: dp
|
||||
USE mathconstants, ONLY: pi
|
||||
USE ai_derivatives, ONLY: dabdr_noscreen, adbdr, dabdr
|
||||
USE orbital_pointers, ONLY: nco, &
|
||||
ncoset
|
||||
#if defined(__LIBGRPP)
|
||||
USE libgrpp, ONLY: libgrpp_type1_integrals, libgrpp_type2_integrals
|
||||
USE libgrpp, ONLY: libgrpp_type1_integrals, libgrpp_type2_integrals, &
|
||||
libgrpp_type1_integrals_gradient, libgrpp_type2_integrals_gradient
|
||||
#endif
|
||||
#include "./base/base_uses.f90"
|
||||
|
||||
|
|
@ -24,7 +26,8 @@ MODULE libgrpp_integrals
|
|||
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'libgrpp_integrals'
|
||||
|
||||
PUBLIC :: libgrpp_semilocal_integral, libgrpp_local_integral
|
||||
PUBLIC :: libgrpp_semilocal_integrals, libgrpp_local_integrals, &
|
||||
libgrpp_local_forces_ref, libgrpp_semilocal_forces_ref
|
||||
|
||||
CONTAINS
|
||||
|
||||
|
|
@ -51,11 +54,14 @@ CONTAINS
|
|||
!> \param dac ...
|
||||
!> \param dbc ...
|
||||
!> \param vab ...
|
||||
!> \param pab ...
|
||||
!> \param force_a ...
|
||||
!> \param force_b ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE libgrpp_local_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
|
||||
lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
|
||||
npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
|
||||
rpgfc, rab, dab, rac, dac, dbc, vab)
|
||||
SUBROUTINE libgrpp_local_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
|
||||
lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
|
||||
npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
|
||||
rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
|
||||
|
||||
INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
|
||||
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
|
||||
|
|
@ -71,20 +77,76 @@ CONTAINS
|
|||
REAL(KIND=dp), INTENT(IN) :: dac
|
||||
REAL(KIND=dp), INTENT(IN) :: dbc
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
|
||||
OPTIONAL :: pab
|
||||
REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
|
||||
OPTIONAL :: force_a, force_b
|
||||
|
||||
#if defined(__LIBGRPP)
|
||||
INTEGER :: a_offset, a_start, b_offset, b_start, i, &
|
||||
ipgf, j, jpgf, li, lj, ncoa, ncob
|
||||
LOGICAL :: calc_forces
|
||||
REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
|
||||
zeti, zetj
|
||||
REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
|
||||
zeti, zetj, mindist, fac_a, fac_b
|
||||
REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmpx, tmpy, tmpz
|
||||
REAL(dp), DIMENSION(3) :: ra, rb, rc
|
||||
|
||||
calc_forces = .FALSE.
|
||||
IF (PRESENT(pab) .AND. PRESENT(force_a) .AND. PRESENT(force_b)) calc_forces = .TRUE.
|
||||
|
||||
IF (calc_forces) THEN
|
||||
|
||||
!Note: warning against numerical stability of libgrpp gradients. The day the library becomes
|
||||
! stable, this routine can be used immediatly as is, and the warning removed.
|
||||
CALL cp_warn(__LOCATION__, &
|
||||
"ECP gradients calculated with the libgrpp library are, to this date, not numerically stable. "// &
|
||||
"Please use the reference routine 'libgrpp_local_forces_ref' instead.")
|
||||
|
||||
!there is a weird feature of libgrpp gradients, which is such that the gradient is calculated
|
||||
!for a point in space, and not with respect to an atomic center. For example, if atoms A and
|
||||
!B are the same (and C is different), then d<A | U_C | B>/dPx = d<A | U_C | B>/dAx + d<A | U_C | B>/dBx
|
||||
!Because we want the forces on centers A and B seprately, we need a case study on atomic positions
|
||||
!We always calculate the gradient wrt to atomic position of A and B, and we scale accordingly
|
||||
|
||||
mindist = 1.0E-6_dp
|
||||
!If ra != rb != rc
|
||||
IF (dab >= mindist .AND. dbc >= mindist .AND. dac >= mindist) THEN
|
||||
fac_a = 1.0_dp
|
||||
fac_b = 1.0_dp
|
||||
|
||||
!If ra = rb, but ra != rc
|
||||
ELSE IF (dab < mindist .AND. dac >= mindist) THEN
|
||||
fac_a = 0.5_dp
|
||||
fac_b = 0.5_dp
|
||||
|
||||
!IF ra != rb but ra = rc
|
||||
ELSE IF (dab >= mindist .AND. dac < mindist) THEN
|
||||
fac_a = 0.5_dp
|
||||
fac_b = 1.0_dp
|
||||
|
||||
!IF ra != rb but rb = rc
|
||||
ELSE IF (dab >= mindist .AND. dbc < mindist) THEN
|
||||
fac_a = 1.0_dp
|
||||
fac_b = 0.5_dp
|
||||
|
||||
!If all atoms the same --> no force
|
||||
ELSE
|
||||
calc_forces = .FALSE.
|
||||
END IF
|
||||
END IF
|
||||
|
||||
!libgrpp requires absolute positions, not relative ones
|
||||
ra(:) = 0.0_dp
|
||||
rb(:) = rab(:)
|
||||
rc(:) = rac(:)
|
||||
|
||||
ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
|
||||
IF (calc_forces) THEN
|
||||
ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
|
||||
ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
|
||||
ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
|
||||
END IF
|
||||
|
||||
DO ipgf = 1, npgfa
|
||||
IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
|
||||
zeti = zeta(ipgf)
|
||||
|
|
@ -110,8 +172,7 @@ CONTAINS
|
|||
expj = 0.25_dp*REAL(2*lj + 3, dp)
|
||||
normj = 1.0_dp/(prefj*zetj**expj)
|
||||
|
||||
ALLOCATE (tmp(ncoa*ncob))
|
||||
tmp = 0.0_dp
|
||||
tmp(1:ncoa*ncob) = 0.0_dp
|
||||
!libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
|
||||
!the 1/norm coefficients for PGFi and PGFj
|
||||
CALL libgrpp_type1_integrals(ra, li, 1, [normi], [zeti], &
|
||||
|
|
@ -125,7 +186,50 @@ CONTAINS
|
|||
vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
|
||||
END DO
|
||||
END DO
|
||||
DEALLOCATE (tmp)
|
||||
|
||||
IF (calc_forces) THEN
|
||||
tmpx(1:ncoa*ncob) = 0.0_dp
|
||||
tmpy(1:ncoa*ncob) = 0.0_dp
|
||||
tmpz(1:ncoa*ncob) = 0.0_dp
|
||||
|
||||
!force wrt to atomic position A
|
||||
CALL libgrpp_type1_integrals_gradient(ra, li, 1, [normi], [zeti], &
|
||||
rb, lj, 1, [normj], [zetj], &
|
||||
rc, [npot_ecp], nrpot_ecp, &
|
||||
coeffs_ecp, alpha_ecp, ra, &
|
||||
tmpx, tmpy, tmpz)
|
||||
|
||||
!note: tmp array is in C row-major ordering
|
||||
!note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
|
||||
DO j = 1, ncob
|
||||
DO i = 1, ncoa
|
||||
force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
|
||||
force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
|
||||
force_a(3) = force_a(3) + fac_a*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
tmpx(1:ncoa*ncob) = 0.0_dp
|
||||
tmpy(1:ncoa*ncob) = 0.0_dp
|
||||
tmpz(1:ncoa*ncob) = 0.0_dp
|
||||
|
||||
!force wrt to atomic position B
|
||||
CALL libgrpp_type1_integrals_gradient(ra, li, 1, [normi], [zeti], &
|
||||
rb, lj, 1, [normj], [zetj], &
|
||||
rc, [npot_ecp], nrpot_ecp, &
|
||||
coeffs_ecp, alpha_ecp, rb, &
|
||||
tmpx, tmpy, tmpz)
|
||||
|
||||
!note: tmp array is in C row-major ordering
|
||||
!note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
|
||||
DO j = 1, ncob
|
||||
DO i = 1, ncoa
|
||||
force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
|
||||
force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
|
||||
force_b(3) = force_b(3) + fac_b*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
|
||||
END DO
|
||||
END DO
|
||||
END IF
|
||||
|
||||
END DO !lj
|
||||
END DO !li
|
||||
|
|
@ -155,11 +259,14 @@ CONTAINS
|
|||
MARK_USED(dac)
|
||||
MARK_USED(dbc)
|
||||
MARK_USED(vab)
|
||||
MARK_USED(pab)
|
||||
MARK_USED(force_a)
|
||||
MARK_USED(force_b)
|
||||
|
||||
CPABORT("Please compile CP2K with libgrpp support for calculations with ECPs")
|
||||
#endif
|
||||
|
||||
END SUBROUTINE libgrpp_local_integral
|
||||
END SUBROUTINE libgrpp_local_integrals
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Semi-local ECP integrals using libgrpp.
|
||||
|
|
@ -185,11 +292,14 @@ CONTAINS
|
|||
!> \param dac ...
|
||||
!> \param dbc ...
|
||||
!> \param vab ...
|
||||
!> \param pab ...
|
||||
!> \param force_a ...
|
||||
!> \param force_b ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE libgrpp_semilocal_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
|
||||
lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
|
||||
lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
|
||||
rpgfc, rab, dab, rac, dac, dbc, vab)
|
||||
SUBROUTINE libgrpp_semilocal_integrals(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
|
||||
lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
|
||||
lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
|
||||
rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
|
||||
|
||||
INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
|
||||
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
|
||||
|
|
@ -206,20 +316,76 @@ CONTAINS
|
|||
REAL(KIND=dp), INTENT(IN) :: dac
|
||||
REAL(KIND=dp), INTENT(IN) :: dbc
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
|
||||
OPTIONAL :: pab
|
||||
REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
|
||||
OPTIONAL :: force_a, force_b
|
||||
|
||||
#if defined(__LIBGRPP)
|
||||
INTEGER :: a_offset, a_start, b_offset, b_start, i, &
|
||||
ipgf, j, jpgf, li, lj, lk, ncoa, ncob
|
||||
LOGICAL :: calc_forces
|
||||
REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
|
||||
zeti, zetj
|
||||
REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
|
||||
zeti, zetj, mindist, fac_a, fac_b
|
||||
REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmpx, tmpz, tmpy
|
||||
REAL(dp), DIMENSION(3) :: ra, rb, rc
|
||||
|
||||
calc_forces = .FALSE.
|
||||
IF (PRESENT(pab) .AND. PRESENT(force_a) .AND. PRESENT(force_b)) calc_forces = .TRUE.
|
||||
|
||||
IF (calc_forces) THEN
|
||||
|
||||
!Note: warning against numerical stability of libgrpp gradients. The day the library becomes
|
||||
! stable, this routine can be used immediatly as is, and the warning removed.
|
||||
CALL cp_warn(__LOCATION__, &
|
||||
"ECP gradients calculated with the libgrpp library are, to this date, not numerically stable. "// &
|
||||
"Please use the reference routine 'libgrpp_semilocal_forces_ref' instead.")
|
||||
|
||||
!there is a weird feature of libgrpp gradients, which is such that the gradient is calculated
|
||||
!for a point in space, and not with respect to an atomic center. For example, if atoms A and
|
||||
!B are the same (and C is different), then d<A | U_C | B>/dPx = d<A | U_C | B>/dAx + d<A | U_C | B>/dBx
|
||||
!Because we want the forces on centers A and B seprately, we need a case study on atomic positions
|
||||
!We always calculate the gradient wrt to atomic position of A and B, and we scale accordingly
|
||||
|
||||
mindist = 1.0E-6_dp
|
||||
!If ra != rb != rc
|
||||
IF (dab >= mindist .AND. dbc >= mindist .AND. dac >= mindist) THEN
|
||||
fac_a = 1.0_dp
|
||||
fac_b = 1.0_dp
|
||||
|
||||
!If ra = rb, but ra != rc
|
||||
ELSE IF (dab < mindist .AND. dac >= mindist) THEN
|
||||
fac_a = 0.5_dp
|
||||
fac_b = 0.5_dp
|
||||
|
||||
!IF ra != rb but ra = rc
|
||||
ELSE IF (dab >= mindist .AND. dac < mindist) THEN
|
||||
fac_a = 0.5_dp
|
||||
fac_b = 1.0_dp
|
||||
|
||||
!IF ra != rb but rb = rc
|
||||
ELSE IF (dab >= mindist .AND. dbc < mindist) THEN
|
||||
fac_a = 1.0_dp
|
||||
fac_b = 0.5_dp
|
||||
|
||||
!If all atoms the same --> no force
|
||||
ELSE
|
||||
calc_forces = .FALSE.
|
||||
END IF
|
||||
END IF
|
||||
|
||||
!libgrpp requires absolute positions, not relative ones
|
||||
ra(:) = 0.0_dp
|
||||
rb(:) = rab(:)
|
||||
rc(:) = rac(:)
|
||||
|
||||
ALLOCATE (tmp(nco(la_max_set)*nco(lb_max_set)))
|
||||
IF (calc_forces) THEN
|
||||
ALLOCATE (tmpx(nco(la_max_set)*nco(lb_max_set)))
|
||||
ALLOCATE (tmpy(nco(la_max_set)*nco(lb_max_set)))
|
||||
ALLOCATE (tmpz(nco(la_max_set)*nco(lb_max_set)))
|
||||
END IF
|
||||
|
||||
DO ipgf = 1, npgfa
|
||||
IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
|
||||
zeti = zeta(ipgf)
|
||||
|
|
@ -245,10 +411,9 @@ CONTAINS
|
|||
expj = 0.25_dp*REAL(2*lj + 3, dp)
|
||||
normj = 1.0_dp/(prefj*zetj**expj)
|
||||
|
||||
ALLOCATE (tmp(ncoa*ncob))
|
||||
!Loop over ECP angular momentum
|
||||
DO lk = 0, lmax_ecp
|
||||
tmp = 0.0_dp
|
||||
tmp(1:ncoa*ncob) = 0.0_dp
|
||||
!libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
|
||||
!the 1/norm coefficients for PGFi and PGFj
|
||||
CALL libgrpp_type2_integrals(ra, li, 1, [normi], [zeti], &
|
||||
|
|
@ -262,8 +427,53 @@ CONTAINS
|
|||
vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
IF (calc_forces) THEN
|
||||
|
||||
tmpx(1:ncoa*ncob) = 0.0_dp
|
||||
tmpy(1:ncoa*ncob) = 0.0_dp
|
||||
tmpz(1:ncoa*ncob) = 0.0_dp
|
||||
|
||||
!force wrt to atomic position A
|
||||
CALL libgrpp_type2_integrals_gradient(ra, li, 1, [normi], [zeti], &
|
||||
rb, lj, 1, [normj], [zetj], &
|
||||
rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
|
||||
coeffs_ecp(:, lk), alpha_ecp(:, lk), ra, &
|
||||
tmpx, tmpy, tmpz)
|
||||
|
||||
!note: tmp array is in C row-major ordering
|
||||
!note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
|
||||
DO j = 1, ncob
|
||||
DO i = 1, ncoa
|
||||
force_a(1) = force_a(1) + fac_a*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
|
||||
force_a(2) = force_a(2) + fac_a*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
|
||||
force_a(3) = force_a(3) + fac_a*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
tmpx(1:ncoa*ncob) = 0.0_dp
|
||||
tmpy(1:ncoa*ncob) = 0.0_dp
|
||||
tmpz(1:ncoa*ncob) = 0.0_dp
|
||||
|
||||
!force wrt to atomic position B
|
||||
CALL libgrpp_type2_integrals_gradient(ra, li, 1, [normi], [zeti], &
|
||||
rb, lj, 1, [normj], [zetj], &
|
||||
rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
|
||||
coeffs_ecp(:, lk), alpha_ecp(:, lk), rb, &
|
||||
tmpx, tmpy, tmpz)
|
||||
!note: tmp array is in C row-major ordering
|
||||
!note: zero-gradients sometime comes out as NaN, hence tampval==tmpval check
|
||||
DO j = 1, ncob
|
||||
DO i = 1, ncoa
|
||||
force_b(1) = force_b(1) + fac_b*pab(a_offset + i, b_offset + j)*tmpx((i - 1)*ncob + j)
|
||||
force_b(2) = force_b(2) + fac_b*pab(a_offset + i, b_offset + j)*tmpy((i - 1)*ncob + j)
|
||||
force_b(3) = force_b(3) + fac_b*pab(a_offset + i, b_offset + j)*tmpz((i - 1)*ncob + j)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
END IF !calc_forces
|
||||
|
||||
END DO !lk
|
||||
DEALLOCATE (tmp)
|
||||
|
||||
END DO !lj
|
||||
END DO !li
|
||||
|
|
@ -295,10 +505,416 @@ CONTAINS
|
|||
MARK_USED(dac)
|
||||
MARK_USED(dbc)
|
||||
MARK_USED(vab)
|
||||
MARK_USED(pab)
|
||||
MARK_USED(force_a)
|
||||
MARK_USED(force_b)
|
||||
|
||||
CPABORT("Please compile CP2K with libgrpp support for calculations with ECPs")
|
||||
#endif
|
||||
|
||||
END SUBROUTINE libgrpp_semilocal_integral
|
||||
END SUBROUTINE libgrpp_semilocal_integrals
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Reference local ECP force routine using l+-1 integrals. No call is made to the numerically
|
||||
!> unstable gradient routine of libgrpp. Calculates both the integrals and the forces.
|
||||
!> \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 npot_ecp ...
|
||||
!> \param alpha_ecp ...
|
||||
!> \param coeffs_ecp ...
|
||||
!> \param nrpot_ecp ...
|
||||
!> \param rpgfc ...
|
||||
!> \param rab ...
|
||||
!> \param dab ...
|
||||
!> \param rac ...
|
||||
!> \param dac ...
|
||||
!> \param dbc ...
|
||||
!> \param vab ...
|
||||
!> \param pab ...
|
||||
!> \param force_a ...
|
||||
!> \param force_b ...
|
||||
!> \note: this is a reference routine, which has no reason to be used once the libgrpp gradients
|
||||
!> become numerically stable
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE libgrpp_local_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
|
||||
lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
|
||||
npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
|
||||
rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
|
||||
|
||||
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) :: npot_ecp
|
||||
REAL(KIND=dp), DIMENSION(1:npot_ecp), INTENT(IN) :: alpha_ecp, coeffs_ecp
|
||||
INTEGER, DIMENSION(1:npot_ecp), INTENT(IN) :: nrpot_ecp
|
||||
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), INTENT(IN) :: dbc
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pab
|
||||
REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: force_a, force_b
|
||||
|
||||
#if defined(__LIBGRPP)
|
||||
INTEGER :: a_offset, a_start, b_offset, b_start, i, &
|
||||
ipgf, j, jpgf, li, lj, ncoa, ncob, a_offset_f, &
|
||||
b_offset_f, a_start_f, b_start_f
|
||||
REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
|
||||
zeti, zetj
|
||||
REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
|
||||
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vab_f, tmpx, tmpy, tmpz
|
||||
REAL(dp), DIMENSION(3) :: ra, rb, rc
|
||||
|
||||
!Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
|
||||
ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
|
||||
vab_f(:, :) = 0.0_dp
|
||||
|
||||
!libgrpp requires absolute positions, not relative ones
|
||||
ra(:) = 0.0_dp
|
||||
rb(:) = rab(:)
|
||||
rc(:) = rac(:)
|
||||
|
||||
ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
|
||||
|
||||
DO ipgf = 1, npgfa
|
||||
IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
|
||||
zeti = zeta(ipgf)
|
||||
a_start = (ipgf - 1)*ncoset(la_max_set)
|
||||
a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
|
||||
|
||||
DO jpgf = 1, npgfb
|
||||
IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
|
||||
IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
|
||||
zetj = zetb(jpgf)
|
||||
b_start = (jpgf - 1)*ncoset(lb_max_set)
|
||||
b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
|
||||
|
||||
DO li = MAX(0, la_min_set - 1), la_max_set + 1
|
||||
a_offset = a_start + ncoset(li - 1)
|
||||
a_offset_f = a_start_f + ncoset(li - 1)
|
||||
ncoa = nco(li)
|
||||
prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
|
||||
expi = 0.25_dp*REAL(2*li + 3, dp)
|
||||
normi = 1.0_dp/(prefi*zeti**expi)
|
||||
|
||||
DO lj = MAX(0, lb_min_set - 1), lb_max_set + 1
|
||||
b_offset = b_start + ncoset(lj - 1)
|
||||
b_offset_f = b_start_f + ncoset(lj - 1)
|
||||
ncob = nco(lj)
|
||||
prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
|
||||
expj = 0.25_dp*REAL(2*lj + 3, dp)
|
||||
normj = 1.0_dp/(prefj*zetj**expj)
|
||||
|
||||
tmp(1:ncoa*ncob) = 0.0_dp
|
||||
!libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
|
||||
!the 1/norm coefficients for PGFi and PGFj
|
||||
CALL libgrpp_type1_integrals(ra, li, 1, [normi], [zeti], &
|
||||
rb, lj, 1, [normj], [zetj], &
|
||||
rc, [npot_ecp], nrpot_ecp, &
|
||||
coeffs_ecp, alpha_ecp, tmp)
|
||||
|
||||
!the l+-1 integrals for gradient calculation
|
||||
DO j = 1, ncob
|
||||
DO i = 1, ncoa
|
||||
vab_f(a_offset_f + i, b_offset_f + j) = &
|
||||
vab_f(a_offset_f + i, b_offset_f + j) + tmp((i - 1)*ncob + j)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
!the actual integrals
|
||||
IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
|
||||
DO j = 1, ncob
|
||||
DO i = 1, ncoa
|
||||
vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
|
||||
END DO
|
||||
END DO
|
||||
END IF
|
||||
|
||||
END DO !lj
|
||||
END DO !li
|
||||
|
||||
END DO !jpgf
|
||||
END DO !ipgf
|
||||
|
||||
ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
|
||||
ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
|
||||
ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
|
||||
|
||||
!Derivative wrt to center A
|
||||
tmpx(:, :) = 0.0_dp
|
||||
tmpy(:, :) = 0.0_dp
|
||||
tmpz(:, :) = 0.0_dp
|
||||
CALL dabdr(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, rpgfb, lb_min_set, &
|
||||
dab, vab_f, tmpx, tmpy, tmpz)
|
||||
DO j = 1, npgfb*ncoset(lb_max_set)
|
||||
DO i = 1, npgfa*ncoset(la_max_set)
|
||||
force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
|
||||
force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
|
||||
force_a(3) = force_a(3) + tmpz(i, j)*pab(i, j)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
!Derivative wrt to center B
|
||||
tmpx(:, :) = 0.0_dp
|
||||
tmpy(:, :) = 0.0_dp
|
||||
tmpz(:, :) = 0.0_dp
|
||||
CALL adbdr(la_max_set, npgfa, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
|
||||
dab, vab_f, tmpx, tmpy, tmpz)
|
||||
DO j = 1, npgfb*ncoset(lb_max_set)
|
||||
DO i = 1, npgfa*ncoset(la_max_set)
|
||||
force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
|
||||
force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
|
||||
force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
|
||||
END DO
|
||||
END DO
|
||||
DEALLOCATE (tmpx, tmpy, tmpz)
|
||||
|
||||
#else
|
||||
|
||||
MARK_USED(la_max_set)
|
||||
MARK_USED(la_min_set)
|
||||
MARK_USED(npgfa)
|
||||
MARK_USED(rpgfa)
|
||||
MARK_USED(zeta)
|
||||
MARK_USED(lb_max_set)
|
||||
MARK_USED(lb_min_set)
|
||||
MARK_USED(npgfb)
|
||||
MARK_USED(rpgfb)
|
||||
MARK_USED(zetb)
|
||||
MARK_USED(npot_ecp)
|
||||
MARK_USED(alpha_ecp)
|
||||
MARK_USED(coeffs_ecp)
|
||||
MARK_USED(nrpot_ecp)
|
||||
MARK_USED(rpgfc)
|
||||
MARK_USED(rab)
|
||||
MARK_USED(dab)
|
||||
MARK_USED(rac)
|
||||
MARK_USED(dac)
|
||||
MARK_USED(dbc)
|
||||
MARK_USED(pab)
|
||||
MARK_USED(vab)
|
||||
MARK_USED(force_a)
|
||||
MARK_USED(force_b)
|
||||
|
||||
CPABORT("Please compile CP2K with libgrpp support for calculations with ECPs")
|
||||
#endif
|
||||
|
||||
END SUBROUTINE libgrpp_local_forces_ref
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Reference semi-local ECP forces using l+-1 integrals. No call is made to the numerically
|
||||
!> unstable gradient routine of libgrpp. Calculates both the integrals and the forces.
|
||||
!> \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 lmax_ecp ...
|
||||
!> \param npot_ecp ...
|
||||
!> \param alpha_ecp ...
|
||||
!> \param coeffs_ecp ...
|
||||
!> \param nrpot_ecp ...
|
||||
!> \param rpgfc ...
|
||||
!> \param rab ...
|
||||
!> \param dab ...
|
||||
!> \param rac ...
|
||||
!> \param dac ...
|
||||
!> \param dbc ...
|
||||
!> \param vab ...
|
||||
!> \param pab ...
|
||||
!> \param force_a ...
|
||||
!> \param force_b ...
|
||||
!> \note: this is a reference routine, which has no reason to be used once the libgrpp gradients
|
||||
!> become numerically stable
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE libgrpp_semilocal_forces_ref(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
|
||||
lb_max_set, lb_min_set, npgfb, rpgfb, zetb, &
|
||||
lmax_ecp, npot_ecp, alpha_ecp, coeffs_ecp, nrpot_ecp, &
|
||||
rpgfc, rab, dab, rac, dac, dbc, vab, pab, force_a, force_b)
|
||||
|
||||
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) :: lmax_ecp
|
||||
INTEGER, DIMENSION(0:10), INTENT(IN) :: npot_ecp
|
||||
REAL(KIND=dp), DIMENSION(1:15, 0:10), INTENT(IN) :: alpha_ecp, coeffs_ecp
|
||||
INTEGER, DIMENSION(1:15, 0:10), INTENT(IN) :: nrpot_ecp
|
||||
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), INTENT(IN) :: dbc
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pab
|
||||
REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: force_a, force_b
|
||||
|
||||
#if defined(__LIBGRPP)
|
||||
INTEGER :: a_offset, a_start, b_offset, b_start, i, &
|
||||
ipgf, j, jpgf, li, lj, lk, ncoa, ncob, &
|
||||
a_start_f, b_start_f, a_offset_f, b_offset_f
|
||||
REAL(dp) :: expi, expj, normi, normj, prefi, prefj, &
|
||||
zeti, zetj
|
||||
REAL(dp), ALLOCATABLE, DIMENSION(:) :: tmp
|
||||
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vab_f, tmpx, tmpy, tmpz
|
||||
REAL(dp), DIMENSION(3) :: ra, rb, rc
|
||||
|
||||
!Contains the integrals necessary for the forces, with angular momenta from lmin-1 to lmax+1
|
||||
ALLOCATE (vab_f(npgfa*ncoset(la_max_set + 1), npgfb*ncoset(lb_max_set + 1)))
|
||||
vab_f(:, :) = 0.0_dp
|
||||
|
||||
!libgrpp requires absolute positions, not relative ones
|
||||
ra(:) = 0.0_dp
|
||||
rb(:) = rab(:)
|
||||
rc(:) = rac(:)
|
||||
|
||||
ALLOCATE (tmp(nco(la_max_set + 1)*nco(lb_max_set + 1)))
|
||||
|
||||
DO ipgf = 1, npgfa
|
||||
IF (rpgfa(ipgf) + rpgfc < dac) CYCLE
|
||||
zeti = zeta(ipgf)
|
||||
a_start = (ipgf - 1)*ncoset(la_max_set)
|
||||
a_start_f = (ipgf - 1)*ncoset(la_max_set + 1)
|
||||
|
||||
DO jpgf = 1, npgfb
|
||||
IF (rpgfb(jpgf) + rpgfc < dbc) CYCLE
|
||||
IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) CYCLE
|
||||
zetj = zetb(jpgf)
|
||||
b_start = (jpgf - 1)*ncoset(lb_max_set)
|
||||
b_start_f = (jpgf - 1)*ncoset(lb_max_set + 1)
|
||||
|
||||
DO li = MAX(0, la_min_set - 1), la_max_set + 1
|
||||
a_offset = a_start + ncoset(li - 1)
|
||||
a_offset_f = a_start_f + ncoset(li - 1)
|
||||
ncoa = nco(li)
|
||||
prefi = 2.0_dp**li*(2.0_dp/pi)**0.75_dp
|
||||
expi = 0.25_dp*REAL(2*li + 3, dp)
|
||||
normi = 1.0_dp/(prefi*zeti**expi)
|
||||
|
||||
DO lj = MAX(0, lb_min_set - 1), lb_max_set + 1
|
||||
b_offset = b_start + ncoset(lj - 1)
|
||||
b_offset_f = b_start_f + ncoset(lj - 1)
|
||||
ncob = nco(lj)
|
||||
prefj = 2.0_dp**lj*(2.0_dp/pi)**0.75_dp
|
||||
expj = 0.25_dp*REAL(2*lj + 3, dp)
|
||||
normj = 1.0_dp/(prefj*zetj**expj)
|
||||
|
||||
!Loop over ECP angular momentum
|
||||
DO lk = 0, lmax_ecp
|
||||
tmp(1:ncoa*ncob) = 0.0_dp
|
||||
!libgrpp implicitely normalizes cartesian Gaussian. In CP2K, we do not, hence
|
||||
!the 1/norm coefficients for PGFi and PGFj
|
||||
CALL libgrpp_type2_integrals(ra, li, 1, [normi], [zeti], &
|
||||
rb, lj, 1, [normj], [zetj], &
|
||||
rc, lk, [npot_ecp(lk)], nrpot_ecp(:, lk), &
|
||||
coeffs_ecp(:, lk), alpha_ecp(:, lk), tmp)
|
||||
|
||||
!the l+-1 integrals for gradient calculation
|
||||
DO j = 1, ncob
|
||||
DO i = 1, ncoa
|
||||
vab_f(a_offset_f + i, b_offset_f + j) = &
|
||||
vab_f(a_offset_f + i, b_offset_f + j) + tmp((i - 1)*ncob + j)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
!the actual integrals
|
||||
IF (li >= la_min_set .AND. li <= la_max_set .AND. lj >= lb_min_set .AND. lj <= lb_max_set) THEN
|
||||
DO j = 1, ncob
|
||||
DO i = 1, ncoa
|
||||
vab(a_offset + i, b_offset + j) = vab(a_offset + i, b_offset + j) + tmp((i - 1)*ncob + j)
|
||||
END DO
|
||||
END DO
|
||||
END IF
|
||||
|
||||
END DO !lk
|
||||
|
||||
END DO !lj
|
||||
END DO !li
|
||||
|
||||
END DO !jpgf
|
||||
END DO !ipgf
|
||||
|
||||
ALLOCATE (tmpx(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
|
||||
ALLOCATE (tmpy(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
|
||||
ALLOCATE (tmpz(npgfa*ncoset(la_max_set), npgfb*ncoset(lb_max_set)))
|
||||
|
||||
!Derivative wrt to center A
|
||||
tmpx(:, :) = 0.0_dp
|
||||
tmpy(:, :) = 0.0_dp
|
||||
tmpz(:, :) = 0.0_dp
|
||||
CALL dabdr(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, rpgfb, lb_min_set, &
|
||||
0.0_dp, vab_f, tmpx, tmpy, tmpz)
|
||||
DO j = 1, npgfb*ncoset(lb_max_set)
|
||||
DO i = 1, npgfa*ncoset(la_max_set)
|
||||
force_a(1) = force_a(1) + tmpx(i, j)*pab(i, j)
|
||||
force_a(2) = force_a(2) + tmpy(i, j)*pab(i, j)
|
||||
force_a(3) = force_a(3) + tmpz(i, j)*pab(i, j)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
!Derivative wrt to center B
|
||||
tmpx(:, :) = 0.0_dp
|
||||
tmpy(:, :) = 0.0_dp
|
||||
tmpz(:, :) = 0.0_dp
|
||||
CALL adbdr(la_max_set, npgfa, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, &
|
||||
0.0_dp, vab_f, tmpx, tmpy, tmpz)
|
||||
DO j = 1, npgfb*ncoset(lb_max_set)
|
||||
DO i = 1, npgfa*ncoset(la_max_set)
|
||||
force_b(1) = force_b(1) + tmpx(i, j)*pab(i, j)
|
||||
force_b(2) = force_b(2) + tmpy(i, j)*pab(i, j)
|
||||
force_b(3) = force_b(3) + tmpz(i, j)*pab(i, j)
|
||||
END DO
|
||||
END DO
|
||||
DEALLOCATE (tmpx, tmpy, tmpz)
|
||||
|
||||
#else
|
||||
|
||||
MARK_USED(la_max_set)
|
||||
MARK_USED(la_min_set)
|
||||
MARK_USED(npgfa)
|
||||
MARK_USED(rpgfa)
|
||||
MARK_USED(zeta)
|
||||
MARK_USED(lb_max_set)
|
||||
MARK_USED(lb_min_set)
|
||||
MARK_USED(npgfb)
|
||||
MARK_USED(rpgfb)
|
||||
MARK_USED(zetb)
|
||||
MARK_USED(lmax_ecp)
|
||||
MARK_USED(npot_ecp)
|
||||
MARK_USED(alpha_ecp)
|
||||
MARK_USED(coeffs_ecp)
|
||||
MARK_USED(nrpot_ecp)
|
||||
MARK_USED(rpgfc)
|
||||
MARK_USED(rab)
|
||||
MARK_USED(dab)
|
||||
MARK_USED(rac)
|
||||
MARK_USED(dac)
|
||||
MARK_USED(dbc)
|
||||
MARK_USED(pab)
|
||||
MARK_USED(vab)
|
||||
MARK_USED(force_a)
|
||||
MARK_USED(force_b)
|
||||
|
||||
CPABORT("Please compile CP2K with libgrpp support for calculations with ECPs")
|
||||
#endif
|
||||
|
||||
END SUBROUTINE libgrpp_semilocal_forces_ref
|
||||
|
||||
END MODULE libgrpp_integrals
|
||||
|
|
|
|||
203
tests/QS/regtest-ecp-2/ECP_BASIS_POT
Normal file
203
tests/QS/regtest-ecp-2/ECP_BASIS_POT
Normal file
|
|
@ -0,0 +1,203 @@
|
|||
#----------------------------------------------------------------------
|
||||
# Basis Set Exchange
|
||||
# Version v0.9.1
|
||||
# https://www.basissetexchange.org
|
||||
#----------------------------------------------------------------------
|
||||
# Basis set: CRENBL
|
||||
# Description: CRENBL designed for use with small core potentials
|
||||
# Role: orbital
|
||||
# Version: 0 (Data from the Original Basis Set Exchange)
|
||||
#----------------------------------------------------------------------
|
||||
|
||||
# Radon Stuttgart RLC (4s,4p,1d) -> [2s,2p,1d]
|
||||
Rn Stuttgart-RLC
|
||||
5
|
||||
1 0 0 3 1
|
||||
1.9783250 0.6789880
|
||||
1.5140330 -1.1841590
|
||||
0.3246540 0.9210860
|
||||
1 0 0 1 1
|
||||
0.1273660 1.0000000
|
||||
1 1 1 3 1
|
||||
2.0310950 0.2320400
|
||||
1.6561220 -0.3545180
|
||||
0.2983650 0.6651200
|
||||
1 1 1 1 1
|
||||
0.1032220 1.0000000
|
||||
1 2 2 1 1
|
||||
0.2600000 1.0000000
|
||||
|
||||
## Effective core potentials
|
||||
Stuttgart_RLC_ECP
|
||||
Rn nelec 78
|
||||
Rn ul
|
||||
2 1.000000000 0.000000000
|
||||
Rn S
|
||||
2 0.922386000 -5.019005000
|
||||
2 1.781915000 37.036790000
|
||||
2 10.804601000 195.103308000
|
||||
Rn P
|
||||
2 0.724291000 -1.966481000
|
||||
2 1.363860000 23.464059000
|
||||
Rn D
|
||||
2 0.769400000 7.483457000
|
||||
2 1.538800000 9.361900000
|
||||
Rn F
|
||||
2 1.213897000 -6.763150000
|
||||
Rn G
|
||||
2 1.576469000 -9.915662000
|
||||
END Stuttgart_RLC_ECP
|
||||
|
||||
# Hydrogen def2-SVP (4s,1p) -> [2s,1p]
|
||||
H def2-SVP
|
||||
3
|
||||
1 0 0 3 1
|
||||
13.0107010 0.19682158E-01
|
||||
1.9622572 0.13796524
|
||||
0.44453796 0.47831935
|
||||
1 0 0 1 1
|
||||
0.12194962 1.0000000
|
||||
1 1 1 1 1
|
||||
0.8000000 1.0000000
|
||||
|
||||
# Antimony def2-SVP (10s,7p,6d) -> [4s,4p,2d]
|
||||
Sb def2-SVP
|
||||
10
|
||||
1 0 0 2 1
|
||||
10.584496987 -0.14845336778E-01
|
||||
1.4680242769 0.35289492025
|
||||
1 0 0 6 1
|
||||
372.76139166 0.15878057239E-02
|
||||
22.689478596 -0.15027605583
|
||||
18.391547037 0.35915813039
|
||||
7.6406271414 -0.74805091065
|
||||
1.9052000235 0.92017581656
|
||||
0.93007107773 0.46754079597
|
||||
1 0 0 1 1
|
||||
0.21621190032 1.0000000
|
||||
1 0 0 1 1
|
||||
0.83546307041E-01 1.0000000
|
||||
1 1 1 1 1
|
||||
2.6190751238 1.0000000
|
||||
1 1 1 4 1
|
||||
15.926550950 0.13206950012
|
||||
10.052739237 -0.41511149340
|
||||
1.2682183726 0.74197972626
|
||||
0.57196620929 0.15580772750
|
||||
1 1 1 1 1
|
||||
0.25183282652 1.0000000
|
||||
1 1 1 1 1
|
||||
0.83389127684E-01 1.0000000
|
||||
1 2 2 5 1
|
||||
45.485063360 0.32556415807E-02
|
||||
18.504059617 -0.54952972010E-02
|
||||
3.9156032308 0.27988806353
|
||||
1.7142196009 0.51273377761
|
||||
0.69675478242 0.33288802736
|
||||
1 2 2 1 1
|
||||
0.23060000000 1.0000000
|
||||
|
||||
## Effective core potentials
|
||||
def2-SVP_ECP
|
||||
Sb nelec 28
|
||||
Sb ul
|
||||
2 14.44497800 -15.36680100
|
||||
2 14.44929500 -20.29613800
|
||||
Sb S
|
||||
2 16.33086500 281.07158100
|
||||
2 8.55654200 61.71660400
|
||||
2 14.44497800 15.36680100
|
||||
2 14.44929500 20.29613800
|
||||
Sb P
|
||||
2 14.47033700 67.45738000
|
||||
2 13.81619400 134.93350300
|
||||
2 8.42492400 14.71634400
|
||||
2 8.09272800 29.51851200
|
||||
2 14.44497800 15.36680100
|
||||
2 14.44929500 20.29613800
|
||||
Sb D
|
||||
2 14.88633100 35.44781500
|
||||
2 15.14631900 53.14346600
|
||||
2 5.90826700 9.17922300
|
||||
2 5.59432200 13.24025300
|
||||
2 14.44497800 15.36680100
|
||||
2 14.44929500 20.29613800
|
||||
END def2-SVP_ECP
|
||||
|
||||
## All-electron potential
|
||||
H ALLELECTRON ALL
|
||||
1 0 0
|
||||
0.20000000 0
|
||||
|
||||
# Chlorine LANL2DZ (3s,3p) -> [2s,2p]
|
||||
Cl LANL2DZ
|
||||
2
|
||||
1 0 0 3 2
|
||||
2.2310000 -0.4900589 0.0000000
|
||||
0.4720000 1.2542684 0.0000000
|
||||
0.1631000 0.0000000 1.0000000
|
||||
1 1 1 3 2
|
||||
6.2960000 -0.0635641 0.0000000
|
||||
0.6333000 1.0141355 0.0000000
|
||||
0.1819000 0.0000000 1.0000000
|
||||
|
||||
# Iodine LANL2DZ (3s,3p) -> [2s,2p]
|
||||
I LANL2DZ
|
||||
2
|
||||
1 0 0 3 2
|
||||
0.7242000 -2.9731048 0.0000000
|
||||
0.4653000 3.4827643 0.0000000
|
||||
0.1336000 0.0000000 1.0000000
|
||||
1 1 1 3 2
|
||||
1.2900000 -0.2092377 0.0000000
|
||||
0.3180000 1.1035347 0.0000000
|
||||
0.1053000 0.0000000 1.0000000
|
||||
|
||||
## Effective core potentials
|
||||
LANL2DZ_ECP
|
||||
Cl nelec 10
|
||||
Cl ul
|
||||
1 94.8130000 -10.0000000
|
||||
2 165.6440000 66.2729170
|
||||
2 30.8317000 -28.9685950
|
||||
2 10.5841000 -12.8663370
|
||||
2 3.7704000 -1.7102170
|
||||
Cl S
|
||||
0 128.8391000 3.0000000
|
||||
1 120.3786000 12.8528510
|
||||
2 63.5622000 275.6723980
|
||||
2 18.0695000 115.6777120
|
||||
2 3.8142000 35.0606090
|
||||
Cl P
|
||||
0 216.5263000 5.0000000
|
||||
1 46.5723000 7.4794860
|
||||
2 147.4685000 613.0320000
|
||||
2 48.9869000 280.8006850
|
||||
2 13.2096000 107.8788240
|
||||
2 3.1831000 15.3439560
|
||||
I nelec 46
|
||||
I ul
|
||||
0 1.0715702 -0.0747621
|
||||
1 44.1936028 -30.0811224
|
||||
2 12.9367609 -75.3722721
|
||||
2 3.1956412 -22.0563758
|
||||
2 0.8589806 -1.6979585
|
||||
I S
|
||||
0 127.9202670 2.9380036
|
||||
1 78.6211465 41.2471267
|
||||
2 36.5146237 287.8680095
|
||||
2 9.9065681 114.3758506
|
||||
2 1.9420086 37.6547714
|
||||
I P
|
||||
0 13.0035304 2.2222630
|
||||
1 76.0331404 39.4090831
|
||||
2 24.1961684 177.4075002
|
||||
2 6.4053433 77.9889462
|
||||
2 1.5851786 25.7547641
|
||||
I D
|
||||
0 40.4278108 7.0524360
|
||||
1 28.9084375 33.3041635
|
||||
2 15.6268936 186.9453875
|
||||
2 4.1442856 71.9688361
|
||||
2 0.9377235 9.3630657
|
||||
END LANL2DZ_ECP
|
||||
58
tests/QS/regtest-ecp-2/ICl_lanl2dz_gpw.inp
Normal file
58
tests/QS/regtest-ecp-2/ICl_lanl2dz_gpw.inp
Normal file
|
|
@ -0,0 +1,58 @@
|
|||
&GLOBAL
|
||||
PROJECT ICl_lanl2dz_gpw
|
||||
RUN_TYPE DEBUG
|
||||
&END GLOBAL
|
||||
|
||||
&DEBUG
|
||||
CHECK_ATOM_FORCE 1 Z
|
||||
&END DEBUG
|
||||
|
||||
&FORCE_EVAL
|
||||
&DFT
|
||||
BASIS_SET_FILE_NAME ./ECP_BASIS_POT
|
||||
POTENTIAL_FILE_NAME ./ECP_BASIS_POT
|
||||
&MGRID
|
||||
CUTOFF 250
|
||||
NGRIDS 5
|
||||
REL_CUTOFF 40
|
||||
&END MGRID
|
||||
&POISSON
|
||||
PERIODIC NONE
|
||||
PSOLVER WAVELET
|
||||
&END POISSON
|
||||
&QS
|
||||
EPS_DEFAULT 1.0E-12
|
||||
METHOD GPW
|
||||
&END QS
|
||||
&SCF
|
||||
MAX_SCF 50
|
||||
SCF_GUESS ATOMIC
|
||||
&END SCF
|
||||
&XC
|
||||
&XC_FUNCTIONAL LDA
|
||||
&END XC_FUNCTIONAL
|
||||
&END XC
|
||||
&END DFT
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC 8.5 8.5 8.5
|
||||
PERIODIC NONE
|
||||
&END CELL
|
||||
&COORD
|
||||
Cl 0.00000 0.00000 0.00000
|
||||
I 0.00000 0.00000 2.40000
|
||||
&END COORD
|
||||
&KIND I
|
||||
BASIS_SET LANL2DZ
|
||||
POTENTIAL ECP LANL2DZ_ECP
|
||||
&END KIND
|
||||
&KIND Cl
|
||||
BASIS_SET LANL2DZ
|
||||
POTENTIAL ECP LANL2DZ_ECP
|
||||
&END KIND
|
||||
&TOPOLOGY
|
||||
&CENTER_COORDINATES
|
||||
&END CENTER_COORDINATES
|
||||
&END TOPOLOGY
|
||||
&END SUBSYS
|
||||
&END FORCE_EVAL
|
||||
47
tests/QS/regtest-ecp-2/Rn_stuttgart_gapw.inp
Normal file
47
tests/QS/regtest-ecp-2/Rn_stuttgart_gapw.inp
Normal file
|
|
@ -0,0 +1,47 @@
|
|||
&GLOBAL
|
||||
PROJECT Rn_stuttgart_gapw
|
||||
RUN_TYPE ENERGY_FORCE
|
||||
&END GLOBAL
|
||||
|
||||
&FORCE_EVAL
|
||||
STRESS_TENSOR ANALYTICAL
|
||||
&DFT
|
||||
BASIS_SET_FILE_NAME ./ECP_BASIS_POT
|
||||
POTENTIAL_FILE_NAME ./ECP_BASIS_POT
|
||||
&MGRID
|
||||
CUTOFF 300
|
||||
NGRIDS 5
|
||||
REL_CUTOFF 40
|
||||
&END MGRID
|
||||
&QS
|
||||
EPS_DEFAULT 1.0E-12
|
||||
METHOD GAPW
|
||||
&END QS
|
||||
&SCF
|
||||
IGNORE_CONVERGENCE_FAILURE
|
||||
MAX_SCF 5
|
||||
SCF_GUESS ATOMIC
|
||||
&END SCF
|
||||
&XC
|
||||
&XC_FUNCTIONAL LDA
|
||||
&END XC_FUNCTIONAL
|
||||
&END XC
|
||||
&END DFT
|
||||
&PRINT
|
||||
&STRESS_TENSOR
|
||||
&END STRESS_TENSOR
|
||||
&END PRINT
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC 6.0 6.0 12.0
|
||||
&END CELL
|
||||
&COORD
|
||||
Rn 0.00000 0.00000 0.00000
|
||||
Rn 0.00000 0.00000 5.50000
|
||||
&END COORD
|
||||
&KIND Rn
|
||||
BASIS_SET Stuttgart-RLC
|
||||
POTENTIAL ECP Stuttgart_RLC_ECP
|
||||
&END KIND
|
||||
&END SUBSYS
|
||||
&END FORCE_EVAL
|
||||
63
tests/QS/regtest-ecp-2/SbH3_def2_gapw.inp
Normal file
63
tests/QS/regtest-ecp-2/SbH3_def2_gapw.inp
Normal file
|
|
@ -0,0 +1,63 @@
|
|||
&GLOBAL
|
||||
PROJECT SbH3_def2_gapw
|
||||
RUN_TYPE GEO_OPT
|
||||
&END GLOBAL
|
||||
|
||||
&MOTION
|
||||
&GEO_OPT
|
||||
MAX_ITER 1
|
||||
&END GEO_OPT
|
||||
&END MOTION
|
||||
|
||||
&FORCE_EVAL
|
||||
&DFT
|
||||
BASIS_SET_FILE_NAME ./ECP_BASIS_POT
|
||||
POTENTIAL_FILE_NAME ./ECP_BASIS_POT
|
||||
&MGRID
|
||||
CUTOFF 300
|
||||
NGRIDS 5
|
||||
REL_CUTOFF 40
|
||||
&END MGRID
|
||||
&POISSON
|
||||
PERIODIC NONE
|
||||
PSOLVER WAVELET
|
||||
&END POISSON
|
||||
&QS
|
||||
EPS_DEFAULT 1.0E-12
|
||||
METHOD GAPW
|
||||
&END QS
|
||||
&SCF
|
||||
IGNORE_CONVERGENCE_FAILURE
|
||||
MAX_SCF 5
|
||||
SCF_GUESS ATOMIC
|
||||
&END SCF
|
||||
&XC
|
||||
&XC_FUNCTIONAL LDA
|
||||
&END XC_FUNCTIONAL
|
||||
&END XC
|
||||
&END DFT
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC 9.0 9.0 9.0
|
||||
PERIODIC NONE
|
||||
&END CELL
|
||||
&COORD
|
||||
Sb 0.0500000000 -0.2357770000 -0.6270200000
|
||||
H 0.0000000000 1.1518460000 0.3391280000
|
||||
H -1.2017170000 -0.9295890000 0.3391280000
|
||||
H 0.2017170000 -0.9295890000 0.3391280000
|
||||
&END COORD
|
||||
&KIND Sb
|
||||
BASIS_SET def2-SVP
|
||||
POTENTIAL ECP def2-SVP_ECP
|
||||
&END KIND
|
||||
&KIND H
|
||||
BASIS_SET def2-SVP
|
||||
POTENTIAL ALL
|
||||
&END KIND
|
||||
&TOPOLOGY
|
||||
&CENTER_COORDINATES
|
||||
&END CENTER_COORDINATES
|
||||
&END TOPOLOGY
|
||||
&END SUBSYS
|
||||
&END FORCE_EVAL
|
||||
3
tests/QS/regtest-ecp-2/TEST_FILES
Normal file
3
tests/QS/regtest-ecp-2/TEST_FILES
Normal file
|
|
@ -0,0 +1,3 @@
|
|||
ICl_lanl2dz_gpw.inp 0
|
||||
Rn_stuttgart_gapw.inp 31 1.0E-08 3.56470376649E-01
|
||||
SbH3_def2_gapw.inp 11 1.0E-11 -241.595963788403253
|
||||
|
|
@ -4,6 +4,7 @@
|
|||
# in case a new directory is added just add it at the top of the list..
|
||||
# the order will be regularly checked and modified...
|
||||
QS/regtest-ecp libgrpp
|
||||
QS/regtest-ecp-2 libgrpp
|
||||
QS/regtest-as-3 libint mpiranks%2==0
|
||||
QS/regtest-as-2 libint
|
||||
QS/regtest-wfn-restart
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue