Skip to content

Commit 79c2aca

Browse files
alongdalongd
authored
Correct over-squaring in internal coordinate objective function for adding an atom to XYZ (#802)
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.
2 parents 72b90ae + ea972cf commit 79c2aca

File tree

2 files changed

+16
-12
lines changed

2 files changed

+16
-12
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

‎arc/species/converter_test.py‎

Lines changed: 8 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1082,7 +1082,6 @@ def test_zmat_from_xyz(self):
10821082
'A_9_0_1': 107.32624053955078, 'D_9_0_1_8': 243.79214484231346,
10831083
'R_6|7_4|4': 1.089653491973877},
10841084
'map': {0: 1, 1: 0, 2: 2, 3: 3, 4: 4, 5: 5, 6: 6, 7: 7, 8: 8, 9: 9}}
1085-
self.assertEqual(zmat_1, expected_zmat_1)
10861085
self.assertTrue(_compare_zmats(zmat_1, expected_zmat_1))
10871086

10881087
zmat_2 = converter.zmat_from_xyz(xyz=self.xyz6['str'],
@@ -4669,7 +4668,7 @@ def test_add_atom_to_xyz_using_internal_coords(self):
46694668
(-3.17488405, 0.21042224, 1.24128129), (-3.53487445, 1.21797729, 1.25310144),
46704669
(-3.56603419, -0.31800758, 2.08547528), (-1.63624745, 0.21425655, 1.30595531),
46714670
(-1.24509731, 0.74268637, 0.46176133), (-1.3183353, 0.6960459, 2.20690531),
4672-
(-1.0413732051724125, -1.4526805980775839, 1.2873769860752802))}
4671+
(-1.0398666203453495, -1.4521333553792743, 1.2864271004043641))}
46734672
self.assertTrue(almost_equal_coords_lists(new_xyz_1, expected_xyz_1))
46744673
self.assertAlmostEqual(calculate_param(coords=new_xyz_1['coords'], atoms=[7, 10]), 1.77, places=2)
46754674
self.assertAlmostEqual(calculate_param(coords=new_xyz_1['coords'], atoms=[4, 7, 10]), 109.5, places=1)
@@ -4685,9 +4684,9 @@ def test_add_atom_to_xyz_using_internal_coords(self):
46854684
d_value=-60.0,
46864685
opt_method='BFGS',
46874686
)
4688-
self.assertAlmostEqual(new_xyz_2['coords'][-1][0], -1.065095539786487, places=3)
4689-
self.assertAlmostEqual(new_xyz_2['coords'][-1][1], -1.4737067126385224, places=3)
4690-
self.assertAlmostEqual(new_xyz_2['coords'][-1][2], 1.2912758496528232, places=3)
4687+
self.assertAlmostEqual(new_xyz_2['coords'][-1][0], -1.0407054376010356, places=2)
4688+
self.assertAlmostEqual(new_xyz_2['coords'][-1][1], -1.4428236405325816, places=2)
4689+
self.assertAlmostEqual(new_xyz_2['coords'][-1][2], 1.2912758496528232, places=2)
46914690
self.assertAlmostEqual(calculate_param(coords=new_xyz_1['coords'], atoms=[4, 10]), 2.70, places=1)
46924691
self.assertAlmostEqual(calculate_param(coords=new_xyz_1['coords'], atoms=[4, 0, 10]), 61.46, places=0)
46934692
self.assertAlmostEqual(calculate_param(coords=new_xyz_1['coords'], atoms=[0, 4, 7, 10]), 300, places=1)
@@ -4715,7 +4714,7 @@ def test_add_atom_to_xyz_using_internal_coords(self):
47154714
(0.35773087, -1.66017412, -0.9786309), (-0.45608483, -1.87500387, -1.86208833),
47164715
(-1.82486467, -0.81522856, 0.14629516), (-1.06962462, 0.60119223, 0.90442455),
47174716
(-1.14968688, 0.45844916, -0.88969505), (1.33643417, -2.15859899, -0.90083808),
4718-
(1.4820209463152687, -2.3772302636106577, 0.30367261748030894))}
4717+
(1.4828150949768655, -2.377060198658591, 0.30306573909787815))}
47194718
self.assertTrue(almost_equal_coords(new_xyz_3, expected_xyz))
47204719
self.assertAlmostEqual(calculate_param(coords=new_xyz_3['coords'], atoms=[2, 8]), 1.85, places=2)
47214720
self.assertAlmostEqual(calculate_param(coords=new_xyz_3['coords'], atoms=[1, 2, 8]), 77.4, places=1)
@@ -4744,9 +4743,9 @@ def test_add_atom_to_xyz_using_internal_coords(self):
47444743
a_value=77.4,
47454744
d_value=140,
47464745
)
4747-
self.assertAlmostEqual(new_xyz_4['coords'][-1][0], 1.8986166253764283)
4748-
self.assertAlmostEqual(new_xyz_4['coords'][-1][1], 0.9987236974936107)
4749-
self.assertAlmostEqual(new_xyz_4['coords'][-1][2], 0.8154061174277444)
4746+
self.assertAlmostEqual(new_xyz_4['coords'][-1][0], 1.8986166253764283, places=2)
4747+
self.assertAlmostEqual(new_xyz_4['coords'][-1][1], 0.9987236974936107, places=2)
4748+
self.assertAlmostEqual(new_xyz_4['coords'][-1][2], 0.8154061174277444, places=2)
47504749
self.assertAlmostEqual(calculate_param(coords=new_xyz_4['coords'], atoms=[1, 14]), 1.85, places=2)
47514750
self.assertAlmostEqual(calculate_param(coords=new_xyz_4['coords'], atoms=[3, 1, 14]), 77.4, places=1)
47524751
self.assertAlmostEqual(calculate_param(coords=new_xyz_4['coords'], atoms=[2, 0, 1, 14]), 140, places=1)

0 commit comments

Comments
 (0)