/source/RMG/jing/chem/QMTP.java

Large files files are truncated, but you can click here to view the full file

// //////////////////////////////////////////////////////////////////////////////
//
// RMG - Reaction Mechanism Generator
//
// Copyright (c) 2002-2011 Prof. William H. Green (whgreen@mit.edu) and the
// RMG Team (rmg_dev@mit.edu)
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// //////////////////////////////////////////////////////////////////////////////
package jing.chem;

import java.util.*;
import jing.chemParser.ChemParser;
import jing.chemUtil.*;
import jing.param.*;
import java.io.File;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.io.BufferedReader;
import java.io.InputStream;
import java.io.InputStreamReader;
import jing.rxnSys.Logger;

// quantum mechanics thermo property estimator; analog of GATP
public class QMTP implements GeneralGAPP {
    final protected static double ENTHALPY_HYDROGEN = 52.1; // needed for HBI
    private static QMTP INSTANCE = new QMTP(); // ## attribute INSTANCE
    protected static PrimaryThermoLibrary primaryLibrary;// Note: may be able to separate this out into GeneralGAPP, as
// this is common to both GATP and QMTP
    public static String qmfolder = System.getProperty("RMG.qmCalculationsDir"); // location for Gaussian/mopac output
    protected LinkedHashMap<String, ThermoData> qmLibrary = new LinkedHashMap<String, ThermoData>();
    // protected static LinkedHashMap library; //as above, may be able to move this and associated functions to GeneralGAPP
// (and possibly change from "x implements y" to "x extends y"), as it is common to both GATP and QMTP
    protected Vector<String> failedQm = new Vector<String>();
    protected ThermoGAGroupLibrary thermoLibrary; // needed for HBI
    public static String qmprogram = "both";// the qmprogram can be "mopac", "gaussian03", "both" (MOPAC and Gaussian),
// or "mm4"/"mm4hr"
    public static boolean usePolar = false; // use polar keyword in MOPAC
    public static boolean useCanTherm = true; // whether to use CanTherm in MM4 cases for interpreting output via
// force-constant matrix; this will hopefully avoid zero frequency issues
    public static boolean useHindRot = false;// whether to use HinderedRotor scans with MM4 (requires useCanTherm=true)
    public static boolean keepQMfiles = true; // whether to keep all the Qmfiles generated by Gaussian/mopac or not
    public static int connectivityCheck = 0; // level of connectivity checking; 0: "off": no checking; 1: "check":
// (print warning if connectivity doesn't seem to match); 2: "confirm": consider the run a failure if the connectivity
// doesn't seem to match
    public static double R = 1.9872; // ideal gas constant in cal/mol-K (does this appear elsewhere in RMG, so I don't
// need to reuse it?)
    public static double Hartree_to_kcal = 627.5095; // conversion from Hartree to kcal/mol taken from Gaussian thermo
// white paper
    public static double Na = 6.02214179E23;// Avagadro's number; cf.
// http://physics.nist.gov/cgi-bin/cuu/Value?na|search_for=physchem_in!
    public static double k = 1.3806504E-23;// Boltzmann's constant in J/K; cf.
// http://physics.nist.gov/cgi-bin/cuu/Value?na|search_for=physchem_in!
    public static double h = 6.62606896E-34;// Planck's constant in J-s; cf.
// http://physics.nist.gov/cgi-bin/cuu/Value?h|search_for=universal_in!
    public static double c = 299792458. * 100;// speed of light in vacuum in cm/s, cf.
// http://physics.nist.gov/cgi-bin/cuu/Value?c|search_for=universal_in!
    public static double deltaTheta = 5.0;// degree increment for rotor scans when using useHindRot

    // Constructors
    // ## operation QMTP()
    private QMTP() {
        initializePrimaryThermoLibrary();
        initalizeQmLibrary();
    }

    public String getQmMethod()
	throws CouldNotDetermineQMMethodException {
        String qmProgram = qmprogram;
        String qmMethod = "";
        if (qmProgram.equals("mm4") || qmProgram.equals("mm4hr")) {
            qmMethod = "mm4";
            if (qmProgram.equals("mm4hr"))
                useHindRot = true;
        } else if (qmProgram.equals("mopac")) {
            qmMethod = "pm7"; // may eventually want to pass this to various functions to choose which "sub-function" to
// call
        } else if (qmProgram.equals("gaussian03") || qmProgram.equals("both")) {
	    qmMethod = "pm3";
	} else {
	    throw new CouldNotDetermineQMMethodException();
        }
        return qmMethod;
    }

    public ThermoData generateFakeThermoData() {
        ThermoData result = null;
        ThermoGAValue impossible = new ThermoGAValue(1000
                , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, null);
        result.plus(impossible);
        return result;
    }

