Skip to content

Commit

Permalink
Merge pull request #160 from ReactionMechanismGenerator/opt_isomorphism
Browse files Browse the repository at this point in the history 
 Check final_xyz isomorphism after opt and composite
  • Loading branch information
@alongd
alongd committed Aug 3, 2019
2 parents e39b3b7 + 6b7425a commit 82debb3
Show file tree
Hide file tree
Showing 3 changed files with 174 additions and 21 deletions.
92 changes: 73 additions & 19 deletions arc/scheduler.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1022,22 +1022,56 @@ def process_conformers(self, label):
logger.info('Only one conformer is available for species {0}, '
'using it as initial xyz'.format(label))
self.species_dict[label].initial_xyz = self.species_dict[label].conformers[0]
if not self.composite_method:
if self.job_types['opt']:
self.run_opt_job(label)
# check whether this conformer is isomorphic to the species 2D graph representation
# (since this won't be checked in determine_most_stable_conformer)
is_isomorphic, spawn_jobs = False, True
try:
b_mol = molecules_from_xyz(self.species_dict[label].initial_xyz,
multiplicity=self.species_dict[label].multiplicity,
charge=self.species_dict[label].charge)[1]
except SanitizationError:
b_mol = None
if b_mol is not None:
try:
is_isomorphic = check_isomorphism(self.species_dict[label].mol, b_mol)
except ValueError as e:
if self.species_dict[label].charge:
logger.error('Could not determine isomorphism for charged species {0}. Got the '
'following error:\n{1}'.format(label, e))
else:
logger.error('Could not determine isomorphism for (non-charged) species {0}. Got the '
'following error:\n{1}'.format(label, e))
if is_isomorphic:
logger.info('The only conformer for species {0} was found to be isomorphic '
'with the 2D graph representation {1}\n'.format(label, b_mol.toSMILES()))
self.output[label]['status'] += 'conformer passed isomorphism check; '
self.species_dict[label].conf_is_isomorphic = True
else:
if self.job_types['freq']:
self.run_freq_job(label)
if self.job_types['sp']:
self.run_sp_job(label)
if self.job_types['1d_rotors']:
self.run_scan_jobs(label)
if self.job_types['onedmin']:
self.run_onedmin_job(label)
if self.job_types['orbitals']:
self.run_orbitals_job(label)
else:
self.run_composite_job(label)
logger.error('The only conformer for species {0} is not isomorphic '
'with the 2D graph representation {1}\n'.format(label, b_mol.toSMILES()))
if self.allow_nonisomorphic_2d:
logger.info('Using this conformer anyway (allow_nonisomorphic_2d was set to True)')
else:
logger.info('Not using this conformer (to change this behavior, set allow_nonisomorphic_2d '
'to True)')
self.species_dict[label].conf_is_isomorphic, spawn_jobs = False, False
if spawn_jobs:
if not self.composite_method:
if self.job_types['opt']:
self.run_opt_job(label)
else:
if self.job_types['freq']:
self.run_freq_job(label)
if self.job_types['sp']:
self.run_sp_job(label)
if self.job_types['1d_rotors']:
self.run_scan_jobs(label)
if self.job_types['onedmin']:
self.run_onedmin_job(label)
if self.job_types['orbitals']:
self.run_orbitals_job(label)
else:
self.run_composite_job(label)

