Skip to content

Commit

Permalink
Merge pull request #1571 from ReactionMechanismGenerator/auto_optical…
Browse files Browse the repository at this point in the history
…_isomer_ex_symmetry_arkane

Auto optical isomer and external symmetry to Arkane
  • Loading branch information
alongd committed Apr 25, 2019
2 parents 88e62a8 + 009a23c commit 4b1c083
Show file tree
Hide file tree
Showing 8 changed files with 197 additions and 70 deletions.
18 changes: 9 additions & 9 deletions arkane/gaussian.py
Expand Up @@ -38,11 +38,12 @@
from rmgpy.exceptions import InputError

from arkane.common import check_conformer_energy, get_element_mass
from arkane.log import Log

################################################################################


class GaussianLog:
class GaussianLog(Log):
"""
Represent a log file from Gaussian. The attribute `path` refers to the
location on disk of the Gaussian log file of interest. Methods are provided
Expand Down Expand Up @@ -156,7 +157,7 @@ def loadGeometry(self):

return coord, number, mass

def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, symfromlog=None, label=''):
def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=None, label=''):
"""
Load the molecular degree of freedom data from a log file created as
the result of a Gaussian "Freq" quantum chemistry calculation. As
Expand All @@ -170,7 +171,12 @@ def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, sym
modes = []
unscaled_frequencies = []
E0 = 0.0

if opticalIsomers is None or symmetry is None:
_opticalIsomers, _symmetry = self.get_optical_isomers_and_symmetry_number()
if opticalIsomers is None:
opticalIsomers = _opticalIsomers
if symmetry is None:
symmetry = _symmetry
f = open(self.path, 'r')
line = f.readline()
while line != '':
Expand All @@ -197,11 +203,6 @@ def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, sym
translation = IdealGasTranslation(mass=(mass,"amu"))
modes.append(translation)

# Read Gaussian's estimate of the external symmetry number
elif 'Rotational symmetry number' in line and symmetry is None:
if symfromlog is True:
symmetry = int(float(line.split()[3]))

# Read moments of inertia for external rotational modes
elif 'Rotational constants (GHZ):' in line:
inertia = [float(d) for d in line.split()[-3:]]
Expand Down Expand Up @@ -251,7 +252,6 @@ def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, sym

# Close file when finished
f.close()

return Conformer(E0=(E0*0.001,"kJ/mol"), modes=modes, spinMultiplicity=spinMultiplicity,
opticalIsomers=opticalIsomers), unscaled_frequencies

Expand Down
34 changes: 31 additions & 3 deletions arkane/gaussianTest.py
Expand Up @@ -37,6 +37,7 @@
from external.wip import work_in_progress

from arkane.gaussian import GaussianLog
from arkane.statmech import determine_qm_software

################################################################################

Expand All @@ -54,7 +55,7 @@ def testLoadEthyleneFromGaussianLog_CBSQB3(self):
"""

log = GaussianLog(os.path.join(os.path.dirname(__file__),'data','ethylene.log'))
conformer, unscaled_frequencies = log.loadConformer(symfromlog=True)
conformer, unscaled_frequencies = log.loadConformer()
E0 = log.loadEnergy()

self.assertTrue(len([mode for mode in conformer.modes if isinstance(mode,IdealGasTranslation)]) == 1)
Expand All @@ -81,7 +82,7 @@ def testLoadOxygenFromGaussianLog(self):
"""

log = GaussianLog(os.path.join(os.path.dirname(__file__),'data','oxygen.log'))
conformer, unscaled_frequencies = log.loadConformer(symfromlog=True)
conformer, unscaled_frequencies = log.loadConformer()
E0 = log.loadEnergy()

