Merge pull request #143 from ReactionMechanismGenerator/conformers

A new way of thinking about conformers
ReactionMechanismGenerator · Jul 25, 2019 · 8aa859b · 8aa859b
2 parents b6eb262 + e2e3abd
commit 8aa859b
Show file tree

Hide file tree

Showing 40 changed files with 10,036 additions and 950 deletions.
diff --git a/ARC.py b/ARC.py
@@ -10,7 +10,8 @@
 import os
 import logging
 
-from arc.main import ARC, read_file
+from arc.main import ARC
+from arc.common import read_yaml_file
 
 ################################################################################
 
@@ -48,10 +49,11 @@ def parse_command_line_arguments(command_line_args=None):
 
 
 def main():
+    """The main ARC executable function"""
     # Parse the command-line arguments (requires the argparse module)
     args = parse_command_line_arguments()
     input_file = args.file
-    input_dict = read_file(input_file)
+    input_dict = read_yaml_file(input_file)
     project_directory = os.path.abspath(os.path.dirname(args.file))
     try:
         input_dict['project']

diff --git a/arc/arc_exceptions.py b/arc/arc_exceptions.py
@@ -6,6 +6,13 @@
 """
 
 
+class ConformerError(Exception):
+    """
+    An exception raised when generating conformers
+    """
+    pass
+
+
 class InputError(Exception):
     """
     An exception raised when parsing an input file for any module

diff --git a/arc/common.py b/arc/common.py
@@ -16,6 +16,11 @@
 import shutil
 import subprocess
 import yaml
+import numpy as np
+
+from rmgpy.molecule.element import getElement
+from rmgpy.qm.qmdata import QMData
+from rmgpy.qm.symmetry import PointGroupCalculator
 
 from arc.settings import arc_path, servers
 from arc.arc_exceptions import InputError, SettingsError
@@ -27,23 +32,6 @@
 VERSION = '1.0.0'
 
 
-def read_file(path):
-    """
-    Read the ARC YAML input file and return the parameters in a dictionary.
-
-    Args:
-        (str, unicode): The input file path.
-
-    Returns:
-        dict: The input dictionary read from the file.
-    """
-    if not os.path.isfile(path):
-        raise InputError('Could not find the input file {0}'.format(path))
-    with open(path, 'r') as f:
-        input_dict = yaml.load(stream=f, Loader=yaml.FullLoader)
-    return input_dict
-
-
 def time_lapse(t0):
     """
     A helper function returning the elapsed time since t0.
@@ -94,7 +82,7 @@ def check_ess_settings(ess_settings=None):
                                 'strings. Got: {0} which is a {1}'.format(server_list, type(server_list)))
     # run checks:
     for ess, server_list in settings.items():
-        if ess.lower() not in ['gaussian', 'qchem', 'molpro', 'onedmin', 'orca']:
+        if ess.lower() not in ['gaussian', 'qchem', 'molpro', 'orca', 'terachem', 'onedmin', 'gromacs']:
             raise SettingsError('Recognized ESS software are Gaussian, QChem, Molpro, Orca or OneDMin. '
                                 'Got: {0}'.format(ess))
         for server in server_list:
@@ -254,14 +242,36 @@ def log_footer(execution_time, level=logging.INFO):
     logger.log(level, 'ARC execution terminated on {0}'.format(time.asctime()))
 
 
-def save_dict_file(path, restart_dict):
+def read_yaml_file(path):
+    """
+    Read a YAML file (usually an input / restart file, but also conformers file)
+    and return the parameters as python variables.
+
+    Args:
+        path (str, unicode): The YAML file path to read.
+
+    Returns:
+        dict, list: The content read from the file.
+    """
+    if not os.path.isfile(path):
+        raise InputError('Could not find the YAML file {0}'.format(path))
+    with open(path, 'r') as f:
+        content = yaml.load(stream=f, Loader=yaml.FullLoader)
+    return content
+
+
+def save_yaml_file(path, content):
     """
