FragmenterTools.java

Go to the documentation of this file.

package denoptim.fragmenter;
 
import java.io.File;
import java.io.FileInputStream;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.logging.Level;
import java.util.logging.Logger;
 
import javax.vecmath.Point3d;
import javax.vecmath.Vector3d;
 
import org.openscience.cdk.Atom;
import org.openscience.cdk.Bond;
import org.openscience.cdk.DefaultChemObjectBuilder;
import org.openscience.cdk.PseudoAtom;
import org.openscience.cdk.config.Isotopes;
import org.openscience.cdk.exception.CDKException;
import org.openscience.cdk.geometry.alignment.KabschAlignment;
import org.openscience.cdk.interfaces.IAtom;
import org.openscience.cdk.interfaces.IAtomContainer;
import org.openscience.cdk.interfaces.IBond;
import org.openscience.cdk.interfaces.IIsotope;
import org.openscience.cdk.io.iterator.IteratingSDFReader;
import org.openscience.cdk.isomorphism.Mappings;
import org.openscience.cdk.isomorphism.Pattern;
import org.openscience.cdk.silent.SilentChemObjectBuilder;
import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
 
import denoptim.constants.DENOPTIMConstants;
import denoptim.exception.DENOPTIMException;
import denoptim.files.FileFormat;
import denoptim.files.UndetectedFileFormatException;
import denoptim.graph.APClass;
import denoptim.graph.AttachmentPoint;
import denoptim.graph.DGraph;
import denoptim.graph.Edge;
import denoptim.graph.Fragment;
import denoptim.graph.Template;
import denoptim.graph.Vertex;
import denoptim.graph.Vertex.BBType;
import denoptim.graph.rings.RingClosingAttractor;
import denoptim.io.DenoptimIO;
import denoptim.io.IteratingAtomContainerReader;
import denoptim.molecularmodeling.ThreeDimTreeBuilder;
import denoptim.programs.RunTimeParameters.ParametersType;
import denoptim.programs.fragmenter.CuttingRule;
import denoptim.programs.fragmenter.FragmenterParameters;
import denoptim.programs.fragmenter.MatchedBond;
import denoptim.utils.CartesianSpaceUtils;
import denoptim.utils.DummyAtomHandler;
import denoptim.utils.FormulaUtils;
import denoptim.utils.ManySMARTSQuery;
import denoptim.utils.MathUtils;
import denoptim.utils.MoleculeUtils;
import denoptim.utils.Randomizer;
 
public class FragmenterTools
{
    
//------------------------------------------------------------------------------
    
    public static void checkElementalAnalysisAgainstFormula(File input,
            File output, Logger logger) 
                    throws DENOPTIMException, IOException
    {
        FileInputStream fis = new FileInputStream(input);
        IteratingSDFReader reader = new IteratingSDFReader(fis, 
                DefaultChemObjectBuilder.getInstance());
 
        int index = -1;
        int maxBufferSize = 2000;
        ArrayList<IAtomContainer> buffer = new ArrayList<IAtomContainer>(500);
        try {
            while (reader.hasNext())
            {
                index++;
                if (logger!=null)
                {
                    logger.log(Level.FINE,"Checking elemental analysis of "
                            + "structure " + index);
                }
                IAtomContainer mol = reader.next();
                if (mol.getProperty(DENOPTIMConstants.FORMULASTR)==null)
                {
                    throw new Error("Property '" + DENOPTIMConstants.FORMULASTR 
                            + "' not found in molecule " + index + " in file "
                            + input + ". Cannot compare formula with elemental"
                            + "analysis.");
                }
                String formula = mol.getProperty(DENOPTIMConstants.FORMULASTR)
                        .toString();
                
                if (FormulaUtils.compareFormulaAndElementalAnalysis(formula, 
                        mol, logger))
                {
                    buffer.add(mol);
                } else {
                    if (logger!=null)
                    {
                        logger.log(Level.INFO,"Inconsistency between elemental "
                                + "analysis of structure and molecular formula."
                                + " Rejecting structure " + index + ": " 
                                + mol.getTitle());
                    }
                }
                
                // If max buffer size is reached, then bump to file
                if (buffer.size() >= maxBufferSize)
                {
                    DenoptimIO.writeSDFFile(output.getAbsolutePath(), buffer, 
                            true);
                    buffer.clear();
                }
            }
        }
        finally {
            reader.close();
        }
        if (buffer.size() < maxBufferSize)
        {
            DenoptimIO.writeSDFFile(output.getAbsolutePath(), buffer, true);
            buffer.clear();
        }
    }
    
//------------------------------------------------------------------------------
    
 
    public static boolean prepareMolToFragmentation(IAtomContainer mol, 
            FragmenterParameters settings, int index)
    {
        try
        {
            if (settings.addExplicitH())
            {
                MoleculeUtils.explicitHydrogens(mol);
            } else {
                MoleculeUtils.setZeroImplicitHydrogensToAllAtoms(mol);
            }
            MoleculeUtils.ensureNoUnsetBondOrders(mol);
        } catch (CDKException e)
        {
            if (e.getMessage().contains("Cannot assign Kekulé structure"))
            {
                if (!settings.acceptUnsetToSingeBO())
                {
                    settings.getLogger().log(Level.WARNING,"Some bond order "
                            + "are unset and attempt to kekulize the "
                            + "system has failed "
                            + "for structure " + index + ". "
                            + "This hampers use of SMARTS queries, which "
                            + "may very well "
                            + "not work as expected. Structure " + index 
                            + " will be rejected. "
                            + "You can avoid rejection by using "
                            + "keyword " 
                            + ParametersType.FRG_PARAMS.getKeywordRoot() 
                            + "UNSETTOSINGLEBO, but you'll "
                            + "still be using a peculiar connectivity "
                            + "table were "
                            + "many bonds are artificially marked as "
                            + "single to "
                            + "avoid use of 'UNSET' bond order. "
                            + "Further details on the problem: " 
                            + e.getMessage());
                    return false;
                } else {
                    settings.getLogger().log(Level.WARNING,"Failed "
                            + "kekulization "
                            + "for structure " + index 
                            + " but UNSETTOSINGLEBO "
                            + "keyword used. Forcing use of single bonds to "
                            + "replace bonds with unset order.");
                    for (IBond bnd : mol.bonds())
                    {
                        if (bnd.getOrder().equals(IBond.Order.UNSET)) 
                        {
                            bnd.setOrder(IBond.Order.SINGLE);
                        }
                    }
                }
            }
        }
        return true;
    }
    
//------------------------------------------------------------------------------
    
    public static void filterStrucutresBySMARTS(File input, Set<String> smarts,
            File output, Logger logger) 
                    throws DENOPTIMException, IOException
    {
        FileInputStream fis = new FileInputStream(input);
        IteratingSDFReader reader = new IteratingSDFReader(fis, 
                DefaultChemObjectBuilder.getInstance());
 
        int i = -1;
        Map<String, String> smartsMap = new HashMap<String, String>();
        for (String s : smarts)
        {
            i++;
            smartsMap.put("prefilter-"+i, s);
        }
        
        int index = -1;
        int maxBufferSize = 2000;
        ArrayList<IAtomContainer> buffer = new ArrayList<IAtomContainer>(500);
        try {
            while (reader.hasNext())
            {
                index++;
                if (logger!=null)
                {
                    logger.log(Level.FINE,"Prefiltering structure " + index);
                }
                IAtomContainer mol = reader.next();
                
                ManySMARTSQuery msq = new ManySMARTSQuery(mol, smartsMap);
                if (msq.hasProblems())
                {
                    String msg = "WARNING! Problems while searching for "
                            + "specific atoms/bonds using SMARTS: " 
                            + msq.getMessage();
                    throw new DENOPTIMException(msg,msq.getProblem());
                }
                Map<String, Mappings> allMatches = msq.getAllMatches();
                
                if (allMatches.size()==0)
                {
                    buffer.add(mol);
                } else {
                    String hits = "";
                    for (String s : allMatches.keySet())
                        hits = hits + DenoptimIO.NL + smartsMap.get(s);
                    if (logger!=null)
                    {
                        logger.log(Level.INFO,"Found match for " + hits
                                + "Rejecting structure " + index + ": " 
                                + mol.getTitle());
                    }
                }
                
                // If max buffer size is reached, then bump to file
                if (buffer.size() >= maxBufferSize)
                {
                    DenoptimIO.writeSDFFile(output.getAbsolutePath(), buffer, 
                            true);
                    buffer.clear();
                }
            }
        } finally {
            reader.close();
        }
        if (buffer.size() < maxBufferSize)
        {
            DenoptimIO.writeSDFFile(output.getAbsolutePath(), buffer, true);
            buffer.clear();
        }
    }
 
//------------------------------------------------------------------------------
 
