diff --git a/src/core_ppl.F b/src/core_ppl.F
index 94dc1fbfaf..29ca85f7e6 100644
--- a/src/core_ppl.F
+++ b/src/core_ppl.F
@@ -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
diff --git a/src/libgrpp_integrals.F b/src/libgrpp_integrals.F
index 3714ad637a..0e82388707 100644
--- a/src/libgrpp_integrals.F
+++ b/src/libgrpp_integrals.F
@@ -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/dPx = d/dAx + d/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/dPx = d/dAx + d/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
diff --git a/tests/QS/regtest-ecp-2/ECP_BASIS_POT b/tests/QS/regtest-ecp-2/ECP_BASIS_POT
new file mode 100644
index 0000000000..de04364c09
--- /dev/null
+++ b/tests/QS/regtest-ecp-2/ECP_BASIS_POT
@@ -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
diff --git a/tests/QS/regtest-ecp-2/ICl_lanl2dz_gpw.inp b/tests/QS/regtest-ecp-2/ICl_lanl2dz_gpw.inp
new file mode 100644
index 0000000000..b9419feac5
--- /dev/null
+++ b/tests/QS/regtest-ecp-2/ICl_lanl2dz_gpw.inp
@@ -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
diff --git a/tests/QS/regtest-ecp-2/Rn_stuttgart_gapw.inp b/tests/QS/regtest-ecp-2/Rn_stuttgart_gapw.inp
new file mode 100644
index 0000000000..4660809934
--- /dev/null
+++ b/tests/QS/regtest-ecp-2/Rn_stuttgart_gapw.inp
@@ -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
diff --git a/tests/QS/regtest-ecp-2/SbH3_def2_gapw.inp b/tests/QS/regtest-ecp-2/SbH3_def2_gapw.inp
new file mode 100644
index 0000000000..1470842d3e
--- /dev/null
+++ b/tests/QS/regtest-ecp-2/SbH3_def2_gapw.inp
@@ -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
diff --git a/tests/QS/regtest-ecp-2/TEST_FILES b/tests/QS/regtest-ecp-2/TEST_FILES
new file mode 100644
index 0000000000..7d82a01f21
--- /dev/null
+++ b/tests/QS/regtest-ecp-2/TEST_FILES
@@ -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
diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS
index 58305d816f..95cfcb6b35 100644
--- a/tests/TEST_DIRS
+++ b/tests/TEST_DIRS
@@ -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