Skip to content

Commit 708e750

Browse files
committed
Correct over-squaring in internal coordinate objective function
The objective function for adding a new atom previously led to over-penalization of angular errors (effectively minimizing the error to the fourth power). This was due to the angular constraint functions (g_theta and g_phi) already returning squared errors, which were then squared again in the final objective function. This commit refactors the constraint functions to return the linear error or linear error magnitude: - `angle_constraint`: Now returns the linear difference (calc_angle - target_angle). - `dihedral_constraint`: Now returns the magnitude of the sine/cosine error vector (square root of the squared sum). The main objective function remains structurally the same (sum of squared, scaled terms), but now correctly implements a standard **relative squared error** minimization.
1 parent da68824 commit 708e750

File tree

1 file changed

+8
-3
lines changed

1 file changed

+8
-3
lines changed

arc/species/converter.py

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,9 @@
4444
ob.obErrorLog.SetOutputLevel(0)
4545
logger = get_logger()
4646

47+
DIST_PRECISION = 0.01 # Angstrom
48+
ANGL_PRECISION = 0.1 # rad (for both bond angle and dihedral)
49+
4750

4851
def str_to_xyz(xyz_str: str,
4952
project_directory: Optional[str] = None,
@@ -2074,7 +2077,7 @@ def calculate_errors(result_coords):
20742077

20752078
def meets_precision(result_coords):
20762079
r_error, a_error, d_error = calculate_errors(result_coords)
2077-
return r_error < 0.01 and a_error < 0.1 and d_error < 0.1
2080+
return r_error < DIST_PRECISION and a_error < ANGL_PRECISION and d_error < ANGL_PRECISION
20782081

20792082
guess_functions = [
20802083
generate_initial_guess_r_a,
@@ -2262,7 +2265,9 @@ def angle_eq(x, y, z):
22622265
cross_product_length = np.linalg.norm(np.cross(BA, BX))
22632266
dot_product = np.dot(BA, BX)
22642267
calc_angle = math.atan2(cross_product_length, dot_product)
2265-
return (calc_angle - target_angle) ** 2
2268+
angle_diff = calc_angle - target_angle
2269+
wrapped_diff = np.arctan2(np.sin(angle_diff), np.cos(angle_diff)) # wrapped to the range (-pi, pi]
2270+
return wrapped_diff
22662271

22672272
return angle_eq
22682273

@@ -2310,7 +2315,7 @@ def dihedral_eq(x, y, z):
23102315
BC_norm = np.linalg.norm(BC)
23112316
cos_calc = np.dot(N1, N2) / (N1_norm * N2_norm)
23122317
sin_calc = np.dot(BC, np.cross(N1, N2)) / (BC_norm * N1_norm * N2_norm)
2313-
return (cos_calc - cos_d) ** 2 + (sin_calc - sin_d) ** 2
2318+
return np.sqrt((cos_calc - cos_d) ** 2 + (sin_calc - sin_d) ** 2)
23142319

23152320
return dihedral_eq
23162321

0 commit comments

Comments
 (0)