    // ## operation generateThermoData(ChemGraph)
    public ThermoData generateThermoData(ChemGraph p_chemGraph) {
        if (!keepQMfiles) {
            cleanQmFiles();
        }
        // #[ operation generateThermoData(ChemGraph)
        // first, check for thermo data in the primary thermo library and library (?); if it is there, use it
        ThermoData result = null;
        ThermoData tmpTherm = primaryLibrary.getThermoData(p_chemGraph
                .getGraph());
        // Logger.info(result);
        if (tmpTherm != null) {
            result = tmpTherm.copyWithExtraInfo();// use a copy of the object!; that way, subsequent modifications of
// this object don't change the primary thermo library
            p_chemGraph.fromprimarythermolibrary = true;
            return result;
        }
        result = new ThermoData();
        int maxRadNumForQM = Global.maxRadNumForQM;
        if (p_chemGraph.getRadicalNumber() > maxRadNumForQM)// use HBI if the molecule has more radicals than
// maxRadNumForQM; this is helpful because ; also MM4 (and MM3) look like they may have issues with radicals
        {// this code is based closely off of GATP saturation (in getGAGroup()), but there are some modifications,
// particularly for symmetry correction
         // find the initial symmetry number
            int sigmaRadical = p_chemGraph.getSymmetryNumber();
            Graph g = p_chemGraph.getGraph();
            LinkedHashMap oldCentralNode = (LinkedHashMap) (p_chemGraph.getCentralNode())
                    .clone();
            // saturate radical site
            int max_radNum_molecule = ChemGraph.getMAX_RADICAL_NUM();
            int max_radNum_atom = Math.min(8, max_radNum_molecule);
            int[] idArray = new int[max_radNum_molecule];
            Atom[] atomArray = new Atom[max_radNum_molecule];
            Node[][] newnode = new Node[max_radNum_molecule][max_radNum_atom];
            int radicalSite = 0;
            Iterator iter = p_chemGraph.getNodeList();
            FreeElectron satuated = FreeElectron.make("0");
            while (iter.hasNext()) {
                Node node = (Node) iter.next();
                Atom atom = (Atom) node.getElement();
                if (atom.isRadical()) {
                    radicalSite++;
                    // save the old radical atom
                    idArray[radicalSite - 1] = node.getID().intValue();
                    atomArray[radicalSite - 1] = atom;
                    // new a satuated atom and replace the old one
                    Atom newAtom = Atom.make(atom.getChemElement(), satuated);
                    node.setElement(newAtom);
                    node.updateFeElement();
                }
            }
            // add H to saturate chem graph
            Atom H = Atom.make(ChemElement.make("H"), satuated);
            Bond S = Bond.make("S");
            for (int i = 0; i < radicalSite; i++) {
                Node node = p_chemGraph.getNodeAt(idArray[i]);
                Atom atom = atomArray[i];
                int HNum = atom.getRadicalNumber();
                for (int j = 0; j < HNum; j++) {
                    newnode[i][j] = g.addNode(H);
                    g.addArcBetween(node, S, newnode[i][j]);
                }
                node.updateFgElement();
            }
            // find the saturated symmetry number
            int sigmaSaturated = p_chemGraph.getSymmetryNumber();
            // result = generateThermoData(g);//I'm not sure what GATP does, but this recursive calling will use HBIs on
// saturated species if it exists in PrimaryThermoLibrary
            // check the primary thermo library for the saturated graph
            tmpTherm = primaryLibrary.getThermoData(p_chemGraph.getGraph());
            // Logger.info(result);
            if (tmpTherm != null) {
                result = tmpTherm.copyWithExtraInfo();// use a copy of the object!; that way, subsequent modifications
// of this object don't change the primary thermo library
                p_chemGraph.fromprimarythermolibrary = false;// we don't want to set fromprimarythermolibrary to true,
// because the result is not directly from the PTL, but comes via PTL + HBI corrections; if true is set here, this would
// affect two things, both in Species.java: 1) the naming; in the HBI case, we don't want to use a name derived from the
// thermo name as weird things can happen 2)findStablestThermoData considers the value to be the final word; however,
// since this is a radical that is not directly in the primaryThermoLibrary, we want to consider alternative thermo for
// all possible resonance isomers
            } else {
                tmpTherm = getQMThermoData(p_chemGraph);
                result = tmpTherm.copyWithExtraInfo();
            }
            // find the BDE for all radical groups
            if (thermoLibrary == null)
                initGAGroupLibrary();
            for (int i = 0; i < radicalSite; i++) {
                int id = idArray[i];
                Node node = g.getNodeAt(id);
                Atom old = (Atom) node.getElement();
                node.setElement(atomArray[i]);
                node.updateFeElement();
                // get rid of the extra H at ith site
                int HNum = atomArray[i].getRadicalNumber();
                for (int j = 0; j < HNum; j++) {
                    g.removeNode(newnode[i][j]);
                }
                node.updateFgElement();
                p_chemGraph.resetThermoSite(node);
                ThermoGAValue thisGAValue = thermoLibrary
                        .findRadicalGroup(p_chemGraph);
                if (thisGAValue == null) {
                    Logger.error("Radical group not found: " + node.getID());
                } else {
                    // Logger.info(node.getID() + " radical correction: " + thisGAValue.getName() +
// "  "+thisGAValue.toString());
                    result.plus(thisGAValue);
                }
                // recover the saturated site for next radical site calculation
                node.setElement(old);
                node.updateFeElement();
                for (int j = 0; j < HNum; j++) {
                    newnode[i][j] = g.addNode(H);
                    g.addArcBetween(node, S, newnode[i][j]);
                }
                node.updateFgElement();
            }
            // recover the chem graph structure
            // recover the radical
            for (int i = 0; i < radicalSite; i++) {
                int id = idArray[i];
                Node node = g.getNodeAt(id);
                node.setElement(atomArray[i]);
                node.updateFeElement();
                int HNum = atomArray[i].getRadicalNumber();
                // get rid of extra H
                for (int j = 0; j < HNum; j++) {
                    g.removeNode(newnode[i][j]);
                }
                node.updateFgElement();
            }
            // subtract the enthalphy of H from the result
            int rad_number = p_chemGraph.getRadicalNumber();
            ThermoGAValue enthalpy_H = new ThermoGAValue(ENTHALPY_HYDROGEN
                    * rad_number, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, null);
            result.minus(enthalpy_H);
            // correct the symmetry number based on the relative radical and saturated symmetry number; this should
// hopefully sidestep potential complications based on the fact that certain symmetry effects could be included in HBI
// value itself, and the fact that the symmetry number correction for saturated molecule has already been implemented,
// and it is likely to be different than symmetry number considered here, since the correction for the saturated
// molecule will have been external symmetry number, whereas RMG's ChemGraph symmetry number estimator includes both
// internal and external symmetry contributions; even so, I don't know if this will handle a change from chiral to
// achiral (or vice versa) properly
            ThermoGAValue symmetryNumberCorrection = new ThermoGAValue(0, -1
                    * GasConstant.getCalMolK()
                    * Math.log((double) (sigmaRadical)
                            / (double) (sigmaSaturated)), 0, 0, 0, 0, 0, 0, 0,
                    0, 0, 0, null);
            result.plus(symmetryNumberCorrection);
            p_chemGraph.setCentralNode(oldCentralNode);
            // display corrected thermo to user
            String[] InChInames = getQMFileName(p_chemGraph);// determine the filename (InChIKey) and InChI with
// appended info for triplets, etc.
            String name = InChInames[0];
            String InChIaug = InChInames[1];
            Logger.info("HBI-based thermo for " + name + "(" + InChIaug + "): "
                    + result.toString());// print result, at least for debugging purposes
        } else {
            // no need for HBI radical corrections, just get the QM result
            tmpTherm = getQMThermoData(p_chemGraph);
            result = tmpTherm.copyWithExtraInfo();
        }
        return result;
        // #]
    }