self.assertTrue(len([mode for mode in conformer.modes if isinstance(mode,IdealGasTranslation)]) == 1)
Expand Down Expand Up @@ -109,7 +110,7 @@ def testLoadEthyleneFromGaussianLog_G3(self):
"""

log = GaussianLog(os.path.join(os.path.dirname(__file__),'data','ethylene_G3.log'))
conformer, unscaled_frequencies = log.loadConformer(symfromlog=True)
conformer, unscaled_frequencies = log.loadConformer()
E0 = log.loadEnergy()

self.assertTrue(len([mode for mode in conformer.modes if isinstance(mode,IdealGasTranslation)]) == 1)
Expand All @@ -130,5 +131,32 @@ def testLoadEthyleneFromGaussianLog_G3(self):
self.assertEqual(conformer.spinMultiplicity, 1)
self.assertEqual(conformer.opticalIsomers, 1)

def testLoadSymmetryAndOptics(self):
"""
Uses a Gaussian03 log file for oxygen (O2) to test that its
molecular degrees of freedom can be properly read.
"""

log = GaussianLog(os.path.join(os.path.dirname(__file__),'data','oxygen.log'))
optical, symmetry = log.get_optical_isomers_and_symmetry_number()
self.assertEqual(optical,1)
self.assertEqual(symmetry,2)

conf = log.loadConformer()[0]
self.assertEqual(conf.opticalIsomers, 1)
found_rotor = False
for mode in conf.modes:
if isinstance(mode,LinearRotor):
self.assertEqual(mode.symmetry,2)
found_rotor = True
self.assertTrue(found_rotor)

def testDetermineQMSoftware(self):
"""
Ensures that determine_qm_software returns a GaussianLog object
"""
log = determine_qm_software(os.path.join(os.path.dirname(__file__),'data','oxygen.log'))
self.assertIsInstance(log,GaussianLog)

if __name__ == '__main__':
unittest.main( testRunner = unittest.TextTestRunner(verbosity=2) )
113 changes: 113 additions & 0 deletions arkane/log.py
@@ -0,0 +1,113 @@
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
A general class for parsing quantum mechanical log files
"""
import os.path
import logging
import shutil

from rmgpy.qm.qmdata import QMData
from rmgpy.qm.symmetry import PointGroupCalculator

class Log(object):
"""
Represent a general log file.
The attribute `path` refers to the location on disk of the log file of interest.
"""
def __init__(self, path):
self.path = path

def getNumberOfAtoms(self):
"""
Return the number of atoms in the molecular configuration used in
the MolPro log file.
"""
raise NotImplementedError("loadGeometry is not implemented for the Log class")

def loadForceConstantMatrix(self):
"""
Return the force constant matrix (in Cartesian coordinates) from the
QChem log file. If multiple such matrices are identified,
only the last is returned. The units of the returned force constants
are J/m^2. If no force constant matrix can be found in the log file,
``None`` is returned.
"""
raise NotImplementedError("loadGeometry is not implemented for the Log class")

def loadGeometry(self):
"""
Return the optimum geometry of the molecular configuration from the
log file. If multiple such geometries are identified, only the
last is returned.
"""
raise NotImplementedError("loadGeometry is not implemented for the Log class")

def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=None, label=''):
"""
Load the molecular degree of freedom data from an output file created as the result of a
QChem "Freq" calculation. As QChem's guess of the external symmetry number is not always correct,
you can use the `symmetry` parameter to substitute your own value;
if not provided, the value in the QChem output file will be adopted.
"""
raise NotImplementedError("loadGeometry is not implemented for the Log class")

def loadEnergy(self, frequencyScaleFactor=1.):
"""
Load the energy in J/mol from a QChem log file. Only the last energy
in the file is returned. The zero-point energy is *not* included in
the returned value.
"""
raise NotImplementedError("loadGeometry is not implemented for the Log class")

def loadZeroPointEnergy(self,frequencyScaleFactor=1.):
"""
Load the unscaled zero-point energy in J/mol from a QChem output file.
"""
raise NotImplementedError("loadGeometry is not implemented for the Log class")

def loadScanEnergies(self):
"""
Extract the optimized energies in J/mol from a QChem log file, e.g. the
result of a QChem "PES Scan" quantum chemistry calculation.
"""
raise NotImplementedError("loadGeometry is not implemented for the Log class")

def loadNegativeFrequency(self):
"""
Return the imaginary frequency from a transition state frequency
calculation in cm^-1.
"""
raise NotImplementedError("loadGeometry is not implemented for the Log class")

