org.openmolecules.chem.conf.so.DistanceRule Maven / Gradle / Ivy
/*
* 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.Molecule;
import com.actelion.research.chem.RingCollection;
import com.actelion.research.chem.StereoMolecule;
import com.actelion.research.chem.conf.BondAngleSet;
import com.actelion.research.chem.conf.BondLengthSet;
import com.actelion.research.chem.conf.Conformer;
import com.actelion.research.chem.conf.VDWRadii;
import com.actelion.research.util.DoubleFormat;
import java.util.ArrayList;
public class DistanceRule extends ConformationRule {
private static final double VDW_RADIUS_CORRECTION = 1.02; // We increase reported VDW radii to increase Lennard-Jones energies
private static final int PRIORITY_ONE_BOND = 10; // We also use priorities to distinguish cases
private static final int PRIORITY_TWO_BONDS = 5;
private static final int PRIORITY_THREE_BONDS = 3;
private static final int PRIORITY_FOUR_AND_MORE_BONDS = 1;
private static final int PRIORITY_DISCONNECTED = 0;
private double[] mDistance;
private int[] mNotList;
private int mPriority;
/**
* Creates a dummy rule to be neglected
*/
public DistanceRule() {
super(null);
}
/**
* Constructor for 2 bonds in between.
* @param atom
* @param notList atoms (direct neighbors of atom[0] and atom[1]) to be excluded from movement
* @param distance
*/
public DistanceRule(int[] atom, int[] notList, double distance, int priority) {
super(atom);
mDistance = new double[1];
mDistance[0] = distance;
mNotList = notList;
mPriority = priority;
}
public DistanceRule(int[] atom, double minDistance, double maxDistance, int priority) {
super(atom);
mDistance = new double[2];
mDistance[0] = minDistance;
mDistance[1] = maxDistance;
mPriority = priority;
}
public DistanceRule(int[] atom, int[] notList, double minDistance, double maxDistance, int priority) {
super(atom);
mDistance = new double[2];
mDistance[0] = minDistance;
mDistance[1] = maxDistance;
mNotList = notList;
mPriority = priority;
}
@Override
public int getRuleType() {
return RULE_TYPE_DISTANCE;
}
public boolean isFixedDistance() {
return mDistance.length == 1;
}
public static void calculateRules(ArrayList ruleList, StereoMolecule mol) {
BondLengthSet bondLengthSet = new BondLengthSet(mol);
BondAngleSet bondAngleSet = new BondAngleSet(mol, bondLengthSet);
DistanceRule[][] rule = new DistanceRule[mol.getAllAtoms()][];
for (int i=1; i 1
&& mol.getAllConnAtoms(atom[1]) > 1) {
// triple bonds
if (mol.getBondOrder(bond) == 3) {
double distance = bondLengthSet.getLength(bond);
int[] outerAtom = new int[2];
for (int i = 0; i < 2; i++) {
for (int j = 0; j < mol.getAllConnAtoms(atom[i]); j++) {
int connBond = mol.getConnBond(atom[i], j);
if (connBond != bond) {
distance += bondLengthSet.getLength(connBond);
outerAtom[i] = mol.getConnAtom(atom[i], j);
break;
}
}
}
int[] notAtom = new int[2];
notAtom[0] = atom[0];
notAtom[1] = atom[1];
setFixedDistance(rule, combineAtoms(outerAtom[0], outerAtom[1]), notAtom, distance, PRIORITY_THREE_BONDS);
}
else if (mol.getBondOrder(bond) == 2
&& mol.getAtomPi(atom[0]) == 1
&& mol.getAtomPi(atom[1]) == 1
&& mol.getBondParity(bond) != Molecule.cBondParityUnknown) {
// strainless double bond with stereo information
// (including symmetrical ones with parityNone)
int[][] connAtom = new int[2][];
int[][] connBond = new int[2][];
double[][] connAngle = new double[2][];
for (int i=0; i<2; i++) {
connAtom[i] = new int[mol.getAllConnAtoms(atom[i])-1];
connBond[i] = new int[mol.getAllConnAtoms(atom[i])-1];
connAngle[i] = new double[mol.getAllConnAtoms(atom[i])-1];
int doubleBondOpponentIndex = -1;
for (int j=0; j connAtom[0][1-i])
isE = !isE;
if (connAtom[1].length == 2 && connAtom[1][j] > connAtom[1][1-j])
isE = !isE;
setDoubleBondDistance(connAtom[0][i], connAtom[1][j],
connBond[0][i], connBond[1][j],
connAngle[0][i], connAngle[1][j],
bond, isE, bondLengthSet, rule, mol);
}
}
}
else if (mol.getBondOrder(bond) != 0) {
int[] opponentIndex = new int[2];
for (int i = 0; i < 2; i++) {
for (int j = 0; j < mol.getAllConnAtoms(atom[i]); j++) {
if (mol.getConnAtom(atom[i], j) == atom[1 - i]) {
opponentIndex[i] = j;
break;
}
}
}
for (int i = 0; i < mol.getAllConnAtoms(atom[0]); i++) {
if (i != opponentIndex[0]) {
for (int j = 0; j < mol.getAllConnAtoms(atom[1]); j++) {
if (j != opponentIndex[1]) {
if (mol.getAtomPi(atom[0]) == 0
&& mol.getAtomPi(atom[1]) == 0)
setSingleBondConnAtomDistance(
mol.getConnAtom(atom[0], i), mol.getConnAtom(atom[1], j),
mol.getConnBond(atom[0], i), mol.getConnBond(atom[1], j),
bondAngleSet.getConnAngle(atom[0], opponentIndex[0], i),
bondAngleSet.getConnAngle(atom[1], opponentIndex[1], j),
bond, bondLengthSet, rule, mol);
else
setAnyBondConnAtomDistance(
mol.getConnAtom(atom[0], i), mol.getConnAtom(atom[1], j),
mol.getConnBond(atom[0], i), mol.getConnBond(atom[1], j),
bondAngleSet.getConnAngle(atom[0], opponentIndex[0], i),
bondAngleSet.getConnAngle(atom[1], opponentIndex[1], j),
bond, bondLengthSet, rule, mol);
}
}
}
}
}
}
}
// distances over 4 bonds in allenes
for (int atom=0; atom 2) {
distanceToRoot[candidate] = distanceToRoot[parent]
+ bondLengthSet.getLength(mol.getConnBond(parent, i));
if (candidate < rootAtom && rule[rootAtom][candidate] == null) {
if (bondCount[candidate] == 3) {
rule[rootAtom][candidate] = new DistanceRule(); // add dummy rule to block this atom combination
}
else {
int[] notList = new int[2];
notList[0] = graphAtom[1];
notList[1] = parent;
rule[rootAtom][candidate] = new DistanceRule(combineAtoms(rootAtom, candidate), notList,
getVDWRadius(rootAtom, mol) + getVDWRadius(candidate, mol), distanceToRoot[candidate], 0);
}
}
}
}
}
current++;
}
}
private static void calculateDisconnectedDistanceRules(DistanceRule[][] rule, StereoMolecule mol) {
for (int atom1=1; atom1 atom2) {
atom[0] = atom1;
atom[1] = atom2;
}
else {
atom[0] = atom2;
atom[1] = atom1;
}
return atom;
}
private static int[] calculateBondRingSizes(StereoMolecule mol) {
int[] bondRingSize = new int[mol.getBonds()];
RingCollection ringSet = mol.getRingSet();
for (int ring=0; ring ringBond.length) {
bondRingSize[ringBond[i]] = ringBond.length;
}
}
}
return bondRingSize;
}
@Override
public boolean apply(Conformer conformer, double cycleFactor) {
double dx = conformer.getX(mAtom[1]) - conformer.getX(mAtom[0]);
double dy = conformer.getY(mAtom[1]) - conformer.getY(mAtom[0]);
double dz = conformer.getZ(mAtom[1]) - conformer.getZ(mAtom[0]);
double distance = Math.sqrt(dx*dx+dy*dy+dz*dz);
double distanceFactor = 0.0f;
if (mDistance.length == 2) { // is min and max
if (distance < mDistance[0]) {
distanceFactor = (distance-mDistance[0]) / distance;
}
else if (distance > mDistance[1]) {
distanceFactor = (distance-mDistance[1]) / distance;
}
}
else { // exact distance
if (distance < mDistance[0]) {
distanceFactor = (distance-mDistance[0]) / distance;
}
else if (distance > mDistance[0]) {
distanceFactor = (distance-mDistance[0]) / distance;
}
}
if (Math.abs(distanceFactor) < 0.001)
return false;
double factor = cycleFactor * distanceFactor;
StereoMolecule mol = conformer.getMolecule();
if (mPriority == PRIORITY_ONE_BOND) {
if (mol.getAllConnAtoms(mAtom[0]) == 1
&& mol.getAllConnAtoms(mAtom[1]) != 1) {
conformer.getCoordinates(mAtom[0]).add(dx*factor, dy*factor, dz*factor);
return true;
}
if (mol.getAllConnAtoms(mAtom[0]) != 1
&& mol.getAllConnAtoms(mAtom[1]) == 1) {
conformer.getCoordinates(mAtom[1]).add(-dx*factor, -dy*factor, -dz*factor);
return true;
}
}
factor /= 2f;
moveGroup(conformer, mAtom[0], mNotList, dx*factor, dy*factor, dz*factor);
moveGroup(conformer, mAtom[1], mNotList, -dx*factor, -dy*factor, -dz*factor);
return true;
}
@Override
public double addStrain(Conformer conformer, double[] atomStrain) {
double strain = getStrain(conformer);
if (atomStrain != null && strain > 0.0) {
atomStrain[mAtom[0]] += strain / 2;
atomStrain[mAtom[1]] += strain / 2;
}
return strain;
}
private double getStrain(Conformer conformer) {
double distance = conformer.getCoordinates(mAtom[1]).distance(conformer.getCoordinates(mAtom[0]));
if (mDistance.length == 2) {
if (distance < VDW_RADIUS_CORRECTION * mDistance[0]) {
return calculateVDWStrain(conformer.getMolecule(), distance);
}
else if (distance > mDistance[1]) {
return calculateRelaxedStrain((distance - mDistance[1]) / distance);
}
}
else {
double strain = (mPriority == PRIORITY_ONE_BOND) ?
calculateDirectConnectionStrain(Math.abs(mDistance[0] - distance) / Math.max(mDistance[0], distance))
: (mPriority == PRIORITY_TWO_BONDS) ?
calculateTwoBondStrain(Math.abs(mDistance[0] - distance) / Math.max(mDistance[0], distance))
: calculateRelaxedStrain(Math.abs(mDistance[0] - distance) / Math.max(mDistance[0], distance));
if (Math.abs(strain) > 0.01f)
return strain;
}
return 0.0;
}
public void printStrain(Conformer conformer) {
double distance = conformer.getCoordinates(mAtom[1]).distance(conformer.getCoordinates(mAtom[0]));
double strain = 0.0;
double must = 0.0;
if (mDistance.length == 2) {
if (distance < VDW_RADIUS_CORRECTION * mDistance[0]) {
must = mDistance[0];
strain = calculateVDWStrain(conformer.getMolecule(), distance);
}
else if (distance > mDistance[1]) {
must = mDistance[1];
strain = calculateRelaxedStrain((distance - mDistance[1]) / distance);
}
}
else {
must = mDistance[0];
strain = (mPriority == PRIORITY_ONE_BOND) ?
calculateDirectConnectionStrain(Math.abs(mDistance[0] - distance) / Math.max(mDistance[0], distance))
: (mPriority == PRIORITY_TWO_BONDS) ?
calculateTwoBondStrain(Math.abs(mDistance[0] - distance) / Math.max(mDistance[0], distance))
: calculateRelaxedStrain(Math.abs(mDistance[0] - distance) / Math.max(mDistance[0], distance));
}
if (strain > 0.001) {
System.out.print("Atoms("+Molecule.cAtomLabel[conformer.getMolecule().getAtomicNo(mAtom[0])]+mAtom[0]
+(mPriority == PRIORITY_ONE_BOND? "-" : mPriority == PRIORITY_TWO_BONDS? "-?-" : mPriority == PRIORITY_THREE_BONDS? "-?-?-" : "-...-")
+Molecule.cAtomLabel[conformer.getMolecule().getAtomicNo(mAtom[1])]+mAtom[1]+") distance:"+DoubleFormat.toString(distance, 4));
System.out.println(" must:"+DoubleFormat.toString(must,4)+" strain:"+DoubleFormat.toString(strain, 4));
}
}
private double calculateDirectConnectionStrain(double deltaDistance) {
// We use a quadratic potential with an energy penalty of 100.0 kcal/mol when 10% off
return 10000.0 * deltaDistance * deltaDistance;
}
private double calculateTwoBondStrain(double deltaDistance) {
// We use a quadratic potential with an energy penalty of 50 kcal/mol when 10% off
return 8000.0 * deltaDistance * deltaDistance;
}
private double calculateRelaxedStrain(double deltaDistance) {
// We use a quadratic potential with an energy penalty of 50 kcal/mol when 10% off
return 4000.0 * deltaDistance * deltaDistance;
}
private double calculateVDWStrain(StereoMolecule mol, double distance) {
// We use the repulsion part of the Lennard-Jones potential
double vdwradii = VDW_RADIUS_CORRECTION * getVDWRadius(mAtom[1], mol) + getVDWRadius(mAtom[0], mol);
double reldist = distance / vdwradii;
double reldist6 = Math.pow(reldist, -6);
double constant = 2.0f; // constant=1.0 causes a return value of 1.66 kcal/mol with reldist=0.9
return (reldist >= 1.0) ? 0.0 : constant * (reldist6 * reldist6 - reldist6);
}
@Override
public String toString() {
StringBuilder sb = new StringBuilder("distance rule:");
super.addAtomList(sb);
if (mDistance.length == 1)
sb.append(" distance:"+DoubleFormat.toString(mDistance[0]));
else
sb.append(" min:"+DoubleFormat.toString(mDistance[0])+" max:"+DoubleFormat.toString(mDistance[1]));
if (mNotList != null) {
sb.append(" not:"+mNotList[0]);
for (int i=1; i © 2015 - 2025 Weber Informatics LLC | Privacy Policy