    public static boolean fragmentationFromGraphs(File input, 
        FragmenterParameters settings, File output, Logger logger) 
        throws DENOPTIMException, Exception
    {
        int totalProd = 0;
        for (DGraph graph : DenoptimIO.readDENOPTIMGraphsFromFile(input))
        {
            List<Vertex> fragments = graph.getVertexList();
 
            // Post-fragmentation processing of fragments
            List<Vertex> keptFragments = new ArrayList<Vertex>();
            int fragCounter = 0;
            for (Vertex frag : fragments)
            {
                // Add metadata
                fragCounter++;
                manageFragmentCollection(frag, fragCounter, settings,
                        keptFragments, logger);
            }
            if (logger!=null)
            {
                logger.log(Level.FINE,"Fragments surviving post-"
                        + "processing: " + keptFragments.size());
            }
 
            if (keptFragments.size()>0)
            {
                totalProd += keptFragments.size();
                DenoptimIO.writeVertexesToFile(output, FileFormat.VRTXSDF, 
                    keptFragments, true);
            }
        }
        return totalProd>0;
    }
 
//------------------------------------------------------------------------------
 
    public static List<Vertex> fragmentation(IAtomContainer mol, 
            FragmenterParameters settings) 
            throws DENOPTIMException
    {
        List<Vertex> fragments = new ArrayList<Vertex>();
        if (settings.getFragmentationTmpls().size()>0)
        {
            fragments = fragmentation(mol, settings.getFragmentationTmpls(), 
                settings.getMaxBufferShellSize(),
                settings.getRandomizer(), settings.getLogger());
        } else {
            fragments = fragmentation(mol, settings.getCuttingRules(), 
                settings.getLogger());
        }
        return fragments;
    }
 
//-----------------------------------------------------------------------------
    
    public static boolean fragmentation(File input, FragmenterParameters settings,
            File output, Logger logger) throws CDKException, IOException, 
    DENOPTIMException, IllegalArgumentException, UndetectedFileFormatException
    {
        IteratingAtomContainerReader iterator = 
                new IteratingAtomContainerReader(input);
 
        int totalProd = 0;
        int totalKept = 0;
        int index = -1;
        try {
            while (iterator.hasNext())
            {
                index++;
                if (logger!=null)
                {
                    logger.log(Level.FINE,"Fragmenting structure " + index);
                }
                IAtomContainer mol = iterator.next();
                String molName = "noname-mol" + index;
                if (mol.getTitle()!=null && !mol.getTitle().isBlank())
                    molName = mol.getTitle();
                
                // Generate the fragments
                List<Vertex> fragments = fragmentation(mol, settings);
                if (logger!=null)
                {
                    logger.log(Level.FINE,"Fragmentation produced " 
                            + fragments.size() + " fragments.");
                }
                totalProd += fragments.size();
                
                // Post-fragmentation processing of fragments
                List<Vertex> keptFragments = new ArrayList<Vertex>();
                int fragCounter = 0;
                for (Vertex frag : fragments)
                {
                    // Add metadata
                    String fragIdStr = "From_" + molName + "_" + fragCounter;
                    frag.setProperty("cdk:Title", fragIdStr);
                    fragCounter++;
                    manageFragmentCollection(frag, fragCounter, settings,
                            keptFragments, logger);
                }
                if (logger!=null)
                {
                    logger.log(Level.FINE,"Fragments surviving post-"
                            + "processing: " + keptFragments.size());
                }
                totalKept += keptFragments.size();
                if (!settings.doManageIsomorphicFamilies() && totalKept>0)
                {
                    DenoptimIO.writeVertexesToFile(output, FileFormat.VRTXSDF, 
                            keptFragments,true);
                }
            }
        } finally {
            iterator.close();
        }
        
        // Did we actually produce anything? We might not...
        if (totalProd==0)
        {
            if (logger!=null)
            {
                logger.log(Level.WARNING,"No fragment produced. Cutting rules "
                        + "were ineffective on the given structures.");
            }
            return false;
        } else if (totalKept==0)
        {
            if (logger!=null)
            {
                logger.log(Level.WARNING,"No fragment kept out of " + totalProd 
                        + " produced fragments. Filtering criteria might be "
                        + "too restrictive.");
            }
            return false;
        }
        return true;
    }
    
//------------------------------------------------------------------------------
    
