com.actelion.research.chem.phesaflex.MetropolisMonteCarloHelper Maven / Gradle / Ivy

Go to download
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