    public ThermoData getQMThermoData(ChemGraph p_chemGraph) {
        // Try to get the QMThermoData from the qmLibrary, if it's not there then call getQMThermoData and save it for
// future.
        ThermoData result = null;
        // First to get from qmLibrary
        String[] InChInames = getQMFileName(p_chemGraph);// determine the filename (InChIKey) and InChI with appended
// info for triplets, etc.
        String name = InChInames[0];
        String InChIaug = InChInames[1];
        ThermoData tempTherm = qmLibrary.get(InChIaug);
        if (tempTherm != null) {
            result = tempTherm.copyWithExtraInfo(); // use a copy of the object!; that way, subsequent modifications of
// this object don't change the QM library
            Logger.info("QM calculation for "
                    + InChIaug
                    + " previously performed. Using Thermo from previous results in QMlibrary.");
            result.setSource(result.comments);
            // Added for debugging
            Logger.info("Thermo Data for " + InChIaug + " is " + result);
        } else { // couldn't find it in qmLibrary
                 // Check to see if QMTP has failed all attempts on this molecule before
                 // If so, generate Thermo using Benson groups
            if (failedQm.contains(InChIaug)) {
                Logger.info("All attempts at QMTP for "
                        + name
                        + "("
                        + InChIaug
                        + ") have previously failed. Falling back to Benson Group Additivity.");
                TDGenerator gen = new BensonTDGenerator();
                return gen.generateFakeThermo();
            //    return gen.generateThermo(p_chemGraph);
            }
            // generate new QM Thermo Data
            try {
                String qmMethod = getQmMethod();
                result = generateQMThermoData(p_chemGraph, qmMethod);
                // now save it for next time
                QMLibraryEditor.addQMTPThermo(p_chemGraph, InChIaug, result,
                        qmMethod, qmprogram);
                qmLibrary.put(InChIaug, result);
            } catch (AllQmtpAttemptsFailedException e) {
                Logger.warning("Falling back to Benson group additivity due to repeated failure in QMTP calculations");
                TDGenerator gen = new BensonTDGenerator();
             //   return gen.generateThermo(p_chemGraph);
		return gen.generateFakeThermo();
            }
        }
        return result;
    }

    private Thread setHoldAndGetClearerThread(final String name) {
        /*
         * Creates a .hold file for the given name, and creates a shutdown hook to delete the hold file. Returns the
         * hook, which is a thread which you can later run and de-register.
         */
        Logger.debug("Creating hold file for " + name);
        final File holdFile = new File(qmfolder, name + ".hold");
        try {
            holdFile.createNewFile();
        } catch (Exception e) {
            Logger.error("Error creating .hold file for " + name + ": "
                    + e.toString());
            System.exit(0);
            throw new RuntimeException("Error creating .hold file for " + name,
                    e);
        }
        Thread hook = new Thread() {
            public void run() {
                Logger.debug("Deleting hold file for " + name);
                holdFile.delete();
            };
        };
        Runtime.getRuntime().addShutdownHook(hook);
        return hook;
    }

    private boolean checkHoldExists(String name) {
        /*
         * Checks for existance of a .hold file for the given name. Returns true if it exists (and you shouldn't run a
         * new job).
         */
        Logger.debug("Checking hold file for " + name);
        final File holdFile = new File(qmfolder, name + ".hold");
        return holdFile.exists();
    }

    private void clearHold(Thread holdClearerThread) {
        /*
         * Clears the hold, by running and unregistering the hold-clearing thread.
         */
        Runtime.getRuntime().removeShutdownHook(holdClearerThread);
        holdClearerThread.start();
        try {
            holdClearerThread.join(); // wait for it to finish
        } catch (InterruptedException e) {
            Thread.currentThread().interrupt();
        }
    }

