/usr/share/RDKit/Contrib/PBF/pbf.py is in rdkit-data 201603.5-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 | #
# Copyright (C) 2015 Greg Landrum
#
# @@ All Rights Reserved @@
# This file is part of the RDKit.
# The contents are covered by the terms of the BSD license
# which is included in the file license.txt, found at the root
# of the RDKit source tree.
from rdkit import Chem
from rdkit.Chem import AllChem
import numpy as np
from numpy import linalg
def GetBestFitPlane(pts,weights=None):
if weights is None:
wSum = len(pts)
origin = np.sum(pts,0)
origin /= wSum
sumXX=0;sumXY=0;sumXZ=0;sumYY=0;sumYZ=0;sumZZ=0;
sums = np.zeros((3,3),np.double)
for pt in pts:
dp = pt-origin
for i in range(3):
sums[i,i] += dp[i]*dp[i]
for j in range(i+1,3):
sums[i,j] += dp[i]*dp[j]
sums[j,i] += dp[i]*dp[j]
sums /= wSum
vals,vects = linalg.eigh(sums)
order = np.argsort(vals)
normal = vects[:,order[0]]
plane = np.zeros((4,),np.double)
plane[:3] = normal
plane[3] = -1 * normal.dot(origin)
return plane
def PBFRD(mol,confId=-1):
conf = mol.GetConformer(confId)
if not conf.Is3D():
return 0
pts = np.array([list(conf.GetAtomPosition(x)) for x in range(mol.GetNumAtoms())])
plane = GetBestFitPlane(pts)
denom = np.dot(plane[:3],plane[:3])
denom = denom**0.5
# add up the distance from the plane for each point:
res = 0.0
for pt in pts:
res += np.abs(pt.dot(plane[:3])+plane[3])
res /= denom
res /= len(pts)
return res
if __name__=='__main__':
suppl = Chem.SDMolSupplier('./testData/egfr.sdf',removeHs=False)
expected = open('./testData/egfr.out','r')
for m in suppl:
res = PBFRD(m);
inl = next(expected).strip().split()
expect = float(inl[1])
assert abs(res-expect)<1e-4
|