TinkerUtils.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.integration.tinker;
 
import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
 
import org.apache.commons.io.input.ReversedLinesFileReader;
import org.openscience.cdk.interfaces.IAtom;
import org.openscience.cdk.interfaces.IAtomContainer;
import org.openscience.cdk.interfaces.IBond;
 
import denoptim.constants.DENOPTIMConstants;
import denoptim.exception.DENOPTIMException;
import denoptim.molecularmodeling.MMBuilderUtils;
import denoptim.utils.ConnectedLigand;
import denoptim.utils.ConnectedLigandComparator;
import denoptim.utils.GeneralUtils;
import denoptim.utils.MathUtils;
import denoptim.utils.MoleculeUtils;
import denoptim.utils.ObjectPair;
import denoptim.utils.RotationalSpaceUtils;
 
 
public class TinkerUtils 
{
    private static final String NL = System.getProperty("line.separator");
    private static boolean debug = false;
    
    
//------------------------------------------------------------------------------    
    public static TinkerMolecule readTinkerIC(String filename)
                                                      throws DENOPTIMException
    {
        BufferedReader br = null;
        String line = null;
        int natom = 0;
 
        TinkerMolecule tinkermol = new TinkerMolecule();
 
        try
        {
            br = new BufferedReader(new FileReader(filename));
            line = br.readLine().trim();
 
            // read the number of atoms and store
            String[] arr = line.split(" +");
            natom = Integer.parseInt(arr[0]);
            if (natom < 1)
            {
                String msg = "Invalid number of atoms in " + filename;
                throw new DENOPTIMException(msg);
            }
 
            if (arr.length >= 2)
            {
                arr = line.split(" +", 2);
                // set the name
                //System.err.println("Name: " + arr[1]);
                tinkermol.setName(arr[1]);
            }
 
 
            ArrayList<int[]> zadd = new ArrayList<>();
            ArrayList<int[]> zdel = new ArrayList<>();
 
            ArrayList<TinkerAtom> lstAtom = new ArrayList<>();
 
            for (int i = 0; i < natom; i++)
            {
                line = br.readLine();
                if (line == null)
                {
                    break;
                }
 
                arr = line.trim().split(" +");
                if (arr == null || arr.length < 3)
                {
                    String msg = "Check atom " + (i + 1) + " in " + filename;
                    throw new DENOPTIMException(msg);
                }
 
                // Atom number, name, type
                String name = arr[1];
                int type = Integer.parseInt(arr[2]);
                double zv[] = new double[3];
                int zi[] = new int[4];
                double d[] =
                {
                    0.0d, 0.0d, 0.0d
                };
 
                // Bond partner and bond value
                if (arr.length >= 5)
                {
                    zi[0] = Integer.parseInt(arr[3]);
                    zv[0] = Double.parseDouble(arr[4]);
                }
                else
                {
                    zi[0] = 0;
                    zv[0] = 0.0d;
                }
                // Angle partner and angle value
                if (arr.length >= 7)
                {
                    zi[1] = Integer.parseInt(arr[5]);
                    zv[1] = Double.parseDouble(arr[6]);
                }
                else
                {
                    zi[1] = 0;
                    zv[1] = 0.0d;
                }
                // Torsion partner and dihedral value
                if (arr.length >= 10)
                {
                    zi[2] = Integer.parseInt(arr[7]);
                    zv[2] = Double.parseDouble(arr[8]);
                    zi[3] = Integer.parseInt(arr[9]);
                }
                else
                {
                    zi[2] = 0;
                    zv[2] = 0.0d;
                    zi[3] = 0;
                }
                TinkerAtom tkatm = new TinkerAtom(i + 1,name,type,d,zi,zv);
                lstAtom.add(tkatm);
            } // end for
 
            if (br.ready())
            {
                line = br.readLine();
                // Check for a first blank line
                if (line.trim().equalsIgnoreCase(""))
                {
                    // Parse bond pairs to add until EOF or a blank line is
                    // reached
                    boolean blank = false;
                    while (br.ready() && !blank)
                    {
                        line = br.readLine();
                        if (line.trim().equalsIgnoreCase(""))
                        {
                            blank = true;
                        }
                        else
                        {
                            arr = line.trim().split(" +");
                            if (arr.length != 2)
                            {
                                String msg = "Check Bond Pair to Remove: " 
                                    + (zadd.size() + 1) + " in " + filename;
                                throw new DENOPTIMException(msg);
                            }
                            int pair[] = new int[2];
                            pair[0] = Integer.parseInt(arr[0]);
                            pair[1] = Integer.parseInt(arr[1]);
                            zadd.add(pair);
                        }
                    }
                    // Parse bond pairs to be removed until EOF
                    while (br.ready())
                    {
                        line = br.readLine();
                        arr = line.trim().split(" +");
                        if (arr.length != 2)
                        {
                            String msg = "Check Bond Pair to Remove: " 
                                    + (zadd.size() + 1) + " in " + filename;
                            throw new DENOPTIMException(msg);
                        }
                        int pair[] = new int[2];
                        pair[0] = Integer.parseInt(arr[0]);
                        pair[1] = Integer.parseInt(arr[1]);
                        zdel.add(pair);
                    }
                }
            }
 
            tinkermol.setAtoms(lstAtom);
            tinkermol.setBondPairs(zdel, zadd);
        }
        catch (NumberFormatException | IOException nfe)
        {
            throw new DENOPTIMException(nfe);
        }
        finally
        {
            try
            {
                if (br != null)
                {
                    br.close();
                }
            }
            catch (IOException ioe)
            {
                throw new DENOPTIMException(ioe);
            }
        }
 
        return tinkermol;
    }
 
//------------------------------------------------------------------------------
    
