From a88238bd0434dd577bd199bfc80684c80da16ad2 Mon Sep 17 00:00:00 2001 From: Jackson Burns Date: Thu, 20 Apr 2023 18:13:03 -0400 Subject: [PATCH 01/14] switch to conda-force mopac in environment.yml, skip install step in CI --- .github/workflows/CI.yml | 10 ---------- environment.yml | 2 +- 2 files changed, 1 insertion(+), 11 deletions(-) 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/environment.yml b/environment.yml index da98a778f1..cbc103154d 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 From 091574522d966d387dccd3e05762bf53e1d9c469 Mon Sep 17 00:00:00 2001 From: Jackson Burns Date: Thu, 20 Apr 2023 18:47:33 -0400 Subject: [PATCH 02/14] add libfortran for mopac, float cantera to potentially unsupported ver. this might end up being blocked by the same PRs as before --- environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/environment.yml b/environment.yml index cbc103154d..0c0f881ce2 100644 --- a/environment.yml +++ b/environment.yml @@ -24,6 +24,7 @@ dependencies: - markupsafe - matplotlib >=1.5 - conda-forge::mopac + - conda-forge::libgfortran-ng>=10.0 - mpmath - rmg::muq2 - networkx From 77be5e502824f6e055b71ebbcd788ff63698a84d Mon Sep 17 00:00:00 2001 From: Jackson Burns <33505528+JacksonBurns@users.noreply.github.com> Date: Tue, 2 May 2023 19:06:47 -0400 Subject: [PATCH 03/14] Use newer pydas and pydqed. I have rebuilt pydas and pydqd in an environment with an updated version of libgfortran-ng in hopes that pydas can be resolved with a modern version of MOPAC. this process will need to be repeated for pydqed --- environment.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/environment.yml b/environment.yml index 0c0f881ce2..50ac8f3394 100644 --- a/environment.yml +++ b/environment.yml @@ -35,9 +35,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 From f38d45aa8e678dd1789b731c597e0ad10bf187cf Mon Sep 17 00:00:00 2001 From: JacksonBurns Date: Mon, 8 May 2023 17:51:52 -0400 Subject: [PATCH 04/14] change mopac unit tests - allowed difference is now 0.01% instead of zero decimal places (we can't expect that between different versions anyway) --- rmgpy/qm/mainTest.py | 15 +++++++++------ rmgpy/qm/mopacTest.py | 25 +++++++++++++------------ 2 files changed, 22 insertions(+), 18 deletions(-) diff --git a/rmgpy/qm/mainTest.py b/rmgpy/qm/mainTest.py index e1b408f28d..3dfd67d748 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.0001 # 0.01% + + class TestQMSettings(unittest.TestCase): """ Contains unit tests for the QMSettings class. @@ -249,12 +252,12 @@ def test_get_thermo_data_mopac(self): 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 + 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): diff --git a/rmgpy/qm/mopacTest.py b/rmgpy/qm/mopacTest.py index 6343724c47..64def404ae 100644 --- a/rmgpy/qm/mopacTest.py +++ b/rmgpy/qm/mopacTest.py @@ -51,6 +51,7 @@ raise mol1 = Molecule().from_smiles('C1=CC=C2C=CC=CC2=C1') +MOPAC_CLOSE_ENOUGH_PERCENT = 0.0001 # 0.01% class TestMopacMolPM3(unittest.TestCase): @@ -94,8 +95,8 @@ def test_generate_thermo_data(self): 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): """ @@ -113,8 +114,8 @@ def test_load_thermo_data(self): 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): @@ -159,8 +160,8 @@ def test_generate_thermo_data(self): 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): """ @@ -178,8 +179,8 @@ def test_load_thermo_data(self): 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, selta=338.099*MOPAC_CLOSE_ENOUGH_PERCENT) class TestMopacMolPM7(unittest.TestCase): @@ -225,8 +226,8 @@ def test_generate_thermo_data(self): 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): """ @@ -244,8 +245,8 @@ def test_load_thermo_data(self): 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) ################################################################################ From 2068326f74b108636f462ff15eed2c30fd77fa79 Mon Sep 17 00:00:00 2001 From: JacksonBurns Date: Mon, 8 May 2023 18:04:21 -0400 Subject: [PATCH 05/14] update documentation regarding use of MOPAC --- .../source/users/rmg/installation/anacondaDeveloper.rst | 4 ---- documentation/source/users/rmg/installation/anacondaUser.rst | 4 ---- documentation/source/users/rmg/installation/dependencies.rst | 3 +-- 3 files changed, 1 insertion(+), 10 deletions(-) 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. From 0fe291cc4a55722f340caf824ec24948098a84f3 Mon Sep 17 00:00:00 2001 From: JacksonBurns Date: Mon, 8 May 2023 18:05:55 -0400 Subject: [PATCH 06/14] delete NO_LICENSE and NO_MOPAC checks since it will now be included always I think the executable finding thing in mopac.py is going to fail --- rmgpy/qm/mainTest.py | 14 -------------- rmgpy/qm/mopac.py | 32 ++++++-------------------------- rmgpy/qm/mopacTest.py | 17 ----------------- 3 files changed, 6 insertions(+), 57 deletions(-) diff --git a/rmgpy/qm/mainTest.py b/rmgpy/qm/mainTest.py index 3dfd67d748..ab1474a7f2 100644 --- a/rmgpy/qm/mainTest.py +++ b/rmgpy/qm/mainTest.py @@ -83,16 +83,6 @@ class TestQMCalculator(unittest.TestCase): """ 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) @@ -222,8 +212,6 @@ def test_get_thermo_data(self): 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. @@ -287,8 +275,6 @@ def test_get_thermo_data_gaussian(self): 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', diff --git a/rmgpy/qm/mopac.py b/rmgpy/qm/mopac.py index a0d7c29184..7bd71005ce 100644 --- a/rmgpy/qm/mopac.py +++ b/rmgpy/qm/mopac.py @@ -100,33 +100,13 @@ class Mopac(object): '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])) - def run(self): - self.test_ready() + # ensure that the executable is present + if not os.path.exists(self.executablePath): + raise DependencyError( + "Couldn't find MOPAC executable at {0}. Please ensure your conda environment is " + "installed correctly.".format(self.executablePath) + ) # submits the input file to mopac dirpath = tempfile.mkdtemp() diff --git a/rmgpy/qm/mopacTest.py b/rmgpy/qm/mopacTest.py index 64def404ae..ec85b45d0b 100644 --- a/rmgpy/qm/mopacTest.py +++ b/rmgpy/qm/mopacTest.py @@ -39,17 +39,6 @@ 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') MOPAC_CLOSE_ENOUGH_PERCENT = 0.0001 # 0.01% @@ -59,8 +48,6 @@ 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. @@ -123,8 +110,6 @@ 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. @@ -188,8 +173,6 @@ 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. From ebea50e0e65e7e47c320668f816a81cb3f2eb8d2 Mon Sep 17 00:00:00 2001 From: JacksonBurns <33505528+JacksonBurns@users.noreply.github.com> Date: Mon, 8 May 2023 21:01:13 -0400 Subject: [PATCH 07/14] format mopac.py with black, remove checks for old unsupported versions --- rmgpy/qm/mopac.py | 261 +++++++++++++++++++++++++++++----------------- 1 file changed, 164 insertions(+), 97 deletions(-) diff --git a/rmgpy/qm/mopac.py b/rmgpy/qm/mopac.py index 7bd71005ce..2f7520cddb 100644 --- a/rmgpy/qm/mopac.py +++ b/rmgpy/qm/mopac.py @@ -45,60 +45,49 @@ 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' - - executablesToTry = ('MOPAC2016.exe', 'MOPAC2012.exe', 'MOPAC2009.exe', 'mopac') + inputFileExtension = ".mop" + outputFileExtension = ".out" - 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)') + try: + executablePath = distutils.spawn.find_executable("mopac") + except: + 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") + executablePath = os.path.join(mopacEnv, "mopac") + if not os.path.exists(executablePath): + executablePath = os.path.join(mopacEnv, "(MOPAC 2009 or 2012 or 2016)") usePolar = False # use polar keyword in MOPAC - "Keywords for the multiplicity" + # 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', + 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', + "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' - ] + successKeys = ["DESCRIPTION OF VIBRATIONS", "MOPAC DONE"] def run(self): # ensure that the executable is present @@ -114,11 +103,14 @@ def run(self): 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! + process = Popen( + [self.executablePath, tempInpFile], 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)) @@ -131,27 +123,31 @@ 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 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 + # Initialize dictionary with "False"s success_keys_found = dict([(key, False) for key in self.successKeys]) with open(self.output_file_path) as outputFile: @@ -160,7 +156,11 @@ def verify_output_file(self): for element in self.failureKeys: # 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 @@ -173,22 +173,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 @@ -201,10 +213,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): @@ -216,18 +236,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): @@ -237,32 +272,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') + mopac_file.write("\n\n\n") mopac_file.write(polar_keys) def input_file_keywords(self, attempt): @@ -281,23 +324,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 @@ -307,17 +368,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 @@ -325,21 +387,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.multiplicityKeywords[ + 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, ) @@ -349,28 +413,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" From 853a9cbbad7c0ed304486bc34ba07368b062ab78 Mon Sep 17 00:00:00 2001 From: JacksonBurns <33505528+JacksonBurns@users.noreply.github.com> Date: Mon, 8 May 2023 22:22:50 -0400 Subject: [PATCH 08/14] fix typo (delta not selta), decrease required close enough, format with black --- rmgpy/qm/mopacTest.py | 176 +++++++++++++++++++++++++++++------------- 1 file changed, 124 insertions(+), 52 deletions(-) diff --git a/rmgpy/qm/mopacTest.py b/rmgpy/qm/mopacTest.py index ec85b45d0b..5e78ed274d 100644 --- a/rmgpy/qm/mopacTest.py +++ b/rmgpy/qm/mopacTest.py @@ -39,8 +39,8 @@ from rmgpy.qm.main import QMCalculator from rmgpy.qm.mopac import Mopac, MopacMolPM3, MopacMolPM6, MopacMolPM7 -mol1 = Molecule().from_smiles('C1=CC=C2C=CC=CC2=C1') -MOPAC_CLOSE_ENOUGH_PERCENT = 0.0001 # 0.01% +mol1 = Molecule().from_smiles("C1=CC=C2C=CC=CC2=C1") +MOPAC_CLOSE_ENOUGH_PERCENT = 0.001 # 0.1% class TestMopacMolPM3(unittest.TestCase): @@ -52,13 +52,14 @@ 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) @@ -70,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, delta=169708*MOPAC_CLOSE_ENOUGH_PERCENT) - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 334.500, delta=334.500*MOPAC_CLOSE_ENOUGH_PERCENT) + 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, delta=169708*MOPAC_CLOSE_ENOUGH_PERCENT) - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 334.500, delta=334.500*MOPAC_CLOSE_ENOUGH_PERCENT) + 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): @@ -114,13 +138,14 @@ 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) @@ -132,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, delta=167704*MOPAC_CLOSE_ENOUGH_PERCENT) - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 338.099, delta=338.099*MOPAC_CLOSE_ENOUGH_PERCENT) + 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, delta=167704*MOPAC_CLOSE_ENOUGH_PERCENT) # to 0 decimal place - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 338.099, selta=338.099*MOPAC_CLOSE_ENOUGH_PERCENT) + 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): @@ -177,18 +225,19 @@ 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): @@ -196,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, 166169, delta=166169*MOPAC_CLOSE_ENOUGH_PERCENT) - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 336.333, delta=336.333*MOPAC_CLOSE_ENOUGH_PERCENT) + 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, 166169, delta=166169*MOPAC_CLOSE_ENOUGH_PERCENT) - self.assertAlmostEqual(self.qmmol1.thermo.S298.value_si, 336.333, delta=336.333*MOPAC_CLOSE_ENOUGH_PERCENT) + 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)) From dd98e07433ea017261becb3085012654de487ad0 Mon Sep 17 00:00:00 2001 From: Jackson Burns Date: Tue, 9 May 2023 08:41:40 -0400 Subject: [PATCH 09/14] set reproduce threshold to half a percent --- rmgpy/qm/mainTest.py | 222 +++++++++++++++++++++++++----------------- rmgpy/qm/mopacTest.py | 2 +- 2 files changed, 131 insertions(+), 93 deletions(-) diff --git a/rmgpy/qm/mainTest.py b/rmgpy/qm/mainTest.py index ab1474a7f2..2a2443ab72 100644 --- a/rmgpy/qm/mainTest.py +++ b/rmgpy/qm/mainTest.py @@ -40,7 +40,7 @@ from rmgpy.species import Species -MOPAC_CLOSE_ENOUGH_PERCENT = 0.0001 # 0.01% +MOPAC_CLOSE_ENOUGH_PERCENT = 0.005 # 0.5% class TestQMSettings(unittest.TestCase): @@ -52,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() @@ -91,46 +92,42 @@ 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) @@ -154,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) @@ -178,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) @@ -194,18 +191,20 @@ 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) @@ -216,78 +215,117 @@ 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, 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.") + 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 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) @@ -295,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/mopacTest.py b/rmgpy/qm/mopacTest.py index 5e78ed274d..7ada11f673 100644 --- a/rmgpy/qm/mopacTest.py +++ b/rmgpy/qm/mopacTest.py @@ -40,7 +40,7 @@ from rmgpy.qm.mopac import Mopac, MopacMolPM3, MopacMolPM6, MopacMolPM7 mol1 = Molecule().from_smiles("C1=CC=C2C=CC=CC2=C1") -MOPAC_CLOSE_ENOUGH_PERCENT = 0.001 # 0.1% +MOPAC_CLOSE_ENOUGH_PERCENT = 0.005 # 0.5% class TestMopacMolPM3(unittest.TestCase): From 795de001e51ef4fc565881ad1fbe82c524093966 Mon Sep 17 00:00:00 2001 From: Jackson Burns <33505528+JacksonBurns@users.noreply.github.com> Date: Tue, 9 May 2023 14:14:11 -0400 Subject: [PATCH 10/14] Update rmgpy/qm/mopac.py Comment from review by @rwest Co-authored-by: Richard West --- rmgpy/qm/mopac.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/qm/mopac.py b/rmgpy/qm/mopac.py index 2f7520cddb..ffdb1a5f4a 100644 --- a/rmgpy/qm/mopac.py +++ b/rmgpy/qm/mopac.py @@ -56,7 +56,7 @@ class Mopac(object): executablePath = distutils.spawn.find_executable("mopac") except: logging.debug( - "Did not find MOPAC on path, checking if it exists in a declared MOPAC_DIR..." + "Did not find mopac on path, checking if it exists in a declared MOPAC_DIR..." ) mopacEnv = os.getenv("MOPAC_DIR", default="/opt/mopac") executablePath = os.path.join(mopacEnv, "mopac") From 101f8b59489c55ade3e2682642e44e4f0779cc68 Mon Sep 17 00:00:00 2001 From: Jackson Burns Date: Tue, 9 May 2023 14:28:45 -0400 Subject: [PATCH 11/14] black formatting, rename variables for PEP8 --- rmgpy/qm/gaussian.py | 207 ++++++++++++++++++++++++------------- rmgpy/qm/gaussianTest.py | 86 ++++++++++------ rmgpy/qm/mainTest.py | 6 +- rmgpy/qm/molecule.py | 215 +++++++++++++++++++++++++-------------- rmgpy/qm/mopac.py | 57 ++++++----- rmgpy/qm/symmetry.py | 102 +++++++++++-------- 6 files changed, 428 insertions(+), 245 deletions(-) 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 2a2443ab72..dcfe492910 100644 --- a/rmgpy/qm/mainTest.py +++ b/rmgpy/qm/mainTest.py @@ -83,10 +83,10 @@ class TestQMCalculator(unittest.TestCase): Contains unit tests for the QMSettings class. """ - mopExecutablePath = Mopac.executablePath + mopexecutable_path = Mopac.executable_path - gaussExecutablePath = Gaussian.executablePath - NO_GAUSSIAN = not os.path.exists(gaussExecutablePath) + gaussexecutable_path = Gaussian.executable_path + NO_GAUSSIAN = not os.path.exists(gaussexecutable_path) def setUp(self): """ diff --git a/rmgpy/qm/molecule.py b/rmgpy/qm/molecule.py index 58abd463bb..67fd30e4d1 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 ffdb1a5f4a..025daa19cf 100644 --- a/rmgpy/qm/mopac.py +++ b/rmgpy/qm/mopac.py @@ -49,24 +49,24 @@ class Mopac(object): Classes such as :class:`MopacMol` will inherit from this class. """ - inputFileExtension = ".mop" - outputFileExtension = ".out" + input_file_extension = ".mop" + output_file_extension = ".out" try: - executablePath = distutils.spawn.find_executable("mopac") + executable_path = distutils.spawn.find_executable("mopac") except: 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") - executablePath = os.path.join(mopacEnv, "mopac") - if not os.path.exists(executablePath): - executablePath = os.path.join(mopacEnv, "(MOPAC 2009 or 2012 or 2016)") + executable_path = os.path.join(mopacEnv, "mopac") + if not os.path.exists(executable_path): + executable_path = os.path.join(mopacEnv, "(MOPAC 2009 or 2012 or 2016)") - usePolar = False # use polar keyword in MOPAC + use_polar = False # use polar keyword in MOPAC # Keywords for the multiplicity - multiplicityKeywords = { + multiplicity_keywords = { 1: "", 2: "uhf doublet", 3: "uhf triplet", @@ -80,31 +80,34 @@ class Mopac(object): #: List of phrases that indicate failure #: NONE of these must be present in a succesful job. - failureKeys = [ + 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"] + success_keys = ["DESCRIPTION OF VIBRATIONS", "MOPAC DONE"] def run(self): # ensure that the executable is present - if not os.path.exists(self.executablePath): + 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.executablePath) + "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) + 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.executablePath, tempInpFile], stdin=PIPE, stdout=PIPE, stderr=PIPE + [self.executable_path, temp_input_file], + stdin=PIPE, + stdout=PIPE, + stderr=PIPE, ) stdout, stderr = process.communicate( input=None @@ -113,11 +116,13 @@ def run(self): 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): @@ -148,13 +153,13 @@ def verify_output_file(self): ) # Initialize dictionary with "False"s - success_keys_found = dict([(key, False) for key in self.successKeys]) + 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( @@ -163,7 +168,7 @@ def verify_output_file(self): ) 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 @@ -304,7 +309,7 @@ def write_input_file(self, attempt): mopac_file.write(input_string) mopac_file.write("\n") mopac_file.write(bottom_keys) - if self.usePolar: + if self.use_polar: mopac_file.write("\n\n\n") mopac_file.write(polar_keys) @@ -387,7 +392,7 @@ def input_file_keywords(self, attempt): if attempt > self.script_attempts: attempt -= self.script_attempts - multiplicity_keys = self.multiplicityKeywords[ + multiplicity_keys = self.multiplicity_keywords[ self.geometry.molecule.multiplicity ] diff --git a/rmgpy/qm/symmetry.py b/rmgpy/qm/symmetry.py index ab4d0c702b..de8fd41e45 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(): """ 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,25 +117,23 @@ 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 @@ -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)) @@ -281,9 +289,15 @@ def calculate(self): self.point_group_found = True 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 From 3b444a92c04eaf8c570ab4deeaa870d630792333 Mon Sep 17 00:00:00 2001 From: Jackson Burns Date: Tue, 9 May 2023 14:41:41 -0400 Subject: [PATCH 12/14] module-level point group dictionary should be all caps (PEP8) --- rmgpy/qm/molecule.py | 2 +- rmgpy/qm/symmetry.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/rmgpy/qm/molecule.py b/rmgpy/qm/molecule.py index 67fd30e4d1..50277de9ef 100644 --- a/rmgpy/qm/molecule.py +++ b/rmgpy/qm/molecule.py @@ -514,7 +514,7 @@ def load_thermo_data(self): assert isinstance(thermo, rmgpy.thermo.ThermoData) self.thermo = thermo - self.point_group = symmetry.point_group_dictionary[ + 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"] diff --git a/rmgpy/qm/symmetry.py b/rmgpy/qm/symmetry.py index de8fd41e45..8fb2d6cdf2 100644 --- a/rmgpy/qm/symmetry.py +++ b/rmgpy/qm/symmetry.py @@ -65,7 +65,7 @@ def __repr__(self): ) -def make_point_group_dictionary(): +def _make_point_group_dictionary(): """ A function to make and fill the point group dictionary. @@ -138,7 +138,7 @@ def make_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): @@ -285,9 +285,9 @@ 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 " From b20c813a06ce47ef960a877fac53bed0da275fe8 Mon Sep 17 00:00:00 2001 From: Jackson Burns Date: Tue, 9 May 2023 14:43:40 -0400 Subject: [PATCH 13/14] remove dead, untested, and undocumented MOPAC_DIR lines --- rmgpy/qm/mopac.py | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/rmgpy/qm/mopac.py b/rmgpy/qm/mopac.py index 025daa19cf..b941a82362 100644 --- a/rmgpy/qm/mopac.py +++ b/rmgpy/qm/mopac.py @@ -52,16 +52,7 @@ class Mopac(object): input_file_extension = ".mop" output_file_extension = ".out" - try: - executable_path = distutils.spawn.find_executable("mopac") - except: - 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") - executable_path = os.path.join(mopacEnv, "mopac") - if not os.path.exists(executable_path): - executable_path = os.path.join(mopacEnv, "(MOPAC 2009 or 2012 or 2016)") + executable_path = distutils.spawn.find_executable("mopac") use_polar = False # use polar keyword in MOPAC From ddf97ee4bed255365f55a5eb05789a1a4d613481 Mon Sep 17 00:00:00 2001 From: Richard West Date: Wed, 10 May 2023 14:00:29 -0400 Subject: [PATCH 14/14] Remove conda-forge::libgfortran-ng>=10.0 from environment.yml for MacOS This package can't be found on MacOS, and apparently doesn't need to be specified explicitly for linux. --- environment.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/environment.yml b/environment.yml index 50ac8f3394..f16342f1c5 100644 --- a/environment.yml +++ b/environment.yml @@ -24,7 +24,6 @@ dependencies: - markupsafe - matplotlib >=1.5 - conda-forge::mopac - - conda-forge::libgfortran-ng>=10.0 - mpmath - rmg::muq2 - networkx