引力指標

CDK descriptor

ニュートンが発見した万有引力の法則、
すなわち、２つの物体の質量をM、m、物体間の距離をr、万有引力定数をGとして、
　物体間にはたらく力 F = G x M x m / r^2
と表せることは、古典力学のもっとも基本的な公式なひとつです。


Katritzkyら*1は、すべての原子間について、
　M x m / r^2
を計算し、これを足し合わせて引力指標（gravitational index）としました。
そして、この値と沸点との相関性を確認したのでした。

  • GRAV = Σ(M x m / r^2)

もうここまでくると直感的に原理を解釈できませんが、
先に式を示しておくとゴールが見えやすくなることもありますので、
これはこれで良いのかもしれません（妥協）。


CDKの実装では、９種類の値が算出されます。

  • heavysum : 結合が存在し、かつ双方とも重原子である原子間について計算
  • heavysumの二乗根
  • heavysumの三乗根
  • sum : 結合が存在する原子間（水素原子を含む）について計算
  • sumの二乗根
  • sumの三乗根
  • allheavysum : すべての原子間で計算
  • allheavysumの二乗根
  • allheavysumの三乗根



/* gravity.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.config.IsotopeFactory;

import org.openscience.cdk.geometry.GeometryTools;

import javax.vecmath.Point3d;

import java.util.ArrayList;
class gravity {
    public static void main(String args[]){
        if(args.length!=1){

            System.err.println("gravity <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");
            IsotopeFactory factory = null;

            double mass1;

            double mass2;

            try {

                factory = IsotopeFactory.getInstance(imol.getBuilder());

            } catch (Exception e1) {

                //logger.debug(e);

                e1.printStackTrace();

            }
            double sum = 0;

            for (int i = 0; i < imol.getBondCount(); i++) {

                IBond bond = imol.getBond(i);
                if (bond.getAtomCount() != 2) {

                    System.err.println("ERROR: Only handles 2 center bonds");

                }
                mass1 = factory.getMajorIsotope( bond.getAtom(0).getSymbol()

                                                ).getMassNumber();

                mass2 = factory.getMajorIsotope( bond.getAtom(1).getSymbol()

                                                 ).getMassNumber();
                Point3d p1 = bond.getAtom(0).getPoint3d();

                Point3d p2 = bond.getAtom(1).getPoint3d();
                double x1 = p1.x;

                double y1 = p1.y;

                double z1 = p1.z;

                double x2 = p2.x;

                double y2 = p2.y;

                double z2 = p2.z;
                double dist = (x1 - x2) * (x1 - x2) +

                    (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);

                sum += (mass1 * mass2) / dist;

            }
            // heavy atoms only

            double heavysum = 0;

            for (int i = 0; i < imol.getBondCount(); i++) {

                IBond b = imol.getBond(i);
                if (b.getAtomCount() != 2) {

                    System.err.println("ERROR: Only handles 2 center bonds");

                }
                if (b.getAtom(0).getSymbol().equals("H") ||

                    b.getAtom(1).getSymbol().equals("H")) continue;



                mass1 = factory.getMajorIsotope(b.getAtom(0).getSymbol()

                                                ).getMassNumber();

                mass2 = factory.getMajorIsotope(b.getAtom(1).getSymbol()

                                                ).getMassNumber();
                Point3d point0 = b.getAtom(0).getPoint3d();

                Point3d point1 = b.getAtom(1).getPoint3d();
                double x1 = point0.x;

                double y1 = point0.y;

                double z1 = point0.z;

                double x2 = point1.x;

                double y2 = point1.y;

                double z2 = point1.z;
                double dist = (x1 - x2) * (x1 - x2) +

                    (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);

                heavysum += (mass1 * mass2) / dist;

            }
            // all pairs

            ArrayList<Integer> x = new ArrayList<Integer>();

            for (int i = 0; i < imol.getAtomCount(); i++) {

                if (!imol.getAtom(i).getSymbol().equals("H")) x.add(i);

            }

            int npair = x.size() * (x.size() - 1) / 2;

            pair[] p = new pair[npair];

            for (int i = 0; i < npair; i++) p[i] = new pair();

            int pcount = 0;

            for (int i = 0; i < x.size() - 1; i++) {

                for (int j = i + 1; j < x.size(); j++) {

                    int present = 0;

                    int a = x.get(i);

                    int b = x.get(j);

                    for (int k = 0; k < pcount; k++) {

                        if ((p[k].x == a && p[k].y == b) ||

                            (p[k].y == a && p[k].x == b)) present = 1;

                    }

                    if (present == 1) continue;

                    p[pcount].x = a;

                    p[pcount].y = b;

                    pcount += 1;

                }

            }

            double allheavysum = 0;

            for (pair aP : p) {

                int atomNumber1 = aP.x;

                int atomNumber2 = aP.y;
                mass1 = factory.getMajorIsotope(imol.getAtom(atomNumber1).getSymbol()

                                                ).getMassNumber();

                mass2 = factory.getMajorIsotope(imol.getAtom(atomNumber2).getSymbol()

                                                ).getMassNumber();
                double x1 = imol.getAtom(atomNumber1).getPoint3d().x;

                double y1 = imol.getAtom(atomNumber1).getPoint3d().y;

                double z1 = imol.getAtom(atomNumber1).getPoint3d().z;

                double x2 = imol.getAtom(atomNumber2).getPoint3d().x;

                double y2 = imol.getAtom(atomNumber2).getPoint3d().y;

                double z2 = imol.getAtom(atomNumber2).getPoint3d().z;
                double dist = (x1 - x2) * (x1 - x2) +

                    (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2);

                allheavysum += (mass1 * mass2) / dist;

            }

            System.out.printf("%s", id);

            System.out.printf("\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\n",

                              heavysum, Math.sqrt(heavysum), Math.pow(heavysum, 1.0/3.0),

                              sum, Math.sqrt(sum), Math.pow(sum, 1.0/3.0),

                              allheavysum, Math.sqrt(allheavysum),

                              Math.pow(allheavysum, 1.0 / 3.0) );

        }

    }

    static private class pair {

        int x, y;

        public pair() {

            x = 0;

            y = 0;

        }

    }

}

*1: Katritzky et al. Correlation of Boiling Points With Molecular Structure. 1. A Training Set of 298 Diverse Organics and a Test Set of 9 Simple Inorganics. J. Phys. Chem. (1996) 100, 10400-10407.