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.calc.ThreadMaster;
import com.actelion.research.chem.*;
import com.actelion.research.chem.conf.Conformer;
import com.actelion.research.chem.conf.TorsionDescriptorHelper;
import com.actelion.research.util.DoubleFormat;
import java.io.BufferedWriter;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.OutputStreamWriter;
import java.nio.charset.StandardCharsets;
import java.util.ArrayList;
import java.util.Random;
public class ConformationSelfOrganizer {
private static final int PHASE_PREPARATION = 0; // without ambiguous rules (i.e. torsion)
private static final int PHASE_PRE_OPTIMIZATION = 1; // with ambiguous rules before breakout
private static final int PHASE_BREAKOUT = 2;
private static final int PHASE_OPTIMIZATION = 3;
private static final int PHASE_MINIMIZATION = 4;
private static final int[] PHASE_CYCLES = { 40, 20, 20, 100, 20 };
private static final double[] PHASE_FACTOR = { 40, 20, 20, 100, 40 };
private static final String[] PHASE_NAME = { "preparation", "pre-optimization", "breakout", "optimization", "minimization" };
private static final int INITIAL_POOL_SIZE = 8;
private static final int MAX_CONFORMER_TRIES = 12;
private static final int MAX_BREAKOUT_ROUNDS = 3;
private static final boolean PREFER_HIGH_STRAIN_RULES = false;
private static final boolean INITIALLY_SKIP_TORSION_RULES = false;
private static final int CONSIDERED_PREVIOUS_CYCLE_COUNT = 10; // must be much smaller than all the ...CYCLES above
private static final double STANDARD_CYCLE_FACTOR = 1.0;
private static final double MINIMIZATION_END_FACTOR = 0.01;
private static final double ATOM_ACCEPTABLE_STRAIN = 2.72; // A conformer is considered acceptable if all atom strains are below this value
private static final double ATOM_FLAT_RING_BREAKOUT_STRAIN = 1000; // TODO check
private static final double ATOM_CAGE_BREAKOUT_STRAIN = 2000; // TODO check Max allowed strain for 4-neighbour atoms
private static final double TWIST_BOAT_ESCAPE_ANGLE = 0.6; // angle to rotate ring member out of plane (from mid point between ring center and atom), after rotating it into the plane from the other side
private static final int TWIST_BOAT_ESCAPE_FREQUENCY = 10; // in every tenth cycle we try escaping trapped twist boats
private static final double BREAKOUT_DISTANCE = 8.0;
private static final double MAX_POOL_STRAIN_DIF = 2.72; // 2 * 1.36 kcal/mol, which is factor 100
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 static BufferedWriter mDWWriter;
private static Conformer mLastDWConformer;
private static int mDWCycle;
private String mDWMode;
private double[] mDWStrain; // TODO get rid of this section
private StereoMolecule mMol;
private Random mRandom;
private int mMaxConformers;
private boolean mPoolIsClosed;
private ArrayList mRuleList;
private ArrayList mConformerList;
private double mMinStrainInPool;
private int[] mRuleCount;
private boolean[] mSkipRule;
private int[] mRotatableBondForDescriptor;
private ThreadMaster mThreadMaster;
/**
* 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;
}
public void setThreadMaster(ThreadMaster tm) {
mThreadMaster = tm;
}
/**
* 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);
SelfOrganizedConformer bestConformer = null;
for (int i=0; i();
mMinStrainInPool = Double.MAX_VALUE;
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 mMinStrainInPool + MAX_POOL_STRAIN_DIF)
return false;
conformer.calculateDescriptor(mRotatableBondForDescriptor);
for (int i=mConformerList.size()-1; i>=0; i--) {
SelfOrganizedConformer soc = mConformerList.get(i);
if (conformer.equals(soc)) {
if (soc.isWorseThan(conformer)) {
mConformerList.remove(i);
mConformerList.add(conformer);
if (mMinStrainInPool > conformer.getTotalStrain())
mMinStrainInPool = conformer.getTotalStrain();
return true;
}
return false;
}
}
mConformerList.add(conformer);
if (mMinStrainInPool > conformer.getTotalStrain()) {
mMinStrainInPool = conformer.getTotalStrain();
for (int i=mConformerList.size()-1; i>=0; i--) {
SelfOrganizedConformer soc = mConformerList.get(i);
if (soc.getTotalStrain() > mMinStrainInPool + MAX_POOL_STRAIN_DIF)
mConformerList.remove(i);
}
}
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.
* Internally, TorsionDescriptors are compared to decide whether a conformer is new.
* These TorsionDescriptors are created on the fly using the default method to determine all
* considered rotatable bonds. Bonds that contain marked atoms are excluded from being considered
* rotatable.
* @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;
}
public void printStrainDetails(SelfOrganizedConformer conformer) {
System.out.println("############ new conformer ###############");
for (ConformationRule rule:mRuleList)
if (rule instanceof DistanceRule)
((DistanceRule)rule).printStrain(conformer);
}
public static void writeDWFileStart() {
mDWCycle = 0;
mLastDWConformer = null;
try {
mDWWriter = new BufferedWriter(new OutputStreamWriter(new FileOutputStream(DATAWARRIOR_DEBUG_FILE), StandardCharsets.UTF_8));
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\tmode\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();
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.close();
mDWWriter = null;
}
catch (IOException ioe) {}
}
/**
* 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;
// Make sure no rule were disabled in previous run.
for (ConformationRule rule:mRuleList)
rule.setEnabled(true);
for (int i=0; i CONSIDERED_PREVIOUS_CYCLE_COUNT) {
double averagePreviousStrain = 0.0;
for (double strain:previousTotalStrains)
averagePreviousStrain += strain;
averagePreviousStrain /= CONSIDERED_PREVIOUS_CYCLE_COUNT;
if (conformer.getTotalStrain() > averagePreviousStrain)
break;
}
if (bestConformer.getTotalStrain() > conformer.getTotalStrain())
bestConformer.copyFrom(conformer);
cycleFactor *= reductionFactor;
}
if (conformer.isWorseThan(bestConformer))
conformer.copyFrom(bestConformer);
}
private boolean isAcceptable(SelfOrganizedConformer conformer) {
for (int atom=0; atom ATOM_ACCEPTABLE_STRAIN)
return false;
return true;
}
private boolean containsTrappedAtom(SelfOrganizedConformer conformer) {
for (int atom=0; atom ATOM_FLAT_RING_BREAKOUT_STRAIN
|| conformer.getAtomStrain(atom) > ATOM_CAGE_BREAKOUT_STRAIN)
return true;
return false;
}
private void writeStrains(SelfOrganizedConformer newConformer, ConformationRule rule, String stepName, double strain1, double strain2) throws Exception {
newConformer.calculateStrain(mRuleList);
double[] strain = new double[ConformationRule.RULE_NAME.length];
double strainSum = 0f;
for (int i=0; i ATOM_FLAT_RING_BREAKOUT_STRAIN
&& tryEscapeFromFlatRingTrap(conformer, atom))
atomCount++;
}
if (atomCount == 0) {
int neighbourCount = 16;
for (int atom=0; atom mMol.getAllConnAtoms(atom)
&& conformer.getAtomStrain(atom) > maxCageBreakoutStrain(atom))
neighbourCount = mMol.getAllConnAtoms(atom);
if (neighbourCount != 16) {
for (int atom=0; atom maxCageBreakoutStrain(atom)) {
//System.out.println("escape "+neighbourCount+" neighbours");
//System.out.print("strains: "); for (int i=0; i
* - 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;
}
private boolean tryEscapeTwistBoats(SelfOrganizedConformer conformer) {
boolean changeDone = false;
RingCollection ringSet = mMol.getRingSet();
for (int r=0; r= 3)) {
int[] notAtom = new int[2];
notAtom[0] = ringAtom[i == 0 ? 5 : i-1];
notAtom[1] = ringAtom[i == 5 ? 0 : i+1];
Coordinates pAtom = conformer.getCoordinates(ringAtom[i]);
Coordinates axis = conformer.getCoordinates(ringAtom[i]).subC(cog).cross(n).unit();
double angle = TWIST_BOAT_ESCAPE_ANGLE + Math.asin(Math.abs(distance[i])/pAtom.distance(cog));
double theta = (distance[i] < 0) ? angle : -angle;
Coordinates center = new Coordinates(pAtom).center(cog);
boolean isBridgedRingAtom = false;
if (mMol.getConnAtoms(ringAtom[i]) > 2) {
for (int j=0; j