diff --git a/man/xtb.1.adoc b/man/xtb.1.adoc
index 4462aaf38..63c53c4fd 100644
--- a/man/xtb.1.adoc
+++ b/man/xtb.1.adoc
@@ -242,7 +242,7 @@ NOTE: You can only select *one* runtyp, only the first runtyp will be used
*--hess*::
perform a numerical hessian calculation on input geometry
-
+
*--ohess* ['LEVEL']::
perform a numerical hessian calculation on an `ancopt(3)` optimized
geometry
@@ -251,6 +251,9 @@ NOTE: You can only select *one* runtyp, only the first runtyp will be used
perform a biased numerical hessian calculation on an `ancopt(3)` optimized
geometry
+ Hessian runs automatically perform a compliance analysis for nonperiodic
+ systems with more than one atom and no frozen atoms.
+
*--o1nh*::
perform the numerical hessian calculation using the ODLR approximation (O1NumHess)
@@ -414,6 +417,11 @@ OUTPUT
*hessian*::
contains the (not mass weighted) cartesian Hessian, *--namespace* compatible.
+*compliance.dat*::
+ compliance analysis from an eligible Hessian run. Lists every internal
+ coordinate with C~ii~ and 1/C~ii~, followed by up to 20 strongest
+ off-diagonal couplings per coordinate in descending order of magnitude.
+
*xtbtopo.mol*::
topology information written in molfile format.
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 92c1b4857..d1069e719 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -23,6 +23,7 @@ add_subdirectory("docking")
add_subdirectory("extern")
add_subdirectory("freq")
add_subdirectory("gfnff")
+add_subdirectory("internals")
add_subdirectory("io")
add_subdirectory("iff")
add_subdirectory("main")
@@ -52,6 +53,7 @@ list(APPEND srcs
"${dir}/broyden.f90"
"${dir}/charge_model.f90"
"${dir}/coffee.f90"
+ "${dir}/compliance.f90"
"${dir}/constrain_param.f90"
"${dir}/constrain_pot.f90"
"${dir}/constr.f90"
diff --git a/src/bmatrix.f90 b/src/bmatrix.f90
index 52d9485c3..5af4cafaf 100644
--- a/src/bmatrix.f90
+++ b/src/bmatrix.f90
@@ -14,23 +14,27 @@
!
! You should have received a copy of the GNU Lesser General Public License
! along with xtb. If not, see .
+
+!> Wilson B rows for internal coordinates and Cartesian Hessian accumulation
+!> helpers shared by internal-coordinate sets and the model Hessians.
module xtb_bmatrix
use xtb_mctc_accuracy, only : wp
use xtb_mctc_constants, only : pi
use xtb_mctc_math, only : crossProd
+ use xtb_internals_type, only : internal_coords_set_type, coord_bond, &
+ & coord_angle, coord_dihedral, coord_linbend
implicit none
private
- public :: bmat_bond, bmat_angle, bmat_linbend, &
+ public :: bmat_bond, bmat_angle, bmat_linbend, linbend_frame, &
& bmat_torsion, bmat_outofplane, oop_angle
public :: bmat_accum_packed, bmat_accum_dense, bmat_accum_pairblock_packed
+ public :: get_bmatrix
- character(len=*), parameter :: source = 'xtb_bmatrix'
+ character(len=*), parameter :: source = "xtb_bmatrix"
!> Tolerance for near-linear/collinear detection (radians)
real(wp), parameter :: tol_collinear = 1.0e-13_wp
- !> Tolerance for near-zero vector norms
- real(wp), parameter :: tol_norm = 1.0e-15_wp
!> Tolerance for cross-product magnitude (collinearity guard)
real(wp), parameter :: tol_cross = 1.0e-10_wp
!> Tolerance for near-linear bend detection
@@ -40,15 +44,20 @@ module xtb_bmatrix
contains
-!> Clamp a cosine to [-1, 1], rounding can push it outside and make acos NaN
+!> Clamp a cosine to [-1, 1]; rounding can push it outside and make acos NaN.
elemental function clamp(co) result(c)
+ !> Cosine to clamp.
real(wp), intent(in) :: co
+
real(wp) :: c
c = min(1.0_wp, max(-1.0_wp, co))
end function clamp
-!> Wilson B matrix for bonds
+!> Wilson B row for a bond length (2 atoms).
+!>
+!> Row layout is [atom1(1:3), atom2(4:6)] for vec = xyz(atom1) - xyz(atom2).
pure function bmat_bond(vec) result(bmat)
+ !> Bond vector, xyz(atom1) - xyz(atom2).
real(wp), intent(in) :: vec(3)
real(wp) :: l, bmat(6)
@@ -60,14 +69,19 @@ pure function bmat_bond(vec) result(bmat)
bmat(4:6) = -vec(:) / l
end function bmat_bond
-!> Wilson B matrix for non-linear angles
+!> Wilson B row for a non-linear bond angle (3 atoms).
+!>
+!> Row layout is [outer1(1:3), centre(4:6), outer2(7:9)].
pure function bmat_angle(vec1, vec2) result(bmat)
- real(wp), intent(in) :: vec1(3), vec2(3)
+ !> Vector from the angle centre to the first outer atom.
+ real(wp), intent(in) :: vec1(3)
+ !> Vector from the angle centre to the second outer atom.
+ real(wp), intent(in) :: vec2(3)
real(wp) :: bmat(9)
real(wp) :: l1, l2, nvec1(3), nvec2(3)
real(wp) :: dl(2, 6), dnvec(2, 3, 6), dinprod(9)
- real(wp) :: dot_n1n2
+ real(wp) :: sinphi
integer :: ii
l1 = norm2(vec1)
@@ -101,51 +115,61 @@ pure function bmat_angle(vec1, vec2) result(bmat)
dinprod(ii+6) = dot_product(dnvec(2, :, ii), nvec1)
end do
- dot_n1n2 = dot_product(nvec1, nvec2)
- bmat = -dinprod / sqrt(max(eps_sqrt, 1.0_wp-dot_n1n2**2))
+ sinphi = norm2(crossProd(nvec1, nvec2))
+ bmat = -dinprod / max(sqrt(eps_sqrt), sinphi)
end function bmat_angle
-!> Wilson B matrix for linear angles
-pure function bmat_linbend(vec1, vec2) result(bmat)
- real(wp), intent(in) :: vec1(3), vec2(3)
+!> Orthonormal pair of directions perpendicular to a unit axis, the fixed frame
+!> a pair of Decius linear bends at a linear centre is measured along.
+pure subroutine linbend_frame(u, e1, e2)
+ !> Unit axis to build the frame around.
+ real(wp), intent(in) :: u(3)
+ !> Orthonormal pair of unit directions perpendicular to u.
+ real(wp), intent(out) :: e1(3), e2(3)
- real(wp) :: bmat(2, 9)
- real(wp) :: l1, l2, nvec1(3), nvec2(3)
- real(wp) :: vn(3), vn2(3), nvn
- real(wp), parameter :: xaxis(3) = [1.0_wp, 0.0_wp, 0.0_wp], yaxis(3) = [0.0_wp, 1.0_wp, 0.0_wp]
+ real(wp) :: ref(3)
- l1 = norm2(vec1)
- l2 = norm2(vec2)
- nvec1 = vec1 / l1
- nvec2 = vec2 / l2
-
- vn = crossProd(vec1, vec2)
- nvn = norm2(vn)
-
- if (nvn < tol_norm) then
- vn = xaxis - dot_product(xaxis, vec1) / l1**2 * vec1
- nvn = norm2(vn)
- if (nvn < tol_norm) then
- vn = yaxis - dot_product(yaxis, vec1) / l1**2 * vec1
- nvn = norm2(vn)
- end if
+ if (abs(u(1)) < 0.9_wp) then
+ ref = [1.0_wp, 0.0_wp, 0.0_wp]
+ else
+ ref = [0.0_wp, 0.0_wp, 1.0_wp]
end if
- vn = vn / nvn
+ e1 = ref - dot_product(ref, u) * u
+ e1 = e1 / sqrt(dot_product(e1, e1))
+ e2 = crossProd(u, e1)
+end subroutine linbend_frame
- vn2 = crossProd(vec1-vec2, vn)
- vn2 = vn2 / norm2(vn2)
+!> Wilson B row of a Decius linear bend along a fixed perpendicular direction.
+!>
+!> Row layout is [outer1(1:3), centre(4:6), outer2(7:9)]; a linear centre
+!> contributes two of these, along the two directions of linbend_frame.
+pure function bmat_linbend(vec1, vec2, evec) result(bmat)
+ !> Vector from the linear centre to the first outer atom.
+ real(wp), intent(in) :: vec1(3)
+ !> Vector from the linear centre to the second outer atom.
+ real(wp), intent(in) :: vec2(3)
+ !> Unit direction perpendicular to the axis the bend is measured along.
+ real(wp), intent(in) :: evec(3)
- bmat = 0.0_wp
- bmat(2, 1:3) = vn / l1
- bmat(2, 7:9) = vn / l2
- bmat(2, 4:6) = -bmat(2, 1:3) - bmat(2, 7:9)
- bmat(1, 1:3) = vn2 / l1
- bmat(1, 7:9) = vn2 / l2
- bmat(1, 4:6) = -bmat(1, 1:3) - bmat(1, 7:9)
+ real(wp) :: bmat(9)
+ real(wp) :: l1, l2, n1(3), n2(3), g1(3), g2(3)
+
+ l1 = norm2(vec1)
+ n1 = vec1 / l1
+ l2 = norm2(vec2)
+ n2 = vec2 / l2
+ g1 = (evec - dot_product(evec, n1)*n1) / l1
+ g2 = (evec - dot_product(evec, n2)*n2) / l2
+ bmat(1:3) = g1
+ bmat(7:9) = g2
+ bmat(4:6) = -g1 - g2
end function bmat_linbend
-!> Wilson B matrix for torsions (4 atoms)
+!> Wilson B rows for a dihedral angle (4 atoms).
+!>
+!> bt(:, a) is the 3-vector of angle derivatives with respect to atom a.
pure function bmat_torsion(xyz) result(bt)
+ !> Cartesian coordinates of the atoms 1-2-3-4 along the torsion.
real(wp), intent(in) :: xyz(3, 4)
real(wp) :: bt(3, 4)
@@ -155,9 +179,9 @@ pure function bmat_torsion(xyz) result(bt)
real(wp) :: fi3, sinfi3, cosfi3
integer :: ix, iy, iz
- rij1 = norm2(xyz(:, 2)-xyz(:, 1))
- rjk1 = norm2(xyz(:, 3)-xyz(:, 2))
- rkl1 = norm2(xyz(:, 4)-xyz(:, 3))
+ rij1 = norm2(xyz(:, 2) - xyz(:, 1))
+ rjk1 = norm2(xyz(:, 3) - xyz(:, 2))
+ rkl1 = norm2(xyz(:, 4) - xyz(:, 3))
brij = bmat_strtch(xyz(:, 1:2))
brjk = bmat_strtch(xyz(:, 2:3))
brkl = bmat_strtch(xyz(:, 3:4))
@@ -174,18 +198,21 @@ pure function bmat_torsion(xyz) result(bt)
if (iy > 3) iy = iy - 3
iz = iy + 1
if (iz > 3) iz = iz - 3
- bt(ix, 1) = (brij(iy, 2)*brjk(iz, 2)-brij(iz, 2)*brjk(iy, 2)) &
+ bt(ix, 1) = (brij(iy, 2)*brjk(iz, 2) - brij(iz, 2)*brjk(iy, 2)) &
& / (rij1*sinfi2**2)
- bt(ix, 4) = (brkl(iy, 1)*brjk(iz, 1)-brkl(iz, 1)*brjk(iy, 1)) &
+ bt(ix, 4) = (brkl(iy, 1)*brjk(iz, 1) - brkl(iz, 1)*brjk(iy, 1)) &
& / (rkl1*sinfi3**2)
- bt(ix, 2) = -((rjk1-rij1*cosfi2)*bt(ix, 1) &
+ bt(ix, 2) = -((rjk1 - rij1*cosfi2)*bt(ix, 1) &
& + rkl1 * cosfi3 * bt(ix, 4)) / rjk1
- bt(ix, 3) = -(bt(ix, 1)+bt(ix, 2)+bt(ix, 4))
+ bt(ix, 3) = -(bt(ix, 1) + bt(ix, 2) + bt(ix, 4))
end do
end function bmat_torsion
-!> Wilson B matrix for out-of-plane angles (4 atoms, atom 4 central)
+!> Wilson B rows for an out-of-plane angle (4 atoms, atom 4 central).
+!>
+!> bt(:, a) is the 3-vector of angle derivatives with respect to atom a.
pure function bmat_outofplane(xyz) result(bt)
+ !> Cartesian coordinates; atom 4 is the out-of-plane centre.
real(wp), intent(in) :: xyz(3, 4)
real(wp) :: bt(3, 4)
@@ -196,35 +223,35 @@ pure function bmat_outofplane(xyz) result(bt)
real(wp) :: r42(3), r43(3)
integer :: ix, iy, iz
-! 4 -> 2 (bond in plane)
+ ! 4 -> 2 (bond in plane)
r2 = xyz(:, 2) - xyz(:, 4)
q42 = norm2(r2)
e42 = r2 / q42
-! 4 -> 3 (bond in plane)
+ ! 4 -> 3 (bond in plane)
r3 = xyz(:, 3) - xyz(:, 4)
q43 = norm2(r3)
e43 = r3 / q43
-! angle between e43 and e42
+ ! angle between e43 and e42
cosfi1 = clamp(dot_product(e43, e42))
fi1 = acos(cosfi1)
-! dirty exit when 2-3-4 collinear
- if (abs(fi1-pi) < tol_collinear) then
+ ! dirty exit when 2-3-4 collinear
+ if (abs(fi1 - pi) < tol_collinear) then
bt = 0.0_wp
return
end if
-! first two centers
+ ! first two centers
c14(:, 1) = xyz(:, 1)
c14(:, 2) = xyz(:, 4)
-! 3rd center: cross product of r42 x r43
+ ! 3rd center: cross product of r42 x r43
r42 = xyz(:, 2) - xyz(:, 4)
r43 = xyz(:, 3) - xyz(:, 4)
c14(:, 3) = crossProd(r42, r43)
-! exit if 2-3-4 collinear
+ ! exit if 2-3-4 collinear
if (sum(c14(:, 3)**2) < tol_cross) then
bt = 0.0_wp
return
@@ -233,15 +260,15 @@ pure function bmat_outofplane(xyz) result(bt)
br14 = bmat_bend(c14)
-! compute the WDC matrix
+ ! compute the WDC matrix
do ix = 1, 3
- iy = mod(ix+1, 4) + (ix+1) / 4
- iz = mod(iy+1, 4) + (iy+1) / 4
+ iy = mod(ix + 1, 4) + (ix + 1) / 4
+ iz = mod(iy + 1, 4) + (iy + 1) / 4
bt(ix, 1) = -br14(ix, 1)
bt(ix, 2) = r43(iz) * br14(iy, 3) - r43(iy) * br14(iz, 3)
bt(ix, 3) = -r42(iz) * br14(iy, 3) + r42(iy) * br14(iz, 3)
- bt(ix, 4) = -(bt(ix, 1)+bt(ix, 2)+bt(ix, 3))
+ bt(ix, 4) = -(bt(ix, 1) + bt(ix, 2) + bt(ix, 3))
end do
bt = -bt
@@ -252,6 +279,7 @@ end function bmat_outofplane
!> between atom 1 and the plane (2,3,4) minus pi/2. Returns 0.0 when
!> atoms 2-3-4 are collinear (ill-defined).
pure function oop_angle(xyz) result(theta)
+ !> Cartesian coordinates; atom 4 is the out-of-plane centre.
real(wp), intent(in) :: xyz(3, 4)
real(wp) :: theta
@@ -269,7 +297,7 @@ pure function oop_angle(xyz) result(theta)
cosfi1 = clamp(dot_product(e43, e42))
fi1 = acos(cosfi1)
- if (abs(fi1-pi) < tol_collinear) then
+ if (abs(fi1 - pi) < tol_collinear) then
theta = 0.0_wp
return
end if
@@ -286,11 +314,58 @@ pure function oop_angle(xyz) result(theta)
end if
c14(:, 3) = c14(:, 3) + xyz(:, 4)
- theta = bend_angle(c14) - 0.5_wp*pi
+ theta = bend_angle(c14) - 0.5_wp * pi
end function oop_angle
-!> Stretch B row for 2-atom fragment (private helper)
+!> Assemble the Wilson B matrix for an internal-coordinate set.
+pure subroutine get_bmatrix(internals, xyz, bmat)
+ !> Internal-coordinate definitions.
+ class(internal_coords_set_type), intent(in) :: internals
+ !> Cartesian coordinates, dimension (3, number of atoms).
+ real(wp), intent(in) :: xyz(:, :)
+ !> Wilson B matrix, dimension (internals%ncoords, 3*number of atoms).
+ real(wp), intent(out) :: bmat(:, :)
+
+ integer :: ic, ia, atom, ncoord
+ real(wp) :: brow(12)
+
+ bmat = 0.0_wp
+ do ic = 1, internals%ncoords
+ select case (internals%kind(ic))
+ case (coord_bond)
+ ncoord = 2
+ brow(1:6) = bmat_bond(xyz(:, internals%atoms(1, ic)) &
+ & - xyz(:, internals%atoms(2, ic)))
+ case (coord_angle)
+ ncoord = 3
+ brow(1:9) = bmat_angle( &
+ & xyz(:, internals%atoms(1, ic)) - xyz(:, internals%atoms(2, ic)), &
+ & xyz(:, internals%atoms(3, ic)) - xyz(:, internals%atoms(2, ic)))
+ case (coord_linbend)
+ ncoord = 3
+ brow(1:9) = bmat_linbend( &
+ & xyz(:, internals%atoms(1, ic)) - xyz(:, internals%atoms(2, ic)), &
+ & xyz(:, internals%atoms(3, ic)) - xyz(:, internals%atoms(2, ic)), &
+ & internals%frame(:, ic))
+ case (coord_dihedral)
+ ncoord = 4
+ brow = reshape(bmat_torsion(xyz(:, internals%atoms(:, ic))), [12])
+ case default
+ error stop "bmatrix: unknown internal coordinate kind"
+ end select
+
+ do ia = 1, ncoord
+ atom = internals%atoms(ia, ic)
+ bmat(ic, 3*atom-2:3*atom) = brow(3*ia-2:3*ia)
+ end do
+ end do
+end subroutine get_bmatrix
+
+!> Stretch B row for a 2-atom fragment (private helper).
+!>
+!> Row layout is [atom1(1:3), atom2(4:6)].
pure function bmat_strtch(xyz) result(b)
+ !> Cartesian coordinates of atoms 1 and 2.
real(wp), intent(in) :: xyz(3, 2)
real(wp) :: b(3, 2)
@@ -302,8 +377,9 @@ pure function bmat_strtch(xyz) result(b)
b(:, 2) = -b(:, 1)
end function bmat_strtch
-!> Bend angle for 3-atom fragment (private helper)
+!> Bend angle for a 3-atom fragment in radians (private helper).
pure function bend_angle(xyz) result(fir)
+ !> Cartesian coordinates of atoms 1-2-3; atom 2 is the angle centre.
real(wp), intent(in) :: xyz(3, 3)
real(wp) :: fir
@@ -317,7 +393,7 @@ pure function bend_angle(xyz) result(fir)
crap = 0.0_wp
do i = 1, 3
co = co + brij(i, 1) * brjk(i, 2)
- crap = crap + (brjk(i, 2)+brij(i, 1))**2
+ crap = crap + (brjk(i, 2) + brij(i, 1))**2
end do
co = clamp(co)
if (sqrt(crap) < tol_bend) then
@@ -325,13 +401,16 @@ pure function bend_angle(xyz) result(fir)
else
fir = acos(co)
end if
- if (abs(fir-pi) < tol_collinear) then
+ if (abs(fir - pi) < tol_collinear) then
fir = pi
end if
end function bend_angle
-!> Bend B row for 3-atom fragment (private helper)
+!> Bend B rows for a 3-atom fragment (private helper).
+!>
+!> bf(:, a) is the 3-vector of angle derivatives with respect to atom a.
pure function bmat_bend(xyz) result(bf)
+ !> Cartesian coordinates of atoms 1-2-3; atom 2 is the angle centre.
real(wp), intent(in) :: xyz(3, 3)
real(wp) :: bf(3, 3)
@@ -342,13 +421,13 @@ pure function bmat_bend(xyz) result(bf)
brij = bmat_strtch(xyz(:, 1:2))
brjk = bmat_strtch(xyz(:, 2:3))
- rij1 = norm2(xyz(:, 2)-xyz(:, 1))
- rjk1 = norm2(xyz(:, 3)-xyz(:, 2))
+ rij1 = norm2(xyz(:, 2) - xyz(:, 1))
+ rjk1 = norm2(xyz(:, 3) - xyz(:, 2))
co = 0.0_wp
crap = 0.0_wp
do i = 1, 3
co = co + brij(i, 1) * brjk(i, 2)
- crap = crap + (brjk(i, 2)+brij(i, 1))**2
+ crap = crap + (brjk(i, 2) + brij(i, 1))**2
end do
co = clamp(co)
if (sqrt(crap) < tol_bend) then
@@ -356,37 +435,35 @@ pure function bmat_bend(xyz) result(bf)
si = sqrt(crap)
else
fir = acos(co)
- si = sqrt(1.0_wp-co**2)
+ si = sqrt(1.0_wp - co**2)
end if
bf = 0.0_wp
- if (abs(fir-pi) < tol_collinear) then
+ if (abs(fir - pi) < tol_collinear) then
return
end if
do i = 1, 3
- bf(i, 1) = (co*brij(i, 1)-brjk(i, 2)) / (si*rij1)
- bf(i, 3) = (co*brjk(i, 2)-brij(i, 1)) / (si*rjk1)
- bf(i, 2) = -(bf(i, 1)+bf(i, 3))
+ bf(i, 1) = (co*brij(i, 1) - brjk(i, 2)) / (si*rij1)
+ bf(i, 3) = (co*brjk(i, 2) - brij(i, 1)) / (si*rjk1)
+ bf(i, 2) = -(bf(i, 1) + bf(i, 3))
end do
end function bmat_bend
!> Accumulate scale * B^T B into a packed lower-triangular Cartesian Hessian.
!>
-!> Layout:
-!> nat - number of atoms in the full molecule
-!> hess - packed array of length (3*nat)*(3*nat+1)/2
-!> atoms - fragment atom indices (1-based global), size nfrag
-!> brow - flattened Wilson B row, length 3*nfrag, ordered as
-!> [a1_x, a1_y, a1_z, a2_x, a2_y, a2_z, ...]
-!> scale - scalar force constant multiplying the outer product
-!>
!> Only lower-packed entries with global Cartesian j <= i are updated.
!> Caller is responsible for ensuring atoms are unique and brow has the
!> expected length; this routine does no validation.
pure subroutine bmat_accum_packed(nat, hess, atoms, brow, scale)
+ !> Number of atoms in the full molecule.
integer, intent(in) :: nat
- real(wp), intent(inout) :: hess((3*nat)*(3*nat+1)/2)
+ !> Packed Hessian of length (3*nat)*(3*nat + 1)/2, updated in place.
+ real(wp), intent(inout) :: hess((3*nat)*(3*nat + 1)/2)
+ !> Fragment atom indices (1-based global), size nfrag.
integer, intent(in) :: atoms(:)
+ !> Flattened Wilson B row, length 3*nfrag, ordered
+ !> [a1_x, a1_y, a1_z, a2_x, a2_y, a2_z, ...].
real(wp), intent(in) :: brow(:)
+ !> Scalar force constant multiplying the outer product.
real(wp), intent(in) :: scale
integer :: nfrag, p, q, cp, cq, gi, gj, imax, imin
@@ -395,15 +472,15 @@ pure subroutine bmat_accum_packed(nat, hess, atoms, brow, scale)
nfrag = size(atoms)
do p = 1, nfrag
do cp = 1, 3
- gi = (atoms(p)-1)*3 + cp
+ gi = (atoms(p) - 1) * 3 + cp
do q = 1, nfrag
do cq = 1, 3
- gj = (atoms(q)-1)*3 + cq
+ gj = (atoms(q) - 1) * 3 + cq
if (gj > gi) cycle
val = scale * brow((p-1)*3+cp) * brow((q-1)*3+cq)
imax = gi
imin = gj
- hess(imax*(imax-1)/2 + imin) = hess(imax*(imax-1)/2 + imin) + val
+ hess(imax*(imax-1)/2+imin) = hess(imax*(imax-1)/2+imin) + val
end do
end do
end do
@@ -412,21 +489,19 @@ end subroutine bmat_accum_packed
!> Accumulate scale * B^T B into a dense Cartesian Hessian.
!>
-!> Layout:
-!> hess - dense (3*nat) x (3*nat) array
-!> atoms - fragment atom indices (1-based global), size nfrag
-!> brow - flattened Wilson B row, length 3*nfrag, ordered as
-!> [a1_x, a1_y, a1_z, a2_x, a2_y, a2_z, ...]
-!> scale - scalar force constant multiplying the outer product
-!>
!> Updates both (i,j) and (j,i) entries. Caller is responsible for
!> ensuring atoms are unique and brow has the expected length; this
!> routine does no validation. Thread-safe when caller passes a
!> thread-local Hessian.
pure subroutine bmat_accum_dense(hess, atoms, brow, scale)
+ !> Dense Cartesian Hessian, updated in place.
real(wp), intent(inout) :: hess(:, :)
+ !> Fragment atom indices (1-based global), size nfrag.
integer, intent(in) :: atoms(:)
+ !> Flattened Wilson B row, length 3*nfrag, ordered
+ !> [a1_x, a1_y, a1_z, a2_x, a2_y, a2_z, ...].
real(wp), intent(in) :: brow(:)
+ !> Scalar force constant multiplying the outer product.
real(wp), intent(in) :: scale
integer :: nfrag, p, q, cp, cq, gi, gj
@@ -435,10 +510,10 @@ pure subroutine bmat_accum_dense(hess, atoms, brow, scale)
nfrag = size(atoms)
do p = 1, nfrag
do cp = 1, 3
- gi = (atoms(p)-1)*3 + cp
+ gi = (atoms(p) - 1) * 3 + cp
do q = 1, nfrag
do cq = 1, 3
- gj = (atoms(q)-1)*3 + cq
+ gj = (atoms(q) - 1) * 3 + cq
val = scale * brow((p-1)*3+cp) * brow((q-1)*3+cq)
hess(gi, gj) = hess(gi, gj) + val
end do
@@ -451,36 +526,39 @@ end subroutine bmat_accum_dense
!> lower-triangular Hessian. Given a symmetric 3x3 block 'mat', stamps:
!> (i,i) -= mat, (j,j) -= mat, (i,j) += mat
!> matching the sign pattern of a pair second derivative where moving
-!> either atom in the same direction cancels and moving them apart adds.
pure subroutine bmat_accum_pairblock_packed(nat, hess, i, j, mat)
+ !> Number of atoms in the full molecule.
integer, intent(in) :: nat
- real(wp), intent(inout) :: hess((3*nat)*(3*nat+1)/2)
+ !> Packed Hessian of length (3*nat)*(3*nat + 1)/2, updated in place.
+ real(wp), intent(inout) :: hess((3*nat)*(3*nat + 1)/2)
+ !> Atom indices (1-based global) of the pair block.
integer, intent(in) :: i, j
- real(wp), intent(in) :: mat(3,3)
+ !> Symmetric 3x3 pair second-derivative block.
+ real(wp), intent(in) :: mat(3, 3)
integer :: ic, jc, gi, gj, imax, imin
do ic = 1, 3
do jc = 1, ic
- gi = (i-1)*3 + ic
- gj = (i-1)*3 + jc
+ gi = (i - 1) * 3 + ic
+ gj = (i - 1) * 3 + jc
imax = gi
imin = gj
- hess(imax*(imax-1)/2 + imin) = hess(imax*(imax-1)/2 + imin) - mat(ic,jc)
- gi = (j-1)*3 + ic
- gj = (j-1)*3 + jc
+ hess(imax*(imax-1)/2+imin) = hess(imax*(imax-1)/2+imin) - mat(ic, jc)
+ gi = (j - 1) * 3 + ic
+ gj = (j - 1) * 3 + jc
imax = gi
imin = gj
- hess(imax*(imax-1)/2 + imin) = hess(imax*(imax-1)/2 + imin) - mat(ic,jc)
+ hess(imax*(imax-1)/2+imin) = hess(imax*(imax-1)/2+imin) - mat(ic, jc)
end do
end do
do ic = 1, 3
do jc = 1, 3
- gi = (i-1)*3 + ic
- gj = (j-1)*3 + jc
+ gi = (i - 1) * 3 + ic
+ gj = (j - 1) * 3 + jc
imax = max(gi, gj)
imin = min(gi, gj)
- hess(imax*(imax-1)/2 + imin) = hess(imax*(imax-1)/2 + imin) + mat(ic,jc)
+ hess(imax*(imax-1)/2+imin) = hess(imax*(imax-1)/2+imin) + mat(ic, jc)
end do
end do
end subroutine bmat_accum_pairblock_packed
diff --git a/src/compliance.f90 b/src/compliance.f90
new file mode 100644
index 000000000..72abf23be
--- /dev/null
+++ b/src/compliance.f90
@@ -0,0 +1,555 @@
+! This file is part of xtb.
+!
+! Copyright (C) 2026 Leopold M. Seidler
+!
+! xtb is free software: you can redistribute it and/or modify it under
+! the terms of the GNU Lesser General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! xtb is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU Lesser General Public License for more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with xtb. If not, see .
+
+!> Compliance matrix driver.
+!>
+!> Coordinates: the redundant internal coordinates of the molecular graph,
+!> built from the reference geometry. Bonds come from a covalent-radius
+!> neighbour list (TNeighbourList%generate_covalent); the graph and
+!> the coordinate set are built by xtb_internals_graph and
+!> xtb_internals_redundant. The set holds bond stretches, ordinary valence
+!> angles, fixed-frame pairs for bends near linearity, and defined dihedrals.
+!> Dihedrals with a near-linear interior angle are omitted because their
+!> torsion is undefined. xtb_bmatrix::get_bmatrix supplies the Wilson B matrix.
+!>
+!> C = B H^+ B^T, H projected out of translations and rotations first;
+!> handles redundant coordinate sets (symmetric tops etc.).
+!>
+!> Equivalent to the projected-force-constant route
+!> F = G^+ Bp Hp Bp^T G^+ , G = Bp Bp^T , C = F^+ ,
+!> where Bp and Hp denote non-rigid projections. Neither G nor F nor their
+!> pseudoinverses are needed: H^+ is formed directly in the non-rigid subspace,
+!> giving the same projection of each B row. Fixed-frame linear bends need not
+!> themselves annihilate rotations.
+!>
+!> Raw numerical Hessians (the ones handed over by the frequency code) carry
+!> residual curvature along those directions, and its reciprocal would
+!> otherwise dominate C.
+!>
+!> Works for any nint and any coordinate set (diatomic, linear, mixed,
+!> general, redundant). nint is passed explicitly -- no hardcoded 3N-6.
+!>
+!> Output: bonds -> angular rows -> dihedrals, with C_ii and 1/C_ii, plus local
+!> mode frequencies for bonds. Each near-linear bend has two consecutive
+!> fixed-frame components. A nonlinear near-linear molecule still has 3N-6
+!> physical modes, and either projected component can be redundant; an exactly
+!> linear molecule has 3N-5 modes. An exactly zero C_ii is reported with
+!> 1/C_ii = +Infinity. compliance.dat reports each diagonal element and the
+!> top-20 off-diagonal couplings per coordinate sorted by |C_ij| descending.
+!>
+!> Units: C in Bohr^2/Hartree, 1/C in Eh/a0^2, 1 Eh/a0^2 = 15.570 N/cm.
+!>
+!> Ref.: K. Brandhorst, J. Grunenberg, Chem. Soc. Rev. 37 (2008), 1558.
+!> J. Grunenberg, Chem. Sci. 6 (2015), 4086.
+!>
+!> SG (with Claude), 05/26
+
+module xtb_compliance
+ use, intrinsic :: ieee_arithmetic, only : ieee_positive_inf, ieee_value
+ use xtb_mctc_accuracy, only : wp
+ use xtb_mctc_math, only : crossProd
+ use xtb_type_environment, only : TEnvironment
+ use xtb_mctc_convert, only : autoamu
+ use xtb_type_molecule, only : TMolecule
+ use xtb_type_neighbourlist, only : TNeighbourList, init
+ use xtb_internals_graph, only : graph_type, init
+ use xtb_internals_type, only : internal_coords_set_type, coord_bond, &
+ & coord_angle, coord_linbend, coord_dihedral
+ use xtb_internals_redundant, only : redundant_type, init
+ use xtb_bmatrix, only : get_bmatrix
+ use xtb_mctc_blas, only : mctc_gemm
+ use xtb_mctc_lapack, only : mctc_syev
+ implicit none
+ private
+
+ public :: compliance_driver, compute_compliance
+
+contains
+
+!> Print compliance constants (bonds, angular rows, dihedrals) for the reference
+!> geometry and write selected compliance couplings to compliance.dat.
+!> Near-linear bends use two consecutive components in a fixed reference frame.
+subroutine compliance_driver(env, mol, hess)
+ !> Calculation environment carrying output unit and error state.
+ type(TEnvironment), intent(inout) :: env
+ !> Molecular structure containing atom count, numbers, geometry, and masses.
+ type(TMolecule), intent(in) :: mol
+ !> Cartesian Hessian in Hartree/Bohr^2, dimension (3*mol%n, 3*mol%n).
+ real(wp), intent(in) :: hess(3*mol%n, 3*mol%n)
+
+ type(TNeighbourList) :: neigh_list
+ type(graph_type) :: graph
+ type(redundant_type) :: internals
+ real(wp), allocatable :: bmat(:, :), compl(:, :)
+ logical :: failed
+
+ write(env%unit, "(a)") ""
+ call generic_header(env%unit, "Compliance Constants", 49, 10)
+ write(env%unit, "(a)") ""
+ write(env%unit, "(a)") "References for compliance constants and relaxed force constants:"
+ write(env%unit, "(a)") &
+ " K. Brandhorst, J. Grunenberg, Chem. Soc. Rev. 37 (2008), 1558."
+ write(env%unit, "(a)") &
+ " J. Grunenberg, Chem. Sci. 6 (2015), 4086."
+
+ call init(neigh_list, mol%n)
+ call neigh_list%generate_covalent(mol%at, mol%xyz)
+ call init(graph, neigh_list)
+ call init(internals, graph, mol%xyz)
+
+ allocate(bmat(internals%ncoords, 3*mol%n), compl(internals%ncoords, internals%ncoords))
+ call get_bmatrix(internals, mol%xyz, bmat)
+ call compute_compliance(env, hess, bmat, mol%xyz, mol%n, internals%ncoords, compl)
+ call env%check(failed)
+ if (failed) return
+ call print_compl(env%unit, mol, internals, compl)
+
+end subroutine compliance_driver
+
+
+!> Print diagonal elements of the compliance matrix in the order of bonds,
+!> angular rows, and dihedrals. Fixed-frame linear-bend pairs are labelled
+!> separately; either component can be redundant after rigid-mode projection.
+!>
+!> Also calls write_compliance_dat to write selected couplings to a file.
+subroutine print_compl(unit, mol, internals, C)
+ !> Molecular structure containing atom count, numbers, and masses.
+ type(TMolecule), intent(in) :: mol
+ !> Redundant coordinate set with the definitions and the reference values.
+ type(redundant_type), intent(in) :: internals
+ !> Formatted output unit.
+ integer, intent(in) :: unit
+ !> Compliance matrix, dimension (ncoords, ncoords), in Bohr^2/Hartree.
+ real(wp), intent(in) :: C(internals%ncoords, internals%ncoords)
+
+ integer :: ic, k, a1, a2, coord_width, atom_width, table_width
+ integer :: idx_ord(internals%ncoords)
+ real(wp) :: cc, invcc, mu, freq
+ character(20) :: s
+ character(32) :: ibuf
+ character(128) :: fmt
+ ! local mode frequency: nu_a = fac * sqrt(k_a[a.u.] / mu[a.u.])
+ ! k_a = 1/C_ii in Eh/a0^2
+ ! mu = m_A*m_B/(m_A+m_B) in electron masses
+ ! fac = 1/(2*pi*c) * sqrt(Eh/(a0^2 * amu))
+ ! = 5140.4869 cm^-1
+ ! Convert mu to amu with autoamu before applying fac.
+ real(wp), parameter :: fac = 5140.4869_wp ! cm^-1, mu converted to amu
+
+ write(unit, "(a)") ""
+ write(unit, "(a)") "Table columns C and 1/C report diagonal compliance elements C_ii and their reciprocals."
+ write(unit, "(a)") "units: C in Bohr^2/Hartree (a0^2/Eh)"
+ write(unit, "(a)") " 1/C in Eh/a0^2 (relaxed force constant)"
+ write(unit, "(a)") " conversion: 1 Eh/a0^2 = 15.570 N/cm"
+ write(unit, "(a)") "coord value: bond stretches in Bohr (a0); angles and dihedrals in radians"
+ write(unit, "(a)") " linear bends are dimensionless fixed-frame projections"
+ write(unit, "(a)") "bond local-mode frequency nu_loc in cm^-1"
+ write(unit, "(a)") "exact zero C is reported as 1/C = +Infinity (zero response)"
+ write(unit, "(a)") ""
+ write(unit, "(a)") "Reference for bond local-mode frequencies:"
+ write(unit, "(a)") " Cremer, Kraka, Zou, J. Chem. Theory Comput. 8 (2012) 2864."
+ write(unit, "(a)") ""
+
+ write(ibuf, "(i0)") internals%ncoords
+ coord_width = max(4, len_trim(ibuf))
+ write(ibuf, "(i0)") mol%n
+ atom_width = max(3, len_trim(ibuf))
+ table_width = coord_width + 4*atom_width + 77
+ k = 0
+
+ ! 1) bond stretches
+ ! nu_a = fac * sqrt(1/(mu_AB * C_ii)) [cm^-1]
+ ! mu_AB is formed from electron-mass inputs and converted to amu below.
+ write(unit, "(a)") &
+ repeat(" ", coord_width + 1) // &
+ "type atoms " // &
+ repeat(" ", 4*(atom_width - 3)) // &
+ "coord value C 1/C nu_loc/cm-1"
+ write(unit, "(a)") repeat("-", table_width)
+ write(fmt, '("(i",i0,",1x,a14,2(a2,i",i0,",3x),",i0,"x,4f10.2)")') &
+ coord_width, atom_width, 2*atom_width + 10
+ do ic = 1, internals%ncoords
+ if (internals%kind(ic) /= coord_bond) cycle
+ k = k + 1
+ idx_ord(k) = ic
+ cc = C(ic, ic)
+ invcc = reciprocal(cc)
+ s = "bond stretch"
+ a1 = internals%atoms(1, ic)
+ a2 = internals%atoms(2, ic)
+ mu = mol%atmass(a1) * mol%atmass(a2) / (mol%atmass(a1) + mol%atmass(a2)) * autoamu
+ freq = fac * sqrt(invcc/mu)
+ write(unit, fmt) &
+ k, s, mol%sym(a1), a1, mol%sym(a2), a2, &
+ internals%q(ic), cc, invcc, freq
+ end do
+
+ ! 2) angular rows
+ write(fmt, '("(i",i0,",1x,a14,3(a2,i",i0,",3x),",i0,"x,3f10.4)")') &
+ coord_width, atom_width, atom_width + 5
+ do ic = 1, internals%ncoords
+ if (internals%kind(ic) /= coord_angle .and. &
+ & internals%kind(ic) /= coord_linbend) cycle
+ k = k + 1
+ idx_ord(k) = ic
+ cc = C(ic, ic)
+ invcc = reciprocal(cc)
+ if (internals%kind(ic) == coord_angle) then
+ s = "angle"
+ else if (is_second_linbend(internals, ic)) then
+ s = "linear bend 2"
+ else
+ s = "linear bend 1"
+ end if
+ write(unit, fmt) &
+ k, s, mol%sym(internals%atoms(1, ic)), internals%atoms(1, ic), &
+ mol%sym(internals%atoms(2, ic)), internals%atoms(2, ic), &
+ mol%sym(internals%atoms(3, ic)), internals%atoms(3, ic), &
+ internals%q(ic), cc, invcc
+ end do
+
+ ! 3) dihedrals
+ write(fmt, '("(i",i0,",1x,a14,4(a2,i",i0,",3x),3f10.4)")') &
+ coord_width, atom_width
+ do ic = 1, internals%ncoords
+ if (internals%kind(ic) /= coord_dihedral) cycle
+ k = k + 1
+ idx_ord(k) = ic
+ cc = C(ic, ic)
+ s = "dihedral"
+ write(unit, fmt) &
+ k, s, mol%sym(internals%atoms(1, ic)), internals%atoms(1, ic), &
+ mol%sym(internals%atoms(2, ic)), internals%atoms(2, ic), &
+ mol%sym(internals%atoms(3, ic)), internals%atoms(3, ic), &
+ mol%sym(internals%atoms(4, ic)), internals%atoms(4, ic), &
+ internals%q(ic), cc, reciprocal(cc)
+ end do
+ write(unit, "(a)") repeat("-", table_width)
+
+ call write_compliance_dat(unit, mol, internals, C, idx_ord, k)
+
+end subroutine print_compl
+
+
+!> Write selected compliance couplings to a file named "compliance.dat".
+!>
+!> For each coordinate it outputs C_ii, its reciprocal (positive infinity when
+!> C_ii is exactly zero), and the top NCOUP off-diagonal couplings |C_ij| sorted
+!> in descending order.
+subroutine write_compliance_dat(unit, mol, internals, C, idx_ord, ncoord)
+ !> Molecular structure containing atom count and atomic numbers.
+ type(TMolecule), intent(in) :: mol
+ !> Internal-coordinate definitions.
+ class(internal_coords_set_type), intent(in) :: internals
+ !> Formatted output unit.
+ integer, intent(in) :: unit
+ !> Number of internal coordinates.
+ integer, intent(in) :: ncoord
+ !> Compliance matrix, dimension (ncoords, ncoords), in Bohr^2/Hartree.
+ real(wp), intent(in) :: C(internals%ncoords, internals%ncoords)
+ !> Coordinate order used for the printed table, dimension (ncoord).
+ integer, intent(in) :: idx_ord(ncoord)
+
+ integer, parameter :: NCOUP = 20
+ integer :: iunit, i, j, p, q, ii, jj, nc_act, coord_width, atom_width
+ integer :: jsort(ncoord)
+ real(wp) :: aval(ncoord), tmp_r
+ integer :: tmp_i
+ character(4) :: coord_type(ncoord)
+ character(:), allocatable :: atoms(:)
+ character(32) :: ibuf
+ character(128) :: fmt_coord, fmt_coupling
+
+ write(ibuf, "(i0)") mol%n
+ atom_width = 11 + 4*len_trim(ibuf)
+ allocate(character(atom_width) :: atoms(ncoord))
+ call build_labels(mol, internals, ncoord, coord_type, atoms)
+ write(ibuf, "(i0)") ncoord
+ coord_width = max(4, len_trim(ibuf))
+ fmt_coord = "(a,1x,a,1x,a,a,es22.12e3,a,es22.12e3)"
+ write(fmt_coupling, '("(1x,i",i0,",2x,a4,1x,a",i0,",2x,es22.12e3)")') &
+ coord_width, atom_width
+
+ open(newunit=iunit, file="compliance.dat", status="replace")
+
+ write(iunit, "(a)") "#"
+ write(iunit, "(a)") "# compliance.dat"
+ write(iunit, "(a)") "#"
+ write(iunit, "(a)") "# units: C in Bohr^2/Hartree (a0^2/Eh)"
+ write(iunit, "(a)") "# 1/C in Eh/a0^2 (relaxed force constant)"
+ write(iunit, "(a)") "# conversion: 1 Eh/a0^2 = 15.570 N/cm"
+ write(iunit, "(a)") "# exact zero C is reported as 1/C = +Infinity (zero response)"
+ write(iunit, "(a)") "#"
+ write(iunit, "(a)") "# Ref.: K. Brandhorst, J. Grunenberg,"
+ write(iunit, "(a)") "# Chem. Soc. Rev. 37 (2008) 1558."
+ write(iunit, "(a)") "#"
+ write(iunit, "(a,i0)") "# number of internal coordinates : ", ncoord
+ write(iunit, "(a,i0)") "# top couplings shown per coord : ", NCOUP
+ write(iunit, "(a)") "#"
+
+ do i = 1, ncoord
+ ii = idx_ord(i)
+
+ write(iunit, "(a)") ""
+ write(iunit, fmt_coord) &
+ "# coord", trim(coord_type(i)), trim(atoms(i)), &
+ " C_ii = ", C(ii, ii), " 1/C_ii = ", reciprocal(C(ii, ii))
+ write(iunit, "(a)") &
+ "#" // repeat(" ", coord_width - 1) // "j type " // &
+ "atoms" // repeat(" ", atom_width - 3) // "C_ij"
+
+ ! collect off-diagonal |C_ij|
+ nc_act = 0
+ do j = 1, ncoord
+ if (j == i) cycle
+ nc_act = nc_act + 1
+ jsort(nc_act) = j
+ aval(nc_act) = abs(C(ii, idx_ord(j)))
+ end do
+
+ ! insertion sort descending
+ do p = 2, nc_act
+ tmp_r = aval(p)
+ tmp_i = jsort(p)
+ q = p - 1
+ do while (q >= 1)
+ if (aval(q) >= tmp_r) exit
+ aval(q+1) = aval(q)
+ jsort(q+1) = jsort(q)
+ q = q - 1
+ end do
+ aval(q+1) = tmp_r
+ jsort(q+1) = tmp_i
+ end do
+
+ ! write top NCOUP couplings
+ do p = 1, min(NCOUP, nc_act)
+ j = jsort(p)
+ jj = idx_ord(j)
+ if (aval(p) < 1.0e-12_wp) exit
+ write(iunit, fmt_coupling) &
+ j, coord_type(j), atoms(j), C(ii, jj)
+ end do
+
+ end do
+
+ write(iunit, "(a)") ""
+ write(iunit, "(a)") "# end of compliance.dat"
+ close(iunit)
+
+ write(unit, *)
+ write(unit, "(a,i0,a)") &
+ "compliance matrix written to compliance.dat (", ncoord, " coordinates)"
+
+end subroutine write_compliance_dat
+
+
+!> Build human-readable labels for the internal coordinates, in the order
+!> bonds, angular rows, dihedrals. Fixed-frame bend pairs use lb1/lb2.
+subroutine build_labels(mol, internals, ncoord, coord_type, atoms)
+ !> Molecular structure containing atomic numbers.
+ type(TMolecule), intent(in) :: mol
+ !> Internal-coordinate definitions.
+ class(internal_coords_set_type), intent(in) :: internals
+ !> Number of internal coordinates.
+ integer, intent(in) :: ncoord
+ !> Coordinate type abbreviations and atom labels, dimension (ncoord).
+ character(4), intent(out) :: coord_type(ncoord)
+ character(*), intent(out) :: atoms(ncoord)
+
+ integer :: ic, k
+
+ k = 0
+ coord_type = "??"
+ atoms = "??"
+
+ do ic = 1, internals%ncoords ! bonds
+ if (internals%kind(ic) /= coord_bond) cycle
+ k = k + 1
+ coord_type(k) = "bond"
+ write(atoms(k), "(a2,i0,a,a2,i0)") &
+ mol%sym(internals%atoms(1, ic)), internals%atoms(1, ic), &
+ "-", mol%sym(internals%atoms(2, ic)), internals%atoms(2, ic)
+ end do
+ do ic = 1, internals%ncoords ! angular rows
+ if (internals%kind(ic) /= coord_angle .and. &
+ & internals%kind(ic) /= coord_linbend) cycle
+ k = k + 1
+ if (internals%kind(ic) == coord_angle) then
+ coord_type(k) = "ang"
+ else if (is_second_linbend(internals, ic)) then
+ coord_type(k) = "lb2"
+ else
+ coord_type(k) = "lb1"
+ end if
+ write(atoms(k), "(a2,i0,a,a2,i0,a,a2,i0)") &
+ mol%sym(internals%atoms(1, ic)), internals%atoms(1, ic), &
+ "-", mol%sym(internals%atoms(2, ic)), internals%atoms(2, ic), &
+ "-", mol%sym(internals%atoms(3, ic)), internals%atoms(3, ic)
+ end do
+ do ic = 1, internals%ncoords ! dihedrals
+ if (internals%kind(ic) /= coord_dihedral) cycle
+ k = k + 1
+ coord_type(k) = "dih"
+ write(atoms(k), "(a2,i0,a,a2,i0,a,a2,i0,a,a2,i0)") &
+ mol%sym(internals%atoms(1, ic)), internals%atoms(1, ic), &
+ "-", mol%sym(internals%atoms(2, ic)), internals%atoms(2, ic), &
+ "-", mol%sym(internals%atoms(3, ic)), internals%atoms(3, ic), &
+ "-", mol%sym(internals%atoms(4, ic)), internals%atoms(4, ic)
+ end do
+
+end subroutine build_labels
+
+!> True for the second row of a consecutive fixed-frame linear-bend pair.
+pure logical function is_second_linbend(internals, ic)
+ class(internal_coords_set_type), intent(in) :: internals
+ integer, intent(in) :: ic
+
+ is_second_linbend = ic > 1
+ if (is_second_linbend) then
+ is_second_linbend = internals%kind(ic - 1) == coord_linbend .and. &
+ & all(internals%atoms(1:3, ic - 1) == internals%atoms(1:3, ic))
+ end if
+
+end function is_second_linbend
+
+
+!> Reciprocal with an explicit positive-infinity result for exact zero.
+pure real(wp) function reciprocal(value)
+ real(wp), intent(in) :: value
+
+ if (value == 0.0_wp) then
+ reciprocal = ieee_value(value, ieee_positive_inf)
+ else
+ reciprocal = 1.0_wp / value
+ end if
+
+end function reciprocal
+
+
+!> Compliance constants of the reference geometry as the pseudoinverse of the
+!> Hessian projected out of translations and rotations.
+!>
+!> C = B H^+ B^T, evaluated as Z (Z D)^T with Z = B V and D = diag(1/w_i) from
+!> the eigen decomposition Hp = V W V^T of the projected, symmetrised Hessian.
+!> Fixed-frame near-linear bend pairs may contain rigid-rotation components;
+!> projection through H^+ removes those components, so one paired row can have
+!> zero response without implying a missing physical bend.
+subroutine compute_compliance(env, H, B, xyz, nat, nint, C)
+ !> Calculation environment carrying output unit and error state.
+ type(TEnvironment), intent(inout) :: env
+ !> Number of atoms.
+ integer, intent(in) :: nat
+ !> Number of internal coordinates, passed explicitly -- no hardcoded 3N-6.
+ integer, intent(in) :: nint
+ !> Cartesian Hessian in Hartree/Bohr^2, dimension (3*natoms, 3*natoms).
+ real(wp), intent(in) :: H(3*nat, 3*nat)
+ !> Wilson B matrix, dimension (nint, 3*natoms).
+ real(wp), intent(in) :: B(nint, 3*nat)
+ !> Cartesian reference coordinates in Bohr, dimension (3, natoms).
+ real(wp), intent(in) :: xyz(3, nat)
+ !> Compliance matrix, dimension (nint, nint), in Bohr^2/Hartree.
+ real(wp), intent(out) :: C(nint, nint)
+
+ integer :: i, j, ndim, nrigid, nvib, rank_h
+ real(wp) :: tol_h, normq, center(3), rotation_scale
+ ! Relative to the largest rotation norm, sqrt(epsilon) makes the rank test size-independent.
+ real(wp), parameter :: rigid_basis_tol = sqrt(epsilon(1.0_wp))
+ real(wp), allocatable :: Hp(:, :), Q(:, :), W(:), Z(:, :), ZD(:, :), &
+ & T1(:, :), T2(:, :)
+ logical :: failed
+
+ if (nint == 0) then
+ C = 0.0_wp
+ return
+ end if
+ ndim = 3 * nat
+ allocate(Hp(ndim, ndim), Q(ndim, 6), W(ndim), Z(nint, ndim), &
+ & ZD(nint, ndim), T1(ndim, 6), T2(6, ndim))
+
+ ! orthonormal basis of the rigid (translation + rotation) space
+ center = sum(xyz, dim=2) / nat
+ Q = 0.0_wp
+ do i = 1, 3
+ Q(i::3, i) = 1.0_wp
+ end do
+ do j = 1, nat
+ Q(3*j-2:3*j, 4) = crossProd([1.0_wp, 0.0_wp, 0.0_wp], xyz(:, j) - center)
+ Q(3*j-2:3*j, 5) = crossProd([0.0_wp, 1.0_wp, 0.0_wp], xyz(:, j) - center)
+ Q(3*j-2:3*j, 6) = crossProd([0.0_wp, 0.0_wp, 1.0_wp], xyz(:, j) - center)
+ end do
+ rotation_scale = max(norm2(Q(:, 4)), norm2(Q(:, 5)), norm2(Q(:, 6)))
+
+ ! modified Gram-Schmidt; the rotation about the molecular axis of a
+ ! linear molecule produces no displacement
+ nrigid = 0
+ do i = 1, 6
+ do j = 1, i - 1
+ normq = dot_product(Q(:, i), Q(:, j))
+ Q(:, i) = Q(:, i) - normq * Q(:, j)
+ end do
+ normq = norm2(Q(:, i))
+ if (i > 3 .and. normq <= rigid_basis_tol*rotation_scale) then
+ Q(:, i) = 0.0_wp
+ cycle
+ end if
+ nrigid = nrigid + 1
+ Q(:, i) = Q(:, i) / normq
+ end do
+ nvib = ndim - nrigid
+
+ ! H averaged into symmetry first: the numerical Hessian handed over by
+ ! the frequency code is only symmetric to within its finite-difference
+ ! noise
+ Hp = 0.5_wp * (H + transpose(H))
+
+ ! Hp = (1 - Q Q^T) Hp (1 - Q Q^T)
+ call mctc_gemm(Hp, Q, T1, alpha=1.0_wp, beta=0.0_wp)
+ call mctc_gemm(T1, Q, Hp, transb="T", alpha=-1.0_wp, beta=1.0_wp)
+ call mctc_gemm(Q, Hp, T2, transa="T", alpha=1.0_wp, beta=0.0_wp)
+ call mctc_gemm(Q, T2, Hp, alpha=-1.0_wp, beta=1.0_wp)
+ Hp = 0.5_wp * (Hp + transpose(Hp))
+
+ ! Hp = V W V^T, overwriting Hp with V
+ call mctc_syev(env, Hp, W, jobz="V", uplo="U")
+ call env%check(failed)
+ if (failed) return
+
+ ! Z = B V
+ ! ZD = Z D, with D(i,i) = 1/W(i) for retained modes and zero otherwise
+ ! C = Z (ZD)^T = B V D V^T B^T = B Hp^+ B^T
+ ! Scale epsilon by matrix size and largest |w| to cover eigensolver roundoff.
+ tol_h = real(ndim, wp) * epsilon(1.0_wp) * maxval(abs(W))
+ rank_h = count(abs(W) > tol_h)
+ if (rank_h /= nvib) then
+ write(env%unit, "(A,I0,A,I0)") &
+ & " Note: Hessian rank=", rank_h, " /= 3N-rigid=", nvib
+ end if
+ call mctc_gemm(B, Hp, Z, alpha=1.0_wp, beta=0.0_wp)
+ ! Apply D: divide retained modal columns by w_i and drop discarded modes.
+ ZD = 0.0_wp
+ do i = 1, ndim
+ if (abs(W(i)) > tol_h) ZD(:, i) = Z(:, i) / W(i)
+ end do
+ call mctc_gemm(Z, ZD, C, transb="T", alpha=1.0_wp, beta=0.0_wp)
+ C = 0.5_wp * (C + transpose(C))
+ deallocate(Hp, Q, W, Z, ZD, T1, T2)
+
+end subroutine compute_compliance
+
+
+end module xtb_compliance
diff --git a/src/hessian.F90 b/src/hessian.F90
index 5440ceb27..e16e7c273 100644
--- a/src/hessian.F90
+++ b/src/hessian.F90
@@ -21,6 +21,7 @@ module xtb_hessian
use xtb_freq_io, only : rdhess, wrhess, writeHessianOut, &
& write_tm_vibspectrum, g98fake, g98fake2, rddipd
use xtb_freq_project, only : trproj
+ use xtb_compliance, only : compliance_driver
implicit none
private
@@ -327,6 +328,12 @@ subroutine numhess( &
enddo
endif
+ if (mol%n > 1 .and. freezeset%n == 0 .and. .not. set%periodic) then
+ call compliance_driver(env, mol, res%hess)
+ call env%check(exitRun)
+ if (exitRun) return
+ end if
+
if (set%pr_dftbp_hessian_out) then
call writeHessianOut('hessian.out', res%hess)
write(env%unit, '(A)') "DFTB+ style hessian.out written"
diff --git a/src/internals/CMakeLists.txt b/src/internals/CMakeLists.txt
new file mode 100644
index 000000000..2c1bcc88d
--- /dev/null
+++ b/src/internals/CMakeLists.txt
@@ -0,0 +1,26 @@
+# This file is part of xtb.
+#
+# Copyright (C) 2026 Leopold M. Seidler
+#
+# xtb is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# xtb is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with xtb. If not, see .
+
+set(dir "${CMAKE_CURRENT_SOURCE_DIR}")
+
+list(APPEND srcs
+ "${dir}/graph.f90"
+ "${dir}/type.f90"
+ "${dir}/redundant.f90"
+)
+
+set(srcs ${srcs} PARENT_SCOPE)
diff --git a/src/internals/graph.f90 b/src/internals/graph.f90
new file mode 100644
index 000000000..3808d9209
--- /dev/null
+++ b/src/internals/graph.f90
@@ -0,0 +1,143 @@
+! This file is part of xtb.
+!
+! Copyright (C) 2026 Leopold M. Seidler
+!
+! xtb is free software: you can redistribute it and/or modify it under
+! the terms of the GNU Lesser General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! xtb is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU Lesser General Public License for more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with xtb. If not, see .
+
+!> For building a molecular graph, starting from a neighbour list.
+!>
+!> Every atom of the central cell is a node. Neighbours in periodic images
+!> are folded onto the atom of the central cell and duplicate entries are
+!> removed, so the adjacency is symmetric and every edge is stored once for
+!> both of its atoms. The BFS then assigns each atom a parent, giving a
+!> spanning tree per connected component.
+module xtb_internals_graph
+ use xtb_type_neighbourlist, only : TNeighbourList
+ implicit none
+ private
+
+ public :: graph_type, init
+
+ !> Molecular graph of a system
+ type :: graph_type
+
+ !> Number of atoms in the system
+ integer :: n = 0
+
+ !> Number of neighbours of each atom, the filled entries of i_neigh(:, iat)
+ integer, allocatable :: neighs(:)
+
+ !> Adjacency list, i_neigh(1:neighs(iat), iat) holds the neighbours of iat
+ integer, allocatable :: i_neigh(:, :)
+
+ !> Parent of each atom in the BFS spanning tree, zero for a component root
+ integer, allocatable :: parent(:)
+
+ !> Atoms in BFS order, order(1:n), each component starting at its root
+ integer, allocatable :: order(:)
+
+ end type graph_type
+
+ !> Build the molecular graph from a neighbour list
+ interface init
+ module procedure :: new_graph
+ end interface init
+
+contains
+
+!> Build the molecular graph from a generated neighbour list.
+!>
+!> The list has to be generated (and, for periodic systems, updated) by the
+!> caller. Only the atoms of the central cell are nodes and only the graph is
+!> built here, the list is not stored.
+pure subroutine new_graph(self, neigh_list)
+ !> Graph to fill, its allocatable components are allocated here
+ type(graph_type), intent(out) :: self
+ !> Neighbour list of the system
+ type(TNeighbourList), intent(in) :: neigh_list
+
+ integer :: nat, iat, jat, j, k, m_neigh, head, tail, root
+ ! Stamp of the last atom whose list was scanned, to drop duplicate entries
+ integer, allocatable :: seen(:)
+ ! Next free entry in i_neigh(:, i_at)
+ integer, allocatable :: pos(:)
+ logical, allocatable :: visited(:)
+
+ nat = size(neigh_list%neighs)
+ self%n = nat
+ allocate(self%neighs(nat), source = 0)
+
+ if (nat == 0) then
+ allocate(self%i_neigh(1, 0), self%parent(0), self%order(0))
+ return
+ end if
+
+ ! count the neighbours of every atom
+ allocate(seen(nat), source = 0)
+ do iat = 1, nat
+ do j = 1, neigh_list%neighs(iat)
+ jat = neigh_list%image(neigh_list%iNeigh(j, iat))
+ if (jat <= iat) cycle
+ if (seen(jat) == iat) cycle
+ seen(jat) = iat
+ self%neighs(iat) = self%neighs(iat) + 1
+ self%neighs(jat) = self%neighs(jat) + 1
+ end do
+ end do
+
+ ! fill the adjacency list
+ m_neigh = max(maxval(self%neighs), 1)
+ allocate(self%i_neigh(m_neigh, nat), source = 0)
+ allocate(pos(nat), source = 0)
+ seen(:) = 0
+ do iat = 1, nat
+ do j = 1, neigh_list%neighs(iat)
+ jat = neigh_list%image(neigh_list%iNeigh(j, iat))
+ if (jat <= iat) cycle
+ if (seen(jat) == iat) cycle
+ seen(jat) = iat
+ pos(iat) = pos(iat) + 1
+ self%i_neigh(pos(iat), iat) = jat
+ pos(jat) = pos(jat) + 1
+ self%i_neigh(pos(jat), jat) = iat
+ end do
+ end do
+
+ ! run BFS spanning forest
+ allocate(self%parent(nat), self%order(nat), source = 0)
+ allocate(visited(nat), source = .false.)
+ tail = 0
+ do root = 1, nat
+ if (visited(root)) cycle
+ visited(root) = .true.
+ tail = tail + 1
+ self%order(tail) = root
+ head = tail
+ do while (head <= tail)
+ iat = self%order(head)
+ head = head + 1
+ do k = 1, self%neighs(iat)
+ jat = self%i_neigh(k, iat)
+ if (visited(jat)) cycle
+ visited(jat) = .true.
+ self%parent(jat) = iat
+ tail = tail + 1
+ self%order(tail) = jat
+ end do
+ end do
+ end do
+
+end subroutine new_graph
+
+end module xtb_internals_graph
diff --git a/src/internals/meson.build b/src/internals/meson.build
new file mode 100644
index 000000000..3cd2c8b48
--- /dev/null
+++ b/src/internals/meson.build
@@ -0,0 +1,22 @@
+# This file is part of xtb.
+#
+# Copyright (C) 2026 Leopold M. Seidler
+#
+# xtb is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# xtb is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with xtb. If not, see .
+
+srcs += files(
+ 'graph.f90',
+ 'type.f90',
+ 'redundant.f90',
+)
diff --git a/src/internals/redundant.f90 b/src/internals/redundant.f90
new file mode 100644
index 000000000..e7de4dd74
--- /dev/null
+++ b/src/internals/redundant.f90
@@ -0,0 +1,227 @@
+! This file is part of xtb.
+!
+! Copyright (C) 2026 Leopold M. Seidler
+!
+! xtb is free software: you can redistribute it and/or modify it under
+! the terms of the GNU Lesser General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! xtb is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU Lesser General Public License for more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with xtb. If not, see .
+
+!> Redundant internal coordinates of a system.
+!>
+!> The coordinate values are evaluated once, at the reference geometry handed
+!> to init, and stored in q. They are a snapshot and go stale when the atoms
+!> move; build the set again from the geometry of interest.
+module xtb_internals_redundant
+ use xtb_mctc_accuracy, only : wp
+ use xtb_mctc_constants, only : pi
+ use xtb_mctc_math, only : crossProd
+ use xtb_basic_geo, only : bangl
+ use xtb_bmatrix, only : linbend_frame
+ use xtb_internals_graph, only : graph_type
+ use xtb_internals_type, only : internal_coords_set_type, coord_bond, &
+ & coord_angle, coord_dihedral, coord_linbend
+ implicit none
+ private
+
+ public :: redundant_type, init
+
+ !> A centre is linear when its reference angle is within five degrees of pi
+ real(wp), parameter :: lin_tol = 5.0_wp*pi/180.0_wp
+
+ !> Redundant internal coordinates of a system. The coordinate definitions
+ !> are the inherited ones, the counts and the reference values are the
+ !> redundant set's own state.
+ type, extends(internal_coords_set_type) :: redundant_type
+
+ !> Number of atoms in the system
+ integer :: n = 0
+
+ !> Number of bonds
+ integer :: nbond = 0
+
+ !> Number of angular coordinate rows; a linear angle contributes two
+ integer :: nangle = 0
+
+ !> Number of dihedrals
+ integer :: ndihedral = 0
+
+ !> Coordinate values at the reference geometry, q(self%ncoords): lengths
+ !> in bohr, regular angles and dihedrals in radians, and linear-bend
+ !> projections dimensionless
+ real(wp), allocatable :: q(:)
+
+ end type redundant_type
+
+ !> Build the redundant coordinate set of a molecular graph
+ interface init
+ module procedure :: new_redundant
+ end interface init
+
+contains
+
+!> Build the coordinate definitions from the molecular graph and evaluate them
+!> at the reference geometry.
+pure subroutine new_redundant(self, graph, xyz)
+ !> Coordinate set to fill, its allocatable components are allocated here
+ type(redundant_type), intent(out) :: self
+ !> Molecular graph of the system
+ type(graph_type), intent(in) :: graph
+ !> Cartesian reference coordinates, dimension (3, graph%n), in the length
+ !> unit the coordinate values are reported in
+ real(wp), intent(in) :: xyz(:, :)
+
+ integer :: i, j, k, a, b, ia, la, ndeg, component
+ real(wp) :: angle, b1(3), b2(3), b3(3), n1(3), n2(3), m(3)
+ real(wp) :: v1(3), v2(3), bend_frame(3, 2)
+
+ self%n = graph%n
+
+ ! count the coordinates of each kind
+ self%nbond = 0
+ do i = 1, graph%n
+ do k = 1, graph%neighs(i)
+ if (graph%i_neigh(k, i) > i) self%nbond = self%nbond + 1
+ end do
+ end do
+
+ self%nangle = 0
+ self%ndihedral = 0
+ do i = 1, graph%n
+ ndeg = graph%neighs(i)
+ do a = 1, ndeg - 1
+ do b = a + 1, ndeg
+ call bangl(xyz, graph%i_neigh(a, i), i, graph%i_neigh(b, i), angle)
+ if (angle > pi - lin_tol) then
+ self%nangle = self%nangle + 2
+ else
+ self%nangle = self%nangle + 1
+ end if
+ end do
+ end do
+ do k = 1, ndeg
+ ! every edge is enumerated once, from its smaller atom
+ j = graph%i_neigh(k, i)
+ if (j < i) cycle
+ ! a dihedral needs a valid outer bond angle at both ends of the edge
+ do a = 1, ndeg
+ ia = graph%i_neigh(a, i)
+ if (ia == j) cycle
+ do b = 1, graph%neighs(j)
+ la = graph%i_neigh(b, j)
+ if (la == i .or. la == ia) cycle
+ if (.not. valid_torsion(xyz, ia, i, j, la)) cycle
+ self%ndihedral = self%ndihedral + 1
+ end do
+ end do
+ end do
+ end do
+
+ self%ncoords = self%nbond + self%nangle + self%ndihedral
+ allocate(self%kind(self%ncoords), self%atoms(4, self%ncoords), source = 0)
+ allocate(self%frame(3, self%ncoords), self%q(self%ncoords), source = 0.0_wp)
+
+ ! fill bonds, then angular rows, then dihedrals
+ ia = 0
+ do i = 1, graph%n
+ do k = 1, graph%neighs(i)
+ j = graph%i_neigh(k, i)
+ if (j < i) cycle
+ ia = ia + 1
+ self%kind(ia) = coord_bond
+ self%atoms(1, ia) = i
+ self%atoms(2, ia) = j
+ self%q(ia) = norm2(xyz(:, i) - xyz(:, j))
+ end do
+ end do
+
+ do i = 1, graph%n
+ ndeg = graph%neighs(i)
+ do a = 1, ndeg - 1
+ do b = a + 1, ndeg
+ call bangl(xyz, graph%i_neigh(a, i), i, graph%i_neigh(b, i), angle)
+ if (angle > pi - lin_tol) then
+ v1 = xyz(:, graph%i_neigh(a, i)) - xyz(:, i)
+ v1 = v1/norm2(v1)
+ v2 = xyz(:, graph%i_neigh(b, i)) - xyz(:, i)
+ v2 = v2/norm2(v2)
+ call linbend_frame(v1, bend_frame(:, 1), bend_frame(:, 2))
+ do component = 1, 2
+ ia = ia + 1
+ self%kind(ia) = coord_linbend
+ self%atoms(1, ia) = graph%i_neigh(a, i)
+ self%atoms(2, ia) = i
+ self%atoms(3, ia) = graph%i_neigh(b, i)
+ self%frame(:, ia) = bend_frame(:, component)
+ self%q(ia) = dot_product(self%frame(:, ia), v1 + v2)
+ end do
+ else
+ ia = ia + 1
+ self%kind(ia) = coord_angle
+ self%atoms(1, ia) = graph%i_neigh(a, i)
+ self%atoms(2, ia) = i
+ self%atoms(3, ia) = graph%i_neigh(b, i)
+ self%q(ia) = angle
+ end if
+ end do
+ end do
+ end do
+
+ do i = 1, graph%n
+ do k = 1, graph%neighs(i)
+ j = graph%i_neigh(k, i)
+ if (j < i) cycle
+ do a = 1, graph%neighs(i)
+ if (graph%i_neigh(a, i) == j) cycle
+ do b = 1, graph%neighs(j)
+ la = graph%i_neigh(b, j)
+ if (la == i .or. la == graph%i_neigh(a, i)) cycle
+ if (.not. valid_torsion(xyz, graph%i_neigh(a, i), i, j, la)) cycle
+ ia = ia + 1
+ self%kind(ia) = coord_dihedral
+ self%atoms(1, ia) = graph%i_neigh(a, i)
+ self%atoms(2, ia) = i
+ self%atoms(3, ia) = j
+ self%atoms(4, ia) = la
+ b1 = xyz(:, i) - xyz(:, self%atoms(1, ia))
+ b2 = xyz(:, j) - xyz(:, i)
+ b3 = xyz(:, la) - xyz(:, j)
+ n1 = crossProd(b1, b2)
+ n2 = crossProd(b2, b3)
+ m = crossProd(n1, b2/norm2(b2))
+ self%q(ia) = atan2(dot_product(m, n2), dot_product(n1, n2))
+ end do
+ end do
+ end do
+ end do
+
+ ! every coordinate of the set was written
+ if (ia /= self%ncoords) error stop "redundant: coordinate count mismatch"
+
+end subroutine new_redundant
+
+!> Whether both interior angles define a nonsingular dihedral.
+pure logical function valid_torsion(xyz, i, j, k, l)
+ !> Cartesian coordinates
+ real(wp), intent(in) :: xyz(:, :)
+ !> Candidate dihedral atoms
+ integer, intent(in) :: i, j, k, l
+
+ real(wp) :: angle1, angle2
+
+ call bangl(xyz, i, j, k, angle1)
+ call bangl(xyz, j, k, l, angle2)
+ valid_torsion = angle1 >= lin_tol .and. angle1 <= pi - lin_tol &
+ & .and. angle2 >= lin_tol .and. angle2 <= pi - lin_tol
+end function valid_torsion
+
+
+end module xtb_internals_redundant
diff --git a/src/internals/type.f90 b/src/internals/type.f90
new file mode 100644
index 000000000..687370682
--- /dev/null
+++ b/src/internals/type.f90
@@ -0,0 +1,61 @@
+! This file is part of xtb.
+!
+! Copyright (C) 2026 Leopold M. Seidler
+!
+! xtb is free software: you can redistribute it and/or modify it under
+! the terms of the GNU Lesser General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! xtb is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU Lesser General Public License for more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with xtb. If not, see .
+!
+!> Internal coordinates of a system.
+!>
+!> A coordinate set is a plain definition of which coordinates exist and
+!> which atoms they involve; the values of the coordinates are not part of
+!> it. Concrete sets extend this type and add their own state.
+module xtb_internals_type
+ use xtb_mctc_accuracy, only : wp
+ implicit none
+ private
+
+ public :: internal_coords_set_type, coord_bond, coord_angle, coord_dihedral, &
+ & coord_linbend
+
+ !> Bond coordinate kind, one coordinate per edge
+ integer, parameter :: coord_bond = 1
+ !> Ordinary-angle coordinate kind, one coordinate per nonlinear neighbour pair
+ integer, parameter :: coord_angle = 2
+ !> Dihedral coordinate kind, one coordinate per pair of outer bonds of an edge
+ integer, parameter :: coord_dihedral = 3
+ !> Linear-bend coordinate kind, two fixed-frame coordinates per linear angle
+ integer, parameter :: coord_linbend = 4
+
+ !> Definition of a set of internal coordinates
+ type :: internal_coords_set_type
+
+ !> Number of coordinates in the set
+ integer :: ncoords = 0
+
+ !> Kind of each coordinate, kind(self%ncoords), one of the coord_* kinds
+ integer, allocatable :: kind(:)
+
+ !> Atoms of each coordinate, allocated as atoms(4, self%ncoords).
+ !> atoms(:, ic) = [i, j, 0, 0] for a bond, [i, j, k, 0] for an angle
+ !> or linear bend with j the centre, and [i, j, k, l] for a dihedral,
+ !> with zero for the atoms the coordinate kind does not use
+ integer, allocatable :: atoms(:, :)
+
+ !> Fixed unit direction of each linear bend, allocated as
+ !> frame(3, self%ncoords), zero for all other coordinate kinds
+ real(wp), allocatable :: frame(:, :)
+
+ end type internal_coords_set_type
+
+end module xtb_internals_type
diff --git a/src/mctc/lapack.f90 b/src/mctc/lapack.f90
index 0a346fe64..38a4602a1 100644
--- a/src/mctc/lapack.f90
+++ b/src/mctc/lapack.f90
@@ -22,7 +22,7 @@ module xtb_mctc_lapack
& lapack_hegvx, lapack_hpgv, lapack_hpgvd, lapack_hpgvx
use xtb_mctc_lapack_gst, only : lapack_sygst, lapack_spgst, lapack_hegst, &
& lapack_hpgst
- use xtb_mctc_lapack_stdeigval, only : lapack_syev, lapack_syevd, lapack_syevx, &
+ use xtb_mctc_lapack_stdeigval, only : mctc_syev, lapack_syev, lapack_syevd, lapack_syevx, &
& lapack_syevr, lapack_spev, lapack_spevd, lapack_spevx, lapack_heev, &
& lapack_heevd, lapack_heevx, lapack_heevr, lapack_hpev, lapack_hpevd, &
& lapack_hpevx
diff --git a/src/mctc/lapack/stdeigval.f90 b/src/mctc/lapack/stdeigval.f90
index 73aaff89f..016ac4d37 100644
--- a/src/mctc/lapack/stdeigval.f90
+++ b/src/mctc/lapack/stdeigval.f90
@@ -18,15 +18,23 @@
!> LAPACK eigenproblem solvers
module xtb_mctc_lapack_stdeigval
use xtb_mctc_accuracy, only : sp, dp
+ use xtb_type_environment, only : TEnvironment
implicit none
private
+ public :: mctc_syev
public :: lapack_syev, lapack_syevd, lapack_syevx, lapack_syevr, lapack_spev
public :: lapack_spevd, lapack_spevx
public :: lapack_heev, lapack_heevd, lapack_heevx, lapack_heevr, lapack_hpev
public :: lapack_hpevd, lapack_hpevx
+ interface mctc_syev
+ module procedure :: mctc_ssyev
+ module procedure :: mctc_dsyev
+ end interface mctc_syev
+
+
interface lapack_syev
pure subroutine ssyev(jobz, uplo, n, a, lda, w, work, lwork, info)
import :: sp
@@ -599,4 +607,144 @@ end subroutine zhpevx
contains
+!> Solve a real symmetric eigenvalue problem in single precision.
+subroutine mctc_ssyev(env, amat, eval, jobz, uplo)
+ character(len=*), parameter :: source = 'mctc_lapack_syev'
+ !> Calculation environment receiving errors.
+ type(TEnvironment), intent(inout) :: env
+ !> Symmetric matrix, overwritten with eigenvectors when requested.
+ real(sp), intent(inout) :: amat(:, :)
+ !> Eigenvalues in ascending order.
+ real(sp), intent(out) :: eval(:)
+ !> Compute eigenvectors ('V') or eigenvalues only ('N').
+ character(len=1), intent(in), optional :: jobz
+ !> Triangle containing matrix data, upper ('U') or lower ('L').
+ character(len=1), intent(in), optional :: uplo
+
+ character(len=1) :: job, ula
+ character(len=80) :: message
+ integer :: info, lda, lwork, n, stat
+ real(sp) :: query(1)
+ real(sp), allocatable :: work(:)
+
+ job = 'V'
+ if (present(jobz)) job = jobz
+ ula = 'U'
+ if (present(uplo)) ula = uplo
+
+ n = size(amat, 1)
+ if (size(amat, 2) /= n) then
+ call env%error("Matrix must be square", source)
+ return
+ end if
+ if (size(eval) < n) then
+ call env%error("Eigenvalue array is too small", source)
+ return
+ end if
+ if (job /= 'V' .and. job /= 'v' .and. job /= 'N' .and. job /= 'n') then
+ call env%error("Invalid eigenvector option", source)
+ return
+ end if
+ if (ula /= 'U' .and. ula /= 'u' .and. ula /= 'L' .and. ula /= 'l') then
+ call env%error("Invalid matrix triangle", source)
+ return
+ end if
+ if (n == 0) return
+
+ lda = max(1, n)
+ call lapack_syev(job, ula, n, amat, lda, eval, query, -1, info)
+ if (info /= 0) then
+ write(message, '(a, i0)') "Workspace query failed with info = ", info
+ call env%error(message, source)
+ return
+ end if
+
+ lwork = max(1, nint(query(1)))
+ allocate(work(lwork), stat=stat)
+ if (stat /= 0) then
+ call env%error("Workspace allocation failed", source)
+ return
+ end if
+
+ call lapack_syev(job, ula, n, amat, lda, eval, work, lwork, info)
+ if (info < 0) then
+ write(message, '(a, i0)') "LAPACK received an illegal argument, info = ", info
+ call env%error(message, source)
+ else if (info > 0) then
+ write(message, '(a, i0)') "Eigenvalue solver did not converge, info = ", info
+ call env%error(message, source)
+ end if
+end subroutine mctc_ssyev
+
+
+!> Solve a real symmetric eigenvalue problem in double precision.
+subroutine mctc_dsyev(env, amat, eval, jobz, uplo)
+ character(len=*), parameter :: source = 'mctc_lapack_syev'
+ !> Calculation environment receiving errors.
+ type(TEnvironment), intent(inout) :: env
+ !> Symmetric matrix, overwritten with eigenvectors when requested.
+ real(dp), intent(inout) :: amat(:, :)
+ !> Eigenvalues in ascending order.
+ real(dp), intent(out) :: eval(:)
+ !> Compute eigenvectors ('V') or eigenvalues only ('N').
+ character(len=1), intent(in), optional :: jobz
+ !> Triangle containing matrix data, upper ('U') or lower ('L').
+ character(len=1), intent(in), optional :: uplo
+
+ character(len=1) :: job, ula
+ character(len=80) :: message
+ integer :: info, lda, lwork, n, stat
+ real(dp) :: query(1)
+ real(dp), allocatable :: work(:)
+
+ job = 'V'
+ if (present(jobz)) job = jobz
+ ula = 'U'
+ if (present(uplo)) ula = uplo
+
+ n = size(amat, 1)
+ if (size(amat, 2) /= n) then
+ call env%error("Matrix must be square", source)
+ return
+ end if
+ if (size(eval) < n) then
+ call env%error("Eigenvalue array is too small", source)
+ return
+ end if
+ if (job /= 'V' .and. job /= 'v' .and. job /= 'N' .and. job /= 'n') then
+ call env%error("Invalid eigenvector option", source)
+ return
+ end if
+ if (ula /= 'U' .and. ula /= 'u' .and. ula /= 'L' .and. ula /= 'l') then
+ call env%error("Invalid matrix triangle", source)
+ return
+ end if
+ if (n == 0) return
+
+ lda = max(1, n)
+ call lapack_syev(job, ula, n, amat, lda, eval, query, -1, info)
+ if (info /= 0) then
+ write(message, '(a, i0)') "Workspace query failed with info = ", info
+ call env%error(message, source)
+ return
+ end if
+
+ lwork = max(1, nint(query(1)))
+ allocate(work(lwork), stat=stat)
+ if (stat /= 0) then
+ call env%error("Workspace allocation failed", source)
+ return
+ end if
+
+ call lapack_syev(job, ula, n, amat, lda, eval, work, lwork, info)
+ if (info < 0) then
+ write(message, '(a, i0)') "LAPACK received an illegal argument, info = ", info
+ call env%error(message, source)
+ else if (info > 0) then
+ write(message, '(a, i0)') "Eigenvalue solver did not converge, info = ", info
+ call env%error(message, source)
+ end if
+end subroutine mctc_dsyev
+
+
end module xtb_mctc_lapack_stdeigval
diff --git a/src/meson.build b/src/meson.build
index 3021fce3d..a95da57a9 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -23,6 +23,7 @@ subdir('docking')
subdir('extern')
subdir('freq')
subdir('gfnff')
+subdir('internals')
subdir('io')
subdir('iff')
subdir('main')
@@ -51,6 +52,7 @@ srcs += files(
'broyden.f90',
'charge_model.f90',
'coffee.f90',
+ 'compliance.f90',
'constrain_param.f90',
'constrain_pot.f90',
'constr.f90',
diff --git a/src/model_hessian/gff.f90 b/src/model_hessian/gff.f90
index 1505e757d..162476d3b 100644
--- a/src/model_hessian/gff.f90
+++ b/src/model_hessian/gff.f90
@@ -20,7 +20,7 @@
module xtb_modelhessian_gff
use xtb_mctc_accuracy, only : wp
use xtb_mctc_math, only : crossProd
- use xtb_bmatrix, only : bmat_bond, bmat_angle, bmat_linbend, &
+ use xtb_bmatrix, only : bmat_bond, bmat_angle, bmat_linbend, linbend_frame, &
& bmat_torsion, bmat_accum_packed, bmat_accum_pairblock_packed
use xtb_gfnff_data, only : TGFFData
use xtb_gfnff_neighbor, only : TNeigh
@@ -32,6 +32,7 @@ module xtb_modelhessian_gff
& distance_threshold => gff_distance_threshold
use xtb_modelhessian_type, only : TModelHessian
use xtb_type_environment, only : TEnvironment
+ use xtb_type_setvar, only : modhess_setvar
implicit none(type, external)
private
@@ -41,12 +42,10 @@ module xtb_modelhessian_gff
type(TGFFData), pointer :: param => null()
type(TGFFTopology), pointer :: topo => null()
type(TNeigh), pointer :: neigh => null()
+ real(wp) :: s6
+ logical :: do_bend, do_torsion
contains
- procedure :: stretch
- procedure :: bend
- procedure :: torsion
- procedure :: outofplane
- procedure :: add_charge
+ procedure, public :: compute_packed
end type TGFFModelHessian
public :: newGFFModelHessian
@@ -54,33 +53,60 @@ module xtb_modelhessian_gff
contains
!> Create a GFN-FF model Hessian bound to calculator-owned data
-function newGFFModelHessian(param, topo, neigh) result(model_hessian)
+function newGFFModelHessian(param, topo, neigh, modh) result(model_hessian)
!> GFN-FF parameters
type(TGFFData), intent(in), target :: param
!> Molecular GFN-FF topology
type(TGFFTopology), intent(in), target :: topo
!> GFN-FF neighbor list
type(TNeigh), intent(in), target :: neigh
+ !> Model Hessian configuration
+ type(modhess_setvar), intent(in) :: modh
type(TGFFModelHessian) :: model_hessian
model_hessian%param => param
model_hessian%topo => topo
model_hessian%neigh => neigh
+ model_hessian%s6 = modh%s6
+ model_hessian%do_bend = modh%kf /= 0.0_wp
+ model_hessian%do_torsion = modh%kt /= 0.0_wp
end function newGFFModelHessian
+!> Compute a packed GFN-FF model Hessian
+subroutine compute_packed(self, env, xyz, n, hess, at)
+ !> GFN-FF model Hessian holding parameters, topology, and neighbour list
+ class(TGFFModelHessian), intent(in) :: self
+ !> Calculation environment
+ type(TEnvironment), intent(inout) :: env
+ !> Number of atoms
+ integer, intent(in) :: n
+ !> Cartesian coordinates
+ real(wp), intent(in) :: xyz(3, n)
+ !> Packed lower-triangle Hessian
+ real(wp), intent(out) :: hess((3*n)*(3*n + 1)/2)
+ !> Atomic numbers
+ integer, intent(in) :: at(n)
+
+ hess = 0.0_wp
+ call stretch(self, xyz, n, hess, at)
+ if (self%do_bend) call bend(self, xyz, n, hess, at)
+ if (self%do_torsion) call torsion(self, xyz, n, hess, at)
+end subroutine compute_packed
+
+
!> Add GFN-FF bond, Coulomb, and dispersion contributions
-subroutine stretch(self, xyz, n, hess, at, kr, kd, s6, lcutoff, rcut)
+subroutine stretch(self, xyz, n, hess, at)
+ !> GFN-FF model Hessian holding parameters, topology, and neighbour list
class(TGFFModelHessian), intent(in) :: self
+ !> Number of atoms
integer, intent(in) :: n
+ !> Cartesian coordinates
real(wp), intent(in) :: xyz(3, n)
+ !> Packed lower-triangle Hessian, updated in place
real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
+ !> Atomic numbers
integer, intent(in) :: at(n)
- real(wp), intent(in) :: kr
- real(wp), intent(in) :: kd
- real(wp), intent(in) :: s6
- logical, intent(inout) :: lcutoff(n, n)
- real(wp), intent(in) :: rcut
integer :: ibond, i, j, ir, jr
integer :: mapped_at(n)
@@ -88,7 +114,6 @@ subroutine stretch(self, xyz, n, hess, at, kr, kd, s6, lcutoff, rcut)
real(wp) :: pair_hessian(3, 3)
mapped_at = gff_atomic_number(at)
- lcutoff = .false.
do ibond = 1, self%neigh%nbond
i = self%neigh%blist(1, ibond)
@@ -106,7 +131,7 @@ subroutine stretch(self, xyz, n, hess, at, kr, kd, s6, lcutoff, rcut)
vec = xyz(:, i) - xyz(:, j)
r2 = dot_product(vec, vec)
if (r2 > 1600.0_wp) cycle
- cdisp = -s6 * sqrt(c6(mapped_at(i))*c6(mapped_at(j)))
+ cdisp = -self%s6 * sqrt(c6(mapped_at(i))*c6(mapped_at(j)))
qq = 2.0_wp * self%topo%qa(i) * self%topo%qa(j)
r0_squared = self%param%d3r0(pair_index(at(i), at(j)))
call get_pair_hessian(vec, qq, cdisp, r0_squared, pair_hessian)
@@ -116,22 +141,24 @@ subroutine stretch(self, xyz, n, hess, at, kr, kd, s6, lcutoff, rcut)
end subroutine stretch
!> Add GFN-FF angle-bending contributions
-subroutine bend(self, xyz, n, hess, at, force_constant, kd, lcutoff)
+subroutine bend(self, xyz, n, hess, at)
+ !> GFN-FF model Hessian holding parameters, topology, and neighbour list
class(TGFFModelHessian), intent(in) :: self
+ !> Number of atoms
integer, intent(in) :: n
+ !> Cartesian coordinates
real(wp), intent(in) :: xyz(3, n)
+ !> Packed lower-triangle Hessian, updated in place
real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
+ !> Atomic numbers
integer, intent(in) :: at(n)
- real(wp), intent(in) :: force_constant
- real(wp), intent(in) :: kd
- logical, intent(in) :: lcutoff(n, n)
- integer :: iangl, i, j, m, ir, jr, mr, ilinear
+ integer :: iangl, i, j, m, ir, jr, mr
integer :: mapped_at(n)
real(wp) :: vec_mi(3), vec_mj(3), vec_ij(3)
real(wp) :: rmi2, rmj2, rmi, rmj, rij
real(wp) :: gij, cross_vec(3), sinphi
- real(wp) :: bmat9(9), linear_bmat(2, 9)
+ real(wp) :: bmat9(9), evec1(3), evec2(3)
mapped_at = gff_atomic_number(at)
do iangl = 1, self%topo%nangl
@@ -162,25 +189,29 @@ subroutine bend(self, xyz, n, hess, at, force_constant, kd, lcutoff)
bmat9 = bmat_angle(vec_mi, vec_mj)
call bmat_accum_packed(n, hess, [i, m, j], bmat9, gij)
else
- linear_bmat = bmat_linbend(vec_mi, vec_mj)
- do ilinear = 1, 2
- call bmat_accum_packed(n, hess, [i, m, j], &
- & linear_bmat(ilinear, :), gij)
- end do
+ ! linear centre: the two Decius bends along a fixed frame
+ ! perpendicular to the axis, k B^T B summed over both
+ call linbend_frame(vec_mi/rmi, evec1, evec2)
+ call bmat_accum_packed(n, hess, [i, m, j], &
+ & bmat_linbend(vec_mi, vec_mj, evec1), gij)
+ call bmat_accum_packed(n, hess, [i, m, j], &
+ & bmat_linbend(vec_mi, vec_mj, evec2), gij)
end if
end do
end subroutine bend
!> Add GFN-FF torsional contributions
-subroutine torsion(self, xyz, n, hess, at, force_constant, kd, lcutoff)
+subroutine torsion(self, xyz, n, hess, at)
+ !> GFN-FF model Hessian holding parameters, topology, and neighbour list
class(TGFFModelHessian), intent(in) :: self
+ !> Number of atoms
integer, intent(in) :: n
+ !> Cartesian coordinates
real(wp), intent(in) :: xyz(3, n)
+ !> Packed lower-triangle Hessian, updated in place
real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
+ !> Atomic numbers
integer, intent(in) :: at(n)
- real(wp), intent(in) :: force_constant
- real(wp), intent(in) :: kd
- logical, intent(in) :: lcutoff(n, n)
integer :: itors, i, j, k, l, ir, jr, kr, lr
integer :: mapped_at(n)
@@ -212,31 +243,9 @@ subroutine torsion(self, xyz, n, hess, at, force_constant, kd, lcutoff)
end do
end subroutine torsion
-!> GFN-FF has no separate out-of-plane model-Hessian term
-subroutine outofplane(self, xyz, n, hess, at, force_constant, kd, lcutoff)
- class(TGFFModelHessian), intent(in) :: self
- integer, intent(in) :: n
- real(wp), intent(in) :: xyz(3, n)
- real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
- integer, intent(in) :: at(n)
- real(wp), intent(in) :: force_constant
- real(wp), intent(in) :: kd
- logical, intent(in) :: lcutoff(n, n)
-end subroutine outofplane
-
-!> GFN-FF charge contributions are included in the pair term
-subroutine add_charge(self, env, xyz, n, hess, at, kq)
- class(TGFFModelHessian), intent(in) :: self
- type(TEnvironment), intent(inout) :: env
- integer, intent(in) :: n
- real(wp), intent(in) :: xyz(3, n)
- real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
- integer, intent(in) :: at(n)
- real(wp), intent(in) :: kq
-end subroutine add_charge
-
!> Map heavy elements to the legacy GFN-FF model-Hessian C6 table
pure elemental function gff_atomic_number(at) result(mapped_at)
+ !> Atomic number
integer, intent(in) :: at
integer :: mapped_at
@@ -250,6 +259,7 @@ end function gff_atomic_number
!> Return the packed pair index for two atomic numbers
pure elemental function pair_index(i, j) result(index)
+ !> Atomic numbers of the pair
integer, intent(in) :: i, j
integer :: index
@@ -259,10 +269,15 @@ end function pair_index
!> Evaluate the damped Coulomb and dispersion Cartesian Hessian block
pure subroutine get_pair_hessian(vec, qq, cdisp, r0_squared, hessian)
+ !> Cartesian pair displacement
real(wp), intent(in) :: vec(3)
+ !> Coulomb pair factor
real(wp), intent(in) :: qq
+ !> Dispersion pair factor
real(wp), intent(in) :: cdisp
+ !> Squared D3 reference distance
real(wp), intent(in) :: r0_squared
+ !> Symmetric Cartesian pair-Hessian block
real(wp), intent(out) :: hessian(3, 3)
real(wp) :: r2, r, r3, damped_r
@@ -284,7 +299,23 @@ end subroutine get_pair_hessian
!> Evaluate a diagonal pair-Hessian element
pure subroutine getqqxx(dx, qq, cdisp, r, r2, r3, damped_r, r0_squared, d2)
- real(wp), intent(in) :: dx, qq, cdisp, r, r2, r3, damped_r, r0_squared
+ !> Cartesian separation component
+ real(wp), intent(in) :: dx
+ !> Coulomb pair factor
+ real(wp), intent(in) :: qq
+ !> Dispersion pair factor
+ real(wp), intent(in) :: cdisp
+ !> Pair distance
+ real(wp), intent(in) :: r
+ !> Squared pair distance
+ real(wp), intent(in) :: r2
+ !> Cubed pair distance
+ real(wp), intent(in) :: r3
+ !> Damped pair distance
+ real(wp), intent(in) :: damped_r
+ !> Squared D3 reference distance
+ real(wp), intent(in) :: r0_squared
+ !> Diagonal second derivative
real(wp), intent(out) :: d2
real(wp) :: damped_r2, dx2, denominator, r6, r8
@@ -301,8 +332,28 @@ pure subroutine getqqxx(dx, qq, cdisp, r, r2, r3, damped_r, r0_squared, d2)
end subroutine getqqxx
!> Evaluate a mixed pair-Hessian element
+!>
+!> `r` is kept for a shared argument list with `getqqxx` and is unused here.
pure subroutine getqqxy(dx, dy, qq, cdisp, r, r2, r3, damped_r, r0_squared, d2)
- real(wp), intent(in) :: dx, dy, qq, cdisp, r, r2, r3, damped_r, r0_squared
+ !> First Cartesian separation component
+ real(wp), intent(in) :: dx
+ !> Second Cartesian separation component
+ real(wp), intent(in) :: dy
+ !> Coulomb pair factor
+ real(wp), intent(in) :: qq
+ !> Dispersion pair factor
+ real(wp), intent(in) :: cdisp
+ !> Pair distance, unused
+ real(wp), intent(in) :: r
+ !> Squared pair distance
+ real(wp), intent(in) :: r2
+ !> Cubed pair distance
+ real(wp), intent(in) :: r3
+ !> Damped pair distance
+ real(wp), intent(in) :: damped_r
+ !> Squared D3 reference distance
+ real(wp), intent(in) :: r0_squared
+ !> Mixed second derivative
real(wp), intent(out) :: d2
real(wp) :: damped_r2, denominator, r6, r8
diff --git a/src/model_hessian/internal.f90 b/src/model_hessian/internal.f90
index 03c69d05b..1c301954f 100644
--- a/src/model_hessian/internal.f90
+++ b/src/model_hessian/internal.f90
@@ -18,8 +18,8 @@
!> Shared internal-coordinate traversal for Swart and Lindh Hessians
module xtb_modelhessian_internal
- use xtb_bmatrix, only : bmat_bond, bmat_angle, bmat_linbend, bmat_torsion, &
- & bmat_outofplane, oop_angle, bmat_accum_packed, &
+ use xtb_bmatrix, only : bmat_bond, bmat_angle, bmat_linbend, linbend_frame, &
+ & bmat_torsion, bmat_outofplane, oop_angle, bmat_accum_packed, &
& bmat_accum_pairblock_packed
use xtb_chargemodel, only : new_charge_model_2019
use xtb_mctc_accuracy, only : wp
@@ -31,28 +31,61 @@ module xtb_modelhessian_internal
use xtb_modelhessian_type, only : TModelHessian
use xtb_type_param, only : chrg_parameter
use xtb_type_environment, only : TEnvironment
+ use xtb_type_setvar, only : modhess_setvar
implicit none(type, external)
private
!> Base for model Hessians defined by redundant internal coordinates
type, public, abstract, extends(TModelHessian) :: TInternalModelHessianBase
+ private
+ !> Bond-stretching force constant
+ real(wp) :: kr
+ !> Angle-bending force constant
+ real(wp) :: kf
+ !> Torsional force constant
+ real(wp) :: kt
+ !> Out-of-plane force constant
+ real(wp) :: ko
+ !> Distance-dependent scaling factor
+ real(wp) :: kd
+ !> Charge-dependent force constant
+ real(wp) :: kq
+ !> Pair distance cutoff
+ real(wp) :: rcut
+ !> Dispersion scaling factor
+ real(wp) :: s6
contains
- procedure :: stretch
- procedure :: bend
- procedure :: torsion
- procedure :: outofplane
- procedure :: add_charge
- procedure(pair_factor_interface), deferred :: pair_factor
+ !> Initialize model-Hessian configuration
+ procedure, public :: init => init_internal_model_hessian
+ !> Compute packed model Hessian
+ procedure, public :: compute_packed
+ !> Add bond-stretching contributions
+ procedure, private :: stretch
+ !> Add angle-bending contributions
+ procedure, private :: bend
+ !> Add torsional contributions
+ procedure, private :: torsion
+ !> Add out-of-plane contributions
+ procedure, private :: outofplane
+ !> Add charge-response contribution
+ procedure, private :: add_charge
+ !> Evaluate pair-distance decay factor
+ procedure(pair_factor_interface), deferred, public :: pair_factor
end type TInternalModelHessianBase
abstract interface
+ !> Pair-distance decay factor for the internal-coordinate force constants
pure function pair_factor_interface(self, at_i, at_j, r2, dispersion_scale, &
& outofplane) result(factor)
import :: TInternalModelHessianBase, wp
implicit none(type, external)
+ !> Model Hessian implementation
class(TInternalModelHessianBase), intent(in) :: self
+ !> Atomic numbers of the pair
integer, intent(in) :: at_i, at_j
+ !> Squared pair distance and dispersion scaling factor
real(wp), intent(in) :: r2, dispersion_scale
+ !> Out-of-plane term indicator
logical, intent(in) :: outofplane
real(wp) :: factor
@@ -61,15 +94,81 @@ end function pair_factor_interface
contains
+!> Copy model-Hessian configuration into an internal-coordinate model
+subroutine init_internal_model_hessian(self, modh)
+ !> Model Hessian implementation
+ class(TInternalModelHessianBase), intent(inout) :: self
+ !> Model Hessian configuration
+ type(modhess_setvar), intent(in) :: modh
+
+ self%kr = modh%kr
+ self%kf = modh%kf
+ self%kt = modh%kt
+ self%ko = modh%ko
+ self%kd = modh%kd
+ self%kq = modh%kq
+ self%rcut = modh%rcut
+ self%s6 = modh%s6
+end subroutine init_internal_model_hessian
+
+!> Compute a packed internal-coordinate model Hessian
+subroutine compute_packed(self, env, xyz, n, hess, at)
+ !> Model Hessian implementation
+ class(TInternalModelHessianBase), intent(in) :: self
+ !> Calculation environment
+ type(TEnvironment), intent(inout) :: env
+ !> Number of atoms
+ integer, intent(in) :: n
+ !> Cartesian coordinates
+ real(wp), intent(in) :: xyz(3, n)
+ !> Packed lower-triangle Hessian
+ real(wp), intent(out) :: hess((3*n)*(3*n + 1)/2)
+ !> Atomic numbers
+ integer, intent(in) :: at(n)
+
+ real(wp) :: kd
+ logical, allocatable :: lcutoff(:, :)
+
+ hess = 0.0_wp
+ allocate(lcutoff(n, n), source=.false.)
+
+ kd = self%kd / self%kr
+ call self%stretch(xyz, n, hess, at, self%kr, kd, self%s6, lcutoff, self%rcut)
+ if (self%kf /= 0.0_wp) then
+ call self%bend(xyz, n, hess, at, self%kf, kd, lcutoff)
+ end if
+ if (self%kt /= 0.0_wp) then
+ call self%torsion(xyz, n, hess, at, self%kt, kd, lcutoff)
+ end if
+ if (self%ko /= 0.0_wp) then
+ call self%outofplane(xyz, n, hess, at, self%ko, kd, lcutoff)
+ end if
+ if (self%kq /= 0.0_wp) then
+ call self%add_charge(env, xyz, n, hess, at, self%kq)
+ end if
+end subroutine compute_packed
+
!> Add bond-stretching and D2 contributions
pure subroutine stretch(self, xyz, n, hess, at, kr, kd, s6, lcutoff, rcut)
+ !> Model Hessian implementation
class(TInternalModelHessianBase), intent(in) :: self
+ !> Number of atoms
integer, intent(in) :: n
+ !> Cartesian coordinates
real(wp), intent(in) :: xyz(3, n)
+ !> Packed Hessian updated in place
real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
+ !> Atomic numbers
integer, intent(in) :: at(n)
- real(wp), intent(in) :: kr, kd, s6
+ !> Bond-stretching force constant
+ real(wp), intent(in) :: kr
+ !> Distance-dependent scaling factor
+ real(wp), intent(in) :: kd
+ !> Dispersion scaling factor
+ real(wp), intent(in) :: s6
+ !> Pair cutoff mask updated in place
logical, intent(inout) :: lcutoff(n, n)
+ !> Distance cutoff
real(wp), intent(in) :: rcut
integer :: i, j
@@ -93,19 +192,31 @@ end subroutine stretch
!> Add angle-bending contributions
pure subroutine bend(self, xyz, n, hess, at, force_constant, kd, lcutoff)
+ !> Model Hessian implementation
class(TInternalModelHessianBase), intent(in) :: self
+ !> Number of atoms
integer, intent(in) :: n
+ !> Cartesian coordinates
real(wp), intent(in) :: xyz(3, n)
+ !> Packed Hessian updated in place
real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
+ !> Atomic numbers
integer, intent(in) :: at(n)
+ !> Internal-coordinate force constant and distance-dependent scaling factor
real(wp), intent(in) :: force_constant, kd
+ !> Pair cutoff mask
logical, intent(in) :: lcutoff(n, n)
- real(wp), parameter :: rzero = 1.0e-10_wp
- integer :: i, j, m, ii
+ ! Minimum arm and outer-pair length (Bohr).
+ real(wp), parameter :: arm_pair_length_tol = 1.0e-10_wp
+ ! Linear-bend sine threshold (dimensionless).
+ real(wp), parameter :: linear_sine_tol = 1.0e-10_wp
+ ! Same-ray cosine proximity threshold (dimensionless).
+ real(wp), parameter :: same_ray_cosine_tol = 1.0e-12_wp
+ integer :: i, j, m
real(wp) :: vec_ij(3), vec_mi(3), vec_mj(3), cross_vec(3)
- real(wp) :: rmi2, rmi, rmj2, rmj, rij2, rrij, gij, rl2, rl
- real(wp) :: sinphi, cosphi, bmat9(9), linear_bmat(2, 9)
+ real(wp) :: rmi2, rmi, rmj2, rmj, rij2, rrij, gij
+ real(wp) :: sinphi, cosphi, bmat9(9), evec1(3), evec2(3)
do m = 1, n
do i = 1, n
@@ -113,38 +224,36 @@ pure subroutine bend(self, xyz, n, hess, at, force_constant, kd, lcutoff)
vec_mi = xyz(:, i) - xyz(:, m)
rmi2 = dot_product(vec_mi, vec_mi)
rmi = sqrt(rmi2)
+ if (rmi <= arm_pair_length_tol) cycle
do j = 1, i - 1
if (j == m) cycle
if (lcutoff(j, i) .or. lcutoff(j, m)) cycle
vec_mj = xyz(:, j) - xyz(:, m)
rmj2 = dot_product(vec_mj, vec_mj)
rmj = sqrt(rmj2)
- cosphi = dot_product(vec_mi, vec_mj) / (rmi*rmj)
- if (abs(cosphi - 1.0_wp) < 1.0e-12_wp) cycle
+ if (rmj <= arm_pair_length_tol) cycle
vec_ij = xyz(:, j) - xyz(:, i)
rij2 = dot_product(vec_ij, vec_ij)
rrij = sqrt(rij2)
+ if (rrij <= arm_pair_length_tol) cycle
+ cosphi = dot_product(vec_mi, vec_mj) / (rmi*rmj)
+ if (abs(cosphi - 1.0_wp) < same_ray_cosine_tol) cycle
gij = force_constant &
* self%pair_factor(at(m), at(i), rmi2, 0.5_wp*kd, .false.) &
* self%pair_factor(at(m), at(j), rmj2, 0.5_wp*kd, .false.)
- cross_vec = crossProd(vec_mi, vec_mj)
- rl2 = dot_product(cross_vec, cross_vec)
- if (rl2 < 1.0e-14_wp) then
- rl = 0.0_wp
- else
- rl = sqrt(rl2)
- end if
- if (rmj <= rzero .or. rmi <= rzero .or. rrij <= rzero) cycle
- sinphi = rl / (rmj*rmi)
- if (sinphi > rzero) then
+ cross_vec = crossProd(vec_mi/rmi, vec_mj/rmj)
+ sinphi = norm2(cross_vec)
+ if (sinphi > linear_sine_tol) then
bmat9 = bmat_angle(vec_mi, vec_mj)
call bmat_accum_packed(n, hess, [i, m, j], bmat9, gij)
else
- linear_bmat = bmat_linbend(vec_mi, vec_mj)
- do ii = 1, 2
- call bmat_accum_packed(n, hess, [i, m, j], &
- & linear_bmat(ii, :), gij)
- end do
+ ! linear centre: the two Decius bends along a fixed frame
+ ! perpendicular to the axis, k B^T B summed over both
+ call linbend_frame(vec_mi/rmi, evec1, evec2)
+ call bmat_accum_packed(n, hess, [i, m, j], &
+ & bmat_linbend(vec_mi, vec_mj, evec1), gij)
+ call bmat_accum_packed(n, hess, [i, m, j], &
+ & bmat_linbend(vec_mi, vec_mj, evec2), gij)
end if
end do
end do
@@ -153,12 +262,19 @@ end subroutine bend
!> Add torsional contributions using one reversal-safe orientation
pure subroutine torsion(self, xyz, n, hess, at, force_constant, kd, lcutoff)
+ !> Model Hessian implementation
class(TInternalModelHessianBase), intent(in) :: self
+ !> Number of atoms
integer, intent(in) :: n
+ !> Cartesian coordinates
real(wp), intent(in) :: xyz(3, n)
+ !> Packed Hessian updated in place
real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
+ !> Atomic numbers
integer, intent(in) :: at(n)
+ !> Internal-coordinate force constant and distance-dependent scaling factor
real(wp), intent(in) :: force_constant, kd
+ !> Pair cutoff mask
logical, intent(in) :: lcutoff(n, n)
real(wp), parameter :: a35 = (35.0_wp/180.0_wp) * pi
@@ -209,12 +325,19 @@ end subroutine torsion
!> Add out-of-plane contributions
pure subroutine outofplane(self, xyz, n, hess, at, force_constant, kd, lcutoff)
+ !> Model Hessian implementation
class(TInternalModelHessianBase), intent(in) :: self
+ !> Number of atoms
integer, intent(in) :: n
+ !> Cartesian coordinates
real(wp), intent(in) :: xyz(3, n)
+ !> Packed Hessian updated in place
real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
+ !> Atomic numbers
integer, intent(in) :: at(n)
+ !> Internal-coordinate force constant and distance-dependent scaling factor
real(wp), intent(in) :: force_constant, kd
+ !> Pair cutoff mask
logical, intent(in) :: lcutoff(n, n)
integer :: i, j, k, l
@@ -264,12 +387,19 @@ end subroutine outofplane
!> Add the neutral-molecule EEQ response contribution
subroutine add_charge(self, env, xyz, n, hess, at, kq)
+ !> Model Hessian implementation
class(TInternalModelHessianBase), intent(in) :: self
+ !> Calculation environment
type(TEnvironment), intent(inout) :: env
+ !> Number of atoms
integer, intent(in) :: n
+ !> Cartesian coordinates
real(wp), intent(in) :: xyz(3, n)
+ !> Packed Hessian updated in place
real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
+ !> Atomic numbers
integer, intent(in) :: at(n)
+ !> Charge-dependent force constant
real(wp), intent(in) :: kq
type(chrg_parameter) :: chrgeq
diff --git a/src/model_hessian/lindh.f90 b/src/model_hessian/lindh.f90
index 0a98aae66..6489dc839 100644
--- a/src/model_hessian/lindh.f90
+++ b/src/model_hessian/lindh.f90
@@ -23,6 +23,7 @@ module xtb_modelhessian_lindh
use xtb_param_model_hessian, only : TLindhParameters, &
& lindh_d2_parameters, lindh_parameters
use xtb_modelhessian_shared, only : itabrow, fk_vdw
+ use xtb_type_setvar, only : modhess_setvar
implicit none(type, external)
private
@@ -37,7 +38,29 @@ module xtb_modelhessian_lindh
type, public, extends(TLindhModelHessianBase) :: TLindhD2ModelHessian
end type TLindhD2ModelHessian
+ public :: newLindhModelHessian, newLindhD2ModelHessian
+
contains
+!> Create a Lindh 2007 model Hessian
+function newLindhModelHessian(modh) result(model_hessian)
+ !> Model Hessian configuration
+ type(modhess_setvar), intent(in) :: modh
+
+ type(TLindhModelHessian) :: model_hessian
+
+ call model_hessian%init(modh)
+end function newLindhModelHessian
+
+!> Create a Lindh-D2 model Hessian
+function newLindhD2ModelHessian(modh) result(model_hessian)
+ !> Model Hessian configuration
+ type(modhess_setvar), intent(in) :: modh
+
+ type(TLindhD2ModelHessian) :: model_hessian
+
+ call model_hessian%init(modh)
+end function newLindhD2ModelHessian
+
!> Evaluate the selected Lindh pair-distance decay factor
pure function pair_factor(self, at_i, at_j, r2, dispersion_scale, outofplane) &
diff --git a/src/model_hessian/swart.f90 b/src/model_hessian/swart.f90
index a32f8b6d4..7bf17e33a 100644
--- a/src/model_hessian/swart.f90
+++ b/src/model_hessian/swart.f90
@@ -22,6 +22,7 @@ module xtb_modelhessian_swart
use xtb_mctc_param, only : covalent_radius_2009
use xtb_modelhessian_internal, only : TInternalModelHessianBase
use xtb_modelhessian_shared, only : fk_vdw
+ use xtb_type_setvar, only : modhess_setvar
implicit none(type, external)
private
@@ -30,7 +31,19 @@ module xtb_modelhessian_swart
procedure :: pair_factor
end type TSwartModelHessian
+ public :: newSwartModelHessian
+
contains
+!> Create a Swart model Hessian
+function newSwartModelHessian(modh) result(model_hessian)
+ !> Model Hessian configuration
+ type(modhess_setvar), intent(in) :: modh
+
+ type(TSwartModelHessian) :: model_hessian
+
+ call model_hessian%init(modh)
+end function newSwartModelHessian
+
!> Evaluate the Swart pair-distance decay factor
pure function pair_factor(self, at_i, at_j, r2, dispersion_scale, outofplane) &
diff --git a/src/model_hessian/type.f90 b/src/model_hessian/type.f90
index d3d408113..6c0836468 100644
--- a/src/model_hessian/type.f90
+++ b/src/model_hessian/type.f90
@@ -17,80 +17,33 @@
!> Common interface for model Hessian implementations
module xtb_modelhessian_type
+ ! TODO:
+ ! Could cut modelhessian module down further by using bmatrix and internal
+ ! coordinate infra more extensively. Coordinate traversal is different from
+ ! generic redundant internals and model specific, so needs special constructors.
+ ! Coordinates may then be assigned model-specific force constants,
+ ! and passed to xtb_bmatrix for shared Wilson B-row construction and
+ ! accumulation of k B^T B. Pairwise Cartesian terms such as dispersion and
+ ! charge contributions remain separate.
use xtb_mctc_accuracy, only : wp
use xtb_type_environment, only : TEnvironment
- use xtb_type_setvar, only : modhess_setvar
implicit none(type, external)
private
!> Abstract model Hessian implementation
type, public, abstract :: TModelHessian
contains
- procedure(model_hessian_stretch), deferred :: stretch
- procedure(model_hessian_mode), deferred :: bend
- procedure(model_hessian_mode), deferred :: torsion
- procedure(model_hessian_mode), deferred :: outofplane
- procedure(model_hessian_charge), deferred :: add_charge
!> Compute Hessian in packed lower-triangle storage
- procedure :: compute_packed => compute_packed
+ procedure(model_hessian_packed), deferred, public :: compute_packed
!> Compute dense symmetric Hessian
- procedure :: compute_dense => compute_dense
- generic :: compute => compute_packed, compute_dense
+ procedure, private :: compute_dense
+ !> Compute packed or dense Hessian
+ generic, public :: compute => compute_packed, compute_dense
end type TModelHessian
abstract interface
- !> Add bond-stretching contributions to a packed Hessian
- subroutine model_hessian_stretch(self, xyz, n, hess, at, kr, kd, s6, &
- & lcutoff, rcut)
- import :: TModelHessian, wp
- implicit none(type, external)
- !> Model Hessian implementation
- class(TModelHessian), intent(in) :: self
- !> Number of atoms
- integer, intent(in) :: n
- !> Cartesian coordinates
- real(wp), intent(in) :: xyz(3, n)
- !> Packed Hessian updated in place
- real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
- !> Atomic numbers
- integer, intent(in) :: at(n)
- !> Stretching force constant
- real(wp), intent(in) :: kr
- !> Distance-dependent scaling factor
- real(wp), intent(in) :: kd
- !> Dispersion scaling factor
- real(wp), intent(in) :: s6
- !> Pair cutoff mask updated in place
- logical, intent(inout) :: lcutoff(n, n)
- !> Distance cutoff
- real(wp), intent(in) :: rcut
- end subroutine model_hessian_stretch
-
- !> Add an internal-coordinate contribution to a packed Hessian
- subroutine model_hessian_mode(self, xyz, n, hess, at, force_constant, &
- & kd, lcutoff)
- import :: TModelHessian, wp
- implicit none(type, external)
- !> Model Hessian implementation
- class(TModelHessian), intent(in) :: self
- !> Number of atoms
- integer, intent(in) :: n
- !> Cartesian coordinates
- real(wp), intent(in) :: xyz(3, n)
- !> Packed Hessian updated in place
- real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
- !> Atomic numbers
- integer, intent(in) :: at(n)
- !> Internal-coordinate force constant
- real(wp), intent(in) :: force_constant
- !> Distance-dependent scaling factor
- real(wp), intent(in) :: kd
- !> Pair cutoff mask
- logical, intent(in) :: lcutoff(n, n)
- end subroutine model_hessian_mode
-
- !> Add charge-dependent contributions to a packed Hessian
- subroutine model_hessian_charge(self, env, xyz, n, hess, at, kq)
+ !> Compute a packed lower-triangle Hessian
+ subroutine model_hessian_packed(self, env, xyz, n, hess, at)
import :: TModelHessian, TEnvironment, wp
implicit none(type, external)
!> Model Hessian implementation
@@ -101,58 +54,17 @@ subroutine model_hessian_charge(self, env, xyz, n, hess, at, kq)
integer, intent(in) :: n
!> Cartesian coordinates
real(wp), intent(in) :: xyz(3, n)
- !> Packed Hessian updated in place
- real(wp), intent(inout) :: hess((3*n)*(3*n + 1)/2)
+ !> Packed lower-triangle Hessian
+ real(wp), intent(out) :: hess((3*n)*(3*n + 1)/2)
!> Atomic numbers
integer, intent(in) :: at(n)
- !> Charge-dependent force constant
- real(wp), intent(in) :: kq
- end subroutine model_hessian_charge
+ end subroutine model_hessian_packed
end interface
contains
-!> Compute Hessian in packed lower-triangle storage
-subroutine compute_packed(self, env, xyz, n, hess, at, modh)
- !> Model Hessian implementation
- class(TModelHessian), intent(in) :: self
- !> Calculation environment
- type(TEnvironment), intent(inout) :: env
- !> Number of atoms
- integer, intent(in) :: n
- !> Cartesian coordinates
- real(wp), intent(in) :: xyz(3, n)
- !> Packed lower-triangle Hessian
- real(wp), intent(out) :: hess((3*n)*(3*n + 1)/2)
- !> Atomic numbers
- integer, intent(in) :: at(n)
- !> Model Hessian parameters
- type(modhess_setvar), intent(in) :: modh
-
- real(wp) :: kd
- logical, allocatable :: lcutoff(:, :)
-
- hess = 0.0_wp
- allocate(lcutoff(n, n), source=.false.)
-
- kd = modh%kd / modh%kr
- call self%stretch(xyz, n, hess, at, modh%kr, kd, modh%s6, lcutoff, modh%rcut)
- if (modh%kf /= 0.0_wp) then
- call self%bend(xyz, n, hess, at, modh%kf, kd, lcutoff)
- end if
- if (modh%kt /= 0.0_wp) then
- call self%torsion(xyz, n, hess, at, modh%kt, kd, lcutoff)
- end if
- if (modh%ko /= 0.0_wp) then
- call self%outofplane(xyz, n, hess, at, modh%ko, kd, lcutoff)
- end if
- if (modh%kq /= 0.0_wp) then
- call self%add_charge(env, xyz, n, hess, at, modh%kq)
- end if
-end subroutine compute_packed
-
!> Compute dense symmetric Hessian from packed implementation
-subroutine compute_dense(self, env, xyz, n, hess, at, modh)
+subroutine compute_dense(self, env, xyz, n, hess, at)
!> Model Hessian implementation
class(TModelHessian), intent(in) :: self
!> Calculation environment
@@ -165,14 +77,12 @@ subroutine compute_dense(self, env, xyz, n, hess, at, modh)
real(wp), intent(out) :: hess(3*n, 3*n)
!> Atomic numbers
integer, intent(in) :: at(n)
- !> Model Hessian parameters
- type(modhess_setvar), intent(in) :: modh
integer :: i, j, ij
real(wp), allocatable :: packed(:)
allocate(packed((3*n)*(3*n + 1)/2))
- call self%compute_packed(env, xyz, n, packed, at, modh)
+ call self%compute_packed(env, xyz, n, packed, at)
ij = 0
do i = 1, 3 * n
diff --git a/src/o1numhess.f90 b/src/o1numhess.f90
index 5fc956e71..f2e5bec4d 100644
--- a/src/o1numhess.f90
+++ b/src/o1numhess.f90
@@ -275,7 +275,7 @@ subroutine lr_loop(env, ndispl, g, hess_out, displdir, final_err)
! Local variables
real(wp), allocatable :: resid(:, :), hcorr(:, :), tmp(:, :)
real(wp), allocatable :: gscale(:), gscaled(:, :), xscaled(:, :)
- real(wp) :: rnorm, rnorm_prev, gnorm, relres
+ real(wp) :: rnorm, rnorm_prev, gnorm
integer :: it, N, j
N = size(g, 1)
@@ -288,10 +288,9 @@ subroutine lr_loop(env, ndispl, g, hess_out, displdir, final_err)
xscaled(:, j) = displdir(:, j) * gscale(j)
end do
- allocate(hcorr(N, N), tmp(N, N))
+ allocate(hcorr(N, N), tmp(N, ndispl))
gnorm = norm2(gscaled)
rnorm_prev = huge(1.0_wp)
- final_err = huge(1.0_wp)
! Iterative Correction Loop
loop_lr: do it = 1, maxiter_LR
@@ -300,12 +299,6 @@ subroutine lr_loop(env, ndispl, g, hess_out, displdir, final_err)
resid = gscaled(:, :ndispl) - tmp
rnorm = norm2(resid)
- if (gnorm > 0.0_wp) then
- relres = rnorm / gnorm
- else
- relres = rnorm
- end if
- final_err = relres
if (rnorm < thresh_LR) then
exit loop_lr
@@ -320,6 +313,16 @@ subroutine lr_loop(env, ndispl, g, hess_out, displdir, final_err)
end do loop_lr
+ call mctc_gemm(hess_out, xscaled(:, :ndispl), tmp)
+ rnorm = norm2(gscaled(:, :ndispl) - tmp)
+ if (gnorm > 0.0_wp) then
+ final_err = rnorm / gnorm
+ else if (rnorm == 0.0_wp) then
+ final_err = 0.0_wp
+ else
+ final_err = huge(1.0_wp)
+ end if
+
if (it > maxiter_LR) then
call env%warning("LR loop failed to converge", source)
end if
diff --git a/src/optimizer.f90 b/src/optimizer.f90
index 287456513..4f9887a7d 100644
--- a/src/optimizer.f90
+++ b/src/optimizer.f90
@@ -1228,8 +1228,9 @@ subroutine modhes(env, calc, modh, natoms, xyz, chg, Hess, pr)
use xtb_type_setvar
use xtb_modelhessian_gff, only : newGFFModelHessian
use xtb_modelhessian_type, only : TModelHessian
- use xtb_modelhessian_lindh, only : TLindhModelHessian, TLindhD2ModelHessian
- use xtb_modelhessian_swart, only : TSwartModelHessian
+ use xtb_modelhessian_lindh, only : newLindhModelHessian, &
+ & newLindhD2ModelHessian
+ use xtb_modelhessian_swart, only : newSwartModelHessian
use xtb_setparam
use xtb_type_calculator
use xtb_gfnff_calculator
@@ -1276,16 +1277,16 @@ subroutine modhes(env, calc, modh, natoms, xyz, chg, Hess, pr)
return
case(p_modh_old)
if (pr) write(env%unit,'(a)') "Using Lindh-Hessian (1995)"
- allocate(TLindhD2ModelHessian :: model_hessian)
+ allocate(model_hessian, source=newLindhD2ModelHessian(modh))
case(p_modh_lindh_d2)
if (pr) write(env%unit,'(a)') "Using Lindh-Hessian"
- allocate(TLindhD2ModelHessian :: model_hessian)
+ allocate(model_hessian, source=newLindhD2ModelHessian(modh))
case(p_modh_lindh)
if (pr) write(env%unit,'(a)') "Using Lindh-Hessian (2007)"
- allocate(TLindhModelHessian :: model_hessian)
+ allocate(model_hessian, source=newLindhModelHessian(modh))
case(p_modh_swart)
if (pr) write(env%unit,'(a)') "Using Swart-Hessian"
- allocate(TSwartModelHessian :: model_hessian)
+ allocate(model_hessian, source=newSwartModelHessian(modh))
end select
type is(TGFFCalculator) ! GFN-FF case
select case(modh%model)
@@ -1295,21 +1296,21 @@ subroutine modhes(env, calc, modh, natoms, xyz, chg, Hess, pr)
case(p_modh_old, p_modh_gff)
if (pr) write(env%unit,'(a)') "Using GFN-FF Model Hessian"
allocate(model_hessian, source=newGFFModelHessian( &
- & calc%param, calc%topo, calc%neigh))
+ & calc%param, calc%topo, calc%neigh, modh))
case(p_modh_lindh_d2)
if (pr) write(env%unit,'(a)') "Using Lindh-Hessian"
- allocate(TLindhD2ModelHessian :: model_hessian)
+ allocate(model_hessian, source=newLindhD2ModelHessian(modh))
case(p_modh_lindh)
if (pr) write(env%unit,'(a)') "Using Lindh-Hessian (2007)"
- allocate(TLindhModelHessian :: model_hessian)
+ allocate(model_hessian, source=newLindhModelHessian(modh))
case(p_modh_swart)
if (pr) write(env%unit,'(a)') "Using Swart-Hessian"
- allocate(TSwartModelHessian :: model_hessian)
+ allocate(model_hessian, source=newSwartModelHessian(modh))
end select
end select
if (allocated(model_hessian)) then
- call model_hessian%compute(env, xyz, natoms, Hess, chg, modh)
+ call model_hessian%compute(env, xyz, natoms, Hess, chg)
end if
! constraints
diff --git a/src/type/calculator.f90 b/src/type/calculator.f90
index cdc69abcb..212065a12 100644
--- a/src/type/calculator.f90
+++ b/src/type/calculator.f90
@@ -29,7 +29,7 @@ module xtb_type_calculator
use xtb_type_restart, only : TRestart
use xtb_o1numhess, only : adj_list, gen_local_hessian, &
& lr_loop, gen_displdir, get_vdw_neighbor_list, find_projected_imag_modes
- use xtb_modelhessian_swart, only : TSwartModelHessian
+ use xtb_modelhessian_swart, only : TSwartModelHessian, newSwartModelHessian
use xtb_type_setvar, only : modhess_setvar
use xtb_setparam, only : set
use xtb_param_uffvdwrad, only : get_rad
@@ -316,7 +316,8 @@ subroutine hessian_odlr(self, env, mol0, chk0, step, hess, final_err, dipgrad, p
! original parameters kept from setparam.f90 for Lindh-D2 model Hessian
modh = modhess_setvar(kr=0.4000_wp, kf=0.1300_wp, kt=0.0075_wp, &
& ko=0.16_wp, kd=0.0_wp, kq=0.0_wp, rcut=70.0_wp, s6=20.0_wp)
- call model_hessian%compute(env, mol%xyz, mol%n, h0, mol%at, modh)
+ model_hessian = newSwartModelHessian(modh)
+ call model_hessian%compute(env, mol%xyz, mol%n, h0, mol%at)
end block
call env%check(terminate_run)
if (terminate_run) then
diff --git a/src/type/neighbourlist.f90 b/src/type/neighbourlist.f90
index 3c8b9edd8..dadad847c 100644
--- a/src/type/neighbourlist.f90
+++ b/src/type/neighbourlist.f90
@@ -23,6 +23,7 @@ module xtb_type_neighbourlist
use xtb_mctc_search, only : bisectSearch
use xtb_mctc_sort, only : indexHeapSort
use xtb_mctc_thresholds, only : tolSameDist2, minNeighDist2
+ use xtb_mctc_param, only : covalent_radius_d3
use xtb_type_environment, only : TEnvironment
implicit none
private
@@ -62,6 +63,8 @@ module xtb_type_neighbourlist
!> Generate neighbour list
procedure :: generate
+ !> Generate neighbour list from covalent radii (nonperiodic)
+ procedure :: generate_covalent
!> Update neighbour list
procedure :: update
@@ -145,6 +148,74 @@ subroutine generate(self, env, coords0, cutoff, latPoint, symmetric)
end subroutine generate
+
+!> Generate neighbour list from covalent radii (nonperiodic).
+subroutine generate_covalent(self, at, coords0)
+ !> Instance of the neighbour list
+ class(TNeighbourList), intent(inout) :: self
+ !> Atomic numbers
+ integer, intent(in) :: at(:)
+ !> Coordinates of the atoms in the central cell
+ real(wp), intent(in) :: coords0(:, :)
+
+ ! Number of atoms in the central cell
+ integer :: nAtom
+
+ ! Capacity of the neighbour related arrays
+ integer :: mNeigh
+
+ integer :: iAt, jAt
+ real(wp) :: r2, dist2
+
+ nAtom = size(at)
+ if (nAtom == 0) then
+ self%cutoff = 0.0_wp
+ return
+ end if
+
+ ! Reset the list to the canonical central-cell half-list
+ self%neighs(:) = 0
+ self%iNeigh(:, :) = 0
+ self%dist2(:, :) = 0.0_wp
+ self%weight(:, :) = 0.0_wp
+ do iAt = 1, nAtom
+ self%image(iAt) = iAt
+ self%trans(iAt) = 1
+ self%coords(:, iAt) = coords0(:, iAt)
+ self%iNeigh(0, iAt) = iAt
+ self%weight(0, iAt) = 1.0_wp
+ end do
+
+ mNeigh = ubound(self%iNeigh, dim=1)
+
+ do iAt = 1, nAtom
+ do jAt = iAt + 1, nAtom
+ r2 = (covalent_radius_d3(at(iAt)) + covalent_radius_d3(at(jAt)))**2
+ dist2 = sum((coords0(:, jAt) - coords0(:, iAt))**2)
+ if (dist2 >= r2) cycle
+
+ self%neighs(iAt) = self%neighs(iAt) + 1
+ if (self%neighs(iAt) > mNeigh) then
+ mNeigh = 2*mNeigh
+ call resizeNeigh(mNeigh, self%iNeigh, self%dist2, self%weight)
+ end if
+ self%iNeigh(self%neighs(iAt), iAt) = jAt
+ self%dist2(self%neighs(iAt), iAt) = dist2
+ self%weight(self%neighs(iAt), iAt) = 1.0_wp
+ end do
+ end do
+ if (maxval(self%neighs) == 0) then
+ self%cutoff = 0.0_wp
+ else
+ self%cutoff = sqrt((2.0_wp*maxval(covalent_radius_d3(at)))**2 &
+ & + 2*tolSameDist2)
+ end if
+
+ call self%sort
+
+end subroutine generate_covalent
+
+
!> Generate neighbour list in serial
subroutine generateSequential(self, env, coords0, cutoff, latPoint, symmetric)
diff --git a/src/xhelp.f90 b/src/xhelp.f90
index b2da88906..1ce372f0f 100644
--- a/src/xhelp.f90
+++ b/src/xhelp.f90
@@ -283,6 +283,10 @@ subroutine help(iunit)
"--bhess [LEVEL]",&
" perform a biased numerical hessian calculation on an ancopt(3) optimized geometry",&
"",&
+ " Molecular Hessian runs with more than one atom, no frozen atoms, and no periodicity",&
+ " automatically print compliance constants, relaxed force constants, and bond local-mode",&
+ " frequencies, and write compliance.dat.",&
+ "",&
"--o1nh",&
" perform the numerical hessian calculation using the ODLR approximation (O1NumHess)",&
"",&
diff --git a/test/unit/test_bmatrix.f90 b/test/unit/test_bmatrix.f90
index 35566b0d2..582a22078 100644
--- a/test/unit/test_bmatrix.f90
+++ b/test/unit/test_bmatrix.f90
@@ -16,18 +16,20 @@
! along with xtb. If not, see .
!> Unit tests for the xtb_bmatrix module (Wilson B-matrix derivatives).
-!> Most tests use central finite differences to validate analytic B rows;
-!> linear-bend uses translational invariance and orthogonality checks instead.
+!> Every analytic B row is validated against central finite differences of the
+!> coordinate it differentiates; the linear-bend pair additionally gets
+!> translational-invariance and orthogonality checks on its two frame rows.
module test_bmatrix
use testdrive, only : new_unittest, unittest_type, error_type, check
use xtb_mctc_accuracy, only : wp
- use xtb_bmatrix, only : bmat_bond, bmat_angle, bmat_linbend, &
+ use xtb_mctc_constants, only : pi
+ use xtb_bmatrix, only : bmat_bond, bmat_angle, bmat_linbend, linbend_frame, &
& bmat_torsion, bmat_outofplane
implicit none
private
real(wp), parameter :: fd_eps = 1.0e-6_wp
- real(wp), parameter :: thr = 100*epsilon(0.0_wp)
+ real(wp), parameter :: thr = 100 * epsilon(0.0_wp)
real(wp), parameter :: thr_loose = 1.0e-9_wp
public :: collect_bmatrix
@@ -36,24 +38,29 @@ module test_bmatrix
!> Collect all exported unit tests
subroutine collect_bmatrix(testsuite)
+ !> Tests exported to the harness.
type(unittest_type), allocatable, intent(out) :: testsuite(:)
- testsuite = [ &
- new_unittest("bond-fd", test_bond_fd), &
- new_unittest("angle-fd", test_angle_fd), &
- new_unittest("linbend-invariance", test_linbend_fd), &
- new_unittest("torsion-fd", test_torsion_fd), &
- new_unittest("torsion-sign", test_torsion_sign), &
- new_unittest("torsion-atom-order", test_torsion_atom_order), &
- new_unittest("outofplane-fd", test_outofplane_fd), &
- new_unittest("outofplane-degenerate", test_outofplane_degenerate) &
- ]
+ testsuite = [ &
+ new_unittest("bond-fd", test_bond_fd), &
+ new_unittest("angle-fd", test_angle_fd), &
+ new_unittest("angle-near-antiparallel", test_angle_near_antiparallel), &
+ new_unittest("linbend-fd", test_linbend_fd), &
+ new_unittest("linbend-invariance", test_linbend_invariance), &
+ new_unittest("torsion-fd", test_torsion_fd), &
+ new_unittest("torsion-sign", test_torsion_sign), &
+ new_unittest("torsion-atom-order", test_torsion_atom_order), &
+ new_unittest("outofplane-fd", test_outofplane_fd), &
+ new_unittest("outofplane-degenerate", test_outofplane_degenerate) &
+ ]
end subroutine collect_bmatrix
-!> Finite-difference check of bmat_bond
+!> Pins the analytic bond row of bmat_bond against central finite differences
+!> of the bond length for every Cartesian displacement.
subroutine test_bond_fd(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
real(wp), parameter :: xyz2(3, 2) = reshape([ &
@@ -66,7 +73,7 @@ subroutine test_bond_fd(error)
fd = 0.0_wp
do k = 1, 2
do c = 1, 3
- fd(c + (k-1)*3) = fd_bond(xyz2, c, k)
+ fd(c+(k-1)*3) = fd_bond(xyz2, c, k)
end do
end do
call check_bmat(error, bmat, fd, thr_loose)
@@ -74,10 +81,12 @@ end subroutine test_bond_fd
!> Central FD for bond B-row entry (c, k)
function fd_bond(xyz0, c, k) result(dval)
+ !> Reference geometry, two atoms.
real(wp), intent(in) :: xyz0(3, 2)
+ !> Displaced Cartesian component and atom index.
integer, intent(in) :: c, k
- real(wp) :: dval
+ real(wp) :: dval
real(wp) :: xyz(3, 2), rp, rm
xyz = xyz0
@@ -86,26 +95,30 @@ function fd_bond(xyz0, c, k) result(dval)
xyz = xyz0
xyz(c, k) = xyz(c, k) - fd_eps
rm = bond_val(xyz)
- dval = (rp - rm) / (2.0_wp * fd_eps)
+ dval = (rp - rm) / (2.0_wp*fd_eps)
end function fd_bond
!> Bond length as scalar function of xyz
function bond_val(xyz) result(val)
+ !> Reference geometry, two atoms.
real(wp), intent(in) :: xyz(3, 2)
+
real(wp) :: val
val = norm2(xyz(:, 2) - xyz(:, 1))
end function bond_val
-!> Finite-difference check of bmat_angle
+!> Pins the analytic angle row of bmat_angle against central finite differences
+!> of the valence angle for every Cartesian displacement.
subroutine test_angle_fd(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
! 3-atom fragment: j at origin, i and k on either side
real(wp), parameter :: xyz3(3, 3) = reshape([ &
- & 1.0_wp, 0.0_wp, 0.0_wp, &
- & 0.0_wp, 0.0_wp, 0.0_wp, &
- & 0.5_wp, 0.8_wp, 0.1_wp], shape(xyz3))
+ & 1.0_wp, 0.0_wp, 0.0_wp, &
+ & 0.0_wp, 0.0_wp, 0.0_wp, &
+ & 0.5_wp, 0.8_wp, 0.1_wp], shape(xyz3))
real(wp) :: bmat(9), fd(9)
integer :: k, c
@@ -114,18 +127,47 @@ subroutine test_angle_fd(error)
fd = 0.0_wp
do k = 1, 3
do c = 1, 3
- fd(c + (k-1)*3) = fd_angle(xyz3, c, k)
+ fd(c+(k-1)*3) = fd_angle(xyz3, c, k)
end do
end do
call check_bmat(error, bmat, fd, thr_loose)
end subroutine test_angle_fd
+!> Near-antiparallel angles retain their finite limiting row and inverse-length
+!> scaling even when rounding makes the normalized dot product exactly -1.
+subroutine test_angle_near_antiparallel(error)
+ use, intrinsic :: ieee_arithmetic, only : ieee_is_finite
+ !> Failure report, allocated when the check fails.
+ type(error_type), allocatable, intent(out) :: error
+
+ real(wp), parameter :: delta = 1.0e-9_wp
+ real(wp), parameter :: scale = 8.0_wp
+ real(wp), parameter :: vec1(3) = [1.0_wp, 0.0_wp, 0.0_wp]
+ real(wp), parameter :: vec2(3) = [-1.0_wp, delta, 0.0_wp]
+ real(wp), parameter :: expected(9) = [ &
+ & 0.0_wp, -1.0_wp, 0.0_wp, &
+ & 0.0_wp, 2.0_wp, 0.0_wp, &
+ & 0.0_wp, -1.0_wp, 0.0_wp]
+ real(wp) :: bmat(9), scaled_bmat(9)
+
+ bmat = bmat_angle(vec1, vec2)
+ scaled_bmat = bmat_angle(scale*vec1, scale*vec2)
+
+ call check(error, all(ieee_is_finite(bmat)))
+ if (allocated(error)) return
+ call check_bmat(error, bmat, expected, thr_loose)
+ if (allocated(error)) return
+ call check_bmat(error, scaled_bmat, expected/scale, thr_loose)
+end subroutine test_angle_near_antiparallel
+
!> Central FD for angle B-row entry (c, k)
function fd_angle(xyz0, c, k) result(dval)
+ !> Reference geometry, three atoms.
real(wp), intent(in) :: xyz0(3, 3)
+ !> Displaced Cartesian component and atom index.
integer, intent(in) :: c, k
- real(wp) :: dval
+ real(wp) :: dval
real(wp) :: xyz(3, 3), rp, rm
xyz = xyz0
@@ -134,48 +176,48 @@ function fd_angle(xyz0, c, k) result(dval)
xyz = xyz0
xyz(c, k) = xyz(c, k) - fd_eps
rm = angle_val(xyz)
- dval = (rp - rm) / (2.0_wp * fd_eps)
+ dval = (rp - rm) / (2.0_wp*fd_eps)
end function fd_angle
!> Valence angle (i-j-k) as scalar function of xyz
function angle_val(xyz) result(val)
+ !> Reference geometry, three atoms.
real(wp), intent(in) :: xyz(3, 3)
+
real(wp) :: val
real(wp) :: v1(3), v2(3), co
v1 = xyz(:, 1) - xyz(:, 2)
v2 = xyz(:, 3) - xyz(:, 2)
- co = dot_product(v1, v2) / (norm2(v1) * norm2(v2))
+ co = dot_product(v1, v2) / (norm2(v1)*norm2(v2))
co = max(-1.0_wp, min(1.0_wp, co))
val = acos(co)
end function angle_val
-!> Linear-bend B-matrix sanity check.
-!> bmat_linbend returns a linearized B-matrix valid at the linear geometry;
-!> finite-difference validation is not meaningful away from that point.
-!> Checks: (1) per-Cartesian-component translational invariance, i.e.
-!> for each row the sum over the x-component entries, y-component entries,
-!> and z-component entries must each be zero independently;
-!> (2) both rows must be nonzero;
-!> (3) the two rows must be mutually orthogonal.
-subroutine test_linbend_fd(error)
+!> Pins translational invariance and mutual orthogonality of the two rows a
+!> linear centre contributes: the linbend_frame directions through bmat_linbend.
+subroutine test_linbend_invariance(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
- ! Near-linear 3-atom fragment: j at origin, i and k nearly collinear
+ ! Near-linear 3-atom fragment: centre at the origin, atoms 1 and 3 collinear
real(wp), parameter :: xyz_lin(3, 3) = reshape([ &
& 1.0_wp, 0.0_wp, 0.0_wp, &
& 0.0_wp, 0.0_wp, 0.0_wp, &
& -0.99_wp, 0.01_wp, 0.0_wp], shape(xyz_lin))
- real(wp) :: bmat(2, 9)
+ real(wp) :: vec1(3), vec2(3), evec1(3), evec2(3), bmat(2, 9)
real(wp) :: row_sum_x, row_sum_y, row_sum_z
real(wp) :: dot_rows
integer :: row, c, nonzero_count
- bmat = bmat_linbend(xyz_lin(:, 1) - xyz_lin(:, 2), &
- & xyz_lin(:, 3) - xyz_lin(:, 2))
+ vec1 = xyz_lin(:, 1) - xyz_lin(:, 2)
+ vec2 = xyz_lin(:, 3) - xyz_lin(:, 2)
+ call linbend_frame(vec1/norm2(vec1), evec1, evec2)
+ bmat(1, :) = bmat_linbend(vec1, vec2, evec1)
+ bmat(2, :) = bmat_linbend(vec1, vec2, evec2)
- ! (1) Per-component translational invariance:
- ! atom k contributes entries 1:3, atom j entries 4:6, atom i entries 7:9
+ ! Per-component translational invariance:
+ ! atom 1 contributes entries 1:3, the centre 4:6, atom 2 entries 7:9
! For each Cartesian component (x=1,4,7; y=2,5,8; z=3,6,9) the row sum must be 0
do row = 1, 2
row_sum_x = bmat(row, 1) + bmat(row, 4) + bmat(row, 7)
@@ -192,7 +234,7 @@ subroutine test_linbend_fd(error)
if (allocated(error)) return
end do
- ! (2) Both rows must be nonzero
+ ! Both rows must be nonzero
do row = 1, 2
nonzero_count = 0
do c = 1, 9
@@ -205,15 +247,89 @@ subroutine test_linbend_fd(error)
end if
end do
- ! (3) The two rows should be mutually orthogonal
+ ! The two rows should be mutually orthogonal
dot_rows = dot_product(bmat(1, :), bmat(2, :))
call check(error, dot_rows, 0.0_wp, thr=thr, &
& message="linbend rows not orthogonal")
+end subroutine test_linbend_invariance
+
+
+!> Pins the analytic linear-bend row of bmat_linbend: the row must be the exact
+!> derivative of q = evec . v1_hat + evec . v2_hat at a near-linear centre that
+!> is not aligned with a coordinate axis.
+subroutine test_linbend_fd(error)
+ !> Failure report, allocated when the check fails.
+ type(error_type), allocatable, intent(out) :: error
+
+ ! unit axis u and a perpendicular unit vector p; the centre sits at the origin
+ real(wp), parameter :: u(3) = [0.6_wp, 0.8_wp, 0.0_wp]
+ real(wp), parameter :: p(3) = [-0.8_wp, 0.6_wp, 0.0_wp]
+ real(wp), parameter :: zax(3) = [0.0_wp, 0.0_wp, 1.0_wp]
+ real(wp) :: xyz(3, 3), evec(3), bmat(9), fd(9)
+ real(wp) :: dev
+ integer :: k, c
+
+ ! bond lengths 1.2 and 1.5, one degree off collinearity along p
+ dev = 1.5_wp * tan(pi/180.0_wp)
+ xyz = 0.0_wp
+ xyz(:, 1) = -1.2_wp * u
+ xyz(:, 3) = 1.5_wp * u + dev * p
+
+ ! fixed frame: exactly perpendicular to vec1, tilted out of the u/p plane
+ evec = p + 0.5_wp * zax
+ evec = evec / norm2(evec)
+
+ ! layout [atom1 (1:3), centre (4:6), atom2 (7:9)]
+ bmat = bmat_linbend(xyz(:, 1) - xyz(:, 2), &
+ & xyz(:, 3) - xyz(:, 2), evec)
+ fd = 0.0_wp
+ do k = 1, 3
+ do c = 1, 3
+ fd(c+(k-1)*3) = fd_linbend(xyz, c, k, evec)
+ end do
+ end do
+ call check_bmat(error, bmat, fd, thr_loose)
end subroutine test_linbend_fd
+!> Central FD for the linear-bend row entry (c, k)
+function fd_linbend(xyz0, c, k, evec) result(dval)
+ !> Reference geometry, three atoms.
+ !> Fixed unit frame direction for the linear-bend coordinate.
+ real(wp), intent(in) :: xyz0(3, 3), evec(3)
+ !> Displaced Cartesian component and atom index.
+ integer, intent(in) :: c, k
+
+ real(wp) :: dval
+ real(wp) :: xyz(3, 3), rp, rm
+
+ xyz = xyz0
+ xyz(c, k) = xyz(c, k) + fd_eps
+ rp = linbend_val(xyz, evec)
+ xyz = xyz0
+ xyz(c, k) = xyz(c, k) - fd_eps
+ rm = linbend_val(xyz, evec)
+ dval = (rp - rm) / (2.0_wp*fd_eps)
+end function fd_linbend
+
+!> Linear-bend coordinate value as scalar function of xyz
+function linbend_val(xyz, evec) result(val)
+ !> Reference geometry, three atoms.
+ !> Fixed unit vector defining the linear-bend coordinate.
+ real(wp), intent(in) :: xyz(3, 3), evec(3)
+
+ real(wp) :: val
+ real(wp) :: v1(3), v2(3)
-!> Finite-difference check of bmat_torsion
+ v1 = xyz(:, 1) - xyz(:, 2)
+ v2 = xyz(:, 3) - xyz(:, 2)
+ val = dot_product(evec, v1) / norm2(v1) &
+ & + dot_product(evec, v2) / norm2(v2)
+end function linbend_val
+
+!> Pins the analytic torsion row of bmat_torsion against central finite
+!> differences of the dihedral angle for every Cartesian displacement.
subroutine test_torsion_fd(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
real(wp), parameter :: xyz4(3, 4) = reshape([ &
@@ -221,7 +337,7 @@ subroutine test_torsion_fd(error)
& 1.0_wp, 0.0_wp, 0.0_wp, &
& 1.5_wp, 0.8_wp, 0.2_wp, &
& 2.0_wp, 0.9_wp, -0.5_wp], shape(xyz4))
- real(wp) :: tau, bt(3, 4), fd(3, 4)
+ real(wp) :: bt(3, 4), fd(3, 4)
integer :: k, c
bt = bmat_torsion(xyz4)
@@ -236,10 +352,12 @@ end subroutine test_torsion_fd
!> Central FD for torsion B-row entry (c, k)
function fd_torsion(xyz0, c, k) result(dval)
+ !> Reference geometry, four atoms.
real(wp), intent(in) :: xyz0(3, 4)
+ !> Displaced Cartesian component and atom index.
integer, intent(in) :: c, k
- real(wp) :: dval
+ real(wp) :: dval
real(wp) :: xyz(3, 4), rp, rm
xyz = xyz0
@@ -248,14 +366,15 @@ function fd_torsion(xyz0, c, k) result(dval)
xyz = xyz0
xyz(c, k) = xyz(c, k) - fd_eps
rm = torsion_val(xyz)
- dval = (rp - rm) / (2.0_wp * fd_eps)
+ dval = (rp - rm) / (2.0_wp*fd_eps)
end function fd_torsion
!> Dihedral angle as scalar function of xyz (4 atoms)
function torsion_val(xyz) result(val)
+ !> Reference geometry, four atoms.
real(wp), intent(in) :: xyz(3, 4)
- real(wp) :: val
+ real(wp) :: val
real(wp) :: v1(3), v2(3), v3(3)
real(wp) :: n1(3), n2(3), s, c
@@ -270,10 +389,10 @@ function torsion_val(xyz) result(val)
end function torsion_val
-!> Torsion sign/branch-wrap check: two geometries that differ only by
-!> flipping atom 1 to the other side of the plane should produce
-!> dihedrals that differ in sign (or wrap from +pi to -pi).
+!> Pins the torsion sign convention: two geometries that differ only by flipping
+!> atom 1 across the i-j-k plane must give dihedrals of opposite sign.
subroutine test_torsion_sign(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
! Geometry A: atom 1 above the i-j-k plane
@@ -299,9 +418,10 @@ subroutine test_torsion_sign(error)
end subroutine test_torsion_sign
-!> Torsion atom-order check: reversing the atom ordering (4-3-2-1 vs 1-2-3-4)
-!> should give the same |tau| and consistent B-matrix relationship.
+!> Pins torsion invariance under atom reversal: the reversed order 4-3-2-1 gives
+!> the same |tau| as 1-2-3-4.
subroutine test_torsion_atom_order(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
real(wp), parameter :: xyz_fwd(3, 4) = reshape([ &
@@ -315,7 +435,7 @@ subroutine test_torsion_atom_order(error)
! Reverse atom order: atom 4 becomes 1, 3 becomes 2, etc.
do k = 1, 4
- xyz_rev(:, k) = xyz_fwd(:, 5 - k)
+ xyz_rev(:, k) = xyz_fwd(:, 5-k)
end do
tau_fwd = torsion_val(xyz_fwd)
@@ -327,14 +447,16 @@ subroutine test_torsion_atom_order(error)
end subroutine test_torsion_atom_order
-!> Finite-difference check of bmat_outofplane
+!> Pins the analytic out-of-plane row of bmat_outofplane against central finite
+!> differences of the out-of-plane angle for every Cartesian displacement.
subroutine test_outofplane_fd(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
! 4-atom fragment: atom 4 central, atom 1 out of plane, 2 and 3 in plane
real(wp), parameter :: xyz4(3, 4) = reshape([ &
- & 0.1_wp, 0.0_wp, 0.5_wp, & ! atom 1 (out of plane, tilted)
- & 1.0_wp, 0.0_wp, 0.0_wp, & ! atom 2 (in plane)
+ & 0.1_wp, 0.0_wp, 0.5_wp, & ! atom 1 (out of plane, tilted)
+ & 1.0_wp, 0.0_wp, 0.0_wp, & ! atom 2 (in plane)
& -0.5_wp, 0.8_wp, 0.0_wp, & ! atom 3 (in plane)
& 0.0_wp, 0.0_wp, 0.0_wp], & ! atom 4 (central)
& shape(xyz4))
@@ -348,20 +470,21 @@ subroutine test_outofplane_fd(error)
fd(c, k) = fd_outofplane(xyz4, c, k)
end do
end do
- call check_bmat_2d_t(error, bt, fd, thr_loose)
+ call check_bmat_2d_t(error, bt, fd, thr_loose)
end subroutine test_outofplane_fd
-!> Degenerate-plane check: when atoms 2, 3, 4 are collinear the
-!> out-of-plane angle is undefined; bmat_outofplane should return
-!> teta=0 and bt=0 without crashing.
+!> Pins the degenerate-plane branch: with atoms 2, 3, 4 collinear the
+!> out-of-plane angle is undefined; bmat_outofplane must return teta=0 and
+!> bt=0 without crashing.
subroutine test_outofplane_degenerate(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
! Atoms 2, 3, 4 collinear along x-axis; atom 1 out of plane
real(wp), parameter :: xyz_deg(3, 4) = reshape([ &
- & 0.0_wp, 0.0_wp, 0.5_wp, & ! atom 1 (out of plane)
- & 1.0_wp, 0.0_wp, 0.0_wp, & ! atom 2 (collinear with 3, 4)
- & 2.0_wp, 0.0_wp, 0.0_wp, & ! atom 3 (collinear with 2, 4)
+ & 0.0_wp, 0.0_wp, 0.5_wp, & ! atom 1 (out of plane)
+ & 1.0_wp, 0.0_wp, 0.0_wp, & ! atom 2 (collinear with 3, 4)
+ & 2.0_wp, 0.0_wp, 0.0_wp, & ! atom 3 (collinear with 2, 4)
& 0.0_wp, 0.0_wp, 0.0_wp], & ! atom 4 (central, collinear)
& shape(xyz_deg))
real(wp) :: bt(3, 4)
@@ -374,10 +497,12 @@ end subroutine test_outofplane_degenerate
!> Central FD for out-of-plane B-row entry (c, k)
function fd_outofplane(xyz0, c, k) result(dval)
+ !> Reference geometry, four atoms with atom 4 central.
real(wp), intent(in) :: xyz0(3, 4)
+ !> Displaced Cartesian component and atom index.
integer, intent(in) :: c, k
- real(wp) :: dval
+ real(wp) :: dval
real(wp) :: xyz(3, 4), rp, rm
xyz = xyz0
@@ -386,14 +511,15 @@ function fd_outofplane(xyz0, c, k) result(dval)
xyz = xyz0
xyz(c, k) = xyz(c, k) - fd_eps
rm = outofplane_val(xyz)
- dval = (rp - rm) / (2.0_wp * fd_eps)
+ dval = (rp - rm) / (2.0_wp*fd_eps)
end function fd_outofplane
!> Out-of-plane angle as scalar function of xyz (4 atoms, atom 4 central)
function outofplane_val(xyz) result(val)
+ !> Reference geometry, four atoms with atom 4 central.
real(wp), intent(in) :: xyz(3, 4)
- real(wp) :: val
+ real(wp) :: val
real(wp) :: r42(3), r43(3), r41(3)
real(wp) :: n(3), nn, h, r41n
@@ -413,24 +539,31 @@ function outofplane_val(xyz) result(val)
val = 0.0_wp
return
end if
- val = -asin(max(-1.0_wp, min(1.0_wp, h / r41n)))
+ val = -asin(max(-1.0_wp, min(1.0_wp, h/r41n)))
end function outofplane_val
!> Cross product of two 3-vectors
function cross(a, b) result(c)
+ !> First vector.
+ !> Second vector.
real(wp), intent(in) :: a(3), b(3)
+
real(wp) :: c(3)
- c(1) = a(2)*b(3) - a(3)*b(2)
- c(2) = a(3)*b(1) - a(1)*b(3)
- c(3) = a(1)*b(2) - a(2)*b(1)
+ c(1) = a(2) * b(3) - a(3) * b(2)
+ c(2) = a(3) * b(1) - a(1) * b(3)
+ c(3) = a(1) * b(2) - a(2) * b(1)
end function cross
!> Compare two 1D B-row arrays element-wise
subroutine check_bmat(error, actual, expected, tol)
+ !> Failure report, allocated when a comparison fails.
type(error_type), allocatable, intent(out) :: error
+ !> Analytic row and finite-difference reference.
real(wp), intent(in) :: actual(:), expected(:)
+ !> Absolute tolerance for the comparison.
real(wp), intent(in) :: tol
+
integer :: i
if (size(actual) /= size(expected)) then
@@ -443,31 +576,15 @@ subroutine check_bmat(error, actual, expected, tol)
end do
end subroutine check_bmat
-!> Compare two 2D B-row arrays (shape 2, N) element-wise
-subroutine check_bmat_2d(error, actual, expected, tol)
- type(error_type), allocatable, intent(out) :: error
- real(wp), intent(in) :: actual(:, :), expected(:, :)
- real(wp), intent(in) :: tol
- integer :: i, j
-
- if (size(actual, 1) /= size(expected, 1) .or. &
- & size(actual, 2) /= size(expected, 2)) then
- call check(error, size(actual, 1), size(expected, 1))
- return
- end if
- do j = 1, size(actual, 2)
- do i = 1, size(actual, 1)
- call check(error, actual(i, j), expected(i, j), thr=tol)
- if (allocated(error)) return
- end do
- end do
-end subroutine check_bmat_2d
-
!> Compare two 2D B-row arrays (shape 3, N) element-wise — transposed view
subroutine check_bmat_2d_t(error, actual, expected, tol)
+ !> Failure report, allocated when a comparison fails.
type(error_type), allocatable, intent(out) :: error
+ !> Analytic array and finite-difference reference, shape (3, N).
real(wp), intent(in) :: actual(:, :), expected(:, :)
+ !> Absolute tolerance for the comparison.
real(wp), intent(in) :: tol
+
integer :: i, j
if (size(actual, 2) /= size(expected, 2)) then
diff --git a/test/unit/test_hessian.f90 b/test/unit/test_hessian.f90
index 6a2a84d2e..492f407c2 100644
--- a/test/unit/test_hessian.f90
+++ b/test/unit/test_hessian.f90
@@ -16,11 +16,19 @@
! You should have received a copy of the GNU Lesser General Public License
! along with xtb. If not, see .
+!> Hessian, compliance and redundant-internal behaviour locks
+!>
+!> Numerical and O1NumHess (ODLR) Hessians of water are pinned against reference
+!> tables and against their eigenvalue spectra, the GFN1 compliance constants of
+!> water are pinned as a reference matrix, and the analytic B-matrix of a
+!> non-collinear redundant coordinate set is validated against central finite
+!> differences.
module test_hessian
use testdrive, only : new_unittest, unittest_type, error_type, check, test_failed
use xtb_mctc_accuracy, only : wp
use xtb_mctc_io, only : stdout
use xtb_mctc_convert, only : autoaa
+ use xtb_mctc_constants, only : pi
use xtb_type_options
use xtb_type_molecule
use xtb_type_restart
@@ -30,12 +38,19 @@ module test_hessian
use xtb_xtb_calculator, only : TxTBCalculator
use xtb_main_setup, only : newXTBCalculator, newWavefunction
+ use xtb_compliance, only : compliance_driver, compute_compliance
+ use xtb_type_neighbourlist, only : TNeighbourList, init
+ use xtb_mctc_param, only : covalent_radius_d3
+ use xtb_internals_graph, only : graph_type, init
+ use xtb_internals_redundant, only : redundant_type, init
+ use xtb_internals_type, only : coord_bond, coord_angle, coord_dihedral, coord_linbend
+ use xtb_bmatrix, only : get_bmatrix
+ use xtb_o1numhess, only : lr_loop
implicit none
private
public :: collect_hessian
-
contains
!> Collect all exported unit tests
@@ -49,60 +64,100 @@ subroutine collect_hessian(testsuite)
new_unittest("gfn1_o1numhess", test_o1numhess_gfn1), &
new_unittest("gfn2_o1numhess", test_o1numhess_gfn2), &
new_unittest("linear_h2o_gfn1_o1numhess", test_o1numhess_linear_h2o_gfn1), &
- new_unittest("linear_h2o_gfn2_o1numhess", test_o1numhess_linear_h2o_gfn2) &
+ new_unittest("linear_h2o_gfn2_o1numhess", test_o1numhess_linear_h2o_gfn2), &
+ new_unittest("rectangular_lr_loop", test_rectangular_lr_loop), &
+ new_unittest("compliance", test_compliance), &
+ new_unittest("compliance_water", test_compliance_water), &
+ new_unittest("covalent_neighbour_list", test_covalent_neighbour_list), &
+ new_unittest("redundant_bmatrix_fd", test_redundant_bmatrix_fd), &
+ new_unittest("linear_chain_compliance", test_linear_chain_compliance), &
+ new_unittest("near_linear_bmatrix_fd", test_near_linear_bmatrix_fd), &
+ new_unittest("compliance_driver_redundant", test_compliance_driver_redundant) &
]
end subroutine collect_hessian
+!> A rectangular displacement set must produce a symmetric correction without
+!> overrunning the GEMM workspace; final_err must describe the returned matrix.
+subroutine test_rectangular_lr_loop(error)
+ type(error_type), allocatable, intent(out) :: error
+
+ real(wp), parameter :: scale = 1.0e-4_wp
+ real(wp), parameter :: displdir(2, 1) = reshape([1.0_wp, 0.0_wp], [2, 1])
+ type(TEnvironment) :: env
+ real(wp) :: gradient(2, 1), hessian(2, 2), final_err, residual
+
+ call init(env)
+ gradient(:, 1) = [scale, 2.0_wp*scale]
+ hessian = 0.0_wp
+ call lr_loop(env, 1, gradient, hessian, displdir, final_err)
+
+ call check(error, hessian(1, 1)/scale, 1.0_wp, thr=1.0e-4_wp)
+ if (allocated(error)) return
+ call check(error, hessian(2, 1)/scale, 2.0_wp, thr=1.0e-4_wp)
+ if (allocated(error)) return
+ call check(error, hessian(1, 2), hessian(2, 1), thr=scale*1.0e-12_wp)
+ if (allocated(error)) return
+ call check(error, hessian(2, 2), 0.0_wp, thr=scale*1.0e-12_wp)
+ if (allocated(error)) return
+
+ residual = norm2(gradient - matmul(hessian, displdir))/norm2(gradient)
+ call check(error, final_err, residual, thr=1.0e-14_wp)
+end subroutine test_rectangular_lr_loop
+
+!> Pins the GFN1 numerical Hessian and dipole gradient of water against the
+!> reference tables below.
subroutine test_gfn1_hessian(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
+
integer, parameter :: nat = 3
- real(wp),parameter :: thr = 1.0e-7_wp
+ real(wp), parameter :: thr = 1.0e-7_wp
character(len=*), parameter :: sym(nat) = ["O", "H", "H"]
real(wp), parameter :: xyz(3, nat) = reshape([&
- & 0.00000000000000_wp, 0.00000000034546_wp, 0.18900383618455_wp, &
- & 0.00000000000000_wp, 1.45674735348811_wp, -0.88650486059828_wp, &
- &-0.00000000000000_wp, -1.45674735383357_wp, -0.88650486086986_wp],&
+ & 0.00000000000000_wp, 0.00000000034546_wp, 0.18900383618455_wp, &
+ & 0.00000000000000_wp, 1.45674735348811_wp, -0.88650486059828_wp, &
+ &-0.00000000000000_wp, -1.45674735383357_wp, -0.88650486086986_wp], &
& shape(xyz))
real(wp), parameter :: dipgrad_ref(3, 3*nat) = reshape([ &
- & -1.013452580143500E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & -1.013452580143500E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
& 3.527708416572683E-11_wp, -0.580105714422270E+00_wp, -4.440892098500626E-10_wp, &
& 1.640601522247707E-11_wp, -1.577916091825909E-09_wp, -0.719923066050399E+00_wp, &
- & 0.506726290141942E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 5.168504173660994E-11_wp, 0.290827888051893E+00_wp, 7.276368407804767E-02_wp, &
- & -3.527728552279606E-11_wp, 0.159475763089199E+00_wp, 0.359943068684032E+00_wp, &
- & 0.506726290062044E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 0.000000000000000E+00_wp, 0.290827888619521E+00_wp, -7.276368707564984E-02_wp, &
- & 0.000000000000000E+00_wp, -0.159475761295180E+00_wp, 0.359943068239943E+00_wp],&
+ & 0.506726290141942E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 5.168504173660994E-11_wp, 0.290827888051893E+00_wp, 7.276368407804767E-02_wp, &
+ & -3.527728552279606E-11_wp, 0.159475763089199E+00_wp, 0.359943068684032E+00_wp, &
+ & 0.506726290062044E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 0.000000000000000E+00_wp, 0.290827888619521E+00_wp, -7.276368707564984E-02_wp, &
+ & 0.000000000000000E+00_wp, -0.159475761295180E+00_wp, 0.359943068239943E+00_wp], &
& shape(dipgrad_ref))
real(wp), parameter :: hessian_ref(3*nat, 3*nat) = reshape([ &
- & 9.642724315151719E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & -4.821349241627352E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & -4.821375073522908E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & -5.183734810089242E-12_wp, 0.654328674106580E+00_wp, -8.411850317864566E-10_wp, &
- & 4.912366132678069E-12_wp, -0.327164337503372E+00_wp, 0.241508066094075E+00_wp, &
+ & 9.642724315151719E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & -4.821349241627352E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & -4.821375073522908E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & -5.183734810089242E-12_wp, 0.654328674106580E+00_wp, -8.411850317864566E-10_wp, &
+ & 4.912366132678069E-12_wp, -0.327164337503372E+00_wp, 0.241508066094075E+00_wp, &
& 2.713686774111726E-13_wp, -0.327164336592530E+00_wp, -0.241508065247398E+00_wp, &
- & 1.269389491007523E-11_wp, -1.003642487147619E-09_wp, 0.396956797942315E+00_wp, &
- & -1.357110357831859E-11_wp, 0.193314810701155E+00_wp, -0.198478399246260E+00_wp, &
+ & 1.269389491007523E-11_wp, -1.003642487147619E-09_wp, 0.396956797942315E+00_wp, &
+ & -1.357110357831859E-11_wp, 0.193314810701155E+00_wp, -0.198478399246260E+00_wp, &
& 8.772086682433594E-13_wp, -0.193314809690279E+00_wp, -0.198478398680934E+00_wp, &
- & -4.821349763250659E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 4.418957886905253E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 4.023918763456385E-06_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 7.510234345234422E-12_wp, -0.327069271103985E+00_wp, 0.193419608180277E+00_wp, &
- & -8.658575901216439E-12_wp, 0.348666988822030E+00_wp, -0.217428743767287E+00_wp, &
- & 1.148341555982017E-12_wp, -2.159771772581407E-02_wp, 2.400913557787203E-02_wp, &
- & 5.183777743453684E-12_wp, 0.241435744915371E+00_wp, -0.198481745303595E+00_wp, &
- & -4.912396810694304E-12_wp, -0.217376911834996E+00_wp, 0.188386308225154E+00_wp, &
- & -2.713809327593799E-13_wp, -2.405883307999757E-02_wp, 1.009543707641887E-02_wp, &
- & -4.821379744720639E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 4.023966094020255E-06_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 4.418983135319131E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & -4.821349763250659E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 4.418957886905253E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 4.023918763456385E-06_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 7.510234345234422E-12_wp, -0.327069271103985E+00_wp, 0.193419608180277E+00_wp, &
+ & -8.658575901216439E-12_wp, 0.348666988822030E+00_wp, -0.217428743767287E+00_wp, &
+ & 1.148341555982017E-12_wp, -2.159771772581407E-02_wp, 2.400913557787203E-02_wp, &
+ & 5.183777743453684E-12_wp, 0.241435744915371E+00_wp, -0.198481745303595E+00_wp, &
+ & -4.912396810694304E-12_wp, -0.217376911834996E+00_wp, 0.188386308225154E+00_wp, &
+ & -2.713809327593799E-13_wp, -2.405883307999757E-02_wp, 1.009543707641887E-02_wp, &
+ & -4.821379744720639E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 4.023966094020255E-06_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 4.418983135319131E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
& 0.000000000000000E+00_wp, -0.327069270111905E+00_wp, -0.193419606975667E+00_wp, &
& 0.000000000000000E+00_wp, -2.159771757181671E-02_wp, -2.400913565342059E-02_wp, &
- & 0.000000000000000E+00_wp, 0.348666987689540E+00_wp, 0.217428742639103E+00_wp, &
+ & 0.000000000000000E+00_wp, 0.348666987689540E+00_wp, 0.217428742639103E+00_wp, &
& 0.000000000000000E+00_wp, -0.241435744067920E+00_wp, -0.198481744759780E+00_wp, &
- & 0.000000000000000E+00_wp, 2.405883312826490E-02_wp, 1.009543692497955E-02_wp, &
- & 0.000000000000000E+00_wp, 0.217376910941702E+00_wp, 0.188386307847977E+00_wp],&
+ & 0.000000000000000E+00_wp, 2.405883312826490E-02_wp, 1.009543692497955E-02_wp, &
+ & 0.000000000000000E+00_wp, 0.217376910941702E+00_wp, 0.188386307847977E+00_wp], &
& shape(hessian_ref))
real(wp), parameter :: step = 1.0e-6_wp
@@ -112,16 +167,16 @@ subroutine test_gfn1_hessian(error)
type(scc_results) :: res
type(TxTBCalculator) :: calc
- integer :: i,j
+ integer :: i, j
real(wp) :: energy, sigma(3, 3)
real(wp) :: hl_gap
- real(wp),allocatable :: gradient(:,:), dipgrad(:,:), hessian(:,:)
+ real(wp), allocatable :: gradient(:, :), dipgrad(:, :), hessian(:, :)
integer, allocatable :: list(:)
call init(env)
call init(mol, sym, xyz)
- allocate(gradient(3,mol%n), dipgrad(3,3*mol%n), hessian(3*mol%n,3*mol%n))
+ allocate(gradient(3, mol%n), dipgrad(3, 3*mol%n), hessian(3*mol%n, 3*mol%n))
energy = 0.0_wp
gradient = 0.0_wp
@@ -139,66 +194,72 @@ subroutine test_gfn1_hessian(error)
do i = 1, size(dipgrad_ref, 2)
do j = 1, size(dipgrad_ref, 1)
call check(error, dipgrad(j, i), dipgrad_ref(j, i), thr=thr)
+ if (allocated(error)) return
end do
end do
do i = 1, size(hessian_ref, 2)
do j = 1, size(hessian_ref, 1)
call check(error, hessian(j, i), hessian_ref(j, i), thr=thr)
+ if (allocated(error)) return
end do
end do
end subroutine test_gfn1_hessian
+!> Pins the GFN2 numerical Hessian and dipole gradient of water against the
+!> reference tables below.
subroutine test_gfn2_hessian(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
+
integer, parameter :: nat = 3
- real(wp),parameter :: thr = 1.0e-7_wp
+ real(wp), parameter :: thr = 1.0e-7_wp
character(len=*), parameter :: sym(nat) = ["O", "H", "H"]
real(wp), parameter :: xyz(3, nat) = reshape([&
- & 0.00000000000000_wp, -0.00000000077760_wp, 0.18829790750029_wp, &
- & 0.00000000000000_wp, 1.45987612440076_wp, -0.88615189669760_wp, &
- &-0.00000000000000_wp, -1.45987612362316_wp, -0.88615189608629_wp],&
+ & 0.00000000000000_wp, -0.00000000077760_wp, 0.18829790750029_wp, &
+ & 0.00000000000000_wp, 1.45987612440076_wp, -0.88615189669760_wp, &
+ &-0.00000000000000_wp, -1.45987612362316_wp, -0.88615189608629_wp], &
& shape(xyz))
real(wp), parameter :: dipgrad_ref(3, 3*nat) = reshape([ &
- & -0.811313733750829E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & -0.811313733750829E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
& -3.379448427312665E-11_wp, -0.412116147045635E+00_wp, -2.220446049250313E-09_wp, &
- & -1.688661825390276E-11_wp, 1.035790186933606E-09_wp, -0.483878709767183E+00_wp, &
- & 0.405648009573818E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 7.146633102892883E-11_wp, 0.205575204598719E+00_wp, 9.503822501200077E-02_wp, &
- & -1.829769461273858E-14_wp, 0.148069825228005E+00_wp, 0.243425242296702E+00_wp, &
- & 0.405648009946286E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & -4.559890805138046E-14_wp, 0.205575207063050E+00_wp, -9.503822462342271E-02_wp, &
- & 3.356458019834560E-11_wp, -0.148069829195679E+00_wp, 0.243425244239592E+00_wp],&
+ & -1.688661825390276E-11_wp, 1.035790186933606E-09_wp, -0.483878709767183E+00_wp, &
+ & 0.405648009573818E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 7.146633102892883E-11_wp, 0.205575204598719E+00_wp, 9.503822501200077E-02_wp, &
+ & -1.829769461273858E-14_wp, 0.148069825228005E+00_wp, 0.243425242296702E+00_wp, &
+ & 0.405648009946286E+00_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & -4.559890805138046E-14_wp, 0.205575207063050E+00_wp, -9.503822462342271E-02_wp, &
+ & 3.356458019834560E-11_wp, -0.148069829195679E+00_wp, 0.243425244239592E+00_wp], &
& shape(dipgrad_ref))
real(wp), parameter :: hessian_ref(3*nat, 3*nat) = reshape([ &
- & -1.939596096290860E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 9.697686009794110E-06_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 9.698274953126006E-06_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & -1.978173339857878E-11_wp, 0.612198809002427E+00_wp, 2.697691516787698E-09_wp, &
- & 2.481358892927353E-11_wp, -0.306099403183025E+00_wp, 0.225852148702132E+00_wp, &
+ & -1.939596096290860E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 9.697686009794110E-06_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 9.698274953126006E-06_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & -1.978173339857878E-11_wp, 0.612198809002427E+00_wp, 2.697691516787698E-09_wp, &
+ & 2.481358892927353E-11_wp, -0.306099403183025E+00_wp, 0.225852148702132E+00_wp, &
& -5.031855530694753E-12_wp, -0.306099405820178E+00_wp, -0.225852151393925E+00_wp, &
- & -1.024802643153815E-11_wp, 2.163320452807998E-09_wp, 0.399810948644134E+00_wp, &
- & 1.255439924419533E-11_wp, 0.184154284413794E+00_wp, -0.199905473616475E+00_wp, &
+ & -1.024802643153815E-11_wp, 2.163320452807998E-09_wp, 0.399810948644134E+00_wp, &
+ & 1.255439924419533E-11_wp, 0.184154284413794E+00_wp, -0.199905473616475E+00_wp, &
& -2.306372812657182E-12_wp, -0.184154286575537E+00_wp, -0.199905475030756E+00_wp, &
- & -6.559552217481455E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 5.182604735260440E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 1.376947482221124E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 3.321776245993920E-11_wp, -0.305714839568990E+00_wp, 0.183399699337614E+00_wp, &
- & -3.597681672480956E-11_wp, 0.339222824678649E+00_wp, -0.204443891042213E+00_wp, &
- & 2.759054264870358E-12_wp, -3.350798510798255E-02_wp, 2.104419170521866E-02_wp, &
- & -2.747668766007676E-13_wp, 0.225612945682292E+00_wp, -0.199989952327133E+00_wp, &
- & 2.050024431283972E-13_wp, -0.204892631012658E+00_wp, 0.183497773686177E+00_wp, &
- & 6.976443347237028E-14_wp, -2.072031467084191E-02_wp, 1.649217862695682E-02_wp, &
- & -6.559491291679194E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 1.376944591474809E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
- & 5.182546700204432E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & -6.559552217481455E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 5.182604735260440E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 1.376947482221124E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 3.321776245993920E-11_wp, -0.305714839568990E+00_wp, 0.183399699337614E+00_wp, &
+ & -3.597681672480956E-11_wp, 0.339222824678649E+00_wp, -0.204443891042213E+00_wp, &
+ & 2.759054264870358E-12_wp, -3.350798510798255E-02_wp, 2.104419170521866E-02_wp, &
+ & -2.747668766007676E-13_wp, 0.225612945682292E+00_wp, -0.199989952327133E+00_wp, &
+ & 2.050024431283972E-13_wp, -0.204892631012658E+00_wp, 0.183497773686177E+00_wp, &
+ & 6.976443347237028E-14_wp, -2.072031467084191E-02_wp, 1.649217862695682E-02_wp, &
+ & -6.559491291679194E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 1.376944591474809E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
+ & 5.182546700204432E-05_wp, 0.000000000000000E+00_wp, 0.000000000000000E+00_wp, &
& -2.703865614441994E-13_wp, -0.305714842513751E+00_wp, -0.183399702402830E+00_wp, &
& 2.719647608323040E-13_wp, -3.350798515432542E-02_wp, -2.104419150402123E-02_wp, &
- & -1.578199388104547E-15_wp, 0.339222827661363E+00_wp, 0.204443893905696E+00_wp, &
+ & -1.578199388104547E-15_wp, 0.339222827661363E+00_wp, 0.204443893905696E+00_wp, &
& 2.029210950346453E-11_wp, -0.225612948054816E+00_wp, -0.199989953443651E+00_wp, &
- & -2.499966571670633E-11_wp, 2.072031476318205E-02_wp, 1.649217857736473E-02_wp, &
- & 4.707556213241799E-12_wp, 0.204892633286792E+00_wp, 0.183497774863884E+00_wp],&
+ & -2.499966571670633E-11_wp, 2.072031476318205E-02_wp, 1.649217857736473E-02_wp, &
+ & 4.707556213241799E-12_wp, 0.204892633286792E+00_wp, 0.183497774863884E+00_wp], &
& shape(hessian_ref))
real(wp), parameter :: step = 1.0e-6_wp
@@ -209,16 +270,16 @@ subroutine test_gfn2_hessian(error)
type(scc_results) :: res
type(TxTBCalculator) :: calc
- integer :: i,j
+ integer :: i, j
real(wp) :: energy, sigma(3, 3)
real(wp) :: hl_gap
- real(wp),allocatable :: gradient(:,:), dipgrad(:,:), hessian(:,:)
+ real(wp), allocatable :: gradient(:, :), dipgrad(:, :), hessian(:, :)
integer, allocatable :: list(:)
call init(env)
call init(mol, sym, xyz)
- allocate(gradient(3,mol%n), dipgrad(3,3*mol%n), hessian(3*mol%n,3*mol%n))
+ allocate(gradient(3, mol%n), dipgrad(3, 3*mol%n), hessian(3*mol%n, 3*mol%n))
energy = 0.0_wp
gradient = 0.0_wp
@@ -236,55 +297,61 @@ subroutine test_gfn2_hessian(error)
do i = 1, size(dipgrad_ref, 2)
do j = 1, size(dipgrad_ref, 1)
call check(error, dipgrad(j, i), dipgrad_ref(j, i), thr=thr)
+ if (allocated(error)) return
end do
end do
do i = 1, size(hessian_ref, 2)
do j = 1, size(hessian_ref, 1)
call check(error, hessian(j, i), hessian_ref(j, i), thr=thr)
+ if (allocated(error)) return
end do
end do
end subroutine test_gfn2_hessian
+!> Pins the GFN1 O1NumHess (ODLR) Hessian of water and its final error
+!> estimate against the reference values below.
subroutine test_o1numhess_gfn1(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
+
integer, parameter :: nat = 3
- real(wp),parameter :: thr1 = 1.4901161193847656e-08_wp, thr2 = 1.4901161193847656e-08_wp
+ real(wp), parameter :: thr1 = 1.4901161193847656e-08_wp, thr2 = 1.4901161193847656e-08_wp
character(len=*), parameter :: sym(nat) = ["O", "H", "H"]
real(wp), parameter :: xyz(3, nat) = reshape([&
- & 0.00000000000000_wp, 0.00000000034546_wp, 0.18900383618455_wp, &
- & 0.00000000000000_wp, 1.45674735348811_wp, -0.88650486059828_wp, &
- &-0.00000000000000_wp, -1.45674735383357_wp, -0.88650486086986_wp],&
+ & 0.00000000000000_wp, 0.00000000034546_wp, 0.18900383618455_wp, &
+ & 0.00000000000000_wp, 1.45674735348811_wp, -0.88650486059828_wp, &
+ &-0.00000000000000_wp, -1.45674735383357_wp, -0.88650486086986_wp], &
& shape(xyz))
real(wp), parameter :: hessian_ref(9, 9) = reshape([&
- & 9.642722735617289E-05_wp, 2.150578923697825E-12_wp, 1.071019141435998E-12_wp, &
- &-4.821347480851584E-05_wp, -4.062528712426575E-12_wp, 2.583698911064602E-13_wp, &
- &-4.821375254766774E-05_wp, 1.911949788728749E-12_wp, -1.329389032542458E-12_wp, &
- & 2.150578923697825E-12_wp, 6.543135697355346E-01_wp, -7.229475687954884E-07_wp, &
- &-1.642628093963891E-12_wp, -3.271571836987667E-01_wp, 2.415028515059434E-01_wp, &
+ & 9.642722735617289E-05_wp, 2.150578923697825E-12_wp, 1.071019141435998E-12_wp, &
+ &-4.821347480851584E-05_wp, -4.062528712426575E-12_wp, 2.583698911064602E-13_wp, &
+ &-4.821375254766774E-05_wp, 1.911949788728749E-12_wp, -1.329389032542458E-12_wp, &
+ & 2.150578923697825E-12_wp, 6.543135697355346E-01_wp, -7.229475687954884E-07_wp, &
+ &-1.642628093963891E-12_wp, -3.271571836987667E-01_wp, 2.415028515059434E-01_wp, &
&-5.079508297339344E-13_wp, -3.271563860451477E-01_wp, -2.415021285496890E-01_wp, &
- & 1.071019141435998E-12_wp, -7.229475687954884E-07_wp, 3.970231512770018E-01_wp, &
- &-6.195795301597118E-13_wp, 1.934311885548390E-01_wp, -1.985118423833666E-01_wp, &
+ & 1.071019141435998E-12_wp, -7.229475687954884E-07_wp, 3.970231512770018E-01_wp, &
+ &-6.195795301597118E-13_wp, 1.934311885548390E-01_wp, -1.985118423833666E-01_wp, &
&-4.514396112762865E-13_wp, -1.934304656060189E-01_wp, -1.985113088943410E-01_wp, &
&-4.821347480851584E-05_wp, -1.642628093963891E-12_wp, -6.195795301597118E-13_wp, &
- & 4.418955528562922E-05_wp, 3.506229450266824E-12_wp, -2.965819713964842E-13_wp, &
- & 4.023919522892954E-06_wp, -1.863601356302934E-12_wp, 9.161615015561960E-13_wp, &
- &-4.062528712426575E-12_wp, -3.271571836987667E-01_wp, 1.934311885548390E-01_wp, &
- & 3.506229450266824E-12_wp, 3.487436663442405E-01_wp, -2.174669866056496E-01_wp, &
- & 5.562992621597500E-13_wp, -2.158648263991621E-02_wp, 2.403579804767703E-02_wp, &
- & 2.583698911064602E-13_wp, 2.415028515059434E-01_wp, -1.985118423833666E-01_wp, &
- &-2.965819713964842E-13_wp, -2.174669866056496E-01_wp, 1.884261289483365E-01_wp, &
- & 3.821208029002418E-14_wp, -2.403586490401287E-02_wp, 1.008571343858931E-02_wp, &
+ & 4.418955528562922E-05_wp, 3.506229450266824E-12_wp, -2.965819713964842E-13_wp, &
+ & 4.023919522892954E-06_wp, -1.863601356302934E-12_wp, 9.161615015561960E-13_wp, &
+ &-4.062528712426575E-12_wp, -3.271571836987667E-01_wp, 1.934311885548390E-01_wp, &
+ & 3.506229450266824E-12_wp, 3.487436663442405E-01_wp, -2.174669866056496E-01_wp, &
+ & 5.562992621597500E-13_wp, -2.158648263991621E-02_wp, 2.403579804767703E-02_wp, &
+ & 2.583698911064602E-13_wp, 2.415028515059434E-01_wp, -1.985118423833666E-01_wp, &
+ &-2.965819713964842E-13_wp, -2.174669866056496E-01_wp, 1.884261289483365E-01_wp, &
+ & 3.821208029002418E-14_wp, -2.403586490401287E-02_wp, 1.008571343858931E-02_wp, &
&-4.821375254766774E-05_wp, -5.079508297339344E-13_wp, -4.514396112762865E-13_wp, &
- & 4.023919522892954E-06_wp, 5.562992621597500E-13_wp, 3.821208029002418E-14_wp, &
- & 4.418983302477912E-05_wp, -4.834843242581555E-14_wp, 4.132275309862622E-13_wp, &
+ & 4.023919522892954E-06_wp, 5.562992621597500E-13_wp, 3.821208029002418E-14_wp, &
+ & 4.418983302477912E-05_wp, -4.834843242581555E-14_wp, 4.132275309862622E-13_wp, &
& 1.911949788728749E-12_wp, -3.271563860451477E-01_wp, -1.934304656060189E-01_wp, &
&-1.863601356302934E-12_wp, -2.158648263991621E-02_wp, -2.403586490401287E-02_wp, &
- &-4.834843242581555E-14_wp, 3.487428686878863E-01_wp, 2.174663305044797E-01_wp, &
+ &-4.834843242581555E-14_wp, 3.487428686878863E-01_wp, 2.174663305044797E-01_wp, &
&-1.329389032542458E-12_wp, -2.415021285496890E-01_wp, -1.985113088943410E-01_wp, &
- & 9.161615015561960E-13_wp, 2.403579804767703E-02_wp, 1.008571343858931E-02_wp, &
- & 4.132275309862622E-13_wp, 2.174663305044797E-01_wp, 1.884255954528983E-01_wp],&
+ & 9.161615015561960E-13_wp, 2.403579804767703E-02_wp, 1.008571343858931E-02_wp, &
+ & 4.132275309862622E-13_wp, 2.174663305044797E-01_wp, 1.884255954528983E-01_wp], &
& shape(hessian_ref))
real(wp), parameter :: final_err_ref = 1.53703515069416350E-06_wp
real(wp), parameter :: step = 1.0e-6_wp
@@ -304,7 +371,7 @@ subroutine test_o1numhess_gfn1(error)
call newXTBCalculator(env, mol, calc, method=1, accuracy=1.0e-8_wp)
call newWavefunction(env, mol, calc, chk)
- allocate(hessian(N, N), dipgrad_dummy(3, N))
+ allocate(hessian(N, N), dipgrad_dummy(3, N), source=0.0_wp)
call calc%hessian(env, mol, chk, [(i, i=1, mol%n)], step, hessian, dipgrad_dummy, odlr=.true., final_err=final_err)
if (any(abs(hessian - hessian_ref) > thr1)) then
@@ -312,62 +379,67 @@ subroutine test_o1numhess_gfn1(error)
print *, "--- hessian ---"
do i = 1, N
- print '(*(F21.14))', hessian(i, :)
+ print "(*(F21.14))", hessian(i, :)
end do
print *, "--- Ref. hessian ---"
do i = 1, N
- print '(*(F21.14))', hessian_ref(i, :)
+ print "(*(F21.14))", hessian_ref(i, :)
end do
+ return
end if
- if (abs(final_err - final_err_ref) > thr2) then
- call test_failed(error, "Final error does not match")
- print *, "--- Final error ---"
- print '(*(F21.14))', final_err
- print *, "--- Ref. final error ---"
- print '(*(F21.14))', final_err_ref
- end if
- end subroutine test_o1numhess_gfn1
+ if (abs(final_err - final_err_ref) > thr2) then
+ call test_failed(error, "Final error does not match")
+ print *, "--- Final error ---"
+ print "(*(F21.14))", final_err
+ print *, "--- Ref. final error ---"
+ print "(*(F21.14))", final_err_ref
+ end if
+end subroutine test_o1numhess_gfn1
+!> Pins the GFN2 O1NumHess (ODLR) Hessian of water and its final error
+!> estimate against the reference values below.
subroutine test_o1numhess_gfn2(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
+
integer, parameter :: nat = 3
- real(wp),parameter :: thr1 = 1.4901161193847656e-08_wp, thr2 = 1.4901161193847656e-08_wp
+ real(wp), parameter :: thr1 = 1.4901161193847656e-08_wp, thr2 = 1.4901161193847656e-08_wp
character(len=*), parameter :: sym(nat) = ["O", "H", "H"]
real(wp), parameter :: xyz(3, nat) = reshape([&
- & 0.00000000000000_wp, 0.00000000034546_wp, 0.18900383618455_wp, &
- & 0.00000000000000_wp, 1.45674735348811_wp, -0.88650486059828_wp, &
- &-0.00000000000000_wp, -1.45674735383357_wp, -0.88650486086986_wp],&
+ & 0.00000000000000_wp, 0.00000000034546_wp, 0.18900383618455_wp, &
+ & 0.00000000000000_wp, 1.45674735348811_wp, -0.88650486059828_wp, &
+ &-0.00000000000000_wp, -1.45674735383357_wp, -0.88650486086986_wp], &
& shape(xyz))
real(wp), parameter :: hessian_ref(9, 9) = reshape([&
- &-7.853049280457453E-04_wp, 1.955977244785751E-12_wp, 1.247080264997499E-12_wp, &
- & 3.926525943422639E-04_wp, -1.210600324495050E-12_wp, 9.850364921710785E-14_wp, &
+ &-7.853049280457453E-04_wp, 1.955977244785751E-12_wp, 1.247080264997499E-12_wp, &
+ & 3.926525943422639E-04_wp, -1.210600324495050E-12_wp, 9.850364921710785E-14_wp, &
& 3.926523337034877E-04_wp, -7.453769202907016E-13_wp, -1.345583914214606E-12_wp, &
- & 1.955977244785751E-12_wp, 6.131511844074383E-01_wp, -7.530713906089167E-07_wp, &
- &-4.508671918596937E-13_wp, -3.065760101726123E-01_wp, 2.266333722977527E-01_wp, &
+ & 1.955977244785751E-12_wp, 6.131511844074383E-01_wp, -7.530713906089167E-07_wp, &
+ &-4.508671918596937E-13_wp, -3.065760101726123E-01_wp, 2.266333722977527E-01_wp, &
&-1.505110052926058E-12_wp, -3.065751742284646E-01_wp, -2.266326192307660E-01_wp, &
- & 1.247080264997499E-12_wp, -7.530713906089167E-07_wp, 4.009666360184508E-01_wp, &
- &-8.063227290784491E-13_wp, 1.851054221071441E-01_wp, -2.004835959201852E-01_wp, &
+ & 1.247080264997499E-12_wp, -7.530713906089167E-07_wp, 4.009666360184508E-01_wp, &
+ &-8.063227290784491E-13_wp, 1.851054221071441E-01_wp, -2.004835959201852E-01_wp, &
&-4.407575359190495E-13_wp, -1.851046690338150E-01_wp, -2.004830401052821E-01_wp, &
& 3.926525943422639E-04_wp, -4.508671918596937E-13_wp, -8.063227290784491E-13_wp, &
- &-5.034357727917302E-04_wp, 9.315986645004943E-13_wp, 2.367249284567081E-13_wp, &
- & 1.107831784494622E-04_wp, -4.807314726408010E-13_wp, 5.695978006217410E-13_wp, &
- &-1.210600324495050E-12_wp, -3.065760101726123E-01_wp, 1.851054221071441E-01_wp, &
- & 9.315986645004943E-13_wp, 3.407104032218639E-01_wp, -2.058693632980307E-01_wp, &
- & 2.790016599945554E-13_wp, -3.413439305260558E-02_wp, 2.076394119763031E-02_wp, &
- & 9.850364921710785E-14_wp, 2.266333722977527E-01_wp, -2.004835959201852E-01_wp, &
- & 2.367249284567081E-13_wp, -2.058693632980307E-01_wp, 1.835957853733738E-01_wp, &
- &-3.352285776738158E-13_wp, -2.076400899834286E-02_wp, 1.688781054869399E-02_wp, &
+ &-5.034357727917302E-04_wp, 9.315986645004943E-13_wp, 2.367249284567081E-13_wp, &
+ & 1.107831784494622E-04_wp, -4.807314726408010E-13_wp, 5.695978006217410E-13_wp, &
+ &-1.210600324495050E-12_wp, -3.065760101726123E-01_wp, 1.851054221071441E-01_wp, &
+ & 9.315986645004943E-13_wp, 3.407104032218639E-01_wp, -2.058693632980307E-01_wp, &
+ & 2.790016599945554E-13_wp, -3.413439305260558E-02_wp, 2.076394119763031E-02_wp, &
+ & 9.850364921710785E-14_wp, 2.266333722977527E-01_wp, -2.004835959201852E-01_wp, &
+ & 2.367249284567081E-13_wp, -2.058693632980307E-01_wp, 1.835957853733738E-01_wp, &
+ &-3.352285776738158E-13_wp, -2.076400899834286E-02_wp, 1.688781054869399E-02_wp, &
& 3.926523337034877E-04_wp, -1.505110052926058E-12_wp, -4.407575359190495E-13_wp, &
- & 1.107831784494622E-04_wp, 2.790016599945554E-13_wp, -3.352285776738158E-13_wp, &
- &-5.034355121529521E-04_wp, 1.226108392931502E-12_wp, 7.759861135928650E-13_wp, &
+ & 1.107831784494622E-04_wp, 2.790016599945554E-13_wp, -3.352285776738158E-13_wp, &
+ &-5.034355121529521E-04_wp, 1.226108392931502E-12_wp, 7.759861135928650E-13_wp, &
&-7.453769202907016E-13_wp, -3.065751742284646E-01_wp, -1.851046690338150E-01_wp, &
&-4.807314726408010E-13_wp, -3.413439305260558E-02_wp, -2.076400899834286E-02_wp, &
- & 1.226108392931502E-12_wp, 3.407095672780631E-01_wp, 2.058686780298182E-01_wp, &
+ & 1.226108392931502E-12_wp, 3.407095672780631E-01_wp, 2.058686780298182E-01_wp, &
&-1.345583914214606E-12_wp, -2.266326192307660E-01_wp, -2.004830401052821E-01_wp, &
- & 5.695978006217410E-13_wp, 2.076394119763031E-02_wp, 1.688781054869399E-02_wp, &
- & 7.759861135928650E-13_wp, 2.058686780298182E-01_wp, 1.835952295617220E-01_wp],&
+ & 5.695978006217410E-13_wp, 2.076394119763031E-02_wp, 1.688781054869399E-02_wp, &
+ & 7.759861135928650E-13_wp, 2.058686780298182E-01_wp, 1.835952295617220E-01_wp], &
& shape(hessian_ref))
real(wp), parameter :: final_err_ref = 1.27700129644838372E-06_wp
real(wp), parameter :: step = 1.0e-6_wp
@@ -387,7 +459,7 @@ subroutine test_o1numhess_gfn2(error)
call newXTBCalculator(env, mol, calc, method=2, accuracy=1.0e-8_wp)
call newWavefunction(env, mol, calc, chk)
- allocate(hessian(N, N), dipgrad_dummy(3, N))
+ allocate(hessian(N, N), dipgrad_dummy(3, N), source=0.0_wp)
call calc%hessian(env, mol, chk, [(i, i=1, mol%n)], step, hessian, dipgrad_dummy, odlr=.true., final_err=final_err)
if (any(abs(hessian - hessian_ref) > thr1)) then
@@ -395,33 +467,38 @@ subroutine test_o1numhess_gfn2(error)
print *, "--- hessian ---"
do i = 1, N
- print '(*(F21.14))', hessian(i, :)
+ print "(*(F21.14))", hessian(i, :)
end do
print *, "--- Ref. hessian ---"
do i = 1, N
- print '(*(F21.14))', hessian_ref(i, :)
+ print "(*(F21.14))", hessian_ref(i, :)
end do
+ return
end if
if (abs(final_err - final_err_ref) > thr2) then
call test_failed(error, "Final error does not match")
print *, "--- Final error ---"
- print '(*(F21.14))', final_err
+ print "(*(F21.14))", final_err
print *, "--- Ref. final error ---"
- print '(*(F21.14))', final_err_ref
+ print "(*(F21.14))", final_err_ref
end if
end subroutine test_o1numhess_gfn2
+!> Pins the GFN1 O1NumHess Hessian of linear water through its eigenvalue
+!> spectrum: three near-zero and two negative frequencies around -0.16.
subroutine test_o1numhess_linear_h2o_gfn1(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
+
integer, parameter :: nat = 3
- real(wp),parameter :: thr1 = 1.0e-9_wp, thr2 = 1.0e-5_wp
+ real(wp), parameter :: thr1 = 1.0e-9_wp, thr2 = 1.0e-5_wp
character(len=*), parameter :: sym(nat) = ["O", "H", "H"]
real(wp), parameter :: xyz(3, nat) = reshape([&
- & 0.00000000000000_wp, 0.00000000000000_wp, 0.00000000000000_wp, &
- & 0.00000000000000_wp, 0.00000000000000_wp, -1.81075448577205_wp, &
- & 0.00000000000000_wp, 0.00000000000000_wp, 1.81075448676713_wp],&
+ & 0.00000000000000_wp, 0.00000000000000_wp, 0.00000000000000_wp, &
+ & 0.00000000000000_wp, 0.00000000000000_wp, -1.81075448577205_wp, &
+ & 0.00000000000000_wp, 0.00000000000000_wp, 1.81075448676713_wp], &
& shape(xyz))
real(wp), parameter :: step = 1.0e-6_wp
@@ -440,36 +517,45 @@ subroutine test_o1numhess_linear_h2o_gfn1(error)
call newXTBCalculator(env, mol, calc, method=1, accuracy=1.0e-8_wp)
call newWavefunction(env, mol, calc, chk)
- allocate(hessian(N, N), dipgrad_dummy(3, N))
+ allocate(hessian(N, N), dipgrad_dummy(3, N), source=0.0_wp)
call calc%hessian(env, mol, chk, [(i, i=1, mol%n)], step, hessian, dipgrad_dummy, odlr=.true., final_err=final_err)
allocate(freq(N))
- lwork = 1 + 6*N + 2*N**2
+ lwork = 1 + 6 * N + 2 * N**2
allocate(aux(lwork))
- call dsyev ('V', 'U', N, hessian, N, freq, aux, lwork, info)
+ call dsyev("V", "U", N, hessian, N, freq, aux, lwork, info)
+ call check(error, info, 0)
+ if (allocated(error)) return
if (count(abs(freq) < 1.0e-4_wp) /= 3) then
call test_failed(error, "Linear H2O should have exactly three ~0 freqs")
+ return
end if
if (count(freq < -1.0e-10_wp) /= 2) then
call test_failed(error, "Linear H2O should have exactly two negative freqs")
+ return
end if
! freqs around -0.16
if (freq(1) > -0.16_wp .or. freq(2) > -0.16_wp) then
call test_failed(error, "First two freqs should be negative")
print *, freq
+ return
end if
end subroutine test_o1numhess_linear_h2o_gfn1
+!> Pins the GFN2 O1NumHess Hessian of linear water through its eigenvalue
+!> spectrum: three near-zero and two negative frequencies around -0.3.
subroutine test_o1numhess_linear_h2o_gfn2(error)
+ !> Failure report, allocated when the check fails.
type(error_type), allocatable, intent(out) :: error
+
integer, parameter :: nat = 3
- real(wp),parameter :: thr1 = 1.0e-9_wp, thr2 = 1.0e-5_wp
+ real(wp), parameter :: thr1 = 1.0e-9_wp, thr2 = 1.0e-5_wp
character(len=*), parameter :: sym(nat) = ["O", "H", "H"]
real(wp), parameter :: xyz(3, nat) = reshape([&
- & 0.00000000000000_wp, 0.00000000000000_wp, 0.00000000000000_wp, &
- & 0.00000000000000_wp, 0.00000000000000_wp, -1.81075448577205_wp, &
- & 0.00000000000000_wp, 0.00000000000000_wp, 1.81075448676713_wp],&
+ & 0.00000000000000_wp, 0.00000000000000_wp, 0.00000000000000_wp, &
+ & 0.00000000000000_wp, 0.00000000000000_wp, -1.81075448577205_wp, &
+ & 0.00000000000000_wp, 0.00000000000000_wp, 1.81075448676713_wp], &
& shape(xyz))
real(wp), parameter :: step = 1.0e-6_wp
@@ -488,25 +574,565 @@ subroutine test_o1numhess_linear_h2o_gfn2(error)
call newXTBCalculator(env, mol, calc, method=2, accuracy=1.0e-8_wp)
call newWavefunction(env, mol, calc, chk)
- allocate(hessian(N, N), dipgrad_dummy(3, N))
+ allocate(hessian(N, N), dipgrad_dummy(3, N), source=0.0_wp)
call calc%hessian(env, mol, chk, [(i, i=1, mol%n)], step, hessian, dipgrad_dummy, odlr=.true., final_err=final_err)
allocate(freq(N))
- lwork = 1 + 6*N + 2*N**2
+ lwork = 1 + 6 * N + 2 * N**2
allocate(aux(lwork))
- call dsyev ('V', 'U', N, hessian, N, freq, aux, lwork, info)
+ call dsyev("V", "U", N, hessian, N, freq, aux, lwork, info)
+ call check(error, info, 0)
+ if (allocated(error)) return
if (count(abs(freq) < 1.0e-4_wp) /= 3) then
call test_failed(error, "Linear H2O should have exactly three ~0 freqs")
+ return
end if
if (count(freq < -1.0e-10_wp) /= 2) then
call test_failed(error, "Linear H2O should have exactly two negative freqs")
+ return
end if
! freqs around -0.3
if (freq(1) > -0.3_wp .or. freq(2) > -0.3_wp) then
call test_failed(error, "First two freqs should be negative")
+ return
end if
end subroutine test_o1numhess_linear_h2o_gfn2
+!> Pins the compliance constants of a two-atom harmonic model against their
+!> analytic values, including the redundant B-matrix with duplicated rows.
+subroutine test_compliance(error)
+ !> Failure report, allocated when the check fails.
+ type(error_type), allocatable, intent(out) :: error
+
+ integer, parameter :: nat = 2
+ integer, parameter :: ndim = 3 * nat
+ real(wp), parameter :: force_constant = 2.0_wp
+ real(wp), parameter :: thr = 1.0e-10_wp
+ integer, parameter :: pairs(2, 1) = reshape([1, 2], [2, 1])
+ real(wp), parameter :: xyz(3, nat) = reshape([ &
+ & 0.0_wp, 0.0_wp, 0.0_wp, &
+ & 2.0_wp, 0.0_wp, 0.0_wp], shape(xyz))
+
+ type(TEnvironment) :: env
+ type(TNeighbourList) :: neigh_list
+ type(graph_type) :: graph
+ type(redundant_type) :: internals
+ integer :: i, j
+ logical :: failed
+ real(wp), allocatable :: bmat(:, :), compliance(:, :), hessian(:, :)
+ real(wp) :: redundant_bmat(2, ndim), redundant_compliance(2, 2)
+
+ call init(env)
+ env%unit = 0
+ neigh_list = build_neigh_list(nat, pairs)
+ call init(graph, neigh_list)
+ call init(internals, graph, xyz)
+ call check(error, internals%ncoords, 1)
+ if (allocated(error)) return
+
+ allocate(bmat(internals%ncoords, ndim), compliance(internals%ncoords, internals%ncoords), &
+ & hessian(ndim, ndim))
+ call get_bmatrix(internals, xyz, bmat)
+
+ do i = 1, ndim
+ do j = 1, ndim
+ hessian(i, j) = force_constant * bmat(1, i) * bmat(1, j)
+ end do
+ end do
+
+ call compute_compliance(env, hessian, bmat, xyz, nat, internals%ncoords, compliance)
+ call env%check(failed)
+ if (failed) then
+ call test_failed(error, "Compliance calculation failed")
+ return
+ end if
+ call check(error, bmat(1, 1), -1.0_wp, thr=thr)
+ if (allocated(error)) return
+ call check(error, bmat(1, 4), 1.0_wp, thr=thr)
+ if (allocated(error)) return
+ call check(error, compliance(1, 1), 1.0_wp/force_constant, thr=thr)
+ if (allocated(error)) return
+
+ redundant_bmat(1, :) = bmat(1, :)
+ redundant_bmat(2, :) = bmat(1, :)
+ call compute_compliance(env, hessian, redundant_bmat, xyz, nat, 2, redundant_compliance)
+ call env%check(failed)
+ if (failed) then
+ call test_failed(error, "Compliance calculation failed")
+ return
+ end if
+ call check(error, redundant_compliance(1, 1), 1.0_wp/force_constant, thr=thr)
+ if (allocated(error)) return
+ call check(error, redundant_compliance(1, 2), 1.0_wp/force_constant, thr=thr)
+ if (allocated(error)) return
+ call check(error, redundant_compliance(2, 1), 1.0_wp/force_constant, thr=thr)
+ if (allocated(error)) return
+ call check(error, redundant_compliance(2, 2), 1.0_wp/force_constant, thr=thr)
+
+end subroutine test_compliance
+
+!> Pins the GFN1 compliance constants of water at the default Hessian step
+!> against the reference matrix below.
+subroutine test_compliance_water(error)
+ !> Failure report, allocated when the check fails.
+ type(error_type), allocatable, intent(out) :: error
+
+ integer, parameter :: nat = 3
+ integer, parameter :: ndim = 3 * nat
+ real(wp), parameter :: step = 5.0e-3_wp
+ integer, parameter :: at(nat) = [8, 1, 1]
+ character(len=*), parameter :: sym(nat) = ["O", "H", "H"]
+ real(wp), parameter :: xyz(3, nat) = reshape([&
+ & 0.00000000000000_wp, 0.00000000034546_wp, 0.18900383618455_wp, &
+ & 0.00000000000000_wp, 1.45674735348811_wp, -0.88650486059828_wp, &
+ &-0.00000000000000_wp, -1.45674735383357_wp, -0.88650486086986_wp], &
+ & shape(xyz))
+ real(wp), parameter :: thr = 1.0e-10_wp
+ !> GFN1 compliance constants of water, in (E_h/Bohr^2)^{-1}, in the
+ !> redundant ordering bond(1,2), bond(1,3), angle(2,1,3):
+ !> the two OH bonds are symmetry equivalent, so C(1,1) == C(2,2)
+ real(wp), parameter :: ref(3, 3) = reshape([&
+ & 2.0178418553745141_wp, 3.9503596809836755e-02_wp, -3.8479905242132134e-01_wp, &
+ & 3.9503596809836755e-02_wp, 2.0178418632138921_wp, -3.8479905539794468e-01_wp, &
+ & -3.8479905242132134e-01_wp, -3.8479905539794468e-01_wp, 8.3576726095245473_wp], &
+ & shape(ref))
+
+ type(TMolecule) :: mol
+ type(TRestart) :: chk
+ type(TEnvironment) :: env
+ type(scc_results) :: res
+ type(TxTBCalculator) :: calc
+
+ type(TNeighbourList) :: neigh_list
+ type(graph_type) :: graph
+ type(redundant_type) :: internals
+ integer :: i, j
+ logical :: failed
+ integer, allocatable :: list(:)
+ real(wp) :: energy, sigma(3, 3), hl_gap
+ real(wp), allocatable :: gradient(:, :), dipgrad(:, :), hessian(:, :)
+ real(wp), allocatable :: bmat(:, :), compliance(:, :)
+
+ call init(env)
+ env%unit = 0
+ call init(mol, sym, xyz)
+
+ allocate(gradient(3, nat), dipgrad(3, ndim), hessian(ndim, ndim))
+ energy = 0.0_wp
+ gradient = 0.0_wp
+
+ call newXTBCalculator(env, mol, calc, method=1)
+ call newWavefunction(env, mol, calc, chk)
+ call calc%singlepoint(env, mol, chk, 2, .false., energy, gradient, sigma, &
+ & hl_gap, res)
+
+ dipgrad = 0.0_wp
+ hessian = 0.0_wp
+ list = [(i, i = 1, nat)]
+ call calc%hessian(env, mol, chk, list, step, hessian, dipgrad)
+
+ call init(neigh_list, nat)
+ call neigh_list%generate_covalent(at, xyz)
+ call init(graph, neigh_list)
+ call init(internals, graph, xyz)
+ call check(error, internals%ncoords, 3)
+ if (allocated(error)) return
+
+ allocate(bmat(internals%ncoords, ndim), compliance(internals%ncoords, internals%ncoords))
+ call get_bmatrix(internals, xyz, bmat)
+ call compute_compliance(env, hessian, bmat, xyz, nat, internals%ncoords, compliance)
+ call env%check(failed)
+ if (failed) then
+ call test_failed(error, "Compliance calculation failed")
+ return
+ end if
+
+ do i = 1, internals%ncoords
+ do j = 1, internals%ncoords
+ call check(error, compliance(i, j), ref(i, j), thr=thr)
+ if (allocated(error)) return
+ end do
+ end do
+
+end subroutine test_compliance_water
+
+!> Pins a redundant-coordinate B-matrix against central finite differences.
+subroutine test_redundant_bmatrix_fd(error)
+ !> Failure report, allocated when the check fails.
+ type(error_type), allocatable, intent(out) :: error
+
+ integer, parameter :: nat = 4
+ integer, parameter :: ncoord = 6
+ integer, parameter :: pairs(2, 3) = reshape([1, 2, 2, 3, 3, 4], [2, 3])
+ integer, parameter :: expected_kind(ncoord) = [ &
+ & coord_bond, coord_bond, coord_bond, coord_angle, coord_angle, coord_dihedral]
+ integer, parameter :: expected_atoms(4, ncoord) = reshape([ &
+ & 1, 2, 0, 0, &
+ & 2, 3, 0, 0, &
+ & 3, 4, 0, 0, &
+ & 1, 2, 3, 0, &
+ & 2, 3, 4, 0, &
+ & 1, 2, 3, 4], [4, ncoord])
+ real(wp), parameter :: xyz(3, nat) = reshape([ &
+ & 0.0_wp, 0.0_wp, 0.0_wp, &
+ & 1.5_wp, 0.0_wp, 0.0_wp, &
+ & 2.5_wp, 1.0_wp, 0.0_wp, &
+ & 3.0_wp, 1.0_wp, 1.0_wp], shape(xyz))
+ real(wp), parameter :: step = 1.0e-6_wp
+ real(wp), parameter :: thr = 1.0e-6_wp
+
+ type(TNeighbourList) :: neigh_list
+ type(graph_type) :: graph
+ type(redundant_type) :: internals, right, left
+ real(wp) :: coord(3, nat), dq
+ real(wp), allocatable :: bmat(:, :)
+ integer :: i, ic, ii, k
+
+ neigh_list = build_neigh_list(nat, pairs)
+ call init(graph, neigh_list)
+ call init(internals, graph, xyz)
+ call check(error, internals%ncoords, ncoord)
+ if (allocated(error)) return
+
+ do k = 1, ncoord
+ call check(error, internals%kind(k), expected_kind(k))
+ if (allocated(error)) return
+ do i = 1, 4
+ call check(error, internals%atoms(i, k), expected_atoms(i, k))
+ if (allocated(error)) return
+ end do
+ end do
+
+ allocate(bmat(ncoord, 3*nat))
+ call get_bmatrix(internals, xyz, bmat)
+
+ do i = 1, nat
+ do ic = 1, 3
+ ii = 3 * (i - 1) + ic
+ coord = xyz
+ coord(ic, i) = coord(ic, i) + step
+ call init(right, graph, coord)
+ coord = xyz
+ coord(ic, i) = coord(ic, i) - step
+ call init(left, graph, coord)
+ do k = 1, ncoord
+ dq = right%q(k) - left%q(k)
+ if (internals%kind(k) == coord_dihedral) then
+ dq = modulo(dq + pi, 2.0_wp*pi) - pi
+ end if
+ call check(error, bmat(k, ii), 0.5_wp*dq/step, thr=thr)
+ if (allocated(error)) return
+ end do
+ end do
+ end do
+
+end subroutine test_redundant_bmatrix_fd
+
+!> A straight four-atom chain has three stretches and both transverse bend
+!> directions at each centre, but no defined torsion.
+subroutine test_linear_chain_compliance(error)
+ !> Failure report, allocated when the check fails.
+ type(error_type), allocatable, intent(out) :: error
+
+ integer, parameter :: nat = 4
+ integer, parameter :: ndim = 3 * nat
+ integer, parameter :: ncoord = 7
+ integer, parameter :: pairs(2, 3) = reshape([1, 2, 2, 3, 3, 4], [2, 3])
+ real(wp), parameter :: xyz(3, nat) = reshape([ &
+ & 0.0_wp, 0.0_wp, 0.0_wp, &
+ & 1.0_wp, 0.0_wp, 0.0_wp, &
+ & 2.0_wp, 0.0_wp, 0.0_wp, &
+ & 3.0_wp, 0.0_wp, 0.0_wp], shape(xyz))
+ real(wp), parameter :: thr = 1.0e-9_wp
+
+ type(TEnvironment) :: env
+ type(TNeighbourList) :: neigh_list
+ type(graph_type) :: graph
+ type(redundant_type) :: internals
+ integer :: i, j
+ logical :: failed
+ real(wp) :: expected_bmat(ncoord, ndim), hessian(ndim, ndim)
+ real(wp), allocatable :: bmat(:, :), compliance(:, :)
+
+ call init(env)
+ env%unit = 0
+ neigh_list = build_neigh_list(nat, pairs)
+ call init(graph, neigh_list)
+ call init(internals, graph, xyz)
+ call check(error, internals%ncoords, ncoord)
+ if (allocated(error)) return
+ call check(error, internals%ndihedral, 0)
+ if (allocated(error)) return
+ do i = 1, 3
+ call check(error, internals%kind(i), coord_bond)
+ if (allocated(error)) return
+ end do
+ do i = 4, ncoord
+ call check(error, internals%kind(i), coord_linbend)
+ if (allocated(error)) return
+ end do
+
+ expected_bmat = 0.0_wp
+ do i = 1, 3
+ expected_bmat(i, 3*(i - 1) + 1) = -1.0_wp
+ expected_bmat(i, 3*i + 1) = 1.0_wp
+ end do
+ expected_bmat(4, [3, 6, 9]) = [1.0_wp, -2.0_wp, 1.0_wp]
+ expected_bmat(5, [2, 5, 8]) = [1.0_wp, -2.0_wp, 1.0_wp]
+ expected_bmat(6, [6, 9, 12]) = [1.0_wp, -2.0_wp, 1.0_wp]
+ expected_bmat(7, [5, 8, 11]) = [1.0_wp, -2.0_wp, 1.0_wp]
+ hessian = matmul(transpose(expected_bmat), expected_bmat)
+
+ allocate(bmat(ncoord, ndim), compliance(ncoord, ncoord))
+ call get_bmatrix(internals, xyz, bmat)
+ call compute_compliance(env, hessian, bmat, xyz, nat, ncoord, compliance)
+ call env%check(failed)
+ if (failed) then
+ call test_failed(error, "Compliance calculation failed")
+ return
+ end if
+ do i = 1, ncoord
+ do j = 1, ncoord
+ if (i == j) then
+ call check(error, compliance(i, j), 1.0_wp, thr=thr)
+ else
+ call check(error, compliance(i, j), 0.0_wp, thr=thr)
+ end if
+ if (allocated(error)) return
+ end do
+ end do
+
+end subroutine test_linear_chain_compliance
+
+
+!> Near-linear coordinates use the reference frame even when their B rows are
+!> evaluated at a displaced geometry.
+subroutine test_near_linear_bmatrix_fd(error)
+ !> Failure report, allocated when the check fails.
+ type(error_type), allocatable, intent(out) :: error
+
+ integer, parameter :: nat = 3
+ integer, parameter :: ncoord = 4
+ integer, parameter :: pairs(2, 2) = reshape([1, 2, 2, 3], [2, 2])
+ real(wp), parameter :: step = 1.0e-6_wp
+ real(wp), parameter :: thr = 1.0e-7_wp
+ real(wp), parameter :: xyz(3, nat) = reshape([ &
+ & -1.0_wp, 0.0_wp, 0.0_wp, &
+ & 0.0_wp, 0.0_wp, 0.0_wp, &
+ & cos(pi/180.0_wp), sin(pi/180.0_wp), 0.0_wp], shape(xyz))
+
+ type(TNeighbourList) :: neigh_list
+ type(graph_type) :: graph
+ type(redundant_type) :: internals
+ real(wp) :: coord(3, nat), right(3, nat), left(3, nat)
+ real(wp) :: v1(3), v2(3), qref, qright, qleft
+ real(wp), allocatable :: bmat(:, :)
+ integer :: i, ic, ii, k
+
+ neigh_list = build_neigh_list(nat, pairs)
+ call init(graph, neigh_list)
+ call init(internals, graph, xyz)
+ call check(error, internals%ncoords, ncoord)
+ if (allocated(error)) return
+ do k = 3, 4
+ call check(error, internals%kind(k), coord_linbend)
+ if (allocated(error)) return
+ v1 = xyz(:, 1) - xyz(:, 2)
+ v2 = xyz(:, 3) - xyz(:, 2)
+ qref = dot_product(internals%frame(:, k), v1/norm2(v1) + v2/norm2(v2))
+ call check(error, internals%q(k), qref, thr=thr)
+ if (allocated(error)) return
+ end do
+
+ coord = xyz
+ coord(2, 1) = coord(2, 1) + 0.02_wp
+ allocate(bmat(ncoord, 3*nat))
+ call get_bmatrix(internals, coord, bmat)
+ do i = 1, nat
+ do ic = 1, 3
+ ii = 3 * (i - 1) + ic
+ right = coord
+ right(ic, i) = right(ic, i) + step
+ left = coord
+ left(ic, i) = left(ic, i) - step
+ do k = 3, 4
+ v1 = right(:, 1) - right(:, 2)
+ v2 = right(:, 3) - right(:, 2)
+ qright = dot_product(internals%frame(:, k), v1/norm2(v1) + v2/norm2(v2))
+ v1 = left(:, 1) - left(:, 2)
+ v2 = left(:, 3) - left(:, 2)
+ qleft = dot_product(internals%frame(:, k), v1/norm2(v1) + v2/norm2(v2))
+ call check(error, bmat(k, ii), 0.5_wp*(qright - qleft)/step, thr=thr)
+ if (allocated(error)) return
+ end do
+ end do
+ end do
+
+end subroutine test_near_linear_bmatrix_fd
+
+
+!> The driver writes two bonds and one angle for water, none for separated atoms.
+subroutine test_compliance_driver_redundant(error)
+ type(error_type), allocatable, intent(out) :: error
+
+ type(TEnvironment) :: env
+ type(TMolecule) :: mol
+ real(wp) :: xyz(3, 3), hess(9, 9)
+ integer :: unit, dat_unit, ios, i, icase, counts(3), expected(3)
+ logical :: failed
+ character(256) :: line
+ character(4) :: kind
+
+ xyz(:, 1) = 0.0_wp
+ xyz(:, 2) = [0.0_wp, 1.8_wp, 0.0_wp]
+ xyz(:, 3) = [1.7_wp, -0.5_wp, 0.0_wp]
+ call init(mol, [8, 1, 1], xyz)
+ hess = 0.0_wp
+ do i = 1, 9
+ hess(i, i) = 1.0_wp
+ end do
+ call init(env)
+ open(newunit=unit, status="scratch")
+ env%unit = unit
+
+ expected = [3, 2, 1] ! Total, bond and angle records.
+ do icase = 1, 2
+ if (icase == 2) then
+ mol%xyz = 10.0_wp * mol%xyz
+ expected = 0
+ end if
+ call compliance_driver(env, mol, hess)
+ call env%check(failed)
+ call check(error, .not. failed)
+ if (allocated(error)) exit
+
+ counts = 0
+ open(newunit=dat_unit, file="compliance.dat", status="old")
+ do
+ read(dat_unit, "(a)", iostat=ios) line
+ if (ios /= 0) exit
+ if (index(line, "# coord ") /= 1) cycle
+ read(line(8:), *) kind
+ counts(1) = counts(1) + 1
+ if (kind == "bond") counts(2) = counts(2) + 1
+ if (kind == "ang") counts(3) = counts(3) + 1
+ end do
+ close(dat_unit, status="delete")
+ do i = 1, 3
+ call check(error, counts(i), expected(i))
+ if (allocated(error)) exit
+ end do
+ if (allocated(error)) exit
+ end do
+ close(unit)
+
+end subroutine test_compliance_driver_redundant
+
+
+!> Locks covalent neighbour-list generation: canonical half-list storage,
+!> replacement semantics, and the supported-cutoff tolerance margin for
+!> getNeighs.
+subroutine test_covalent_neighbour_list(error)
+ !> Failure report, allocated when the check fails.
+ type(error_type), allocatable, intent(out) :: error
+
+ integer, parameter :: nat = 4
+ integer, parameter :: at(nat) = [8, 1, 1, 1]
+ real(wp), parameter :: xyz(3, nat) = reshape([ &
+ & 0.0_wp, 0.0_wp, 0.0_wp, &
+ & 0.0_wp, 1.8_wp, 0.0_wp, &
+ & 1.7_wp, -0.5_wp, 0.0_wp, &
+ & 10.0_wp, 0.0_wp, 0.0_wp], shape(xyz))
+ integer, parameter :: deg_ref(nat) = [2, 1, 1, 0]
+ real(wp), parameter :: dist2_ref(2) = [3.14_wp, 3.24_wp]
+ real(wp), parameter :: thr = 1.0e-10_wp
+
+ type(TNeighbourList) :: neigh_list
+ type(graph_type) :: graph
+ integer :: neighs(nat), i
+
+ ! Empty list must have zero supported cutoff.
+ call init(neigh_list, 0)
+ call neigh_list%generate_covalent(at(:0), xyz(:, :0))
+ call check(error, neigh_list%cutoff, 0.0_wp)
+ if (allocated(error)) return
+
+ ! O/H/H plus one isolated H: two O-H bonds, no H-H bond
+ call init(neigh_list, nat)
+ call neigh_list%generate_covalent(at, xyz)
+ call check(error, neigh_list%neighs(1), 2)
+ if (allocated(error)) return
+ call check(error, all(neigh_list%neighs(2:4) == 0), .true.)
+ if (allocated(error)) return
+ call check(error, all(neigh_list%image == [(i, i = 1, nat)]), .true.)
+ if (allocated(error)) return
+ call check(error, all(neigh_list%iNeigh(0, :) == [(i, i = 1, nat)]), .true.)
+ if (allocated(error)) return
+ call check(error, all(abs(neigh_list%dist2(1:2, 1) - dist2_ref) < thr), .true.)
+ if (allocated(error)) return
+
+ call init(graph, neigh_list)
+ do i = 1, nat
+ call check(error, graph%neighs(i), deg_ref(i))
+ if (allocated(error)) return
+ end do
+ call check(error, graph%parent(4), 0)
+ if (allocated(error)) return
+
+ ! Replacement semantics: regenerating the same geometry must replace
+ ! (not append to) the previous list
+ call init(neigh_list, nat)
+ call neigh_list%generate_covalent(at, xyz)
+ call neigh_list%generate_covalent(at, xyz)
+ call check(error, neigh_list%neighs(1), 2)
+ if (allocated(error)) return
+ call check(error, all(neigh_list%neighs(2:4) == 0), .true.)
+ if (allocated(error)) return
+
+ ! Default getNeighs must return every stored neighbour, also when the
+ ! bond sits infinitesimally below the largest pair threshold
+ call init(neigh_list, 3)
+ call neigh_list%generate_covalent([1, 1, 1], &
+ & reshape([0.0_wp, 0.0_wp, 0.0_wp, &
+ & 0.5_wp, 0.0_wp, 0.0_wp, &
+ & 2.0_wp*covalent_radius_d3(1) - 1.0e-12_wp, 0.0_wp, 0.0_wp], &
+ & [3, 3]))
+ call neigh_list%getNeighs(neighs(:3))
+ call check(error, neighs(1), 2)
+ if (allocated(error)) return
+ call check(error, neighs(2), 1)
+ if (allocated(error)) return
+ call check(error, neighs(3), 0)
+
+end subroutine test_covalent_neighbour_list
+
+
+!> Build a symmetric neighbour list from undirected atom pairs.
+function build_neigh_list(nat, pairs) result(neigh_list)
+ !> Number of atoms.
+ integer, intent(in) :: nat
+ !> Undirected atom pairs, dimension (2, number of pairs).
+ integer, intent(in) :: pairs(:, :)
+ !> Neighbour list with both adjacency columns filled.
+ type(TNeighbourList) :: neigh_list
+
+ integer :: i, j, k
+
+ call init(neigh_list, nat)
+ neigh_list%neighs = 0
+ neigh_list%iNeigh = 0
+ neigh_list%image = [(i, i = 1, nat)]
+ do k = 1, size(pairs, 2)
+ i = pairs(1, k)
+ j = pairs(2, k)
+ neigh_list%neighs(i) = neigh_list%neighs(i) + 1
+ neigh_list%iNeigh(neigh_list%neighs(i), i) = j
+ neigh_list%neighs(j) = neigh_list%neighs(j) + 1
+ neigh_list%iNeigh(neigh_list%neighs(j), j) = i
+ end do
+
+end function build_neigh_list
+
end module test_hessian
diff --git a/test/unit/test_model_hessian.f90 b/test/unit/test_model_hessian.f90
index b37024035..d6d875c93 100644
--- a/test/unit/test_model_hessian.f90
+++ b/test/unit/test_model_hessian.f90
@@ -30,8 +30,9 @@ module test_model_hessian
use xtb_type_molecule, only : TMolecule
use xtb_modelhessian_eeq, only : add_eeq_hessian
use xtb_modelhessian_type, only : TModelHessian
- use xtb_modelhessian_lindh, only : TLindhModelHessian, TLindhD2ModelHessian
- use xtb_modelhessian_swart, only : TSwartModelHessian
+ use xtb_modelhessian_lindh, only : TLindhModelHessian, &
+ & newLindhModelHessian, newLindhD2ModelHessian
+ use xtb_modelhessian_swart, only : newSwartModelHessian
use xtb_type_param, only : chrg_parameter
use xtb_type_setvar, only : modhess_setvar
use xtb_test_molstock, only : getMolecule
@@ -72,6 +73,7 @@ subroutine collect_model_hessian(testsuite)
! Misc
new_unittest("model_hessian_dense", test_model_hessian_dense), &
new_unittest("model_hessian_charge", test_model_hessian_charge), &
+ new_unittest("lindh_bend_scale_invariance", test_lindh_bend_scale_invariance), &
new_unittest("eeq_addition", test_eeq_addition), &
new_unittest("gff_h2o", test_gff_h2o), &
new_unittest("torsion_reversal_key", test_torsion_reversal_key), &
@@ -130,23 +132,24 @@ subroutine compute_mh_packed(mol, variant, modh, hess_packed)
call init(env)
n3 = 3 * mol%n
allocate(hess_packed(n3*(n3 + 1)/2))
- call new_model_hessian(variant, model_hessian)
- call model_hessian%compute(env, mol%xyz, mol%n, hess_packed, mol%at, modh)
+ call new_model_hessian(variant, modh, model_hessian)
+ call model_hessian%compute(env, mol%xyz, mol%n, hess_packed, mol%at)
end subroutine compute_mh_packed
!> Allocate model Hessian implementation for a test variant
-subroutine new_model_hessian(variant, model_hessian)
+subroutine new_model_hessian(variant, modh, model_hessian)
integer, intent(in) :: variant
+ type(modhess_setvar), intent(in) :: modh
class(TModelHessian), allocatable, intent(out) :: model_hessian
select case (variant)
case(VAR_LINDH_D2)
- allocate(TLindhD2ModelHessian :: model_hessian)
+ allocate(model_hessian, source=newLindhD2ModelHessian(modh))
case(VAR_LINDH)
- allocate(TLindhModelHessian :: model_hessian)
+ allocate(model_hessian, source=newLindhModelHessian(modh))
case(VAR_SWART)
- allocate(TSwartModelHessian :: model_hessian)
+ allocate(model_hessian, source=newSwartModelHessian(modh))
case default
error stop "Unknown model Hessian variant"
end select
@@ -169,9 +172,9 @@ subroutine test_model_hessian_dense(error)
allocate(hess_packed(n3*(n3 + 1)/2), hess_dense(n3, n3))
do variant = VAR_LINDH_D2, VAR_SWART
- call new_model_hessian(variant, model_hessian)
- call model_hessian%compute(env, mol%xyz, mol%n, hess_packed, mol%at, default_modh())
- call model_hessian%compute(env, mol%xyz, mol%n, hess_dense, mol%at, default_modh())
+ call new_model_hessian(variant, default_modh(), model_hessian)
+ call model_hessian%compute(env, mol%xyz, mol%n, hess_packed, mol%at)
+ call model_hessian%compute(env, mol%xyz, mol%n, hess_dense, mol%at)
ij = 0
do i = 1, n3
@@ -209,11 +212,12 @@ subroutine test_model_hessian_charge(error)
call add_eeq_hessian(env, mol%n, mol%at, mol%xyz, 0.0_wp, chrgeq, 0.1_wp, contribution)
do variant = VAR_LINDH_D2, VAR_SWART
- call new_model_hessian(variant, model_hessian)
modh = default_modh()
- call model_hessian%compute(env, mol%xyz, mol%n, base, mol%at, modh)
+ call new_model_hessian(variant, modh, model_hessian)
+ call model_hessian%compute(env, mol%xyz, mol%n, base, mol%at)
modh%kq = 0.1_wp
- call model_hessian%compute(env, mol%xyz, mol%n, charged, mol%at, modh)
+ call new_model_hessian(variant, modh, model_hessian)
+ call model_hessian%compute(env, mol%xyz, mol%n, charged, mol%at)
do i = 1, size(base)
call check(error, charged(i), base(i) + contribution(i), &
@@ -224,6 +228,55 @@ subroutine test_model_hessian_charge(error)
end subroutine test_model_hessian_charge
+!> Near-linear bend selection is independent of uniformly scaled arm lengths
+subroutine test_lindh_bend_scale_invariance(error)
+ type(error_type), allocatable, intent(out) :: error
+
+ real(wp), parameter :: delta = 5.0e-8_wp
+ integer, parameter :: at(3) = 1
+ real(wp) :: xyz(3, 3), hess(45), without_bend(45), normalized(45, 2)
+ real(wp) :: scale, pair_amplitude, bend_force, bend_trace
+ type(TEnvironment) :: env
+ type(TLindhModelHessian) :: model_hessian
+ type(modhess_setvar) :: modh
+ integer :: i, diagonal
+
+ call init(env)
+ do i = 1, 2
+ scale = real(i, wp)
+ xyz = 0.0_wp
+ xyz(:, 1) = scale * [1.0_wp, 0.0_wp, 0.0_wp]
+ xyz(:, 3) = scale * [-1.0_wp, delta, 0.0_wp]
+
+ modh = default_modh()
+ modh%s6 = 0.0_wp
+ bend_force = modh%kf
+ model_hessian = newLindhModelHessian(modh)
+ call model_hessian%compute(env, xyz, 3, hess, at)
+ pair_amplitude = model_hessian%pair_factor(1, 1, scale**2, 0.0_wp, .false.) &
+ & * model_hessian%pair_factor(1, 1, scale**2*(1.0_wp + delta**2), &
+ & 0.0_wp, .false.)
+
+ modh%kf = 0.0_wp
+ model_hessian = newLindhModelHessian(modh)
+ call model_hessian%compute(env, xyz, 3, without_bend, at)
+ normalized(:, i) = (hess - without_bend) * scale**2 &
+ & / (bend_force * pair_amplitude)
+ end do
+
+ ! The old dimensional cross-product cutoff chose linear at 1 Bohr and
+ ! ordinary at 2 Bohr for this same dimensionless angle.
+ call compare(error, normalized(:, 1), normalized(:, 2))
+ if (allocated(error)) return
+
+ bend_trace = 0.0_wp
+ do diagonal = 1, 9
+ bend_trace = bend_trace + normalized(diagonal*(diagonal + 1)/2, 1)
+ end do
+ call check(error, bend_trace, 6.0_wp, thr=1.0e-6_wp)
+end subroutine test_lindh_bend_scale_invariance
+
+
!> Check EEQ utility adds to, rather than replaces, packed Hessian values
subroutine test_eeq_addition(error)
type(error_type), allocatable, intent(out) :: error
@@ -319,8 +372,8 @@ subroutine test_gff_h2o(error)
n3 = 3 * mol%n
allocate(hessian(n3*(n3 + 1)/2))
modh = default_modh()
- model_hessian = newGFFModelHessian(calc%param, calc%topo, calc%neigh)
- call model_hessian%compute(env, mol%xyz, mol%n, hessian, mol%at, modh)
+ model_hessian = newGFFModelHessian(calc%param, calc%topo, calc%neigh, modh)
+ call model_hessian%compute(env, mol%xyz, mol%n, hessian, mol%at)
call compare(error, ref_gff_h2o, hessian)
end subroutine test_gff_h2o