org.openmolecules.chem.conf.so.DistanceRule 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.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