def get_optical_isomers_and_symmetry_number(self):
"""
This method uses the symmetry package from RMG's QM module
and returns a tuple where the first element is the number
of optical isomers and the second element is the symmetry number.
"""
coordinates, atom_numbers, _ = self.loadGeometry()
unique_id = '0' # Just some name that the SYMMETRY code gives to one of its jobs
scr_dir = os.path.join(os.path.abspath('.'), str('scratch')) # Scratch directory that the SYMMETRY code writes its files in
if not os.path.exists(scr_dir):
os.makedirs(scr_dir)
try:
qmdata = QMData(
groundStateDegeneracy=1, # Only needed to check if valid QMData
numberOfAtoms=len(atom_numbers),
atomicNumbers=atom_numbers,
atomCoords=(coordinates, str('angstrom')),
energy=(0.0, str('kcal/mol')) # Only needed to avoid error
)
settings = type(str(''), (), dict(symmetryPath=str('symmetry'), scratchDirectory=scr_dir))() # Creates anonymous class
pgc = PointGroupCalculator(settings, unique_id, qmdata)
pg = pgc.calculate()
if pg is not None:
optical_isomers = 2 if pg.chiral else 1
symmetry = pg.symmetryNumber
logging.debug("Symmetry algorithm found {0} optical isomers and a symmetry number of {1}".format(optical_isomers,symmetry))
else:
logging.error("Symmetry algorithm errored when computing point group\nfor log file located at{0}.\nManually provide values in Arkane input.".format(self.path))
return optical_isomers, symmetry
finally:
shutil.rmtree(scr_dir)
17 changes: 9 additions & 8 deletions arkane/molpro.py
Expand Up @@ -37,11 +37,11 @@
from rmgpy.statmech import IdealGasTranslation, NonlinearRotor, LinearRotor, HarmonicOscillator, Conformer

from arkane.common import get_element_mass

from arkane.log import Log
################################################################################


class MolproLog:
class MolproLog(Log):
"""
Represents a Molpro log file. The attribute `path` refers to the
location on disk of the Molpro log file of interest. Methods are provided
Expand Down Expand Up @@ -167,7 +167,7 @@ def loadGeometry(self):

return coord, number, mass

def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, symfromlog=None, label=''):
def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=None, label=''):
"""
Load the molecular degree of freedom data from a log file created as
the result of a MolPro "Freq" quantum chemistry calculation with the thermo printed.
Expand All @@ -176,7 +176,12 @@ def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, sym
modes = []
unscaled_frequencies = []
E0 = 0.0

if opticalIsomers is None or symmetry is None:
_opticalIsomers, _symmetry = self.get_optical_isomers_and_symmetry_number()
if opticalIsomers is None:
opticalIsomers = _opticalIsomers
if symmetry is None:
symmetry = _symmetry
f = open(self.path, 'r')
line = f.readline()
while line != '':
Expand Down Expand Up @@ -224,10 +229,6 @@ def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, sym
mass = float(line.split()[2])
translation = IdealGasTranslation(mass=(mass,"amu"))
modes.append(translation)
# Read MolPro's estimate of the external symmetry number
elif 'Rotational Symmetry factor' in line and symmetry is None:
if symfromlog is True:
symmetry = int(float(line.split()[3]))

# Read moments of inertia for external rotational modes
elif 'Rotational Constants' in line and line.split()[-1]=='[GHz]':
Expand Down
2 changes: 1 addition & 1 deletion arkane/molproTest.py
Expand Up @@ -86,7 +86,7 @@ def testLoadHOSIFromMolpro_log(self):
"""

log = MolproLog(os.path.join(os.path.dirname(__file__),'data','HOSI_ccsd_t1.out'))
conformer, unscaled_frequencies = log.loadConformer(symfromlog=True, spinMultiplicity=1)
conformer, unscaled_frequencies = log.loadConformer(spinMultiplicity=1)
E0 = log.loadEnergy()

self.assertTrue(len([mode for mode in conformer.modes if isinstance(mode,IdealGasTranslation)]) == 1)
Expand Down
18 changes: 9 additions & 9 deletions arkane/qchem.py
Expand Up @@ -38,11 +38,12 @@
from rmgpy.statmech import IdealGasTranslation, NonlinearRotor, LinearRotor, HarmonicOscillator, Conformer

from arkane.common import check_conformer_energy, get_element_mass
from arkane.log import Log

################################################################################


class QChemLog:
class QChemLog(Log):
"""
Represent an output file from QChem. The attribute `path` refers to the
location on disk of the QChem output file of interest. Methods are provided
Expand Down Expand Up @@ -166,7 +167,7 @@ def loadGeometry(self):

return coord, number, mass

