Skip to content

Commit

Permalink
Fix H_Abstraction TS search heuristics for CH3 + RH <=> CH4 + R (#595)
Browse files Browse the repository at this point in the history 
In the case of the reaction `CH3 + RH <=> CH4 + R`, we get that the ZMat
constructed for CH4 has parameters which rely on the abstracted H atom.
This is an edge case, usually higher angles and dihedral angles are not
defined w.r.t H atoms where possible, but here there's no choice. The
issue reported by @naddeu was that removing the H atom from the ZMat
resulted in illegal ZMat parameters (a parameter that refers twice to a
certain atom).
This PR fixes this issue and adds a bunch of tests

fixes #638
  • Loading branch information
@alongd
alongd committed Apr 15, 2023
2 parents fac0c4e + dcf3c61 commit ebe2642
Show file tree
Hide file tree
Showing 5 changed files with 292 additions and 92 deletions.
17 changes: 17 additions & 0 deletions arc/job/adapters/ts/heuristics_test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -901,6 +901,23 @@ def test_heuristics_for_h_abstraction(self):
heuristics_1.execute_incore()
self.assertEqual(len(rxn1.ts_species.ts_guesses), 6)

ch3 = ARCSpecies(label='CH3', smiles='[CH3]')
cdcdc = ARCSpecies(label='C=C=C', smiles='C=C=C')
ch4 = ARCSpecies(label='CH4', smiles='C')
chdcdc = ARCSpecies(label='CH=C=C', smiles='[CH]=C=C')
rxn1 = ARCReaction(r_species=[ch3, cdcdc], p_species=[ch4, chdcdc])
rxn1.determine_family(rmg_database=self.rmgdb)
self.assertEqual(rxn1.family.label, 'H_Abstraction')
heuristics_1 = HeuristicsAdapter(job_type='tsg',
reactions=[rxn1],
testing=True,
project='test',
project_directory=os.path.join(ARC_PATH, 'arc', 'testing', 'heuristics'),
dihedral_increment=120,
)
heuristics_1.execute_incore()
self.assertEqual(len(rxn1.ts_species.ts_guesses), 3)

def test_keeping_atom_order_in_ts(self):
"""Test that the generated TS has the same atom order as in the reactants"""
# reactant_reversed, products_reversed = False, False
Expand Down
26 changes: 26 additions & 0 deletions arc/species/vectors.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -232,6 +232,32 @@ def calculate_dihedral_angle(coords: Union[list, tuple, dict],
return get_dihedral(v1, v2, v3, units=units)


def calculate_param(coords: Union[list, tuple, dict],
atoms: list,
index: int = 0,
) -> float:
"""
Calculate a distance / angle / dihedral angle parameter.
Default units (deg for angles) are used.
Args:
coords (list, tuple, dict): The array-format or tuple-format coordinates, or the xyz dict.
atoms (list): The 2 atoms to calculate defining the vector for which the length will be calculated.
index (int, optional): Whether ``atoms`` is 0-indexed or 1-indexed (values are 0 or 1).
Returns: float
The calculated parameter.
"""
if len(atoms) == 2:
return calculate_distance(coords=coords, atoms=atoms, index=index)
elif len(atoms) == 3:
return calculate_angle(coords=coords, atoms=atoms, index=index)
elif len(atoms) == 4:
return calculate_dihedral_angle(coords=coords, torsion=atoms, index=index)
else:
raise ValueError(f'Expected 2, 3, or 4 indices in coords, got {len(coords)}: {coords}')


def unit_vector(vector: List[float]) -> List[float]:
"""
Calculate a unit vector in the same direction as the input vector.
Expand Down
75 changes: 37 additions & 38 deletions arc/species/vectors_test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,15 @@ def setUpClass(cls):
A method that is run before all unit tests in this class.
"""
cls.maxDiff = None
cls.propene = converter.str_to_xyz("""C 1.22905000 -0.16449200 0.00000000
C -0.13529200 0.45314000 0.00000000
C -1.27957200 -0.21983000 0.00000000
H 1.17363000 -1.25551200 0.00000000
H 1.79909600 0.15138400 0.87934300
H 1.79909600 0.15138400 -0.87934300
H -0.16831500 1.54137600 0.00000000
H -2.23664600 0.28960500 0.00000000
H -1.29848800 -1.30626200 0.00000000""")

def test_get_normal(self):
"""Test calculating a normal vector"""
Expand Down Expand Up @@ -64,6 +73,11 @@ def test_get_angle(self):
theta = vectors.get_angle(v1, v2)
self.assertAlmostEqual(theta, math.pi)

v1 = [-0.9110001674469675, 0.41240600991486, 0.0]
v2 = [0.8619789448920263, 0.5069440956654351, 0.0]
theta = vectors.get_angle(v1, v2)
self.assertAlmostEqual(theta, 2.1848632513137702, 5)

def test_get_dihedral(self):
"""Test calculating a dihedral angle from vectors"""
v1, v2, v3 = [1., 1., 0.], [0., 1., 1.], [1., 0., 1.] # test floats
Expand All @@ -89,20 +103,11 @@ def test_get_dihedral(self):

def test_calculate_distance(self):
"""Test calculating a distance between two atoms"""
propene = converter.str_to_xyz("""C 1.22905000 -0.16449200 0.00000000
C -0.13529200 0.45314000 0.00000000
C -1.27957200 -0.21983000 0.00000000
H 1.17363000 -1.25551200 0.00000000
H 1.79909600 0.15138400 0.87934300
H 1.79909600 0.15138400 -0.87934300
H -0.16831500 1.54137600 0.00000000
H -2.23664600 0.28960500 0.00000000
H -1.29848800 -1.30626200 0.00000000""")
distance = vectors.calculate_distance(coords=propene['coords'], atoms=[1, 4], index=1)
distance = vectors.calculate_distance(coords=self.propene['coords'], atoms=[1, 4], index=1)
self.assertAlmostEqual(distance, 1.092426698)
distance = vectors.calculate_distance(coords=propene['coords'], atoms=[1, 2], index=1)
distance = vectors.calculate_distance(coords=self.propene['coords'], atoms=[1, 2], index=1)
self.assertAlmostEqual(distance, 1.49763087)
distance = vectors.calculate_distance(coords=propene['coords'], atoms=[2, 3], index=1)
distance = vectors.calculate_distance(coords=self.propene['coords'], atoms=[2, 3], index=1)
self.assertAlmostEqual(distance, 1.32750337)

def test_calculate_angle(self):
Expand All @@ -124,36 +129,17 @@ def test_calculate_angle(self):
O -1.40465894 -0.03095532 0.00000000""")
angle = vectors.calculate_angle(coords=fake_co2['coords'], atoms=[0, 1, 2], index=0, units='degs')
self.assertEqual(angle, 0.0)

propene = converter.str_to_xyz("""C 1.22905000 -0.16449200 0.00000000
C -0.13529200 0.45314000 0.00000000
C -1.27957200 -0.21983000 0.00000000
H 1.17363000 -1.25551200 0.00000000
H 1.79909600 0.15138400 0.87934300
H 1.79909600 0.15138400 -0.87934300
H -0.16831500 1.54137600 0.00000000
H -2.23664600 0.28960500 0.00000000
H -1.29848800 -1.30626200 0.00000000""")
angle = vectors.calculate_angle(coords=propene['coords'], atoms=[8, 3, 9], index=1, units='degs')
angle = vectors.calculate_angle(coords=self.propene['coords'], atoms=[8, 3, 9], index=1, units='degs')
self.assertAlmostEqual(angle, 117.02817, 4)
angle = vectors.calculate_angle(coords=propene['coords'], atoms=[9, 3, 8], index=1, units='degs')
angle = vectors.calculate_angle(coords=self.propene['coords'], atoms=[9, 3, 8], index=1, units='degs')
self.assertAlmostEqual(angle, 117.02817, 4)
angle = vectors.calculate_angle(coords=propene['coords'], atoms=[1, 2, 3], index=1, units='degs')
angle = vectors.calculate_angle(coords=self.propene['coords'], atoms=[1, 2, 3], index=1, units='degs')
self.assertAlmostEqual(angle, 125.18344, 4)
angle = vectors.calculate_angle(coords=propene['coords'], atoms=[5, 1, 2], index=1, units='degs')
angle = vectors.calculate_angle(coords=self.propene['coords'], atoms=[5, 1, 2], index=1, units='degs')
self.assertAlmostEqual(angle, 110.82078, 4)

def test_calculate_dihedral_angle(self):
"""Test calculating a dihedral angle from xyz coordinates"""
propene = converter.str_to_xyz("""C 1.22905000 -0.16449200 0.00000000
C -0.13529200 0.45314000 0.00000000
C -1.27957200 -0.21983000 0.00000000
H 1.17363000 -1.25551200 0.00000000
H 1.79909600 0.15138400 0.87934300
H 1.79909600 0.15138400 -0.87934300
H -0.16831500 1.54137600 0.00000000
H -2.23664600 0.28960500 0.00000000
H -1.29848800 -1.30626200 0.00000000""")
hydrazine = converter.str_to_xyz("""N 0.70683700 -0.07371000 -0.21400700
N -0.70683700 0.07371000 -0.21400700
H 1.11984200 0.81113900 -0.47587600
Expand Down Expand Up @@ -233,11 +219,11 @@ def test_calculate_dihedral_angle(self):
H 2.134 -1.659 -3.429
H 2.951 -3.078 -4.102""")

dihedral0 = vectors.calculate_dihedral_angle(coords=propene['coords'], torsion=[9, 3, 2, 7], index=1)
dihedral1 = vectors.calculate_dihedral_angle(coords=propene['coords'], torsion=[5, 1, 2, 7], index=1)
dihedral0 = vectors.calculate_dihedral_angle(coords=self.propene['coords'], torsion=[9, 3, 2, 7], index=1)
dihedral1 = vectors.calculate_dihedral_angle(coords=self.propene['coords'], torsion=[5, 1, 2, 7], index=1)
self.assertAlmostEqual(dihedral0, 180, 2)
self.assertAlmostEqual(dihedral1, 59.26447, 2)
dihedral2 = vectors.calculate_dihedral_angle(coords=propene['coords'], torsion=[8, 2, 1, 6], index=0)
dihedral2 = vectors.calculate_dihedral_angle(coords=self.propene['coords'], torsion=[8, 2, 1, 6], index=0)
self.assertEqual(dihedral0, dihedral2)

dihedral3 = vectors.calculate_dihedral_angle(coords=hydrazine['coords'], torsion=[3, 1, 2, 5], index=1)
Expand All @@ -246,6 +232,19 @@ def test_calculate_dihedral_angle(self):
dihedral2 = vectors.calculate_dihedral_angle(coords=cj_11974['coords'], torsion=[15, 18, 19, 20], index=1)
self.assertAlmostEqual(dihedral2, 308.04758, 2)

def test_calculate_param(self):
"""Test calculating a parameter"""
distance = vectors.calculate_param(coords=self.propene['coords'], atoms=[0, 1])
self.assertAlmostEqual(distance, 1.497630871004034)
angle = vectors.calculate_param(coords=self.propene['coords'], atoms=[0, 1, 2])
self.assertAlmostEqual(angle, 125.18344391469404)
dihedral = vectors.calculate_param(coords=self.propene['coords'], atoms=[0, 1, 2, 6])
self.assertAlmostEqual(dihedral, 180.0)
with self.assertRaises(ValueError):
vectors.calculate_param(coords=self.propene['coords'], atoms=[0])
with self.assertRaises(ValueError):
vectors.calculate_param(coords=self.propene['coords'], atoms=[0, 1, 2, 3, 4])

def test_unit_vector(self):
"""Test calculating a unit vector"""
v1 = [1, 0, 0]
Expand Down

0 comments on commit ebe2642

Please sign in to comment.