    public static ArrayList<double[]> readTinkerXYZ(String filename)
                                                        throws DENOPTIMException
    {
        ArrayList<double[]> coords = new ArrayList<>();
 
        BufferedReader br = null;
        String line = null;
 
        int natom = 0;
 
        try
        {
            br = new BufferedReader(new FileReader(filename));
            line = br.readLine().trim();
 
            // read the number of atoms and store
            String[] arr = line.split("\\s+");
            natom = Integer.parseInt(arr[0]);
            if (natom < 1)
            {
               String msg = "Invalid number of atoms in " + filename;
               throw new DENOPTIMException(msg);
            }
 
            while ((line = br.readLine()) != null)
            {
                if (line.trim().length() == 0)
                {
                    continue;
                }
 
                arr = line.trim().split("\\s+");
                if (arr.length >= 5)
                {
                    // Atom number, name, type
                    //String name = arr[1];
                    double[] f = new double[3];
                    f[0] = Double.parseDouble(arr[2]);
                    f[1] = Double.parseDouble(arr[3]);
                    f[2] = Double.parseDouble(arr[4]);
 
                    coords.add(f);
                }
            }
        }
        catch (NumberFormatException | IOException nfe)
        {
            throw new DENOPTIMException(nfe);
        }
        finally
        {
            try
            {
                if (br != null)
                {
                    br.close();
                }
            }
            catch (IOException ioe)
            {
                throw new DENOPTIMException(ioe);
            }
        }
 
        if (coords.size() != natom)
        {
            throw new DENOPTIMException("Incorrect number of atoms perceived "
                    + "in " + filename);
        }
 
        return coords;
 
    }
 
//------------------------------------------------------------------------------   
    
