MoleculeUtils.java

Go to the documentation of this file.

/*
 *   DENOPTIM
 *   Copyright (C) 2019 Vishwesh Venkatraman <vishwesh.venkatraman@ntnu.no> and
 *   Marco Foscato <marco.foscato@uib.no>
 * 
 *   This program is free software: you can redistribute it and/or modify
 *   it under the terms of the GNU Affero General Public License as published
 *   by the Free Software Foundation, either version 3 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU Affero General Public License for more details.
 *
 *   You should have received a copy of the GNU Affero General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
package denoptim.utils;
 
import java.io.File;
import java.util.ArrayList;
import java.util.Collections;
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.Point2d;
import javax.vecmath.Point3d;
 
import org.openscience.cdk.Atom;
import org.openscience.cdk.AtomContainer;
import org.openscience.cdk.AtomRef;
import org.openscience.cdk.CDKConstants;
import org.openscience.cdk.PseudoAtom;
import org.openscience.cdk.aromaticity.Kekulization;
import org.openscience.cdk.config.IsotopeFactory;
import org.openscience.cdk.config.Isotopes;
import org.openscience.cdk.depict.Depiction;
import org.openscience.cdk.depict.DepictionGenerator;
import org.openscience.cdk.exception.CDKException;
import org.openscience.cdk.geometry.GeometryUtil;
import org.openscience.cdk.graph.ConnectivityChecker;
import org.openscience.cdk.inchi.InChIGenerator;
import org.openscience.cdk.inchi.InChIGeneratorFactory;
import org.openscience.cdk.interfaces.IAtom;
import org.openscience.cdk.interfaces.IAtomContainer;
import org.openscience.cdk.interfaces.IAtomContainerSet;
import org.openscience.cdk.interfaces.IAtomType.Hybridization;
import org.openscience.cdk.interfaces.IBond;
import org.openscience.cdk.interfaces.IChemObjectBuilder;
import org.openscience.cdk.interfaces.IElement;
import org.openscience.cdk.layout.StructureDiagramGenerator;
import org.openscience.cdk.qsar.DescriptorValue;
import org.openscience.cdk.qsar.IMolecularDescriptor;
import org.openscience.cdk.qsar.descriptors.molecular.RotatableBondsCountDescriptor;
import org.openscience.cdk.qsar.descriptors.molecular.WeightDescriptor;
import org.openscience.cdk.qsar.result.DoubleResult;
import org.openscience.cdk.qsar.result.IntegerResult;
import org.openscience.cdk.silent.SilentChemObjectBuilder;
import org.openscience.cdk.smiles.SmiFlavor;
import org.openscience.cdk.smiles.SmilesGenerator;
import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
 
import denoptim.constants.DENOPTIMConstants;
import denoptim.exception.DENOPTIMException;
import denoptim.graph.APClass;
import denoptim.graph.AttachmentPoint;
import denoptim.graph.DGraph;
import denoptim.graph.Edge.BondType;
import denoptim.graph.Fragment;
import denoptim.graph.Ring;
import denoptim.graph.Vertex;
import denoptim.graph.Vertex.BBType;
import denoptim.graph.rings.RingClosingAttractor;
import denoptim.io.DenoptimIO;
import denoptim.logging.StaticLogger;
import io.github.dan2097.jnainchi.InchiFlag;
import io.github.dan2097.jnainchi.InchiOptions;
import io.github.dan2097.jnainchi.InchiStatus;
 
 
public class MoleculeUtils
{
    private static final StructureDiagramGenerator SDG = 
            new StructureDiagramGenerator();
    private static final SmilesGenerator SMGEN = new SmilesGenerator(
            SmiFlavor.Generic);
    private static final  IChemObjectBuilder builder = 
            SilentChemObjectBuilder.getInstance();
 
//------------------------------------------------------------------------------
 
    public static double getMolecularWeight(IAtomContainer mol) 
            throws DENOPTIMException
    {
        double ret_wd = 0;
        try
        {
            WeightDescriptor wd = new WeightDescriptor();
            Object[] pars = {"*"};
            wd.setParameters(pars);
            ret_wd = ((DoubleResult) wd.calculate(mol).getValue())
                    .doubleValue();
        }
        catch (Exception ex)
        {
            throw new DENOPTIMException(ex);
        }
        return ret_wd;
    }
    
//------------------------------------------------------------------------------
    
    public static boolean isDummy(IAtom atm)
    {
        String el = getSymbolOrLabel(atm);
        return DENOPTIMConstants.DUMMYATMSYMBOL.equals(el);
    }
    
//------------------------------------------------------------------------------
 
    public static boolean isElement(IAtom atom)
    {
        String symbol = atom.getSymbol();
        return isElement(symbol);
    }
    
//------------------------------------------------------------------------------
 
    public static boolean isElement(String symbol)
    {
        boolean res = false;
        IsotopeFactory ifact = null;
        try {
            ifact = Isotopes.getInstance();
            if (ifact.isElement(symbol))
            {
                @SuppressWarnings("unused")
                IElement el = ifact.getElement(symbol);
                res = true;
            }
        } catch (Throwable t) {
            throw new Error("ERROR! Unable to create Isotope.");
        }
        return res;
    }
 
//------------------------------------------------------------------------------
 
    public static void removeRCA(IAtomContainer mol) {
 
        // convert PseudoAtoms to H
        for (IAtom a : mol.atoms())
        {
            boolean isRca = false;
            Set<String> rcaElSymbols = RingClosingAttractor.RCATYPEMAP.keySet();
            for (String rcaEl : rcaElSymbols)
            {
                if (MoleculeUtils.getSymbolOrLabel(a).equals(rcaEl))
                {
                    isRca = true;
                    break;
                }
            }
            if (isRca)
            {
                IAtom newAtm = new Atom("H",new Point3d(a.getPoint3d()));
                newAtm.setProperties(a.getProperties());
                AtomContainerManipulator.replaceAtomByAtom(mol,a,newAtm);
            }
        }
    }
 
//------------------------------------------------------------------------------
 
    public static void removeUsedRCA(IAtomContainer mol, DGraph graph, 
            Logger logger) throws DENOPTIMException {
 
        // add ring-closing bonds
        ArrayList<Vertex> usedRcvs = graph.getUsedRCVertices();
        Map<Vertex,ArrayList<Integer>> vIdToAtmId =
                MoleculeUtils.getVertexToAtomIdMap(usedRcvs,mol);
        if (vIdToAtmId.size() == 0)
        {
            // No used RCV to remove.
            return;
        }
        ArrayList<IAtom> atmsToRemove = new ArrayList<>();
        ArrayList<Boolean> doneVertices =
                new ArrayList<>(Collections.nCopies(usedRcvs.size(),false));
 
        for (Vertex v : usedRcvs)
        {
            if (doneVertices.get(usedRcvs.indexOf(v)))
            {
                continue;
            }
            ArrayList<Ring> rings = graph.getRingsInvolvingVertex(v);
            if (rings.size() != 1)
            {
                String s = "Unexpected inconsistency between used RCV list "
                       + v + " in {" + usedRcvs + "}"
                       + "and list of DENOPTIMRings "
                       + "{" + rings + "}. Check Code!";
                throw new DENOPTIMException(s);
            }
            Vertex vH = rings.get(0).getHeadVertex();
            Vertex vT = rings.get(0).getTailVertex();
            IAtom aH = mol.getAtom(vIdToAtmId.get(vH).get(0));
            IAtom aT = mol.getAtom(vIdToAtmId.get(vT).get(0));
            if (mol.getConnectedAtomsList(aH).size() == 0
                    || mol.getConnectedAtomsList(aT).size() == 0)
            {
                // This can happen when building a graph with empty vertexes
                continue;
            }
            int iSrcH = mol.indexOf(mol.getConnectedAtomsList(aH).get(0));
            int iSrcT = mol.indexOf(mol.getConnectedAtomsList(aT).get(0));
            atmsToRemove.add(aH);
            atmsToRemove.add(aT);
 
            BondType bndTyp = rings.get(0).getBondType();
            if (bndTyp.hasCDKAnalogue())
            {
                mol.addBond(iSrcH, iSrcT, bndTyp.getCDKOrder());
            } else {
                logger.log(Level.WARNING, "WARNING! "
                        + "Attempt to add ring closing bond "
                        + "did not add any actual chemical bond because the "
                        + "bond type of the chord is '" + bndTyp +"'.");
            }
 
            doneVertices.set(usedRcvs.indexOf(vH),true);
            doneVertices.set(usedRcvs.indexOf(vT),true);
        }
        
        // Adapt atom indexes in APs to the upcoming change of atom list
        ArrayList<Integer> removedIds = new ArrayList<Integer>();
        for (IAtom a : atmsToRemove)
        {
            removedIds.add(mol.indexOf(a));
        }
        Collections.sort(removedIds);
        for (Vertex v : graph.getVertexList())
        {
            for (AttachmentPoint ap : v.getAttachmentPoints())
            {
                int apSrcId = ap.getAtomPositionNumberInMol();
                int countOfAtmsBEforeSrc = 0;
                for (Integer removingId : removedIds)
                {
                    if (removingId < apSrcId)
                    {
                        countOfAtmsBEforeSrc++; 
                    } else if (removingId > apSrcId)
                    {
                        break;
                    }
                }
                ap.setAtomPositionNumberInMol(apSrcId-countOfAtmsBEforeSrc);
            }
        }
 
        // remove used RCAs
        for (IAtom a : atmsToRemove)
        {
            mol.removeAtom(a);
        }
    }
 
//------------------------------------------------------------------------------
 
    public static String getSMILESForMolecule(IAtomContainer mol, Logger logger)
            throws DENOPTIMException {
        IAtomContainer fmol = builder.newAtomContainer();
        try 
        {
            fmol = mol.clone();
        }  
        catch (CloneNotSupportedException e)
        {
            throw new DENOPTIMException(e);
        }
 
        // remove Dummy atoms
        DummyAtomHandler dan = new DummyAtomHandler(
                DENOPTIMConstants.DUMMYATMSYMBOL, logger);
        fmol = dan.removeDummyInHapto(fmol);
        fmol = dan.removeDummy(fmol);
 
        // convert PseudoAtoms to H
        removeRCA(fmol);
        
        // WARNING: assumptions on implicit H count and bond orders!
        MoleculeUtils.setZeroImplicitHydrogensToAllAtoms(fmol);
        MoleculeUtils.ensureNoUnsetBondOrdersSilent(fmol);
 
        String smiles = "";
        try
        {
            smiles = SMGEN.create(fmol);
        }
        catch (Throwable t)
        {
            t.printStackTrace();
            String fileName = "failed_generation_of_SMILES.sdf";
            logger.log(Level.WARNING, "WARNING: Skipping calculation of SMILES. "
                    + "See file '" + fileName + "'");
            DenoptimIO.writeSDFFile(fileName,fmol,false);
            smiles = "calculation_or_SMILES_crashed";
        }
        return smiles;
    }
 
//------------------------------------------------------------------------------
 
    public static IAtomContainer generate2DCoordinates(IAtomContainer ac, 
            Logger logger) throws DENOPTIMException
    {
        IAtomContainer fmol = builder.newAtomContainer();
        try 
        { 
            fmol = ac.clone();
        }  
        catch (CloneNotSupportedException e)
        {
            throw new DENOPTIMException(e);
        }
 
        // remove Dummy atoms
        DummyAtomHandler dan = new DummyAtomHandler(
                DENOPTIMConstants.DUMMYATMSYMBOL, logger);
        fmol = dan.removeDummyInHapto(fmol);
 
        // remove ring-closing attractors
        removeRCA(fmol);
 
        // Generate 2D structure diagram (for each connected component).
        IAtomContainer ac2d = builder.newAtomContainer();
        IAtomContainerSet som = ConnectivityChecker.partitionIntoMolecules(
                fmol);
 
        for (int n = 0; n < som.getAtomContainerCount(); n++)
        {
            synchronized (SDG)
            {
                IAtomContainer mol = som.getAtomContainer(n);
                SDG.setMolecule(mol, true);
                try
                {
                    // Generate 2D coordinates for this molecule.
                    SDG.generateCoordinates();
                    mol = SDG.getMolecule();
                }
                catch (Exception e)
                {
                    throw new DENOPTIMException(e);
                }
 
                ac2d.add(mol);  // add 2D molecule.
            }
        }
        
        return GeometryUtil.has2DCoordinates(ac2d) ? ac2d : null;
    }
 
//------------------------------------------------------------------------------
 
    public static String getInChIKeyForMolecule(IAtomContainer mol, 
            Logger logger) throws DENOPTIMException {
        InchiOptions options = new InchiOptions.InchiOptionsBuilder()
                .withFlag(InchiFlag.AuxNone)
                .withFlag(InchiFlag.RecMet)
                .withFlag(InchiFlag.SUU)
                .build();
        return getInChIKeyForMolecule(mol, options, logger);
    }
    
//------------------------------------------------------------------------------
 
    public static String getInChIKeyForMolecule(IAtomContainer mol, 
            InchiOptions options, Logger logger) throws DENOPTIMException {
        IAtomContainer fmol = builder.newAtomContainer();
        try 
        { 
            fmol = mol.clone();
        }  
        catch (Throwable t) 
        {
            throw new DENOPTIMException(t);
        }
 
        // remove Dummy atoms before generating the inchi
        DummyAtomHandler dan = new DummyAtomHandler(
                DENOPTIMConstants.DUMMYATMSYMBOL, logger);
        fmol = dan.removeDummyInHapto(fmol);
 
        // remove PseudoAtoms
        removeRCA(fmol);
 
        String inchikey;
        try
        {
            InChIGeneratorFactory factory = InChIGeneratorFactory.getInstance();
            InChIGenerator gen = factory.getInChIGenerator(fmol, options);
            InchiStatus ret = gen.getStatus();
            if (ret == InchiStatus.WARNING)
            {
                //String error = gen.getMessage();
                //InChI generated, but with warning message
                //return new ObjectPair(null, error);
            }
            else if (ret != InchiStatus.SUCCESS)
            {
                return null;
            }
            inchikey = gen.getInchiKey();
        }
        catch (CDKException cdke)
        {
            throw new DENOPTIMException(cdke);
        }
        if (inchikey.length() > 0)
            return inchikey;
        else
            return null;
    }
 
//------------------------------------------------------------------------------
 
    public static int getNumberOfRotatableBonds(IAtomContainer mol)
            throws DENOPTIMException {
        int value;
        try
        {
            IMolecularDescriptor descriptor = 
                    new RotatableBondsCountDescriptor();
            descriptor.setParameters(new Object[]{Boolean.FALSE,Boolean.FALSE});
            DescriptorValue result = descriptor.calculate(mol);
            value = ((IntegerResult)result.getValue()).intValue();
        }
        catch (CDKException cdke)
        {
            throw new DENOPTIMException(cdke);
        }
        return value;
    }
 
//------------------------------------------------------------------------------
 
    public static int getHeavyAtomCount(IAtomContainer mol)
    {
        int n = 0;
        for (IAtom atm : mol.atoms())
        {
            if (MoleculeUtils.isElement(atm)
                    && !MoleculeUtils.getSymbolOrLabel(atm).equals("H"))
                n++;
        }
        return n;
    }
    
//------------------------------------------------------------------------------
 
    public static int countAtomsOfElement(IAtomContainer mol, String symbol)
    {
        int n = 0;
        for (IAtom atm : mol.atoms())
        {
            if (atm.getSymbol().equals(symbol))
                n++;
        }
        return n;
    }
 
//------------------------------------------------------------------------------
 
    public static Map<Vertex,ArrayList<Integer>> getVertexToAtomIdMap(
            ArrayList<Vertex> vertLst,
            IAtomContainer mol
    ) {
 
        ArrayList<Long> vertIDs = new ArrayList<>();
        for (Vertex v : vertLst) {
            vertIDs.add(v.getVertexId());
        }
 
        Map<Vertex,ArrayList<Integer>> map = new HashMap<>();
        for (IAtom atm : mol.atoms())
        {
            long vID = Long.parseLong(atm.getProperty(
                                 DENOPTIMConstants.ATMPROPVERTEXID).toString());
            if (vertIDs.contains(vID))
            {
                Vertex v = vertLst.get(vertIDs.indexOf(vID));
                int atmID = mol.indexOf(atm);
                if (map.containsKey(v))
                {
                    map.get(v).add(atmID);
                }
                else
                {
                    ArrayList<Integer> atmLst = new ArrayList<>();
                    atmLst.add(atmID);
                    map.put(v, atmLst);
                }
            }
        }
        return map;
    }
 
//------------------------------------------------------------------------------
    
    public static void moleculeToPNG(IAtomContainer mol, String filename, 
            Logger logger) throws DENOPTIMException
    {
        IAtomContainer iac = mol;
 
        if (!GeometryUtil.has2DCoordinates(mol))
        {
            iac = generate2DCoordinates(mol, logger);
        }
        
        if (iac == null)
        {
            throw new DENOPTIMException("Failed to generate 2D coordinates.");
        }
        
        try
        {
            Depiction depiction = new DepictionGenerator().depict(iac);
            depiction.writeTo(filename);
        } catch (Exception e)
        {
            throw new DENOPTIMException("Failed to write image to "+filename);
        }
    }
    
//------------------------------------------------------------------------------
 
    public static Point3d getPoint3d(IAtom atm)
    {
        Point3d p = atm.getPoint3d();
 
        if (p == null)
        {
            Point2d p2d = atm.getPoint2d();
            if (p2d == null)
            {
                p = new Point3d(0.0, 0.0, 0.0);
            } else {
                p = new Point3d(p2d.x, p2d.y, 0.0);
            }
        }
        return p;
    }
    
//------------------------------------------------------------------------------
 
    //TODO: should we set only values that would otherwise be null?
    
    public static void setZeroImplicitHydrogensToAllAtoms(IAtomContainer iac)
    {
        for (IAtom atm : iac.atoms()) {
            atm.setImplicitHydrogenCount(0);
        }
    }
 
//------------------------------------------------------------------------------
    
    public static void explicitHydrogens(IAtomContainer mol)
    {
        AtomContainerManipulator.convertImplicitToExplicitHydrogens(mol);
    
    }
    
//------------------------------------------------------------------------------    
 
    public static void ensureNoUnsetBondOrdersSilent(IAtomContainer iac)
    {
        try {
            ensureNoUnsetBondOrders(iac);
        } catch (CDKException e) {
            StaticLogger.appLogger.log(Level.WARNING, "Kekulization failed. "
                    + "Bond orders will be unreliable as all unset bonds are"
                    + "now converted to single-order bonds.");
        }
        for (IBond bnd : iac.bonds())
        {
            if (bnd.getOrder().equals(IBond.Order.UNSET)) 
            {
                bnd.setOrder(IBond.Order.SINGLE);
            }
        }
    }
    
//------------------------------------------------------------------------------
 
    public static void ensureNoUnsetBondOrders(IAtomContainer iac) throws CDKException
    {
        Kekulization.kekulize(iac);
        
        for (IBond bnd : iac.bonds())
        {
            if (bnd.getOrder().equals(IBond.Order.UNSET)) 
            {
                bnd.setOrder(IBond.Order.SINGLE);
            }
        }
    }
    
//------------------------------------------------------------------------------
    
    public static String missmatchingAromaticity(IAtomContainer mol)
    {
        String cause = "";
        for (IAtom atm : mol.atoms())
        {
            //Check carbons with or without aromatic flags
            if (atm.getSymbol().equals("C") && atm.getFormalCharge() == 0
                    && mol.getConnectedBondsCount(atm) == 3)
            {
                if (atm.getFlag(CDKConstants.ISAROMATIC))
                {
                    int n = numOfBondsWithBO(atm, mol, IBond.Order.DOUBLE);
                    if (n == 0)
                    {
                        cause = "Aromatic atom " + getAtomRef(atm,mol) 
                        + " has 3 connected atoms but no double bonds";
                        return cause;
                    }
                } else {
                    for (IAtom nbr : mol.getConnectedAtomsList(atm))
                    {
                        if (nbr.getSymbol().equals("C"))
                        {
                            if (nbr.getFormalCharge() == 0)
                            {
                                if (mol.getConnectedBondsCount(nbr) == 3)
                                {
                                    int nNbr = numOfBondsWithBO(nbr, mol, 
                                            IBond.Order.SINGLE);
                                    int nAtm = numOfBondsWithBO(atm, mol, 
                                            IBond.Order.SINGLE);
                                    if ((nNbr == 3) && (nAtm == 3))
                                    {
                                        cause = "Connected atoms "
                                            + getAtomRef(atm, mol) + " " 
                                            + getAtomRef(nbr, mol) 
                                            + " have 3 connected atoms "
                                            + "but no double bond. They are "
                                            + "likely to be aromatic but no "
                                            + "aromaticity has been reported.";
                                        return cause;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        return cause;
    }
    
//------------------------------------------------------------------------------
    
    public static int numOfBondsWithBO(IAtom atm, IAtomContainer mol, 
            IBond.Order order)
    {
        int n = 0;
        for (IBond bnd : mol.getConnectedBondsList(atm))
        {
            if (bnd.getOrder() == order)
                n++;
        }
        return n;
    }
    
//------------------------------------------------------------------------------
    
    public static IAtomContainer makeSameAs(IAtomContainer mol) throws DENOPTIMException
    {
        IChemObjectBuilder builder = SilentChemObjectBuilder.getInstance();
        IAtomContainer iac = builder.newAtomContainer();
        
        for (IAtom oAtm : mol.atoms())
        {
            IAtom nAtm = MoleculeUtils.makeSameAtomAs(oAtm,true,true);
            iac.addAtom(nAtm);
        }
        
        for (IBond oBnd : mol.bonds())
        {
            if (oBnd.getAtomCount() != 2)
            {
                throw new DENOPTIMException("Unable to deal with bonds "
                        + "involving more than two atoms.");
            }
            int ia = mol.indexOf(oBnd.getAtom(0));
            int ib = mol.indexOf(oBnd.getAtom(1));
            iac.addBond(ia,ib,oBnd.getOrder());
        }
        return iac;
    }
    
//------------------------------------------------------------------------------
 
    public static IAtom makeSameAtomAs(IAtom oAtm)
    {
        return makeSameAtomAs(oAtm, false, false);
    }
    
//------------------------------------------------------------------------------
 
    public static IAtom makeSameAtomAs(IAtom oAtm, boolean ignoreValence,
            boolean ignoreImplicitH)
    {
        IAtom nAtm = null;
        String s = getSymbolOrLabel(oAtm);
        if (MoleculeUtils.isElement(oAtm))
        {
            nAtm = new Atom(s);
        } else {
            nAtm = new PseudoAtom(s);
        }
        if (oAtm.getPoint3d() != null)
        {
            Point3d p3d = oAtm.getPoint3d();
            nAtm.setPoint3d(new Point3d(p3d.x, p3d.y, p3d.z));
        } else if (oAtm.getPoint2d() != null)
        {
            Point2d p2d = oAtm.getPoint2d();
            nAtm.setPoint3d(new Point3d(p2d.x, p2d.y, 0.00001));
        }
        if (oAtm.getFormalCharge() != null)
            nAtm.setFormalCharge(oAtm.getFormalCharge());
        if (oAtm.getBondOrderSum() != null)
            nAtm.setBondOrderSum(oAtm.getBondOrderSum());
        if (oAtm.getCharge() != null)
            nAtm.setCharge(oAtm.getCharge());
        if (oAtm.getValency() != null && !ignoreValence)
            nAtm.setValency(oAtm.getValency());
        if (oAtm.getExactMass() != null)
            nAtm.setExactMass(oAtm.getExactMass());
        if (oAtm.getMassNumber() != null)
            nAtm.setMassNumber(oAtm.getMassNumber());
        if (oAtm.getFormalNeighbourCount() != null)
            nAtm.setFormalNeighbourCount(oAtm.getFormalNeighbourCount());
        if (oAtm.getFractionalPoint3d() != null)
            nAtm.setFractionalPoint3d(new Point3d(oAtm.getFractionalPoint3d()));
        if (oAtm.getHybridization() != null)
            nAtm.setHybridization(Hybridization.valueOf(
                    oAtm.getHybridization().toString()));
        if (oAtm.getImplicitHydrogenCount() != null && !ignoreImplicitH)
            nAtm.setImplicitHydrogenCount(oAtm.getImplicitHydrogenCount());
        if (oAtm.getMaxBondOrder() != null)
            nAtm.setMaxBondOrder(oAtm.getMaxBondOrder());
        if (oAtm.getNaturalAbundance() != null)
            nAtm.setNaturalAbundance(oAtm.getNaturalAbundance());
        
        return nAtm;
    }
 
//------------------------------------------------------------------------------
    
    public static String getSymbolOrLabel(IAtom atm)
    {
        String s = "none";
        if (MoleculeUtils.isElement(atm))
        {
            s = atm.getSymbol();
        } else {
            //TODO: one day will account for the possibility of having any 
            // other implementation of IAtom
            IAtom a = null;
            if (atm instanceof AtomRef) 
            {
                a = ((AtomRef) atm).deref();
                if (a instanceof PseudoAtom)
                {
                    s = ((PseudoAtom) a).getLabel();
                } else {
                    // WARNING: we fall back to standard behavior, but there
                    // could still be cases where this is not good...
                    s = a.getSymbol();
                }
            } else if (atm instanceof PseudoAtom)
            {
                s = ((PseudoAtom) atm).getLabel();
            }
        }      
        return s;
    }
    
//------------------------------------------------------------------------------
 
    public static String getAtomRef(IAtom atm, IAtomContainer mol)
    {
        return atm.getSymbol() + (mol.indexOf(atm) +1);
    }
    
//------------------------------------------------------------------------------
 
    public static Point3d calculateCentroid(IAtomContainer mol)
    {
        Point3d c = new Point3d(0,0,0);
        for (IAtom atm : mol.atoms())
        {
            c.add(getPoint3d(atm));
        }
        c.scale(1.0 / ((double) mol.getAtomCount()));
        return c;
    }
    
//------------------------------------------------------------------------------
    
    public static IAtomContainer extractIACForSubgraph(IAtomContainer wholeIAC, 
            DGraph subGraph, DGraph wholeGraph, Logger logger, 
            Randomizer randomizer) throws DENOPTIMException
    {
        IAtomContainer iac = makeSameAs(wholeIAC);
        
        Set<Long> wantedVIDs = new HashSet<Long>();
        Map<Long,Vertex> wantedVertexesMap = new HashMap<>();
        for (Vertex v : subGraph.getVertexList())
        {
            Object o = v.getProperty(DENOPTIMConstants.STOREDVID);
            if (o == null)
            {
                throw new DENOPTIMException("Property '" 
                        + DENOPTIMConstants.STOREDVID + "' not defined in "
                        + "vertex " + v + ", but is needed to extract "
                                + "substructure.");
            }
            wantedVIDs.add(((Long) o).longValue());
            wantedVertexesMap.put(((Long) o).longValue(), v);
        }
        
        // Identify the destiny of each atom: keep, remove, or make AP from it.
        IAtomContainer toRemove = new AtomContainer();
        Map<IAtom,IAtom> toAP = new HashMap<IAtom,IAtom>();
        Map<IAtom,AttachmentPoint> mapAtmToAPInG = 
                new HashMap<IAtom,AttachmentPoint>();
        Map<IAtom,APClass> apcMap = new HashMap<IAtom,APClass>();
        for (int i=0; i<wholeIAC.getAtomCount(); i++)
        {
            IAtom cpAtm = iac.getAtom(i);
            IAtom oriAtm = wholeIAC.getAtom(i);
            Object o = oriAtm.getProperty(DENOPTIMConstants.ATMPROPVERTEXID);
            if (o == null)
            {
                throw new DENOPTIMException("Property '" 
                        + DENOPTIMConstants.ATMPROPVERTEXID 
                        + "' not defined in atom "
                        + oriAtm.getSymbol() + wholeIAC.indexOf(oriAtm) 
                        + ", but is needed to extract substructure.");
            }
            long vid = ((Long) o).longValue();
            if (wantedVIDs.contains(vid))
            {
                continue; //keep this atom cpAtm
            }
            // Now, decide if the current atom should become an AP
            boolean willBecomeAP = false;
            for (IAtom nbr : wholeIAC.getConnectedAtomsList(oriAtm))
            {
                Object oNbr = nbr.getProperty(DENOPTIMConstants.ATMPROPVERTEXID);
                if (oNbr == null)
                {
                    throw new DENOPTIMException("Property '" 
                            + DENOPTIMConstants.ATMPROPVERTEXID 
                            + "' not defined in atom "
                            + nbr.getSymbol() + wholeIAC.indexOf(nbr) 
                            + ", but is needed to extract substructure.");
                }
                long nbrVid = ((Long) oNbr).longValue();
                if (wantedVIDs.contains(nbrVid))
                {
                    // cpAtm is connected to an atom to keep and will therefore
                    // become an AP. Note that the connection may go through
                    // a chord in the graph, i.e., a pairs of RCVs!
                    toAP.put(cpAtm,iac.getAtom(wholeIAC.indexOf(nbr)));
                    AttachmentPoint apInWholeGraph = 
                            wholeGraph.getAPOnLeftVertexID(nbrVid,vid);
                    if (apInWholeGraph == null)
                    {
                        String debugFile = "failedAPIdentificationIACSubGraph" 
                                + wholeGraph.getGraphId() + ".sdf";
                        DenoptimIO.writeGraphToSDF(new File(debugFile), 
                                wholeGraph, willBecomeAP, logger, randomizer);
                        throw new DENOPTIMException("Unmexpected null AP from "
                                + nbrVid + " " + vid +" on " + wholeGraph 
                                + " See " + debugFile);
                    }
                    AttachmentPoint apInSubGraph = 
                            wantedVertexesMap.get(nbrVid).getAP(apInWholeGraph.getIndexInOwner());
                    mapAtmToAPInG.put(cpAtm, apInSubGraph);
                    APClass apc = apInWholeGraph.getAPClass();
                    apcMap.put(cpAtm, apc);
                    willBecomeAP = true;
                    break;
                }
            }
            if (!willBecomeAP)
                toRemove.addAtom(cpAtm);
        }
        
        iac.remove(toRemove);
        
        // NB: the molecular representation in frag is NOT iac! It's a clone of it
        Fragment frag = new Fragment(iac,BBType.FRAGMENT);
        
        List<IAtom> atmosToRemove = new ArrayList<>();
        List<IBond> bondsToRemove = new ArrayList<>();
        for (IAtom trgAtmInIAC : toAP.keySet())
        {
            IAtom trgAtm = frag.getAtom(iac.indexOf(trgAtmInIAC));
            IAtom srcAtmInISC = toAP.get(trgAtmInIAC);
            IAtom srcAtm = frag.getAtom(iac.indexOf(srcAtmInISC));
            
            AttachmentPoint apInG = mapAtmToAPInG.get(trgAtmInIAC);
            
            // Make Attachment point
            Point3d srcP3d = MoleculeUtils.getPoint3d(srcAtm);
            Point3d trgP3d = MoleculeUtils.getPoint3d(trgAtm);
            double currentLength = srcP3d.distance(trgP3d);
            //TODO-V3+? change hard-coded value with property of AP, when such
            // property will be available, i. e., one refactoring of AP and 
            // atom coordinates is done.
            double idealLength = 1.53;
            /*
            double idealLength = apInG.getProperty(
                    DENOPTIMConstants.APORIGINALLENGTH);
                    */
            Point3d vector = new Point3d();
            vector.x = srcP3d.x + (trgP3d.x - srcP3d.x)*(idealLength/currentLength);
            vector.y = srcP3d.y + (trgP3d.y - srcP3d.y)*(idealLength/currentLength);
            vector.z = srcP3d.z + (trgP3d.z - srcP3d.z)*(idealLength/currentLength);
            
            AttachmentPoint createdAP = frag.addAPOnAtom(srcAtm, 
                    apcMap.get(trgAtmInIAC), vector);
            createdAP.setProperty(DENOPTIMConstants.LINKAPS, apInG);
            
            for (IBond bnd : frag.bonds())
            {
                if (bnd.contains(trgAtm))
                {
                    bondsToRemove.add(bnd);
                }
            }
            atmosToRemove.add(trgAtm);
        }
        
        // Remove atom that has become an AP and associated bonds
        for (IBond bnd : bondsToRemove)
        {
            frag.removeBond(bnd);
        }
        for (IAtom a : atmosToRemove)
        {
            frag.removeAtom(a);
        }
        
        frag.updateAPs();
        
        return frag.getIAtomContainer();
    }
    
