[Git][java-team/openchemlib][master] 7 commits: New upstream version 2021.10.1+dfsg

Andrius Merkys (@merkys) gitlab at salsa.debian.org
Tue Oct 26 07:20:57 BST 2021



Andrius Merkys pushed to branch master at Debian Java Maintainers / openchemlib


Commits:
5c2ee3a2 by Andrius Merkys at 2021-10-26T01:43:00-04:00
New upstream version 2021.10.1+dfsg
- - - - -
b4239435 by Andrius Merkys at 2021-10-26T01:43:06-04:00
Update upstream source from tag 'upstream/2021.10.1+dfsg'

Update to upstream version '2021.10.1+dfsg'
with Debian dir 7581eca60d78616c904ceb27449c0df09c3e0712
- - - - -
642ec862 by Andrius Merkys at 2021-10-26T01:43:45-04:00
Removing empty debian/maven.cleanIgnoreRules and debian/maven.publishedRules.

- - - - -
91b5e386 by Andrius Merkys at 2021-10-26T01:44:06-04:00
Removing empty debian/maven.properties.

- - - - -
7bd2286b by Andrius Merkys at 2021-10-26T01:44:33-04:00
Removing boilerplate debian/README.source.

- - - - -
dcac7ad2 by Andrius Merkys at 2021-10-26T01:47:26-04:00
Describing patches.

- - - - -
0856cc21 by Andrius Merkys at 2021-10-26T01:48:00-04:00
Update changelog for 2021.10.1+dfsg-1 release

- - - - -


25 changed files:

- − debian/README.source
- debian/changelog
- − debian/maven.cleanIgnoreRules
- − debian/maven.properties
- − debian/maven.publishedRules
- debian/patches/javafx.patch
- debian/patches/libtablelayout-java.patch
- pom.xml
- src/main/java/com/actelion/research/chem/AromaticityResolver.java
- src/main/java/com/actelion/research/chem/Molecule.java
- src/main/java/com/actelion/research/chem/SmilesParser.java
- src/main/java/com/actelion/research/chem/alignment3d/KabschAlignment.java
- src/main/java/com/actelion/research/chem/docking/DockingEngine.java
- src/main/java/com/actelion/research/chem/docking/LigandPose.java
- src/main/java/com/actelion/research/chem/docking/receptorpharmacophore/NegativeReceptorImage.java
- src/main/java/com/actelion/research/chem/docking/scoring/AbstractScoringEngine.java
- + src/main/java/com/actelion/research/chem/forcefield/mmff/MMFFExternalPositionConstraint.java
- src/main/java/com/actelion/research/chem/forcefield/mmff/PositionConstraint.java → src/main/java/com/actelion/research/chem/forcefield/mmff/MMFFPositionConstraint.java
- src/main/java/com/actelion/research/chem/io/CompoundFileHelper.java
- src/main/java/com/actelion/research/chem/mcs/MCS.java
- − src/main/java/com/actelion/research/chem/mcs/SubStructSearchBondIndex.java
- src/main/java/com/actelion/research/chem/mcs/SubStructSearchExhaustiveTreeWalker.java
- src/main/java/com/actelion/research/chem/phesa/pharmacophore/PharmacophoreCalculator.java
- src/main/java/com/actelion/research/chem/phesaflex/FlexibleShapeAlignment.java
- src/main/java/com/actelion/research/chem/potentialenergy/PositionConstraint.java


Changes:

=====================================
debian/README.source deleted
=====================================
@@ -1,9 +0,0 @@
-Information about openchemlib
------------------------------
-
-This package was debianized using the mh_make command
-from the maven-debian-helper package.
-
-The build system uses Maven but prevents it from downloading
-anything from the Internet, making the build compliant with
-the Debian policy.


=====================================
debian/changelog
=====================================
@@ -1,3 +1,13 @@
+openchemlib (2021.10.1+dfsg-1) unstable; urgency=medium
+
+  * New upstream version 2021.10.1+dfsg
+  * Removing empty debian/maven.cleanIgnoreRules, debian/maven.publishedRules
+    and debian/maven.properties.
+  * Removing boilerplate debian/README.source.
+  * Describing patches.
+
+ -- Andrius Merkys <merkys at debian.org>  Tue, 26 Oct 2021 01:47:35 -0400
+
 openchemlib (2021.10.0+dfsg-1) unstable; urgency=medium
 
   * New upstream version 2021.10.0+dfsg


=====================================
debian/maven.cleanIgnoreRules deleted
=====================================
@@ -1 +0,0 @@
-


=====================================
debian/maven.properties deleted
=====================================
@@ -1,5 +0,0 @@
-# Include here properties to pass to Maven during the build.
-# For example:
-# maven.test.skip=true
-# project.build.sourceEncoding=UTF-8
-


=====================================
debian/maven.publishedRules deleted
=====================================
@@ -1 +0,0 @@
-


=====================================
debian/patches/javafx.patch
=====================================
@@ -1,3 +1,8 @@
+Description: Adding dependencies required to build the package in Debian..
+ Forwarded initially upstream (https://github.com/Actelion/openchemlib/pull/47),
+ but retracted due to patch being Debian-specific.
+Author: Andrius Merkys <merkys at debian.org>
+Forwarded: not-needed
 --- a/pom.xml
 +++ b/pom.xml
 @@ -107,6 +107,21 @@


=====================================
debian/patches/libtablelayout-java.patch
=====================================
@@ -1,3 +1,7 @@
+Description: Using org.debian.tablelayout.TableLayout instead of
+ info.clearthought.layout.TableLayout.
+Author: Andrius Merkys <merkys at debian.org>
+Forwarded: not-needed
 --- a/pom.xml
 +++ b/pom.xml
 @@ -122,6 +122,10 @@


=====================================
pom.xml
=====================================
@@ -8,7 +8,7 @@
     Please follow the naming scheme YEAR.MONTH.RELEASE_NO_OF_MONTH
     (eg. 2016.4.1 for second release in Apr 2016)
     -->
-    <version>2021.10.0</version>
+    <version>2021.10.1</version>
 
     <name>OpenChemLib</name>
     <description>Open Source Chemistry Library</description>
@@ -195,7 +195,7 @@
         <connection>scm:git:git at github.com:Actelion/openchemlib.git</connection>
         <developerConnection>scm:git:git at github.com:Actelion/openchemlib.git</developerConnection>
         <url>https://github.com/Actelion/openchemlib</url>
-      <tag>openchemlib-2021.10.0</tag>
+      <tag>openchemlib-2021.10.1</tag>
   </scm>
 
     <distributionManagement>


=====================================
src/main/java/com/actelion/research/chem/AromaticityResolver.java
=====================================
@@ -332,8 +332,9 @@ public class AromaticityResolver {
 			if (mIsDelocalizedAtom[atom]
 			 && mMol.getLowestFreeValence(atom) == 0
 			 && (!mayChangeAtomCharges
-			  || !mMol.isElectronegative(atom)
-			  || mMol.getAtomCharge(atom) > 0))
+			  || (mMol.getAtomicNo(atom) == 5 && mMol.getAtomCharge(atom) < 0)
+			  || (mMol.getAtomicNo(atom) == 6 || mMol.getAtomicNo(atom) == 14)
+			  || (mMol.isElectronegative(atom) && mMol.getAtomCharge(atom) > 0)))
 				protectAtom(atom);
 		}
 
