Skip to content

Commit

Permalink
Merge pull request #134 from ReactionMechanismGenerator/model_chemistry
Browse files Browse the repository at this point in the history 
Modifications to how the model chemistry is treated to be consistent with Arkane
  • Loading branch information
@alongd
alongd committed Jun 9, 2019
2 parents f893b48 + d90e195 commit faed1a0
Show file tree
Hide file tree
Showing 13 changed files with 166 additions and 296 deletions.
6 changes: 2 additions & 4 deletions arc/job/inputs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,15 +116,13 @@
'arkane_input_species': """#!/usr/bin/env python
# -*- coding: utf-8 -*-
linear = {linear}{bonds}
externalSymmetry = {symmetry}
{bonds}externalSymmetry = {symmetry}
spinMultiplicity = {multiplicity}
opticalIsomers = {optical}
energy = {{'{model_chemistry}': Log('{sp_path}')}}
energy = {{'{sp_level}': Log('{sp_path}')}}
geometry = Log('{opt_path}')
Expand Down
99 changes: 47 additions & 52 deletions arc/main.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,8 +67,8 @@ class ARC(object):
`output` ``dict`` Output dictionary with status and final QM file paths for all species
Only used for restarting, the actual object used is in the Scheduler class
`use_bac` ``bool`` Whether or not to use bond additivity corrections for thermo calculations
`model_chemistry` ``list`` The model chemistry in Arkane for energy corrections (AE, BAC).
This can be usually determined automatically.
`model_chemistry` ``str`` The model chemistry in Arkane for energy corrections (AE, BAC)
and frequencies/ZPE scaling factor. Can usually be determined automatically.
`ess_settings` ``dict`` A dictionary of available ESS (keys) and a corresponding server list (values)
`initial_trsh` ``dict`` Troubleshooting methods to try by default. Keys are ESS software, values are trshs
't0' ``float`` Initial time when the project was spawned
Expand Down Expand Up @@ -772,66 +772,61 @@ def log_footer(self, level=logging.INFO):
logging.log(level, 'ARC execution terminated on {0}'.format(time.asctime()))

def determine_model_chemistry(self):
"""Determine the model_chemistry used in Arkane"""
"""Determine the model_chemistry to be used in Arkane
Todo:
* Determine whether the model chemistry exists in Arkane automaticaly instead of hard coding
"""
if self.model_chemistry:
self.model_chemistry = self.model_chemistry.lower()
if self.model_chemistry not in ['cbs-qb3', 'cbs-qb3-paraskevas', 'ccsd(t)-f12/cc-pvdz-f12',
'ccsd(t)-f12/cc-pvtz-f12', 'ccsd(t)-f12/cc-pvqz-f12',
'b3lyp/cbsb7', 'b3lyp/6-311g(2d,d,p)', 'b3lyp/6-311+g(3df,2p)',
'b3lyp/6-31g**']:
if self.model_chemistry.split('//')[0] not in [
'cbs-qb3', 'cbs-qb3-paraskevas', 'ccsd(t)-f12/cc-pvdz-f12', 'ccsd(t)-f12/cc-pvtz-f12',
'ccsd(t)-f12/cc-pvqz-f12', 'b3lyp/cbsb7', 'b3lyp/6-311g(2d,d,p)', 'b3lyp/6-311+g(3df,2p)',
'b3lyp/6-31g**']:
logging.warn('No bond additivity corrections (BAC) are available in Arkane for "model chemistry"'
' {0}. As a result, thermodynamic parameters are expected to be inaccurate. Make sure that'
' atom energy corrections (AEC) were supplied or are available in Arkane to avoid'
' error.'.format(self.model_chemistry))
else:
# model chemistry was not given, try to determine it from the sp_level
# model chemistry was not given, try to determine it from the sp_level and freq_level
model_chemistry = ''
if not self.composite_method:
sp_level = self.sp_level.lower()
else:
sp_level = self.composite_method
sp_level = sp_level.replace('f12a', 'f12').replace('f12b', 'f12')
if sp_level in ['ccsd(t)-f12/cc-pvdz', 'ccsd(t)-f12/cc-pvtz', 'ccsd(t)-f12/cc-pvqz']:
logging.warning('Using model chemistry {0} based on sp level {1}.'.format(
sp_level + '-f12', sp_level))
model_chemistry = sp_level + '-f12'
elif not model_chemistry and sp_level in ['cbs-qb3', 'cbs-qb3-paraskevas', 'ccsd(t)-f12/cc-pvdz-f12',
'ccsd(t)-f12/cc-pvtz-f12', 'ccsd(t)-f12/cc-pvqz-f12',
'b3lyp/cbsb7', 'b3lyp/6-311g(2d,d,p)', 'b3lyp/6-311+g(3df,2p)',
'b3lyp/6-31g**']:
model_chemistry = sp_level
elif self.use_bac:
logging.info('\n\n')
logging.warning('Could not determine appropriate Model Chemistry to be used in Arkane for'
' thermochemical parameter calculations. Not using atom energy corrections and '
'bond additivity corrections!\n\n')
self.use_bac = False
if self.composite_method:
self.model_chemistry = self.composite_method.lower()
else:
# use_bac is False, and no model chemistry was specified
if sp_level in ['m06-2x/cc-pvtz', 'g3', 'm08so/mg3s*', 'klip_1', 'klip_2', 'klip_3', 'klip_2_cc',
'ccsd(t)-f12/cc-pvdz-f12_h-tz', 'ccsd(t)-f12/cc-pvdz-f12_h-qz',
'ccsd(t)-f12/cc-pvdz-f12', 'ccsd(t)-f12/cc-pvtz-f12', 'ccsd(t)-f12/cc-pvqz-f12',
'ccsd(t)-f12/cc-pcvdz-f12', 'ccsd(t)-f12/cc-pcvtz-f12', 'ccsd(t)-f12/cc-pcvqz-f12',
'ccsd(t)-f12/cc-pvtz-f12(-pp)', 'ccsd(t)/aug-cc-pvtz(-pp)', 'ccsd(t)-f12/aug-cc-pvdz',
'ccsd(t)-f12/aug-cc-pvtz', 'ccsd(t)-f12/aug-cc-pvqz', 'b-ccsd(t)-f12/cc-pvdz-f12',
'b-ccsd(t)-f12/cc-pvtz-f12', 'b-ccsd(t)-f12/cc-pvqz-f12', 'b-ccsd(t)-f12/cc-pcvdz-f12',
'b-ccsd(t)-f12/cc-pcvtz-f12', 'b-ccsd(t)-f12/cc-pcvqz-f12', 'b-ccsd(t)-f12/aug-cc-pvdz',
'b-ccsd(t)-f12/aug-cc-pvtz', 'b-ccsd(t)-f12/aug-cc-pvqz', 'mp2_rmp2_pvdz',
'mp2_rmp2_pvtz', 'mp2_rmp2_pvqz', 'ccsd-f12/cc-pvdz-f12',
'ccsd(t)-f12/cc-pvdz-f12_noscale', 'g03_pbepbe_6-311++g_d_p', 'fci/cc-pvdz',
'fci/cc-pvtz', 'fci/cc-pvqz', 'bmk/cbsb7', 'bmk/6-311g(2d,d,p)', 'b3lyp/6-31g**',
'b3lyp/6-311+g(3df,2p)', 'MRCI+Davidson/aug-cc-pV(T+d)Z']:
model_chemistry = sp_level
self.model_chemistry = model_chemistry
logging.debug('Using {0} as model chemistry for energy corrections in Arkane'.format(
self.model_chemistry))
if not self.model_chemistry:
logging.warn('Could not determine a Model Chemistry to be used in Arkane, NOT calculating thermodata')
for spc in self.arc_species_list:
spc.generate_thermo = False
sp_level = self.sp_level.replace('f12a', 'f12').replace('f12b', 'f12').lower()
freq_level = self.freq_level.replace('f12a', 'f12').replace('f12b', 'f12').lower()
if sp_level in ['ccsd(t)-f12/cc-pvdz', 'ccsd(t)-f12/cc-pvtz', 'ccsd(t)-f12/cc-pvqz']:
logging.warning('Using model chemistry {0} based on sp level {1}.'.format(
sp_level + '-f12', sp_level))
sp_level += '-f12'
if sp_level not in ['ccsd(t)-f12/cc-pvdz-f12', 'ccsd(t)-f12/cc-pvtz-f12', 'ccsd(t)-f12/cc-pvqz-f12',
'b3lyp/cbsb7', 'b3lyp/6-311g(2d,d,p)', 'b3lyp/6-311+g(3df,2p)', 'b3lyp/6-31g**']\
and self.use_bac:
logging.info('\n\n')
logging.warning('Could not determine appropriate Model Chemistry to be used in Arkane for '
'thermochemical parameter calculations. Not using atom energy corrections and '
'bond additivity corrections!\n\n')
self.use_bac = False
elif sp_level not in ['m06-2x/cc-pvtz', 'g3', 'm08so/mg3s*', 'klip_1', 'klip_2', 'klip_3', 'klip_2_cc',
'ccsd(t)-f12/cc-pvdz-f12_h-tz', 'ccsd(t)-f12/cc-pvdz-f12_h-qz',
'ccsd(t)-f12/cc-pvdz-f12', 'ccsd(t)-f12/cc-pvtz-f12', 'ccsd(t)-f12/cc-pvqz-f12',
'ccsd(t)-f12/cc-pcvdz-f12', 'ccsd(t)-f12/cc-pcvtz-f12', 'ccsd(t)-f12/cc-pcvqz-f12',
'ccsd(t)-f12/cc-pvtz-f12(-pp)', 'ccsd(t)/aug-cc-pvtz(-pp)', 'ccsd(t)-f12/aug-cc-pvdz',
'ccsd(t)-f12/aug-cc-pvtz', 'ccsd(t)-f12/aug-cc-pvqz', 'b-ccsd(t)-f12/cc-pvdz-f12',
'b-ccsd(t)-f12/cc-pvtz-f12', 'b-ccsd(t)-f12/cc-pvqz-f12', 'b-ccsd(t)-f12/cc-pcvdz-f12',
'b-ccsd(t)-f12/cc-pcvtz-f12', 'b-ccsd(t)-f12/cc-pcvqz-f12', 'b-ccsd(t)-f12/aug-cc-pvdz',
'b-ccsd(t)-f12/aug-cc-pvtz', 'b-ccsd(t)-f12/aug-cc-pvqz', 'mp2_rmp2_pvdz',
'mp2_rmp2_pvtz', 'mp2_rmp2_pvqz', 'ccsd-f12/cc-pvdz-f12',
'ccsd(t)-f12/cc-pvdz-f12_noscale', 'g03_pbepbe_6-311++g_d_p', 'fci/cc-pvdz',
'fci/cc-pvtz', 'fci/cc-pvqz', 'bmk/cbsb7', 'bmk/6-311g(2d,d,p)', 'b3lyp/6-31g**',
'b3lyp/6-311+g(3df,2p)', 'MRCI+Davidson/aug-cc-pV(T+d)Z']:
logging.warn('Could not determine a Model Chemistry to be used in Arkane, '
'NOT calculating thermodata')
for spc in self.arc_species_list:
spc.generate_thermo = False
self.model_chemistry = sp_level + '//' + freq_level
if self.model_chemistry:
logging.info('Using {0} as model chemistry for energy corrections in Arkane'.format(
self.model_chemistry))
logging.info('Using {0} as a model chemistry in Arkane'.format(self.model_chemistry))

def determine_ess_settings(self, diagnostics=False):
"""
Expand Down
33 changes: 27 additions & 6 deletions arc/mainTest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,13 +53,13 @@ def test_as_dict(self):
arc_species_list=[spc1])
restart_dict = arc0.as_dict()
expected_dict = {'composite_method': '',
'conformer_level': 'b3lyp/6-31+g(d,p)',
'ts_guess_level': 'b3lyp/6-31+g(d,p)',
'conformer_level': 'b97d3/6-31+g(d,p)',
'ts_guess_level': 'b97d3/6-31+g(d,p)',
'opt_level': 'wb97xd/6-311++g(d,p)',
'freq_level': 'wb97xd/6-311++g(d,p)',
'initial_trsh': 'scf=(NDump=30)',
'max_job_time': 120,
'model_chemistry': 'ccsd(t)-f12/cc-pvtz-f12',
'model_chemistry': 'ccsd(t)-f12/cc-pvtz-f12//wb97xd/6-311++g(d,p)',
'output': {},
'project': 'arc_test',
'running_jobs': {},
Expand Down Expand Up @@ -266,9 +266,9 @@ def test_restart(self):
spc2 = Species().fromSMILES(str('CC([O])=O'))
spc2.generate_resonance_structures()
spc2.thermo = db.thermo.getThermoData(spc2)
self.assertAlmostEqual(spc2.getEnthalpy(298), -178003.44650359568, 1)
self.assertAlmostEqual(spc2.getEntropy(298), 283.5983103176096, 1)
self.assertAlmostEqual(spc2.getHeatCapacity(1000), 118.99753808225603, 1)
self.assertAlmostEqual(spc2.getEnthalpy(298), -176074.01886272896, 1)
self.assertAlmostEqual(spc2.getEntropy(298), 283.2225158405262, 1)
self.assertAlmostEqual(spc2.getHeatCapacity(1000), 118.28356605714401, 1)
self.assertTrue('arc_project_for_testing_delete_after_usage2' in spc2.thermo.comment)

# delete the generated library from RMG-database
Expand Down Expand Up @@ -332,6 +332,27 @@ def test_read_file(self):
with self.assertRaises(InputError):
read_file('nopath')

def test_determine_model_chemistry(self):
"""Test determining the model chemistry"""
arc0 = ARC(project='arc_model_chemistry_test_0', level_of_theory='CBS-QB3')
# arc0.determine_model_chemistry()
self.assertEqual(arc0.model_chemistry, 'cbs-qb3')

arc1 = ARC(project='arc_model_chemistry_test_1', model_chemistry='CBS-QB3-Paraskevas')
# arc1.determine_model_chemistry()
self.assertEqual(arc1.model_chemistry, 'cbs-qb3-paraskevas')

arc2 = ARC(project='arc_model_chemistry_test_2',
level_of_theory='ccsd(t)-f12/cc-pvtz-f12//wb97xd/6-311++g(d,p)')
# arc2.determine_model_chemistry()
self.assertEqual(arc2.model_chemistry, 'ccsd(t)-f12/cc-pvtz-f12//wb97xd/6-311++g(d,p)')

arc3 = ARC(project='arc_model_chemistry_test_2',
sp_level='ccsd(t)-f12/cc-pvtz-f12', opt_level='wb97xd/6-311++g(d,p)')
# arc3.determine_model_chemistry()
self.assertEqual(arc3.model_chemistry, 'ccsd(t)-f12/cc-pvtz-f12//wb97xd/6-311++g(d,p)')


@classmethod
def tearDownClass(cls):
"""
Expand Down
10 changes: 5 additions & 5 deletions arc/parser.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,19 +77,19 @@ def parse_t1(path):
return t1


def parse_e0(path):
def parse_e_elect(path, zpe_scale_factor=1.):
"""
Parse the zero K energy, E0, from an sp job output file
"""
if not os.path.isfile(path):
raise InputError('Could not find file {0}'.format(path))
log = determine_qm_software(fullpath=path)
try:
e0 = log.loadEnergy(frequencyScaleFactor=1.) * 0.001 # convert to kJ/mol
e_elect = log.loadEnergy(zpe_scale_factor) * 0.001 # convert to kJ/mol
except Exception:
logging.warning('Could not read E0 from {0}'.format(path))
e0 = None
return e0
logging.warning('Could not read e_elect from {0}'.format(path))
e_elect = None
return e_elect


def parse_xyz_from_file(path):
Expand Down
10 changes: 5 additions & 5 deletions arc/parserTest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,15 +75,15 @@ def test_parse_t1(self):
t1 = parser.parse_t1(path)
self.assertEqual(t1, 0.0086766)

def test_parse_e0(self):
def test_parse_e_elect(self):
"""Test parsing E0 from an sp job output file"""
path1 = os.path.join(arc_path, 'arc', 'testing', 'mehylamine_CCSD(T).out')
e0 = parser.parse_e0(path1)
self.assertEqual(e0, -251377.49160993524)
e_elect = parser.parse_e_elect(path1)
self.assertEqual(e_elect, -251377.49160993524)

path2 = os.path.join(arc_path, 'arc', 'testing', 'SO2OO_CBS-QB3.log')
e0 = parser.parse_e0(path2)
self.assertEqual(e0, -1833127.0939478774)
e_elect = parser.parse_e_elect(path2, zpe_scale_factor=0.99)
self.assertEqual(e_elect, -1833127.0939478774)

def test_parse_dipole_moment(self):
"""Test parsing the dipole moment from an opt job output file"""
Expand Down
2 changes: 1 addition & 1 deletion arc/plotter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@
def draw_3d(xyz=None, species=None, project_directory=None, save_only=False):
"""
Draws the molecule in a "3D-balls" style
If xyz ig given, it will be used, otherwise the function looks for species.final_xyz
If xyz is given, it will be used, otherwise the function looks for species.final_xyz
Input coordinates are in string format
Saves an image if a species and `project_directory` are provided
If `save_only` is ``True``, then don't plot, only save the image
Expand Down

0 comments on commit faed1a0

Please sign in to comment.