Skip to content

Commit

Permalink
make the input rdkit
Browse files Browse the repository at this point in the history
  • Loading branch information
JintaoWu98 committed Dec 3, 2024
1 parent 7d46e62 commit 2a26437
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 34 deletions.
49 changes: 15 additions & 34 deletions arc/species/converter.py
Original file line number Diff line number Diff line change
Expand Up @@ -1852,43 +1852,34 @@ def set_rdkit_dihedrals(conf, rd_mol, torsion, deg_increment=None, deg_abs=None)
return new_xyz


def set_rdkit_ring_dihedrals(xyz, ring_head, ring_tail, torsions, dihedrals):
def set_rdkit_ring_dihedrals(rd_mol, ring_head, ring_tail, torsions, dihedrals):
"""
A helper function for setting dihedral angles using RDKit.
Either ``deg_increment`` or ``deg_abs`` must be specified.
A helper function for setting dihedral angles in a ring using RDKit.
Args:
conf: The RDKit conformer with the current xyz information.
rd_mol: The respective RDKit molecule.
torsion (list, tuple): The 0-indexed atom indices of the four atoms defining the torsion.
deg_increment (float, optional): The required dihedral increment in degrees.
deg_abs (float, optional): The required dihedral in degrees.
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.
Raises:
ConverterError: If the dihedral cannot be set.
"""

# xyz = conformers[22]['xyz']
s_mol, b_mol = molecules_from_xyz(xyz)
mol = b_mol if b_mol is not None else s_mol
conf, rd_mol = rdkit_conf_from_mol(mol, xyz)


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()
# torsions = [(0, 1, 2, 3), (1, 2, 3, 4), (2, 3, 4, 5)]
# dihedrals = [305.736462037467,
# 54.263604783044116,
# 305.7364054482706]
for torsion, dihedral in zip(torsions, dihedrals):
torsion_0_indexed = [tor - 0 for tor in torsion]
xyz_dihedrals = set_rdkit_dihedrals(conf_mod, rd_mol_mod, torsion_0_indexed, deg_abs=dihedral)
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)
Expand All @@ -1901,16 +1892,6 @@ def set_rdkit_ring_dihedrals(xyz, ring_head, ring_tail, torsions, dihedrals):
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)



# rdMT.SetDihedralDeg(conf, torsion[0], torsion[1], torsion[2], torsion[3], deg_abs)
# coords = list()
# symbols = list()
# for i, atom in enumerate(list(rd_mol.GetAtoms())):
# coords.append([conf.GetAtomPosition(i).x, conf.GetAtomPosition(i).y, conf.GetAtomPosition(i).z])
# symbols.append(atom.GetSymbol())
# new_xyz = xyz_from_data(coords=coords, symbols=symbols)
return new_xyz


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 2a26437

Please sign in to comment.