com.actelion.research.chem.phesaflex.EvaluableFlexibleOverlap 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.phesaflex;
import java.util.Map;
import com.actelion.research.chem.Coordinates;
import com.actelion.research.chem.StereoMolecule;
import com.actelion.research.chem.alignment3d.PheSAAlignmentOptimizer.PheSASetting;
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.RotationDerivatives;
import com.actelion.research.chem.alignment3d.transformation.TransformationSequence;
import com.actelion.research.chem.alignment3d.transformation.Translation;
import com.actelion.research.chem.conf.BondRotationHelper;
import com.actelion.research.chem.conf.Conformer;
import com.actelion.research.chem.conf.TorsionDB;
import com.actelion.research.chem.forcefield.mmff.ForceFieldMMFF94;
import com.actelion.research.chem.optimization.Evaluable;
import com.actelion.research.chem.phesa.AtomicGaussian;
import com.actelion.research.chem.phesa.Gaussian3D;
import com.actelion.research.chem.phesa.MolecularVolume;
import com.actelion.research.chem.phesa.QuickMathCalculator;
import com.actelion.research.chem.phesa.ShapeVolume;
import com.actelion.research.chem.phesa.VolumeGaussian;
import com.actelion.research.chem.phesa.pharmacophore.pp.PPGaussian;
import com.actelion.research.chem.phesa.PheSAAlignment;
/**
* @author JW, Oktober 2019
* functionality for optimizing PheSA overlap (Pharmacophore+Shape) allowing for molecular flexibility
*/
public class EvaluableFlexibleOverlap implements Evaluable {
//private static final double SCALE = -250;
//private static final double DELTA = -0.01;
private static final double LAMBDA = 0.0625;
private double e0 = 0.0;
private Conformer fitConf;
private StereoMolecule refMol;
private PheSAAlignment shapeAlign;
private boolean[] isHydrogen;
private double[] v; // internal coordinates of the atoms
private double[][] precalcPow;
private double[] precalcExp;
private double oAA;
private double oAApp;
private ForceFieldMMFF94 ff;
private Map ffOptions;
private PheSASetting settings;
private Coordinates[] origCoords;
private Coordinates[] cachedCoords; //used for gradient calculation: ligand coordinates with adjusted dihedral angles, but before rotation and translation
private double[][] dRdvi1;
private double[][] dRdvi2;
private double[][] dRdvi3;
private Coordinates origCOM;
private BondRotationHelper torsionHelper;
public EvaluableFlexibleOverlap(PheSAAlignment shapeAlign, StereoMolecule refMol, Conformer fitConf, BondRotationHelper torsionHelper,
PheSASetting settings, boolean[] isHydrogen, Map ffOptions) {
ForceFieldMMFF94.initialize(ForceFieldMMFF94.MMFF94SPLUS);
this.ffOptions = ffOptions;
this.shapeAlign = shapeAlign;
this.fitConf = fitConf;
this.isHydrogen = isHydrogen;
this.settings = settings;
this.refMol = refMol;
this.torsionHelper = torsionHelper;
init();
}
private void init() {
ff = new ForceFieldMMFF94(fitConf.getMolecule(), ForceFieldMMFF94.MMFF94SPLUS, this.ffOptions);
setInitialState();
origCoords = new Coordinates[fitConf.getMolecule().getAllAtoms()];
cachedCoords = new Coordinates[fitConf.getMolecule().getAllAtoms()];
origCOM = new Coordinates();
for(int a=0;a5) { //torsions
if(v[i]>Math.PI) {
v[i] -= 2*Math.PI;
}
}
this.v[i] = v[i];
}
updateLigandCoordinates();
shapeAlign.getMolGauss().update(fitConf);
}
public double[] getState(double[] v){
for(int i=0;i0.0) {
totalOverlap += atomOverlap;
gradientPrefactor = atomOverlap*-2*refAt.getWidth()*fitAt.getWidth()/(refAt.getWidth()+fitAt.getWidth());
}
}
grad[3*a] += (xj-xi)*gradientPrefactor;
grad[3*a+1] += (yj-yi)*gradientPrefactor;
grad[3*a+2] += (zj-zi)*gradientPrefactor;
}
}
if(refMolGauss instanceof MolecularVolume) {
for(VolumeGaussian refVG:((MolecularVolume)refMolGauss).getVolumeGaussians()){
double xi = refVG.getCenter().x;
double yi = refVG.getCenter().y;
double zi = refVG.getCenter().z;
for(AtomicGaussian fitAt:molGauss.getAtomicGaussians()){
int a = fitAt.getAtomId();
double atomOverlap = 0.0;
double xj = v[3*a];
double yj = v[3*a+1];
double zj = v[3*a+2];
double dx = xi-xj;
double dy = yi-yj;
double dz = zi-zj;
double Rij2 = dx*dx + dy*dy + dz*dz;
double alphaSum = refVG.getWidth() + fitAt.getWidth();
double gradientPrefactor=0.0;
if(Rij20.0) {
totalOverlap += atomOverlap;
gradientPrefactor = atomOverlap*-2*refVG.getWidth()*fitAt.getWidth()/(refVG.getWidth()+fitAt.getWidth());
}
}
grad[3*a] += (xj-xi)*gradientPrefactor;
grad[3*a+1] += (yj-yi)*gradientPrefactor;
grad[3*a+2] += (zj-zi)*gradientPrefactor;
}
}
}
return totalOverlap;
}
public double getFGValuePP() {
ShapeVolume molGauss = shapeAlign.getMolGauss();
ShapeVolume refMolGauss = shapeAlign.getRefMolGauss();
/**
* derivative of ShapeOverlap with respect to the four elements of the quaternion and three elements of translation
*
*/
double totalOverlap = 0.0;
for(PPGaussian refPP:refMolGauss.getPPGaussians()){
double xi = refPP.getCenter().x;
double yi = refPP.getCenter().y;
double zi = refPP.getCenter().z;
for(PPGaussian fitPP:molGauss.getPPGaussians()){
double atomOverlap = 0.0;
Coordinates fitCenterModCoord = fitPP.getCenter();
double xj = fitCenterModCoord.x;
double yj = fitCenterModCoord.y;
double zj = fitCenterModCoord.z;
double dx = xi-xj;
double dy = yi-yj;
double dz = zi-zj;
double Rij2 = dx*dx + dy*dy + dz*dz;
double alphaSum = refPP.getWidth() + fitPP.getWidth();
if(Rij20.0) {
double sim = refPP.getInteractionSimilarity(fitPP);
atomOverlap *= sim;
totalOverlap += atomOverlap;
}
}
}
}
return totalOverlap;
}
public double getFGValueShapeSelf(double[] grad, ShapeVolume molGauss,boolean rigid) {
for(int i=0;i0.0) {
gradientPrefactor = atomOverlap*-2*refAt.getWidth()*fitAt.getWidth()/(refAt.getWidth()+fitAt.getWidth());
}
}
grad[3*b] += (2*xj-2*xi)*gradientPrefactor;
grad[3*b+1] += (2*yj-2*yi)*gradientPrefactor;
grad[3*b+2] += (2*zj-2*zi)*gradientPrefactor;
return atomOverlap;
}
public double getFGValueSelfPP(ShapeVolume molVol,boolean rigid) {
double xi,yi,zi,xj,yj,zj;
/**
* derivative of ShapeOverlap with respect to Cartesian coordinates
*/
double totalOverlap = 0.0;
for(PPGaussian refPP:molVol.getPPGaussians()){
if(rigid) {
xi = refPP.getCenter().x;
yi = refPP.getCenter().y;
zi = refPP.getCenter().z;
}
else {
xi = refPP.getCenter().x;
yi = refPP.getCenter().y;
zi = refPP.getCenter().z;
}
for(PPGaussian fitPP:molVol.getPPGaussians()){
double atomOverlap = 0.0;
if(rigid) {
xj = fitPP.getCenter().x;
yj = fitPP.getCenter().y;
zj = fitPP.getCenter().z;
}
else {
xj = fitPP.getCenter().x;
yj = fitPP.getCenter().y;
zj = fitPP.getCenter().z;
}
double dx = xi-xj;
double dy = yi-yj;
double dz = zi-zj;
double Rij2 = dx*dx + dy*dy + dz*dz;
double alphaSum = fitPP.getWidth() + fitPP.getWidth();
if(Rij20.0) {
double sim = refPP.getInteractionSimilarity(fitPP);
atomOverlap *= sim;
totalOverlap += atomOverlap;
}
}
}
}
return totalOverlap;
}
public double[] getCartState(){
double[] cartState = new double[3*fitConf.getMolecule().getAllAtoms()];
for(int a=0;a