Petitjean形状指標
「直径 = 2 x 半径」
と小学生のときに習いましたが、これは円や球に通用する公式であって、
その他の物体では成立しないことがあります。化学構造も然りです。
まずは、二次元記述子計算のときと同じく、化学構造をグラフと考えてみます。
すべての原子について、その原子から最も遠い原子までの距離（パス長）を計算したときに、
その最大値を「直径」、最小値を「半径」と定義します。
つまり、いつか解説したeccentricity ηを使って、
- 半径R = min( η )
- 直径D = max( η )
と定式化できます。
(注: eccentricity = その原子から最も遠い原子までの距離。See this page.)
このとき
- Petitjean形状指標 I = (D - R) / D
により、分子の形状を数値化することができます*1。
この指標は、三次元座標空間にも拡張できます。
つまり、パス長の代わりにユークリッド距離を使えば、
上と同様にして、Petitjean形状指標の三次元版が計算できるというわけです。
以下は、CDKによる実装です。上記の二次元版（topoShape）と三次元版（geomShape）を計算します。後者のみ水素込みで計算しています。
/* petit_shape.java */import java.io.*;
import org.openscience.cdk.interfaces.*;
import org.openscience.cdk.io.iterator.IteratingMDLReader;
import org.openscience.cdk.DefaultChemObjectBuilder;import org.openscience.cdk.geometry.GeometryTools;
import org.openscience.cdk.graph.PathTools;
import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
import javax.vecmath.Point3d;class petit_shape {
public static void main(String args[]){
if(args.length!=1){
System.err.println("petit_shape <sd-file>");
System.exit(1);
}
FileInputStream fis = null;
IteratingMDLReader isr = null;
try{
fis = new FileInputStream(new File(args[0]) );
isr = new IteratingMDLReader(fis, DefaultChemObjectBuilder.getInstance() );
} catch (Exception e) {
e.printStackTrace();
}while( isr.hasNext() ){
IMolecule imol = (IMolecule)isr.next();
String id = (String)imol.getProperty("cdk:Title");
IAtomContainer container = (IAtomContainer) imol;
IAtomContainer local = AtomContainerManipulator.removeHydrogens(container);int tradius = PathTools.getMolecularGraphRadius(local);
int tdiameter = PathTools.getMolecularGraphDiameter(local);
double topoShape = (double) (tdiameter - tradius) / (double) tradius;double geomShape = 0.0;
// get the 3D distance matrix
if (GeometryTools.has3DCoordinates(container)) {
int natom = container.getAtomCount();
double[][] distanceMatrix = new double[natom][natom];
for (int i = 0; i < natom; i++) {
for (int j = 0; j < natom; j++) {
if (i == j) {
distanceMatrix[i][j] = 0.0;
continue;
}Point3d a = container.getAtom(i).getPoint3d();
Point3d b = container.getAtom(j).getPoint3d();
distanceMatrix[i][j] = Math.sqrt((a.x - b.x) * (a.x - b.x) +
(a.y - b.y) * (a.y - b.y) +
(a.z - b.z) * (a.z - b.z));
}
}double gradius = 999999;
double gdiameter = -999999;
double[] geta = new double[natom];
for (int i = 0; i < natom; i++) {
double max = -99999;
for (int j = 0; j < natom; j++) {
if (distanceMatrix[i][j] > max) max = distanceMatrix[i][j];
}
geta[i] = max;
}
for (int i = 0; i < natom; i++) {
if (geta[i] < gradius) gradius = geta[i];
if (geta[i] > gdiameter) gdiameter = geta[i];
}
geomShape = (gdiameter - gradius) / gradius;
}
System.out.printf("%s\t%.3f\t%.3f\n", id, topoShape, geomShape);
}
}
}