    public static List<Vertex> fragmentation(IAtomContainer mol, 
            List<DGraph> templates, int maxBufferShellSize, 
            Randomizer randomizer, Logger logger) 
            throws DENOPTIMException
    { 
        // The best set of fragments is the one with as many fragments that are isomorfic 
        // // with the template
        List<Vertex> fragments = new ArrayList<Vertex>();
        int maxNumIsomorfFrgs = -1;
        for (DGraph templateGraph : templates)
        {
            // Check if the template is connected
            if (!templateGraph.isConnected())
            {
                // This limitation derives from the fact that to speed up substructure
                // search we reduce the template graph to the smallest set of atoms that
                // identify the topology of APs. The 
                // can identify the matching topology. This is not possible for disconnected
                // graphs.
                throw new DENOPTIMException("Template graph is not connected. "
                        + "Cannot use it for fragmentation. Please check the "
                        + "template graph and make sure it is connected.");
            }
 
            ThreeDimTreeBuilder tb = new ThreeDimTreeBuilder(logger, randomizer);
            IAtomContainer templateMol = tb.convertGraphTo3DAtomContainer(templateGraph, 
                true, true, true);
 
            //TODO: consider splitting the template into connected components and
            // reduce each component independently.
            
            // Optimized creation of a reduced templateMol for faster isomorphism search
            // Keep full templateMol for exploreDGraphForMappings
            TopoTemplateProducer ttp = new TopoTemplateProducer(templateMol);
            
            // The best set of fragments from this template may not be perfect, so we
            // consider a score that consists of how many vertexes are isomorphic
            // to the template graph vertexes.
            List<Vertex> fragsFromThisTemplate = new ArrayList<>();
            int maxNumIsomorfFrgsFromThisTmpl = -1;
 
            // Define the maximum theoretical score considering that the template may
            // be matched multiple times in the molecule
            List<Vertex> templateVertexes = templateGraph.getVertexAnyLevel();
            int maxTheoreticalScorePerMapping = 0;
            for (int iva=0; iva<templateVertexes.size(); iva++)
            {
                Vertex va = templateVertexes.get(iva);
                if (!(va instanceof Fragment))
                {
                    continue;
                }
                for (int ivb=0; ivb<templateVertexes.size(); ivb++)
                {
                    Vertex vb = templateVertexes.get(ivb);
                    if (!(vb instanceof Fragment))
                    {
                        continue;
                    }
                    if (iva == ivb)
                    {
                        // we know that the same
                        maxTheoreticalScorePerMapping++;
                    } else if (((Fragment) va).isIsomorphicTo(vb))
                    {
                        maxTheoreticalScorePerMapping++;
                    }
                }
            }
 
            // Use the smallest possible template needed to identify the graph topology
            // in the molecular structure.
            for (int bufferShellSize = 0; bufferShellSize < maxBufferShellSize; bufferShellSize++)
            {
                // Use a reduced templateMol for isomorphism search (faster)
                IAtomContainer reducedTemplateMol = ttp.getTemplateWithBufferShell(bufferShellSize);
                List<Map<IAtom,IAtom>> reducedAtomMappings = MoleculeUtils.findUniqueAtomMappings(
                        reducedTemplateMol, mol, logger);
                
                // Convert mappings from reducedTemplateMol:mol to templateMol:mol using stored indices
                List<Map<IAtom,IAtom>> atomMappings = new ArrayList<>();
                for (Map<IAtom,IAtom> reducedMapping : reducedAtomMappings)
                {
                    Map<IAtom,IAtom> fullMapping = new HashMap<>();
                    for (Map.Entry<IAtom,IAtom> entry : reducedMapping.entrySet())
                    {
                        IAtom reducedAtom = entry.getKey();
                        IAtom molAtom = entry.getValue();
                        
                        // Get original atom index from reduced atom property
                        Object indexObj = reducedAtom.getProperty("DENOPTIM_ORIGINAL_ATOM_INDEX");
                        if (indexObj != null)
                        {
                            int originalIndex = ((Number) indexObj).intValue();
                            IAtom originalAtom = templateMol.getAtom(originalIndex);
                            if (originalAtom != null)
                            {
                                fullMapping.put(originalAtom, molAtom);
                            }
                        }
                    }
                    if (!fullMapping.isEmpty())
                    {
                        atomMappings.add(fullMapping);
                    }
                }
                
                Fragment masterFrag = new Fragment(mol, BBType.UNDEFINED);
                IAtomContainer masterFragIAC = masterFrag.getIAtomContainer();
 
                // For each atom mapping, replace bonds corresponding to edges 
                // (at any embedding level) with attachment points and annotate embedding level.
                // In case of symmetry, multiple mappings will effectivly cut the same bonds,
                // and produce the same fragments of the template graph only for
                // symmetry-redundant mappings.
                for (Map<IAtom,IAtom> atomMapping : atomMappings)
                {
                    try {
                        exploreDGraphForMappings(templateGraph, templateMol,
                            masterFrag, masterFragIAC, mol, atomMapping);
                    } catch (Throwable e) {
                        e.printStackTrace();
                        logger.log(Level.WARNING, "Error while exploring the template graph: " + e.getMessage());
                        continue;
                    }
                }
 
                // Extract isolated fragments
                List<Vertex> locfragments = new ArrayList<Vertex>();
                Set<Integer> doneAlready = new HashSet<Integer>();
                for (int idx=0 ; idx<masterFrag.getAtomCount(); idx++)
                {
                    if (doneAlready.contains(idx))
                        continue;
                    
                    Fragment cloneOfMaster = masterFrag.clone();
                    IAtomContainer iac = cloneOfMaster.getIAtomContainer();
                    Set<IAtom> atmsToKeep = exploreConnectivity(iac.getAtom(idx), iac);
                    atmsToKeep.stream().forEach(atm -> doneAlready.add(iac.indexOf(atm)));
                    
                    Set<IAtom> atmsToRemove = new HashSet<IAtom>();
                    for (IAtom atm : cloneOfMaster.atoms())
                    {
                        if (!atmsToKeep.contains(atm))
                        {
                            atmsToRemove.add(atm);
                        }
                    }
                    cloneOfMaster.removeAtoms(atmsToRemove);
                    if (cloneOfMaster.getAttachmentPoints().size()>0)
                        locfragments.add(cloneOfMaster);
                }
            
                int localScore = 0;
                for (Vertex frag : locfragments)
                {
                    if (!(frag instanceof Fragment))
                    {
                        continue;
                    }
                    for (Vertex templateVertex : templateGraph.getVertexList())
                    {
                        if (((Fragment) frag).isIsomorphicTo(templateVertex))
                        {
                            localScore++;
                        }
                    }
                }
 
                if (localScore > maxNumIsomorfFrgsFromThisTmpl)
                {
                    maxNumIsomorfFrgsFromThisTmpl = localScore;
                    fragsFromThisTemplate = locfragments;
                    if (maxNumIsomorfFrgsFromThisTmpl == 
                        maxTheoreticalScorePerMapping*atomMappings.size())
                    {
                        break;
                    }
                }
            }
 
            if (maxNumIsomorfFrgsFromThisTmpl > maxNumIsomorfFrgs)
            {
                maxNumIsomorfFrgs = maxNumIsomorfFrgsFromThisTmpl;
                fragments = fragsFromThisTemplate;
            }
        }
        return fragments;
    }
 
//------------------------------------------------------------------------------
 
    private static void exploreDGraphForMappings(DGraph graph, IAtomContainer graphIAC, 
        Fragment masterFrag, IAtomContainer masterFragIAC, IAtomContainer mol,
        Map<IAtom,IAtom> graphToMolMapping) 
        throws DENOPTIMException
    {
        int cutId = -1;
        for (Edge edge : graph.getEdgeList())
        {
            cutId++;
            int srcAPIdx = edge.getSrcAP().getAtomPositionNumberInMol();
            int trgAPIdx = edge.getTrgAP().getAtomPositionNumberInMol();
            IAtom graphAtmSrc = graphIAC.getAtom(srcAPIdx);
            IAtom graphAtmTrg = graphIAC.getAtom(trgAPIdx);
 
            IAtom masterFragAtmSrc = masterFragIAC.getAtom(mol.indexOf(graphToMolMapping.get(graphAtmSrc)));
            IAtom masterFragAtmTrg = masterFragIAC.getAtom(mol.indexOf(graphToMolMapping.get(graphAtmTrg)));
            
            IBond bnd = masterFragIAC.getBond(masterFragAtmSrc, masterFragAtmTrg);
            if (bnd != null)
            {
                masterFragIAC.removeBond(bnd);
            }
 
            AttachmentPoint srcAP = masterFrag.addAPOnAtom(masterFragAtmSrc, 
                edge.getSrcAP().getAPClass(), 
                MoleculeUtils.getPoint3d(masterFragAtmTrg));
            srcAP.setCutId(cutId);
            AttachmentPoint trgAP = masterFrag.addAPOnAtom(masterFragAtmTrg, 
                edge.getTrgAP().getAPClass(), 
                MoleculeUtils.getPoint3d(masterFragAtmSrc));
            trgAP.setCutId(cutId);
        }
 
        for (Vertex v : graph.getVertexList())
        {
            if (v instanceof Template)
            {
                DGraph innerGraph = ((Template) v).getInnerGraph();
                exploreDGraphForMappings(innerGraph, graphIAC, masterFrag, masterFragIAC, mol, graphToMolMapping);
                continue;
            }
 
            // Deal with vertexes without edges, i.e., free APs
            for (AttachmentPoint ap : v.getAttachmentPoints())
            {
                if (!ap.isAvailable())
                {
                    continue;
                }
                IAtom graphAtmSrc = graphIAC.getAtom(ap.getAtomPositionNumberInMol());
                IAtom masterFragAtmSrc = masterFragIAC.getAtom(mol.indexOf(graphToMolMapping.get(graphAtmSrc)));
 
                // find the position of the AP head in 3D space by aligning geoetry to template geometry
                Point3d apHead = findPointAlignedWithTmpl(ap, graphAtmSrc, 
                    graphIAC, masterFragIAC, mol, graphToMolMapping);
 
                // Make the free AP on the master
                masterFrag.addAPOnAtom(masterFragAtmSrc, ap.getAPClass(), apHead);
            }
        }
    }
 
//------------------------------------------------------------------------------
 
