org.openmolecules.chem.conf.gen.ConformerGenerator 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
/*
* Copyright 2013-2020 Thomas Sander, openmolecules.org
*
* Redistribution and use in source and binary forms, with or without modification,
* are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
* 3. Neither the name of the copyright holder nor the names of its contributors
* may be used to endorse or promote products derived from this software without
* specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
* EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
* SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
* @author Thomas Sander
*/
package org.openmolecules.chem.conf.gen;
import com.actelion.research.calc.ThreadMaster;
import com.actelion.research.chem.*;
import com.actelion.research.chem.conf.Conformer;
import com.actelion.research.chem.conf.TorsionDB;
import com.actelion.research.chem.conf.VDWRadii;
import com.actelion.research.util.DoubleFormat;
import com.actelion.research.util.IntArrayComparator;
import org.openmolecules.chem.conf.so.ConformationSelfOrganizer;
import org.openmolecules.chem.conf.so.SelfOrganizedConformer;
import java.util.Arrays;
import java.util.Random;
import java.util.TreeMap;
/**
* This class generates 3D-conformers of a given molecule using the following strategy:
* All rotatable, non-ring bonds are determined.
* Fragments separated by rotatable bonds are considered rigid, but there may be more than one possible
* fragment conformer, e.g. chair- and boat conformers of a saturated 6-membered ring.
* These atom coordinate sets of rigid fragments are handed out by a dedicated RigidFragmentProvider instance,
* which either generates them using a self organization algorithm, or which takes it from a cache.
* For every rotatable bond a list of preferred torsion angles is determined based on from a COD statistics
* of similar bond environments.
* For individual torsion values likelihoods are estimated based on frequency and atom collisions of vicinal fragments.
* A dedicated (systematic, biased or random) torsion set strategy delivers collision-free torsion sets, i.e. conformers.
*
* For generating conformers in multiple threads, every thread needs its own ConformerGenerator instance.
* If they use a RigidFragmentCache, then the cache is shared among all ConformerGenerators.
* Important: Input molecules should contain absolute stereo centers. If they contain undefined or ESR type '&' or 'or'
* stereo centers, then a ConformerGenerator randomly takes one of the possible stereo isomers and generates conformers
* for that. If you want conformers for all possible stereo isomers of a molecules with non-absolute stereo centers,
* you should use a StereoIsomerEnumerator to produce all possible stereo isomers and then produce conformers for every
* one of them. If half of a set of stereo isomers consists of the enantiomers of the other half, then it is advisable
* to generate conformes for one half only and to generate the second half by just mirroring the first halfs coordinates.
* To do that use option skipEnantiomers==true create a mirrored set of conformers, if isSkippingEnantiomers() of the
* StereoIsomerEnumerator returns true.
*/
public class ConformerGenerator {
public static final int STRATEGY_LIKELY_SYSTEMATIC = 1;
public static final int STRATEGY_PURE_RANDOM = 2;
public static final int STRATEGY_LIKELY_RANDOM = 3;
public static final int STRATEGY_ADAPTIVE_RANDOM = 4;
protected static final double VDW_TOLERANCE_HYDROGEN = 0.90; // factor on VDW radii for minimum tolerated non bound atom distances
protected static final double VDW_TOLERANCE_OTHER = 0.90; // factor on VDW radii for minimum tolerated non bound atom distances
private static final int ESCAPE_ANGLE = 8; // degrees to rotate two adjacent rotatable bonds to escape collisions
private static final int ESCAPE_STEPS = 3; // how often we apply this rotation trying to solve the collision
private static final double MIN_ESCAPE_GAIN_PER_STEP = 1.0;
// We try to translate arbitrary collision values into a kcal/mol energy scale
public static final double COLLISION_STRAIN_TO_ENERGY_FACTOR = 20;
private StereoMolecule mMolecule;
private TreeMap mBaseConformerMap;
private RotatableBond[] mRotatableBond;
private RigidFragment[] mRigidFragment;
private ConformationSelfOrganizer mSelfOrganizer;
private RigidFragmentProvider mRigidFragmentProvider;
private TorsionSetStrategy mTorsionSetStrategy;
private TorsionSet mTorsionSet;
private long mRandomSeed;
private int mDisconnectedFragmentCount,mAllConformerCount,mReturnedConformerCount;
private boolean mUseSelfOrganizerIfAllFails,mIsDiagnosticsMode,mIsFinished;
private int[] mFragmentNo,mDisconnectedFragmentNo,mDisconnectedFragmentSize;
private boolean[][] mSkipCollisionCheck;
private Random mRandom;
private ThreadMaster mThreadMaster;
private ConformerSetDiagnostics mDiagnostics;
/**
* Assuming that the given molecule has 2D-coordinates, this method
* converts all implicit hydrogen atoms into explicit ones by filling valences
* and adapting for atom charges. New hydrogen atoms receive new 2D-coordinates
* by equally locating them between those two neighbors with the widest angle between
* their bonds. Any stereo configurations deducible from 2D-coordinates are retained.
* @param mol
*/
public static void addHydrogenAtoms(StereoMolecule mol) {
// We may have parities but empty coordinates. In this case we need to protect parities.
int oldStereoHelperBits = mol.getHelperArrayStatus() & Molecule.cHelperBitsStereo;
mol.ensureHelperArrays(Molecule.cHelperNeighbours);
int[] implicitHydrogen = new int[mol.getAtoms()];
for (int atom=0; atom(new IntArrayComparator());
return getNextConformer();
}
else {
ConformationSelfOrganizer sampler = new ConformationSelfOrganizer(mol, true);
Conformer conformer = sampler.generateOneConformer(mRandomSeed);
separateDisconnectedFragments(conformer);
conformer.setName("SO#1");
return conformer;
}
}
/**
* Don't call this method directly. Rather call initializeConformers() !!!
* Adds implicit hydrogens to the molecule and determines all rotatable bonds,
* which are not part of a ring. Generates rigid fragments between rotatable bonds.
* The base conformer is constructed by connecting all rigid fragments using the
* most frequent torsions. Atoms of different fragments in the base conformer may
* collide. In order to obtain collision free conformers choose a TorsionStrategy
* and call getNextConformer() at least once.
* @param mol
* @param use60degreeSteps use 60 degree steps for every rotatable bond instead of torsion DB
*/
protected boolean initialize(StereoMolecule mol, boolean use60degreeSteps) {
mSelfOrganizer = null;
mol.ensureHelperArrays(Molecule.cHelperNeighbours);
for (int atom=0; atom mol.getMaxValence(atom))
return false;
addHydrogenAtoms(mol);
// we need to protect explicit hydrogens in the fragments that we create from mol
mol.setHydrogenProtection(true);
int[] atomParity = null;
int[] bondParity = null;
boolean[] atomParityIsPseudo = null;
boolean[] bondParityIsPseudo = null;
// if we have parities and no atom coordinates, parities cannot be correctly recreated
// by the Canonizer and, thus, need to be cached and copied back after the symmetry detection.
if ((mol.getHelperArrayStatus() & Molecule.cHelperBitParities) != 0) {
atomParity = new int[mol.getAtoms()];
atomParityIsPseudo = new boolean[mol.getAtoms()];
for (int atom=0; atom
Integer.compare(b2.getSmallerSideAtoms().length, b1.getSmallerSideAtoms().length) );
if (mIsDiagnosticsMode)
mDiagnostics = new ConformerSetDiagnostics(this);
// TODO this is actually only used with random conformers. Using a conformer mode
// may save some time, if systematic conformers are created.
initializeCollisionCheck();
return true;
}
/**
* After calling initializeConformers() this method returns the number of rotatable bonds,
* which are used to separate the molecule into rigid fragments.
* @return
*/
public int getRotatableBondCount() {
return mRotatableBond == null ? 0 : mRotatableBond.length;
}
public TreeMap getBaseConformerMap() {
return mBaseConformerMap;
}
public TorsionSetStrategy getTorsionSetStrategy() {
return mTorsionSetStrategy;
}
public BaseConformer getBaseConformer(int[] fragmentPermutation) {
BaseConformer baseConformer = mBaseConformerMap.get(fragmentPermutation);
if (baseConformer != null)
return baseConformer;
baseConformer = new BaseConformer(mMolecule, mRigidFragment, mRotatableBond, fragmentPermutation, mRandom);
//printDebugConformers(baseConformer);
mBaseConformerMap.put(fragmentPermutation.clone(), baseConformer);
return baseConformer;
}
/**
* Creates the next random, likely or systematic new(!) conformer of the molecule
* that was passed when calling initializeConformers(). A new conformer is one,
* whose combination of torsion angles was not used in a previous conformer
* created by this function since the last call of initializeConformers().
* If parameter mol is null, then a compact copy of the original molecule with the new
* coordinates is returned. You may pass the original molecule or a copy of it to
* recycle the original molecule to receive new 3D coordinates.
* Every call of this method creates a new collision-free conformer until the employed torsion set
* strategy decides that it cannot generate any more suitable torsion sets.
* @param mol null or molecule used during initialization or a copy of it
* @return conformer or null, if all/maximum torsion permutations have been tried
*/
public StereoMolecule getNextConformerAsMolecule(StereoMolecule mol) {
Conformer conformer = getNextConformer();
return (conformer == null) ? null : conformer.toMolecule(mol);
}
public Conformer getNextConformer() {
return getNextConformer(null);
}
/**
* Creates the next random, likely or systematic new(!) conformer of the molecule
* that was passed when calling initializeConformers(). A new conformer is one,
* whose combination of torsion angles was not used in a previous conformer
* created by this function since the last call of initializeConformers().
* Every call of this method creates a new collision-free conformer until the employed torsion set
* strategy decides that it cannot generate any more suitable torsion sets.
* @param torsionSetHolder : will contain the TorsionSet that gave rise to the conformer.
* @return conformer or null, if all/maximum torsion permutations have been tried
*/
public Conformer getNextConformer(TorsionSet[] torsionSetHolder) {
if (mRotatableBond == null && mSelfOrganizer == null)
return null;
if (mSelfOrganizer != null) {
SelfOrganizedConformer conformer = mSelfOrganizer.getNextConformer();
if (conformer != null) {
separateDisconnectedFragments(conformer);
mReturnedConformerCount++;
conformer.setName("SO#"+(++mAllConformerCount));
return conformer;
}
return null;
}
if (mIsFinished)
return null;
// create a base conformer from first set of fragments and calculate torsion likelihoods
if (mBaseConformerMap.size() == 0)
getBaseConformer(new int[mRigidFragment.length]);
mTorsionSet = mTorsionSetStrategy.getNextTorsionSet(mTorsionSet, mDiagnostics);
while (mTorsionSet != null && (mThreadMaster == null || !mThreadMaster.threadMustDie())) {
BaseConformer baseConformer = getBaseConformer(mTorsionSet.getConformerIndexes());
if (mTorsionSet.getConformer() == null) {
mTorsionSet.setConformer(baseConformer.deriveConformer(mTorsionSet.getTorsionIndexes(), "#" + (++mAllConformerCount)));
if (mIsDiagnosticsMode)
mDiagnostics.addNew(mTorsionSet);
}
// If the torsionSet has already a collision value, then it is a second choice torsion set
// with a collision value below the tolerance. No need to check or fix again.
if (mTorsionSet.getCollisionStrainSum() == 0.0) {
calculateCollision(mTorsionSet, mTorsionSet.getConformer());
if (mTorsionSet.getCollisionStrainSum() != 0.0)
tryFixCollisions(baseConformer, mTorsionSet);
if (mIsDiagnosticsMode) {
mDiagnostics.get(mTorsionSet).setConformer(baseConformer, mTorsionSet.getConformer());
mDiagnostics.get(mTorsionSet).setCollisionStrain(mTorsionSet.getCollisionStrainSum());
}
}
if (mTorsionSet.getCollisionStrainSum() > mTorsionSetStrategy.calculateCollisionTolerance()) {
mTorsionSet = mTorsionSetStrategy.getNextTorsionSet(mTorsionSet, mDiagnostics);
if (mTorsionSet != null || mReturnedConformerCount != 0)
continue;
if (mTorsionSet == null && mUseSelfOrganizerIfAllFails) {
// We couldn't create torsion strategy based conformers: switch to self organizer!
mSelfOrganizer = new ConformationSelfOrganizer(mMolecule, true);
mSelfOrganizer.setThreadMaster(mThreadMaster);
mSelfOrganizer.initializeConformers(mRandomSeed, -1);
SelfOrganizedConformer conformer = mSelfOrganizer.getNextConformer();
if (conformer != null) {
separateDisconnectedFragments(conformer);
mReturnedConformerCount++;
conformer.setName("SO#"+(++mAllConformerCount));
return conformer;
}
}
// we didn't get any torsion set that didn't collide; take the best we had
mTorsionSet = mTorsionSetStrategy.getBestCollidingTorsionIndexes();
mIsFinished = true; // we are finished with conformers
}
if (mTorsionSet != null) {
separateDisconnectedFragments(mTorsionSet.getConformer());
mTorsionSet.setUsed();
mReturnedConformerCount++;
if(torsionSetHolder != null)
torsionSetHolder[0] = new TorsionSet(mTorsionSet);
if (mIsDiagnosticsMode)
mDiagnostics.get(mTorsionSet).setSuccess(true);
return mTorsionSet.getConformer();
}
// This poses the danger to produce highly strained unrealistic conformers, e.g. for impossible compounds...
// else if (mReturnedConformerCount == 0) {
// // We couldn't create torsion strategy based conformers: try creating one conformer with self organizer!
// ConformationSelfOrganizer sampler = new ConformationSelfOrganizer(mMolecule, true);
// Conformer conformer = sampler.generateOneConformer(mRandomSeed);
// separateDisconnectedFragments(conformer);
// mReturnedConformerCount++;
// conformer.setName("SO#1");
// return conformer;
// }
}
return null;
}
/*private void printDebugConformers(BaseConformer bc) {
int degreesA = bc.getBondTorsion(mRotatableBond[0].getBond());
int degreesB = bc.getBondTorsion(mRotatableBond[1].getBond());
System.out.println("Startangles A:"+degreesA+" B:"+degreesB);
System.out.println("idcode\tidcoords\tdegrees A\tdegrees B");
for (int degrees=0; degrees<100; degrees+=20) {
Conformer co = new Conformer(bc);
bc.rotateTo(co, mRotatableBond[0], (short)(degreesA + degrees));
Canonizer ca = new Canonizer(co.toMolecule(null));
System.out.println(ca.getIDCode() + "\t" + ca.getEncodedCoordinates() + "\t" + degrees +"\t0");
}
for (int degrees=0; degrees<100; degrees+=20) {
Conformer co = new Conformer(bc);
bc.rotateTo(co, mRotatableBond[1], (short)(degreesB + degrees));
Canonizer ca = new Canonizer(co.toMolecule(null));
System.out.println(ca.getIDCode() + "\t" + ca.getEncodedCoordinates() + "\t0\t" + degrees);
}
}*/
/**
* @return count of valid delivered conformers
*/
public int getConformerCount() {
return mReturnedConformerCount;
}
/**
* Calculates the potential count of conformers by multiplying degrees of freedom
* (torsions per rotatable bond & rigid fragment multiplicities).
* Cannot be called before calling initializeConformers().
* @return
*/
public int getPotentialConformerCount() {
return (mTorsionSetStrategy == null) ? 1 : mTorsionSetStrategy.getPermutationCount();
}
/**
* If a molecule has at least one rotatable bond if all permutations
* of torsions collide beyond a tolerated strain, then the standard
* behaviour of this class is to return that clashing conformer with
* the lowest strain.
* If passing true to this method, the ConformerGenerator will use
* the ConformerSelfOrganizer in these cases to generate conformers.
* getNextConformer will then deliver conformers until the self organization
* fails to create new conformers.
* @param b
*/
public void setUseSelfOrganizerIfAllFails(boolean b) {
mUseSelfOrganizerIfAllFails = b;
}
/**
* One of the initializeConformers() methods needs to be called, before getting individual
* conformers of the same molecule by getNextConformer().
* Open valences of the passed molecule are filled with hydrogen atoms.
* The passed molecule may repeatedly be used as container for a new conformer's atom
* coordinates, if it is passed as parameter to getNextConformer().
* This method uses the STRATEGY_LIKELY_RANDOM strategy with a maximum of 100.000 distinct
* torsion sets and uses torsions from crystallographic data.
* @param mol will be saturated with hydrogen atoms
* @return false if there is a structure problem
*/
public boolean initializeConformers(StereoMolecule mol) {
return initializeConformers(mol, STRATEGY_LIKELY_RANDOM, 100000, false);
}
/**
* One of the initializeConformers() methods needs to be called, before getting individual
* conformers of the same molecule by getNextConformer().
* Open valences of the passed molecule are filled with hydrogen atoms.
* The passed molecule may repeatedly be used as container for a new conformer's atom
* coordinates, if it is passed as parameter to getNextConformer().
* @param mol will be saturated with hydrogen atoms
* @param strategy one of the STRATEGY_ constants
* @param maxTorsionSets maximum number of distinct torsion sets the strategy will try (default 100000)
* @param use60degreeSteps use 60 degree steps for every rotatable bond instead of torsion DB
* @return false if there is a structure problem
*/
public boolean initializeConformers(StereoMolecule mol, int strategy, int maxTorsionSets, boolean use60degreeSteps) {
if (!initialize(mol, use60degreeSteps))
return false;
if (mRotatableBond == null) {
mSelfOrganizer = new ConformationSelfOrganizer(mol, true);
mSelfOrganizer.setThreadMaster(mThreadMaster);
mSelfOrganizer.initializeConformers(mRandomSeed, -1);
}
else {
mBaseConformerMap = new TreeMap<>(new IntArrayComparator());
switch(strategy) {
case STRATEGY_PURE_RANDOM:
mTorsionSetStrategy = new TorsionSetStrategyRandom(this, false, mRandomSeed);
break;
case STRATEGY_LIKELY_RANDOM:
mTorsionSetStrategy = new TorsionSetStrategyRandom(this, true, mRandomSeed);
break;
case STRATEGY_ADAPTIVE_RANDOM:
mTorsionSetStrategy = new TorsionSetStrategyAdaptiveRandom(this, true, true, mRandomSeed);
break;
case STRATEGY_LIKELY_SYSTEMATIC:
mTorsionSetStrategy = new TorsionSetStrategyLikelySystematic(this);
break;
}
mTorsionSetStrategy.setMaxTotalCount(maxTorsionSets);
}
return true;
}
public RotatableBond[] getRotatableBonds() {
return mRotatableBond;
}
public RigidFragment[] getRigidFragments() {
return mRigidFragment;
}
/**
* Moves disconnected fragments along Z-axis such that there is an
* empty z-space SEPARATION_DISTANCE thick between the fragments.
* @param conformer
*/
private void separateDisconnectedFragments(Conformer conformer) {
final double SEPARATION_DISTANCE = 3.0;
if (mDisconnectedFragmentCount > 1) {
double[] meanX = new double[mDisconnectedFragmentCount];
double[] meanY = new double[mDisconnectedFragmentCount];
double[] minZ = new double[mDisconnectedFragmentCount];
double[] maxZ = new double[mDisconnectedFragmentCount];
for (int i=0; i conformer.getZ(atom))
minZ[mDisconnectedFragmentNo[atom]] = conformer.getZ(atom);
if (maxZ[mDisconnectedFragmentNo[atom]] < conformer.getZ(atom))
maxZ[mDisconnectedFragmentNo[atom]] = conformer.getZ(atom);
}
for (int i=0; i TorsionSetStrategy.MAX_COLLISION_STRAIN)
return false;
Conformer conformer = null;
Conformer localBest = null;
Conformer globalBest = null;
double bestCollisionStrain = Double.MAX_VALUE;
double origCollisionStrainSum = torsionSet.getCollisionStrainSum();
for (int f1=1; f1 MIN_ESCAPE_GAIN_PER_STEP
&& mTorsionSetStrategy.getBondsBetweenFragments(f1, f2).length == 2) {
int[] rbIndex = mTorsionSetStrategy.getBondsBetweenFragments(f1, f2);
int[] startTorsion = new int[2];
for (int i=0; i<2; i++) {
RotatableBond rotatableBond = mRotatableBond[rbIndex[i]];
startTorsion[i] = torsionSet.getConformer().getBondTorsion(rotatableBond.getBond());
// short[] torsionRange = rotatableBond.getDefaultTorsionRanges()[torsionSet.getTorsionIndexes()[rbIndex[i]]];
}
int[] torsionDif = new int[2];
for (int r1=-1; r1<=1; r1+=2) {
torsionDif[0] = r1 * ESCAPE_ANGLE;
for (int r2=-1; r2<=1; r2+=2) {
torsionDif[1] = r2 * ESCAPE_ANGLE;
double collisionStrain = origCollisionStrainMatrix[f1][f2];
for (int step=1; step<=ESCAPE_STEPS; step++) {
if (conformer == null)
conformer = new Conformer(torsionSet.getConformer());
else
conformer.copyFrom(torsionSet.getConformer());
for (int i=0; i<2; i++)
baseConformer.rotateTo(conformer, mRotatableBond[rbIndex[i]], (short)(startTorsion[i] + step * torsionDif[i]));
double strain = calculateCollisionStrain(conformer, mRigidFragment[f1], mRigidFragment[f2]);
if (strain < collisionStrain - MIN_ESCAPE_GAIN_PER_STEP) {
collisionStrain = strain;
if (localBest == null)
localBest = new Conformer(conformer);
else
localBest.copyFrom(conformer);
// not enough collision strain left for correction
if (collisionStrain < MIN_ESCAPE_GAIN_PER_STEP)
break;
}
else {
break;
}
}
if (collisionStrain < origCollisionStrainMatrix[f1][f2]
&& collisionStrain < bestCollisionStrain) {
bestCollisionStrain = collisionStrain;
if (globalBest == null)
globalBest = new Conformer(localBest);
else
globalBest.copyFrom(localBest);
}
}
}
}
}
}
}
if (globalBest == null)
return false;
calculateCollision(torsionSet, globalBest);
if (torsionSet.getCollisionStrainSum() >= origCollisionStrainSum) {
// local improvement causes worsening somewhere else
torsionSet.setCollisionStrain(origCollisionStrainSum, origCollisionStrainMatrix);
return false;
}
torsionSet.getConformer().copyFrom(globalBest);
return true;
}
private double calculateCollisionStrain(Conformer conformer, RigidFragment f1, RigidFragment f2) {
double collisionStrainSum = 0;
for (int i=0; i