    public static void writeIC(String filename, TinkerMolecule tmol)
            throws DENOPTIMException
    {
        FileWriter fw = null;
        try
        {
            fw = new FileWriter(new File(filename));
            int numatoms = tmol.getAtoms().size();
            // write out the number of atoms and the title
            String header = String.format("%6d  %s%n", numatoms, tmol.getName());
            fw.write(header);
            fw.flush();
 
            String line = "";
 
            for (int i = 0; i < numatoms; i++)
            {
                TinkerAtom atom = tmol.getAtoms().get(i);
                int[] d1 = atom.getAtomNeighbours();
                double[] d2 = atom.getDistAngle();
                // output of first three atoms is handled separately
                if (i == 0)
                {
                    line = String.format("%6d  %-3s%6d%n",
                            atom.getXYZIndex(), atom.getAtomString(),
                            atom.getAtomType());
                    fw.write(line);
                    fw.flush();
                }
                else
                {
                    if (i == 1)
                    {
                        line = String.format("%6d  %-3s%6d%6d%10.5f%n",
                                atom.getXYZIndex(), atom.getAtomString(),
                                atom.getAtomType(), d1[0], d2[0]);
                        fw.write(line);
                        fw.flush();
                    }
                    else
                    {
                        if (i == 2)
                        {
                            line = String.format("%6d  %-3s%6d%6d%10.5f%6d%10.4f%n",
                                    atom.getXYZIndex(), atom.getAtomString(),
                                    atom.getAtomType(), d1[0], d2[0],
                                    d1[1], d2[1]);
                            fw.write(line);
                            fw.flush();
                        } // output the fourth through final atoms
                        else
                        {
                            line = String.format("%6d  %-3s%6d%6d%10.5f%6d%10.4f%6d%10.4f%6d%n",
                                    atom.getXYZIndex(), atom.getAtomString(),
                                    atom.getAtomType(), d1[0], d2[0],
                                    d1[1], d2[1], d1[2], d2[2], d1[3]);
                            fw.write(line);
                            fw.flush();
                        }
                    }
                }
            }
 
            // addition and deletion of bonds as required
            ArrayList<int[]> zadd = tmol.getBondAdd();
            ArrayList<int[]> zdel = tmol.getBondDel();
 
            if (zadd.size() > 0 || zdel.size() > 0)
            {
                fw.write(NL);
                fw.flush();
 
                for (int i = 0; i < zadd.size(); i++)
                {
                    int[] z = zadd.get(i);
                    line = String.format("%6d%6d%n", z[0], z[1]);
                    fw.write(line);
                    fw.flush();
                }
 
                if (zdel.size() > 0)
                {
                    fw.write(NL);
                    fw.flush();
 
                    for (int i = 0; i < zdel.size(); i++)
                    {
                        int[] z = zdel.get(i);
                        line = String.format("%6d%6d%n", z[0], z[1]);
                        fw.write(line);
                        fw.flush();
                    }
                }
            }
 
        }
        catch (IOException ioe)
        {
            throw new DENOPTIMException(ioe);
        }
        finally
        {
            try
            {
                if (fw != null)
                {
                    fw.close();
                }
            }
            catch (IOException ioe)
            {
                throw new DENOPTIMException(ioe);
            }
        }
    }
 
//------------------------------------------------------------------------------
 
    public static ArrayList<Double> readPSSROTOutput(String filename)
            throws DENOPTIMException
    {
        ArrayList<Double> energies = new ArrayList<>();
 
        BufferedReader br = null;
        String line;
        try
        {
            br = new BufferedReader(new FileReader(filename));
 
            while ((line = br.readLine()) != null)
            {
                line = line.trim();
                if (line.contains("Final Function Value and Deformation"))
                {
                    String str = line.substring(38);
                    // read the number of atoms and store
                    String[] arr = str.split("\\s+");
                    energies.add(Double.parseDouble(arr[1]));
                }
            }
        }
        catch (NumberFormatException | IOException nfe)
        {
            throw new DENOPTIMException(nfe);
        }
        finally
        {
            try
            {
                if (br != null)
                {
                    br.close();
                }
            }
            catch (IOException ioe)
            {
                throw new DENOPTIMException(ioe);
            }
        }
 
        if (energies.isEmpty())
        {
            String msg = "No data found in file: " + filename;
            throw new DENOPTIMException(msg);
        }
        return energies;
    }
 
//------------------------------------------------------------------------------
 
    public static void readPSSROTParams(String filename, 
            ArrayList<String> initPars, ArrayList<String> restPars)
                    throws DENOPTIMException
    {
        BufferedReader br = null;
        String line;
        int fnd = -1;
 
        try
        {
            br = new BufferedReader(new FileReader(filename));
            while ((line = br.readLine()) != null)
            {
                if (line.trim().length() == 0)
                    continue;
                if (line.contains("INIT"))
                {
                    fnd = 0;
                    continue;
                }
                if (line.contains("REST"))
                {
                    fnd = 1;
                    continue;
                }
                if (fnd == 0)
                    initPars.add(line.trim());
                if (fnd == 1)
                    restPars.add(line.trim());
            }
        }
        catch (IOException nfe)
        {
            String msg = "File '" + filename + "' not found.";
            throw new DENOPTIMException(msg, nfe);
        }
        finally
        {
            try
            {
                if (br != null)
                {
                    br.close();
                }
            }
            catch (IOException ioe)
            {
                throw new DENOPTIMException(ioe);
            }
        }
 
        if (initPars.isEmpty())
        {
            String msg = "No data found in file: " + filename;
            throw new DENOPTIMException(msg);
        }
        if (restPars.isEmpty())
        {
            String msg = "No data found in file: " + filename;
            throw new DENOPTIMException(msg);
        }
    }
 
//------------------------------------------------------------------------------
 
