com.actelion.research.chem.phesa.PheSAAlignment Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of openchemlib Show documentation
Show all versions of openchemlib Show documentation
Open Source Chemistry Library
package com.actelion.research.chem.phesa;
import com.actelion.research.chem.Coordinates;
import com.actelion.research.chem.StereoMolecule;
import com.actelion.research.chem.conf.Conformer;
import com.actelion.research.chem.phesa.pharmacophore.PPGaussian;
import com.actelion.research.calc.Matrix;
import com.actelion.research.calc.SingularValueDecomposition;
import java.util.Arrays;
import java.util.stream.DoubleStream;
import java.util.stream.IntStream;
import java.util.ArrayList;
/**
* @version: 1.0, February 2018
* Author: J. Wahl
* this class provides functionalities to calculate the overlap between two molecules
*
*/
public class PheSAAlignment {
private MolecularVolume refMolGauss;
private MolecularVolume molGauss;
private double ppWeight;
public enum axis {X,Y,Z};
public PheSAAlignment(StereoMolecule refMol, StereoMolecule mol,double ppWeight) {
this.ppWeight = ppWeight;
this.refMolGauss = new MolecularVolume(refMol);
this.molGauss = new MolecularVolume(mol);
}
public PheSAAlignment(StereoMolecule refMol, StereoMolecule mol) {
this(refMol,mol,0.5);
}
public PheSAAlignment(MolecularVolume refMolGauss, MolecularVolume molGauss) {
this(refMolGauss,molGauss,0.5);
}
public PheSAAlignment(MolecularVolume refMolGauss, MolecularVolume molGauss,double ppWeight) {
this.ppWeight = ppWeight;
this.refMolGauss= refMolGauss;
this.molGauss = molGauss;
}
public MolecularVolume getRefMolGauss() {
return refMolGauss;
}
public MolecularVolume getMolGauss() {
return molGauss;
}
/**
* Move COM of the molecular volume to the origin of the lab-frame and orient molecules so that their principal moments
* of inertia coincide with the 3 axis of the coordinate system
* @param mol
* @param molVol
*/
public static Matrix preProcess(Conformer conf, MolecularVolume molVol) {
Coordinates COM = molVol.getCOM();
int nrOfAtoms = conf.getSize();
for (int i=0;i {
Coordinates coords = conf.getCoordinates(i);
coords.y = -coords.y;
coords.z = -coords.z;
});
}
else if (a == axis.Y) {
IntStream.range(0,conf.getSize()).forEach(i -> {
Coordinates coords = conf.getCoordinates(i);
coords.x = -coords.x;
coords.z = -coords.z;
});
}
else {
IntStream.range(0,conf.getSize()).forEach(i -> {
Coordinates coords = conf.getCoordinates(i);
coords.x = -coords.x;
coords.y = -coords.y;
});
}
}
private static Matrix getCovarianceMatrix(MolecularVolume molGauss) {
Matrix massMatrix = new Matrix(3,3);
double volume = 0.0;
for (AtomicGaussian ag : molGauss.getAtomicGaussians()){
volume += ag.getVolume();
double value = ag.getVolume()*ag.getCenter().x*ag.getCenter().x;
massMatrix.addToElement(0,0,value);
value = ag.getVolume()*ag.getCenter().x*ag.getCenter().y;
massMatrix.addToElement(0,1,value);
value = ag.getVolume()*ag.getCenter().x*ag.getCenter().z;
massMatrix.addToElement(0,2,value);
value = ag.getVolume()*ag.getCenter().y*ag.getCenter().y;
massMatrix.addToElement(1,1,value);
value = ag.getVolume()*ag.getCenter().y*ag.getCenter().z;
massMatrix.addToElement(1,2,value);
value = ag.getVolume()*ag.getCenter().z*ag.getCenter().z;
massMatrix.addToElement(2,2,value);
}
for (VolumeGaussian vg : molGauss.getVolumeGaussians()){
volume += vg.getRole()*vg.getVolume();
double value = vg.getRole()*vg.getVolume()*vg.getCenter().x*vg.getCenter().x;
massMatrix.addToElement(0,0,value);
value = vg.getRole()*vg.getVolume()*vg.getCenter().x*vg.getCenter().y;
massMatrix.addToElement(0,1,value);
value = vg.getRole()*vg.getVolume()*vg.getCenter().x*vg.getCenter().z;
massMatrix.addToElement(0,2,value);
value = vg.getRole()*vg.getVolume()*vg.getCenter().y*vg.getCenter().y;
massMatrix.addToElement(1,1,value);
value = vg.getRole()*vg.getVolume()*vg.getCenter().y*vg.getCenter().z;
massMatrix.addToElement(1,2,value);
value = vg.getRole()*vg.getVolume()*vg.getCenter().z*vg.getCenter().z;
massMatrix.addToElement(2,2,value);
}
massMatrix.set(0,0,massMatrix.get(0,0)/volume);
massMatrix.set(0,1,massMatrix.get(0,1)/volume);
massMatrix.set(0,2,massMatrix.get(0,2)/volume);
massMatrix.set(1,1,massMatrix.get(1,1)/volume);
massMatrix.set(1,2,massMatrix.get(1,2)/volume);
massMatrix.set(2,2,massMatrix.get(2,2)/volume);
massMatrix.set(1,0,massMatrix.get(0,1));
massMatrix.set(2,0,massMatrix.get(0,2));
massMatrix.set(2,1,massMatrix.get(1,2));
return massMatrix;
}
private static boolean isCanonicalOrientation(MolecularVolume molGauss) {
double xxPos = 0;
double xxNeg = 0;
double yyPos = 0;
double yyNeg = 0;
int nXPos = 0;
int nXNeg = 0;
int nYPos = 0;
int nYNeg = 0;
for (AtomicGaussian ag : molGauss.getAtomicGaussians()){
double x = ag.center.x;
double y = ag.center.y;
if(x>0) {
xxPos += x*x;
nXPos++;
}
else {
xxNeg += x*x;
nXNeg++;
}
if(y>0) {
yyPos += y*y;
nYPos++;
}
else {
yyNeg += y*y;
nYNeg++;
}
}
xxPos/=nXPos;
yyPos/=nYPos;
xxNeg/=nXNeg;
yyNeg/=nYNeg;
if(xxPos>xxNeg && yyPos>yyNeg)
return true;
else
return false;
}
/**
* .
* generate initial orientations of the molecule:
* mode1: 4 orientations: initial orientation and 180 degree rotation about each axis
* mode2: mode1 and 90 degree rotations about each axis
* a transformation vector consists of 7 elements: the first 4 elements form a Quaternion and describe the rotation
* the last three elements are the translation vector
* @param mode
* @return
*/
public static double[][] initialTransform(int mode) {
double c = 0.707106781;
switch(mode){
case 1:
double[][] transforms1 = {{1.0,0.0,0.0,0.0,0.0,0.0,0.0},{0.0,1.0,0.0,0.0,0.0,0.0,0.0},{0.0,0.0,1.0,0.0,0.0,0.0,0.0},
{0.0,0.0,0.0,1.0,0.0,0.0,0.0}};
return transforms1;
case 2:
double[][] transforms2 = {{1,0,0,0,0,0,0},{0,1,0,0,0,0,0},{0,0,1,0,0,0,0},
{0,0,0,1,0,0,0},
{c,c,0,0,0,0,0},
{c,0,c,0,0,0,0},
{c,0,0,c,0,0,0},
{-0.5,0.5,0.5,-0.5,0,0,0},
{0.5,-0.5,0.5,-0.5,0,0,0},
{0.5,0.5,0.5,-0.5,0,0,0},
{0.5,-0.5,-0.5,-0.5,0,0,0},
{0.5,0.5,-0.5,-0.5,0,0,0}
};
return transforms2;
default:
double [][] transform = {{1.0,0.0,0.0,0.0,0.0,0.0,0.0}};
return transform;
}
}
/**
* calculate the Overlap of the two molecular volumes as a function a transform vector that is applied to the query molecule
* overlap Volume of two molecular Volumes formulated as a summed overlap of atomic Gaussians
* taken from Grant, Gallardo, Pickup, Journal of Computational Chemistry, 17, 1653-1666, 1996
* returns a double[2]: the first double is the total overlap, whereas the second value is the specific
* contribution of additional volume gaussians (inclusion, exclusion)
* @param transform
* @return
*/
public double[] getTotalAtomOverlap(double[] transform){
double[] result = new double[2];
Quaternion quat = new Quaternion(transform[0],transform[1],transform[2],transform[3]);
double Vtot = 0.0;
double Vvol = 0.0;
double[][] rotMatrix = quat.getRotMatrix().getArray();
ArrayList atomicGaussians = molGauss.getAtomicGaussians();
Coordinates[] fitCenterModCoords = new Coordinates[atomicGaussians.size()];
double normFactor = 1/(transform[0]*transform[0]+transform[1]*transform[1]+transform[2]*transform[2]+transform[3]*transform[3]);
for(int k=0;k ppGaussians = molGauss.getPPGaussians();
Coordinates[] fitCenterModCoords = new Coordinates[ppGaussians.size()];
Coordinates[] fitDirectionalityMod = new Coordinates[ppGaussians.size()];
double normFactor = 1/(transform[0]*transform[0]+transform[1]*transform[1]+transform[2]*transform[2]+transform[3]*transform[3]);
for(int k=0;k g.getWeight()).sum();
ppSimilarity*=correctionFactor;
if(ppSimilarity>1.0) //can happen because of weights
ppSimilarity = 1.0f;
double[] result = getTotalAtomOverlap(bestTransform);
atomOverlap = result[0];
double additionalVolOverlap = result[1];
double atomSimilarity = atomOverlap/(Oaa+Obb-atomOverlap);
if(atomSimilarity>1.0) //can happen because of weights
atomSimilarity = 1.0f;
volSimilarity = (additionalVolOverlap/atomOverlap);
similarity = (1.0-ppWeight)*atomSimilarity + ppWeight*ppSimilarity;
if (similarity>maxSimilarity) {
maxSimilarity = similarity;
maxVolSimilarity = volSimilarity;
maxPPSimilarity = ppSimilarity;
alignment = bestTransform;
}
}
if(maxSimilarity>1.0) // can happen because of manually placed inclusion spheres
maxSimilarity = 1.0;
return DoubleStream.concat(Arrays.stream(new double[] {maxSimilarity,maxPPSimilarity,maxVolSimilarity}), Arrays.stream(alignment)).toArray();
}
public static void rotateMol(Conformer conf, Matrix rotMat) {
int nrOfAtoms = conf.getSize();
for (int i=0;i © 2015 - 2025 Weber Informatics LLC | Privacy Policy