org.openmolecules.chem.conf.so.ConformationSelfOrganizer 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.so;
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.*;
import com.actelion.research.util.DoubleFormat;
import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Random;
public class ConformationSelfOrganizer {
private static final int INITIAL_POOL_SIZE = 4;
private static final int MAX_CONFORMER_TRIES = 6;
private static final int MAX_BREAKOUT_ROUNDS = 3;
private static final int PREPARATION_CYCLES = 40; // without ambiguous rules (i.e. torsion)
private static final int PRE_OPTIMIZATION_CYCLES = 20; // with ambiguous rules before breakout
private static final int BREAKOUT_CYCLES = 20;
private static final int OPTIMIZATION_CYCLES = 100;
private static final int MINIMIZATION_CYCLES = 20;
private static final double STANDARD_CYCLE_FACTOR = 1.0;
private static final double MINIMIZATION_REDUCTION = 20.0;
private static final double ATOM_FLAT_RING_BREAKOUT_STRAIN = 0.25;
private static final double ATOM_CAGE_BREAKOUT_STRAIN = 2.0;
private static final double BREAKOUT_DISTANCE = 8.0;
private static final double MAX_AVERAGE_ATOM_STRAIN = 0.025;
private static final double MAX_HIGHEST_ATOM_STRAIN = 0.05;
private static final double MAX_STRAIN_TOLERANCE = 1.5;
public static boolean KEEP_INITIAL_COORDINATES = false;
public static boolean WRITE_DW_FILE = false;
private static final String DATAWARRIOR_DEBUG_FILE = "/home/thomas/data/debug/conformationSampler.dwar";
private BufferedWriter mDWWriter;
private Conformer mLastDWConformer;
private int mDWCycle;
private double[] mDWStrain; // TODO get rid of this
private StereoMolecule mMol;
private Random mRandom;
private int mMaxConformers;
private boolean mPoolIsClosed;
private ArrayList mRuleList;
private ArrayList mConformerList;
private double mMinAverageAtomStrainInPool,mMinHighestAtomStrainInPool;
private int[] mRuleCount;
private boolean[] mSkipRule;
private int[] mRotatableBondForDescriptor;
/**
* Generates a new ConformationSelfOrganizer from the given molecule.
* Explicit hydrogens are removed from the molecule, unless the
* keepHydrogen flag is set.
* One conformer can be generated with the getOneConformer()
* or getOneConformerInPlace().
* Multiple different conformers can be generated with initializeConformers()
* and getNextConformer(). In this case conformers are considered different
* if at least one dihedral angle of a rotatable bond is substantially different.
* If atoms of the molecule are marked, these are not considered part of the molecule,
* when the rotatable bonds for the difference check are located.
* @param mol
*/
public ConformationSelfOrganizer(final StereoMolecule mol, boolean keepHydrogen) {
/*// winkel zwischen zwei vektoren:
final double[] v1 = { Math.sqrt(3.0)/2.0, 0.5, 0 };
final double[] v2 = { -1, 0, 1 };
double cosa = (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
/ (Math.sqrt(v1[0]*v1[0]+v1[1]*v1[1]+v1[2]*v1[2])
* Math.sqrt(v2[0]*v2[0]+v2[1]*v2[1]+v2[2]*v2[2]));
double a = Math.acos(cosa);
System.out.println("angle:"+a+" in degrees:"+(a*180/Math.PI));
*/
mMol = mol;
if (!keepHydrogen)
mMol.removeExplicitHydrogens();
mMol.ensureHelperArrays(Molecule.cHelperParities);
mRuleList = new ArrayList();
mSkipRule = new boolean[ConformationRule.RULE_NAME.length];
mRuleCount = new int[ConformationRule.RULE_NAME.length];
DistanceRule.calculateRules(mRuleList, mol);
mRuleCount[ConformationRule.RULE_TYPE_DISTANCE] = mRuleList.size();
mRuleCount[ConformationRule.RULE_TYPE_PLANE] = -mRuleList.size();
PlaneRule.calculateRules(mRuleList, mol);
mRuleCount[ConformationRule.RULE_TYPE_PLANE] += mRuleList.size();
mRuleCount[ConformationRule.RULE_TYPE_LINE] = -mRuleList.size();
StraightLineRule.calculateRules(mRuleList, mol);
mRuleCount[ConformationRule.RULE_TYPE_LINE] += mRuleList.size();
mRuleCount[ConformationRule.RULE_TYPE_STEREO] = -mRuleList.size();
TetrahedralStereoRule.calculateRules(mRuleList, mol);
mRuleCount[ConformationRule.RULE_TYPE_STEREO] += mRuleList.size();
mRuleCount[ConformationRule.RULE_TYPE_BINAP] = -mRuleList.size();
AxialStereoRule.calculateRules(mRuleList, mol);
mRuleCount[ConformationRule.RULE_TYPE_BINAP] += mRuleList.size();
mRuleCount[ConformationRule.RULE_TYPE_TORSION] = -mRuleList.size();
TorsionRule.calculateRules(mRuleList, mol);
mRuleCount[ConformationRule.RULE_TYPE_TORSION] += mRuleList.size();
// listRules();
}
/* private void listRules() {
System.out.println("---------------------------------------------------------------------");
for (int i=0; i getRuleList() {
return mRuleList;
}
/**
* @return returns the molecule that was passed to the constructor.
*/
public StereoMolecule getMolecule() {
return mMol;
}
/**
* This convenience method returns the StereoMolecule that has been passed
* to the constructor after modifying its atom coordinates
* to reflect the conformer internally created by generateOneConformer().
* @return
*/
public StereoMolecule generateOneConformerInPlace(long randomSeed) {
SelfOrganizedConformer conformer = generateOneConformer(randomSeed);
return (conformer == null) ? null : conformer.toMolecule(mMol);
}
/**
* Generates the coordinates for one conformer in the calling thread.
* This is done by trying MAX_CONFORMER_TRIES times to create a random
* conformer that meets MAX_ATOM_STRAIN and MAX_TOTAL_STRAIN criteria.
* If one is found it is returned. Otherwise the conformer with the lowest
* total strain is returned.
* Note: If randomSeed is different from 0, then only one conformer
* is done produced and returned, no matter whether it meets any strain criteria.
* @param randomSeed 0 or specific seed
*/
public SelfOrganizedConformer generateOneConformer(long randomSeed) {
mRandom = (randomSeed == 0) ? new Random() : new Random(randomSeed);
SelfOrganizedConformer conformer = new SelfOrganizedConformer(mMol);
if (WRITE_DW_FILE) {
try {
writeDWFileStart();
mDWCycle = 0;
mLastDWConformer = null;
tryGenerateConformer(conformer);
writeDWFileEnd();
mDWWriter.close();
return conformer;
}
catch (IOException e) {
e.printStackTrace();
return null;
}
}
SelfOrganizedConformer bestConformer = null;
for (int i=0; i();
mMinHighestAtomStrainInPool = MAX_STRAIN_TOLERANCE * MAX_HIGHEST_ATOM_STRAIN;
mMinAverageAtomStrainInPool = MAX_STRAIN_TOLERANCE * MAX_AVERAGE_ATOM_STRAIN;
mPoolIsClosed = false;
int freeBondCount = 0;
int ringBondCount = 0;
for (int bond=0; bond 1
&& mMol.getAllConnAtoms(mMol.getBondAtom(1, bond)) > 1) {
if (!mMol.isRingBond(bond))
freeBondCount++;
else if (mMol.getBondRingSize(bond) > 4)
ringBondCount++;
}
}
mMaxConformers = (maxConformers == -1) ? (1 << (1+freeBondCount+ringBondCount/2)) : maxConformers;
increaseConformerPool(Math.min(INITIAL_POOL_SIZE, mMaxConformers));
}
private void increaseConformerPool(int newConformerCount) {
SelfOrganizedConformer bestRefusedConformer = null;
SelfOrganizedConformer conformer = null;
int finalPoolSize = mConformerList.size() + newConformerCount;
int tryCount = newConformerCount * MAX_CONFORMER_TRIES;
for (int i=0; i averageStrain)
|| (mMinHighestAtomStrainInPool > highestStrain)) {
if (mMinAverageAtomStrainInPool > averageStrain)
mMinAverageAtomStrainInPool = averageStrain;
if (mMinHighestAtomStrainInPool > highestStrain)
mMinHighestAtomStrainInPool = highestStrain;
for (int j=mConformerList.size()-1; j>=0; j--) {
SelfOrganizedConformer soc = mConformerList.get(j);
if (!soc.isAcceptable(mRuleList)
&& (soc.getTotalStrain() / soc.getSize() > MAX_STRAIN_TOLERANCE * mMinAverageAtomStrainInPool
|| soc.getHighestAtomStrain() > MAX_STRAIN_TOLERANCE * mMinHighestAtomStrainInPool))
mConformerList.remove(j);
}
}
return true;
}
/**
* Picks a new conformer from the conformer pool created by initializeConformers().
* Low strain conformers and conformers being most different from already selected ones
* are selected first. If the pool is getting short on conformers, new conformers are
* created as long as molecule flexibility allows. If a representative set of low strain
* molecules have been picked, this method returns null, provided that at least one conformer
* was returned.
* @return
*/
public SelfOrganizedConformer getNextConformer() {
if (mConformerList == null)
return null;
if (!mPoolIsClosed)
increaseConformerPool(1);
SelfOrganizedConformer bestConformer = null;
for (SelfOrganizedConformer conformer:mConformerList) {
if (!conformer.isUsed()
&& (bestConformer == null || bestConformer.isWorseThan(conformer))) {
bestConformer = conformer;
}
}
if (bestConformer != null)
bestConformer.setUsed(true);
else
mConformerList = null;
return bestConformer;
}
private void writeDWFileStart() throws IOException {
mDWWriter = new BufferedWriter(new FileWriter(DATAWARRIOR_DEBUG_FILE));
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write(" ");
mDWWriter.newLine();
mDWWriter.write("Structure\tbefore\tafter\tcycle\truleName\truleAtoms\truleDetail");
for (int i=0; i");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("\">");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
mDWWriter.write("");
mDWWriter.newLine();
}
/**
* Tries to create one random conformer in current thread by once going through
* sequence of preparation, break-out, optimization and fine-tuning phases.
* Stops successfully early, if none of the atoms' strain exceeds MAX_ATOM_STRAIN
* and the conformer's total strain is below MAX_TOTAL_STRAIN.
* Returns false, if after going though all phases still one of these two conditions
* is not met.
* @param conformer receives coordinates of the conformer
* @return true if conformer with satisfactory low strain could be generated
*/
private boolean tryGenerateConformer(SelfOrganizedConformer conformer) {
if (mMol.getAllAtoms() < 2)
return true;
if (!KEEP_INITIAL_COORDINATES)
jumbleAtoms(conformer);
mSkipRule[ConformationRule.RULE_TYPE_TORSION] = true;
optimize(conformer, PREPARATION_CYCLES, STANDARD_CYCLE_FACTOR, 1.0);
boolean done = false;
if (mRuleCount[ConformationRule.RULE_TYPE_TORSION] != 0) {
mSkipRule[ConformationRule.RULE_TYPE_TORSION] = false;
done = optimize(conformer, PRE_OPTIMIZATION_CYCLES, STANDARD_CYCLE_FACTOR, 1.0);
}
for (int i=0; !done && i ATOM_FLAT_RING_BREAKOUT_STRAIN) {
if (tryEscapeFromFlatRingTrap(conformer, atom)) {
atomCount++;
continue;
}
}
if (atomStrain > ATOM_CAGE_BREAKOUT_STRAIN) {
if (mDWWriter != null) {
try {
writeStrains(conformer, null, "escapeCage", atomStrain, Double.NaN);
}
catch (Exception e) { e.printStackTrace(); }
}
Coordinates c = conformer.getCoordinates(atom);
c.add(BREAKOUT_DISTANCE * mRandom.nextDouble() - BREAKOUT_DISTANCE / 2,
BREAKOUT_DISTANCE * mRandom.nextDouble() - BREAKOUT_DISTANCE / 2,
BREAKOUT_DISTANCE * mRandom.nextDouble() - BREAKOUT_DISTANCE / 2);
atomCount++;
}
}
if (atomCount != 0)
conformer.invalidateStrain();
return atomCount;
}
/**
* Sometimes individual exocyclic atoms end up trapped inside a flat ring,
* because a plane rule combined with a distance rule stabilize the situation.
* This method checks, whether atom
* - has only one neighbour
* - is connected to a small ring atom
* - all angles atom-neighbour-otherRingAtom are below 90 degrees
* If all conditions are true then atom is moved to the opposite side of the neighbour atom.
* @param conformer
* @param atom
* @return whether atom was moved from trapped state
*/
private boolean tryEscapeFromFlatRingTrap(SelfOrganizedConformer conformer, int atom) {
if (mMol.getAllConnAtoms(atom) == 1) {
int connAtom = mMol.getConnAtom(atom, 0);
if (mMol.isSmallRingAtom(connAtom)) {
Coordinates ca = conformer.getCoordinates(atom);
Coordinates cc = conformer.getCoordinates(connAtom);
Coordinates va = cc.subC(ca);
for (int i=0; i Math.PI / 2)
return false;
}
}
ca.add(va).add(va);
if (mDWWriter != null) {
try {
writeStrains(conformer, null, "escapeFlatRing", conformer.getAtomStrain(atom), Double.NaN);
}
catch (Exception e) { e.printStackTrace(); }
}
return true;
}
}
return false;
}
public boolean disableCollidingTorsionRules(SelfOrganizedConformer conformer) {
boolean found = false;
conformer.calculateStrain(mRuleList);
StereoMolecule mol = mMol;
boolean[] isInvolvedAtom = new boolean[mol.getAllAtoms()];
for (ConformationRule rule:mRuleList) {
if (rule instanceof TorsionRule) {
if (((TorsionRule)rule).disableIfColliding(conformer)) {
int[] atom = rule.getAtomList();
for (int i=1; i<=2; i++) {
for (int j=0; j © 2015 - 2025 Weber Informatics LLC | Privacy Policy