    private static Point3d findPointAlignedWithTmpl(AttachmentPoint ap, 
        IAtom graphAtmSrc, IAtomContainer graphIAC, IAtomContainer masterFragIAC, 
        IAtomContainer mol,
        Map<IAtom,IAtom> graphToMolMapping) throws DENOPTIMException
    {
        // Ensure assumptions
        if (!graphIAC.contains(graphAtmSrc))
        {
            throw new DENOPTIMException("The atom " + graphAtmSrc.getSymbol() + " is not in the graphIAC.");
        }
        if (mol.getAtomCount() != masterFragIAC.getAtomCount())
        {
            throw new DENOPTIMException("The number of atoms in the molecule and the master fragment are not the same.");
        }
 
        // Work with cloned molecules
        IAtomContainer molA = MoleculeUtils.makeSameAs(masterFragIAC);
        IAtomContainer molB = MoleculeUtils.makeSameAs(graphIAC);
 
        // Analogue of graphAtmSrc (belongs to graphIAC) on molA, which is a clone of masterFragIAC
        // which is consistent with mol.
        IAtom graphAtmSrcOnMolA = molA.getAtom(mol.indexOf(graphToMolMapping.get(graphAtmSrc)));
 
        // Add the AP as a dummy atom to the molecule that will be rototranslated
        IAtom dummyAtm = new Atom("He"); //element is irrelevant
        dummyAtm.setPoint3d(ap.getDirectionVector());
        molB.addAtom(dummyAtm);
 
        List<IAtom> closestNeighborsOnGrphIAC = graphIAC.getConnectedAtomsList(graphAtmSrc);
        List<IAtom> closestNeighborsOnMolA = molA.getConnectedAtomsList(graphAtmSrcOnMolA);
        int minNumClosestNeighbors = Math.min(closestNeighborsOnGrphIAC.size(), closestNeighborsOnMolA.size());
        if (minNumClosestNeighbors < 1)
        {
            // We do not have enough atoms to do an alignement: use only distance
            // NB: the case of a biatomic fragment is handles below.
            double dist = graphAtmSrc.getPoint3d().distance(ap.getDirectionVector());
            Point3d apHead = new Point3d(
                graphAtmSrcOnMolA.getPoint3d().x + dist,
                graphAtmSrcOnMolA.getPoint3d().y,
                graphAtmSrcOnMolA.getPoint3d().z);
            return apHead;
        } else if (minNumClosestNeighbors == 1) {
            // For biatomic fragments
            double dist = graphAtmSrc.getPoint3d().distance(ap.getDirectionVector());
            double angle = MathUtils.angle(
                closestNeighborsOnGrphIAC.get(0).getPoint3d(), 
                graphAtmSrc.getPoint3d(), 
                ap.getDirectionVector());
            IAtom nbrAtmOnMolA = molA.getAtom(mol.indexOf(graphToMolMapping.get(closestNeighborsOnGrphIAC.get(0))));
            Point3d pSrc = graphAtmSrcOnMolA.getPoint3d();
            Point3d pNbr = nbrAtmOnMolA.getPoint3d();
            Point3d pSrcTmpl = graphAtmSrc.getPoint3d();
            Point3d pNbrTmpl = closestNeighborsOnGrphIAC.get(0).getPoint3d();
            Point3d pApTmpl = ap.getDirectionVector();
 
            Vector3d uMol = CartesianSpaceUtils.getVectorFromTo(pSrc, pNbr);
            uMol.normalize();
            Vector3d uTmpl = CartesianSpaceUtils.getVectorFromTo(pSrcTmpl, pNbrTmpl);
            uTmpl.normalize();
            Vector3d wTmpl = CartesianSpaceUtils.getVectorFromTo(pSrcTmpl, pApTmpl);
            wTmpl.normalize();
 
            double angleRad = Math.toRadians(angle);
            Vector3d apDir = new Vector3d();
            apDir.scale(Math.cos(angleRad), uMol);
            double sinAngle = Math.sin(angleRad);
            if (sinAngle > DENOPTIMConstants.FLOATCOMPARISONTOLERANCE)
            {
                Vector3d vRefTmpl = new Vector3d(wTmpl);
                vRefTmpl.scaleAdd(-uTmpl.dot(wTmpl), uTmpl, vRefTmpl);
                vRefTmpl.normalize();
                Vector3d rotAxis = new Vector3d();
                rotAxis.cross(uTmpl, uMol);
                double sinRot = rotAxis.length();
                double cosRot = uTmpl.dot(uMol);
                if (sinRot > DENOPTIMConstants.FLOATCOMPARISONTOLERANCE)
                {
                    rotAxis.normalize();
                    double rotAng = Math.toDegrees(Math.atan2(sinRot, cosRot));
                    CartesianSpaceUtils.rotatedVectorWAxisAngle(vRefTmpl, rotAxis, rotAng);
                }
                else if (cosRot < 0.0)
                {
                    Vector3d flipAxis = CartesianSpaceUtils.getNormalDirection(uMol);
                    CartesianSpaceUtils.rotatedVectorWAxisAngle(vRefTmpl, flipAxis, 180.0);
                }
                vRefTmpl.scale(sinAngle);
                apDir.add(vRefTmpl);
            }
 
            Point3d apHead = new Point3d(
                pSrc.x + dist * apDir.x,
                pSrc.y + dist * apDir.y,
                pSrc.z + dist * apDir.z);
            return apHead;
        } else {
            // Try to alogn atoms around AP sourc to get the position of the AP head
            Map<IAtom,IAtom> molToGraphMappingAroundAPSrc = new HashMap<>();
            Set<IAtom> atomsAlreadyUsed = new HashSet<>();
            for (Map.Entry<IAtom,IAtom> pair : graphToMolMapping.entrySet())
            {
                IAtom atmGraphIAC = pair.getKey();
                if (graphIAC.getConnectedAtomsList(graphAtmSrc).contains(atmGraphIAC)
                    || atmGraphIAC == graphAtmSrc)
                {
                    atomsAlreadyUsed.add(atmGraphIAC);
                    molToGraphMappingAroundAPSrc.put(
                        molA.getAtom(mol.indexOf(pair.getValue())),
                        molB.getAtom(graphIAC.indexOf(atmGraphIAC)));
                }
            }
            // When we have few atoms we take also the second shell of neighbors
            if (minNumClosestNeighbors<3)
            {
                for (Map.Entry<IAtom,IAtom> pair : graphToMolMapping.entrySet())
                {
                    IAtom atmGraphIAC = pair.getKey();
                    if (closestNeighborsOnGrphIAC.contains(atmGraphIAC) 
                        && !atomsAlreadyUsed.contains(atmGraphIAC))
                    {
                        atomsAlreadyUsed.add(atmGraphIAC);
                        molToGraphMappingAroundAPSrc.put(
                            molA.getAtom(mol.indexOf(pair.getValue())),
                            molB.getAtom(graphIAC.indexOf(atmGraphIAC)));
                    }
                }
            }
 
            // Make mapping in a format suitable for KabschAlignment
            IAtom[] lstA = new IAtom[molToGraphMappingAroundAPSrc.size()];
            IAtom[] lstB = new IAtom[molToGraphMappingAroundAPSrc.size()];
            int k = 0;
            for (Map.Entry<IAtom,IAtom> pair : molToGraphMappingAroundAPSrc.entrySet())
            {
                lstA[k] = pair.getKey();
                lstB[k] = pair.getValue();
                k++;
            }
 
            //Calculate rotational matrix and align
            KabschAlignment sa;
            try
            {
                sa = new KabschAlignment(lstA, lstB);
                sa.align();
            } catch (Throwable t) {
                throw new DENOPTIMException("KabschAlignment failed.", t);
            }
 
            //Rototranslation of molecule B (origin is in A's center of mass!!!)
            sa.rotateAtomContainer(molB);
 
            // Translate B to to actual coordinates of A (instead of c.o.m.)
            Point3d cm = sa.getCenterOfMass();
            for (int ib = 0; ib < molB.getAtomCount(); ib++)
            {
                IAtom a = molB.getAtom(ib);
                Point3d oldPlace = MoleculeUtils.getPoint3d(a);
                Point3d newPlace = new Point3d(oldPlace.x + cm.x,
                                                    oldPlace.y + cm.y,
                                                    oldPlace.z + cm.z);
                a.setPoint3d(newPlace);
            }
        }
        return MoleculeUtils.getPoint3d(dummyAtm);
    }
 
//------------------------------------------------------------------------------
    
