ChemicalObjectModel.java

Go to the documentation of this file.

/*
 *   DENOPTIM
 *   Copyright (C) 2019 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.molecularmodeling;
 
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.graph.PathTools;
import org.openscience.cdk.graph.SpanningTree;
import org.openscience.cdk.interfaces.IAtom;
import org.openscience.cdk.interfaces.IAtomContainer;
import org.openscience.cdk.interfaces.IBond;
import org.openscience.cdk.interfaces.IChemObjectBuilder;
import org.openscience.cdk.interfaces.IRingSet;
import org.openscience.cdk.silent.RingSet;
import org.openscience.cdk.silent.SilentChemObjectBuilder;
 
import denoptim.constants.DENOPTIMConstants;
import denoptim.exception.DENOPTIMException;
import denoptim.graph.DGraph;
import denoptim.graph.Edge.BondType;
import denoptim.graph.Ring;
import denoptim.graph.rings.RingClosingAttractor;
import denoptim.graph.rings.RingClosure;
import denoptim.integration.tinker.TinkerAtom;
import denoptim.integration.tinker.TinkerMolecule;
import denoptim.integration.tinker.TinkerUtils;
import denoptim.io.DenoptimIO;
import denoptim.utils.MathUtils;
import denoptim.utils.ObjectPair;
 
public class ChemicalObjectModel
{
    private DGraph molGraph;
 
    private IAtomContainer fmol;
 
    private TinkerMolecule tmol;
 
    private String molName;
 
    private ArrayList<RingClosingAttractor> attractors;
 
    private Map<RingClosingAttractor,Integer> attToAtmID;
 
    private ArrayList<Set<ObjectPair>> allRCACombs;
 
    private ArrayList<ObjectPair> rotatableBnds;
 
    private ArrayList<RingClosure> newRingClosures;
 
    private double overalRCScore = Double.NaN;
 
    private double atmOveralScore = Double.NaN;
 
    private Logger logger;
    
    private ArrayList<Integer> oldToNewOrder;
    private ArrayList<Integer> newToOldOrder;
 
//------------------------------------------------------------------------------
 
    public ChemicalObjectModel()
    {
        this.molGraph = new DGraph();
        IChemObjectBuilder builder = SilentChemObjectBuilder.getInstance();
        this.fmol = builder.newAtomContainer();
        this.tmol = new TinkerMolecule();
        this.attractors = new ArrayList<RingClosingAttractor>();
        this.attToAtmID = new HashMap<RingClosingAttractor,Integer>();
        this.allRCACombs = new ArrayList<Set<ObjectPair>>();
        this.rotatableBnds = new ArrayList<ObjectPair>();
        this.newRingClosures = new ArrayList<RingClosure>();
        this.molName = "none";
        this.oldToNewOrder = new ArrayList<Integer>();
        this.newToOldOrder =  new ArrayList<Integer>();
    }
 
//------------------------------------------------------------------------------
 
    public ChemicalObjectModel(DGraph molGraph, 
            IAtomContainer fmol, 
            TinkerMolecule tmol, 
            String molName, 
            ArrayList<ObjectPair> rotatableBnds,
            ArrayList<RingClosingAttractor> attractors,
            Map<RingClosingAttractor,Integer> attToAtmID,
            ArrayList<Set<ObjectPair>> allRCACombs,
            ArrayList<RingClosure> ringClosures,
            ArrayList<Integer> oldToNewOrder,
            ArrayList<Integer> newToOldOrder,
            Logger logger)
    {
        this.molGraph = molGraph;
        this.fmol = fmol;
        this.tmol = tmol;
        this.attractors = attractors;
        this.attToAtmID = attToAtmID;
        this.allRCACombs = allRCACombs;
        this.rotatableBnds = rotatableBnds;
        this.molName = molName;
        this.newRingClosures = ringClosures;
        this.overalRCScore = Double.NaN;
        this.atmOveralScore = Double.NaN;
        if (this.attractors.size() != 0)
        {
            if (molGraph.hasOrEmbedsRings())
            {
                // ring closures defined in the input
                convertRingsToRCACombinations();
            }
        }
        this.oldToNewOrder = oldToNewOrder;
        this.newToOldOrder = newToOldOrder;
        this.logger = logger;
    }
 
//------------------------------------------------------------------------------
 
    public ChemicalObjectModel(DGraph molGraph, IAtomContainer fmol, 
            TinkerMolecule tmol, String molName,
            ArrayList<ObjectPair> rotatableBnds,
            ArrayList<Integer> oldToNewOrder,
            ArrayList<Integer> newToOldOrder,
            Logger logger) throws DENOPTIMException
    {
        this.molGraph = molGraph;
        this.fmol = fmol;
        this.tmol = tmol;
        updateXYZFromINT();
        this.molName = molName;
        this.rotatableBnds = rotatableBnds;
        this.molName = molName;
        this.attractors = new ArrayList<RingClosingAttractor>();
        this.attToAtmID = new HashMap<RingClosingAttractor,Integer>();
        findAttractors();
        this.allRCACombs = new ArrayList<Set<ObjectPair>>();
        if (this.attractors.size() != 0)
        {
            if (molGraph.hasOrEmbedsRings())
            {
                // ring closures defined in the input
                convertRingsToRCACombinations();
            }
        }
        this.newRingClosures = new ArrayList<RingClosure>();
        this.overalRCScore = Double.NaN;
        this.atmOveralScore = Double.NaN;
        this.oldToNewOrder = oldToNewOrder;
        this.newToOldOrder = newToOldOrder;
        this.logger = logger;
    }
    
//------------------------------------------------------------------------------
    
    public List<Integer> getOldToNewOrder()
    {
        return oldToNewOrder;
    }
    
//------------------------------------------------------------------------------
    
    public List<Integer> getNewToOldOrder()
    {
        return newToOldOrder;
    }
    
//------------------------------------------------------------------------------    
 
    private void findAttractors()
    {
        // Identify all RingClosingAttractors
        for (IAtom atm : fmol.atoms())
        {
            RingClosingAttractor rca = new RingClosingAttractor(atm,fmol);
            if (rca.isAttractor())
            {
                attractors.add(rca);
                attToAtmID.put(rca, fmol.indexOf(atm));
            }
        }
    }
 
//------------------------------------------------------------------------------
 
    private void convertRingsToRCACombinations()
    {
        //When we give the a graph with rings, we are defining the RCA 
        //combinations, so there will be only one RCA combination
        Set<ObjectPair> singleRCAcomb = new HashSet<ObjectPair>();
        
        for (int i=0; i<attractors.size(); i++)
        {
            RingClosingAttractor rcaI = attractors.get(i);
            Ring ringOwnerI = rcaI.getRingUser();
            if (ringOwnerI==null)
                continue;
            
            for (int j=i+1; j<attractors.size(); j++)
            {
                RingClosingAttractor rcaJ = attractors.get(j);
                Ring ringOwnerJ = rcaJ.getRingUser();
                if (ringOwnerJ==null)
                    continue;
                if (ringOwnerI==ringOwnerJ)
                {
                    singleRCAcomb.add(new ObjectPair(rcaI, rcaJ));
                }
            }
        }
        allRCACombs.add(singleRCAcomb);
    }
 
//------------------------------------------------------------------------------
 
    public void updateXYZFromINT() throws DENOPTIMException
    {
        ArrayList<Point3d> newCoords = new ArrayList<Point3d>();
        for (int i=0; i<fmol.getAtomCount(); i++)
        {
            int iTnk = i + 1;
            TinkerAtom tAtm = tmol.getAtom(iTnk);
            int ia = tAtm.getAtomNeighbours()[0];
            int ib = tAtm.getAtomNeighbours()[1];
            int ic = tAtm.getAtomNeighbours()[2];
            if (ia >  i || ib > i || ic > i)
            {
                String msg = "ERROR! Cannot convert internal coordinates of "
                                   + "atom " + i + " " + tAtm + ". Reference "
                                   + "atoms have higher ID. The atoms list in "
                                   + "Molecule3DBuilder needs to be reordered.";
                throw new DENOPTIMException(msg);
            }
            int angleFlag = tAtm.getAtomNeighbours()[3];
            double bond = tAtm.getDistAngle()[0];
            double angle1 = tAtm.getDistAngle()[1] * 2 * Math.PI / 360.0;
            double angle2 = tAtm.getDistAngle()[2] * 2 * Math.PI / 360.0;
            double[] newXYZ = {0.0, 0.0, 0.0};
            double smallDble = DENOPTIMConstants.FLOATCOMPARISONTOLERANCE;
            double s1 = Math.sin(angle1); 
            double s2 = Math.sin(angle2);
            double c1 = Math.cos(angle1);
            double c2 = Math.cos(angle2);
            if (ia == 0)
            {
                //first atom remains on the origin
            }
            else if (ib == 0)
            {
                newXYZ[2] = bond;
            }
            else if (ic == 0)
            {
                Vector3d nab = MathUtils.normDist(
                                       newCoords.get(ia-1),newCoords.get(ib-1));
                double rab = MathUtils.distance(
                                       newCoords.get(ia-1),newCoords.get(ib-1));
                newXYZ[0] = bond * s1;
                newXYZ[2] = newCoords.get(ib-1).z + (rab - bond*c1)*nab.z;
            }
            else if (angleFlag == 0)
            {
                Vector3d nab = MathUtils.normDist(
                                       newCoords.get(ia-1),newCoords.get(ib-1));
                Vector3d nbc = MathUtils.normDist(
                                       newCoords.get(ib-1),newCoords.get(ic-1));
                Vector3d t = new Vector3d();
                t.cross(nbc,nab);
                double c = nab.x*nbc.x + nab.y*nbc.y + nab.z*nbc.z;
                if (Math.abs(c) > (1.0 - smallDble))
                {
                    String msg = "ERROR! Linearity does not allow definition "
                        + "of the dihedral angle for atom " + i + " " + tAtm 
                        + " (Value of c="+c+" too close to unity; threshold is "
                        + (1.0 - smallDble) + "). "
                        + "You better use dummy atoms to avoid linearities.";
                    throw new DENOPTIMException(msg);
                }
                t.scale(1/Math.sqrt(Math.max(1.00 - c*c, smallDble)));
                Vector3d u = new Vector3d();
                u.cross(t,nab);
                newXYZ[0] = newCoords.get(ia-1).x
                                      + bond*(u.x*s1*c2 + t.x*s1*s2 - nab.x*c1);
                newXYZ[1] = newCoords.get(ia-1).y
                                      + bond*(u.y*s1*c2 + t.y*s1*s2 - nab.y*c1);
                newXYZ[2] = newCoords.get(ia-1).z
                                      + bond*(u.z*s1*c2 + t.z*s1*s2 - nab.z*c1);
            }
            else if (Math.abs(angleFlag) == 1)
            {
                Vector3d nba = MathUtils.normDist(
                                       newCoords.get(ib-1),newCoords.get(ia-1));
                Vector3d nac = MathUtils.normDist(
                                       newCoords.get(ia-1),newCoords.get(ic-1));
                Vector3d t = new Vector3d();
                t.cross(nac,nba);
                double c = nba.x*nac.x + nba.y*nac.y + nba.z*nac.z;
                if (Math.abs(c) > (1.0 - smallDble))
                {
                    logger.log(Level.WARNING, "WARNING! close-to-linear system "
                        + "in the definition of atom " + i + " " + tAtm + ". "
                        + "You better use dummy atoms to avoid linearities.");
                }
                double s = Math.max(1.00 - c*c, smallDble);
                double a = (-c2 - c*c1) / s;
                double b = (c1 + c*c2) / s;
                double ci = (1.00 + a*c2 - b*c1) / s;
                if (ci > smallDble)
                {
                    ci = angleFlag * Math.sqrt(ci);
                } 
                else if (ci < -smallDble)
                {
                    ci = Math.sqrt((a*nac.x+b*nba.x)*(a*nac.x+b*nba.x) 
                                + (a*nac.y+b*nba.y)*(a*nac.y+b*nba.y)
                                + (a*nac.z+b*nba.z)*(a*nac.z+b*nba.z));
                    a = a / ci;
                    b = b / ci;
                    ci = 0.0;
                    logger.log(Level.WARNING, "WARNING! close-to-linear system "
                        + "in the definition of atom " + i + " " + tAtm + ". "
                        + "You better use dummy atoms to avoid linearities.");
                }
                else 
                {
                   ci = 0.0;
                }
                newXYZ[0] = newCoords.get(ia-1).x 
                                            + bond*(a*nac.x + b*nba.x + ci*t.x);
                newXYZ[1] = newCoords.get(ia-1).y 
                                            + bond*(a*nac.y + b*nba.y + ci*t.y);
                newXYZ[2] = newCoords.get(ia-1).z 
                                            + bond*(a*nac.z + b*nba.z + ci*t.z);
            }
            Point3d p3d = new Point3d(newXYZ[0],newXYZ[1],newXYZ[2]);
            newCoords.add(p3d);
        }
 
        // Update Cartesian coordinates
        for (int i=0; i<fmol.getAtomCount(); i++)
        {
            fmol.getAtom(i).setPoint3d(newCoords.get(i));
            tmol.getAtom(i+1).moveTo(newCoords.get(i).x, 
                                     newCoords.get(i).y,
                                     newCoords.get(i).z);
        }
    }
 
//------------------------------------------------------------------------------
 
    public DGraph getGraph()
    {
        return molGraph;
    }
 
//------------------------------------------------------------------------------
 
    public IAtomContainer getIAtomContainer()
    {
        return fmol;
    }
 
//------------------------------------------------------------------------------
 
    public TinkerMolecule getTinkerMolecule()
    {
        return tmol;
    }
 
//------------------------------------------------------------------------------
 
    public ArrayList<RingClosingAttractor> getAttractorsList()
    {
        return attractors;
    }
 
//------------------------------------------------------------------------------
 
    public RingClosingAttractor getAttractor(int i)
    {
        return attractors.get(i);
    }
 
//------------------------------------------------------------------------------
 
    public int getAtmIdOfRCA(RingClosingAttractor rca)
    {
        return attToAtmID.get(rca);
    }
 
//------------------------------------------------------------------------------
 
    public int getTnkAtmIdOfRCA(RingClosingAttractor rca)
    {
        return attToAtmID.get(rca) + 1;
    }
 
//------------------------------------------------------------------------------
 
    public ArrayList<Set<ObjectPair>> getRCACombinations()
    {
        return allRCACombs;
    }
 
//------------------------------------------------------------------------------
 
    public ArrayList<ObjectPair> getRotatableBonds()
    {
        return rotatableBnds;
    }
 
//------------------------------------------------------------------------------
 
    public int getNumberRotatableBonds()
    {
        return rotatableBnds.size();
    }
 
//------------------------------------------------------------------------------
 
    public ArrayList<RingClosure> getNewRingClosures()
    {
        return newRingClosures;
    }
 
//------------------------------------------------------------------------------
 
    public double getNewRingClosuresQuality()
    {
        if (Double.isNaN(overalRCScore))
        {
            double score = 0.0;
            for (RingClosure rc : newRingClosures)
            {
                score = score + rc.getRingClosureQuality();
            }
            overalRCScore = score;
        }
        return overalRCScore;
    }
 
//------------------------------------------------------------------------------
 
    public double getAtomOverlapScore()
    {
        if (Double.isNaN(atmOveralScore))
        {
            double score = 0.0;
            for (IAtom atmA : fmol.atoms())
            {
                String elA = atmA.getSymbol();
 
                if (DENOPTIMConstants.DUMMYATMSYMBOL.equals(elA))
                    continue;
 
                IAtom[] toExclude = PathTools.findClosestByBond(fmol,atmA,4);
                for (IAtom atmB : fmol.atoms())
                {
                    if (atmA == atmB)
                        continue;
 
                    if (Arrays.asList(toExclude).contains(atmB))
                        continue;
 
                    String elB = atmB.getSymbol();
                    if (DENOPTIMConstants.DUMMYATMSYMBOL.equals(elB))
                        continue;
 
                    double dist = atmA.getPoint3d().distance(atmB.getPoint3d());
                    score = score + dist;
                }
            }
            atmOveralScore = score;
        }
        return atmOveralScore;
    }
 
//------------------------------------------------------------------------------
 
    public String getName()
    {
        return molName;
    }
 
//------------------------------------------------------------------------------
 
    public void addBond(IAtom atmA, IAtom atmB, RingClosure nRc, 
            BondType bndTyp)
    {
        this.newRingClosures.add(nRc);
        this.overalRCScore = Double.NaN;
 
        int iA = fmol.indexOf(atmA);
        int iB = fmol.indexOf(atmB);
       
        if (bndTyp.hasCDKAnalogue())
        {
            this.fmol.addBond(iA, iB, bndTyp.getCDKOrder());
            logger.log(Level.FINE, "ADDING BOND: "+iA+" "+iB);
        } else {
            logger.log(Level.FINE, "WARNING! Attempt to add ring closing bond "
                    + "did not add any actual chemical bond because the "
                    + "bond type of the chord is '" + bndTyp +"'.");
        }
 
        if (iA < iB)
        {
            this.tmol.addBond(iA+1, iB+1);
        }
        else
        {
            this.tmol.addBond(iB+1, iA+1);
        }
    }
 
//------------------------------------------------------------------------------
 
    public void purgeListRotatableBonds() throws DENOPTIMException
    {
        //Get all rings
        SpanningTree st = new SpanningTree(fmol);
        IRingSet allRings = new RingSet();
        try {
            allRings = st.getAllRings();
        } catch (Exception ex) {
            throw new DENOPTIMException(ex);
        }
 
        //Identify bonds to remove (cyclic bonds)
        ArrayList<ObjectPair> toRemove = new ArrayList<ObjectPair>();
        for (ObjectPair op : rotatableBnds)
        {
            int i1 = ((Integer)op.getFirst()).intValue();
            int i2 = ((Integer)op.getSecond()).intValue();
 
            IBond bnd = fmol.getBond(fmol.getAtom(i1), fmol.getAtom(i2));
            IRingSet rs = allRings.getRings(bnd);
            if (!rs.isEmpty())
            {
                toRemove.add(op);
                logger.log(Level.FINE, "Bond " + i1 + "-" + i2 + " (TnkID:"
                    + (i1+1) + "-" + (i2+1) + ") is not rotatable anymore");
            }
        }
 
        //Remove bonds
        for (ObjectPair op : toRemove)
        {
            rotatableBnds.remove(op);
        }        
    }
 
//------------------------------------------------------------------------------
 
    @SuppressWarnings("unchecked")
    public ChemicalObjectModel deepcopy() throws DENOPTIMException
    {
        String nMolName = this.molName;
        DGraph nMolGraph = this.molGraph.clone();
        IAtomContainer nFMol;
        TinkerMolecule nTMol;
        ArrayList<ObjectPair> nRotBnds;
 
        try 
        {
            nFMol = this.fmol.clone();
            nTMol = (TinkerMolecule) this.tmol.deepCopy();
            nRotBnds = (ArrayList<ObjectPair>) this.rotatableBnds.clone();
        }
        catch (CloneNotSupportedException cns) 
        {
            throw new DENOPTIMException(cns);
        }
 
        ArrayList<RingClosingAttractor> nAttractors = 
                new ArrayList<RingClosingAttractor>();
        Map<RingClosingAttractor,Integer> nAttToAtmID = 
                new HashMap<RingClosingAttractor,Integer>();
        Map<RingClosingAttractor,RingClosingAttractor> oldToNewRCA =
                new HashMap<RingClosingAttractor,RingClosingAttractor>();
        for (int iorca=0; iorca<attractors.size(); iorca++)
        {
            RingClosingAttractor oRca = attractors.get(iorca);
            int ioatm = this.getAtmIdOfRCA(oRca);
            IAtom atm = nFMol.getAtom(ioatm);
            RingClosingAttractor nRca = new RingClosingAttractor(atm,nFMol);
            nAttractors.add(nRca);
            nAttToAtmID.put(nRca, ioatm);
            oldToNewRCA.put(oRca, nRca);
        }
 
        ArrayList<Set<ObjectPair>> nAllRCACombs = 
                new ArrayList<Set<ObjectPair>>();
        for (Set<ObjectPair> sop : this.allRCACombs)
        {
            Set<ObjectPair> nSop = new HashSet<ObjectPair>();
            for (ObjectPair op : sop)
            {
                ObjectPair nOp = new ObjectPair(
                        oldToNewRCA.get(op.getFirst()), 
                        oldToNewRCA.get(op.getSecond()));
                nSop.add(nOp);
            }
            nAllRCACombs.add(nSop);
        }
 
        ArrayList<RingClosure> nNewRingClosures = new ArrayList<RingClosure>();
        for (RingClosure rc : this.newRingClosures)
        {
            nNewRingClosures.add(rc.deepCopy());
        }
        
        ArrayList<Integer> newOldToNewOrder = new ArrayList<Integer>();
        for (Integer i : oldToNewOrder)
            newOldToNewOrder.add(i.intValue());
        
        ArrayList<Integer> newNewToOldOrder = new ArrayList<Integer>();
        for (Integer i : newToOldOrder)
            newNewToOldOrder.add(i.intValue());
        
        ChemicalObjectModel molClone = new ChemicalObjectModel(nMolGraph,
                nFMol,
                nTMol,
                nMolName,
                nRotBnds,
                nAttractors,
                nAttToAtmID,
                nAllRCACombs,
                nNewRingClosures,
                newOldToNewOrder,
                newNewToOldOrder, logger);
 
        return molClone;
    }
 
//------------------------------------------------------------------------------
}