Skip to content

Commit

Permalink
Merge pull request #104 from ReactionMechanismGenerator/confs
Browse files Browse the repository at this point in the history 
Allow ARCSpecies to be initialized with a conformer list file
  • Loading branch information
@alongd
alongd committed Mar 28, 2019
2 parents 8b5db7a + 1a8b806 commit f96500f
Show file tree
Hide file tree
Showing 9 changed files with 215 additions and 77 deletions.
5 changes: 4 additions & 1 deletion arc/job/job.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -234,7 +234,10 @@ def __init__(self, project, settings, species_name, xyz, job_type, level_of_theo
self.submit = ''
self.input = ''
self.server_nodes = list()
self._write_initiated_job_to_csv_file()
if job_num is None:
# this checks jon_num and not self.job_num on purpose
# if job_num was given, then don't save as initiated jobs, this is a restarted job
self._write_initiated_job_to_csv_file()

def as_dict(self):
"""A helper function for dumping this object as a dictionary in a YAML file for restarting ARC"""
Expand Down
6 changes: 3 additions & 3 deletions arc/job/submit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,12 +31,12 @@
SubmitDir=`pwd`
GAUSS_SCRDIR=/scratch/users/{un}/g09/$SLURM_JOB_NAME-$SLURM_JOB_ID
export GAUSS_SCRDIR
export GAUSS_SCRDIR
mkdir -p $GAUSS_SCRDIR
mkdir -p $WorkDir
cd $WorkDir
cd $WorkDir
. $g09root/g09/bsd/g09.profile
cp $SubmitDir/input.gjf .
Expand Down Expand Up @@ -79,7 +79,7 @@
SubmitDir=`pwd`
mkdir -p $WorkDir
cd $WorkDir
cd $WorkDir
cp $SubmitDir/input.inp .
Expand Down
3 changes: 3 additions & 0 deletions arc/main.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

from __future__ import (absolute_import, division, print_function, unicode_literals)
import logging
import warnings
import sys
import os
import time
Expand Down Expand Up @@ -776,6 +777,8 @@ def log_header(self, level=logging.INFO):
logging.log(level, 'The current git HEAD for ARC is:')
logging.log(level, ' {0}\n {1}\n'.format(head, date))
logging.info('Starting project {0}'.format(self.project))
# ignore Paramiko warnings:
warnings.filterwarnings(action='ignore', module='.*paramiko.*')

def log_footer(self, level=logging.INFO):
"""
Expand Down
1 change: 0 additions & 1 deletion arc/mainTest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,6 @@ def test_as_dict(self):
'multiplicity': 1,
'neg_freqs_trshed': [],
'number_of_rotors': 0,
'opt_level': '',
'rotors_dict': {},
't1': None}],
}
Expand Down
130 changes: 78 additions & 52 deletions arc/scheduler.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -308,7 +308,7 @@ def __init__(self, project, settings, species_list, composite_method, conformer_
# restart-related check are performed in run_scan_jobs()
self.run_scan_jobs(species.label)
elif not self.species_dict[species.label].is_ts and self.generate_conformers\
and 'geo' not in self.output[species.label]:
and 'geo' not in self.output[species.label] and not species.conformers:
self.species_dict[species.label].generate_conformers()
else:
# Species is loaded from a YAML file
Expand All @@ -325,6 +325,13 @@ def schedule_jobs(self):
"""
The main job scheduling block
"""
for species in self.species_dict.values():
if not species.initial_xyz and not species.final_xyz and species.conformers and species.conformer_energies:
self.determine_most_stable_conformer(species.label)
if self.composite_method:
self.run_composite_job(species.label)
else:
self.run_opt_job(species.label)
if self.generate_conformers:
self.run_conformer_jobs()
while self.running_jobs != {}: # loop while jobs are still running
Expand Down Expand Up @@ -442,7 +449,8 @@ def schedule_jobs(self):
folder_name = 'rxns' if self.species_dict[label].is_ts else 'Species'
orbitals_path = os.path.join(self.project_directory, 'output', folder_name, label, 'geometry',
'orbitals.fchk')
shutil.copyfile(job.local_path_to_orbitals_file, orbitals_path)
if os.path.isfile(job.local_path_to_orbitals_file):
shutil.copyfile(job.local_path_to_orbitals_file, orbitals_path)
self.timer = False
break

Expand Down Expand Up @@ -549,15 +557,8 @@ def run_conformer_jobs(self):
"""
for label in self.unique_species_labels:
if not self.species_dict[label].is_ts and 'opt converged' not in self.output[label]['status']\
and 'opt' not in self.job_dict[label]:
geo_dir = os.path.join(self.project_directory, 'output', 'Species', label, 'geometry')
if not os.path.exists(geo_dir):
os.makedirs(geo_dir)
conf_path = os.path.join(geo_dir, 'conformers_before_optimization.txt')
with open(conf_path, 'w') as f:
for conf in self.species_dict[label].conformers:
f.write(conf)
f.write('\n\n')
and 'opt' not in self.job_dict[label] and not self.species_dict[label].conformer_energies:
self.save_conformers_file(label)
if not self.testing:
if len(self.species_dict[label].conformers) > 1:
self.job_dict[label]['conformers'] = dict()
Expand Down Expand Up @@ -714,13 +715,13 @@ def run_orbitals_job(self, label):

def parse_conformer_energy(self, job, label, i):
"""
Parse E0 (Hartree) from the conformer opt output file, and save it in the 'conformer_energies' attribute.
Parse E0 (J/mol) from the conformer opt output file, and save it in the 'conformer_energies' attribute.
"""
if job.job_status[1] == 'done':
log = Log(path='')
log.determine_qm_software(fullpath=job.local_path_to_output_file)
e0 = log.software_log.loadEnergy()
self.species_dict[label].conformer_energies[i] = e0 # in J/mol
self.species_dict[label].conformer_energies[i] = log.software_log.loadEnergy() # in J/mol
logging.debug('energy for {0} is {1}'.format(i, self.species_dict[label].conformer_energies[i]))
else:
logging.warn('Conformer {i} for {label} did not converge!'.format(i=i, label=label))

Expand All @@ -733,7 +734,7 @@ def determine_most_stable_conformer(self, label):
if self.species_dict[label].is_ts:
raise SchedulerError('The determine_most_stable_conformer() method does not deal with transition'
' state guesses.')
if all(e == 0.0 for e in self.species_dict[label].conformer_energies):
if 'conformers' in self.job_dict[label] and all(e is None for e in self.species_dict[label].conformer_energies):
logging.error('No conformer converged for species {0}! Trying to troubleshoot conformer jobs...'.format(
label))
for i, job in self.job_dict[label]['conformers'].items():
Expand All @@ -743,34 +744,19 @@ def determine_most_stable_conformer(self, label):
conformer_xyz = None
xyzs = list()
log = Log(path='')
for job in self.job_dict[label]['conformers'].values():
log.determine_qm_software(fullpath=job.local_path_to_output_file)
try:
coord, number, _ = log.software_log.loadGeometry()
except RMGInputError:
xyzs.append(None)
else:
xyzs.append(get_xyz_string(xyz=coord, number=number))
energies, xyzs = (list(t) for t in zip(*sorted(zip(self.species_dict[label].conformer_energies, xyzs))))
smiles_list = list()
for xyz in xyzs:
b_mol = molecules_from_xyz(xyz, multiplicity=self.species_dict[label].multiplicity,
charge=self.species_dict[label].charge)[1]
smiles = b_mol.toSMILES() if b_mol is not None else 'no 2D structure'
smiles_list.append(smiles)
geo_dir = os.path.join(self.project_directory, 'output', 'Species', label, 'geometry')
if not os.path.exists(geo_dir):
os.makedirs(geo_dir)
conf_path = os.path.join(geo_dir, 'conformers_after_optimization.txt')
with open(conf_path, 'w') as f:
for i, xyz in enumerate(xyzs):
f.write('conformer {0}:\n'.format(i))
if xyz is not None:
f.write(xyz + '\n')
f.write('SMILES: ' + smiles_list[i] + '\n')
f.write('Relative Energy: {0} kJ/mol\n\n\n'.format((energies[i] - min(energies)) * 0.001))
if self.species_dict[label].conformer_energies:
xyzs = self.species_dict[label].conformers
else:
for job in self.job_dict[label]['conformers'].values():
log.determine_qm_software(fullpath=job.local_path_to_output_file)
try:
coord, number, _ = log.software_log.loadGeometry()
except RMGInputError:
xyzs.append(None)
else:
f.write('Failed to converge')
xyzs.append(get_xyz_string(xyz=coord, number=number))
energies, xyzs = (list(t) for t in zip(*sorted(zip(self.species_dict[label].conformer_energies, xyzs))))
self.save_conformers_file(label, xyzs=xyzs, energies=energies)
# Run isomorphism checks if a 2D representation is available
if self.species_dict[label].mol is not None:
for i, xyz in enumerate(xyzs):
Expand All @@ -781,11 +767,14 @@ def determine_most_stable_conformer(self, label):
is_isomorphic = check_isomorphism(self.species_dict[label].mol, b_mol)
except ValueError as e:
if self.species_dict[label].charge:
logging.error('Could not determine isomorphism for charged species. Got the '
'following error:\n{0}'.format(e.message))
logging.error('Could not determine isomorphism for charged species {0}. '
'Optimizing the most stable conformer anyway. Got the '
'following error:\n{1}'.format(label, e))
else:
logging.error('Could not determine isomorphism for (non-charged) species. Got the '
'following error:\n{0}'.format(e.message))
logging.error('Could not determine isomorphism for (non-charged) species {0}. '
'Optimizing the most stable conformer anyway. Got the '
'following error:\n{1}'.format(label, e))
conformer_xyz = xyzs[0]
break
if is_isomorphic:
if i == 0:
Expand Down Expand Up @@ -818,7 +807,7 @@ def determine_most_stable_conformer(self, label):
smiles_list.append(molecules_from_xyz(xyz, multiplicity=self.species_dict[label].multiplicity,
charge=self.species_dict[label].charge)[1])
if self.allow_nonisomorphic_2d or self.species_dict[label].charge:
# we'll optimize the most stable conformer even if it not isomorphic to the 2D graph
# we'll optimize the most stable conformer even if it is not isomorphic to the 2D graph
logging.error('No conformer for {0} was found to be isomorphic with the 2D graph representation'
' {1} (got: {2}). Optimizing the most stable conformer anyway.'.format(
label, self.species_dict[label].mol.toSMILES(), smiles_list))
Expand Down Expand Up @@ -846,7 +835,7 @@ def determine_most_likely_ts_conformer(self, label):
if not self.species_dict[label].is_ts:
raise SchedulerError('The determine_most_likely_ts_conformer() method only deals with transition'
' state guesses.')
if all(e == 0.0 for e in self.species_dict[label].conformer_energies):
if all(e is None for e in self.species_dict[label].conformer_energies):
logging.error('No guess converged for TS {0}!')
# for i, job in self.job_dict[label]['conformers'].items():
# self.troubleshoot_ess(label, job, level_of_theory=job.level_of_theory, job_type='conformer',
Expand Down Expand Up @@ -893,6 +882,7 @@ def parse_composite_geo(self, label, job):
self.species_dict[label].final_xyz = get_xyz_string(xyz=coord, number=number)
self.output[label]['status'] += 'composite converged; '
self.output[label]['composite'] = os.path.join(job.local_path, 'output.out')
self.species_dict[label].opt_level = self.composite_method
rxn_str = ''
if self.species_dict[label].is_ts:
rxn_str = ' of reaction {0}'.format(self.species_dict[label].rxn_label)
Expand All @@ -914,7 +904,6 @@ def parse_composite_geo(self, label, job):
return True # run freq / scan jobs on this optimized geometry
elif not self.species_dict[label].is_ts:
self.troubleshoot_negative_freq(label=label, job=job)
self.species_dict[label].opt_level = self.composite_method
if job.job_status[1] != 'done' or not freq_ok:
self.troubleshoot_ess(label=label, job=job, level_of_theory=job.level_of_theory, job_type='composite')
return False # return ``False``, so no freq / scan jobs are initiated for this unoptimized geometry
Expand Down Expand Up @@ -1314,7 +1303,7 @@ def troubleshoot_negative_freq(self, label, job):
xyz2 = atomcoords - factor * displacement
self.species_dict[label].conformers.append(get_xyz_string(xyz=xyz1, number=atomnos))
self.species_dict[label].conformers.append(get_xyz_string(xyz=xyz2, number=atomnos))
self.species_dict[label].conformer_energies.extend([0.0, 0.0]) # a placeholder (lists are synced)
self.species_dict[label].conformer_energies.extend([None, None]) # a placeholder (lists are synced)
self.job_dict[label]['conformers'] = dict() # initialize the conformer job dictionary
for i, xyz in enumerate(self.species_dict[label].conformers):
self.run_job(label=label, xyz=xyz, level_of_theory=self.conformer_level, job_type='conformer', conformer=i)
Expand Down Expand Up @@ -1736,16 +1725,53 @@ def make_reaction_labels_info_file(self):
f.write('Reaction labels and respective TS labels:\n\n')
return rxn_info_path

def save_conformers_file(self, label, xyzs=None, energies=None):
"""
A helper function for saving the conformers for species `label` before and after optimization
If `xyzs` is given, then it is used as the conformers xyz list, otherwise it is taken from the species.conformer
attribute. If `energies`, a list of respective conformer energies in J/mol, is given, it will also be reported.
"""
geo_dir = os.path.join(self.project_directory, 'output', 'Species', label, 'geometry')
if not os.path.exists(geo_dir):
os.makedirs(geo_dir)
smiles_list = list()
xyzs = xyzs or self.species_dict[label].conformers
for xyz in xyzs:
b_mol = molecules_from_xyz(xyz, multiplicity=self.species_dict[label].multiplicity,
charge=self.species_dict[label].charge)[1]
smiles = b_mol.toSMILES() if b_mol is not None else 'no 2D structure'
smiles_list.append(smiles)
if energies is not None:
conf_path = os.path.join(geo_dir, 'conformers_after_optimization.txt')
else:
conf_path = os.path.join(geo_dir, 'conformers_before_optimization.txt')
with open(conf_path, 'w') as f:
if energies is not None:
f.write('conformers optimized at {0}\n\n'.format(self.conformer_level))
for i, conf in enumerate(xyzs):
f.write('conformer {0}:\n'.format(i))
if conf is not None:
f.write(conf + '\n')
f.write('SMILES: ' + smiles_list[i] + '\n')
if energies is not None:
if energies[i] == min(energies):
f.write('Relative Energy: 0 kJ/mol (lowest)')
else:
f.write('Relative Energy: {0:.3f} kJ/mol'.format((energies[i] - min(energies)) * 0.001))
else:
f.write('Failed to converge')
f.write('\n\n\n')


# Add a custom string representer to use block literals for multiline strings
def string_representer(dumper, data):
"""Add a custom string representer to use block literals for multiline strings"""
if len(data.splitlines()) > 1:
return dumper.represent_scalar(tag='tag:yaml.org,2002:str', value=data, style='|')
return dumper.represent_scalar(tag='tag:yaml.org,2002:str', value=data)


# Add a custom unicode representer to use block literals for multiline strings
def unicode_representer(dumper, data):
"""Add a custom unicode representer to use block literals for multiline strings"""
if len(data.splitlines()) > 1:
return yaml.ScalarNode(tag='tag:yaml.org,2002:str', value=data, style='|')
return yaml.ScalarNode(tag='tag:yaml.org,2002:str', value=data)
17 changes: 11 additions & 6 deletions arc/schedulerTest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,10 +63,14 @@ def test_conformers(self):
self.sched1.job_dict[label]['conformers'] = dict()
self.sched1.job_dict[label]['conformers'][0] = self.job1
self.sched1.job_dict[label]['conformers'][1] = self.job2
self.sched1.species_dict[label].conformer_energies = [None, None]
self.sched1.species_dict[label].conformers = [None, None]
self.sched1.parse_conformer_energy(job=self.job1, label=label, i=0)
self.sched1.parse_conformer_energy(job=self.job2, label=label, i=1)
expecting = [-251596443.5088726, -254221943.3698632]
self.assertEqual(self.sched1.species_dict[label].conformer_energies, expecting)
self.sched1.species_dict[label].conformers[0] = parser.parse_xyz_from_file(self.job1.local_path_to_output_file)
self.sched1.species_dict[label].conformers[1] = parser.parse_xyz_from_file(self.job2.local_path_to_output_file)

self.sched1.determine_most_stable_conformer(label=label)
expecting = """N -0.75555952 -0.12937106 0.00000000
Expand All @@ -83,20 +87,21 @@ def test_conformers(self):
self.assertTrue(os.path.isfile(methylamine_conf_path))
with open(methylamine_conf_path, 'r') as f:
lines = f.readlines()
self.assertEqual(lines[0], 'conformer 0:\n')
self.assertEqual(lines[9], 'SMILES: CN\n')
self.assertTrue('Relative Energy:' in lines[10])
self.assertEqual(lines[14][0], 'N')
self.assertTrue('conformers optimized at' in lines[0])
self.assertEqual(lines[11], 'SMILES: CN\n')
self.assertTrue('Relative Energy:' in lines[12])
self.assertEqual(lines[16][0], 'N')

self.sched1.run_conformer_jobs()
self.sched1.save_conformers_file(label='C2H6')
c2h6_conf_path = os.path.join(self.sched1.project_directory, 'output', 'Species', 'C2H6', 'geometry',
'conformers_before_optimization.txt')
self.assertTrue(os.path.isfile(c2h6_conf_path))
with open(c2h6_conf_path, 'r') as f:
lines = f.readlines()
self.assertEqual(lines[0][0], 'C')
self.assertEqual(lines[1][0], 'C')
self.assertEqual(lines[9], '\n')
self.assertEqual(lines[12][0], 'H')
self.assertEqual(lines[17][0], 'H')

def test_check_negative_freq(self):
"""Test the check_negative_freq() method"""
Expand Down
Loading

0 comments on commit f96500f

Please sign in to comment.