pub_soft/rdkit
3
0
Fork 0
mirror of https://github.com/rdkit/rdkit.git synced 2026-06-04 21:54:27 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
master
rdkit/Code/GraphMol/MolTransforms/Wrap/testMolTransforms.py
Ric 880a8e5725 Reformat Python code for 2023.03 release (#6294) 
* run yapf

* run isort

---------

Co-authored-by: Greg Landrum <greg.landrum@gmail.com>
2023-04-28 06:53:56 +02:00

291 lines
11 KiB
Python
Raw Permalink Blame History

import math
import os
import sys
import unittest
import numpy as np
from rdkit import Chem, DataStructs, RDConfig
from rdkit.Chem import rdMolTransforms as rdmt
from rdkit.Geometry import rdGeometry as geom
def feq(v1, v2, tol=1.0e-4):
return abs(v1 - v2) < tol
def ptEq(pt, tpt, tol=0.0001):
pt -= tpt
return feq(pt.Length(), 0.0)
class TestCase(unittest.TestCase):
def setUp(self):
pass
def test1Canonicalization(self):
mol = Chem.MolFromSmiles("C")
conf = Chem.Conformer(1)
conf.SetAtomPosition(0, (4.0, 5.0, 6.0))
mol.AddConformer(conf, 1)
conf = mol.GetConformer()
pt = rdmt.ComputeCentroid(conf)
self.assertTrue(ptEq(pt, geom.Point3D(4.0, 5.0, 6.0)))
fileN = os.path.join(RDConfig.RDBaseDir, 'Code', 'GraphMol', 'MolTransforms', 'test_data',
'1oir.mol')
m = Chem.MolFromMolFile(fileN)
cpt = rdmt.ComputeCentroid(m.GetConformer())
trans = rdmt.ComputeCanonicalTransform(m.GetConformer(), cpt)
trans2 = rdmt.ComputeCanonicalTransform(m.GetConformer())
for i in range(4):
for j in range(4):
self.assertTrue(feq(trans[i, j], trans2[i, j]))
rdmt.TransformConformer(m.GetConformer(), trans2)
m2 = Chem.MolFromMolFile(fileN)
rdmt.CanonicalizeConformer(m2.GetConformer())
nats = m.GetNumAtoms()
cnf1 = m.GetConformer()
cnf2 = m2.GetConformer()
for i in range(nats):
p1 = list(cnf1.GetAtomPosition(i))
p2 = list(cnf2.GetAtomPosition(i))
self.assertTrue(feq(p1[0], p2[0]))
self.assertTrue(feq(p1[1], p2[1]))
self.assertTrue(feq(p1[2], p2[2]))
m3 = Chem.MolFromMolFile(fileN)
rdmt.CanonicalizeMol(m3)
cnf1 = m.GetConformer()
cnf2 = m3.GetConformer()
for i in range(nats):
p1 = list(cnf1.GetAtomPosition(i))
p2 = list(cnf2.GetAtomPosition(i))
self.assertTrue(feq(p1[0], p2[0]))
self.assertTrue(feq(p1[1], p2[1]))
self.assertTrue(feq(p1[2], p2[2]))
def testComputePrincipalAxesAndMoments(self):
if (not hasattr(rdmt, 'ComputePrincipalAxesAndMoments')):
return
molBlock = '''\
RDKit 3D
4 3 0 0 0 0 0 0 0 0999 V2000
-0.1888 1.3224 -0.2048 F 0 0 0 0 0 0 0 0 0 0 0 0
-0.0790 0.1671 0.6070 C 0 0 1 0 0 0 0 0 0 0 0 0
-1.3083 -0.8939 -0.2236 Cl 0 0 0 0 0 0 0 0 0 0 0 0
1.5761 -0.5956 -0.1786 Br 0 0 0 0 0 0 0 0 0 0 0 0
2 1 1 6
2 3 1 0
2 4 1 0
M END
'''
axesRef = ((-0.9997, -0.0246, 0.0009), (0.0246, -0.9981, 0.0559), (0.0004, -0.0559, -0.9984))
momentsRef = (3.4220, 4.7230, 7.1757)
m = Chem.MolFromMolBlock(molBlock)
axes, moments = rdmt.ComputePrincipalAxesAndMoments(m.GetConformer())
self.assertIsNotNone(axes)
self.assertIsNotNone(moments)
for y in range(3):
for x in range(3):
self.assertAlmostEqual(axes[y][x], axesRef[y][x], 3)
self.assertAlmostEqual(moments[y], momentsRef[y], 3)
failed = False
try:
axes, moments = rdmt.ComputePrincipalAxesAndMoments(m.GetConformer(), weights=(0.5, 0.5))
except Exception:
failed = True
self.assertTrue(failed)
axesWeightedRef = ((-0.9998, -0.0114, -0.0189), (-0.0153, 0.9744, 0.2245), (0.0158, 0.2247,
-0.9743))
momentsWeightedRef = (0.5496, 1.5559, 1.9361)
axesWeighted, momentsWeighted = rdmt.ComputePrincipalAxesAndMoments(
m.GetConformer(), weights=(0.1, 0.2, 0.3, 0.4))
self.assertIsNotNone(axesWeighted)
self.assertIsNotNone(momentsWeighted)
for y in range(3):
for x in range(3):
self.assertAlmostEqual(axesWeighted[y][x], axesWeightedRef[y][x], 3)
self.assertAlmostEqual(momentsWeighted[y], momentsWeightedRef[y], 3)
def testComputePrincipalAxesAndMomentsFromGyrationMatrix(self):
if (not hasattr(rdmt, 'ComputePrincipalAxesAndMomentsFromGyrationMatrix')):
return
molBlock = '''\
RDKit 3D
4 3 0 0 0 0 0 0 0 0999 V2000
-0.1888 1.3224 -0.2048 F 0 0 0 0 0 0 0 0 0 0 0 0
-0.0790 0.1671 0.6070 C 0 0 1 0 0 0 0 0 0 0 0 0
-1.3083 -0.8939 -0.2236 Cl 0 0 0 0 0 0 0 0 0 0 0 0
1.5761 -0.5956 -0.1786 Br 0 0 0 0 0 0 0 0 0 0 0 0
2 1 1 6
2 3 1 0
2 4 1 0
M END
'''
axesRef = ((-0.0009, -0.0246, 0.9997), (-0.0559, -0.9981, -0.0246), (0.9984, -0.0559, -0.0004))
momentsRef = (0.1212, 0.7343, 1.0596)
m = Chem.MolFromMolBlock(molBlock)
axes, moments = rdmt.ComputePrincipalAxesAndMomentsFromGyrationMatrix(m.GetConformer())
self.assertIsNotNone(axes)
self.assertIsNotNone(moments)
for y in range(3):
for x in range(3):
self.assertAlmostEqual(axes[y][x], axesRef[y][x], 3)
self.assertAlmostEqual(moments[y], momentsRef[y], 3)
failed = False
try:
axes, moments = rdmt.ComputePrincipalAxesAndMomentsFromGyrationMatrix(
m.GetConformer(), weights=(0.5, 0.5))
except Exception:
failed = True
self.assertTrue(failed)
axesWeightedRef = ((0.0189, -0.0114, 0.9998), (-0.2245, 0.9744, 0.0153), (0.9743, 0.2247,
-0.0158))
momentsWeightedRef = (0.0847, 0.4649, 1.4712)
axesWeighted, momentsWeighted = rdmt.ComputePrincipalAxesAndMomentsFromGyrationMatrix(
m.GetConformer(), weights=(0.1, 0.2, 0.3, 0.4))
self.assertIsNotNone(axesWeighted)
self.assertIsNotNone(momentsWeighted)
for y in range(3):
for x in range(3):
self.assertAlmostEqual(abs(axesWeighted[y][x]), abs(axesWeightedRef[y][x]), 3)
self.assertAlmostEqual(abs(momentsWeighted[y]), abs(momentsWeightedRef[y]), 3)
def testGetSetBondLength(self):
file = os.path.join(RDConfig.RDBaseDir, 'Code', 'GraphMol', 'MolTransforms', 'test_data',
'3-cyclohexylpyridine.mol')
m = Chem.MolFromMolFile(file, True, False)
conf = m.GetConformer()
dist = rdmt.GetBondLength(conf, 0, 19)
self.assertAlmostEqual(dist, 1.36, 2)
rdmt.SetBondLength(conf, 0, 19, 2.5)
dist = rdmt.GetBondLength(conf, 0, 19)
self.assertAlmostEqual(dist, 2.5, 1)
rdmt.SetBondLength(conf, 19, 0, 3.0)
dist = rdmt.GetBondLength(conf, 0, 19)
self.assertAlmostEqual(dist, 3.0, 1)
def testGetSetAngle(self):
file = os.path.join(RDConfig.RDBaseDir, 'Code', 'GraphMol', 'MolTransforms', 'test_data',
'3-cyclohexylpyridine.mol')
m = Chem.MolFromMolFile(file, True, False)
conf = m.GetConformer()
angle = rdmt.GetAngleDeg(conf, 0, 19, 21)
self.assertAlmostEqual(angle, 109.7, 1)
rdmt.SetAngleDeg(conf, 0, 19, 21, 125.0)
angle = rdmt.GetAngleDeg(conf, 0, 19, 21)
self.assertAlmostEqual(angle, 125.0, 1)
rdmt.SetAngleRad(conf, 21, 19, 0, math.pi / 2.)
angle = rdmt.GetAngleRad(conf, 0, 19, 21)
self.assertAlmostEqual(angle, math.pi / 2., 1)
angle = rdmt.GetAngleDeg(conf, 0, 19, 21)
self.assertAlmostEqual(angle, 90.0, 1)
def testGetSetDihedral(self):
file = os.path.join(RDConfig.RDBaseDir, 'Code', 'GraphMol', 'MolTransforms', 'test_data',
'3-cyclohexylpyridine.mol')
m = Chem.MolFromMolFile(file, True, False)
conf = m.GetConformer()
dihedral = rdmt.GetDihedralDeg(conf, 0, 19, 21, 24)
self.assertAlmostEqual(dihedral, 176.05, 2)
rdmt.SetDihedralDeg(conf, 8, 0, 19, 21, 65.0)
dihedral = rdmt.GetDihedralDeg(conf, 8, 0, 19, 21)
self.assertAlmostEqual(dihedral, 65.0, 1)
rdmt.SetDihedralDeg(conf, 8, 0, 19, 21, -130.0)
dihedral = rdmt.GetDihedralDeg(conf, 8, 0, 19, 21)
self.assertAlmostEqual(dihedral, -130.0, 1)
rdmt.SetDihedralRad(conf, 21, 19, 0, 8, -2. / 3. * math.pi)
dihedral = rdmt.GetDihedralRad(conf, 8, 0, 19, 21)
self.assertAlmostEqual(dihedral, -2. / 3. * math.pi, 1)
dihedral = rdmt.GetDihedralDeg(conf, 8, 0, 19, 21)
self.assertAlmostEqual(dihedral, -120.0, 1)
def testGetSetDihedralThroughTripleBond(self):
file = os.path.join(RDConfig.RDBaseDir, 'Code', 'GraphMol', 'MolTransforms', 'test_data',
'github1262_2.mol')
m = Chem.MolFromMolFile(file, True, False)
conf = m.GetConformer()
rdmt.SetDihedralDeg(conf, 6, 1, 2, 9, 0.0)
dihedral = rdmt.GetDihedralDeg(conf, 6, 1, 2, 9)
self.assertAlmostEqual(dihedral, 0.0, 1)
dist = rdmt.GetBondLength(conf, 6, 9)
rdmt.SetDihedralDeg(conf, 6, 1, 2, 9, 120.0)
dihedral = rdmt.GetDihedralDeg(conf, 6, 1, 2, 9)
self.assertAlmostEqual(dihedral, 120.0, 1)
dist2 = rdmt.GetBondLength(conf, 6, 7)
self.assertAlmostEqual(dist, dist2, 1)
rdmt.SetDihedralDeg(conf, 6, 1, 2, 9, 180.0)
dihedral = rdmt.GetDihedralDeg(conf, 6, 1, 2, 9)
self.assertAlmostEqual(dihedral, 180.0, 1)
dist3 = rdmt.GetBondLength(conf, 6, 9)
self.assertNotAlmostEqual(dist, dist3, 1)
exceptionRaised = False
try:
rdmt.SetDihedralDeg(conf, 6, 0, 3, 9, 0.0)
except ValueError:
exceptionRaised = True
self.assertTrue(exceptionRaised)
def testEigen3CanonicalTransformAgainstNumpy(self):
def canonicalize_conf_rdkit(mol, conf_id=-1):
mol = Chem.Mol(mol)
conf = mol.GetConformer(conf_id)
ctd = rdmt.ComputeCentroid(conf)
canon_trans = rdmt.ComputeCanonicalTransform(conf, ctd)
rdmt.TransformConformer(conf, canon_trans)
return mol
def canonicalize_conf_numpy(mol, conf_id=-1):
mol = Chem.Mol(mol)
conf = mol.GetConformer(conf_id)
pos = conf.GetPositions()
ctd = pos.mean(axis=0)
trans_pos = pos - ctd
cov_mat = np.cov(trans_pos, bias=1, rowvar=False) * conf.GetNumAtoms()
eigval, eigvect = np.linalg.eig(cov_mat)
eigval_sorted = sorted(enumerate(eigval), key=lambda x: x[1], reverse=True)
eigvect_sorted = [
eigvect[:, i] * (1.0 if eigvect[:, i].sum() > 0.0 else -1.0) for i, _ in eigval_sorted
]
canon_trans = np.array([
[*eigvect_sorted[0], 0.],
[*eigvect_sorted[1], 0.],
[*eigvect_sorted[2], 0.],
[0., 0., 0., 1.],
], dtype=np.double)
canon_trans[0:3, 3] = np.array([*ctd, 1.]).dot(canon_trans)[:3]
rdmt.TransformConformer(conf, canon_trans)
return mol
file = os.path.join(RDConfig.RDBaseDir, 'Code', 'GraphMol', 'MolTransforms', 'test_data',
'github4302.sdf')
built_against_eigen3 = hasattr(rdmt, 'ComputePrincipalAxesAndMomentsFromGyrationMatrix')
with Chem.SDMolSupplier(file) as suppl:
for mol in suppl:
mol_orig = Chem.Mol(mol)
rdkit_canon = canonicalize_conf_rdkit(mol)
numpy_canon = canonicalize_conf_numpy(mol)
for i in range(mol.GetNumAtoms()):
orig_coord = mol_orig.GetConformer().GetAtomPosition(i)
rdkit_coord = rdkit_canon.GetConformer().GetAtomPosition(i)
numpy_coord = numpy_canon.GetConformer().GetAtomPosition(i)
self.assertLess(orig_coord.Distance(numpy_coord), 1.e-4)
if built_against_eigen3:
self.assertLess(orig_coord.Distance(rdkit_coord), 1.e-4)
if __name__ == "__main__":
unittest.main()
Reference in New Issue View Git Blame Copy Permalink