From 31b9a200d467d15decd377fb1c27b3cd35b14f68 Mon Sep 17 00:00:00 2001 From: Teodoro Laino Date: Thu, 1 May 2008 18:52:48 +0000 Subject: [PATCH] Analytical derivatives for d-orbitals! added few more regtests together with a "real" chemistry example. svn-origin-rev: 7331 --- src/semi_empirical_int_ana.F | 194 +++++++++++++++++++--- src/semi_empirical_int_debug.F | 6 +- src/semi_empirical_int_num.F | 14 +- src/semi_empirical_int_utils.F | 207 +++++++++++++++++++++++- tests/SE/regtest-2/O2-ROKS-PM6.inp | 61 +++++++ tests/SE/regtest-2/Pt-cis-2xpet3Cl2.inp | 103 ++++++++++++ tests/SE/regtest-2/TEST_FILES | 6 + tests/SE/regtest-2/hcl.inp | 66 ++++++++ tests/SE/regtest-2/hcl_ot.inp | 68 ++++++++ tests/SE/regtest-2/hcl_z.inp | 66 ++++++++ 10 files changed, 758 insertions(+), 33 deletions(-) create mode 100644 tests/SE/regtest-2/O2-ROKS-PM6.inp create mode 100644 tests/SE/regtest-2/Pt-cis-2xpet3Cl2.inp create mode 100644 tests/SE/regtest-2/hcl.inp create mode 100644 tests/SE/regtest-2/hcl_ot.inp create mode 100644 tests/SE/regtest-2/hcl_z.inp diff --git a/src/semi_empirical_int_ana.F b/src/semi_empirical_int_ana.F index fc5c95e597..7bdeed9c68 100644 --- a/src/semi_empirical_int_ana.F +++ b/src/semi_empirical_int_ana.F @@ -27,9 +27,12 @@ MODULE semi_empirical_int_ana USE semi_empirical_int_arrays, ONLY: & fac_x_to_z, ijkl_ind, ijkl_sym, inddd, inddp, indexa, indexb, indpp, & int2c_type, l_index, map_x_to_z, rij_threshold - USE semi_empirical_int_num, ONLY: nucint_sp_num,& + USE semi_empirical_int_num, ONLY: nucint_d_num,& + nucint_sp_num,& + terep_d_num,& terep_sp_num - USE semi_empirical_int_utils, ONLY: d_ijkl_sp,& + USE semi_empirical_int_utils, ONLY: d_ijkl_d,& + d_ijkl_sp,& dtaper_eval,& rot_2el_2c_first,& rotmat,& @@ -622,7 +625,7 @@ CONTAINS DO k = 1, sepj%natorb DO l = 1, k ndim = ndim + 1 - + ! Get the integral in the original frame (along z) DO ind = 1, 2 imap = map_x_to_z(ind,i) @@ -646,7 +649,7 @@ CONTAINS gknd = indexb(kmap,lmap) tind = (gind-1)*tdim+gknd - tmp2el(ndim) = tmp2el(ndim) + ifac*jfac*int2el(tind) + tmp2el(ndim) = tmp2el(ndim) + ifac*jfac*lfac*kfac*int2el(tind) END DO END DO @@ -773,14 +776,16 @@ CONTAINS dssss = 0.0_dp END IF - ! Value of the integrals + ! Value of the integrals for sp shell CALL nucint_sp_num(sepi, sepj, rij, ssss, core, itype, longrange, error=error) - ! Integrals derivatives + ! Integrals derivatives for sp shell CALL dnucint_sp_ana(sepi, sepj, rij,dssss, dcore, itype, longrange, error=error) - IF (sepi%dorb.OR.sepj%dorb) THEN - ! compute the contribution from d-orbitals - STOP "D-ORBITALS: dnucint_ana not working with analytical derivatives" + IF ((sepi%dorb.OR.sepj%dorb).AND.l_core) THEN + ! Compute the contribution from d-orbitals + CALL nucint_d_num(sepi, sepj, rij, core, itype, longrange, error=error) + ! Compute the contribution from d-orbitals + CALL dnucint_d_ana(sepi, sepj, rij,dcore, itype, longrange, error=error) END IF IF (.NOT.longrange) THEN @@ -847,7 +852,9 @@ CONTAINS !> ssss = derivative of (ssss) integral !> derivatives are intended w.r.t. rij !> The storage of the nuclear attraction integrals core(kl/ij) iS -!> (SS/)=1, (SO/)=2, (OO/)=3, (PP/)=4 +!> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5, +!> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10 +!> !> where ij=1 if the orbitals centred on atom i, =2 if on atom j. !> \note !> Analytical version - Analytical evaluation of gradients @@ -926,6 +933,94 @@ CONTAINS END SUBROUTINE dnucint_sp_ana +! ***************************************************************************** +!> \brief Calculates the analytical derivative of the nuclear attraction +!> integrals involving d orbitals +!> \param on input rij = interatomic distance +!> sepi = paramters of atom i +!> sepj = paramters of atom j +!> itype = type of semi_empirical model +!> extension to the original routine to compute qm/mm integrals +!> short_range = requests only the evaluation of the short range part +!> longrange = requests only the evaluation of the long range part +!> on output core = 4 X 2 array of electron-core attraction integrals +!> +!> The storage of the nuclear attraction integrals core(kl/ij) iS +!> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5, +!> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10 +!> +!> where ij=1 if the orbitals centred on atom i, =2 if on atom j. +!> \author +!> Teodoro Laino (05.2008) [tlaino] - University of Zurich: created +! ***************************************************************************** + SUBROUTINE dnucint_d_ana( sepi, sepj, rij, dcore, itype, longrange, error) + TYPE(semi_empirical_type), POINTER :: sepi, sepj + REAL(dp), INTENT(IN) :: rij + REAL(dp), DIMENSION(10, 2), & + INTENT(INOUT) :: dcore + INTEGER, INTENT(IN) :: itype + LOGICAL, INTENT(IN) :: longrange + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'dnucint_d_ana', & + routineP = moduleN//':'//routineN + + INTEGER :: i, ij, kl + LOGICAL :: failure + + failure = .FALSE. + ! Check if d-orbitals are present + IF (sepi%dorb.OR.sepj%dorb) THEN + DO i = 5, 10 + dcore(i, 1) = 0.0_dp + dcore(i, 2) = 0.0_dp + END DO + ! + ij = indexa(1, 1) + IF (sepj%dorb) THEN + ! + kl = indexa(5, 1) + dcore(5, 2) = d_ijkl_d (sepi, sepj, ij, kl, 0, 0, 2, 0, 1, rij, longrange, itype, error) * sepi%zeff + ! + kl = indexa(5, 2) + dcore(6, 2) = d_ijkl_d (sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, longrange, itype, error) * sepi%zeff + ! + kl = indexa(5, 5) + dcore(7, 2) = d_ijkl_d (sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, longrange, itype, error) * sepi%zeff + ! + kl = indexa(6, 3) + dcore(8, 2) = d_ijkl_d (sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, longrange, itype, error) * sepi%zeff + ! + kl = indexa(6, 6) + dcore(9, 2) = d_ijkl_d (sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, longrange, itype, error) * sepi%zeff + ! + kl = indexa(8, 8) + dcore(10, 2)= d_ijkl_d (sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, longrange, itype, error) * sepi%zeff + END IF + IF (sepi%dorb) THEN + ! + kl = indexa(5, 1) + dcore(5, 1) = d_ijkl_d (sepi, sepj, kl, ij, 2, 0, 0, 0, 2, rij, longrange, itype, error) * sepj%zeff + ! + kl = indexa(5, 2) + dcore(6, 1) = d_ijkl_d (sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, longrange, itype, error) * sepj%zeff + ! + kl = indexa(5, 5) + dcore(7, 1) = d_ijkl_d (sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, longrange, itype, error) * sepj%zeff + ! + kl = indexa(6, 3) + dcore(8, 1) = d_ijkl_d (sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, longrange, itype, error) * sepj%zeff + ! + kl = indexa(6, 6) + dcore(9, 1) = d_ijkl_d (sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, longrange, itype, error) * sepj%zeff + ! + kl = indexa(8, 8) + dcore(10, 1)= d_ijkl_d (sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, longrange, itype, error) * sepj%zeff + END IF + + END IF + END SUBROUTINE dnucint_d_ana + ! ***************************************************************************** !> \brief calculates the derivative of the two-particle interactions !> \param on input sepi = Atomic parameters of first atom @@ -1268,9 +1363,11 @@ CONTAINS ! Integral Derivatives CALL dterep_sp_ana(sepi, sepj, rij, rep_d(1:34), longrange, error=error) - IF (sepi%dorb.OR.sepi%dorb) THEN - ! compute the contribution from d-orbitals - STOP "D-ORBITALS: dterep_ana not working with analytical derivatives" + IF (sepi%dorb.OR.sepj%dorb) THEN + ! Compute the contribution from d-orbitals + CALL terep_d_num (sepi, sepj, rij, rep, longrange, error=error ) + ! Compute the derivatives from d-orbitals + CALL dterep_d_ana (sepi, sepj, rij, rep_d, longrange, error=error ) END IF ! Tapering function @@ -1298,19 +1395,12 @@ CONTAINS ! ***************************************************************************** !> \brief Calculates the derivative and the value of two-electron repulsion -!> integrals and the nuclear attraction integrals w.r.t. |r| - sp core only +!> integrals and the nuclear attraction integrals w.r.t. |r| - sp shells only !> \param on input rij = interatomic distance !> sepi = paramters of atom i !> sepj = paramters of atom j !> on output ri = array of two-electron repulsion integrals !> dri = array of derivatives of two-electron repulsion integrals -!> The two-centre repulsion integrals (over local coordinates) are -!> stored as follows (where p-sigma = O, and p-pi = P and P* ) -!> (SS/SS)=1, (SO/SS)=2, (OO/SS)=3, (PP/SS)=4, (SS/OS)=5, -!> (SO/SO)=6, (SP/SP)=7, (OO/SO)=8, (PP/SO)=9, (PO/SP)=10, -!> (SS/OO)=11, (SS/PP)=12, (SO/OO)=13, (SO/PP)=14, (SP/OP)=15, -!> (OO/OO)=16, (PP/OO)=17, (OO/PP)=18, (PP/PP)=19, (PO/PO)=20, -!> (PP/P*P*)=21, (P*P/P*P)=22. !> \note !> Analytical version - Analytical evaluation of gradients !> Teodoro Laino - Zurich University 04.2007 @@ -1371,4 +1461,66 @@ CONTAINS END SUBROUTINE dterep_sp_ana +! ***************************************************************************** +!> \brief Calculates the derivatives of the two-electron repulsion integrals - d shell only +!> \param on input rij = interatomic distance +!> sepi = paramters of atom i +!> sepj = paramters of atom j +!> longrange = requests only the evaluation of the long range part +!> +!> on output rep = array of two-electron repulsion integrals +!> +!> \par History +!> Teodoro Laino (05.2008) [tlaino] - University of Zurich : new driver +!> for computing integral derivatives for d-orbitals +! ***************************************************************************** + SUBROUTINE dterep_d_ana (sepi, sepj, rij, drep, longrange, error ) + TYPE(semi_empirical_type), POINTER :: sepi, sepj + REAL(dp), INTENT(IN) :: rij + REAL(dp), DIMENSION(491), INTENT(INOUT) :: drep + LOGICAL, INTENT(IN) :: longrange + TYPE(cp_error_type), INTENT(inout) :: error + + CHARACTER(len=*), PARAMETER :: routineN = 'dterep_d_ana', & + routineP = moduleN//':'//routineN + + INTEGER :: i, ij, j, k, kl, l, lasti, & + lastj, li, lj, lk, ll, nold, & + numb + LOGICAL :: failure + + failure = .FALSE. + lasti = sepi%natorb + lastj = sepj%natorb + DO i = 1, lasti + li = l_index(i) + DO j = 1, i + lj = l_index(j) + ij = indexa(i, j) + DO k = 1, lastj + lk = l_index(k) + DO l = 1, k + ll = l_index(l) + kl = indexa(k, l) + + numb = ijkl_ind(ij, kl) + ! From 1 to 34 we store integrals involving sp shells + IF (numb > 34) THEN + nold = ijkl_sym(numb) + IF (nold > 34) THEN + drep(numb) = drep( nold) + ELSE IF (nold <-34) THEN + drep(numb) = -drep(-nold) + ELSE IF (nold == 0) THEN + drep(numb) = d_ijkl_d (sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij,& + longrange, do_method_undef, error) + END IF + END IF + END DO + END DO + END DO + END DO + + END SUBROUTINE dterep_d_ana + END MODULE semi_empirical_int_ana diff --git a/src/semi_empirical_int_debug.F b/src/semi_empirical_int_debug.F index 8557e95ee0..8a18cce1fd 100644 --- a/src/semi_empirical_int_debug.F +++ b/src/semi_empirical_int_debug.F @@ -498,7 +498,7 @@ SUBROUTINE check_dterep_ana (sepi,sepj,r,ri,dri,shortrange,taper,error) routineN = 'check_dterep_ana', routineP = moduleN//':'//routineN LOGICAL :: check_value - INTEGER :: j + INTEGER :: j, i REAL(dp) :: delta, od, rn REAL(dp), DIMENSION(491) :: nri, ri0, rim, rip @@ -520,6 +520,10 @@ SUBROUTINE check_dterep_ana (sepi,sepj,r,ri,dri,shortrange,taper,error) END IF IF (.NOT.check_value(dri(j), nri(j), delta*100.0_dp, 0.1_dp)) THEN WRITE(*,*)"ERROR for derivative of RI integral, RI(j), j::",j + WRITE(*,*)"FULL SET OF INTEGRALS: INDX ANAL NUM DIFF" + DO i = 1, 491 + WRITE(*,'(I5,3F15.9)')i, dri(i), nri(i), nri(i)-dri(i) + END DO STOP END IF END DO diff --git a/src/semi_empirical_int_num.F b/src/semi_empirical_int_num.F index 631c687cb1..8896327618 100644 --- a/src/semi_empirical_int_num.F +++ b/src/semi_empirical_int_num.F @@ -48,7 +48,7 @@ MODULE semi_empirical_int_num PUBLIC :: rotint_num, rotnuc_num, corecore_num,& drotint_num, drotnuc_num, dcorecore_num,& nucint_num, terep_num, nucint_sp_num,& - terep_sp_num + terep_sp_num, terep_d_num, nucint_d_num CONTAINS @@ -202,7 +202,7 @@ CONTAINS limij = (ii*(ii+1))/2 limkl = (kk*(kk+1))/2 CALL store_2el_2c_diag(limij, limkl, ww, w, error=error) - ENDIF + ENDIF END SUBROUTINE rotint_num ! ***************************************************************************** @@ -782,7 +782,9 @@ CONTAINS !> ssss = (ssss) integral !> derivatives are intended w.r.t. rij !> The storage of the nuclear attraction integrals core(kl/ij) iS -!> (SS/)=1, (SO/)=2, (OO/)=3, (PP/)=4 +!> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5, +!> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10 +!> !> where ij=1 if the orbitals centred on atom i, =2 if on atom j. !> \par History !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver @@ -864,11 +866,11 @@ CONTAINS !> short_range = requests only the evaluation of the short range part !> longrange = requests only the evaluation of the long range part !> on output core = 4 X 2 array of electron-core attraction integrals -!> ssss = (ssss) integral -!> derivatives are intended w.r.t. rij !> !> The storage of the nuclear attraction integrals core(kl/ij) iS -!> (SS/)=1, (SO/)=2, (OO/)=3, (PP/)=4 +!> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5, +!> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10 +!> !> where ij=1 if the orbitals centred on atom i, =2 if on atom j. !> \author !> Teodoro Laino (03.2008) [tlaino] - University of Zurich diff --git a/src/semi_empirical_int_utils.F b/src/semi_empirical_int_utils.F index d86340e38a..1682c9b4b0 100644 --- a/src/semi_empirical_int_utils.F +++ b/src/semi_empirical_int_utils.F @@ -29,7 +29,7 @@ MODULE semi_empirical_int_utils CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_utils' PUBLIC :: taper_eval, dtaper_eval, ijkl_d, ijkl_sp, rotmat, rot_2el_2c_first,& - store_2el_2c_diag, d_ijkl_sp + store_2el_2c_diag, d_ijkl_sp, d_ijkl_d CONTAINS @@ -626,8 +626,9 @@ CONTAINS END FUNCTION dcharg_int_nri ! ***************************************************************************** -!> \brief General driver for computing semi-empirical integrals involving -!> d-orbitals. The choice of the linear quadrupole was REALLY unhappy +!> \brief General driver for computing semi-empirical integrals +!> involving d-orbitals. +!> The choice of the linear quadrupole was REALLY unhappy !> in the first development of the NDDO codes. That choice makes !> impossible the merging of the integral code with or without d-orbitals !> unless a reparametrization of all NDDO codes for s and p orbitals will @@ -652,6 +653,71 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'ijkl_d', & routineP = moduleN//':'//routineN + res = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, longrange, itype, charg_int_ri, error) + + END FUNCTION ijkl_d + +! ***************************************************************************** +!> \brief General driver for computing the derivatives of semi-empirical integrals +!> involving d-orbitals. +!> The choice of the linear quadrupole was REALLY unhappy +!> in the first development of the NDDO codes. That choice makes +!> impossible the merging of the integral code with or without d-orbitals +!> unless a reparametrization of all NDDO codes for s and p orbitals will +!> be performed.. more over the choice of the linear quadrupole does not make +!> calculations rotational invariants (of course the rotational invariant +!> can be forced). The definitions of quadrupoles for d-orbitals is the +!> correct one in order to have the rotational invariant property by +!> construction.. +!> +!> \author Teodoro Laino [tlaino] - University of Zurich +!> \date 03.2008 [tlaino] +! ***************************************************************************** + FUNCTION d_ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, longrange, itype, error) RESULT(res) + TYPE(semi_empirical_type), POINTER :: sepi, sepj + INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic + REAL(KIND=dp), INTENT(IN) :: r + LOGICAL, INTENT(IN) :: longrange + INTEGER, INTENT(IN) :: itype + TYPE(cp_error_type), INTENT(inout) :: error + REAL(KIND=dp) :: res + + CHARACTER(len=*), PARAMETER :: routineN = 'd_ijkl_d', & + routineP = moduleN//':'//routineN + + res = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, longrange, itype, dcharg_int_ri, error) + + END FUNCTION d_ijkl_d + +! ***************************************************************************** +!> \brief Low level general driver for computing semi-empirical integrals +!> and their derivatives involving d-orbitals. +!> The choice of the linear quadrupole was REALLY unhappy +!> in the first development of the NDDO codes. That choice makes +!> impossible the merging of the integral code with or without d-orbitals +!> unless a reparametrization of all NDDO codes for s and p orbitals will +!> be performed.. more over the choice of the linear quadrupole does not make +!> calculations rotational invariants (of course the rotational invariant +!> can be forced). The definitions of quadrupoles for d-orbitals is the +!> correct one in order to have the rotational invariant property by +!> construction.. +!> +!> \author Teodoro Laino [tlaino] - University of Zurich +!> \date 03.2008 [tlaino] +! ***************************************************************************** + FUNCTION ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, longrange, itype, eval, error) RESULT(res) + TYPE(semi_empirical_type), POINTER :: sepi, sepj + INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic + REAL(KIND=dp), INTENT(IN) :: r + LOGICAL, INTENT(IN) :: longrange + INTEGER, INTENT(IN) :: itype + REAL(KIND=dp) :: eval + TYPE(cp_error_type), INTENT(inout) :: error + REAL(KIND=dp) :: res + + CHARACTER(len=*), PARAMETER :: routineN = 'ijkl_d_low', & + routineP = moduleN//':'//routineN + INTEGER :: l1, l1max, l1min, l2, l2max, & l2min, lij, lkl, lmin, m, mm REAL(KIND=dp) :: add, ccc, chrg, dij, dkl, & @@ -715,7 +781,7 @@ CONTAINS ccc = clm_d(ij, l1, m) * clm_d(kl, l2, m) IF (ABS(ccc) > EPSILON(0.0_dp)) THEN mm = ABS (m) - chrg = charg_int_ri(r, l1, l2, mm, dij, dkl, add, error) + chrg = eval(r, l1, l2, mm, dij, dkl, add, error) s1 = s1 + chrg * ccc END IF END DO @@ -723,7 +789,7 @@ CONTAINS END DO END DO res = sum - END FUNCTION ijkl_d + END FUNCTION ijkl_d_low ! ***************************************************************************** !> \brief Interaction function between two point-charge configurations (MNDO/d) @@ -854,6 +920,137 @@ CONTAINS CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) END FUNCTION charg_int_ri +! ***************************************************************************** +!> \brief Derivatives of the interaction function between two point-charge +!> configurations (MNDO/d) +!> Rotational invariant property built-in in the quadrupole definition +!> r - Distance r12 +!> l1,m - Quantum numbers for multipole of configuration 1 +!> l2,m - Quantum numbers for multipole of configuration 2 +!> da - charge separation of configuration 1 +!> db - charge separation of configuration 2 +!> add - additive term +!> +!> \author Teodoro Laino [tlaino] - University of Zurich +!> \date 03.2008 [tlaino] +! ***************************************************************************** + FUNCTION dcharg_int_ri(r, l1_i, l2_i, m, da_i, db_i, add, error) RESULT(charg) + REAL(KIND=dp), INTENT(in) :: r + INTEGER, INTENT(in) :: l1_i, l2_i, m + REAL(KIND=dp), INTENT(in) :: da_i, db_i, add + TYPE(cp_error_type), INTENT(inout) :: error + REAL(KIND=dp) :: charg + + CHARACTER(len=*), PARAMETER :: routineN = 'dcharg_int_ri', & + routineP = moduleN//':'//routineN + + INTEGER :: l1, l2 + LOGICAL :: failure + REAL(KIND=dp) :: aa, ab, da, db, dxdx, dxqxz, & + dzdz, dzqzz, fact, qqzz, & + qxzqxz, xyxy, zzzz + + failure = .FALSE. + IF ( l1_i <= l2_i) THEN + l1 = l1_i + l2 = l2_i + da = da_i + db = db_i + fact = 1.0_dp + ELSE IF ( l1_i > l2_i ) THEN + l1 = l2_i + l2 = l1_i + da = db_i + db = da_i + fact = (-1.0_dp)**(l1+l2) + END IF + charg = 0.0_dp + ! Q - Q. + IF (l1 == 0 .AND. l2 == 0) THEN + charg = r/SQRT(r**2+add)**3 + charg = -fact*charg + RETURN + END IF + ! Q - Z. + IF (l1 == 0 .AND. l2 == 1) THEN + charg = (r+db)/SQRT((r+db)**2+add)**3 - (r-db)/SQRT((r-db)**2+add)**3 + charg = -charg*0.5_dp*fact + RETURN + END IF + ! Z - Z. + IF (l1 == 1 .AND. l2 == 1 .AND. m == 0) THEN + dzdz = & + +(r+da-db)/SQRT((r+da-db)**2+add)**3 + (r-da+db)/SQRT((r-da+db)**2+add)**3 & + -(r-da-db)/SQRT((r-da-db)**2+add)**3 - (r+da+db)/SQRT((r+da+db)**2+add)**3 + charg = -dzdz*0.25_dp*fact + RETURN + END IF + ! X - X + IF (l1 == 1 .AND. l2 == 1 .AND. m == 1) THEN + dxdx = 2.0_dp*r/SQRT(r**2+(da-db)**2+add)**3 - 2.0_dp*r/SQRT(r**2+(da+db)**2+add)**3 + charg = -dxdx*0.25_dp*fact + RETURN + END IF + ! Q - ZZ + IF (l1 == 0 .AND. l2 == 2) THEN + qqzz = (r-db)/SQRT((r-db)**2+add)**3 - 2.0_dp*r/SQRT(r**2+db**2+add)**3 + (r+db)/SQRT((r+db)**2+add)**3 + charg = -qqzz*0.25_dp*fact + RETURN + END IF + ! Z - ZZ + IF (l1 == 1 .AND. l2 == 2 .AND. m == 0) THEN + dzqzz = & + +(r-da-db)/SQRT((r-da-db)**2+add)**3 - 2.0_dp*(r-da)/SQRT((r-da)**2+db**2+add)**3 & + +(r+db-da)/SQRT((r+db-da)**2+add)**3 - (r-db+da)/SQRT((r-db+da)**2+add)**3 & + +2.0_dp*(r+da)/SQRT((r+da)**2+db**2+add)**3 - (r+da+db)/SQRT((r+da+db)**2+add)**3 + charg = -dzqzz*0.125_dp*fact + RETURN + END IF + ! ZZ - ZZ + IF (l1 == 2 .AND. l2 == 2 .AND. m == 0) THEN + zzzz = & + +(r-da-db)/SQRT((r-da-db)**2+add)**3 + (r+da+db)/SQRT((r+da+db)**2+add)**3 & + +(r-da+db)/SQRT((r-da+db)**2+add)**3 + (r+da-db)/SQRT((r+da-db)**2+add)**3 & + -2.0_dp*(r-da)/SQRT((r-da)**2+db**2+add)**3 - 2.0_dp*(r-db)/SQRT((r-db)**2+da**2+add)**3 & + -2.0_dp*(r+da)/SQRT((r+da)**2+db**2+add)**3 - 2.0_dp*(r+db)/SQRT((r+db)**2+da**2+add)**3 & + +2.0_dp*r/SQRT(r**2+(da-db)**2+add)**3 + 2.0_dp*r/SQRT(r**2+(da+db)**2+add)**3 + xyxy = & + +4.0_dp*r/SQRT(r**2+(da-db)**2+add)**3 + 4.0_dp*r/SQRT(r**2+(da+db)**2+add)**3 & + -8.0_dp*r/SQRT(r**2+da**2+db**2+add)**3 + charg = -(zzzz*0.0625_dp - xyxy*0.015625_dp)*fact + RETURN + END IF + ! X - ZX + IF (l1 == 1 .AND. l2 == 2 .AND. m == 1) THEN + ab = db/SQRT(2.0_dp) + dxqxz = & + -2.0_dp*(r-ab)/SQRT((r-ab)**2+(da-ab)**2+add)**3 + 2.0_dp*(r+ab)/SQRT((r+ab)**2+(da-ab)**2+add)**3 & + +2.0_dp*(r-ab)/SQRT((r-ab)**2+(da+ab)**2+add)**3 - 2.0_dp*(r+ab)/SQRT((r+ab)**2+(da+ab)**2+add)**3 + charg = -dxqxz*0.125_dp*fact + RETURN + END IF + ! ZX - ZX + IF (l1 == 2 .AND. l2 == 2 .AND. m == 1) THEN + aa = da/SQRT (2.0_dp) + ab = db/SQRT (2.0_dp) + qxzqxz = & + +2.0_dp*(r+aa-ab)/SQRT((r+aa-ab)**2+(aa-ab)**2+add)**3 - 2.0_dp*(r+aa+ab)/SQRT((r+aa+ab)**2+(aa-ab)**2+add)**3 & + -2.0_dp*(r-aa-ab)/SQRT((r-aa-ab)**2+(aa-ab)**2+add)**3 + 2.0_dp*(r-aa+ab)/SQRT((r-aa+ab)**2+(aa-ab)**2+add)**3 & + -2.0_dp*(r+aa-ab)/SQRT((r+aa-ab)**2+(aa+ab)**2+add)**3 + 2.0_dp*(r+aa+ab)/SQRT((r+aa+ab)**2+(aa+ab)**2+add)**3 & + +2.0_dp*(r-aa-ab)/SQRT((r-aa-ab)**2+(aa+ab)**2+add)**3 - 2.0_dp*(r-aa+ab)/SQRT((r-aa+ab)**2+(aa+ab)**2+add)**3 + charg = -qxzqxz*0.0625_dp*fact + RETURN + END IF + ! XX - XX + IF (l1 == 2 .AND. l2 == 2 .AND. m == 2) THEN + xyxy = 4.0_dp*r/SQRT(r**2+(da-db)**2+add)**3 + 4.0_dp*r/SQRT(r**2+(da+db)**2+add)**3 - 8.0_dp*r/SQRT(r**2+da**2+db**2+add)**3 + charg = -xyxy*0.0625_dp*fact + RETURN + END IF + ! We should NEVER reach this point + CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) + END FUNCTION dcharg_int_ri + ! ***************************************************************************** !> \brief Computes the general rotation matrices for the integrals !> between I and J (J>=I) diff --git a/tests/SE/regtest-2/O2-ROKS-PM6.inp b/tests/SE/regtest-2/O2-ROKS-PM6.inp new file mode 100644 index 0000000000..3eb18d00da --- /dev/null +++ b/tests/SE/regtest-2/O2-ROKS-PM6.inp @@ -0,0 +1,61 @@ +&GLOBAL + + PROJECT O2 + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL + +&MOTION + &GEO_OPT + OPTIMIZER BFGS + MAX_ITER 30 + MAX_DR 0.001 + RMS_DR 0.0005 + MAX_FORCE 0.00015 + RMS_FORCE 0.0001 + &END GEO_OPT + &PRINT + &STRUCTURE_DATA + DISTANCE 1 2 + &END STRUCTURE_DATA + &END +&END MOTION + +&FORCE_EVAL + METHOD Quickstep + &DFT + ROKS + MULTIP 3 + &QS + METHOD PM6 + &SE + + &END + + &END QS + &SCF + EPS_SCF 1.0E-6 + SCF_GUESS atomic + &END SCF + &END DFT + &SUBSYS + &CELL + ABC 10.0 10.0 10.0 + UNIT ANGSTROM + PERIODIC NONE + &END CELL + &COORD + O 0.000000 0.000000 -0.600000 + O 0.000000 0.000000 0.600000 + &END COORD + &PRINT + &STRUCTURE_DATA + DISTANCE 1 2 + &END STRUCTURE_DATA + &KINDS + BASIS_SET + POTENTIAL + &END KINDS + &END PRINT + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/SE/regtest-2/Pt-cis-2xpet3Cl2.inp b/tests/SE/regtest-2/Pt-cis-2xpet3Cl2.inp new file mode 100644 index 0000000000..b6e2d3bf6d --- /dev/null +++ b/tests/SE/regtest-2/Pt-cis-2xpet3Cl2.inp @@ -0,0 +1,103 @@ +&GLOBAL + PROJECT Pt-cmplx + RUN_TYPE GEO_OPT + PRINT_LEVEL LOW +&END GLOBAL + +&MOTION + &GEO_OPT + OPTIMIZER BFGS + MAX_ITER 1 + MAX_DR 0.01 + RMS_DR 0.005 + MAX_FORCE 0.0015 + RMS_FORCE 0.001 + &END GEO_OPT +&END MOTION + +&FORCE_EVAL + METHOD Quickstep + &DFT + &QS + METHOD PM6 + &SE + &END + &END QS + &SCF + MAX_SCF 20 + EPS_SCF 1.0E-4 + SCF_GUESS atomic + &OT + MINIMIZER DIIS + N_DIIS 7 + PRECONDITIONER NONE + &END + &OUTER_SCF + MAX_SCF 5 + EPS_SCF 1.0E-4 + &END + &END SCF + &END DFT + &SUBSYS + &CELL + ABC 30.0 30.0 30.0 + UNIT ANGSTROM + PERIODIC NONE + &END CELL + &COORD +Pt 5.664194 4.348452 6.089697 +Cl 5.664194 4.348452 8.432092 +Cl 7.963412 4.348452 5.646152 + P 3.475479 4.361633 6.698833 + P 5.863001 4.289642 3.829353 + C 3.079497 2.768742 7.500998 + C 2.146824 4.433308 5.436883 + C 2.947475 5.675961 7.848487 + C 6.606138 5.844605 3.230467 + C 4.333126 4.193223 2.823153 + C 6.882213 2.963555 3.111233 + C 3.872823 1.730794 6.749773 + H 3.395411 2.809009 8.564407 + H 1.999780 2.536631 7.504078 + H 2.367198 3.695945 4.640643 + C 2.008353 5.821781 4.874901 + H 1.178219 4.100756 5.857728 + C 1.580673 5.457600 8.429502 + H 2.997252 6.645945 7.314322 + H 3.695390 5.752483 8.669544 + H 7.710505 5.775320 3.313801 + H 6.388757 6.052665 2.168119 + C 6.070582 6.928852 4.131128 + H 3.777599 5.141109 2.972847 + C 4.530970 3.969292 1.350200 + H 3.690088 3.388023 3.232376 + C 6.299818 1.633085 3.501303 + H 6.960889 3.064331 2.013199 + H 7.923146 3.056503 3.490605 + H 3.592315 4.142662 0.793766 + H 4.846572 2.934673 1.132732 + H 5.288590 4.651249 0.928709 + H 6.744835 7.802927 4.148583 + H 5.981224 6.552987 5.178905 + H 5.079755 7.281473 3.796705 + H 6.919516 0.805474 3.112487 + H 5.280316 1.494462 3.102453 + H 6.260516 1.526321 4.600993 + H 1.323767 6.269682 9.133831 + H 1.527561 4.510683 8.993283 + H 0.793766 5.438218 7.657200 + H 3.954499 0.793766 7.327787 + H 4.907143 2.104697 6.559379 + H 3.404264 1.484302 5.781142 + H 1.511300 5.810548 3.888687 + H 2.990791 6.311146 4.747730 + H 1.399317 6.465175 5.533297 + &END COORD + &PRINT + &KINDS + BASIS_SET + POTENTIAL + &END KINDS + &END PRINT + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/SE/regtest-2/TEST_FILES b/tests/SE/regtest-2/TEST_FILES index 65d1e87082..d0a552fe37 100644 --- a/tests/SE/regtest-2/TEST_FILES +++ b/tests/SE/regtest-2/TEST_FILES @@ -31,3 +31,9 @@ b2h6_pm6.inp 3 ch2o_pm6.inp 3 hcn_pm6.inp 3 md.inp 11 +# d-orbitals analytical der +O2-ROKS-PM6.inp 3 +Pt-cis-2xpet3Cl2.inp 3 +hcl.inp 3 +hcl_ot.inp 3 +hcl_z.inp 3 diff --git a/tests/SE/regtest-2/hcl.inp b/tests/SE/regtest-2/hcl.inp new file mode 100644 index 0000000000..1b460d70e0 --- /dev/null +++ b/tests/SE/regtest-2/hcl.inp @@ -0,0 +1,66 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &QS + METHOD MNDOD + &SE + STO_NG 6 + &END + + &END QS + &SCF + SCF_GUESS MOPAC + &END SCF + &PRINT + &MO + EIGENVECTORS + EIGENVALUES + &END + &AO_MATRICES + DENSITY + CORE_HAMILTONIAN + OVERLAP + KOHN_SHAM_MATRIX + &END + &END + &END DFT + &SUBSYS + &CELL + ABC 30.0 30.0 30.0 + UNIT ANGSTROM + &END CELL + &COORD + H .00000000 .00000000 .00000000 + Cl .00000000 1.40000000 .00000000 + &END COORD + &PRINT + &KINDS + BASIS_SET + SE_PARAMETERS + &END + &END + &END SUBSYS + &PRINT + &FORCES + &END + &END +&END FORCE_EVAL +&GLOBAL + PROJECT HCl + RUN_TYPE MD +&END GLOBAL +&MOTION + &MD + ENSEMBLE NVE + STEPS 100 + TIMESTEP 0.5 + TEMPERATURE 300 + TEMP_TOL 100 + &END MD + &PRINT + &VELOCITIES OFF + &END + &TRAJECTORY OFF + &END + &END +&END diff --git a/tests/SE/regtest-2/hcl_ot.inp b/tests/SE/regtest-2/hcl_ot.inp new file mode 100644 index 0000000000..26af8f6eed --- /dev/null +++ b/tests/SE/regtest-2/hcl_ot.inp @@ -0,0 +1,68 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &QS + METHOD MNDOD + &SE + STO_NG 6 + &END + + &END QS + &SCF + SCF_GUESS MOPAC + &OT + &END + &END SCF + &PRINT + &MO + EIGENVECTORS + EIGENVALUES + &END + &AO_MATRICES + DENSITY + CORE_HAMILTONIAN + OVERLAP + KOHN_SHAM_MATRIX + &END + &END + &END DFT + &SUBSYS + &CELL + ABC 30.0 30.0 30.0 + UNIT ANGSTROM + &END CELL + &COORD + H .00000000 .00000000 .00000000 + Cl .00000000 1.40000000 .00000000 + &END COORD + &PRINT + &KINDS + BASIS_SET + SE_PARAMETERS + &END + &END + &END SUBSYS + &PRINT + &FORCES + &END + &END +&END FORCE_EVAL +&GLOBAL + PROJECT HCl + RUN_TYPE MD +&END GLOBAL +&MOTION + &MD + ENSEMBLE NVE + STEPS 100 + TIMESTEP 0.5 + TEMPERATURE 300 + TEMP_TOL 100 + &END MD + &PRINT + &VELOCITIES OFF + &END + &TRAJECTORY OFF + &END + &END +&END diff --git a/tests/SE/regtest-2/hcl_z.inp b/tests/SE/regtest-2/hcl_z.inp new file mode 100644 index 0000000000..4a901da2b6 --- /dev/null +++ b/tests/SE/regtest-2/hcl_z.inp @@ -0,0 +1,66 @@ +&FORCE_EVAL + METHOD Quickstep + &DFT + &QS + METHOD MNDOD + &SE + STO_NG 6 + &END + + &END QS + &SCF + SCF_GUESS MOPAC + &END SCF + &PRINT + &MO + EIGENVECTORS + EIGENVALUES + &END + &AO_MATRICES + DENSITY + CORE_HAMILTONIAN + OVERLAP + KOHN_SHAM_MATRIX + &END + &END + &END DFT + &SUBSYS + &CELL + ABC 30.0 30.0 30.0 + UNIT ANGSTROM + &END CELL + &COORD + H .00000000 .00000000 .00000000 + Cl .00000000 .00000000 1.40000000 + &END COORD + &PRINT + &KINDS + BASIS_SET + SE_PARAMETERS + &END + &END + &END SUBSYS + &PRINT + &FORCES + &END + &END +&END FORCE_EVAL +&GLOBAL + PROJECT HCl + RUN_TYPE MD +&END GLOBAL +&MOTION + &MD + ENSEMBLE NVE + STEPS 100 + TIMESTEP 0.5 + TEMPERATURE 300 + TEMP_TOL 100 + &END MD + &PRINT + &VELOCITIES OFF + &END + &TRAJECTORY OFF + &END + &END +&END