-    Save an input / restart YAML file
+    Save a YAML file (usually an input / restart file, but also conformers file)
+
+    Args:
+        path (str, unicode): The YAML file path to save.
+        content (str, unicode): The content to save.
     """
     yaml.add_representer(str, string_representer)
     yaml.add_representer(unicode, unicode_representer)
     logger.debug('Creating a restart file...')
-    content = yaml.dump(data=restart_dict, encoding='utf-8', allow_unicode=True)
+    content = yaml.dump(data=content, encoding='utf-8', allow_unicode=True)
     with open(path, 'w') as f:
         f.write(content)
 
@@ -278,3 +288,103 @@ def unicode_representer(dumper, data):
     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)
+
+
+def get_ordinal_indicator(number):
+    """
+    Returns the ordinal indicator for an integer.
+
+    Args:
+        number (int): An integer for which the ordinal indicator will be determined.
+
+    Returns:
+        str, unicode: The integer's ordinal indicator.
+    """
+    ordinal_dict = {1: 'st', 2: 'nd', 3: 'rd'}
+    if number > 13:
+        number %= 10
+    if number in ordinal_dict.keys():
+        return ordinal_dict[number]
+    return 'th'
+
+
+# ATOM_RADII data taken from DOI: 10.1039/b801115j
+ATOM_RADII = {'H': 0.31, 'He': 0.28,
+              'Li': 1.28, 'Be': 0.96, 'B': 0.84, 'C': 0.76, 'N': 0.71, 'O': 0.66, 'F': 0.57, 'Ne': 0.58,
+              'Na': 1.66, 'Mg': 1.41, 'Al': 1.21, 'Si': 1.11, 'P': 1.07, 'S': 1.05, 'Cl': 1.02, 'Ar': 1.06,
+              'K': 2.03, 'Ca': 1.76, 'Sc': 1.70, 'Ti': 1.60, 'V': 1.53, 'Cr': 1.39, 'Mn': 1.39, 'Fe': 1.32,
+              'Co': 1.26, 'Ni': 1.24, 'Cu': 1.32, 'Zn': 1.22, 'Ga': 1.22, 'Ge': 1.20, 'As': 1.19,
+              'Se': 1.20, 'Br': 1.20, 'Kr': 1.16, 'Rb': 2.20, 'Sr': 1.95, 'Y': 1.90, 'Zr': 1.75,
+              'Nb': 1.64, 'Mo': 1.54, 'Tc': 1.47, 'Ru': 1.46, 'Rh': 1.42, 'Pd': 1.39, 'Ag': 1.45,
+              'Cd': 1.44, 'In': 1.42, 'Sn': 1.39, 'Sb': 1.39, 'Te': 1.38, 'I': 1.39, 'Xe': 1.40,
+              'Cs': 2.44, 'Ba': 2.15, 'Pt': 1.36, 'Au': 1.36, 'Hg': 1.32, 'Tl': 1.45, 'Pb': 1.46,
+              'Bi': 1.48, 'Po': 1.40, 'At': 1.50, 'Rn': 1.50, 'Fr': 2.60, 'Ra': 2.21, 'U': 1.96}
+
+
+def get_atom_radius(symbol):
+    """
+    Get the atom covalent radius in Angstroms, data in the ATOM_RADII dict.
+    (Change to QCElemental after transitioning to Py3)
+
+    Args:
+        symbol (str, unicode): The atomic symbol.
+
+    Returns:
+        float: The atomic covalent radius (None if not found).
+    """
+    if not isinstance(symbol, (str, unicode)):
+        raise InputError('the symbol argument must be string, got {0} which is a {1}'.format(symbol, type(symbol)))
+
+    if symbol in ATOM_RADII:
+        return ATOM_RADII[symbol]
+    else:
+        return None
+
+
+def determine_symmetry(coords, symbols):
+    """
+    Determine external symmetry and chirality (optical isomers) of the species.
+
+    Args:
+        coords (list): The 3D coordinates of a molecule in array-format.
+        symbols (list): Entries are atomic symbols correcponding to coords.
+
+    Returns:
+        int: The external symmetry number.
+        int: 1 if no chiral centers are present, 2 if chiral centers are present.
+    """
+    atom_numbers = list()  # List of atomic numbers
+    for symbol in symbols:
+        atom_numbers.append(getElement(str(symbol)).number)
+    coords = np.array(coords, np.float64)  # N x 3 numpy.ndarray of atomic coordinates
+    #  in the same order as `atom_numbers`
+    unique_id = '0'  # Just some name that the SYMMETRY code gives to one of its jobs
+    scr_dir = os.path.join(arc_path, str('scratch'))  # Scratch directory that the SYMMETRY code writes its files in
+    if not os.path.exists(scr_dir):
+        os.makedirs(scr_dir)
+    symmetry = optical_isomers = 1
+    qmdata = QMData(
+        groundStateDegeneracy=1,  # Only needed to check if valid QMData
+        numberOfAtoms=len(atom_numbers),
+        atomicNumbers=atom_numbers,
+        atomCoords=(coords, str('angstrom')),
+        energy=(0.0, str('kcal/mol'))  # Only needed to avoid error
+    )
+    settings = type(str(''), (), dict(symmetryPath=str('symmetry'), scratchDirectory=scr_dir))()
+    pgc = PointGroupCalculator(settings, unique_id, qmdata)
+    pg = pgc.calculate()
+    if pg is not None:
+        symmetry = pg.symmetryNumber
+        optical_isomers = 2 if pg.chiral else optical_isomers
+    return symmetry, optical_isomers
+
+
+def min_list(lst):
+    """A helper function for finding the minimum of a list of integers where some of the entries might be None"""
+    if len(lst) == 0:
+        return None
+    elif len(lst) == 1:
+        return lst[0]
+    elif all([entry is None for entry in lst]):
+        return None
+    return min([entry for entry in lst if entry is not None])
diff --git a/arc/commonTest.py b/arc/commonTest.py
@@ -10,7 +10,7 @@
 import os
 import time
 
-from arc.common import read_file, get_git_commit, time_lapse, check_ess_settings
+import arc.common as common
 from arc.settings import arc_path, servers
 from arc.arc_exceptions import InputError, SettingsError
 
@@ -22,10 +22,10 @@ class TestARC(unittest.TestCase):
     Contains unit tests for ARC's common module
     """
 