    public static HashMap<String, Integer> readTinkerAtomTypes(String filename)
            throws DENOPTIMException
    {
        HashMap<String, Integer> atomTypesMap = new HashMap<>();
 
        BufferedReader br = null;
        String line;
 
        try
        {
            br = new BufferedReader(new FileReader(filename));
            while ((line = br.readLine()) != null)
            {
                //Read only lines starting with keyword "atom"
                line = line.trim();
                if (!line.startsWith("atom"))
                {
                    continue;
                }
 
                //Format:
                //key, class, symbol, "label", Z, atomic weight, connectivity
                try
                {
                    //extract atom type (or 'class' according to Tinker's 
                    // nomenclature) and atom symbol
                    String[] dq = line.split("\"");
                    String fp = dq[0];
                    String[] str1 = fp.split("\\s+");
 
                    //Check the format by reading all parts of atom type def.
                    int atomType = Integer.parseInt(str1[1]);
                    String symbol = str1[2];
                    /* Not needed so far...
                    String sp = dq[2];
                    String[] str2 = sp.split("\\s+");
                    String label = dq[1];
                    int z = Integer.parseInt(str2[1]);
                    double atmWeight = Double.parseDouble(str2[2]);
                    int cn = Integer.parseInt(str2[3]);
                    */
 
                    //Store
                    atomTypesMap.put(symbol,atomType);                    
                }
                catch (Throwable t)
                {
                    String msg = "Format of Tinker's atom type definition not "
                            + "recognized. " + NL + "Details: " + NL 
                            + t.getMessage();
                    throw new DENOPTIMException(msg);
                }
            }
        }
        catch (NumberFormatException | IOException nfe)
        {
            throw new DENOPTIMException(nfe);
        }
        finally
        {
            try
            {
                if (br != null)
                {
                    br.close();
                }
            }
            catch (IOException ioe)
            {
                throw new DENOPTIMException(ioe);
            }
        }
 
        if (atomTypesMap.isEmpty())
        {
            String msg = "No data found in file: " + filename;
            throw new DENOPTIMException(msg);
        }
 
        return atomTypesMap;
    }
 
//------------------------------------------------------------------------------
 
