com.actelion.research.chem.io.pdb.converter.BondsCalculator Maven / Gradle / Ivy
package com.actelion.research.chem.io.pdb.converter;
/*
* Copyright 2016 Actelion Pharmaceuticals Ltd., Gewerbestrasse 16, CH-4123 Allschwil, Switzerland
*
* This file is part of cif2molecule.
*
* cif2molecule is free software: you can redistribute it and/or modify it under the terms of the
* GNU General Public License as published by the Free Software Foundation, either version 3 of
* the License, or (at your option) any later version.
*
* cif2molecule is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
* See the GNU General Public License for more details.
* You should have received a copy of the GNU General Public License along with DataWarrior.
* If not, see http://www.gnu.org/licenses/.
*
* @author Antanas Vaitkus
*/
import com.actelion.research.chem.*;
import com.actelion.research.chem.conf.BondLengthSet;
import com.actelion.research.util.IntQueue;
import java.util.ArrayList;
import java.util.HashSet;
import java.util.LinkedList;
import java.util.List;
import java.util.Set;
import static com.actelion.research.chem.conf.VDWRadii.COVALENT_RADIUS;
/**
* BondsCalculator is used to recreate the bonds and / or calculate the bonds orders
* based on the 3D coordinates of the atoms
*/
public class BondsCalculator {
private static final boolean DEBUG = false;
private static final double BUMP_FACTOR = 0.5;
private static final double SPECIAL_POSITION_CUTOFF = 0.01;
/**
* Calculates the bonds of a molecule by checking the distance between
* all atoms. The bond order is not set with this function.
*
* Complexity O(nAtoms)
* Memory O(nAtoms)
*
* @param mol
* The molecule for which the bonds should be calculated.
* @param lenient
* Suppresses exceptions that would normally be raised
* due to chemical discrepancies found in the molecule.
* @throws Exception
* Exception containing information about the exceeded
* maximum valence of an atom.
*/
public static void createBonds(Molecule3D mol, boolean lenient) throws Exception {
if ( mol.getAllAtoms()==0 ) return;
// Processing organic atoms first:
// 1. Create a grid
MoleculeGrid grid = new MoleculeGrid(mol);
int[] neighborCount = new int[mol.getAllAtoms()];
// 2. For each atom, check the neighbours and create
// a connection if the distance is close to the sum
// of covalent radii.
for (int i = 0; i < mol.getAllAtoms(); i++) {
if (!mol.isOrganicAtom(i)) continue;
// Get the neighbours
Set set = grid.getNeighbours(mol, i, 3.2);
for (int j : set) {
if (i > j) continue;
if (!mol.isOrganicAtom(j)) continue;
// Two hydrogen (H) atoms are very unlikely to form bonds in crystal structures
if (mol.getAtomicNo(i) == 1 && mol.getAtomicNo(j) == 1) continue;
double dist = getBondLength(mol, i, j);
double idealDist = COVALENT_RADIUS[mol.getAtomicNo(i)] + COVALENT_RADIUS[mol.getAtomicNo(j)];
if (dist > idealDist + .45) continue;
if (neighborCount[i] >= maxNeighborCount(mol, i) ||
neighborCount[j] >= maxNeighborCount(mol, j)) {
int hypervalentAtom = (neighborCount[i] >= maxNeighborCount(mol, i)) ? i : j;
String error = "maximum valence of " + maxNeighborCount(mol, hypervalentAtom) +
" exceeded for atom '" + mol.getAtomName(hypervalentAtom) + "'";
if (!lenient) {
throw new BondCalculatorException(error);
} else {
System.err.println();
}
continue;
}
try {
mol.addBond(i, j, 1);
neighborCount[i]++;
neighborCount[j]++;
} catch (Exception e) {
if (!lenient) throw e;
}
}
}
mol.ensureHelperArrays(Molecule.cHelperRings);
// Deleting incorrectly assigned bonds
for (int i = 0; i < mol.getRingSet().getSize(); i++) {
int[] ringAtoms = mol.getRingSet().getRingAtoms(i);
if (ringAtoms.length == 3) {
if (markImprobableDiphosphorusBondForDeletionInThreeAtomRing(mol, i)) continue;
markImprobableBondWithHydrogenForDeletionInThreeAtomRing(mol, i);
} else if (ringAtoms.length == 4) {
// Small rings sometimes get assigned incorrect interring bonds.
// For example, a 4-membered ring with atom connectivity a1-a2,
// a2-a3, a3-a4, a4-a1 sometimes get assigned an additional
// a1-a3 (or a2-a4) bond.
// The ring must be more or less flat to avoid deleting bonds
// in molecules that have a closely packed tetrahedral geometry
// (e. g. tetrahedrane).
if (StrictMath.abs(getDihedral(mol, ringAtoms[0], ringAtoms[1],
ringAtoms[2], ringAtoms[3])) > 20 * StrictMath.PI / 180) continue;
for (int ringAtom : ringAtoms) {
int inRingNeighbours = 0;
for (int k = 0; k < mol.getConnAtoms(ringAtom); k++) {
if (isRingAtom(mol, i, mol.getConnAtom(ringAtom, k))) {
inRingNeighbours++;
}
}
if (inRingNeighbours > 2) {
int longestBondIndex = -1;
double longestBondLength = -Double.MAX_VALUE;
for (int k = 0; k < mol.getConnAtoms(ringAtom); k++) {
if (!isRingAtom(mol, i, mol.getConnAtom(ringAtom, k))) continue;
if (longestBondLength < getBondLength(mol, mol.getConnBond(ringAtom, k))) {
longestBondIndex = k;
longestBondLength = getBondLength(mol, mol.getConnBond(ringAtom, k));
}
}
// if (isPlanar(mol, ringAtom, mol.getConnAtom(ringAtom, longestBondIndex))) {
mol.markBondForDeletion(mol.getConnBond(ringAtom, longestBondIndex));
// }
}
}
}
}
mol.deleteMarkedAtomsAndBonds();
// Processing non-organic atoms
grid = new MoleculeGrid(mol);
for( int i = 0; i < mol.getAllAtoms(); i++ ) {
if (mol.isOrganicAtom(i)) continue;
// TODO: since these are not covalent bonds, maybe none of the metals should be excluded?
// Not Lithium or Magnesium
// TODO: also Be in organometallic compounds?
// TODO: allow Alkali metals to coordinate water?
// if ( (isAlkaliMetalAtom(mol, i) && mol.getAtomicNo(i) != 3 ) ||
// (isAlkalineEarthMetalAtom(mol, i) && mol.getAtomicNo(i) != 12 ) ) continue;
Set set = grid.getNeighbours(mol, i, 3.2);
for (int j : set) {
double dist = getBondLength( mol, i, j );
double idealDist = COVALENT_RADIUS[mol.getAtomicNo(i)] + COVALENT_RADIUS[mol.getAtomicNo(j)];
if (dist > idealDist + .45) continue;
mol.addBond(i, j, Molecule3D.cBondTypeMetalLigand);
}
}
int ringNumber = mol.getRingSet().getSize();
Coordinates[] ringCenters = new Coordinates[ringNumber];
for (int i = 0; i < ringNumber; i++) {
ringCenters[i] = getAtomRingCenter(mol, i);
}
ArrayList[] atomsToRings = getAtomToRings(mol);
for ( int i = 0; i < mol.getAllAtoms(); i++ ) {
if ( !mol.isMetalAtom(i) ) continue;
boolean[] ringProcessed = new boolean[ringNumber];
boolean[] ringCoordinated = new boolean[ringNumber];
for ( int j = mol.getAllConnAtoms(i); j < mol.getAllConnAtomsPlusMetalBonds(i); j++ ) {
int ligand = mol.getConnAtom(i, j);
if ( mol.isMetalAtom(ligand) ) continue;
if ( atomsToRings[ligand] == null ) continue;
int closestRing = -1;
double closestDist = 0;
for ( int ringIndex: atomsToRings[ligand] ) {
if ( closestRing < 0 ||
ringCenters[ringIndex].distance(mol.getCoordinates(i)) < closestDist ) {
closestDist = ringCenters[ringIndex].distance(mol.getCoordinates(i));
closestRing = ringIndex;
}
}
if ( ringProcessed[closestRing] ) continue;
int bondCount = 0;
for ( int k = 0; k < mol.getAllConnAtomsPlusMetalBonds(i); k++ ) {
if ( isRingAtom(mol, closestRing, mol.getConnAtom(i, k) ) ) bondCount++;
}
// TODO: check is ring is (potentially) aromatic?
ringProcessed[closestRing] = true;
if ( bondCount > mol.getRingSet().getRingSize(closestRing) / 2 ) {
ringCoordinated[closestRing] = true;
}
}
for ( int j = 0; j < ringNumber; j++ ) {
if ( ringCoordinated[j] ) {
int[] ringAtoms = mol.getRingSet().getRingAtoms(j);
for (int ringAtom : ringAtoms) {
mol.addBond(i, ringAtom, Molecule.cBondTypeMetalLigand);
}
}
}
}
mol.deleteMarkedAtomsAndBonds();
mol.ensureHelperArrays(Molecule.cHelperNeighbours);
// TODO: ask about the applicability of this part
// remove erroneously assigned metal-ligand bonds
ArrayList bondsToDelete = new ArrayList<>();
for (int i = 0; i < mol.getAllBonds(); i++) {
if ( mol.getBondType(i) != Molecule3D.cBondTypeMetalLigand ) continue;
int metal;
int ligand;
if ( mol.isMetalAtom( mol.getBondAtom(0, i) ) ) {
metal = mol.getBondAtom(0, i);
ligand = mol.getBondAtom(1, i);
} else {
metal = mol.getBondAtom(1, i);
ligand = mol.getBondAtom(0, i);
}
for (int j = 0; j < mol.getAllConnAtomsPlusMetalBonds(ligand); j++ ) {
int ligandNeighbour = mol.getConnAtom(ligand, j);
if ( ligandNeighbour == metal ) continue;
double sigma = 0.2;
if ( getBondLength( mol, metal, ligand ) - getBondLength(mol, metal, ligandNeighbour) > sigma) {
bondsToDelete.add(i);
}
}
}
bondsToDelete.forEach(mol::markBondForDeletion);
mol.deleteMarkedAtomsAndBonds();
}
/**
* Evaluates if the given three-membered ring contains an unlikely bond
* between two phosphorus atoms. The bond is considered unlikely if it
* occurs in a three-membered ring and is above the specified bond length
* threshold.
*
* Such incorrect bonds are likely to be placed in molecules where two
* phosphorus (P) atoms are bridged by a third atom forming an acute angle.
* For example, COD entry 4071806 contains the PCP pattern that is incorrectly
* fused into a P1CP1 ring due to the distance between the two P atoms (2.426 A)
* falling into the range of previously observed long P-P bonds. It should be noted,
* however, that not all diphosporus bonds appearing in three-membered rings are
* incorrect as illustrated by COD entry 4062100. As such, additional bond length
* threshold is used to determine is the bond should be marked for deletion.
*
* @param mol
* Molecule that contains the ring.
* @param ringIndex
* Index of the ring in the molecule.
* @return
* Boolean value denoting if an improbable diphosphorus
* bond was identified and marked for deletion.
*/
private static boolean markImprobableDiphosphorusBondForDeletionInThreeAtomRing(Molecule3D mol, int ringIndex) {
if (mol.getRingSet().getRingSize(ringIndex) != 3) return false;
int[] ringBonds = mol.getRingSet().getRingBonds(ringIndex);
int longestDiphosphoroBondIndex = -1;
double longestDiphosphoroBondLength = 0;
for (int ringBond : ringBonds) {
int a1 = mol.getBondAtom(0, ringBond);
if (mol.getAtomicNo(a1) != 15) continue;
int a2 = mol.getBondAtom(1, ringBond);
if (mol.getAtomicNo(a2) != 15) continue;
if (getBondLength(mol, ringBond) > longestDiphosphoroBondLength) {
longestDiphosphoroBondIndex = ringBond;
longestDiphosphoroBondLength = getBondLength(mol, ringBond);
}
}
if (longestDiphosphoroBondIndex == -1) return false;
// TODO: eventually replace the hardcoded value with one based on statistics
// Explanation: the current cutoff values of 2.42 A was selected
// after analysing about a hundred structures that were assigned
// extremely long diphosphorus bonds. The longest correct
// disphophorus bond in a three-membered ring with the bond
// length of 2.409 was observed in COD entry 7107333.
if (longestDiphosphoroBondLength < 2.42) return false;
mol.markBondForDeletion(longestDiphosphoroBondIndex);
return true;
}
/**
* Evaluates if the given three-membered ring contains an unlikely bond
* that involves a hydrogen atom.
*
* Such incorrect bonds are placed in molecules where a hydrogen atom
* is positioned withing a bonding distance of two atoms. For example,
* COD entry 7027071 contains a hydrogen atom ('H4C') that is within
* a 0.8 A distance of the nitrogen atom ('N4') and a 1.46 A distance
* of a carbon atom ('C8'). As a result, an incorrect bond between
* the 'H4C' and 'C8' atoms gets placed forming a three-membered ring
* (C1HN1).
*
* @param mol
* Molecule that contains the ring.
* @param ringIndex
* Index of the ring in the molecule.
* @return
* Boolean value denoting if an improbable bond that involves
* a hydrogen atom was identified and marked for deletion.
*/
private static boolean markImprobableBondWithHydrogenForDeletionInThreeAtomRing(Molecule3D mol, int ringIndex){
if (mol.getRingSet().getRingSize(ringIndex) != 3) return false;
// TODO: allow two boron atoms to share a hydrogen
boolean bondMarkedForDeletion = false;
int[] ringAtoms = mol.getRingSet().getRingAtoms(ringIndex);
for (int j = 0; j < ringAtoms.length; j++) {
if (mol.getAtomicNo(ringAtoms[j]) != 1) continue;
int n1 = ringAtoms[((j - 1) + ringAtoms.length) % ringAtoms.length];
int n2 = ringAtoms[((j + 1) + ringAtoms.length) % ringAtoms.length];
int boundAtom = n1;
int unboundAtom = n2;
if (getBondLength(mol, ringAtoms[j], n1) > getBondLength(mol, ringAtoms[j], n2)) {
boundAtom = n2;
unboundAtom = n1;
}
mol.markBondForDeletion(mol.getBond(ringAtoms[j], unboundAtom));
System.err.printf(
"WARNING",
"a bond between atoms '" + mol.getAtomName(ringAtoms[j]) + "' and '" +
mol.getAtomName(unboundAtom) + "' was initially assigned based on the interatomic " +
"distance (" + "%.3f" + " nm), but was " +
"subsequently removed based on the molecular geometry -- atoms '" +
mol.getAtomName(ringAtoms[j]) + "', '" + mol.getAtomName(boundAtom) + "' and '" +
mol.getAtomName(unboundAtom) + "' are very unlikely to form a three-membered atom ring" + "\n",
getBondLength(mol, ringAtoms[j], unboundAtom)/10
);
bondMarkedForDeletion = true;
}
return bondMarkedForDeletion;
}
/**
* Calculates the geometrical center of the atom ring.
* @param mol
* The molecule that contains the ring.
* @param ringIndex
* The ring index in the molecule.
* @return
* Coordinates of the geometrical center of the atom ring.
*/
private static Coordinates getAtomRingCenter(Molecule3D mol, int ringIndex) {
Coordinates c = new Coordinates();
for ( int i: mol.getRingSet().getRingAtoms(ringIndex) ) {
c.x += mol.getAtomX(i);
c.y += mol.getAtomY(i);
c.z += mol.getAtomZ(i);
}
int ringAtomCount = mol.getRingSet().getRingAtoms(ringIndex).length;
c.x /= ringAtomCount;
c.y /= ringAtomCount;
c.z /= ringAtomCount;
return c;
}
/**
* Calculates for the organic subset of atoms the maximum number of neighbors
* @param mol
* The molecule that contains the atom.
* @param atom
* The index of the atom in the molecule.
* @return
* The maximum neighbour count for the atom.
*/
private static int maxNeighborCount(Molecule3D mol, int atom) {
int atomicNo = mol.getAtomicNo(atom);
return atomicNo == 5 ? 6 // B
: atomicNo <= 7 ? 4 // N, C
: atomicNo == 8 ? 3 // O
: atomicNo == 9 ? 1 // F
: atomicNo == 14 ? 6 // Si
: atomicNo == 15 ? 6 // P
: atomicNo == 16 ? 6 // S
: atomicNo == 17 ? 4 // Cl
: atomicNo == 33 ? 6 // As
: atomicNo == 34 ? 6 // Se
: atomicNo == 35 ? 6 // Br
: atomicNo == 52 ? 6
: atomicNo == 53 ? 6 : 8;
}
/**
* Calculate the bond orders of the molecule (without knowing the hydrogens).
* The calculation is based on the bond distance between each atoms.
*
* http://www.ccp14.ac.uk/ccp/web-mirrors/i_d_brown/valence.txt
* s = exp((Ro - R)/B)
*
* The implementation of this method is also most similar to the algorithm
* described in:
* http://www.daylight.com/meetings/mug01/Sayle/m4xbondage.html
*
* @param mol
* The molecule for which the bond orders should be calculated.
*/
public static void calculateBondOrders(Molecule3D mol) {
mol.ensureHelperArrays(Molecule.cHelperRings);
/*
* Pass 1: assign higher order bonds in unambiguous locations
*/
assignObviousHigherOrderBonds(mol);
// Place disilene bonds
assignDisileneBonds(mol);
//////////////////////////////////
// Functional Group Recognition
// Hybridization State Determination
resolveKnownIons(mol);
int[] spOrder = calculateHybridizationStates(mol);
recogniseFunctionalGroups(mol, spOrder);
// FIXME: this method should eventually be moved to the start of the method
resolveSharedHydrogenAtoms(mol);
// hybridization states need to be recalculated since
// the resolveSharedHydrogenAtoms() method might affect
// the order of atoms in the array
spOrder = calculateHybridizationStates(mol);
boolean[] visited = new boolean[mol.getAllBonds()];
if (DEBUG) {
System.err.println("> DEBUG: entering preliminary pass phase");
}
//Preliminary pass: process obvious bonds outside rings
for (int i = 0; i < mol.getBonds(); i++) {
if (mol.getBondType(i) == Molecule.cBondTypeMetalLigand) continue;
if (mol.isRingBond(i)) continue;
int a1 = mol.getBondAtom(0, i);
int a2 = mol.getBondAtom(1, i);
if (mol.getLowestFreeValence(a1) < 2) continue;
if (mol.getLowestFreeValence(a2) < 2) continue;
// if (!isPlanar(mol, a1, a2)) continue;
if (spOrder[a1] != 1 || spOrder[a2] != 1) continue;
int index = BondLengthSet.getBondIndex(3, false, false, mol.getAtomicNo(a1), mol.getAtomicNo(a2), 2, 2);
if (index == -1) continue;
double bondLength = BondLengthSet.getBondLength(index);
double bondLengthStd = BondLengthSet.getBondStdDev(index);
// Used to be 3 order > TRIPLE_BOND_CUTOFF
if (bondLength + bondLengthStd > getBondLength(mol, i)) {
mol.setBondOrder(i, 3);
visited[i] = true;
if (DEBUG) {
System.err.println(">> DEBUG: placing a triple bond between atoms '" +
mol.getAtomName(mol.getBondAtom(0, i)) + "' and '" +
mol.getAtomName(mol.getBondAtom(1, i)) + "' based on " +
"the geometry of the molecule and bond length statistics.");
}
}
}
if (DEBUG) {
System.err.println("> DEBUG: exiting preliminary pass phase");
}
identifySquarates(mol, spOrder);
if (DEBUG) {
System.err.println("> DEBUG: entering exocyclic double bond detection phase.");
}
// go through all exocyclic atoms and set the most likely double bonds
boolean[] isNonAromatic = assignExocyclicDoubleBonds(mol, spOrder);
// check exocyclic atoms that are part of non-aromatic rings
for (int i = 0; i < mol.getRingSet().getSize(); i++) {
if (isNonAromatic[i]) continue;
int[] ringAtoms = mol.getRingSet().getRingAtoms(i);
for (int ringAtom : ringAtoms) {
if (mol.getLowestFreeValence(ringAtom) < 1) continue;
for (int j = 0; j < mol.getConnAtoms(ringAtom); j++) {
// check if atom is exocyclic and capable of a double bond
int exoAtom = mol.getConnAtom(ringAtom, j);
if (mol.getAtomCharge(exoAtom) != 0) continue;
if (mol.getLowestFreeValence(exoAtom) < 1) continue;
if (spOrder[exoAtom] != 2) continue;
// if ( mol.getAtomicNo(exoAtom) == 8 ) continue;
// Check if the neighbour atom doesn't belong to the same ring
if (isRingAtom(mol, i, exoAtom)) continue;
boolean allowsDoubleBonds = true;
if (mol.isRingAtom(exoAtom)) {
for (int k = 0; k < mol.getRingSet().getSize(); k++) {
if (isRingAtom(mol, k, exoAtom) && !isNonAromatic[k]) {
allowsDoubleBonds = false;
break;
}
}
}
if (!allowsDoubleBonds) continue;
if (!isPlanar(mol, ringAtom, exoAtom)) continue;
int minPiCount = (mol.getConnAtoms(exoAtom) > 1) ? 1 : 0;
double averageSingleBondLength = getIdealisedBondLength(1, mol.getAtomicNo(ringAtom), mol.getAtomicNo(exoAtom), 1, minPiCount);
double averageDoubleBondLength = getIdealisedBondLength(2, mol.getAtomicNo(ringAtom), mol.getAtomicNo(exoAtom), 1, 1);
double averageDoubleBondStd = BondLengthSet.getBondStdDev(BondLengthSet.getBondIndex(2, false, false, mol.getAtomicNo(ringAtom), mol.getAtomicNo(exoAtom), 1, 1));
double calculatedLength = getBondLength(mol, ringAtom, exoAtom);
double doubleBondDiff = StrictMath.abs(averageDoubleBondLength - calculatedLength);
if (averageSingleBondLength < averageDoubleBondLength) {
if (DEBUG) {
System.err.println(">>> DEBUG: single '" + mol.getAtomicNo(ringAtom) + "-" +
mol.getAtomicNo(exoAtom) + "' bond seems to be shorter than a the " +
"double bond of the same type -- skipping the estimation of this bond.");
}
continue;
}
if (doubleBondDiff < 0 || doubleBondDiff <= averageDoubleBondStd * 1.5) {
mol.setBondOrder(mol.getBond(ringAtom, exoAtom), 2);
if (DEBUG) {
System.err.println(">>> DEBUG: placing a double bond between atoms '" +
mol.getAtomName(ringAtom) + "' and '" + mol.getAtomName(exoAtom) + "'.");
}
}
}
}
}
if (DEBUG) {
System.err.println("> DEBUG: exiting exocyclic ring atom double bond detection phase.");
}
/////////////////////////////////////////////////////////
// Aromatic Ring Perception
// This procedure calculates a normal to the ring and check that all
// atoms in the ring and their neighbours are within the plane defined by the normal
if ( DEBUG ) {
System.err.println("> DEBUG: entering aromatic ring perception phase");
}
RingCollection ringSet = mol.getRingSet();
boolean[] aromaticRing = findAromaticRings(mol, spOrder);
for (int i = 0; i < ringSet.getSize(); i++) {
if ( aromaticRing[i] ) {
boolean isResolved = false;
if (resolveTetraatomicAromaticRingIons(mol,i)) {
isResolved = true;
// Cyclothiaselenazenium
} else if ( ringSet.getRingSize(i) == 7 ) {
int NSBondCount = 0;
int SSBondCount = 0;
for (int bond : ringSet.getRingBonds(i)) {
if ( mol.getAtomicNo(mol.getBondAtom(0, bond) ) == 16 &&
mol.getAtomicNo(mol.getBondAtom(1, bond) ) == 16 ) {
SSBondCount++;
} else if (
( mol.getAtomicNo(mol.getBondAtom(0, bond)) == 7 &&
mol.getAtomicNo(mol.getBondAtom(1, bond)) == 16 ) ||
( mol.getAtomicNo(mol.getBondAtom(0, bond)) == 16 &&
mol.getAtomicNo(mol.getBondAtom(1, bond)) == 7 ) ) {
NSBondCount++;
} else {
break;
}
}
if ( NSBondCount == 6 && SSBondCount == 1 ) {
int SSBond = -1;
int ringSize = ringSet.getRingSize(i);
for (int j = 0; j < ringSet.getRingBonds(i).length; j++) {
int bond = ringSet.getRingBonds(i)[j];
if ( mol.getAtomicNo(mol.getBondAtom(0, bond) ) == 16 &&
mol.getAtomicNo(mol.getBondAtom(1, bond) ) == 16 ) {
SSBond = j;
break;
}
}
if (SSBond == -1) {
throw new RuntimeException("the S-S bond could not be found even though " +
"it was previously detected.");
} else {
int doubleBond1 = ringSet.getRingBonds(i)[(SSBond - 2 + ringSize) % ringSize];
int doubleBond2 = ringSet.getRingBonds(i)[(SSBond + 2 + ringSize) % ringSize];
int chargedNBond = ringSet.getRingBonds(i)[(SSBond + 3 + ringSize) % ringSize];
mol.setBondOrder(doubleBond1, 2);
int sulphurAtom = mol.getAtomicNo(mol.getBondAtom(0,doubleBond1)) == 16 ? 0 : 1;
mol.setAtomCharge(mol.getBondAtom(sulphurAtom, doubleBond1), 1);
mol.setBondOrder(doubleBond2, 2);
sulphurAtom = mol.getAtomicNo(mol.getBondAtom(0,doubleBond2)) == 16 ? 0 : 1;
mol.setAtomCharge(mol.getBondAtom(sulphurAtom, doubleBond2), 1);
int nitrogenAtom = mol.getAtomicNo(mol.getBondAtom(0,doubleBond1)) == 7 ? 0 : 1;
mol.setAtomCharge(mol.getBondAtom(nitrogenAtom, chargedNBond), 1);
}
isResolved = true;
}
}
if (!isResolved) {
for (int bond : ringSet.getRingBonds(i)) {
mol.setBondType(bond, Molecule3D.cBondTypeDelocalized);
}
}
}
}
//Aromatizer
AromaticityResolver resolver = new AromaticityResolver(mol);
resolver.locateDelocalizedDoubleBonds(null, true, true);
mol.ensureHelperArrays(Molecule3D.cHelperRings);
if (DEBUG) {
for (int i = 0; i < ringSet.getSize(); i++) {
if ( aromaticRing[i] ) {
System.err.println(">> DEBUG: entering ring number " + i );
for (int bond : ringSet.getRingBonds(i)) {
System.err.println(">>> DEBUG: aromatic bond between atoms " +
mol.getAtomName(mol.getBondAtom(0, bond) )+ "' and '" +
mol.getAtomName(mol.getBondAtom(1, bond) ) + "' was marked as a " +
"bond of order " + mol.getBondOrder(bond));
}
}
}
System.err.println("> DEBUG: exiting aromatic ring perception phase");
}
assignObviousHigherOrderBonds(mol);
// Find pi bonds
if ( DEBUG ) {
System.err.println("> DEBUG: entering pi bond detection phase");
}
boolean[] visitedAtoms = new boolean[mol.getAtoms()];
List piBonds = new LinkedList<>();
for (int i = 0; i < mol.getAtoms(); i++) {
// skip atoms that cannot accommodate an additional double bond
if (mol.getConnAtoms(i) < 2) continue;
if (spOrder[i] != 2) continue;
if (mol.getAtomPi(i) != 0) continue;
if (mol.isAromaticAtom(i)) continue;
if (mol.getAtomCharge(i) != 0) continue;
if (mol.getLowestFreeValence(i) < 1) continue;
visitedAtoms[i] = true;
for (int j = 0; j < mol.getConnAtoms(i); j++) {
int a1 = mol.getConnAtom(i, (j + 1) % mol.getConnAtoms(i));
int a2 = mol.getConnAtom(i, j);
if (visitedAtoms[a2]) continue;
// skip atoms that cannot accommodate an additional double bond
if (spOrder[a2] != 2) continue;
if (mol.getAtomPi(a2) != 0) continue;
if (mol.getAtomCharge(a2) != 0) continue;
if (mol.getLowestFreeValence(a2) < 1) continue;
if (!fitsAromaticLengthConstraints(mol, i, a2)) continue;
Coordinates origin = mol.getCoordinates(i);
Coordinates normal = calculatePlaneNormal(origin, mol.getCoordinates(a1), mol.getCoordinates(a2));
if (normal.distSq() == 0) continue;
boolean isFlat = true;
// TODO: join them in a list?
for (int k = 0; k < mol.getConnAtoms(i); k++) {
if (calculatePointToPlaneDistance(normal, origin, mol.getCoordinates(mol.getConnAtom(i, k))) > 0.4) {
isFlat = false;
break;
}
}
for (int k = 0; k < mol.getConnAtoms(a2); k++) {
if (calculatePointToPlaneDistance(normal, origin, mol.getCoordinates(mol.getConnAtom(a2, k))) > 0.4) {
isFlat = false;
break;
}
}
if (isFlat) {
piBonds.add(mol.getConnBond(i, j));
if (DEBUG) {
System.err.println(">> DEBUG: bond between atoms '" + mol.getAtomName(i) +
"' and '" + mol.getAtomName(mol.getConnAtom(i, j)) + "' has been marked " +
"as delocalized");
}
}
}
}
if ( DEBUG ) {
System.err.println("> DEBUG: exiting pi bond detection phase");
}
for (int piBond : piBonds) {
mol.setBondType(piBond, Molecule3D.cBondTypeDelocalized);
}
mol.ensureHelperArrays(Molecule3D.cHelperRings);
List> fragments = getPiBondFragments(mol);
for (List fragment : fragments) {
if (!(mol.getLowestFreeValence(fragment.get(0)) > 0 &&
mol.getLowestFreeValence(fragment.get(fragment.size() - 1)) > 0)) {
continue;
}
// if the fragment is capable of interchanging single-double bonds
// then the more likely pattern should be selected
double[] singleBonds = new double[fragment.size() - 1];
double[] doubleBonds = new double[fragment.size() - 1];
for (int j = 0; j < fragment.size() - 1; j++) {
singleBonds[j] = estimatePiBondFitness(mol, fragment.get(j), fragment.get(j + 1), 1);
doubleBonds[j] = estimatePiBondFitness(mol, fragment.get(j), fragment.get(j + 1), 2);
}
if (fragmentCanResonate(mol, fragment)) {
double leftToRight = 0;
double rightToLeft = 0;
for (int j = 0; j < fragment.size() - 1; j++) {
// start with a single bond
leftToRight += (j % 2) == 0 ? singleBonds[j] : doubleBonds[j];
// start with a double bond
rightToLeft += (j % 2) == 0 ? doubleBonds[j] : singleBonds[j];
}
if (fragment.size() % 2 == 1) {
for (int j = 0; j < fragment.size() - 1; j++) {
int bondOrder;
if (leftToRight < rightToLeft) {
bondOrder = (j % 2 == 0 ? 1 : 2);
} else {
bondOrder = (j % 2 == 0 ? 2 : 1);
}
mol.setBondOrder(mol.getBond(fragment.get(j), fragment.get(j + 1)), bondOrder);
if (DEBUG) {
System.err.println(">>> DEBUG: a bond of order " + bondOrder + " was placed between " +
"atoms '" + mol.getAtomName(fragment.get(j)) + "' and '" +
mol.getAtomName(fragment.get(j + 1)) + "' outside of the aromatizer");
}
}
}
} else {
for (int j = 0; j < fragment.size() - 1; j++) {
if (doubleBonds[j] < 0.03) {
mol.setBondOrder(mol.getBond(fragment.get(j), fragment.get(j + 1)), 2);
}
}
}
}
mol.ensureHelperArrays(Molecule3D.cHelperRings);
resolver = new AromaticityResolver(mol);
resolver.locateDelocalizedDoubleBonds(null, true, true);
mol.ensureHelperArrays(Molecule3D.cHelperRings);
if ( DEBUG ) {
for (int piBond : piBonds) {
int a1 = mol.getBondAtom(0, piBond);
int a2 = mol.getBondAtom(1, piBond);
System.err.println(">>>> DEBUG: a bond of order " + mol.getBondOrder(piBond) + " was " +
"placed between atoms '" + mol.getAtomName(a1) + "' and '" + mol.getAtomName(a2) +
"'.");
}
}
// Try to recognise the N-C-N pattern and set the charges accordingly
int index = BondLengthSet.getBondIndex(1, true, false, 6, 7, 1 , 1);
double bondLength = BondLengthSet.getBondLength(index);
double bondStd = BondLengthSet.getBondStdDev(index);
for ( int i = 0; i < mol.getAtoms(); i++ ) {
if ( mol.getAtomicNo(i) == 6 && spOrder[i] == 2 && mol.getConnAtoms(i) > 1 ) {
for (int j = 0; j < mol.getConnAtoms(i); j++) {
int a1 = mol.getConnAtom(i, j);
int a2 = mol.getConnAtom(i, (j + 1) % mol.getConnAtoms(i));
if (spOrder[a1] != 2 || spOrder[a2] != 2) continue;
if (mol.getAtomicNo(a1) != 7 || mol.getAtomicNo(a2) != 7) continue;
double a1Distance = getBondLength(mol, i, a1);
double a2Distance = getBondLength(mol, i, a2);
if ( ( a1Distance > bondLength + 3 * bondStd ) ||
mol.getBondOrder(mol.getBond(i, a1)) == 2 ) continue;
if ( ( a2Distance > bondLength + 3 * bondStd ) ||
mol.getBondOrder(mol.getBond(i, a1)) == 2 ) continue;
int shorterBondNitrogen = a1Distance < a2Distance ? a1 : a2;
// TODO: add an angle constraint
if ( mol.getAtomPi(i) == 0 ) {
if ( mol.getAtomCharge(a1) != 0 || mol.getAtomCharge(a2) != 0 ) continue;
if ( mol.getAtomPi(a1) != 0 || mol.getAtomPi(a2) != 0 ) continue;
if ( mol.getBondOrder(mol.getBond(i, a1)) < 2 &&
mol.getBondOrder(mol.getBond(i, a2)) < 2) {
mol.setBondOrder(mol.getBond(i, shorterBondNitrogen), 2);
mol.setAtomCharge(shorterBondNitrogen, +1);
mol.ensureHelperArrays(Molecule.cHelperRings);
}
} else {
// Switch N-C=N to N=C-N in if the latter resonance form seems more likely
// exocyclic aromatic resonance should not be considered
if ( mol.isAromaticAtom(a1) != mol.isAromaticAtom(a2) ) continue;
int longerBondNitrogen = ( shorterBondNitrogen == a1 ) ? a2 : a1;
if ( mol.getBondOrder(mol.getBond(i, longerBondNitrogen)) == 2 &&
mol.getAtomCharge(longerBondNitrogen) == 1 &&
mol.getAtomPi(shorterBondNitrogen) == 0 ) {
mol.setBondOrder(mol.getBond(i, longerBondNitrogen), 1);
mol.setBondOrder(mol.getBond(i, shorterBondNitrogen), 2);
mol.setAtomCharge(longerBondNitrogen, 0);
mol.setAtomCharge(longerBondNitrogen, +1);
mol.ensureHelperArrays(Molecule.cHelperRings);
}
}
}
}
}
mol.ensureHelperArrays(Molecule.cHelperRings);
/*
* 2nd pass: find obvious double bonds on sp and sp2 B, C and N outside aromatic rings
*/
assignObviousDoubleBonds(mol, spOrder);
/*
* 3rd pass: double and triple bonds that correspond to the bond length
* (but might deviate from the expected geometry)
*/
IntQueue queue = new IntQueue();
for (int i = 0; i < mol.getAllBonds(); i++) {
if ( mol.getBondType(i) == Molecule3D.cBondTypeMetalLigand ) continue;
if (!visited[i] && !mol.isAromaticBond(i)) queue.push(i);
while (!queue.isEmpty()) {
int bnd = queue.pop();
if (visited[bnd]) continue;
visited[bnd] = true;
int a1 = mol.getBondAtom(0, bnd);
int a2 = mol.getBondAtom(1, bnd);
if ( a1 >= spOrder.length || a2 >= spOrder.length) continue;
// Push the neighbour bonds into the queue, but
// avoid processing hydrogen and metal atoms
for (int j = 0; j < mol.getConnAtoms(a1); j++) {
queue.push(mol.getConnBond(a1, j));
}
for (int j = 0; j < mol.getConnAtoms(a2); j++) {
queue.push(mol.getConnBond(a2, j));
}
//Compute the free valence and increase the bond order if needed
double realBondLength = getBondLength(mol, a1, a2);
int doubleBondIndex = BondLengthSet.getBondIndex(2, false, false, mol.getAtomicNo(a1), mol.getAtomicNo(a2), 1, 1);
if ( doubleBondIndex == -1 ) continue;
double idealisedDoubleBondLength = BondLengthSet.getBondLength(doubleBondIndex);
double idealisedDoubleBondStdDev = BondLengthSet.getBondStdDev(doubleBondIndex);
if ( idealisedDoubleBondStdDev == BondLengthSet.DEFAULT_BOND_STDDEV ) {
idealisedDoubleBondStdDev = 0.1;
}
if ( realBondLength < idealisedDoubleBondLength + 3 * idealisedDoubleBondStdDev &&
!mol.isAromaticAtom(a1) && !mol.isAromaticAtom(a2) ) {
//Special case CS
if (mol.getAtomicNo(a1) == 16 && mol.getAllConnAtoms(a1) <= 2) continue;
if (mol.getAtomicNo(a2) == 16 && mol.getAllConnAtoms(a2) <= 2) continue;
if ( mol.getAtomCharge(a1) != 0 || mol.getAtomCharge(a2) != 0 ) continue;
int freeValence1 = StrictMath.max( mol.getLowestFreeValence(a1),
StrictMath.max( ( mol.getAtomicNo(a1) == 7 && spOrder[a1] < 3 ) ? 1 : 0,
( mol.getAtomicNo(a1) == 5 && spOrder[a1] < 3 ) ? 1 : 0 ) );
int freeValence2 = StrictMath.max( mol.getLowestFreeValence(a2),
StrictMath.max( ( mol.getAtomicNo(a2) == 7 && spOrder[a2] < 3 ) ? 1 : 0,
( mol.getAtomicNo(a2) == 5 && spOrder[a2] < 2 ) ? 1 : 0 ) );
if ( freeValence1 < 1 || freeValence2 < 1 ) continue;
boolean aligned = ( spOrder[a1] == 1 && spOrder[a2] == 1 );
boolean fitsTripleBondLength;
int tripleBondIndex = BondLengthSet.getBondIndex(3, false, false, mol.getAtomicNo(a1), mol.getAtomicNo(a2), 2, 2);
if ( tripleBondIndex != -1 ) {
double idealisedTripleBondLength = BondLengthSet.getBondLength(tripleBondIndex);
double idealisedTripleBondStdDev = BondLengthSet.getBondStdDev(tripleBondIndex);
if (idealisedTripleBondStdDev == BondLengthSet.DEFAULT_BOND_STDDEV) {
idealisedTripleBondStdDev = 0.1;
}
fitsTripleBondLength = realBondLength > idealisedTripleBondLength + 3 * idealisedTripleBondStdDev;
} else {
fitsTripleBondLength = realBondLength < idealisedDoubleBondLength;
}
if ( aligned && mol.getAtomPi(a1) < 2 && mol.getAtomPi(a2) < 2 && fitsTripleBondLength ) {
mol.setBondOrder(bnd, 3);
} else if ( mol.getAtomPi(a1) < ( spOrder[a1] == 1 ? 2 : 1 ) &&
mol.getAtomPi(a2) < ( spOrder[a2] == 1 ? 2 : 1 ) ) {
mol.setBondOrder(bnd, 2);
}
}
}
}
mol.ensureHelperArrays(Molecule3D.cHelperRings);
/*
* 4th pass: increase bond order if the geometry requires it
* (even though the bond length might seem too great)
*/
int[] localAtomPi = getAtomPiCount(mol);
for (int i = 0; i < mol.getAtoms(); i++) {
if (!canAccommodateDoubleBond(mol, i, localAtomPi[i], spOrder[i])) continue;
// sp3 atoms smaller than Si should not have double bonds
if (spOrder[i] == 3 && mol.getAtomicNo(i) < 14) continue;
for (int j = 0; j < mol.getConnAtoms(i); j++) {
int a1 = mol.getConnAtom(i, j);
if (canAccommodateDoubleBond(mol, a1, localAtomPi[a1], spOrder[a1])) {
mol.setBondOrder(mol.getBond(i, a1), 2);
localAtomPi[i] += 1;
localAtomPi[a1] += 1;
}
if (localAtomPi[i] >= (spOrder[i] == 1 ? 2 : 1)) break;
}
}
// 5th pass: redistributing charges
/*
* v0 is the start atom
* v1 is the 1st order neighbour of atom v0
* v2 is the 2nd order neighbour of atom v0
* v3 is the 3rd order neighbour of atom v0
*/
localAtomPi = getAtomPiCount(mol);
for (int v0 = 0; v0 < mol.getAllAtoms(); v0++) {
if ( mol.getAtomicNo(v0) == 6 ) {
// C[-]-
if ( spOrder[v0] != 2 ) continue;
if ( mol.getLowestFreeValence(v0) != 1 ) continue;
if ( localAtomPi[v0] != 0 ) continue;
if ( hasMetalLigandBonds(mol, v0) ) continue;
// O=C-C(-)-C=O -> O=C-C=C-O(-)
List terminalOxygens = new LinkedList<>();
// C(-)-N -> C->N(+)
List sp2NitrogenNeighbours = new LinkedList<>();
for (int j = 0; j < mol.getConnAtoms(v0); j++) {
int v1 = mol.getConnAtom(v0, j);
if ( mol.getAtomicNo(v1) == 6 ) {
// C[-]-C(R)=
if ( spOrder[v1] != 2 ) continue;
if ( localAtomPi[v1] == 0 ) continue;
for (int k = 0; k < mol.getConnAtoms(v1); k++) {
int v2 = mol.getConnAtom(v1, k);
// select the second order neighbour atom that is double bonded
// to the first order neighbour atom
if ( v2 == v0 ) continue;
if ( mol.getBondOrder(mol.getConnBond(v1, k)) != 2 ) continue;
if ( mol.getAtomicNo(v2) == 8 &&
mol.getConnAtoms(v2) == 1 ) {
terminalOxygens.add(v2);
} else if ( mol.getAtomicNo(v2) == 6 ) {
// C[-]-C(R)=C(R)(R)
for ( int l = 0; l < mol.getConnAtoms(v2); l++ ) {
if ( localAtomPi[v0] > 0 ) break;
int v3 = mol.getConnAtom(v2, l);
if ( v1 == v3 ) continue;
if ( mol.getAtomicNo(v3) != 6 ) continue;
if ( spOrder[v3] != 2 ) continue;
if ( mol.getLowestFreeValence(v3) != 1 ) continue;
if ( hasMetalLigandBonds(mol, v3) ) continue;
mol.setBondOrder(mol.getBond(v0, v1), 2);
mol.setBondOrder(mol.getBond(v1, v2), 1);
mol.setBondOrder(mol.getBond(v2, v3), 2);
localAtomPi[v0]++;
localAtomPi[v3]++;
break;
}
}
}
} else if (mol.getAtomicNo(v1) == 7 && spOrder[v1] == 2 &&
mol.getAtomPi(v1) == 0 && mol.getImplicitHydrogens(v1) == 0) {
sp2NitrogenNeighbours.add(j);
}
}
// FIXME: the valency of associated C atoms might change due to previous loop
// This needs serious refactoring
// moving charges from carbon atoms to oxygen atoms where possible
if ( mol.getLowestFreeValence(v0) > 0 ) {
if (terminalOxygens.size() > 0 ) {
mol.setBondOrder(mol.getConnBond(terminalOxygens.get(0), 0), 1);
mol.setBondOrder(mol.getBond(mol.getConnAtom(terminalOxygens.get(0), 0), v0), 2);
} else if (sp2NitrogenNeighbours.size() > 0) {
int[] sortedAtom = sortAtomsByBondLength(mol, v0, sp2NitrogenNeighbours);
mol.setBondOrder(mol.getConnBond(v0, sortedAtom[0]), 2);
// nitrogen atoms may contain a negative charge that was assigned during previous processing
int nAtom = mol.getConnAtom(v0, sortedAtom[0]);
if (mol.getAtomCharge(nAtom) < 0) {
mol.setAtomCharge(nAtom, mol.getAtomCharge(nAtom) + 1);
}
}
}
} else if ( mol.getAtomicNo(v0) == 16 ) {
// Sp2 sulphur with +1 charge with no metal bonds
if ( spOrder[v0] != 2 ) continue;
if ( mol.getAtomPi(v0) == 0 ) continue;
if ( mol.getAtomCharge(v0) != 1 ) continue;
if ( hasMetalLigandBonds(mol, v0) ) continue;
// (S+)=C-C=(S+) -> S-C=C-S
// (S+)=C-C -> (S+)-C=C
for (int j = 0; j < mol.getConnAtoms(v0); j++) {
// Sp2 carbon that is double bonded to the v0 sulphur
int v1 = mol.getConnAtom(v0, j);
if ( mol.getAtomicNo(v1) != 6 ) continue;
if ( mol.getBondOrder( mol.getBond(v0, v1) ) != 2 ) continue;
for (int k = 0; k < mol.getConnAtoms(v1); k++) {
// Sp2 carbon that is single bonded to the v1 carbon
int v2 = mol.getConnAtom(v1, k);
if ( mol.getAtomicNo(v2) != 6 ) continue;
if ( mol.getBondOrder(mol.getBond(v1,v2)) != 1 ) continue;
if ( mol.getLowestFreeValence(v2) == 1 ) {
mol.setBondOrder(mol.getBond(v1, v2), 2);
mol.setBondOrder(mol.getBond(v0, v1), 1);
mol.setAtomCharge(v0, 0);
} else if (mol.getAtomPi(v2) != 0){
for (int l = 0; l < mol.getConnAtoms(v2); l++) {
int v3 = mol.getConnAtom(v2, l);
if ( mol.getAtomicNo(v3) != 16 ) continue;
if ( mol.getAtomCharge(v3) != 1 ) continue;
if ( mol.getBondOrder(mol.getBond(v2, v3)) != 2 ) continue;
mol.setBondOrder(mol.getBond(v0, v1), 1);
mol.setBondOrder(mol.getBond(v1, v2), 2);
mol.setBondOrder(mol.getBond(v2, v3), 1);
mol.setAtomCharge(v0, 0);
mol.setAtomCharge(v3, 0);
}
}
}
}
}
}
mol.ensureHelperArrays(Molecule.cHelperRings);
// Set charges for atoms other than metals
assignAtomChargesBasedOnValence(mol, spOrder);
/*
* Transform the molecule to acquire a more reasonable resonance structure
*/
// R1(+)=X-R2(-) -> R1-X=R2
for (int i = 0; i < mol.getAtoms(); i++) {
if ( mol.getAtomCharge(i) != 1 ) continue;
if ( mol.getAtomPi(i) == 0 ) continue;
for ( int j = 0; j < mol.getConnAtoms(i); j++ ) {
int a1 = mol.getConnAtom(i, j);
if ( mol.getAtomCharge(a1) == 0 && mol.getConnBondOrder(i, j) == 2 ) {
for (int k = 0; k < mol.getConnAtoms(a1); k++) {
int a2 = mol.getConnAtom(a1, k);
double R1XBondLengthDiff = BondLengthSet.getBondLength(
BondLengthSet.getBondIndex(2, false, false,
mol.getAtomicNo(i),
mol.getAtomicNo(a1), 1, 1) ) - getBondLength(mol, i, a1);
if ( mol.getAtomCharge(i) < 0 && spOrder[a2] == 2 &&
mol.getAtomPi(a2) == 0 &&
( mol.getAtomCharge(i) + mol.getAtomCharge(a2) ) == 0 ) {
double R2XBondLengthDiff = BondLengthSet.getBondLength(
BondLengthSet.getBondIndex(2, false, false,
mol.getAtomicNo(a1),
mol.getAtomicNo(a2), 1, 1)) - getBondLength(mol, a1, a2);
if ( R2XBondLengthDiff < R1XBondLengthDiff ) {
mol.setBondOrder(mol.getBond(a1, a2), 2);
mol.setBondOrder(mol.getBond(i, a1), 1);
mol.setAtomCharge(i, 0);
mol.setAtomCharge(a2, 0);
break;
}
}
}
}
}
}
processKnownIonsInvolvingMetals(mol);
assignMetalCharges(mol);
mol.ensureHelperArrays(Molecule.cHelperNeighbours);
}
/**
* Gets the count of used pi electrons for each atom in the molecule.
* @param mol
* The molecule that contains the bonded atoms.
* @return
* Counts of used pi electron.
*/
private static int[] getAtomPiCount(Molecule3D mol) {
int[] atomPiCount = new int[mol.getAllAtoms()];
for ( int i = 0; i < atomPiCount.length; i++ ) {
if (mol.getAtomicNo(i) == 1) {
atomPiCount[i] = 0;
} else {
atomPiCount[i] = mol.getAtomPi(i);
}
}
return atomPiCount;
}
/**
* Assigns obvious higher order bonds where the geometry explicitly requires it.
*
* @param mol
* Molecule that contains the bonded atoms.
*/
private static void assignObviousHigherOrderBonds(Molecule3D mol) {
// Place obvious triple bonds for sp atoms
boolean changeOccurred = true;
int tripleBondIterationNo = 1;
while (changeOccurred) {
if (DEBUG) {
System.err.println(">> DEBUG: entered obvious triple bond assignment iteration " +
tripleBondIterationNo + ".");
}
changeOccurred = assignObviousTripleBonds(mol);
if (DEBUG) {
System.err.println(">> DEBUG: exited obvious triple bond assignment iteration " +
tripleBondIterationNo + ".");
}
tripleBondIterationNo++;
}
}
/**
* Assigns obvious triple bonds where the geometry explicitly requires it.
*
* @param mol
* Molecule that contains the bonded atoms.
* @return
* Boolean value denoting if any triple bonds were assigned.
*/
private static boolean assignObviousTripleBonds(Molecule3D mol) {
// the local pi count array is used to avoid calling ensureHelperArrays()
// each time a triple bond is assigned
boolean changeOccured = false;
int[] localAtomPi = getAtomPiCount(mol);
for (int i = 0; i < mol.getAtoms(); i++) {
if (mol.isMetalAtom(i)) continue;
if (localAtomPi[i] != 0) continue;
if (calculateHybridizationState(mol, i) != 1) continue;
if (mol.getAllConnAtoms(i) != 2) continue; // skip terminal atoms
boolean allowNeighbourCharge = true;
if (!canAccommodateTripleBond(mol, i, false)) {
if (!canAccommodateTripleBond(mol, i, true)) continue;
allowNeighbourCharge = false;
}
int a1 = mol.getConnAtom(i, 0);
int a2 = mol.getConnAtom(i, 1);
// skip if both atoms are capable of at least a single double bond
boolean a1CanCarryDoubleBond = canAccommodateDoubleBond(mol, a1, localAtomPi[a1], calculateHybridizationState(mol, a1));
boolean a2CanCarryDoubleBond = canAccommodateDoubleBond(mol, a2, localAtomPi[a2], calculateHybridizationState(mol, a2));
boolean a1CanCarryTripleBond = canAccommodateTripleBond(mol, a1, allowNeighbourCharge);
boolean a2CanCarryTripleBond = canAccommodateTripleBond(mol, a2, allowNeighbourCharge);
if ((a1CanCarryDoubleBond || a1CanCarryTripleBond) &&
(a2CanCarryDoubleBond || a2CanCarryTripleBond)) {
continue;
}
if (!a1CanCarryTripleBond && !a2CanCarryTripleBond) continue;
int tripleBondAtom = a1CanCarryTripleBond ? a1 : a2;
mol.setBondOrder(mol.getBond(i, tripleBondAtom), 3);
localAtomPi[i] = 2;
localAtomPi[tripleBondAtom] = 2;
changeOccured = true;
if (DEBUG) {
System.err.println(">>> DEBUG: placed a triple bond between atoms '" +
mol.getAtomName(i) + "' and '" + mol.getAtomName(tripleBondAtom) +
"' without relying on the bond length statistics.");
}
}
mol.ensureHelperArrays(Molecule.cHelperRings);
return changeOccured;
}
/**
* Evaluates is the given atom is capable of carrying a triple bond.
*
* @param mol
* Molecule that contains the atom.
* @param atom
* Index of the atom that should be evaluated.
* @param allowCharge
* Boolean value denoting if the atom is allowed to gain
* a formal charge in order to accommodate the triple bond
* (i.e. N+, B-).
* @return
* Boolean value denoting if the atom is capable of carrying a triple bond.
*/
private static boolean canAccommodateTripleBond(Molecule3D mol, int atom, boolean allowCharge) {
if (mol.getAtomicNo(atom) == 1) return false;
if (mol.isHalogene(atom)) return false;
int requiredValence = 2;
if (allowCharge) {
// boron (B) with the charge of -1 may carry a triple bond
if (mol.getAtomicNo(atom) == 5 ) requiredValence = 1;
// nitrogen (N) with the charge of +1 may carry a triple bond
if (mol.getAtomicNo(atom) == 7 ) requiredValence = 1;
}
if (mol.getLowestFreeValence(atom) < requiredValence) return false;
int spOrder = calculateHybridizationState(mol, atom);
return spOrder == 1;
}
/**
* Assigns obvious double bonds to boron, nitrogen and carbon atoms.
*
* @param mol
* The molecule that contains the bonded atoms.
* @param spOrder
* The sp hybridization type for each atom.
*/
private static void assignObviousDoubleBonds(Molecule3D mol, int[] spOrder) {
// using the local pi count array to avoid calling ensureHelperArrays
// each time a double bond is assigned
int[] atomPi = getAtomPiCount(mol);
for (int i = 0; i < mol.getAtoms(); i++) {
if ( mol.isAromaticAtom(i) ) continue;
if ( mol.getAtomCharge(i) != 0 ) continue;
// B, N, C with valence lower than 4
if ( mol.getAtomicNo(i) > 3 && mol.getAtomicNo(i) < 8 && mol.getOccupiedValence(i) < 4 ) {
if ( spOrder[i] == 1 ) {
for (int j = 0; j < mol.getConnAtoms(i); j++) {
int a1 = mol.getConnAtom(i, j);
if ( mol.getLowestFreeValence(a1) > 0 &&
mol.getAtomCharge(a1) == 0 &&
spOrder[a1] < 3 && atomPi[a1] <
( spOrder[a1] == 1 ? 2 :
spOrder[a1] == 2 ? 1 : 0 ) ) {
mol.setBondOrder(mol.getConnBond(i, j), 2);
atomPi[i]++;
atomPi[a1]++;
}
}
} else if (spOrder[i] == 2) {
if (atomPi[i] != 0) continue;
List potentialDoubleBondNeighbours = new LinkedList<>();
// nitrogen atoms with filled valencies are handled a little differently
boolean isFilledNitrogenAtom = mol.getAtomicNo(i) == 7 && mol.getLowestFreeValence(i) == 0;
for (int j = 0; j < mol.getConnAtoms(i); j++) {
int neighbour = mol.getConnAtom(i, j);
if ( mol.getAtomCharge(neighbour) != 0 ) continue;
int freeValence = StrictMath.max( mol.getLowestFreeValence(neighbour),
mol.getAtomicNo(neighbour) == 5 && spOrder[neighbour] < 3 ? 1 : 0 );
if ( freeValence < 1 ) continue;
if ( isFilledNitrogenAtom && mol.getAtomicNo(neighbour) == 7 && spOrder[neighbour] == 2 ) continue;
if ( atomPi[neighbour] < ( spOrder[neighbour] == 1 ? 2 :
spOrder[neighbour] == 2 ? 1 : 0 ) ||
mol.getAtomicNo(neighbour) == 15 ) {
potentialDoubleBondNeighbours.add(j);
}
}
if ( potentialDoubleBondNeighbours.size() > 1 ) {
double minDiff = -Double.MAX_VALUE;
int minDiffNeighbour = -1;
for (Integer neighbourIndex : potentialDoubleBondNeighbours) {
int neighbour = mol.getConnAtom(i, neighbourIndex);
double singleBondLength =BondLengthSet.getBondLength(
BondLengthSet.getBondIndex(1, false, false, mol.getAtomicNo(i), mol.getAtomicNo(neighbour), 1, 1));
double doubleBondLength = BondLengthSet.getBondLength(
BondLengthSet.getBondIndex(2, false, false, mol.getAtomicNo(i), mol.getAtomicNo(neighbour), 1, 1));
// TODO: think of a more legitimate way of detecting statistical anomalies
// safeguard against statistical anomalies
if ( singleBondLength < doubleBondLength || singleBondLength - doubleBondLength < 0.1 ) {
continue;
}
double bondLengthDifference = doubleBondLength - getBondLength(mol, i, neighbour);
if (bondLengthDifference > minDiff) {
minDiff = bondLengthDifference;
minDiffNeighbour = neighbour;
}
}
if (minDiffNeighbour > -1) {
mol.setBondOrder(mol.getBond(i, minDiffNeighbour), 2);
atomPi[i]++;
atomPi[minDiffNeighbour]++;
}
}
}
}
}
mol.ensureHelperArrays(Molecule.cHelperRings);
}
/**
* Assigns obvious disilene bonds.
*
* @param mol
* The molecule that contains the bonded atoms.
*/
private static void assignDisileneBonds(Molecule3D mol) {
double doubleBondCutoff = 2.29;
// FIXME: properly check if a normal index is returned
int bondIndex = BondLengthSet.getBondIndex(2, false, false, 14, 14, 1, 1 );
double siliconDoubleBondLength = BondLengthSet.getBondLength( bondIndex );
double siliconDoubleBondStdDev = BondLengthSet.getBondStdDev( bondIndex );
if ( siliconDoubleBondLength != -1 && siliconDoubleBondStdDev != -1 ) {
if ( ( siliconDoubleBondLength + 3 * siliconDoubleBondStdDev ) < doubleBondCutoff ) {
doubleBondCutoff = siliconDoubleBondLength + 3 * siliconDoubleBondStdDev;
}
}
int[] atomPi = getAtomPiCount(mol);
for (int i = 0; i < mol.getAtoms(); i++) {
if ( mol.getAtomicNo(i) != 14 ) continue;
if ( atomPi[i] != 0 ) continue;
if ( mol.getOccupiedValence(i) > 3 ) continue;
if ( mol.isAromaticAtom(i) ) continue;
for (int j = 0; j < mol.getConnAtoms(i); j++) {
int a1 = mol.getConnAtom(i, j);
// ring bonds are handled in aromaticity perception stage
if ( mol.isRingBond( mol.getBond( i, a1 ) ) ) continue;
if ( mol.getAtomicNo(a1) != 14 ) continue;
if ( mol.getOccupiedValence(a1) > 3 ) continue;
if ( atomPi[a1] != 0 ) continue;
if ( mol.isAromaticAtom(a1) ) continue;
if ( getBondLength( mol, i, a1 ) < doubleBondCutoff ) {
mol.setBondOrder( mol.getBond(i, a1), 2 );
atomPi[i]++;
atomPi[a1]++;
break;
}
}
}
}
/**
* Detects and assigns exocyclic double bonds that are more likely than
* double bonds inside the rings.
* @param mol
* Molecule that contains the atoms.
* @param spOrder
* The sp hybridization type for each atom.
* @return
* Boolean values denoting is an exocyclic bond has been placed
* for the corresponding atom ring.
*/
private static boolean[] assignExocyclicDoubleBonds(Molecule3D mol, int[] spOrder) {
boolean[] hasExocyclicDoubleBond = new boolean[mol.getRingSet().getSize()];
int[] localAtomPi = getAtomPiCount(mol);
for (int i = 0; i < mol.getRingSet().getSize(); i++) {
int[] ringAtoms = mol.getRingSet().getRingAtoms(i);
for (int ringAtom : ringAtoms) {
// skip the atom if the only exocyclic neighbours are hydrogen atoms
if (mol.getConnAtoms(ringAtom) == 2) continue;
// additional exclusion rules for certain elements
if (mol.getAtomicNo(ringAtom) == 14 || // Silicon (Si)
mol.getAtomicNo(ringAtom) == 15 || // Phosphorus (P)
mol.getAtomicNo(ringAtom) == 32) { // Germanium (Ge)
// skip atoms with full valency
if (mol.getFreeValence(ringAtom) < 1) continue;
// skip atoms that already contain a double bond
if (localAtomPi[ringAtom] > 0) continue;
}
// sp3d2 P and Te
if (spOrder[ringAtom] == 5 &&
(mol.getAtomicNo(ringAtom) == 15 || mol.getAtomicNo(ringAtom) == 52)) continue;
// skip the atom if it is not capable of a double bond
if (mol.getLowestFreeValence(ringAtom) <= 0 && mol.getAtomicNo(ringAtom) < 14) {
continue;
}
// find the best fit for an aromatic double bond inside the ring
double inRingMinDeviation = Double.MAX_VALUE;
for (int j = 0; j < mol.getConnAtoms(ringAtom); j++) {
int ringAtomNeighbour = mol.getConnAtom(ringAtom, j);
if (isRingAtom(mol, i, ringAtomNeighbour)) {
int bondTypeIndex = BondLengthSet.getBondIndex(2, true, false, mol.getAtomicNo(ringAtom),
mol.getAtomicNo(ringAtomNeighbour), 1, 1);
/*
* FIXME: the logic used for exocyclic bonds involving silicon (Si) should be generalised.
* In the future the logic should be generalised for all elements to initially try
* to retrieve the aromatic double bond distance and if it is not available fall back
* to the simple double bond one.
*/
if (mol.getAtomicNo(ringAtom) == 14 || mol.getAtomicNo(ringAtomNeighbour) == 14) {
if (bondTypeIndex == -1) {
bondTypeIndex = BondLengthSet.getBondIndex(2, false, false, mol.getAtomicNo(ringAtom),
mol.getAtomicNo(ringAtomNeighbour), 1, 1);
}
}
double currentDeviation = determineExpectedBondDeviation(mol,
mol.getBond(ringAtom, ringAtomNeighbour),
bondTypeIndex);
if (currentDeviation < inRingMinDeviation) {
inRingMinDeviation = currentDeviation;
}
}
}
for (int k = 0; k < mol.getConnAtoms(ringAtom); k++) {
// check if the ring atom is still capable of
// carrying an additional double bond
if (mol.getAtomicNo(ringAtom) == 15) { // phosphorus (P)
if (localAtomPi[ringAtom] > 0) break;
}
// check if atom is exocyclic and capable of a double bond
int exoAtom = mol.getConnAtom(ringAtom, k);
if (mol.isRingBond(mol.getBond(ringAtom, exoAtom))) continue;
if (mol.getAtomCharge(exoAtom) != 0) continue;
if (mol.getLowestFreeValence(exoAtom) < 1) continue;
// additional exclusion rules for certain elements
if (mol.getAtomicNo(exoAtom) == 14 || // Silicon (Si)
mol.getAtomicNo(exoAtom) == 15 || // Phosphorus (P)
mol.getAtomicNo(exoAtom) == 32) { // Germanium (Ge)
// skip atoms with full valency
if (mol.getFreeValence(exoAtom) < 1) continue;
// skip atoms that already contains a double bond
if (localAtomPi[exoAtom] > 0) continue;
}
// exocyclic double-bonded atoms should normally be planar,
// however, some large atoms may disrupt the planarity
if (!isPlanar(mol, ringAtom, exoAtom) &&
mol.getAtomicNo(ringAtom) < 14 &&
mol.getAtomicNo(exoAtom) < 14) {
continue;
}
// FIXME: pentavalent phosphorus (P) should be considered conjungated,
// but not strictly aromatic as stated in 10.1016/j.ccr.2015.02.010
if ((spOrder[exoAtom] == 2 || getAllConnAtoms(mol, exoAtom, true) == 1)) {
// check exocyclic atom bond length
int minPiCount = (getAllConnAtoms(mol, exoAtom, true) > 1) ? 1 : 0;
double averageSingleBondLength = getIdealisedBondLength(1, mol.getAtomicNo(ringAtom), mol.getAtomicNo(exoAtom), 1, minPiCount);
double averageSingleBondStd = BondLengthSet.getBondStdDev(BondLengthSet.getBondIndex(1, false, false, mol.getAtomicNo(ringAtom), mol.getAtomicNo(exoAtom), 1, minPiCount));
double averageDoubleBondLength = getIdealisedBondLength(2, mol.getAtomicNo(ringAtom), mol.getAtomicNo(exoAtom), 1, 1);
double averageDoubleBondStd = BondLengthSet.getBondStdDev(BondLengthSet.getBondIndex(2, false, false, mol.getAtomicNo(ringAtom), mol.getAtomicNo(exoAtom), 1, 1));
double calculatedLength = getBondLength(mol, ringAtom, exoAtom);
double singleBondDiff = calculatedLength - averageSingleBondLength;
double doubleBondDiff = calculatedLength - averageDoubleBondLength;
if (averageSingleBondLength < averageDoubleBondLength) {
if (DEBUG) {
System.err.println(">>> DEBUG: single '" + mol.getAtomicNo(ringAtom) + "-" +
mol.getAtomicNo(exoAtom) + "' bond seems to be shorter than a the " +
"double bond of the same type -- skipping the estimation of this bond.");
}
continue;
}
// add an exocyclic double bond is the exocyclic atom
// has no alternative way to fill its valency
if (mol.getConnAtoms(exoAtom) > 1 && !hasMetalLigandBonds(mol, exoAtom)) {
int doubleBondNeighbours = 0;
for (int j = 0; j < mol.getConnAtoms(exoAtom); j++) {
if (mol.getLowestFreeValence(mol.getConnAtom(exoAtom, j)) != 0) {
doubleBondNeighbours++;
}
}
if (doubleBondNeighbours == 1 && doubleBondDiff <= 3 * averageDoubleBondStd) {
mol.setBondOrder(mol.getBond(ringAtom, exoAtom), 2);
localAtomPi[ringAtom] += 1;
localAtomPi[exoAtom] += 1;
continue;
}
}
// skip the atom if we are more likely to find a double bond inside the ring
if (doubleBondDiff > inRingMinDeviation) continue;
if (doubleBondDiff < 0 || doubleBondDiff <= 1.5 * averageDoubleBondStd
|| StrictMath.abs(singleBondDiff) > 2 * averageSingleBondStd) {
mol.setBondOrder(mol.getBond(ringAtom, exoAtom), 2);
localAtomPi[ringAtom] += 1;
localAtomPi[exoAtom] += 1;
hasExocyclicDoubleBond[i] = true;
if (DEBUG) {
System.err.println(">>> DEBUG: placing a double bond between atoms '" +
mol.getAtomName(ringAtom) + "' and '" + mol.getAtomName(exoAtom) + "'.");
}
}
}
}
}
}
mol.ensureHelperArrays(Molecule3D.cHelperRings);
return hasExocyclicDoubleBond;
}
/**
* Assigns positive charges to metal atoms in well-known ions.
* @param mol
* The molecule that contains the metal atoms.
*/
private static void processKnownIonsInvolvingMetals(Molecule3D mol) {
for (int i = 0; i < mol.getAllAtoms(); i++) {
if (!mol.isMetalAtom(i)) continue;
if (resolveMetalHexafluoride(mol, i)) continue;
if (resolveGoldInLinearCoordination(mol, i)) continue;
// Aluminium (Al)
if (mol.getAtomicNo(i) == 13) {
// AlCl4(-)
if (mol.getAllConnAtomsPlusMetalBonds(i) == 4) {
boolean isIon = true;
for (int j = 0; j < mol.getAllConnAtomsPlusMetalBonds(i); j++) {
int connAtom = mol.getConnAtom(i, j);
// is terminal Cl atom
if (!(mol.getAtomicNo(connAtom) == 17 && mol.getAllConnAtoms(connAtom) == 0)) {
isIon = false;
break;
}
}
if (isIon) {
mol.setAtomCharge(i, 3);
}
}
}
}
}
/**
* Evaluates if the given atom serves as the central metal atom of a hexafluoride complex
* and assigns appropriate charges if needed. Only hexafluorides that can be resolved
* unambiguously are modified.
*
* @param mol
* Molecule that contains the metal atom.
* @param metalAtom
* Metal atom to be evaluated.
* @return
* Boolean value denoting if the charge of the metal atom was modified.
*/
private static boolean resolveMetalHexafluoride(Molecule3D mol, int metalAtom) {
int charge;
switch (mol.getAtomicNo(metalAtom)) {
case 13: // AlF6(3-)
case 26: // FeF6(3-)
charge = +3;
break;
case 51: // SbF6(-1)
case 79: // AuF6(-1)
charge = +5;
break;
// -1 is treated as a special value
default:
charge = -1;
}
if (charge == -1) return false;
if (!isCentralHexafluoriteIonAtom(mol, metalAtom)) return false;
mol.setAtomCharge(metalAtom, charge);
return true;
}
/**
* Evaluates if the given atom is a gold atom in a linear coordination
* and assigns appropriate charges if needed.
*
* NOTE: gold atoms appearing in linear coordination most likely have
* the oxidation state of +1.
*
* @param mol
* Molecule that contains the metal atom.
* @param metalAtom
* Metal atom to be evaluated.
* @return
* Boolean value denoting if the charge of the metal atom was modified.
*/
private static boolean resolveGoldInLinearCoordination(Molecule3D mol, int metalAtom) {
if (mol.getAtomicNo(metalAtom) != 79) return false; // Gold (Au)
if (mol.getAllConnAtomsPlusMetalBonds(metalAtom) != 2) return false;
int n1 = mol.getConnAtom(metalAtom, 0);
int n2 = mol.getConnAtom(metalAtom, 1);
// 2.80 rad ~ 160 deg
if (getAngle(mol, n1, metalAtom, n2) < 2.80) return false;
mol.setAtomCharge(metalAtom, +1);
return true;
}
/**
* Evaluates if the given atom serves as the central atom of a hexafluoride complex.
* @param mol
* Molecule that contains the atom.
* @param atom
* Atom to be evaluated.
* @return
* Boolean value denoting if given atom serves as the central atom of
* a hexafluoride complex
*/
static boolean isCentralHexafluoriteIonAtom(Molecule3D mol, int atom) {
boolean isIon = true;
if (mol.getAllConnAtomsPlusMetalBonds(atom) != 6) return false;
for (int j = 0; j < mol.getAllConnAtomsPlusMetalBonds(atom); j++) {
int connAtom = mol.getConnAtom(atom, j);
// is terminal F atom
if (!(mol.getAtomicNo(connAtom) == 9 && mol.getAllConnAtoms(connAtom) == 0)) {
isIon = false;
break;
}
}
return isIon;
}
private static void assignAtomChargesBasedOnValence(Molecule3D mol, int[] spOrder) {
// Set charges for non-metal atoms
for (int i = 0; i < mol.getAllAtoms(); i++) {
if (mol.isMetalAtom(i)) continue;
// FIXME: check and properly test if the ( mol.getAtomCharge(i) != 0 )
// can be refactored to the start of this loop
// only apply this set of rules to atoms that do not yet have a charge
if (mol.getAtomCharge(i) == 0) {
if (mol.getAtomicNo(i) == 5) { // boron (B)
if (mol.getOccupiedValence(i) == 4) {
mol.setAtomCharge(i, -1);
}
// } else if ( getAllConnAtoms(mol, i, true) == 6 ) {
// mol.setAtomAbnormalValence(i, 6);
// }
}
// Phoshorus (P)
if (mol.getAtomicNo(i) == 15) {
// Phoshorus (P) atoms with 4 neighbours and no higher order bonds have a +1 formal charge
if (getAllConnAtoms(mol, i, true) == 4 && mol.getAtomPi(i) == 0) {
mol.setAtomCharge(i, +1);
}
// Phosphorus (P) atoms with the valence of 6 have a -1 formal charge
if (mol.getOccupiedValence(i) == 6) {
mol.setAtomCharge(i, -1);
}
}
// Sulphur (S)
if (mol.getAtomicNo(i) == 16) {
// Sulphur (S) sp3 atoms with 3 neighbours have a +1 formal charge
if (spOrder[i] == 3 && getAllConnAtoms(mol, i, true) == 3) {
mol.setAtomCharge(i, +1);
}
}
// Arsernic (As) atoms with 4 neighbours and no double bonds have a +1 formal charge
if (mol.getAtomicNo(i) == 33) {
if (getAllConnAtoms(mol, i, true) == 4 && mol.getAtomPi(i) == 0) {
mol.setAtomCharge(i, +1);
}
}
}
// Hydrogen (H) ions
if (mol.getAtomicNo(i) == 1 && mol.getAllConnAtoms(i) == 0) {
if (mol.getAllConnAtomsPlusMetalBonds(i) > 0 &&
mol.isMetalAtom(mol.getConnAtom(i, 0))) {
// Hydride anion H- (usually bonded to metals)
mol.setAtomCharge(i, -1);
} else {
// Hydron cation H+
mol.setAtomCharge(i, +1);
}
}
// Removing implicit hydrogens by setting a negative charge.
// '_atom_site_attached_hydrogens' data item values are taken
// into consideration
if (mol.getImplicitHydrogens(i) > 0 || mol.getAttachedHydrogenCount(i) > 0) {
// sp2 carbon with three non-hydrogen neighbours should have a positive charge
if (mol.getAtomicNo(i) == 6 && spOrder[i] == 2 &&
mol.getConnAtoms(i) == 3 && mol.getAtomPi(i) == 0
&& mol.getAtomCharge(i) == 0 && !mol.isAromaticAtom(i)) {
mol.setAtomCharge(i, +1);
} else {
mol.setAtomCharge(i, mol.getAtomCharge(i) - mol.getImplicitHydrogens(i)
+ mol.getAttachedHydrogenCount(i));
}
}
// Adding a positive charge to the elements that have a high valency
if (mol.getOccupiedValence(i) > mol.getMaxValence(i) && mol.getAtomCharge(i) == 0) {
mol.setAtomCharge(i, mol.getAtomCharge(i) + mol.getOccupiedValence(i) - mol.getMaxValence(i));
}
}
}
/**
* Assigns positive charges to metal atoms if needed to balance out the overall
* charge of the molecule. It is assumed that the metal atoms were not assigned
* a charge up to this point.
* @param mol
* The molecule that contains the metal atoms.
*/
private static void assignMetalCharges(Molecule3D mol) {
int netCharge = 0;
for (int i = 0; i < mol.getAllAtoms(); i++) {
netCharge += mol.getAtomCharge(i);
}
if ( netCharge >= 0 ) {
return;
}
List metalAtoms = new ArrayList<>();
for (int i = 0; i < mol.getAtoms(); i++) {
if ( !mol.isMetalAtom(i) ) continue;
if ( mol.getAtomCharge(i) != 0 ) continue;
metalAtoms.add(i);
}
if ( metalAtoms.size() == 0 ) {
return;
}
// Put all of the remaining negative charge on the single metal atom
if ( metalAtoms.size() == 1 ) {
mol.setAtomCharge(metalAtoms.get(0), -Math.max(netCharge, -8));
} else {
/*
// Process the metal coordination sphere surroundings first
// All negatively charged ligands are assumed to increase the
// oxidation state by one despite their actual charge
if ( netCharge < 0 ) {
for (Integer metal : metalAtoms) {
if (maximumIonChargeReached(mol, metal)) continue;
if (!hasMetalLigandBonds(mol, metal)) continue;
for (int k = mol.getAllConnAtoms(metal); k < mol.getAllConnAtomsPlusMetalBonds(metal); k++) {
int ligandCharge = mol.getAtomCharge(mol.getConnAtom(metal, k)) < 0 ? 1 : 0;
if (ligandCharge < 0) {
int oldMetalCharge = mol.getAtomCharge(metal);
int newMetalCharge =
StrictMath.min( oldMetalCharge + StrictMath.abs(ligandCharge),
getMaximumIonCharge(mol, metal) );
mol.setAtomCharge(metal, newMetalCharge);
netCharge += newMetalCharge - oldMetalCharge;
}
}
}
}
*/
// TODO: select only atoms that do not have maximum charge
// If all metals are of the same type, valence and charge
// and the net charge is divisible by the number of metals
// atoms, distribute the net charge equally
boolean similarMetals = false;
if ( StrictMath.abs(netCharge) % metalAtoms.size() == 0 ) {
similarMetals = true;
for (Integer metalAtom : metalAtoms) {
similarMetals = similarMetals &&
mol.getAtomicNo(metalAtoms.get(0)) == mol.getAtomicNo(metalAtom) &&
mol.getAtomCharge(metalAtoms.get(0)) == mol.getAtomCharge(metalAtom) &&
mol.getAllConnAtoms(metalAtoms.get(0)) == mol.getAllConnAtoms(metalAtom);
}
}
if (similarMetals) {
int increment = Math.min( StrictMath.abs(netCharge) / metalAtoms.size(), 8 );
for (Integer metalAtom : metalAtoms) {
mol.setAtomCharge(metalAtom, increment);
}
} else {
// Otherwise try to distribute it over all of the metal atoms
boolean changeOccurred = true;
while (netCharge < 0 && changeOccurred) {
changeOccurred = false;
for (Integer metal : metalAtoms) {
int newMetalCharge = mol.getAtomCharge(metal) + 1;
if (newMetalCharge <= getMaximumIonCharge(mol, metal)) {
netCharge += 1;
mol.setAtomCharge(metal, newMetalCharge);
changeOccurred = true;
}
if (netCharge == 0) break;
}
}
}
}
// Set abnormal valence to metals that exhibit uncommon
// charges in order to generate a proper id codes
for (Integer metal : metalAtoms) {
// TODO: extend the logic to other metals as well
// Currently, the logic is only applied to Aluminium (Al).
// The code should work is other metals as well, but
// more extensive testing is needed
if (mol.getAtomicNo(metal) != 13) continue;
if (Molecule.cAtomValence[mol.getAtomicNo(metal)] == null) continue;
if (mol.getMaxValence(metal) > 0) {
mol.setAtomAbnormalValence(metal, mol.getDefaultMaxValenceUncharged(metal) - mol.getMaxValence(metal));
}
}
}
/**
* Creates a list of atom ring indexes that an atom is a
* part of for each atom of the molecule.
*
* @param mol
* The molecule that contains the atoms and rings.
* @return
* Array of atom ring index lists. The array indexes
* match the indexes of the molecule atoms. If an
* atom is not a part of any rings, null value is
* placed in the array instead of an empty list.
*/
private static ArrayList[] getAtomToRings(Molecule3D mol) {
ArrayList[] atomToRings = new ArrayList[mol.getAllAtoms()];
RingCollection ringSet = mol.getRingSet();
for ( int r = 0; r < ringSet.getSize(); r++ ) {
int[] ringAtom = ringSet.getRingAtoms(r);
for (int atom: ringAtom) {
if (atomToRings[atom] == null) {
atomToRings[atom] = new ArrayList<>();
}
atomToRings[atom].add(r);
}
}
return atomToRings;
}
/**
* Determines which atom rings in the molecule are aromatic.
* The aromaticity is determined by checking the ring geometry
* against various constraints like atom chemical type, atom
* planarity, atom bond lengths, atom bond angles.
*
* @param mol
* The molecule that contains the atom rings.
* @param spOrder
* The sp hybridization type for each atom.
* @return
* An array of boolean values marking which atom rings could
* potentially accommodate aromaticity. The indexes of the
* returned array correspond to the atom ring indexes in
* the molecule.
*/
private static boolean[] findAromaticRings(Molecule3D mol, int[] spOrder) {
RingCollection ringSet = mol.getRingSet();
boolean[] looksLikeAromaticRing = new boolean[ringSet.getSize()];
for (int ringNo = 0; ringNo < ringSet.getSize(); ringNo++) {
looksLikeAromaticRing[ringNo] = looksLikeAromaticRing(mol, ringNo, spOrder);
}
return looksLikeAromaticRing;
}
/**
* Determines is the given atom ring could potentially accommodate aromaticity.
* False positive values may be returned.
* @param mol
* Molecule that contains the atom ring.
* @param ringNo
* Index of the ring in the molecule.
* @param spOrder
* The sp hybridization type for each atom.
* @return
* Boolean values denoting is the ring could potentially accommodate aromaticity.
*/
private static boolean looksLikeAromaticRing(Molecule3D mol, int ringNo, int[] spOrder) {
if (mol.getRingSet().getRingSize(ringNo) < 3) return false;
if (mol.getRingSet().getRingSize(ringNo) > 7) return false;
int[] ringAtoms = mol.getRingSet().getRingAtoms(ringNo);
// Rings that are potentially interacting with a metal using
// pi bonds should be allowed to be a little more deformed
boolean isCoordinatedRing = isRingCoordinated(mol, ringNo);
for (int i = 0; i < ringAtoms.length; i++) {
int ringAtom = ringAtoms[i];
if (!looksLikeAromaticAtom(mol, ringAtom, spOrder[ringAtom])) return false;
int a0 = ringAtoms[(i - 1 + ringAtoms.length) % ringAtoms.length];
int a1 = ringAtoms[(i) % ringAtoms.length];
int a2 = ringAtoms[(i + 1) % ringAtoms.length];
int a3 = ringAtoms[(i + 2) % ringAtoms.length];
// Make sure that all boron, carbon, nitrogen and silicon
// atoms in the ring are planar, but allow some deviation
// due to neighbouring sulphur atoms
if ((mol.getAtomicNo(a1) == 5 ||
mol.getAtomicNo(a1) == 6 ||
mol.getAtomicNo(a1) == 7 ||
mol.getAtomicNo(a1) == 14) &&
spOrder[ringAtom] != 2) {
if (!((mol.getAtomicNo(a0) == 16 || mol.getAtomicNo(a2) == 16) &&
(getAllConnAtoms(mol, a1, true)) < 4)) {
return false;
}
}
if (!(fitsAromaticLengthConstraints(mol, a1, a2) ||
hasMetalLigandBonds(mol, a1) ||
hasMetalLigandBonds(mol, a2) ||
isCoordinatedRing)) {
if (DEBUG) {
System.err.println(">>> DEBUG: bond length between atoms '" + mol.getAtomName(a1) + "' and '" +
mol.getAtomName(a2) + "' caused ring " + ringNo + " to be marked as non aromatic.");
}
return false;
}
// check if all of the ring atoms lie in a plane
double dihedral = getDihedral(mol, a0, a1, a2, a3);
if (Math.abs(Math.toDegrees(dihedral)) > 20) {
return false;
}
if (mol.getAtomicNo(a1) == 5 || mol.getAtomicNo(a1) == 6 || mol.getAtomicNo(a1) == 7) {
if (isCoordinatedRing) continue;
for (int j = 0; j < mol.getAllConnAtoms(a1); j++) {
if (isRingAtom(mol, ringNo, mol.getConnAtom(a1, j))) continue;
dihedral = getDihedral(mol, mol.getConnAtom(a1, j), a1, a2, a3);
if (StrictMath.min(Math.abs(Math.toDegrees(dihedral)), 180 - Math.abs(Math.toDegrees(dihedral))) > 25) {
if (DEBUG) {
System.err.println(">>> DEBUG: the torsion angle involving atoms '"
+ mol.getAtomName(mol.getConnAtom(a1, j)) + "', '"
+ mol.getAtomName(a1) + "' , "
+ mol.getAtomName(a2) + "' and '"
+ mol.getAtomName(a3) + "' caused ring " + ringNo +
" to be marked as non aromatic.");
}
return false;
}
}
}
}
// excluding 3-membered rings that are part of a tetrahedron
if (ringAtoms.length == 3) {
int ringAtom = ringAtoms[0];
int exocyclicAtom = -1;
if (mol.getConnAtoms(ringAtom) == 4) {
for (int j = 0; j < mol.getConnAtoms(ringAtom); j++) {
if (!isRingAtom(mol, ringNo, mol.getConnAtom(ringAtom, j))) {
exocyclicAtom = mol.getConnAtom(ringAtom, j);
break;
}
}
}
return (mol.getBond(ringAtoms[1], exocyclicAtom) != -1) &&
(mol.getBond(ringAtoms[2], exocyclicAtom) != -1);
}
return true;
}
/**
* Evaluates if the given atom could potentially participate in aromaticity.
* False positive values may be returned.
*
* @param mol
* Molecule that contains the atom.
* @param atom
* Index of the atom.
* @param spOrder
* The sp hybridization type for each atom.
* @return
* Boolean value denoting if the atom could participate in aromaticity.
*/
private static boolean looksLikeAromaticAtom(Molecule3D mol, int atom, int spOrder) {
// boron (B), carbon (C), nitrogen (N), oxygen (O), silicon (Si)
if (mol.getAtomicNo(atom) < 15) {
if (getAllConnAtoms(mol, atom, true) > 3) return false;
}
// germanium (Ge)
if (mol.getAtomicNo(atom) == 32) {
if (getAllConnAtoms(mol, atom, true) > 3) return false;
if (spOrder != 2) return false;
}
if (getAllConnAtoms(mol, atom, true) > 4) return false;
// Atoms that already have double bonds cannot participate in aromaticity
// TODO: check if it is only true in some cases, sulphur (S) for example
if (mol.getAtomPi(atom) != 0) return false;
return true;
}
/**
* Calculates the hybridization state for each atom in a molecule.
*
* @param mol
* Molecule that contains the atoms.
* @return
* Array of integers with values {1,2,3,4,5} that correspond to
* the the {sp1, sp2, sp3, sp3d1, sp3d2} hybridization states
* accordingly.
*/
private static int[] calculateHybridizationStates(Molecule3D mol) {
int[] spOrder = new int[mol.getAllAtoms()];
for (int i = 0; i < mol.getAllAtoms(); i++) {
spOrder[i] = calculateHybridizationState(mol, i);
if (DEBUG) {
System.err.println("Atom '" + mol.getAtomName(i) + "' was assigned hybridization state " + spOrder[i]);
}
}
return spOrder;
}
/**
* Calculates the hybridization state for the given atom.
*
* @param mol
* Molecule that contains the atom.
* @return
* Integers from the set of {1,2,3,4,5} that correspond to
* the the {sp1, sp2, sp3, sp3d1, sp3d2} hybridization states
* accordingly.
*/
private static int calculateHybridizationState(Molecule3D mol, int i) {
int spOrder = 0; // special value denoting an unrecognised hybridization state
if (getAllConnAtoms(mol, i, true) <= 1) {
// TODO: ask if this can always be assumed
if (mol.getAtomicNo(i) == 8) {
spOrder = 2;
} else {
spOrder = 1;
}
} else if (getAllConnAtoms(mol, i, true) == 2) {
if (getAllConnAtoms(mol, i, false) <= 1) {
spOrder = 2;
} else {
double angle = getAngle(mol, mol.getConnAtom(i, 0), i, mol.getConnAtom(i, 1));
if (Math.abs(angle - Math.PI) < Math.PI / 6) {
spOrder = 1;
} else {
spOrder = 2;
}
}
} else if (getAllConnAtoms(mol, i, true) == 3) {
if (getAllConnAtoms(mol, i, false) < 2) {
spOrder = 0; // TODO: to be changed to 2 or 3
} else if (getAllConnAtoms(mol, i, false) < 3) {
spOrder = 2;
} else {
Coordinates origin = mol.getCoordinates(i);
Coordinates normal = calculatePlaneNormal(origin,
mol.getCoordinates(mol.getConnAtom(i, 0)),
mol.getCoordinates(mol.getConnAtom(i, 1)));
if (normal.distSq() > 0) {
if (calculatePointToPlaneDistance(normal, origin, mol.getCoordinates(mol.getConnAtom(i, 2))) < 0.4) {
spOrder = 2;
} else {
spOrder = 3;
}
} else {
spOrder = 3;
}
}
} else if (getAllConnAtoms(mol, i, true) == 4) {
spOrder = 3;
} else if (getAllConnAtoms(mol, i, true) == 5) {
spOrder = 4;
} else if (getAllConnAtoms(mol, i, true) == 6) {
spOrder = 5;
}
return spOrder;
}
/**
* Assigns bond orders and charges to hydrogen (H) atoms that
* are bonded to two non-metal(lic) atoms.
*
* @param mol
* The molecule that contains the shared hydrogen atoms.
*/
private static void resolveSharedHydrogenAtoms(Molecule3D mol) {
// Shared hydrogen atoms are handled in one of the two ways
// depending on the bond lengths to the neighbouring atoms:
// - in case one of the bonds is significantly shorter,
// the longer bond gets removed;
// - in case neither of the bonds is significantly shorter,
// both bonds are marked as metal-ligand (coordinate) bonds
// and the hydrogen atom is assigned a +1 charge.
for (int i = 0; i < mol.getAllAtoms(); i++) {
// R-H-R bridging hydrogen atom
if (mol.getAtomicNo(i) != 1) continue;
if (mol.getConnAtoms(i) != 2) continue;
// boron compounds must be handled separately
if (mol.getAtomicNo(mol.getConnAtom(i, 0)) == 5) continue;
if (mol.getAtomicNo(mol.getConnAtom(i, 1)) == 5) continue;
int bond1 = mol.getConnBond(i, 0);
int bond2 = mol.getConnBond(i, 1);
double bondLength1 = getBondLength(mol, bond1);
double bondLength2 = getBondLength(mol, bond2);
if (StrictMath.abs(bondLength1 - bondLength2) < 0.1) {
mol.setBondType(bond1, Molecule3D.cBondTypeMetalLigand);
mol.setBondType(bond2, Molecule3D.cBondTypeMetalLigand);
mol.setAtomCharge(i, +1);
} else {
// delete the longer single bond
if (bondLength1 > bondLength2) {
// mol.setBondType(bond1, Molecule3D.cBondTypeMetalLigand);
mol.markBondForDeletion(bond1);
} else {
// mol.setBondType(bond2, Molecule3D.cBondTypeMetalLigand);
mol.markBondForDeletion(bond2);
}
}
}
mol.deleteMarkedAtomsAndBonds();
mol.ensureHelperArrays(Molecule.cHelperNeighbours);
}
/**
* Recognises known ions in a molecule and assigns the appropriate
* bond orders and atom charges if needed.
*
* @param mol
* Molecule that contains the ions.
*/
private static void resolveKnownIons(Molecule3D mol) {
for (int i = 0; i < mol.getAtoms(); i++) {
if (mol.getAtomicNo(i) != 6) continue;
if (calculateHybridizationState(mol, i) != 2) continue;
if (mol.getAllConnAtoms(i) != 3) continue;
if (mol.getAttachedHydrogenCount(i) != 0) continue;
resolveTricyanomethanideIon(mol, i);
}
}
/**
* Evaluates if the given atom acts as the central atom of the tricyanomethanide
* ion ([C-](C#N)(C#N)(C#N)) and assigns appropriate bond orders and atom charges
* if needed.
*
* @param mol
* Molecule that contains the atoms.
* @param cCentralAtom
* Index of the central carbon atom to be evaluated.
* @return
* Boolean value denoting if any bond order or atom charge
* modifications were applied.
*/
private static boolean resolveTricyanomethanideIon(Molecule3D mol, int cCentralAtom) {
if (mol.getAtomicNo(cCentralAtom) != 6) return false;
if (calculateHybridizationState(mol, cCentralAtom) != 2) return false;
if (mol.getAllConnAtoms(cCentralAtom) != 3) return false;
if (mol.getAttachedHydrogenCount(cCentralAtom) != 0) return false;
int[] cAtoms = new int[3];
int[] nAtoms = new int[3];
for (int i = 0; i < mol.getAllConnAtoms(cCentralAtom); i++) {
int cAtom = mol.getConnAtom(cCentralAtom, i);
if (mol.getAllConnAtoms(cAtom) != 2) return false;
for (int j = 0; j < mol.getAllConnAtoms(cAtom); j++) {
int nAtom = mol.getConnAtom(cAtom, j);
if (cCentralAtom == nAtom) continue;
if (!looksLikeCyanoGroup(mol, cAtom, nAtom)) return false;
cAtoms[i] = cAtom;
nAtoms[i] = nAtom;
}
}
// [C-](C#N)(C#N)(C#N)
for (int j = 0; j < 3; j++) {
mol.setBondOrder(mol.getBond(cAtoms[j], nAtoms[j]), 3);
mol.setAtomCharge(cCentralAtom, -1);
}
mol.ensureHelperArrays(Molecule3D.cHelperNeighbours);
return true;
}
/**
* Evaluates if the given atom ring looks like a tetraatomic aromatic ring ion
* and assigns appropriate bond orders and atom charges if needed.
*
* @param mol
* Molecule that contains the atoms ring.
* @param ringNo
* Index of the ring in the molecule.
* @return
* Boolean value denoting if any bond order or atom charge
* modifications were applied.
*/
private static boolean resolveTetraatomicAromaticRingIons(Molecule3D mol, int ringNo) {
if (!looksLikeTetraatomicAromaticRingIon(mol, ringNo)) return false;
boolean isChalcogen = true;
int atomicNo = mol.getAtomicNo(mol.getRingSet().getRingAtoms(ringNo)[0]);
switch (atomicNo) {
case 7: // nitrogen (N)
case 15: // phosphorus (P);
case 33: // arsenic (As);
isChalcogen = false;
break;
case 8: // oxygen (O)
case 16: // sulphur (S)
case 34: // selenium (Se); tetraselenium ion ([Se]1[Se+]=[Se+][Se]1)
case 52: // tellurium (Te)
break;
default:
return false;
}
int[] ringBonds = mol.getRingSet().getRingBonds(ringNo);
int shortestBondIndex = 0;
double shortestBondLength = getBondLength(mol, ringBonds[shortestBondIndex]);
for (int i = 0; i < ringBonds.length; i++) {
if (getBondLength(mol, ringBonds[i]) < shortestBondLength) {
shortestBondLength = getBondLength(mol, ringBonds[i]);
shortestBondIndex = i;
}
}
mol.setBondOrder(ringBonds[shortestBondIndex], 2);
int dbAtom1 = mol.getBondAtom(0, ringBonds[shortestBondIndex]);
int dbAtom2 = mol.getBondAtom(1, ringBonds[shortestBondIndex]);
if (isChalcogen) {
// chalcogen ring, i.e. [Se]1[Se][Se+]=[Se+]1
mol.setAtomCharge(dbAtom1, +1);
mol.setAtomCharge(dbAtom2, +1);
} else {
// pnictogen ring, i.e. [P-]1[P-][P]=[P]1
for (int atom : mol.getRingSet().getRingAtoms(ringNo)) {
if (atom != dbAtom1 && atom != dbAtom2) {
mol.setAtomCharge(atom, -1);
}
}
}
return true;
}
/**
* Evaluates if the given atom ring looks like a tetraatomic aromatic ring ion,
* i.e. the tetraselenium ion ([Se]1[Se+]=[Se+][Se]1).
*
* @param mol
* Molecule that contains the atoms ring.
* @param ringNo
* Index of the ring in the molecule.
* @return
* Boolean value denoting if the ring looks like a tetraatomic aromatic
* ring ion.
*/
private static boolean looksLikeTetraatomicAromaticRingIon(Molecule3D mol, int ringNo) {
int[] ringAtoms = mol.getRingSet().getRingAtoms(ringNo);
if (ringAtoms.length != 4) return false;
// rings must be homocyclic
int atomicNo = mol.getAtomicNo(ringAtoms[0]);
for (int atom : ringAtoms) {
if (mol.getAtomicNo(atom) != atomicNo) return false;
if (getAllConnAtoms(mol, atom, true) != 2) return false;
if (mol.getAtomRingCount(atom, 7) != 1) return false;
}
return true;
}
/**
* Recognises functional groups in a molecule and sets the bonds and charges accordingly.
* @param mol
* The molecule that contains the functional groups.
* @param spOrder
* Array that contains the orbital hybridization type for each atom.
*/
private static void recogniseFunctionalGroups(Molecule3D mol, int[] spOrder) {
// TODO: consider the atom type of terminal bonds when sorting
// TODO: ask about the best way to keep track of new double bonds
// calling molecule helpers to ensure pi count after each double
// bond assignment is inefficient so a local array is used instead
int[] localAtomPiCount = getAtomPiCount(mol);
for(int i = 0; i < mol.getAtoms(); i++) {
if ( mol.getAtomCharge(i) != 0 ) continue;
if ( mol.getAtomRadical(i) != Molecule.cAtomRadicalStateNone ) continue;
if ( spOrder[i] == 2 && localAtomPiCount[i] != 0 && mol.getAtomicNo(i) < 7 ) continue;
// Terminal oxygen atoms capable of double bonds
List terminalOxygenAtoms = new LinkedList<>();
// All terminal atoms capable of double bonds (including the oxygen atoms)
List terminalDoubleBondAtoms = new LinkedList<>();
int oxygenCount = 0;
int terminalAtomCount = 0;
for (int j = 0; j < mol.getConnAtoms(i); j++ ) {
int neighbour = mol.getConnAtom(i, j);
boolean isTerminal = ( getAllConnAtoms(mol, neighbour, true) == 1 );
boolean isOxygen = ( mol.getAtomicNo(neighbour) == 8 );
if (isOxygen) oxygenCount++;
if ( isTerminal && isOxygen ) {
terminalOxygenAtoms.add(j);
}
if ( mol.getAtomicNo(neighbour) == 1 ) continue;
if ( isTerminal ) {
terminalAtomCount++;
// for now only consider terminal N, O, P and S for double bonds
if ( localAtomPiCount[neighbour] == 0 &&
( mol.getAtomicNo(neighbour) == 7 ||
mol.getAtomicNo(neighbour) == 8 ||
mol.getAtomicNo(neighbour) == 15 ||
mol.getAtomicNo(neighbour) == 16 ) ) {
terminalDoubleBondAtoms.add(j);
}
}
}
if ( getAllConnAtoms(mol, i, true) == 1 ) {
if (mol.getConnAtoms(i) == 0) continue;
int a1 = mol.getConnAtom(i, 0);
if (resolveDiatomicPnictogenIons(mol, i)) continue;
if (resolveDiatomicChalcogenIons(mol, i)) continue;
if ( mol.getAtomicNo(i) == 6 ) {
// R-N(+)=-=C(-)
if (mol.getAtomicNo(a1) == 7 && spOrder[a1] == 1) {
mol.setBondOrder(mol.getBond(i, a1), 3);
mol.setAtomCharge(i, -1);
mol.setAtomCharge(a1, +1);
localAtomPiCount[i] += 2;
localAtomPiCount[a1] += 2;
}
} else if ( mol.getAtomicNo(i) == 7 ) {
if (mol.getAtomicNo(a1) == 16 ) {
if ( getAllConnAtoms(mol, a1, true) == 1 ) {
if ( hasMetalLigandBonds(mol, i) ) {
mol.setBondOrder(mol.getBond(i, a1), 2);
localAtomPiCount[i]++;
localAtomPiCount[a1]++;
mol.setAtomCharge(i, -1);
} else {
mol.setBondOrder(mol.getBond(i, a1), 3);
localAtomPiCount[i] = localAtomPiCount[i] + 2;
localAtomPiCount[a1] = localAtomPiCount[a1] + 2;
mol.setAtomCharge(a1, 1);
}
} else if (spOrder[a1] == 1) {
mol.setBondOrder(mol.getBond(i, a1), 2);
localAtomPiCount[i]++;
localAtomPiCount[a1]++;
} else if (mol.getAtomicNo(a1) == 16 && spOrder[a1] == 3) {
// R-S#N
// FIXME: eventually replace the hardcoded S#N bond length value once enough statistics
// are gathered
if (getBondLength(mol, mol.getBond(i, a1)) < 1.47) {
mol.setBondOrder(mol.getBond(i, a1), 3);
localAtomPiCount[i] = localAtomPiCount[i] + 2;
localAtomPiCount[a1] = localAtomPiCount[a1] + 2;
}
}
}
} else if ( mol.getAtomicNo(i) == 8 ) {
// TODO: change O-C distances with proper values
if ( mol.getAtomicNo(a1) == 6 ) {
// C-=-O
if ( getAllConnAtoms(mol, a1, true) == 1 ) {
mol.setBondOrder(mol.getBond(i, a1), 3);
mol.setAtomCharge(a1, -1);
mol.setAtomCharge(i, +1);
localAtomPiCount[i] += 2;
localAtomPiCount[a1] += 2;
}
} else if ( mol.getAtomicNo(a1) == 7 ) {
// TODO: add the correct processing of the nitrosonium ion (N-=-O+)
// N=O ion
if ( getAllConnAtoms(mol, a1, true) == 1 ) {
mol.setBondOrder(mol.getBond(i, a1), 2);
mol.setAtomRadical(a1, Molecule3D.cAtomRadicalStateD);
localAtomPiCount[i]++;
localAtomPiCount[a1]++;
}
} else if ( mol.getAtomicNo(a1) == 8 ) {
// O-O ions
double sigma = 0.05;
if ( getAllConnAtoms(mol, a1, true) == 1 ) {
if (getBondLength(mol, i, a1) < 1.21 + sigma) {
mol.setBondOrder(mol.getBond(i, a1), 2);
localAtomPiCount[i]++;
localAtomPiCount[a1]++;
} else if (getBondLength(mol, i, a1) < 1.28 + sigma) {
mol.setBondOrder(mol.getConnBond(i, 0), 1);
mol.setAtomCharge(i, -1);
mol.setAtomRadical(a1, Molecule3D.cAtomRadicalStateD);
} else {
mol.setBondOrder(mol.getConnBond(i, 0), 1);
mol.setAtomCharge(i, -1);
mol.setAtomCharge(a1, -1);
}
}
}
}
} else if ( mol.getAllConnAtoms(i) == 2 ) {
if ( localAtomPiCount[i] == 2 ) continue;
// assigning neighbour atoms ordered by chemical number
int a1 = mol.getAtomicNo(mol.getConnAtom(i, 0)) <=
mol.getAtomicNo(mol.getConnAtom(i, 1)) ?
mol.getConnAtom(i, 0) : mol.getConnAtom(i, 1);
int a2 = a1 != mol.getConnAtom(i, 0) ?
mol.getConnAtom(i, 0) : mol.getConnAtom(i, 1);
if ( mol.getAtomicNo(i) == 5 ) {
if ( spOrder[i] == 1 ) {
if ( terminalDoubleBondAtoms.size() == 2 ) {
mol.setBondOrder(mol.getBond(i, a1), 2);
mol.setBondOrder(mol.getBond(i, a2), 2);
localAtomPiCount[i] += 2;
localAtomPiCount[a1] += 1;
localAtomPiCount[a2] += 1;
mol.setAtomCharge(i, -1);
}
}
} else if ( mol.getAtomicNo(i) == 6 ) {
if ( spOrder[i] == 1 ) {
if (terminalDoubleBondAtoms.size() == 1) {
int terminalAtom = mol.getConnAtom(i, terminalDoubleBondAtoms.get(0));
int nonTerminalAtom = (terminalAtom != a1) ? a1 : a2;
// cyano group (R-C#N)
if (resolveCyanoGroup(mol, i, terminalAtom)) {
localAtomPiCount[i] += 2;
localAtomPiCount[terminalAtom] += 2;
continue;
}
// isocyanate group (R-N=C=O)
if (resolveIsocyanateGroup(mol, i, nonTerminalAtom, terminalAtom)) {
localAtomPiCount[i] += 2;
localAtomPiCount[terminalAtom] += 1;
localAtomPiCount[nonTerminalAtom] += 1;
continue;
}
// isothiocyanate group (R-N=C=S)
if (resolveIsothiocyanateGroup(mol, i, nonTerminalAtom, terminalAtom)) {
localAtomPiCount[i] += 2;
localAtomPiCount[terminalAtom] += 1;
localAtomPiCount[nonTerminalAtom] += 1;
continue;
}
} else if (terminalDoubleBondAtoms.size() == 2) {
// thiocyanate ion (S=C=N[-])
if (resolverThiocyanateIon(mol, i, a2, a1)) {
localAtomPiCount[i] += 2;
localAtomPiCount[a1] += 1;
localAtomPiCount[a2] += 1;
continue;
}
// atoms are of the same species
if (mol.getAtomicNo(a1) == mol.getAtomicNo(a2)) {
mol.setBondOrder(mol.getBond(i, a1), 2);
mol.setBondOrder(mol.getBond(i, a2), 2);
localAtomPiCount[i] += 2;
localAtomPiCount[a1] += 1;
localAtomPiCount[a2] += 1;
}
}
}
} else if ( mol.getAtomicNo(i) == 7 ) {
if ( spOrder[i] == 1 ) {
if (terminalDoubleBondAtoms.size() > 0) {
// azide anion ([N-]=[N+]=[N-])
if (resolveAzideIon(mol, i)) {
localAtomPiCount[i] += 2;
localAtomPiCount[a1] += 1;
localAtomPiCount[a2] += 1;
continue;
}
// azide group (R[N]=[N+]=[N-])
if (resolveAzideGroup(mol, i)) {
localAtomPiCount[i] = mol.getAtomPi(i);
localAtomPiCount[a1] = mol.getAtomPi(a1);
localAtomPiCount[a2] = mol.getAtomPi(a2);
continue;
}
if ((mol.getAtomicNo(a1) == 6 && mol.getAtomicNo(a2) == 7) ||
(mol.getAtomicNo(a1) == 7 && mol.getAtomicNo(a2) == 6)) {
int carbonAtom = (mol.getAtomicNo(a1) == 6) ? a1 : a2;
int nitrogenAtom = (mol.getAtomicNo(a1) == 7) ? a1 : a2;
if ( getAllConnAtoms(mol, nitrogenAtom, true) == 1 &&
spOrder[carbonAtom] == 2) {
double doubleNNBond = getIdealisedBondLength(2, 7, 7, 1, 2);
double tripleNNBond = getIdealisedBondLength(3, 7, 7, 2, 2);
double singleCNBond = getIdealisedBondLength(1, 6, 7, 1, 2);
double doubleCNBond = getIdealisedBondLength(2, 6, 7, 1, 2);
double NNBondLength = getBondLength(mol, i, nitrogenAtom);
double CNBondLength = getBondLength(mol, i, carbonAtom);
double C2N2N = StrictMath.max( doubleNNBond - NNBondLength, 0 ) +
StrictMath.max( doubleCNBond - CNBondLength, 0 );
double C1N3N = StrictMath.max( tripleNNBond - NNBondLength, 0 ) +
StrictMath.max( singleCNBond - CNBondLength, 0 );
// RC-[N+]#N
if ( C1N3N < C2N2N ) {
mol.setBondOrder(mol.getBond(i, nitrogenAtom), 3);
mol.setAtomCharge(i, +1);
localAtomPiCount[i] += 2;
localAtomPiCount[nitrogenAtom] += 1;
} else {
// RC=[N+]=[N-]
mol.setBondOrder(mol.getBond(i, nitrogenAtom), 2);
mol.setBondOrder(mol.getBond(i, carbonAtom), 2);
mol.setAtomCharge(i, +1);
mol.setAtomCharge(nitrogenAtom, -1);
localAtomPiCount[i] += 2;
localAtomPiCount[nitrogenAtom] += 1;
localAtomPiCount[carbonAtom] += 1;
}
}
// O=N=O, S=N=S
} else if ( terminalOxygenAtoms.size() == 2 ||
mol.getAtomicNo(a1) == 16 && mol.getAtomicNo(a2) == 16 ) {
mol.setAtomCharge(i, +1);
for (Integer terminalDoubleBondAtom : terminalDoubleBondAtoms) {
mol.setBondOrder(mol.getConnBond(i, terminalDoubleBondAtom), 2);
localAtomPiCount[i] += 1;
localAtomPiCount[mol.getConnAtom(i, terminalDoubleBondAtom)] += 1;
}
}
} else {
// R-P=[N+]=P-R
if ( mol.getAtomicNo(a1) == 15 && spOrder[a1] == 3 &&
mol.getAtomicNo(a2) == 15 && spOrder[a2] == 3 ) {
// cyclophosphaneses
if ( mol.isRingAtom(i) ) {
double NPBondLength1 = getBondLength(mol, i, mol.getConnAtom(i, 0));
double NPBondLength2 = getBondLength(mol, i, mol.getConnAtom(i, 1));
if ( NPBondLength1 < NPBondLength2 && localAtomPiCount[mol.getConnAtom(i, 0)] == 0) {
mol.setBondOrder(mol.getConnBond(i, 0), 2);
localAtomPiCount[i] += 1;
localAtomPiCount[mol.getConnAtom(i, 0)] += 1;
} else if ( localAtomPiCount[mol.getConnAtom(i, 1)] == 0 ) {
mol.setBondOrder(mol.getConnBond(i, 1), 2);
localAtomPiCount[i] += 1;
localAtomPiCount[mol.getConnAtom(i, 1)] += 1;
}
// R-P=N(+)=P-R (bridging)
} else {
mol.setBondOrder(mol.getConnBond(i, 0), 2);
mol.setBondOrder(mol.getConnBond(i, 1), 2);
mol.setAtomCharge(i, +1);
localAtomPiCount[i] += 2;
localAtomPiCount[mol.getConnAtom(i, 0)] += 1;
localAtomPiCount[mol.getConnAtom(i, 1)] += 1;
}
}
}
} else if ( spOrder[i] == 2 ) {
if (terminalOxygenAtoms.size() == 1) {
if (resolveNitrosoGroup(mol, i)) {
localAtomPiCount[i] += 1;
localAtomPiCount[mol.getConnAtom(i, terminalOxygenAtoms.get(0))] += 1;
}
} else {
// R-P=N(+)-P-R
if (mol.getAtomicNo(a1) != 15) continue;
if (spOrder[a1] != 3) continue;
if (mol.getAtomicNo(a2) != 15) continue;
if (spOrder[a2] != 3) continue;
// skip if both phosphorus atoms already have a double bond
if (localAtomPiCount[a1] != 0 && localAtomPiCount[a2] != 0) continue;
// select the shortest bond out of the available ones
int doubleBondAtom = a1;
if (localAtomPiCount[a1] != 0) {
doubleBondAtom = a2;
} else if (localAtomPiCount[a2] == 0) {
if (getBondLength(mol, mol.getBond(i, a1)) > getBondLength(mol, mol.getBond(i, a2))) {
doubleBondAtom = a2;
}
}
mol.setBondOrder(mol.getBond(i, doubleBondAtom), 2);
localAtomPiCount[i] += 1;
localAtomPiCount[doubleBondAtom] += 1;
}
}
// S
} else if ( mol.getAtomicNo(i) == 16 ) {
if (spOrder[i] == 2) {
// SO2
if (terminalOxygenAtoms.size() > 0) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalOxygenAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 1);
mol.setAtomCharge(i, +1);
}
}
// F, Cl, Br, I
} else if ( mol.getAtomicNo(i) == 7 ||
mol.getAtomicNo(i) == 17 ||
mol.getAtomicNo(i) == 35 ||
mol.getAtomicNo(i) == 53 ) {
if (terminalDoubleBondAtoms.size() > 0) {
// F is not capable of a double bond
if (mol.getAtomicNo(i) == 7) {
mol.setAtomCharge(i, +1);
} else {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalDoubleBondAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 1);
}
} else {
// R-X-R (linear)
if (spOrder[i] == 1) {
mol.setAtomCharge(i, -1);
// R-X-R (bent)
} else if (spOrder[i] == 2) {
mol.setAtomCharge(i, +1);
}
}
}
} else if ( mol.getAllConnAtoms(i) == 3 ) {
int a, b;
int a1 = mol.getConnAtom(i, 0);
int a2 = mol.getConnAtom(i, 1);
int a3 = mol.getConnAtom(i, 2);
if (mol.getAtomicNo(i) == 6 && spOrder[i] == 2) {
if (mol.getAtomRingSize(i) > 0) {
if (terminalDoubleBondAtoms.isEmpty()) continue;
// thioketone in TTF: C(R)(R)(=S)
if (mol.getAtomRingSize(i) != 5) continue;
if (mol.getAtomicNo(a1) != 16) continue;
if (mol.getAtomicNo(a2) != 16) continue;
if (mol.getAtomicNo(a3) != 16) continue;
mol.setBondOrder(mol.getConnBond(i, terminalDoubleBondAtoms.get(0)), 2);
} else {
//C(N)(N)(=O)
//C(R)(O)(=O)
if (terminalOxygenAtoms.size() >= 2) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalOxygenAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 1);
continue;
}
//C(R)(OR)(=O)
a = connectedAtom(mol, i, 8, 2, 0, 0);
b = connectedBond(mol, i, 8, 1);
if (a >= 0 && b >= 0) {
mol.setBondOrder(b, 2);
continue;
}
//C(R)(SR)(=O)
a = connectedAtom(mol, i, 16, 2, 0, 0);
b = connectedBond(mol, i, 8, 1);
if (a >= 0 && b >= 0) {
mol.setBondOrder(b, 2);
continue;
}
//C(R)(NR)(=O)
a = connectedAtom(mol, i, 7, 3, 0, 0);
b = connectedBond(mol, i, 8, 1);
if (a >= 0 && b >= 0) {
mol.setBondOrder(b, 2);
continue;
}
//C(R)(SR)(=S)
a = connectedAtom(mol, i, 16, 2, 0, 0);
b = connectedBond(mol, i, 16, 1);
if (a >= 0 && b >= 0) {
mol.setBondOrder(b, 2);
continue;
}
//C(R)(NR)(=S)
a = connectedAtom(mol, i, 7, 3, 0, 0);
b = connectedBond(mol, i, 16, 1);
if (a >= 0 && b >= 0) {
mol.setBondOrder(b, 2);
continue;
}
//C(CR)(N)(=N)
if (terminalDoubleBondAtoms.size() == 2 &&
mol.getAtomicNo(terminalDoubleBondAtoms.get(0)) == 7 &&
mol.getAtomicNo(terminalDoubleBondAtoms.get(1)) == 7 &&
(mol.getAtomicNo(a1) == 6 ||
mol.getAtomicNo(a2) == 6 ||
mol.getAtomicNo(a3) == 6)) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalDoubleBondAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 1);
}
}
} else if( mol.getAtomicNo(i) == 7 ) {
// TODO: evaluate the probability of a certain SP order from the geometry
// slightly disordered nitrate ion
if ( spOrder[i] == 2 ) {
//N(R)(R)C=O -> Amide
a = connectedAtom(mol, i, 6, 2, 8, 1);
b = connectedBond(mol, a, 8, 1);
if (a >= 0 && b >= 0) {
mol.setBondOrder(b, 2);
localAtomPiCount[i] += 1;
localAtomPiCount[b] += 1;
} else
// N(O)(=O) -> Nitro
// nitrate ion
// NF2O(1+)
if ( terminalAtomCount > 1 && terminalOxygenAtoms.size() > 0) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalOxygenAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 1);
mol.setAtomCharge(i, +1);
}
}
} else if (mol.getAtomicNo(i) == 16) {
if (spOrder[i] == 2) {
// sulphur trioxide
if (terminalOxygenAtoms.size() > 0) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalOxygenAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 3);
}
} else if (spOrder[i] == 3) {
// sulphite ion
if (terminalOxygenAtoms.size() > 0) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalOxygenAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 1);
}
}
// Cl, Br, I
} else if ( mol.getAtomicNo(i) == 17 ||
mol.getAtomicNo(i) == 35 ||
mol.getAtomicNo(i) == 53 ) {
if ( spOrder[i] == 3 && terminalDoubleBondAtoms.size() > 0 ) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalDoubleBondAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 2);
}
} else if ( mol.getAtomicNo(i) == 34 ) {
if (oxygenCount == 3) {
if (spOrder[i] == 2) {
// selenium trioxide
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalOxygenAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 2);
} else if (spOrder[i] == 3) {
// selenite ion
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalOxygenAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 1);
}
} else {
mol.setAtomCharge(i, +1);
}
// Te
} else if( mol.getAtomicNo(i) == 52 ) {
mol.setAtomCharge(i, +1);
}
} else if( mol.getAllConnAtoms(i) == 4 ) {
// phosphorus (P) and arsenic (As)
if( mol.getAtomicNo(i) == 15 || mol.getAtomicNo(i) == 33 ) {
if (terminalDoubleBondAtoms.isEmpty()) continue;
// skip is the atom already has a double bond
if (localAtomPiCount[i] != 0) continue;
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalDoubleBondAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 1);
localAtomPiCount[i] += 1;
localAtomPiCount[mol.getConnAtom(i, sortedAtoms[0])] += 1;
} else if( mol.getAtomicNo(i)==16 ) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalDoubleBondAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 2);
// Cl, Br, I
} else if( mol.getAtomicNo(i) == 17 ||
mol.getAtomicNo(i) == 35 ||
mol.getAtomicNo(i) == 53 ) {
if ( terminalDoubleBondAtoms.size() > 0 ) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalDoubleBondAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 3);
} else {
// Set a negative charge to hypervalent atom
mol.setAtomCharge(i, -1);
}
} else if( mol.getAtomicNo(i) == 34 ) {
if ( oxygenCount == 4 && terminalOxygenAtoms.size() > 1 ) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalOxygenAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 2);
}
}
} else if( mol.getAllConnAtoms(i) == 5 ) {
if ( mol.getAtomicNo(i) == 14 ) { // Silicon (Si)
mol.setAtomCharge(i, -1);
} else if ( mol.getAtomicNo(i) == 53 ) { // Iodine (I)
if ( terminalOxygenAtoms.size() > 0 ) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalOxygenAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 2);
}
}
} else if ( mol.getAllConnAtoms(i) == 6) {
if (mol.getAtomicNo(i) == 14) { // Silicon (Si)
mol.setAtomCharge(i, -2);
} else if (mol.getAtomicNo(i) == 32) { // Germanium (Ge)
mol.setAtomCharge(i, -2);
} else if (mol.getAtomicNo(i) == 33) { // Arsenic (As)
mol.setAtomCharge(i, -1);
} else if (mol.getAtomicNo(i) == 34) { // Selenium (Se)
if ( terminalDoubleBondAtoms.size() > 0 ) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalDoubleBondAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 1);
} else {
mol.setAtomCharge(i, -2);
}
} else if( mol.getAtomicNo(i) == 52 ) { // Tellurium (Te)
// telluric acid has a net charge of 0
// even though it has 6 neighbours
if ( terminalDoubleBondAtoms.size() == 0 && oxygenCount == 0 ) {
mol.setAtomCharge(i, -2);
}
} else if( mol.getAtomicNo(i) == 53 ) { // Iodine (I)
if ( terminalDoubleBondAtoms.size() > 0 ) {
int[] sortedAtoms = sortAtomsByBondLength(mol, i, terminalDoubleBondAtoms);
setTerminalDoubleBonds(mol, i, sortedAtoms, 1);
} else {
mol.setAtomCharge(i, +1);
}
}
}
}
mol.ensureHelperArrays(Molecule3D.cHelperRings);
}
/**
* Evaluates if the given nitrogen atom serves as a part of a nitroso (R-N=O)
* group and assigns appropriate bond orders and atom charges if needed.
*
* Current implementation uses the following assumptions that are not explicitly checked:
* - the provided nitrogen atom has only two neighbour atoms;
* - one of the neighbour atoms is an oxygen atom capable of a double bond;
*
* @param mol
* Molecule that contains the atoms.
* @param nAtom
* Index of the nitrogen atom to be evaluated.
* @return
* Boolean value denoting if any bond order or atom charge
* modifications were applied.
*/
private static boolean resolveNitrosoGroup(Molecule3D mol, int nAtom) {
int oAtom;
int rAtom;
if (mol.getAtomicNo(mol.getConnAtom(nAtom, 0)) == 8) {
oAtom = mol.getConnAtom(nAtom, 0);
rAtom = mol.getConnAtom(nAtom, 1);
} else {
oAtom = mol.getConnAtom(nAtom, 1);
rAtom = mol.getConnAtom(nAtom, 0);
}
// c-nitroso group (RC-N=O)
if (mol.getAtomicNo(rAtom) == 6) {
mol.setBondOrder(mol.getBond(nAtom, oAtom), 2);
return true;
}
return false;
}
/**
* Evaluates if the given atoms could potentially comprise a cyano group (R-C#N).
*
* @param mol
* Molecule that contains the atoms.
* @param cAtom
* Index of the carbon atom to be evaluated.
* @param nAtom
* Index of the nitrogen atom to be evaluated.
* @return
* Boolean value denoting if the provided atoms could potentially
* comprise the group.
*/
private static boolean looksLikeCyanoGroup(Molecule3D mol, int cAtom, int nAtom) {
if (mol.getAtomicNo(cAtom) != 6) return false;
if (calculateHybridizationState(mol, cAtom) != 1) return false;
if (mol.getAtomicNo(nAtom) != 7) return false;
if (mol.getAllConnAtoms(nAtom) != 1) return false;
if (mol.getAttachedHydrogenCount(nAtom) != 0) return false;
return true;
}
/**
* Evaluates if the given atoms comprise a cyano group (R-C#N)
* and assigns appropriate bond orders and atom charges if needed.
*
* @param mol
* Molecule that contains the atoms.
* @param cAtom
* Index of the carbon atom to be evaluated.
* @param nAtom
* Index of the nitrogen atom to be evaluated.
* @return
* Boolean value denoting if any bond order or atom charge
* modifications were applied.
*/
private static boolean resolveCyanoGroup(Molecule3D mol, int cAtom, int nAtom) {
if (!looksLikeCyanoGroup(mol, cAtom, nAtom)) return false;
// cyano group (R-C#N)
mol.setBondOrder(mol.getBond(cAtom, nAtom), 3);
return true;
}
/**
* Evaluates if the given atoms comprise an isocyanate group (R-N=C=O)
* and assigns appropriate bond orders and atom charges if needed.
*
* @param mol
* Molecule that contains the atoms.
* @param cAtom
* Index of the carbon atom to be evaluated.
* @param nAtom
* Index of the nitrogen atom to be evaluated.
* @param oAtom
* Index of the oxygen atom to be evaluated.
* @return
* Boolean value denoting if any bond order or atom charge
* modifications were applied.
*/
private static boolean resolveIsocyanateGroup(Molecule3D mol, int cAtom, int nAtom, int oAtom) {
if (mol.getAtomicNo(nAtom) != 7) return false;
if (mol.getAtomicNo(cAtom) != 6) return false;
if (calculateHybridizationState(mol, cAtom) != 1) return false;
if (mol.getAtomicNo(oAtom) != 8) return false;
if (mol.getAllConnAtoms(oAtom) != 1) return false;
if (mol.getAttachedHydrogenCount(oAtom) != 0) return false;
// isocyanate group (N=C=O)
mol.setBondOrder(mol.getBond(nAtom, cAtom), 2);
mol.setBondOrder(mol.getBond(cAtom, oAtom), 2);
return true;
}
/**
* Evaluates if the given atoms comprise an isothiocyanate group (R-N=C=S)
* and assigns appropriate bond orders and atom charges if needed.
*
* @param mol
* Molecule that contains the atoms.
* @param cAtom
* Index of the carbon atom to be evaluated.
* @param nAtom
* Index of the nitrogen atom to be evaluated.
* @param sAtom
* Index of the sulphur atom to be evaluated.
* @return
* Boolean value denoting if any bond order or atom charge
* modifications were applied.
*/
private static boolean resolveIsothiocyanateGroup(Molecule3D mol, int cAtom, int nAtom, int sAtom) {
if (mol.getAtomicNo(nAtom) != 7) return false;
if (mol.getAtomicNo(cAtom) != 6) return false;
if (calculateHybridizationState(mol, cAtom) != 1) return false;
if (mol.getAtomicNo(sAtom) != 16) return false;
if (mol.getAllConnAtoms(sAtom) != 1) return false;
if (mol.getAttachedHydrogenCount(sAtom) != 0) return false;
// isothiocyanate group (N=C=S)
mol.setBondOrder(mol.getBond(nAtom, cAtom), 2);
mol.setBondOrder(mol.getBond(cAtom, sAtom), 2);
return true;
}
/**
* Evaluates if the given atoms comprise a thiocyanate ion (S=C=[N-])
* and assigns appropriate bond orders and atom charges if needed.
*
* @param mol
* Molecule that contains the atoms.
* @param cAtom
* Index of the carbon atom to be evaluated.
* @param sAtom
* Index of the sulphur atom to be evaluated.
* @param nAtom
* Index of the nitrogen atom to be evaluated.
* @return
* Boolean value denoting if any bond order or atom charge
* modifications were applied.
*/
private static boolean resolverThiocyanateIon(Molecule3D mol, int cAtom, int sAtom, int nAtom) {
if (mol.getAtomicNo(cAtom) != 6) return false;
if (calculateHybridizationState(mol, cAtom) != 1) return false;
if (mol.getAtomicNo(nAtom) != 7) return false;
if (mol.getAtomicNo(sAtom) != 16) return false;
// TODO: consider the resonance structure of SCN ([S-]-C=N <=> S=C=N[-])?
// thiocyanate ion (S=C=N[-])
mol.setBondOrder(mol.getBond(sAtom, cAtom), 2);
mol.setBondOrder(mol.getBond(cAtom, nAtom), 2);
mol.setAtomCharge(nAtom, -1);
return true;
}
/**
* Evaluate if the given atom could potentially act as the central
* atom of a azide group (R[N]=[N+]=[N-]).
*
* @param mol
* Molecule that contains the atoms.
* @param nCentralAtom
* Index of the nitrogen atom to be evaluated.
* @return
* Boolean value denoting if the given atom could potentially
* act as the central atom of an azide group.
*/
private static boolean looksLikeUnresolvedAzideGroup(Molecule3D mol, int nCentralAtom) {
if (mol.getAtomicNo(nCentralAtom) != 7) return false;
if (calculateHybridizationState(mol, nCentralAtom) != 1) return false;
if (mol.getConnAtoms(nCentralAtom) != 2) return false;
if (getAllConnAtoms(mol, nCentralAtom, true) != 2) return false;
int nAtom1 = mol.getConnAtom(nCentralAtom, 0);
int nAtom2 = mol.getConnAtom(nCentralAtom, 1);
if (mol.getAtomicNo(nAtom1) != 7) return false;
if (getAllConnAtoms(mol, nAtom1, true) > 2) return false;
if (mol.getAtomicNo(nAtom2) != 7) return false;
if (getAllConnAtoms(mol, nAtom2, true) > 2) return false;
// check if there is at least one terminal atom
if (getAllConnAtoms(mol, nAtom1, true) != 1 &&
getAllConnAtoms(mol, nAtom2, true) != 1) {
return false;
}
return true;
}
/**
* Evaluate if the given atom could potentially act as the central
* atom of an azide group (R[N]=[N+]=[N-]) and assigns appropriate bond
* orders and atom charges for the entire group if needed.
*
* @param mol
* Molecule that contains the atoms.
* @param nCentralAtom
* Index of the nitrogen atom to be evaluated.
* @return
* Boolean value denoting if any bond order or atom charge
* modifications were applied.
*/
private static boolean resolveAzideGroup(Molecule3D mol, int nCentralAtom) {
if (!looksLikeUnresolvedAzideGroup(mol, nCentralAtom)) return false;
int nAtom1 = mol.getConnAtom(nCentralAtom, 0);
int nAtom2 = mol.getConnAtom(nCentralAtom, 1);
int shorterBondAtom;
int longerBondAtom;
if (getBondLength(mol, nCentralAtom, nAtom1) < getBondLength(mol, nCentralAtom, nAtom2)) {
shorterBondAtom = nAtom1;
longerBondAtom = nAtom2;
} else {
shorterBondAtom = nAtom2;
longerBondAtom = nAtom1;
}
double shorterBondLength = getBondLength(mol, nCentralAtom, shorterBondAtom);
double longerBondLength = getBondLength(mol, nCentralAtom, longerBondAtom);
boolean accommodatesTripleBond = calculateHybridizationState(mol, shorterBondAtom) == 1 &&
StrictMath.abs(shorterBondAtom - longerBondAtom) > 0.05;
if (accommodatesTripleBond) {
double NNSingleBondLength = getIdealisedBondLength(2, 7, 7, 1, 2);
double NNDoubleBondLength = getIdealisedBondLength(2, 7, 7, 1, 2);
double NNTripleBondLength = getIdealisedBondLength(3, 7, 7, 2, 2);
double N2N2N = StrictMath.abs(NNDoubleBondLength - shorterBondLength) +
StrictMath.abs(NNDoubleBondLength - longerBondLength);
double N1N3N = StrictMath.abs(NNTripleBondLength - shorterBondLength) +
StrictMath.abs(NNSingleBondLength - longerBondLength);
accommodatesTripleBond = N1N3N < N2N2N;
}
// choose between the resonance structures
// RN=[N+]=[N-] and R[N-]-[N+]#N
if (accommodatesTripleBond) {
mol.setBondOrder(mol.getBond(nCentralAtom, shorterBondAtom), 3);
mol.setBondOrder(mol.getBond(nCentralAtom, longerBondAtom), 1);
mol.setAtomCharge(nCentralAtom, +1);
// in case of NH2-[N+]#N
// TODO: check if this conditional is needed
if (mol.getFreeValence(longerBondAtom) > 0) {
mol.setAtomCharge(longerBondAtom, -1);
}
} else {
mol.setBondOrder(mol.getBond(nCentralAtom, shorterBondAtom), 2);
mol.setBondOrder(mol.getBond(nCentralAtom, longerBondAtom), 2);
mol.setAtomCharge(nCentralAtom, +1);
if (mol.getFreeValence(shorterBondAtom) > 0 || calculateHybridizationState(mol, shorterBondAtom) == 1) {
mol.setAtomCharge(shorterBondAtom, -1);
}
if (mol.getFreeValence(longerBondAtom) > 0 || calculateHybridizationState(mol, longerBondAtom) == 1) {
mol.setAtomCharge(longerBondAtom, -1);
}
}
mol.ensureHelperArrays(Molecule.cHelperRings);
return true;
}
/**
* Evaluate if the given atom could potentially act as the central
* atom of a azide ion ([N-]=[N+]=[N-]).
*
* @param mol
* Molecule that contains the atoms.
* @param nCentralAtom
* Index of the nitrogen atom to be evaluated.
* @return
* Boolean value denoting if the given atom could potentially
* act as the central atom of an azide ion.
*/
private static boolean looksLikeUnresolvedAzideIon(Molecule3D mol, int nCentralAtom) {
if (!looksLikeUnresolvedAzideGroup(mol, nCentralAtom)) return false;
int nAtom1 = mol.getConnAtom(nCentralAtom, 0);
int nAtom2 = mol.getConnAtom(nCentralAtom, 1);
if (getAllConnAtoms(mol, nAtom1, true) != 1) return false;
if (getAllConnAtoms(mol, nAtom2, true) != 1) return false;
return true;
}
/**
* Evaluate if the given atom could potentially act as the central
* atom of an azide ion ([N-]=[N+]=[N-]) and assigns appropriate bond
* orders and atom charges for the entire ion if needed.
*
* @param mol
* Molecule that contains the atoms.
* @param nCentralAtom
* Index of the nitrogen atom to be evaluated.
* @return
* Boolean value denoting if any bond order or atom charge
* modifications were applied.
*/
private static boolean resolveAzideIon(Molecule3D mol, int nCentralAtom) {
if (!looksLikeUnresolvedAzideIon(mol, nCentralAtom)) return false;
int a1 = mol.getConnAtom(nCentralAtom, 0);
int a2 = mol.getConnAtom(nCentralAtom, 1);
mol.setBondOrder(mol.getBond(nCentralAtom, a1), 2);
mol.setBondOrder(mol.getBond(nCentralAtom, a2), 2);
mol.setAtomCharge(nCentralAtom, +1);
mol.setAtomCharge(a1, -1);
mol.setAtomCharge(a2, -1);
mol.ensureHelperArrays(Molecule.cHelperRings);
return true;
}
/**
* Evaluates if the given atom is part of a diatomic pnictogen ion of
* the forms [E-2][E-2], [E-1]=[E-1] where E = P, As (i.e. [P-2][P-2],
* [P-1]=[P-1]) and assigns appropriate bond orders and atom charges
* for the entire ion if needed.
*
* @param mol
* Molecule that contains the atoms.
* @param atom
* Index of the atom to be evaluated.
* @return
* Boolean value denoting if any atom charge modifications
* were applied.
*/
private static boolean resolveDiatomicPnictogenIons(Molecule3D mol, int atom) {
if (!mol.isNitrogenFamily(atom)) return false;
// skip nitrogen (N) for now
if (mol.getAtomicNo(atom) == 7) return false;
if (getAllConnAtoms(mol, atom, true) != 1) return false;
if (mol.getConnAtoms(atom) != 1) return false;
int nAtom = mol.getConnAtom(atom, 0);
if (!mol.isNitrogenFamily(nAtom)) return false;
// skip nitrogen (N) for now
if (mol.getAtomicNo(nAtom) == 7) return false;
if (mol.getConnAtoms(nAtom) != 1) return false;
// double bonded form ([E-]=[E-])
if (mol.getAtomicNo(atom) == 15 && mol.getAtomicNo(nAtom) == 15) {
// TODO: eventually replace the hardcoded value with one based on statistics
// Explanation: the provided cutoff value was chosen based
// on bond distances from 6 crystal structures that were
// show to contain the P=P ligand. The hardcoded value may
// eventually be removed if sufficiently granular and accurate
// statistics become available.
double phosphorusDoubleBondCuttoff = 2.05;
if (phosphorusDoubleBondCuttoff > getBondLength(mol, atom, nAtom)) {
mol.setBondOrder(mol.getBond(atom, nAtom), 2);
mol.setAtomCharge(atom, -1);
mol.setAtomCharge(nAtom, -1);
return true;
}
}
// TODO: investigate if a cutoff value for As is also needed
// single bonded form ([E-2][E-2])
mol.setAtomCharge(atom, -2);
mol.setAtomCharge(nAtom, -2);
return true;
}
/**
* Evaluates if the given atom is part of a diatomic chalcogen ion of
* the form [E-][E-] where E = S, Se, Te (i.e. [S-][S-]) and assigns
* appropriate atom charges for the entire ion if needed.
*
* @param mol
* Molecule that contains the atoms.
* @param atom
* Index of the atom to be evaluated.
* @return
* Boolean value denoting if any atom charge modifications
* were applied.
*/
private static boolean resolveDiatomicChalcogenIons(Molecule3D mol, int atom) {
if (!mol.isChalcogene(atom)) return false;
// skip oxygen (O) for now
if (mol.getAtomicNo(atom) == 8) return false;
if (getAllConnAtoms(mol, atom, true) != 1) return false;
if (mol.getConnAtoms(atom) != 1) return false;
int nAtom = mol.getConnAtom(atom, 0);
if (!mol.isChalcogene(nAtom)) return false;
// skip oxygen (O) for now
if (mol.getAtomicNo(nAtom) == 8) return false;
if (mol.getConnAtoms(nAtom) != 1) return false;
mol.setAtomCharge(atom, -1);
mol.setAtomCharge(nAtom, -1);
return true;
}
static void coordinationBondsToCovalent(Molecule3D mol) {
for (int i = 0; i < mol.getRingSet().getSize(); i++ ) {
if ( isRingCoordinated(mol, i) ) {
int ringAtom = mol.getRingSet().getRingAtoms(i)[0];
int metalAtom = -1;
for ( int j = mol.getAllConnAtoms(ringAtom); j < mol.getAllConnAtomsPlusMetalBonds(ringAtom); j++) {
if ( mol.isMetalAtom( mol.getConnAtom(ringAtom, j) ) ) {
metalAtom = mol.getConnAtom(ringAtom, j);
break;
}
}
for (int j = 0; j < mol.getRingSet().getRingSize(i); j++) {
mol.markBondForDeletion(mol.getBond(mol.getRingSet().getRingAtoms(i)[j], metalAtom));
}
}
}
mol.deleteMarkedAtomsAndBonds();
mol.ensureHelperArrays(Molecule3D.cHelperNeighbours);
for (int i = 0; i < mol.getAllBonds(); i++ ) {
if ( mol.getBondType(i) == Molecule3D.cBondTypeMetalLigand ) {
mol.setBondOrder(i, 1);
int a1 = mol.getBondAtom(0, i);
int a2 = mol.getBondAtom(1, i);
if ( mol.isMetalAtom(a1) && mol.isMetalAtom(a2) ) continue;
if ( mol.isMetalAtom(a1) || mol.getAtomicNo(a1) == 1 ) {
mol.setAtomCharge( a1, mol.getAtomCharge(a1) - 1 );
mol.setAtomCharge( a2, mol.getAtomCharge(a2) + 1 );
} else if ( mol.isMetalAtom(a2) || mol.getAtomicNo(a2) == 1 ) {
mol.setAtomCharge( a1, mol.getAtomCharge(a1) + 1 );
mol.setAtomCharge( a2, mol.getAtomCharge(a2) - 1 );
} else {
throw new RuntimeException( "a metal-ligand bond was assigned to a non-metal atoms " +
"'" + mol.getAtomName(a1) + "' and '" + mol.getAtomName(a2) + "'" );
}
}
}
}
private static boolean isPlanar(Molecule3D mol, int a1, int a2) {
Coordinates ci = mol.getCoordinates(a1);
Coordinates u = null, v = null;
for (int i = 0; v == null && i < mol.getAllConnAtoms(a1); i++) {
if ( u == null ) {
u = mol.getCoordinates(mol.getConnAtom(a1, i)).subC(ci);
} else {
v = mol.getCoordinates(mol.getConnAtom(a1, i)).subC(ci);
}
}
for (int i = 0; v == null && i < mol.getAllConnAtoms(a2); i++) {
if ( u == null ) {
u = mol.getCoordinates(mol.getConnAtom(a2, i)).subC(ci);
} else {
v = mol.getCoordinates(mol.getConnAtom(a2, i)).subC(ci);
}
}
if ( u == null ) return false;
Coordinates normal = u.cross(v);
if ( normal.distSq() == 0 ) return false; //what to do?
normal = normal.unitC();
Coordinates cj = mol.getCoordinates(a2);
for(int k = 0; k < mol.getAllConnAtoms(a2); k++) {
Coordinates ck = mol.getCoordinates(mol.getConnAtom(a2, k));
if(Math.abs(ck.subC(cj).dot(normal))>0.25) return false;
}
for(int k=0; k0.25) return false;
}
return true;
}
/**
* Util function for substructure searches. It searches for the substructure
* of form a-b-c where a, b, c are consecutively connected atoms and a != c.
* If a negative value is provided as the atomic number or valency of atom c
* this atom will be ignored and only the a-b substructure search will be
* carried out.
*
* @param mol
* The molecule that contains the atoms.
* @param aIndex
* The index of the atom 'a' in the molecule. The atom 'a' is used
* as the starting position of the substructure search.
* @param bAtomicNo
* The atomic number of the substructure atom 'b' that should be
* connected directly to the atom 'a'.
* @param bValence
* The valence of the substructure atom 'b' that should be connected
* directly to the atom 'a'.
* @param cAtomicNo
* The atomic number of the substructure atom 'c' that should be
* connected directly to the atom 'b'.
* @param cValence
* The valence of the substructure atom 'c' that should be connected
* directly to the atom 'b'.
* @return
* Index of the atom 'b' in the molecule is the search is successful
* and -1 otherwise.
*/
private static int connectedAtom( Molecule3D mol, int aIndex, int bAtomicNo,
int bValence, int cAtomicNo, int cValence ) {
loop: for(int i=0; i 0 && mol.getAtomicNo(atm) != bAtomicNo ) continue;
if ( bValence > 0 && mol.getAllConnAtoms(atm) != bValence ) continue;
// Checking if the atom 'c' that is connected to atom 'b'
// is of the the requested chemical type and valence
if ( cAtomicNo > 0 || cValence > 0 ) {
for ( int j = 0; j < mol.getAllConnAtoms(atm); j++ ) {
int otherAtm = mol.getConnAtom(atm, j);
if(otherAtm==aIndex) continue loop;
if ( cAtomicNo > 0 && mol.getAtomicNo(otherAtm) != cAtomicNo ) continue loop;
if ( cValence > 0 && mol.getAllConnAtoms(otherAtm) != cValence ) continue loop;
}
}
return atm;
}
return -1;
}
private static double getBondLength(Molecule3D mol, int a1, int a2) {
return mol.getCoordinates(a1).distance(mol.getCoordinates(a2));
}
private static double getBondLength(Molecule3D mol, int bond) {
return getBondLength(mol, mol.getBondAtom(0, bond), mol.getBondAtom(1, bond));
}
private static int connectedBond(Molecule3D mol, int a, int atomicNo, int valence) {
if(a<0) return -1;
for(int i=0; i< mol.getAllConnAtoms(a); i++) {
int atm = mol.getConnAtom(a, i);
if(atomicNo>0 && mol.getAtomicNo(atm)!=atomicNo) continue;
if(valence>0 && mol.getAllConnAtoms(atm)!=valence) continue;
return mol.getConnBond(a, i);
}
return -1;
}
/**
* Returns the average bond length value based on bond order
* and the type and pi electron count of the bonded atoms.
* The values are taken from the statistical model
* constructed from previously observed bond parameter
* values.
* @param bondOrder
* Bond order expressed as a positive integer.
* @param atomicNo1
* Atomic number of the first atom.
* @param atomicNo2
* Atomic number of the second atom.
* @param atomPi1
* Pi electron count of the first atom.
* @param atomPi2
* Pi electron count of the second atom.
* @return
* The average bond length.
*/
private static double getIdealisedBondLength(int bondOrder, int atomicNo1, int atomicNo2, int atomPi1, int atomPi2) {
return BondLengthSet.getBondLength(BondLengthSet.getBondIndex(bondOrder, false, false, atomicNo1, atomicNo2, atomPi1, atomPi2));
}
/**
* Produce a logical estimate of the bond being aromatic based
* strictly on length constraints. The bond length is compared
* to the average single bond.
* @param mol
* The molecule that contains the bond.
* @param a1
* Index of the first atom that forms a bond.
* @param a2
* Index of the second atom that forms a bond.
* @return
* Estimate deciding if the bond could be aromatic
* based strictly on the bond length.
*/
private static boolean fitsAromaticLengthConstraints(Molecule3D mol, int a1, int a2) {
int atomType1 = Math.min(mol.getAtomicNo(a1), mol.getAtomicNo(a2));
int atomType2 = Math.max(mol.getAtomicNo(a1), mol.getAtomicNo(a2));
int pi1 = BondLengthSet.isPiConsidered(atomType1) ? 1 : -1;
int pi2 = BondLengthSet.isPiConsidered(atomType2) ? 1 : -1;
// Try to get statistics for the aromatic single bond
boolean aromaticConstraintsUsed = true;
int index = BondLengthSet.getBondIndex(1, true, false, atomType1, atomType2, pi1, pi2);
// If no such bond exist, try to get statistics for the
// regular single bond
if ( index == -1 ) {
aromaticConstraintsUsed = false;
index = BondLengthSet.getBondIndex(1, false, false, atomType1, atomType2, pi1, pi2);
if (index == -1) return false;
}
double bondLength = BondLengthSet.getBondLength(index);
double sigma = BondLengthSet.getBondStdDev(index);
if ( sigma < 0 ) sigma = 0.1;
double distance = getBondLength(mol, a1, a2);
// allow 3 standard deviations for aromatic bonds and 1 for non-aromatic
return ( distance < bondLength + ( aromaticConstraintsUsed ? 3 : 1 ) * sigma);
}
/**
* Produce a numerical estimate of how well the bond length fits the statistical model
* constructed from previously observed bond parameter values. The estimate is provided
* as a deviation value from the idealised value (smaller values mean better fitness).
* Currently it only compares the bond length with the average bond length of the same
* type, but a more intricate model might be implemented in the future.
* @param mol
* The molecule that contains the bond.
* @param a1
* Index of the first atom that forms a bond.
* @param a2
* Index of the second atom that forms a bond.
* @param bondOrder
* Order of the bond. Currently only values '1' (single bond)
* and '2' (double bond) are supported.
* @return
* Bond length deviation from the previously observed bond parameter value.
*/
private static double estimatePiBondFitness(Molecule3D mol, int a1, int a2, int bondOrder) {
int atomType1 = Math.min(mol.getAtomicNo(a1), mol.getAtomicNo(a2));
int atomType2 = Math.max(mol.getAtomicNo(a1), mol.getAtomicNo(a2));
double singleBond = getIdealisedBondLength(1, atomType1, atomType2, 1, 1);
double doubleBond = getIdealisedBondLength(2, atomType1, atomType2, 1, 1);
boolean statisticsAvailable = singleBond != -1 || doubleBond != -1;
if (!statisticsAvailable) {
throw new RuntimeException("bond length statistics for bond '" + atomType1
+ "-" + atomType2 + "' is not available");
}
double idealLength = ( bondOrder == 1 ) ? singleBond : doubleBond;
if ( bondOrder == 2 && getBondLength(mol, a1, a2) < idealLength ) {
return 0;
}
return StrictMath.abs( getBondLength(mol, a1, a2) - idealLength );
}
private static boolean isAlkaliMetalAtom (Molecule3D mol, int atom) {
return ( mol.getAtomicNo(atom) == 3 ) ||
( mol.getAtomicNo(atom) == 11 ) ||
( mol.getAtomicNo(atom) == 19 ) ||
( mol.getAtomicNo(atom) == 37 ) ||
( mol.getAtomicNo(atom) == 55 ) ||
( mol.getAtomicNo(atom) == 87 );
}
private static boolean isAlkalineEarthMetalAtom(Molecule3D mol, int atom) {
return ( mol.getAtomicNo(atom) == 4 ) ||
( mol.getAtomicNo(atom) == 12 ) ||
( mol.getAtomicNo(atom) == 20 ) ||
( mol.getAtomicNo(atom) == 38 ) ||
( mol.getAtomicNo(atom) == 56 ) ||
( mol.getAtomicNo(atom) == 88 );
}
List checkMoleculeForBumps(Molecule3D mol) {
List bumps = new LinkedList<>();
double distance;
for (int i = 0; i < mol.getAllAtoms(); i++) {
for (int j = i+1; j < mol.getAllAtoms(); j++) {
if ( i < j ) {
distance = mol.getCoordinates(i).distance(mol.getCoordinates(j));
if ( distance < SPECIAL_POSITION_CUTOFF ||
distance < BUMP_FACTOR * ( COVALENT_RADIUS[mol.getAtomicNo(i)] +
COVALENT_RADIUS[mol.getAtomicNo(j)] ) ) {
bumps.add( "atoms '" + mol.getAtomName(i) + "' and '" + mol.getAtomName(j) + "' " +
"are too close (distance = " + String.format("%6.4f", distance) + ") and " +
"are considered a bump" );
}
}
}
}
return bumps;
}
private static boolean hasMetalLigandBonds(Molecule3D mol, int atom) {
return (mol.getMetalBondedConnAtoms(atom) > 0);
}
private static boolean isCommonIonCharge(int atomicNo, int charge) {
for ( int i = 0; i < Molecule.cCommonOxidationState[atomicNo].length; i++ ) {
if ( Molecule.cCommonOxidationState[atomicNo][i] == charge ) {
return true;
} else if ( Molecule.cCommonOxidationState[atomicNo][i] > charge ) {
return false;
}
}
return false;
}
private static boolean isCommonIonCharge(Molecule3D mol, int atomicNo, int charge) {
return isCommonIonCharge(mol.getAtomicNo(atomicNo), charge);
}
public static int getMaximumIonCharge(Molecule3D mol, int atom) {
if ( Molecule.cCommonOxidationState[mol.getAtomicNo(atom)] == null ) {
return 0;
} else {
return Molecule.cCommonOxidationState[mol.getAtomicNo(atom)]
[Molecule.cCommonOxidationState[mol.getAtomicNo(atom)].length - 1];
}
}
private static boolean maximumIonChargeReached(Molecule3D mol, int atom) {
return mol.getAtomCharge(atom) >= getMaximumIonCharge(mol, atom);
}
private static int getDelocalizedBondCount(Molecule3D mol, int atom) {
return getDelocalizedBonds(mol, atom).size();
}
private static List getDelocalizedBonds(Molecule3D mol, int atom) {
List delocalizedBonds = new LinkedList<>();
for( int i = 0; i < mol.getAllConnAtoms(atom); i++) {
if ( mol.getBondType(mol.getConnBond(atom, i)) == Molecule3D.cBondTypeDelocalized ) {
delocalizedBonds.add(mol.getConnBond(atom, i));
}
}
return delocalizedBonds;
}
private static List> getPiBondFragments (Molecule3D mol) {
List> fragments = new LinkedList<>();
mol.ensureHelperArrays(Molecule3D.cHelperRings);
for (int i = 0; i < mol.getAllAtoms(); i++) {
if (getDelocalizedBondCount(mol, i) == 1) {
// Check if the selected atom is not already at the end of analyzed fragments
boolean used = false;
for (List fragment1 : fragments) {
if (fragment1.get(fragment1.size() - 1) == i) {
used = true;
break;
}
}
if ( used ) { continue; }
List fragment = new LinkedList<>();
fragment.add(i);
int prevBond = getDelocalizedBonds(mol, i).get(0);
int atom;
if ( mol.getBondAtom( 0, prevBond ) != i ) {
atom = mol.getBondAtom( 0, prevBond );
} else {
atom = mol.getBondAtom( 1, prevBond );
}
boolean finished = false;
while (!finished) {
// Currently only process linear pi-bond systems
// TODO: extend the logic to process branched systems
if (getDelocalizedBondCount(mol, atom) > 2) {
break;
}
fragment.add(atom);
if (getDelocalizedBondCount(mol, atom) == 1) {
fragments.add(fragment);
finished = true;
} else {
int currentBond = ( getDelocalizedBonds(mol, atom).get(0) != prevBond ) ?
getDelocalizedBonds(mol, atom).get(0) :
getDelocalizedBonds(mol, atom).get(1);
if ( mol.getBondAtom( 0, currentBond ) != atom ) {
atom = mol.getBondAtom( 0, currentBond );
} else {
atom = mol.getBondAtom( 1, currentBond );
}
prevBond = currentBond;
}
}
}
}
return fragments;
}
private static boolean fragmentCanResonate(Molecule3D mol, List fragment) {
for (Integer atom : fragment) {
if (mol.getFreeValence(atom) - 1 < 0) {
return false;
}
}
return true;
}
private static int[] sortAtomsByBondLength(Molecule3D mol, int a, List bondAtoms) {
double[] bondLength = new double[bondAtoms.size()];
for (int i = 0; i < bondLength.length; i++ ) {
bondLength[i] = getBondLength(mol, a, mol.getConnAtom(a, bondAtoms.get(i)));
}
// using primitive bubble sort for now to sort the indexes
int[] indexes = new int[bondAtoms.size()];
for (int i = 0; i < bondLength.length; i++) {
indexes[i] = i;
}
for (int i = 0; i < bondLength.length; i++) {
for (int j = i+1; j < bondLength.length; j++ ) {
if( bondLength[indexes[i]] > bondLength[indexes[j]] ) {
int tmp = indexes[i];
indexes[i] = indexes[j];
indexes[j] = tmp;
}
}
}
// replacing sort indexes with atom indexes
for (int i = 0; i < bondLength.length; i++) {
indexes[i] = bondAtoms.get(indexes[i]);
}
return indexes;
}
private static Coordinates calculatePlaneNormal (Coordinates p1, Coordinates p2, Coordinates p3) {
Coordinates u = p1.subC(p2);
Coordinates v = p1.subC(p3);
return u.cross(v);
}
private static double calculatePointToPlaneDistance ( Coordinates normal, Coordinates p0, Coordinates p1 ) {
Coordinates v1 = p0.subC(p1);
return StrictMath.abs(normal.unitC().dot(v1) / v1.dist());
}
private static void setTerminalDoubleBonds( Molecule3D mol, int a, int[] sortedAtoms,
int maxDoubleBondCount ) {
int doubleBondCount = 0;
for (int sortedAtom : sortedAtoms) {
if (doubleBondCount < maxDoubleBondCount) {
mol.setBondOrder(mol.getConnBond(a, sortedAtom), 2);
doubleBondCount++;
} else {
mol.setBondOrder(mol.getConnBond(a, sortedAtom), 1);
mol.setAtomCharge(mol.getConnAtom(a, sortedAtom), -1);
}
}
}
// TODO: ask to add to the molecule
private static boolean isRingAtom(Molecule3D mol, int ring, int a) {
int[] ringAtoms = mol.getRingSet().getRingAtoms(ring);
for (int ringAtom : ringAtoms) {
if (ringAtom == a) return true;
}
return false;
}
private static void identifySquarates(Molecule3D mol, int[] spOrder) {
for (int i = 0; i < mol.getAtoms(); i++) {
if (mol.getAtomRingSize(i) != 4) continue;
if (getRingCount(mol, i) != 1) continue;
int ring = getSmallestRing(mol, i);
int[] ringAtoms = mol.getRingSet().getRingAtoms(ring);
boolean[] allowsExocyclicDoubleBond = new boolean[ringAtoms.length];
int[] exocyclicAtoms = new int[ringAtoms.length];
boolean fitsConstraints = true;
for (int j = 0; j < ringAtoms.length; j++) {
int a = ringAtoms[j];
if (mol.getAtomicNo(a) == 6 && mol.getConnAtoms(a) == 3 &&
spOrder[a] == 2 && mol.getLowestFreeValence(a) > 0) {
int a1 = ringAtoms[(j - 1 + ringAtoms.length) % ringAtoms.length];
int a2 = ringAtoms[(j + 1 + ringAtoms.length) % ringAtoms.length];
if (StrictMath.abs(90 - StrictMath.toDegrees(
getAngle(mol, a1, a, a2))) > 5) {
fitsConstraints = false;
break;
}
for (int k = 0; k < mol.getConnAtoms(a); k++) {
if (!isRingAtom(mol, ring, mol.getConnAtom(a, k))) {
allowsExocyclicDoubleBond[j] = (mol.getLowestFreeValence(mol.getConnAtom(a, k)) > 0);
exocyclicAtoms[j] = mol.getConnAtom(a, k);
break;
}
}
} else {
fitsConstraints = false;
break;
}
}
if (fitsConstraints) {
for (int a1 = 0; a1 < ringAtoms.length; a1++) {
int a2 = (a1 + 1) % ringAtoms.length;
int a3 = (a1 + 2) % ringAtoms.length;
int a4 = (a1 + 3) % ringAtoms.length;
if (allowsExocyclicDoubleBond[a1] && allowsExocyclicDoubleBond[a2]) {
mol.setBondOrder(mol.getBond(ringAtoms[a1], exocyclicAtoms[a1]), 2);
mol.setBondOrder(mol.getBond(ringAtoms[a2], exocyclicAtoms[a2]), 2);
mol.setBondOrder(mol.getBond(ringAtoms[a3], ringAtoms[a4]), 2);
if (allowsExocyclicDoubleBond[a3]) mol.setAtomCharge(exocyclicAtoms[a3], -1);
if (allowsExocyclicDoubleBond[a4]) mol.setAtomCharge(exocyclicAtoms[a4], -1);
break;
}
}
}
}
}
private static int getRingCount(Molecule3D mol, int a) {
int ringCount = 0;
RingCollection rings = mol.getRingSet();
for (int i = 0; i < mol.getRingSet().getSize(); i++) {
int[] ringAtoms = rings.getRingAtoms(i);
for (int ringAtom : ringAtoms) {
if (ringAtom == a) {
ringCount++;
break;
}
}
}
return ringCount;
}
private static int getSmallestRing(Molecule3D mol, int a) {
int smallestRing = -1;
int ringSize = -1;
RingCollection rings = mol.getRingSet();
for (int i = 0; i < mol.getRingSet().getSize(); i++) {
int[] ringAtoms = rings.getRingAtoms(i);
for (int ringAtom : ringAtoms) {
if (ringAtom == a) {
if (ringSize == -1 || ringAtoms.length < ringSize) {
ringSize = ringAtoms.length;
smallestRing = i;
break;
}
}
}
}
return smallestRing;
}
private static int getAllConnAtoms(Molecule3D mol, int a, boolean includeAttachedHydrogens) {
if ( includeAttachedHydrogens ) {
return mol.getAllConnAtoms(a) + mol.getAttachedHydrogenCount(a);
} else {
return mol.getAllConnAtoms(a);
}
}
private static boolean isRingCoordinated(Molecule3D mol, int ringNo) {
int[] ringAtoms = mol.getRingSet().getRingAtoms(ringNo);
boolean atomsAreCoordinated = true;
for (int ringAtom : ringAtoms) {
atomsAreCoordinated = atomsAreCoordinated && hasMetalLigandBonds(mol, ringAtom);
}
return atomsAreCoordinated;
}
/**
* Determines if the given atom is capable of accommodating an additional double bond.
* @param mol
* Molecule that contains the atom.
* @param atom
* Index of the atom inside the molecule data structure.
* @param atomPi
* The number of pi electrons that the atom is currently using.
* @param spOrder
* The sp order (1, 2, 3) of the atom.
* @return
* Logical value denoting is an atom can accommodate an additional double bond.
*/
private static boolean canAccommodateDoubleBond(Molecule3D mol, int atom, int atomPi, int spOrder) {
// Tetravalent silicon (Si) cannot accommodate any double bonds
if (mol.getAtomicNo(atom) == 14 && mol.getOccupiedValence(atom) > 3) return false;
// Germanium (Ge)
if (mol.getAtomicNo(atom) == 32) {
// only sp3 germanium (Ge) can potentially accommodate a double bond
if (mol.getOccupiedValence(atom) > 3) return false;
if (spOrder != 2) return false;
}
// Phosphorus (P)
if (mol.getAtomicNo(atom) == 15) {
if (atomPi > 0 && spOrder != 1) return false;
}
if ( mol.getLowestFreeValence(atom) < 1 ) return false;
if ( mol.getAtomCharge(atom) != 0 ) return false;
// Sp1-like atoms are not allowed to have more than 2 pi electrons
if ( spOrder == 1 && atomPi > 1 ) return false;
// Sp2-like atoms are not allowed to have more than 1 pi electron
if ( spOrder == 2 && atomPi > 0 ) return false;
// Atoms smaller than phosphorus are not allowed a double bond in Sp3-like configuration
if ( spOrder == 3 && mol.getAtomicNo(atom) < 14 ) return false;
return true;
}
/**
* Calculates the bond length deviation from the statistical
* average for the bond assuming it is of the given bond type.
* @param mol
* Molecule that contains the bond.
* @param bond
* The bond that is being evaluated.
* @param bondTypeIndex
* The bond type index that is used to retrieve
* the average bond length and standard deviation
* for the given bond type.
* @return
* The difference between the actual and the average
* bond lengths.
*/
private static double determineExpectedBondDeviation(Molecule3D mol, int bond, int bondTypeIndex) {
return ( getBondLength(mol, bond ) - BondLengthSet.getBondLength(bondTypeIndex));
}
/**
* Calculates the angle a1-a2-a3 between 3 atoms.
* @param mol
* The molecule that contains the atoms.
* @param a1
* Index of the first atom in the molecule.
* @param a2
* Index of the second atom in the molecule.
* @param a3
* Index of the third atom in the molecule.
* @return
* The angle a1-a2-a3 in radians.
*/
private static double getAngle(Molecule mol, int a1, int a2, int a3) {
Coordinates c1 = mol.getCoordinates(a1);
Coordinates c2 = mol.getCoordinates(a2);
Coordinates c3 = mol.getCoordinates(a3);
return c1.subC(c2).getAngle(c3.subC(c2));
}
/**
* Calculates the dihedral angle a1-a2-a3-a4 between 4 atoms.
* @param mol
* The molecule that contains the atoms.
* @param a1
* Index of the first atom in the molecule.
* @param a2
* Index of the second atom in the molecule.
* @param a3
* Index of the third atom in the molecule.
* @param a4
* Index of the fourth atom in the molecule.
* @return
* The dihedral angle a1-a2-a3-a4 in radians.
*/
private static double getDihedral(Molecule mol, int a1, int a2, int a3, int a4) {
Coordinates c1 = mol.getCoordinates(a1);
Coordinates c2 = mol.getCoordinates(a2);
Coordinates c3 = mol.getCoordinates(a3);
Coordinates c4 = mol.getCoordinates(a4);
return c1.getDihedral(c2, c3, c4);
}
public static int connected(StereoMolecule mol, int a, int atomicNo, int bondOrder) {
for(int i=0; i=0 && mol.getAtomicNo(atm)!=atomicNo) continue;
if(bondOrder>0 && mol.getConnBondOrder(a, i)!=bondOrder) continue;
return atm;
}
return -1;
}
public static boolean aromatize(Molecule3D mol, Set aromaticAtoms, Set aromaticBonds) {
ArrayList[] atomToRings = getAtomToRings(mol);
RingCollection ringSet = mol.getRingSetSimple();
return aromatize(mol,atomToRings,ringSet,aromaticAtoms,aromaticBonds);
}
public static boolean aromatize(Molecule3D mol, ArrayList[] atomToRings, RingCollection ringSet, Set aromaticAtoms, Set aromaticBonds) {
//RingCollection ringSet = mol.getRingSetSimple();
//Flag the aromatic rings
boolean[] aromaticRings = new boolean[ringSet.getSize()];
for (int i = 0; i < ringSet.getSize(); i++) {
//Is is an aromatic ring
boolean isAromatic = true;
int ringSize = ringSet.getRingSize(i);
for (int j = -1; j < ringSize; j++) {
int[] ringAtom = ringSet.getRingAtoms(i);
int a1 = j==-1? ringAtom[ringSize-1]: ringAtom[j];
int a2 = j==ringSize-1? ringAtom[0]: ringAtom[j+1];
int b = mol.getBond(a1, a2);
if(!aromaticBonds.contains(b)) {
isAromatic = false;
}
}
aromaticRings[i] = isAromatic;
}
Set nonAromaticAtoms = new HashSet();
for (int i=0;i(), true);
if(!success) success = aromatize(mol, atomToRings, ringSet, i, aromaticRings, nonAromaticAtoms, new boolean[mol.getAllAtoms()], 0, ringSet.getRingSize(i)%2, new ArrayList(), false);
if(!success) {
System.out.println("Could not aromatize ring "+i);
aromaticRings[i] = false;
ok = false;
}
}
}
return ok;
}
private static boolean aromatize(Molecule3D mol, ArrayList[] atomToRings, RingCollection ringSet, int index, boolean[] aromatic, Set nonAromaticAtoms, boolean[] visited, int seed, int left, List bondsMade, boolean easy) {
//RingCollection ringSet = mol.getRingSetSimple();
//Ends if the ring has been fully visited
int[] ring = (int[]) ringSet.getRingAtoms(index);
boolean allVisited = true;
int bestSeed = -1;
for (int i = 0; i < ring.length; i++) {
if(!visited[ring[(seed + i)%ring.length]]) {
if(bestSeed<0) bestSeed = seed + i;
allVisited = false;
}
}
if(allVisited) {
return true;
} else {
seed = bestSeed;
}
int a = ring[seed%ring.length];
int ap = ring[(seed+1)%ring.length];
if(visited[a]) { //already treated
return aromatize(mol, atomToRings, ringSet, index, aromatic, nonAromaticAtoms, visited, seed+1, left, bondsMade, easy);
} else {
//Try to create a double bond from the atom a to a connected one
for (int j = -1; j < mol.getAllConnAtoms(a); j++) {
int a2 = j==-1? ap: mol.getConnAtom(a, j);
if(visited[a2]) continue;
if(j>=0 && a2==ap) continue;
if(nonAromaticAtoms.contains(a2)) continue;
if(nonAromaticAtoms.contains(a)) continue;
if(mol.getAtomicNo(a)==8 || mol.getAtomicNo(a)==16) continue;
if(mol.getAtomicNo(a2)==8 || mol.getAtomicNo(a2)==16) continue;
if(easy && mol.getFreeValence(a) <= 0) continue;
if(easy && mol.getFreeValence(a2) <= 0) continue;
if(connected(mol, a, -1, 2)>=0) continue;
if(connected(mol, a2, -1, 2)>=0) continue;
visited[a] = visited[a2] = true;
int b = mol.getBond(a, a2);
mol.setBondOrder(b, 2);
//Test whole ring
List trackBondsMade = new ArrayList();
boolean success = aromatize(mol, atomToRings, ringSet, index, aromatic, nonAromaticAtoms, visited, seed+1, left, trackBondsMade, easy);
//Test connecting rings
if(success) {
List rings = atomToRings[a2];
if(rings.size()>1) {
for (int r : rings) {
if(r!=index && r>=0 && r " +success+" "+trackBondsMade.size());
if(!success) break;
}
}
}
}
if(success) {
//It works!!!
bondsMade.add(b);
bondsMade.addAll(trackBondsMade);
return true;
} else {
//Backtrack changes
visited[a] = visited[a2] = false;
mol.setBondOrder(b, 1);
for (int b2 : trackBondsMade) {
//System.out.println("retrack "+mol.getBondAtom(0, b2)+"-"+mol.getBondAtom(1, b2));
mol.setBondOrder(b2, 1);
visited[mol.getBondAtom(0, b2)] = visited[mol.getBondAtom(1, b2)] = false;
}
}
}
//Try to skip this atom
if(left>0 && (mol.getAtomicNo(a)!=6)) {
visited[a] = true;
List trackBondsMade = new ArrayList();
boolean success = aromatize(mol, atomToRings, ringSet, index, aromatic, nonAromaticAtoms, visited, seed+1, left-1, trackBondsMade, easy);
if(success) {
bondsMade.addAll(trackBondsMade);
return true;
} else {
visited[a] = false;
for (int b2 : trackBondsMade) {
mol.setBondOrder(b2, 1);
visited[mol.getBondAtom(0, b2)] = visited[mol.getBondAtom(1, b2)] = false;
}
}
}
return false;
}
}
} © 2015 - 2025 Weber Informatics LLC | Privacy Policy