comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Analytical derivatives for d-orbitals! added few more regtests together with

Browse source 
a "real" chemistry example.


svn-origin-rev: 7331
This commit is contained in:
Teodoro Laino 2008-05-01 18:52:48 +00:00
parent 6435094b8e
commit 31b9a200d4
10 changed files with 758 additions and 33 deletions
Download patch file Download diff file Expand all files Collapse all files

194
src/semi_empirical_int_ana.F
View file

@ -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
! <S S | D S>
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
! <S S | D P >
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
! <S S | D D >
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
! <S S | D+P+>
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
! <S S | D+D+>
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
! <S S | D#D#>
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
! <D S | S S>
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
! <D P | S S >
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
! <D D | S S >
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
! <D+P+| S S >
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
! <D+D+| S S >
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
! <D#D#| S S >
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

6
src/semi_empirical_int_debug.F
View file

@ -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

14
src/semi_empirical_int_num.F
View file

@ -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

207
src/semi_empirical_int_utils.F
View file

@ -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 <ij|kl> involving
!> d-orbitals. The choice of the linear quadrupole was REALLY unhappy
!> \brief General driver for computing semi-empirical integrals <ij|kl>
!> 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 <ij|kl>
!> 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 <ij|kl>
!> 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)

61
tests/SE/regtest-2/O2-ROKS-PM6.inp Normal file
View file

@ -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

103
tests/SE/regtest-2/Pt-cis-2xpet3Cl2.inp Normal file
View file

@ -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

6
tests/SE/regtest-2/TEST_FILES
View file

@ -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

66
tests/SE/regtest-2/hcl.inp Normal file
View file

@ -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

68
tests/SE/regtest-2/hcl_ot.inp Normal file
View file

@ -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

66
tests/SE/regtest-2/hcl_z.inp Normal file
View file

@ -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