    public static List<Vertex> fragmentation(IAtomContainer mol, 
            List<CuttingRule> rules, Logger logger) throws DENOPTIMException
    {   
        Fragment masterFrag = new Fragment(mol,BBType.UNDEFINED);
        IAtomContainer fragsMol = masterFrag.getIAtomContainer();
        
        // Identify bonds
        Map<String, List<MatchedBond>> matchingbonds = 
                FragmenterTools.getMatchingBondsAllInOne(fragsMol,rules,logger);
        
        // Select bonds to cut and what rule to use for cutting them
        int cutId = -1;
        for (CuttingRule rule : rules) // NB: iterator follows rule's priority
        {
            String ruleName = rule.getName();
 
            // Skip unmatched rules
            if (!matchingbonds.keySet().contains(ruleName))
                continue;
 
            for (MatchedBond tb: matchingbonds.get(ruleName)) 
            {
                IAtom atmA = tb.getAtmSubClass0();
                IAtom atmB = tb.getAtmSubClass1();
 
                //ignore if bond already broken
                if (!fragsMol.getConnectedAtomsList(atmA).contains(atmB))
                { 
                    continue;
                }
 
                // Initialize the source atoms for the attachment points
                // but may be changed to handle hapticity
                IAtom srcAPA = atmA;
                IAtom srcAPB = atmB;
 
                //treatment of n-hapto ligands
                if (rule.isHAPTO())
                {
                    // Get central atom (i.e., the "mono-hapto" side, 
                    // typically the metal)
                    // As a convention the central atom has subclass '0'
                    IAtom centralAtm = atmA;
 
                    // Get list of candidates for hapto-system: 
                    // they have same cutting Rule and central metal
                    ArrayList<IAtom> candidatesForHapto = new ArrayList<IAtom>();
                    for (MatchedBond tbForHapto : matchingbonds.get(ruleName))
                    {
                        //Consider only bond involving same central atom
                        if (tbForHapto.getAtmSubClass0() == centralAtm)
                            candidatesForHapto.add(tbForHapto.getAtmSubClass1());
                    }
 
                    // Select atoms in n-hapto system: contiguous neighbors with  
                    // same type of bond with the same central atom.
                    Set<IAtom> atmsInHapto = new HashSet<IAtom>();
                    atmsInHapto.add(tb.getAtmSubClass1());
                    atmsInHapto = exploreHapticity(tb.getAtmSubClass1(), 
                            centralAtm, candidatesForHapto, fragsMol);
                    if (atmsInHapto.size() == 1)
                    {
                        logger.log(Level.WARNING,"Unable to find more than one "
                                + "bond involved in high-hapticity ligand! "
                                + "Bond ignored.");
                        continue;
                    }
 
                    // Check existence of all bonds involved in multi-hapto system
                    boolean isSystemIntact = true;
                    for (IAtom ligAtm : atmsInHapto)
                    {
                        List<IAtom> nbrsOfLigAtm = 
                                fragsMol.getConnectedAtomsList(ligAtm);
                        if (!nbrsOfLigAtm.contains(centralAtm))
                        {
                            isSystemIntact = false;
                            break;
                        }
                    } 
 
                    // If not, it means that another rule already acted on the  
                    // system thus kill this attempt without generating du-atom
                    if (!isSystemIntact)
                        continue;
 
                    srcAPA = centralAtm;
                    srcAPB = getAPSourceAtom(masterFrag, new ArrayList<>(atmsInHapto), Arrays.asList(centralAtm));
                }
 
                // Storage of info for AP creation
                List<List<List<IAtom>>> listOfAtomPairs = new ArrayList<>();
                List<List<APClass>> listOfAPClasses = new ArrayList<>();
                listOfAtomPairs.add(Arrays.asList(Arrays.asList(srcAPA), Arrays.asList(srcAPB)));
                listOfAPClasses.add(Arrays.asList(rule.getAPClass0(), rule.getAPClass1()));
                cutId++;
                makeAPPairs(masterFrag, listOfAtomPairs, listOfAPClasses, cutId);
            } //end of loop over matching bonds
        } //end of loop over rules
        
        return isolateFragments(masterFrag);
    }
 
//------------------------------------------------------------------------------
 
    public static List<Vertex> fragmentation(IAtomContainer mol, 
        List<List<List<IAtom>>> atomPairs, List<List<APClass>> apClasses) 
        throws DENOPTIMException
    {   
        Fragment masterFrag = new Fragment(mol,BBType.UNDEFINED);
        IAtomContainer fragsMol = masterFrag.getIAtomContainer();
 
        // Project the list of atoms to the new atom container
        List<List<List<IAtom>>> projectedAtomPairs = new ArrayList<>();
        for (List<List<IAtom>> atomPair : atomPairs)
        {
            List<List<IAtom>> projectedAtomPair = new ArrayList<>();
            for (List<IAtom> leftOrRightMembers : atomPair)
            {
                List<IAtom> projectedLeftOrRightMembers = new ArrayList<>();
                for (IAtom atm : leftOrRightMembers)
                {
                    projectedLeftOrRightMembers.add(fragsMol.getAtom(mol.indexOf(atm)));
                }
                projectedAtomPair.add(projectedLeftOrRightMembers);
            }
            projectedAtomPairs.add(projectedAtomPair);
        }
        makeAPPairs(masterFrag, projectedAtomPairs, apClasses, 1);
        return isolateFragments(masterFrag);
    }
 
//------------------------------------------------------------------------------
 
    private static List<Vertex> isolateFragments(Fragment masterFrag) 
    throws DENOPTIMException
    {
        List<Vertex>  fragments = new ArrayList<Vertex>();
        Set<Integer> doneAlready = new HashSet<Integer>();
        for (int idx=0 ; idx<masterFrag.getAtomCount(); idx++)
        {
            if (doneAlready.contains(idx))
                continue;
            
            Fragment cloneOfMaster = masterFrag.clone();
            IAtomContainer iac = cloneOfMaster.getIAtomContainer();
            Set<IAtom> atmsToKeep = exploreConnectivity(iac.getAtom(idx), iac);
            atmsToKeep.stream().forEach(atm -> doneAlready.add(iac.indexOf(atm)));
            
            Set<IAtom> atmsToRemove = new HashSet<IAtom>();
            for (IAtom atm : cloneOfMaster.atoms())
            {
                if (!atmsToKeep.contains(atm))
                {
                    atmsToRemove.add(atm);
                }
            }
            cloneOfMaster.removeAtoms(atmsToRemove);
            if (cloneOfMaster.getAttachmentPoints().size()>0)
                fragments.add(cloneOfMaster);
        }
        return fragments;
    }
 
//------------------------------------------------------------------------------
 
    private static void  makeAPPairs(Fragment masterFrag, 
        List<List<List<IAtom>>> atomPairs, List<List<APClass>> apClasses, Integer cutId) 
        throws DENOPTIMException
    {
        for (int i = 0; i < atomPairs.size(); i++)
        {
            // "pair" in the sense that the list must have two items, but each 
            // item may actually be a list of atoms in case of hapticity>1.
            List<List<IAtom>> atomPair = atomPairs.get(i);
            APClass apClassA = apClasses.get(i).get(0);
            APClass apClassB = apClasses.get(i).get(1);
 
            IAtom srcAPA = getAPSourceAtom(masterFrag, atomPair.get(0), atomPair.get(1));
            IAtom srcAPB =  getAPSourceAtom(masterFrag, atomPair.get(1), atomPair.get(0));
            
            //treatment of mono-hapto ligands
            IBond bnd = masterFrag.getIAtomContainer().getBond(srcAPA,srcAPB);
            if (bnd != null)
            {
                masterFrag.removeBond(bnd);
            }
 
            AttachmentPoint apA = masterFrag.addAPOnAtom(srcAPA, apClassA, 
                MoleculeUtils.getPoint3d(srcAPB));
            AttachmentPoint apB = masterFrag.addAPOnAtom(srcAPB, apClassB, 
                MoleculeUtils.getPoint3d(srcAPA));
 
            cutId++;
            apA.setCutId(cutId);
            apB.setCutId(cutId);
        }
    }
 
//------------------------------------------------------------------------------
 