    public static TinkerMolecule getICFromIAC(IAtomContainer mol, 
                       HashMap<String,Integer> tMap )  throws DENOPTIMException
    {
        TinkerMolecule tm = new TinkerMolecule();
        String doneBnd = "visitedBond";
        for (int i=0; i<mol.getAtomCount(); i++)
        {
            int i2 = 0;
            int i3 = 0;
            int i4 = 0;
            int i5 = 0;
            double d = 0.0;
            double a = 0.0;
            double t = 0.0;
 
            int[] nbrs = new int[] {0, 0, 0, 0};
 
            IAtom atmI = mol.getAtom(i);
            if (debug)
            {
                System.err.println("Atom to IC: "
                        +MoleculeUtils.getSymbolOrLabel(atmI)+" "+i);
            }
 
            // define the bond length
            if (i>0)
            {
                i2 = getFirstRefAtomId(i,mol);
                d = atmI.getPoint3d().distance(mol.getAtom(i2).getPoint3d());
                mol.getBond(atmI,mol.getAtom(i2)).setProperty(doneBnd,"T");
                if (debug)
                {
                    System.err.println(" i2 = " + i2 + " d: " + d);
                }
                nbrs[0] = i2+1;
            }
 
            // define the bond angle
            if (i>1)
            {
                i3 = getSecondRefAtomId(i,i2,mol);
                a = MathUtils.angle(atmI.getPoint3d(),
                                            mol.getAtom(i2).getPoint3d(),
                                            mol.getAtom(i3).getPoint3d());
                if (debug)
                {
                    System.err.println(" i3 = "+ i3 + " a: " + a);
                }
                nbrs[0] = i2+1;
                nbrs[1] = i3+1;
            }
 
            // decide on dihedral or second angle
            if (i>2)
            {
                ObjectPair op = getThirdRefAtomId(i,i2,i3,mol,tm);
                i4 = (int) op.getFirst();
                i5 = (int) op.getSecond();
                if (i5==1)
                {
                    t = MathUtils.angle(atmI.getPoint3d(),
                                                mol.getAtom(i2).getPoint3d(),
                                                mol.getAtom(i4).getPoint3d());
                    double sign = MathUtils.torsion(
                                                atmI.getPoint3d(),
                                                mol.getAtom(i2).getPoint3d(),
                                                mol.getAtom(i3).getPoint3d(),
                                                mol.getAtom(i4).getPoint3d());
                    if (sign > 0.0)
                    {
                        i5 = -1;
                    } 
                }
                else
                {
                    IAtom atmJ = mol.getAtom(i2);
                    IAtom atmK = mol.getAtom(i3);
                    IAtom atmL = mol.getAtom(i4);
                    double valueIJKL = MathUtils.torsion(
                            atmI.getPoint3d(), atmJ.getPoint3d(),
                            atmK.getPoint3d(), atmL.getPoint3d());
                    
                    t = valueIJKL;
                    
                    IBond bnd = mol.getBond(atmJ, atmK);
                    if (bnd!=null)
                    {
                        Object cnstrDefObj = bnd.getProperty(
                                RotationalSpaceUtils.PROPERTY_ROTDBDCSTR_DEF);
                        if (cnstrDefObj!=null)
                        {
                            // For clarity, IJKL are the atoms identifying the
                            // dihedral used in the internal coordinates,
                            // while ABCD are the atoms used to define the 
                            // constrained dihedral angle.
                            IAtom[] atomsABCD = (IAtom[]) cnstrDefObj;
                            double cnstrABCD = (double) bnd.getProperty(
                                    RotationalSpaceUtils.PROPERTY_ROTDBDCSTR_VALUE);
                            if (atmJ==atomsABCD[2] && atmK==atomsABCD[1])
                            {
                                // We ensure consistent orfer in the definitions
                                // of the two dihedral angles
                                atomsABCD = new IAtom[]{
                                        atomsABCD[3], atomsABCD[2],
                                        atomsABCD[1], atomsABCD[0]};
                                if (cnstrABCD>0)
                                {
                                    cnstrABCD = 360.0 - cnstrABCD;
                                } else {
                                    cnstrABCD = cnstrABCD - 360.0;
                                }
                            }
 
                            double valueABCD = MathUtils.torsion(
                                    atomsABCD[0].getPoint3d(), 
                                    atomsABCD[1].getPoint3d(),
                                    atomsABCD[2].getPoint3d(), 
                                    atomsABCD[3].getPoint3d());
                            double correctionABCD = cnstrABCD - valueABCD;
                            
                            t = valueIJKL + correctionABCD;
                            
                            if (t > 360)
                            {
                                t = t - 360;
                            } else if (t < -360) {
                                t = t + 360;
                            }
 
                            if (debug)
                            {
                                System.err.println(" dihedral constrain along "
                                        + i2 + "-" + i3 + ": " 
                                        + valueIJKL + " -> " + t 
                                        + " (changed by " + correctionABCD + ")");
                            }
                        }
                    }
                }
                if (debug)
                {
                    System.err.println(" i4 = "+ i4 + " t: " + t + " " + i5);
                }
                nbrs[0] = i2+1;
                nbrs[1] = i3+1;
                nbrs[2] = i4+1;
                nbrs[3] = i5;
            }
 
            String symb = MoleculeUtils.getSymbolOrLabel(atmI);
            int atyp = 0; // Tinker types are assigned later
 
            TinkerAtom ta = new TinkerAtom(i+1, symb, atyp,
                                           new double[] {atmI.getPoint3d().x,
                                                         atmI.getPoint3d().y, 
                                                         atmI.getPoint3d().z},
                                                                         nbrs,
                                                        new double[] {d,a,t});
 
            long vidx = atmI.getProperty(DENOPTIMConstants.ATMPROPVERTEXID);
            ta.setVertexId(vidx);
 
            if (debug)
            {
                System.err.println(" TinkerAtom: "+ta.toString());
            }
 
            tm.addAtom(ta);
 
            if (debug)
            {
                System.err.println("TinkerMolecule: ");
                tm.printIC();
            }
        }
 
        // Add bonds not visited
        for (IBond b : mol.bonds())
        {
            if (b.getProperty(doneBnd) == null)
            {
                tm.addBond(mol.indexOf(b.getAtom(0))+1,
                           mol.indexOf(b.getAtom(1))+1);
            }
        }
 
        // Due to the assumption that all atoms are part of the same 
        // connected network, no bond has to be deleted
 
        // Fix TinkerAtom types
        setTinkerTypes(tm,tMap);
 
        return tm;
    }
 
//----------------------------------------------------------------------------
   