    public ThermoData generateQMThermoData(ChemGraph p_chemGraph, String qmMethod)
            throws AllQmtpAttemptsFailedException {
        // if there is no data in the libraries, calculate the result based on QM or MM calculations; the below steps
// will be generalized later to allow for other quantum mechanics packages, etc.
        String qmProgram = qmprogram;
        ThermoData result = new ThermoData();
        double[] dihedralMinima = null;
        String[] InChInames = getQMFileName(p_chemGraph);// determine the filename (InChIKey) and InChI with appended
// info for triplets, etc.
        String name = InChInames[0];
        String InChIaug = InChInames[1];
        String directory = qmfolder;
        File dir = new File(directory);
        directory = dir.getAbsolutePath();// this and previous three lines get the absolute path for the directory
        // check for existing hold file before starting calculations (to avoid the possibility of interference with
// other jobs using the same QMfiles folder)
        try {
            while (checkHoldExists(name)) {
                Logger.info("Existence of hold file for "
                        + name
                        + " suggests that another RMG process is currently running calculations on this molecule; waiting for other RMG process to finish; will check again in 60 seconds...");
                Thread.sleep(60000);
            }
        } catch (Exception e) {
            Logger.error("Unexpected error while waiting for hold to be lifted: "
                    + e.toString());
            System.exit(0);
        }
        if (qmMethod.equals("pm3") || qmMethod.equals("pm7")) { 
            // first, check to see if the result already exists and the job terminated successfully
            boolean gaussianResultExists = successfulGaussianResultExistsQ(
                    name, directory, InChIaug);
            boolean mopacResultExists = successfulMopacResultExistsQ(name,
                    directory, InChIaug);
            if (!gaussianResultExists && !mopacResultExists) {// if a successful result doesn't exist from previous run
// (or from this run), run the calculation; if a successful result exists, we will skip directly to parsing the file
                // step 0: create a .hold file to prevent other jobs running with the same directory from interfering by
// trying to run with the same molecule
                Thread holdClearer = setHoldAndGetClearerThread(name);
                // steps 1 and 2: create 2D and 3D mole files
                molFile p_3dfile = create3Dmolfile(name, p_chemGraph);
                // 3. create the Gaussian or MOPAC input file
                directory = qmfolder;
                dir = new File(directory);
                directory = dir.getAbsolutePath();// this and previous three lines get the absolute path for the
// directory
                int attemptNumber = 1;// counter for attempts using different keywords
                int successFlag = 0;// flag for success of Gaussian run; 0 means it failed, 1 means it succeeded
                int maxAttemptNumber = 1;
                int multiplicity = p_chemGraph.getRadicalNumber() + 1; // multiplicity = radical number + 1
                while (successFlag == 0 && attemptNumber <= maxAttemptNumber) {
                    // IF block to check which program to use
                    if (qmProgram.equals("gaussian03")) {
                        if (p_chemGraph.getAtomNumber() > 1) {
                            maxAttemptNumber = createGaussianInput(name,
                                    directory, p_3dfile, attemptNumber,
                                    InChIaug, multiplicity, qmMethod);
                        } else {
                            maxAttemptNumber = createGaussianInput(name,
                                    directory, p_3dfile, -1, InChIaug,
                                    multiplicity, qmMethod);// use -1 for attemptNumber for monoatomic case
                        }
                        // 4. run Gaussian
                        successFlag = runGaussian(name, directory, InChIaug);
                    } else if (qmProgram.equals("mopac")
                            || qmProgram.equals("both")) {
			// Modified by Aaron
                        maxAttemptNumber = createMopacInput(name, directory,
                                p_3dfile, attemptNumber, InChIaug, multiplicity, qmMethod);
                        successFlag = runMOPAC(name, directory, InChIaug);
                    } else {
                        Logger.critical("Unsupported quantum chemistry program");
                        System.exit(0);
                    }
                    // new IF block to check success
                    if (successFlag == 1) {
                        clearHold(holdClearer);
                        Logger.info("Attempt #" + attemptNumber
                                + " on species " + name + " (" + InChIaug
                                + ") succeeded.");
                    } else if (successFlag == 0) {
                        if (attemptNumber == maxAttemptNumber) {// if this is the last possible attempt, and the
// calculation fails, exit with an error message
                            if (qmProgram.equals("both")) { // if we are running with "both" option and all keywords
// fail, try with Gaussian
                                qmProgram = "gaussian03";
                                Logger.info("*****Final MOPAC attempt (#"
                                        + maxAttemptNumber + ") on species "
                                        + name + " (" + InChIaug
                                        + ") failed. Trying to use Gaussian.");
                                attemptNumber = 0;// this needs to be 0 so that when we increment attemptNumber below,
// it becomes 1 when returning to the beginning of the for loop
                                maxAttemptNumber = 1;
                            } else {
                                Logger.info("*****Final attempt (#"
                                        + maxAttemptNumber + ") on species "
                                        + name + " (" + InChIaug + ") failed.");
                                Logger.info(p_chemGraph.toString());
                                clearHold(holdClearer);
                                // Add to augmented InChI failedQm so that RMG will not try to run this molecule in QMTP
// again
                                failedQm.add(InChIaug);
                                throw new AllQmtpAttemptsFailedException();
                            }
                        }
                        Logger.info("*****Attempt #" + attemptNumber
                                + " on species " + name + " (" + InChIaug
                                + ") failed. Will attempt a new keyword.");
                        attemptNumber++;// try again with new keyword
                    }
                }
            }
            // 5. parse QM output and record as thermo data (function includes symmetry/point group calcs, etc.); if
// both Gaussian and MOPAC results exist, Gaussian result is used
            if (gaussianResultExists
                    || (qmProgram.equals("gaussian03") && !mopacResultExists)) {
                result = parseGaussianPM3(name, directory, p_chemGraph);
            } else if (mopacResultExists || qmProgram.equals("mopac")
                    || qmProgram.equals("both")) {
                result = parseMopac(name, directory, p_chemGraph, qmMethod);
            } else {
                Logger.critical("Unexpected situation in QMTP thermo estimation");
                System.exit(0);
            }
        } else {// mm4 case
                // first, check to see if the result already exists and the job terminated successfully
            boolean mm4ResultExists = successfulMM4ResultExistsQ(name,
                    directory, InChIaug);
            if (!mm4ResultExists) {// if a successful result doesn't exist from previous run (or from this run), run the
// calculation; if a successful result exists, we will skip directly to parsing the file
                // step 0: create a .hold file to prevent other jobs running with the same directory from interfering by
// trying to run with the same molecule
                Thread holdClearer = setHoldAndGetClearerThread(name);
                // steps 1 and 2: create 2D and 3D mole files
                molFile p_3dfile = create3Dmolfile(name, p_chemGraph);
                // 3. create the MM4 input file
                directory = qmfolder;
                dir = new File(directory);
                directory = dir.getAbsolutePath();// this and previous three lines get the absolute path for the
// directory
                int attemptNumber = 1;// counter for attempts using different keywords
                int successFlag = 0;// flag for success of MM4 run; 0 means it failed, 1 means it succeeded
                int maxAttemptNumber = 1;
                int multiplicity = p_chemGraph.getRadicalNumber() + 1; // multiplicity = radical number + 1
                while (successFlag == 0 && attemptNumber <= maxAttemptNumber) {
                    maxAttemptNumber = createMM4Input(name, directory,
                            p_3dfile, attemptNumber, InChIaug, multiplicity);
                    // 4. run MM4
                    successFlag = runMM4(name, directory, InChIaug);
                    // new IF block to check success
                    if (successFlag == 1) {
                        clearHold(holdClearer);
                        Logger.info("Attempt #" + attemptNumber
                                + " on species " + name + " (" + InChIaug
                                + ") succeeded.");
                        // run rotor calculations if necessary
                        int rotors = p_chemGraph.getInternalRotor();
                        if (useHindRot && rotors > 0) {
                            // we should re-run scans even if pre-existing scans exist because atom numbering may change
// from case to case; a better solution would be to check for stored CanTherm output and use that if available
                            Logger.info("Running rotor scans on " + name
                                    + "...");
                            dihedralMinima = createMM4RotorInput(name,
                                    directory, p_chemGraph, rotors);// we don't worry about checking InChI here; if
// there is an InChI mismatch it should be caught
                            runMM4Rotor(name, directory, rotors);
                        }
                    } else if (successFlag == 0) {
                        if (attemptNumber == maxAttemptNumber) {// if this is the last possible attempt, and the
// calculation fails, exit with an error message
                            Logger.info("*****Final attempt (#"
                                    + maxAttemptNumber + ") on species " + name
                                    + " (" + InChIaug + ") failed.");
                            Logger.info(p_chemGraph.toString());
                            // delete the hold file
                            clearHold(holdClearer);
                            throw new AllQmtpAttemptsFailedException();
                        }
                        Logger.info("*****Attempt #" + attemptNumber
                                + " on species " + name + " (" + InChIaug
                                + ") failed. Will attempt a new keyword.");
                        attemptNumber++;// try again with new keyword
                    }
                }
                if (useCanTherm) {
                    performCanThermCalcs(name, directory, p_chemGraph,
                            dihedralMinima, false);
                    if (p_chemGraph.getInternalRotor() > 0)
                        performCanThermCalcs(name, directory, p_chemGraph,
                                dihedralMinima, true);// calculate RRHO case for comparison
                }
            }
            // 5. parse MM4 output and record as thermo data (function includes symmetry/point group calcs, etc.)
            if (!useCanTherm)
                result = parseMM4(name, directory, p_chemGraph);
            else {
                // if (qmdata==null) qmdata = getQMDataWithCClib(name, directory, p_chemGraph, true);//get qmdata if it
// is null (i.e. if a pre-existing successful result exists and it wasn't read in above)
                result = parseCanThermFile(name, directory, p_chemGraph);
                if (p_chemGraph.getInternalRotor() > 0)
                    parseCanThermFile(name + "_RRHO", directory, p_chemGraph);// print the RRHO result for comparison
            }
        }
        return result;
    }