    private static IAtom getAPSourceAtom(Fragment masterFrag, List<IAtom> atmsInHapto, 
        List<IAtom> atmsOutsideHapto) throws DENOPTIMException
    {
        // Not really a multihapto system: return the single atoms
        if (atmsInHapto.size() == 1)
        {
            return atmsInHapto.get(0);
        }
 
        IAtomContainer fragsMol = masterFrag.getIAtomContainer();
        // A dummy atom will be used to define attachment point of
        // ligand with high hapticity
        Point3d dummyP3d = new Point3d(); //Used also for 2D
        for (IAtom ligAtm : atmsInHapto)
        {
            Point3d ligP3d = MoleculeUtils.getPoint3d(ligAtm);
            dummyP3d.x = dummyP3d.x + ligP3d.x;
            dummyP3d.y = dummyP3d.y + ligP3d.y;
            dummyP3d.z = dummyP3d.z + ligP3d.z;
        }
 
        dummyP3d.x = dummyP3d.x / (double) atmsInHapto.size();
        dummyP3d.y = dummyP3d.y / (double) atmsInHapto.size();
        dummyP3d.z = dummyP3d.z / (double) atmsInHapto.size();
 
        //Add Dummy atom to molecular object
        //if no other Du is already in the same position
        IAtom dummyAtm = null;
        for (IAtom oldDu : fragsMol.atoms())
        {
            if (MoleculeUtils.getSymbolOrLabel(oldDu) 
                    == DENOPTIMConstants.DUMMYATMSYMBOL)
            {
                Point3d oldDuP3d = oldDu.getPoint3d();
                if (oldDuP3d.distance(dummyP3d) < 0.002)
                {
                    dummyAtm = oldDu;
                    break;
                }
            } 
        }
    
        // Make the dummy, since it is not there already
        if (dummyAtm==null)
        {
            dummyAtm = new PseudoAtom(DENOPTIMConstants.DUMMYATMSYMBOL);
            dummyAtm.setPoint3d(dummyP3d);
            fragsMol.addAtom(dummyAtm);
        }
 
        // Modify connectivity of atoms involved in high-hapticity 
        // coordination creation of Du-to-ATM bonds 
        // By internal convention the bond order is "SINGLE".
        IBond.Order border = IBond.Order.valueOf("SINGLE");
        
        for (IAtom ligAtm : atmsInHapto)
        {
            List<IAtom> nbrsOfDu = fragsMol.getConnectedAtomsList(
                    dummyAtm);
            if (!nbrsOfDu.contains(ligAtm))
            {
                // Add bond with dummy
                Bond bnd = new Bond(dummyAtm,ligAtm,border);
                fragsMol.addBond(bnd);
            }
            // Remove bonds between central and coordinating atoms
            for (IAtom atmOutsideHapto : atmsOutsideHapto)
            {
                IBond oldBnd = fragsMol.getBond(atmOutsideHapto,ligAtm);
                fragsMol.removeBond(oldBnd);
            }
        }
        return dummyAtm;
    }
    
//------------------------------------------------------------------------------
    static Set<IAtom> exploreHapticity(IAtom seed, IAtom centralAtom, 
            ArrayList<IAtom> candidates, IAtomContainer mol)
    {
        Set<IAtom> atmsInHapto = new HashSet<IAtom>();
        atmsInHapto.add(seed);
        ArrayList<IAtom> toVisitAtoms = new ArrayList<IAtom>();
        toVisitAtoms.add(seed);
        ArrayList<IAtom> visitedAtoms = new ArrayList<IAtom>();
        while (toVisitAtoms.size()>0)
        {
            ArrayList<IAtom> toVisitLater = new ArrayList<IAtom>();
            for (IAtom atomInFocus : toVisitAtoms)
            {
                if (visitedAtoms.contains(atomInFocus) 
                        || atomInFocus==centralAtom)
                    continue;
                else 
                    visitedAtoms.add(atomInFocus);
                
                if (candidates.contains(atomInFocus))
                {
                    atmsInHapto.add(atomInFocus);
                    toVisitLater.addAll(mol.getConnectedAtomsList(atomInFocus));
                }
            }
            toVisitAtoms.clear();
            toVisitAtoms.addAll(toVisitLater);
        }
        return atmsInHapto;
    }
    
//------------------------------------------------------------------------------
    static Set<IAtom> exploreConnectivity(IAtom seed, IAtomContainer mol)
    {
        Set<IAtom> atmsReachableFromSeed = new HashSet<IAtom>();
        ArrayList<IAtom> toVisitAtoms = new ArrayList<IAtom>();
        toVisitAtoms.add(seed);
        ArrayList<IAtom> visitedAtoms = new ArrayList<IAtom>();
        while (toVisitAtoms.size()>0)
        {
            ArrayList<IAtom> toVisitLater = new ArrayList<IAtom>();
            for (IAtom atomInFocus : toVisitAtoms)
            {
                if (visitedAtoms.contains(atomInFocus))
                    continue;
                else 
                    visitedAtoms.add(atomInFocus);
                
                atmsReachableFromSeed.add(atomInFocus);
                toVisitLater.addAll(mol.getConnectedAtomsList(atomInFocus));
            }
            toVisitAtoms.clear();
            toVisitAtoms.addAll(toVisitLater);
        }
        return atmsReachableFromSeed;
    }
    
//-----------------------------------------------------------------------------
 
    static Map<String, List<MatchedBond>> getMatchingBondsAllInOne(
            IAtomContainer mol, List<CuttingRule> rules, Logger logger)
    {
        // Collect all SMARTS queries
        Map<String,String> smarts = new HashMap<String,String>();
        for (CuttingRule rule : rules)
        {
            smarts.put(rule.getName(),rule.getWholeSMARTSRule());
        }
 
        // Prepare a data structure for the return value
        Map<String, List<MatchedBond>> bondsMatchingRules = 
                new HashMap<String, List<MatchedBond>>();
 
        // Get all the matches to the SMARTS queries
        ManySMARTSQuery msq = new ManySMARTSQuery(mol, smarts);
        if (msq.hasProblems())
        {
            if (logger!=null)
            {
                logger.log(Level.WARNING, "Problem matching SMARTS: " 
                        + msq.getMessage());
            }
            return bondsMatchingRules;
        }
 
        for (CuttingRule rule : rules)
        {
            String ruleName = rule.getName();
 
            if (msq.getNumMatchesOfQuery(ruleName) == 0)
            {
                continue;
            }
           
            // Get atoms matching cutting rule queries
            Mappings purgedPairs = msq.getMatchesOfSMARTS(ruleName);
            
            // Evaluate subclass membership and eventually store target bonds
            ArrayList<MatchedBond> bondsMatched = new ArrayList<MatchedBond>();
            for (int[] pair : purgedPairs) 
            {
                if (pair.length!=2)
                {
                    throw new Error("Cutting rule: " + ruleName 
                            + " has identified " + pair.length + " atoms "
                            + "instead of 2. Modify rule to make it find a "
                            + "pair of atoms.");
                }
                MatchedBond tb = new MatchedBond(mol.getAtom(pair[0]),
                        mol.getAtom(pair[1]), rule);
                
                // Apply any further option of the cutting rule
                if (tb.satisfiesRuleOptions(logger))
                    bondsMatched.add(tb);
            }
 
            if (!bondsMatched.isEmpty())
                bondsMatchingRules.put(ruleName, bondsMatched);
        }
 
        return bondsMatchingRules;
    }
    
//------------------------------------------------------------------------------
    