//------------------------------------------------------------------------------
 
    public static int getDimensions(IAtomContainer mol)
    {
        final int is2D = 2;
        final int is3D = 3;
        final int not2or3D = -1;
 
        int numOf2D = 0;
        int numOf3D = 0;
 
        for (IAtom atm : mol.atoms())
        {
            Point2d p2d = new Point2d();
            Point3d p3d = new Point3d();
            p2d = atm.getPoint2d();
            boolean have2D = true;
            if (p2d == null)
            {
                have2D = false;
                p3d = atm.getPoint3d();
                if (p3d == null)
                {
                    return not2or3D;
                }
            }
            ArrayList<Double> pointer = new ArrayList<Double>();
            try {
                if (have2D)
                {
                    pointer.add(p2d.x);
                    pointer.add(p2d.y);
                    numOf2D++;
                } else {
                    pointer.add(p3d.x);
                    pointer.add(p3d.y);
                    pointer.add(p3d.z);
                    numOf3D++;
                }
            } catch (Throwable t) {
                return not2or3D;
            }
        }
 
        if (numOf2D == mol.getAtomCount())
            return is2D;
        else if (numOf3D == mol.getAtomCount())
            return is3D;
        else
            return not2or3D;
    }
 
//------------------------------------------------------------------------------
 
}