Skip to content

Commit c3a1693

Browse files
ss0832ss0832
authored
Add files via upload
1 parent 7b29661 commit c3a1693

3 files changed

Lines changed: 462 additions & 38 deletions

File tree

‎multioptpy/ModelFunction/binary_image_ts_search_model_function.py‎

Lines changed: 111 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -1,47 +1,140 @@
11
import numpy as np
22

3-
#J. Chem. Phys. 157, 124107 (2022)
4-
#https://doi.org/10.1063/5.0102145
3+
# J. Chem. Phys. 157, 124107 (2022)
4+
# https://doi.org/10.1063/5.0102145
55

66
class BITSSModelFunction:
77
def __init__(self, geom_num_list_1, geom_num_list_2):
88
self.f = 0.5
99
self.alpha = 10.0
10-
self.beta = 0.1
11-
self.d = np.linalg.norm(geom_num_list_1 - geom_num_list_2)
10+
# self.beta controls the strength of the distance constraint.
11+
# Smaller beta -> Larger kappa_d -> Stronger attractive force.
12+
# Default: 0.1 -> Modified: 0.06 for stronger constraint.
13+
self.beta = 0.06
1214

13-
return
15+
# Initial distance
16+
diff = geom_num_list_1 - geom_num_list_2
17+
self.d = np.linalg.norm(diff)
1418

19+
# Initialize variables to avoid AttributeError if calc_hess/grad called before iter % 500 == 0
20+
self.kappa_e = 0.0
21+
self.kappa_d = 0.0
22+
self.E_B = 0.0
23+
24+
# Store vector info for Hessian calculation
25+
self.diff_vec = diff.reshape(-1, 1) # (3N, 1)
26+
self.current_dist = self.d
27+
1528
def calc_energy(self, energy_1, energy_2, geom_num_list_1, geom_num_list_2, gradient_1, gradient_2, iter):
16-
current_distance = np.linalg.norm(geom_num_list_1 - geom_num_list_2)
29+
# Update distances
30+
diff_vec = geom_num_list_1 - geom_num_list_2
31+
current_distance = np.linalg.norm(diff_vec)
32+
33+
# Update parameters periodically
1734
if iter % 500 == 0:
18-
1935
self.E_B = abs(energy_1 - energy_2)
36+
# Avoid division by zero
2037
self.kappa_e = self.alpha / (2.0 * self.E_B + 1e-10)
2138

22-
unit_vec = (geom_num_list_1 - geom_num_list_2) / current_distance
23-
24-
39+
unit_vec = diff_vec / (current_distance + 1e-10)
2540

26-
proj_grad_1 = gradient_1 * unit_vec * (-1)
27-
proj_grad_2 = gradient_2 * unit_vec
41+
# Project gradients onto the distance vector direction
42+
# grad_1 is (N, 3), unit_vec is (N, 3). Element-wise mult then sum gives dot product.
43+
proj_grad_1 = np.sum(gradient_1 * (-1) * unit_vec)
44+
proj_grad_2 = np.sum(gradient_2 * unit_vec)
2845

29-
a = np.sqrt(np.linalg.norm(proj_grad_1) + np.linalg.norm(proj_grad_2)) / (2 ** 1.5 * self.beta * self.d + 1e-10)
46+
# Eq. (5) logic
47+
grad_norm_term = np.sqrt(proj_grad_1**2 + proj_grad_2**2)
48+
a = grad_norm_term / (2 ** 1.5 * self.beta * self.d + 1e-10)
3049
b = self.E_B / (self.beta * self.d ** 2 + 1e-10)
3150
self.kappa_d = max(a, b)
32-
self.d = np.linalg.norm(geom_num_list_1 - geom_num_list_2)
51+
52+
# Reset target distance d to current distance at update step
53+
self.d = current_distance
3354

55+
# Reduce target distance
3456
self.d = max((1.0 - self.f) * self.d, 1e-10)
35-
energy = energy_1 + energy_2 + self.kappa_e * (energy_1 + energy_2) ** 2 + self.kappa_d * (current_distance - self.d) ** 2
57+
58+
# Calculate BITSS Energy
59+
# Formula: E1 + E2 + ke * (E1 - E2)^2 + kd * (d - d0)^2
60+
energy = energy_1 + energy_2 + self.kappa_e * (energy_1 - energy_2) ** 2 + self.kappa_d * (current_distance - self.d) ** 2
3661

3762
return energy
3863

3964
def calc_grad(self, energy_1, energy_2, geom_num_list_1, geom_num_list_2, gradient_1, gradient_2):
65+
# Calculate vector r = x1 - x2 and distance d
4066
current_vec = geom_num_list_1 - geom_num_list_2
4167
current_dist = np.linalg.norm(current_vec) + 1e-10
4268

