From 48d3695c3ced91faad61a1fd2c67ef2ebe356476 Mon Sep 17 00:00:00 2001 From: abussy Date: Fri, 15 Dec 2023 12:46:12 +0100 Subject: [PATCH] ECP nuclear gradients --- src/core_ppl.F | 59 +- src/libgrpp_integrals.F | 660 ++++++++++++++++++- tests/QS/regtest-ecp-2/ECP_BASIS_POT | 203 ++++++ tests/QS/regtest-ecp-2/ICl_lanl2dz_gpw.inp | 58 ++ tests/QS/regtest-ecp-2/Rn_stuttgart_gapw.inp | 47 ++ tests/QS/regtest-ecp-2/SbH3_def2_gapw.inp | 63 ++ tests/QS/regtest-ecp-2/TEST_FILES | 3 + tests/TEST_DIRS | 1 + 8 files changed, 1053 insertions(+), 41 deletions(-) create mode 100644 tests/QS/regtest-ecp-2/ECP_BASIS_POT create mode 100644 tests/QS/regtest-ecp-2/ICl_lanl2dz_gpw.inp create mode 100644 tests/QS/regtest-ecp-2/Rn_stuttgart_gapw.inp create mode 100644 tests/QS/regtest-ecp-2/SbH3_def2_gapw.inp create mode 100644 tests/QS/regtest-ecp-2/TEST_FILES 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