org.openmolecules.chem.conf.gen.BaseConformer 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 org.openmolecules.chem.conf.gen;
import com.actelion.research.chem.Coordinates;
import com.actelion.research.chem.StereoMolecule;
import com.actelion.research.chem.conf.Conformer;
import com.actelion.research.chem.conf.TorsionDB;
import java.util.ArrayList;
import java.util.Random;
/**
* A BaseConformer is an extended Conformer that uses a specific set of rigid fragment conformations.
* A BaseConformer knows all rotatable bonds, which are those bonds connecting the rigid fragments.
* For every rotatable bond, the BaseConformer calculates a set of possible torsion angles and their
* likelyhoods considering the specific conformations of the pair of connected rigid fragments.
* For every rotatable bond the number and order of torsion angles matches the one of the TorsionDB for
* matching bond environments, however, torsion values and their likelyhoods may be adjusted to reflect
* the particular sterical situation of the fragments being connected.
* Typically, multiple Conformer instances are derived from the same BaseConformer by creating random or
* likely torsion sets using the deriveConformer() method.
* A BaseConformer also keeps track of EliminationRules, which is a dictionary of bond torsion combinations
* that cause severe atom clashes. The dictionary grows when trying to assemble new Conformers based on
* torsion sets, which cause clashes.
*/
public class BaseConformer extends Conformer {
private static final double ANGLE_TOLERANCE = 0.001f; // limit for considering bonds as parallel
// For every optimum torsion we check for collisions and in case we try, whether left or right
// range ends are substantially better.
private static final double ACCEPTABLE_CENTER_STRAIN = 0.05f; // if we have less strain than this we don't check whether range edges improve strain
private static final double NECESSARY_EDGE_STRAIN_IMPROVEMENT_FACTOR = 0.8f; // necessary strain improvement to use edge torsion rather than center torsion
private static final double MAXIMUM_CENTER_STRAIN = 0.2f; // if center torsion is above this limit we refuse that torsion
// If no acceptable torsion remains after checking and attempting to use range edge values,
// we try to modify the least bad torsion stepwise until it is acceptable.
private static final int ESCAPE_ANGLE = 8; // degrees to rotate the rotatable bond to escape collisions
private static final int ESCAPE_STEPS = 4; // how often we apply this rotation trying to solve the collision
private static final double MIN_ESCAPE_GAIN_PER_STEP = 0.05;
private short[][] mTorsion;
private short[][] mFrequency;
private double[][] mLikelyhood; // considering directly connected rigid fragments (frequency and collision strain)
private int[] mBestTorsionIndex;
private RotatableBond[] mRotatableBond;
private ArrayList mEliminationRuleList;
private Random mRandom;
private int mConformerCount;
public BaseConformer(StereoMolecule mol, RigidFragment[] rigidFragment, RotatableBond[] rotatableBond, int[] fragmentPermutation, Random random) {
super(mol);
mRandom = random;
mRotatableBond = rotatableBond;
mTorsion = new short[rotatableBond.length][];
mFrequency = new short[rotatableBond.length][];
mLikelyhood = new double[rotatableBond.length][];
short[][][] defaultTorsionRange = new short[rotatableBond.length][][];
for (int bondIndex=0; bondIndex();
}
public RotatableBond[] getRotatableBonds() {
return mRotatableBond;
}
public int getTorsion(int bondIndex, int torsionIndex) {
return mTorsion[bondIndex][torsionIndex];
}
public double getTorsionLikelyhood(int bondIndex, int torsionIndex) {
return mLikelyhood[bondIndex][torsionIndex];
}
public int[] getMostLikelyTorsionIndexes() {
return mBestTorsionIndex;
}
public Conformer deriveConformer(int[] torsionIndex, String name) {
mConformerCount++;
Conformer conformer = new Conformer(this);
conformer.setName(name);
for (int rb=mRotatableBond.length-1; rb>=0; rb--)
rotateToIndex(conformer, mRotatableBond[rb], rb, torsionIndex[rb]);
return conformer;
}
/**
* @return number of derived (not necessarily successful) conformers
*/
public int getDerivedConformerCount() {
return mConformerCount;
}
public ArrayList getEliminationRules() {
return mEliminationRuleList;
}
/**
* Calculates a random torsion index giving torsions with higher likelyhoods
* (i.e. frequencies and collision strains) a higher chance to be selected.
* With a progress value of 0.0 selection likelyhoods are proportional to
* the torsion frequencies. With increasing progress value the higher frequent
* torsions are less and less preferred until 1.0 without any preference.
* @param random
* @param progress 0...1
*/
public int getLikelyRandomTorsionIndex(int bondIndex, double random, double progress) {
double sum = 0;
for (int t=0; t fragment2.getCoreSize()) ? fragment1 : fragment2;
int largerFragmentNo = (fragment1.getCoreSize() > fragment2.getCoreSize()) ? fragmentNo1 : fragmentNo2;
isAttached[largerFragmentNo] = true;
int fragmentConformer = (fragmentPermutation == null) ? 0 : fragmentPermutation[largerFragmentNo];
for (int i=0; i Math.PI / 2) ?
froot.subC(fragment.getExtendedCoordinates(fragmentConformer, i))
: fragment.getExtendedCoordinates(fragmentConformer, i).subC(froot);
}
}
else {
Coordinates rotationAxis;
if (alpha < Math.PI - ANGLE_TOLERANCE) { // normal case, just rotate around the plane orthogonal
rotationAxis = uv.cross(fuv);
}
else { // special case both axes anti-parallel: for cross-product we need any vector being different from uv
if (Math.abs(uv.x) >= Math.abs(uv.y) && Math.abs(uv.x) >= Math.abs(uv.z))
rotationAxis = new Coordinates(-(uv.y + uv.z) / uv.x, 1.0, 1.0);
else if (Math.abs(uv.y) >= Math.abs(uv.x) && Math.abs(uv.y) >= Math.abs(uv.z))
rotationAxis = new Coordinates(1.0, -(uv.x + uv.z) / uv.y, 1.0);
else
rotationAxis = new Coordinates(1.0, 1.0, -(uv.x + uv.y) / uv.z);
}
double[][] m = getRotationMatrix(rotationAxis.unit(), alpha);
for (int i=0; i 0f)
totalLikelyhood += mLikelyhood[bondIndex][t];
// make sure, we have at least one torsion with positive likelyhood, because only those are considered later
if (mLikelyhood[bondIndex][mBestTorsionIndex[bondIndex]] <= 0f) {
mLikelyhood[bondIndex][mBestTorsionIndex[bondIndex]] = 1.0f;
int angle = bestTorsionEdgeUsed == 1 ? -ESCAPE_ANGLE
: bestTorsionEdgeUsed == 2 ? ESCAPE_ANGLE
: (mRandom.nextDouble() < 0.5) ? -ESCAPE_ANGLE : ESCAPE_ANGLE;
for (int step=1; step<=ESCAPE_STEPS; step++) {
currentTorsion = (short)(mTorsion[bondIndex][mBestTorsionIndex[bondIndex]]+angle*step);
double[][] m = getRotationMatrix(uv, Math.PI * currentTorsion / 180 - startTorsion);
for (int i=0; i= 360)
torsion -= 360;
int bond = rotatableBond.getBond();
if (torsion != conformer.getBondTorsion(bond)) {
int deltaTorsion = torsion - conformer.getBondTorsion(bond);
rotateSmallerSide(conformer, rotatableBond, Math.PI * deltaTorsion / 180.0);
conformer.setBondTorsion(bond, torsion);
}
}
/**
* Rotates all atoms in atomList around the axis leading from atom1 through atom2
* by angle. The coordinates of atom1 and atom2 are not touched.
* @param alpha
*/
private void rotateSmallerSide(Conformer conformer, RotatableBond rotatableBond, double alpha) {
int[] torsionAtom = rotatableBond.getTorsionAtoms();
int rotationCenter = rotatableBond.getRotationCenter();
Coordinates t2 = conformer.getCoordinates(torsionAtom[2]);
Coordinates unit = t2.subC(conformer.getCoordinates(torsionAtom[1])).unit();
double[][] m = getRotationMatrix(unit, (rotationCenter == torsionAtom[1]) ? alpha : -alpha);
for (int atom:rotatableBond.getSmallerSideAtoms()) {
if (atom != rotationCenter) {
double x = conformer.getX(atom) - t2.x;
double y = conformer.getY(atom) - t2.y;
double z = conformer.getZ(atom) - t2.z;
conformer.setX(atom, x*m[0][0]+y*m[0][1]+z*m[0][2] + t2.x);
conformer.setY(atom, x*m[1][0]+y*m[1][1]+z*m[1][2] + t2.y);
conformer.setZ(atom, x*m[2][0]+y*m[2][1]+z*m[2][2] + t2.z);
}
}
}
}