    private static int getFirstRefAtomId(int i1, IAtomContainer mol)
    {
        List<ConnectedLigand> candidates = new ArrayList<ConnectedLigand>();
        for (IAtom nbr : mol.getConnectedAtomsList(mol.getAtom(i1)))
        {
            if (mol.indexOf(nbr) < i1)
            {
                ConnectedLigand cl = new ConnectedLigand(nbr,1);
                candidates.add(cl);
            }
        }
        Collections.sort(candidates, new ConnectedLigandComparator());
        int i2 = mol.indexOf(candidates.get(0).getAtom());
        return i2;
    }
 
//----------------------------------------------------------------------------
 
    private static int getSecondRefAtomId(int i1, int i2, IAtomContainer mol)
    {
        List<ConnectedLigand> candidates = new ArrayList<ConnectedLigand>();
        for (IAtom nbr : mol.getConnectedAtomsList(mol.getAtom(i2)))
        {
            if ((mol.indexOf(nbr) < i1) && (nbr != mol.getAtom(i1)))
            {
                ConnectedLigand cl = new ConnectedLigand(nbr,1);
                candidates.add(cl);
            }
        }
        Collections.sort(candidates, new ConnectedLigandComparator());
        int i3 = mol.indexOf(candidates.get(0).getAtom());
        return i3;
    }
 
//----------------------------------------------------------------------------
 
    private static ObjectPair getThirdRefAtomId(int i1, int i2, int i3, 
               IAtomContainer mol, TinkerMolecule tm) throws DENOPTIMException
    {
        int i5 = 0;
        IAtom atmI1 = mol.getAtom(i1);
        IAtom atmI2 = mol.getAtom(i2);
        IAtom atmI3 = mol.getAtom(i3);
        List<ConnectedLigand> candidates = new ArrayList<ConnectedLigand>();
        if (tm.isTorsionUsed(i2+1, i3+1) || 
            countPredefinedNeighbours(i1,atmI3,mol)==1)
        {
            i5 = 1;
            for (IAtom nbr : mol.getConnectedAtomsList(atmI2))
            {
                if (debug)
                {
                   System.err.println("  Eval. 3rd (ANG): " + 
                           MoleculeUtils.getSymbolOrLabel(nbr)
                   + mol.indexOf(nbr) + " " 
                   + (mol.indexOf(nbr) < i1) + " "
                   + (nbr != atmI1) + " " 
                   + (nbr != atmI3));
                }
                if ((mol.indexOf(nbr) < i1) && (nbr != atmI1) && 
                                                        (nbr != atmI3))
                {
                    double dbcAng = MathUtils.angle(nbr.getPoint3d(),
                                                          atmI2.getPoint3d(),
                                                          atmI3.getPoint3d());
                    if(dbcAng > 1.0)
                    {
                        ConnectedLigand cl = new ConnectedLigand(nbr,1);
                        candidates.add(cl);
                    }
                    else
                    {
                        if (debug)
                        {
                            System.err.println("  ...but collinear with "
                                    + MoleculeUtils.getSymbolOrLabel(atmI3) + i3
                                               + " (i4-i2-i3: " + dbcAng 
                                               + ")");
                        }
                    }
                }
            }
        }
        else
        {
            i5 = 0;
            for (IAtom nbr : mol.getConnectedAtomsList(atmI3))
            {
                if (debug)
                {
                   System.err.println("  Eval. 3rd (TOR): " 
                            + MoleculeUtils.getSymbolOrLabel(nbr)
                   + mol.indexOf(nbr) + " "
                   + (mol.indexOf(nbr) < i1) + " "
                   + (nbr != atmI1) + " "
                   + (nbr != atmI2));
                }
                if ((mol.indexOf(nbr) < i1) && (nbr != atmI1) &&
                                                        (nbr != atmI2))
                {
                    ConnectedLigand cl = new ConnectedLigand(nbr,1);
                    candidates.add(cl);
                }
            }
        }
        Collections.sort(candidates, new ConnectedLigandComparator());
        if (candidates.size() == 0)
        {
            String msg = "Unable to make internal coordinates. Please, "
                         + "consider the use of dummy atoms in proximity "
                         + "of atom " + tm.getAtom(i1+1);
            throw new DENOPTIMException(msg);
        }
        int i4 = mol.indexOf(candidates.get(0).getAtom());
 
        ObjectPair op = new ObjectPair(i4,i5);
        
        return op;
    }
 
//------------------------------------------------------------------------------
 
