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.Canonizer;
import com.actelion.research.chem.Coordinates;
import com.actelion.research.chem.Molecule;
import com.actelion.research.chem.StereoMolecule;
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.
*/
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.80; // factor on VDW radii for minimum tolerated non bound atom distances
protected static final double VDW_TOLERANCE_OTHER = 0.80; // 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 = 4; // how often we apply this rotation trying to solve the collision
private static final double MIN_ESCAPE_GAIN_PER_STEP = 0.05;
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,mConformerCount;
private boolean mUseSelfOrganizerIfAllFails;
private double mContribution;
private int[] mFragmentNo,mDisconnectedFragmentNo,mDisconnectedFragmentSize;
private boolean[][] mSkipCollisionCheck;
private Random mRandom;
private ThreadMaster mThreadMaster;
public static final boolean PRINT_TORSION_AND_FRAGMENT_LIKELYHOODS = false;
public static final boolean PRINT_DEBUG_INDEXES = false;
public static final boolean PRINT_EXIT_REASON = false;
public static final boolean PRINT_ELIMINATION_RULES_WITH_STRUCTURES = false;
public static final String DW_FRAGMENTS_FILE = "/home/thomas/data/debug/conformationGeneratorFragments.dwar";
public String mDiagnosticCollisionString,mDiagnosticTorsionString; // TODO get rid of this
public int[] mDiagnosticCollisionAtoms; // TODO get rid of this
public static boolean WRITE_DW_FRAGMENT_FILE = false;
/**
* Adds explicit hydrogen atoms where they are implicit 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.
boolean paritiesValid = (mol.getHelperArrayStatus() & Molecule.cHelperBitParities) != 0;
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);
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");
// writer.newLine();
// writer.write("");
// writer.newLine();
// writer.write("");
// writer.newLine();
// writer.write("");
// writer.newLine();
// writer.write("");
// writer.newLine();
// writer.write("");
// writer.newLine();
// writer.write("");
// writer.newLine();
// writer.write("Structure\tcoords");
// writer.newLine();
// for (Rigid3DFragment f:mRigidFragment) {
// StereoMolecule[] fragment = f.getFragment();
// for (int i=0; i
Integer.compare(b2.getSmallerSideAtomCount(), b1.getSmallerSideAtomCount()) );
// 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;
}
private Conformer getBaseConformer(int[] fragmentPermutation) {
Conformer baseConformer = mBaseConformerMap.get(fragmentPermutation);
if (baseConformer != null)
return baseConformer;
baseConformer = new Conformer(mMolecule);
boolean[] isAttached = new boolean[mRigidFragment.length];
for (RotatableBond rb:mRotatableBond)
rb.connectFragments(baseConformer, isAttached, fragmentPermutation);
// for separated fragments without connection points we need to get coordinates
for (int i=0; i=0; j--)
mRotatableBond[j].rotateToIndex(conformer, mTorsionSet.getTorsionIndexes()[j]);
if (WRITE_DW_FRAGMENT_FILE) {
mDiagnosticTorsionString = ""+mTorsionSet.getTorsionIndexes()[0];
for (int i=1; i" + mTorsionSet.getConformerIndexes()[0];
for (int i=1; i mTorsionSetStrategy.calculateCollisionTolerance()) {
//System.out.println("COLLIDES!");
// TODO remove
String idcode,coords;
int elimRules;
if (PRINT_ELIMINATION_RULES_WITH_STRUCTURES) {
Canonizer can = new Canonizer(conformer.toMolecule(null));
idcode = can.getIDCode();
coords = can.getEncodedCoordinates();
elimRules = mTorsionSetStrategy.getEliminationRuleList().size();
}
mTorsionSet = mTorsionSetStrategy.getNextTorsionSet(mTorsionSet);
// TODO remove
if (PRINT_ELIMINATION_RULES_WITH_STRUCTURES) {
if (elimRules != mTorsionSetStrategy.getEliminationRuleList().size()) {
if (elimRules == 0)
System.out.println("idcode\tidcoords\telimRules\tstrategy");
StringBuilder sb = new StringBuilder();
for (int i = elimRules; i < mTorsionSetStrategy.getEliminationRuleList().size(); i++) {
TorsionSetEliminationRule rule = mTorsionSetStrategy.getEliminationRuleList().get(i);
sb.append(mTorsionSetStrategy.eliminationRuleString(rule));
sb.append(" m:" + Long.toHexString(rule.getMask()[0]));
sb.append(" d:" + Long.toHexString(rule.getData()[0]));
sb.append("");
}
System.out.println(idcode + "\t" + coords + "\t" + sb.toString() + "\tlikelyRandom");
}}
if (mTorsionSet != null || mConformerCount != 0)
continue;
if (mUseSelfOrganizerIfAllFails) {
mSelfOrganizer = new ConformationSelfOrganizer(mMolecule, true);
mSelfOrganizer.setThreadMaster(mThreadMaster);
mSelfOrganizer.initializeConformers(mRandomSeed, -1);
conformer = mSelfOrganizer.getNextConformer();
if (conformer != null) {
separateDisconnectedFragments(conformer);
mConformerCount++;
return conformer;
}
}
// we didn't get any torsion set that didn't collide; take the best we had
mTorsionSet = mTorsionSetStrategy.getBestCollidingTorsionIndexes();
conformer = new Conformer(getBaseConformer(mTorsionSet.getConformerIndexes()));
for (int j=mRotatableBond.length-1; j>=0; j--)
mRotatableBond[j].rotateToIndex(conformer, mTorsionSet.getTorsionIndexes()[j]);
mBaseConformerMap = null; // we are finished with conformers
}
//System.out.println("passed collision check! "+mTorsionSet.toString());
separateDisconnectedFragments(conformer);
mContribution = mTorsionSetStrategy.getContribution(mTorsionSet);
mTorsionSet.setUsed();
mConformerCount++;
if(torsion_set_out != null) {
torsion_set_out[0] = new TorsionSet(mTorsionSet);
}
return conformer;
}
return null;
}
/**
* @return count of valid delivered conformers
*/
public int getConformerCount() {
return mConformerCount;
}
/**
* 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();
}
/**
* With best current knowledge about colliding torsion combinations
* and based on the individual frequencies of currently active torsions
* this method returns the conformers's overall contribution to the
* total set of non colliding conformers.
* @return this conformer's contribution to all conformers
*/
public double getPreviousConformerContribution() {
return mRotatableBond == null ? 1f : mContribution;
}
/**
* 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 {
switch(strategy) {
case STRATEGY_PURE_RANDOM:
mTorsionSetStrategy = new TorsionSetStrategyRandom(mRotatableBond, mRigidFragment, false, mRandomSeed);
break;
case STRATEGY_LIKELY_RANDOM:
mTorsionSetStrategy = new TorsionSetStrategyRandom(mRotatableBond, mRigidFragment, true, mRandomSeed);
break;
case STRATEGY_ADAPTIVE_RANDOM:
mTorsionSetStrategy = new TorsionSetStrategyAdaptiveRandom(mRotatableBond, mRigidFragment, true, true, mRandomSeed);
break;
case STRATEGY_LIKELY_SYSTEMATIC:
mTorsionSetStrategy = new TorsionSetStrategyLikelySystematic(mRotatableBond, mRigidFragment);
break;
}
mTorsionSetStrategy.setMaxTotalCount(maxTorsionSets);
mBaseConformerMap = new TreeMap(new IntArrayComparator());
}
return true;
}
/**
* 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";
mDiagnosticCollisionString = mDiagnosticCollisionString+"a1:"+atom1+" f1:"+mFragmentNo[atom1]+" a2:"+atom2+" f2:"+mFragmentNo[atom2]+" distance:"+distance+" min:"+minDistance;
if (mDiagnosticCollisionAtoms == null) {
mDiagnosticCollisionAtoms = new int[2];
mDiagnosticCollisionAtoms[0] = atom1;
mDiagnosticCollisionAtoms[1] = atom2;
}
}
if (collisionIntensityMatrix == null)
collisionIntensityMatrix = new double[mRigidFragment.length][];
int f1 = mFragmentNo[atom1];
int f2 = mFragmentNo[atom2];
if (f1 < f2) {
if (collisionIntensityMatrix[f2] == null)
collisionIntensityMatrix[f2] = new double[f2];
collisionIntensityMatrix[f2][f1] += collisionIntensity;
}
else {
if (collisionIntensityMatrix[f1] == null)
collisionIntensityMatrix[f1] = new double[f1];
collisionIntensityMatrix[f1][f2] += collisionIntensity;
}
}
}
}
}
}
}
}
torsionSet.setCollisionIntensity(collisionIntensitySum, collisionIntensityMatrix);
return (collisionIntensitySum != 0);
}
/**
* If we have collisions between two fragments that have two rotatable bonds in between,
* then we try to rotate both rotatable bonds in any of these situations stepwise
* by a small angle to reduce the collision intensity sum. If we achieve a sufficient
* improvement without hampering the rest of the molecule, the passed conformer is modified.
* @param origConformer
* @param torsionSet
* @return true if successful; false if no fix possible that reduces total collision intensity below tolerance
*/
private boolean tryFixCollisions(Conformer origConformer, TorsionSet torsionSet) {
if (mThreadMaster != null && mThreadMaster.threadMustDie())
return false;
double[][] origCollisionIntensityMatrix = torsionSet.getCollisionIntensityMatrix();
double remainingCollisionIntensity = 0;
for (int f1=1; f1 TorsionSetStrategy.MAX_ALLOWED_COLLISION_INTENSITY)
return false;
boolean changeDone = false;
Conformer conformer = null;
Conformer backup = null;
double origCollisionIntensitySum = torsionSet.getCollisionIntensitySum();
for (int f1=1; f1 MIN_ESCAPE_GAIN_PER_STEP
&& mTorsionSetStrategy.getBondsBetweenFragments(f1, f2).length == 2) {
int[] rotatableBond = mTorsionSetStrategy.getBondsBetweenFragments(f1, f2);
int angle = (mRandom.nextDouble() < 0.5) ? -ESCAPE_ANGLE : ESCAPE_ANGLE;
double origCollisionIntensity = origCollisionIntensityMatrix[f1][f2];
for (int i=1; i<=ESCAPE_STEPS; i++) {
// make sure we keep the original coordinates
if (conformer == null)
conformer = new Conformer(origConformer);
else if (backup == null)
backup = new Conformer(conformer);
else
backup.copyFrom(conformer);
for (int bond : rotatableBond)
mRotatableBond[bond].rotateTo(conformer, (short) (conformer.getBondTorsion(bond) + i*angle));
double localCollisionIntensity = calculateCollisionIntensity(conformer, mRigidFragment[f1], mRigidFragment[f2]);
if (localCollisionIntensity < origCollisionIntensity - MIN_ESCAPE_GAIN_PER_STEP) {
origCollisionIntensity = localCollisionIntensity;
changeDone = true;
// not enough collision intensity left for correction
if (localCollisionIntensity < MIN_ESCAPE_GAIN_PER_STEP)
break;
}
else {
if (backup != null)
conformer.copyFrom(backup);
else
conformer.copyFrom(origConformer);
break;
}
}
}
}
}
}
if (!changeDone)
return false;
checkCollision(conformer, torsionSet);
if (torsionSet.getCollisionIntensitySum() >= origCollisionIntensitySum) {
// local improvement causes worsening somewhere else
torsionSet.setCollisionIntensity(origCollisionIntensitySum, origCollisionIntensityMatrix);
return false;
}
origConformer.copyFrom(conformer);
return true;
}
private double calculateCollisionIntensity(Conformer conformer, RigidFragment f1, RigidFragment f2) {
double collisionIntensitySum = 0;
for (int i=0; i=0; j--) {
Conformer ci = new Conformer(conformer);
mRotatableBond[j].rotateToIndex(ci, torsionIndexes[j]);
// System.out.print(torsionIndexes[j]+ " ");
}
//System.out.println("passed collision check! "+mTorsionSet.toString());
separateDisconnectedFragments(conformer);
return conformer;
}
} © 2015 - 2025 Weber Informatics LLC | Privacy Policy