    public static void manageFragmentCollection(File input, 
            FragmenterParameters settings,
            File output, Logger logger) throws DENOPTIMException, IOException, 
                IllegalArgumentException, UndetectedFileFormatException
    {
        FileInputStream fis = new FileInputStream(input);
        IteratingSDFReader reader = new IteratingSDFReader(fis, 
                DefaultChemObjectBuilder.getInstance());
 
        int index = -1;
        int maxBufferSize = 2000;
        ArrayList<Vertex> buffer = new ArrayList<Vertex>(500);
        try {
            while (reader.hasNext())
            {
                index++;
                if (logger!=null)
                {
                    logger.log(Level.FINE,"Processing fragment " + index);
                }
                Vertex frag = new Fragment(reader.next(), BBType.UNDEFINED);
                manageFragmentCollection(frag, index, settings,
                        buffer, logger);
                
                // If max buffer size is reached, then bump to file
                if (buffer.size() >= maxBufferSize)
                {
                    DenoptimIO.writeVertexesToFile(output, FileFormat.VRTXSDF, 
                            buffer, true);
                    buffer.clear();
                }
            }
        } finally {
            reader.close();
        }
        if (buffer.size() < maxBufferSize)
        {
            DenoptimIO.writeVertexesToFile(output, FileFormat.VRTXSDF, 
                    buffer, true);
            buffer.clear();
        }
    }
    
//------------------------------------------------------------------------------
    
    public static void manageFragmentCollection(Vertex frag, int fragCounter, 
            FragmenterParameters settings, 
            List<Vertex> collector, Logger logger) 
                    throws DENOPTIMException, IllegalArgumentException, 
                    UndetectedFileFormatException, IOException
    {
 
        if (!filterFragment((Fragment) frag, settings, logger))
        {
            return;
        }
        
        //Compare with list of fragments to ignore
        if (settings.getIgnorableFragments().size() > 0)
        {
            if (settings.getIgnorableFragments().stream()
                    .anyMatch(ignorable -> ((Fragment)frag)
                            .isIsomorphicTo(ignorable)))
            {
                if (logger!=null)
                {
                    logger.log(Level.FINE,"Fragment " + fragCounter 
                            + " is ignorable.");
                }
                return;
            }
        }
        
        //Compare with list of fragments to retain
        if (settings.getTargetFragments().size() > 0)
        {
            if (!settings.getTargetFragments().stream()
                    .anyMatch(ignorable -> ((Fragment)frag)
                            .isIsomorphicTo(ignorable)))
            {
                if (logger!=null)
                {
                    logger.log(Level.FINE,"Fragment " + fragCounter 
                          + " doesn't match any target: rejected.");
                }
                return;
            }
        }
        
        // Add dummy atoms on linearities
        if (MoleculeUtils.getDimensions(frag.getIAtomContainer())==3 
                && settings.doAddDuOnLinearity())
        {
            DummyAtomHandler.addDummiesOnLinearities((Fragment) frag,
                    settings.getLinearAngleLimit());
        }
        
         // Management of duplicate fragments:
         // -> identify duplicates (isomorphic fragments), 
         // -> keep one (or more, if we want to sample the isomorphs),
         // -> reject the rest.
        if (settings.doManageIsomorphicFamilies())
        {
            synchronized (settings.MANAGEMWSLOTSSLOCK)
            {
                String mwSlotID = getMWSlotIdentifier(frag, 
                        settings.getMWSlotSize());
 
                File mwFileUnq = settings.getMWSlotFileNameUnqFrags(
                        mwSlotID);
                File mwFileAll = settings.getMWSlotFileNameAllFrags(
                        mwSlotID);
                
                // Compare this fragment with previously seen ones
                Vertex unqVersion = null;
                if (mwFileUnq.exists())
                {
                    ArrayList<Vertex> knownFrags = 
                            DenoptimIO.readVertexes(mwFileUnq,BBType.UNDEFINED);
                    unqVersion = knownFrags.stream()
                        .filter(knownFrag -> 
                            ((Fragment)frag).isIsomorphicTo(knownFrag))
                        .findAny()
                        .orElse(null);
                }
                if (unqVersion!=null)
                {
                    // Identify this unique fragment
                    String isoFamID = unqVersion.getProperty(
                            DENOPTIMConstants.ISOMORPHICFAMILYID)
                            .toString();
                    
                    // Do we already have enough isomorphic family members 
                    // for this fragment?
                    int sampleSize = settings.getIsomorphsCount()
                            .get(isoFamID);
                    if (sampleSize < settings.getIsomorphicSampleSize())
                    {
                        // Add this isomorphic version to the sample
                        frag.setProperty(
                                DENOPTIMConstants.ISOMORPHICFAMILYID,
                                isoFamID);
                        settings.getIsomorphsCount().put(isoFamID,
                                sampleSize+1);
                        DenoptimIO.writeVertexToFile(mwFileAll, 
                                FileFormat.VRTXSDF, frag, true);
                        collector.add(frag);
                    } else {
                        // This would be inefficient in the long run
                        // because it by-passes the splitting by MW. 
                        // Do not do it!
                        /*
                        if (logger!=null)
                        {
                            logger.log(Level.FINE,"Fragment " 
                                  + fragCounter 
                                  + " is isomorphic to unique fragment " 
                                  + unqVersionID + ", but we already "
                                  + "have a sample of " + sampleSize
                                  + ": ignoring this fragment from now "
                                  + "on.");
                        }
                        settings.getIgnorableFragments().add(frag);
                        */
                    }
                } else {
                    // This is a never-seen fragment
                    String isoFamID = settings.newIsomorphicFamilyID();
                    frag.setProperty(
                            DENOPTIMConstants.ISOMORPHICFAMILYID,
                            isoFamID);
                    settings.getIsomorphsCount().put(isoFamID, 1);
                    DenoptimIO.writeVertexToFile(mwFileUnq, 
                            FileFormat.VRTXSDF, frag, true);
                    DenoptimIO.writeVertexToFile(mwFileAll, 
                            FileFormat.VRTXSDF, frag, true);
                    collector.add(frag);
                }
            } // end synchronized block
        } else {
            //If we are here, we did not ask to remove duplicates
            collector.add(frag);
        }
    }
 
//------------------------------------------------------------------------------
    
    public static boolean filterFragment(Fragment frag, 
            FragmenterParameters settings)
    {
       return filterFragment(frag, settings, settings.getLogger());
    }
 
//------------------------------------------------------------------------------
    