    protected static QMTP getINSTANCE() {
        return INSTANCE;
    }

    public void initializePrimaryThermoLibrary() {
        primaryLibrary = PrimaryThermoLibrary.getINSTANCE();
    }

    // Checks to see if QmLibrary already exists. If not, creates the necessary directorys and files
    public void initalizeQmLibrary() {
        try {
            qmLibrary = QMLibraryEditor
                    .readLibrary("QMThermoLibrary/Library.txt");
            QMLibraryEditor.initialize();
        } catch (IOException e) {
            Logger.info("No QM Thermo Library detected. New QM Library being created");
            QMLibraryEditor.initialize();
        }
    }

    // Deletes the contents of the QmFiles directory
    public void cleanQmFiles() {
        File qmFolder = new File(qmfolder);
        // Deletes folder and everything inside
        ChemParser.deleteDir(qmFolder);
        // Remake empty directory
        qmFolder.mkdir();
    }

    // creates a 3D molFile; for monoatomic species, it just returns the 2D molFile
    public molFile create3Dmolfile(String name, ChemGraph p_chemGraph) {
        // 1. create a 2D file
        // use the absolute path for directory, so we can easily reference from other directories in command-line paths
        String directory = System.getProperty("RMG.2DmolfilesDir");
        directory = new File(directory).getAbsolutePath();
        molFile p_2dfile = new molFile(name, directory, p_chemGraph);
        molFile p_3dfile = new molFile();// it seems this must be initialized, so we initialize to empty object
        // 2. convert from 2D to 3D using RDKit if the 2D molfile is for a molecule with 2 or more atoms
        int atoms = p_chemGraph.getAtomNumber();
        if (atoms > 1) {
            int distGeomAttempts = 1;
            if (atoms > 3) {// this check prevents the number of attempts from being negative
                distGeomAttempts = 5 * (p_chemGraph.getAtomNumber() - 3); // number of conformer attempts is just a
// linear scaling with molecule size, due to time considerations; in practice, it is probably more like 3^(n-3) or
// something like that
            }
            p_3dfile = embed3D(p_2dfile, distGeomAttempts);
            return p_3dfile;
        } else {
            return p_2dfile;
        }
    }