def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, symfromlog=None, label=''):
def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=None, label=''):
"""
Load the molecular degree of freedom data from an output file created as the result of a
QChem "Freq" calculation. As QChem's guess of the external symmetry number is not always correct,
Expand All @@ -176,6 +177,12 @@ def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, sym
modes = []; freq = []; mmass = []; rot = []; inertia = []
unscaled_frequencies = []
E0 = 0.0
if opticalIsomers is None or symmetry is None:
_opticalIsomers, _symmetry = self.get_optical_isomers_and_symmetry_number()
if opticalIsomers is None:
opticalIsomers = _opticalIsomers
if symmetry is None:
symmetry = _symmetry
f = open(self.path, 'r')
line = f.readline()
while line != '':
Expand Down Expand Up @@ -227,20 +234,13 @@ def loadConformer(self, symmetry=None, spinMultiplicity=0, opticalIsomers=1, sym
elif 'Eigenvalues --' in line:
inertia = [float(d) for d in line.split()[-3:]]

# Read QChem's estimate of the external rotational symmetry number, which may very well be incorrect
elif 'Rotational Symmetry Number is' in line and symfromlog:
symmetry = int(float(line.split()[-1]))
logging.debug('Rotational Symmetry read from QChem is {}'.format(str(symmetry)))

# Read the next line in the file
line = f.readline()

# Read the next line in the file
line = f.readline()

if len(inertia):
if symmetry is None:
symmetry = 1
if inertia[0] == 0.0:
# If the first eigenvalue is 0, the rotor is linear
inertia.remove(0.0)
Expand Down
12 changes: 6 additions & 6 deletions arkane/qchemTest.py
Expand Up @@ -74,11 +74,11 @@ def testLoadVibrationsFromQChemLog(self):
molecular energies can be properly read.
"""
log = QChemLog(os.path.join(os.path.dirname(__file__),'data','npropyl.out'))
conformer, unscaled_frequencies = log.loadConformer(symfromlog=True)
conformer, unscaled_frequencies = log.loadConformer()
self.assertEqual(len(conformer.modes[2]._frequencies.getValue()), 24)
self.assertEqual(conformer.modes[2]._frequencies.getValue()[5], 881.79)
log = QChemLog(os.path.join(os.path.dirname(__file__),'data','co.out'))
conformer, unscaled_frequencies = log.loadConformer(symfromlog=True)
conformer, unscaled_frequencies = log.loadConformer()
self.assertEqual(len(conformer.modes[2]._frequencies.getValue()), 1)
self.assertEqual(conformer.modes[2]._frequencies.getValue(), 2253.16)

Expand All @@ -88,7 +88,7 @@ def testLoadNpropylModesFromQChemLog(self):
molecular modes can be properly read.
"""
log = QChemLog(os.path.join(os.path.dirname(__file__),'data','npropyl.out'))
conformer, unscaled_frequencies = log.loadConformer(symfromlog=True)
conformer, unscaled_frequencies = log.loadConformer()

self.assertTrue(len([mode for mode in conformer.modes if isinstance(mode,IdealGasTranslation)]) == 1)
self.assertTrue(len([mode for mode in conformer.modes if isinstance(mode,NonlinearRotor)]) == 1)
Expand All @@ -101,10 +101,10 @@ def testSpinMultiplicityFromQChemLog(self):
molecular degrees of freedom can be properly read.
"""
log = QChemLog(os.path.join(os.path.dirname(__file__),'data','npropyl.out'))
conformer, unscaled_frequencies = log.loadConformer(symfromlog=True)
conformer, unscaled_frequencies = log.loadConformer()
self.assertEqual(conformer.spinMultiplicity, 2)
log = QChemLog(os.path.join(os.path.dirname(__file__),'data','co.out'))
conformer, unscaled_frequencies = log.loadConformer(symfromlog=True)
conformer, unscaled_frequencies = log.loadConformer()
self.assertEqual(conformer.spinMultiplicity, 1)

def testLoadCOModesFromQChemLog(self):
Expand All @@ -113,7 +113,7 @@ def testLoadCOModesFromQChemLog(self):
molecular degrees of freedom can be properly read.
"""
log = QChemLog(os.path.join(os.path.dirname(__file__),'data','co.out'))
conformer, unscaled_frequencies = log.loadConformer(symfromlog=True)
conformer, unscaled_frequencies = log.loadConformer()
E0 = log.loadEnergy()

self.assertTrue(len([mode for mode in conformer.modes if isinstance(mode,IdealGasTranslation)]) == 1)
Expand Down

0 comments on commit 4b1c083

Please sign in to comment.