-    def test_read_file(self):
-        """Test the read_file() function"""
+    def test_read_yaml_file(self):
+        """Test the read_yaml_file() function"""
         restart_path = os.path.join(arc_path, 'arc', 'testing', 'restart(H,H2O2,N2H3,CH3CO2).yml')
-        input_dict = read_file(restart_path)
+        input_dict = common.read_yaml_file(restart_path)
         self.assertIsInstance(input_dict, dict)
         self.assertTrue('reactions' in input_dict)
         self.assertTrue('freq_level' in input_dict)
@@ -34,11 +34,11 @@ def test_read_file(self):
         self.assertTrue('running_jobs' in input_dict)
 
         with self.assertRaises(InputError):
-            read_file('nopath')
+            common.read_yaml_file('nopath')
 
     def test_get_git_commit(self):
         """Test the get_git_commit() function"""
-        git_commit = get_git_commit()
+        git_commit = common.get_git_commit()
         # output format: ['fafdb957049917ede565cebc58b29899f597fb5a', 'Fri Mar 29 11:09:50 2019 -0400']
         self.assertEqual(len(git_commit[0]), 40)
         self.assertEqual(len(git_commit[1].split()), 6)
@@ -47,7 +47,7 @@ def test_time_lapse(self):
         """Test the time_lapse() function"""
         t0 = time.time()
         time.sleep(2)
-        lap = time_lapse(t0)
+        lap = common.time_lapse(t0)
         self.assertEqual(lap, '00:00:02')
 
     def test_check_ess_settings(self):
@@ -61,11 +61,11 @@ def test_check_ess_settings(self):
         ess_settings4 = {'gaussian': server_names[0], 'molpro': server_names[1], 'qchem': server_names[0]}
         ess_settings5 = {'gaussian': 'local', 'molpro': server_names[1], 'qchem': server_names[0]}
 
-        ess_settings1 = check_ess_settings(ess_settings1)
-        ess_settings2 = check_ess_settings(ess_settings2)
-        ess_settings3 = check_ess_settings(ess_settings3)
-        ess_settings4 = check_ess_settings(ess_settings4)
-        ess_settings5 = check_ess_settings(ess_settings5)
+        ess_settings1 = common.check_ess_settings(ess_settings1)
+        ess_settings2 = common.check_ess_settings(ess_settings2)
+        ess_settings3 = common.check_ess_settings(ess_settings3)
+        ess_settings4 = common.check_ess_settings(ess_settings4)
+        ess_settings5 = common.check_ess_settings(ess_settings5)
 
         ess_list = [ess_settings1, ess_settings2, ess_settings3, ess_settings4, ess_settings5]
 
@@ -76,10 +76,48 @@ def test_check_ess_settings(self):
 
         with self.assertRaises(SettingsError):
             ess_settings6 = {'nosoft': ['server1']}
-            check_ess_settings(ess_settings6)
+            common.check_ess_settings(ess_settings6)
         with self.assertRaises(SettingsError):
             ess_settings7 = {'gaussian': ['noserver']}
-            check_ess_settings(ess_settings7)
+            common.check_ess_settings(ess_settings7)
+
+    def test_min_list(self):
+        """Test the min_list() function"""
+        lst = []
+        min_lst = common.min_list(lst)
+        self.assertEqual(min_lst, None)
+
+        lst = [None]
+        min_lst = common.min_list(lst)
+        self.assertEqual(min_lst, None)
+
+        lst = [None, None]
+        min_lst = common.min_list(lst)
+        self.assertEqual(min_lst, None)
+
+        lst = [0]
+        min_lst = common.min_list(lst)
+        self.assertEqual(min_lst, 0)
+
+        lst = [-8]
+        min_lst = common.min_list(lst)
+        self.assertEqual(min_lst, -8)
+
+        lst = [-8, -80]
+        min_lst = common.min_list(lst)
+        self.assertEqual(min_lst, -80)
+
+        lst = [-8, None]
+        min_lst = common.min_list(lst)
+        self.assertEqual(min_lst, -8)
+
+        lst = [-8, -8, -8, -8]
+        min_lst = common.min_list(lst)
+        self.assertEqual(min_lst, -8)
+
+        lst = [-8, None, None, 100, -79, None]
+        min_lst = common.min_list(lst)
+        self.assertEqual(min_lst, -79)
 
 ################################################################################