43-
bitss_grad_1 = gradient_1 + gradient_2 + 2.0 * self.kappa_e * (energy_1 - energy_2) * (gradient_1 - gradient_2) + current_vec * 2.0 * self.kappa_d * (current_dist - self.d) / current_dist
69+
# Store for calc_hess (flattened for matrix ops)
70+
self.diff_vec = current_vec.reshape(-1, 1)
71+
self.current_dist = current_dist
72+
73+
# Common terms
74+
delta_E = energy_1 - energy_2
75+
dist_diff = current_dist - self.d
76+
77+
# Gradient term for distance: 2 * kd * (d - d0) * (r / d)
78+
grad_dist_term = current_vec * 2.0 * self.kappa_d * dist_diff / current_dist
79+
80+
# Gradient term for energy: 2 * ke * (E1 - E2) * (g1 - g2)
81+
# Total Gradient 1: g1 + 2*ke*dE*g1 + dist_term
82+
bitss_grad_1 = gradient_1 * (1.0 + 2.0 * self.kappa_e * delta_E) + grad_dist_term
83+
84+
# Total Gradient 2: g2 + 2*ke*dE*(-g2) - dist_term (since d(r)/dx2 = -r/d)
85+
bitss_grad_2 = gradient_2 * (1.0 - 2.0 * self.kappa_e * delta_E) - grad_dist_term
86+
87+
return bitss_grad_1, bitss_grad_2
4488

45-
bitss_grad_2 = gradient_1 + gradient_2 + 2.0 * self.kappa_e * (energy_1 - energy_2) * (gradient_1 - gradient_2) - current_vec * 2.0 * self.kappa_d * (current_dist - self.d) / current_dist
89+
def calc_hess(self, energy_1, energy_2, grad_1, grad_2, hess_1, hess_2):
90+
"""
91+
Calculate the 6N x 6N Hessian matrix for BITSS.
92+
H = [ H11 H12 ]
93+
[ H21 H22 ]
94+
"""
95+
# Ensure inputs are flattened (3N, 1) or (3N, 3N)
96+
N3 = self.diff_vec.shape[0]
97+
g1 = grad_1.reshape(N3, 1)
98+
g2 = grad_2.reshape(N3, 1)
99+
100+
delta_E = energy_1 - energy_2
101+
dist_diff = self.current_dist - self.d
102+
103+
# --- Distance Constraint Hessian Terms ---
104+
# Vd = kd * (d - d0)^2
105+
# P = r * r.T / d^2 (Projection onto bond axis)
106+
r = self.diff_vec
107+
d = self.current_dist
108+
P = np.dot(r, r.T) / (d**2)
109+
I = np.eye(N3)
110+
# H_dist_block = 2*kd * [ P + (d-d0)/d * (I - P) ]
111+
term_d = P + (dist_diff / d) * (I - P)
112+
H_dist = 2.0 * self.kappa_d * term_d
113+
114+
# --- Total Hessian Blocks ---
115+
116+
# Block 11: d^2 E / dx1^2
117+
# = H1 * (1 + 2*ke*dE) + 2*ke * g1 * g1.T + H_dist
118+
H11 = hess_1 * (1.0 + 2.0 * self.kappa_e * delta_E) + \
119+
2.0 * self.kappa_e * np.dot(g1, g1.T) + \
120+
H_dist
121+
122+
# Block 22: d^2 E / dx2^2
123+
# = H2 * (1 - 2*ke*dE) + 2*ke * g2 * g2.T + H_dist
124+
H22 = hess_2 * (1.0 - 2.0 * self.kappa_e * delta_E) + \
125+
2.0 * self.kappa_e * np.dot(g2, g2.T) + \
126+
H_dist
127+
128+
# Block 12: d^2 E / dx1 dx2
129+
# = -2*ke * g1 * g2.T - H_dist
130+
H12 = -2.0 * self.kappa_e * np.dot(g1, g2.T) - H_dist
131+
132+
# Block 21: d^2 E / dx2 dx1
133+
H21 = H12.T # Symmetric
134+
135+
# Construct Full Matrix
136+
H_top = np.hstack((H11, H12))
137+
H_bot = np.hstack((H21, H22))
138+
H_total = np.vstack((H_top, H_bot))
46139

47-
return bitss_grad_1, bitss_grad_2
140+
return H_total

‎multioptpy/fileio.py‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -406,9 +406,12 @@ def print_geometry_list(self, new_geometry, element_list, electric_charge_and_mu
406406
return geometry_list
407407

408408

409-
def make_psi4_input_file(self, geometry_list, iter):#geometry_list: ang.
409+
def make_psi4_input_file(self, geometry_list, iter, path=None):#geometry_list: ang.
410410
"""structure updated geometry is saved."""
411-
file_directory = self.work_directory+"samples_"+self.NOEXT_START_FILE+"_"+str(iter)
411+
if path is not None:
412+
file_directory = os.path.join(path, f"samples_{self.NOEXT_START_FILE}_{iter}")
413+
else:
414+
file_directory = self.work_directory+"samples_"+self.NOEXT_START_FILE+"_"+str(iter)
412415
tmp_cs = ["SAMPLE"+str(iter), ""]
413416

414417

0 commit comments

Comments
 (0)