    // embed a molecule in 3D, using RDKit
    public molFile embed3D(molFile twoDmolFile, int numConfAttempts) {
        // convert to 3D MOL file using RDKit script
        int flag = 0;
        String directory = System.getProperty("RMG.3DmolfilesDir");
        directory = new File(directory).getAbsolutePath(); // get the absolute path for the directory
        String name = twoDmolFile.getName();
        
        String rdbase = System.getenv("RDBASE");
        if (rdbase == null) {
            Logger.critical("Please set your RDBASE environment variable to the directory containing RDKit.");
            System.exit(0);
        }

        try {
            File runningdir = new File(directory);
            String command = "";
            if (System.getProperty("os.name").toLowerCase().contains("windows")) {// special windows case where paths
// can have spaces and are allowed to be surrounded by quotes
                command = "python \""
                        + System.getProperty("RMG.workingDirectory")
                        + "/scripts/distGeomScriptMolLowestEnergyConf.py\" ";
                String twoDmolpath = twoDmolFile.getPath();
                command = command.concat("\"" + twoDmolpath + "\" ");
                command = command.concat("\"" + name + ".mol\" ");// this is the target file name; use the same name as
// the twoDmolFile (but it will be in he 3Dmolfiles folder
                command = command.concat("\"" + name + ".cmol\" ");// this is the target file name for crude coordinates
// (corresponding to the minimum energy conformation based on UFF refinement); use the same name as the twoDmolFile (but
// it will be in he 3Dmolfiles folder) and have suffix .cmol
                command = command.concat(numConfAttempts + " ");
                command = command.concat("\"" + rdbase + "\"");// pass the $RDBASE environment variable
// to the script so it can use the approprate directory when importing rdkit
            } else {// non-Windows case
                command = "python "
                        + System.getProperty("RMG.workingDirectory")
                        + "/scripts/distGeomScriptMolLowestEnergyConf.py ";
                String twoDmolpath = twoDmolFile.getPath();
                command = command.concat("" + twoDmolpath + " ");
                command = command.concat(name + ".mol ");// this is the target file name; use the same name as the
// twoDmolFile (but it will be in he 3Dmolfiles folder
                command = command.concat(name + ".cmol ");// this is the target file name for crude coordinates
// (corresponding to the minimum energy conformation based on UFF refinement); use the same name as the twoDmolFile (but
// it will be in he 3Dmolfiles folder) and have suffix .cmol
                command = command.concat(numConfAttempts + " ");
                command = command.concat(rdbase);// pass the $RDBASE environment variable to the script
// so it can use the approprate directory when importing rdkit
            }
            Process pythonProc = Runtime.getRuntime().exec(command, null,
                    runningdir);
            String killmsg = "Python process for "
                    + twoDmolFile.getName()
                    + " did not complete within 120 seconds, and the process was killed. File was probably not written.";// message
// to print if the process times out
            // check for errors and display the error if there is one
            InputStream is = pythonProc.getErrorStream();
            InputStreamReader isr = new InputStreamReader(is);
            BufferedReader br = new BufferedReader(isr);
            String line = null;
            while ((line = br.readLine()) != null) {
                line = line.trim();
                Logger.error(line);
                flag = 1;
            }
            // if there was an error, indicate the file and InChI
            if (flag == 1) {
                Logger.info("RDKit received error (see above) on "
                        + twoDmolFile.getName()
                        + ". File was probably not written.");
            }
            int exitValue = pythonProc.waitFor();
            pythonProc.getInputStream().close();
            pythonProc.getOutputStream().close();
            br.close();
            isr.close();
            is.close();
        } catch (Exception e) {
            Logger.logStackTrace(e);
            String err = "Error in running RDKit Python process \n";
            err += e.toString();
            Logger.critical(err);
            System.exit(0);
        }
        // construct molFile pointer to new file (name will be same as 2D mol file
        return new molFile(name, directory);
    }

