All Downloads are FREE. Search and download functionalities are using the official Maven repository.

com.actelion.research.chem.CoordinateInventor Maven / Gradle / Ivy

There is a newer version: 2024.12.1
Show newest version
/*
* Copyright (c) 1997 - 2016
* Actelion Pharmaceuticals Ltd.
* Gewerbestrasse 16
* CH-4123 Allschwil, Switzerland
*
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice, this
*    list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
*    this list of conditions and the following disclaimer in the documentation
*    and/or other materials provided with the distribution.
* 3. Neither the name of the the copyright holder nor the
*    names of its contributors may be used to endorse or promote products
*    derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*/

package com.actelion.research.chem;

import java.util.*;

public class CoordinateInventor {
	public static final int MODE_REMOVE_HYDROGEN = 1;
	public static final int MODE_KEEP_MARKED_ATOM_COORDS = 2;
	public static final int MODE_PREFER_MARKED_ATOM_COORDS = 4;
	private static final int MODE_CONSIDER_MARKED_ATOMS = MODE_KEEP_MARKED_ATOM_COORDS | MODE_PREFER_MARKED_ATOM_COORDS;

	private static final byte FLIP_AS_LAST_RESORT = 1;
	private static final byte FLIP_POSSIBLE = 2;
	private static final byte FLIP_PREFERRED = 3;
	private static final int  PREFERRED_FLIPS = 32;
	private static final int  POSSIBLE_FLIPS = 64;
	private static final int  LAST_RESORT_FLIPS = 128;
	private static final int  TOTAL_FLIPS = PREFERRED_FLIPS + POSSIBLE_FLIPS + LAST_RESORT_FLIPS;

	private StereoMolecule mMol;
	private ArrayList   mFragmentList;
	private Random		mRandom;
	private boolean[]	mAtomHandled;
	private boolean[]	mBondHandled;
	private int			mMode,mAtoms,mBonds;
	private int[]		mConnAtoms;

	/**
	 * Creates an CoordinateInventor, which removes unneeded hydrogen atoms
	 * and creates new atom coordinates for all(!) atoms.
	 */
	public CoordinateInventor () {
		mMode = MODE_REMOVE_HYDROGEN;
		}


	/**
	 * Creates an CoordinateInventor, which removes unneeded hydrogens, if mode flags include
	 * MODE_REMOVE_HYDROGEN. If mode includes MODE_KEEP_MARKED_ATOM_COORDS, then marked atoms
	 * keep their coordinates. If mode includes MODE_PREFER_MARKED_ATOM_COORDS, then coordinates
	 * of marked atoms are changed only, if perfect coordinates are not possible without.
	 * @param mode
	 */
	public CoordinateInventor (int mode) {
		mMode = mode;
		}


	public void setRandomSeed(long seed) {
		mRandom = new Random(seed);
		}