    public static boolean filterFragment(Fragment frag, 
            FragmenterParameters settings, Logger logger)
    {
        // Default filtering criteria: get ring of R/*/X/Xx
        for (IAtom atm : frag.atoms())
        {
            if (MoleculeUtils.isElement(atm))
            {
                continue;
            }
            String smb = MoleculeUtils.getSymbolOrLabel(atm);
            if (DENOPTIMConstants.DUMMYATMSYMBOL.equals(smb))
            {
                continue;
            }
            logger.log(Level.FINE,"Removing fragment contains non-element '"
                    + smb + "'");
            return false;
        }
 
        if (settings.isWorkingIn3D()) 
        {   
            // Incomplete 3D fragmentation: an atom has the same coords of an AP.
            for (AttachmentPoint ap : frag.getAttachmentPoints())
            {
                Point3d ap3d = ap.getDirectionVector();
                if (ap3d!=null)
                {
                    for (IAtom atm : frag.atoms())
                    {
                        Point3d atm3d = MoleculeUtils.getPoint3d(atm);
                        double dist = ap3d.distance(atm3d);
                        if (dist < 0.0002)
                        {
                            logger.log(Level.FINE,"Removing fragment with AP"
                            + frag.getIAtomContainer().indexOf(atm)
                            + " and atom " + MoleculeUtils.getSymbolOrLabel(atm)
                            + " coincide.");
                            return false;
                        }   
                    }
                }
            }
        }
        if (settings.doRejectWeirdIsotopes())
        {
            for (IAtom atm : frag.atoms())
            {
                if (MoleculeUtils.isElement(atm))
                {
                    // Unconfigured isotope has null mass number
                    if (atm.getMassNumber() == null)
                        continue;
                    
                    String symb =  MoleculeUtils.getSymbolOrLabel(atm);
                    int a = atm.getMassNumber();
                    try {
                        IIsotope major = Isotopes.getInstance().getMajorIsotope(symb);
                        if (a != major.getMassNumber())
                        {
                            logger.log(Level.FINE,"Removing fragment containing "
                                    + "isotope "+symb+a+".");
                            return false;
                        }
                    } catch (Throwable t) {
                        logger.log(Level.WARNING,"Not able to perform Isotope"
                                + "detection.");
                    }
                }
                
            }
        }
        
        // User-controlled filtering criteria
        
        if (settings.getRejectedElements().size() > 0)
        {
            for (IAtom atm : frag.atoms())
            {
                String symb = MoleculeUtils.getSymbolOrLabel(atm);
                if (settings.getRejectedElements().contains(symb))
                {
                    logger.log(Level.FINE,"Removing fragment containing '" 
                            + symb + "'.");
                    return false;
                }
            }
        }
        
        if (settings.getRejectedFormulaLessThan().size() > 0
                || settings.getRejectedFormulaMoreThan().size() > 0)
        {
            Map<String,Double> eaMol = FormulaUtils.getElementalanalysis(
                    frag.getIAtomContainer());
            
            for (Map<String,Double> criterion : 
                settings.getRejectedFormulaMoreThan())
            {
                for (String el : criterion.keySet())
                {
                    if (eaMol.containsKey(el))
                    {
                        // -0.5 to make it strictly less-than
                        if (eaMol.get(el) - criterion.get(el) > 0.5)
                        {
                            logger.log(Level.FINE,"Removing fragment that "
                                    + "contains too much '" + el + "' "
                                    + "as requested by formula"
                                    + "-based (more-than) settings (" + el
                                    + eaMol.get(el)  + " > " + criterion + ").");
                            return false;
                        }
                    }
                }
            }
            
            Map<String,Double> criterion = settings.getRejectedFormulaLessThan();
            for (String el : criterion.keySet())
            {
                if (!eaMol.containsKey(el))
                {
                    logger.log(Level.FINE,"Removing fragment that does not "
                            + "contain '" + el + "' as requested by formula"
                            + "-based (less-than) settings.");
                    return false;
                } else {
                    // 0.5 to make it strictly more-than
                    if (eaMol.get(el) - criterion.get(el) < -0.5)
                    {
                        logger.log(Level.FINE,"Removing fragment that "
                                + "contains too little '" + el + "' "
                                + "as requested by formula"
                                + "-based settings (" + el
                                + eaMol.get(el)  + " < " + criterion + ").");
                        return false;
                    }
                }
            }
            
        }
        
        if (settings.getRejectedAPClasses().size() > 0)
        {
            for (APClass apc : frag.getAllAPClasses())
            {
                for (String s : settings.getRejectedAPClasses())
                {
                    if (apc.toString().startsWith(s))
                    {
                        logger.log(Level.FINE,"Removing fragment with APClass "
                                + apc);
                        return false;
                    }    
                }
            }
        }
        
        if (settings.getRejectedAPClassCombinations().size() > 0)
        {
            loopOverCombinations:
            for (String[] conditions : settings.getRejectedAPClassCombinations())
            {
                for (int ip=0; ip<conditions.length; ip++)
                {   
                    String condition = conditions[ip];
                    boolean found = false;
                    for (APClass apc : frag.getAllAPClasses())
                    {
                        if (apc.toString().startsWith(condition))
                        {
                            found = true;
                            continue;
                        }
                    }
                    if (!found)
                        continue loopOverCombinations;
                    // Here we do have at least one AP satisfying the condition.
                }
                // Here we manage or satisfy all conditions. Therefore, we can
                // reject this fragment
                
                String allCondsAsString = "";
                for (int i=0; i<conditions.length; i++)
                    allCondsAsString = allCondsAsString + " " + conditions[i];
                
                logger.log(Level.FINE,"Removing fragment with combination of "
                        + "APClasses matching '" + allCondsAsString + "'.");
                return false;
            }
        }
        
        if (settings.getMaxFragHeavyAtomCount()>0 
                || settings.getMinFragHeavyAtomCount()>0)
        {
            int totHeavyAtm = 0;
            for (IAtom atm : frag.atoms())
            {
                if (MoleculeUtils.isElement(atm))
                {
                    String symb = MoleculeUtils.getSymbolOrLabel(atm);
                    if ((!symb.equals("H")) && (!symb.equals(
                            DENOPTIMConstants.DUMMYATMSYMBOL)))
                        totHeavyAtm++;
                }
            }
            if (settings.getMaxFragHeavyAtomCount() > 0 
                    && totHeavyAtm > settings.getMaxFragHeavyAtomCount())
            {
                logger.log(Level.FINE,"Removing fragment with too many atoms (" 
                        + totHeavyAtm + " < " 
                        + settings.getMaxFragHeavyAtomCount() 
                        + ")");
                return false;
            }
            if (settings.getMinFragHeavyAtomCount() > 0 
                    && totHeavyAtm < settings.getMinFragHeavyAtomCount())
            {
                logger.log(Level.FINE,"Removing fragment with too few atoms (" 
                        + totHeavyAtm + " < " 
                        + settings.getMinFragHeavyAtomCount() 
                        + ")");
                return false;
            }
        }
        
        if (settings.getFragRejectionSMARTS().size() > 0)
        {
            ManySMARTSQuery msq = new ManySMARTSQuery(frag.getIAtomContainer(),
                    settings.getFragRejectionSMARTS());
            if (msq.hasProblems())
            {
                logger.log(Level.WARNING,"Problems evaluating SMARTS-based "
                        + "rejection criteria. " + msq.getMessage());
            }
        
            for (String criterion : settings.getFragRejectionSMARTS().keySet())
            {
                if (msq.getNumMatchesOfQuery(criterion)>0)
                {
                    logger.log(Level.FINE,"Removing fragment that matches "
                            + "SMARTS-based rejection criteria '" + criterion 
                            + "'.");
                    return false;
                }
            }
        }
        
        if (settings.getFragRetentionSMARTS().size() > 0)
        {
            ManySMARTSQuery msq = new ManySMARTSQuery(frag.getIAtomContainer(),
                    settings.getFragRetentionSMARTS());
            if (msq.hasProblems())
            {
                logger.log(Level.WARNING,"Problems evaluating SMARTS-based "
                        + "rejection criteria. " + msq.getMessage());
            }
        
            boolean matchesAny = false;
            for (String criterion : settings.getFragRetentionSMARTS().keySet())
            {
                if (msq.getNumMatchesOfQuery(criterion) > 0)
                {
                    matchesAny = true;
                    break;
                }
            }
            if (!matchesAny)
            {
                logger.log(Level.FINE,"Removing fragment that does not "
                        + "match any SMARTS-based retention criteria.");
                return false;
            }
        }
        return true;
    }
    
//------------------------------------------------------------------------------
  
    public static String getMWSlotIdentifier(Vertex frag, int slotSize)
    {
        for (IAtom a : frag.getIAtomContainer().atoms())
        {
            if (a.getImplicitHydrogenCount()==null)
                a.setImplicitHydrogenCount(0);
        }
        double mw = AtomContainerManipulator.getMass(frag.getIAtomContainer());
        int slotNum = (int) (mw / (Double.valueOf(slotSize)));
        return slotNum*slotSize + "-" + (slotNum+1)*slotSize;
    }
    
//------------------------------------------------------------------------------
    
    public static Vertex getRCVForAP(AttachmentPoint ap, APClass rcvApClass) 
            throws DENOPTIMException
    {
        IAtomContainer mol = SilentChemObjectBuilder.getInstance()
                .newAtomContainer();
        Point3d apv = ap.getDirectionVector();
        mol.addAtom(new PseudoAtom(RingClosingAttractor.RCALABELPERAPCLASS.get(
                    rcvApClass), 
                new Point3d(
                    Double.valueOf(apv.x),
                    Double.valueOf(apv.y),
                    Double.valueOf(apv.z))));
        
        Fragment rcv = new Fragment(mol, BBType.FRAGMENT);
        rcv.setAsRCV(true);
 
        Point3d aps = MoleculeUtils.getPoint3d(
                ap.getOwner().getIAtomContainer().getAtom(
                        ap.getAtomPositionNumber()));
        rcv.addAP(0, rcvApClass, new Point3d(
                Double.valueOf(aps.x),
                Double.valueOf(aps.y),
                Double.valueOf(aps.z)));
        return rcv;
    }
    
//------------------------------------------------------------------------------
 
}