com.actelion.research.chem.phesa.ShapeVolume 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 java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.stream.Collectors;
import com.actelion.research.calc.Matrix;
import com.actelion.research.calc.SingularValueDecomposition;
import com.actelion.research.chem.Coordinates;
import com.actelion.research.chem.StereoMolecule;
import com.actelion.research.chem.alignment3d.transformation.ExponentialMap;
import com.actelion.research.chem.alignment3d.transformation.Quaternion;
import com.actelion.research.chem.alignment3d.transformation.Rotation;
import com.actelion.research.chem.alignment3d.transformation.Transformation;
import com.actelion.research.chem.alignment3d.transformation.TransformationSequence;
import com.actelion.research.chem.alignment3d.transformation.Translation;
import com.actelion.research.chem.conf.Conformer;
import com.actelion.research.chem.phesa.PheSAAlignment.axis;
import com.actelion.research.chem.phesa.pharmacophore.PharmacophoreCalculator;
import com.actelion.research.chem.phesa.pharmacophore.pp.ExitVectorPoint;
import com.actelion.research.chem.phesa.pharmacophore.pp.PPGaussian;
public class ShapeVolume {
protected List ppGaussians;
protected List atomicGaussians;
protected Coordinates com;
public ShapeVolume() {
ppGaussians = new ArrayList<>();
atomicGaussians = new ArrayList<>();
}
public ShapeVolume(ShapeVolume original) {
this.atomicGaussians = new ArrayList();
this.ppGaussians = new ArrayList();
for(AtomicGaussian ag : original.getAtomicGaussians()) {
this.atomicGaussians.add(new AtomicGaussian(ag));
}
for(PPGaussian pg : original.getPPGaussians()) {
this.ppGaussians.add(new PPGaussian(pg));
}
this.com = new Coordinates(original.com);
}
public void addPharmacophorePoint(PPGaussian ppGaussian) {
ppGaussians.add(ppGaussian);
}
public void addAtomVolume(AtomicGaussian atomGaussian) {
atomicGaussians.add(atomGaussian);
}
public void update(StereoMolecule mol) {
updateCoordinates(getAtomicGaussians(),mol.getAtomCoordinates());
updateCoordinates(getPPGaussians(),mol.getAtomCoordinates());
}
public void update(Conformer conf) {
updateCoordinates(getAtomicGaussians(),conf.getCoordinates());
updateCoordinates(getPPGaussians(),conf.getCoordinates());
}
public void transform(Transformation transform) {
transformGaussians(transform);
}
/**
* 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 Rotation preProcess(Conformer conf) {
Coordinates COM = getCOM();
if(conf!=null) {
int nrOfAtoms = conf.getSize();
for (int i=0;i toRemove = new ArrayList<>();
int i=0;
for(PPGaussian ppg : ppGaussians) {
if (ppg.getPharmacophorePoint().getFunctionalityIndex()==PharmacophoreCalculator.AROM_ID)
toRemove.add(i);
i++;
}
Collections.reverse(toRemove);
for(int index : toRemove) {
ppGaussians.remove(index);
}
}
public Matrix createCanonicalOrientation(Conformer conf) {
Matrix m = getCovarianceMatrix();
SingularValueDecomposition svd = new SingularValueDecomposition(m.getArray(),null,null);
Matrix u = new Matrix(svd.getU());
double det = u.det();
if(det<0) {
u.set(0,1,-u.get(0, 1));
u.set(1,1,-u.get(1, 1));
u.set(2,1,-u.get(2, 1));
}
Rotation rot = new Rotation(u.getArray());
if(conf!=null) {
rot.apply(conf);
update(conf);
}
else {
transformGaussians(rot);
}
if(!isCanonicalOrientation()) {
if(conf!=null) {
PheSAAlignment.rotateMolAroundAxis180(conf,axis.X);
update(conf);
}
else {
rotate180DegreeAroundAxis(axis.X);
}
if(isCanonicalOrientation()) {
u.set(0,1,-u.get(0, 1));
u.set(1,1,-u.get(1, 1));
u.set(2,1,-u.get(2, 1));
u.set(0,2,-u.get(0, 2));
u.set(1,2,-u.get(1, 2));
u.set(2,2,-u.get(2, 2));
}
else {
if(conf!=null) {
PheSAAlignment.rotateMolAroundAxis180(conf,axis.X); // rotate back
update(conf);
PheSAAlignment.rotateMolAroundAxis180(conf,axis.Y);
update(conf);
}
else {
rotate180DegreeAroundAxis(axis.X);
rotate180DegreeAroundAxis(axis.Y);
}
if(isCanonicalOrientation()) {
u.set(0,0,-u.get(0, 0));
u.set(1,0,-u.get(1, 0));
u.set(2,0,-u.get(2, 0));
u.set(0,2,-u.get(0, 2));
u.set(1,2,-u.get(1, 2));
u.set(2,2,-u.get(2, 2));
}
else {
if(conf!=null) {
PheSAAlignment.rotateMolAroundAxis180(conf,axis.Y);
update(conf);
PheSAAlignment.rotateMolAroundAxis180(conf,axis.Z);
update(conf);
}
else {
rotate180DegreeAroundAxis(axis.Y);
rotate180DegreeAroundAxis(axis.Z);
}
if(isCanonicalOrientation()) {
u.set(0,0,-u.get(0, 0));
u.set(1,0,-u.get(1, 0));
u.set(2,0,-u.get(2, 0));
u.set(0,1,-u.get(0, 1));
u.set(1,1,-u.get(1, 1));
u.set(2,1,-u.get(2, 1));
}
}
}
}
return u;
}
protected boolean isCanonicalOrientation() {
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 : atomicGaussians){
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;
}
public Coordinates getCOM() {
this.calcCOM();
return this.com;
}
public void calcCOM(){
double volume = 0.0;
double comX = 0.0;
double comY = 0.0;
double comZ = 0.0;
for(AtomicGaussian atGauss : this.atomicGaussians){
volume += atGauss.getVolume();
comX += atGauss.getCenter().x*atGauss.getVolume();
comY += atGauss.getCenter().y*atGauss.getVolume();
comZ += atGauss.getCenter().z*atGauss.getVolume();
}
comX = comX/volume;
comY = comY/volume;
comZ = comZ/volume;
this.com = new Coordinates(comX,comY,comZ);
}
protected void updateCoordinates(List gaussians, Coordinates[] coords) {
for(Gaussian3D gaussian : gaussians) {
gaussian.updateCoordinates(coords);
}
}
protected void transformGaussians(Transformation transform) {
transformGaussians(atomicGaussians,transform);
transformGaussians(ppGaussians,transform);
}
protected static void transformGaussians(List gaussians, Transformation transform) {
for (Gaussian3D gaussian : gaussians) {
gaussian.transform(transform);
}
}
protected static void rotateGaussian(List gaussians, double[][] rot) {
for (Gaussian3D gaussian : gaussians) {
Coordinates coords1 = new Coordinates(gaussian.getCenter());
coords1.rotate(rot);
gaussian.setCenter(coords1);
}
}
protected static void rotateGaussians180DegreeAroundAxis(List gaussians ,PheSAAlignment.axis a) {
for (Gaussian3D gaussian : gaussians) {
Coordinates coords1 = new Coordinates(gaussian.getCenter());
PheSAAlignment.rotateCoordsAroundAxis180(coords1, a);
gaussian.setCenter(coords1);
}
}
protected void rotate180DegreeAroundAxis(PheSAAlignment.axis a) {
rotateGaussians180DegreeAroundAxis(atomicGaussians,a);
rotateGaussians180DegreeAroundAxis(ppGaussians,a);
}
public String encode() {
StringBuilder sb = new StringBuilder();
sb.append(Integer.toString(atomicGaussians.size()));
sb.append(" ");
atomicGaussians.forEach(e -> {
sb.append(e.encode());
sb.append(" ");
});
sb.append(Integer.toString(ppGaussians.size()));
sb.append(" ");
ppGaussians.forEach(e -> {
sb.append(e.encode());
sb.append(" ");
});
return sb.toString();
}
public static ShapeVolume decode(String s) {
String[] splitString = s.split(" ");
int nrOfAtomicGaussians = Integer.decode(splitString[0].trim());
int firstIndex = 1;
int lastIndex = 1+nrOfAtomicGaussians;
List atomicGaussians = new ArrayList();
List ppGaussians = new ArrayList();
for(int i=firstIndex;i receptorVol.addAtomVolume(e));
ppGaussians.forEach(e -> receptorVol.addPharmacophorePoint(e));
return receptorVol;
}
public List getExitVectorGaussians() {
return ppGaussians.stream().filter(e -> (e.getPharmacophorePoint() instanceof ExitVectorPoint)).collect(Collectors.toList());
}
public List getPPGaussians() {
return ppGaussians;
}
public void setPPGaussians(List ppGaussians) {
this.ppGaussians = ppGaussians;
}
public List getAtomicGaussians() {
return atomicGaussians;
}
public void setAtomicGaussians(List atomicGaussians) {
this.atomicGaussians = atomicGaussians;
}
public Matrix getCovarianceMatrix() {
Matrix massMatrix = new Matrix(3,3);
double volume = 0.0;
for (AtomicGaussian ag : atomicGaussians){
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);
}
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;
}
/**
* calculate the self-overlap of the base molecule
* @return
*/
public double getSelfAtomOverlap(){
double Vtot = 0.0;
for(AtomicGaussian at:atomicGaussians){
for(AtomicGaussian at2:atomicGaussians){
Vtot += at.getVolumeOverlap(at2);
}
}
return Vtot;
}
public double getSelfPPOverlap(){
double Vtot = 0.0;
for(PPGaussian pp:ppGaussians){
for(PPGaussian pp2:ppGaussians){
Vtot+=pp.getWeight()*pp.getSimilarity(pp2)* pp.getVolumeOverlap(pp2);
}
}
return Vtot;
}
/**
* 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, ShapeVolume fitVol){
double[] result = new double[2];
ExponentialMap eMap = new ExponentialMap(transform[0],transform[1],transform[2]);
double Vtot = 0.0;
double Vvol = 0.0;
Coordinates com = fitVol.getCOM();
double[][] rotMatrix = eMap.toQuaternion().getRotMatrix().getArray();
List fitGaussians = fitVol.atomicGaussians;
Coordinates[] fitCenterModCoords = new Coordinates[fitGaussians.size()];
for(int k=0;k fitPPGaussians = fitVol.ppGaussians;
Coordinates[] fitCenterModCoords = new Coordinates[fitPPGaussians.size()];
Coordinates[] fitDirectionalityMod = new Coordinates[fitPPGaussians.size()];
for(int k=0;k