	/**
	 * Creates new atom 2D-coordinates for a molecule or a part of a molecule.
	 * Coordinates will correctly reflect E/Z double bond parities, unless the double bond is in a ring.
	 * If atom parities are available, this call is typically followed by a calling mol.setStereoBondsFromParity();
	 * Unneeded explicit hydrogens are removed, if mode includes MODE_REMOVE_HYDROGEN.
	 * The relative orientation of all marked atoms is retained, if mode includes MODE_KEEP_MARKED_ATOM_COORDS.
	 * The relative orientation of all marked atoms is changed as last resort only, if mode includes MODE_PREFER_MARKED_ATOM_COORDS.
	 * @param mol the molecule that gets new 2D coordinates in place
	 */
	public void invent(StereoMolecule mol) {
		if (mRandom == null)
			mRandom = new Random();

		mMol = mol;
		mMol.ensureHelperArrays(Molecule.cHelperRings);

		if ((mMode & MODE_REMOVE_HYDROGEN) == 0) {
			mAtoms = mMol.getAllAtoms();
			mBonds = mMol.getAllBonds();
			mConnAtoms = new int[mAtoms];
			for (int atom=0; atom();
		mAtomHandled = new boolean[mAtoms];
		mBondHandled = new boolean[mBonds];

		if ((mMode & MODE_CONSIDER_MARKED_ATOMS) != 0)
			locateCoreFragment();

		locateInitialFragments();
		joinOverlappingFragments();

		locateChainFragments();
		joinOverlappingFragments();

		locateFragmentBonds();
		correctChainEZParities();
		optimizeFragments();

		locateSingleAtoms();

		arrangeAllFragments();

		for (int i=0; i 4) {
				InventorFragment f = new InventorFragment(mMol, 1+mConnAtoms[atom]);

				f.mAtomX[mConnAtoms[atom]] = 0.0;
				f.mAtomY[mConnAtoms[atom]] = 0.0;
				f.mPriority[mConnAtoms[atom]] = 32;
				f.mAtom[mConnAtoms[atom]] = atom;
				mAtomHandled[atom] = true;

				for (int i=0; i 2) {
					InventorFragment f = new InventorFragment(mMol, members);
					int count = 0;
					for (int i=0; i 2)
								break;

							mAtomHandled[alleneAtom[last+1]] = true;
							mBondHandled[mMol.getConnBond(alleneAtom[last], nextIndex)] = true;
							last++;
							} while (mMol.getAtomPi(alleneAtom[last]) == 2
								  && mConnAtoms[alleneAtom[last]] == 2);

						int members = mConnAtoms[alleneAtom[0]]
									+ mConnAtoms[alleneAtom[last]]
									+ last - 1;
						InventorFragment f = new InventorFragment(mMol, members);
						for (int j=0; j<=last; j++) {
							f.mAtomX[j] = (double)j;
							f.mAtomY[j] = 0.0;
							f.mPriority[j] = 64;
							f.mAtom[j] = alleneAtom[j];
							}

						int current = last+1;
						boolean found = false;
						for (int j=0; j connAtom)
								secondAtomCounts[i] = true;
							f.mAtomX[count] = (i==0) ? -1.0 : 1.0;
							f.mAtomY[count] = (neighbours==0) ? -Math.sin(Math.PI/3) : Math.sin(Math.PI/3);
							f.mAtom[count++] = connAtom;
							mAtomHandled[connAtom] = true;
							mBondHandled[mMol.getConnBond(bondAtom[i], j)] = true;
							neighbours++;
							}
						}
					}

				if ((bondParity == Molecule.cBondParityE) ^ (secondAtomCounts[0] ^ secondAtomCounts[1]))
					for (int i=1; i longestChain.getChainLength())
						longestChain = theChain;
					}
				}

			if (longestChain == null)
				break;

			InventorFragment f = new InventorFragment(mMol, longestChain.getChainLength());
			for (int i=0; i= 10 && f.size() <= 24 && cBondZList[ringIndex] != null) {
			int maxBit = (1 << f.size());
			int bondEConstraint = 0;
			int bondZConstraint = 0;
			for (int i=0; i>>= 1;
				bondZConstraint >>>= 1;
				}

			for (int zList=0; zList>>= 1;
								}
							return;
							}
						if ((bondZList & 1) != 0)
							bondZList |= maxBit;
						bondZList >>>= 1;
						}
					}
				}
			}

			// if not successful so far
		createRegularRingFragment(f);
		}


	private InventorChain getSmallestRingFromBond(int bond) {
			// find smallest ring of given bond
		int atom1 = mMol.getBondAtom(0, bond);
		int atom2 = mMol.getBondAtom(1, bond);
		int graphAtom[] = new int[mAtoms];
		int graphBond[] = new int[mAtoms];
		int graphLevel[] = new int[mAtoms];
		int graphParent[] = new int[mAtoms];
		graphAtom[0] = atom1;
		graphAtom[1] = atom2;
		graphBond[1] = bond;
		graphLevel[atom1] = 1;
		graphLevel[atom2] = 2;
		graphParent[0] = -1;
		graphParent[1] = 0;
		int current = 1;
		int highest = 1;
		while (current <= highest) {
//			if (graphLevel[graphAtom[current]] > RingCollection.MAX_LARGE_RING_SIZE)
//				return null;		// disabled ring size limit;  TLS 20130613
			for (int i=0; i 1) && candidate == atom1) {
					InventorChain theRing = new InventorChain(graphLevel[graphAtom[current]]);
					graphBond[0] = mMol.getConnBond(graphAtom[current], i);
					int index = current;
					for (int j=0; j RingCollection.MAX_LARGE_RING_SIZE)
//				return 0;		// disabled ring size limit;  TLS 20130613
			for (int i=0; i 0) {
						int handlePreferred = (commonAtoms == 1
											&& getConnAtoms(f1, commonAtom) == 1
											&& getConnAtoms(f2, commonAtom) == 1) ? 0 : 1;

						int joinPriority;
						if (maxF1Priority > maxF2Priority)
							joinPriority = (handlePreferred << 24)
										 + (maxF1Priority << 16)
										 + (maxF2Priority << 8)
										 + commonAtoms;
						else
							joinPriority = (handlePreferred << 24)
										 + (maxF2Priority << 16)
										 + (maxF1Priority << 8)
										 + commonAtoms;

						if (maxJoinPriority < joinPriority) {
							maxJoinPriority = joinPriority;
							maxCommonAtoms = commonAtoms;

							// retain coordinates of fragment with highest priority atom
							maxF1Priority = 0;
							maxF2Priority = 0;
							for (int k=0; k maxF2Priority) {
								maxFragment1 = f1;
								maxFragment2 = f2;
								}
							else {
								maxFragment1 = f2;
								maxFragment2 = f1;
								}
							}
						}
					}
				}

			if (maxJoinPriority == 0)
				break;

			if (maxCommonAtoms == maxFragment1.size())
				mFragmentList.remove(maxFragment1);
			else if (maxCommonAtoms == maxFragment2.size())
				mFragmentList.remove(maxFragment2);
			else
				joinFragments(maxFragment1, maxFragment2, maxCommonAtoms);
			}
		}


	private void joinFragments(InventorFragment f1, InventorFragment f2, int commonAtoms) {
		int[] commonAtom = new int[commonAtoms];
		int count = 0;
		for (int i=0; i Math.abs(getAngleDif(meanNeighbourAngleF1.mAngle, meanNeighbourAngleF2Flip.mAngle))) {
			f2.rotate(meanX1, meanY1, meanAngleDif.mAngle);
			}
		else {
			f2.flip(meanX1, meanY1, 0.0);
			f2.rotate(meanX1, meanY1, meanAngleDifFlip.mAngle);
			}

		return getMergedFragment(f1, f2, commonAtoms);
		}


	private InventorFragment getMergedFragment(InventorFragment f1, InventorFragment f2, int commonAtoms) {
			// merges all atoms of two fragments into a new one retaining original coordinates
		InventorFragment f = new InventorFragment(mMol, f1.mAtom.length + f2.mAtom.length - commonAtoms);
		int count = 0;
		for (int i=0; i0; i--) {	// bubble sort
			for (int j=0; j connAngle[j+1]) {
					double tempAngle = connAngle[j];
					connAngle[j] = connAngle[j+1];
					connAngle[j+1] = tempAngle;
					int tempAtom = connAtom[j];
					connAtom[j] = connAtom[j+1];
					connAtom[j+1] = tempAtom;
					int tempBond = connBond[j];
					connBond[j] = connBond[j+1];
					connBond[j+1] = tempBond;
					}
				}
			}

		connAngle[connAngles] = connAngle[0] + 2 * Math.PI;
		connAtom[connAngles] = connAtom[0];
		connBond[connAngles] = connBond[0];

		double maxAngleDif = -100.0;
		int maxIndex = 0;
		for (int i=0; i 2
			 && mMol.isRingBond(connBond[i])
			 && mMol.isRingBond(connBond[i+1])) {
				int ringSize = getSmallestRingSize(connAtom[i], atom, connAtom[i+1]);
				if (ringSize != 0)
					angleDif -= 100.0 - ringSize;
				}

			if (maxAngleDif < angleDif) {
				maxAngleDif = angleDif;
				maxIndex = i;
				}
			}

		return (connAngle[maxIndex] + connAngle[maxIndex+1]) / 2;
		}


	private double getAngleDif(double angle1, double angle2) {
		double angleDif = angle1 - angle2;
		while (angleDif < -Math.PI)
			angleDif += 2 * Math.PI;
		while (angleDif > Math.PI)
			angleDif -= 2 * Math.PI;
		return angleDif;
		}


	protected static InventorAngle getMeanAngle(InventorAngle[] angle, int noOfAngles) {
		// adds noOfAngles vectors of length=1 with angles angle[i]
		// and returns angle of sum-vector
		// length of sum-vector is criteria for deviation
		double sinSum = 0;
		double cosSum = 0;
		for (int i=0; i 0) ? Math.PI/2 : -Math.PI/2;
		else {
			meanAngle = Math.atan(sinSum/cosSum);
			if (cosSum < 0)
				meanAngle += Math.PI;
			}

		double length = Math.sqrt(sinSum * sinSum + cosSum * cosSum) / noOfAngles;

		return new InventorAngle(meanAngle, length);
		}


	private int getConnAtoms(InventorFragment f, int atom) {
		int connAtoms = 0;
		for (int i=0; i 1)
					 && (mMol.getConnAtoms(mMol.getBondAtom(1, bond)) > 1)
					 && (mMol.getBondParity(bond) == Molecule.cBondParityEor1
					  || mMol.getBondParity(bond) == Molecule.cBondParityZor2)) {
						int[] minConnAtom = new int[2];
						int[] bondAtom = new int[2];
						for (int j=0; j<2; j++) {
							minConnAtom[j] = mMol.getMaxAtoms();
							bondAtom[j] = mMol.getBondAtom(j, bond);
							for (int k=0; k connAtom)
									minConnAtom[j] = connAtom;
								}
							}
	
						double dbAngle = InventorAngle.getAngle(f.mAtomX[f.mAtomIndex[bondAtom[0]]],
																f.mAtomY[f.mAtomIndex[bondAtom[0]]],
																f.mAtomX[f.mAtomIndex[bondAtom[1]]],
																f.mAtomY[f.mAtomIndex[bondAtom[1]]]);
						double angle1  = InventorAngle.getAngle(f.mAtomX[f.mAtomIndex[minConnAtom[0]]],
																f.mAtomY[f.mAtomIndex[minConnAtom[0]]],
																f.mAtomX[f.mAtomIndex[bondAtom[0]]],
																f.mAtomY[f.mAtomIndex[bondAtom[0]]]);
						double angle2  = InventorAngle.getAngle(f.mAtomX[f.mAtomIndex[bondAtom[1]]],
																f.mAtomY[f.mAtomIndex[bondAtom[1]]],
																f.mAtomX[f.mAtomIndex[minConnAtom[1]]],
																f.mAtomY[f.mAtomIndex[minConnAtom[1]]]);
	
						if (((getAngleDif(dbAngle, angle1) < 0)
						   ^ (getAngleDif(dbAngle, angle2) < 0))
						  ^ (mMol.getBondParity(bond) == Molecule.cBondParityZor2)) {
							f.flipOneSide(bond);
							}
						}
					}
				}
			}
		}


	private void optimizeFragments() {
		int[] atomSymRank = calculateAtomSymmetries();

		byte[] bondFlipPriority = new byte[mBonds];
		locateFlipBonds(bondFlipPriority, atomSymRank);
		for (int bond=0; bond collisionList = f.getCollisionList();
			double minCollisionPanalty = f.getCollisionPanalty();
			InventorFragment minCollisionFragment = new InventorFragment(f);

			int lastBond = -1;
			for (int flip=0; flip= FLIP_POSSIBLE)
							availableBond[availableBonds++] = bondSequence[i];
					}
				else {
					for (int i=1; i= FLIP_AS_LAST_RESORT)
							availableBond[availableBonds++] = bondSequence[i];
					}

				if (availableBonds != 0) {
					int theBond = availableBond[0];
					if (availableBonds > 1) {
						do {	// don't rotate 2 times around same bond
							theBond = availableBond[mRandom.nextInt(availableBonds)];
							} while (theBond == lastBond);
						}

					if (theBond != lastBond) {
						lastBond = theBond;
		
						f.flipOneSide(theBond);

							// getCollisionList() is necessary to update collision panalty
						collisionList = f.getCollisionList();
						if (minCollisionPanalty > f.getCollisionPanalty()) {
							minCollisionPanalty = f.getCollisionPanalty();
							minCollisionFragment = new InventorFragment(f);
							}
						}
					}
				}

			mFragmentList.set(fragmentNo, minCollisionFragment);
				// finished optimization by rotating around single bonds

				// starting optimization by moving individual atoms
