com.actelion.research.chem.phesaflex.MetropolisMonteCarloHelper 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.ArrayList;
import java.util.List;
import java.util.Random;
import java.util.stream.IntStream;
import com.actelion.research.chem.StereoMolecule;
import com.actelion.research.chem.conf.BondRotationHelper;
import com.actelion.research.chem.conf.TorsionDB;
import com.actelion.research.chem.conf.TorsionRelevanceHelper;
/**
* @author JW
* Provides functionality to perform random dihedral angle perturbations on the 3D conformation of the molecule.
* Central bonds are perturbed by smaller values (given by torsion relevance), whereas terminal bonds can be perturbed by
* 60 degrees
*
*/
public class MetropolisMonteCarloHelper {
private StereoMolecule mol;
private static double MAX_ANGLE = 60.0/180.0*Math.PI;
private static double MIN_ANGLE = 5.0/180.0*Math.PI;
private static double TEMPERATURE = 0.0043; //move that reduces Tanimoto by 0.01 has 10% chance of acceptance
private long seed = 1234L;
private double rmax;
private double rmin;
private float[] torsionRelevance;
private int[] rotatableBonds;
private Random random;
private int previousBond;
private double previousAngle;
private double slope;
private BondRotationHelper bondRotationHelper;
public MetropolisMonteCarloHelper(StereoMolecule mol) {
this.mol = mol;
}
/**
*
* @return: boolean, indicating if MMC Helper could successfully be initialized
*/
public boolean init() {
boolean success = true;
bondRotationHelper = new BondRotationHelper(mol);
random = new Random(seed);
boolean[] isRotatableBond = new boolean[mol.getBonds()];
TorsionDB.findRotatableBonds(mol,true, isRotatableBond);
List rotBonds = new ArrayList();
IntStream.range(0, isRotatableBond.length).forEach(e -> {
if(isRotatableBond[e])
rotBonds.add(e);
});
rotatableBonds = rotBonds.stream().mapToInt(i->i).toArray();
torsionRelevance = TorsionRelevanceHelper.getRelevance(mol, isRotatableBond);
rmin = Float.MAX_VALUE;
rmax = 0.0f;
for(float relevance : torsionRelevance) {
if (relevancermax)
rmax = relevance;
}
slope = (MIN_ANGLE-MAX_ANGLE)/(rmax-rmin);
if(rotatableBonds.length==0)
success = false;
return success;
}
public void step() {
int prefactor = random.nextInt(2)<1 ? -1 : 1;
previousBond = rotatableBonds[random.nextInt(rotatableBonds.length)];
previousAngle = MAX_ANGLE+(torsionRelevance[previousBond]-rmin)*slope;
previousAngle*=prefactor;
bondRotationHelper.rotateSmallerSide(previousBond, previousAngle);
}
public void undoStep() {
bondRotationHelper.rotateSmallerSide(previousBond, -previousAngle);
}
public boolean accept(double oldScore, double newScore) {
boolean accept = false;
if(newScore>oldScore)
accept = true;
else {
double delta = -(newScore-oldScore);
double p = Math.exp(-delta/TEMPERATURE);
double rnd = random.nextDouble();
if(rnd