com.actelion.research.chem.conf.TorsionDB 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 (c) 1997 - 2016
* Actelion Pharmaceuticals Ltd.
* Gewerbestrasse 16
* CH-4123 Allschwil, Switzerland
*
* All rights reserved.
*
* 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 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 OWNER 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 com.actelion.research.chem.conf;
import com.actelion.research.chem.Molecule;
import com.actelion.research.chem.StereoMolecule;
import java.io.*;
import java.nio.charset.StandardCharsets;
import java.util.TreeMap;
public class TorsionDB {
public static final int MODE_ANGLES = 1;
public static final int MODE_BINS = 2;
public static final int DEFAULT_MERGE_SPAN = 15;
private static final String DATABASE_COD = "cod/";
private static final String DATABASE_CSD = "csd/";
public static final int TORSION_NOT_FOUND = -1;
public static final int TORSION_GREEN = 0;
public static final int TORSION_YELLOW = 1;
public static final int TORSION_RED = 2;
private static final String cBasePath = "/resources/";
private static final String cTorsionIDFile = "torsionID.txt";
private static final String cTorsionAngleFile = "torsionAngle.txt";
private static final String cTorsionRangeFile = "torsionRange.txt";
private static final String cTorsionFrequencyFile = "torsionFrequency.txt";
private static final String cTorsionBinsFile = "torsionBins.txt";
private static final int MAX_SP_CHAIN_LENGTH = 15;
private static TorsionDB sInstance;
private static String sDatabase;
private static String sExternalResourcePath;
private int mSupportedModes;
private int mMergeSpan;
private TreeMap mTreeMap;
/**
* Use this method before initializing the TorsionDB, if you want to open the resource
* files from an external location in the file system rather than from the project resources.
* @param path absolute path of resource file directory ending with a file separator
*/
public static void setResourcePath(String path) {
sExternalResourcePath = path;
}
/**
* Initializes the in memory torsion database by reading torsion list,
* associated torsion frequencies, torsion peak widths (ranges), and/or
* normalized 3-degree bins into memory.
* @param mode combination of MODE_ flags
*/
public static void initialize(int mode) {
initialize(mode, DEFAULT_MERGE_SPAN);
}
/**
* Initializes the in memory torsion database by reading torsion list,
* associated torsion frequencies, torsion peak widths (ranges), and/or
* normalized 3-degree bins into memory.
* If mergeSpan is given, then torsions closer than mergeSpan are merged
* using their frequencies as weights. Frequencies and ranges are adapted
* to new merged and reduced torsion lists.
* @param mode combination of MODE_ flags
* @param mergeSpan 0 or merge span in degrees
*/
public static void initialize(int mode, int mergeSpan) {
synchronized(TorsionDB.class) {
if (sInstance == null)
sInstance = new TorsionDB();
// if the following happens for good reason then we need to support multiple torsionDB instances
if (mode == MODE_ANGLES
&& (sInstance.mSupportedModes & MODE_ANGLES) != 0
&& sInstance.mMergeSpan != mergeSpan)
System.out.println("WARNING: TorsionDB already initialized in ANGLE mode with different merge span.");
int newMode = mode & ~sInstance.mSupportedModes;
if (newMode != 0) {
sInstance.init(newMode);
if (newMode == MODE_ANGLES && mergeSpan != 0)
sInstance.mergeTorsions(mergeSpan);
}
}
}
/**
* Locates all relevant rotatable bonds, i.e. single bonds, which are not in a ring with less than 6 members,
* where both atoms are sp2 or sp3 and carry at least one more non-hydrogen neighbor,
* which are not redundant, and those where a torsion change modifies the relative location of atoms.
* For chains of conjugated triple bonds the following applies:
* If at least one terminal sp2/sp3 atom has no external neighbor, then no single bond is considered rotatable.
* Otherwise that terminal single bond connecting the smaller substituent is considered the only rotatable bond of the linear atom strand.
* Bonds considered rotatable by this method are guaranteed to produce a valid torsionID with getTorsionID().
* @param mol
* @param skipAllRingBonds if true, then considers bonds in small and large(!) rings as non-rotatable
* @param isRotatableBond null or empty(!) array that will be filled by this method
* @return number of rotatable bonds
*/
public static int findRotatableBonds(StereoMolecule mol, boolean skipAllRingBonds, boolean[] isRotatableBond) {
mol.ensureHelperArrays(Molecule.cHelperRings);
if (isRotatableBond == null)
isRotatableBond = new boolean[mol.getBonds()];
int count = 0;
for (int bond=0; bond 1
&& mol.getNonHydrogenNeighbourCount(mol.getBondAtom(1, bond)) > 1
&& !(skipAllRingBonds && mol.isRingBond(bond))
&& !(mol.isSmallRingBond(bond) && mol.getBondRingSize(bond) <= 5)) {
isRotatableBond[bond] = true;
count++;
}
}
int[] centralAtom = new int[2];
int[] rearAtom = new int[2];
boolean[] bondHandled = null;
// For every triple bond (or conjugated chain of triple bonds) keep
// only one terminal single bond as rotatable, more precisely that bond
// that is connected to the smaller molecule part.
for (int bond=0; bond 1
&& mol.getNonHydrogenNeighbourCount(centralAtom[1]) > 1) {
int substituentSize0 = mol.getSubstituentSize(rearAtom[0], centralAtom[0]);
int substituentSize1 = mol.getSubstituentSize(rearAtom[1], centralAtom[1]);
int i = (substituentSize0 < substituentSize1) ? 0 : 1;
int relevantBond = mol.getBond(rearAtom[i], centralAtom[i]);
if (bondHandled == null)
bondHandled = new boolean[mol.getBonds()];
bondHandled[relevantBond] = true;
isRotatableBond[relevantBond] = true;
count++;
}
}
}
}
return count;
}
private TorsionDB() {
mTreeMap = new TreeMap();
mSupportedModes = 0;
}
private void mergeTorsions(int mergeSpan) {
mMergeSpan = mergeSpan;
for (TorsionInfo ti:mTreeMap.values())
ti.mergeTorsions(mergeSpan);
}
private void init(int mode) {
mSupportedModes |= mode;
try {
BufferedReader tr = openReader(cTorsionIDFile);
BufferedReader ar = ((mode & MODE_ANGLES) == 0) ? null : openReader(cTorsionAngleFile);
BufferedReader rr = ((mode & MODE_ANGLES) == 0) ? null : openReader(cTorsionRangeFile);
BufferedReader fr = ((mode & MODE_ANGLES) == 0) ? null : openReader(cTorsionFrequencyFile);
BufferedReader br = ((mode & MODE_BINS) == 0) ? null : openReader(cTorsionBinsFile);
String type = tr.readLine();
while (type != null) {
TorsionInfo torsionInfo = mTreeMap.get(type);
if (torsionInfo == null) {
torsionInfo = new TorsionInfo(getSymmetryType(type));
mTreeMap.put(type, torsionInfo);
}
if (ar != null) {
String[] angle = ar.readLine().split(",");
torsionInfo.angle = new short[angle.length];
for (int i=0; i' indicates
* a non-inverted stereo center.
* If a TorsionDetail is given, then this will be filled with all classification information.
* If an torsionAtom array is given instead, then this is filled with the four torsion defining atoms.
* Terminal torsion atoms in torsionAtom (indexes 0,3) may be -1, which indicates a virtual atom
* assumed to lie on the third terminal sp3 position, when we have two real terminal atoms with
* equal symmetry ranks. If bond is part of a consecutive sp-sp atom chain, then
* the classifying fragment covers all linear atoms plus two end atoms not being member
* of the linear sp-atom strand. In this case torsionAtom (indexes 1,2) contains the two
* end-atoms of the linear atom strand, which are sp2- or sp3-hybridized.
* @param mol source molecule
* @param bond of which to determine torsion statistics
* @param torsionAtom null or int[4] to receive torsion atoms
* @param detail null or empty TorsionDetail to receive torsion fragment, atom map, etc.
* @return null or idcode identifying the bond and its vicinity
*/
public static String getTorsionID(StereoMolecule mol, int bond, int[] torsionAtom, TorsionDetail detail) {
if (detail == null)
detail = new TorsionDetail();
detail.classify(mol, bond);
if (torsionAtom != null) {
torsionAtom[0] = detail.getReferenceAtom(0);
torsionAtom[1] = detail.getCentralAtom(0);
torsionAtom[2] = detail.getCentralAtom(1);
torsionAtom[3] = detail.getReferenceAtom(1);
}
return detail.getID();
}
/**
* A fragment's symmetry type defines whether one part of a full 360 degree
* matches symmetrically another range.
* SYMMETRY_D1D1 (0 -> 180 match 0 -> -180): Two times D1 symmetry, i.e. one
* terminal neighbor or two symmetrical sp3 neighbors.
* SYMMETRY_C1D2 (0 -> 180 match -180 -> 0): One stereo center on one side
* and two symmetrical sp2 neighbors on the other.
* SYMMETRY_D1D2_OR_D2D2 (0 -> 90 match 180 -> 90, 0 -> -90, -180 -> -90):
* Two symmetrical sp2 neighbors on one side and at least D1 (mirror plane) on the other.
* SYMMETRY_C1C1_OR_C1D1 (no symmetry): one stereo center on one side and
* not more than D1 (mirror plane) on the other.
* This is a quick lookup, because this information is encoded in the torsionID itself.
* @param torsionID
* @return
*/
public static int getSymmetryType(String torsionID) {
return torsionID.endsWith("<") || torsionID.endsWith(">") ? TorsionDetail.SYMMETRY_C1C1_OR_C1D1
: torsionID.endsWith("-") || torsionID.endsWith("+") ? TorsionDetail.SYMMETRY_C1D2
: torsionID.endsWith("=") ? TorsionDetail.SYMMETRY_D1D2_OR_D2D2
: TorsionDetail.SYMMETRY_D1D1;
}
public static boolean isInverted(String torsionID) {
return torsionID.endsWith("<") || torsionID.endsWith("-");
}
public static String normalizeID(String torsionID) {
return torsionID.endsWith("<") ? torsionID.substring(0, torsionID.length()-1)+">"
: torsionID.endsWith("-") ? torsionID.substring(0, torsionID.length()-1)+"+"
: torsionID;
}
/**
* Calculates a signed torsion like calculateTorsion(). However, terminal atoms
* of the defined atom sequence may be set to -1, referring to a hydrogen
* or electron pair in the 3rd terminal sp2 position, when the other two
* positions are occupied by two atoms that share the same symmetry rank.
* @param mol
* @param atom valid 4-atom sequence with terminal atoms optionally set to -1
* @return torsion in the range: -pi <= torsion <= pi
*/
public static double calculateTorsionExtended(StereoMolecule mol, int[] atom) {
if (atom[0] != -1 && atom[3] != -1)
return mol.calculateTorsion(atom);
return calculateTorsionExtended(mol, null, atom);
}
/**
* Calculates a signed torsion like calculateTorsion(). However, terminal atoms
* of the defined atom sequence may be set to -1, referring to a hydrogen
* or electron pair in the 3rd terminal sp2 position, when the other two
* positions are occupied by two atoms that share the same symmetry rank.
* @param conformer
* @param atom valid 4-atom sequence with terminal atoms optionally set to -1
* @return torsion in the range: -pi <= torsion <= pi (NaN, e.g. if 3 atoms on straight line)
*/
public static double calculateTorsionExtended(Conformer conformer, int[] atom) {
if (atom[0] != -1 && atom[3] != -1)
return conformer.calculateTorsion(atom);
return calculateTorsionExtended(conformer.getMolecule(), conformer, atom);
}
/**
* @param mol used for connectivity info, used for coords if conformer==null
* @param conformer provides coords or is null
* @param atom
* @return
*/
private static double calculateTorsionExtended(StereoMolecule mol, Conformer conformer, int[] atom) {
double[] angle1 = new double[2];
int[] rearAtom = findRearAtoms(mol, atom);
for (int i=0; i<2; i++) {
if (atom[3*i] != -1) {
int central = 2-i;
int terminal = 3-3*i;
assert(mol.getNonHydrogenNeighbourCount(atom[central]) == 3);
int index = 0;
for (int j=0; j<3; j++) {
int connAtom = mol.getConnAtom(atom[central], j);
if (connAtom != rearAtom[1-i] && mol.getAtomicNo(connAtom) != 1) {
// if (connAtom != atom[3-central]) {
atom[terminal] = connAtom;
if (conformer != null)
angle1[index++] = conformer.calculateTorsion(atom);
else
angle1[index++] = mol.calculateTorsion(atom);
}
}
atom[terminal] = -1; // restore value for virtual atom
return calculateVirtualTorsion(angle1);
}
}
double[] angle2 = new double[2];
int index1 = 0;
assert(mol.getNonHydrogenNeighbourCount(atom[1]) == 3);
for (int i=0; i<3; i++) {
int terminal1 = mol.getConnAtom(atom[1], i);
if (terminal1 != rearAtom[0] && mol.getAtomicNo(terminal1) != 1) {
// if (terminal1 != atom[2]) {
assert(mol.getNonHydrogenNeighbourCount(atom[2]) == 3);
atom[0] = terminal1;
int index2 = 0;
for (int j=0; j<3; j++) {
int terminal2 = mol.getConnAtom(atom[2], j);
if (terminal2 != rearAtom[1] && mol.getAtomicNo(terminal2) != 1) {
// if (terminal2 != atom[1]) {
atom[3] = terminal2;
if (conformer != null)
angle2[index2++] = conformer.calculateTorsion(atom);
else
angle2[index2++] = mol.calculateTorsion(atom);
}
}
angle1[index1++] = calculateVirtualTorsion(angle2);
}
}
atom[0] = -1;
atom[3] = -1;
return calculateVirtualTorsion(angle1);
}
/**
* If the atom sequence contains a straight chain of sp-hybridized atoms,
* then the atom array contain the first two and last two atoms of a sequence
* of 2n+4 atoms (n: number of conjugated triple bonds). This method finds
* the 3rd and 3rd-last atoms of the chain. rearAtom[0] is connected to atom[1]
* and rearAtom[1] is connected to atom[2].
* @param mol
* @param torsionAtom
* @return
*/
public static int[] findRearAtoms(StereoMolecule mol, int[] torsionAtom) {
int[] rearAtom = new int[2];
if (mol.getBond(torsionAtom[1], torsionAtom[2]) != -1) {
rearAtom[0] = torsionAtom[2];
rearAtom[1] = torsionAtom[1];
}
else {
int[] pathAtom = new int[MAX_SP_CHAIN_LENGTH+1];
int pathLength = mol.getPath(pathAtom, torsionAtom[1], torsionAtom[2], MAX_SP_CHAIN_LENGTH, null);
rearAtom[0] = pathAtom[1];
rearAtom[1] = pathAtom[pathLength-1];
}
return rearAtom;
}
/**
* If the atom sequence contains a straight chain of sp-hybridized atoms,
* then the atom array contain the first two and last two atoms of a sequence
* of 2n+4 atoms (n: number of conjugated triple bonds). This method finds
* extends the chain by inserting all sp-hybridised atoms into the strand.
* @param mol
* @param atom torsion defining (1st two and last two) atoms of atom sequence
* @return 2n+4 atom indexes (n: number of conjugated triple bonds)
*/
public static int[] getExtendedAtomSequence(StereoMolecule mol, int[] atom) {
if (mol.getBond(atom[1], atom[2]) != -1)
return atom;
int[] pathAtom = new int[MAX_SP_CHAIN_LENGTH+1];
int pathLength = mol.getPath(pathAtom, atom[1], atom[2], MAX_SP_CHAIN_LENGTH, null);
int[] realAtom = new int[3+pathLength];
realAtom[0] = atom[0];
for (int i=0; i<=pathLength; i++)
realAtom[i+1] = pathAtom[i];
realAtom[pathLength+2] = atom[3];
return realAtom;
}
private static double calculateVirtualTorsion(double[] angle) {
double meanAngle = (angle[1] + angle[0]) / 2;
double angleDif = Math.abs(angle[1] - angle[0]);
if (angleDif > Math.PI)
return meanAngle;
return (meanAngle < 0) ? meanAngle + Math.PI : meanAngle - Math.PI;
}
/**
* Checks, whether the torsion angle lies within statistical limits.
* The torsion angle used should be obtained from a call to calculateTorsionExtended()
* in order to be canonical concerning to the chosen terminal (pseudo)atoms.
* Requires TorsionDB to be initialized with MODE_RANGES.
* @param torsionID as obtained by getTorsionID()
* @param angle obtained from a call to calculateTorsionExtended()
* @return TORSION_GREEN,TORSION_YELLOW,TORSION_RED or TORSION_NOT_FOUND
*/
public static int getTorsionStrainClass(String torsionID, double angle) {
short[][] rangeList = getTorsionRanges(normalizeID(torsionID));
if (rangeList == null)
return TORSION_NOT_FOUND;
int strain = TORSION_RED;
int angleIndex = getNormalizedTorsionIndex(torsionID, angle);
for (short[] range:rangeList) {
if (angleIndex >= range[0] && angleIndex <= range[1])
return TORSION_GREEN;
if (angleIndex >= range[0]-5 && angleIndex <= range[1]+5)
strain = TORSION_YELLOW;
}
return strain;
}
/**
* Normalizes a torsion angle considering the fragments symmetry type
* by returning the lowest symmetrically equivalent torsion that
* is >= 0.
* @param torsionID
* @param angle
* @return angle within native range of symmetry type
*/
public static int getNormalizedTorsionIndex(String torsionID, double angle) {
if (isInverted(torsionID))
angle = -angle;
// indices range from -180 to 179 indicating an angle = index
angle = 180.0 * angle / Math.PI;
int index = (int)Math.floor(0.5 + angle);
if (index == 180)
index = -180;
int symmetryClass = getSymmetryType(torsionID);
switch (symmetryClass) {
case TorsionDetail.SYMMETRY_C1C1_OR_C1D1:
if (index < 0)
index += 360;
break;
case TorsionDetail.SYMMETRY_C1D2:
if (index < 0)
index += 180;
break;
case TorsionDetail.SYMMETRY_D1D1:
index = Math.abs(index);
break;
case TorsionDetail.SYMMETRY_D1D2_OR_D2D2:
if (index < 0)
index += 180;
if (index > 90)
index = 180 - index;
break;
}
return index;
}
}
class TorsionInfo {
short[] angle;
short[][] range;
short[] frequency;
byte[] binSize;
int symmetryClass;
public TorsionInfo(int symmetryClass) {
this.symmetryClass = symmetryClass;
}
/**
* Create TorsionInfo describing a fragment with inverted stereo situation.
* @param ti TorsionInfo of non-symmetrical fragment
*/
public TorsionInfo(TorsionInfo ti) {
symmetryClass = ti.symmetryClass;
if (ti.angle != null) {
angle = new short[ti.angle.length];
for (int i=0; i= 0) ? angle[i1] : -angle[0];
int angle2 = (i2 < angle.length) ? angle[i2]
: (symmetryClass == TorsionDetail.SYMMETRY_C1D2) ? angle[0] + 180
: (symmetryClass == TorsionDetail.SYMMETRY_D1D1) ? 360 - angle[angle[i1] == 180 ? i1-1 : i1]
: (symmetryClass == TorsionDetail.SYMMETRY_D1D2_OR_D2D2) ? 180 - angle[angle[i1] == 90 ? i1-1 : i1] : angle[0] + 360;
if (minAngleDif > angle2 - angle1) {
minAngleDif = angle2 - angle1;
minIndex = i1;
}
}
if (minAngleDif > mergeSpan)
break;
int i1 = minIndex;
int i2 = i1+1;
if (i1 == -1) { // angle[0] != 0: we have a symmetrical peak at -angle[0]
additionalFrequency += frequency[0];
angle[0] = 0;
frequency[0] *= 2;
range[0][0] = (short)-range[0][1];
}
else if (i2 < angle.length) {
int mergedAngle = mergedAngle(i1, i2, angle[i2]);
mergeAngles(mergedAngle, i1, i2, range[i1][0], range[i2][1]);
}
else if (symmetryClass == TorsionDetail.SYMMETRY_C1D2) {
i2 = 0;
int mergedAngle = mergedAngle(i1, i2, angle[i2] + 180);
if (mergedAngle<=180)
mergeAngles(mergedAngle, i1, i2, range[i1][0], range[i2][1] + 180);
else
mergeAngles(mergedAngle - 180, i2, i1, range[i2][0], range[i1][1] - 180);
}
else if (symmetryClass == TorsionDetail.SYMMETRY_D1D1) {
i2 = angle[i1] == 180 ? i1-1 : i1;
if (angle[i1] < 180) { // we merge with a symmetrical peak
additionalFrequency += frequency[i1];
angle[i1] = 180;
frequency[i1] *= 2;
range[i1][1] = (short)(360 - range[i1][0]);
}
else { // should not happen, because we should have caught the symmetrical situation one peak earlier first
System.out.println("WARNING: TorsionDB - Should have caught the symmetrical situation one peak earlier.");
int mergedAngle = mergedAngle(i1, i2, angle[i2]);
mergeAngles(mergedAngle, i2, i1, range[i2][0], range[i1][1]);
}
}
else if (symmetryClass == TorsionDetail.SYMMETRY_D1D2_OR_D2D2) {
i2 = angle[i1] == 90 ? i1-1 : i1;
if (angle[i1] < 90) { // we merge with a symmetrical peak
additionalFrequency += frequency[i1];
angle[i1] = 90;
frequency[i1] *= 2;
range[i1][1] = (short)(180 - range[i1][0]);
}
else { // should not happen, because we should have caught the symmetrical situation one peak earlier first
System.out.println("WARNING: TorsionDB - Should have caught the symmetrical situation one peak earlier.");
int mergedAngle = mergedAngle(i1, i2, angle[i2]);
mergeAngles(mergedAngle, i2, i1, range[i2][0], range[i1][1]);
}
}
else { // TorsionDetail.SYMMETRY_C1C1_OR_C1D1
i2 = 0;
int mergedAngle = mergedAngle(i1, i2, angle[i2] + 360);
if (mergedAngle<=360)
mergeAngles(mergedAngle, i1, i2, range[i1][0], range[i2][1] + 360);
else
mergeAngles(mergedAngle - 360, i2, i1, range[i2][0], range[i1][1] - 360);
}
}
if (additionalFrequency != 0)
for (int i=0; i