/source/RMG/jing/chem/QMTP.java
Java | 3717 lines | 2641 code | 57 blank | 1019 comment | 766 complexity | 7d2a8e4aaaf0e894bcd0c17001f4c8b5 MD5 | raw file
Possible License(s): LGPL-2.1
- // //////////////////////////////////////////////////////////////////////////////
- //
- // 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(multiplicity == 3) inpKeyStr+= " guess=mix"; //assumed to be triplet biradical...use guess=mix to perform
- // unrestricted ; nevermind...I think this would only be for singlet biradicals based on
- // http://www.gaussian.com/g_tech/g_ur/k_guess.htm
- if (usePolar)
- inpKeyStr += " polar";
- // call the OpenBabel process (note that this requires OpenBabel environment variable)
- try {
- File runningdir = new File(directory);
- String molPath = null;
- if (attemptNumber <= scriptAttempts) {// use UFF-refined coordinates
- molPath = p_molfile.getPath();
- } else {// use crude coordinates
- molPath = p_molfile.getCrudePath();
- }
- String[] command = { "babel", "-imol", molPath, "-ogjf",
- name + ".gjf", "-xk", inpKeyStr, "--title", InChIaug };
- Process babelProc = Runtime.getRuntime().exec(command, null,
- runningdir);
- // read in output
- InputStream is = babelProc.getInputStream();
- InputStreamReader isr = new InputStreamReader(is);
- BufferedReader br = new BufferedReader(isr);
- String line = null;
- while ((line = br.readLine()) != null) {
- // do nothing
- }
- int exitValue = babelProc.waitFor();
- babelProc.getErrorStream().close();
- babelProc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- String err = "Error in running OpenBabel MOL to GJF process \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- return maxAttemptNumber;
- }
- // creates MM4 input file and MM4 batch file in directory with filenames name.mm4 and name.com, respectively using
- // MoleCoor
- // 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
- public int createMM4Input(String name, String directory, molFile p_molfile,
- int attemptNumber, String InChIaug, int multiplicity) {
- // Step 1: write the script for MM4 batch operation
- // Example script file:
- // #! /bin/csh
- // cp testEthylene.mm4 CPD.MM4
- // cp $MM4_DATDIR/BLANK.DAT PARA.MM4
- // cp $MM4_DATDIR/CONST.MM4 .
- // $MM4_EXEDIR/mm4 <<%
- // 1
- // 2
- // 0
- // %
- // mv TAPE4.MM4 testEthyleneBatch.out
- // mv TAPE9.MM4 testEthyleneBatch.opt
- // exit
- int scriptAttempts = 2;// the number of script 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
- try {
- // create batch file with executable permissions: cf.
- // http://java.sun.com/docs/books/tutorial/essential/io/fileAttr.html#posix
- File inpKey = new File(directory + "/" + name + ".com");
- String inpKeyStr = "#! /bin/csh\n";
- inpKeyStr += "cp " + name + ".mm4 CPD.MM4\n";
- inpKeyStr += "cp $MM4_DATDIR/BLANK.DAT PARA.MM4\n";
- inpKeyStr += "cp $MM4_DATDIR/CONST.MM4 .\n";
- inpKeyStr += "$MM4_EXEDIR/mm4 <<%\n";
- inpKeyStr += "1\n";// read from first line of .mm4 file
- if (!useCanTherm) {
- if (attemptNumber % scriptAttempts == 1)
- inpKeyStr += "2\n"; // Block-Diagonal Method then Full-Matrix Method
- else if (attemptNumber % scriptAttempts == 0)
- inpKeyStr += "3\n"; // Full-Matrix Method only
- else
- throw new Exception();// this point should not be reached
- } else {// CanTherm case: write the FORCE.MAT file
- if (attemptNumber % scriptAttempts == 1)
- inpKeyStr += "4\n"; // Block-Diagonal Method then Full-Matrix Method
- else if (attemptNumber % scriptAttempts == 0)
- inpKeyStr += "5\n"; // Full-Matrix Method only
- else
- throw new Exception();// this point should not be reached
- inpKeyStr += "\n";// <RETURN> for temperature
- inpKeyStr += "4\n";// unofficial option 4 for vibrational eigenvector printout to generate Cartesian
- // force constant matrix in FORCE.MAT file
- inpKeyStr += "0\n";// no vibrational amplitude printout
- }
- inpKeyStr += "0\n";// terminate the job
- inpKeyStr += "%\n";
- inpKeyStr += "mv TAPE4.MM4 " + name + ".mm4out\n";
- inpKeyStr += "mv TAPE9.MM4 " + name + ".mm4opt\n";
- if (useCanTherm) {
- inpKeyStr += "mv FORCE.MAT " + name + ".fmat\n";
- }
- inpKeyStr += "exit\n";
- FileWriter fw = new FileWriter(inpKey);
- fw.write(inpKeyStr);
- fw.close();
- } catch (Exception e) {
- String err = "Error in writing MM4 script file\n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // Step 2: call the MoleCoor process to create the MM4 input file from the mole file
- try {
- File runningdir = new File(directory);
- // this will only be run on Linux so we don't have to worry about Linux vs. Windows issues
- String command = "python " + System.getenv("RMG")
- + "/scripts/MM4InputFileMaker.py ";
- // first argument: input file path; for the first attempts, we will use UFF refined coordinates; if that
- // doesn't work, (e.g. case of cyclopropene, InChI=1/C3H4/c1-2-3-1/h1-2H,3H2 OOXWYYGXTJLWHA-UHFFFAOYAJ) we will try
- // crude coordinates (.cmol suffix)
- if (attemptNumber <= scriptAttempts) {
- command = command.concat(p_molfile.getPath() + " ");
- } else {
- command = command.concat(p_molfile.getCrudePath() + " ");
- }
- // second argument: output path
- String inpfilepath = directory + "/" + name + ".mm4";
- command = command.concat(inpfilepath + " ");
- // third argument: molecule name (the augmented InChI)
- command = command.concat(InChIaug + " ");
- // fourth argument: PYTHONPATH
- command = command.concat(System.getenv("RMG") + "/source/MoleCoor");// this will pass $RMG/source/MoleCoor
- // to the script (in order to get the appropriate path for importing
- Process molecoorProc = Runtime.getRuntime().exec(command, null,
- runningdir);
- // read in output
- InputStream is = molecoorProc.getInputStream();
- InputStreamReader isr = new InputStreamReader(is);
- BufferedReader br = new BufferedReader(isr);
- String line = null;
- while ((line = br.readLine()) != null) {
- // do nothing
- }
- int exitValue = molecoorProc.waitFor();
- molecoorProc.getErrorStream().close();
- molecoorProc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- String err = "Error in running MoleCoor MOL to .MM4 process \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- return maxAttemptNumber;
- }
- // creates MM4 rotor input file and MM4 batch file in directory with filenames name.mm4roti and name.comi,
- // respectively
- // the function returns the set of rotor dihedral angles for the minimum energy conformation
- public double[] createMM4RotorInput(String name, String directory,
- ChemGraph p_chemgraph, int rotors) {
- // read in the optimized coordinates from the completed "normal" MM4 job; this will be used as a template for
- // the rotor input files
- String mm4optContents = "";
- int indexForFirstAtom = -1;
- int lineIndex = 0;
- double[] dihedralMinima = new double[rotors];
- try {
- FileReader mm4opt = new FileReader(directory + "/" + name
- + ".mm4opt");
- BufferedReader reader = new BufferedReader(mm4opt);
- String line = reader.readLine();
- while (line != null) {
- if (indexForFirstAtom < 0 && line.length() >= 40
- && line.substring(35, 40).equals("( 1)"))
- indexForFirstAtom = lineIndex;// store the location of the first atom coordinate
- mm4optContents += line + "\n";
- line = reader.readLine();
- lineIndex++;
- }
- reader.close();
- mm4opt.close();
- } catch (Exception e) {
- String err = "Error in reading .mm4opt file\n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- String[] lines = mm4optContents.split("[\\r?\\n]+");// split by newlines, excluding blank lines; cf.
- // http://stackoverflow.com/questions/454908/split-java-string-by-new-line
- int flag = 0;// flag to indicate whether there is a "line 1a" with pi system information
- if (lines[1].startsWith("T") || lines[1].startsWith("F"))
- flag = 1;// pi system information is indicated by the second line beginning with T's and potentially F's
- lines[1 + flag] = lines[1 + flag].substring(0, 78) + " 2";// take the first 78 characters of line 2 (or 3 if it
- // is a pi-system compound), and append the option number for the NDRIVE option; in other words, we are replacing the
- // NDRIVE=0 option with the desired option number
- // reconstruct mm4optContents
- mm4optContents = "";
- for (int j = 0; j < lines.length; j++) {
- mm4optContents += lines[j] + "\n";
- }
- // iterate over all the rotors in the molecule
- int i = 0;// rotor index
- LinkedHashMap rotorInfo = p_chemgraph.getInternalRotorInformation();
- Iterator iter = rotorInfo.keySet().iterator();
- while (iter.hasNext()) {
- i++;
- int[] rotorAtoms = (int[]) iter.next();
- try {
- // write one script file for each rotor
- // Step 1: write the script for MM4 batch operation
- // Example script file:
- // #! /bin/csh
- // cp testEthylene.mm4 CPD.MM4
- // cp $MM4_DATDIR/BLANK.DAT PARA.MM4
- // cp $MM4_DATDIR/CONST.MM4 .
- // $MM4_EXEDIR/mm4 <<%
- // 1
- // 2
- // 0
- // %
- // mv TAPE4.MM4 testEthyleneBatch.out
- // mv TAPE9.MM4 testEthyleneBatch.opt
- // exit
- // create batch file with executable permissions: cf.
- // http://java.sun.com/docs/books/tutorial/essential/io/fileAttr.html#posix
- File inpKey = new File(directory + "/" + name + ".com" + i);
- String inpKeyStr = "#! /bin/csh\n";
- inpKeyStr += "cp " + name + ".mm4rot" + i + " CPD.MM4\n";
- inpKeyStr += "cp $MM4_DATDIR/BLANK.DAT PARA.MM4\n";
- inpKeyStr += "cp $MM4_DATDIR/CONST.MM4 .\n";
- inpKeyStr += "$MM4_EXEDIR/mm4 <<%\n";
- inpKeyStr += "1\n";// read from first line of .mm4 file
- inpKeyStr += "3\n"; // Full-Matrix Method only
- // inpKeyStr+="2\n"; //Block-Diagonal Method then Full-Matrix Method
- inpKeyStr += "0\n";// terminate the job
- inpKeyStr += "%\n";
- inpKeyStr += "mv TAPE4.MM4 " + name + ".mm4rotout" + i + "\n";
- inpKeyStr += "mv TAPE9.MM4 " + name + ".mm4rotopt" + i + "\n";
- inpKeyStr += "exit\n";
- FileWriter fw = new FileWriter(inpKey);
- fw.write(inpKeyStr);
- fw.close();
- } catch (Exception e) {
- String err = "Error in writing MM4 script files\n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // Step 2: create the MM4 rotor job input file for rotor i, using the output file from the "normal" MM4 job
- // as a template
- // Step 2a: find the dihedral angle for the minimized conformation (we need to start at the minimum energy
- // conformation because CanTherm assumes that theta=0 is the minimum
- // extract the lines for dihedral atoms
- String dihedral1s = lines[indexForFirstAtom + rotorAtoms[0] - 1];
- String atom1s = lines[indexForFirstAtom + rotorAtoms[1] - 1];
- String atom2s = lines[indexForFirstAtom + rotorAtoms[2] - 1];
- String dihedral2s = lines[indexForFirstAtom + rotorAtoms[3] - 1];
- // extract the x,y,z coordinates for each atom
- double[] dihedral1 = {
- Double.parseDouble(dihedral1s.substring(0, 10)),
- Double.parseDouble(dihedral1s.substring(10, 20)),
- Double.parseDouble(dihedral1s.substring(20, 30)) };
- double[] atom1 = { Double.parseDouble(atom1s.substring(0, 10)),
- Double.parseDouble(atom1s.substring(10, 20)),
- Double.parseDouble(atom1s.substring(20, 30)) };
- double[] atom2 = { Double.parseDouble(atom2s.substring(0, 10)),
- Double.parseDouble(atom2s.substring(10, 20)),
- Double.parseDouble(atom2s.substring(20, 30)) };
- double[] dihedral2 = {
- Double.parseDouble(dihedral2s.substring(0, 10)),
- Double.parseDouble(dihedral2s.substring(10, 20)),
- Double.parseDouble(dihedral2s.substring(20, 30)) };
- // determine the dihedral angle
- dihedralMinima[i - 1] = calculateDihedral(dihedral1, atom1, atom2,
- dihedral2);
- // double dihedral = calculateDihedral(dihedral1,atom1,atom2,dihedral2);
- // if (dihedral < 0) dihedral = dihedral + 360;//make sure dihedral is positive; this will save an extra
- // character in the limited space for specifying starting and ending angles
- // eventually, when problems arise due to collinear atoms (both arguments to atan2 are zero) we can iterate
- // over other atom combinations (with atoms in each piece determined by the corresponding value in rotorInfo for atom2
- // and the complement of these (full set = p_chemgraph.getNodeIDs()) for atom1) until they are not collinear
- // Step 2b: write the file for rotor i
- try {
- FileWriter mm4roti = new FileWriter(directory + "/" + name
- + ".mm4rot" + i);
- // mm4roti.write(mm4optContents+"\n"+String.format(" %3d %3d %3d %3d %5.1f%5.1f%5.1f",
- // rotorAtoms[0],rotorAtoms[1],rotorAtoms[2],rotorAtoms[3], dihedral, dihedral + 360.0 - deltaTheta,
- // deltaTheta)+"\n");//deltaTheta should be less than 100 degrees so that dihedral-deltaTheta still fits
- // mm4roti.write(mm4optContents+"\n"+String.format(" %3d %3d %3d %3d %5.1f%5.1f%5.1f",
- // rotorAtoms[0],rotorAtoms[1],rotorAtoms[2],rotorAtoms[3], dihedral, dihedral - deltaTheta,
- // deltaTheta)+"\n");//deltaTheta should be less than 100 degrees so that dihedral-deltaTheta still fits
- mm4roti.write(mm4optContents
- + String.format(
- " %3d %3d %3d %3d %5.1f%5.1f%5.1f",
- rotorAtoms[0], rotorAtoms[1], rotorAtoms[2],
- rotorAtoms[3], 0.0, 360.0 - deltaTheta,
- deltaTheta) + "\n");// M1, M2, M3, M4, START, FINISH, DIFF (as described in MM4 manual);
- // //this would require miminum to be stored and used to adjust actual angles before sending to CanTherm
- mm4roti.close();
- } catch (Exception e) {
- String err = "Error in writing MM4 rotor input file\n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- }
- return dihedralMinima;
- }
- // given x,y,z (cartesian) coordinates for dihedral1 atom, (rotor) atom 1, (rotor) atom 2, and dihedral2 atom,
- // calculates the dihedral angle (in degrees, between +180 and -180) using the atan2 formula at
- // http://en.wikipedia.org/w/index.php?title=Dihedral_angle&oldid=373614697
- public double calculateDihedral(double[] dihedral1, double[] atom1,
- double[] atom2, double[] dihedral2) {
- // calculate the vectors between the atoms
- double[] b1 = { atom1[0] - dihedral1[0], atom1[1] - dihedral1[1],
- atom1[2] - dihedral1[2] };
- double[] b2 = { atom2[0] - atom1[0], atom2[1] - atom1[1],
- atom2[2] - atom1[2] };
- double[] b3 = { dihedral2[0] - atom2[0], dihedral2[1] - atom2[1],
- dihedral2[2] - atom2[2] };
- // calculate norm of b2
- double normb2 = Math
- .sqrt(b2[0] * b2[0] + b2[1] * b2[1] + b2[2] * b2[2]);
- // calculate necessary cross products
- double[] b1xb2 = { b1[1] * b2[2] - b1[2] * b2[1],
- b2[0] * b1[2] - b1[0] * b2[2], b1[0] * b2[1] - b1[1] * b2[0] };
- double[] b2xb3 = { b2[1] * b3[2] - b2[2] * b3[1],
- b3[0] * b2[2] - b2[0] * b3[2], b2[0] * b3[1] - b2[1] * b3[0] };
- // compute arguments for atan2 function (includes dot products)
- double y = normb2
- * (b1[0] * b2xb3[0] + b1[1] * b2xb3[1] + b1[2] * b2xb3[2]);// |b2|*b1.(b2xb3)
- double x = b1xb2[0] * b2xb3[0] + b1xb2[1] * b2xb3[1] + b1xb2[2]
- * b2xb3[2];// (b1xb2).(b2xb3)
- double dihedral = Math.atan2(y, x);
- // return dihedral*180.0/Math.PI;//converts from radians to degrees
- return Math.toDegrees(dihedral);// converts from radians to degrees
- }
- // returns the extra Mopac keywords to use for radical species, given the spin multiplicity (radical number + 1)
- public String getMopacRadicalString(int multiplicity) {
- if (multiplicity == 1)
- return "";
- else if (multiplicity == 2)
- return "uhf doublet";
- else if (multiplicity == 3)
- return "uhf triplet";
- else if (multiplicity == 4)
- return "uhf quartet";
- else if (multiplicity == 5)
- return "uhf quintet";
- else if (multiplicity == 6)
- return "uhf sextet";
- else if (multiplicity == 7)
- return "uhf septet";
- else if (multiplicity == 8)
- return "uhf octet";
- else if (multiplicity == 9)
- return "uhf nonet";
- else {
- Logger.critical("Invalid multiplicity encountered: " + multiplicity);
- System.exit(0);
- }
- return "this should not be returned: error associated with getMopacRadicalString()";
- }
- // creates MOPAC input file in directory with filename name.mop 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
- // unlike createGaussianPM3 input, this requires an additional input specifying the spin multiplicity (radical
- // number + 1) for the species
- public int createMopacInput(String name, String directory,
- molFile p_molfile, int attemptNumber, String InChIaug,
- int multiplicity, String qmMethod) {
- // write a file with the input keywords
- int scriptAttempts = 5;// 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 inpKeyStrBoth = "";// this string will be written at both the top (for optimization) and the bottom (for
- // thermo/force calc)
- String inpKeyStrTop = "";// this string will be written only at the top
- String inpKeyStrBottom = "";// this string will be written at the bottom
- String radicalString = getMopacRadicalString(multiplicity);
- if (attemptNumber % scriptAttempts == 1) {
- inpKeyStrBoth = qmMethod + " " + radicalString;
- inpKeyStrTop = " precise nosym";
- inpKeyStrBottom = "oldgeo thermo nosym precise ";// 7/10/09: based on a quick review of recent results,
- // keyword combo #1 rarely works, and when it did (CJAINEUZFLXGFA-UHFFFAOYAUmult3
- // (InChI=1/C8H16O5Si/c1-4-11-14(9,12-5-2)13-8-6-10-7(8)3/h7-8H,3-6H2,1-2H3/mult3)), the grad. norm on the force step
- // was about 1.7 (too large); I manually removed this result and re-ran...the entropy was increased by nearly 20
- // cal/mol-K...perhaps we should add a check for the "WARNING" that MOPAC prints out when the gradient is high; 7/22/09:
- // for the case of FUGDBSHZYPTWLG-UHFFFAOYADmult3 (InChI=1/C5H8/c1-4-3-5(4)2/h4-5H,1-3H2/mult3), adding nosym seemed to
- // resolve 1. large discrepancies from Gaussian and 2. negative frequencies in mass-weighted coordinates and possibly
- // related issue in discrepancies between regular and mass-weighted coordinate frequencies
- } else if (attemptNumber % scriptAttempts == 2) {// 7/9/09: used for VCSJVABXVCFDRA-UHFFFAOYAI
- // (InChI=1/C8H19O5Si/c1-5-10-8(4)13-14(9,11-6-2)12-7-3/h8H,5-7H2,1-4H3); all existing Gaussian keywords also failed;
- // the Gaussian result was also rectified, but the resulting molecule was over 70 kcal/mol less stable, probably due to
- // a large amount of spin contamination (~1.75 in fixed Gaussian result vs. 0.754 for MOPAC)
- inpKeyStrBoth = qmMethod + " " + radicalString;
- inpKeyStrTop = " precise nosym gnorm=0.0 nonr";
- inpKeyStrBottom = "oldgeo thermo nosym precise ";
- } else if (attemptNumber % scriptAttempts == 3) {// 7/8/09: 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); all existing Gaussian keywords also failed;
- // however, the Gaussian result was also rectified, and the resulting conformation was about 1.0 kcal/mol more stable
- // than the one resulting from this, so fixed Gaussian result was manually copied to QMFiles folder
- inpKeyStrBoth = qmMethod + " " + radicalString;
- inpKeyStrTop = " precise nosym gnorm=0.0";
- inpKeyStrBottom = "oldgeo thermo nosym precise "; // precise appeared to be necessary for the problematic
- // case (to avoid negative frequencies);
- } else if (attemptNumber % scriptAttempts == 4) {// 7/8/09: used for 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 (negative frequency issues in MOPAC)
- // (also, none of the existing Gaussian combinations worked with it); note that the Gaussian result appears to be a
- // different conformation as it is about 0.85 kcal/mol more stable, so the Gaussian result was manually copied to
- // QMFiles directory; note that the MOPAC output included a very low frequency (4-5 cm^-1)
- inpKeyStrBoth = qmMethod + " " + radicalString;
- inpKeyStrTop = " precise nosym gnorm=0.0 bfgs";
- inpKeyStrBottom = "oldgeo thermo nosym precise "; // precise appeared to be necessary for the problematic
- // case (to avoid negative frequencies)
- }
- // else if(attemptNumber==5){
- // inpKeyStrBoth="pm3 "+radicalString;
- // inpKeyStrTop=" precise nosym gnorm=0.0 ddmin=0.0";
- // inpKeyStrBottom="oldgeo thermo nosym precise ";
- // }
- // else if(attemptNumber==6){
- // inpKeyStrBoth="pm3 "+radicalString;
- // inpKeyStrTop=" precise nosym gnorm=0.0 nonr ddmin=0.0";
- // inpKeyStrBottom="oldgeo thermo nosym precise ";
- // }
- // else if(attemptNumber==7){
- // inpKeyStrBoth="pm3 "+radicalString;
- // inpKeyStrTop=" precise nosym bfgs gnorm=0.0 ddmin=0.0";
- // inpKeyStrBottom="oldgeo thermo nosym precise ";
- // }
- else if (attemptNumber % scriptAttempts == 0) {// used for troublesome HGRZRPHFLAXXBT-UHFFFAOYAVmult3
- // (InChI=1/C3H2O4/c4-2(5)1-3(6)7/h1H2/mult3) case (negative frequency and large gradient issues)
- inpKeyStrBoth = qmMethod + " " + radicalString;
- inpKeyStrTop = " precise nosym recalc=10 dmax=0.10 nonr cycles=2000 t=2000";
- inpKeyStrBottom = "oldgeo thermo nosym precise ";
- }
- // else if(attemptNumber==9){//used for troublesome CMARQPBQDRXBTN-UHFFFAOYAAmult3
- // (InChI=1/C3H2O4/c1-3(5)7-6-2-4/h1H2/mult3) case (negative frequency issues)
- // inpKeyStrBoth="pm3 "+radicalString;
- // inpKeyStrTop=" precise nosym recalc=1 dmax=0.05 gnorm=0.0 cycles=1000 t=1000";
- // inpKeyStrBottom="oldgeo thermo nosym precise ";
- // }
- // else if(attemptNumber==10){//used for ATCYLHQLTOSVFK-UHFFFAOYAMmult4
- // (InChI=1/C4H5O5/c1-3(5)8-9-4(2)7-6/h6H,1-2H2/mult4) case (timeout issue; also, negative frequency issues); note that
- // this is very similar to the keyword below, so we may want to consolidate
- // inpKeyStrBoth="pm3 "+radicalString;
- // inpKeyStrTop=" precise nosym recalc=1 dmax=0.05 gnorm=0.2 cycles=1000 t=1000";
- // inpKeyStrBottom="oldgeo thermo nosym precise ";
- // }
- else {// this point should not be reached
- Logger.error("Unexpected error in determining which MOPAC keyword set to use");
- System.exit(0);
- }
- String polarString = "";
- if (usePolar) {
- if (multiplicity == 1)
- polarString = System.getProperty("line.separator")
- + System.getProperty("line.separator")
- + System.getProperty("line.separator")
- + "oldgeo polar nosym precise " + inpKeyStrBoth;
- else
- polarString = System.getProperty("line.separator")
- + System.getProperty("line.separator")
- + System.getProperty("line.separator")
- + "oldgeo static nosym precise " + inpKeyStrBoth;
- }
- // call the OpenBabel process (note that this requires OpenBabel environment variable)
- try {
- File runningdir = new File(directory);
- String inpKeyStrTopCombined = inpKeyStrBoth + inpKeyStrTop;
- Process babelProc = null;
- if (attemptNumber <= scriptAttempts) {// use UFF-refined coordinates
- String[] command = { "babel", "-imol", p_molfile.getPath(),
- "-xk", inpKeyStrTopCombined, "--title", InChIaug,
- "-omop", name + ".mop" };
- babelProc = Runtime.getRuntime()
- .exec(command, null, runningdir);
- } else {// use the crude coordinates
- String[] command = { "babel", "-imol",
- p_molfile.getCrudePath(), "-xk", inpKeyStrTopCombined,
- "--title", InChIaug, "-omop", name + ".mop" };
- babelProc = Runtime.getRuntime()
- .exec(command, null, runningdir);
- }
- // read in output
- InputStream is = babelProc.getInputStream();
- InputStreamReader isr = new InputStreamReader(is);
- BufferedReader br = new BufferedReader(isr);
- String line = null;
- while ((line = br.readLine()) != null) {
- // do nothing
- }
- int exitValue = babelProc.waitFor();
- babelProc.getErrorStream().close();
- babelProc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- // append the final keywords to the end of the file just written
- // File mopacInpFile = new File(directory+"/"+name+".mop");
- FileWriter fw = new FileWriter(directory + "/" + name + ".mop",
- true);// filewriter with append = true
- fw.write(System.getProperty("line.separator") + inpKeyStrBottom
- + inpKeyStrBoth + polarString);// on Windows Vista, "\n" appears correctly in WordPad view, but not
- // Notepad view (cf. http://forums.sun.com/thread.jspa?threadID=5386822)
- fw.close();
- } catch (Exception e) {
- String err = "Error in running OpenBabel MOL to MOP process \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- return maxAttemptNumber;
- }
- // name and directory are the name and directory for the input (and output) file;
- // input is assumed to be preexisting and have the .gjf suffix
- // returns an integer indicating success or failure of the Gaussian calculation: 1 for success, 0 for failure;
- public int runGaussian(String name, String directory, String InChIaug) {
- int flag = 0;
- int successFlag = 0;
- try {
- String command = "g03 ";
- 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 = command.concat("\"" + qmfolder + "/" + name
- + ".gjf\" ");// specify the input file; space is important
- command = command.concat("\"" + qmfolder + "/" + name
- + ".log\"");// specify the output file
- } else {// non-Windows case
- command = command.concat(qmfolder + "/" + name + ".gjf ");// specify the input file; space is important
- command = command.concat(qmfolder + "/" + name + ".log");// specify the output file
- }
- Process gaussianProc = Runtime.getRuntime().exec(command);
- // check for errors and display the error if there is one
- InputStream is = gaussianProc.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 that an error was obtained
- if (flag == 1) {
- Logger.info("Gaussian process received error (see above) on "
- + name);
- }
- int exitValue = gaussianProc.waitFor();
- gaussianProc.getInputStream().close();
- gaussianProc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- String err = "Error in running Gaussian process \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // look in the output file to check for the successful termination of the Gaussian calculation
- // failed jobs will contain the a line beginning with " Error termination" near the end of the file
- int failureFlag = 0;
- int completeFlag = 0;
- String errorLine = "";// string to store the error
- try {
- FileReader in = new FileReader(directory + "/" + name + ".log");
- BufferedReader reader = new BufferedReader(in);
- String line = reader.readLine();
- while (line != null) {
- if (line.startsWith(" Error termination ")) {
- failureFlag = 1;
- errorLine = line.trim();
- Logger.info("*****Error in Gaussian log file: " + errorLine);// print the error (note that in
- // general, I think two lines will be printed)
- } else if (line.startsWith(" ******")) {// also look for imaginary frequencies
- if (line.contains("imaginary frequencies")) {
- Logger.info("*****Imaginary freqencies found:");
- failureFlag = 1;
- }
- } else if (line
- .startsWith(" Normal termination of Gaussian 03 at")) {// check for completion notice at end of
- // file
- if (reader.readLine() == null) {// the above line will occur once in the middle of most files
- // (except for monatomic species), but the following line is non-critical (e.g.
- // "Link1: Proceeding to internal job step number 2.") so reading over this line should not cause a problem
- completeFlag = 1;
- }
- }
- line = reader.readLine();
- }
- reader.close();
- in.close();
- } catch (Exception e) {
- String err = "Error in reading Gaussian log file \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // if the complete flag is still 0, the process did not complete and is a failure
- if (completeFlag == 0)
- failureFlag = 1;
- if (failureFlag == 0 && connectivityCheck > 0) {// if the process was successful, check the connectivity (if
- // requested by the user)
- boolean connectivityMatch = connectivityMatchInPM3ResultQ(name,
- directory, InChIaug, ".log", "g03");
- if (!connectivityMatch)
- if (connectivityCheck > 1) {// connectivityCheck=2
- failureFlag = 1;
- } else {// connectivityCheck=1; print warning and procede with "return true"
- Logger.warning("Connectivity in quantum result for "
- + name
- + " ("
- + InChIaug
- + ") does not appear to match desired connectivity.");
- }
- }
- // if the failure flag is still 0, the process should have been successful
- if (failureFlag == 0)
- successFlag = 1;
- return successFlag;
- }
- // name and directory are the name and directory for the input (and output) file;
- // input script is assumed to be preexisting and have the .com suffix
- // returns an integer indicating success or failure of the calculation: 1 for success, 0 for failure
- public int runMM4(String name, String directory, String InChIaug) {
- int successFlag = 0;
- // int flag = 0;
- try {
- File runningDirectory = new File(qmfolder);
- String command = name + ".com";
- // File script = new File(qmfolder+command);
- // command = "./"+command;
- // script.setExecutable(true);
- command = "csh " + command;
- Process mm4Proc = Runtime.getRuntime().exec(command, null,
- runningDirectory);
- // check for errors and display the error if there is one
- // InputStream is = mm4Proc.getErrorStream();
- // InputStreamReader isr = new InputStreamReader(is);
- // BufferedReader br = new BufferedReader(isr);
- // String line=null;
- // while ( (line = br.readLine()) != null) {
- // line = line.trim();
- // if(!line.equals("STOP statement executed")){//string listed here seems to be typical
- // Logger.error(line);
- // flag=1;
- // }
- // }
- InputStream is = mm4Proc.getInputStream();
- InputStreamReader isr = new InputStreamReader(is);
- BufferedReader br = new BufferedReader(isr);
- String line = null;
- while ((line = br.readLine()) != null) {
- // do nothing
- }
- // if there was an error, indicate that an error was obtained
- // if(flag==1){
- // Logger.info("MM4 process received error (see above) on " + name);
- // }
- int exitValue = mm4Proc.waitFor();
- mm4Proc.getErrorStream().close();
- mm4Proc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- String err = "Error in running MM4 process \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // look in the output file to check for the successful termination of the MM4 calculation (cf.
- // successfulMM4ResultExistsQ)
- File file = new File(directory + "/" + name + ".mm4out");
- int failureFlag = 1;// flag (1 or 0) indicating whether the MM4 job failed
- int failureOverrideFlag = 0;// flag (1 or 0) to override success as measured by failureFlag
- if (file.exists()) {// if the file exists, do further checks; otherwise, we will skip to final statement and
- // return false
- try {
- FileReader in = new FileReader(file);
- BufferedReader reader = new BufferedReader(in);
- String line = reader.readLine();
- while (line != null) {
- String trimLine = line.trim();
- if (trimLine.equals("STATISTICAL THERMODYNAMICS ANALYSIS")) {
- failureFlag = 0;
- } else if (trimLine.endsWith("imaginary frequencies,")) {// read the number of imaginary frequencies
- // and make sure it is zero
- String[] split = trimLine.split("\\s+");
- if (Integer.parseInt(split[3]) > 0) {
- Logger.info("*****Imaginary freqencies found:");
- failureOverrideFlag = 1;
- }
- } else if (trimLine.contains(" 0.0 (fir )")) {
- if (useCanTherm) {// zero frequencies are only acceptable when CanTherm is used
- Logger.info("*****Warning: zero freqencies found (values lower than 7.7 cm^-1 are rounded to zero in MM4 output); CanTherm should hopefully correct this:");
- } else {
- Logger.info("*****Zero freqencies found:");
- failureOverrideFlag = 1;
- }
- }
- line = reader.readLine();
- }
- reader.close();
- in.close();
- } catch (Exception e) {
- String err = "Error in reading MM4 output file \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- }
- // if the failure flag is still 0, the process should have been successful
- if (failureOverrideFlag == 1)
- failureFlag = 1; // job will be considered a failure if there are imaginary frequencies or if job terminates
- // to to excess time/cycles
- if (failureFlag == 0 && connectivityCheck > 0) {// if the process was successful, check the connectivity (if
- // requested by the user)
- boolean connectivityMatch = connectivityMatchInMM4ResultQ(name,
- directory, InChIaug);
- if (!connectivityMatch)
- if (connectivityCheck > 1) {// connectivityCheck=2
- failureFlag = 1;
- } else {// connectivityCheck=1; print warning and procede with "return true"
- Logger.warning("Connectivity in MM4 result for "
- + name
- + " ("
- + InChIaug
- + ") does not appear to match desired connectivity.");
- }
- }
- // if the failure flag is 0 and there are no negative frequencies (and connectivity matches, if requested by the
- // user), the process should have been successful
- if (failureFlag == 0)
- successFlag = 1;
- return successFlag;
- }
- // name and directory are the name and directory for the input (and output) file;
- // input script is assumed to be preexisting and have the .comi suffix where i is a number between 1 and rotors
- public void runMM4Rotor(String name, String directory, int rotors) {
- for (int j = 1; j <= rotors; j++) {
- try {
- File runningDirectory = new File(qmfolder);
- String command = name + ".com" + j;
- // File script = new File(qmfolder+command);
- // command = "./"+command;
- // script.setExecutable(true);
- command = "csh " + command;
- Process mm4Proc = Runtime.getRuntime().exec(command, null,
- runningDirectory);
- // check for errors and display the error if there is one
- // InputStream is = mm4Proc.getErrorStream();
- // InputStreamReader isr = new InputStreamReader(is);
- // BufferedReader br = new BufferedReader(isr);
- // String line=null;
- // while ( (line = br.readLine()) != null) {
- // line = line.trim();
- // if(!line.equals("STOP statement executed")){//string listed here seems to be typical
- // Logger.error(line);
- // flag=1;
- // }
- // }
- InputStream is = mm4Proc.getInputStream();
- InputStreamReader isr = new InputStreamReader(is);
- BufferedReader br = new BufferedReader(isr);
- String line = null;
- while ((line = br.readLine()) != null) {
- // do nothing
- }
- // if there was an error, indicate that an error was obtained
- // if(flag==1){
- // Logger.info("MM4 process received error (see above) on " + name);
- // }
- int exitValue = mm4Proc.waitFor();
- mm4Proc.getErrorStream().close();
- mm4Proc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- String err = "Error in running MM4 rotor process \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- }
- return;
- }
- // name and directory are the name and directory for the input (and output) file;
- // input is assumed to be preexisting and have the .mop suffix
- // returns an integer indicating success or failure of the MOPAC calculation: 1 for success, 0 for failure;
- // this function is based on the Gaussian analogue
- public int runMOPAC(String name, String directory, String InChIaug) {
- if (!new File(System.getenv("MOPAC_LICENSE"), "MOPAC2009.exe").exists()) {
- Logger.error("Please set your MOPAC_LICENSE environment variable to the directory containing MOPAC2009.exe");
- System.exit(1);
- }
- int flag = 0;
- int successFlag = 0;
- try {
- String command = new File(System.getenv("MOPAC_LICENSE"),
- "MOPAC2009.exe").getAbsolutePath() + " ";
- 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 = command.concat("\"" + directory + "/" + name
- + ".mop\" ");// specify the input file; space is important
- command = command.concat("\"" + directory + "/" + name
- + ".out\"");// specify the output file
- } else {// non-Windows case
- command = command.concat(directory + "/" + name + ".mop ");// specify the input file; space is important
- command = command.concat(directory + "/" + name + ".out");// specify the output file
- }
- Process mopacProc = Runtime.getRuntime().exec(command);
- // check for errors and display the error if there is one
- InputStream is = mopacProc.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 that an error was obtained
- if (flag == 1) {
- Logger.info("MOPAC process received error (see above) on "
- + name);
- }
- int exitValue = mopacProc.waitFor();
- mopacProc.getInputStream().close();
- mopacProc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- String err = "Error in running MOPAC process \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // look in the output file to check for the successful termination of the calculation (this is a trimmed down
- // version of what appears in successfulMOPACResultExistsQ (it doesn't have the InChI check)
- File file = new File(directory + "/" + name + ".out");
- int failureFlag = 1;// flag (1 or 0) indicating whether the MOPAC job failed
- int failureOverrideFlag = 0;// flag (1 or 0) to override success as measured by failureFlag
- if (file.exists()) {// if the file exists, do further checks; otherwise, we will skip to final statement and
- // return false
- try {
- FileReader in = new FileReader(file);
- BufferedReader reader = new BufferedReader(in);
- String line = reader.readLine();
- while (line != null) {
- String trimLine = line.trim();
- if (trimLine.equals("DESCRIPTION OF VIBRATIONS")) {// check for this line; if it is here, check for
- // negative frequencies
- // if(!MopacFileContainsNegativeFreqsQ(name, directory)) failureFlag=0;
- failureFlag = 0;
- }
- // negative frequencies notice example:
- // NOTE: SYSTEM IS NOT A GROUND STATE, THEREFORE ZERO POINT
- // ENERGY IS NOT MEANINGFULL. ZERO POINT ENERGY PRINTED
- // DOES NOT INCLUDE THE 2 IMAGINARY FREQUENCIES
- else if (trimLine.endsWith("IMAGINARY FREQUENCIES")) {
- Logger.info("*****Imaginary freqencies found:");
- failureOverrideFlag = 1;
- } else if (trimLine
- .equals("EXCESS NUMBER OF OPTIMIZATION CYCLES")) {// exceeding max cycles error
- failureOverrideFlag = 1;
- } else if (trimLine
- .equals("NOT ENOUGH TIME FOR ANOTHER CYCLE")) {// timeout error
- failureOverrideFlag = 1;
- }
- line = reader.readLine();
- }
- reader.close();
- in.close();
- } catch (Exception e) {
- String err = "Error in reading MOPAC output file \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- }
- if (failureOverrideFlag == 1)
- failureFlag = 1; // job will be considered a failure if there are imaginary frequencies or if job terminates
- // to to excess time/cycles
- if (failureFlag == 0 && connectivityCheck > 0) {// if the process was successful, check the connectivity (if
- // requested by the user)
- boolean connectivityMatch = connectivityMatchInPM3ResultQ(name,
- directory, InChIaug, ".out", "moo");
- if (!connectivityMatch)
- if (connectivityCheck > 1) {// connectivityCheck=2
- failureFlag = 1;
- } else {// connectivityCheck=1; print warning and procede with "return true"
- Logger.warning("Connectivity in quantum result for "
- + name
- + " ("
- + InChIaug
- + ") does not appear to match desired connectivity.");
- }
- }
- // if the failure flag is 0 (completed, there are no negative frequencies, and connectivity matches (if
- // requested by user)), the process should have been successful
- if (failureFlag == 0)
- successFlag = 1;
- return successFlag;
- }
- public String getGaussianPM3ParseCommand(String name, String directory) {
- String command = null;
- 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/GaussianPM3ParsingScript.py\" ";
- String logfilepath = "\"" + directory + "/" + name + ".log\"";
- command = command.concat(logfilepath);
- command = command
- .concat(" \"" + System.getenv("RMG") + "/source\"");// this will pass $RMG/source to the script (in
- // order to get the appropriate path for importing
- } else {// non-Windows case
- command = "python " + System.getProperty("RMG.workingDirectory")
- + "/scripts/GaussianPM3ParsingScript.py ";
- String logfilepath = directory + "/" + name + ".log";
- command = command.concat(logfilepath);
- command = command.concat(" " + System.getenv("RMG") + "/source");// this will pass $RMG/source to the script
- // (in order to get the appropriate path for importing
- }
- return command;
- }
- // parse the results using cclib and return a ThermoData object; name and directory indicate the location of the
- // Gaussian .log file
- public ThermoData parseGaussianPM3(String name, String directory,
- ChemGraph p_chemGraph) {
- // //parse the Gaussian file using cclib
- String command = getGaussianPM3ParseCommand(name, directory);
- ThermoData result = getPM3MM4ThermoDataUsingCCLib(name, directory,
- p_chemGraph, command);
- result.setSource("Gaussian PM3 calculation");
- Logger.info("Thermo for " + name + ": " + result.toString());// print result, at least for debugging purposes
- return result;
- }
- public String getMM4ParseCommand(String name, String directory) {
- String command = "python " + System.getProperty("RMG.workingDirectory")
- + "/scripts/MM4ParsingScript.py ";
- String logfilepath = directory + "/" + name + ".mm4out";
- command = command.concat(logfilepath);
- command = command.concat(" " + System.getenv("RMG") + "/source");// this will pass $RMG/source to the script (in
- // order to get the appropriate path for importing
- return command;
- }
- // parse the results using cclib and return a ThermoData object; name and directory indicate the location of the MM4
- // .mm4out file
- public ThermoData parseMM4(String name, String directory,
- ChemGraph p_chemGraph) {
- String command = getMM4ParseCommand(name, directory);
- ThermoData result = getPM3MM4ThermoDataUsingCCLib(name, directory,
- p_chemGraph, command);
- result.setSource("MM4 calculation");
- Logger.info("Thermo for " + name + ": " + result.toString());// print result, at least for debugging purposes
- return result;
- }
- // parse the results using cclib and CanTherm and return a ThermoData object; name and directory indicate the
- // location of the MM4 .mm4out file
- // formerly known as parseMM4withForceMat
- public QMData performCanThermCalcs(String name, String directory,
- ChemGraph p_chemGraph, double[] dihedralMinima, boolean forceRRHO) {
- // 1. parse the MM4 file with cclib to get atomic number vector and geometry
- QMData qmdata = getMM4QMDataWithCClib(name, directory, true);
- // unpack the needed results
- double energy = qmdata.energy;
- double stericEnergy = qmdata.stericEnergy;
- ArrayList freqs = qmdata.freqs;
- // 2. compute H0/E0; note that we will compute H0 for CanTherm by H0=Hf298(harmonicMM4)-(H298-H0)harmonicMM4,
- // where harmonicMM4 values come from cclib parsing; 298.16 K is the standard temperature used by MM4; also note that
- // Hthermal should include the R*T contribution (R*298.16 (enthalpy vs. energy difference)) and H0=E0 (T=0)
- double T_MM4 = 298.16;
- energy *= Hartree_to_kcal;// convert from Hartree to kcal/mol
- stericEnergy *= Hartree_to_kcal;// convert from Hartree to kcal/mol
- double Hthermal = 5. / 2. * R * T_MM4 / 1000.;// enthalpy vs. energy difference(RT)+translation(3/2*R*T)
- // contribution to thermal energy
- // rotational contribution
- if (p_chemGraph.getAtomNumber() == 1)
- Hthermal += 0.0;
- else if (p_chemGraph.isLinear())
- Hthermal += R * T_MM4 / 1000.;
- else
- Hthermal += 3. / 2. * R * T_MM4 / 1000;
- // vibrational contribution
- if (p_chemGraph.getAtomNumber() != 1)
- Hthermal += R * calcVibH(freqs, T_MM4, h, k, c) / 1000.;
- energy = energy - Hthermal;
- // 3. write CanTherm input file
- // determine point group using the SYMMETRY Program
- String geom = qmdata.getSYMMETRYinput();
- String pointGroup = determinePointGroupUsingSYMMETRYProgram(geom, name);
- double sigmaCorr = getSigmaCorr(pointGroup);
- String canInp = "Calculation: Thermo\n";
- canInp += "Trange: 300 100 13\n";// temperatures from 300 to 1500 in increments of 100
- canInp += "Scale: 1.0\n";// scale factor of 1
- canInp += "Mol 1\n";
- if (p_chemGraph.getAtomNumber() == 1)
- canInp += "ATOM\n";
- else if (p_chemGraph.isLinear())
- canInp += "LINEAR\n";
- else
- canInp += "NONLINEAR\n";
- canInp += "GEOM MM4File " + name + ".mm4out\n";// geometry file; ***special MM4 treatment in CanTherm; another
- // option would be to use mm4opt file, but CanTherm code would need to be modified accordingly
- canInp += "FORCEC MM4File " + name + ".fmat\n";// force constant file; ***special MM4 treatment in CanTherm
- if (forceRRHO)
- name = name + "_RRHO"; // "_RRHO" will be appended to InChIKey for RRHO calcs (though we want to use
- // unmodified name in getQMDataWithCClib above, and in GEOM and FORCEC sections
- canInp += "ENERGY " + energy + " MM4\n";// ***special MM4 treatment in CanTherm
- canInp += "EXTSYM " + Math.exp(-sigmaCorr) + "\n";// ***modified treatment in CanTherm; traditional EXTSYM
- // integer replaced by EXTSYM double, to allow fractional values that take chirality into account
- canInp += "NELEC 1\n";// multiplicity = 1; all cases we consider will be non-radicals
- int rotors = p_chemGraph.getInternalRotor();
- String rotInput = null;
- if (!useHindRot || rotors == 0 || forceRRHO)
- canInp += "ROTORS 0\n";// do not consider hindered rotors
- else {
- int rotorCount = 0;
- canInp += "ROTORS " + rotors + " " + name + ".rotinfo\n";
- canInp += "POTENTIAL separable mm4files_inertia\n";// ***special MM4 treatment in canTherm;
- String rotNumbersLine = "" + stericEnergy;// the line that will contain V0 (kcal/mol), and all the dihedral
- // minima (degrees)
- rotInput = "L1: 1 2 3\n";
- LinkedHashMap rotorInfo = p_chemGraph.getInternalRotorInformation();
- Iterator iter = rotorInfo.keySet().iterator();
- while (iter.hasNext()) {
- rotorCount++;
- int[] rotorAtoms = (int[]) iter.next();
- LinkedHashSet rotatingGroup = (LinkedHashSet) rotorInfo
- .get(rotorAtoms);
- Iterator iterGroup = rotatingGroup.iterator();
- rotInput += "L2: " + rotorAtoms[4] + " " + rotorAtoms[1] + " "
- + rotorAtoms[2];// note: rotorAtoms[4] is the rotorSymmetry number as estimated by
- // calculateRotorSymmetryNumber; this will be taken into account elsewhere
- while (iterGroup.hasNext()) {// print the atoms associated with rotorAtom 2
- Integer id = (Integer) iterGroup.next();
- if (id != rotorAtoms[2])
- rotInput += " " + id;// don't print atom2
- }
- rotInput += "\n";
- canInp += name + ".mm4rotopt" + rotorCount + " ";// potential files will be named as name.mm4rotopti
- rotNumbersLine += " " + dihedralMinima[rotorCount - 1];
- }
- canInp += "\n" + rotNumbersLine + "\n";
- }
- canInp += "0\n";// no bond-additivity corrections
- try {
- File canFile = new File(directory + "/" + name + ".can");
- FileWriter fw = new FileWriter(canFile);
- fw.write(canInp);
- fw.close();
- if (rotInput != null) {// write the rotor information
- File rotFile = new File(directory + "/" + name + ".rotinfo");
- FileWriter fwr = new FileWriter(rotFile);
- fwr.write(rotInput);
- fwr.close();
- }
- } catch (Exception e) {
- String err = "Error in writing CanTherm input \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // 4 call CanTherm
- try {
- File runningDirectory = new File(qmfolder);
- String canCommand = "python " + System.getenv("RMG")
- + "/source/CanTherm/source/CanTherm.py " + name + ".can";
- Process canProc = Runtime.getRuntime().exec(canCommand, null,
- runningDirectory);
- InputStream is = canProc.getInputStream();
- InputStreamReader isr = new InputStreamReader(is);
- BufferedReader br = new BufferedReader(isr);
- String line = null;
- while ((line = br.readLine()) != null) {
- // do nothing
- }
- int exitValue = canProc.waitFor();
- canProc.getErrorStream().close();
- canProc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- String err = "Error in running CanTherm process \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- return qmdata;
- }
- public ThermoData parseCanThermFile(String name, String directory,
- ChemGraph p_chemGraph) {
- // 5. read CanTherm output
- Double Hf298 = null;
- Double S298 = null;
- Double Cp300 = null;
- Double Cp400 = null;
- Double Cp500 = null;
- Double Cp600 = null;
- Double Cp800 = null;
- Double Cp1000 = null;
- Double Cp1500 = null;
- File file = new File(directory + "/" + name + ".canout");
- try {
- FileReader in = new FileReader(file);
- BufferedReader reader = new BufferedReader(in);
- String line = reader.readLine();
- while (!line.startsWith("Hf298 S298 Cps:")) {// get to the end of the file with the data we want
- line = reader.readLine();
- }
- String[] split = reader.readLine().trim().split("\\s+");// read the next line, which contains the
- // information we want
- Hf298 = Double.parseDouble(split[0]);
- S298 = Double.parseDouble(split[1]);
- Cp300 = Double.parseDouble(split[2]);
- Cp400 = Double.parseDouble(split[3]);
- Cp500 = Double.parseDouble(split[4]);
- Cp600 = Double.parseDouble(split[5]);
- Cp800 = Double.parseDouble(split[7]);
- Cp1000 = Double.parseDouble(split[9]);
- Cp1500 = Double.parseDouble(split[14]);
- reader.close();
- in.close();
- } catch (Exception e) {
- String err = "Error in reading CanTherm .canout file \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- ThermoData result = new ThermoData(Hf298, S298, Cp300, Cp400, Cp500,
- Cp600, Cp800, Cp1000, Cp1500, 3, 1, 1,
- "MM4 calculation; includes CanTherm analysis of force-constant matrix");// this includes rough estimates
- // of uncertainty
- result.setSource("MM4 calculation with CanTherm analysis");
- Logger.info("Thermo for " + name + ": " + result.toString());// print result, at least for debugging purposes
- return result;
- }
- public String getMopacPM3ParseCommand(String name, String directory) {
- String command = null;
- 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/MopacPM3ParsingScript.py\" ";
- String logfilepath = "\"" + directory + "/" + name + ".out\"";
- command = command.concat(logfilepath);
- command = command
- .concat(" \"" + System.getenv("RMG") + "/source\"");// this will pass $RMG/source to the script (in
- // order to get the appropriate path for importing
- } else {// non-Windows case
- command = "python " + System.getProperty("RMG.workingDirectory")
- + "/scripts/MopacPM3ParsingScript.py ";
- String logfilepath = directory + "/" + name + ".out";
- command = command.concat(logfilepath);
- command = command.concat(" " + System.getenv("RMG") + "/source");// this will pass $RMG/source to the script
- // (in order to get the appropriate path for importing
- }
- return command;
- }
- // parse the results using cclib and return a ThermoData object; name and directory indicate the location of the
- // MOPAC .out file
- public ThermoData parseMopac(String name, String directory,
- ChemGraph p_chemGraph, String qmMethod) {
- String command = getMopacPM3ParseCommand(name, directory);
- ThermoData result = getPM3MM4ThermoDataUsingCCLib(name, directory,
- p_chemGraph, command);
- result.setSource("MOPAC "+qmMethod+" calculation");
- Logger.info("Thermo for " + name + ": " + result.toString());// print result, at least for debugging purposes
- return result;
- }
- // separated from parseMopacPM3, since the function was originally based off of parseGaussianPM3 and was very
- // similar (differences being command and logfilepath variables);
- public ThermoData getPM3MM4ThermoDataUsingCCLib(String name,
- String directory, ChemGraph p_chemGraph, String command) {
- // parse the file using cclib
- int gdStateDegen = p_chemGraph.getRadicalNumber() + 1;// calculate ground state degeneracy from the number of
- // radicals; this should give the same result as spin multiplicity in Gaussian input file (and output file), but we do
- // not explicitly check this (we could use "mult" which cclib reads in if we wanted to do so); also, note that this is
- // not always correct, as there can apparently be additional spatial degeneracy for non-symmetric linear molecules like
- // OH radical (cf. http://cccbdb.nist.gov/thermo.asp)
- QMData qmdata = getQMDataWithCClib(name, directory, command, false);
- ThermoData result = calculateThermoFromPM3MM4Calc(qmdata, gdStateDegen,
- name);
- return result;
- }
- // returns a thermo result, given results from quantum PM3 calculation or MM4 calculation (originally, this was in
- // parseGaussianPM3 function
- public ThermoData calculateThermoFromPM3MM4Calc(QMData qmdata,
- int gdStateDegen, String name) {
- // unpack qmdata
- int natoms = qmdata.natoms;
- ArrayList atomicNumber = qmdata.atomicNumber; // vector of atomic numbers (integers) (apparently Vector is
- // thread-safe; cf. http://answers.yahoo.com/question/index?qid=20081214065127AArZDT3; ...should I be using this
- // instead?)
- ArrayList x_coor = qmdata.x_coor; // vectors of x-, y-, and z-coordinates (doubles) (Angstroms) (in order
- // corresponding to above atomic numbers)
- ArrayList y_coor = qmdata.y_coor;
- ArrayList z_coor = qmdata.z_coor;
- double energy = qmdata.energy; // PM3 energy (Hf298) in Hartree (***note: in the case of MOPAC, the MOPAC file
- // will contain in units of kcal/mol, but modified ccLib will return in Hartree)
- double molmass = qmdata.molmass; // molecular mass in amu
- ArrayList freqs = qmdata.freqs; // list of frequencies in units of cm^-1
- double rotCons_1 = qmdata.rotCons_1;// rotational constants in (1/s)
- double rotCons_2 = qmdata.rotCons_2;
- double rotCons_3 = qmdata.rotCons_3;
- // determine point group using the SYMMETRY Program
- String geom = natoms + "\n";
- for (int i = 0; i < natoms; i++) {
- geom += atomicNumber.get(i) + " " + x_coor.get(i) + " "
- + y_coor.get(i) + " " + z_coor.get(i) + "\n";
- }
- // String pointGroup = determinePointGroupUsingSYMMETRYProgram(geom, 0.01);
- String pointGroup = determinePointGroupUsingSYMMETRYProgram(geom, name);
- // calculate thermo quantities using stat. mech. equations
- // boolean linearity = p_chemGraph.isLinear();//determine linearity (perhaps it would be more appropriate to
- // determine this from point group?)
- boolean linearity = false;
- if (pointGroup.equals("Cinfv") || pointGroup.equals("Dinfh"))
- linearity = true;// determine linearity from 3D-geometry; changed to correctly consider linear ketene
- // radical case
- // we will use number of atoms from above (alternatively, we could use the chemGraph); this is needed to test
- // whether the species is monoatomic
- double Hf298, S298, Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500;
- double sigmaCorr = getSigmaCorr(pointGroup);
- Hf298 = energy * Hartree_to_kcal;
- S298 = R
- * Math.log(gdStateDegen)
- + R
- * (3.
- / 2.
- * Math.log(2. * Math.PI * molmass
- / (1000. * Na * Math.pow(h, 2.))) + 5. / 2.
- * Math.log(k * 298.15) - Math.log(100000.) + 5. / 2.);// electronic + translation; note use of
- // 10^5 Pa for standard pressure; also note that molecular mass needs to be divided by 1000 for kg units
- Cp300 = 5. / 2. * R;
- Cp400 = 5. / 2. * R;
- Cp500 = 5. / 2. * R;
- Cp600 = 5. / 2. * R;
- Cp800 = 5. / 2. * R;
- Cp1000 = 5. / 2. * R;
- Cp1500 = 5. / 2. * R;
- if (natoms > 1) {// include statistical correction and rotational (without symmetry number, vibrational
- // contributions if species is polyatomic
- if (linearity) {// linear case
- // determine the rotational constant (note that one of the rotcons will be zero)
- double rotCons;
- if (rotCons_1 > 0.0001)
- rotCons = rotCons_1;
- else
- rotCons = rotCons_2;
- S298 += R * sigmaCorr + R
- * (Math.log(k * 298.15 / (h * rotCons)) + 1) + R
- * calcVibS(freqs, 298.15, h, k, c);
- Cp300 += R + R * calcVibCp(freqs, 300., h, k, c);
- Cp400 += R + R * calcVibCp(freqs, 400., h, k, c);
- Cp500 += R + R * calcVibCp(freqs, 500., h, k, c);
- Cp600 += R + R * calcVibCp(freqs, 600., h, k, c);
- Cp800 += R + R * calcVibCp(freqs, 800., h, k, c);
- Cp1000 += R + R * calcVibCp(freqs, 1000., h, k, c);
- Cp1500 += R + R * calcVibCp(freqs, 1500., h, k, c);
- } else {// nonlinear case
- S298 += R
- * sigmaCorr
- + R
- * (3.
- / 2.
- * Math.log(k * 298.15 / h)
- - 1.
- / 2.
- * Math.log(rotCons_1 * rotCons_2 * rotCons_3
- / Math.PI) + 3. / 2.) + R
- * calcVibS(freqs, 298.15, h, k, c);
- Cp300 += 3. / 2. * R + R * calcVibCp(freqs, 300., h, k, c);
- Cp400 += 3. / 2. * R + R * calcVibCp(freqs, 400., h, k, c);
- Cp500 += 3. / 2. * R + R * calcVibCp(freqs, 500., h, k, c);
- Cp600 += 3. / 2. * R + R * calcVibCp(freqs, 600., h, k, c);
- Cp800 += 3. / 2. * R + R * calcVibCp(freqs, 800., h, k, c);
- Cp1000 += 3. / 2. * R + R * calcVibCp(freqs, 1000., h, k, c);
- Cp1500 += 3. / 2. * R + R * calcVibCp(freqs, 1500., h, k, c);
- }
- }
- ThermoData result = new ThermoData(Hf298, S298, Cp300, Cp400, Cp500,
- Cp600, Cp800, Cp1000, Cp1500, 5, 1, 1, "PM3 or MM4 calculation");// this includes rough estimates of
- // uncertainty
- return result;
- }
- // gets the statistical correction for S in dimensionless units (divided by R)
- public double getSigmaCorr(String pointGroup) {
- double sigmaCorr = 0;
- // determine statistical correction factor for 1. external rotational symmetry (affects rotational partition
- // function) and 2. chirality (will add R*ln2 to entropy) based on point group
- // ref: http://cccbdb.nist.gov/thermo.asp
- // assumptions below for Sn, T, Th, O, I seem to be in line with expectations based on order reported at:
- // http://en.wikipedia.org/w/index.php?title=List_of_character_tables_for_chemically_important_3D_point_groups&oldid=287261611
- // (assuming order = symmetry number * 2 (/2 if chiral))...this appears to be true for all point groups I "know" to be
- // correct
- // minor concern: does SYMMETRY appropriately calculate all Sn groups considering 2007 discovery of previous
- // errors in character tables (cf. Wikipedia article above)
- if (pointGroup.equals("C1"))
- sigmaCorr = +Math.log(2.);// rot. sym. = 1, chiral
- else if (pointGroup.equals("Cs"))
- sigmaCorr = 0; // rot. sym. = 1
- else if (pointGroup.equals("Ci"))
- sigmaCorr = 0; // rot. sym. = 1
- else if (pointGroup.equals("C2"))
- sigmaCorr = 0;// rot. sym. = 2, chiral (corrections cancel)
- else if (pointGroup.equals("C3"))
- sigmaCorr = +Math.log(2.) - Math.log(3.);// rot. sym. = 3, chiral
- else if (pointGroup.equals("C4"))
- sigmaCorr = +Math.log(2.) - Math.log(4.);// rot. sym. = 4, chiral
- else if (pointGroup.equals("C5"))
- sigmaCorr = +Math.log(2.) - Math.log(5.);// rot. sym. = 5, chiral
- else if (pointGroup.equals("C6"))
- sigmaCorr = +Math.log(2.) - Math.log(6.);// rot. sym. = 6, chiral
- else if (pointGroup.equals("C7"))
- sigmaCorr = +Math.log(2.) - Math.log(7.);// rot. sym. = 7, chiral
- else if (pointGroup.equals("C8"))
- sigmaCorr = +Math.log(2.) - Math.log(8.);// rot. sym. = 8, chiral
- else if (pointGroup.equals("D2"))
- sigmaCorr = +Math.log(2.) - Math.log(4.);// rot. sym. = 4, chiral
- else if (pointGroup.equals("D3"))
- sigmaCorr = +Math.log(2.) - Math.log(6.);// rot. sym. = 6, chiral
- else if (pointGroup.equals("D4"))
- sigmaCorr = +Math.log(2.) - Math.log(8.);// rot. sym. = 8, chiral
- else if (pointGroup.equals("D5"))
- sigmaCorr = +Math.log(2.) - Math.log(10.);// rot. sym. = 10, chiral
- else if (pointGroup.equals("D6"))
- sigmaCorr = +Math.log(2.) - Math.log(12.);// rot. sym. = 12, chiral
- else if (pointGroup.equals("D7"))
- sigmaCorr = +Math.log(2.) - Math.log(14.);// rot. sym. = 14, chiral
- else if (pointGroup.equals("D8"))
- sigmaCorr = +Math.log(2.) - Math.log(16.);// rot. sym. = 16, chiral
- else if (pointGroup.equals("C2v"))
- sigmaCorr = -Math.log(2.);// rot. sym. = 2
- else if (pointGroup.equals("C3v"))
- sigmaCorr = -Math.log(3.);// rot. sym. = 3
- else if (pointGroup.equals("C4v"))
- sigmaCorr = -Math.log(4.);// rot. sym. = 4
- else if (pointGroup.equals("C5v"))
- sigmaCorr = -Math.log(5.);// rot. sym. = 5
- else if (pointGroup.equals("C6v"))
- sigmaCorr = -Math.log(6.);// rot. sym. = 6
- else if (pointGroup.equals("C7v"))
- sigmaCorr = -Math.log(7.);// rot. sym. = 7
- else if (pointGroup.equals("C8v"))
- sigmaCorr = -Math.log(8.);// rot. sym. = 8
- else if (pointGroup.equals("C2h"))
- sigmaCorr = -Math.log(2.);// rot. sym. = 2
- else if (pointGroup.equals("C3h"))
- sigmaCorr = -Math.log(3.);// rot. sym. = 3
- else if (pointGroup.equals("C4h"))
- sigmaCorr = -Math.log(4.);// rot. sym. = 4
- else if (pointGroup.equals("C5h"))
- sigmaCorr = -Math.log(5.);// rot. sym. = 5
- else if (pointGroup.equals("C6h"))
- sigmaCorr = -Math.log(6.);// rot. sym. = 6
- else if (pointGroup.equals("C7h"))
- sigmaCorr = -Math.log(7.);// rot. sym. = 7
- else if (pointGroup.equals("C8h"))
- sigmaCorr = -Math.log(8.);// rot. sym. = 8
- else if (pointGroup.equals("D2h"))
- sigmaCorr = -Math.log(4.);// rot. sym. = 4
- else if (pointGroup.equals("D3h"))
- sigmaCorr = -Math.log(6.);// rot. sym. = 6
- else if (pointGroup.equals("D4h"))
- sigmaCorr = -Math.log(8.);// rot. sym. = 8
- else if (pointGroup.equals("D5h"))
- sigmaCorr = -Math.log(10.);// rot. sym. = 10
- else if (pointGroup.equals("D6h"))
- sigmaCorr = -Math.log(12.);// rot. sym. = 12
- else if (pointGroup.equals("D7h"))
- sigmaCorr = -Math.log(14.);// rot. sym. = 14
- else if (pointGroup.equals("D8h"))
- sigmaCorr = -Math.log(16.);// rot. sym. = 16
- else if (pointGroup.equals("D2d"))
- sigmaCorr = -Math.log(4.);// rot. sym. = 4
- else if (pointGroup.equals("D3d"))
- sigmaCorr = -Math.log(6.);// rot. sym. = 6
- else if (pointGroup.equals("D4d"))
- sigmaCorr = -Math.log(8.);// rot. sym. = 8
- else if (pointGroup.equals("D5d"))
- sigmaCorr = -Math.log(10.);// rot. sym. = 10
- else if (pointGroup.equals("D6d"))
- sigmaCorr = -Math.log(12.);// rot. sym. = 12
- else if (pointGroup.equals("D7d"))
- sigmaCorr = -Math.log(14.);// rot. sym. = 14
- else if (pointGroup.equals("D8d"))
- sigmaCorr = -Math.log(16.);// rot. sym. = 16
- else if (pointGroup.equals("S4"))
- sigmaCorr = -Math.log(2.);// rot. sym. = 2 ;*** assumed achiral
- else if (pointGroup.equals("S6"))
- sigmaCorr = -Math.log(3.);// rot. sym. = 3 ;*** assumed achiral
- else if (pointGroup.equals("S8"))
- sigmaCorr = -Math.log(4.);// rot. sym. = 4 ;*** assumed achiral
- else if (pointGroup.equals("T"))
- sigmaCorr = +Math.log(2.) - Math.log(12.);// rot. sym. = 12, *** assumed chiral
- else if (pointGroup.equals("Th"))
- sigmaCorr = -Math.log(12.);// ***assumed rot. sym. = 12
- else if (pointGroup.equals("Td"))
- sigmaCorr = -Math.log(12.);// rot. sym. = 12
- else if (pointGroup.equals("O"))
- sigmaCorr = +Math.log(2.) - Math.log(24.);// ***assumed rot. sym. = 24, chiral
- else if (pointGroup.equals("Oh"))
- sigmaCorr = -Math.log(24.);// rot. sym. = 24
- else if (pointGroup.equals("Cinfv"))
- sigmaCorr = 0;// rot. sym. = 1
- else if (pointGroup.equals("Dinfh"))
- sigmaCorr = -Math.log(2.);// rot. sym. = 2
- else if (pointGroup.equals("I"))
- sigmaCorr = +Math.log(2.) - Math.log(60.);// ***assumed rot. sym. = 60, chiral
- else if (pointGroup.equals("Ih"))
- sigmaCorr = -Math.log(60.);// rot. sym. = 60
- else if (pointGroup.equals("Kh"))
- sigmaCorr = 0;// arbitrarily set to zero...one could argue that it is infinite; apparently this is the point
- // group of a single atom (cf. http://www.cobalt.chem.ucalgary.ca/ps/symmetry/tests/G_Kh); this should not have a
- // rotational partition function, and we should not use the symmetry correction in this case
- else {// this point should not be reached, based on checks performed in determinePointGroupUsingSYMMETRYProgram
- Logger.critical("Unrecognized point group: " + pointGroup);
- System.exit(0);
- }
- return sigmaCorr;
- }
- // determine the point group using the SYMMETRY program (http://www.cobalt.chem.ucalgary.ca/ps/symmetry/)
- // required input is a line with number of atoms followed by lines for each atom including atom number and x,y,z
- // coordinates
- // finalTol determines how loose the point group criteria are; values are comparable to those specifed in the
- // GaussView point group interface
- // public String determinePointGroupUsingSYMMETRYProgram(String geom, double finalTol){
- public String determinePointGroupUsingSYMMETRYProgram(String geom,
- String name) {
- int attemptNumber = 1;
- int maxAttemptNumber = 4;
- boolean pointGroupFound = false;
- // write the input file
- File inputFile = new File(qmfolder, name + ".symm");
- try {
- FileWriter fw = new FileWriter(inputFile);
- fw.write(geom);
- fw.close();
- } catch (IOException e) {
- String err = "Error writing input file for point group calculation";
- err += e.toString();
- Logger.critical(err);
- System.exit(0);
- }
- String result = "";
- String command = "";
- while (attemptNumber <= maxAttemptNumber && !pointGroupFound) {
- // call the program and read the result
- result = "";
- String[] lineArray;
- try {
- if (System.getProperty("os.name").toLowerCase()
- .contains("windows")) {// the Windows case where the precompiled executable seems to need to be
- // called from a batch script
- if (attemptNumber == 1)
- command = "\""
- + System.getProperty("RMG.workingDirectory")
- + "/scripts/symmetryDefault2.bat\" \""
- + inputFile.getAbsolutePath() + "\"";// 12/1/09 gmagoon: switched to use slightly looser
- // criteria of 0.02 rather than 0.01 to handle methylperoxyl radical result from MOPAC
- else if (attemptNumber == 2)
- command = "\""
- + System.getProperty("RMG.workingDirectory")
- + "/scripts/symmetryLoose.bat\" \""
- + inputFile.getAbsolutePath() + "\"";// looser criteria (0.1 instead of 0.01) to
- // properly identify C2v group in VBURLMBUVWIEMQ-UHFFFAOYAVmult5 (InChI=1/C3H4O2/c1-3(2,4)5/h1-2H2/mult5) MOPAC result;
- // C2 and sigma were identified with default, but it should be C2 and sigma*2
- else if (attemptNumber == 3)
- command = "\""
- + System.getProperty("RMG.workingDirectory")
- + "/scripts/symmetryLoose2.bat\" \""
- + inputFile.getAbsolutePath() + "\"";// looser criteria to properly identify D2d group
- // in XXHDHKZTASMVSX-UHFFFAOYAM
- else if (attemptNumber == 4) {
- Logger.error("*****WARNING****: Using last-resort symmetry estimation options; symmetry may be underestimated");
- command = "\""
- + System.getProperty("RMG.workingDirectory")
- + "/scripts/symmetryLastResort.bat\" \""
- + inputFile.getAbsolutePath() + "\"";// last resort criteria to avoid crashing (this
- // will likely result in identification of C1 point group)
- } else {
- Logger.critical("Invalid attemptNumber: "
- + attemptNumber);
- System.exit(0);
- }
- } else {// in other (non-Windows) cases, where it is compiled from scratch, we should be able to run
- // this directly
- if (attemptNumber == 1)
- command = System.getProperty("RMG.workingDirectory")
- + "/bin/SYMMETRY.EXE -final 0.02 "
- + inputFile.getAbsolutePath();// 12/1/09 gmagoon: switched to use slightly looser
- // criteria of 0.02 rather than 0.01 to handle methylperoxyl radical result from MOPAC
- else if (attemptNumber == 2)
- command = System.getProperty("RMG.workingDirectory")
- + "/bin/SYMMETRY.EXE -final 0.1 "
- + inputFile.getAbsolutePath();// looser criteria (0.1 instead of 0.01) to properly
- // identify C2v group in VBURLMBUVWIEMQ-UHFFFAOYAVmult5 (InChI=1/C3H4O2/c1-3(2,4)5/h1-2H2/mult5) MOPAC result; C2 and
- // sigma were identified with default, but it should be C2 and sigma*2
- else if (attemptNumber == 3)
- command = System.getProperty("RMG.workingDirectory")
- + "/bin/SYMMETRY.EXE -primary 0.2 -final 0.1 "
- + inputFile.getAbsolutePath();// looser criteria to identify D2d group in
- // XXHDHKZTASMVSX-UHFFFAOYAM (InChI=1/C12H16/c1-5-9-10(6-2)12(8-4)11(9)7-3/h5-12H,1-4H2)
- else if (attemptNumber == 4) {
- Logger.warning("*****WARNING****: Using last-resort symmetry estimation options; symmetry may be underestimated");
- command = System.getProperty("RMG.workingDirectory")
- + "/bin/SYMMETRY.EXE -final 0.0 "
- + inputFile.getAbsolutePath();// last resort criteria to avoid crashing (this will
- // likely result in identification of C1 point group)
- } else {
- Logger.critical("Invalid attemptNumber: "
- + attemptNumber);
- System.exit(0);
- }
- }
- Process symmProc = Runtime.getRuntime().exec(command);
- // check for errors and display the error if there is one
- InputStream is = symmProc.getInputStream();
- InputStreamReader isr = new InputStreamReader(is);
- BufferedReader br = new BufferedReader(isr);
- String line = null;
- while ((line = br.readLine()) != null) {
- if (line.startsWith("It seems to be the ")) {// last line, ("It seems to be the [x] point group")
- // indicates point group
- lineArray = line.split(" ");// split the line around spaces
- result = lineArray[5];// point group string should be the 6th word
- }
- }
- int exitValue = symmProc.waitFor();
- symmProc.getErrorStream().close();
- symmProc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- String err = "Error in running point group calculation process using SYMMETRY \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // check for a recognized point group
- if (result.equals("C1") || result.equals("Cs")
- || result.equals("Ci") || result.equals("C2")
- || result.equals("C3") || result.equals("C4")
- || result.equals("C5") || result.equals("C6")
- || result.equals("C7") || result.equals("C8")
- || result.equals("D2") || result.equals("D3")
- || result.equals("D4") || result.equals("D5")
- || result.equals("D6") || result.equals("D7")
- || result.equals("D8") || result.equals("C2v")
- || result.equals("C3v") || result.equals("C4v")
- || result.equals("C5v") || result.equals("C6v")
- || result.equals("C7v") || result.equals("C8v")
- || result.equals("C2h") || result.equals("C3h")
- || result.equals("C4h") || result.equals("C5h")
- || result.equals("C6h") || result.equals("C7h")
- || result.equals("C8h") || result.equals("D2h")
- || result.equals("D3h") || result.equals("D4h")
- || result.equals("D5h") || result.equals("D6h")
- || result.equals("D7h") || result.equals("D8h")
- || result.equals("D2d") || result.equals("D3d")
- || result.equals("D4d") || result.equals("D5d")
- || result.equals("D6d") || result.equals("D7d")
- || result.equals("D8d") || result.equals("S4")
- || result.equals("S6") || result.equals("S8")
- || result.equals("T") || result.equals("Th")
- || result.equals("Td") || result.equals("O")
- || result.equals("Oh") || result.equals("Cinfv")
- || result.equals("Dinfh") || result.equals("I")
- || result.equals("Ih") || result.equals("Kh"))
- pointGroupFound = true;
- else {
- if (attemptNumber < maxAttemptNumber)
- Logger.info("Attempt number " + attemptNumber
- + " did not identify a recognized point group ("
- + result
- + "). Will retry with looser point group criteria.");
- else {
- Logger.critical("Final attempt number " + attemptNumber
- + " did not identify a recognized point group ("
- + result + "). Exiting.");
- System.exit(0);
- }
- attemptNumber++;
- }
- }
- Logger.info("Point group: " + result);// print result, at least for debugging purposes
- return result;
- }
- // gmagoon 6/8/09
- // calculate the vibrational contribution (divided by R, dimensionless) at temperature, T, in Kelvin to entropy
- // p_freqs in cm^-1; c in cm/s; k in J/K; h in J-s
- // ref.: http://cccbdb.nist.gov/thermo.asp
- public double calcVibS(ArrayList p_freqs, double p_T, double h, double k,
- double c) {
- double Scontrib = 0;
- double dr;
- for (int i = 0; i < p_freqs.size(); i++) {
- double freq = (Double) p_freqs.get(i);
- dr = h * c * freq / (k * p_T); // frequently used dimensionless ratio
- Scontrib = Scontrib - Math.log(1. - Math.exp(-dr)) + dr
- * Math.exp(-dr) / (1. - Math.exp(-dr));
- }
- return Scontrib;
- }
- // gmagoon 6/8/09
- // calculate the vibrational contribution (divided by R, dimensionless) at temperature, T, in Kelvin to heat
- // capacity, Cp
- // p_freqs in cm^-1; c in cm/s; k in J/K; h in J-s
- // ref.: http://cccbdb.nist.gov/thermo.asp
- public double calcVibCp(ArrayList p_freqs, double p_T, double h, double k,
- double c) {
- double Cpcontrib = 0;
- double dr;
- for (int i = 0; i < p_freqs.size(); i++) {
- double freq = (Double) p_freqs.get(i);
- dr = h * c * freq / (k * p_T); // frequently used dimensionless ratio
- Cpcontrib = Cpcontrib + Math.pow(dr, 2.) * Math.exp(-dr)
- / Math.pow(1. - Math.exp(-dr), 2.);
- }
- return Cpcontrib;
- }
- // gmagoon 6/23/10
- // calculate the vibrational contribution (divided by R, units of K) at temperature, T, in Kelvin to Hthermal (ZPE
- // not included)
- // p_freqs in cm^-1; c in cm/s; k in J/K; h in J-s
- // we need to ignore zero frequencies, as MM4 does, to avoid dividing by zero; based on equation for Ev in
- // http://www.gaussian.com/g_whitepap/thermo.htm; however, we ignore zero point contribution to be consistent with the
- // input that CanTherm currently takes; note, however, that this could introduce some small inaccuracies, as the
- // frequencies may be slightly different in CanTherm vs. MM4, particularly for a frequency < 7.7 cm^-1 (treated as zero
- // in MM4)
- public double calcVibH(ArrayList p_freqs, double p_T, double h, double k,
- double c) {
- double Hcontrib = 0;
- double dr;
- for (int i = 0; i < p_freqs.size(); i++) {
- double freq = (Double) p_freqs.get(i);
- if (freq > 0.0) {// ignore zero frequencies as MM4 does
- dr = h * c * freq / (k * p_T); // frequently used dimensionless ratio
- Hcontrib = Hcontrib + dr * p_T / (Math.exp(dr) - 1.);
- }
- }
- return Hcontrib;
- }
- // determine the QM filename (element 0) and augmented InChI (element 1) for a ChemGraph
- // QM filename is InChIKey appended with mult3, mult4, mult5, or mult6 for multiplicities of 3 or higher
- // augmented InChI is InChI appended with /mult3, /mult4, /mult5, or /mult6 for multiplicities of 3 or higher
- public String[] getQMFileName(ChemGraph p_chemGraph) {
- String[] result = new String[2];
- result[0] = p_chemGraph.getModifiedInChIKeyAnew();// need to generate InChI and key anew because ChemGraph may
- // have changed (in particular, adding/removing hydrogens in HBI process)
- result[1] = p_chemGraph.getModifiedInChIAnew();
- return result;
- }
- public String getInChIFromModifiedInChI(String modInChI) {
- int pos = modInChI.indexOf("/mult");
- if (pos < 0)
- return modInChI;// no mult component
- else
- return modInChI.substring(0, pos); // return the component before the mult section
- }
- // removes b, t, m, and s layers from InChI
- public String stripStereochemLayersFromInChI(String InChI) {
- // String result1=InChI.replaceFirst("/b[^/]*", "");//replaces forward slash followed by b followed by anything
- // except forward slash with an empty string
- // String result2=result1.replaceFirst("/t[^/]*", "");//replaces forward slash followed by t followed by
- // anything except forward slash with an empty string
- // String result3=result2.replaceFirst("/m[^/]*", "");//replaces forward slash followed by m followed by
- // anything except forward slash with an empty string
- // String result4=result3.replaceFirst("/s[^/]*", "");//replaces forward slash followed by s followed by
- // anything except forward slash with an empty string
- return InChI.replaceFirst("/b[^/]*", "").replaceFirst("/t[^/]*", "")
- .replaceFirst("/m[^/]*", "").replaceFirst("/s[^/]*", "");// a shorted version of the commented out code
- // above
- }
- // returns true if a Gaussian file for the given name and directory (.log suffix) exists and indicates successful
- // completion (same criteria as used after calculation runs); terminates if the InChI doesn't match the InChI in the
- // file or if there is no InChI in the file; returns false otherwise
- public boolean successfulGaussianResultExistsQ(String name,
- String directory, String InChIaug) {
- // part of the code is taken from runGaussian code above
- // look in the output file to check for the successful termination of the Gaussian calculation
- // failed jobs will contain the a line beginning with " Error termination" near the end of the file
- File file = new File(directory + "/" + name + ".log");
- if (file.exists()) {// if the file exists, do further checks; otherwise, we will skip to final statement and
- // return false
- int failureFlag = 0;// flag (1 or 0) indicating whether the Gaussian job failed
- int InChIMatch = 0;// flag (1 or 0) indicating whether the InChI in the file matches InChIaug; this can only
- // be 1 if InChIFound is also 1;
- int InChIFound = 0;// flag (1 or 0) indicating whether an InChI was found in the log file
- int InChIPartialMatch = 0;// flag (1 or 0) indicating whether the InChI in the log file is a substring of
- // the InChI in RMG's memory
- int completeFlag = 0;
- String logFileInChI = "";
- try {
- FileReader in = new FileReader(file);
- BufferedReader reader = new BufferedReader(in);
- String line = reader.readLine();
- while (line != null) {
- if (line.startsWith(" Error termination "))
- failureFlag = 1;
- else if (line.startsWith(" ******")) {// also look for imaginary frequencies
- if (line.contains("imaginary frequencies"))
- failureFlag = 1;
- } else if (line.startsWith(" InChI=")) {
- logFileInChI = line.trim();
- // continue reading lines until a line of dashes is found (in this way, we can read InChIs that
- // span multiple lines)
- line = reader.readLine();
- while (!line.startsWith(" --------")) {
- logFileInChI += line.trim();
- line = reader.readLine();
- }
- InChIFound = 1;
- if (logFileInChI.equals(InChIaug))
- InChIMatch = 1;
- else if (InChIaug.startsWith(logFileInChI))
- InChIPartialMatch = 1;
- } else if (line
- .startsWith(" Normal termination of Gaussian 03 at")) {// check for completion notice at end
- // of file
- if (reader.readLine() == null) {// the above line will occur once in the middle of most files
- // (except for monatomic species), but the following line is non-critical (e.g.
- // "Link1: Proceeding to internal job step number 2.") so reading over this line should not cause a problem
- completeFlag = 1;
- }
- }
- line = reader.readLine();
- }
- reader.close();
- in.close();
- } catch (Exception e) {
- String err = "Error in reading preexisting Gaussian log file \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // if the complete flag is still 0, the process did not complete and is a failure
- if (completeFlag == 0)
- failureFlag = 1;
- if (failureFlag == 0 && connectivityCheck > 0) {// if the process was successful, check the connectivity (if
- // requested by the user)
- boolean connectivityMatch = connectivityMatchInPM3ResultQ(name,
- directory, InChIaug, ".log", "g03");
- if (!connectivityMatch)
- if (connectivityCheck > 1) {// connectivityCheck=2
- failureFlag = 1;
- } else {// connectivityCheck=1; print warning and procede with "return true"
- Logger.warning("Connectivity in quantum result for "
- + name
- + " ("
- + InChIaug
- + ") does not appear to match desired connectivity.");
- }
- }
- // if the failure flag is still 0, the process should have been successful
- if (failureFlag == 0 && InChIMatch == 1) {
- Logger.info("Pre-existing successful quantum result for "
- + name + " (" + InChIaug
- + ") has been found. This log file will be used.");
- return true;
- } else if (InChIFound == 1 && InChIMatch == 0) {// InChIs do not match (most likely due to limited name
- // length mirrored in log file (79 characters), but possibly due to a collision)
- if (InChIPartialMatch == 1) {// case where the InChI in memory begins with the InChI in the log file; we
- // will continue and check the input file, printing a warning if there is no match
- File inputFile = new File(directory + "/" + name + ".gjf");
- if (inputFile.exists()) {// read the Gaussian inputFile
- String inputFileInChI = "";
- try {
- FileReader inI = new FileReader(inputFile);
- BufferedReader readerI = new BufferedReader(inI);
- String lineI = readerI.readLine();
- while (lineI != null) {
- if (lineI.startsWith(" InChI=")) {
- inputFileInChI = lineI.trim();
- }
- lineI = readerI.readLine();
- }
- readerI.close();
- inI.close();
- } catch (Exception e) {
- String err = "Error in reading preexisting Gaussian gjf file \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- if (inputFileInChI.equals(InChIaug)) {
- if (failureFlag == 0) {
- Logger.info("Pre-existing successful quantum result for "
- + name
- + " ("
- + InChIaug
- + ") has been found. This log file will be used. *Note that input file was read to confirm lack of InChIKey collision (InChI probably more than 79 characters)");
- return true;
- } else {// otherwise, failureFlag==1
- Logger.info("Pre-existing quantum result for "
- + name
- + " ("
- + InChIaug
- + ") has been found, but the result was apparently unsuccessful. The file will be overwritten with a new calculation. *Note that input file was read to confirm lack of InChIKey collision (InChI probably more than 79 characters)");
- return false;
- }
- } else {
- if (inputFileInChI.equals("")) {// InChI was not found in input file
- Logger.info("*****Warning: potential InChIKey collision: InChIKey(augmented) = "
- + name
- + " RMG Augmented InChI = "
- + InChIaug
- + " Log file Augmented InChI = "
- + logFileInChI
- + " . InChI could not be found in the Gaussian input file. You should manually check that the log file contains the intended species.");
- return true;
- } else {// InChI was found but doesn't match
- Logger.critical("Congratulations! You may have discovered one of the first recorded instances of an InChIKey collision: InChIKey(augmented) = "
- + name
- + " RMG Augmented InChI = "
- + InChIaug
- + " Gaussian input file Augmented InChI = "
- + inputFileInChI);
- System.exit(0);
- }
- }
- } else {
- Logger.info("*****Warning: potential InChIKey collision: InChIKey(augmented) = "
- + name
- + " RMG Augmented InChI = "
- + InChIaug
- + " Log file Augmented InChI = "
- + logFileInChI
- + " . Gaussian input file could not be found to check full InChI. You should manually check that the log file contains the intended species.");
- return true;
- }
- } else {
- Logger.critical("Congratulations! You may have discovered one of the first recorded instances of an InChIKey collision: InChIKey(augmented) = "
- + name
- + " RMG Augmented InChI = "
- + InChIaug
- + " Log file Augmented InChI = " + logFileInChI);
- System.exit(0);
- }
- } else if (InChIFound == 0) {
- Logger.critical("An InChI was not found in file: " + name
- + ".log");
- System.exit(0);
- } else if (failureFlag == 1) {// note these should cover all possible results for this block, and if the
- // file.exists block is entered, it should return from within the block and should not reach the return statement below
- Logger.info("Pre-existing quantum result for "
- + name
- + " ("
- + InChIaug
- + ") has been found, but the result was apparently unsuccessful. The file will be overwritten with a new calculation.");
- return false;
- }
- }
- // we could print a line here for cases where the file doesn't exist, but this would probably be too verbose
- return false;
- }
- // returns true if a successful result exists (either Gaussian or MOPAC)
- // public boolean [] successfulResultExistsQ(String name, String directory, String InChIaug){
- // boolean gaussianResult=successfulGaussianResultExistsQ(name, directory, InChIaug);
- // boolean mopacResult=successfulMOPACResultExistsQ(name, directory, InChIaug);
- // return (gaussianResult || mopacResult);// returns true if either a successful Gaussian or MOPAC result exists
- // }
- // returns true if a MOPAC output file for the given name and directory (.out suffix) exists and indicates
- // successful completion (same criteria as used after calculation runs); terminates if the InChI doesn't match the InChI
- // in the file or if there is no InChI in the file; returns false otherwise
- public boolean successfulMopacResultExistsQ(String name, String directory,
- String InChIaug) {
- // part of the code is taken from analogous code for Gaussian
- // look in the output file to check for the successful termination of the calculation (assumed to be successful
- // if "description of vibrations appears)
- File file = new File(directory + "/" + name + ".out");
- if (file.exists()) {// if the file exists, do further checks; otherwise, we will skip to final statement and
- // return false
- int failureFlag = 1;// flag (1 or 0) indicating whether the MOPAC job failed
- int failureOverrideFlag = 0;// flag (1 or 0) to override success as measured by failureFlag
- int InChIMatch = 0;// flag (1 or 0) indicating whether the InChI in the file matches InChIaug; this can only
- // be 1 if InChIFound is also 1;
- int InChIFound = 0;// flag (1 or 0) indicating whether an InChI was found in the log file
- int InChIPartialMatch = 0;// flag (1 or 0) indicating whether the InChI in the log file is a substring of
- // the InChI in RMG's memory
- String logFileInChI = "";
- try {
- FileReader in = new FileReader(file);
- BufferedReader reader = new BufferedReader(in);
- String line = reader.readLine();
- while (line != null) {
- String trimLine = line.trim();
- if (trimLine.equals("DESCRIPTION OF VIBRATIONS")) {
- // if(!MopacFileContainsNegativeFreqsQ(name, directory)) failureFlag=0;//check for this line; if
- // it is here, check for negative frequencies
- failureFlag = 0;
- }
- // negative frequencies notice example:
- // NOTE: SYSTEM IS NOT A GROUND STATE, THEREFORE ZERO POINT
- // ENERGY IS NOT MEANINGFULL. ZERO POINT ENERGY PRINTED
- // DOES NOT INCLUDE THE 2 IMAGINARY FREQUENCIES
- else if (trimLine.endsWith("IMAGINARY FREQUENCIES")) {
- // Logger.info("*****Imaginary freqencies found:");
- failureOverrideFlag = 1;
- } else if (trimLine
- .equals("EXCESS NUMBER OF OPTIMIZATION CYCLES")) {// exceeding max cycles error
- failureOverrideFlag = 1;
- } else if (trimLine
- .equals("NOT ENOUGH TIME FOR ANOTHER CYCLE")) {// timeout error
- failureOverrideFlag = 1;
- } else if (line.startsWith(" InChI=")) {
- logFileInChI = line.trim();// output files should take up to 240 characters of the name in the
- // input file
- InChIFound = 1;
- if (logFileInChI.equals(InChIaug))
- InChIMatch = 1;
- else if (InChIaug.startsWith(logFileInChI))
- InChIPartialMatch = 1;
- }
- line = reader.readLine();
- }
- reader.close();
- in.close();
- } catch (Exception e) {
- String err = "Error in reading preexisting MOPAC output file \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- if (failureOverrideFlag == 1)
- failureFlag = 1; // job will be considered a failure if there are imaginary frequencies or if job
- // terminates to to excess time/cycles
- if (failureFlag == 0 && connectivityCheck > 0) {// if the process was successful, check the connectivity (if
- // requested by the user)
- boolean connectivityMatch = connectivityMatchInPM3ResultQ(name,
- directory, InChIaug, ".out", "moo");
- if (!connectivityMatch)
- if (connectivityCheck > 1) {// connectivityCheck=2
- failureFlag = 1;
- } else {// connectivityCheck=1; print warning and procede with "return true"
- Logger.warning("Connectivity in quantum result for "
- + name
- + " ("
- + InChIaug
- + ") does not appear to match desired connectivity.");
- }
- }
- // if the failure flag is still 0, the process should have been successful
- if (failureFlag == 0 && InChIMatch == 1) {
- Logger.info("Pre-existing successful MOPAC quantum result for "
- + name + " (" + InChIaug
- + ") has been found. This log file will be used.");
- return true;
- } else if (InChIFound == 1 && InChIMatch == 0) {// InChIs do not match (most likely due to limited name
- // length mirrored in log file (240 characters), but possibly due to a collision)
- // if(InChIPartialMatch == 1){//case where the InChI in memory begins with the InChI in the log file; we
- // will continue and check the input file, printing a warning if there is no match
- // look in the input file if the InChI doesn't match (apparently, certain characters can be deleted in
- // MOPAC output file for long InChIs)
- File inputFile = new File(directory + "/" + name + ".mop");
- if (inputFile.exists()) {// read the MOPAC inputFile
- String inputFileInChI = "";
- try {
- FileReader inI = new FileReader(inputFile);
- BufferedReader readerI = new BufferedReader(inI);
- String lineI = readerI.readLine();
- while (lineI != null) {
- if (lineI.startsWith("InChI=")) {
- inputFileInChI = lineI.trim();
- }
- lineI = readerI.readLine();
- }
- readerI.close();
- inI.close();
- } catch (Exception e) {
- String err = "Error in reading preexisting MOPAC input file \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- if (inputFileInChI.equals(InChIaug)) {
- if (failureFlag == 0) {
- Logger.info("Pre-existing successful MOPAC quantum result for "
- + name
- + " ("
- + InChIaug
- + ") has been found. This log file will be used. *Note that input file was read to confirm lack of InChIKey collision (InChI probably more than 240 characters or characters probably deleted from InChI in .out file)");
- return true;
- } else {// otherwise, failureFlag==1
- Logger.info("Pre-existing MOPAC quantum result for "
- + name
- + " ("
- + InChIaug
- + ") has been found, but the result was apparently unsuccessful. The file will be overwritten with a new calculation or Gaussian result (if available) will be used. *Note that input file was read to confirm lack of InChIKey collision (InChI probably more than 240 characters or characters probably deleted from InChI in .out file)");
- return false;
- }
- } else {
- if (inputFileInChI.equals("")) {// InChI was not found in input file
- Logger.info("*****Warning: potential InChIKey collision: InChIKey(augmented) = "
- + name
- + " RMG Augmented InChI = "
- + InChIaug
- + " Log file Augmented InChI = "
- + logFileInChI
- + " . InChI could not be found in the MOPAC input file. You should manually check that the output file contains the intended species.");
- return true;
- } else {// InChI was found but doesn't match
- Logger.critical("Congratulations! You may have discovered one of the first recorded instances of an InChIKey collision: InChIKey(augmented) = "
- + name
- + " RMG Augmented InChI = "
- + InChIaug
- + " MOPAC input file Augmented InChI = "
- + inputFileInChI
- + " Log file Augmented InChI = "
- + logFileInChI);
- System.exit(0);
- }
- }
- } else {
- Logger.info("*****Warning: potential InChIKey collision: InChIKey(augmented) = "
- + name
- + " RMG Augmented InChI = "
- + InChIaug
- + " Log file Augmented InChI = "
- + logFileInChI
- + " . MOPAC input file could not be found to check full InChI. You should manually check that the log file contains the intended species.");
- return true;
- }
- // }
- // else{
- // Logger.critical("Congratulations! You may have discovered one of the first recorded instances of an InChIKey collision: InChIKey(augmented) = "
- // + name + " RMG Augmented InChI = "+ InChIaug + " MOPAC output file Augmented InChI = "+logFileInChI);
- // System.exit(0);
- // }
- } else if (InChIFound == 0) {
- Logger.critical("An InChI was not found in file: " + name
- + ".out");
- System.exit(0);
- } else if (failureFlag == 1) {// note these should cover all possible results for this block, and if the
- // file.exists block is entered, it should return from within the block and should not reach the return statement below
- Logger.info("Pre-existing MOPAC quantum result for "
- + name
- + " ("
- + InChIaug
- + ") has been found, but the result was apparently unsuccessful. The file will be overwritten with a new calculation or Gaussian result (if available) will be used.");
- return false;
- }
- }
- // we could print a line here for cases where the file doesn't exist, but this would probably be too verbose
- return false;
- }
- // returns true if an MM4 output file for the given name and directory (.mm4out suffix) exists and indicates
- // successful completion (same criteria as used after calculation runs); terminates if the InChI doesn't match the InChI
- // in the file or if there is no InChI in the file; returns false otherwise
- public boolean successfulMM4ResultExistsQ(String name, String directory,
- String InChIaug) {
- // part of the code is taken from analogous code for MOPAC (first ~half) and Gaussian (second ~half)
- // look in the output file to check for the successful termination of the calculation (assumed to be successful
- // if "description of vibrations appears)
- int failureFlag = 1;// flag (1 or 0) indicating whether the MM4 job failed
- int failureOverrideFlag = 0;// flag (1 or 0) to override success as measured by failureFlag
- File file = new File(directory + "/" + name + ".mm4out");
- File canFile = new File(directory + "/" + name + ".canout");
- int InChIMatch = 0;// flag (1 or 0) indicating whether the InChI in the file matches InChIaug; this can only be
- // 1 if InChIFound is also 1;
- int InChIFound = 0;// flag (1 or 0) indicating whether an InChI was found in the log file
- int InChIPartialMatch = 0;// flag (1 or 0) indicating whether the InChI in the log file is a substring of the
- // InChI in RMG's memory
- if (useCanTherm) {// if we are using CanTherm, check whether a CanTherm output file exists...if it does, we will
- // continue on, otherwise, we will rerun calculations (including MM4 calculation) from scratch to ensure our atom
- // numbering is consistent; note: if the .canout file exists, we still want to check for the actual MM4 file even if we
- // are using CanTherm and reading CanTherm output because 1. it ensures we have the correct species and don't have an
- // InChI collision 2. it is needed for getting the geometry which is needed for the symmetry number corrections applied
- // to CanTherm output (which doesn't include symmetry number considerations)
- if (!canFile.exists())
- return false;
- }
- if (file.exists()) {// if the file exists, do further checks; otherwise, we will skip to final statement and
- // return false
- String logFileInChI = "";
- try {
- FileReader in = new FileReader(file);
- BufferedReader reader = new BufferedReader(in);
- String line = reader.readLine();
- while (line != null) {
- String trimLine = line.trim();
- if (trimLine.equals("STATISTICAL THERMODYNAMICS ANALYSIS")) {
- failureFlag = 0;
- } else if (trimLine.endsWith("imaginary frequencies,")) {// read the number of imaginary frequencies
- // and make sure it is zero
- String[] split = trimLine.split("\\s+");
- if (Integer.parseInt(split[3]) > 0) {
- Logger.info("*****Imaginary freqencies found:");
- failureOverrideFlag = 1;
- }
- } else if (trimLine.contains(" 0.0 (fir )")) {
- if (useCanTherm) {// zero frequencies are only acceptable when CanTherm is used
- Logger.info("*****Warning: zero freqencies found (values lower than 7.7 cm^-1 are rounded to zero in MM4 output); CanTherm should hopefully correct this:");
- } else {
- Logger.info("*****Zero freqencies found:");
- failureOverrideFlag = 1;
- }
- } else if (trimLine.startsWith("InChI=")) {
- logFileInChI = line.trim();// output files should take up to about 60 (?) characters of the name
- // in the input file
- InChIFound = 1;
- if (logFileInChI.equals(InChIaug))
- InChIMatch = 1;
- else if (InChIaug.startsWith(logFileInChI))
- InChIPartialMatch = 1;
- }
- line = reader.readLine();
- }
- reader.close();
- in.close();
- } catch (Exception e) {
- String err = "Error in reading preexisting MM4 output file \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- if (failureOverrideFlag == 1)
- failureFlag = 1; // job will be considered a failure if there are imaginary frequencies or if job
- // terminates to to excess time/cycles
- if (failureFlag == 0 && connectivityCheck > 0) {// if the process was successful, check the connectivity (if
- // requested by the user)
- boolean connectivityMatch = connectivityMatchInMM4ResultQ(name,
- directory, InChIaug);
- if (!connectivityMatch)
- if (connectivityCheck > 1) {// connectivityCheck=2
- failureFlag = 1;
- } else {// connectivityCheck=1; print warning and procede with "return true"
- Logger.warning("Connectivity in MM4 result for "
- + name
- + " ("
- + InChIaug
- + ") does not appear to match desired connectivity.");
- }
- }
- // if the failure flag is still 0, the process should have been successful
- if (failureFlag == 0 && InChIMatch == 1) {
- Logger.info("Pre-existing successful MM4 result for " + name
- + " (" + InChIaug
- + ") has been found. This log file will be used.");
- return true;
- } else if (InChIFound == 1 && InChIMatch == 0) {// InChIs do not match (most likely due to limited name
- // length mirrored in log file (79 characters), but possibly due to a collision)
- if (InChIPartialMatch == 1) {// case where the InChI in memory begins with the InChI in the log file; we
- // will continue and check the input file, printing a warning if there is no match
- File inputFile = new File(directory + "/" + name + ".mm4");
- if (inputFile.exists()) {// read the MM4 inputFile
- String inputFileInChI = "";
- try {
- FileReader inI = new FileReader(inputFile);
- BufferedReader readerI = new BufferedReader(inI);
- String lineI = readerI.readLine();
- // InChI should be repeated after in the first line of the input file
- inputFileInChI = lineI.trim().substring(80);// extract the string starting with character 81
- readerI.close();
- inI.close();
- } catch (Exception e) {
- String err = "Error in reading preexisting MM4 .mm4 file \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- if (inputFileInChI.equals(InChIaug)) {
- if (failureFlag == 0) {
- Logger.info("Pre-existing successful MM4 result for "
- + name
- + " ("
- + InChIaug
- + ") has been found. This log file will be used. *Note that input file was read to confirm lack of InChIKey collision (InChI probably more than ~60 characters)");
- return true;
- } else {// otherwise, failureFlag==1
- Logger.info("Pre-existing MM4 result for "
- + name
- + " ("
- + InChIaug
- + ") has been found, but the result was apparently unsuccessful. The file will be overwritten with a new calculation. *Note that input file was read to confirm lack of InChIKey collision (InChI probably more than ~60 characters)");
- return false;
- }
- } else {
- if (inputFileInChI.equals("")) {// InChI was not found in input file
- Logger.info("*****Warning: potential InChIKey collision: InChIKey(augmented) = "
- + name
- + " RMG Augmented InChI = "
- + InChIaug
- + " Log file Augmented InChI = "
- + logFileInChI
- + " . InChI could not be found in the MM4 input file. You should manually check that the log file contains the intended species.");
- return true;
- } else {// InChI was found but doesn't match
- Logger.critical("Congratulations! You may have discovered one of the first recorded instances of an InChIKey collision: InChIKey(augmented) = "
- + name
- + " RMG Augmented InChI = "
- + InChIaug
- + " MM4 input file Augmented InChI = "
- + inputFileInChI);
- System.exit(0);
- }
- }
- } else {
- Logger.info("*****Warning: potential InChIKey collision: InChIKey(augmented) = "
- + name
- + " RMG Augmented InChI = "
- + InChIaug
- + " Log file Augmented InChI = "
- + logFileInChI
- + " . MM4 input file could not be found to check full InChI. You should manually check that the log file contains the intended species.");
- return true;
- }
- } else {
- Logger.critical("Congratulations! You may have discovered one of the first recorded instances of an InChIKey collision: InChIKey(augmented) = "
- + name
- + " RMG Augmented InChI = "
- + InChIaug
- + " Log file Augmented InChI = " + logFileInChI);
- System.exit(0);
- }
- } else if (InChIFound == 0) {
- Logger.critical("An InChI was not found in file: " + name
- + ".mm4out");
- System.exit(0);
- } else if (failureFlag == 1) {// note these should cover all possible results for this block, and if the
- // file.exists block is entered, it should return from within the block and should not reach the return statement below
- Logger.info("Pre-existing MM4 result for "
- + name
- + " ("
- + InChIaug
- + ") has been found, but the result was apparently unsuccessful. The file will be overwritten with a new calculation.");
- return false;
- }
- }
- // we could print a line here for cases where the file doesn't exist, but this would probably be too verbose
- return false;
- }
- // //checks the MOPAC file for negative frequencies
- // public boolean MopacFileContainsNegativeFreqsQ(String name, String directory){
- // boolean negativeFreq=false;
- //
- // //code below copied from parseMopacPM3()
- // String command = "c:/Python25/python.exe c:/Python25/MopacPM3ParsingScript.py ";//this should eventually be modified
- // for added generality
- // String logfilepath=directory+"/"+name+".out";
- // command=command.concat(logfilepath);
- //
- // //much of code below is copied from calculateThermoFromPM3Calc()
- // //parse the Mopac file using cclib
- // int natoms = 0; //number of atoms from Mopac file; in principle, this should agree with number of chemGraph atoms
- // ArrayList atomicNumber = new ArrayList(); //vector of atomic numbers (integers) (apparently Vector is thread-safe;
- // cf. http://answers.yahoo.com/question/index?qid=20081214065127AArZDT3; ...should I be using this instead?)
- // ArrayList x_coor = new ArrayList(); //vectors of x-, y-, and z-coordinates (doubles) (Angstroms) (in order
- // corresponding to above atomic numbers)
- // ArrayList y_coor = new ArrayList();
- // ArrayList z_coor = new ArrayList();
- // double energy = 0; //PM3 energy (Hf298) in Hartree (***note: in the case of MOPAC, the MOPAC file will contain in
- // units of kcal/mol, but modified ccLib will return in Hartree)
- // double molmass = 0; //molecular mass in amu
- // ArrayList freqs = new ArrayList(); //list of frequencies in units of cm^-1
- // double rotCons_1 = 0;//rotational constants in (1/s)
- // double rotCons_2 = 0;
- // double rotCons_3 = 0;
- // //int gdStateDegen = p_chemGraph.getRadicalNumber()+1;//calculate ground state degeneracy from the number of
- // radicals; this should give the same result as spin multiplicity in Gaussian input file (and output file), but we do
- // not explicitly check this (we could use "mult" which cclib reads in if we wanted to do so); also, note that this is
- // not always correct, as there can apparently be additional spatial degeneracy for non-symmetric linear molecules like
- // OH radical (cf. http://cccbdb.nist.gov/thermo.asp)
- // try{
- // File runningdir=new File(directory);
- // Process cclibProc = Runtime.getRuntime().exec(command, null, runningdir);
- // //read the stdout of the process, which should contain the desired information in a particular format
- // InputStream is = cclibProc.getInputStream();
- // InputStreamReader isr = new InputStreamReader(is);
- // BufferedReader br = new BufferedReader(isr);
- // String line=null;
- // //example output:
- // // C:\Python25>python.exe GaussianPM3ParsingScript.py TEOS.out
- // // 33
- // // [ 6 6 8 14 8 6 6 8 6 6 8 6 6 1 1 1 1 1 1 1 1 1 1 1 1
- // // 1 1 1 1 1 1 1 1]
- // // [[ 2.049061 -0.210375 3.133106]
- // // [ 1.654646 0.321749 1.762752]
- // // [ 0.359284 -0.110429 1.471465]
- // // [-0.201871 -0.013365 -0.12819 ]
- // // [ 0.086307 1.504918 -0.82893 ]
- // // [-0.559186 2.619928 -0.284003]
- // // [-0.180246 3.839463 -1.113029]
- // // [ 0.523347 -1.188305 -1.112765]
- // // [ 1.857584 -1.018167 -1.495088]
- // // [ 2.375559 -2.344392 -2.033403]
- // // [-1.870397 -0.297297 -0.075427]
- // // [-2.313824 -1.571765 0.300245]
- // // [-3.83427 -1.535927 0.372171]
- // // [ 1.360346 0.128852 3.917699]
- // // [ 2.053945 -1.307678 3.160474]
- // // [ 3.055397 0.133647 3.403037]
- // // [ 1.677262 1.430072 1.750899]
- // // [ 2.372265 -0.029237 0.985204]
- // // [-0.245956 2.754188 0.771433]
- // // [-1.656897 2.472855 -0.287156]
- // // [-0.664186 4.739148 -0.712606]
- // // [-0.489413 3.734366 -2.161038]
- // // [ 0.903055 4.016867 -1.112198]
- // // [ 1.919521 -0.229395 -2.269681]
- // // [ 2.474031 -0.680069 -0.629949]
- // // [ 2.344478 -3.136247 -1.273862]
- // // [ 1.786854 -2.695974 -2.890647]
- // // [ 3.41648 -2.242409 -2.365094]
- // // [-1.884889 -1.858617 1.28054 ]
- // // [-1.976206 -2.322432 -0.440995]
- // // [-4.284706 -1.26469 -0.591463]
- // // [-4.225999 -2.520759 0.656131]
- // // [-4.193468 -0.809557 1.112677]]
- // // -14.1664924726
- // // [ 9.9615 18.102 27.0569 31.8459 39.0096 55.0091
- // // 66.4992 80.4552 86.4912 123.3551 141.6058 155.5448
- // // 159.4747 167.0013 178.5676 207.3738 237.3201 255.3487
- // // 264.5649 292.867 309.4248 344.6503 434.8231 470.2074
- // // 488.9717 749.1722 834.257 834.6594 837.7292 839.6352
- // // 887.9767 892.9538 899.5374 992.1851 1020.6164 1020.8671
- // // 1028.3897 1046.7945 1049.1768 1059.4704 1065.1505 1107.4001
- // // 1108.1567 1109.0466 1112.6677 1122.7785 1124.4315 1128.4163
- // // 1153.3438 1167.6705 1170.9627 1174.9613 1232.1826 1331.8459
- // // 1335.3932 1335.8677 1343.9556 1371.37 1372.8127 1375.5428
- // // 1396.0344 1402.4082 1402.7554 1403.2463 1403.396 1411.6946
- // // 1412.2456 1412.3519 1414.5982 1415.3613 1415.5698 1415.7993
- // // 1418.5409 2870.7446 2905.3132 2907.0361 2914.1662 2949.2646
- // // 2965.825 2967.7667 2971.5223 3086.3849 3086.3878 3086.6448
- // // 3086.687 3089.2274 3089.4105 3089.4743 3089.5841 3186.0753
- // // 3186.1375 3186.3511 3186.365 ]
- // // [ 0.52729 0.49992 0.42466]
- // //note: above example has since been updated to print molecular mass; also frequency and atomic number format has
- // been updated
- // String [] stringArray;
- // natoms = Integer.parseInt(br.readLine());//read line 1: number of atoms
- // stringArray = br.readLine().replace("[", "").replace("]","").trim().split(",\\s+");//read line 2: the atomic numbers
- // (first removing braces)
- // // line = br.readLine().replace("[", "").replace("]","");//read line 2: the atomic numbers (first removing braces)
- // // StringTokenizer st = new StringTokenizer(line); //apprently the stringTokenizer class is deprecated, but I am
- // having trouble getting the regular expressions to work properly
- // for(int i=0; i < natoms; i++){
- // // atomicNumber.add(i,Integer.parseInt(stringArray[i]));
- // atomicNumber.add(i,Integer.parseInt(stringArray[i]));
- // }
- // for(int i=0; i < natoms; i++){
- // stringArray = br.readLine().replace("[", "").replace("]","").trim().split("\\s+");//read line 3+: coordinates for
- // atom i; used /s+ for split; using spaces with default limit of 0 was giving empty string
- // x_coor.add(i,Double.parseDouble(stringArray[0]));
- // y_coor.add(i,Double.parseDouble(stringArray[1]));
- // z_coor.add(i,Double.parseDouble(stringArray[2]));
- // }
- // energy = Double.parseDouble(br.readLine());//read next line: energy
- // molmass = Double.parseDouble(br.readLine());//read next line: molecular mass (in amu)
- // if (natoms>1){//read additional info for non-monoatomic species
- // stringArray = br.readLine().replace("[", "").replace("]","").trim().split(",\\s+");//read next line: frequencies
- // for(int i=0; i < stringArray.length; i++){
- // freqs.add(i,Double.parseDouble(stringArray[i]));
- // }
- // stringArray = br.readLine().replace("[", "").replace("]","").trim().split("\\s+");//read next line rotational
- // constants (converting from GHz to Hz in the process)
- // rotCons_1 = Double.parseDouble(stringArray[0])*1000000000;
- // rotCons_2 = Double.parseDouble(stringArray[1])*1000000000;
- // rotCons_3 = Double.parseDouble(stringArray[2])*1000000000;
- // }
- // while ( (line = br.readLine()) != null) {
- // //do nothing (there shouldn't be any more information, but this is included to get all the output)
- // }
- // int exitValue = cclibProc.waitFor();
- // }
- // catch (Exception e) {
- // Logger.logStackTrace(e);
- // String err = "Error in running ccLib Python process \n";
- // err += e.toString();
- // Logger.critical(err);
- // System.exit(0);
- // }
- //
- // //start of code "new" to this function (aside from initialization of negativeFreq)
- // if(natoms > 0){
- // for (int i=0; i<freqs.size(); i++){
- // if((Double)freqs.get(i) < 0) negativeFreq = true;
- // }
- // }
- // return negativeFreq;
- // }
- // check the connectivity in a Gaussian/MOPAC result (or XYZ file); returns true if the there appears to be a match,
- // and returns false otherwise
- // the connectivity is assumed to match if the InChI produced by processing the result/coordinates through OpenBabel
- // into MOL file and then using InChI utility produces an InChI that is equivalent (following stereochemical layer
- // removal) to the InChI stored in memory
- // we do not need/want to compare the augmented multN part of the InChI, so InChI will be converted within to the
- // non-modified version
- // for gaussian03, outfileExtension should be ".log" and babelType should be "g03"
- // for MOPAC, outfileExtension should be ".out" and babelType should be "moo"
- // for XYZ, outfileExtension should be ".xyz" and babelType should be "xyz"
- public boolean connectivityMatchInPM3ResultQ(String name, String directory,
- String augInChI, String outfileExtension, String babelType) {
- // Step 0: get the un-modified InChI
- String InChI = getInChIFromModifiedInChI(augInChI);
- // Step 1: call the OpenBabel process (note that this requires OpenBabel environment variable) to convert to MOL
- // file
- String molPath = directory + "/" + name + ".mol";// the path to the MOL file to be created (with connectivity);
- try {
- File runningdir = new File(directory);
- String command = null;
- String inpPath = directory + "/" + name + outfileExtension;// the path to the QM output file with
- // coordinates
- if (System.getProperty("os.name").toLowerCase().contains("windows")) {// special windows case
- command = "babel -i" + babelType + " \"" + inpPath
- + "\" -omol \"" + molPath + "\"";
- } else {
- command = "babel -i" + babelType + " " + inpPath + " -omol "
- + molPath;
- }
- Process babelProc = Runtime.getRuntime().exec(command, null,
- runningdir);
- // read in output
- InputStream is = babelProc.getInputStream();
- InputStreamReader isr = new InputStreamReader(is);
- BufferedReader br = new BufferedReader(isr);
- String line = null;
- while ((line = br.readLine()) != null) {
- // do nothing
- // Logger.info(line);
- }
- int exitValue = babelProc.waitFor();
- babelProc.getErrorStream().close();
- babelProc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- String err = "Error in running OpenBabel G03/MOO/XYZ to MOL process \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // Step 2. convert the MOL file to InChI (with stereochem layers removed)
- String[] result = Species.runInChIProcess(new File(molPath), new File(
- System.getProperty("RMG.InChI_running_directory")
- + "/species.txt"), true);
- String InChI3D = stripStereochemLayersFromInChI(result[0]);
- // Step 3. check whether there is a match (i.e. InChI equals InChI3D)
- if (InChI3D.equals(InChI)) {
- return true;
- } else {// if there is a mismatch using OpenBabel (either due to OpenBabel segmentation fault (e.g.
- // http://sourceforge.net/tracker/?func=detail&aid=3417992&group_id=40728&atid=428740) or tolerances that are too tight
- // (e.g. InChI=1/C10H12O2/c1-2-6-9(5-1)11-12-10-7-3-4-8-10/h1,7H,2-4,6,8H2)), print a warning and go to the backup check
- // using MoleCoor
- Logger.info("***For species "
- + name
- + " an OpenBabel-based check suggests the optimized three-dimensional InChI ("
- + InChI3D
- + ") does not match the intended (unmodified) InChI ("
- + InChI + "). Will retry connectivity check with MoleCoor");
- return connectivityMatchInPM3ResultQSecondary(name, directory,
- augInChI, babelType);
- }
- }
- // a backup/fallback connectivity check using MoleCoor for when the primary approach using OpenBabel fails
- public boolean connectivityMatchInPM3ResultQSecondary(String name,
- String directory, String augInChI, String babelType) {
- // Step 0: get the un-modified InChI
- String InChI = getInChIFromModifiedInChI(augInChI);
- String molPath = directory + "/" + name + ".mol";// the path to the MOL file to be created (with connectivity);
- // Step 1. parse the file according to the program used, as determined by babelType (we could just use
- // coordinates from a pre-existing MOL file, but that wouldn't handle cases where OpenBabel crashed due to segmentation
- // fault
- String command = null;
- if (babelType.equals("moo"))
- command = getMopacPM3ParseCommand(name, directory); // MOPAC case
- else if (babelType.equals("g03"))
- command = getGaussianPM3ParseCommand(name, directory); // G03 case
- else if (babelType.equals("xyz"))
- command = getMM4ParseCommand(name, directory); // MM4 case
- else {
- Logger.critical("Unrecognized babelType during backup connectivity check");
- System.exit(0);
- }
- QMData qmdata = getQMDataWithCClib(name, directory, command, false);
- // Step 2a. write the coordinates as an xyz file
- String inpFilePath = directory + "/" + name + ".xyz";
- // construct the string to write
- int natoms = qmdata.natoms;
- ArrayList atomicNumber = qmdata.atomicNumber; // vector of atomic numbers (integers) (apparently Vector is
- // thread-safe; cf. http://answers.yahoo.com/question/index?qid=20081214065127AArZDT3; ...should I be using this
- // instead?)
- ArrayList x_coor = qmdata.x_coor; // vectors of x-, y-, and z-coordinates (doubles) (Angstroms) (in order
- // corresponding to above atomic numbers)
- ArrayList y_coor = qmdata.y_coor;
- ArrayList z_coor = qmdata.z_coor;
- String geomXYZ = natoms + "\n";
- geomXYZ += augInChI + "\n"; // use modified/augmented InChI as molecule name
- for (int i = 0; i < natoms; i++) {
- geomXYZ += getAtomicSymbol((Integer) atomicNumber.get(i)) + " "
- + x_coor.get(i) + " " + y_coor.get(i) + " " + z_coor.get(i)
- + "\n";
- }
- // write the string to file
- try {
- // File inputFilePath=new File(qmfolder,"symminput.txt");//SYMMETRY program directory
- FileWriter fw = new FileWriter(inpFilePath);
- fw.write(geomXYZ);
- fw.close();
- } catch (IOException e) {
- String err = "Error writing XYZ file for backup MoleCoor connectivity perception";
- err += e.toString();
- Logger.critical(err);
- System.exit(0);
- }
- // Step 2b. call the MoleCoor script to perceive connectivity and write MOL file (with connectivity)
- try {
- File runningdir = new File(directory);
- String command2 = null;
- if (System.getProperty("os.name").toLowerCase().contains("windows")) {// special windows case where paths
- // can have spaces and are allowed to be surrounded by quotes
- command2 = "python \"" + System.getenv("RMG")
- + "/scripts/XYZtoMOLconn.py\" ";
- // first argument: input file path
- command2 = command2.concat("\"" + inpFilePath + "\" ");
- // second argument: output path
- command2 = command2.concat("\"" + molPath + "\" ");
- // third argument: molecule name (the modified InChI)
- command2 = command2.concat(augInChI + " ");
- // fourth argument: tolerance
- command2 = command2.concat("0.50 ");// InChI=1/C10H12O2/c1-2-6-9(5-1)11-12-10-7-3-4-8-10/h1,7H,2-4,6,8H2
- // case has very strained with long O-O bond length; O-O bond 1.83873 Ang.; would need tol = 0.47873 or higher with
- // MoleCoor scheme; 0.50 is used here
- // fifth argument: PYTHONPATH
- command2 = command2.concat("\"" + System.getenv("RMG")
- + "/source/MoleCoor\"");// this will pass $RMG/source/MoleCoor to the script (in order to get
- // the appropriate path for importing
- } else {// non-Windows case
- command2 = "python " + System.getenv("RMG")
- + "/scripts/XYZtoMOLconn.py ";
- // first argument: input file path
- command2 = command2.concat(inpFilePath + " ");
- // second argument: output path
- command2 = command2.concat(molPath + " ");
- // third argument: molecule name (the modified InChI)
- command2 = command2.concat(augInChI + " ");
- // fourth argument: tolerance
- command2 = command2.concat("0.50 ");// InChI=1/C10H12O2/c1-2-6-9(5-1)11-12-10-7-3-4-8-10/h1,7H,2-4,6,8H2
- // case has very strained with long O-O bond length; O-O bond 1.83873 Ang.; would need tol = 0.47873 or higher with
- // MoleCoor scheme; 0.50 is used here
- // fifth argument: PYTHONPATH
- command2 = command2.concat(System.getenv("RMG")
- + "/source/MoleCoor");// this will pass $RMG/source/MoleCoor to the script (in order to get the
- // appropriate path for importing
- }
- Process molecoorProc = Runtime.getRuntime().exec(command2, null,
- runningdir);
- // read in output
- InputStream is = molecoorProc.getInputStream();
- InputStreamReader isr = new InputStreamReader(is);
- BufferedReader br = new BufferedReader(isr);
- String line = null;
- while ((line = br.readLine()) != null) {
- // do nothing
- }
- int exitValue = molecoorProc.waitFor();
- molecoorProc.getErrorStream().close();
- molecoorProc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- String err = "Error in running MoleCoor .XYZ to .MOL process \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // Step 3. convert the MOL file to InChI (with stereochem layers removed)
- String[] result = Species.runInChIProcess(new File(molPath), new File(
- System.getProperty("RMG.InChI_running_directory")
- + "/species.txt"), true);
- String InChI3D = stripStereochemLayersFromInChI(result[0]);
- // Step 4. check whether there is a match (i.e. InChI equals InChI3D)
- if (InChI3D.equals(InChI)) {
- Logger.info("For species "
- + name
- + " a MoleCoor-based check suggests the optimized three-dimensional InChI ("
- + InChI3D
- + ") actually DOES match the intended (unmodified) InChI ("
- + InChI
- + "). Therefore, the result will be assumed to be successful.");
- return true;
- } else {
- Logger.info("***For species "
- + name
- + " a MoleCoor-based check suggests the optimized three-dimensional InChI ("
- + InChI3D
- + ") does not match the intended (unmodified) InChI ("
- + InChI + ").");
- return false;
- }
- }
- // check the connectivity in a MM4 result; returns true if the there appears to be a match, and returns false
- // otherwise
- // the connectivity is assumed to match if the InChI produced by processing the result/coordinates through OpenBabel
- // produces in InChI that is equivalent or superstring to the InChI stored in memory
- // we do not need/want to compare the augmented multN part of the InChI, so InChI will be converted within to the
- // non-modified version
- public boolean connectivityMatchInMM4ResultQ(String name, String directory,
- String augInChI) {
- // first, we convert the .mm4opt file into an xyz file that can be processed by OpenBabel; we use a python
- // script that calls MoleCoor for this
- // call the MoleCoor script to create the XYZ file from the .mm4opt file
- try {
- File runningdir = new File(directory);
- // this will only be run on Linux so we don't have to worry about Linux vs. Windows issues
- String command = "python " + System.getenv("RMG")
- + "/scripts/XYZfromMM4.py ";
- // first argument: input file path
- command = command.concat(directory + "/" + name + ".mm4opt ");
- // second argument: output path
- String inpfilepath = directory + "/" + name + ".xyz";
- command = command.concat(inpfilepath + " ");
- // third argument: molecule name (the modified InChI)
- command = command.concat(augInChI + " ");
- // fourth argument: PYTHONPATH
- command = command.concat(System.getenv("RMG") + "/source/MoleCoor");// this will pass $RMG/source/MoleCoor
- // to the script (in order to get the appropriate path for importing
- Process molecoorProc = Runtime.getRuntime().exec(command, null,
- runningdir);
- // read in output
- InputStream is = molecoorProc.getInputStream();
- InputStreamReader isr = new InputStreamReader(is);
- BufferedReader br = new BufferedReader(isr);
- String line = null;
- while ((line = br.readLine()) != null) {
- // do nothing
- }
- int exitValue = molecoorProc.waitFor();
- molecoorProc.getErrorStream().close();
- molecoorProc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- String err = "Error in running MoleCoor .MM4OPT to .XYZ process \n";
- err += e.toString();
- Logger.logStackTrace(e);
- System.exit(0);
- }
- // check whether there is a match (i.e. InChI is a substring of InChI3D); note that this makes use of the
- // connectivityMatchInPM3ResultQ function, which is sufficiently general to process XYZ files
- return connectivityMatchInPM3ResultQ(name, directory, augInChI, ".xyz",
- "xyz");
- }
- // ## operation initGAGroupLibrary()
- protected void initGAGroupLibrary() {
- // #[ operation initGAGroupLibrary()
- thermoLibrary = ThermoGAGroupLibrary.getINSTANCE();
- // #]
- }
- public QMData getMM4QMDataWithCClib(String name, String directory,
- boolean getStericEnergy) {
- String command = "python " + System.getProperty("RMG.workingDirectory")
- + "/scripts/MM4ParsingScript.py ";
- String logfilepath = directory + "/" + name + ".mm4out";
- command = command.concat(logfilepath);
- command = command.concat(" " + System.getenv("RMG") + "/source");// this will pass $RMG/source to the script (in
- // order to get the appropriate path for importing
- if (getStericEnergy)
- command = command.concat(" 1");// option to print stericenergy before molar mass (this will only be used in
- // useHindRot cases, but it is always read in with this function
- return getQMDataWithCClib(name, directory, command, getStericEnergy);
- }
- public QMData getQMDataWithCClib(String name, String directory,
- String command, boolean getStericEnergy) {
- // parse the file using cclib
- int natoms = 0; // number of atoms from Mopac file; in principle, this should agree with number of chemGraph
- // atoms
- ArrayList atomicNumber = new ArrayList(); // vector of atomic numbers (integers) (apparently Vector is
- // thread-safe; cf. http://answers.yahoo.com/question/index?qid=20081214065127AArZDT3; ...should I be using this
- // instead?)
- ArrayList x_coor = new ArrayList(); // vectors of x-, y-, and z-coordinates (doubles) (Angstroms) (in order
- // corresponding to above atomic numbers)
- ArrayList y_coor = new ArrayList();
- ArrayList z_coor = new ArrayList();
- double energy = 0; // energy (Hf298) in Hartree
- double molmass = 0; // molecular mass in amu
- double stericEnergy = 99999;// steric energy in Hartree
- ArrayList freqs = new ArrayList(); // list of frequencies in units of cm^-1
- double rotCons_1 = 0;// rotational constants in (1/s)
- double rotCons_2 = 0;
- double rotCons_3 = 0;
- try {
- Process cclibProc = Runtime.getRuntime().exec(command);
- // read the stdout of the process, which should contain the desired information in a particular format
- InputStream is = cclibProc.getInputStream();
- InputStreamReader isr = new InputStreamReader(is);
- BufferedReader br = new BufferedReader(isr);
- String line = null;
- // example output:
- // C:\Python25>python.exe GaussianPM3ParsingScript.py TEOS.out
- // 33
- // [ 6 6 8 14 8 6 6 8 6 6 8 6 6 1 1 1 1 1 1 1 1 1 1 1 1
- // 1 1 1 1 1 1 1 1]
- // [[ 2.049061 -0.210375 3.133106]
- // [ 1.654646 0.321749 1.762752]
- // [ 0.359284 -0.110429 1.471465]
- // [-0.201871 -0.013365 -0.12819 ]
- // [ 0.086307 1.504918 -0.82893 ]
- // [-0.559186 2.619928 -0.284003]
- // [-0.180246 3.839463 -1.113029]
- // [ 0.523347 -1.188305 -1.112765]
- // [ 1.857584 -1.018167 -1.495088]
- // [ 2.375559 -2.344392 -2.033403]
- // [-1.870397 -0.297297 -0.075427]
- // [-2.313824 -1.571765 0.300245]
- // [-3.83427 -1.535927 0.372171]
- // [ 1.360346 0.128852 3.917699]
- // [ 2.053945 -1.307678 3.160474]
- // [ 3.055397 0.133647 3.403037]
- // [ 1.677262 1.430072 1.750899]
- // [ 2.372265 -0.029237 0.985204]
- // [-0.245956 2.754188 0.771433]
- // [-1.656897 2.472855 -0.287156]
- // [-0.664186 4.739148 -0.712606]
- // [-0.489413 3.734366 -2.161038]
- // [ 0.903055 4.016867 -1.112198]
- // [ 1.919521 -0.229395 -2.269681]
- // [ 2.474031 -0.680069 -0.629949]
- // [ 2.344478 -3.136247 -1.273862]
- // [ 1.786854 -2.695974 -2.890647]
- // [ 3.41648 -2.242409 -2.365094]
- // [-1.884889 -1.858617 1.28054 ]
- // [-1.976206 -2.322432 -0.440995]
- // [-4.284706 -1.26469 -0.591463]
- // [-4.225999 -2.520759 0.656131]
- // [-4.193468 -0.809557 1.112677]]
- // -14.1664924726
- // [ 9.9615 18.102 27.0569 31.8459 39.0096 55.0091
- // 66.4992 80.4552 86.4912 123.3551 141.6058 155.5448
- // 159.4747 167.0013 178.5676 207.3738 237.3201 255.3487
- // 264.5649 292.867 309.4248 344.6503 434.8231 470.2074
- // 488.9717 749.1722 834.257 834.6594 837.7292 839.6352
- // 887.9767 892.9538 899.5374 992.1851 1020.6164 1020.8671
- // 1028.3897 1046.7945 1049.1768 1059.4704 1065.1505 1107.4001
- // 1108.1567 1109.0466 1112.6677 1122.7785 1124.4315 1128.4163
- // 1153.3438 1167.6705 1170.9627 1174.9613 1232.1826 1331.8459
- // 1335.3932 1335.8677 1343.9556 1371.37 1372.8127 1375.5428
- // 1396.0344 1402.4082 1402.7554 1403.2463 1403.396 1411.6946
- // 1412.2456 1412.3519 1414.5982 1415.3613 1415.5698 1415.7993
- // 1418.5409 2870.7446 2905.3132 2907.0361 2914.1662 2949.2646
- // 2965.825 2967.7667 2971.5223 3086.3849 3086.3878 3086.6448
- // 3086.687 3089.2274 3089.4105 3089.4743 3089.5841 3186.0753
- // 3186.1375 3186.3511 3186.365 ]
- // [ 0.52729 0.49992 0.42466]
- // note: above example has since been updated to print molecular mass and steric energy; also frequency and atomic
- // number format has been updated
- String[] stringArray;
- natoms = Integer.parseInt(br.readLine());// read line 1: number of atoms
- stringArray = br.readLine().replace("[", "").replace("]", "")
- .trim().split(",\\s+");// read line 2: the atomic numbers (first removing braces)
- // line = br.readLine().replace("[", "").replace("]","");//read line 2: the atomic numbers (first removing
- // braces)
- // StringTokenizer st = new StringTokenizer(line); //apprently the stringTokenizer class is deprecated, but
- // I am having trouble getting the regular expressions to work properly
- for (int i = 0; i < natoms; i++) {
- // atomicNumber.add(i,Integer.parseInt(stringArray[i]));
- atomicNumber.add(i, Integer.parseInt(stringArray[i]));
- }
- for (int i = 0; i < natoms; i++) {
- stringArray = br.readLine().replace("[", "").replace("]", "")
- .trim().split("\\s+");// read line 3+: coordinates for atom i; used /s+ for split; using spaces
- // with default limit of 0 was giving empty string
- x_coor.add(i, Double.parseDouble(stringArray[0]));
- y_coor.add(i, Double.parseDouble(stringArray[1]));
- z_coor.add(i, Double.parseDouble(stringArray[2]));
- }
- energy = Double.parseDouble(br.readLine());// read next line: energy
- if (getStericEnergy)
- stericEnergy = Double.parseDouble(br.readLine());// read next line: steric energy (in Hartree)
- molmass = Double.parseDouble(br.readLine());// read next line: molecular mass (in amu)
- if (natoms > 1) {// read additional info for non-monoatomic species
- stringArray = br.readLine().replace("[", "").replace("]", "")
- .trim().split(",\\s+");// read next line: frequencies
- for (int i = 0; i < stringArray.length; i++) {
- freqs.add(i, Double.parseDouble(stringArray[i]));
- }
- stringArray = br.readLine().replace("[", "").replace("]", "")
- .trim().split("\\s+");// read next line rotational constants (converting from GHz to Hz in the
- // process)
- rotCons_1 = Double.parseDouble(stringArray[0]) * 1000000000;
- rotCons_2 = Double.parseDouble(stringArray[1]) * 1000000000;
- rotCons_3 = Double.parseDouble(stringArray[2]) * 1000000000;
- }
- while ((line = br.readLine()) != null) {
- // do nothing (there shouldn't be any more information, but this is included to get all the output)
- }
- int exitValue = cclibProc.waitFor();
- cclibProc.getErrorStream().close();
- cclibProc.getOutputStream().close();
- br.close();
- isr.close();
- is.close();
- } catch (Exception e) {
- Logger.logStackTrace(e);
- String err = "Error in running ccLib Python process \n";
- err += e.toString();
- Logger.critical(err);
- System.exit(0);
- }
- // package up the result
- return new QMData(natoms, atomicNumber, x_coor, y_coor, z_coor, energy,
- stericEnergy, molmass, freqs, rotCons_1, rotCons_2, rotCons_3);
- }
- // converts atomicnumber to atomic symbol
- // returns Err if the atomic number is not recognized; currently, H, C, O, N, S, and Si are recognized
- public String getAtomicSymbol(Integer atomicNumber) {
- if (atomicNumber == 1)
- return "H";
- else if (atomicNumber == 6)
- return "C";
- else if (atomicNumber == 8)
- return "O";
- else if (atomicNumber == 7)
- return "N";
- else if (atomicNumber == 16)
- return "S";
- else if (atomicNumber == 14)
- return "Si";
- else
- return "Err";//
- }
- }
- /*********************************************************************
- * File Path : RMG\RMG\jing\chem\QMTP.java
- *********************************************************************/