diff --git a/arc/species/converter.py b/arc/species/converter.py index 6db652dae0..b8f844b995 100644 --- a/arc/species/converter.py +++ b/arc/species/converter.py @@ -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) @@ -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 diff --git a/arc/species/converter_test.py b/arc/species/converter_test.py index 9b302e7272..b0f32306e1 100644 --- a/arc/species/converter_test.py +++ b/arc/species/converter_test.py @@ -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