Skip to content

Commit

Permalink
define set_rdkit_ring_dihedrals and make test
Browse files Browse the repository at this point in the history
  • Loading branch information
JintaoWu98 committed Sep 29, 2024
1 parent fb16f4a commit ecb1262
Show file tree
Hide file tree
Showing 2 changed files with 80 additions and 0 deletions.
43 changes: 43 additions & 0 deletions arc/species/converter.py
Original file line number Diff line number Diff line change
Expand Up @@ -1852,6 +1852,49 @@ def set_rdkit_dihedrals(conf, rd_mol, torsion, deg_increment=None, deg_abs=None)
return new_xyz


def set_rdkit_ring_dihedrals(rd_mol, ring_head, ring_tail, torsions, dihedrals):
"""
A helper function for setting dihedral angles in a ring using RDKit.
Args:
rd_mol: The respective RDKit molecule.
ring_head: The first atom index of the ring(0-indexed).
ring_tail: The last atom index of the ring(0-indexed).
torsions: A list of torsions, each corresponding to a dihedral.
dihedrals: A list of dihedral angles in degrees, each corresponding to a torsion.
Example of a 6-membered ring:
ring_head = 0
ring_tail = 5
torsions = [(0, 1, 2, 3), (1, 2, 3, 4), (2, 3, 4, 5)]
dihedrals = [30, 300, 30]
Returns:
dict: The xyz with the new dihedral, ordered according to the map.
"""

rd_mol_mod = Chem.RWMol(rd_mol)
rd_mol_mod.RemoveBond(ring_head, ring_tail)
Chem.SanitizeMol(rd_mol_mod)
conf_mod = rd_mol_mod.GetConformer()
for torsion, dihedral in zip(torsions, dihedrals):
torsion_0_indexed = [tor - 0 for tor in torsion]
set_rdkit_dihedrals(conf_mod, rd_mol_mod, torsion_0_indexed, deg_abs=dihedral)

rd_mol_mod.AddBond(ring_head, ring_tail, Chem.BondType.SINGLE)
Chem.SanitizeMol(rd_mol_mod)
Chem.AllChem.MMFFOptimizeMolecule(rd_mol_mod)
conf_mod = rd_mol_mod.GetConformer()

coords = list()
symbols = list()
for i, atom in enumerate(list(rd_mol_mod.GetAtoms())):
coords.append([conf_mod.GetAtomPosition(i).x, conf_mod.GetAtomPosition(i).y, conf_mod.GetAtomPosition(i).z])
symbols.append(atom.GetSymbol())
new_xyz = xyz_from_data(coords=coords, symbols=symbols)
return new_xyz


def check_isomorphism(mol1, mol2, filter_structures=True, convert_to_single_bonds=False):
"""
Convert ``mol1`` and ``mol2`` to RMG Species objects, and generate resonance structures.
Expand Down
37 changes: 37 additions & 0 deletions arc/species/converter_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -4066,6 +4066,43 @@ def test_set_rdkit_dihedrals(self):
H 2.16336803 0.09985803 0.03295192"""
self.assertEqual(converter.xyz_to_str(new_xyz4), expected_xyz4)

def test_set_rdkit_ring_dihedrals(self):
"""Test setting the dihedral angles of an RDKit ring molecule"""
xyz_original = """C 1.17528959 0.88689342 -0.09425887
C -0.23165323 1.40815606 -0.37444021
C -1.28915380 0.60789983 0.38119602
C -1.17528947 -0.88689346 0.09425817
C 0.23165279 -1.40815571 0.37444068
C 1.28915350 -0.60789979 -0.38119592
H 1.90063595 1.43610053 -0.70501194
H 1.43190556 1.07695419 0.95523181
H -0.29672067 2.46483469 -0.09164586
H -0.43309289 1.35229514 -1.45133707
H -2.28822258 0.96189799 0.10312701
H -1.17664390 0.78164432 1.45848873
H -1.43190253 -1.07695588 -0.95523291
H -1.90063264 -1.43610606 0.70501042
H 0.29671416 -2.46483479 0.09164748
H 0.43309139 -1.35229454 1.45133785
H 1.17664469 -0.78164459 -1.45848883
H 2.28822409 -0.96189136 -0.10312655"""
xyz_original = converter.str_to_xyz(xyz_original)

ring_head = 0
ring_tail = 5
torsions = [(0, 1, 2, 3), (1, 2, 3, 4), (2, 3, 4, 5)]
dihedrals = [29.167577928701704, 299.8936870462789, 29.167577208303104]

s_mol, b_mol = converter.molecules_from_xyz(xyz_original)
mol = b_mol if b_mol is not None else s_mol
_, rd_mol = converter.rdkit_conf_from_mol(mol, xyz_original)

xyz_final = converter.set_rdkit_ring_dihedrals(rd_mol, ring_head, ring_tail, torsions, dihedrals)

self.assertAlmostEqual(calculate_dihedral_angle(xyz_final,[0,1,2,3]), 29.167577928701704, 2)
self.assertAlmostEqual(calculate_dihedral_angle(xyz_final,[1,2,3,4]), 299.8936870462789, 2)
self.assertAlmostEqual(calculate_dihedral_angle(xyz_final,[2,3,4,5]), 29.167577208303104, 2)

def test_get_center_of_mass(self):
"""Test calculating the center of mass for coordinates"""
xyz = """O 1.28706525 0.52121353 0.04219198
Expand Down

0 comments on commit ecb1262

Please sign in to comment.