diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 6a30a980d4..7a4179be0a 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -86,16 +86,6 @@ jobs: julia -e 'using Pkg; Pkg.add(PackageSpec(name="ReactionMechanismSimulator",rev="main")); using ReactionMechanismSimulator' || true ln -sfn $(which python-jl) $(which python) - # Attempt to install MOPAC - - name: Install MOPAC - env: - MOPACKEY: ${{ secrets.MOPACKEY }} - timeout-minutes: 1 - continue-on-error: true # allowed to fail on pull request from a forked repository - run: | - set +o pipefail - yes 'Yes' | ${CONDA_PREFIX}/bin/mopac "$MOPACKEY" - # non-regression testing - name: Unit tests run: make test-unittests diff --git a/documentation/source/users/rmg/installation/anacondaDeveloper.rst b/documentation/source/users/rmg/installation/anacondaDeveloper.rst index 5a28d0632e..63d1716e74 100644 --- a/documentation/source/users/rmg/installation/anacondaDeveloper.rst +++ b/documentation/source/users/rmg/installation/anacondaDeveloper.rst @@ -137,10 +137,6 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux python-jl replace/with/path/to/rmg.py input.py -#. Optional: If you wish to use the :ref:`QMTP interface ` with `MOPAC `_ to run quantum mechanical calculations for improved thermochemistry estimates of cyclic species, please obtain a legal license through the `MOPAC License Request Form `_. Once you have it, type the following into your terminal :: - - mopac password_string_here - You may now use RMG-Py, Arkane, as well as any of the :ref:`Standalone Modules ` included in the RMG-Py package. diff --git a/documentation/source/users/rmg/installation/anacondaUser.rst b/documentation/source/users/rmg/installation/anacondaUser.rst index 513394b621..0d2bc4ba03 100644 --- a/documentation/source/users/rmg/installation/anacondaUser.rst +++ b/documentation/source/users/rmg/installation/anacondaUser.rst @@ -22,10 +22,6 @@ Binary Installation Using Anaconda for Unix-Based Systems: Linux and Mac OSX source activate rmg_env -#. Optional: If you wish to use the :ref:`QMTP interface ` with `MOPAC `_ to run quantum mechanical calculations for improved thermochemistry estimates of cyclic species, please obtain a legal license through the `MOPAC License Request Form `_. Once you have it, type the following into your terminal :: - - mopac password_string_here - #. You may now run an RMG test job. Save the `Minimal Example Input File `_ to a local directory. Use the terminal to run your RMG job inside that folder using the following command :: diff --git a/documentation/source/users/rmg/installation/dependencies.rst b/documentation/source/users/rmg/installation/dependencies.rst index 70709a3fec..1e257148f9 100644 --- a/documentation/source/users/rmg/installation/dependencies.rst +++ b/documentation/source/users/rmg/installation/dependencies.rst @@ -63,5 +63,4 @@ All of RMG's dependencies except the ones listed below are freely available and If you wish to do on-the-fly quantum chemistry calculations of thermochemistry (advisable for fused cyclic species in particular, where the ring corrections to group additive estimates are lacking), the then you will need the third-party software for the QM calculations: -* **gaussian**: Gaussian03 and Gaussian09 are currently supported and commercially available. See `https://gaussian.com `_ for more details. -* **mopac** MOPAC can be found at `http://openmopac.net/ `_. Though it is not free for industrial use, it is free for non-profit and academic research use. +* **gaussian**: Gaussian03 and Gaussian09 are currently supported and commercially available. See `https://gaussian.com `_ for more details. diff --git a/environment.yml b/environment.yml index da98a778f1..f16342f1c5 100644 --- a/environment.yml +++ b/environment.yml @@ -23,7 +23,7 @@ dependencies: - rmg::lpsolve55 - markupsafe - matplotlib >=1.5 - - rmg::mopac + - conda-forge::mopac - mpmath - rmg::muq2 - networkx @@ -34,9 +34,9 @@ dependencies: - conda-forge::openbabel >= 3 - pandas - psutil - - rmg::pydas >=1.0.2 + - rmg::pydas >=1.0.3 - pydot - - rmg::pydqed >=1.0.1 + - rmg::pydqed >=1.0.3 - rmg::pyjulia - pymongo - pyparsing diff --git a/rmgpy/qm/gaussian.py b/rmgpy/qm/gaussian.py index a1d33b7845..993f341ded 100644 --- a/rmgpy/qm/gaussian.py +++ b/rmgpy/qm/gaussian.py @@ -44,57 +44,67 @@ class Gaussian: """ A base class for all QM calculations that use Gaussian. - + Classes such as :class:`GaussianMol` will inherit from this class. """ - inputFileExtension = '.gjf' - outputFileExtension = '.log' + input_file_extension = ".gjf" + output_file_extension = ".log" - executablesToTry = ('g16', 'g09', 'g03') + executablesToTry = ("g16", "g09", "g03") for exe in executablesToTry: try: - executablePath = distutils.spawn.find_executable(exe) + executable_path = distutils.spawn.find_executable(exe) except: - executablePath = None - if executablePath is not None: + executable_path = None + if executable_path is not None: break else: # didn't break - logging.debug("Did not find Gaussian on path, checking if it exists in a declared GAUSS_EXEDIR, g09root or g03root...") - gaussEnv = os.getenv('GAUSS_EXEDIR') or os.getenv('g09root') or os.getenv('g03root') or "" - possibleDirs = gaussEnv.split(':') # GAUSS_EXEDIR may be a list like "path1:path2:path3" + logging.debug( + "Did not find Gaussian on path, checking if it exists in a declared GAUSS_EXEDIR, g09root or g03root..." + ) + gaussEnv = ( + os.getenv("GAUSS_EXEDIR") + or os.getenv("g09root") + or os.getenv("g03root") + or "" + ) + possibleDirs = gaussEnv.split( + ":" + ) # GAUSS_EXEDIR may be a list like "path1:path2:path3" for exe, possibleDir in itertools.product(executablesToTry, possibleDirs): - executablePath = os.path.join(possibleDir, exe) - if os.path.exists(executablePath): + executable_path = os.path.join(possibleDir, exe) + if os.path.exists(executable_path): break else: # didn't break - executablePath = os.path.join(gaussEnv, '(Gaussian 2003 or 2009)') + executable_path = os.path.join(gaussEnv, "(Gaussian 2003 or 2009)") - usePolar = False + use_polar = False #: List of phrases that indicate failure #: NONE of these must be present in a succesful job. - failureKeys = [ - 'ERROR TERMINATION', - 'IMAGINARY FREQUENCIES' - ] + failure_keys = ["ERROR TERMINATION", "IMAGINARY FREQUENCIES"] #: List of phrases to indicate success. #: ALL of these must be present in a successful job. - successKeys = [ - 'Normal termination of Gaussian' - ] + success_keys = ["Normal termination of Gaussian"] def test_ready(self): - if not os.path.exists(self.executablePath): - raise Exception("Couldn't find Gaussian executable at {0}. " - "Try setting your GAUSS_EXEDIR environment variable.".format(self.executablePath)) + if not os.path.exists(self.executable_path): + raise Exception( + "Couldn't find Gaussian executable at {0}. " + "Try setting your GAUSS_EXEDIR environment variable.".format( + self.executable_path + ) + ) def run(self): self.test_ready() # submits the input file to Gaussian - process = Popen([self.executablePath, self.input_file_path, self.output_file_path]) + process = Popen( + [self.executable_path, self.input_file_path, self.output_file_path] + ) process.communicate() # necessary to wait for executable termination! return self.verify_output_file() @@ -102,10 +112,10 @@ def run(self): def verify_output_file(self): """ Check's that an output file exists and was successful. - - Returns a boolean flag that states whether a successful GAUSSIAN simulation already exists for the molecule with the + + Returns a boolean flag that states whether a successful GAUSSIAN simulation already exists for the molecule with the given (augmented) InChI Key. - + The definition of finding a successful simulation is based on these criteria: 1) finding an output file with the file name equal to the InChI Key 2) NOT finding any of the keywords that are denote a calculation failure @@ -117,25 +127,33 @@ def verify_output_file(self): If all are satisfied, it will return True. """ if not os.path.exists(self.output_file_path): - logging.info("Output file {0} does not exist.".format(self.output_file_path)) + logging.info( + "Output file {0} does not exist.".format(self.output_file_path) + ) return False inchi_match = False # flag (1 or 0) indicating whether the InChI in the file matches InChIaug this can only be 1 if inchi_found is also 1 - inchi_found = False # flag (1 or 0) indicating whether an InChI was found in the log file + inchi_found = ( + False # flag (1 or 0) indicating whether an InChI was found in the log file + ) - # Initialize dictionary with "False"s - success_keys_found = dict([(key, False) for key in self.successKeys]) + # Initialize dictionary with "False"s + success_keys_found = dict([(key, False) for key in self.success_keys]) with open(self.output_file_path) as outputFile: for line in outputFile: line = line.strip() - for element in self.failureKeys: # search for failure keywords + for element in self.failure_keys: # search for failure keywords if element in line: - logging.error("Gaussian output file contains the following error: {0}".format(element)) + logging.error( + "Gaussian output file contains the following error: {0}".format( + element + ) + ) return False - for element in self.successKeys: # search for success keywords + for element in self.success_keys: # search for success keywords if element in line: success_keys_found[element] = True @@ -145,26 +163,42 @@ def verify_output_file(self): if self.unique_id_long in log_file_inchi: inchi_match = True elif self.unique_id_long.startswith(log_file_inchi): - logging.info("InChI too long to check, but beginning matches so assuming OK.") + logging.info( + "InChI too long to check, but beginning matches so assuming OK." + ) inchi_match = True else: - logging.warning("InChI in log file ({0}) didn't match that in geometry " - "({1}).".format(log_file_inchi, self.geometry.unique_id_long)) + logging.warning( + "InChI in log file ({0}) didn't match that in geometry " + "({1}).".format( + log_file_inchi, self.geometry.unique_id_long + ) + ) if self.geometry.unique_id_long.startswith(log_file_inchi): - logging.warning("but the beginning matches so it's probably just a truncation problem.") + logging.warning( + "but the beginning matches so it's probably just a truncation problem." + ) inchi_match = True # Check that ALL 'success' keywords were found in the file. if not all(success_keys_found.values()): - logging.error('Not all of the required keywords for success were found in the output file!') + logging.error( + "Not all of the required keywords for success were found in the output file!" + ) return False if not inchi_found: - logging.error("No InChI was found in the Gaussian output file {0}".format(self.output_file_path)) + logging.error( + "No InChI was found in the Gaussian output file {0}".format( + self.output_file_path + ) + ) return False if not inchi_match: # InChIs do not match (most likely due to limited name length mirrored in log file (240 characters), but possibly due to a collision) - return self.checkForInChiKeyCollision(log_file_inchi) # Not yet implemented! + return self.checkForInChiKeyCollision( + log_file_inchi + ) # Not yet implemented! # Compare the optimized geometry to the original molecule qm_data = self.parse() @@ -172,10 +206,18 @@ def verify_output_file(self): cclib_mol.from_xyz(qm_data.atomicNumbers, qm_data.atomCoords.value) test_mol = self.molecule.to_single_bonds() if not cclib_mol.is_isomorphic(test_mol): - logging.info("Incorrect connectivity for optimized geometry in file {0}".format(self.output_file_path)) + logging.info( + "Incorrect connectivity for optimized geometry in file {0}".format( + self.output_file_path + ) + ) return False - logging.info("Successful {1} quantum result in {0}".format(self.output_file_path, self.__class__.__name__)) + logging.info( + "Successful {1} quantum result in {0}".format( + self.output_file_path, self.__class__.__name__ + ) + ) return True def parse(self): @@ -183,7 +225,9 @@ def parse(self): Parses the results of the Gaussian calculation, and returns a QMData object. """ parser = cclib.parser.Gaussian(self.output_file_path) - parser.logger.setLevel(logging.ERROR) # cf. http://cclib.sourceforge.net/wiki/index.php/Using_cclib#Additional_information + parser.logger.setLevel( + logging.ERROR + ) # cf. http://cclib.sourceforge.net/wiki/index.php/Using_cclib#Additional_information cclib_data = parser.parse() radical_number = sum([i.radical_electrons for i in self.molecule.atoms]) qm_data = parse_cclib_data(cclib_data, radical_number + 1) @@ -192,15 +236,15 @@ def parse(self): class GaussianMol(QMMolecule, Gaussian): """ - A base Class for calculations of molecules using Gaussian. - + A base Class for calculations of molecules using Gaussian. + Inherits from both :class:`QMMolecule` and :class:`Gaussian`. """ def input_file_keywords(self, attempt): """ Return the top keywords for attempt number `attempt`. - + NB. `attempt` begins at 1, not 0. """ assert attempt <= self.max_attempts @@ -214,10 +258,16 @@ def write_input_file(self, attempt): for the `attempt`. """ molfile = self.get_mol_file_path_for_calculation(attempt) - atomline = re.compile(r'\s*([\- ][0-9.]+\s+[\-0-9.]+\s+[\-0-9.]+)\s+([A-Za-z]+)') + atomline = re.compile( + r"\s*([\- ][0-9.]+\s+[\-0-9.]+\s+[\-0-9.]+)\s+([A-Za-z]+)" + ) - output = ['', self.geometry.unique_id_long, ''] - output.append("{charge} {mult}".format(charge=0, mult=(self.molecule.get_radical_count() + 1))) + output = ["", self.geometry.unique_id_long, ""] + output.append( + "{charge} {mult}".format( + charge=0, mult=(self.molecule.get_radical_count() + 1) + ) + ) atom_count = 0 with open(molfile) as molinput: @@ -228,17 +278,17 @@ def write_input_file(self, attempt): atom_count += 1 assert atom_count == len(self.molecule.atoms) - output.append('') - input_string = '\n'.join(output) + output.append("") + input_string = "\n".join(output) top_keys = self.input_file_keywords(attempt) - with open(self.input_file_path, 'w') as gaussian_file: + with open(self.input_file_path, "w") as gaussian_file: gaussian_file.write(top_keys) - gaussian_file.write('\n') + gaussian_file.write("\n") gaussian_file.write(input_string) - gaussian_file.write('\n') - if self.usePolar: - gaussian_file.write('\n\n\n') + gaussian_file.write("\n") + if self.use_polar: + gaussian_file.write("\n\n\n") raise NotImplementedError("Not sure what should be here, if anything.") # gaussian_file.write(polar_keys) @@ -253,23 +303,41 @@ def generate_qm_data(self): if self.verify_output_file(): logging.info("Found a successful output file already; using that.") - source = "QM {0} calculation found from previous run.".format(self.__class__.__name__) + source = "QM {0} calculation found from previous run.".format( + self.__class__.__name__ + ) else: self.create_geometry() success = False for attempt in range(1, self.max_attempts + 1): self.write_input_file(attempt) - logging.info('Trying {3} attempt {0} of {1} on molecule {2}.'.format(attempt, self.max_attempts, - self.molecule.to_smiles(), - self.__class__.__name__)) + logging.info( + "Trying {3} attempt {0} of {1} on molecule {2}.".format( + attempt, + self.max_attempts, + self.molecule.to_smiles(), + self.__class__.__name__, + ) + ) success = self.run() if success: - logging.info('Attempt {0} of {1} on species {2} succeeded.'.format(attempt, self.max_attempts, - self.molecule.to_augmented_inchi())) - source = "QM {0} calculation attempt {1}".format(self.__class__.__name__, attempt) + logging.info( + "Attempt {0} of {1} on species {2} succeeded.".format( + attempt, + self.max_attempts, + self.molecule.to_augmented_inchi(), + ) + ) + source = "QM {0} calculation attempt {1}".format( + self.__class__.__name__, attempt + ) break else: - logging.error('QM thermo calculation failed for {0}.'.format(self.molecule.to_augmented_inchi())) + logging.error( + "QM thermo calculation failed for {0}.".format( + self.molecule.to_augmented_inchi() + ) + ) return None result = self.parse() # parsed in cclib result.source = source @@ -285,10 +353,11 @@ def get_parser(self, output_file): class GaussianMolPM3(GaussianMol): """ Gaussian PM3 calculations for molecules - + This is a class of its own in case you wish to do anything differently, but for now it's only the 'pm3' in the keywords that differs. """ + #: Keywords that will be added at the top of the qm input file keywords = [ # The combinations of keywords were derived by Greg Magoon for pm3 in Gaussian. His comments are attached to each combination. @@ -312,13 +381,15 @@ class GaussianMolPM3(GaussianMol): "# pm3 opt=(calcall,small,maxcyc=100) IOP(2/16=3)", # used to address troublesome FILUFGAZMJGNEN-UHFFFAOYAImult3 case (InChI=1/C5H6/c1-3-5-4-2/h3H,1H2,2H3/mult3) ] + class GaussianMolPM6(GaussianMol): """ Gaussian PM6 calculations for molecules - + This is a class of its own in case you wish to do anything differently, but for now it's only the 'pm6' in the keywords that differs. """ + #: Keywords that will be added at the top of the qm input file keywords = [ # The combinations of keywords were derived by Greg Magoon for pm3. For now, we assume similar ones will work for pm6: diff --git a/rmgpy/qm/gaussianTest.py b/rmgpy/qm/gaussianTest.py index 8ba8723b1e..713d9c5365 100644 --- a/rmgpy/qm/gaussianTest.py +++ b/rmgpy/qm/gaussianTest.py @@ -38,10 +38,10 @@ from rmgpy.qm.gaussian import Gaussian, GaussianMolPM3, GaussianMolPM6 from rmgpy.qm.main import QMCalculator -executablePath = Gaussian.executablePath -NO_GAUSSIAN = not os.path.exists(executablePath) +executable_path = Gaussian.executable_path +NO_GAUSSIAN = not os.path.exists(executable_path) -mol1 = Molecule().from_smiles('C1=CC=C2C=CC=CC2=C1') +mol1 = Molecule().from_smiles("C1=CC=C2C=CC=CC2=C1") class TestGaussianMolPM3(unittest.TestCase): @@ -49,18 +49,22 @@ class TestGaussianMolPM3(unittest.TestCase): Contains unit tests for the Geometry class. """ - @unittest.skipIf(NO_GAUSSIAN, "Gaussian not found. Try resetting your environment variables if you want to use it.") + @unittest.skipIf( + NO_GAUSSIAN, + "Gaussian not found. Try resetting your environment variables if you want to use it.", + ) def setUp(self): """ A function run before each unit test in this class. """ - rmgpy_path = os.path.normpath(os.path.join(get_path(), '..')) + rmgpy_path = os.path.normpath(os.path.join(get_path(), "..")) - qm = QMCalculator(software='gaussian', - method='pm3', - fileStore=os.path.join(rmgpy_path, 'testing', 'qm', 'QMfiles'), - scratchDirectory=os.path.join(rmgpy_path, 'testing', 'qm', 'QMscratch'), - ) + qm = QMCalculator( + software="gaussian", + method="pm3", + fileStore=os.path.join(rmgpy_path, "testing", "qm", "QMfiles"), + scratchDirectory=os.path.join(rmgpy_path, "testing", "qm", "QMscratch"), + ) if not os.path.exists(qm.settings.fileStore): os.makedirs(qm.settings.fileStore) @@ -72,16 +76,21 @@ def test_generate_thermo_data(self): Test that generate_thermo_data() works correctly on gaussian PM3. """ # First ensure any old data are removed, or else they'll be reused! - for directory in (self.qmmol1.settings.fileStore, self.qmmol1.settings.scratchDirectory): + for directory in ( + self.qmmol1.settings.fileStore, + self.qmmol1.settings.scratchDirectory, + ): shutil.rmtree(directory, ignore_errors=True) self.qmmol1.generate_thermo_data() result = self.qmmol1.qm_data - self.assertTrue(self.qmmol1.thermo.comment.startswith('QM GaussianMolPM3 calculation')) + self.assertTrue( + self.qmmol1.thermo.comment.startswith("QM GaussianMolPM3 calculation") + ) self.assertEqual(result.numberOfAtoms, 18) self.assertIsInstance(result.atomicNumbers, np.ndarray) - if result.molecularMass.units == 'amu': + if result.molecularMass.units == "amu": self.assertAlmostEqual(result.molecularMass.value, 128.0626, 3) def test_load_thermo_data(self): @@ -94,10 +103,12 @@ def test_load_thermo_data(self): self.qmmol1.generate_thermo_data() result = self.qmmol1.qm_data - self.assertTrue(self.qmmol1.thermo.comment.startswith('QM GaussianMolPM3 calculation')) + self.assertTrue( + self.qmmol1.thermo.comment.startswith("QM GaussianMolPM3 calculation") + ) self.assertEqual(result.numberOfAtoms, 18) self.assertIsInstance(result.atomicNumbers, np.ndarray) - if result.molecularMass.units == 'amu': + if result.molecularMass.units == "amu": self.assertAlmostEqual(result.molecularMass.value, 128.0626, 3) @@ -106,43 +117,56 @@ class TestGaussianMolPM6(unittest.TestCase): Contains unit tests for the Geometry class. """ - @unittest.skipIf(NO_GAUSSIAN, "Gaussian not found. Try resetting your environment variables if you want to use it.") + @unittest.skipIf( + NO_GAUSSIAN, + "Gaussian not found. Try resetting your environment variables if you want to use it.", + ) def setUp(self): """ A function run before each unit test in this class. """ - rmgpy_path = os.path.normpath(os.path.join(get_path(), '..')) + rmgpy_path = os.path.normpath(os.path.join(get_path(), "..")) - qm = QMCalculator(software='gaussian', - method='pm6', - fileStore=os.path.join(rmgpy_path, 'testing', 'qm', 'QMfiles'), - scratchDirectory=os.path.join(rmgpy_path, 'testing', 'qm', 'QMscratch'), - ) + qm = QMCalculator( + software="gaussian", + method="pm6", + fileStore=os.path.join(rmgpy_path, "testing", "qm", "QMfiles"), + scratchDirectory=os.path.join(rmgpy_path, "testing", "qm", "QMscratch"), + ) if not os.path.exists(qm.settings.fileStore): os.makedirs(qm.settings.fileStore) self.qmmol1 = GaussianMolPM6(mol1, qm.settings) - @unittest.skipIf('g03' in executablePath, "This test was shown not to work on g03.") + @unittest.skipIf( + "g03" in executable_path, "This test was shown not to work on g03." + ) def test_generate_thermo_data(self): """ Test that generate_thermo_data() works correctly for gaussian PM6. """ # First ensure any old data are removed, or else they'll be reused! - for directory in (self.qmmol1.settings.fileStore, self.qmmol1.settings.scratchDirectory): + for directory in ( + self.qmmol1.settings.fileStore, + self.qmmol1.settings.scratchDirectory, + ): shutil.rmtree(directory, ignore_errors=True) self.qmmol1.generate_thermo_data() result = self.qmmol1.qm_data - self.assertTrue(self.qmmol1.thermo.comment.startswith('QM GaussianMolPM6 calculation')) + self.assertTrue( + self.qmmol1.thermo.comment.startswith("QM GaussianMolPM6 calculation") + ) self.assertEqual(result.numberOfAtoms, 18) self.assertIsInstance(result.atomicNumbers, np.ndarray) - if result.molecularMass.units == 'amu': + if result.molecularMass.units == "amu": self.assertAlmostEqual(result.molecularMass.value, 128.0626, 3) - @unittest.skipIf('g03' in executablePath, "This test was shown not to work on g03.") + @unittest.skipIf( + "g03" in executable_path, "This test was shown not to work on g03." + ) def test_load_thermo_data(self): """ Test that generate_thermo_data() can load thermo from the previous gaussian PM6 run. @@ -153,14 +177,16 @@ def test_load_thermo_data(self): self.qmmol1.generate_thermo_data() result = self.qmmol1.qm_data - self.assertTrue(self.qmmol1.thermo.comment.startswith('QM GaussianMolPM6 calculation')) + self.assertTrue( + self.qmmol1.thermo.comment.startswith("QM GaussianMolPM6 calculation") + ) self.assertEqual(result.numberOfAtoms, 18) self.assertIsInstance(result.atomicNumbers, np.ndarray) - if result.molecularMass.units == 'amu': + if result.molecularMass.units == "amu": self.assertAlmostEqual(result.molecularMass.value, 128.0626, 3) ################################################################################ -if __name__ == '__main__': +if __name__ == "__main__": unittest.main(testRunner=unittest.TextTestRunner(verbosity=2)) diff --git a/rmgpy/qm/mainTest.py b/rmgpy/qm/mainTest.py index e1b408f28d..dcfe492910 100644 --- a/rmgpy/qm/mainTest.py +++ b/rmgpy/qm/mainTest.py @@ -40,6 +40,9 @@ from rmgpy.species import Species +MOPAC_CLOSE_ENOUGH_PERCENT = 0.005 # 0.5% + + class TestQMSettings(unittest.TestCase): """ Contains unit tests for the QMSettings class. @@ -49,15 +52,16 @@ def setUp(self): """ A function run before each unit test in this class. """ - rmg_path = os.path.normpath(os.path.join(get_path(), '..')) + rmg_path = os.path.normpath(os.path.join(get_path(), "..")) - self.settings1 = QMSettings(software='mopac', - method='pm3', - fileStore=os.path.join(rmg_path, 'testing', 'qm', 'QMfiles'), - scratchDirectory=None, - onlyCyclics=False, - maxRadicalNumber=0, - ) + self.settings1 = QMSettings( + software="mopac", + method="pm3", + fileStore=os.path.join(rmg_path, "testing", "qm", "QMfiles"), + scratchDirectory=None, + onlyCyclics=False, + maxRadicalNumber=0, + ) self.settings2 = QMSettings() @@ -79,65 +83,51 @@ class TestQMCalculator(unittest.TestCase): Contains unit tests for the QMSettings class. """ - mopExecutablePath = Mopac.executablePath - if not os.path.exists(mopExecutablePath): - NO_MOPAC = NO_LICENCE = True - else: - NO_MOPAC = False - process = subprocess.Popen(mopExecutablePath, - stdin=subprocess.PIPE, - stdout=subprocess.PIPE, - stderr=subprocess.PIPE) - stdut, stderr = process.communicate(b'\n') - NO_LICENCE = b'To install the MOPAC license' in stderr - - gaussExecutablePath = Gaussian.executablePath - NO_GAUSSIAN = not os.path.exists(gaussExecutablePath) + mopexecutable_path = Mopac.executable_path + + gaussexecutable_path = Gaussian.executable_path + NO_GAUSSIAN = not os.path.exists(gaussexecutable_path) def setUp(self): """ A function run before each unit test in this class. """ - rmg_path = os.path.normpath(os.path.join(get_path(), '..')) - self.fileStore = os.path.join(rmg_path, 'testing', 'qm', 'QMfiles') - - self.mop1 = QMCalculator(software='mopac', - method='pm3', - fileStore=self.fileStore - ) - - self.mop2 = QMCalculator(software='mopac', - method='pm6', - ) - - self.mop3 = QMCalculator(software='mopac', - method='pm7', - fileStore=self.fileStore - ) - - self.mop4 = QMCalculator(software='mopac', - method='pm8', - fileStore=self.fileStore - ) - - self.gauss1 = QMCalculator(software='gaussian', - method='pm3', - ) - - self.gauss2 = QMCalculator(software='gaussian', - method='pm6', - fileStore=self.fileStore - ) - - self.gauss3 = QMCalculator(software='gaussian', - method='pm7', - fileStore=self.fileStore - ) - - self.molpro1 = QMCalculator(software='molpro', - method='mp2', - fileStore=self.fileStore - ) + rmg_path = os.path.normpath(os.path.join(get_path(), "..")) + self.fileStore = os.path.join(rmg_path, "testing", "qm", "QMfiles") + + self.mop1 = QMCalculator( + software="mopac", method="pm3", fileStore=self.fileStore + ) + + self.mop2 = QMCalculator( + software="mopac", + method="pm6", + ) + + self.mop3 = QMCalculator( + software="mopac", method="pm7", fileStore=self.fileStore + ) + + self.mop4 = QMCalculator( + software="mopac", method="pm8", fileStore=self.fileStore + ) + + self.gauss1 = QMCalculator( + software="gaussian", + method="pm3", + ) + + self.gauss2 = QMCalculator( + software="gaussian", method="pm6", fileStore=self.fileStore + ) + + self.gauss3 = QMCalculator( + software="gaussian", method="pm7", fileStore=self.fileStore + ) + + self.molpro1 = QMCalculator( + software="molpro", method="mp2", fileStore=self.fileStore + ) self.qmmol1 = QMCalculator(fileStore=self.fileStore) @@ -161,7 +151,7 @@ def test_set_default_output_directory(self): self.assertIsNone(self.gauss2.settings.scratchDirectory) # Now set the default directories for those not set - outputDirectory = os.path.join(self.mop1.settings.fileStore, '..', '..') + outputDirectory = os.path.join(self.mop1.settings.fileStore, "..", "..") self.mop1.set_default_output_directory(outputDirectory) self.mop2.set_default_output_directory(outputDirectory) self.mop3.set_default_output_directory(outputDirectory) @@ -185,7 +175,7 @@ def test_initialize(self): """ # Now set the default directories for those not set - outputDirectory = os.path.join(self.mop1.settings.fileStore, '..', '..') + outputDirectory = os.path.join(self.mop1.settings.fileStore, "..", "..") self.mop1.set_default_output_directory(outputDirectory) self.mop2.set_default_output_directory(outputDirectory) self.mop3.set_default_output_directory(outputDirectory) @@ -201,104 +191,141 @@ def test_initialize(self): except AssertionError: self.fail("initialize() raised unexpected AssertionError.") except Exception: - self.fail("initialize() raised Exception. Output file paths not correctly set.") + self.fail( + "initialize() raised Exception. Output file paths not correctly set." + ) def test_get_thermo_data(self): """ Test that get_thermo_data() fails when expected. """ - output_directory = os.path.join(self.mop4.settings.fileStore, '..', '..') + output_directory = os.path.join(self.mop4.settings.fileStore, "..", "..") self.mop4.set_default_output_directory(output_directory) self.gauss3.set_default_output_directory(output_directory) self.molpro1.set_default_output_directory(output_directory) - mol = Molecule().from_smiles('C1=CC=C2C=CC=CC2=C1') + mol = Molecule().from_smiles("C1=CC=C2C=CC=CC2=C1") with self.assertRaises(Exception): self.mop4.get_thermo_data(mol) self.gauss3.get_thermo_data(mol) self.molpro1.get_thermo_data(mol) - @unittest.skipIf(NO_MOPAC, "MOPAC not found. Try resetting your environment variables if you want to use it.") - @unittest.skipIf(NO_LICENCE, "MOPAC license not installed. Run mopac for instructions") def test_get_thermo_data_mopac(self): """ Test that Mocpac get_thermo_data() works correctly. """ - output_directory = os.path.join(self.mop1.settings.fileStore, '..', '..') + output_directory = os.path.join(self.mop1.settings.fileStore, "..", "..") self.mop1.set_default_output_directory(output_directory) self.mop2.set_default_output_directory(output_directory) self.mop3.set_default_output_directory(output_directory) - mol = Molecule().from_smiles('C1=CC=C2C=CC=CC2=C1') + mol = Molecule().from_smiles("C1=CC=C2C=CC=CC2=C1") - for directory in (self.mop1.settings.fileStore, self.mop1.settings.scratchDirectory): + for directory in ( + self.mop1.settings.fileStore, + self.mop1.settings.scratchDirectory, + ): shutil.rmtree(directory, ignore_errors=True) - for directory in (self.mop2.settings.fileStore, self.mop2.settings.scratchDirectory): + for directory in ( + self.mop2.settings.fileStore, + self.mop2.settings.scratchDirectory, + ): shutil.rmtree(directory, ignore_errors=True) - for directory in (self.mop3.settings.fileStore, self.mop3.settings.scratchDirectory): + for directory in ( + self.mop3.settings.fileStore, + self.mop3.settings.scratchDirectory, + ): shutil.rmtree(directory, ignore_errors=True) thermo1 = self.mop1.get_thermo_data(mol) thermo2 = self.mop2.get_thermo_data(mol) thermo3 = self.mop3.get_thermo_data(mol) - self.assertTrue(thermo1.comment.startswith('QM MopacMolPM3')) - self.assertTrue(thermo2.comment.startswith('QM MopacMolPM6')) - self.assertTrue(thermo3.comment.startswith('QM MopacMolPM7')) - - self.assertAlmostEqual(thermo1.H298.value_si, 169708.0608, 1) # to 1 decimal place - self.assertAlmostEqual(thermo1.S298.value_si, 334.5007584, 1) # to 1 decimal place - self.assertAlmostEqual(thermo2.H298.value_si, 167704.4270, 1) # to 1 decimal place - self.assertAlmostEqual(thermo2.S298.value_si, 338.0999241, 1) # to 1 decimal place - self.assertAlmostEqual(thermo3.H298.value_si, 166168.8571, 1) # to 1 decimal place - self.assertAlmostEqual(thermo3.S298.value_si, 336.3330406, 1) # to 1 decimal place - - @unittest.skipIf(NO_GAUSSIAN, "Gaussian not found. Try resetting your environment variables if you want to use it.") + self.assertTrue(thermo1.comment.startswith("QM MopacMolPM3")) + self.assertTrue(thermo2.comment.startswith("QM MopacMolPM6")) + self.assertTrue(thermo3.comment.startswith("QM MopacMolPM7")) + + self.assertAlmostEqual( + thermo1.H298.value_si, 169708, delta=169708 * MOPAC_CLOSE_ENOUGH_PERCENT + ) + self.assertAlmostEqual( + thermo1.S298.value_si, 334.500, delta=334.500 * MOPAC_CLOSE_ENOUGH_PERCENT + ) + self.assertAlmostEqual( + thermo2.H298.value_si, 167704, delta=167704 * MOPAC_CLOSE_ENOUGH_PERCENT + ) + self.assertAlmostEqual( + thermo2.S298.value_si, 338.099, delta=338.099 * MOPAC_CLOSE_ENOUGH_PERCENT + ) + self.assertAlmostEqual( + thermo3.H298.value_si, 166168, delta=166168 * MOPAC_CLOSE_ENOUGH_PERCENT + ) + self.assertAlmostEqual( + thermo3.S298.value_si, 336.333, delta=336.333 * MOPAC_CLOSE_ENOUGH_PERCENT + ) + + @unittest.skipIf( + NO_GAUSSIAN, + "Gaussian not found. Try resetting your environment variables if you want to use it.", + ) def test_get_thermo_data_gaussian(self): """ Test that Gaussian get_thermo_data() works correctly. """ - output_directory = os.path.join(self.mop1.settings.fileStore, '..', '..') + output_directory = os.path.join(self.mop1.settings.fileStore, "..", "..") self.gauss1.set_default_output_directory(output_directory) self.gauss2.set_default_output_directory(output_directory) - mol = Molecule().from_smiles('C1=CC=C2C=CC=CC2=C1') + mol = Molecule().from_smiles("C1=CC=C2C=CC=CC2=C1") - for directory in (self.gauss1.settings.fileStore, self.gauss1.settings.scratchDirectory): + for directory in ( + self.gauss1.settings.fileStore, + self.gauss1.settings.scratchDirectory, + ): shutil.rmtree(directory, ignore_errors=True) - for directory in (self.gauss1.settings.fileStore, self.gauss2.settings.scratchDirectory): + for directory in ( + self.gauss1.settings.fileStore, + self.gauss2.settings.scratchDirectory, + ): shutil.rmtree(directory, ignore_errors=True) thermo1 = self.gauss1.get_thermo_data(mol) thermo2 = self.gauss2.get_thermo_data(mol) - self.assertTrue(thermo1.comment.startswith('QM GaussianMolPM3')) - self.assertTrue(thermo2.comment.startswith('QM GaussianMolPM6')) - - self.assertAlmostEqual(thermo1.H298.value_si, 169908.3376, 0) # to 1 decimal place - self.assertAlmostEqual(thermo1.S298.value_si, 335.5438748, 0) # to 1 decimal place - self.assertAlmostEqual(thermo2.H298.value_si, 169326.2504, 0) # to 1 decimal place - self.assertAlmostEqual(thermo2.S298.value_si, 338.2696063, 0) # to 1 decimal place + self.assertTrue(thermo1.comment.startswith("QM GaussianMolPM3")) + self.assertTrue(thermo2.comment.startswith("QM GaussianMolPM6")) + + self.assertAlmostEqual( + thermo1.H298.value_si, 169908.3376, 0 + ) # to 1 decimal place + self.assertAlmostEqual( + thermo1.S298.value_si, 335.5438748, 0 + ) # to 1 decimal place + self.assertAlmostEqual( + thermo2.H298.value_si, 169326.2504, 0 + ) # to 1 decimal place + self.assertAlmostEqual( + thermo2.S298.value_si, 338.2696063, 0 + ) # to 1 decimal place - @unittest.skipIf(NO_MOPAC, "MOPAC not found. Try resetting your environment variables if you want to use it.") - @unittest.skipIf(NO_LICENCE, "MOPAC license not installed. Run mopac for instructions") def test_run_jobs(self): """Test that run_jobs() works properly.""" - qm = QMCalculator(software='mopac', - method='pm3', - fileStore=self.fileStore, - onlyCyclics=True, - maxRadicalNumber=0, - ) - output_directory = os.path.join(qm.settings.fileStore, '..', '..') + qm = QMCalculator( + software="mopac", + method="pm3", + fileStore=self.fileStore, + onlyCyclics=True, + maxRadicalNumber=0, + ) + output_directory = os.path.join(qm.settings.fileStore, "..", "..") qm.set_default_output_directory(output_directory) - spc1 = Species().from_smiles('c1ccccc1') - spc2 = Species().from_smiles('CC1C=CC=CC=1') + spc1 = Species().from_smiles("c1ccccc1") + spc2 = Species().from_smiles("CC1C=CC=CC=1") spc_list = [spc1, spc2] qm.run_jobs(spc_list, procnum=1) @@ -306,5 +333,5 @@ def test_run_jobs(self): ################################################################################ -if __name__ == '__main__': +if __name__ == "__main__": unittest.main(testRunner=unittest.TextTestRunner(verbosity=2)) diff --git a/rmgpy/qm/molecule.py b/rmgpy/qm/molecule.py index 58abd463bb..50277de9ef 100644 --- a/rmgpy/qm/molecule.py +++ b/rmgpy/qm/molecule.py @@ -32,6 +32,7 @@ import os import numpy as np + try: from rdkit import Chem from rdkit.Chem import AllChem @@ -89,11 +90,17 @@ def __init__(self, settings, unique_id, molecule, unique_id_long=None): self.scratchDirectory = None if self.fileStore and not os.path.exists(self.fileStore): - logging.info("Creating permanent directory %s for qm files." % os.path.abspath(self.fileStore)) + logging.info( + "Creating permanent directory %s for qm files." + % os.path.abspath(self.fileStore) + ) os.makedirs(self.fileStore) if self.scratchDirectory and not os.path.exists(self.scratchDirectory): - logging.info("Creating scratch directory %s for qm files." % os.path.abspath(self.scratchDirectory)) + logging.info( + "Creating scratch directory %s for qm files." + % os.path.abspath(self.scratchDirectory) + ) os.makedirs(self.scratchDirectory) def get_file_path(self, extension, scratch=True): @@ -106,16 +113,16 @@ def get_file_path(self, extension, scratch=True): """ return os.path.join( self.settings.scratchDirectory if scratch else self.settings.fileStore, - self.unique_id + extension + self.unique_id + extension, ) def get_crude_mol_file_path(self): """Returns the path of the crude mol file.""" - return self.get_file_path('.crude.mol') + return self.get_file_path(".crude.mol") def get_refined_mol_file_path(self): """Returns the path the the refined mol file.""" - return self.get_file_path('.refined.mol') + return self.get_file_path(".refined.mol") def generate_rdkit_geometries(self): """ @@ -129,7 +136,9 @@ def generate_rdkit_geometries(self): atoms = len(self.molecule.atoms) dist_geom_attempts = 1 if atoms > 3: # this check prevents the number of attempts from being negative - dist_geom_attempts = 5 * (atoms - 3) # number of conformer attempts is just a linear scaling with molecule size, due to time considerations in practice, it is probably more like 3^(n-3) or something like that + dist_geom_attempts = 5 * ( + atoms - 3 + ) # number of conformer attempts is just a linear scaling with molecule size, due to time considerations in practice, it is probably more like 3^(n-3) or something like that rdmol, min_e_id = self.rd_embed(rdmol, dist_geom_attempts) self.save_coordinates_from_rdmol(rdmol, min_e_id, rdatom_idx) @@ -155,7 +164,9 @@ def rd_embed(self, rdmol, num_conf_attempts): # Try to embed using random coordinates. However, there # are still cases that it can fails: # https://github.com/rdkit/rdkit/issues/2996 - AllChem.EmbedMultipleConfs(rdmol, num_conf_attempts, useRandomCoords=True, randomSeed=1) + AllChem.EmbedMultipleConfs( + rdmol, num_conf_attempts, useRandomCoords=True, randomSeed=1 + ) if rdmol.GetNumConformers() == 0: # As a workaround, one can also build a molecule using its 2D Geometry @@ -163,19 +174,25 @@ def rd_embed(self, rdmol, num_conf_attempts): for _ in range(num_conf_attempts): # Since it can only add one conformer at a time, clearConfs is set to False # to keep the previously generated 2D geometries - AllChem.Compute2DCoords(rdmol, nSample=10, sampleSeed=1, clearConfs=False) + AllChem.Compute2DCoords( + rdmol, nSample=10, sampleSeed=1, clearConfs=False + ) if rdmol.GetNumConformers() == 0: # EDKTG + random coords + 2D should be able to cover almost all # the cases. This is to check and report if further workaround # needs to be implemented - raise RuntimeError(f'RDKit has issue in embedding the conformer of {self.molecule}.' - f' Please raise an issue on the RMG Github page mentioning the' - f' error and the molecule.') + raise RuntimeError( + f"RDKit has issue in embedding the conformer of {self.molecule}." + f" Please raise an issue on the RMG Github page mentioning the" + f" error and the molecule." + ) energy = 0.0 min_e_id = 0 - min_e = 9.999999e99 # start with a very high number, which would never be reached + min_e = ( + 9.999999e99 # start with a very high number, which would never be reached + ) crude = Chem.Mol(rdmol.ToBinary()) @@ -196,13 +213,15 @@ def rd_embed(self, rdmol, num_conf_attempts): if not good_embed and not good_opt: # TODO: Probably add something to prevent the following QM calculation if a molecule # cannot embed and force field opt correctly. - logging.debug(f"Encounted a molecule {AllChem.MolToInchi(rdmol)} that cannot be embedded" - f" and optimized correctly.") + logging.debug( + f"Encounted a molecule {AllChem.MolToInchi(rdmol)} that cannot be embedded" + f" and optimized correctly." + ) - with open(self.get_crude_mol_file_path(), 'w') as out_3d_crude: + with open(self.get_crude_mol_file_path(), "w") as out_3d_crude: out_3d_crude.write(Chem.MolToMolBlock(crude, confId=min_e_id)) - with open(self.get_refined_mol_file_path(), 'w') as out_3d: + with open(self.get_refined_mol_file_path(), "w") as out_3d: out_3d.write(Chem.MolToMolBlock(rdmol, confId=min_e_id)) return rdmol, min_e_id @@ -232,17 +251,17 @@ def load_thermo_data_file(file_path): return None try: with open(file_path) as result_file: - logging.info('Reading existing thermo file {0}'.format(file_path)) - global_context = {'__builtins__': None} + logging.info("Reading existing thermo file {0}".format(file_path)) + global_context = {"__builtins__": None} local_context = { - '__builtins__': None, - 'True': True, - 'False': False, - 'ThermoData': rmgpy.thermo.ThermoData, - 'PointGroup': symmetry.PointGroup, - 'QMData': qmdata.QMData, - 'array': np.array, - 'int32': np.int32, + "__builtins__": None, + "True": True, + "False": False, + "ThermoData": rmgpy.thermo.ThermoData, + "PointGroup": symmetry.PointGroup, + "QMData": qmdata.QMData, + "array": np.array, + "int32": np.int32, } exec(result_file.read(), global_context, local_context) except IOError: @@ -252,14 +271,20 @@ def load_thermo_data_file(file_path): logging.error('The thermo file "{0}" was invalid:'.format(file_path)) logging.exception(e) return None - if 'InChI' not in local_context: - logging.error('The thermo file "{0}" did not contain an InChI.'.format(file_path)) + if "InChI" not in local_context: + logging.error( + 'The thermo file "{0}" did not contain an InChI.'.format(file_path) + ) return None - if 'adjacencyList' not in local_context: - logging.error('The thermo file "{0}" did not contain adjacencyList.'.format(file_path)) + if "adjacencyList" not in local_context: + logging.error( + 'The thermo file "{0}" did not contain adjacencyList.'.format(file_path) + ) return None - if 'thermoData' not in local_context: - logging.error('The thermo file "{0}" did not contain thermoData.'.format(file_path)) + if "thermoData" not in local_context: + logging.error( + 'The thermo file "{0}" did not contain thermoData.'.format(file_path) + ) return None return local_context @@ -271,8 +296,8 @@ class QMMolecule(object): Specific programs and methods should inherit from this and define some extra attributes and methods: - * outputFileExtension - * inputFileExtension + * output_file_extension + * input_file_extension * generate_qm_data() ...and whatever else is needed to make this method work. The attributes are: @@ -306,22 +331,22 @@ def get_file_path(self, extension, scratch=True): # ToDo: this is duplicated in Geometry class. Should be refactored. return os.path.join( self.settings.scratchDirectory if scratch else self.settings.fileStore, - self.unique_id + extension + self.unique_id + extension, ) @property def output_file_path(self): """Get the output file name.""" - return self.get_file_path(self.outputFileExtension) + return self.get_file_path(self.output_file_extension) @property def input_file_path(self): """Get the input file name.""" - return self.get_file_path(self.inputFileExtension) + return self.get_file_path(self.input_file_extension) def get_thermo_file_path(self): """Returns the path the thermo data file.""" - return self.get_file_path('.thermo', scratch=False) + return self.get_file_path(".thermo", scratch=False) @property def script_attempts(self): @@ -350,18 +375,29 @@ def check_paths(self): """ Check the paths in the settings are OK. Make folders as necessary. """ - self.settings.fileStore = os.path.expandvars(self.settings.fileStore) # to allow things like $HOME or $RMGpy - self.settings.scratchDirectory = os.path.expandvars(self.settings.scratchDirectory) + self.settings.fileStore = os.path.expandvars( + self.settings.fileStore + ) # to allow things like $HOME or $RMGpy + self.settings.scratchDirectory = os.path.expandvars( + self.settings.scratchDirectory + ) for path in [self.settings.fileStore, self.settings.scratchDirectory]: if not os.path.exists(path): - logging.info("Creating directory %s for QM files." % os.path.abspath(path)) + logging.info( + "Creating directory %s for QM files." % os.path.abspath(path) + ) os.makedirs(path) def create_geometry(self): """ Creates self.geometry with RDKit geometries """ - self.geometry = Geometry(self.settings, self.unique_id, self.molecule, unique_id_long=self.unique_id_long) + self.geometry = Geometry( + self.settings, + self.unique_id, + self.molecule, + unique_id_long=self.unique_id_long, + ) self.geometry.generate_rdkit_geometries() return self.geometry @@ -370,21 +406,33 @@ def parse(self): Parses the results of the Mopac calculation, and returns a QMData object. """ parser = self.get_parser(self.output_file_path) - parser.logger.setLevel(logging.ERROR) # cf. http://cclib.sourceforge.net/wiki/index.php/Using_cclib#Additional_information - parser.rotcons = [] # give it an attribute and it won't delete it, leaving it on the parser object - parser.molmass = None # give it an attribute and it won't delete it, leaving it on the parser object + parser.logger.setLevel( + logging.ERROR + ) # cf. http://cclib.sourceforge.net/wiki/index.php/Using_cclib#Additional_information + parser.rotcons = ( + [] + ) # give it an attribute and it won't delete it, leaving it on the parser object + parser.molmass = None # give it an attribute and it won't delete it, leaving it on the parser object cclib_data = parser.parse() radical_number = self.molecule.get_radical_count() - cclib_data.rotcons = parser.rotcons # this hack required because rotcons not part of a default cclib data object - cclib_data.molmass = parser.molmass # this hack required because rotcons not part of a default cclib data object - qm_data = parse_cclib_data(cclib_data, radical_number + 1) # Should `radical_number+1` be `self.molecule.multiplicity` in the next line of code? It's the electronic ground state degeneracy. + cclib_data.rotcons = ( + parser.rotcons + ) # this hack required because rotcons not part of a default cclib data object + cclib_data.molmass = ( + parser.molmass + ) # this hack required because rotcons not part of a default cclib data object + qm_data = parse_cclib_data( + cclib_data, radical_number + 1 + ) # Should `radical_number+1` be `self.molecule.multiplicity` in the next line of code? It's the electronic ground state degeneracy. return qm_data def generate_qm_data(self): """ Calculate the QM data somehow and return a CCLibData object, or None if it fails. """ - raise NotImplementedError("This should be defined in a subclass that inherits from QMMolecule") + raise NotImplementedError( + "This should be defined in a subclass that inherits from QMMolecule" + ) return qmdata.QMData() or None def generate_thermo_data(self): @@ -424,15 +472,19 @@ def save_thermo_data(self): """ Save the generated thermo data. """ - self.thermo.H298.units = 'kcal/mol' - self.thermo.S298.units = 'cal/mol/K' - self.thermo.Cpdata.units = 'cal/mol/K' - with open(self.get_thermo_file_path(), 'w') as result_file: + self.thermo.H298.units = "kcal/mol" + self.thermo.S298.units = "cal/mol/K" + self.thermo.Cpdata.units = "cal/mol/K" + with open(self.get_thermo_file_path(), "w") as result_file: result_file.write('InChI = "{0!s}"\n'.format(self.unique_id_long)) result_file.write("thermoData = {0!r}\n".format(self.thermo)) result_file.write("pointGroup = {0!r}\n".format(self.point_group)) result_file.write("qmData = {0!r}\n".format(self.qm_data)) - result_file.write('adjacencyList = """\n{0!s}"""\n'.format(self.molecule.to_adjacency_list(remove_h=False))) + result_file.write( + 'adjacencyList = """\n{0!s}"""\n'.format( + self.molecule.to_adjacency_list(remove_h=False) + ) + ) def load_thermo_data(self): """ @@ -443,26 +495,34 @@ def load_thermo_data(self): if local_context is None: # file does not exist or is invalid return None - if local_context['InChI'] != self.unique_id_long: - logging.error('The InChI in the thermo file {0} did not match the ' - 'current molecule {1}'.format(file_path, self.unique_id_long)) + if local_context["InChI"] != self.unique_id_long: + logging.error( + "The InChI in the thermo file {0} did not match the " + "current molecule {1}".format(file_path, self.unique_id_long) + ) return None - loaded_molecule = rmgpy.molecule.Molecule().from_adjacency_list(local_context['adjacencyList']) + loaded_molecule = rmgpy.molecule.Molecule().from_adjacency_list( + local_context["adjacencyList"] + ) if not loaded_molecule.is_isomorphic(self.molecule): - logging.error('The adjacencyList in thermo file {0} did not match the ' - 'current molecule {1}'.format(file_path, self.unique_id_long)) + logging.error( + "The adjacencyList in thermo file {0} did not match the " + "current molecule {1}".format(file_path, self.unique_id_long) + ) return None - thermo = local_context['thermoData'] + thermo = local_context["thermoData"] assert isinstance(thermo, rmgpy.thermo.ThermoData) self.thermo = thermo - self.point_group = symmetry.point_group_dictionary[local_context['pointGroup'].point_group] # point to the one in the module level dictionary with the same name - self.qm_data = local_context['qmData'] + self.point_group = symmetry.POINT_GROUP_DICTIONARY[ + local_context["pointGroup"].point_group + ] # point to the one in the module level dictionary with the same name + self.qm_data = local_context["qmData"] return thermo def get_augmented_inchi_key(self): """ - Returns the augmented InChI from self.molecule + Returns the augmented InChI from self.molecule """ return self.molecule.to_augmented_inchi_key() @@ -495,7 +555,7 @@ def calculate_chirality_correction(self): if self.point_group.chiral: return rmgpy.quantity.constants.R * math.log(2) else: - return 0. + return 0.0 def calculate_thermo_data(self): """ @@ -507,16 +567,21 @@ def calculate_thermo_data(self): assert self.qm_data, "Need QM Data first in order to calculate thermo." assert self.point_group, "Need Point Group first in order to calculate thermo." - mass = getattr(self.qm_data, 'molecularMass', None) + mass = getattr(self.qm_data, "molecularMass", None) if mass is None: # If using a cclib that doesn't read molecular mass, for example - mass = sum(rmgpy.molecule.element.get_element(int(a)).mass for a in self.qm_data.atomicNumbers) - mass = rmgpy.quantity.Mass(mass, 'kg/mol') + mass = sum( + rmgpy.molecule.element.get_element(int(a)).mass + for a in self.qm_data.atomicNumbers + ) + mass = rmgpy.quantity.Mass(mass, "kg/mol") trans = rmgpy.statmech.IdealGasTranslation(mass=mass) if self.point_group.linear: # there should only be one rotational constant for a linear rotor - rotational_constant = rmgpy.quantity.Frequency(max(self.qm_data.rotationalConstants.value), - self.qm_data.rotationalConstants.units) + rotational_constant = rmgpy.quantity.Frequency( + max(self.qm_data.rotationalConstants.value), + self.qm_data.rotationalConstants.units, + ) rot = rmgpy.statmech.LinearRotor( rotationalConstant=rotational_constant, symmetry=self.point_group.symmetry_number, @@ -531,9 +596,11 @@ def calculate_thermo_data(self): # @todo: We need to extract or calculate E0 somehow from the qmdata E0 = (0, "kJ/mol") - self.statesmodel = rmgpy.statmech.Conformer(E0=E0, - modes=[trans, rot, vib], - spin_multiplicity=self.qm_data.groundStateDegeneracy) + self.statesmodel = rmgpy.statmech.Conformer( + E0=E0, + modes=[trans, rot, vib], + spin_multiplicity=self.qm_data.groundStateDegeneracy, + ) # we will use number of atoms from above (alternatively, we could use the chemGraph); this is needed to test whether the species is monoatomic # SI units are J/mol, but converted to kJ/mol for generating the thermo. @@ -552,7 +619,7 @@ def calculate_thermo_data(self): S298=(S298, "J/(mol*K)"), Tmin=(300.0, "K"), Tmax=(2000.0, "K"), - comment=comment + comment=comment, ) self.thermo = thermo return thermo diff --git a/rmgpy/qm/mopac.py b/rmgpy/qm/mopac.py index a0d7c29184..b941a82362 100644 --- a/rmgpy/qm/mopac.py +++ b/rmgpy/qm/mopac.py @@ -45,145 +45,121 @@ class Mopac(object): """ A base class for all QM calculations that use MOPAC. - + Classes such as :class:`MopacMol` will inherit from this class. """ - inputFileExtension = '.mop' - outputFileExtension = '.out' + input_file_extension = ".mop" + output_file_extension = ".out" - executablesToTry = ('MOPAC2016.exe', 'MOPAC2012.exe', 'MOPAC2009.exe', 'mopac') + executable_path = distutils.spawn.find_executable("mopac") - for exe in executablesToTry: - try: - executablePath = distutils.spawn.find_executable(exe) - except: - executablePath = None - if executablePath is not None: - break - else: # didn't break - logging.debug("Did not find MOPAC on path, checking if it exists in a declared MOPAC_DIR...") - mopacEnv = os.getenv('MOPAC_DIR', default="/opt/mopac") - for exe in executablesToTry: - executablePath = os.path.join(mopacEnv, exe) - if os.path.exists(executablePath): - break - else: # didn't break - executablePath = os.path.join(mopacEnv, '(MOPAC 2009 or 2012 or 2016)') - - usePolar = False # use polar keyword in MOPAC - - "Keywords for the multiplicity" - multiplicityKeywords = { - 1: '', - 2: 'uhf doublet', - 3: 'uhf triplet', - 4: 'uhf quartet', - 5: 'uhf quintet', - 6: 'uhf sextet', - 7: 'uhf septet', - 8: 'uhf octet', - 9: 'uhf nonet', + use_polar = False # use polar keyword in MOPAC + + # Keywords for the multiplicity + multiplicity_keywords = { + 1: "", + 2: "uhf doublet", + 3: "uhf triplet", + 4: "uhf quartet", + 5: "uhf quintet", + 6: "uhf sextet", + 7: "uhf septet", + 8: "uhf octet", + 9: "uhf nonet", } #: List of phrases that indicate failure #: NONE of these must be present in a succesful job. - failureKeys = [ - 'IMAGINARY FREQUENCIES', - 'EXCESS NUMBER OF OPTIMIZATION CYCLES', - 'NOT ENOUGH TIME FOR ANOTHER CYCLE', + failure_keys = [ + "IMAGINARY FREQUENCIES", + "EXCESS NUMBER OF OPTIMIZATION CYCLES", + "NOT ENOUGH TIME FOR ANOTHER CYCLE", ] #: List of phrases to indicate success. #: ALL of these must be present in a successful job. - successKeys = [ - 'DESCRIPTION OF VIBRATIONS', - 'MOPAC DONE' - ] - - def test_ready(self): - if not os.path.exists(self.executablePath): - raise DependencyError("Couldn't find MOPAC executable at {0}. Try setting your MOPAC_DIR " - "environment variable.".format(self.executablePath)) - - # Check if MOPAC executable works properly - process = Popen(self.executablePath, - stdin=PIPE, - stdout=PIPE, - stderr=PIPE) - stdout, stderr = process.communicate() - - self.expired = False - stderr = stderr.decode('utf-8') - if 'has expired' in stderr: - # The MOPAC executable is expired - logging.warning('\n'.join(stderr.split('\n')[2:7])) - self.expired = True - elif 'To install the MOPAC license' in stderr: - # The MOPAC executable exists, but the license has not been installed - raise DependencyError('\n'.join(stderr.split('\n')[0:9])) - elif 'MOPAC_LICENSE' in stderr: - # The MOPAC executable is in the wrong location on Windows; MOPAC_LICENSE must be set - raise DependencyError('\n'.join(stderr.split('\n')[0:11])) + success_keys = ["DESCRIPTION OF VIBRATIONS", "MOPAC DONE"] def run(self): - self.test_ready() + # ensure that the executable is present + if not os.path.exists(self.executable_path): + raise DependencyError( + "Couldn't find MOPAC executable at {0}. Please ensure your conda environment is " + "installed correctly.".format(self.executable_path) + ) # submits the input file to mopac - dirpath = tempfile.mkdtemp() + dir_path = tempfile.mkdtemp() # copy input file to temp dir: - tempInpFile = os.path.join(dirpath, os.path.basename(self.input_file_path)) - shutil.copy(self.input_file_path, dirpath) - - process = Popen([self.executablePath, tempInpFile], stdin=PIPE, stdout=PIPE, stderr=PIPE) - command = b'\n' if self.expired else None # press enter to pass expiration notice - stdout, stderr = process.communicate(input=command) # necessary to wait for executable termination! + temp_input_file = os.path.join(dir_path, os.path.basename(self.input_file_path)) + shutil.copy(self.input_file_path, dir_path) + + process = Popen( + [self.executable_path, temp_input_file], + stdin=PIPE, + stdout=PIPE, + stderr=PIPE, + ) + stdout, stderr = process.communicate( + input=None + ) # necessary to wait for executable termination! if b"ended normally" not in stderr.strip(): - logging.warning("Mopac error message:" + stderr.decode('utf-8')) + logging.warning("Mopac error message:" + stderr.decode("utf-8")) # copy output file from temp dir to output dir: - tempOutFile = os.path.join(dirpath, os.path.basename(self.output_file_path)) - shutil.copy(tempOutFile, self.output_file_path) + temp_output_file = os.path.join( + dir_path, os.path.basename(self.output_file_path) + ) + shutil.copy(temp_output_file, self.output_file_path) # delete temp folder: - shutil.rmtree(dirpath) + shutil.rmtree(dir_path) return self.verify_output_file() def verify_output_file(self): """ Check's that an output file exists and was successful. - - Returns a boolean flag that states whether a successful MOPAC simulation already exists for the molecule with the + + Returns a boolean flag that states whether a successful MOPAC simulation already exists for the molecule with the given (augmented) InChI Key. - + The definition of finding a successful simulation is based on these criteria: 1) finding an output file with the file name equal to the InChI Key 2) NOT finding any of the keywords that are denote a calculation failure 3) finding all the keywords that denote a calculation success. 4) finding a match between the InChI of the given molecule and the InchI found in the calculation files 5) checking that the optimized geometry, when connected by single bonds, is isomorphic with self.molecule (converted to single bonds) - + If any of the above criteria is not matched, False will be returned. If all succeed, then it will return True. """ if not os.path.exists(self.output_file_path): - logging.debug("Output file {0} does not (yet) exist.".format(self.output_file_path)) + logging.debug( + "Output file {0} does not (yet) exist.".format(self.output_file_path) + ) return False - inchi_found = False # flag (1 or 0) indicating whether an InChI was found in the log file + inchi_found = ( + False # flag (1 or 0) indicating whether an InChI was found in the log file + ) - # Initialize dictionary with "False"s - success_keys_found = dict([(key, False) for key in self.successKeys]) + # Initialize dictionary with "False"s + success_keys_found = dict([(key, False) for key in self.success_keys]) - with open(self.output_file_path) as outputFile: - for line in outputFile: + with open(self.output_file_path) as output_file: + for line in output_file: line = line.strip() - for element in self.failureKeys: # search for failure keywords + for element in self.failure_keys: # search for failure keywords if element in line: - logging.error("MOPAC output file contains the following error: {0}".format(element)) + logging.error( + "MOPAC output file contains the following error: {0}".format( + element + ) + ) return False - for element in self.successKeys: # search for success keywords + for element in self.success_keys: # search for success keywords if element in line: success_keys_found[element] = True @@ -193,22 +169,34 @@ def verify_output_file(self): if self.unique_id_long in log_file_inchi: pass elif self.unique_id_long.startswith(log_file_inchi): - logging.info("InChI too long to check, but beginning matches so assuming OK.") + logging.info( + "InChI too long to check, but beginning matches so assuming OK." + ) else: - logging.warning("InChI in log file ({0}) didn't match that in geometry " - "({1}).".format(log_file_inchi, self.unique_id_long)) + logging.warning( + "InChI in log file ({0}) didn't match that in geometry " + "({1}).".format(log_file_inchi, self.unique_id_long) + ) # Use only up to first 80 characters to match due to MOPAC bug which deletes 81st character of InChI string if self.unique_id_long.startswith(log_file_inchi[:80]): - logging.warning("but the beginning matches so it's probably just a truncation problem.") + logging.warning( + "but the beginning matches so it's probably just a truncation problem." + ) # Check that ALL 'success' keywords were found in the file. if not all(success_keys_found.values()): - logging.error('Not all of the required keywords for success were found in the output file!') + logging.error( + "Not all of the required keywords for success were found in the output file!" + ) return False if not inchi_found: - logging.error("No InChI was found in the MOPAC output file {0}".format(self.output_file_path)) + logging.error( + "No InChI was found in the MOPAC output file {0}".format( + self.output_file_path + ) + ) return False # Compare the optimized geometry to the original molecule @@ -221,10 +209,18 @@ def verify_output_file(self): return False test_mol = self.molecule.to_single_bonds() if not cclib_mol.is_isomorphic(test_mol): - logging.info("Incorrect connectivity for optimized geometry in file {0}".format(self.output_file_path)) + logging.info( + "Incorrect connectivity for optimized geometry in file {0}".format( + self.output_file_path + ) + ) return False - logging.info("Successful {1} quantum result in {0}".format(self.output_file_path, self.__class__.__name__)) + logging.info( + "Successful {1} quantum result in {0}".format( + self.output_file_path, self.__class__.__name__ + ) + ) return True def get_parser(self, output_file): @@ -236,18 +232,33 @@ def get_parser(self, output_file): class MopacMol(QMMolecule, Mopac): """ - A base Class for calculations of molecules using MOPAC. - + A base Class for calculations of molecules using MOPAC. + Inherits from both :class:`QMMolecule` and :class:`Mopac`. """ #: Keywords that will be added at the top and bottom of the qm input file keywords = [ - {'top': "precise nosym THREADS=1", 'bottom': "oldgeo thermo nosym precise THREADS=1 "}, - {'top': "precise nosym gnorm=0.0 nonr THREADS=1", 'bottom': "oldgeo thermo nosym precise THREADS=1 "}, - {'top': "precise nosym gnorm=0.0 THREADS=1", 'bottom': "oldgeo thermo nosym precise THREADS=1 "}, - {'top': "precise nosym gnorm=0.0 bfgs THREADS=1", 'bottom': "oldgeo thermo nosym precise THREADS=1 "}, - {'top': "precise nosym recalc=10 dmax=0.10 nonr cycles=2000 t=2000 THREADS=1", 'bottom': "oldgeo thermo nosym precise THREADS=1 "}, + { + "top": "precise nosym THREADS=1", + "bottom": "oldgeo thermo nosym precise THREADS=1 ", + }, + { + "top": "precise nosym gnorm=0.0 nonr THREADS=1", + "bottom": "oldgeo thermo nosym precise THREADS=1 ", + }, + { + "top": "precise nosym gnorm=0.0 THREADS=1", + "bottom": "oldgeo thermo nosym precise THREADS=1 ", + }, + { + "top": "precise nosym gnorm=0.0 bfgs THREADS=1", + "bottom": "oldgeo thermo nosym precise THREADS=1 ", + }, + { + "top": "precise nosym recalc=10 dmax=0.10 nonr cycles=2000 t=2000 THREADS=1", + "bottom": "oldgeo thermo nosym precise THREADS=1 ", + }, ] def write_input_file(self, attempt): @@ -257,32 +268,40 @@ def write_input_file(self, attempt): """ molfile = self.get_mol_file_path_for_calculation(attempt) - atomline = re.compile(r'\s*([\- ][0-9.]+)\s+([\- ][0-9.]+)+\s+([\- ][0-9.]+)\s+([A-Za-z]+)') + atomline = re.compile( + r"\s*([\- ][0-9.]+)\s+([\- ][0-9.]+)+\s+([\- ][0-9.]+)\s+([A-Za-z]+)" + ) - output = [self.geometry.unique_id_long, ''] + output = [self.geometry.unique_id_long, ""] atom_count = 0 with open(molfile) as molinput: for line in molinput: match = atomline.match(line) if match: - output.append("{0:4s} {1} 1 {2} 1 {3} 1".format(match.group(4), match.group(1), - match.group(2), match.group(3))) + output.append( + "{0:4s} {1} 1 {2} 1 {3} 1".format( + match.group(4), + match.group(1), + match.group(2), + match.group(3), + ) + ) atom_count += 1 assert atom_count == len(self.molecule.atoms) - output.append('') - input_string = '\n'.join(output) + output.append("") + input_string = "\n".join(output) top_keys, bottom_keys, polar_keys = self.input_file_keywords(attempt) - with open(self.input_file_path, 'w') as mopac_file: + with open(self.input_file_path, "w") as mopac_file: mopac_file.write(top_keys) - mopac_file.write('\n') + mopac_file.write("\n") mopac_file.write(input_string) - mopac_file.write('\n') + mopac_file.write("\n") mopac_file.write(bottom_keys) - if self.usePolar: - mopac_file.write('\n\n\n') + if self.use_polar: + mopac_file.write("\n\n\n") mopac_file.write(polar_keys) def input_file_keywords(self, attempt): @@ -301,23 +320,41 @@ def generate_qm_data(self): if self.verify_output_file(): logging.info("Found a successful output file already; using that.") - source = "QM {0} calculation found from previous run.".format(self.__class__.__name__) + source = "QM {0} calculation found from previous run.".format( + self.__class__.__name__ + ) else: self.create_geometry() success = False for attempt in range(1, self.max_attempts + 1): self.write_input_file(attempt) - logging.info('Trying {3} attempt {0} of {1} on molecule {2}.'.format(attempt, self.max_attempts, - self.molecule.to_smiles(), - self.__class__.__name__)) + logging.info( + "Trying {3} attempt {0} of {1} on molecule {2}.".format( + attempt, + self.max_attempts, + self.molecule.to_smiles(), + self.__class__.__name__, + ) + ) success = self.run() if success: - logging.info('Attempt {0} of {1} on species {2} succeeded.'.format(attempt, self.max_attempts, - self.molecule.to_augmented_inchi())) - source = "QM {0} calculation attempt {1}".format(self.__class__.__name__, attempt) + logging.info( + "Attempt {0} of {1} on species {2} succeeded.".format( + attempt, + self.max_attempts, + self.molecule.to_augmented_inchi(), + ) + ) + source = "QM {0} calculation attempt {1}".format( + self.__class__.__name__, attempt + ) break else: - logging.error('QM thermo calculation failed for {0}.'.format(self.molecule.to_augmented_inchi())) + logging.error( + "QM thermo calculation failed for {0}.".format( + self.molecule.to_augmented_inchi() + ) + ) return None result = self.parse() # parsed in cclib result.source = source @@ -327,17 +364,18 @@ def generate_qm_data(self): class MopacMolPMn(MopacMol): """ Mopac PMn calculations for molecules (n undefined here) - + This is a parent class for MOPAC PMn calculations. - Inherit it, and define the pm_method, then redefine + Inherit it, and define the pm_method, then redefine anything you wish to do differently. """ - pm_method = '(should be defined by sub class)' + + pm_method = "(should be defined by sub class)" def input_file_keywords(self, attempt): """ Return the top, bottom, and polar keywords for attempt number `attempt`. - + NB. `attempt` begins at 1, not 0. """ assert attempt <= self.max_attempts @@ -345,21 +383,23 @@ def input_file_keywords(self, attempt): if attempt > self.script_attempts: attempt -= self.script_attempts - multiplicity_keys = self.multiplicityKeywords[self.geometry.molecule.multiplicity] + multiplicity_keys = self.multiplicity_keywords[ + self.geometry.molecule.multiplicity + ] top_keys = "{method} {mult} {top}".format( method=self.pm_method, mult=multiplicity_keys, - top=self.keywords[attempt - 1]['top'], + top=self.keywords[attempt - 1]["top"], ) bottom_keys = "{bottom} {method} {mult}".format( method=self.pm_method, - bottom=self.keywords[attempt - 1]['bottom'], + bottom=self.keywords[attempt - 1]["bottom"], mult=multiplicity_keys, ) polar_keys = "oldgeo {polar} nosym precise {method} {mult}".format( method=self.pm_method, - polar=('polar' if self.geometry.molecule.multiplicity == 1 else 'static'), + polar=("polar" if self.geometry.molecule.multiplicity == 1 else "static"), mult=multiplicity_keys, ) @@ -369,28 +409,31 @@ def input_file_keywords(self, attempt): class MopacMolPM3(MopacMolPMn): """ Mopac PM3 calculations for molecules - + This is a class of its own in case you wish to do anything differently, but for now it's the same as all the MOPAC PMn calculations, only pm3 """ - pm_method = 'pm3' + + pm_method = "pm3" class MopacMolPM6(MopacMolPMn): """ Mopac PM6 calculations for molecules - + This is a class of its own in case you wish to do anything differently, but for now it's the same as all the MOPAC PMn calculations, only pm6 """ - pm_method = 'pm6' + + pm_method = "pm6" class MopacMolPM7(MopacMolPMn): """ Mopac PM7 calculations for molecules - + This is a class of its own in case you wish to do anything differently, but for now it's the same as all the MOPAC PMn calculations, only pm7 """ - pm_method = 'pm7' + + pm_method = "pm7" diff --git a/rmgpy/qm/mopacTest.py b/rmgpy/qm/mopacTest.py index 6343724c47..7ada11f673 100644 --- a/rmgpy/qm/mopacTest.py +++ b/rmgpy/qm/mopacTest.py @@ -39,18 +39,8 @@ from rmgpy.qm.main import QMCalculator from rmgpy.qm.mopac import Mopac, MopacMolPM3, MopacMolPM6, MopacMolPM7 -NO_MOPAC = NO_LICENCE = False -try: - Mopac().test_ready() -except DependencyError as e: - if "Couldn't find MOPAC executable" in str(e): - NO_MOPAC = NO_LICENCE = True - elif 'To install the MOPAC license' in str(e) or 'MOPAC_LICENSE' in str(e): - NO_LICENCE = True - else: - raise - -mol1 = Molecule().from_smiles('C1=CC=C2C=CC=CC2=C1') +mol1 = Molecule().from_smiles("C1=CC=C2C=CC=CC2=C1") +MOPAC_CLOSE_ENOUGH_PERCENT = 0.005 # 0.5% class TestMopacMolPM3(unittest.TestCase): @@ -58,19 +48,18 @@ class TestMopacMolPM3(unittest.TestCase): Contains unit tests for the Geometry class. """ - @unittest.skipIf(NO_MOPAC, "MOPAC not found. Try resetting your environment variables if you want to use it.") - @unittest.skipIf(NO_LICENCE, "MOPAC license not installed. Run mopac for instructions") def setUp(self): """ A function run before each unit test in this class. """ - rmgpy_path = os.path.normpath(os.path.join(get_path(), '..')) + rmgpy_path = os.path.normpath(os.path.join(get_path(), "..")) - qm = QMCalculator(software='mopac', - method='pm3', - fileStore=os.path.join(rmgpy_path, 'testing', 'qm', 'QMfiles'), - scratchDirectory=os.path.join(rmgpy_path, 'testing', 'qm', 'QMscratch'), - ) + qm = QMCalculator( + software="mopac", + method="pm3", + fileStore=os.path.join(rmgpy_path, "testing", "qm", "QMfiles"), + scratchDirectory=os.path.join(rmgpy_path, "testing", "qm", "QMscratch"), + ) if not os.path.exists(qm.settings.fileStore): os.makedirs(qm.settings.fileStore) @@ -82,39 +71,62 @@ def test_generate_thermo_data(self): Test that generate_thermo_data() works correctly for MOPAC PM3 """ # First ensure any old data are removed, or else they'll be reused! - for directory in (self.qmmol1.settings.fileStore, self.qmmol1.settings.scratchDirectory): + for directory in ( + self.qmmol1.settings.fileStore, + self.qmmol1.settings.scratchDirectory, + ): shutil.rmtree(directory, ignore_errors=True) self.qmmol1.generate_thermo_data() result = self.qmmol1.qm_data self.assertTrue(self.qmmol1.verify_output_file()) - self.assertTrue(self.qmmol1.thermo.comment.startswith('QM MopacMolPM3 calculation')) + self.assertTrue( + self.qmmol1.thermo.comment.startswith("QM MopacMolPM3 calculation") + ) self.assertEqual(result.numberOfAtoms, 18) self.assertIsInstance(result.atomicNumbers, np.ndarray) - if result.molecularMass.units == 'amu': + if result.molecularMass.units == "amu": self.assertAlmostEqual(result.molecularMass.value, 128.173, 2) - self.assertAlmostEqual(self.qmmol1.thermo.H298.value_si, 169708.0608, 0) # to 1 decimal place - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 334.5007584, 1) # to 1 decimal place + self.assertAlmostEqual( + self.qmmol1.thermo.H298.value_si, + 169708, + delta=169708 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) + self.assertAlmostEqual( + self.qmmol1.thermo.S298.value_si, + 334.500, + delta=334.500 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) def test_load_thermo_data(self): """ Test that generate_thermo_data() can load thermo from the previous MOPAC PM3 run. - + Check that it loaded, and the values are the same as above. """ self.qmmol1.generate_thermo_data() result = self.qmmol1.qm_data - self.assertTrue(self.qmmol1.thermo.comment.startswith('QM MopacMolPM3 calculation')) + self.assertTrue( + self.qmmol1.thermo.comment.startswith("QM MopacMolPM3 calculation") + ) self.assertEqual(result.numberOfAtoms, 18) self.assertIsInstance(result.atomicNumbers, np.ndarray) - if result.molecularMass.units == 'amu': + if result.molecularMass.units == "amu": self.assertAlmostEqual(result.molecularMass.value, 128.173, 2) - self.assertAlmostEqual(self.qmmol1.thermo.H298.value_si, 169708.0608, 0) # to 1 decimal place - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 334.5007584, 1) # to 1 decimal place + self.assertAlmostEqual( + self.qmmol1.thermo.H298.value_si, + 169708, + delta=169708 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) + self.assertAlmostEqual( + self.qmmol1.thermo.S298.value_si, + 334.500, + delta=334.500 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) class TestMopacMolPM6(unittest.TestCase): @@ -122,19 +134,18 @@ class TestMopacMolPM6(unittest.TestCase): Contains unit tests for the Geometry class. """ - @unittest.skipIf(NO_MOPAC, "MOPAC not found. Try resetting your environment variables if you want to use it.") - @unittest.skipIf(NO_LICENCE, "MOPAC license not installed. Run mopac for instructions") def setUp(self): """ A function run before each unit test in this class. """ - rmgpy_path = os.path.normpath(os.path.join(get_path(), '..')) + rmgpy_path = os.path.normpath(os.path.join(get_path(), "..")) - qm = QMCalculator(software='mopac', - method='pm6', - fileStore=os.path.join(rmgpy_path, 'testing', 'qm', 'QMfiles'), - scratchDirectory=os.path.join(rmgpy_path, 'testing', 'qm', 'QMscratch'), - ) + qm = QMCalculator( + software="mopac", + method="pm6", + fileStore=os.path.join(rmgpy_path, "testing", "qm", "QMfiles"), + scratchDirectory=os.path.join(rmgpy_path, "testing", "qm", "QMscratch"), + ) if not os.path.exists(qm.settings.fileStore): os.makedirs(qm.settings.fileStore) @@ -146,40 +157,63 @@ def test_generate_thermo_data(self): Test that generate_thermo_data() works correctly for MOPAC PM6 """ # First ensure any old data are removed, or else they'll be reused! - for directory in (self.qmmol1.settings.fileStore, self.qmmol1.settings.scratchDirectory): + for directory in ( + self.qmmol1.settings.fileStore, + self.qmmol1.settings.scratchDirectory, + ): shutil.rmtree(directory, ignore_errors=True) self.qmmol1.generate_thermo_data() result = self.qmmol1.qm_data self.assertTrue(self.qmmol1.verify_output_file()) - self.assertTrue(self.qmmol1.thermo.comment.startswith('QM MopacMolPM6 calculation')) + self.assertTrue( + self.qmmol1.thermo.comment.startswith("QM MopacMolPM6 calculation") + ) self.assertEqual(result.numberOfAtoms, 18) self.assertIsInstance(result.atomicNumbers, np.ndarray) - if result.molecularMass.units == 'amu': + if result.molecularMass.units == "amu": self.assertAlmostEqual(result.molecularMass.value, 128.173, 2) - self.assertAlmostEqual(self.qmmol1.thermo.H298.value_si, 167704.4270, 0) # to 1 decimal place - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 338.0999241, 1) # to 1 decimal place + self.assertAlmostEqual( + self.qmmol1.thermo.H298.value_si, + 167704, + delta=167704 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) + self.assertAlmostEqual( + self.qmmol1.thermo.S298.value_si, + 338.099, + delta=338.099 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) def test_load_thermo_data(self): """ Test that generate_thermo_data() can load thermo from the previous MOPAC PM6 run. - + Check that it loaded, and the values are the same as above. """ self.qmmol1.generate_thermo_data() result = self.qmmol1.qm_data - self.assertTrue(self.qmmol1.thermo.comment.startswith('QM MopacMolPM6 calculation')) + self.assertTrue( + self.qmmol1.thermo.comment.startswith("QM MopacMolPM6 calculation") + ) self.assertEqual(result.numberOfAtoms, 18) self.assertIsInstance(result.atomicNumbers, np.ndarray) - if result.molecularMass.units == 'amu': + if result.molecularMass.units == "amu": self.assertEqual(result.molecularMass.value, 128.173) - self.assertAlmostEqual(self.qmmol1.thermo.H298.value_si, 167704.0681, 0) # to 0 decimal place - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 338.0999241, 1) # to 1 decimal place + self.assertAlmostEqual( + self.qmmol1.thermo.H298.value_si, + 167704, + delta=167704 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) # to 0 decimal place + self.assertAlmostEqual( + self.qmmol1.thermo.S298.value_si, + 338.099, + delta=338.099 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) class TestMopacMolPM7(unittest.TestCase): @@ -187,24 +221,23 @@ class TestMopacMolPM7(unittest.TestCase): Contains unit tests for the Geometry class. """ - @unittest.skipIf(NO_MOPAC, "MOPAC not found. Try resetting your environment variables if you want to use it.") - @unittest.skipIf(NO_LICENCE, "MOPAC license not installed. Run mopac for instructions") def setUp(self): """ A function run before each unit test in this class. """ - rmgpy_path = os.path.normpath(os.path.join(get_path(), '..')) + rmgpy_path = os.path.normpath(os.path.join(get_path(), "..")) - qm = QMCalculator(software='mopac', - method='pm7', - fileStore=os.path.join(rmgpy_path, 'testing', 'qm', 'QMfiles'), - scratchDirectory=os.path.join(rmgpy_path, 'testing', 'qm', 'QMscratch'), - ) + qm = QMCalculator( + software="mopac", + method="pm7", + fileStore=os.path.join(rmgpy_path, "testing", "qm", "QMfiles"), + scratchDirectory=os.path.join(rmgpy_path, "testing", "qm", "QMscratch"), + ) if not os.path.exists(qm.settings.fileStore): os.makedirs(qm.settings.fileStore) - mol1 = Molecule().from_smiles('C1=CC=C2C=CC=CC2=C1') + mol1 = Molecule().from_smiles("C1=CC=C2C=CC=CC2=C1") self.qmmol1 = MopacMolPM7(mol1, qm.settings) def test_generate_thermo_data(self): @@ -212,43 +245,66 @@ def test_generate_thermo_data(self): Test that generate_thermo_data() works correctly for MOPAC PM7 """ # First ensure any old data are removed, or else they'll be reused! - for directory in (self.qmmol1.settings.fileStore, self.qmmol1.settings.scratchDirectory): + for directory in ( + self.qmmol1.settings.fileStore, + self.qmmol1.settings.scratchDirectory, + ): shutil.rmtree(directory, ignore_errors=True) self.qmmol1.generate_thermo_data() result = self.qmmol1.qm_data self.assertTrue(self.qmmol1.verify_output_file()) - self.assertTrue(self.qmmol1.thermo.comment.startswith('QM MopacMolPM7 calculation')) + self.assertTrue( + self.qmmol1.thermo.comment.startswith("QM MopacMolPM7 calculation") + ) self.assertEqual(result.numberOfAtoms, 18) self.assertIsInstance(result.atomicNumbers, np.ndarray) - if result.molecularMass.units == 'amu': + if result.molecularMass.units == "amu": self.assertAlmostEqual(result.molecularMass.value, 128.173, 2) - self.assertAlmostEqual(self.qmmol1.thermo.H298.value_si, 166168.9863, 0) # to 1 decimal place - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 336.3330406, 1) # to 1 decimal place + self.assertAlmostEqual( + self.qmmol1.thermo.H298.value_si, + 166169, + delta=166169 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) + self.assertAlmostEqual( + self.qmmol1.thermo.S298.value_si, + 336.333, + delta=336.333 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) def test_load_thermo_data(self): """ Test that generate_thermo_data() can load thermo from the previous MOPAC PM7 run. - + Check that it loaded, and the values are the same as above. """ self.qmmol1.generate_thermo_data() result = self.qmmol1.qm_data - self.assertTrue(self.qmmol1.thermo.comment.startswith('QM MopacMolPM7 calculation')) + self.assertTrue( + self.qmmol1.thermo.comment.startswith("QM MopacMolPM7 calculation") + ) self.assertEqual(result.numberOfAtoms, 18) self.assertIsInstance(result.atomicNumbers, np.ndarray) - if result.molecularMass.units == 'amu': + if result.molecularMass.units == "amu": self.assertAlmostEqual(result.molecularMass.value, 128.173, 2) - self.assertAlmostEqual(self.qmmol1.thermo.H298.value_si, 166168.8571, 0) # to 1 decimal place - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 336.3330406, 1) # to 1 decimal place + self.assertAlmostEqual( + self.qmmol1.thermo.H298.value_si, + 166169, + delta=166169 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) + self.assertAlmostEqual( + self.qmmol1.thermo.S298.value_si, + 336.333, + delta=336.333 * MOPAC_CLOSE_ENOUGH_PERCENT, + ) ################################################################################ -if __name__ == '__main__': +if __name__ == "__main__": unittest.main(testRunner=unittest.TextTestRunner(verbosity=2)) diff --git a/rmgpy/qm/symmetry.py b/rmgpy/qm/symmetry.py index ab4d0c702b..8fb2d6cdf2 100644 --- a/rmgpy/qm/symmetry.py +++ b/rmgpy/qm/symmetry.py @@ -39,17 +39,16 @@ class PointGroup(object): """ A symmetry Point Group. - + Attributes are: - + * point_group * symmetry_number * chiral - * linear + * linear """ def __init__(self, point_group, symmetry_number, chiral): - self.point_group = point_group self.symmetry_number = symmetry_number self.chiral = chiral @@ -61,18 +60,19 @@ def __init__(self, point_group, symmetry_number, chiral): self.linear = False def __repr__(self): - return 'PointGroup("{0}", symmetry_number={1}, chiral={2})'.format(self.point_group, self.symmetry_number, - self.chiral) + return 'PointGroup("{0}", symmetry_number={1}, chiral={2})'.format( + self.point_group, self.symmetry_number, self.chiral + ) -def make_point_group_dictionary(): +def _make_point_group_dictionary(): """ A function to make and fill the point group dictionary. - + This will be stored once as a module level variable. """ point_group_dictionary = dict() - for (point_group, symmetry_number, chiral) in [ + for point_group, symmetry_number, chiral in [ ("C1", 1, True), ("Cs", 1, False), ("Ci", 1, False), @@ -83,7 +83,6 @@ def make_point_group_dictionary(): ("C6", 6, True), ("C7", 7, True), ("C8", 8, True), - ("D2", 4, True), ("D3", 6, True), ("D4", 8, True), @@ -91,7 +90,6 @@ def make_point_group_dictionary(): ("D6", 12, True), ("D7", 14, True), ("D8", 16, True), - ("C2v", 2, False), ("C3v", 3, False), ("C4v", 4, False), @@ -99,14 +97,12 @@ def make_point_group_dictionary(): ("C6v", 6, False), ("C7v", 7, False), ("C8v", 8, False), - ("C2h", 2, False), ("C3h", 3, False), ("C4h", 4, False), ("C5h", 5, False), ("C6h", 6, False), ("C8h", 8, False), - ("D2h", 4, False), ("D3h", 6, False), ("D4h", 8, False), @@ -114,7 +110,6 @@ def make_point_group_dictionary(): ("D6h", 12, False), ("D7h", 14, False), ("D8h", 16, False), - ("D2d", 4, False), ("D3d", 6, False), ("D4d", 8, False), @@ -122,30 +117,28 @@ def make_point_group_dictionary(): ("D6d", 12, False), ("D7d", 14, False), ("D8d", 16, False), - ("S4", 2, True), ("S6", 3, True), ("S8", 4, True), - ("T", 12, True), ("Th", 12, False), ("Td", 12, False), - ("O", 24, True), ("Oh", 24, False), - ("Cinfv", 1, False), ("Dinfh", 2, False), ("I", 60, True), ("Ih", 60, False), ("Kh", 1, False), ]: - point_group_dictionary[point_group] = PointGroup(point_group, symmetry_number, chiral) + point_group_dictionary[point_group] = PointGroup( + point_group, symmetry_number, chiral + ) return point_group_dictionary #: A dictionary of PointGroup objects, stored as a module level variable. -point_group_dictionary = make_point_group_dictionary() +POINT_GROUP_DICTIONARY = _make_point_group_dictionary() class PointGroupCalculator(object): @@ -166,12 +159,12 @@ def calculate(self): class SymmetryJob(object): """ - Determine the point group using the SYMMETRY program - + Determine the point group using the SYMMETRY program + (Originally ``http://www.cobalt.chem.ucalgary.ca/ps/symmetry/`` now mirrored at https://github.com/nquesada/symmetry). - Required input is a line with number of atoms followed by lines for each atom + Required input is a line with number of atoms followed by lines for each atom including: 1) atom number 2) x,y,z coordinates @@ -181,15 +174,15 @@ class SymmetryJob(object): """ - 'Arguments that will be passed as an argument for the consecutive attempts' + "Arguments that will be passed as an argument for the consecutive attempts" argumentsList = [ - ['-final', '0.02'], - ['-final', '0.1'], - ['-primary', '0.2', '-final', '0.1'], - ['-final', '0.0'], + ["-final", "0.02"], + ["-final", "0.1"], + ["-primary", "0.2", "-final", "0.1"], + ["-final", "0.0"], ] - inputFileExtension = '.symm' + input_file_extension = ".symm" def __init__(self, settings, unique_id, qm_data): self.settings = settings @@ -203,18 +196,24 @@ def __init__(self, settings, unique_id, qm_data): @property def input_file_path(self): """The input file's path""" - return os.path.join(self.settings.scratchDirectory, self.unique_id + self.inputFileExtension) + return os.path.join( + self.settings.scratchDirectory, self.unique_id + self.input_file_extension + ) def parse(self, output): """ Check the `output` string and extract the resulting point group, which is returned. """ - for line in output.split('\n'): - if line.startswith("It seems to be the "): # "It seems to be the [x] point group" indicates point group. + for line in output.split("\n"): + if line.startswith( + "It seems to be the " + ): # "It seems to be the [x] point group" indicates point group. result = line.split(" ")[5] break else: - logging.exception("Couldn't find point group from symmetry output:\n{0}".format(output)) + logging.exception( + "Couldn't find point group from symmetry output:\n{0}".format(output) + ) return "Not found" logging.info("Point group: " + result) @@ -228,13 +227,15 @@ def run(self, command): pp = Popen(command, stdout=PIPE, stderr=PIPE) except OSError as e: logging.error(e) - raise Exception('Running symmetry on the point group calculation has failed. Please check if symmetry ' - 'program is installed on your system in RMG-Py/bin/symmetry or on your path.') + raise Exception( + "Running symmetry on the point group calculation has failed. Please check if symmetry " + "program is installed on your system in RMG-Py/bin/symmetry or on your path." + ) stdout, stderr = pp.communicate() if stderr: logging.error("Error message from SYMMETRY calculation:") logging.error(stderr) - return stdout.decode('utf-8') + return stdout.decode("utf-8") def write_input_file(self): """ @@ -243,12 +244,19 @@ def write_input_file(self): geom = str(self.qm_data.numberOfAtoms) + "\n" coords_in_angstrom = self.qm_data.atomCoords.value_si * 1e10 for i in range(self.qm_data.numberOfAtoms): - geom = geom + " ".join((str(self.qm_data.atomicNumbers[i]), - str(coords_in_angstrom[i][0]), - str(coords_in_angstrom[i][1]), - str(coords_in_angstrom[i][2]) - )) + "\n" - with open(self.input_file_path, 'w') as input_file: + geom = ( + geom + + " ".join( + ( + str(self.qm_data.atomicNumbers[i]), + str(coords_in_angstrom[i][0]), + str(coords_in_angstrom[i][1]), + str(coords_in_angstrom[i][2]), + ) + ) + + "\n" + ) + with open(self.input_file_path, "w") as input_file: input_file.write(geom) input_file.close() logging.info("Symmetry input file written to {0}".format(self.input_file_path)) @@ -277,13 +285,19 @@ def calculate(self): # parse the output to get a point group name point_group_name = self.parse(output) - if point_group_name in point_group_dictionary: + if point_group_name in POINT_GROUP_DICTIONARY: self.point_group_found = True - return point_group_dictionary[point_group_name] + return POINT_GROUP_DICTIONARY[point_group_name] else: - logging.info("Attempt number {0} did not identify a recognized " - "point group ({1}).".format(attempt, point_group_name)) + logging.info( + "Attempt number {0} did not identify a recognized " + "point group ({1}).".format(attempt, point_group_name) + ) if attempt + 2 == len(self.argumentsList): - logging.warning("Using last-resort symmetry estimation options; symmetry may be underestimated.") - logging.critical("Final attempt did not identify a recognized point group. Exiting.") + logging.warning( + "Using last-resort symmetry estimation options; symmetry may be underestimated." + ) + logging.critical( + "Final attempt did not identify a recognized point group. Exiting." + ) return None