    private static int countPredefinedNeighbours(int i, IAtom a, IAtomContainer mol)
    {
        int tot = 0;
        for (IAtom nbr : mol.getConnectedAtomsList(a))
        {
            if (mol.indexOf(nbr) < i)
                tot++;
        }
        return tot;
    }
 
//------------------------------------------------------------------------------
 
    public static void setTinkerTypes(TinkerMolecule tmol, 
                        HashMap<String,Integer> tMap) throws DENOPTIMException
    {
        ArrayList<TinkerAtom> lstAtoms = tmol.getAtoms();
        int numberOfAtoms = lstAtoms.size();
        for (int i = 0; i < numberOfAtoms; i++)
        {
            TinkerAtom tatom = lstAtoms.get(i);
 
            String st = tatom.getAtomString().trim();
            if (!tMap.containsKey(st))
            {
                String msg = "Unable to assign atom type to atom '" + st +"'. ";
                if (st.equals("R"))
                {
                    msg = msg + "Unusual dummy atom symbols get atom symbol "
                      + "'R'. Please, add atom type 'R' in your atom type map.";
                }
                throw new DENOPTIMException(msg);
            }
            Integer val = tMap.get(st);
            if (val != null)
            {
                tatom.setAtomType(val.intValue());
                if (debug)
                    System.err.println("Set parameter for " + st + " " + val);
            }
            else
            {
                String msg = "No valid Tinker atom type assigned for atom "
                            + st + ".\n";
                throw new DENOPTIMException(msg);
            }
        }
    }
    
//------------------------------------------------------------------------------
    
    public static void ensureOutputExistsOrRelayError(String outputPathName, 
            String logPathName, String taskName) throws TinkerException
    {
        File output = new File(outputPathName);
        if (output.exists() && output.canRead())
        {
            return;
        }
        
        String errMsg = "TINKER ERROR: ";
        ReversedLinesFileReader fr = null;
        try
        {
            fr = new ReversedLinesFileReader(new File(
                    logPathName), null);
 
            int numEmpty = 0;
            for (int i=0; i<100; i++) //at most 100 lines are read
            {
                String line = fr.readLine();
                if (line==null)
                    break;
                
                if (line.trim().isEmpty())
                    numEmpty++;
                
                if (numEmpty==2)
                    break;
                errMsg = errMsg + NL + line;
            }
            fr.close();
        } catch (IOException e)
        {
            throw new TinkerException("Missing Tinker log '" + logPathName + "'",
                    taskName);
        }
        
        throw new TinkerException(errMsg, taskName);
    }
 
//------------------------------------------------------------------------------
    
    public static String getNameLastCycleFile(String workDir, String fname, 
            String tinkerLog, String pattern) throws DENOPTIMException
    {
        int lastI = MMBuilderUtils.countLinesWKeywordInFile(tinkerLog, pattern);
        String xyzfile = workDir + System.getProperty("file.separator") + fname 
                + "." + GeneralUtils.getPaddedString(3, lastI - 1);
        return xyzfile;
    }
 
    
//------------------------------------------------------------------------------
 
}