@@ -679,39 +680,23 @@ public class AromaticityResolver {
 	 */
 	private boolean checkAtomTypePi1(int atom, boolean correctCharge) {
 		int atomicNo = mMol.getAtomicNo(atom);
-		if ((atomicNo >=5 && atomicNo <= 8)
+		if ((atomicNo >= 5 && atomicNo <= 8)
 		 || atomicNo == 15 || atomicNo == 16 || atomicNo == 33 || atomicNo == 34 || atomicNo == 52) {	// P,S,As,Se,Te
 
-// Old logic seems fishy to me; TLS 10Dec2020
-//			int freeValence = mMol.getFreeValence(atom);
-//			if (freeValence == 1 || freeValence == 2)	// we allow one more free valence, because the atom may have a missing charge
-//  			return true;
-
 			int freeValence = mMol.getLowestFreeValence(atom);
 			if (freeValence != 0)
 				return true;
 
-			if (mMol.getAtomCharge(atom) == 0) {
-				if ((atomicNo == 15 || atomicNo == 33) /* && freeValence == 3 */) {
-					if (correctCharge)
-						mMol.setAtomCharge(atom, 1);
-					return true;
-					}
-				if ((atomicNo == 16 || atomicNo == 34 || atomicNo == 52) /* && freeValence == 4 */) {
-					if (correctCharge)
-						mMol.setAtomCharge(atom, 1);
-					return true;
-					}
-				if (atomicNo == 5 /* && freeValence == 0 */) {
-					if (correctCharge)
-						mMol.setAtomCharge(atom, -1);
-					return true;
-					}
-				if ((atomicNo == 7 || atomicNo == 8) /* && freeValence == 0 */) {
-					if (correctCharge)
-						mMol.setAtomCharge(atom, 1);
-					return true;
-					}
+			int charge = mMol.getAtomCharge(atom);
+			if (atomicNo == 5 && charge >= 0) {
+				if (correctCharge)
+					mMol.setAtomCharge(atom, charge-1);
+				return true;
+				}
+			if (atomicNo != 5 && charge <= 0) {
+				if (correctCharge)
+					mMol.setAtomCharge(atom, charge+1);
+				return true;
 				}
 			}
 


=====================================
src/main/java/com/actelion/research/chem/Molecule.java
=====================================
@@ -370,7 +370,7 @@ public class Molecule implements Serializable {
 	public static final int cDefaultAtomValence = 6;
 	private static final byte[] cDefaultAtomValences = { cDefaultAtomValence };
 	private static final byte[] cAminoAcidValences = { 2 };
-	private static final byte[][] cAtomValence = {null,
+	public static final byte[][] cAtomValence = {null,
 			{1}, {0}, {1}, {2}, {3}, {4}, {3}, {2}, {1}, {0},			// H to Ne
 			{1}, {2}, {3}, {4}, {3, 5}, {2, 4, 6}, {1, 3, 5, 7}, {0},	// Na to Ar
 			{1}, {2}, null, null, null, null, null, null, null, null,	// K to Ni


=====================================
src/main/java/com/actelion/research/chem/SmilesParser.java
=====================================
@@ -50,6 +50,7 @@ public class SmilesParser {
 	public static final int SMARTS_MODE_IS_SMARTS = 2;
 
 	public static final int MODE_SKIP_COORDINATE_TEMPLATES = 4;
+	public static final int MODE_MAKE_HYDROGEN_EXPLICIT = 8;
 
 	private static final int INITIAL_CONNECTIONS = 16;
 	private static final int MAX_CONNECTIONS = 100; // largest allowed one in SMILES is 99
@@ -65,7 +66,7 @@ public class SmilesParser {
 	private StereoMolecule mMol;
 	private boolean[] mIsAromaticBond;
 	private int mAromaticAtoms,mAromaticBonds,mSmartsMode,mCoordinateMode;
-	private boolean mCreateSmartsWarnings,mSkipTemplates;
+	private boolean mCreateSmartsWarnings, mMakeHydrogenExplicit;
 	private StringBuilder mSmartsWarningBuffer;
 
 	/**
@@ -91,11 +92,12 @@ public class SmilesParser {
 	public SmilesParser(int mode, boolean createSmartsWarnings) {
 		mSmartsMode = mode & SMARTS_MODE_MASK;
 		mCreateSmartsWarnings = createSmartsWarnings;
+		mMakeHydrogenExplicit = ((mode & MODE_MAKE_HYDROGEN_EXPLICIT) != 0);
 		mCoordinateMode = CoordinateInventor.MODE_DEFAULT;
 		if ((mode & MODE_SKIP_COORDINATE_TEMPLATES) != 0)
 			mCoordinateMode |= CoordinateInventor.MODE_SKIP_DEFAULT_TEMPLATES;
-
-		mSkipTemplates = ((mode & MODE_SKIP_COORDINATE_TEMPLATES) != 0);
+		if (mMakeHydrogenExplicit)
+			mCoordinateMode &= ~CoordinateInventor.MODE_REMOVE_HYDROGEN;
 		}
 
 	public StereoMolecule parseMolecule(String smiles) {
@@ -624,7 +626,7 @@ public class SmilesParser {
 
 				// mark aromatic atoms
 				if (Character.isLowerCase(theChar)) {
-					if (atomicNo != 5 && atomicNo != 6 && atomicNo != 7 && atomicNo != 8 && atomicNo != 15 &&atomicNo != 16)
+					if (atomicNo != 5 && atomicNo != 6 && atomicNo != 7 && atomicNo != 8 && atomicNo != 15 &&atomicNo != 16 && atomicNo != 33  && atomicNo != 34)
 						throw new Exception("SmilesParser: atomicNo "+atomicNo+" must not be aromatic");
 
 					mMol.setAtomMarker(atom, true);
@@ -637,7 +639,7 @@ public class SmilesParser {
 				// put explicitHydrogen into atomCustomLabel to keep atom-relation when hydrogens move to end of atom list in handleHydrogen()
 				if (explicitHydrogens != HYDROGEN_ANY && atomicNo != 1) {	// no custom labels for hydrogen to get useful results in getHandleHydrogenMap()
 					byte[] bytes = new byte[1];
-					bytes[0] = (byte)explicitHydrogens;
+					bytes[0] = (byte)(explicitHydrogens == HYDROGEN_IMPLICIT_ZERO ? 0 : explicitHydrogens);
 					mMol.setAtomCustomLabel(atom, bytes);
 					}
 
@@ -908,7 +910,11 @@ public class SmilesParser {
 			if (mMol.getAtomCustomLabel(atom) != null) {	// if we have the exact number of hydrogens
 				int explicitHydrogen = mMol.getAtomCustomLabelBytes(atom)[0];
 
-				if (smartsFeatureFound || mSmartsMode == SMARTS_MODE_IS_SMARTS) {
+				if (mMakeHydrogenExplicit) {
+					for (int i=0; i<explicitHydrogen; i++)
+						mMol.addBond(atom, mMol.addAtom(1), 1);
+					}
+				else if (smartsFeatureFound || mSmartsMode == SMARTS_MODE_IS_SMARTS) {
 					if (explicitHydrogen == 0)
 						mMol.setAtomQueryFeature(atom, Molecule.cAtomQFHydrogen & ~Molecule.cAtomQFNot0Hydrogen, true);
 					if (explicitHydrogen == 1)
@@ -919,9 +925,6 @@ public class SmilesParser {
 						mMol.setAtomQueryFeature(atom, Molecule.cAtomQFHydrogen & ~Molecule.cAtomQFNot3Hydrogen, true);
 					}
 				else {
-					if (explicitHydrogen == HYDROGEN_IMPLICIT_ZERO)
-						explicitHydrogen = 0;
-
 					if (!mMol.isMarkedAtom(atom)) {
 						// We don't correct aromatic atoms, because for aromatic atoms the number of
 						// explicit hydrogens encodes whether a pi-bond needs to be placed at the atom
@@ -954,8 +957,24 @@ public class SmilesParser {
 						}
 					}
 				}
+			else if (!mMakeHydrogenExplicit && (smartsFeatureFound || mSmartsMode == SMARTS_MODE_IS_SMARTS)) {
+				// if we don't have a hydrogen count on the atom, but we have explicit hydrogen atoms
+				// and if we decode a SMARTS, then we convert explicit hydrogens into an 'at least n hydrogen'
+				int explicitHydrogen = mMol.getExplicitHydrogens(atom);
+				if (explicitHydrogen >= 1)
+					mMol.setAtomQueryFeature(atom, Molecule.cAtomQFNot0Hydrogen, true);
+				if (explicitHydrogen >= 2)
+					mMol.setAtomQueryFeature(atom, Molecule.cAtomQFNot1Hydrogen, true);
+				if (explicitHydrogen >= 3)
+					mMol.setAtomQueryFeature(atom, Molecule.cAtomQFNot2Hydrogen, true);
+				if (explicitHydrogen >= 4)
+					mMol.setAtomQueryFeature(atom, Molecule.cAtomQFNot3Hydrogen, true);
+				}
 			}
 
+		if (!mMakeHydrogenExplicit && (smartsFeatureFound || mSmartsMode == SMARTS_MODE_IS_SMARTS))
+			mMol.removeExplicitHydrogens();
+
 		mMol.ensureHelperArrays(Molecule.cHelperNeighbours);
 
 		correctValenceExceededNitrogen();	// convert pyridine oxides and nitro into polar structures with valid nitrogen valences
@@ -1343,9 +1362,8 @@ public class SmilesParser {
 		 || !mMol.isMarkedAtom(atom))	// already marked as hetero-atom of another ring
 			return false;
 
-		int explicitHydrogens = (mMol.getAtomCustomLabel(atom) == null
-							  || mMol.getAtomCustomLabelBytes(atom)[0] == HYDROGEN_IMPLICIT_ZERO) ?
-											0 : mMol.getAtomCustomLabelBytes(atom)[0];
+		int explicitHydrogens = (mMol.getAtomCustomLabel(atom) == null) ?
+								0 : mMol.getAtomCustomLabelBytes(atom)[0];
 		int freeValence = mMol.getFreeValence(atom) - explicitHydrogens;
 		if (freeValence < 1)
 			return false;


=====================================
src/main/java/com/actelion/research/chem/alignment3d/KabschAlignment.java
=====================================
@@ -66,9 +66,14 @@ public class KabschAlignment {
 		for(Coordinates c2 : coords2)
 			c2.sub(com2);
 		
-		trans1 = com2.scaleC(-1.0);
-		trans2 = new Coordinates(com1);
-
+		Coordinates t1 = com2.scaleC(-1.0);
+		trans1.x = t1.x;
+		trans1.y = t1.y;
+		trans1.z = t1.z;
+		
+		trans2.x = com1.x;
+		trans2.y = com1.y;
+		trans2.z = com1.z;
 		
 		Matrix m = new Matrix(3,3);
 		double [][] c1 = Arrays.stream(coords1).map(e -> new double[] {e.x,e.y,e.z}).toArray(double[][]::new);
@@ -97,15 +102,15 @@ public class KabschAlignment {
 		ma.set(2,2,det);
 		
 		
-		rot = ma.multiply(ut);
-		rot = v.multiply(rot);
-		assert(rot.det()>0.0);
-		rot = rot.getTranspose();
+		Matrix r = ma.multiply(ut);
+		r = v.multiply(r);
+		assert(r.det()>0.0);
+		r = r.getTranspose();
 		
 		
 		
 	    for(Coordinates c : coords2) {
-	    	c.rotate(rot.getArray());
+	    	c.rotate(r.getArray());
 	    	c.add(com1);
 	    }
 	    
@@ -113,6 +118,8 @@ public class KabschAlignment {
 	    	c.add(com1);
 	    }
 	    
+	    rot.set(r.getArray());
+	    
 
 		
 	}


=====================================
src/main/java/com/actelion/research/chem/docking/DockingEngine.java
=====================================
@@ -5,15 +5,20 @@ import java.util.ArrayList;
 import java.util.Base64;
 import java.util.HashMap;
 import java.util.HashSet;
+import java.util.Iterator;
+import java.util.LinkedHashMap;
 import java.util.List;
 import java.util.Map;
 import java.util.Random;
 import java.util.Set;
 import java.util.Base64.Decoder;
 import java.util.Base64.Encoder;
+import java.util.Map.Entry;
+import java.util.stream.Collectors;
 
 import org.openmolecules.chem.conf.gen.ConformerGenerator;
 
+import com.actelion.research.calc.Matrix;
 import com.actelion.research.calc.ThreadMaster;
 import com.actelion.research.chem.Canonizer;
 import com.actelion.research.chem.Coordinates;
@@ -21,7 +26,9 @@ import com.actelion.research.chem.IDCodeParser;
 import com.actelion.research.chem.IDCodeParserWithoutCoordinateInvention;
 import com.actelion.research.chem.Molecule;
 import com.actelion.research.chem.Molecule3D;
+import com.actelion.research.chem.SSSearcher;
 import com.actelion.research.chem.StereoMolecule;
+import com.actelion.research.chem.alignment3d.KabschAlignment;
 import com.actelion.research.chem.alignment3d.transformation.ExponentialMap;
 import com.actelion.research.chem.alignment3d.transformation.Quaternion;
 import com.actelion.research.chem.alignment3d.transformation.Rotation;
@@ -35,15 +42,18 @@ import com.actelion.research.chem.docking.scoring.AbstractScoringEngine;
 import com.actelion.research.chem.docking.scoring.ChemPLP;
 import com.actelion.research.chem.docking.scoring.IdoScore;
 import com.actelion.research.chem.docking.shape.ShapeDocking;
+import com.actelion.research.chem.forcefield.mmff.MMFFExternalPositionConstraint;
 import com.actelion.research.chem.forcefield.mmff.ForceFieldMMFF94;
-import com.actelion.research.chem.forcefield.mmff.PositionConstraint;
+import com.actelion.research.chem.forcefield.mmff.MMFFPositionConstraint;
 import com.actelion.research.chem.interactionstatistics.InteractionAtomTypeCalculator;
 import com.actelion.research.chem.io.pdb.converter.MoleculeGrid;
+import com.actelion.research.chem.mcs.MCS;
 import com.actelion.research.chem.optimization.OptimizerLBFGS;
 import com.actelion.research.chem.phesa.EncodeFunctions;
 import com.actelion.research.chem.phesa.MolecularVolume;
 import com.actelion.research.chem.phesa.PheSAAlignment;
 import com.actelion.research.chem.phesa.ShapeVolume;
+import com.actelion.research.chem.potentialenergy.PositionConstraint;
 import com.actelion.research.chem.phesa.PheSAAlignment.PheSAResult;
 
 public class DockingEngine {
@@ -62,6 +72,7 @@ public class DockingEngine {
 	public static final double GRID_DIMENSION = 6.0;
 	public static final double GRID_RESOLUTION = 0.5;
 	public static final double MINI_CUTOFF = 100; //if energy higher after MC step, don't minimize
+	public static final int MCS_EXHAUSTIVENESS = 3; //MCS docking requires more starting positions, but less MC sampling (reduced DOF)
 	
 
 	private Rotation rotation; //for initial prealignment to native ligand
@@ -73,6 +84,9 @@ public class DockingEngine {
 	private StereoMolecule nativeLigand;
 	private ShapeDocking shapeDocking;
 	private ThreadMaster threadMaster;
+	private StereoMolecule mcsRef;
+	private List<Integer> mcsConstrainedBonds;
+	private List<Integer> mcsConstrainedAtoms;
 
 	
 	public DockingEngine(StereoMolecule rec, StereoMolecule nativeLig, int mcSteps, int startPositions,
@@ -122,20 +136,24 @@ public class DockingEngine {
 		threadMaster = tm;
 	}
 	
+	
+	/**
+	 * generate initial poses: 
+	 * 1) shape docking into the negative receptor image
+	 * 2) constrained optimization of initial poses to reduce strain energy
+	 * @param mol
+	 * @param initialPos
+	 * @return
+	 * @throws DockingFailedException
+	 */
 	private double getStartingPositions(StereoMolecule mol, List<Conformer> initialPos) throws DockingFailedException {
 
 		double eMin = Double.MAX_VALUE;
 		Map<String, Object> ffOptions = new HashMap<String, Object>();
 		ffOptions.put("dielectric constant", 80.0);
-		//ffOptions.put("angle bend", false);
-		//ffOptions.put("stretch bend", false);
-		//ffOptions.put("bond stretch", false);
-		//ffOptions.put("out of plane", false);
 
 		ConformerSet confSet = new ConformerSet();
-		long t0 = System.currentTimeMillis();
 		List<StereoMolecule> alignedMol = shapeDocking.dock(mol);
-		long t1 = System.currentTimeMillis();
 		alignedMol.stream().forEach(e -> confSet.add(new Conformer(e)));
 		for(Conformer c : confSet) {
 			if(c!=null) {
@@ -149,7 +167,7 @@ public class DockingEngine {
 				catch(Exception e) {
 					throw new DockingFailedException("could not assess atom types");
 				}
-				PositionConstraint constraint = new PositionConstraint(conf,50,0.2);
+				MMFFPositionConstraint constraint = new MMFFPositionConstraint(conf,50,0.2);
 				mmff.addEnergyTerm(constraint);
 				mmff.minimise();
 				Conformer ligConf = new Conformer(conf);
@@ -168,7 +186,136 @@ public class DockingEngine {
 		return eMin;
 		
 	}
+	/**
+	 * 1) find MCS between reference molecule and candidate molecule
+	 * 2) 
+	 * @param mol
+	 */
+	private double getStartingPositionsMCS(StereoMolecule mol, List<Conformer> initialPos) {
+		double eMin = Double.MAX_VALUE;
+		int rotBondsMol = mol.getRotatableBondCount();
+		int rotBondsRef = nativeLigand.getRotatableBondCount();
+		MCS mcs = new MCS();
+		//in MCS, the first molecule should have more rotatable bonds than the second 
+		boolean[] bondMCSMol = null;
+		boolean[] bondMCSFrag = null;
+		boolean isNativeLigFrag;//if native ligand is smaller than the candidate molecule
+		if(rotBondsMol > rotBondsRef) {
+			mcs.set(mol,nativeLigand);
+			bondMCSMol = new boolean[mol.getAllBonds()];
+			bondMCSFrag = new boolean[nativeLigand.getAllBonds()];
+			isNativeLigFrag = true;
+		}
+		else {
+			mcs.set(mol,nativeLigand);
+			bondMCSMol = new boolean[mol.getAllBonds()];
+			bondMCSFrag = new boolean[nativeLigand.getAllBonds()];
+			isNativeLigFrag = false;
+		}
+
+		StereoMolecule mcsMol = mcs.getMCS();
+		SSSearcher sss = new SSSearcher();
+		sss.setFragment(mcsMol);
+		sss.setMolecule(nativeLigand);
+		sss.findFragmentInMolecule(SSSearcher.cCountModeOverlapping, SSSearcher.cMatchDBondToDelocalized);
+		int[] map1 = sss.getMatchList().get(0);
+		sss.setMolecule(mol);
+		sss.findFragmentInMolecule(SSSearcher.cCountModeOverlapping, SSSearcher.cMatchDBondToDelocalized);
+		int[] map2 = sss.getMatchList().get(0);
+		mcsConstrainedAtoms = new ArrayList<>();
+		for(int a: map2) {
+			mcsConstrainedAtoms.add(a);
+		}
+		//map from reference to mol
+		Map<Integer,Integer> map = new HashMap<Integer,Integer>();
+		for(int i=0;i<map1.length;i++) {
+			map.put(map1[i], map2[i]);
+		}
+		//find bonds not part of MCS
+		mcs.getMCSBondArray(bondMCSMol, bondMCSFrag);
+		mcsConstrainedBonds = new ArrayList<>();
+		boolean[] bondArray = null;
+		if(isNativeLigFrag) 
+			bondArray = bondMCSMol;
+		else 
+			bondArray = bondMCSFrag;
+		for(int b=0;b<bondArray.length;b++) {
+				if(bondArray[b])
+					mcsConstrainedBonds.add(b);
+		}
+		ForceFieldMMFF94.initialize(ForceFieldMMFF94.MMFF94SPLUS);
+		Map<String, Object> ffOptions = new HashMap<String, Object>();
+		ffOptions.put("dielectric constant", 80.0);
+		ConformerSetGenerator csGen = new ConformerSetGenerator();
+		ConformerSet confSet = csGen.generateConformerSet(mol);
+		Map<Conformer,Double> confsWithRMSDs = new HashMap<Conformer,Double>();
+		for(Conformer conf : confSet) {
+			//align MCS using Kabsch algorithm
+			Coordinates[] coords1 = new Coordinates[map1.length];
+			Coordinates[] coords2 = new Coordinates[map1.length];
+			int counter = 0;
+			for(int key : map.keySet()) {
+				coords1[counter] = new Coordinates(nativeLigand.getCoordinates(key));
+				coords2[counter] = new Coordinates(conf.getCoordinates(map.get(key)));
+				counter++;
+			}
+			int[][] mapping = new int[coords1.length][2];
+			counter = 0;
+			for(int[] m : mapping) {
+				m[0] = counter;
+				m[1] = counter;
+				counter++;
+			}
+			Coordinates trans1 = new Coordinates();
+			Matrix rot = new Matrix(3,3); 
+			Coordinates trans2 = new Coordinates();
+			KabschAlignment alignment = new KabschAlignment(coords1,coords2,mapping);
+			alignment.align(trans1,rot,trans2);
+			double rmsd = getCoreRMSD(coords1,coords2);
+			for(Coordinates coord : conf.getCoordinates()) {
+				coord.add(trans1);
+				coord.rotate(rot.getArray());
+				coord.add(trans2);
+			}
+			confsWithRMSDs.put(conf, rmsd);
+		}
+		
+		LinkedHashMap<Conformer,Double> sortedConfs = confsWithRMSDs.entrySet().stream().sorted(Map.Entry.<Conformer,Double>comparingByValue()).collect(Collectors.toMap(Map.Entry::getKey,Map.Entry::getValue,
+				(e1, e2) -> e1, LinkedHashMap::new));
+		Iterator<Entry<Conformer,Double>> iterator = sortedConfs.entrySet().iterator();
+		boolean done = false;
+		Entry<Conformer,Double> entry;
+		int c=0;
+		while(!done && c<startPositions*MCS_EXHAUSTIVENESS) {
+			c++;
+			try	{
+				entry = iterator.next();
+			}
+			catch(Exception e) {
+				done = true;
+				continue;
+			}
+			Conformer conf = entry.getKey();
+			StereoMolecule aligned = new StereoMolecule();
+			aligned = conf.toMolecule(null);
+			aligned.ensureHelperArrays(Molecule.cHelperParities);
+			ForceFieldMMFF94 mmff;
+			List<MMFFExternalPositionConstraint> constraints = new ArrayList<>();
+			MMFFPositionConstraint constraint = new MMFFPositionConstraint(aligned,50,0.2);
+			mmff = new ForceFieldMMFF94(aligned, ForceFieldMMFF94.MMFF94SPLUS,ffOptions);
+			mmff.addEnergyTerm(constraint);
+			mmff.minimise();
+			double e = mmff.getTotalEnergy();
+			for(int a=0;a<aligned.getAllAtoms();a++) {
+				conf.setCoordinates(a, new Coordinates(aligned.getCoordinates(a)));
+			}
+			initialPos.add(conf);
+			if(e<eMin)
+				eMin = e;
+		}
+		return eMin;
 	
+	}
 	
 	
 	public DockingResult dockMolecule(StereoMolecule mol) throws DockingFailedException {
@@ -177,35 +324,38 @@ public class DockingEngine {
 		if(ForceFieldMMFF94.table(ForceFieldMMFF94.MMFF94SPLUS)==null)
 			ForceFieldMMFF94.initialize(ForceFieldMMFF94.MMFF94SPLUS);
 		List<Conformer> startPoints = new ArrayList<>();
-		double eMin = getStartingPositions(mol, startPoints);
+		double eMin = 0.0;
+		int steps = mcSteps;
+		if(mcsRef!=null) {
+			eMin = getStartingPositionsMCS(mol, startPoints);
+			steps=steps/MCS_EXHAUSTIVENESS;
+		}
+		else {
+			eMin = getStartingPositions(mol, startPoints);
+		}
+
 		Map<String,Double> contributions = null;
 		for(Conformer ligConf : startPoints) {
-			for(double[] transform : PheSAAlignment.initialTransform(0)) {
-				Conformer newLigConf = new Conformer(ligConf);
-				ExponentialMap eMap = new ExponentialMap(transform[0],transform[1],transform[2]);
-				Quaternion q = eMap.toQuaternion();
-				Coordinates com = DockingUtils.getCOM(newLigConf);
-				Rotation rot = new Rotation(q.getRotMatrix().getArray());
-				Translation trans = new Translation(new double[] {transform[3],transform[4],transform[5]});
-				TransformationSequence t = new TransformationSequence();
-				Translation t1 = new Translation(-com.x,-com.y,-com.z);
-				Translation t2 = new Translation(com.x,com.y,com.z);
-				t.addTransformation(t1);
-				t.addTransformation(rot);
-				t.addTransformation(t2);
-				t.addTransformation(trans);
-				t.apply(newLigConf);
-				LigandPose pose = initiate(newLigConf,eMin);
-				double energy = mcSearch(pose);
-				if(energy<bestEnergy) {
-					bestEnergy = energy;
-					bestPose = pose.getLigConf();
-					contributions = pose.getContributions();
+			Conformer newLigConf = new Conformer(ligConf);		
+			LigandPose pose = initiate(newLigConf,eMin);
+			if(mcsRef!=null) {
+				pose.setMCSBondConstraints(mcsConstrainedBonds);
+				for(int a : mcsConstrainedAtoms) {
+					PositionConstraint constr = new PositionConstraint(newLigConf,a,50,1.0);
+					pose.addConstraint(constr);
 				}
-				if(threadMaster!=null && threadMaster.threadMustDie())
-					break;
+				
 			}
+			double energy = mcSearch(pose,steps);
+			if(energy<bestEnergy) {
+				bestEnergy = energy;
+				bestPose = pose.getLigConf();
+				contributions = pose.getContributions();
+			}
+			if(threadMaster!=null && threadMaster.threadMustDie())
+				break;
 		}
+		
 		if(bestPose!=null) {
 			StereoMolecule best = bestPose.toMolecule();
 			Rotation rot = rotation.getInvert();
@@ -222,8 +372,13 @@ public class DockingEngine {
 	
 	
 
-	
-	private double mcSearch(LigandPose pose) {
+	/**
+	 * use monte carlo steps to permute molecular rotation, translation, torsion angles
+	 * promising poses (below a certain cutoff) are optimized
+	 * @param pose
+	 * @return
+	 */
+	private double mcSearch(LigandPose pose, int steps) {
 		double[] bestState = new double[pose.getState().length];
 		double[] oldState = new double[pose.getState().length];
 		double[] state = new double[pose.getState().length];
@@ -236,7 +391,7 @@ public class DockingEngine {
 		oldEnergy = pose.getFGValue(new double[bestState.length]);
 		bestEnergy = oldEnergy;
 	
-		for(int i=0;i<mcSteps;i++) {
+		for(int i=0;i<steps;i++) {
 			pose.randomPerturbation(random);
 			double energyMC = pose.getFGValue(new double[bestState.length]);
 			if(energyMC<MINI_CUTOFF) {
@@ -274,6 +429,8 @@ public class DockingEngine {
 		}
 
 		pose.setState(bestState);
+		pose.removeConstraints();
+		bestEnergy = pose.getFGValue(new double[bestState.length]);
 		return bestEnergy;
 		
 	}
@@ -286,6 +443,10 @@ public class DockingEngine {
 		
 	}
 	
+	public void setMCSReference(StereoMolecule referencePose) {
+		mcsRef = referencePose;
+	}
+	
 	public static void getBindingSiteAtoms(StereoMolecule receptor, Set<Integer> bindingSiteAtoms, MoleculeGrid grid,
 			boolean includeHydrogens) {
 		int[] gridSize = grid.getGridSize();
@@ -341,15 +502,12 @@ public class DockingEngine {
 	public double refineNativePose(double d, double[] coords) {
 		Map<String, Object> ffOptions = new HashMap<String, Object>();
 		ffOptions.put("dielectric constant", 80.0);
-		//ffOptions.put("angle bend", false);
-		//ffOptions.put("stretch bend", false);
-		//ffOptions.put("bond stretch", false);
 		if(ForceFieldMMFF94.table(ForceFieldMMFF94.MMFF94SPLUS)==null)
 			ForceFieldMMFF94.initialize(ForceFieldMMFF94.MMFF94SPLUS);
 		Molecule3D nativePose = new Molecule3D(nativeLigand);
 		new Canonizer(nativePose);
 		ForceFieldMMFF94 mmff = new ForceFieldMMFF94(nativePose, ForceFieldMMFF94.MMFF94SPLUS, ffOptions);
-		PositionConstraint constraint = new PositionConstraint(nativePose,50,0.2);
+		MMFFPositionConstraint constraint = new MMFFPositionConstraint(nativePose,50,0.2);
 		mmff.addEnergyTerm(constraint);
 		mmff.minimise();
 		ConformerSetGenerator confSetGen = new ConformerSetGenerator(100,ConformerGenerator.STRATEGY_LIKELY_RANDOM, false,
@@ -363,7 +521,7 @@ public class DockingEngine {
 				StereoMolecule conf = conformer.toMolecule();
 				conf.ensureHelperArrays(Molecule.cHelperParities);
 				mmff = new ForceFieldMMFF94(conf, ForceFieldMMFF94.MMFF94SPLUS, ffOptions);
-				constraint = new PositionConstraint(conf,50,0.2);
+				constraint = new MMFFPositionConstraint(conf,50,0.2);
 				mmff.addEnergyTerm(constraint);
 				mmff.minimise();
 				Conformer ligConf = new Conformer(conf);
@@ -396,6 +554,18 @@ public class DockingEngine {
 		
 	}
 	
+	private static double getCoreRMSD(Coordinates[] coords1, Coordinates[] coords2) {
+		double rmsd = 0.0;
+		for(int i=0;i<coords1.length;i++) { 
+			Coordinates c1 = coords1[i];
+			Coordinates c2 = coords2[i];
+			rmsd+=c1.distanceSquared(c2);
+		}
+		rmsd/=coords1.length;
+		rmsd = Math.sqrt(rmsd);
+		return rmsd;
+	}
+	
 	public static class DockingResult implements Comparable<DockingResult>  {
 		private double score;
 		private StereoMolecule pose;
@@ -486,6 +656,7 @@ public class DockingEngine {
             return Double.compare( this.score, o.score);
            
 		}
+		
 	}
 	
 


=====================================
src/main/java/com/actelion/research/chem/docking/LigandPose.java
=====================================
@@ -1,7 +1,9 @@
 package com.actelion.research.chem.docking;
 
+import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.HashMap;
+import java.util.List;
 import java.util.Map;
 import java.util.Random;
 
@@ -21,6 +23,7 @@ import com.actelion.research.chem.conf.TorsionDB;
 import com.actelion.research.chem.docking.scoring.AbstractScoringEngine;
 import com.actelion.research.chem.optimization.Evaluable;
 import com.actelion.research.chem.potentialenergy.PositionConstraint;
+import com.actelion.research.chem.potentialenergy.PotentialEnergyTerm;
 
 
 public class LigandPose implements Evaluable{
@@ -39,6 +42,8 @@ public class LigandPose implements Evaluable{
 	private AbstractScoringEngine engine;
 	public static long SEED = 12345L;
 	private Coordinates origCOM;
+	private int[] mcsRotBondIndeces; //for MCS docking, only bonds not part of the MCS are sampled
+	
 	
 	public LigandPose(Conformer ligConf, AbstractScoringEngine engine, double e0) {
 		
@@ -47,6 +52,25 @@ public class LigandPose implements Evaluable{
 		init(e0);
 	}
 	
+	/**
+	 * for MCS docking:create array of bond indices that are allowed to be permuted
+	 * @param constraints
+	 */
+	public void setMCSBondConstraints(List<Integer> constraints) {
+		int[] rotBonds = torsionHelper.getRotatableBonds();
+		List<Integer> allowedIndices = new ArrayList<>();
+		for(int rbIndex=0;rbIndex<rotBonds.length;rbIndex++) {
+			int rb = rotBonds[rbIndex];
+			if(!constraints.contains(rb))
+				allowedIndices.add(rbIndex);
+		}
+		mcsRotBondIndeces = new int[allowedIndices.size()];
+		for(int i=0;i<allowedIndices.size();i++) {
+			mcsRotBondIndeces[i]=allowedIndices.get(i);
+		}
+		
+	}
+	
 	private void init(double e0) {
 		mol = ligConf.getMolecule();
 		//rotationCenter = calcCOM();
@@ -255,46 +279,67 @@ public class LigandPose implements Evaluable{
 	}
 	
 	public void randomPerturbation(Random random) {
-		int num = (int) (3*random.nextDouble());
-		if(num==0) { //translation
-			Coordinates shift = DockingUtils.randomVectorInSphere(random).scale(MOVE_AMPLITUDE);
-			state[0]+=shift.x;
-			state[1]+=shift.y;
-			state[2]+=shift.z;
-			/*
-			for(int a=0;a<ligConf.getMolecule().getAllAtoms();a++) { 
-				Coordinates c = ligConf.getCoordinates(a);
-				c.add(shift);
+		if(mcsRotBondIndeces==null) {
+			int num = (int) (3*random.nextDouble());
+			if(num==0) { //translation
+				Coordinates shift = DockingUtils.randomVectorInSphere(random).scale(MOVE_AMPLITUDE);
+				state[0]+=shift.x;
+				state[1]+=shift.y;
+				state[2]+=shift.z;
+				/*
+				for(int a=0;a<ligConf.getMolecule().getAllAtoms();a++) { 
+					Coordinates c = ligConf.getCoordinates(a);
+					c.add(shift);
+				}
+				*/
 			}
-			*/
-		}
-		else if(num==1) {
-			double r = getGyrationRadius();
-			Coordinates rot = DockingUtils.randomVectorInSphere(random).scale(MOVE_AMPLITUDE/r);
-			double angle = rot.dist();
-			Quaternion q = new Quaternion(1.0,0.0,0.0,0.0);
-			if(angle>0.0001) {
-				Coordinates axis = rot.scale(1.0/angle);
-				q = new Quaternion(axis,angle);
+			else if(num==1) {
+				double r = getGyrationRadius();
+				Coordinates rot = DockingUtils.randomVectorInSphere(random).scale(MOVE_AMPLITUDE/r);
+				double angle = rot.dist();
+				Quaternion q = new Quaternion(1.0,0.0,0.0,0.0);
+				if(angle>0.0001) {
+					Coordinates axis = rot.scale(1.0/angle);
+					q = new Quaternion(axis,angle);
+				}
+				ExponentialMap em = new ExponentialMap(state[3],state[4],state[5]);
+				Quaternion qOrig = em.toQuaternion();
+				q.multiply(qOrig);
+				ExponentialMap emNew = new ExponentialMap(q);
+				state[3] = emNew.getP().x;
+				state[4] = emNew.getP().y;
+				state[5] = emNew.getP().z;
+	
+				
+			}
+			else  {
+				if(torsionHelper.getRotatableBonds().length==0)
+					return;
+				double rnd = random.nextDouble();
+				rnd*=180.0;
+				double rotateBy = random.nextBoolean() ? rnd : -rnd;
+				rotateBy = rotateBy*Math.PI/180.0;
+				int bond = random.nextInt(torsionHelper.getRotatableBonds().length);
+				double previousAngle = state[6+bond];
+				double newAngle = previousAngle+rotateBy; 
+				if(newAngle>Math.PI) {
+					newAngle -= 2*Math.PI;
+				}
+				state[6+bond] = newAngle;
+				
 			}
-			ExponentialMap em = new ExponentialMap(state[3],state[4],state[5]);
-			Quaternion qOrig = em.toQuaternion();
-			q.multiply(qOrig);
-			ExponentialMap emNew = new ExponentialMap(q);
-			state[3] = emNew.getP().x;
-			state[4] = emNew.getP().y;
-			state[5] = emNew.getP().z;
-
-			
 		}
-		else  {
+		//mcs docking
+		else {
 			if(torsionHelper.getRotatableBonds().length==0)
 				return;
+			else if(mcsRotBondIndeces.length==0)
+				return;
 			double rnd = random.nextDouble();
 			rnd*=180.0;
 			double rotateBy = random.nextBoolean() ? rnd : -rnd;
 			rotateBy = rotateBy*Math.PI/180.0;
-			int bond = random.nextInt(torsionHelper.getRotatableBonds().length);
+			int bond = mcsRotBondIndeces[random.nextInt(mcsRotBondIndeces.length)];
 			double previousAngle = state[6+bond];
 			double newAngle = previousAngle+rotateBy; 
 			if(newAngle>Math.PI) {
@@ -302,7 +347,9 @@ public class LigandPose implements Evaluable{
 			}
 			state[6+bond] = newAngle;
 			
+			
 		}
+		
 		updateLigandCoordinates();
 		
 	}
@@ -314,20 +361,20 @@ public class LigandPose implements Evaluable{
 		}
 	}
 	
+	public void addConstraint(PositionConstraint constraint) {
+		engine.addConstraint(constraint);
+	}
+	
+	public void removeConstraints() {
+		engine.removeConstraints();
+	}
+	
+	
 	public Conformer getLigConf() {
 		return ligConf;
 	}
 	
-	private Coordinates calcCOM(){ 
-		Coordinates com = new Coordinates();
-		for(int a=0;a<ligConf.getMolecule().getAtoms();a++){
-			com.add(ligConf.getCoordinates(a));
-		}
-		com.scale(1.0/ligConf.getMolecule().getAtoms());
-		return com;
-		
 
-	}
 	
 
 }


=====================================
src/main/java/com/actelion/research/chem/docking/receptorpharmacophore/NegativeReceptorImage.java
=====================================
@@ -1,37 +1,23 @@
 package com.actelion.research.chem.docking.receptorpharmacophore;
 
-import java.util.ArrayList;
-import java.util.Arrays;
-import java.util.HashMap;
-import java.util.HashSet;
-import java.util.LinkedList;
-import java.util.List;
-import java.util.Map;
-import java.util.Queue;
-import java.util.Random;
-import java.util.Set;
-import java.util.stream.IntStream;
-
 import com.actelion.research.chem.Coordinates;
 import com.actelion.research.chem.StereoMolecule;
-import com.actelion.research.chem.descriptor.pharmacophoretree.IonizableGroupDetector2D;
 import com.actelion.research.chem.docking.DockingUtils;
-import com.actelion.research.chem.docking.ScoringTask;
 import com.actelion.research.chem.docking.scoring.ProbeScanning;
-import com.actelion.research.chem.interactionstatistics.InteractionAtomTypeCalculator;
 import com.actelion.research.chem.io.pdb.converter.MoleculeGrid;
-import com.actelion.research.chem.phesa.MolecularVolume;
 import com.actelion.research.chem.phesa.AtomicGaussian;
-import com.actelion.research.chem.phesa.ShapeVolume;
 import com.actelion.research.chem.phesa.Gaussian3D;
+import com.actelion.research.chem.phesa.MolecularVolume;
+import com.actelion.research.chem.phesa.ShapeVolume;
 import com.actelion.research.chem.phesa.pharmacophore.PharmacophoreCalculator;
-import com.actelion.research.chem.phesa.pharmacophore.pp.ChargePoint;
 import com.actelion.research.chem.phesa.pharmacophore.pp.IPharmacophorePoint;
 import com.actelion.research.chem.phesa.pharmacophore.pp.PPGaussian;
 import com.actelion.research.chem.phesa.pharmacophore.pp.SimplePharmacophorePoint;
-
 import smile.clustering.KMeans;
 
+import java.util.*;
+import java.util.stream.IntStream;
+
 
 public class NegativeReceptorImage extends MoleculeGrid {
 	
@@ -324,9 +310,19 @@ public class NegativeReceptorImage extends MoleculeGrid {
 	
 	private void analyzeBumps() {
 		double radiusSq = BUMP_RADIUS2*BUMP_RADIUS2;
-		for(int x=0;x<this.gridSize[0];x++) {
-			for(int y=0;y<this.gridSize[1];y++) {
-				for(int z=0;z<this.gridSize[2];z++) {
+
+// Replaced because of out-of-bounds exceptions; TLS 17Oct2021
+//		for(int x=0;x<this.gridSize[0];x++) {
+//			for(int y=0;y<this.gridSize[1];y++) {
+//				for(int z=0;z<this.gridSize[2];z++) {
+
+		int xmax = Math.min(this.gridSize[0], bumpGrid.length);
+		int ymax = Math.min(this.gridSize[1], bumpGrid[0].length);
+		int zmax = Math.min(this.gridSize[2], bumpGrid[0][0].length);
+		for(int x=0;x<xmax;x++) {
+			for(int y=0;y<ymax;y++) {
+				for(int z=0;z<zmax;z++) {
+
 					Coordinates probeCoords = this.getCartCoordinates(new int[] {x,y,z});
 					for(int atom=0;atom<receptor.getAtoms();atom++) {
 						Coordinates c = receptor.getCoordinates(atom);


=====================================
src/main/java/com/actelion/research/chem/docking/scoring/AbstractScoringEngine.java
=====================================
@@ -72,6 +72,10 @@ public abstract class AbstractScoringEngine  {
 	public void addConstraint(PotentialEnergyTerm constraint) {
 		this.constraints.add(constraint);
 	}
+	
+	public void removeConstraints() {
+		this.constraints = new ArrayList<>();
+	}
 
 	public abstract void init(LigandPose candidatePose, double e0);
 	


=====================================
src/main/java/com/actelion/research/chem/forcefield/mmff/MMFFExternalPositionConstraint.java
=====================================
@@ -0,0 +1,56 @@
+package com.actelion.research.chem.forcefield.mmff;
+
+import com.actelion.research.chem.StereoMolecule;
+
+public class MMFFExternalPositionConstraint implements EnergyTerm {
+	private double[] refPos;
+	private int constrainedAtom;
+	private double k;
+	private double d;
+	
+	
+	public MMFFExternalPositionConstraint (int atom, double[] refPos, double k, double d) {
+		this.refPos = refPos;
+		this.constrainedAtom = atom;
+		this.k = k;
+		this.d = d;
+	}
+	
+
+	@Override
+	public double getEnergy(double[] pos) {
+		double energy = 0.0;
+		double dx = pos[3*constrainedAtom]-refPos[0];
+		double dy = pos[3*constrainedAtom+1]-refPos[1];
+		double dz = pos[3*constrainedAtom+2]-refPos[2];
+		double dist = Math.sqrt(dx*dx + dy*dy + dz*dz);
+		double prefactor = 0.0;
+		if(dist>d)
+			prefactor = dist-d;
+		else 
+			prefactor = 0.0;
+		energy+=0.5*k*prefactor*prefactor;
+		
+		return energy;
+	}
+
+	@Override
+	public void getGradient(double[] pos, double[] grad) {
+			int a = constrainedAtom;
+			double dx = pos[3*a]-refPos[0];
+			double dy = pos[3*a+1]-refPos[1];
+			double dz = pos[3*a+2]-refPos[2];
+			double dist = Math.sqrt(dx*dx + dy*dy + dz*dz);
+			double prefactor = 0.0;
+			if(dist>d)
+				prefactor = dist-d;
+			else 
+				prefactor = 0.0;
+			grad[a] += prefactor*dx/Math.max(dist, 1E-08);
+			grad[a+1] += prefactor*dy/Math.max(dist, 1E-08);
+			grad[a+2] += prefactor*dz/Math.max(dist, 1E-08);
+	}
+		
+	
+
+}


=====================================
src/main/java/com/actelion/research/chem/forcefield/mmff/PositionConstraint.java → src/main/java/com/actelion/research/chem/forcefield/mmff/MMFFPositionConstraint.java
=====================================
@@ -2,14 +2,14 @@ package com.actelion.research.chem.forcefield.mmff;
 
 import com.actelion.research.chem.StereoMolecule;
 
-public class PositionConstraint implements EnergyTerm {
+public class MMFFPositionConstraint implements EnergyTerm {
 	private double[] refPos;
 	private boolean[] constrained;
 	private double k;
 	private double d;
 	
 	
-	public PositionConstraint (StereoMolecule mol, double k, double d) {
+	public MMFFPositionConstraint (StereoMolecule mol, double k, double d) {
 		refPos = new double[3*mol.getAllAtoms()];
 		constrained = new boolean[mol.getAllAtoms()];
 		for(int a=0;a<mol.getAllAtoms();a++) {


=====================================
src/main/java/com/actelion/research/chem/io/CompoundFileHelper.java
=====================================
@@ -69,7 +69,9 @@ public abstract class CompoundFileHelper {
 	public static final int cFileTypeMOL = 0x00010000;
 	public static final int cFileTypeMOL2 = 0x00020000;
 	public static final int cFileTypePDB = 0x00040000;
-	public static final int cFileTypeSDGZ = 0x00080000;
+	public static final int cFileTypeMMTF = 0x00080000;
+	public static final int cFileTypeProtein = 0x000C0000;
+	public static final int cFileTypeSDGZ = 0x00100000;
     public static final int cFileTypeUnknown = -1;
 	public static final int cFileTypeDirectory = -2;
 
@@ -379,6 +381,10 @@ public abstract class CompoundFileHelper {
 			filter.addExtension("pdb");
 			filter.addDescription("Protein Data Bank files");
 			}
+		if ((filetypes & cFileTypeMMTF) != 0) {
+			filter.addExtension("mmtf");
+			filter.addDescription("Binary Protein Data Bank files");
+			}
 		if ((filetypes & cFileTypeMOL) != 0) {
 			filter.addExtension("mol");
 			filter.addDescription("MDL Molfiles");
@@ -490,6 +496,8 @@ public abstract class CompoundFileHelper {
 			return cFileTypeMOL2;
 		if (extension.equals(".pdb"))
 			return cFileTypePDB;
+		if (extension.equals(".mmtf"))
+			return cFileTypeMMTF;
 
         return cFileTypeUnknown;
         }
@@ -574,7 +582,10 @@ public abstract class CompoundFileHelper {
 		case cFileTypePDB:
 			extension = ".pdb";
 			break;
-			case cFileTypeSDGZ:
+		case cFileTypeMMTF:
+			extension = ".mmtf";
+			break;
+		case cFileTypeSDGZ:
 			extension = ".sdf.gz";
 			break;
 			}


=====================================
src/main/java/com/actelion/research/chem/mcs/MCS.java
=====================================
@@ -37,13 +37,12 @@ import java.util.*;
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 */
-public class MCS
-{
-	public static final int PAR_CLEAVE_RINGS=0; 
-	public static final int PAR_KEEP_RINGS=1;
-	public static final int PAR_KEEP_AROMATIC_RINGS=2; 
+public class MCS {
+	public static final int PAR_CLEAVE_RINGS = 0;
+	public static final int PAR_KEEP_RINGS = 1;
+	public static final int PAR_KEEP_AROMATIC_RINGS = 2;
 	private static final boolean DEBUG = false;
-	private static final int CAPACITY = (60+60)*10; 
+	private static final int CAPACITY = (60 + 60) * 10;
 	private StereoMolecule mol;
 	private StereoMolecule frag;
 	private HashSet<IntVec> hsIndexFragCandidates;
@@ -56,74 +55,71 @@ public class MCS
 	private boolean considerRings;
 	// The largest valid solutions.
 	private List<IntVec> liMCSSolutions;
-	private HashMap<Integer, List<int []>> hmRingBnd_ListRingBnds;
-	private HashMap<Integer, List<int []>> hmAromaticRingBnd_ListRingBnds;
+	private HashMap<Integer, List<int[]>> hmRingBnd_ListRingBnds;
+	private HashMap<Integer, List<int[]>> hmAromaticRingBnd_ListRingBnds;
 	private RingCollection ringCollection;
-    boolean excluded[] = null;
+	boolean excluded[] = null;
+	private int [] arrMatchListFrag2Mol;
 
-    public MCS()
-    {
-        this(PAR_CLEAVE_RINGS,null);
+	public MCS() {
+		this(PAR_CLEAVE_RINGS, null);
 	}
-	
-    public MCS(int ringStatus)
-    {
-        this(ringStatus,null);
-    }
-		
-    public MCS(int ringStatus,SSSearcher searcher)
-    {
-		considerRings=false;
-		considerAromaticRings=false;
+
+	public MCS(int ringStatus) {
+		this(ringStatus, null);
+	}
+
+	public MCS(int ringStatus, SSSearcher searcher) {
+		considerRings = false;
+		considerAromaticRings = false;
 		switch (ringStatus) {
-		case PAR_CLEAVE_RINGS:
-			break;
-		case PAR_KEEP_RINGS:
-			considerRings=true;
-			break;
-		case PAR_KEEP_AROMATIC_RINGS:
-			considerAromaticRings=true;
-			break;
-		default:
-			break;
+			case PAR_CLEAVE_RINGS:
+				break;
+			case PAR_KEEP_RINGS:
+				considerRings = true;
+				break;
+			case PAR_KEEP_AROMATIC_RINGS:
+				considerAromaticRings = true;
+				break;
+			default:
+				break;
 		}
 		hsIndexFragCandidates = new HashSet<IntVec>(CAPACITY);
 		hsIndexFragGarbage = new HashSet<IntVec>(CAPACITY);
 		hsIndexFragSolution = new HashSet<IntVec>(CAPACITY);
 		liMCSSolutions = new ArrayList<IntVec>(CAPACITY);
-        if (searcher == null)
-		sss = new SSSearcher();
-        else
-            sss = searcher;
+		if (searcher == null)
+			sss = new SSSearcher();
+		else
+			sss = searcher;
 		comparatorBitsSet = new ComparatorBitsSet();
-		hmRingBnd_ListRingBnds = new HashMap<Integer, List<int []>>();
-		hmAromaticRingBnd_ListRingBnds = new HashMap<Integer, List<int []>>();
+		hmRingBnd_ListRingBnds = new HashMap<Integer, List<int[]>>();
+		hmAromaticRingBnd_ListRingBnds = new HashMap<Integer, List<int[]>>();
+	}
+
+	public void setSSSearcher(SSSearcher sss) {
+		this.sss = sss;
 	}
-	
-    public void setSSSearcher(SSSearcher sss)
-    {
-        this.sss  = sss;
-    }
 
 	/**
 	 * mol should contain equal or more bonds than frag.
+	 *
 	 * @param mol
 	 * @param frag
 	 */
-    public void set(StereoMolecule mol, StereoMolecule frag)
-    {
-        set(mol,frag,null);
+	public void set(StereoMolecule mol, StereoMolecule frag) {
+		set(mol, frag, null);
 	}
 
 	/**
 	 * mol should contain equal or more bonds than frag.
 	 * If frag contains more than one molecule only the biggest one is considered.
+	 *
 	 * @param mol
 	 * @param frag
 	 * @param excluded
 	 */
-    public void set(StereoMolecule mol, StereoMolecule frag, boolean excluded[])
-    {
+	public void set(StereoMolecule mol, StereoMolecule frag, boolean excluded[]) {
 
 		StereoMolecule fragBiggestSub = new StereoMolecule(frag);
 
@@ -135,15 +131,14 @@ public class MCS
 
 		this.mol = mol;
 
-        this.frag = fragBiggestSub;
+		this.frag = fragBiggestSub;
 
 		this.excluded = excluded;
 
-        init();
-    }
+		init();
+	}
 
-    private void init()
-    {
+	private void init() {
 		hsIndexFragCandidates.clear();
 		hsIndexFragGarbage.clear();
 		hsIndexFragSolution.clear();
@@ -152,32 +147,31 @@ public class MCS
 		hmAromaticRingBnd_ListRingBnds.clear();
 		initCandidates();
 	}
-	
-    private void initCandidates()
-    {
+
+	private void initCandidates() {
 		ringCollection = frag.getRingSet();
 		int rings = ringCollection.getSize();
 		for (int i = 0; i < rings; i++) {
-			int [] arrIndexBnd = ringCollection.getRingBonds(i);
+			int[] arrIndexBnd = ringCollection.getRingBonds(i);
 			for (int j = 0; j < arrIndexBnd.length; j++) {
-				if(!hmRingBnd_ListRingBnds.containsKey(arrIndexBnd[j])){
-					hmRingBnd_ListRingBnds.put(arrIndexBnd[j], new ArrayList<int []>());
+				if (!hmRingBnd_ListRingBnds.containsKey(arrIndexBnd[j])) {
+					hmRingBnd_ListRingBnds.put(arrIndexBnd[j], new ArrayList<int[]>());
 				}
-				List<int []> li = hmRingBnd_ListRingBnds.get(arrIndexBnd[j]);
+				List<int[]> li = hmRingBnd_ListRingBnds.get(arrIndexBnd[j]);
 				li.add(arrIndexBnd);
 			}
-			if(ringCollection.isAromatic(i)){
+			if (ringCollection.isAromatic(i)) {
 				for (int j = 0; j < arrIndexBnd.length; j++) {
-					if(!hmAromaticRingBnd_ListRingBnds.containsKey(arrIndexBnd[j])){
-						hmAromaticRingBnd_ListRingBnds.put(arrIndexBnd[j], new ArrayList<int []>());
+					if (!hmAromaticRingBnd_ListRingBnds.containsKey(arrIndexBnd[j])) {
+						hmAromaticRingBnd_ListRingBnds.put(arrIndexBnd[j], new ArrayList<int[]>());
 					}
-					List<int []> li = hmAromaticRingBnd_ListRingBnds.get(arrIndexBnd[j]);
+					List<int[]> li = hmAromaticRingBnd_ListRingBnds.get(arrIndexBnd[j]);
 					li.add(arrIndexBnd);
 				}
 			}
 		}
 		// Array length for bit list.
-		int nInts = (int)(((double)frag.getBonds() / Integer.SIZE) + ((double)(Integer.SIZE-1)/Integer.SIZE));
+		int nInts = (int) (((double) frag.getBonds() / Integer.SIZE) + ((double) (Integer.SIZE - 1) / Integer.SIZE));
 		// The vector is used as bit vector.
 		// For each bond one bit.
 		for (int i = 0; i < frag.getBonds(); i++) {
@@ -186,28 +180,27 @@ public class MCS
 			hsIndexFragCandidates.add(iv);
 		}
 	}
-	
-    private void setBitAndAddRelatedRingBonds(int bit, IntVec iv)
-    {
-		if(!considerAromaticRings && !considerRings) {
+
+	private void setBitAndAddRelatedRingBonds(int bit, IntVec iv) {
+		if (!considerAromaticRings && !considerRings) {
 			iv.setBit(bit);
-		} else if(considerRings) {
+		} else if (considerRings) {
 			iv.setBit(bit);
-			if(hmRingBnd_ListRingBnds.containsKey(bit)){
-				List<int []> li = hmRingBnd_ListRingBnds.get(bit);
+			if (hmRingBnd_ListRingBnds.containsKey(bit)) {
+				List<int[]> li = hmRingBnd_ListRingBnds.get(bit);
 				for (int i = 0; i < li.size(); i++) {
-					int [] a = li.get(i);
+					int[] a = li.get(i);
 					for (int j = 0; j < a.length; j++) {
 						iv.setBit(a[j]);
 					}
 				}
 			}
-		} else if(considerAromaticRings) {
+		} else if (considerAromaticRings) {
 			iv.setBit(bit);
-			if(hmAromaticRingBnd_ListRingBnds.containsKey(bit)){
-				List<int []> li = hmAromaticRingBnd_ListRingBnds.get(bit);
+			if (hmAromaticRingBnd_ListRingBnds.containsKey(bit)) {
+				List<int[]> li = hmAromaticRingBnd_ListRingBnds.get(bit);
 				for (int i = 0; i < li.size(); i++) {
-					int [] a = li.get(i);
+					int[] a = li.get(i);
 					for (int j = 0; j < a.length; j++) {
 						iv.setBit(a[j]);
 					}
@@ -216,17 +209,16 @@ public class MCS
 		}
 		iv.calculateHashCode();
 	}
-	
+
 	/**
 	 * Checks first fragment for being sub structure of molecule. If true, frag is returned.
 	 * The returned IntVec contains the bond information of frag. The IntVec is used bit wise. Each bit corresponds to
 	 * one bond in the fragment. The information for the matching substructure in the molecule is not contained in the
 	 * IntVec.
-     *
+	 *
 	 * @return
 	 */
-    private List<IntVec> getAllSolutionsForCommonSubstructures()
-    {
+	private List<IntVec> getAllSolutionsForCommonSubstructures() {
 		sss.setMolecule(mol);
 		frag.setFragment(true);
 		sss.setFragment(frag);
@@ -234,10 +226,10 @@ public class MCS
 		// The MCS is the complete fragment.
 		//
 		try {
-            if (sss.findFragmentInMolecule(SSSearcher.cCountModeOverlapping, SSSearcher.cMatchDBondToDelocalized, excluded) > 0) {
+			if (sss.findFragmentInMolecule(SSSearcher.cCountModeOverlapping, SSSearcher.cMatchDBondToDelocalized, excluded) > 0) {
 				molMCS = frag;
 				List<IntVec> liIndexFragCandidates = new ArrayList<IntVec>(hsIndexFragCandidates);
-                if (liIndexFragCandidates != null && !liIndexFragCandidates.isEmpty()) {
+				if (liIndexFragCandidates != null && !liIndexFragCandidates.isEmpty()) {
 					IntVec iv = liIndexFragCandidates.get(0);
 					iv.setBits(0, iv.sizeBits());
 					iv.calculateHashCode();
@@ -245,70 +237,69 @@ public class MCS
 					li.add(iv);
 					return li;
 				}
-            }
+			}
 		} catch (Exception e) {
 			e.printStackTrace();
 		}
 		int maxSizeCandidates = 0;
-		while(!hsIndexFragCandidates.isEmpty()){
+		while (!hsIndexFragCandidates.isEmpty()) {
 			List<IntVec> liIndexFragCandidates = new ArrayList<IntVec>(hsIndexFragCandidates);
 			Collections.sort(liIndexFragCandidates, comparatorBitsSet);
 			// Get largest mcs.
-			IntVec iv = liIndexFragCandidates.get(liIndexFragCandidates.size()-1);
+			IntVec iv = liIndexFragCandidates.get(liIndexFragCandidates.size() - 1);
 			hsIndexFragCandidates.remove(iv);
-			if(DEBUG) {
+			if (DEBUG) {
 				System.out.println("Bits set " + iv.getBitsSet());
-				if(iv.getBitsSet() == frag.getBonds()){
+				if (iv.getBitsSet() == frag.getBonds()) {
 					System.out.println("Full structure in iv.");
 				}
 			}
 			StereoMolecule fragSub = getSubFrag(frag, iv);
 			sss.setFragment(fragSub);
-			if(sss.findFragmentInMolecule(SSSearcher.cCountModeOverlapping, SSSearcher.cDefaultMatchMode,excluded)>0){
+			if (sss.findFragmentInMolecule(SSSearcher.cCountModeOverlapping, SSSearcher.cDefaultMatchMode, excluded) > 0) {
 				hsIndexFragSolution.add(iv);
 				removeAllSubSolutions(iv);
-				if(iv.getBitsSet() != frag.getBonds()) {
+				if (iv.getBitsSet() != frag.getBonds()) {
 					List<IntVec> liIV = getAllPlusOneAtomCombinations(iv, frag);
 					for (IntVec ivPlus : liIV) {
-						if(DEBUG) {
-							if(ivPlus.getBitsSet() == frag.getBonds()){
+						if (DEBUG) {
+							if (ivPlus.getBitsSet() == frag.getBonds()) {
 								System.out.println("Full structure in ivPlus.");
 							}
 						}
-						if((!hsIndexFragGarbage.contains(ivPlus)) && (!hsIndexFragSolution.contains(ivPlus))){
-							hsIndexFragCandidates.add(ivPlus);	
+						if ((!hsIndexFragGarbage.contains(ivPlus)) && (!hsIndexFragSolution.contains(ivPlus))) {
+							hsIndexFragCandidates.add(ivPlus);
 						}
 					}
-                    if (DEBUG) {
+					if (DEBUG) {
 						System.out.println("tsIndexFragCandidates " + hsIndexFragCandidates.size());
+					}
 				}
-                }
-                if (maxSizeCandidates < hsIndexFragCandidates.size()) {
+				if (maxSizeCandidates < hsIndexFragCandidates.size()) {
 					maxSizeCandidates = hsIndexFragCandidates.size();
-                }
+				}
 			} else {
 				hsIndexFragGarbage.add(iv);
 			}
 		}
-		if(hsIndexFragSolution.size()==0){
+		if (hsIndexFragSolution.size() == 0) {
 			return null;
 		}
 		return getFinalSolutionSet(hsIndexFragSolution);
 	}
-	
+
 	/**
 	 * All molecules which are sub structures of an other molecule in the list are removed.
-     *
+	 *
 	 * @return
 	 */
-    public LinkedList<StereoMolecule> getAllCommonSubstructures()
-    {
+	public LinkedList<StereoMolecule> getAllCommonSubstructures() {
 		List<IntVec> liIndexFragSolution = getAllSolutionsForCommonSubstructures();
-		if(liIndexFragSolution==null){
+		if (liIndexFragSolution == null) {
 			return null;
 		}
 		Collections.sort(liIndexFragSolution, comparatorBitsSet);
-		IntVec ivMCS = liIndexFragSolution.get(liIndexFragSolution.size()-1);
+		IntVec ivMCS = liIndexFragSolution.get(liIndexFragSolution.size() - 1);
 		molMCS = getSubFrag(frag, ivMCS);
 		List<StereoMolecule> li = new ArrayList<StereoMolecule>();
 		for (IntVec iv : liIndexFragSolution) {
@@ -320,28 +311,27 @@ public class MCS
 	/**
 	 * @return maximum common substructure at top of list or null of none common MCS was found.
 	 */
-    public StereoMolecule getMCS()
-    {
-		List<IntVec> liIndexFragSolution = getAllSolutionsForCommonSubstructures ();
-		if(liIndexFragSolution==null){
+	public StereoMolecule getMCS() {
+		List<IntVec> liIndexFragSolution = getAllSolutionsForCommonSubstructures();
+		if (liIndexFragSolution == null) {
 			return null;
 		}
 		Collections.sort(liIndexFragSolution, comparatorBitsSet);
-		IntVec ivMCS = liIndexFragSolution.get(liIndexFragSolution.size()-1);
+		IntVec ivMCS = liIndexFragSolution.get(liIndexFragSolution.size() - 1);
 		molMCS = getSubFrag(frag, ivMCS);
 		return molMCS;
 	}
-	
+
 	/**
 	 * Calculates the bond arrays for molecule and fragment for their maximum common substructure.
-     *
-	 * @param arrBondMCSMol result array. Null or initialized with the length number of bonds in molecule. 
-	 * @param arrBondFrag result array. Null or initialized with the length number of bonds in fragment.
-	 * @return two arrays, in the first array the bits are set 'true' which corresponds to the bonds for maximum common substructure in the molecule. 
+	 *
+	 * @param arrBondMCSMol result array. Null or initialized with the length number of bonds in molecule.
+	 * @param arrBondFrag   result array. Null or initialized with the length number of bonds in fragment.
+	 * @return two arrays, in the first array the bits are set 'true' which corresponds to the bonds for maximum common substructure in the molecule..
 	 * In the second are the bits set 'true' which corresponds to the bonds for maximum common substructure in the fragment.
 	 */
-    public boolean[][] getMCSBondArray(boolean[] arrBondMCSMol, boolean[] arrBondFrag)
-    {
+
+    public boolean[][] getMCSBondArray(boolean[] arrBondMCSMol, boolean[] arrBondFrag) {
 		boolean [][] arrBondMol_Result = new boolean [2][];
 		List<IntVec> liIndexFragSolution = getAllSolutionsForCommonSubstructures();
 		if(liIndexFragSolution==null){
@@ -366,9 +356,10 @@ public class MCS
 		//
 		// The substructure has to be searched in the molecule because the mapping indices are not contained in the
 		// IntVec that is the solution from the MCS search.
+		arrMatchListFrag2Mol = null;
 		if(sss.findFragmentInMolecule(SSSearcher.cCountModeOverlapping, SSSearcher.cDefaultMatchMode,excluded)>0){
 			ArrayList<int[]> liMatchSubFrag2Mol = sss.getMatchList();
-			int [] arrMatchListFrag2Mol = getMappedMatchListFrag2Mol(fragSub, liMatchSubFrag2Mol.get(0));
+			arrMatchListFrag2Mol = getMappedMatchListFrag2Mol(fragSub, liMatchSubFrag2Mol.get(0));
 			int bondsMol = mol.getBonds();
 			if(arrBondMCSMol==null) {
 				arrBondMCSMol = new boolean [bondsMol];
@@ -382,7 +373,16 @@ public class MCS
 		return arrBondMol_Result;
 	}
 
-    private int[] getMappedMatchListFrag2Mol(StereoMolecule fragSub, int[] arrMatchListSubFrag2Mol)
+	/**
+	 * delivers the match list indices from getMCSBondArray(...)
+	 * index is the Fragment atom index, value is the matched atom index in molecule.
+	 * @return null if no match was found.
+	 */
+	public int[] getArrMatchListFrag2Mol() {
+		return arrMatchListFrag2Mol;
+	}
+
+	private int[] getMappedMatchListFrag2Mol(StereoMolecule fragSub, int[] arrMatchListSubFrag2Mol)
     {
 		int [] arrMatchListFrag2Mol = new int [frag.getAtoms()];
 //        SSSearcher sss = new SSSearcher();


=====================================
src/main/java/com/actelion/research/chem/mcs/SubStructSearchBondIndex.java deleted
=====================================
@@ -1,68 +0,0 @@
-package com.actelion.research.chem.mcs;
-
-import com.actelion.research.chem.StereoMolecule;
-import com.actelion.research.util.datamodel.IntVec;
-
-/*
-* 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.
-*
-*/
-public class SubStructSearchBondIndex {
-	
-	
-	private StereoMolecule mol;
-	
-	private StereoMolecule frag;
-	
-	private IntVec ivBondIndex;
-	
-	public void setMolecule(StereoMolecule mol){
-		this.mol=mol;
-	}
-	
-	public void setFragment(StereoMolecule frag){
-		this.frag=frag;
-	}
-	
-	public void setBondIndex(IntVec ivBondIndex){
-		this.ivBondIndex=ivBondIndex;
-		
-		
-	}
-	
-	public int findFragmentInMolecule() {
-		
-		
-		
-		return -1;
-	}
-
-}


=====================================
src/main/java/com/actelion/research/chem/mcs/SubStructSearchExhaustiveTreeWalker.java
=====================================
@@ -281,7 +281,7 @@ public class SubStructSearchExhaustiveTreeWalker {
 	
 	/**
 	 * Generates a tree from the fragment
-	 * @param possMapFrag2MolStart
+	 * @param
 	 * @return an array with one field for each fragment atom. 
 	 */
 	private void determineProcessingArray(){


=====================================
src/main/java/com/actelion/research/chem/phesa/pharmacophore/PharmacophoreCalculator.java
=====================================
@@ -1,10 +1,5 @@
 package com.actelion.research.chem.phesa.pharmacophore;
 
-import java.util.ArrayList;
-import java.util.Arrays;
-import java.util.List;
-import java.util.stream.Collectors;
-
 import com.actelion.research.chem.AtomFunctionAnalyzer;
 import com.actelion.research.chem.RingCollection;
 import com.actelion.research.chem.StereoMolecule;
@@ -12,9 +7,13 @@ import com.actelion.research.chem.interactionstatistics.InteractionAtomTypeCalcu
 import com.actelion.research.chem.phesa.pharmacophore.pp.AcceptorPoint;
 import com.actelion.research.chem.phesa.pharmacophore.pp.AromRingPoint;
 import com.actelion.research.chem.phesa.pharmacophore.pp.DonorPoint;
-import com.actelion.research.chem.phesa.pharmacophore.pp.ExitVectorPoint;
 import com.actelion.research.chem.phesa.pharmacophore.pp.IPharmacophorePoint;
 
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.List;
+import java.util.stream.Collectors;
+
 public class PharmacophoreCalculator {
 	
 	public static final int ACCEPTOR_ID = 0;
@@ -79,36 +78,36 @@ public class PharmacophoreCalculator {
 					}
 					if(mol.getAtomicNo(i)==8 && neighbours==1 && (mol.getConnBondOrder(i, 0)==2 || AtomFunctionAnalyzer.isAcidicOxygen(mol, i) || mol.getAtomCharge(i)==-1 )) {
 						int a1 = mol.getConnAtom(i,0);
-						if(!(mol.getAtomicNo(a1)==16 || mol.getAtomicNo(a1)==15)) {		
-							int aa1 = mol.getConnAtom(a1,0);
-							if(aa1==i) 
-								aa1 = mol.getConnAtom(a1,1);
-							neighbourList.add(aa1);
-							AcceptorPoint ap = new AcceptorPoint(mol,i,neighbourList,interactionClass,1);
+						if (mol.getConnAtoms(a1) > 1) { // added this check to prevent OOB exceptions with cropped proteins; TLS 17Oct2021
+							if(!(mol.getAtomicNo(a1)==16 || mol.getAtomicNo(a1)==15)) {
+								int aa1 = mol.getConnAtom(a1,0);
+								if(aa1==i)
+									aa1 = mol.getConnAtom(a1,1);
+								neighbourList.add(aa1);
+								AcceptorPoint ap = new AcceptorPoint(mol,i,neighbourList,interactionClass,1);
 
-							ppPoints.add(ap);
-							List<Integer> neighbourList2 = new ArrayList<Integer>();
-							for(int neighbour : neighbourList) {
-								neighbourList2.add(neighbour);
-							}
-							ap = new AcceptorPoint(mol,i,neighbourList2,interactionClass,2);
-							ppPoints.add(ap);
+								ppPoints.add(ap);
+								List<Integer> neighbourList2 = new ArrayList<Integer>();
+								for(int neighbour : neighbourList) {
+									neighbourList2.add(neighbour);
+								}
+								ap = new AcceptorPoint(mol,i,neighbourList2,interactionClass,2);
+								ppPoints.add(ap);
 
+								}
+							else {
+								AcceptorPoint ap = new AcceptorPoint(mol,i,neighbourList,interactionClass);
+								ppPoints.add(ap);
 							}
-						else { 
-							AcceptorPoint ap = new AcceptorPoint(mol,i,neighbourList,interactionClass);
-							ppPoints.add(ap);
 						}
-	
 					}
-					else {
-					AcceptorPoint ap = new AcceptorPoint(mol,i,neighbourList,interactionClass);
-					ppPoints.add(ap);
-
+					else if (neighbourList.size() != 0) { // added this check to prevent OOB exceptions with cropped proteins; TLS 17Oct2021
+						AcceptorPoint ap = new AcceptorPoint(mol,i,neighbourList,interactionClass);
+						ppPoints.add(ap);
 					}
+				}
 			}
 		}
-		}
 		return ppPoints;
 	}
 	
@@ -124,7 +123,7 @@ public class PharmacophoreCalculator {
 				}
 			}
 			else if (mol.getAtomicNo(a)==7){ // atom is not aromatic
-				if(mol.getConnBondOrder(a, 0)==3) //nitrile
+				if(mol.getConnAtoms(a) == 1 && mol.getConnBondOrder(a, 0)==3) //nitrile
 					return true;
 				if (mol.isFlatNitrogen(a)) { 
 					for(int b=0;b<mol.getConnAtoms(a);b++) {


=====================================
src/main/java/com/actelion/research/chem/phesaflex/FlexibleShapeAlignment.java
=====================================
@@ -7,7 +7,7 @@ import java.util.Map;
 import com.actelion.research.chem.StereoMolecule;
 import com.actelion.research.chem.alignment3d.transformation.TransformationSequence;
 import com.actelion.research.chem.forcefield.mmff.ForceFieldMMFF94;
-import com.actelion.research.chem.forcefield.mmff.PositionConstraint;
+import com.actelion.research.chem.forcefield.mmff.MMFFPositionConstraint;
 import com.actelion.research.chem.optimization.OptimizerLBFGS;
 import com.actelion.research.chem.phesa.MolecularVolume;
 import com.actelion.research.chem.phesa.PheSAAlignment;
@@ -166,7 +166,7 @@ public class FlexibleShapeAlignment {
 		int cycles = 0;
 		while(notRelaxed && cycles<maxCycles) {
 			ForceFieldMMFF94 forceField = new ForceFieldMMFF94(fitMol, ForceFieldMMFF94.MMFF94SPLUS, ffOptions);
-			PositionConstraint constraint = new PositionConstraint(fitMol,50,init);
+			MMFFPositionConstraint constraint = new MMFFPositionConstraint(fitMol,50,init);
 			forceField.addEnergyTerm(constraint);
 			forceField.minimise();
 			double e = forceField.getTotalEnergy();


=====================================
src/main/java/com/actelion/research/chem/potentialenergy/PositionConstraint.java
=====================================
@@ -15,9 +15,9 @@ public class PositionConstraint implements PotentialEnergyTerm {
 		refPos = new double[3];
 		this.atom = atom;
 		this.conf = conf;
-		refPos[0] = conf.getMolecule().getAtomX(atom);
-		refPos[1] = conf.getMolecule().getAtomY(atom);
-		refPos[2] = conf.getMolecule().getAtomZ(atom);
+		refPos[0] = conf.getX(atom);
+		refPos[1] = conf.getY(atom);
+		refPos[2] = conf.getZ(atom);
 		
 		this.k = k;
 		this.d = d;



View it on GitLab: https://salsa.debian.org/java-team/openchemlib/-/compare/c5ed69cf9116b44e98bb5bead30d862756a07422...0856cc213369aa602d497dc60712407421cedc2d

-- 
View it on GitLab: https://salsa.debian.org/java-team/openchemlib/-/compare/c5ed69cf9116b44e98bb5bead30d862756a07422...0856cc213369aa602d497dc60712407421cedc2d
You're receiving this email because of your account on salsa.debian.org.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-java-commits/attachments/20211026/887e4424/attachment.htm>


More information about the pkg-java-commits mailing list