/*
double avbl = mMol.getAverageBondLength();
for (int i=0; i currentRank && theRank < nextAvailableRank)
						nextAvailableRank = theRank;
					}
				currentRank = nextAvailableRank;
				} while (nextAvailableRank != 9999);
			}
		}


	private int[] getShortestConnection(int atom1, int atom2) {
		int graphAtom[] = new int[mAtoms];
		int graphBond[] = new int[mAtoms];
		int graphLevel[] = new int[mAtoms];
		int graphParent[] = new int[mAtoms];
		graphAtom[0] = atom2;
		graphLevel[atom2] = 1;
		graphParent[0] = -1;
		int current = 0;
		int highest = 0;
		while (current <= highest) {
			for (int i=0; i 2) {
					boolean symmetricEndFound = true;
					int connSymRank = -1;
					for (int j=0; jj; k--)
					connRank[k] = connRank[k-1];
				connRank[j] = rank;
				}

			int neighbours = Math.min(6, mConnAtoms[atom]);
			baseValue[atom].init(atom);
			baseValue[atom].add(Canonizer.ATOM_BITS, symRank[atom]);
			baseValue[atom].add((6 - neighbours)*(Canonizer.ATOM_BITS + 1), 0);
			for (int i=0; i 1) {
			double[] largeFragmentSize = new double[2];
			InventorFragment[] largeFragment = new InventorFragment[2];

				// find largest two Fragments (in terms of display size)
			InventorFragment f0 = mFragmentList.get(0);
			InventorFragment f1 = mFragmentList.get(1);
			double s0 = f0.getWidth() + f0.getHeight();
			double s1 = f1.getWidth() + f1.getHeight();
			if (s0 > s1) {
				largeFragment[0] = f0;
				largeFragmentSize[0] = s0;
				largeFragment[1] = f1;
				largeFragmentSize[1] = s1;
				}
			else {
				largeFragment[0] = f1;
				largeFragmentSize[0] = s1;
				largeFragment[1] = f0;
				largeFragmentSize[1] = s0;
				}

			for (int i=2; ithe bond atom that lies on the larger side of the bond
				// [1]->the bond atom on the smaller side of the bond
				// [2...n]->all other atoms on the smaller side of the bond.
				//		  These are the ones getting flipped on the mirror
				//		  line defined by the bond.
			if (mFlipList == null)
				mFlipList = new int[mBonds][];

			if (mFlipList[bond] == null) {
				int[] graphAtom = new int[mAtom.length];
				boolean[] isOnSide = new boolean[mAtoms];
				int atom1 = mMol.getBondAtom(0, bond);
				int atom2 = mMol.getBondAtom(1, bond);
				graphAtom[0] = atom1;
				isOnSide[atom1] = true;
				int current = 0;
				int highest = 0;
				while (current <= highest) {
					for (int i=0; i mAtom.length/2);

				// if we retain core atoms and the smaller side contains core atoms, then flip the larger side
				if ((mMode & MODE_CONSIDER_MARKED_ATOMS) != 0) {
					boolean coreOnSide = false;
					boolean coreOffSide = false;
					for (int i=0; i=2) ? corner-2 : corner+2);
				if (maxGain < gain) {
					maxGain = gain;
					maxCorner = corner;
					}
				}

			double sumHeight = getHeight() + f.getHeight();
			double sumWidth = 0.75 * (getWidth() + f.getWidth());
			double maxHeight = Math.max(getHeight(), f.getHeight());
			double maxWidth = 0.75 * Math.max(getWidth(), f.getWidth());

			double bestCornerSize = Math.sqrt((sumHeight - maxGain) * (sumHeight - maxGain)
											+ (sumWidth - 0.75 * maxGain) * (sumWidth - 0.75 * maxGain));
			double toppedSize = Math.max(maxWidth, sumHeight);
			double besideSize = Math.max(maxHeight, sumWidth);

			if (bestCornerSize < toppedSize && bestCornerSize < besideSize) {
				switch(maxCorner) {
				case 0:
					f.translate(mMaxX - f.mMinX - maxGain + 1.0, mMinY - f.mMaxY + maxGain - 1.0);
					break;
				case 1:
					f.translate(mMaxX - f.mMinX - maxGain + 1.0, mMaxY - f.mMinY - maxGain + 1.0);
					break;
				case 2:
					f.translate(mMinX - f.mMaxX + maxGain - 1.0, mMaxY - f.mMinY - maxGain + 1.0);
					break;
				case 3:
					f.translate(mMinX - f.mMaxX + maxGain - 1.0, mMinY - f.mMaxY + maxGain - 1.0);
					break;
					}
				}
			else if (besideSize < toppedSize) {
				f.translate(mMaxX - f.mMinX + 1.0, (mMaxY + mMinY - f.mMaxY - f.mMinY) / 2);
				}
			else {
				f.translate((mMaxX + mMinX - f.mMaxX - f.mMinX) / 2, mMaxY - f.mMinY + 1.0);
				return;
				}
			}

		private void calculateMinMax() {
			if (mMinMaxAvail)
				return;

			mMinX = mAtomX[0];
			mMaxX = mAtomX[0];
			mMinY = mAtomY[0];
			mMaxY = mAtomY[0];
			for (int i=0; i mAtomX[i] - surplus)
					mMinX = mAtomX[i] - surplus;
				if (mMaxX < mAtomX[i] + surplus)
					mMaxX = mAtomX[i] + surplus;
				if (mMinY > mAtomY[i] - surplus)
					mMinY = mAtomY[i] - surplus;
				if (mMaxY < mAtomY[i] + surplus)
					mMaxY = mAtomY[i] + surplus;
				}

			mMinMaxAvail = true;
			}

		private double getCornerDistance(int corner) {
			double minDistance = 9999.0;
			for (int atom=0; atom d - surplus)
					minDistance = d - surplus;
				}

			return minDistance;
			}

		private double getAtomSurplus(int atom) {
			return (mMol.getAtomQueryFeatures(mAtom[atom]) != 0) ? 0.6
						  : (mMol.getAtomicNo(mAtom[atom]) != 6) ? 0.25 : 0.0;
			}

		protected void generateAtomListsOfFlipBonds(byte[] flipBondType) {
			
			}

		protected ArrayList getCollisionList() {
			mCollisionPanalty = 0.0;
			ArrayList collisionList = new ArrayList();
			for (int i=1; i mAtom[i])
						fragmentBonds++;

			mBond = new int[fragmentBonds];
			mAtomIndex = new int[mAtoms];

			fragmentBonds = 0;
			for (int i=0; i mAtom[i]) {
						mBond[fragmentBonds] = mMol.getConnBond(mAtom[i], j);
						fragmentBonds++;
						}
					}
				}
			}

		protected void optimizeAtomCoordinates(int atomIndex) {
			double x = mAtomX[atomIndex];
			double y = mAtomY[atomIndex];

			InventorAngle[] collisionForce = new InventorAngle[4];

			int forces = 0;
			for (int i=0; i= 4)
					break;

				if (atomIndex == mAtomIndex[mMol.getBondAtom(0, mBond[i])]
				 || atomIndex == mAtomIndex[mMol.getBondAtom(1, mBond[i])])
					continue;

				double x1 = mAtomX[mAtomIndex[mMol.getBondAtom(0, mBond[i])]];
				double y1 = mAtomY[mAtomIndex[mMol.getBondAtom(0, mBond[i])]];
				double x2 = mAtomX[mAtomIndex[mMol.getBondAtom(1, mBond[i])]];
				double y2 = mAtomY[mAtomIndex[mMol.getBondAtom(1, mBond[i])]];
				double d1 = Math.sqrt((x1-x)*(x1-x)+(y1-y)*(y1-y));
				double d2 = Math.sqrt((x2-x)*(x2-x)+(y2-y)*(y2-y));
				double bondLength = Math.sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));

				if (d1 0) {
				InventorAngle force = CoordinateInventor.getMeanAngle(collisionForce, forces);
				mAtomX[atomIndex] += force.mLength * Math.sin(force.mAngle);
				mAtomY[atomIndex] += force.mLength * Math.cos(force.mAngle);
				}
			}
		}
	}

class InventorAngle {
	double mAngle;
	double mLength;

	protected static double getAngle(double x1, double y1, double x2, double y2) {
		double angle;
		double xdif = x2 - x1;
		double ydif = y2 - y1;

		if (ydif != 0) {
			angle = Math.atan(xdif/ydif);
			if (ydif < 0) {
				if (xdif < 0)
					angle -= Math.PI;
				else
					angle += Math.PI;
				}
			}
		else
			angle = (xdif >0) ? Math.PI/2 : -Math.PI/2;

		return angle;
		}

	protected InventorAngle(double angle, double length) {
		mAngle = angle;
		mLength = length;
		}

	protected InventorAngle(double x1, double y1, double x2, double y2) {
		mAngle = getAngle(x1, y1, x2, y2);
		double xdif = x2 - x1;
		double ydif = y2 - y1;
		mLength = Math.sqrt(xdif * xdif + ydif * ydif);
		}
	}


class InventorChain {
	protected int[] mAtom;
	protected int[] mBond;

	public InventorChain(int chainLength) {
		mAtom = new int[chainLength];
		mBond = new int[chainLength];
			// last mBond[] value isn't needed if chain isn't a ring
		}

	protected int getChainLength() {
		return mAtom.length;
		}

	protected int[] getRingAtoms() {
		return mAtom;
		}

	protected int[] getRingBonds() {
		return mBond;
		}
	}


//class FlipBondIterator




© 2015 - 2025 Weber Informatics LLC | Privacy Policy