def parse_conformer(self, job, label, i):
"""
Expand Down Expand Up @@ -1125,7 +1159,9 @@ def determine_most_stable_conformer(self, label):
logger.info('Most stable conformer for species {0} was found to be isomorphic '
'with the 2D graph representation {1}\n'.format(label, b_mol.toSMILES()))
conformer_xyz = xyz
self.output[label]['status'] += 'passed isomorphism check; '
self.output[label]['status'] += 'most stable conformer ({0}) passed isomorphism ' \
'check; '.format(i)
self.species_dict[label].conf_is_isomorphic = True
else:
if energies[i] is not None:
logger.info('A conformer for species {0} was found to be isomorphic '
Expand All @@ -1139,11 +1175,12 @@ def determine_most_stable_conformer(self, label):
multiplicity=self.species_dict[label].multiplicity,
charge=self.species_dict[label].charge)[1]))
conformer_xyz = xyz
self.output[label]['status'] += 'passed isomorphism check but not for the most stable' \
' conformer; '
break
else:
if i == 0:
self.output[label]['status'] += 'most stable conformer ({0}) did not pass isomorphism ' \
'check; '.format(i)
self.species_dict[label].conf_is_isomorphic = False
logger.warning('Most stable conformer for species {0} with structure {1} was found to '
'be NON-isomorphic with the 2D graph representation {2}. Searching for '
'a different conformer that is isomorphic...'.format(
Expand Down Expand Up @@ -1281,7 +1318,16 @@ def parse_composite_geo(self, label, job):
if freq_ok:
# Update restart dictionary and save the yaml restart file:
self.save_restart_dict()
return True # run freq / scan jobs on this optimized geometry
success = True # run freq / scan jobs on this optimized geometry
if not self.species_dict[label].is_ts:
is_isomorphic = self.species_dict[label].check_final_xyz_isomorphism(
allow_nonisomorphic_2d=self.allow_nonisomorphic_2d)
if is_isomorphic:
self.output[label]['status'] += 'composite passed isomorphism check; '
else:
self.output[label]['status'] += 'composite did not pass isomorphism check; '
success &= is_isomorphic
return success
elif not self.species_dict[label].is_ts:
self.troubleshoot_negative_freq(label=label, job=job)
if job.job_status[1] != 'done' or not freq_ok:
Expand Down Expand Up @@ -1347,6 +1393,14 @@ def parse_opt_geo(self, label, job):
# for TSs, only use `draw_3d()`, not `show_sticks()` which gets connectivity wrong:
plotter.draw_structure(species=self.species_dict[label], project_directory=self.project_directory,
method='draw_3d')
if not self.species_dict[label].is_ts:
is_isomorphic = self.species_dict[label].check_final_xyz_isomorphism(
allow_nonisomorphic_2d=self.allow_nonisomorphic_2d)
if is_isomorphic:
self.output[label]['status'] += 'opt passed isomorphism check; '
else:
self.output[label]['status'] += 'opt did not pass isomorphism check; '
success &= is_isomorphic
else:
self.troubleshoot_opt_jobs(label=label)
if success:
Expand Down
56 changes: 55 additions & 1 deletion arc/species/species.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
from rmgpy.exceptions import InvalidAdjacencyListError

from arc.common import get_logger, get_atom_radius, determine_symmetry
from arc.arc_exceptions import SpeciesError, RotorError, InputError, TSError
from arc.arc_exceptions import SpeciesError, RotorError, InputError, TSError, SanitizationError
from arc.settings import default_ts_methods, valid_chars, minimum_barrier
from arc.parser import parse_xyz_from_file, parse_dipole_moment, parse_polarizability, process_conformers_file
from arc.species.converter import get_xyz_string, get_xyz_matrix, rdkit_conf_from_mol, standardize_xyz_string,\
Expand Down Expand Up @@ -172,6 +172,9 @@ class ARCSpecies(object):
with the fitted force field (recommended for drug-like species and species with many
heteroatoms). Another option is specifying 'cheap', and the "old" RDKit embedding method
will be used.
conf_is_isomorphic (bool): Whether the lowest conformer is isomorphic with the 2D graph representation
of the species. `True` if it is. Defaults to `None`. If `True`, an isomorphism check
will be strictly enforced for the final optimized coordinates.
"""
def __init__(self, label=None, is_ts=False, rmg_species=None, mol=None, xyz=None, multiplicity=None, charge=None,
Expand Down Expand Up @@ -218,6 +221,7 @@ def __init__(self, label=None, is_ts=False, rmg_species=None, mol=None, xyz=None
self.e_elect = None
self.arkane_file = None
self.svpfit_output_file = svpfit_output_file
self.conf_is_isomorphic = None
if self.is_ts:
if ts_methods is None:
self.ts_methods = default_ts_methods
Expand Down Expand Up @@ -450,6 +454,8 @@ def as_dict(self):
species_dict['external_symmetry'] = self.external_symmetry
species_dict['optical_isomers'] = self.optical_isomers
species_dict['neg_freqs_trshed'] = self.neg_freqs_trshed
if self.conf_is_isomorphic is not None:
species_dict['conf_is_isomorphic'] = self.conf_is_isomorphic
if self.bond_corrections is not None:
species_dict['bond_corrections'] = self.bond_corrections
if self.mol is not None:
Expand Down Expand Up @@ -502,6 +508,7 @@ def from_dict(self, species_dict):
self.initial_xyz = standardize_xyz_string(species_dict['initial_xyz']) if 'initial_xyz' in species_dict\
else None
self.final_xyz = standardize_xyz_string(species_dict['final_xyz']) if 'final_xyz' in species_dict else None
self.conf_is_isomorphic = species_dict['conf_is_isomorphic'] if 'conf_is_isomorphic' in species_dict else None
self.is_ts = species_dict['is_ts'] if 'is_ts' in species_dict else False
if self.is_ts:
self.ts_conf_spawned = species_dict['ts_conf_spawned'] if 'ts_conf_spawned' in species_dict else False
Expand Down Expand Up @@ -986,6 +993,53 @@ def set_transport_data(self, lj_path, opt_path, bath_gas, opt_level, freq_path='
comment=str(comment)
)

def check_final_xyz_isomorphism(self, allow_nonisomorphic_2d=False):
"""
Check whether the perception of self.final_xyz is isomorphic with self.mol.
Args:
allow_nonisomorphic_2d (bool): Whether to continue spawning jobs for the species even if this test fails.
`True` to allow (default is `False`).
Returns:
bool: Whether the perception of self.final_xyz is isomorphic with self.mol, `True` if it is.
"""
passed_test, return_value = False, False
if self.mol is not None:
try:
b_mol = molecules_from_xyz(self.final_xyz, multiplicity=self.multiplicity, charge=self.charge)[1]
except SanitizationError:
b_mol = None
if b_mol is not None:
is_isomorphic = check_isomorphism(self.mol, b_mol)
else:
is_isomorphic = False
if is_isomorphic:
passed_test, return_value = True, True
else:
# isomorphism test failed
passed_test = False
if self.conf_is_isomorphic:
if allow_nonisomorphic_2d:
# conformer was isomorphic, we **do** allow nonisomorphism, and the optimized structure isn't
return_value = True
else:
# conformer was isomorphic, we don't allow nonisomorphism, but the optimized structure isn't
return_value = False
else:
# conformer was not isomorphic, don't strictly enforce isomorphism here
return_value = True
if not passed_test:
logger.error('The optimized geometry of species {0} is not isomorphic with the 2D structure {1}'.format(
self.label, self.mol.toSMILES()))
if not return_value:
logger.error('Not spawning additional jobs for this species!')
else:
logger.info('Species {0} was found to be isomorphic with the perception of its optimized coordinates.')
else:
logger.error('Cannot check isomorphism for species {0}'.format(self.label))
return return_value


class TSGuess(object):
"""
Expand Down
47 changes: 46 additions & 1 deletion arc/species/speciesTest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ def test_conformers(self):

self.spc4.conformers = list()
self.spc4.generate_conformers()
self.assertIn(len(self.spc4.conformers), [4, 5])
self.assertIn(len(self.spc4.conformers), [3, 4, 5])

self.spc5.conformers = list()
self.spc5.generate_conformers()
Expand Down Expand Up @@ -839,6 +839,51 @@ def test_check_xyz(self):
self.assertFalse(check_xyz(xyz1, multiplicity=2, charge=1))
self.assertTrue(check_xyz(xyz1, multiplicity=1, charge=-1))

def test_check_final_xyz_isomorphism(self):
"""Test the check_final_xyz_isomorphism() method"""
xyz1 = """C -1.9681540 0.0333440 -0.0059220
C -0.6684360 -0.7562450 0.0092140
C 0.5595480 0.1456260 -0.0036480
O 0.4958540 1.3585920 0.0068500
N 1.7440770 -0.5331650 -0.0224050
H -2.8220500 -0.6418490 0.0045680
H -2.0324190 0.6893210 0.8584580
H -2.0300690 0.6574940 -0.8939330
H -0.6121640 -1.4252590 -0.8518010
H -0.6152180 -1.3931710 0.8949150
H 1.7901420 -1.5328370 0.0516350
H 2.6086580 -0.0266360 0.0403330"""
spc1 = ARCSpecies(label='propanamide1', smiles='CCC(=O)N', xyz=xyz1)
spc1.final_xyz = spc1.conformers[0]
is_isomorphic1 = spc1.check_final_xyz_isomorphism()
self.assertTrue(is_isomorphic1)

xyz2 = """C 0.6937910 -0.8316510 0.0000000
O 0.2043990 -1.9431180 0.0000000
N 0.0000000 0.3473430 0.0000000
C -1.4529360 0.3420060 0.0000000
C 0.6724520 1.6302410 0.0000000
H 1.7967050 -0.6569720 0.0000000
H -1.7909050 -0.7040620 0.0000000
H -1.8465330 0.8552430 0.8979440
H -1.8465330 0.8552430 -0.8979440
H 1.7641260 1.4761740 0.0000000
H 0.4040540 2.2221690 -0.8962980
H 0.4040540 2.2221690 0.8962980""" # dimethylformamide, O=CN(C)C
spc2 = ARCSpecies(label='propanamide2', smiles='CCC(=O)N', xyz=xyz1) # define w/ the correct xyz
spc2.final_xyz = xyz2 # set .final_xyz to the incorrect isomer

spc2.conf_is_isomorphic = True # set to True so that isomorphism is strictly enforced
is_isomorphic2 = spc2.check_final_xyz_isomorphism()
self.assertFalse(is_isomorphic2)

is_isomorphic3 = spc2.check_final_xyz_isomorphism(allow_nonisomorphic_2d=True)
self.assertTrue(is_isomorphic3)

spc2.conf_is_isomorphic = False # set to False so that isomorphism is not strictly enforced
is_isomorphic4 = spc2.check_final_xyz_isomorphism()
self.assertTrue(is_isomorphic4)

@classmethod
def tearDownClass(cls):
"""
Expand Down

0 comments on commit 82debb3

Please sign in to comment.