Skip to content
This repository was archived by the owner on May 5, 2026. It is now read-only.

Commit a0bbb11

Browse files
ss0832ss0832
authored
Add files via upload
1 parent ef6eb06 commit a0bbb11

3 files changed

Lines changed: 335 additions & 3 deletions

File tree

‎multioptpy/ModelHessian/approx_hessian.py‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
from multioptpy.ModelHessian.schlegel import SchlegelApproxHessian
1515
from multioptpy.ModelHessian.schlegeld3 import SchlegelD3ApproxHessian
1616
from multioptpy.ModelHessian.schlegeld4 import SchlegelD4ApproxHessian
17+
from multioptpy.ModelHessian.swart import SwartApproxHessian
1718
from multioptpy.ModelHessian.swartd2 import SwartD2ApproxHessian
1819
from multioptpy.ModelHessian.swartd3 import SwartD3ApproxHessian
1920
from multioptpy.ModelHessian.swartd4 import SwartD4ApproxHessian
@@ -62,9 +63,12 @@ def main(self, coord, element_list, cart_gradient, approx_hess_type="lindh2007d3
6263
elif "swartd4" in approx_hess_type.lower():
6364
SWHD4 = SwartD4ApproxHessian()
6465
hess_proj = SWHD4.main(coord, element_list, cart_gradient)
65-
elif "swart" in approx_hess_type.lower():
66+
elif "swartd2" in approx_hess_type.lower():
6667
SWH = SwartD2ApproxHessian()
6768
hess_proj = SWH.main(coord, element_list, cart_gradient)
69+
elif "swart" in approx_hess_type.lower():
70+
SW = SwartApproxHessian()
71+
hess_proj = SW.main(coord, element_list, cart_gradient)
6872
elif "lindh2007d3" in approx_hess_type.lower():
6973
LH2007D3 = Lindh2007D3ApproxHessian()
7074
hess_proj = LH2007D3.main(coord, element_list, cart_gradient)

‎multioptpy/ModelHessian/o1numhess.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
from scipy.spatial.distance import cdist
66

77

8-
from multioptpy.ModelHessian.swartd3 import SwartD3ApproxHessian
8+
from multioptpy.ModelHessian.swart import SwartApproxHessian
99
from multioptpy.Parameters.parameter import UnitValueLib
1010

1111

@@ -57,7 +57,7 @@ def __init__(self, calculation_engine, element_list, charge_mult, method,
5757
self.bohr2ang = self.uvl.bohr2angstroms
5858

5959
# Fischer Hessian Calculator (for initial guess)
60-
self.fischer_calc = SwartD3ApproxHessian()
60+
self.fischer_calc = SwartApproxHessian()
6161

6262
# Paper parameters
6363
self.lam = 1.0e-2 # Regularization parameter λ (Eq. 28)

‎multioptpy/ModelHessian/swart.py‎

Lines changed: 328 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,328 @@
1+
import numpy as np
2+
import math
3+
4+
from multioptpy.Utils.calc_tools import Calculationtools
5+
6+
# ============================================================================
7+
# Library / Helper Functions (Acting as multioptpy replacements for this logic)
8+
# ============================================================================
9+
10+
# Periodic table (element symbols)
11+
PERIODIC_TABLE = [
12+
'X', # Dummy element at index 0
13+
'H', 'He',
14+
'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne',
15+
'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar',
16+
'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn',
17+
'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr',
18+
'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd',
19+
'In', 'Sn', 'Sb', 'Te', 'I', 'Xe',
20+
'Cs', 'Ba',
21+
'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu',
22+
'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg',
23+
'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn'
24+
]
25+
26+
# Covalent radii in Bohr (from Pyykkö & Atsumi, Chem. Eur. J. 2009, 15, 186)
27+
COVALENT_RADII = {
28+
'H': 0.59, 'He': 0.54,
29+
'Li': 2.43, 'Be': 1.72, 'B': 1.53, 'C': 1.40, 'N': 1.34, 'O': 1.25, 'F': 1.18, 'Ne': 1.14,
30+
'Na': 2.89, 'Mg': 2.53, 'Al': 2.19, 'Si': 2.10, 'P': 2.04, 'S': 1.97, 'Cl': 1.87, 'Ar': 1.82,
31+
'K': 3.42, 'Ca': 3.06, 'Sc': 2.85, 'Ti': 2.70, 'V': 2.55, 'Cr': 2.49, 'Mn': 2.49,
32+
'Fe': 2.44, 'Co': 2.38, 'Ni': 2.32, 'Cu': 2.42, 'Zn': 2.40,
33+
'Ga': 2.27, 'Ge': 2.19, 'As': 2.17, 'Se': 2.10, 'Br': 2.04, 'Kr': 2.06,
34+
'Rb': 3.70, 'Sr': 3.40, 'Y': 3.21, 'Zr': 2.98, 'Nb': 2.85, 'Mo': 2.72, 'Tc': 2.61,
35+
'Ru': 2.55, 'Rh': 2.51, 'Pd': 2.55, 'Ag': 2.68, 'Cd': 2.72,
36+
'In': 2.61, 'Sn': 2.55, 'Sb': 2.51, 'Te': 2.48, 'I': 2.44, 'Xe': 2.48,
37+
'Cs': 4.03, 'Ba': 3.59,
38+
'La': 3.34, 'Ce': 3.25, 'Pr': 3.23, 'Nd': 3.21, 'Pm': 3.19, 'Sm': 3.17, 'Eu': 3.17,
39+
'Gd': 3.15, 'Tb': 3.13, 'Dy': 3.13, 'Ho': 3.11, 'Er': 3.11, 'Tm': 3.09, 'Yb': 3.09, 'Lu': 3.06,
40+
'Hf': 2.89, 'Ta': 2.76, 'W': 2.61, 'Re': 2.49, 'Os': 2.46, 'Ir': 2.42, 'Pt': 2.42, 'Au': 2.55, 'Hg': 2.72,
41+
'Tl': 2.68, 'Pb': 2.68, 'Bi': 2.68, 'Po': 2.61, 'At': 2.57, 'Rn': 2.63,
42+
}
43+
44+
def covalent_radii_lib(element):
45+
"""Mimics multioptpy covalent_radii_lib"""
46+
elem = element.capitalize()
47+
if elem not in COVALENT_RADII:
48+
raise ValueError(f"Unknown element: {element}")
49+
return COVALENT_RADII[elem]
50+
51+
def bmat_bond(xyz: np.ndarray, i: int, j: int) -> np.ndarray:
52+
vec = xyz[i, :] - xyz[j, :]
53+
l = np.linalg.norm(vec)
54+
B = np.zeros(6)
55+
B[0:3] = vec / l
56+
B[3:6] = -vec / l
57+
return B
58+
59+
def bmat_angle(xyz: np.ndarray, i: int, j: int, k: int) -> np.ndarray:
60+
vec1 = xyz[i, :] - xyz[j, :]
61+
vec2 = xyz[k, :] - xyz[j, :]
62+
l1 = np.linalg.norm(vec1)
63+
l2 = np.linalg.norm(vec2)
64+
65+
if l1 < 1e-10 or l2 < 1e-10:
66+
return np.zeros(9)
67+
68+
nvec1 = vec1 / l1
69+
nvec2 = vec2 / l2
70+
71+
dot_prod = np.dot(nvec1, nvec2)
72+
sin_theta_sq = max(1e-15, 1.0 - dot_prod**2)
73+
74+
dl = np.zeros((2, 6))
75+
dl[0, 0:3] = nvec1
76+
dl[0, 3:6] = -nvec1
77+
dl[1, 0:3] = nvec2
78+
dl[1, 3:6] = -nvec2
79+
80+
dnvec = np.zeros((2, 3, 6))
81+
for ii in range(6):
82+
dnvec[0, :, ii] = -nvec1 * dl[0, ii] / l1
83+
dnvec[1, :, ii] = -nvec2 * dl[1, ii] / l2
84+
85+
for ii in range(3):
86+
dnvec[0, ii, ii] += 1.0 / l1
87+
dnvec[1, ii, ii] += 1.0 / l2
88+
dnvec[0, ii, ii + 3] -= 1.0 / l1
89+
dnvec[1, ii, ii + 3] -= 1.0 / l2
90+
91+
dinprod = np.zeros(9)
92+
for ii in range(3):
93+
dinprod[ii] = np.dot(dnvec[0, :, ii], nvec2)
94+
dinprod[ii + 3] = np.dot(dnvec[0, :, ii + 3], nvec2) + np.dot(dnvec[1, :, ii + 3], nvec1)
95+
dinprod[ii + 6] = np.dot(dnvec[1, :, ii], nvec1)
96+
97+
B = -dinprod / math.sqrt(sin_theta_sq)
98+
return B
99+
100+
def bmat_linear_angle(xyz: np.ndarray, i: int, j: int, k: int) -> np.ndarray:
101+
vec1 = xyz[i, :] - xyz[j, :]
102+
vec2 = xyz[k, :] - xyz[j, :]
103+
l1 = np.linalg.norm(vec1)
104+
l2 = np.linalg.norm(vec2)
105+
106+
if l1 < 1e-10 or l2 < 1e-10:
107+
return np.zeros((2, 9))
108+
109+
nvec1 = vec1 / l1
110+
nvec2 = vec2 / l2
111+
112+
vn = np.cross(vec1, vec2)
113+
nvn = np.linalg.norm(vn)
114+
115+
if nvn < 1e-15:
116+
vn = np.array([1.0, 0.0, 0.0])
117+
vn = vn - np.dot(vn, vec1) / l1**2 * vec1
118+
nvn = np.linalg.norm(vn)
119+
120+
if nvn < 1e-15:
121+
vn = np.array([0.0, 1.0, 0.0])
122+
vn = vn - np.dot(vn, vec1) / l1**2 * vec1
123+
nvn = np.linalg.norm(vn)
124+
125+
vn = vn / nvn
126+
vn2 = np.cross(vec1 - vec2, vn)
127+
vn2 = vn2 / np.linalg.norm(vn2)
128+
129+
B = np.zeros((2, 9))
130+
131+
B[1, 0:3] = vn / l1
132+
B[1, 6:9] = vn / l2
133+
B[1, 3:6] = -B[1, 0:3] - B[1, 6:9]
134+
135+
B[0, 0:3] = vn2 / l1
136+
B[0, 6:9] = vn2 / l2
137+
B[0, 3:6] = -B[0, 0:3] - B[0, 6:9]
138+
139+
return B
140+
141+
def bond_length(xyz: np.ndarray, i: int, j: int) -> float:
142+
return np.linalg.norm(xyz[i, :] - xyz[j, :])
143+
144+
145+
# ============================================================================
146+
# Main Class
147+
# ============================================================================
148+
149+
class SwartApproxHessian:
150+
def __init__(self):
151+
# Swart's model Hessian parameters from swart.py
152+
# ref.: M. Swart, F. M. Bickelhaupt, Int. J. Quantum Chem., 2006, 106, 2536
153+
# ref.: O1NumHess paper (arXiv:2508.07544v1) Appendix B
154+
155+
self.wthr = 0.3
156+
self.f = 0.12
157+
self.tolth = 0.2
158+
159+
# Derived parameters
160+
self.eps1 = self.wthr**2
161+
self.eps2 = self.wthr**2 / math.exp(1)
162+
163+
self.cart_hess = None
164+
return
165+
166+
def screening_function(self, distance, cov_radius_sum):
167+
"""Screening function ρ_AB = exp(1 - R_AB / (r_A + r_B))"""
168+
return math.exp(1.0 - distance / cov_radius_sum)
169+
170+
def _cos_angle(self, xyz, i, j, k):
171+
"""Calculate cos(angle i-j-k)."""
172+
vec1 = xyz[i, :] - xyz[j, :]
173+
vec2 = xyz[k, :] - xyz[j, :]
174+
175+
l1 = np.linalg.norm(vec1)
176+
l2 = np.linalg.norm(vec2)
177+
178+
if l1 < 1e-10 or l2 < 1e-10:
179+
return 1.0
180+
181+
cos_theta = np.dot(vec1, vec2) / (l1 * l2)
182+
return np.clip(cos_theta, -1.0, 1.0)
183+
184+
def swart_bond(self, coord, element_list):
185+
"""
186+
Calculate bond stretching contributions.
187+
Paper: All atom pairs are treated as bonds, with force constant H_int = 0.35 * ρ_AB^3
188+
"""
189+
N = len(coord)
190+
# Re-calculate screening matrix locally as per original flow, or reuse if optimized.
191+
# To strictly follow algorithmic configuration of swart.py, we perform calculations as needed.
192+
193+
for i in range(N):
194+
for j in range(i + 1, N):
195+
cov_sum = covalent_radii_lib(element_list[i]) + covalent_radii_lib(element_list[j])
196+
dist = bond_length(coord, i, j)
197+
198+
screen = self.screening_function(dist, cov_sum)
199+
200+
# Force constant (Paper: 0.35 * ρ^3)
201+
H_int = 0.35 * screen**3
202+
203+
# Wilson B matrix for bond
204+
B = bmat_bond(coord, i, j)
205+
206+
# Indices for atoms i and j
207+
range_i = list(range(3 * i, 3 * (i + 1)))
208+
range_j = list(range(3 * j, 3 * (j + 1)))
209+
range_ij = range_i + range_j
210+
211+
# Add contribution: H += H_int * B * B^T
212+
# Using np.ix_ to match logic of swart.py efficient update
213+
self.cart_hess[np.ix_(range_ij, range_ij)] += H_int * np.outer(B, B)
214+
return
215+
216+
def swart_angle(self, coord, element_list):
217+
"""
218+
Calculate angle bending contributions.
219+
Paper: Force constant H_int = 0.075 * s_ij,jk^2 * (f + (1-f)*sin(θ))^2
220+
"""
221+
N = len(coord)
222+
223+
# Precompute screening for efficiency inside the angle loop (consistent with swart.py logic)
224+
screen_matrix = np.zeros((N, N))
225+
for i in range(N):
226+
for j in range(N):
227+
if i == j: continue
228+
cov_sum = covalent_radii_lib(element_list[i]) + covalent_radii_lib(element_list[j])
229+
dist = bond_length(coord, i, j)
230+
screen_matrix[i, j] = self.screening_function(dist, cov_sum)
231+
232+
for i in range(N):
233+
for j in range(N):
234+
if i == j:
235+
continue
236+
237+
# Check if i-j screening is sufficient
238+
if screen_matrix[i, j] < self.eps2:
239+
continue
240+
241+
for k in range(i + 1, N):
242+
if k == j:
243+
continue
244+
245+
# Combined screening
246+
s_ij_jk = screen_matrix[i, j] * screen_matrix[j, k]
247+
if s_ij_jk < self.eps1:
248+
continue
249+
250+
# Calculate angle
251+
cos_theta = self._cos_angle(coord, i, j, k)
252+
sin_theta = math.sqrt(max(0.0, 1.0 - cos_theta**2))
253+
254+
# Force constant (Paper Eq. with modified coefficient 0.075)
255+
H_int = 0.075 * s_ij_jk**2 * (self.f + (1 - self.f) * sin_theta)**2
256+
257+
# Check for linear or zero angle
258+
if cos_theta > 1 - self.tolth:
259+
th1 = 1.0 - cos_theta
260+
else:
261+
th1 = 1.0 + cos_theta
262+
263+
range_i = list(range(3 * i, 3 * (i + 1)))
264+
range_j = list(range(3 * j, 3 * (j + 1)))
265+
range_k = list(range(3 * k, 3 * (k + 1)))
266+
range_ijk = range_i + range_j + range_k
267+
268+
if th1 < self.tolth:
269+
# Near-linear or near-zero angle
270+
scale_lin = (1.0 - (th1 / self.tolth)**2)**2
271+
272+
if cos_theta > 1 - self.tolth:
273+
# Near 180 degrees (linear)
274+
B_lin = bmat_linear_angle(coord, i, j, k)
275+
B = bmat_angle(coord, i, j, k)
276+
277+
# Scale between linear and normal
278+
B_combined = scale_lin * B_lin[0, :] + (1.0 - scale_lin) * B
279+
280+
# Add linear bending mode
281+
self.cart_hess[np.ix_(range_ijk, range_ijk)] += H_int * np.outer(B_lin[1, :], B_lin[1, :])
282+
# Add combined mode
283+
self.cart_hess[np.ix_(range_ijk, range_ijk)] += H_int * np.outer(B_combined, B_combined)
284+
else:
285+
# Near 0 degrees
286+
B = bmat_angle(coord, i, j, k)
287+
B_scaled = (1.0 - scale_lin) * B
288+
self.cart_hess[np.ix_(range_ijk, range_ijk)] += H_int * np.outer(B_scaled, B_scaled)
289+
else:
290+
# Normal angle
291+
B = bmat_angle(coord, i, j, k)
292+
self.cart_hess[np.ix_(range_ijk, range_ijk)] += H_int * np.outer(B, B)
293+
return
294+
295+
def swart_dihedral_angle(self, coord, element_list):
296+
"""
297+
Calculate dihedral angle contributions.
298+
Note: swart.py logic states 'No dihedral terms (paper notes they are not necessary)'.
299+
Method exists to maintain interface compatibility with SwartD3ApproxHessian.
300+
"""
301+
pass
302+
303+
def swart_out_of_plane(self, coord, element_list):
304+
"""
305+
Calculate out-of-plane bending contributions.
306+
Note: swart.py logic states 'No out-of-plane terms'.
307+
Method exists to maintain interface compatibility with SwartD3ApproxHessian.
308+
"""
309+
pass
310+
311+
def main(self, coord, element_list, cart_gradient=None):
312+
"""
313+
Main method to calculate the approximate Hessian using Swart's model.
314+
Note: logic expects 'coord' in Bohr.
315+
"""
316+
print("Generating Swart's approximate hessian (O1NumHess variant)...")
317+
self.cart_hess = np.zeros((len(coord)*3, len(coord)*3), dtype="float64")
318+
319+
self.swart_bond(coord, element_list)
320+
self.swart_angle(coord, element_list)
321+
self.swart_dihedral_angle(coord, element_list)
322+
self.swart_out_of_plane(coord, element_list)
323+
324+
# Note: swart.py does not implement projection of translational/rotational modes
325+
# in the main logic, returning the raw internal Hessian matrix directly.
326+
hess_proj = Calculationtools().project_out_hess_tr_and_rot_for_coord(self.cart_hess, element_list, coord)
327+
return hess_proj
328+

0 commit comments

Comments
 (0)