    // creates Gaussian PM3 input file in directory with filename name.gjf by using OpenBabel to convert p_molfile
    // attemptNumber determines which keywords to try
    // the function returns the maximum number of keywords that can be attempted; this will be the same throughout the
// evaluation of the code, so it may be more appropriate to have this as a "constant" attribute of some sort
    // attemptNumber=-1 will call a special set of keywords for the monoatomic case
    public int createGaussianInput(String name, String directory,
            molFile p_molfile, int attemptNumber, String InChIaug,
            int multiplicity, String qmMethod) {
        // write a file with the input keywords
        int scriptAttempts = 18;// the number of keyword permutations available; update as additional options are added
        int maxAttemptNumber = 2 * scriptAttempts;// we will try a second time with crude coordinates if the UFF refined
// coordinates do not work
        // String inpKeyStr="%chk="+directory+"/RMGrunCHKfile.chk\n";//by not specifying a filename, Gaussian will
// create its own pseudo-random filename, and it will be deleted at the conclusion of a successful calculation; an
// alternative would be to specify a checkpoint filename based on the modified InChIKey... however, it wouldn't be
// automatically deleted after the completion of the calculation
        String inpKeyStr = "%mem=6MW\n";
        inpKeyStr += "%nproc=1\n";
        if (attemptNumber == -1)
            inpKeyStr += "# " + qmMethod + " freq";// keywords for monoatomic case (still need to use freq keyword to get molecular
// mass)
        else if (attemptNumber % scriptAttempts == 1)
            inpKeyStr += "# " + qmMethod + " opt=(verytight,gdiis) freq IOP(2/16=3)";// added IOP option to avoid aborting when
// symmetry changes; 3 is supposed to be default according to documentation, but it seems that 0 (the default) is the
// only option that doesn't work from 0-4; also, it is interesting to note that all 4 options seem to work for test case
// with z-matrix input rather than xyz coords; cf. http://www.ccl.net/cgi-bin/ccl/message-new?2006+10+17+005 for
// original idea for solution
        else if (attemptNumber % scriptAttempts == 2)
            inpKeyStr += "# " + qmMethod + " opt=(verytight,gdiis) freq IOP(2/16=3) IOP(4/21=2)";// use different SCF method; this
// addresses at least one case of failure for a C4H7J species
        else if (attemptNumber % scriptAttempts == 3)
            inpKeyStr += "# " + qmMethod + " opt=(verytight,calcfc,maxcyc=200) freq IOP(2/16=3) nosymm";// try multiple different
// options (no gdiis, use calcfc, nosymm); 7/21/09: added maxcyc option to fix case of MPTBUKVAJYJXDE-UHFFFAOYAPmult3
// (InChI=1/C4H10O5Si/c1-3-7-9-10(5,6)8-4-2/h4-5H,3H2,1-2H3/mult3) (file manually copied to speed things along)
        else if (attemptNumber % scriptAttempts == 4)
            inpKeyStr += "# " + qmMethod + " opt=(verytight,calcfc,maxcyc=200) freq=numerical IOP(2/16=3) nosymm";// 7/8/09:
// numerical frequency keyword version of keyword #3; used to address GYFVJYRUZAKGFA-UHFFFAOYALmult3
// (InChI=1/C6H14O6Si/c1-3-10-13(8,11-4-2)12-6-5-9-7/h6-7H,3-5H2,1-2H3/mult3) case; (none of the existing Gaussian or
// MOPAC combinations worked with it)
        else if (attemptNumber % scriptAttempts == 5)
            inpKeyStr += "# " + qmMethod + " opt=(verytight,gdiis,small) freq IOP(2/16=3)";// 7/10/09: somehow, this worked for
// problematic case of ZGAWAHRALACNPM-UHFFFAOYAF
// (InChI=1/C8H17O5Si/c1-3-11-14(10,12-4-2)13-8-5-7(9)6-8/h7-9H,3-6H2,1-2H3); (was otherwise giving l402 errors); even
// though I had a keyword that worked for this case, I manually copied the fixed log file to QMfiles folder to speed
// things along; note that there are a couple of very low frequencies (~5-6 cm^-1 for this case)
        else if (attemptNumber % scriptAttempts == 6)
            inpKeyStr += "# " + qmMethod + " opt=(verytight,nolinear,calcfc,small) freq IOP(2/16=3)";// used for troublesome C5H7J2
// case (similar error to C5H7J below); calcfc is not necessary for this particular species, but it speeds convergence
// and probably makes it more robust for other species
        else if (attemptNumber % scriptAttempts == 7)
            inpKeyStr += "# " + qmMethod + " opt=(verytight,gdiis,maxcyc=200) freq=numerical IOP(2/16=3)"; // use numerical
// frequencies; this takes a relatively long time, so should only be used as one of the last resorts; this seemed to
// address at least one case of failure for a C6H10JJ species; 7/15/09: maxcyc=200 added to address
// GVCMURUDAUQXEY-UHFFFAOYAVmult3 (InChI=1/C3H4O7Si/c1-2(9-6)10-11(7,8)3(4)5/h6-7H,1H2/mult3)...however, result was
// manually pasted in QMfiles folder to speed things along
        else if (attemptNumber % scriptAttempts == 8)
            inpKeyStr += "# " + qmMethod + " opt=tight freq IOP(2/16=3)";// 7/10/09: this worked for problematic case of
// SZSSHFMXPBKYPR-UHFFFAOYAF (InChI=1/C7H15O5Si/c1-3-10-13(8,11-4-2)12-7-5-6-9-7/h7H,3-6H2,1-2H3) (otherwise, it had
// l402.exe errors); corrected log file was manually copied to QMfiles to speed things along; we could also add a
// freq=numerical version of this keyword combination for added robustness; UPDATE: see below
        else if (attemptNumber % scriptAttempts == 9)
            inpKeyStr += "# " + qmMethod + " opt=tight freq=numerical IOP(2/16=3)";// 7/10/09: used for problematic case of
// CIKDVMUGTARZCK-UHFFFAOYAImult4 (InChI=1/C8H15O6Si/c1-4-12-15(10,13-5-2)14-7-6-11-8(7,3)9/h7H,3-6H2,1-2H3/mult4 (most
// other cases had l402.exe errors); corrected log file was manually copied to QMfiles to speed things along
        else if (attemptNumber % scriptAttempts == 10)
            inpKeyStr += "# " + qmMethod + " opt=(tight,nolinear,calcfc,small,maxcyc=200) freq IOP(2/16=3)";// 7/8/09: similar to
// existing #5, but uses tight rather than verytight; used for ADMPQLGIEMRGAT-UHFFFAOYAUmult3
// (InChI=1/C6H14O5Si/c1-4-9-12(8,10-5-2)11-6(3)7/h6-7H,3-5H2,1-2H3/mult3)
        else if (attemptNumber % scriptAttempts == 11)
            inpKeyStr += "# " + qmMethod + " opt freq IOP(2/16=3)"; // use default (not verytight) convergence criteria; use this as
// last resort
        else if (attemptNumber % scriptAttempts == 12)
            inpKeyStr += "# " + qmMethod + " opt=(verytight,gdiis) freq=numerical IOP(2/16=3) IOP(4/21=200)";// to address
// problematic C10H14JJ case
        else if (attemptNumber % scriptAttempts == 13)
            inpKeyStr += "# " + qmMethod + " opt=(calcfc,verytight,newton,notrustupdate,small,maxcyc=100,maxstep=100) freq=(numerical,step=10) IOP(2/16=3) nosymm";// added
// 6/10/09 for very troublesome RRMZRNPRCUANER-UHFFFAOYAQ (InChI=1/C5H7/c1-3-5-4-2/h3H,1-2H3) case...there were troubles
// with negative frequencies, where I don't think they should have been; step size of numerical frequency was adjusted
// to give positive result; accuracy of result is questionable; it is possible that not all of these keywords are
// needed; note that for this and other nearly free rotor cases, I think heat capacity will be overestimated by R/2 (R
// vs. R/2) (but this is a separate issue)
        else if (attemptNumber % scriptAttempts == 14)
            inpKeyStr += "# " + qmMethod + " opt=(tight,gdiis,small,maxcyc=200,maxstep=100) freq=numerical IOP(2/16=3) nosymm";// added
// 6/22/09 for troublesome
// QDERTVAGQZYPHT-UHFFFAOYAHmult3(InChI=1/C6H14O4Si/c1-4-8-11(7,9-5-2)10-6-3/h4H,5-6H2,1-3H3/mult3); key aspects appear
// to be tight (rather than verytight) convergence criteria, no calculation of frequencies during optimization, use of
// numerical frequencies, and probably also the use of opt=small
        // gmagoon 7/9/09: commented out since although this produces a "reasonable" result for the problematic case,
// there is a large amount of spin contamination, apparently introducing 70+ kcal/mol of instability else
// if(attemptNumber==12)
// inpKeyStr+="# pm3 opt=(verytight,gdiis,small) freq=numerical IOP(2/16=3) IOP(4/21=200)";//7/9/09: similar to current
// number 9 with keyword small; this addresses case of VCSJVABXVCFDRA-UHFFFAOYAI
// (InChI=1/C8H19O5Si/c1-5-10-8(4)13-14(9,11-6-2)12-7-3/h8H,5-7H2,1-4H3)
        else if (attemptNumber % scriptAttempts == 15)
            inpKeyStr += "# " + qmMethod + " opt=(verytight,gdiis,calcall) IOP(2/16=3)";// used for troublesome C5H7J case; note that
// before fixing, I got errors like the following:
// "Incomplete coordinate system.  Try restarting with Geom=Check Guess=Read Opt=(ReadFC,NewRedundant) Incomplete coordinate system. Error termination via Lnk1e in l103.exe";
// we could try to restart, but it is probably preferrable to have each keyword combination standalone; another keyword
// that may be helpful if additional problematic cases are encountered is opt=small; 6/9/09 note: originally, this had #
// pm3 opt=(verytight,gdiis,calcall) freq IOP(2/16=3)" (with freq keyword), but I discovered that in this case, there
// are two thermochemistry sections and cclib parses frequencies twice, giving twice the number of desired frequencies
// and hence produces incorrect thermo; this turned up on C5H6JJ isomer
// gmagoon 7/3/09: it is probably best to retire this keyword combination in light of the similar combination below
// //else if(attemptNumber==6) inpKeyStr+="# pm3 opt=(verytight,gdiis,calcall,small) IOP(2/16=3) IOP(4/21=2)";//6/10/09:
// worked for OJZYSFFHCAPVGA-UHFFFAOYAK (InChI=1/C5H7/c1-3-5-4-2/h1,4H2,2H3) case; IOP(4/21) keyword was key
        else if (attemptNumber % scriptAttempts == 16)
            inpKeyStr += "# " + qmMethod + " opt=(verytight,gdiis,calcall,small,maxcyc=200) IOP(2/16=3) IOP(4/21=2) nosymm";// 6/29/09:
// worked for troublesome ketene case: CCGKOQOJPYTBIH-UHFFFAOYAO (InChI=1/C2H2O/c1-2-3/h1H2) (could just increase number
// of iterations for similar keyword combination above (#6 at the time of this writing), allowing symmetry, but nosymm
// seemed to reduce # of iterations; I think one of nosymm or higher number of iterations would allow the similar
// keyword combination to converge; both are included here for robustness)
        else if (attemptNumber % scriptAttempts == 17)
            inpKeyStr += "# " + qmMethod + " opt=(verytight,gdiis,calcall,small) IOP(2/16=3) nosymm";// 7/1/09: added for case of
// ZWMVZWMBTVHPBS-UHFFFAOYAEmult3 (InChI=1/C4H4O2/c1-3-5-6-4-2/h1-2H2/mult3)
        else if (attemptNumber % scriptAttempts == 0)
            inpKeyStr += "# " + qmMethod + " opt=(calcall,small,maxcyc=100) IOP(2/16=3)"; // 6/10/09: used to address troublesome
// FILUFGAZMJGNEN-UHFFFAOYAImult3 case (InChI=1/C5H6/c1-3-5-4-2/h3H,1H2,2H3/mult3)
        else {// this point should not be reached
            Logger.error("Unexpected error in determining which G03 PM3 keyword set to use");
            System.exit(0);
        }
        // if(m…

Large files files are truncated, but you can click here to view the full file