com.actelion.research.chem.phesaflex.FlexibleShapeAlignment 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.Arrays;
import java.util.HashMap;
import java.util.Map;
import java.util.Random;
import com.actelion.research.chem.StereoMolecule;
import com.actelion.research.chem.alignment3d.PheSAAlignmentOptimizer.PheSASetting;
import com.actelion.research.chem.alignment3d.PheSAAlignmentOptimizer.SimilarityMode;
import com.actelion.research.chem.alignment3d.transformation.TransformationSequence;
import com.actelion.research.chem.conf.BondRotationHelper;
import com.actelion.research.chem.conf.Conformer;
import com.actelion.research.chem.forcefield.mmff.ForceFieldMMFF94;
import com.actelion.research.chem.forcefield.mmff.MMFFPositionConstraint;
import com.actelion.research.chem.optimization.MCHelper;
import com.actelion.research.chem.optimization.OptimizerLBFGS;
import com.actelion.research.chem.phesa.MolecularVolume;
import com.actelion.research.chem.phesa.PheSAAlignment;
/**
* Performs flexible Alignment of two Molecules that are prealigned. A fit molecule is thereby
* aligned to a rigid template. The internal degrees of freedom of the fit molecule are flexible
* in order to not create molecular conformations with too high strain energy, the force field potential
* energy is part of the objective function of the alignment optimization. The procedure is inspired by
* doi: 10.1021/acs.jcim.7b00618, but uses analytical gradients and optimizes the directionality overlap of
* pharmacophore features. Furthermore, in addition to local optimization, Monte Carlo steps are performed to randomly
* perturb dihedral angles of the base molecule.
* The first step is to perform constrained minimizations of the rigidly prealigned fit molecule. The well-depth of the positional
* constraint is gradually increased until to strain energy is lower than 10 kcal/mol.
* This yields a good initial guess of a energetically low-lying conformation of the fit molecule that still
* has a reasonable PheSA overlap. From this first conformation, a Monte Carlo optimization with local minimization is started.
* The objective function is:
* O = -Tphesa + lambda*Estrain*Estrain
* Estrain = E-E0, E0=10 kcal/mol
* the strain penalty is only applied if the strain energy is higher than 10 kcal/mol, otherwise the
* objective function is solely depending on the PheSA similarity Tphesa
* @author JW
*
*/
public class FlexibleShapeAlignment {
private static final int MC_STEPS = 50;
public static final double ENERGY_CUTOFF = 10.0;
private StereoMolecule refMol;
private StereoMolecule fitMol;
private MolecularVolume refVol;
private MolecularVolume fitVol;
private Map ffOptions;
private PheSASetting settings;
public PheSASetting getSettings() {
return settings;
}
public void setSettings(PheSASetting settings) {
this.settings = settings;
}
public FlexibleShapeAlignment(StereoMolecule refMol, StereoMolecule fitMol) {
this(refMol, fitMol, new MolecularVolume(refMol), new MolecularVolume(fitMol));
}
public FlexibleShapeAlignment(StereoMolecule refMol,StereoMolecule fitMol, MolecularVolume refVol, MolecularVolume fitVol) {
this.refMol = refMol;
this.fitMol = fitMol;
this.refVol = refVol;
this.fitVol = fitVol;
ffOptions = new HashMap();
ffOptions.put("dielectric constant", 4.0);
settings = new PheSASetting();
}
public double[] align() {
double[] result = new double[3];
PheSAAlignment shapeAlign = new PheSAAlignment(refVol,fitVol);
double e0 = calcMin(fitMol);
if(Double.isNaN(e0)) {
System.err.print("no force field parameters for this structure");
return result;
}
restrainedRelaxation(fitMol,e0);
boolean[] isHydrogen = new boolean[fitMol.getAllAtoms()];
for(int at=0;at1.0) //can happen because of weights
Tshape = 1.0f;
double correctionFactor = shapeAlign.getRefMolGauss().getPPGaussians().size()/shapeAlign.getRefMolGauss().getPPGaussians().stream().mapToDouble(g -> g.getWeight()).sum();
double ObbPP = eval.getFGValueSelfPP(shapeAlign.getMolGauss(),false);
double OaaPP = eval.getFGValueSelfPP(shapeAlign.getRefMolGauss(),true);
double OabPP = eval.getFGValuePP();
double Tpp = 0.0;
if(shapeAlign.getRefMolGauss().getPPGaussians().size()==0 && shapeAlign.getMolGauss().getPPGaussians().size()==0 )
Tpp = 1.0;
else {
if(tversky)
Tpp = OabPP/(tverskyCoeff*ObbPP+(1.0-tverskyCoeff)*OaaPP);
else
Tpp = (OabPP/(OaaPP+ObbPP-OabPP));
}
Tpp*=correctionFactor;
if(!tversky && Tshape>1.0) //can happen because of weights
Tshape = 1.0f;
if(!tversky && Tpp>1.0) //can happen because of weights
Tpp = 1.0f;
double T = (1.0f-(float)ppWeight)*Tshape + (float)ppWeight*Tpp;
return T;
}
private double[] getResult() {
TransformationSequence sequence = new TransformationSequence();
PheSAAlignment pa = new PheSAAlignment(fitMol,refMol, settings.getPpWeight());
double[] r = pa.findAlignment(new double[][] {{0.00, 0.00, 0.00,0.0,0.0,0.0}},sequence,false,settings.getSimMode());
return new double[] {r[0],r[1],r[2], r[3]};
}
public double calcMin(StereoMolecule fitMol) {
ForceFieldMMFF94.initialize(ForceFieldMMFF94.MMFF94SPLUS);
ForceFieldMMFF94 forceField = new ForceFieldMMFF94(new StereoMolecule(fitMol), ForceFieldMMFF94.MMFF94SPLUS, ffOptions);
forceField.minimise();
double e0 = forceField.getTotalEnergy();
return e0;
}
public void restrainedRelaxation(StereoMolecule fitMol, double e0) {
ForceFieldMMFF94.initialize(ForceFieldMMFF94.MMFF94SPLUS);
double init = 0.2;
boolean notRelaxed = true;
int maxCycles = 10;
int cycles = 0;
while(notRelaxed && cyclese0) && (e-e0>ENERGY_CUTOFF);
init += 0.2;
cycles++;
}
}
}