Skip to content

Commit 800af3f

Browse files
ss0832ss0832
authored
Add files via upload
1 parent 10662f9 commit 800af3f

2 files changed

Lines changed: 539 additions & 0 deletions

File tree

‎multioptpy/Optimizer/mf_rsirfo.py‎

Lines changed: 266 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,266 @@
1+
import numpy as np
2+
from multioptpy.Optimizer.hessian_update import ModelHessianUpdate
3+
from multioptpy.Optimizer.block_hessian_update import BlockHessianUpdate
4+
from multioptpy.Utils.calc_tools import Calculationtools
5+
6+
# Import the original RSIRFO class for inheritance
7+
from multioptpy.Optimizer.rsirfo import RSIRFO
8+
from multioptpy.Optimizer.mode_following import ModeFollowing
9+
10+
11+
class MF_RSIRFO(RSIRFO):
12+
"""
13+
Mode-Following RS-I-RFO Optimizer.
14+
15+
References:
16+
[1] Banerjee et al., Phys. Chem., 89, 52-57 (1985)
17+
[2] Heyden et al., J. Chem. Phys., 123, 224101 (2005)
18+
[3] Baker, J. Comput. Chem., 7, 385-395 (1986)
19+
[4] Besalú and Bofill, Theor. Chem. Acc., 100, 265-274 (1998)
20+
21+
This code is made based on the below codes.
22+
1, https://github.com/eljost/pysisyphus/blob/master/pysisyphus/tsoptimizers/TSHessianOptimizer.py
23+
2, https://github.com/eljost/pysisyphus/blob/master/pysisyphus/tsoptimizers/RSIRFOptimizer.py
24+
25+
26+
Extended 'method' string support:
27+
"method_name : target_index : ema<val> : grad<val>"
28+
29+
Examples:
30+
"block_fsb" -> Default (Index 0, EMA=1.0 if adaptive, Grad=0.0)
31+
"block_fsb:1" -> Track Mode 1
32+
"block_fsb:0:ema0.5" -> Track Mode 0 with EMA alpha=0.5
33+
"block_fsb:ema0.1:grad0.3" -> Track Mode 0, EMA=0.1, Gradient Bias=0.3
34+
"""
35+
def __init__(self, **config):
36+
# 1. Parse 'method' string for advanced configs
37+
raw_method_str = config.get("method", "auto")
38+
39+
# Initial Defaults
40+
update_method = raw_method_str
41+
target_mode_index = 0
42+
43+
# Check config for fallback defaults
44+
is_adaptive_config = config.get("adaptive_mode_following", True)
45+
46+
# Placeholders for parsed values
47+
parsed_update_rate = None
48+
parsed_gradient_weight = 0.0
49+
50+
# Parse logic
51+
if ":" in raw_method_str:
52+
parts = raw_method_str.split(":")
53+
update_method = parts[0].strip() # First part is always method name
54+
55+
for part in parts[1:]:
56+
part = part.strip().lower()
57+
if not part: continue
58+
59+
if part.isdigit():
60+
# "1", "2" -> Target Index
61+
target_mode_index = int(part)
62+
63+
elif part.startswith("ema"):
64+
# "ema0.5" -> Update Rate
65+
try:
66+
val = float(part[3:])
67+
parsed_update_rate = val
68+
except ValueError:
69+
print(f"Warning: Invalid ema value in '{part}'. Ignoring.")
70+
71+
elif part.startswith("grad"):
72+
# "grad0.5" -> Gradient Weight
73+
try:
74+
val = float(part[4:])
75+
parsed_gradient_weight = val
76+
except ValueError:
77+
print(f"Warning: Invalid grad value in '{part}'. Ignoring.")
78+
79+
# Resolve Update Rate (EMA) and Adaptive Flag
80+
if parsed_update_rate is not None:
81+
# Explicit string config overrides config dict
82+
update_rate = parsed_update_rate
83+
# If ema > 0, we must enable adaptive mode
84+
adaptive = (update_rate > 1e-12)
85+
else:
86+
# Fallback to config dict or defaults
87+
# "If ema not specified: 0 if static, 1 if adaptive"
88+
adaptive = is_adaptive_config
89+
update_rate = 1.0 if adaptive else 0.0
90+
91+
# Resolve Gradient Weight
92+
gradient_weight = parsed_gradient_weight
93+
94+
# Update config for parent class
95+
config['method'] = update_method
96+
97+
# Initialize parent RSIRFO
98+
super().__init__(**config)
99+
self.hessian_update_method = update_method
100+
101+
self.use_mode_following = config.get("use_mode_following", True)
102+
103+
# Other configs
104+
use_hungarian = config.get("use_hungarian", True)
105+
element_list = config.get("element_list", None)
106+
107+
# Initialize Mode Following with resolved parameters
108+
self.mode_follower = ModeFollowing(
109+
self.saddle_order,
110+
atoms=element_list,
111+
initial_target_index=target_mode_index,
112+
adaptive=adaptive,
113+
update_rate=update_rate,
114+
use_hungarian=use_hungarian,
115+
gradient_weight=gradient_weight,
116+
debug_mode=config.get("debug_mode", False)
117+
)
118+
119+
if self.display_flag:
120+
print(f"MF-RS-I-RFO Initialized:")
121+
print(f" - Update Method: {self.hessian_update_method}")
122+
print(f" - Target Index: {target_mode_index}")
123+
print(f" - Mode Following: Adaptive={adaptive} (EMA Rate={update_rate})")
124+
print(f" - Gradient Bias: {gradient_weight}")
125+
print(f" - Matching: {'Hungarian' if use_hungarian else 'Greedy'}")
126+
print(f" - Mass-Weighted: {'Yes' if element_list else 'No'}")
127+
128+
def run(self, geom_num_list, B_g, pre_B_g=[], pre_geom=[], B_e=0.0, pre_B_e=0.0, pre_move_vector=[], initial_geom_num_list=[], g=[], pre_g=[]):
129+
"""
130+
Execute one step of RS-I-RFO with Advanced Mode Following.
131+
"""
132+
self.log(f"\n{'='*50}\nMF-RS-I-RFO Iteration {self.iteration}\n{'='*50}", force=True)
133+
134+
if self.Initialization:
135+
self.prev_eigvec_min = None
136+
self.prev_eigvec_size = None
137+
self.predicted_energy_changes = []
138+
self.actual_energy_changes = []
139+
self.prev_geometry = None
140+
self.prev_gradient = None
141+
self.prev_energy = None
142+
self.converged = False
143+
self.iteration = 0
144+
self.Initialization = False
145+
146+
if self.hessian is None:
147+
raise ValueError("Hessian matrix must be set")
148+
149+
if self.prev_geometry is not None and self.prev_gradient is not None and len(pre_g) > 0 and len(pre_geom) > 0:
150+
self.update_hessian(geom_num_list, g, pre_geom, pre_g)
151+
152+
gradient_norm = np.linalg.norm(B_g)
153+
self.log(f"Gradient norm: {gradient_norm:.6f}", force=True)
154+
155+
if gradient_norm < self.gradient_norm_threshold:
156+
self.log(f"Converged: Gradient norm", force=True)
157+
self.converged = True
158+
159+
if self.actual_energy_changes and abs(self.actual_energy_changes[-1]) < self.energy_change_threshold:
160+
self.log(f"Converged: Energy change", force=True)
161+
self.converged = True
162+
163+
current_energy = B_e
164+
gradient = np.asarray(B_g).ravel()
165+
166+
tmp_hess = self.hessian
167+
if self.bias_hessian is not None:
168+
H_base = tmp_hess + self.bias_hessian
169+
else:
170+
H_base = tmp_hess
171+
172+
H = Calculationtools().project_out_hess_tr_and_rot_for_coord(
173+
H_base, geom_num_list.reshape(-1, 3), geom_num_list.reshape(-1, 3), False
174+
)
175+
H = 0.5 * (H + H.T)
176+
177+
eigvals, eigvecs = self.compute_eigendecomposition_with_shift(H)
178+
self.check_hessian_conditioning(eigvals)
179+
180+
# =========================================================================
181+
# Mode Following: Identify Targets
182+
# =========================================================================
183+
target_indices = []
184+
185+
if self.saddle_order > 0:
186+
if self.use_mode_following:
187+
if self.iteration == 0:
188+
self.log(f"Init Mode Following...")
189+
self.mode_follower.set_references(eigvecs, eigvals)
190+
# For iter 0, use start indices directly
191+
start = self.mode_follower.target_offset
192+
target_indices = list(range(start, start + self.saddle_order))
193+
else:
194+
# Find matching modes (pass gradient for optional bias)
195+
target_indices = self.mode_follower.get_matched_indices(
196+
eigvecs, eigvals, current_gradient=gradient
197+
)
198+
else:
199+
target_indices = list(range(self.saddle_order))
200+
201+
# =========================================================================
202+
# Trust Radius
203+
# =========================================================================
204+
if self.iteration > 0 and self.prev_energy is not None:
205+
actual_change = B_e - self.prev_energy
206+
if len(self.actual_energy_changes) >= 3:
207+
self.actual_energy_changes.pop(0)
208+
self.actual_energy_changes.append(actual_change)
209+
210+
if self.predicted_energy_changes:
211+
# Use curvature of the TRACKED mode
212+
min_eigval_for_tr = eigvals[0]
213+
if target_indices and target_indices[0] < len(eigvals):
214+
min_eigval_for_tr = eigvals[target_indices[0]]
215+
216+
self.adjust_trust_radius(
217+
actual_change,
218+
self.predicted_energy_changes[-1],
219+
min_eigval_for_tr,
220+
gradient_norm
221+
)
222+
223+
# =========================================================================
224+
# Image Surface Construction
225+
# =========================================================================
226+
P = np.eye(gradient.size)
227+
228+
for idx in target_indices:
229+
if idx < len(eigvals) and np.abs(eigvals[idx]) > 1e-10:
230+
trans_vec = eigvecs[:, idx]
231+
if self.NEB_mode:
232+
P -= np.outer(trans_vec, trans_vec)
233+
else:
234+
P -= 2 * np.outer(trans_vec, trans_vec)
235+
236+
H_star = np.dot(P, H)
237+
H_star = 0.5 * (H_star + H_star.T)
238+
grad_star = np.dot(P, gradient)
239+
240+
eigvals_star, eigvecs_star = self.compute_eigendecomposition_with_shift(H_star)
241+
eigvals_star, eigvecs_star = self.filter_small_eigvals(eigvals_star, eigvecs_star)
242+
243+
current_eigvec_size = eigvecs_star.shape[1]
244+
if self.prev_eigvec_size is not None and self.prev_eigvec_size != current_eigvec_size:
245+
self.prev_eigvec_min = None
246+
self.prev_eigvec_size = current_eigvec_size
247+
248+
move_vector = self.get_rs_step(eigvals_star, eigvecs_star, grad_star)
249+
250+
predicted_change = self.rfo_model(gradient, H, move_vector)
251+
252+
if len(self.predicted_energy_changes) >= 3:
253+
self.predicted_energy_changes.pop(0)
254+
self.predicted_energy_changes.append(predicted_change)
255+
256+
self.log(f"Predicted energy change: {predicted_change:.6f}", force=True)
257+
258+
if self.actual_energy_changes and len(self.predicted_energy_changes) > 1:
259+
self.evaluate_step_quality()
260+
261+
self.prev_geometry = geom_num_list
262+
self.prev_gradient = B_g
263+
self.prev_energy = current_energy
264+
self.iteration += 1
265+
266+
return -1 * move_vector.reshape(-1, 1)

0 commit comments

Comments
 (0)