com.actelion.research.chem.phesa.EvaluableOverlap 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.List;
import com.actelion.research.chem.Coordinates;
import com.actelion.research.chem.alignment3d.transformation.ExponentialMap;
import com.actelion.research.chem.alignment3d.transformation.Quaternion;
import com.actelion.research.chem.alignment3d.transformation.RotationDerivatives;
import com.actelion.research.chem.optimization.Evaluable;
import com.actelion.research.chem.phesa.pharmacophore.pp.PPGaussian;
/**
* @author J.Wahl, February 2018
* describes the overlap function, the state (relative orientation, translation of the two molecules)
* returns the objective function and the gradient of the overlap with respect to translation and rotation
* accessed by the optimization algorithm
*/
public class EvaluableOverlap implements Evaluable {
private double ppWeight;
private PheSAAlignment shapeAlign;
private double[] transform;
private Coordinates[] cachedCoords; //coords before rotation + translation is applied, but COM at coordinate origin: Xi-p
private Coordinates[] cachedCoordsPP;
private Coordinates origCOM;
private double [][] dv0At;
private double [][] dv1At;
private double [][] dv2At;
private double [][] dv0PP;
private double [][] dv1PP;
private double [][] dv2PP;
private double[][] results;
private Coordinates[] fitAtGaussModCoords;
private Coordinates[] fitPPGaussModCoords;
public EvaluableOverlap(PheSAAlignment shapeAlign, double[] transform) {
this(shapeAlign, transform, 0.5);
}
public EvaluableOverlap(PheSAAlignment shapeAlign, double[] transform, double ppWeight) {
this.ppWeight = ppWeight;
this.shapeAlign = shapeAlign;
this.transform = transform;
this.fitAtGaussModCoords = new Coordinates[shapeAlign.getMolGauss().getAtomicGaussians().size()];
this.fitPPGaussModCoords = new Coordinates[shapeAlign.getMolGauss().getPPGaussians().size()];
this.dv0At = new double[fitAtGaussModCoords.length][3];
this.dv1At = new double[fitAtGaussModCoords.length][3];
this.dv2At = new double[fitAtGaussModCoords.length][3];
this.dv0PP = new double[fitPPGaussModCoords.length][3];
this.dv1PP = new double[fitPPGaussModCoords.length][3];
this.dv2PP = new double[fitPPGaussModCoords.length][3];
this.results = new double[shapeAlign.getRefMolGauss().getAtomicGaussians().size()][shapeAlign.getRefMolGauss().getAtomicGaussians().size()];
cachedCoords = new Coordinates[shapeAlign.getMolGauss().getAtomicGaussians().size()];
cachedCoordsPP = new Coordinates[shapeAlign.getMolGauss().getPPGaussians().size()];
origCOM = new Coordinates();
for(int i=0;i fitMolGauss) {
ExponentialMap em = new ExponentialMap(transform[0],transform[1],transform[2]);
Quaternion q = em.toQuaternion();
double[][] rotMatrix = q.getRotMatrix().getArray();
for(int k=0;k volumeGaussians = new ArrayList<>();
if(refMolGauss instanceof MolecularVolume)
volumeGaussians = ((MolecularVolume)refMolGauss).getVolumeGaussians();
value += (1.0-ppWeight)*this.getFGValueOverlap(atomGrad,refMolGauss.getAtomicGaussians(),fitMolGauss.getAtomicGaussians(),volumeGaussians,
dv0At,dv1At,dv2At,fitAtGaussModCoords);
double[] ppGrad = new double[grad.length];
value += ppWeight*this.getFGValueOverlapPP(ppGrad,refMolGauss.getPPGaussians(),fitMolGauss.getPPGaussians(),
dv0PP,dv1PP,dv2PP,fitPPGaussModCoords);
for(int i=0;i refMolGauss,List fitMolGauss,List volGaussians,
double[][] dRdv0, double[][] dRdv1, double[][] dRdv2, Coordinates[] fitGaussModCoords) {
/**
* we first calculate the partial derivatives with respect to the four elements of the quaternion q,r,s,u
* the final gradient has 7 elements, the first four elements are the gradients for the quaternion (rotation),
* the last three elements are for the translation
*/
getTransformedCoordinates(fitGaussModCoords, fitMolGauss);
this.getEMapGradient(dRdv0, dRdv1, dRdv2,cachedCoords);
/**
* derivative of ShapeOverlap with respect to the four elements of the quaternion and three elements of translation
*
*/
double totalOverlap = 0.0;
Coordinates fitCenterModCoord;
for(int i=0; i=Gaussian3D.DIST_CUTOFF)
continue;
atomOverlap = refAt.getHeight()*fitAt.getHeight()*QuickMathCalculator.getInstance().quickExp(-( refAt.getWidth() * fitAt.getWidth()* Rij2)/alphaSum) *
QuickMathCalculator.getInstance().getPrefactor(refAt.getAtomicNo(),fitAt.getAtomicNo());
if (atomOverlap>0.0) {
totalOverlap += atomOverlap;
double gradientPrefactor = atomOverlap*-2*refAt.getWidth()*fitAt.getWidth()/(refAt.getWidth()+fitAt.getWidth());
double dv0 = dRdv0[j][0]*dx+dRdv0[j][1]*dy+dRdv0[j][2]*dz;
double dv1 = dRdv1[j][0]*dx+dRdv1[j][1]*dy+dRdv1[j][2]*dz;
double dv2 = dRdv2[j][0]*dx+dRdv2[j][1]*dy+dRdv2[j][2]*dz;
grad[0] += gradientPrefactor*dv0;
grad[1] += gradientPrefactor*dv1;
grad[2] += gradientPrefactor*dv2;
grad[3] += gradientPrefactor*dx;
grad[4] += gradientPrefactor*dy;
grad[5] += gradientPrefactor*dz;
}
}
}
for(int k=0; k=Gaussian3D.DIST_CUTOFF)
continue;
atomOverlap = refVol.getRole()*refVol.getHeight()*fitAt.getHeight()*QuickMathCalculator.getInstance().quickExp(-( refVol.getWidth() * fitAt.getWidth()* Rij2)/alphaSum) *
QuickMathCalculator.getInstance().getPrefactor(refVol.getAtomicNo(),fitAt.getAtomicNo());
if (Math.abs(atomOverlap)>0.0) {
totalOverlap += atomOverlap;
double gradientPrefactor = atomOverlap*-2*refVol.getWidth()*fitAt.getWidth()/(refVol.getWidth()+fitAt.getWidth());
double dv0 = dRdv0[j][0]*dx+dRdv0[j][1]*dy+dRdv0[j][2]*dz;
double dv1 = dRdv1[j][0]*dx+dRdv1[j][1]*dy+dRdv1[j][2]*dz;
double dv2 = dRdv2[j][0]*dx+dRdv2[j][1]*dy+dRdv2[j][2]*dz;
grad[0] += gradientPrefactor*dv0;
grad[1] += gradientPrefactor*dv1;
grad[2] += gradientPrefactor*dv2;
grad[3] += gradientPrefactor*dx;
grad[4] += gradientPrefactor*dy;
grad[5] += gradientPrefactor*dz;
}
}
}
return (-1.0*totalOverlap); //the negative overlap is returned as the objective, since we minimize the objective in the optimization algorithm
}
private double getFGValueOverlapPP(double[] grad, List refMolGauss,List fitMolGauss, double[][] dRdv0, double[][] dRdv1, double[][] dRdv2, Coordinates[] fitGaussModCoords) {
getTransformedCoordinates(fitGaussModCoords,fitMolGauss);
this.getEMapGradient(dRdv0, dRdv1, dRdv2,cachedCoordsPP);
/**
* derivative of ShapeOverlap with respect to the four elements of the quaternion and three elements of translation
*
*/
double totalOverlap = 0.0;
Coordinates fitCenterModCoord;
for(int i=0; i=Gaussian3D.DIST_CUTOFF) {
continue;
}
atomOverlap = refAt.getWeight()*refAt.getHeight()*fitAt.getHeight()*QuickMathCalculator.getInstance().quickExp(-( refAt.getWidth() * fitAt.getWidth()* Rij2)/alphaSum) *
QuickMathCalculator.getInstance().getPrefactor(refAt.getAtomicNo(),fitAt.getAtomicNo());
if (atomOverlap>0.0) {
double sim = refAt.getInteractionSimilarity(fitAt);
if(sim==0.0)
continue;
atomOverlap *= sim;
totalOverlap += atomOverlap;
double gradientPrefactor = atomOverlap*-2*refAt.getWidth()*fitAt.getWidth()/(refAt.getWidth()+fitAt.getWidth());
double dv0 = dRdv0[j][0]*dx+dRdv0[j][1]*dy+dRdv0[j][2]*dz;
double dv1 = dRdv1[j][0]*dx+dRdv1[j][1]*dy+dRdv1[j][2]*dz;
double dv2 = dRdv2[j][0]*dx+dRdv2[j][1]*dy+dRdv2[j][2]*dz;
grad[0] += gradientPrefactor*dv0;
grad[1] += gradientPrefactor*dv1;
grad[2] += gradientPrefactor*dv2;
grad[3] += gradientPrefactor*dx;
grad[4] += gradientPrefactor*dy;
grad[5] += gradientPrefactor*dz;
}
}
}
return (-1.0*totalOverlap); //the negative overlap is returned as the objective, since we minimize the objective in the optimization algorithm
}
public EvaluableOverlap clone() {
return new EvaluableOverlap(this);
}
}