diff --git a/.conda/meta.yaml b/.conda/meta.yaml index 7cd02969fe7..e781b460496 100644 --- a/.conda/meta.yaml +++ b/.conda/meta.yaml @@ -14,7 +14,6 @@ requirements: - {{ compiler('c') }} # [unix] host: - cython >=0.25.2 - - lpsolve55 - numpy - openbabel >=3 - pydas >=1.0.2 @@ -39,7 +38,6 @@ requirements: - h5py - jinja2 - jupyter - - lpsolve55 - markupsafe - matplotlib >=1.5 - mopac diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 7230afa5b7e..3d79209c3b6 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -25,6 +25,7 @@ # 2023-07-17 - made it pass by default # 2023-07-21 - upload the regression results summary as artifact (for use as a comment on PRs) # 2023-07-31 - removed option to run from RMG-database with GitHub resuable workflows +# 2024-10-01 - deprecated Mambaforge with Miniforge3 for environment creation name: Continuous Integration @@ -32,8 +33,8 @@ on: schedule: # * is a special character in YAML so you have to quote this string - cron: "0 8 * * *" - # runs on all branches on both RMG-Py and forks - push: + # allow running on RMG-Py on a pushed branch, only if triggered manually + workflow_dispatch: # runs on PRs against RMG-Py (and anywhere else, but we add this for RMG-Py) pull_request: @@ -56,8 +57,9 @@ jobs: matrix: python-version: ["3.9"] os: [macos-13, macos-latest, ubuntu-latest] + test_julia: [yes, no] runs-on: ${{ matrix.os }} - name: ${{ matrix.os }} Build and Test Python ${{ matrix.python-version }} + name: ${{ matrix.os }} Build and Test Python ${{ matrix.python-version }} (Julia? ${{ matrix.test_julia }}) # skip scheduled runs from forks if: ${{ !( github.repository != 'ReactionMechanismGenerator/RMG-Py' && github.event_name == 'schedule' ) }} env: @@ -70,27 +72,28 @@ jobs: - name: Checkout RMG-Py uses: actions/checkout@v4 + # Step to create a custom condarc.yml before setting up conda + - name: Create custom condarc.yml + run: | + RUNNER_CWD=$(pwd) + echo "channels:" > $RUNNER_CWD/condarc.yml + echo " - conda-forge" >> $RUNNER_CWD/condarc.yml + echo " - rmg" >> $RUNNER_CWD/condarc.yml + echo " - cantera" >> $RUNNER_CWD/condarc.yml + echo "show_channel_urls: true" >> $RUNNER_CWD/condarc.yml + # configures the mamba environment manager and builds the environment - name: Setup Mambaforge Python ${{ matrix.python-version }} uses: conda-incubator/setup-miniconda@v3 with: environment-file: environment.yml - miniforge-variant: Mambaforge + miniforge-variant: Miniforge3 miniforge-version: latest python-version: ${{ matrix.python-version }} + condarc-file: condarc.yml activate-environment: rmg_env use-mamba: true - - name: Fix ARM Mac environment - if: matrix.os == 'macos-latest' - run: | - conda deactivate - conda env remove -n rmg_env - conda create -n rmg_env - conda activate rmg_env - conda config --env --set subdir osx-64 - conda env update -f environment.yml - # list the environment for debugging purposes - name: mamba info run: | @@ -118,31 +121,36 @@ jobs: # Setup Juliaup - name: Set Julia paths + if: matrix.test_julia == 'yes' run: | - echo "JULIAUP_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_ENV - echo "JULIAUP_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_PATH - echo "JULIA_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_ENV - echo "JULIA_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_PATH - echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_ENV - echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_PATH - echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_ENV - echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_PATH + # echo "JULIAUP_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_ENV + # echo "JULIAUP_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_PATH + # echo "JULIA_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_ENV + # echo "JULIA_DEPOT_PATH=$CONDA/envs/rmg_env/.julia" >> $GITHUB_PATH + # echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_ENV + # echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_PATH + # echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_ENV + # echo "JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba" >> $GITHUB_PATH echo "JULIA_CONDAPKG_BACKEND=Current" >> $GITHUB_ENV - echo "JULIA_CONDAPKG_BACKEND=Current" >> $GITHUB_PATH + # echo "JULIA_CONDAPKG_BACKEND=Current" >> $GITHUB_PATH - name: Setup Juliaup + if: matrix.test_julia == 'yes' uses: julia-actions/install-juliaup@v2 with: - channel: '1.10' + channel: '1.9' - name: Check Julia version + if: matrix.test_julia == 'yes' run: julia --version # RMS installation and linking to Julia - name: Install and link Julia dependencies + if: matrix.test_julia == 'yes' timeout-minutes: 120 # this usually takes 20-45 minutes (or hangs for 6+ hours). # JULIA_CONDAPKG_EXE points to the existing conda/mamba to avoid JuliaCall from installing their own. See https://juliapy.github.io/PythonCall.jl/stable/pythoncall/#If-you-already-have-a-Conda-environment. run: | + mamba install conda-forge::pyjuliacall julia -e 'using Pkg; Pkg.add(Pkg.PackageSpec(name="ReactionMechanismSimulator", url="https://github.com/hwpang/ReactionMechanismSimulator.jl.git", rev="fix_installation")); using ReactionMechanismSimulator' - name: Install Q2DTor @@ -158,7 +166,7 @@ jobs: id: regression-execution timeout-minutes: 60 run: | - for regr_test in aromatics liquid_oxidation nitrogen oxidation sulfur superminimal RMS_constantVIdealGasReactor_superminimal RMS_CSTR_liquid_oxidation RMS_liquidSurface_ch4o2cat fragment RMS_constantVIdealGasReactor_fragment; + for regr_test in aromatics liquid_oxidation nitrogen oxidation sulfur superminimal RMS_constantVIdealGasReactor_superminimal RMS_CSTR_liquid_oxidation RMS_liquidSurface_ch4o2cat fragment RMS_constantVIdealGasReactor_fragment minimal_surface; do if python rmg.py test/regression/"$regr_test"/input.py; then echo "$regr_test" "Executed Successfully" @@ -209,16 +217,16 @@ jobs: - name : Find ID of Reference Results env: GH_TOKEN: ${{ secrets.GITHUB_TOKEN }} + # this will search for the last successful execution of CI on main run: | - run_id=$(gh run list -R ReactionMechanismGenerator/RMG-Py --workflow="Continuous Integration" --branch main --limit 1 --json databaseId --jq '.[0].databaseId') + run_id=$(gh run list -R ReactionMechanismGenerator/RMG-Py --workflow="Continuous Integration" --branch main --limit 15 --json databaseId,conclusion --jq 'map(select(.conclusion == "success")) | .[0].databaseId') echo "CI_RUN_ID=$run_id" >> $GITHUB_ENV - name: Retrieve Stable Regression Results if: ${{ env.REFERENCE_JOB == 'false' }} uses: actions/download-artifact@v4 with: - # this will search for the last successful execution of CI on main and download - # the stable regression results + # download stable regression results run-id: ${{ env.CI_RUN_ID }} repository: ReactionMechanismGenerator/RMG-Py github-token: ${{ secrets.GITHUB_TOKEN }} @@ -234,7 +242,7 @@ jobs: run: | exec 2> >(tee -a regression.stderr >&2) 1> >(tee -a regression.stdout) mkdir -p "test/regression-diff" - for regr_test in aromatics liquid_oxidation nitrogen oxidation sulfur superminimal RMS_constantVIdealGasReactor_superminimal RMS_CSTR_liquid_oxidation fragment RMS_constantVIdealGasReactor_fragment; + for regr_test in aromatics liquid_oxidation nitrogen oxidation sulfur superminimal RMS_constantVIdealGasReactor_superminimal RMS_CSTR_liquid_oxidation fragment RMS_constantVIdealGasReactor_fragment minimal_surface; do echo "" echo "### Regression test $regr_test:" @@ -369,7 +377,7 @@ jobs: password: ${{ secrets.DOCKERHUB_TOKEN }} - name: Build and Push - uses: docker/build-push-action@v4 + uses: docker/build-push-action@v6 with: push: true tags: reactionmechanismgenerator/rmg:latest diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 7adffb55fa9..b500afd43d0 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -30,7 +30,7 @@ jobs: uses: conda-incubator/setup-miniconda@v2 with: environment-file: environment.yml - miniforge-variant: Mambaforge + miniforge-variant: Miniforge3 miniforge-version: latest python-version: 3.9 activate-environment: rmg_env @@ -60,13 +60,6 @@ jobs: make clean make - - name: Install and link Julia dependencies - timeout-minutes: 120 # this usually takes 20-45 minutes (or hangs for 6+ hours). - run: | - which julia - export JULIA_CONDAPKG_EXE=$CONDA/condabin/mamba - julia -e 'ENV["JULIA_CONDAPKG_BACKEND"] = "Current"; using Pkg; Pkg.add(PackageSpec(name="ReactionMechanismSimulator", url="https://github.com/hwpang/ReactionMechanismSimulator.jl.git", rev="fix_installation")); using ReactionMechanismSimulator' - - name: Checkout gh-pages Branch uses: actions/checkout@v2 with: diff --git a/.gitignore b/.gitignore index 946eb94d4ac..38d73e2bc18 100644 --- a/.gitignore +++ b/.gitignore @@ -103,6 +103,7 @@ test/rmgpy/test_data/temp_dir_for_testing/cantera/chem001.yaml rmgpy/test_data/copied_kinetic_lib/ testing/qm/* test_log.txt +rmgpy/tools/data/flux/flux/1/*.dot # example directory - save the inputs but not the outputs # cantera input files diff --git a/Dockerfile b/Dockerfile index d237a659c7d..974235c21bb 100644 --- a/Dockerfile +++ b/Dockerfile @@ -23,29 +23,30 @@ RUN apt-get update && \ apt-get clean -y # Install conda -RUN wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh && \ - bash Miniconda3-latest-Linux-x86_64.sh -b -p /miniconda && \ - rm Miniconda3-latest-Linux-x86_64.sh -ENV PATH="$PATH:/miniconda/bin" - -# Set solver backend to mamba for speed -RUN conda install -n base conda-libmamba-solver && \ - conda config --set solver libmamba +RUN wget "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-x86_64.sh" && \ + bash Miniforge3-Linux-x86_64.sh -b -p /miniforge && \ + rm Miniforge3-Linux-x86_64.sh +ENV PATH="$PATH:/miniforge/bin" # Set Bash as the default shell for following commands SHELL ["/bin/bash", "-c"] +# Add build arguments for RMG-Py, RMG-database, and RMS branches. +ARG RMG_Py_Branch=main +ARG RMG_Database_Branch=main +ARG RMS_Branch=main +ENV rmsbranch=${RMS_Branch} + # cd WORKDIR /rmg # Clone the RMG base and database repositories -RUN git clone --single-branch --branch main --depth 1 https://github.com/ReactionMechanismGenerator/RMG-Py.git && \ - git clone --single-branch --branch main --depth 1 https://github.com/ReactionMechanismGenerator/RMG-database.git +RUN git clone --single-branch --branch ${RMG_Py_Branch} --depth 1 https://github.com/ReactionMechanismGenerator/RMG-Py.git && \ + git clone --single-branch --branch ${RMG_Database_Branch} --depth 1 https://github.com/ReactionMechanismGenerator/RMG-database.git WORKDIR /rmg/RMG-Py # build the conda environment -RUN conda env create --file environment.yml && \ - conda clean --all --yes +RUN conda env create --file environment.yml # This runs all subsequent commands inside the rmg_env conda environment # @@ -54,6 +55,10 @@ RUN conda env create --file environment.yml && \ # in a Dockerfile build script) SHELL ["conda", "run", "--no-capture-output", "-n", "rmg_env", "/bin/bash", "-c"] +RUN conda install -c conda-forge julia=1.9.1 pyjulia>=0.6 && \ + conda install -c rmg pyrms diffeqpy && \ + conda clean --all --yes + # Set environment variables as directed in the RMG installation instructions ENV RUNNER_CWD=/rmg ENV PYTHONPATH="$RUNNER_CWD/RMG-Py:$PYTHONPATH" diff --git a/arkane/data/molpro/HO2_RS2C.out b/arkane/data/molpro/HO2_RS2C.out new file mode 100644 index 00000000000..9142f71520d --- /dev/null +++ b/arkane/data/molpro/HO2_RS2C.out @@ -0,0 +1,531 @@ + + Working directory : /gtmp/molpro.T66wyGqmSD/ + Global scratch directory : /gtmp/molpro.T66wyGqmSD/ + Wavefunction directory : /home/alon/wfu/ + Main file repository : /gtmp/molpro.T66wyGqmSD/ + + id : dana + + Nodes nprocs + n136 16 + GA implementation: MPI file + GA implementation (serial work in mppx): MPI file + + Using customized tuning parameters: mindgm=1; mindgv=20; mindgc=4; mindgr=1; noblas=0; minvec=7 + default implementation of scratch files=sf + + + Variables initialized (1025), CPU time= 0.01 sec + ***,HO2 + memory,Total=15625,m; + + geometry={angstrom; + O -0.16513900 0.44815400 0.00000000 + O 1.00112100 -0.18975400 -0.00000000 + H -0.83598100 -0.25840000 0.00000000} + + gprint,orbitals; + + basis=aug-cc-pvtz + + + + int; + + {hf; + maxit,999; + wf,spin=1,charge=0; + } + + + {mp2; + wf,spin=1,charge=0; + } + + {casscf; + maxit,999; + wf,spin=1,charge=0; + occ,8,2; + closed,5,1; + state,1; + } + + {rs2c; + maxit,999; + wf,spin=1,charge=0; + } + + + + + Commands initialized (840), CPU time= 0.00 sec, 684 directives. + Default parameters read. Elapsed time= 0.32 sec + + Checking input... + Passed +1 + + + *** PROGRAM SYSTEM MOLPRO *** + Copyright, TTI GmbH Stuttgart, 2015 + Version 2022.3 linked Wed Nov 30 13:40:34 2022 + + + ********************************************************************************************************************************** + LABEL * HO2 + 64 bit mpp version DATE: 08-Sep-24 TIME: 09:43:40 + ********************************************************************************************************************************** + + SHA1: e31ec9a5ea85254ab76f59d122cbdd51c71cf98b + ********************************************************************************************************************************** + + Memory per process: 732 MW + Total memory per node: 15625 MW + Total GA space: 3906 MW + + GA preallocation enabled + GA check enabled + + Variable memory set to 732.4 MW + + + Geometry recognized as XYZ + + + SETTING BASIS = AUG-CC-PVTZ + + + Using spherical harmonics + + Library entry O S aug-cc-pVTZ selected for orbital group 1 + Library entry O P aug-cc-pVTZ selected for orbital group 1 + Library entry O D aug-cc-pVTZ selected for orbital group 1 + Library entry O F aug-cc-pVTZ selected for orbital group 1 + Library entry H S aug-cc-pVTZ selected for orbital group 2 + Library entry H P aug-cc-pVTZ selected for orbital group 2 + Library entry H D aug-cc-pVTZ selected for orbital group 2 + + + PROGRAM * SEWARD (Integral evaluation for generally contracted gaussian basis sets) Author: Roland Lindh, 1990 + + Geometry written to block 1 of record 700 + + + Point group Cs + + + + ATOMIC COORDINATES + + NR ATOM CHARGE X Y Z + + 1 O 8.00 -0.312067482 0.846888322 0.000000000 + 2 O 8.00 1.891844508 -0.358583091 0.000000000 + 3 H 1.00 -1.579775135 -0.488305231 0.000000000 + + Bond lengths in Bohr (Angstrom) + + 1-2 2.512048842 1-3 1.841147608 + ( 1.329319000) ( 0.974293356) + + Bond angles + + 2-1-3 104.83752478 + + NUCLEAR CHARGE: 17 + NUMBER OF PRIMITIVE AOS: 161 + NUMBER OF SYMMETRY AOS: 141 + NUMBER OF CONTRACTIONS: 115 ( 76A' + 39A" ) + NUMBER OF INNER CORE ORBITALS: 0 ( 0A' + 0A" ) + NUMBER OF OUTER CORE ORBITALS: 2 ( 2A' + 0A" ) + NUMBER OF VALENCE ORBITALS: 9 ( 7A' + 2A" ) + + + NUCLEAR REPULSION ENERGY 32.12512051 + + + Allocated 245 MW GA space on the current processor, total: 3907 MW. Time: 0.00 sec + + Eigenvalues of metric + + 1 0.240E-03 0.415E-03 0.122E-02 0.267E-02 0.374E-02 0.493E-02 0.520E-02 0.693E-02 + 2 0.587E-02 0.114E-01 0.380E-01 0.419E-01 0.523E-01 0.813E-01 0.118E+00 0.128E+00 + + + Contracted 2-electron integrals neglected if value below 1.0D-11 + AO integral compression algorithm 1 Integral accuracy 1.0D-11 + + 65.798 MB (compressed) written to integral file ( 55.4%) + + Node minimum: 3.408 MB, node maximum: 5.243 MB + + + NUMBER OF SORTED TWO-ELECTRON INTEGRALS: 706608. BUFFER LENGTH: 32768 + NUMBER OF SEGMENTS: 1 SEGMENT LENGTH: 706608 RECORD LENGTH: 524288 + + Memory used in sort: 1.26 MW + + SORT1 READ 14525730. AND WROTE 675498. INTEGRALS IN 2 RECORDS. CPU TIME: 0.14 SEC, REAL TIME: 0.15 SEC + SORT2 READ 10930363. AND WROTE 11433631. INTEGRALS IN 272 RECORDS. CPU TIME: 0.01 SEC, REAL TIME: 0.01 SEC + + Node minimum: 706378. Node maximum: 718170. integrals + + OPERATOR DM FOR CENTER 0 COORDINATES: 0.000000 0.000000 0.000000 + + + ********************************************************************************************************************************** + DATASETS * FILE NREC LENGTH (MB) RECORD NAMES + 1 18 28.86 500 610 700 900 950 970 1000 129 960 1100 + VAR BASINP GEOM SYMINP ZMAT AOBASIS BASIS P2S ABASIS S + 1400 1410 1200 1210 1080 1600 1650 1700 + T V H0 H01 AOSYM SMH MOLCAS OPER + + PROGRAMS * TOTAL INT + CPU TIMES * 0.60 0.47 + REAL TIME * 1.03 SEC + DISK USED * 29.18 MB (local), 644.53 MB (total) + GA USED * 0.00 MB (max) 0.00 MB (current) + ********************************************************************************************************************************** + + + Program * Restricted Hartree-Fock + + Orbital guess generated from atomic densities. Full valence occupancy: 9 2 + + Initial alpha occupancy: 7 2 + Initial beta occupancy: 7 1 + + NELEC= 17 SYM=2 MS2= 1 THRE=1.0D-08 THRD=3.2D-06 THRG=3.2D-06 HFMA2=F DIIS_START=2 DIIS_MAX=10 DIIS_INCORE=F + + Level shifts: 0.00 (CLOSED) 0.00 (OPEN) 0.30 (GAP_MIN) + + ITER ETOT DE GRAD DDIFF DIIS NEXP TIME(IT) TIME(TOT) DIAG + 1 -150.17011817 -150.17011817 0.00D+00 0.88D-01 0 0 0.00 0.00 start + 2 -150.20465019 -0.03453202 0.90D-02 0.92D-02 1 0 0.01 0.01 diag2 + 3 -150.22033330 -0.01568311 0.64D-02 0.47D-02 2 0 0.00 0.01 diag2 + 4 -150.22490807 -0.00457477 0.13D-02 0.14D-02 3 0 0.01 0.02 diag2 + 5 -150.22856663 -0.00365856 0.12D-02 0.11D-02 4 0 0.01 0.03 diag2 + 6 -150.23199458 -0.00342794 0.80D-03 0.16D-02 5 0 0.00 0.03 diag2 + 7 -150.23257275 -0.00057817 0.35D-03 0.88D-03 6 0 0.00 0.03 diag2 + 8 -150.23258496 -0.00001221 0.98D-04 0.16D-03 7 0 0.01 0.04 fixocc + 9 -150.23259040 -0.00000544 0.69D-04 0.13D-03 8 0 0.00 0.04 diag2 + 10 -150.23259103 -0.00000063 0.28D-04 0.43D-04 9 0 0.01 0.05 diag2/orth + 11 -150.23259126 -0.00000022 0.13D-04 0.19D-04 9 0 0.00 0.05 diag2 + 12 -150.23259126 -0.00000000 0.95D-06 0.19D-05 9 0 0.01 0.06 diag2 + 13 -150.23259126 -0.00000000 0.51D-06 0.64D-06 0 0 0.00 0.06 diag + + Final alpha occupancy: 7 2 + Final beta occupancy: 7 1 + + !RHF STATE 1.2 Energy -150.232591258293 + RHF One-electron energy -270.928629502666 + RHF Two-electron energy 88.570917737861 + RHF Kinetic energy 150.046375775003 + RHF Nuclear energy 32.125120506513 + RHF Virial quotient -1.001241052857 + + !RHF STATE 1.2 Dipole moment -0.62605931 -0.43865566 0.00000000 + Dipole moment /Debye -1.59128405 -1.11495147 0.00000000 + + ELECTRON ORBITALS + ================= + + Orbital Occupation Energy Cen Mu Typ Coefficients + 1.1 2.00000 -20.67277 2 1 s 0.99737 + 2.1 2.00000 -20.66835 1 1 s 0.99740 + 3.1 2.00000 -1.54217 1 2 s 0.68214 2 2 s 0.56715 + 4.1 2.00000 -1.15777 1 2 s -0.57468 1 1 px 0.25089 2 2 s 0.74537 3 1 s -0.36249 + 5.1 2.00000 -0.76227 1 1 px 0.46397 1 1 py 0.46612 2 1 py 0.47384 3 1 s -0.53346 + 6.1 2.00000 -0.66270 1 1 px -0.39599 1 1 py 0.48511 2 1 px 0.64562 + 7.1 2.00000 -0.50672 1 1 py -0.47553 2 1 px 0.32427 2 1 py 0.79533 3 1 s 0.36930 + 1.2 2.00000 -0.62448 1 1 pz 0.88505 + 2.2 1.00000 -0.59233 1 1 pz -0.36862 2 1 pz 0.99301 + + + HOMO 7.1 -0.506721 = -13.7886eV + LUMO 8.1 0.030590 = 0.8324eV + LUMO-HOMO 0.537311 = 14.6210eV + + Orbitals saved in record 2100.2 + + + ********************************************************************************************************************************** + DATASETS * FILE NREC LENGTH (MB) RECORD NAMES + 1 18 28.87 500 610 700 900 950 970 1000 129 960 1100 + VAR BASINP GEOM SYMINP ZMAT AOBASIS BASIS P2S ABASIS S + 1400 1410 1200 1210 1080 1600 1650 1700 + T V H0 H01 AOSYM SMH MOLCAS OPER + + 2 4 0.62 700 1000 520 2100 + GEOM BASIS MCVARS RHF + + PROGRAMS * TOTAL HF-SCF INT + CPU TIMES * 0.74 0.14 0.47 + REAL TIME * 5.11 SEC + DISK USED * 30.56 MB (local), 666.61 MB (total) + GA USED * 0.00 MB (max) 0.00 MB (current) + ********************************************************************************************************************************** + + + PROGRAM * RMP2 (Restricted open-shell) Authors: H.-J. Werner, 2001 + + + Convergence thresholds: THRVAR = 1.00D-08 THRDEN = 1.00D-06 + + Number of core orbitals: 2 ( 2 0 ) + Number of closed-shell orbitals: 6 ( 5 1 ) + Number of active orbitals: 1 ( 0 1 ) + Number of external orbitals: 106 ( 69 37 ) + + Reading RHF information from record 2100.2 + + TOTAL A-A B-B A-B + Pair energies: -0.45576889 -0.06238867 -0.05285220 -0.34052802 + Singles energies: -0.00441624 -0.00220300 -0.00221324 + + !Reference STATE 1.2 energy -150.232591258293 + !RMP2 STATE 1.2 correl. energy -0.460185122779 + !RMP2 STATE 1.2 total energy -150.692776381072 + + !SCS-RMP2 STATE 1.2 correl. energy -0.447047243835 + !SCS-RMP2 STATE 1.2 total energy -150.679638502128 + + + ********************************************************************************************************************************** + DATASETS * FILE NREC LENGTH (MB) RECORD NAMES + 1 18 28.87 500 610 700 900 950 970 1000 129 960 1100 + VAR BASINP GEOM SYMINP ZMAT AOBASIS BASIS P2S ABASIS S + 1400 1410 1200 1210 1080 1600 1650 1700 + T V H0 H01 AOSYM SMH MOLCAS OPER + + 2 4 0.62 700 1000 520 2100 + GEOM BASIS MCVARS RHF + + PROGRAMS * TOTAL MP2 HF-SCF INT + CPU TIMES * 0.81 0.06 0.14 0.47 + REAL TIME * 5.19 SEC + DISK USED * 30.56 MB (local), 666.61 MB (total) + SF USED * 8.12 MB + GA USED * 0.00 MB (max) 0.00 MB (current) + ********************************************************************************************************************************** + + PROGRAM * MULTI (Direct Multiconfiguration SCF) Authors: P.J. Knowles, H.-J. Werner (1984) + D.A. Kreplin, P.J. Knowles, H.-J. Werner (2019) + + Number of closed-shell orbitals: 6 ( 5 1) + Number of active orbitals: 4 ( 3 1) + Number of external orbitals: 105 ( 68 37) + + State symmetry 1 + + Number of active electrons: 5 Spin symmetry=Doublet Space symmetry=2 + Number of states: 1 + Number of CSFs: 9 (12 determinants, 24 intermediate states) + + Molecular orbitals read from record 2100.2 Type=RHF/CANONICAL + + Wavefunction dump at record 2140.2 + + Convergence thresholds 0.10E-04 (gradient) 0.10E-05 (energy) + + Number of orbital rotations: 634 ( 16 closed/active, 377 closed/virtual, 0 active/active, 241 active/virtual ) + Total number of variables: 646 + + Second-order MCSCF: L-BFGS accelerated + + ITER MIC NCI NEG ENERGY(VAR) ENERGY(PROJ) ENERGY CHANGE GRAD(0) GRAD(ORB) GRAD(CI) STEP TIME + 1 55 20 28 -150.23260182 -150.25752980 -0.02492798 0.00026496 0.00000027 0.00000000 0.94E+00 2.24 + 2 26 13 0 -150.26507911 -150.26990559 -0.00482648 0.06074391 0.00034661 0.00000000 0.23E+00 2.53 + 3 31 12 0 -150.27012384 -150.27015191 -0.00002807 0.00393155 0.00000450 0.00000000 0.19E-01 2.83 + 4 23 9 0 -150.27015263 -150.27015272 -0.00000009 0.00025798 0.00000025 0.00000000 0.98E-03 3.07 + + CONVERGENCE REACHED! Final gradient: 0.00000052 ( 0.52E-06) + Final energy: -150.27015272 + + First order charge density matrix for state 1.2 saved on record 2140.2 (density set 1) + + Results for state 1.2 + ===================== + !MCSCF STATE 1.2 Energy -150.270152720089 + Nuclear energy 32.12512051 + Kinetic energy 150.20743631 + One electron energy -270.97568640 + Two electron energy 88.58041317 + Virial ratio 2.00041753 + + !MCSCF STATE 1.2 Dipole moment -0.62063164 -0.43831921 0.00000000 + Dipole moment /Debye -1.57748829 -1.11409631 0.00000000 + + No non-zero expectation values + + NATURAL ORBITALS + ================ + + Orbital Occupation Energy Cen Mu Typ Coefficients + 1.1 2.00000 -20.65876 1 1 s 0.91254 2 1 s 0.40637 + 2.1 2.00000 -20.64159 1 1 s -0.40589 2 1 s 0.91190 + 3.1 2.00000 -1.40648 1 2 s 0.81184 2 2 s 0.37367 3 1 s 0.37720 + 4.1 2.00000 -0.75836 1 1 px 0.31043 1 1 py 0.64168 2 1 py 0.35375 3 1 s -0.54763 + 5.1 2.00000 -0.52211 1 1 py 0.34480 2 1 px -0.43218 2 1 py -0.78756 3 1 s -0.38960 + 6.1 1.99791 -1.14889 1 2 s -0.42578 2 2 s 0.89357 + 7.1 1.95773 -0.78780 1 1 px -0.54740 1 1 py 0.37155 2 1 px 0.55821 2 1 py -0.33506 + 8.1 0.04436 0.48087 1 2 s -0.40162 1 4 s -0.35196 1 1 px -0.82535 1 1 py 0.49487 + 2 2 s 0.35945 2 4 s 0.38540 2 5 s 0.26437 2 1 px -0.86572 + 2 1 py 0.48426 + 1.2 2.00000 -0.62077 1 1 pz 0.88536 + 2.2 1.00000 -0.21314 1 1 pz -0.36073 2 1 pz 0.99357 + + Natural orbital dump at molpro section 2140.2 (Orbital set 2) + + + ********************************************************************************************************************************** + DATASETS * FILE NREC LENGTH (MB) RECORD NAMES + 1 19 33.34 500 610 700 900 950 970 1000 129 960 1100 + VAR BASINP GEOM SYMINP ZMAT AOBASIS BASIS P2S ABASIS S + 1400 1410 1200 1210 1080 1600 1650 1700 1380 + T V H0 H01 AOSYM SMH MOLCAS OPER JKOP + + 2 5 0.97 700 1000 520 2100 2140 + GEOM BASIS MCVARS RHF MCSCF + + PROGRAMS * TOTAL CASSCF MP2 HF-SCF INT + CPU TIMES * 2.08 1.27 0.06 0.14 0.47 + REAL TIME * 9.10 SEC + DISK USED * 34.37 MB (local), 727.57 MB (total) + SF USED * 12.99 MB + GA USED * 0.00 MB (max) 0.00 MB (current) + ********************************************************************************************************************************** + + + PROGRAM * RS2C (Multireference RS Perturbation Theory) Authors: H.-J. Werner (1993), P. Celani (1998) + + + Convergence thresholds: THRVAR = 1.00D-06 THRDEN = 1.00D-06 + + Number of optimized states: 1 Roots: 1 + Number of reference states: 1 Roots: 1 + + Using rs2c + + Reference symmetry: 2 Doublet + Number of electrons: 17 + Maximum number of shells: 4 + Maximum number of spin couplings: 14 + + Reference space: 6 conf 9 CSFs + N elec internal: 16 conf 20 CSFs + N-1 el internal: 19 conf 35 CSFs + N-2 el internal: 13 conf 38 CSFs + + Number of electrons in valence space: 13 + Maximum number of open shell orbitals in reference space: 3 + Maximum number of open shell orbitals in internal spaces: 7 + + + Number of core orbitals: 2 ( 2 0 ) + Number of closed-shell orbitals: 4 ( 3 1 ) + Number of active orbitals: 4 ( 3 1 ) + Number of external orbitals: 105 ( 68 37 ) + + Molecular orbitals read from record 2140.2 Type=MCSCF/NATURAL (state 1.2) + + Integral transformation finished. Total CPU: 0.02 sec, npass= 1 Memory used: 0.36 MW + + Number of p-space configurations: 1 + + Reference wavefunction optimized for reference space (refopt=1) + + State Reference Energy + 1 -150.27015272 + + Number of blocks in overlap matrix: 4 Smallest eigenvalue: 0.28D-03 + Number of N-2 electron functions: 63 + Number of N-1 electron functions: 35 + + Number of internal configurations: 9 + Number of singly external configurations: 1822 + Number of doubly external configurations: 177615 + Total number of contracted configurations: 179446 + Total number of uncontracted configurations: 107372 + + Weight factors for SA-density in H0: 1.000000 + + FIMAX= 0.28D+00 FXMAX= 0.11D-08 DIAG= F F NOREF=1 NOINT=0 IHPPD=2 + + Nuclear energy: 32.12512051 + Core energy: -130.99803678 + Zeroth-order valence energy: -10.64492796 + Zeroth-order total energy: -109.51784424 + First-order energy: -40.75230848 + + + Using contracted singles + + Number of contracted N-1 electron functions: 180 + Number of contracted internal configurations: 8 + + + Diagonal Coupling coefficients finished. Storage: 175576 words, CPU-Time: 0.00 seconds. + Energy denominators for pairs finished in 0 passes. Storage: 189560 words, CPU-time: 0.00 seconds. + + ITER. STATE ROOT SQ.NORM CORR.ENERGY TOTAL ENERGY ENERGY CHANGE DEN1 VAR(S) VAR(P) TIME + 1 1 1 1.03247754 0.00000000 -150.27015272 0.00000000 -0.41711899 0.32D-01 0.63D-01 0.04 + 2 1 1 1.10304126 -0.42756068 -150.69771340 -0.42756068 -0.00025015 0.45D-03 0.16D-03 0.05 + 3 1 1 1.10435585 -0.42946425 -150.69961697 -0.00190357 -0.00181114 0.19D-03 0.17D-05 0.06 + 4 1 1 1.10498583 -0.42995826 -150.70011098 -0.00049401 0.00003012 0.64D-05 0.57D-06 0.07 + 5 1 1 1.10513054 -0.42999346 -150.70014618 -0.00003519 -0.00008892 0.17D-05 0.18D-07 0.08 + 6 1 1 1.10517684 -0.42999665 -150.70014937 -0.00000319 0.00000193 0.38D-06 0.52D-08 0.09 + 7 1 1 1.10519789 -0.42999869 -150.70015141 -0.00000204 -0.00000860 0.15D-06 0.10D-08 0.10 + 8 1 1 1.10520600 -0.42999909 -150.70015181 -0.00000040 -0.00000018 0.62D-07 0.38D-09 0.10 + + Energy contributions for state 1.2: + =================================== + + Energy contr. SQ.Norm of FOWF + Space I0 0.00000000 0.00000000 + Space I1 -0.01058738 0.00712445 + Space I2 -0.00124151 0.00064953 + Space S0 -0.02885894 0.00870840 + Space S1 -0.05688864 0.01910554 + Space S2 -0.01800045 0.00649351 + Space P0 -0.04231239 0.00657076 + Space P1 -0.15768459 0.03130729 + Space P2 -0.11442520 0.02524654 + + + RESULTS FOR STATE 1.2 + ===================== + + + Reference energy -150.27015272 + Nuclear energy 32.12512051 + Correlation energy -0.42999909 + !RSPT2 STATE 1.2 Energy -150.700151809954 + !RS2C STATE 1.2 Dipole moment -0.65532698 -0.41398806 0.00000000 + Dipole moment /Debye -1.66567504 -1.05225269 0.00000000 + + + ********************************************************************************************************************************** + DATASETS * FILE NREC LENGTH (MB) RECORD NAMES + 1 19 33.34 500 610 700 900 950 970 1000 129 960 1100 + VAR BASINP GEOM SYMINP ZMAT AOBASIS BASIS P2S ABASIS S + 1400 1410 1200 1210 1080 1600 1650 1700 1380 + T V H0 H01 AOSYM SMH MOLCAS OPER JKOP + + 2 5 0.97 700 1000 520 2100 2140 + GEOM BASIS MCVARS RHF MCSCF + + PROGRAMS * TOTAL RS2C CASSCF MP2 HF-SCF INT + CPU TIMES * 2.22 0.14 1.27 0.06 0.14 0.47 + REAL TIME * 9.81 SEC + DISK USED * 34.37 MB (local), 727.57 MB (total) + SF USED * 18.18 MB + GA USED * 0.00 MB (max) 0.00 MB (current) + ********************************************************************************************************************************** + + RS2C/aug-cc-pVTZ energy= -150.700151809954 + + RS2C CASSCF MP2 HF-SCF + -150.70015181 -150.27015272 -150.69277638 -150.23259126 + ********************************************************************************************************************************** + Molpro calculation terminated diff --git a/arkane/data/orca/N_MRCI.log b/arkane/data/orca/N_MRCI.log new file mode 100644 index 00000000000..80430936782 --- /dev/null +++ b/arkane/data/orca/N_MRCI.log @@ -0,0 +1,1662 @@ + + ***************** + * O R C A * + ***************** + + #, + ### + #### + ##### + ###### + ########, + ,,################,,,,, + ,,#################################,, + ,,##########################################,, + ,#########################################, ''#####, + ,#############################################,, '####, + ,##################################################,,,,####, + ,###########'''' ''''############################### + ,#####'' ,,,,##########,,,, '''####''' '#### + ,##' ,,,,###########################,,, '## + ' ,,###'''' '''############,,, + ,,##'' '''############,,,, ,,,,,,###'' + ,#'' '''#######################''' + ' ''''####'''' + ,#######, #######, ,#######, ## + ,#' '#, ## ## ,#' '#, #''# ###### ,####, + ## ## ## ,#' ## #' '# # #' '# + ## ## ####### ## ,######, #####, # # + '#, ,#' ## ## '#, ,#' ,# #, ## #, ,# + '#######' ## ## '#######' #' '# #####' # '####' + + + + ####################################################### + # -***- # + # Department of theory and spectroscopy # + # Directorship and core code : Frank Neese # + # Max Planck Institute fuer Kohlenforschung # + # Kaiser Wilhelm Platz 1 # + # D-45470 Muelheim/Ruhr # + # Germany # + # # + # All rights reserved # + # -***- # + ####################################################### + + + Program Version 5.0.4 - RELEASE - + + + With contributions from (in alphabetic order): + Daniel Aravena : Magnetic Suceptibility + Michael Atanasov : Ab Initio Ligand Field Theory (pilot matlab implementation) + Alexander A. Auer : GIAO ZORA, VPT2 properties, NMR spectrum + Ute Becker : Parallelization + Giovanni Bistoni : ED, misc. LED, open-shell LED, HFLD + Martin Brehm : Molecular dynamics + Dmytro Bykov : SCF Hessian + Vijay G. Chilkuri : MRCI spin determinant printing, contributions to CSF-ICE + Dipayan Datta : RHF DLPNO-CCSD density + Achintya Kumar Dutta : EOM-CC, STEOM-CC + Dmitry Ganyushin : Spin-Orbit,Spin-Spin,Magnetic field MRCI + Miquel Garcia : C-PCM and meta-GGA Hessian, CC/C-PCM, Gaussian charge scheme + Yang Guo : DLPNO-NEVPT2, F12-NEVPT2, CIM, IAO-localization + Andreas Hansen : Spin unrestricted coupled pair/coupled cluster methods + Benjamin Helmich-Paris : MC-RPA, TRAH-SCF, COSX integrals + Lee Huntington : MR-EOM, pCC + Robert Izsak : Overlap fitted RIJCOSX, COSX-SCS-MP3, EOM + Marcus Kettner : VPT2 + Christian Kollmar : KDIIS, OOCD, Brueckner-CCSD(T), CCSD density, CASPT2, CASPT2-K + Simone Kossmann : Meta GGA functionals, TD-DFT gradient, OOMP2, MP2 Hessian + Martin Krupicka : Initial AUTO-CI + Lucas Lang : DCDCAS + Marvin Lechner : AUTO-CI (C++ implementation), FIC-MRCC + Dagmar Lenk : GEPOL surface, SMD + Dimitrios Liakos : Extrapolation schemes; Compound Job, initial MDCI parallelization + Dimitrios Manganas : Further ROCIS development; embedding schemes + Dimitrios Pantazis : SARC Basis sets + Anastasios Papadopoulos: AUTO-CI, single reference methods and gradients + Taras Petrenko : DFT Hessian,TD-DFT gradient, ASA, ECA, R-Raman, ABS, FL, XAS/XES, NRVS + Peter Pinski : DLPNO-MP2, DLPNO-MP2 Gradient + Christoph Reimann : Effective Core Potentials + Marius Retegan : Local ZFS, SOC + Christoph Riplinger : Optimizer, TS searches, QM/MM, DLPNO-CCSD(T), (RO)-DLPNO pert. Triples + Tobias Risthaus : Range-separated hybrids, TD-DFT gradient, RPA, STAB + Michael Roemelt : Original ROCIS implementation + Masaaki Saitow : Open-shell DLPNO-CCSD energy and density + Barbara Sandhoefer : DKH picture change effects + Avijit Sen : IP-ROCIS + Kantharuban Sivalingam : CASSCF convergence, NEVPT2, FIC-MRCI + Bernardo de Souza : ESD, SOC TD-DFT + Georgi Stoychev : AutoAux, RI-MP2 NMR, DLPNO-MP2 response + Willem Van den Heuvel : Paramagnetic NMR + Boris Wezisla : Elementary symmetry handling + Frank Wennmohs : Technical directorship + + + We gratefully acknowledge several colleagues who have allowed us to + interface, adapt or use parts of their codes: + Stefan Grimme, W. Hujo, H. Kruse, P. Pracht, : VdW corrections, initial TS optimization, + C. Bannwarth, S. Ehlert DFT functionals, gCP, sTDA/sTD-DF + Ed Valeev, F. Pavosevic, A. Kumar : LibInt (2-el integral package), F12 methods + Garnet Chan, S. Sharma, J. Yang, R. Olivares : DMRG + Ulf Ekstrom : XCFun DFT Library + Mihaly Kallay : mrcc (arbitrary order and MRCC methods) + Jiri Pittner, Ondrej Demel : Mk-CCSD + Frank Weinhold : gennbo (NPA and NBO analysis) + Christopher J. Cramer and Donald G. Truhlar : smd solvation model + Lars Goerigk : TD-DFT with DH, B97 family of functionals + V. Asgeirsson, H. Jonsson : NEB implementation + FAccTs GmbH : IRC, NEB, NEB-TS, DLPNO-Multilevel, CI-OPT + MM, QMMM, 2- and 3-layer-ONIOM, Crystal-QMMM, + LR-CPCM, SF, NACMEs, symmetry and pop. for TD-DFT, + nearIR, NL-DFT gradient (VV10), updates on ESD, + ML-optimized integration grids + S Lehtola, MJT Oliveira, MAL Marques : LibXC Library + Liviu Ungur et al : ANISO software + + + Your calculation uses the libint2 library for the computation of 2-el integrals + For citations please refer to: http://libint.valeyev.net + + Your ORCA version has been built with support for libXC version: 5.1.0 + For citations please refer to: https://tddft.org/programs/libxc/ + + This ORCA versions uses: + CBLAS interface : Fast vector & matrix operations + LAPACKE interface : Fast linear algebra routines + SCALAPACK package : Parallel linear algebra routines + Shared memory : Shared parallel matrices + BLAS/LAPACK : OpenBLAS 0.3.15 USE64BITINT DYNAMIC_ARCH NO_AFFINITY Zen SINGLE_THREADED + Core in use : Zen + Copyright (c) 2011-2014, The OpenBLAS Project + + +================================================================================ + +----- Orbital basis set information ----- +Your calculation utilizes the basis: aug-cc-pVTZ + H, B-Ne : Obtained from the ccRepo (grant-hill.group.shef.ac.uk/ccrepo) Feb. 2017 + R. A. Kendall, T. H. Dunning, Jr., R. J. Harrison, J. Chem. Phys. 96, 6796 (1992) + He : Obtained from the ccRepo (grant-hill.group.shef.ac.uk/ccrepo) Feb. 2017 + D. E. Woon, T. H. Dunning, Jr., J. Chem. Phys. 100, 2975 (1994) + Li-Be, Na : Obtained from the ccRepo (grant-hill.group.shef.ac.uk/ccrepo) Feb. 2017 + B. P. Prascher, D. E. Woon, K. A. Peterson, T. H. Dunning, Jr., A. K. Wilson, Theor. Chem. Acc. 128, 69 (2011) + Mg : Obtained from the Peterson Research Group Website (tyr0.chem.wsu.edu/~kipeters) Feb. 2017 + B. P. Prascher, D. E. Woon, K. A. Peterson, T. H. Dunning, Jr., A. K. Wilson, Theor. Chem. Acc. 128, 69 (2011) + Al-Ar : Obtained from the ccRepo (grant-hill.group.shef.ac.uk/ccrepo) Feb. 2017 + D. E. Woon, T. H. Dunning, Jr., J. Chem. Phys. 98, 1358 (1993) + Sc-Zn : Obtained from the ccRepo (grant-hill.group.shef.ac.uk/ccrepo) Feb. 2017 + N. B. Balabanov, K. A. Peterson, J. Chem. Phys. 123, 064107 (2005) + N. B. Balabanov, K. A. Peterson, J. Chem. Phys. 125, 074110 (2006) + Ga-Kr : Obtained from the ccRepo (grant-hill.group.shef.ac.uk/ccrepo) Feb. 2017 + A. K. Wilson, D. E. Woon, K. A. Peterson, T. H. Dunning, Jr., J. Chem. Phys. 110, 7667 (1999) + Ag, Au : Obtained from the Peterson Research Group Website (tyr0.chem.wsu.edu/~kipeters) Feb. 2017 + K. A. Peterson, C. Puzzarini, Theor. Chem. Acc. 114, 283 (2005) + +================================================================================ + WARNINGS + Please study these warnings very carefully! +================================================================================ + + +INFO : the flag for use of the SHARK integral package has been found! + +================================================================================ + INPUT FILE +================================================================================ +NAME = input.in +| 1> !uHF aug-cc-pvtz tightscf +| 2> !sp +| 3> +| 4> %maxcore 3200 +| 5> %pal nprocs 16 end +| 6> +| 7> * xyz 0 4 +| 8> N 0.00000000 0.00000000 0.00000000 +| 9> * +| 10> +| 11> %scf +| 12> MaxIter 999 +| 13> end +| 14> +| 15> +| 16> %mp2 +| 17> RI true +| 18> end +| 19> +| 20> %casscf +| 21> nel 5 +| 22> norb 4 +| 23> nroots 1 +| 24> maxiter 999 +| 25> end +| 26> +| 27> %mrci +| 28> citype MRCI +| 29> davidsonopt true +| 30> maxiter 999 +| 31> end +| 32> +| 33> +| 34> ****END OF INPUT**** +================================================================================ + + **************************** + * Single Point Calculation * + **************************** + +--------------------------------- +CARTESIAN COORDINATES (ANGSTROEM) +--------------------------------- + N 0.000000 0.000000 0.000000 + +---------------------------- +CARTESIAN COORDINATES (A.U.) +---------------------------- + NO LB ZA FRAG MASS X Y Z + 0 N 7.0000 0 14.007 0.000000 0.000000 0.000000 + +-------------------------------- +INTERNAL COORDINATES (ANGSTROEM) +-------------------------------- + N 0 0 0 0.000000000000 0.00000000 0.00000000 + +--------------------------- +INTERNAL COORDINATES (A.U.) +--------------------------- + N 0 0 0 0.000000000000 0.00000000 0.00000000 + +--------------------- +BASIS SET INFORMATION +--------------------- +There are 1 groups of distinct atoms + + Group 1 Type N : 19s6p3d2f contracted to 5s4p3d2f pattern {88111/3111/111/11} + +Atom 0N basis set group => 1 + + + ************************************************************ + * Program running with 16 parallel MPI-processes * + * working on a common directory * + ************************************************************ +------------------------------------------------------------------------------ + ORCA GTO INTEGRAL CALCULATION +------------------------------------------------------------------------------ +------------------------------------------------------------------------------ + ___ + / \ - P O W E R E D B Y - + / \ + | | | _ _ __ _____ __ __ + | | | | | | | / \ | _ \ | | / | + \ \/ | | | | / \ | | | | | | / / + / \ \ | |__| | / /\ \ | |_| | | |/ / + | | | | __ | / /__\ \ | / | \ + | | | | | | | | __ | | \ | |\ \ + \ / | | | | | | | | | |\ \ | | \ \ + \___/ |_| |_| |__| |__| |_| \__\ |__| \__/ + + - O R C A' S B I G F R I E N D - + & + - I N T E G R A L F E E D E R - + + v1 FN, 2020, v2 2021 +------------------------------------------------------------------------------ + + +Reading SHARK input file input.SHARKINP.tmp ... ok +---------------------- +SHARK INTEGRAL PACKAGE +---------------------- + +Number of atoms ... 1 +Number of basis functions ... 52 +Number of shells ... 20 +Maximum angular momentum ... 3 +Integral batch strategy ... SHARK/LIBINT Hybrid +RI-J (if used) integral strategy ... SPLIT-RIJ (Revised 2003 algorithm where possible) +Printlevel ... 1 +Contraction scheme used ... PARTIAL GENERAL contraction +Coulomb Range Separation ... NOT USED +Exchange Range Separation ... NOT USED +Finite Nucleus Model ... NOT USED +Auxiliary Coulomb fitting basis ... NOT available +Auxiliary J/K fitting basis ... NOT available +Auxiliary Correlation fitting basis ... NOT available +Auxiliary 'external' fitting basis ... NOT available +Integral threshold ... 2.500000e-11 +Primitive cut-off ... 2.500000e-12 +Primitive pair pre-selection threshold ... 2.500000e-12 + +Calculating pre-screening integrals ... done ( 0.0 sec) Dimension = 20 +Calculating pre-screening integrals (ORCA) ... done ( 0.1 sec) Dimension = 14 +Organizing shell pair data ... done ( 0.0 sec) +Shell pair information +Total number of shell pairs ... 210 +Shell pairs after pre-screening ... 210 +Total number of primitive shell pairs ... 256 +Primitive shell pairs kept ... 256 + la=0 lb=0: 66 shell pairs + la=1 lb=0: 44 shell pairs + la=1 lb=1: 10 shell pairs + la=2 lb=0: 33 shell pairs + la=2 lb=1: 12 shell pairs + la=2 lb=2: 6 shell pairs + la=3 lb=0: 22 shell pairs + la=3 lb=1: 8 shell pairs + la=3 lb=2: 6 shell pairs + la=3 lb=3: 3 shell pairs + +Calculating one electron integrals ... done ( 0.0 sec) +Calculating Nuclear repulsion ... done ( 0.0 sec) ENN= 0.000000000000 Eh + +SHARK setup successfully completed in 0.2 seconds + +Maximum memory used throughout the entire GTOINT-calculation: 16.0 MB + + + ************************************************************ + * Program running with 16 parallel MPI-processes * + * working on a common directory * + ************************************************************ +------------------------------------------------------------------------------- + ORCA SCF +------------------------------------------------------------------------------- + +------------ +SCF SETTINGS +------------ +Hamiltonian: + Ab initio Hamiltonian Method .... Hartree-Fock(GTOs) + + +General Settings: + Integral files IntName .... input + Hartree-Fock type HFTyp .... CASSCF + Total Charge Charge .... 0 + Multiplicity Mult .... 4 + Number of Electrons NEL .... 7 + Basis Dimension Dim .... 46 + Nuclear Repulsion ENuc .... 0.0000000000 Eh + + +Diagonalization of the overlap matrix: +Smallest eigenvalue ... 2.690e-02 +Time for diagonalization ... 0.000 sec +Threshold for overlap eigenvalues ... 1.000e-08 +Number of eigenvalues below threshold ... 0 +Time for construction of square roots ... 0.001 sec +Total time needed ... 0.001 sec + +Time for model grid setup = 0.007 sec + +------------------------------ +INITIAL GUESS: MODEL POTENTIAL +------------------------------ +Loading Hartree-Fock densities ... done +Calculating cut-offs ... done +Initializing the effective Hamiltonian ... done +Setting up the integral package (SHARK) ... done +Starting the Coulomb interaction ... done ( 0.0 sec) +Reading the grid ... done +Mapping shells ... done +Starting the XC term evaluation ... done ( 0.0 sec) +Transforming the Hamiltonian ... done ( 0.0 sec) +Diagonalizing the Hamiltonian ... done ( 0.0 sec) +Back transforming the eigenvectors ... done ( 0.0 sec) +Now organizing SCF variables ... done + ------------------ + INITIAL GUESS DONE ( 0.0 sec) + ------------------ + + + ... the calculation is a CASSCF calculation -I'm leaving here GOOD LUCK!!! + + + + ************************************************************ + * Program running with 16 parallel MPI-processes * + * working on a common directory * + ************************************************************ +------------------------------------------------------------------------------- + ORCA-CASSCF +------------------------------------------------------------------------------- + +Setting up the integral package ... done +Building the CAS space ... done (4 configurations for Mult=4) + +SYSTEM-SPECIFIC SETTINGS: +Number of active electrons ... 5 +Number of active orbitals ... 4 +Total number of electrons ... 7 +Total number of orbitals ... 46 + +Determined orbital ranges: + Internal 0 - 0 ( 1 orbitals) + Active 1 - 4 ( 4 orbitals) + External 5 - 45 ( 41 orbitals) +Number of rotation parameters ... 209 + +CI-STEP: +CI strategy ... General CI +Number of multiplicity blocks ... 1 +BLOCK 1 WEIGHT= 1.0000 + Multiplicity ... 4 + #(Configurations) ... 4 + #(CSFs) ... 4 + #(Roots) ... 1 + ROOT=0 WEIGHT= 1.000000 + + PrintLevel ... 1 + N(GuessMat) ... 512 + MaxDim(CI) ... 10 + MaxIter(CI) ... 64 + Energy Tolerance CI ... 2.50e-09 + Residual Tolerance CI ... 2.50e-09 + Shift(CI) ... 1.00e-04 + +INTEGRAL-TRANSFORMATION-STEP: + Algorithm ... EXACT + +ORBITAL-IMPROVEMENT-STEP: + Algorithm ... SuperCI(PT) + Default Parametrization ... CAYLEY + Act-Act rotations ... depends on algorithm used + + Note: SuperCI(PT) will ignore FreezeIE, FreezeAct and FreezeGrad. In general Default settings are encouraged. + In conjunction with ShiftUp, ShiftDn or GradScaling the performance of SuperCI(PT) is less optimal. + + MaxRot ... 2.00e-01 + Max. no of vectors (DIIS) ... 15 + DThresh (cut-off) metric ... 1.00e-06 + Switch step at gradient ... 3.00e-02 + Switch step at iteration ... 50 + Switch step to ... SuperCI(PT) + +SCF-SETTINGS: + Incremental ... on + RIJCOSX approximation ... off + RI-JK approximation ... off + AO integral handling ... DIRECT + Integral Neglect Thresh ... 2.50e-11 + Primitive cutoff TCut ... 2.50e-12 + Energy convergence tolerance ... 2.50e-08 + Orbital gradient convergence ... 2.50e-04 + Max. number of iterations ... 999 + + +FINAL ORBITALS: + Active Orbitals ... natural + Internal Orbitals ... canonical + External Orbitals ... canonical + +------------------ +CAS-SCF ITERATIONS +------------------ + + +MACRO-ITERATION 1: + --- Inactive Energy E0 = -44.72092697 Eh +CI-ITERATION 0: + -54.381848153 0.000000000000 ( 0.00) + CI-PROBLEM SOLVED + DENSITIES MADE + + <<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>> + +BLOCK 1 MULT= 4 NROOTS= 1 +ROOT 0: E= -54.3818481526 Eh + 1.00000 [ 3]: 2111 + + <<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>> + + E(CAS)= -54.381848153 Eh DE= 0.000000e+00 + --- Energy gap subspaces: Ext-Act = 0.358 Act-Int = 14.743 + N(occ)= 2.00000 1.00000 1.00000 1.00000 + ||g|| = 8.377747e-01 Max(G)= -7.011926e-01 Rot=45,0 + --- Orbital Update [SuperCI(PT)] + --- Canonicalize Internal Space + --- Canonicalize External Space + --- SX_PT (Skipped TA=0 IT=1): ||X|| = 0.097273197 Max(X)(6,2) = -0.044359639 + --- SFit(Active Orbitals) + +MACRO-ITERATION 2: + --- Inactive Energy E0 = -44.72879418 Eh +CI-ITERATION 0: + -54.396929340 0.000000000000 ( 0.00) + CI-PROBLEM SOLVED + DENSITIES MADE + E(CAS)= -54.396929340 Eh DE= -1.508119e-02 + --- Energy gap subspaces: Ext-Act = 0.317 Act-Int = 14.676 + N(occ)= 2.00000 1.00000 1.00000 1.00000 + ||g|| = 8.661958e-02 Max(G)= -4.686820e-02 Rot=45,1 + --- Orbital Update [SuperCI(PT)] + --- Canonicalize Internal Space + --- Canonicalize External Space + --- SX_PT (Skipped TA=0 IT=1): ||X|| = 0.018251693 Max(X)(6,2) = -0.008989526 + --- SFit(Active Orbitals) + +MACRO-ITERATION 3: + --- Inactive Energy E0 = -44.72888355 Eh +CI-ITERATION 0: + -54.397580561 0.000000000000 ( 0.00) + CI-PROBLEM SOLVED + DENSITIES MADE + E(CAS)= -54.397580561 Eh DE= -6.512206e-04 + --- Energy gap subspaces: Ext-Act = 0.318 Act-Int = 14.686 + N(occ)= 2.00000 1.00000 1.00000 1.00000 + ||g|| = 1.169199e-02 Max(G)= -3.787582e-03 Rot=17,0 + --- Orbital Update [SuperCI(PT)] + --- Canonicalize Internal Space + --- Canonicalize External Space + --- SX_PT (Skipped TA=0 IT=1): ||X|| = 0.004765813 Max(X)(8,2) = 0.002361881 + --- SFit(Active Orbitals) + +MACRO-ITERATION 4: + --- Inactive Energy E0 = -44.72891366 Eh +CI-ITERATION 0: + -54.397608742 0.000000000000 ( 0.00) + CI-PROBLEM SOLVED + DENSITIES MADE + E(CAS)= -54.397608742 Eh DE= -2.818090e-05 + --- Energy gap subspaces: Ext-Act = 0.316 Act-Int = 14.684 + N(occ)= 2.00000 1.00000 1.00000 1.00000 + ||g|| = 3.401689e-03 Max(G)= 1.928798e-03 Rot=45,0 + --- Orbital Update [SuperCI(PT)] + --- Canonicalize Internal Space + --- Canonicalize External Space + --- SX_PT (Skipped TA=0 IT=1): ||X|| = 0.000448352 Max(X)(5,1) = -0.000198017 + --- SFit(Active Orbitals) + +MACRO-ITERATION 5: + --- Inactive Energy E0 = -44.72890218 Eh +CI-ITERATION 0: + -54.397609514 0.000000000000 ( 0.00) + CI-PROBLEM SOLVED + DENSITIES MADE + E(CAS)= -54.397609514 Eh DE= -7.726311e-07 + --- Energy gap subspaces: Ext-Act = 0.316 Act-Int = 14.684 + N(occ)= 2.00000 1.00000 1.00000 1.00000 + ||g|| = 3.651368e-04 Max(G)= -2.086419e-04 Rot=45,1 + --- Orbital Update [SuperCI(PT)] + --- Canonicalize Internal Space + --- Canonicalize External Space + --- SX_PT (Skipped TA=0 IT=1): ||X|| = 0.000048550 Max(X)(5,1) = -0.000029015 + --- SFit(Active Orbitals) + +MACRO-ITERATION 6: + --- Inactive Energy E0 = -44.72890337 Eh +CI-ITERATION 0: + -54.397609522 0.000000000000 ( 0.00) + CI-PROBLEM SOLVED + DENSITIES MADE + E(CAS)= -54.397609522 Eh DE= -7.825491e-09 + --- Energy gap subspaces: Ext-Act = 0.316 Act-Int = 14.684 + N(occ)= 2.00000 1.00000 1.00000 1.00000 + ||g|| = 5.710511e-05 Max(G)= 3.612727e-05 Rot=45,1 + ---- THE CAS-SCF ENERGY HAS CONVERGED ---- + ---- THE CAS-SCF GRADIENT HAS CONVERGED ---- + --- FINALIZING ORBITALS --- + ---- DOING ONE FINAL ITERATION FOR PRINTING ---- + --- Forming Natural Orbitals + --- Canonicalize Internal Space + --- Canonicalize External Space + +MACRO-ITERATION 7: + --- Inactive Energy E0 = -44.72890337 Eh + --- All densities will be recomputed +CI-ITERATION 0: + -54.397609522 0.000000000000 ( 0.00) + CI-PROBLEM SOLVED + DENSITIES MADE + E(CAS)= -54.397609522 Eh DE= -1.004565e-10 + --- Energy gap subspaces: Ext-Act = 0.311 Act-Int = 14.684 + N(occ)= 2.00000 1.00000 1.00000 1.00000 + ||g|| = 5.710508e-05 Max(G)= 3.643662e-05 Rot=45,1 +-------------- +CASSCF RESULTS +-------------- + +Final CASSCF energy : -54.397609522 Eh -1480.2342 eV + +---------------- +ORBITAL ENERGIES +---------------- + + NO OCC E(Eh) E(eV) + 0 2.0000 -15.630414 -425.3252 + 1 2.0000 -0.946888 -25.7661 + 2 1.0000 -0.180692 -4.9169 + 3 1.0000 -0.180692 -4.9169 + 4 1.0000 -0.180692 -4.9169 + 5 0.0000 0.129931 3.5356 + 6 0.0000 0.137198 3.7333 + 7 0.0000 0.137198 3.7333 + 8 0.0000 0.137198 3.7333 + 9 0.0000 0.464734 12.6461 + 10 0.0000 0.464734 12.6461 + 11 0.0000 0.464734 12.6461 + 12 0.0000 0.464734 12.6461 + 13 0.0000 0.464734 12.6461 + 14 0.0000 0.779518 21.2118 + 15 0.0000 0.779518 21.2118 + 16 0.0000 0.779518 21.2118 + 17 0.0000 1.009605 27.4727 + 18 0.0000 1.496268 40.7155 + 19 0.0000 1.496268 40.7155 + 20 0.0000 1.496268 40.7155 + 21 0.0000 1.496268 40.7155 + 22 0.0000 1.496268 40.7155 + 23 0.0000 1.496268 40.7155 + 24 0.0000 1.496268 40.7155 + 25 0.0000 1.559115 42.4257 + 26 0.0000 1.559115 42.4257 + 27 0.0000 1.559115 42.4257 + 28 0.0000 1.559115 42.4257 + 29 0.0000 1.559115 42.4257 + 30 0.0000 3.350863 91.1816 + 31 0.0000 3.350863 91.1816 + 32 0.0000 3.350863 91.1816 + 33 0.0000 4.877040 132.7110 + 34 0.0000 4.877040 132.7110 + 35 0.0000 4.877040 132.7110 + 36 0.0000 4.877040 132.7110 + 37 0.0000 4.877040 132.7110 + 38 0.0000 4.877040 132.7110 + 39 0.0000 4.877040 132.7110 + 40 0.0000 5.059015 137.6628 + 41 0.0000 5.059015 137.6628 + 42 0.0000 5.059015 137.6628 + 43 0.0000 5.059015 137.6628 + 44 0.0000 5.059015 137.6628 + 45 0.0000 6.334698 172.3759 + + +--------------------------------------------- +CAS-SCF STATES FOR BLOCK 1 MULT= 4 NROOTS= 1 +--------------------------------------------- + +ROOT 0: E= -54.3976095224 Eh + 1.00000 [ 3]: 2111 + + +-------------- +DENSITY MATRIX +-------------- + + 0 1 2 3 + 0 2.000000 0.000000 0.000000 0.000000 + 1 0.000000 1.000000 0.000000 0.000000 + 2 0.000000 0.000000 1.000000 0.000000 + 3 0.000000 0.000000 0.000000 1.000000 +Trace of the electron density: 5.000000 +Extracting Spin-Density from 2-RDM (MULT=4) ... done + +------------------- +SPIN-DENSITY MATRIX +------------------- + + 0 1 2 3 + 0 0.000000 0.000000 0.000000 0.000000 + 1 0.000000 1.000000 0.000000 0.000000 + 2 0.000000 0.000000 1.000000 0.000000 + 3 0.000000 0.000000 0.000000 1.000000 +Trace of the spin density: 3.000000 + +----------------- +ENERGY COMPONENTS +----------------- + +One electron energy : -73.937937606 Eh -2011.9536 eV +Two electron energy : 19.540328083 Eh 531.7194 eV +Nuclear repulsion energy : 0.000000000 Eh 0.0000 eV + ---------------- + -54.397609522 + +Kinetic energy : 54.383096271 Eh 1479.8393 eV +Potential energy : -108.780705794 Eh -2960.0735 eV +Virial ratio : -2.000266871 + ---------------- + -54.397609522 + +Core energy : -44.728903374 Eh -1217.1353 eV + + +---------------------------- +LOEWDIN ORBITAL-COMPOSITIONS +---------------------------- + + 0 1 2 3 4 5 + -15.63041 -0.94689 -0.18069 -0.18069 -0.18069 0.12993 + 2.00000 2.00000 1.00000 1.00000 1.00000 0.00000 + -------- -------- -------- -------- -------- -------- + 0 N s 100.0 100.0 0.0 0.0 0.0 100.0 + 0 N pz 0.0 0.0 98.7 1.3 0.0 0.0 + 0 N px 0.0 0.0 1.3 98.7 0.0 0.0 + 0 N py 0.0 0.0 0.0 0.0 100.0 0.0 + + 6 7 8 9 10 11 + 0.13720 0.13720 0.13720 0.46473 0.46473 0.46473 + 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + -------- -------- -------- -------- -------- -------- + 0 N pz 0.6 28.8 70.6 0.0 0.0 0.0 + 0 N px 0.2 70.6 29.2 0.0 0.0 0.0 + 0 N py 99.2 0.6 0.2 0.0 0.0 0.0 + 0 N dz2 0.0 0.0 0.0 0.0 1.2 86.4 + 0 N dxz 0.0 0.0 0.0 69.8 23.3 0.0 + 0 N dyz 0.0 0.0 0.0 3.5 0.0 11.6 + 0 N dx2y2 0.0 0.0 0.0 12.1 66.6 0.9 + 0 N dxy 0.0 0.0 0.0 14.6 8.9 1.2 + + 12 13 14 15 16 17 + 0.46473 0.46473 0.77952 0.77952 0.77952 1.00960 + 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + -------- -------- -------- -------- -------- -------- + 0 N s 0.0 0.0 0.0 0.0 0.0 100.0 + 0 N pz 0.0 0.0 0.0 27.3 72.7 0.0 + 0 N px 0.0 0.0 0.1 72.6 27.3 0.0 + 0 N py 0.0 0.0 99.9 0.0 0.0 0.0 + 0 N dz2 0.9 11.5 0.0 0.0 0.0 0.0 + 0 N dxz 3.6 3.3 0.0 0.0 0.0 0.0 + 0 N dyz 1.2 83.7 0.0 0.0 0.0 0.0 + 0 N dx2y2 19.3 1.1 0.0 0.0 0.0 0.0 + 0 N dxy 75.0 0.3 0.0 0.0 0.0 0.0 + + 18 19 20 21 22 23 + 1.49627 1.49627 1.49627 1.49627 1.49627 1.49627 + 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + -------- -------- -------- -------- -------- -------- + 0 N f0 5.6 24.2 37.4 21.3 4.3 2.9 + 0 N f+1 46.0 0.9 0.5 24.9 1.5 4.7 + 0 N f-1 18.8 40.5 9.7 5.5 0.0 0.5 + 0 N f+2 2.5 2.5 11.3 4.0 70.8 8.9 + 0 N f-2 1.7 9.1 3.0 18.4 20.1 29.2 + 0 N f+3 22.5 3.1 31.4 24.5 0.5 17.9 + 0 N f-3 2.8 19.7 6.7 1.3 2.8 35.9 + + 24 25 26 27 28 29 + 1.49627 1.55912 1.55912 1.55912 1.55912 1.55912 + 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + -------- -------- -------- -------- -------- -------- + 0 N dz2 0.0 5.6 15.9 74.3 1.7 2.5 + 0 N dxz 0.0 12.7 2.5 0.9 6.7 77.2 + 0 N dyz 0.0 5.3 37.5 20.9 35.1 1.2 + 0 N dx2y2 0.0 65.9 18.1 0.9 0.0 15.1 + 0 N dxy 0.0 10.6 26.0 3.0 56.4 4.0 + 0 N f0 4.2 0.0 0.0 0.0 0.0 0.0 + 0 N f+1 21.5 0.0 0.0 0.0 0.0 0.0 + 0 N f-1 25.0 0.0 0.0 0.0 0.0 0.0 + 0 N f-2 18.5 0.0 0.0 0.0 0.0 0.0 + 0 N f-3 30.8 0.0 0.0 0.0 0.0 0.0 + + 30 31 32 33 34 35 + 3.35086 3.35086 3.35086 4.87704 4.87704 4.87704 + 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + -------- -------- -------- -------- -------- -------- + 0 N pz 23.6 12.5 63.8 0.0 0.0 0.0 + 0 N px 56.8 8.5 34.6 0.0 0.0 0.0 + 0 N py 19.5 78.9 1.6 0.0 0.0 0.0 + 0 N f0 0.0 0.0 0.0 5.7 0.4 27.9 + 0 N f+1 0.0 0.0 0.0 5.5 0.0 0.0 + 0 N f-1 0.0 0.0 0.0 0.5 2.4 0.5 + 0 N f+2 0.0 0.0 0.0 36.3 53.6 2.9 + 0 N f-2 0.0 0.0 0.0 5.6 1.5 55.0 + 0 N f+3 0.0 0.0 0.0 12.8 0.1 0.3 + 0 N f-3 0.0 0.0 0.0 33.6 42.0 13.5 + + 36 37 38 39 40 41 + 4.87704 4.87704 4.87704 4.87704 5.05902 5.05902 + 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 + -------- -------- -------- -------- -------- -------- + 0 N dz2 0.0 0.0 0.0 0.0 44.5 6.2 + 0 N dxz 0.0 0.0 0.0 0.0 1.2 0.0 + 0 N dyz 0.0 0.0 0.0 0.0 3.8 33.6 + 0 N dx2y2 0.0 0.0 0.0 0.0 47.3 7.1 + 0 N dxy 0.0 0.0 0.0 0.0 3.2 53.1 + 0 N f0 8.9 0.1 3.3 53.6 0.0 0.0 + 0 N f+1 10.6 64.8 19.1 0.0 0.0 0.0 + 0 N f-1 1.1 26.1 66.5 3.0 0.0 0.0 + 0 N f+2 5.1 0.8 0.1 1.3 0.0 0.0 + 0 N f-2 0.1 0.0 1.3 36.6 0.0 0.0 + 0 N f+3 66.4 8.1 7.1 5.2 0.0 0.0 + 0 N f-3 7.9 0.1 2.6 0.3 0.0 0.0 + + 42 43 44 45 + 5.05902 5.05902 5.05902 6.33470 + 0.00000 0.00000 0.00000 0.00000 + -------- -------- -------- -------- + 0 N s 0.0 0.0 0.0 100.0 + 0 N dz2 32.3 16.9 0.0 0.0 + 0 N dxz 1.5 0.1 97.2 0.0 + 0 N dyz 19.7 41.9 1.1 0.0 + 0 N dx2y2 5.9 39.0 0.8 0.0 + 0 N dxy 40.7 2.1 0.9 0.0 + +---------------------------- +LOEWDIN REDUCED ACTIVE MOs +---------------------------- + + 0 1 2 3 4 5 + -15.63041 -0.94689 -0.18069 -0.18069 -0.18069 0.12993 + 2.00000 2.00000 1.00000 1.00000 1.00000 0.00000 + -------- -------- -------- -------- -------- -------- + 0 N s 100.0 100.0 0.0 0.0 0.0 100.0 + 0 N pz 0.0 0.0 98.7 1.3 0.0 0.0 + 0 N px 0.0 0.0 1.3 98.7 0.0 0.0 + 0 N py 0.0 0.0 0.0 0.0 100.0 0.0 + +------------------------------------------------------------------------------ + ORCA POPULATION ANALYSIS +------------------------------------------------------------------------------ +Input electron density ... input.scfp +Input spin density ... input.scfr +BaseName (.gbw .S,...) ... input + + ******************************** + * MULLIKEN POPULATION ANALYSIS * + ******************************** + +------------------------------------------ +MULLIKEN ATOMIC CHARGES AND SPIN DENSITIES +------------------------------------------ + 0 N : -0.000000 3.000000 +Sum of atomic charges : -0.0000000 +Sum of atomic spin densities: 3.0000000 + +--------------------------------------------------- +MULLIKEN REDUCED ORBITAL CHARGES AND SPIN DENSITIES +--------------------------------------------------- +CHARGE + 0 N s : 4.000000 s : 4.000000 + pz : 1.000000 p : 3.000000 + px : 1.000000 + py : 1.000000 + dz2 : 0.000000 d : 0.000000 + dxz : 0.000000 + dyz : 0.000000 + dx2y2 : 0.000000 + dxy : 0.000000 + f0 : 0.000000 f : 0.000000 + f+1 : 0.000000 + f-1 : 0.000000 + f+2 : 0.000000 + f-2 : 0.000000 + f+3 : 0.000000 + f-3 : 0.000000 + +SPIN + 0 N s : 0.000000 s : 0.000000 + pz : 1.000000 p : 3.000000 + px : 1.000000 + py : 1.000000 + dz2 : 0.000000 d : 0.000000 + dxz : 0.000000 + dyz : 0.000000 + dx2y2 : 0.000000 + dxy : 0.000000 + f0 : 0.000000 f : 0.000000 + f+1 : 0.000000 + f-1 : 0.000000 + f+2 : 0.000000 + f-2 : 0.000000 + f+3 : 0.000000 + f-3 : 0.000000 + + + ******************************* + * LOEWDIN POPULATION ANALYSIS * + ******************************* + +----------------------------------------- +LOEWDIN ATOMIC CHARGES AND SPIN DENSITIES +----------------------------------------- + 0 N : -0.000000 3.000000 + +-------------------------------------------------- +LOEWDIN REDUCED ORBITAL CHARGES AND SPIN DENSITIES +-------------------------------------------------- +CHARGE + 0 N s : 4.000000 s : 4.000000 + pz : 1.000000 p : 3.000000 + px : 1.000000 + py : 1.000000 + dz2 : 0.000000 d : 0.000000 + dxz : 0.000000 + dyz : 0.000000 + dx2y2 : 0.000000 + dxy : 0.000000 + f0 : 0.000000 f : 0.000000 + f+1 : 0.000000 + f-1 : 0.000000 + f+2 : 0.000000 + f-2 : 0.000000 + f+3 : 0.000000 + f-3 : 0.000000 + +SPIN + 0 N s : 0.000000 s : 0.000000 + pz : 1.000000 p : 3.000000 + px : 1.000000 + py : 1.000000 + dz2 : 0.000000 d : 0.000000 + dxz : 0.000000 + dyz : 0.000000 + dx2y2 : 0.000000 + dxy : 0.000000 + f0 : 0.000000 f : 0.000000 + f+1 : 0.000000 + f-1 : 0.000000 + f+2 : 0.000000 + f-2 : 0.000000 + f+3 : 0.000000 + f-3 : 0.000000 + + + ***************************** + * MAYER POPULATION ANALYSIS * + ***************************** + + NA - Mulliken gross atomic population + ZA - Total nuclear charge + QA - Mulliken gross atomic charge + VA - Mayer's total valence + BVA - Mayer's bonded valence + FA - Mayer's free valence + + ATOM NA ZA QA VA BVA FA + 0 N 7.0000 7.0000 -0.0000 3.0000 0.0000 3.0000 + + + +Transition Dipole Moments (a.u.) with input orbitals: +MO 0: + MO 2: 0.074753310 ( 0.008585038, 0.000026009, -0.074258695) + MO 3: 0.074753310 ( -0.074258616, -0.000109345, -0.008585067) + MO 4: 0.074753310 ( -0.000111609, 0.074753225, 0.000013280) +MO 1: + MO 2: 0.777510748 ( 0.089293158, 0.000270525, -0.772366248) + MO 3: 0.777510748 ( -0.772365423, -0.001137300, -0.089293461) + MO 4: 0.777510748 ( -0.001160843, 0.777509869, 0.000138122) +MO 2: +MO 3: +MO 4: + newgto N + 0 18 + 1 11420.000000000000 0.000522768355 + 2 1712.000000000000 0.004043208403 + 3 389.300000000000 0.020765798413 + 4 110.000000000000 0.080691244693 + 5 35.570000000000 0.232970767719 + 6 12.540000000000 0.433308995327 + 7 4.644000000000 0.347318098976 + 8 0.511800000000 -0.008504231668 + 9 11420.000000000000 -0.000000607139 + 10 1712.000000000000 -0.000004725125 + 11 389.300000000000 -0.000024412268 + 12 110.000000000000 -0.000097818014 + 13 35.570000000000 -0.000302719511 + 14 12.540000000000 -0.000697291235 + 15 4.644000000000 -0.000910761311 + 16 0.511800000000 0.003167386147 + 17 1.293000000000 0.042437946502 + 18 0.178700000000 0.006072346930 + 0 18 + 1 11420.000000000000 0.000002720308 + 2 1712.000000000000 0.000021039475 + 3 389.300000000000 0.000108058121 + 4 110.000000000000 0.000419889672 + 5 35.570000000000 0.001212300289 + 6 12.540000000000 0.002254791986 + 7 4.644000000000 0.001807324738 + 8 11420.000000000000 0.000115506572 + 9 1712.000000000000 0.000898942453 + 10 389.300000000000 0.004644368606 + 11 110.000000000000 0.018609615382 + 12 35.570000000000 0.057591576877 + 13 12.540000000000 0.132657791513 + 14 4.644000000000 0.173269902282 + 15 0.511800000000 -0.602586738479 + 16 1.293000000000 -0.151583562870 + 17 0.178700000000 -0.380619025879 + 18 0.057600000000 -0.006167374158 + 1 6 + 1 26.630000000000 0.014646107832 + 2 5.948000000000 0.091614549360 + 3 1.742000000000 0.298196552532 + 4 0.555000000000 0.500982554894 + 5 0.172500000000 0.327306035784 + 6 0.049100000000 0.015517087267 + end + +------------------------------------------------------------- + Forming the transition density ... done in 0.000164 sec +------------------------------------------------------------- + + + +========================================== +CASSCF UV, CD spectra and dipole moments +========================================== +------------------- +ABSORPTION SPECTRUM +------------------- + +Center of mass = ( 0.0000, 0.0000, 0.0000) +Nuclear contribution to the dipole moment = 0.000000, 0.000000, 0.000000 au + +Calculating the Dipole integrals ... done +Transforming integrals ... done +Calculating the Linear Momentum integrals ... done +Transforming integrals ... done +Calculating the Angular Momentum integrals ... done +Transforming integrals ... done + +------------------------------------------------------------------------------ + DIPOLE MOMENTS +------------------------------------------------------------------------------ + Root Block TX TY TZ |T| + (Debye) (Debye) (Debye) (Debye) +------------------------------------------------------------------------------ + 0 0 0.00000 0.00000 0.00000 0.00000 + +-------------- +CASSCF TIMINGS +-------------- + +Total time ... 7.2 sec +Sum of individual times ... 3.1 sec ( 42.6%) + +Calculation of AO operators + F(Core) operator ... 1.4 sec ( 18.9%) + G(Act) operator ... 0.3 sec ( 4.1%) + J(AO) operators ... 0.0 sec ( 0.0%) +Calculation of MO transformed quantities + J(MO) operators ... 0.1 sec ( 1.2%) + (pq|rs) integrals ... 0.0 sec ( 0.0%) + AO->MO one electron integrals ... 0.0 sec ( 0.0%) +Configuration interaction steps + CI-setup phase ... 0.0 sec ( 0.0%) + CI-solution phase ... 0.4 sec ( 6.0%) + Generation of densities ... 0.0 sec ( 0.0%) +Orbital improvement steps + Orbital gradient ... 0.9 sec ( 12.4%) + O(1) converger ... 0.0 sec ( 0.0%) +Properties ... 0.0 sec ( 0.0%) + SOC integral calculation ... 0.0 sec ( 0.0%) + SSC RMEs (incl. integrals) ... 0.0 sec ( 0.0%) + SOC RMEs ... 0.0 sec ( 0.0%) + +Maximum memory used throughout the entire CASSCF-calculation: 117.8 MB +Warning: no block have been defined for MRCI - copying CASSCF information! + + + ************************************************************ + * Program running with 16 parallel MPI-processes * + * working on a common directory * + ************************************************************ + +------------------------------------------------------------------------------ + ORCA MRCI +------------------------------------------------------------------------------ + +Transformed one-electron integrals ... input.cih.tmp +Transformed RI-integrals ... input.rijt.tmp +Output vectors ... input.mrci + +---------------- +MRCI-SETUP PHASE +---------------- + +GBW-Name ... input.gbw +IVO-Name ... input.ivo +HName ... input.H.tmp +SName ... input.S.tmp +CIHName ... input.cih.tmp +CIFName ... input.cif.tmp +MRCIName ... input.mrci +IntFiles ... input.rijt.tmp +LocName ... input.loc +MRCIInput ... input.mrciinp +Basis dimension ... 46 +Improved Virtual Orbitals ... OFF +Orbital localization ... OFF + +Figured out Orbital ranges +Set data to blocks + +----------------------- +FROZEN CORE FOCK MATRIX +----------------------- + +Improved virtual orbitals (IVOs) WILL NOT be constructed +Orbital Range is before PREP ... Int= 1- 0 Act= 1- 4 Ext= 5- 45 +Calculating Frozen Core Matrices ... (look at input.lastciprep) +Warning: FirstInternal=1 LastInternal=0 NInternals=0 - setting LastInternal for FI=-1 + +------------------------------------------------------------------------------ + ORCA CI-PREPARATION +------------------------------------------------------------------------------ +Reading the GBW file ... done + + +One-Electron Matrix ... input.H.tmp +GBW-File ... input.gbw +First MO in the CI ... 1 +Integral package used ... LIBINT +Reading the GBW file ... done + +Reading the one-electron matrix ... done +Forming inactive density ... done +Forming Fock matrix/matrices ... +Nuclear repulsion ... 0.000000 +Core repulsion ... 0.000000 +One-electron energy ... -48.843003 +Fock-energy ... -40.614803 +Final value ... -44.728903 +done +Transforming integrals ... done +Storing passive energy ... done ( -44.72890337 Eh) +The internal FI matrix is equal to the CIH matrix; storing it as such! + .... done with the Frozen Core Fock matrices +Final Orbital Range is now ... Int= 1- 0 Act= 1- 4 Ext= 5- 45 + +------------------------------ +PARTIAL COULOMB TRANSFORMATION +------------------------------ + +Dimension of the basis ... 46 +Number of internal MOs ... 45 (1-45) +Pair cutoff ... 2.500e-12 Eh +Number of AO pairs included in the trafo ... 1081 +Total Number of distinct AO pairs ... 1081 +Memory devoted for trafo ... 3200 MB +Memory needed per MO pair ... 0 MB +Max. Number of MO pairs treated together ... 388002 +Memory needed per MO ... 0 MB +Max. Number of MOs treated per batch ... 45 +Number Format for Storage ... Double (8 Byte) +Integral package used ... LIBINT + + --->>> The Coulomb operators (i,j|mue,nue) will be calculated + +Starting integral evaluation: +: : 1575 b 0 skpd 0.002 s ( 0.001 ms/b) +: 2100 b 0 skpd 0.002 s ( 0.001 ms/b) +: 1575 b 0 skpd 0.001 s ( 0.001 ms/b) +: 1050 b 0 skpd 0.002 s ( 0.002 ms/b) +: 1050 b 0 skpd 0.001 s ( 0.001 ms/b) +: 1260 b 0 skpd 0.001 s ( 0.001 ms/b) +: 840 b 0 skpd 0.002 s ( 0.003 ms/b) + 630 b 0 skpd 0.002 s ( 0.003 ms/b) +: 630 b 0 skpd 0.003 s ( 0.005 ms/b) +: 315 b 0 skpd 0.004 s ( 0.014 ms/b) +Collecting buffer AOJ + ... done with AO integral generation +Closing buffer AOJ ( 0.00 GB; CompressionRatio= 1.77) +Number of MO pairs included in the trafo ... 1035 + ... Now sorting integrals +IBATCH = 1 of 2 +IBATCH = 2 of 2 +Closing buffer JAO ( 0.01 GB; CompressionRatio= 1.00) +TOTAL TIME for half transformation ... 0.030 sec +AO-integral generation ... 0.020 sec +Half transformation ... 0.003 sec +J-integral sorting ... 0.007 sec +Collecting buffer JAO + +------------------- +FULL TRANSFORMATION +------------------- + +Processing MO 10 +Processing MO 20 +Processing MO 30 +Processing MO 40 +Full transformation done +Number of integrals made ... 536130 +Number of integrals stored ... 261508 +Timings: +Time for first half transformation ... 0.032 sec +Time for second half transformation ... 0.002 sec +Total time ... 0.037 sec + +------------------ +CI-BLOCK STRUCTURE +------------------ + +Number of CI-blocks ... 1 + +=========== +CI BLOCK 1 +=========== +Multiplicity ... 4 +Irrep ... -1 +Number of reference defs ... 1 + Reference 1: CAS(5,4) + +Excitation type ... CISD +Excitation flags for singles: + 1 1 1 1 +Excitation flags for doubles: + 1 1 1 / 1 1 1 / 1 1 1 + + + + -------------------------------------------------------------------- + -------------------- ALL SETUP TASKS ACCOMPLISHED ------------------ + -------------------- ( 0.258 sec) ------------------ + -------------------------------------------------------------------- + + + + ################################################### + # # + # M R C I # + # # + # TSel = 1.000e-06 Eh # + # TPre = 1.000e-04 # + # TIntCut = 1.000e-10 Eh # + # Extrapolation to unselected MR-CI by full MP2 # + # DAVIDSON-1 Correction to full CI # + # # + ################################################### + + +--------------------- +INTEGRAL ORGANIZATION +--------------------- + +Reading the one-Electron matrix ... done +E0 read was -44.728903374198 +Reading the internal Fock matrix ... done +Preparing the integral list ... done +Loading the full integral list ... done +Making the simple integrals ... done + + *************************************** + * CI-BLOCK 1 * + *************************************** + +Configurations with insufficient # of SOMOs WILL be rejected +Building a CAS(5,4) for multiplicity 4 +Reference Space: + Initial Number of Configurations : 4 + Internal Orbitals : 1 - 0 + Active Orbitals : 1 - 4 + External Orbitals : 5 - 45 +The number of CSFs in the reference is 4 +Calling MRPT_Selection with N(ref)=4 +------------------ +REFERENCE SPACE CI +------------------ + + Pre-diagonalization threshold : 1.000e-04 +Warning: Setting NGuessMat to 512 +N(ref-CFG)=4 N(ref-CSF)=4 + + ****Iteration 0**** + Lowest Energy : -54.397609522415 + Maximum Energy change : 54.397609522415 (vector 0) + Maximum residual norm : 0.000000000000 + + *** CONVERGENCE OF RESIDUAL NORM REACHED *** +Reference space selection using TPre= 1.00e-04 + + ... found 1 reference configurations (1 CSFs) + ... now redoing the reference space CI ... + +Warning: Setting NGuessMat to 512 +N(ref-CFG)=1 N(ref-CSF)=1 + + ****Iteration 0**** + Lowest Energy : -54.397609522415 + Maximum Energy change : 54.397609522415 (vector 0) + Maximum residual norm : 0.000000000000 + + *** CONVERGENCE OF RESIDUAL NORM REACHED *** + +---------- +CI-RESULTS +---------- + +The threshold for printing is 0.30 percent +The weights of configurations will be printed. The weights are summed over +all CSFs that belong to a given configuration before printing + +STATE 0: Energy= -54.397609522 Eh RefWeight= 1.0000 0.00 eV 0.0 cm**-1 + 1.0000 : h---h---[2111] + +------------------------------ +MR-PT SELECTION TSel= 1.00e-06 +------------------------------ + + +Setting reference configurations WITHOUT use of symmetry +Building active patterns WITHOUT use of symmetry + +Selection will be done from 1 spatial configurations + ( 0) Refs : Sel: 1CFGs/ 1CSFs Gen: 1CFGs/ 1CSFs ( 0.000 sec) +Building active space densities ... 0.002 sec +Building active space Fock operators ... 0.000 sec + ( 1) (p,q)->(r,s): Sel: 3CFGs/ 3CSFs Gen: 3CFGs/ 3CSFs ( 0.001 sec) + ( 5) (p,-)->(a,-): Sel: 164CFGs/ 287CSFs Gen: 164CFGs/ 287CSFs ( 0.001 sec) + ( 6) (i,-)->(a,-): Sel: 0CFGs/ 0CSFs Gen: 0CFGs/ 0CSFs ( 0.001 sec) + ( 7) (i,j)->(p,a): Sel: 0CFGs/ 0CSFs Gen: 0CFGs/ 0CSFs ( 0.000 sec) + ( 8) (i,p)->(q,a): Sel: 0CFGs/ 0CSFs Gen: 0CFGs/ 0CSFs ( 0.001 sec) + ( 9) (p,q)->(r,a): Sel: 108CFGs/ 108CSFs Gen: 369CFGs/ 369CSFs ( 0.001 sec) + (10) (i,p)->(a,b): Sel: 0CFGs/ 0CSFs Gen: 0CFGs/ 0CSFs ( 0.001 sec) + (11) (p,q)->(a,b): Sel: 1279CFGs/ 3493CSFs Gen: 5904CFGs/ 15744CSFs ( 0.002 sec) + (12) (i,j)->(a,b): Sel: 0CFGs/ 0CSFs Gen: 0CFGs/ 0CSFs ( 2.945 sec) + +Selection results: +Total number of generated configurations: 6441 +Number of selected configurations : 1555 ( 24.1%) +Total number of generated CSFs : 16404 +Number of selected CSFS : 3892 ( 23.7%) + +The selected tree structure: +Number of selected Internal Portions : 4 +Number of selected Singly External Portions: 13 + average number of VMOs/Portion : 19.43 + percentage of selected singly externals : 50.94 +Number of selected Doubly External Portions: 7 + average number of VMOs/Portion : 159.88 + percentage of selected doubly externals : 22.28 + +Diagonal second order perturbation results: +State E(tot) E(0)+E(1) E2(sel) E2(unsel) + Eh Eh Eh Eh +---------------------------------------------------------------- + 0 -54.503566265 -54.397609522 -0.105685 -0.000272 + +Computing the reference space Hamiltonian ... done (DIM=1) +Storing the reference space Hamiltonian ... done +Finding start and stop indices ... done +Collecting additional information ... done +Entering the DIIS solver section +------------------------ +MULTIROOT DIIS CI SOLVER +------------------------ + +Number of CSFs ... 3892 +Number of configurations ... 1555 +Maximum number of DIIS vectors stored ... 5 +Level shift ... 0.20 +Convergence tolerance on energies ... 2.500e-07 +Convergence tolerance on residual ... 2.500e-07 +Partitioning used ... MOELLER-PLESSET + + + ****Iteration 0**** + State 0: E= -54.466179124 Ec=-0.068569602 R= 0.299648465 W0= 0.96521 + Max energy change = 6.8570e-02 Eh + Max Residual Norm = 2.9965e-01 + + ****Iteration 1**** + State 0: E= -54.512269425 Ec=-0.114659903 R= 0.002512143 W0= 0.95744 + Max energy change = 4.6090e-02 Eh + Max Residual Norm = 2.5121e-03 + + ****Iteration 2**** + State 0: E= -54.512978160 Ec=-0.115368637 R= 0.000325639 W0= 0.96224 + Max energy change = 7.0873e-04 Eh + Max Residual Norm = 3.2564e-04 + + ****Iteration 3**** + State 0: E= -54.513035558 Ec=-0.115426036 R= 0.000008152 W0= 0.96143 + Max energy change = 5.7399e-05 Eh + Max Residual Norm = 8.1520e-06 + + ****Iteration 4**** + State 0: E= -54.513039684 Ec=-0.115430162 R= 0.000001364 W0= 0.96186 + Max energy change = 4.1259e-06 Eh + Max Residual Norm = 1.3636e-06 + + ****Iteration 5**** + State 0: residual converged + State 0: E= -54.513040333 Ec=-0.115430811 R= 0.000000068 W0= 0.96184 + Max energy change = 0.0000e+00 Eh + Max Residual Norm = 6.7607e-08 + *** Convergence of energies reached *** + *** Convergence of residual reached *** + *** All vectors converged *** +Returned from DIIS section + +---------------------------------------------- +MULTI-REFERENCE/MULTI-ROOT DAVIDSON CORRECTION +---------------------------------------------- + + +Summary of multireference corrections: + +Root W(ref) E(MR-CI) E(ref) Delta-E None Davidson-1 Davidson-2 Siegbahn Pople +------------------------------------------------------------------------------------------------------------ + 0 0.962 -54.513040333 -54.397609522 0.115430811 0.000000 -0.004405 -0.004769 -0.004580 -0.002769 +------------------------------------------------------------------------------------------------------------ +Active option = Davidson-1 + +Unselected CSF estimate: +Full relaxed MR-MP2 calculation ... + +Selection will be done from 1 spatial configurations + +Selection will be done from 1 spatial configurations + +Selection will be done from 1 spatial configurations +done +Selected MR-MP2 energies ... + + Root= 0 E(unsel)= -0.000272010 + +---------- +CI-RESULTS +---------- + +The threshold for printing is 0.30 percent +The weights of configurations will be printed. The weights are summed over +all CSFs that belong to a given configuration before printing + +STATE 0: Energy= -54.517717250 Eh RefWeight= 0.9618 0.00 eV 0.0 cm**-1 + 0.9618 : h---h---[2111] +Now choosing densities with flags StateDens=3 and TransDens=1 NStates(total)=1 +State density of the lowest state in each block NStateDens= 2 +GS to excited state transition electron densities NTransDens=0 +NDens(total)=2 +All lowest density information prepared Cnt(Dens)=2 +GS to ES state electron density information prepared Cnt(Dens)=2 + +------------------ +DENSITY GENERATION +------------------ + + ... generating densities (input.mrci.vec0) ... o.k. +MRCI-Population analysis: looping over 2 densities +Found state electron-density state=0 block=0 +Found state spin-density state=0 block=0 +------------------------------------------------------------------------------ + ORCA POPULATION ANALYSIS +------------------------------------------------------------------------------ +Input electron density ... input.state_0_block_0.el.tmp +Input spin density ... input.state_0_block_0.spin.tmp +BaseName (.gbw .S,...) ... input + + ******************************** + * MULLIKEN POPULATION ANALYSIS * + ******************************** + +------------------------------------------ +MULLIKEN ATOMIC CHARGES AND SPIN DENSITIES +------------------------------------------ + 0 N : -0.000000 3.000000 +Sum of atomic charges : -0.0000000 +Sum of atomic spin densities: 3.0000000 + +--------------------------------------------------- +MULLIKEN REDUCED ORBITAL CHARGES AND SPIN DENSITIES +--------------------------------------------------- +CHARGE + 0 N s : 3.975014 s : 3.975014 + pz : 0.998985 p : 2.996924 + px : 0.998969 + py : 0.998970 + dz2 : 0.005285 d : 0.026369 + dxz : 0.005270 + dyz : 0.005265 + dx2y2 : 0.005292 + dxy : 0.005257 + f0 : 0.000241 f : 0.001693 + f+1 : 0.000239 + f-1 : 0.000253 + f+2 : 0.000241 + f-2 : 0.000236 + f+3 : 0.000234 + f-3 : 0.000251 + +SPIN + 0 N s : 0.021881 s : 0.021881 + pz : 0.984931 p : 2.954801 + px : 0.984937 + py : 0.984933 + dz2 : 0.004517 d : 0.022525 + dxz : 0.004496 + dyz : 0.004491 + dx2y2 : 0.004526 + dxy : 0.004496 + f0 : 0.000113 f : 0.000793 + f+1 : 0.000112 + f-1 : 0.000119 + f+2 : 0.000112 + f-2 : 0.000109 + f+3 : 0.000110 + f-3 : 0.000117 + + + ******************************* + * LOEWDIN POPULATION ANALYSIS * + ******************************* + +----------------------------------------- +LOEWDIN ATOMIC CHARGES AND SPIN DENSITIES +----------------------------------------- + 0 N : -0.000000 3.000000 + +-------------------------------------------------- +LOEWDIN REDUCED ORBITAL CHARGES AND SPIN DENSITIES +-------------------------------------------------- +CHARGE + 0 N s : 3.975014 s : 3.975014 + pz : 0.998985 p : 2.996924 + px : 0.998969 + py : 0.998970 + dz2 : 0.005285 d : 0.026369 + dxz : 0.005270 + dyz : 0.005265 + dx2y2 : 0.005292 + dxy : 0.005257 + f0 : 0.000241 f : 0.001693 + f+1 : 0.000239 + f-1 : 0.000253 + f+2 : 0.000241 + f-2 : 0.000236 + f+3 : 0.000234 + f-3 : 0.000251 + +SPIN + 0 N s : 0.021881 s : 0.021881 + pz : 0.984931 p : 2.954801 + px : 0.984937 + py : 0.984933 + dz2 : 0.004517 d : 0.022525 + dxz : 0.004496 + dyz : 0.004491 + dx2y2 : 0.004526 + dxy : 0.004496 + f0 : 0.000113 f : 0.000793 + f+1 : 0.000112 + f-1 : 0.000119 + f+2 : 0.000112 + f-2 : 0.000109 + f+3 : 0.000110 + f-3 : 0.000117 + + + ***************************** + * MAYER POPULATION ANALYSIS * + ***************************** + + NA - Mulliken gross atomic population + ZA - Total nuclear charge + QA - Mulliken gross atomic charge + VA - Mayer's total valence + BVA - Mayer's bonded valence + FA - Mayer's free valence + + ATOM NA ZA QA VA BVA FA + 0 N 7.0000 7.0000 -0.0000 3.1783 0.0000 3.1783 + + + + +--------------------- +CI-EXCITATION SPECTRA +--------------------- + +Center of mass = ( 0.0000, 0.0000, 0.0000) + +Nuclear contribution to the dipole moment= 0.000000, 0.000000, 0.000000 au + +Calculating the Dipole integrals ... done +Transforming integrals ... done +Calculating the Linear Momentum integrals ... done +Transforming integrals ... done +Calculating the Angular momentum integrals ... done +Transforming integrals ... done + +------------------------------------------------------------------------------------------ + ABSORPTION SPECTRUM +------------------------------------------------------------------------------------------ + States Energy Wavelength fosc T2 TX TY TZ + (cm-1) (nm) (D**2) (D) (D) (D) +------------------------------------------------------------------------------------------ + +------------------------------------------------------------------------------ + CD SPECTRUM +------------------------------------------------------------------------------ + States Energy Wavelength R*T RX RY RZ + (cm-1) (nm) (1e40*sgs) (au) (au) (au) +------------------------------------------------------------------------------ + +------------------------------------------------------------------------------ + STATE DIPOLE MOMENTS +------------------------------------------------------------------------------ + Root Block TX TY TZ |T| + (Debye) (Debye) (Debye) (Debye) +------------------------------------------------------------------------------ + 0 0 0.00000 0.00000 0.00000 0.00000 + +Maximum memory used throughout the entire MRCI-calculation: 76.9 MB + +------------------------- -------------------- +FINAL SINGLE POINT ENERGY -54.517717249721 +------------------------- -------------------- + + + *************************************** + * ORCA property calculations * + *************************************** + + --------------------- + Active property flags + --------------------- + (+) Dipole Moment + + +------------------------------------------------------------------------------ + ORCA ELECTRIC PROPERTIES CALCULATION +------------------------------------------------------------------------------ + +Dipole Moment Calculation ... on +Quadrupole Moment Calculation ... off +Polarizability Calculation ... off +GBWName ... input.gbw +Electron density ... input.scfp +The origin for moment calculation is the CENTER OF MASS = ( 0.000000, 0.000000 0.000000) + +------------- +DIPOLE MOMENT +------------- + X Y Z +Electronic contribution: 0.00000 0.00000 0.00000 +Nuclear contribution : 0.00000 0.00000 0.00000 + ----------------------------------------- +Total Dipole Moment : 0.00000 0.00000 0.00000 + ----------------------------------------- +Magnitude (a.u.) : 0.00000 +Magnitude (Debye) : 0.00000 + + + +-------------------- +Rotational spectrum +-------------------- + +Rotational constants in cm-1: 0.000000 0.000000 0.000000 +Rotational constants in MHz : 0.000000 0.000000 0.000000 + + Dipole components along the rotational axes: +x,y,z [a.u.] : 0.000000 0.000000 0.000000 +x,y,z [Debye]: 0.000000 0.000000 0.000000 + + + +Timings for individual modules: + +Sum of individual times ... 17.846 sec (= 0.297 min) +GTO integral calculation ... 0.534 sec (= 0.009 min) 3.0 % +SCF iterations ... 0.410 sec (= 0.007 min) 2.3 % +CASSCF iterations ... 9.882 sec (= 0.165 min) 55.4 % +Multireference CI module ... 7.020 sec (= 0.117 min) 39.3 % + ****ORCA TERMINATED NORMALLY**** +TOTAL RUN TIME: 0 days 0 hours 0 minutes 18 seconds 202 msec diff --git a/arkane/encorr/isodesmic.py b/arkane/encorr/isodesmic.py index 83a9164e175..f4f39d0d85d 100644 --- a/arkane/encorr/isodesmic.py +++ b/arkane/encorr/isodesmic.py @@ -4,7 +4,7 @@ # # # RMG - Reaction Mechanism Generator # # # -# Copyright (c) 2002-2023 Prof. William H. Green (whgreen@mit.edu), # +# Copyright (c) 2002-2024 Prof. William H. Green (whgreen@mit.edu), # # Prof. Richard H. West (r.west@neu.edu) and the RMG Team (rmg_dev@mit.edu) # # # # Permission is hereby granted, free of charge, to any person obtaining a # @@ -41,28 +41,29 @@ https://doi.org/10.1021/jp404158v """ -import signal -from collections import deque +import logging +from copy import deepcopy +from typing import List, Union -from lpsolve55 import lpsolve, EQ, LE import numpy as np - -from rmgpy.molecule import Molecule -from rmgpy.quantity import ScalarQuantity +from scipy.optimize import LinearConstraint, milp from arkane.modelchem import LOT - -# Optional Imports -try: - import pyomo.environ as pyo -except ImportError: - pyo = None +from rmgpy.molecule import Bond, Molecule +from rmgpy.quantity import ScalarQuantity class ErrorCancelingSpecies: """Class for target and known (reference) species participating in an error canceling reaction""" - def __init__(self, molecule, low_level_hf298, level_of_theory, high_level_hf298=None, source=None): + def __init__( + self, + molecule, + low_level_hf298, + level_of_theory, + high_level_hf298=None, + source=None, + ): """ Args: @@ -76,35 +77,45 @@ def __init__(self, molecule, low_level_hf298, level_of_theory, high_level_hf298= if isinstance(molecule, Molecule): self.molecule = molecule else: - raise ValueError(f'ErrorCancelingSpecies molecule attribute must be an rmgpy Molecule object. Instead a ' - f'{type(molecule)} object was given') + raise ValueError( + f"ErrorCancelingSpecies molecule attribute must be an rmgpy Molecule object. Instead a " + f"{type(molecule)} object was given" + ) if isinstance(level_of_theory, LOT): self.level_of_theory = level_of_theory else: - raise ValueError(f'The level of theory used to calculate the low level Hf298 must be provided ' - f'for consistency checks. Instead, a {type(level_of_theory)} object was given') + raise ValueError( + f"The level of theory used to calculate the low level Hf298 must be provided " + f"for consistency checks. Instead, a {type(level_of_theory)} object was given" + ) if not isinstance(low_level_hf298, ScalarQuantity): if isinstance(low_level_hf298, tuple): low_level_hf298 = ScalarQuantity(*low_level_hf298) else: - raise TypeError(f'Low level Hf298 should be a ScalarQuantity object or its tuple representation, but ' - f'received {low_level_hf298} instead.') + raise TypeError( + f"Low level Hf298 should be a ScalarQuantity object or its tuple representation, but " + f"received {low_level_hf298} instead." + ) self.low_level_hf298 = low_level_hf298 # If the species is a reference species, then the high level data is already known - if high_level_hf298 is not None and not isinstance(high_level_hf298, ScalarQuantity): + if high_level_hf298 is not None and not isinstance( + high_level_hf298, ScalarQuantity + ): if isinstance(high_level_hf298, tuple): high_level_hf298 = ScalarQuantity(*high_level_hf298) else: - raise TypeError(f'High level Hf298 should be a ScalarQuantity object or its tuple representation, but ' - f'received {high_level_hf298} instead.') + raise TypeError( + f"High level Hf298 should be a ScalarQuantity object or its tuple representation, but " + f"received {high_level_hf298} instead." + ) self.high_level_hf298 = high_level_hf298 self.source = source def __repr__(self): - return f'' + return f"" class ErrorCancelingReaction: @@ -131,22 +142,24 @@ def __init__(self, target, species): # Perform a consistency check that all species are using the same level of theory for spcs in species.keys(): if spcs.level_of_theory != self.level_of_theory: - raise ValueError(f'Species {spcs} has level of theory {spcs.level_of_theory}, which does not match the ' - f'level of theory of the reaction of {self.level_of_theory}') + raise ValueError( + f"Species {spcs} has level of theory {spcs.level_of_theory}, which does not match the " + f"level of theory of the reaction of {self.level_of_theory}" + ) # Does not include the target, which is handled separately. self.species = species def __repr__(self): - reactant_string = f'1*{self.target.molecule.to_smiles()}' - product_string = '' + reactant_string = f"1*{self.target.molecule.to_smiles()}" + product_string = "" for spcs, coeff in self.species.items(): if coeff > 0: - product_string += f' + {int(coeff)}*{spcs.molecule.to_smiles()}' + product_string += f" + {int(coeff)}*{spcs.molecule.to_smiles()}" else: - reactant_string += f' + {-1*int(coeff)}*{spcs.molecule.to_smiles()}' + reactant_string += f" + {-1*int(coeff)}*{spcs.molecule.to_smiles()}" - return f' {product_string[3:]} >' + return f" {product_string[3:]} >" def calculate_target_thermo(self): """ @@ -155,12 +168,191 @@ def calculate_target_thermo(self): Returns: rmgpy.quantity.ScalarQuantity: Hf298 in 'J/mol' estimated for the target species """ - low_level_h_rxn = sum(spec[0].low_level_hf298.value_si*spec[1] for spec in self.species.items()) - \ - self.target.low_level_hf298.value_si + low_level_h_rxn = ( + sum( + spec.low_level_hf298.value_si * coeff + for spec, coeff in self.species.items() + ) + - self.target.low_level_hf298.value_si + ) + + target_thermo = ( + sum( + spec.high_level_hf298.value_si * coeff + for spec, coeff in self.species.items() + ) + - low_level_h_rxn + ) + return ScalarQuantity(target_thermo, "J/mol") + + +class AtomConstraint: + def __init__(self, label, connections=None): + self.label = label + self.connections = connections if connections is not None else [] + + def __eq__(self, other): + if isinstance(other, AtomConstraint): + if self.label == other.label: + return self.match_connections(other) + + return False + + else: + raise NotImplementedError( + f"AtomConstraint object has no __eq__ defined for other object of type " + f"{type(other)}" + ) + + def match_connections(self, other): + if len(self.connections) != len(other.connections): + return False + + connections = deepcopy(other.connections) + for c in self.connections: + for i, c_other in enumerate(connections): + if c == c_other: + break + else: + return False + connections.pop(i) + + return True + + def __repr__(self): + return f"{self.label}" + "".join([f"({c})" for c in self.connections]) + + +class BondConstraint: + def __init__(self, atom1, atom2, bond_order): + self.atom1 = atom1 + self.atom2 = atom2 + self.bond_order = bond_order + + def __eq__(self, other): + if isinstance(other, BondConstraint): + if self.bond_order == other.bond_order: + if (self.atom1 == other.atom1 and self.atom2 == other.atom2) or ( + self.atom1 == other.atom2 and self.atom2 == other.atom1 + ): + return True + return False + + if isinstance(other, GenericConstraint): + return False + + else: + raise NotImplementedError( + f"BondConstraint object has no __eq__ defined for other object of type " + f"{type(other)}" + ) + + def __repr__(self): + symbols = ["", "-", "=", "#"] + return f"{self.atom1}{symbols[self.bond_order]}{self.atom2}" + + +class Connection: + def __init__(self, atom, bond_order): + self.atom = atom + self.bond_order = bond_order + + def __eq__(self, other): + if isinstance(other, Connection): + if self.bond_order == other.bond_order: + if self.atom == other.atom: + return True + return False + + else: + raise NotImplementedError( + f"Connection object has no __eq__ defined for other object of type {type(other)}" + ) + + def __repr__(self): + symbols = ["", "-", "=", "#"] + return f"{symbols[self.bond_order]}{self.atom}" + + +class GenericConstraint: + def __init__(self, constraint_str: str): + self.constraint_str = constraint_str + + def __eq__(self, other): + if isinstance(other, BondConstraint): + return False + elif isinstance(other, GenericConstraint): + return self.constraint_str == other.constraint_str + else: + raise NotImplementedError( + f"GenericConstraint object has no __eq__ defined for other object of " + f"type {type(other)}" + ) + def __repr__(self): + return self.constraint_str + + +def bond_centric_constraints( + species: ErrorCancelingSpecies, constraint_class: str +) -> List[BondConstraint]: + constraints = [] + contraint_func = CONSTRAINT_CLASSES[constraint_class] + molecule = species.molecule + + for bond in molecule.get_all_edges(): + constraints.append(contraint_func(bond)) + + return constraints + + +def _buerger_rc2(bond: Bond) -> BondConstraint: + atom1 = bond.atom1 + atom2 = bond.atom2 + + if atom1.symbol > atom2.symbol: + atom1, atom2 = atom2, atom1 + + atom1 = AtomConstraint(label=atom1.symbol) + atom2 = AtomConstraint(label=atom2.symbol) + + return BondConstraint(atom1=atom1, atom2=atom2, bond_order=int(bond.order)) + + +def _buerger_rc3(bond: Bond) -> BondConstraint: + atom1 = bond.atom1 + atom2 = bond.atom2 + + if atom1.symbol > atom2.symbol: + atom1, atom2 = atom2, atom1 - target_thermo = sum(spec[0].high_level_hf298.value_si*spec[1] for spec in self.species.items()) - \ - low_level_h_rxn - return ScalarQuantity(target_thermo, 'J/mol') + atom1 = AtomConstraint(label=f"{atom1.symbol}{atom1.connectivity1}") + atom2 = AtomConstraint(label=f"{atom2.symbol}{atom2.connectivity1}") + + return BondConstraint(atom1=atom1, atom2=atom2, bond_order=int(bond.order)) + + +def _buerger_rc4(bond: Bond) -> BondConstraint: + atom1 = bond.atom1 + atom2 = bond.atom2 + + if atom1.symbol > atom2.symbol: + atom1, atom2 = atom2, atom1 + + atoms = [] + + for atom in [atom1, atom2]: + connections = [] + for a, b in atom.bonds.items(): + ac = AtomConstraint(label=f"{a.symbol}{a.connectivity1}") + bond_order = b.order + connections.append(Connection(atom=ac, bond_order=bond_order)) + atoms.append( + AtomConstraint( + label=f"{atom.symbol}{atom.connectivity1}", connections=connections + ) + ) + + return BondConstraint(atom1=atoms[0], atom2=atoms[1], bond_order=int(bond.order)) class SpeciesConstraints: @@ -168,76 +360,181 @@ class SpeciesConstraints: A class for defining and enumerating constraints to ReferenceSpecies objects for error canceling reactions """ - def __init__(self, target, reference_list, conserve_bonds=True, conserve_ring_size=True): + def __init__( + self, + target, + reference_list, + isodesmic_class="rc2", + conserve_ring_size=True, + limit_charges=True, + limit_scope=True, + ): """ Define the constraints that will be enforced, and determine the mapping of indices in the constraint vector to the labels for these constraints. - To reduce the size of the linear programming problem that will try to find error canceling reactions of the - target and subsets of the reference species, the `reference_species` list is automatically pruned to remove - species that have additional atom, bond, and/or ring attributes not found in the target molecule. + Notes: + To reduce the size of the linear programming problem that will try to find error canceling reactions of the + target and subsets of the reference species, the `reference_species` list is automatically pruned to remove + species that have additional atom, bond, and/or ring attributes not found in the target molecule. + + Charge is also explicitly conserved, as there are charged species in the reference database Args: target (ErrorCancelingSpecies): The target species of the error canceling reaction scheme reference_list(list): A list of ErrorCancelingSpecies objects for the reference species that can participate in the error canceling reaction scheme - conserve_bonds (bool, optional): Enforce the number of each bond type be conserved - conserve_ring_size (bool, optional): Enforce that the number of each ring size be conserved + isodesmic_class (str, optional): Reaction classes as defined by Buerger et al. that determine how specific + the constraints are. + conserve_ring_size (bool, optional): Enforce that the number of each ring size be conserved. + limit_charges (bool, optional): Only allow species in the reaction that are within the range [C, 0] for + anions or [0, C] for cations where "C" is the charge of the target + limit_scope (bool, optional): Exclude any molecules from the reference set that have features not contained + in the target molecule. This will reduce the size of the linear programing problem being solved to yield + faster, possibly more accurate results """ self.target = target self.all_reference_species = reference_list self.reference_species = [] - self.conserve_bonds = conserve_bonds + self.isodesmic_class = isodesmic_class self.conserve_ring_size = conserve_ring_size - self.constraint_map = self._get_constraint_map() - - def _get_constraint_map(self): - # Enumerate all of the constraints in the target molecule to initialize the constraint mapping - constraint_map = {label: i for i, label in enumerate(self.target.molecule.get_element_count().keys())} - if self.conserve_bonds: - j = len(constraint_map) - constraint_map.update( - {label: j + i for i, label in enumerate(self.target.molecule.enumerate_bonds().keys())}) + self.limit_charges = limit_charges + self.limit_scope = limit_scope + + def _get_constraint_lists(self): + full_constraints_list = [self._get_all_constraints(self.target)] + for ref_spcs in self.all_reference_species: + full_constraints_list.append(self._get_all_constraints(ref_spcs)) + + return full_constraints_list + + def _get_ring_constraints( + self, species: ErrorCancelingSpecies + ) -> List[GenericConstraint]: + ring_features = [] + rings = species.molecule.get_smallest_set_of_smallest_rings() + for ring in rings: + ring_features.append(GenericConstraint(constraint_str=f"{len(ring)}_ring")) + + return ring_features + + def _get_all_constraints( + self, species: ErrorCancelingSpecies + ) -> List[Union[BondConstraint, GenericConstraint]]: + features = bond_centric_constraints(species, self.isodesmic_class) if self.conserve_ring_size: - j = len(constraint_map) - possible_rings_sizes = set(f'{len(sssr)}_ring' for sssr in - self.target.molecule.get_smallest_set_of_smallest_rings()) - constraint_map.update({label: j + i for i, label in enumerate(possible_rings_sizes)}) + features += self._get_ring_constraints(species) - return constraint_map + features.sort(key=lambda x: x.__repr__()) + return features - def _enumerate_constraints(self, species): + def _enumerate_constraints(self, full_constraints_list): + """ + Return the target constraint counts and the reference constraint counts. """ - Determine the constraint vector for a molecule given the enforced constraints - Args: - species (ErrorCancelingSpecies): Species whose constraints are to be enumerated + target_constraints = full_constraints_list[0] + reference_constraintss = full_constraints_list[1:] - Returns: - np.ndarray: vector of the number of instances of each constraining feature e.g. number of carbon atoms - """ - constraint_vector = np.zeros(len(self.constraint_map)) - molecule = species.molecule + # Enumerate through the constraints of reference species and keep only those that are present in the target + enumerated_reference_constraintss = [] + + if self.limit_scope: + + for i, ref_spc in enumerate(self.all_reference_species): - try: - atoms = molecule.get_element_count() - for atom_label, count in atoms.items(): - constraint_vector[self.constraint_map[atom_label]] += count + if not all(constraint in target_constraints for constraint in reference_constraintss[i]): + continue - if self.conserve_bonds: - bonds = molecule.enumerate_bonds() - for bond_label, count in bonds.items(): - constraint_vector[self.constraint_map[bond_label]] += count + self.reference_species.append(ref_spc) + enumerated_reference_constraintss.append(reference_constraintss[i]) - if self.conserve_ring_size: - rings = molecule.get_smallest_set_of_smallest_rings() - for ring in rings: - constraint_vector[self.constraint_map[f'{len(ring)}_ring']] += 1 - except KeyError: # This molecule has a feature not found in the target molecule. Return None to exclude this - return None + else: + self.reference_species = self.all_reference_species + enumerated_reference_constraintss = reference_constraintss + + # Get a list of the unique constraints sorted by their string representation + if self.limit_scope: - return constraint_vector + # The scope of constraints to consider is the target constraints + unique_constraints = self._get_unique_constraints(target_constraints) + unique_constraints.sort(key=lambda x: x.__repr__()) + + else: + all_constraints = target_constraints + [constraint for constraints in enumerated_reference_constraintss for constraint in constraints] + unique_constraints = self._get_unique_constraints(all_constraints) + unique_constraints.sort(key=lambda x: x.__repr__()) + + # Get the counts of each unique constraint in the target and reference constraints + target_constraint_counts = [target_constraints.count(c) for c in unique_constraints] + reference_constraint_counts = [] + + for i, ref_constraints in enumerate(enumerated_reference_constraintss): + reference_constraint_counts.append([ref_constraints.count(c) for c in unique_constraints]) + + return target_constraint_counts, reference_constraint_counts + + def _get_unique_constraints(self, constraints): + # Constraints are unhashable, so we need to use some workaround to get unique constraints + constraint_dict = {constraint.__repr__(): constraint for constraint in constraints} + return list(constraint_dict.values()) + + def _enumerate_charge_constraints(self, target_constraints, reference_constraints): + charge = self.target.molecule.get_net_charge() + target_constraints.append(charge) + + for i, spcs in enumerate(self.reference_species): + reference_constraints[i].append(spcs.molecule.get_net_charge()) + + if self.limit_charges: + allowed_reference_species = [] + new_constraints = [] + + if charge < 0: + allowable_charges = list(range(charge, 0)) + elif charge > 0: + allowable_charges = list(range(1, charge + 1)) + else: + allowable_charges = [0] + for i, spcs in enumerate(self.reference_species): + if reference_constraints[i][-1] in allowable_charges: + allowed_reference_species.append(spcs) + new_constraints.append(reference_constraints[i]) + + reference_constraints = new_constraints + self.reference_species = allowed_reference_species + + return target_constraints, reference_constraints + + def _enumerate_element_constraints(self, target_constraints, reference_constraints): + all_elements = set() + for spc in self.reference_species: + all_elements.update(spc.molecule.get_element_count().keys()) + + # Check that the target and reference species have the same elements to be able to satisfy mass conservation + if set(self.target.molecule.get_element_count().keys()) != all_elements: + logging.warning( + "Target species and reference species do not have the same elements:" + f"Target: {' '.join(self.target.molecule.get_element_count().keys())}" + f"Reference: {all_elements}" + ) + + all_elements.update(self.target.molecule.get_element_count().keys()) + all_elements = sorted(list(all_elements)) + + element_count = self.target.molecule.get_element_count() + new_constraints = [element_count.get(element, 0) for element in all_elements] + target_constraints.extend(new_constraints) + + for i, spc in enumerate(self.reference_species): + element_count = spc.molecule.get_element_count() + new_constraints = [ + element_count.get(element, 0) for element in all_elements + ] + reference_constraints[i].extend(new_constraints) + + return target_constraints, reference_constraints def calculate_constraints(self): """ @@ -248,16 +545,41 @@ def calculate_constraints(self): - target constraint vector (1 x len(constraints)) - constraint matrix for allowable reference species (len(self.reference_species) x len(constraints)) """ - target_constraints = self._enumerate_constraints(self.target) - constraint_matrix = [] - - for spcs in self.all_reference_species: - spcs_constraints = self._enumerate_constraints(spcs) - if spcs_constraints is not None: # This species is allowed - self.reference_species.append(spcs) - constraint_matrix.append(spcs_constraints) + full_constraint_list = self._get_constraint_lists() + target_constraints, reference_constraints = self._enumerate_constraints( + full_constraint_list + ) + target_constraints, reference_constraints = self._enumerate_charge_constraints( + target_constraints, reference_constraints + ) + target_constraints, reference_constraints = self._enumerate_element_constraints( + target_constraints, reference_constraints + ) + + target_constraints = np.array(target_constraints, dtype=int) + constraint_matrix = np.array(reference_constraints, dtype=int) + + return target_constraints, constraint_matrix + + +def _clean_up_constraints(target_constraints, constraint_matrix): + # make sure that the constraint matrix is 2d + if len(constraint_matrix.shape) == 1: + constraint_matrix = np.array([constraint_matrix], dtype=int) + + # Remove any columns that are all zeros + zero_indices = np.where(~constraint_matrix.any(axis=0))[0] + # Check that this wouldn't eliminate a non-zero target entry + for z in zero_indices: + if ( + target_constraints[z] != 0 + ): # This problem is not solvable. Return None to signal this + return None, None + indices = [i for i in range(constraint_matrix.shape[1]) if i not in zero_indices] + constraint_matrix = np.take(constraint_matrix, indices=indices, axis=1) + target_constraints = np.take(target_constraints, indices=indices) - return target_constraints, np.array(constraint_matrix, dtype=int) + return target_constraints, constraint_matrix class ErrorCancelingScheme: @@ -265,27 +587,42 @@ class ErrorCancelingScheme: A Base class for calculating target species thermochemistry using error canceling reactions """ - def __init__(self, target, reference_set, conserve_bonds, conserve_ring_size): + def __init__( + self, + target, + reference_set, + isodesmic_class, + conserve_ring_size, + limit_charges, + limit_scope, + ): """ Args: target (ErrorCancelingSpecies): Species whose Hf298 will be calculated using error canceling reactions reference_set (list): list of reference species (as ErrorCancelingSpecies) that can participate in error canceling reactions to help calculate the thermochemistry of the target - conserve_bonds (bool): Flag to determine if the number and type of each bond must be conserved in each error - canceling reaction conserve_ring_size (bool): Flag to determine if the number of each ring size must be conserved in each error canceling reaction """ self.target = target - self.constraints = SpeciesConstraints(target, reference_set, conserve_bonds=conserve_bonds, - conserve_ring_size=conserve_ring_size) - - self.target_constraint, self.constraint_matrix = self.constraints.calculate_constraints() + self.constraints = SpeciesConstraints( + target, + reference_set, + isodesmic_class=isodesmic_class, + conserve_ring_size=conserve_ring_size, + limit_charges=limit_charges, + limit_scope=limit_scope, + ) + + ( + self.target_constraint, + self.constraint_matrix, + ) = self.constraints.calculate_constraints() self.reference_species = self.constraints.reference_species - def _find_error_canceling_reaction(self, reference_subset, milp_software=None): + def _find_error_canceling_reaction(self, reference_subset): """ Automatically find a valid error canceling reaction given a subset of the available benchmark species. This is done by solving a mixed integer linear programming (MILP) problem similiar to @@ -293,107 +630,45 @@ def _find_error_canceling_reaction(self, reference_subset, milp_software=None): Args: reference_subset (list): A list of indices from self.reference_species that can participate in the reaction - milp_software (list, optional): Solvers to try in order. Defaults to ['lpsolve'] or if pyomo is available - defaults to ['lpsolve', 'pyomo']. lpsolve is usually faster. Returns: tuple(ErrorCancelingReaction, np.ndarray) - Reaction with the target species (if a valid reaction is found, else ``None``) - Indices (of the subset) for the species that participated in the return reaction """ - if milp_software is None: - milp_software = ['lpsolve'] - if pyo is not None: - milp_software.append('pyomo') # Define the constraints based on the provided subset c_matrix = np.take(self.constraint_matrix, reference_subset, axis=0) + + # Remove unnecessary constraints + target_constraint, c_matrix = _clean_up_constraints( + self.target_constraint, c_matrix + ) + if target_constraint is None or c_matrix is None: # The problem is not solvable + return None, None + + # Setup MILP problem c_matrix = np.tile(c_matrix, (2, 1)) sum_constraints = np.sum(c_matrix, 1, dtype=int) - targets = -1*self.target_constraint + targets = -1 * target_constraint m = c_matrix.shape[0] n = c_matrix.shape[1] - split = int(m/2) - - for solver in milp_software: - if solver == 'pyomo': - # Check that pyomo is available - if pyo is None: - raise ImportError('Cannot import optional package pyomo. Either install this dependency with ' - '`conda install -c conda-forge pyomo glpk` or set milp_software to `lpsolve`') - - # Setup the MILP problem using pyomo - lp_model = pyo.ConcreteModel() - lp_model.i = pyo.RangeSet(0, m - 1) - lp_model.j = pyo.RangeSet(0, n - 1) - lp_model.r = pyo.RangeSet(0, split-1) # indices before the split correspond to reactants - lp_model.p = pyo.RangeSet(split, m - 1) # indices after the split correspond to products - lp_model.v = pyo.Var(lp_model.i, domain=pyo.NonNegativeIntegers) # The stoich. coef. we are solving for - lp_model.c = pyo.Param(lp_model.i, lp_model.j, initialize=lambda _, i_ind, j_ind: c_matrix[i_ind, - j_ind]) - lp_model.s = pyo.Param(lp_model.i, initialize=lambda _, i_ind: sum_constraints[i_ind]) - lp_model.t = pyo.Param(lp_model.j, initialize=lambda _, j_ind: targets[j_ind]) - - lp_model.obj = pyo.Objective(rule=_pyo_obj_expression) - lp_model.constraints = pyo.Constraint(lp_model.j, rule=_pyo_constraint_rule) - - # Solve the MILP problem using the GLPK MILP solver (https://www.gnu.org/software/glpk/) - opt = pyo.SolverFactory('glpk') - results = opt.solve(lp_model, timelimit=1) - - # Return the solution if a valid reaction is found. Otherwise continue to next solver - if results.solver.termination_condition == pyo.TerminationCondition.optimal: - # Extract the solution and find the species with non-zero stoichiometric coefficients - solution = lp_model.v.extract_values().values() - break + split = int(m / 2) - elif solver == 'lpsolve': - # Save the current signal handler - sig = signal.getsignal(signal.SIGINT) - - # Setup the MILP problem using lpsolve - lp = lpsolve('make_lp', 0, m) - lpsolve('set_verbose', lp, 2) # Reduce the logging from lpsolve - lpsolve('set_obj_fn', lp, sum_constraints) - lpsolve('set_minim', lp) - - for j in range(n): - lpsolve('add_constraint', lp, np.concatenate((c_matrix[:split, j], -1*c_matrix[split:, j])), EQ, - targets[j]) - - lpsolve('add_constraint', lp, np.ones(m), LE, 20) # Use at most 20 species (including replicates) - lpsolve('set_timeout', lp, 1) # Move on if lpsolve can't find a solution quickly - - # Constrain v_i to be 4 or less - for i in range(m): - lpsolve('set_upbo', lp, i, 4) - - # All v_i must be integers - lpsolve('set_int', lp, [True]*m) - - status = lpsolve('solve', lp) - - # Reset signal handling since lpsolve changed it - try: - signal.signal(signal.SIGINT, sig) - except ValueError: - # This is not being run in the main thread, so we cannot reset signal - pass - except TypeError: - print( - "Failed to reset signal handling in LPSolve - are you running pytest?" - ) + constraints = tuple((LinearConstraint(A=np.concatenate((c_matrix[:split, j], -1 * c_matrix[split:, j])), lb=targets[j], ub=targets[j]) for j in range(n))) - # Return the solution if a valid reaction is found. Otherwise continue to next solver - if status == 0: - _, solution = lpsolve('get_solution', lp)[:2] - break + result = milp( + sum_constraints, + integrality=1, + constraints=constraints, + options={"time_limit": 10}, + ) - else: - raise ValueError(f'Unrecognized MILP solver {solver} for isodesmic reaction generation') - - else: + if result.status != 0: + logging.warning("Optimization could not find a valid solution.") return None, None + + solution = result.x reaction = ErrorCancelingReaction(self.target, dict()) subset_indices = [] @@ -401,13 +676,19 @@ def _find_error_canceling_reaction(self, reference_subset, milp_software=None): if v > 0: subset_indices.append(index % split) if index < split: - reaction.species.update({self.reference_species[reference_subset[index]]: -v}) + reaction.species.update( + {self.reference_species[reference_subset[index]]: -v} + ) else: - reaction.species.update({self.reference_species[reference_subset[index % split]]: v}) + reaction.species.update( + {self.reference_species[reference_subset[index % split]]: v} + ) return reaction, np.array(subset_indices) - def multiple_error_canceling_reaction_search(self, n_reactions_max=20, milp_software=None): + def multiple_error_canceling_reaction_search( + self, n_reactions_max=20, + ): """ Generate multiple error canceling reactions involving the target and a subset of the reference species. @@ -418,21 +699,20 @@ def multiple_error_canceling_reaction_search(self, n_reactions_max=20, milp_soft Args: n_reactions_max (int, optional): The maximum number of found reactions that will be returned, after which no further searching will occur even if there are possible subsets left in the queue. - milp_software (list, optional): Solvers to try in order. Defaults to ['lpsolve'] or if pyomo is available - defaults to ['lpsolve', 'pyomo']. lpsolve is usually faster. Returns: list: A list of the found error canceling reactions """ - subset_queue = deque() - subset_queue.append(np.arange(0, len(self.reference_species))) + subset_queue = [np.arange(0, len(self.reference_species))] reaction_list = [] while (len(subset_queue) != 0) and (len(reaction_list) < n_reactions_max): - subset = subset_queue.popleft() + subset = subset_queue.pop() if len(subset) == 0: continue - reaction, subset_indices = self._find_error_canceling_reaction(subset, milp_software=milp_software) + reaction, subset_indices = self._find_error_canceling_reaction( + subset + ) if reaction is None: continue else: @@ -441,9 +721,18 @@ def multiple_error_canceling_reaction_search(self, n_reactions_max=20, milp_soft for index in subset_indices: subset_queue.append(np.delete(subset, index)) + # Clean up the queue to remove subsets that would allow for already found reactions + new_queue = [] + reaction_indices = [subset[i] for i in subset_indices] + for s in subset_queue: + if not all([i in s for i in reaction_indices]): + new_queue.append(s) + + subset_queue = new_queue + return reaction_list - def calculate_target_enthalpy(self, n_reactions_max=20, milp_software=None): + def calculate_target_enthalpy(self, n_reactions_max=5): """ Perform a multiple error canceling reactions search and calculate hf298 for the target species by taking the median hf298 value from among the error canceling reactions found @@ -451,47 +740,57 @@ def calculate_target_enthalpy(self, n_reactions_max=20, milp_software=None): Args: n_reactions_max (int, optional): The maximum number of found reactions that will returned, after which no further searching will occur even if there are possible subsets left in the queue. - milp_software (list, optional): Solvers to try in order. Defaults to ['lpsolve'] or if pyomo is available - defaults to ['lpsolve', 'pyomo']. lpsolve is usually faster. Returns: tuple(ScalarQuantity, list) - Standard heat of formation at 298 K calculated for the target species - reaction list containing all error canceling reactions found """ - reaction_list = self.multiple_error_canceling_reaction_search(n_reactions_max, milp_software) + reaction_list = self.multiple_error_canceling_reaction_search(n_reactions_max) + if len(reaction_list) == 0: # No reactions found + return None, reaction_list h298_list = np.zeros(len(reaction_list)) for i, rxn in enumerate(reaction_list): h298_list[i] = rxn.calculate_target_thermo().value_si - return ScalarQuantity(np.median(h298_list), 'J/mol'), reaction_list - - -def _pyo_obj_expression(model): - return pyo.summation(model.v, model.s, index=model.i) - - -def _pyo_constraint_rule(model, col): - return sum(model.v[row] * model.c[row, col] for row in model.r) - \ - sum(model.v[row] * model.c[row, col] for row in model.p) == model.t[col] + return ScalarQuantity(np.median(h298_list), "J/mol"), reaction_list class IsodesmicScheme(ErrorCancelingScheme): """ An error canceling reaction where the number and type of both atoms and bonds are conserved """ + def __init__(self, target, reference_set): - super().__init__(target, reference_set, conserve_bonds=True, conserve_ring_size=False) + super().__init__( + target, + reference_set, + isodesmic_class="rc2", + conserve_ring_size=False, + limit_charges=True, + limit_scope=True, + ) class IsodesmicRingScheme(ErrorCancelingScheme): """ A stricter form of the traditional isodesmic reaction scheme where the number of each ring size is also conserved """ + def __init__(self, target, reference_set): - super().__init__(target, reference_set, conserve_bonds=True, conserve_ring_size=True) + super().__init__( + target, + reference_set, + isodesmic_class="rc2", + conserve_ring_size=True, + limit_charges=True, + limit_scope=True, + ) + + +CONSTRAINT_CLASSES = {"rc2": _buerger_rc2, "rc3": _buerger_rc3, "rc4": _buerger_rc4} -if __name__ == '__main__': +if __name__ == "__main__": pass diff --git a/arkane/ess/adapter.py b/arkane/ess/adapter.py index 9b25b1e2c1c..f0a30c932a2 100644 --- a/arkane/ess/adapter.py +++ b/arkane/ess/adapter.py @@ -45,10 +45,11 @@ class ESSAdapter(ABC): An abstract ESS Adapter class """ - def __init__(self, path, check_for_errors=True): + def __init__(self, path, check_for_errors=True, scratch_directory=None): self.path = path if check_for_errors: self.check_for_errors() + self.scratch_directory = scratch_directory if scratch_directory is not None else os.path.join(os.path.abspath('.'), str('scratch')) @abstractmethod def check_for_errors(self): @@ -174,7 +175,7 @@ def get_symmetry_properties(self): coordinates, atom_numbers, _ = self.load_geometry() unique_id = '0' # Just some name that the SYMMETRY code gives to one of its jobs # Scratch directory that the SYMMETRY code writes its files in: - scr_dir = os.path.join(os.path.abspath('.'), str('scratch')) + scr_dir = self.scratch_directory if not os.path.exists(scr_dir): os.makedirs(scr_dir) try: diff --git a/arkane/ess/factory.py b/arkane/ess/factory.py index 8e503900938..ca0baee26c6 100644 --- a/arkane/ess/factory.py +++ b/arkane/ess/factory.py @@ -61,6 +61,7 @@ def register_ess_adapter(ess: str, def ess_factory(fullpath: str, check_for_errors: bool = True, + scratch_directory: str = None, ) -> Type[ESSAdapter]: """ A factory generating the ESS adapter corresponding to ``ess_adapter``. @@ -106,4 +107,4 @@ def ess_factory(fullpath: str, raise InputError(f'The file at {fullpath} could not be identified as a ' f'Gaussian, Molpro, Orca, Psi4, QChem, or TeraChem log file.') - return _registered_ess_adapters[ess_name](path=fullpath, check_for_errors=check_for_errors) + return _registered_ess_adapters[ess_name](path=fullpath, check_for_errors=check_for_errors, scratch_directory=scratch_directory) diff --git a/arkane/ess/molpro.py b/arkane/ess/molpro.py index 35fcfbb00c0..01be43adac9 100644 --- a/arkane/ess/molpro.py +++ b/arkane/ess/molpro.py @@ -371,6 +371,9 @@ def load_energy(self, zpe_scale_factor=1.): if 'CCSD' in line and 'energy=' in line: e_elect = float(line.split()[-1]) break + if 'RS2C' in line and 'energy' in line: + e_elect = float(line.split()[-1]) + break if e_elect is None and mrci: # No Davidson correction is given, search for MRCI energy read_e_elect = False @@ -435,10 +438,10 @@ def load_negative_frequency(self): if len(freqs) == 1: return -float(freqs[0]) elif len(freqs) > 1: - logging.info('More than one imaginary frequency in Molpro output file {0}.'.format(self.path)) + logging.info(f'More than one imaginary frequency in Molpro output file {self.path}.') return -float(freqs[0]) else: - raise LogError('Unable to find imaginary frequency in Molpro output file {0}'.format(self.path)) + raise LogError(f'Unable to find imaginary frequency in Molpro output file {self.path}') def load_scan_energies(self): """ @@ -458,7 +461,8 @@ def get_T1_diagnostic(self): if 'T1 diagnostic: ' in line: items = line.split() return float(items[-1]) - raise LogError('Unable to find T1 diagnostic in energy file: {0}'.format(self.path)) + logging.warning(f'Unable to find T1 diagnostic in energy file: {self.path}') + return None def get_D1_diagnostic(self): """ @@ -472,7 +476,8 @@ def get_D1_diagnostic(self): if 'D1 diagnostic: ' in line: items = line.split() return float(items[-1]) - raise LogError('Unable to find D1 diagnostic in energy file: {0}'.format(self.path)) + logging.warning(f'Unable to find D1 diagnostic in energy file: {self.path}') + return None def load_scan_pivot_atoms(self): """Not implemented for Molpro""" diff --git a/arkane/ess/orca.py b/arkane/ess/orca.py index e80d5ad180b..a538f4dea31 100644 --- a/arkane/ess/orca.py +++ b/arkane/ess/orca.py @@ -109,6 +109,9 @@ def load_force_constant_matrix(self): The units of the returned force constants are J/m^2. If no force constant matrix can be found, ``None`` is returned. """ + n_atoms = self.get_number_of_atoms() + if n_atoms == 1: + return None hess_files = list() for (_, _, files) in os.walk(os.path.dirname(self.path)): for file_ in files: @@ -196,7 +199,7 @@ def load_conformer(self, symmetry=None, spin_multiplicity=0, optical_isomers=Non you can use the `symmetry` parameter to substitute your own value; if not provided, the value in the Orca log file will be adopted. """ - freq, mmass, rot, unscaled_frequencies = [], [], [], [] + freq, mmass, rot, unscaled_frequencies, modes = [], [], [], [], [] e0 = 0.0 if optical_isomers is None or symmetry is None: @@ -230,7 +233,7 @@ def load_conformer(self, symmetry=None, spin_multiplicity=0, optical_isomers=Non symbols = [symbol_by_number[i] for i in number] inertia = get_principal_moments_of_inertia(coord, numbers=number, symbols=symbols) inertia = list(inertia[0]) - if len(inertia): + if len(inertia) and not all(i == 0.0 for i in inertia): if any(i == 0.0 for i in inertia): inertia.remove(0.0) rot.append(LinearRotor(inertia=(inertia, "amu*angstrom^2"), symmetry=symmetry)) @@ -241,7 +244,11 @@ def load_conformer(self, symmetry=None, spin_multiplicity=0, optical_isomers=Non mmass.append(translation) # Take only the last modes found (in the event of multiple jobs). - modes = [mmass[-1], rot[-1], freq[-1]] + modes = [mmass[-1]] + if len(rot): + modes.append(rot[-1]) + if len(freq): + modes.append(freq[-1]) return Conformer(E0=(e0 * 0.001, "kJ/mol"), modes=modes, diff --git a/arkane/statmech.py b/arkane/statmech.py index 4380bad121f..fdfac0b5122 100644 --- a/arkane/statmech.py +++ b/arkane/statmech.py @@ -37,6 +37,7 @@ import math import os import pathlib +import warnings import matplotlib.pyplot as plt import numpy as np @@ -1132,8 +1133,8 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_ logging.debug('Frequencies from internal Hessian') for i in range(3 * n_atoms - external): - with np.warnings.catch_warnings(): - np.warnings.filterwarnings('ignore', r'invalid value encountered in sqrt') + with warnings.catch_warnings(): + warnings.filterwarnings('ignore', r'invalid value encountered in sqrt') logging.debug(np.sqrt(eig[i]) / (2 * math.pi * constants.c * 100)) # Now we can start thinking about projecting out the internal rotations @@ -1243,8 +1244,8 @@ def project_rotors(conformer, hessian, rotors, linear, is_ts, get_projected_out_ logging.debug('Frequencies from projected Hessian') for i in range(3 * n_atoms): - with np.warnings.catch_warnings(): - np.warnings.filterwarnings('ignore', r'invalid value encountered in sqrt') + with warnings.catch_warnings(): + warnings.filterwarnings('ignore', r'invalid value encountered in sqrt') logging.debug(np.sqrt(eig[i]) / (2 * math.pi * constants.c * 100)) return np.sqrt(eig[-n_vib:]) / (2 * math.pi * constants.c * 100) diff --git a/documentation/source/users/rmg/input.rst b/documentation/source/users/rmg/input.rst index 38dead45171..fbbff659405 100644 --- a/documentation/source/users/rmg/input.rst +++ b/documentation/source/users/rmg/input.rst @@ -59,8 +59,8 @@ ThermoLibrary field must be with respect to the :file:`$RMG/RMG-database/input/t directory. .. note:: - Checks during the initialization are maid to avoid users to use "liquid thermo librairies" in gas phase simulations or to use - "liquid phase libraries" obtained in another solvent that the one defined in the input file in liquid phase simulations. + Checks during the initialization are made to avoid users using "liquid thermo libraries" in gas phase simulations or using + "liquid phase libraries" obtained in another solvent than the one defined in the input file in liquid phase simulations. .. _reactionlibraries: diff --git a/documentation/source/users/rmg/installation/anacondaDeveloper.rst b/documentation/source/users/rmg/installation/anacondaDeveloper.rst index 9eb74479fd9..1926e75f2c4 100644 --- a/documentation/source/users/rmg/installation/anacondaDeveloper.rst +++ b/documentation/source/users/rmg/installation/anacondaDeveloper.rst @@ -4,26 +4,12 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux and Mac OSX ******************************************************************************************* -#. Install the `conda` package manager, if you do not already have it (or Anaconda). - Select one of the following options: +#. Install the `conda` package manager via `miniforge`, if you do not already have it (or Anaconda), by following the `Miniforge installation instructions `_. - a. Users of Fedora Linux and Red Hat derivatives (RHEL, CentOS Stream) may install from the official repositories and EPEL, respectively, with the command :: +#. If your `conda` version is older than 23.10.0, switch the solver backend to `libmamba` :: - sudo dnf install conda - - b. All other users, download and install `Miniconda `_. - - The download will be a .sh file with a name like ``Miniconda3-latest-Linux-x86_64.sh``. - Open a terminal in the same directory as this file, and type the following to install Conda - (replace the name of your .sh file below). :: - - bash Miniconda3-latest-Linux-x86_64.sh - - **When prompted to append Anaconda to your PATH, select or type Yes**. - Install the Conda folder inside your home directory - (typically ``/home/YourUsername/`` in Linux and ``/Users/YourUsername`` in Mac). - - Note that you should reinitialize or restart your terminal in order for the changes to take effect, as the installer will tell you. + conda install -n base conda-libmamba-solver + conda config --set solver libmamba #. There are a few system-level dependencies which are required and should not be installed via Conda. These include `Git `_ for version control, `GNU Make `_, and the C and C++ compilers from the `GNU Compiler Collection (GCC) `_ for compiling RMG. @@ -71,11 +57,6 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux For information on using ``ssh`` with GitHub see the `Connecting to GitHub with SSH `_ -#. Switch the conda solver backend to speed up creation of the RMG environment :: - - conda install -n base conda-libmamba-solver - conda config --set solver libmamba - #. Navigate to the RMG-Py directory :: cd RMG-Py @@ -110,16 +91,11 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux conda activate rmg_env -#. Switch the conda solver to libmamba again, to accelerate any changes you might make to this conda environment in the future:: - - conda config --set solver libmamba - #. Compile RMG-Py after activating the conda environment :: make #. Modify environment variables. Add RMG-Py to the PYTHONPATH to ensure that you can access RMG modules from any folder. - *This is important before the next step in which julia dependencies are installed.* Also, add your RMG-Py folder to PATH to launch ``rmg.py`` from any folder. In general, these commands should be placed in the appropriate shell initialization file. @@ -134,16 +110,17 @@ Installation by Source Using Anaconda Environment for Unix-based Systems: Linux Be sure to either close and reopen your terminal to refresh your environment variables (``source ~/.bashrc`` or ``source ~/.zshrc``). -#. Install and Link Julia dependencies: :: +#. **Optional (Recommended)**: Install and Link Julia dependencies. Ensure that you have modified your environment variables as described above, and then run the following: :: + conda install conda-forge::pyjuliacall + julia -e 'using Pkg; Pkg.add("PyCall");Pkg.build("PyCall");Pkg.add(PackageSpec(name="ReactionMechanismSimulator",rev="main")); using ReactionMechanismSimulator;' - python -c "import julia; julia.install(); import diffeqpy; diffeqpy.install()" - + Installing these dependencies will allow using ``method='ode'`` when solving the Master Equation with Arkane and using ``ReactionMechanismSimulator.jl``-based reactors in RMG. #. Finally, you can run RMG from any location by typing the following (given that you have prepared the input file as ``input.py`` in the current folder). :: - python replace/with/path/to/rmg.py input.py + python rmg.py input.py 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 1c42d02e283..286f77db143 100644 --- a/documentation/source/users/rmg/installation/anacondaUser.rst +++ b/documentation/source/users/rmg/installation/anacondaUser.rst @@ -5,26 +5,7 @@ Binary Installation Using Anaconda for Unix-Based Systems: Linux and Mac OSX **************************************************************************** -#. Install the `conda` package manager, if you do not already have it (or Anaconda). - Select one of the following options: - - a. Users of Fedora Linux and Red Hat derivatives (RHEL, CentOS Stream) may install from the official repositories and EPEL, respectively, with the command :: - - sudo dnf install conda - - b. All other users, download and install `Miniconda `_. - - The download will be a .sh file with a name like ``Miniconda3-latest-Linux-x86_64.sh``. - Open a terminal in the same directory as this file, and type the following to install Conda - (replace the name of your .sh file below). :: - - bash Miniconda3-latest-Linux-x86_64.sh - - **When prompted to append Anaconda to your PATH, select or type Yes**. - Install the Conda folder inside your home directory - (typically ``/home/YourUsername/`` in Linux and ``/Users/YourUsername`` in Mac). - - Note that you should reinitialize or restart your terminal in order for the changes to take effect, as the installer will tell you. +#. Install the `conda` package manager via `miniforge`, if you do not already have it (or Anaconda), by following the `Miniforge installation instructions `_. #. Install both RMG and the RMG-database binaries through the terminal. Dependencies will be installed automatically. It is safest to make a new conda environment for RMG and its dependencies. Type the following command into the terminal to create the new environment named 'rmg_env' containing the latest stable version of the RMG program and its database. :: diff --git a/documentation/source/users/rmg/installation/dependencies.rst b/documentation/source/users/rmg/installation/dependencies.rst index ba8a3ad66c7..f91b2807612 100644 --- a/documentation/source/users/rmg/installation/dependencies.rst +++ b/documentation/source/users/rmg/installation/dependencies.rst @@ -24,7 +24,6 @@ Briefly, RMG depends on the following packages, almost all of which can be found * **graphviz:** generating flux diagrams * **jinja2:** Python templating language for html rendering * **jupyter:** (optional) for using IPython notebooks -* **lpsolve:** mixed integer linear programming solver, used for resonance structure generation. Must also install Python extension. * **markupsafe:** implements XML/HTML/XHTML markup safe strings for Python * **matplotlib:** library for making plots * **mock:** for unit-testing diff --git a/environment.yml b/environment.yml index fd792336c8c..9080a480e5e 100644 --- a/environment.yml +++ b/environment.yml @@ -18,65 +18,66 @@ # moved diffeqpy to pip and (temporarily) removed chemprop # - May 14, 2024 Removed diffeqpy by switching to call SciMLBase and Sundials using JuliaCall # - March 15, 2024 - started migration to Python 3.9 -# +# - Mar 11, 2024 Removed Julia dependencies, now considered optional +# - April 17, 2024 Limit versions of cclib at advice of maintainers. +# - August 4, 2024 Restricted pyrms to <2 name: rmg_env channels: - - defaults - - rmg - conda-forge - cantera + - rmg dependencies: # System-level dependencies - we could install these at the OS level # but by installing them in the conda environment we get better control - - cairo - - cairocffi - - ffmpeg - - xlrd - - xlwt - - h5py - - graphviz - - markupsafe - - psutil + - conda-forge::cairo + - conda-forge::cairocffi + - conda-forge::ffmpeg >= 7 + - conda-forge::xlrd + - conda-forge::xlwt + - conda-forge::h5py + - conda-forge::graphviz >=12 + - conda-forge::markupsafe + - conda-forge::psutil # conda-forge not default, since default has a version information bug # (see https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2421) - conda-forge::ncurses - conda-forge::suitesparse # external software tools for chemistry - - coolprop - - cantera::cantera >=3 + - conda-forge::coolprop + - cantera::cantera =2.6 - conda-forge::mopac - - conda-forge::cclib >=1.6.3 + # see https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2639#issuecomment-2050292972 + - conda-forge::cclib >=1.6.3,<1.9 - conda-forge::openbabel >= 3 - conda-forge::rdkit >=2022.09.1 -# general-purpose external software tools - - conda-forge::pyjuliacall # for calling julia packages - # Python tools - - python >=3.9 # leave as GEQ so that GitHub actions can add EQ w/o breaking (contradictory deps) - - coverage - - cython >=0.25.2 - - scikit-learn - - scipy - - numpy >=1.10.0 - - pydot - - jinja2 - - jupyter - - pymongo - - pyparsing - - pyyaml - - networkx - - pytest - - pytest-cov + - conda-forge::python >=3.9 # leave as GEQ so that GitHub actions can add EQ w/o breaking (contradictory deps) + - conda-forge::coverage + - conda-forge::cython >=0.25.2 + - conda-forge::scikit-learn + - conda-forge::scipy >=1.9 + - conda-forge::numpy >=1.10.0 + - conda-forge::pydot + - conda-forge::jinja2 + - conda-forge::jupyter + - conda-forge::pymongo + - conda-forge::pyparsing + - conda-forge::pyyaml + - conda-forge::networkx + - conda-forge::pytest + - conda-forge::pytest-cov - conda-forge::pytest-check - - matplotlib >=1.5 - - mpmath - - pandas + - conda-forge::pyutilib + - conda-forge::matplotlib >=1.5 + - conda-forge::mpmath + - conda-forge::pandas - conda-forge::gprof2dot - conda-forge::numdifftools - - conda-forge::quantities - - conda-forge::lpsolve55 + # bug in quantities, see: + # https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2694#issuecomment-2489286263 + - conda-forge::quantities !=0.16.0,!=0.16.1 - conda-forge::ringdecomposerlib-python # packages we maintain @@ -84,8 +85,12 @@ dependencies: - rmg::pydqed >=1.0.3 - rmg::symmetry -# conda mutex metapackage - - nomkl +# configure packages to use OpenBLAS instead of Intel MKL + - blas=*=openblas + +# optional dependencies for using ReactionMechanismSimulator +# remove the leading '#' to install the required dependencies + # - conda-forge::pyjuliacall # additional packages that are required, but not specified here (and why) # pydqed, pydas, mopac, and likely others require a fortran compiler (specifically gfortran) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index f8818c78acb..c4084147d5c 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -3379,8 +3379,8 @@ def make_bm_rules_from_template_rxn_map(self, template_rxn_map, nprocs=1, Tref=1 entries = list(self.groups.entries.values()) rxnlists = [(template_rxn_map[entry.label], entry.label) if entry.label in template_rxn_map.keys() else [] for entry in entries] - inputs = np.array([(self.forward_recipe.actions, rxns, Tref, fmax, label, [r.rank for r in rxns]) - for rxns, label in rxnlists]) + inputs = [(self.forward_recipe.actions, rxns, Tref, fmax, label, [r.rank for r in rxns]) + for rxns, label in rxnlists] inds = np.arange(len(inputs)) np.random.shuffle(inds) # want to parallelize in random order @@ -3389,9 +3389,9 @@ def make_bm_rules_from_template_rxn_map(self, template_rxn_map, nprocs=1, Tref=1 if nprocs > 1: pool = mp.Pool(nprocs) - kinetics_list = np.array(pool.map(_make_rule, inputs[inds])) + kinetics_list = np.array(pool.map(_make_rule, list(inputs[i] for i in inds))) else: - kinetics_list = np.array(list(map(_make_rule, inputs[inds]))) + kinetics_list = np.array(list(map(_make_rule, list(inputs[i] for i in inds)))) kinetics_list = kinetics_list[revinds] # fix order diff --git a/rmgpy/exceptions.py b/rmgpy/exceptions.py index 7b4420cf66a..f2e7dae4602 100644 --- a/rmgpy/exceptions.py +++ b/rmgpy/exceptions.py @@ -109,15 +109,6 @@ class ForbiddenStructureException(Exception): pass -class ILPSolutionError(Exception): - """ - An exception to be raised when solving an integer linear programming problem if a solution - could not be found or the solution is not valid. Can pass a string to indicate the reason - that the solution is invalid. - """ - pass - - class ImplicitBenzeneError(Exception): """ An exception class when encountering a group with too many implicit benzene diff --git a/rmgpy/molecule/draw.py b/rmgpy/molecule/draw.py index 6dd6c8b46b1..3792772af32 100644 --- a/rmgpy/molecule/draw.py +++ b/rmgpy/molecule/draw.py @@ -49,6 +49,7 @@ import math import os.path import re +import itertools try: import cairocffi as cairo @@ -60,7 +61,8 @@ import numpy as np from rdkit.Chem import AllChem -from rmgpy.molecule.molecule import Atom, Molecule +from rmgpy.molecule.molecule import Atom, Molecule, Bond +from rmgpy.molecule.pathfinder import find_shortest_path from rmgpy.qm.molecule import Geometry @@ -96,6 +98,12 @@ def create_new_surface(file_format, target=None, width=1024, height=768): ################################################################################ +class AdsorbateDrawingError(Exception): + """ + When something goes wrong trying to draw an adsorbate. + """ + pass + class MoleculeDrawer(object): """ This class provides functionality for drawing the skeletal formula of @@ -207,7 +215,16 @@ def draw(self, molecule, file_format, target=None): # replace the bonds after generating coordinates. This avoids # bugs with RDKit old_bond_dictionary = self._make_single_bonds() - self._generate_coordinates() + if molecule.contains_surface_site(): + try: + self._connect_surface_sites() + self._generate_coordinates() + self._disconnect_surface_sites() + except AdsorbateDrawingError as e: + self._disconnect_surface_sites() + self._generate_coordinates(fix_surface_sites=False) + else: + self._generate_coordinates() self._replace_bonds(old_bond_dictionary) # Generate labels to use @@ -323,11 +340,13 @@ def _find_ring_groups(self): if not found: self.ringSystems.append([cycle]) - def _generate_coordinates(self): + def _generate_coordinates(self, fix_surface_sites=True): """ Generate the 2D coordinates to be used when drawing the current molecule. The function uses rdKits 2D coordinate generation. Updates the self.coordinates Array in place. + If `fix_surface_sites` is True, then the surface sites are placed + at the bottom of the molecule. """ atoms = self.molecule.atoms natoms = len(atoms) @@ -390,7 +409,6 @@ def _generate_coordinates(self): # If two atoms lie on top of each other, push them apart a bit # This is ugly, but at least the mess you end up with isn't as misleading # as leaving everything piled on top of each other at the origin - import itertools for atom1, atom2 in itertools.combinations(backbone, 2): i1, i2 = atoms.index(atom1), atoms.index(atom2) if np.linalg.norm(coordinates[i1, :] - coordinates[i2, :]) < 0.5: @@ -402,7 +420,6 @@ def _generate_coordinates(self): # If two atoms lie on top of each other, push them apart a bit # This is ugly, but at least the mess you end up with isn't as misleading # as leaving everything piled on top of each other at the origin - import itertools for atom1, atom2 in itertools.combinations(backbone, 2): i1, i2 = atoms.index(atom1), atoms.index(atom2) if np.linalg.norm(coordinates[i1, :] - coordinates[i2, :]) < 0.5: @@ -457,26 +474,59 @@ def _generate_coordinates(self): coordinates[:, 0] = temp[:, 1] coordinates[:, 1] = temp[:, 0] - # For surface species, rotate them so the site is at the bottom. - if self.molecule.contains_surface_site(): + # For surface species + if fix_surface_sites and self.molecule.contains_surface_site(): if len(self.molecule.atoms) == 1: return coordinates - for site in self.molecule.atoms: - if site.is_surface_site(): - break - else: - raise Exception("Can't find surface site") - if site.bonds: - adsorbate = next(iter(site.bonds)) - vector0 = coordinates[atoms.index(site), :] - coordinates[atoms.index(adsorbate), :] - angle = math.atan2(vector0[0], vector0[1]) - math.pi + sites = [atom for atom in self.molecule.atoms if atom.is_surface_site()] + if len(sites) == 1: + # rotate them so the site is at the bottom. + site = sites[0] + if site.bonds: + adatom = next(iter(site.bonds)) + vector0 = coordinates[atoms.index(site), :] - coordinates[atoms.index(adatom), :] + angle = math.atan2(vector0[0], vector0[1]) - math.pi + rot = np.array([[math.cos(angle), math.sin(angle)], [-math.sin(angle), math.cos(angle)]], float) + self.coordinates = coordinates = np.dot(coordinates, rot) + else: + # van der Waals + index = atoms.index(site) + coordinates[index, 1] = min(coordinates[:, 1]) - 0.8 # just move the site down a bit + coordinates[index, 0] = coordinates[:, 0].mean() # and center it + elif len(sites) <= 4: + # Rotate so the line of best fit through the adatoms is horizontal. + # find atoms bonded to sites + adatoms = [next(iter(site.bonds)) for site in sites] + adatom_indices = [atoms.index(a) for a in adatoms] + # find the best fit line through the bonded atoms + x = coordinates[adatom_indices, 0] + y = coordinates[adatom_indices, 1] + A = np.vstack([x, np.ones(len(x))]).T + m, c = np.linalg.lstsq(A, y, rcond=None)[0] + # rotate so the line is horizontal + angle = -math.atan(m) rot = np.array([[math.cos(angle), math.sin(angle)], [-math.sin(angle), math.cos(angle)]], float) self.coordinates = coordinates = np.dot(coordinates, rot) + # if the line is above the middle, flip it + not_site_indices = [atoms.index(a) for a in atoms if not a.is_surface_site()] + if coordinates[adatom_indices, 1].mean() > coordinates[not_site_indices, 1].mean(): + coordinates[:, 1] *= -1 + x = coordinates[adatom_indices, 0] + y = coordinates[adatom_indices, 1] + site_y_pos = min(min(y) - 0.8, min(coordinates[not_site_indices, 1]) - 0.5) + if max(y) - site_y_pos > 1.5: + raise AdsorbateDrawingError("Adsorbate bond too long") + for x1, x2 in itertools.combinations(x, 2): + if abs(x1 - x2) < 0.2: + raise AdsorbateDrawingError("Sites overlapping") + for site, x_pos in zip(sites, x): + index = atoms.index(site) + coordinates[index, 1] = site_y_pos + coordinates[index, 0] = x_pos + else: - # van der waals - index = atoms.index(site) - coordinates[index, 1] = min(coordinates[:, 1]) - 0.8 # just move the site down a bit - coordinates[index, 0] = coordinates[:, 0].mean() # and center it + # more than 4 surface sites? leave them alone + pass def _find_cyclic_backbone(self): """ @@ -854,7 +904,7 @@ def _generate_functional_group_coordinates(self, atom0, atom1): # Check to see if atom1 is in any cycles in the molecule ring_system = None for ring_sys in self.ringSystems: - if any([atom1 in ring for ring in ring_sys]): + if any(atom1 in ring for ring in ring_sys): ring_system = ring_sys if ring_system is not None: @@ -1624,6 +1674,40 @@ def _replace_bonds(self, bond_order_dictionary): for bond, order in bond_order_dictionary.items(): bond.set_order_num(order) + def _connect_surface_sites(self): + """ + Creates single bonds between atoms that are surface sites. + This is to help make multidentate adsorbates look better. + """ + sites = [a for a in self.molecule.atoms if a.is_surface_site()] + if len(sites) > 4: + return + for site1 in sites: + other_sites = [a for a in sites if a != site1] + if not other_sites: break + # connect to the nearest site + site2 = min(other_sites, key=lambda a: len(find_shortest_path(site1, a))) + if len(find_shortest_path(site1, site2)) > 2 and len(sites) > 3: + # if there are more than 3 sites, don't connect sites that aren't neighbors + continue + + bond = site1.bonds.get(site2) + if bond is None: + bond = Bond(site1, site2, 1) + site1.bonds[site2] = bond + site2.bonds[site1] = bond + + def _disconnect_surface_sites(self): + """ + Removes all bonds between atoms that are surface sites. + """ + for site1 in self.molecule.atoms: + if site1.is_surface_site(): + for site2 in list(site1.bonds.keys()): # make a list copy so we can delete from the dict + if site2.is_surface_site(): + del site1.bonds[site2] + del site2.bonds[site1] + ################################################################################ @@ -1644,6 +1728,7 @@ def __init__(self, options=None): self.options = MoleculeDrawer().options.copy() self.options.update({ 'arrowLength': 36, + 'drawReversibleArrow': True }) if options: self.options.update(options) @@ -1744,11 +1829,30 @@ def draw(self, reaction, file_format, path=None): rxn_cr.save() rxn_cr.set_source_rgba(0.0, 0.0, 0.0, 1.0) rxn_cr.set_line_width(1.0) - rxn_cr.move_to(rxn_x + 8, rxn_top + 0.5 * rxn_height) - rxn_cr.line_to(rxn_x + arrow_width - 8, rxn_top + 0.5 * rxn_height) - rxn_cr.move_to(rxn_x + arrow_width - 14, rxn_top + 0.5 * rxn_height - 3.0) - rxn_cr.line_to(rxn_x + arrow_width - 8, rxn_top + 0.5 * rxn_height) - rxn_cr.line_to(rxn_x + arrow_width - 14, rxn_top + 0.5 * rxn_height + 3.0) + if self.options['drawReversibleArrow'] and reaction.reversible: # draw double harpoons + TOP_HARPOON_Y = rxn_top + (0.5 * rxn_height - 1.75) + BOTTOM_HARPOON_Y = rxn_top + (0.5 * rxn_height + 1.75) + + # Draw top harpoon + rxn_cr.move_to(rxn_x + 8, TOP_HARPOON_Y) + rxn_cr.line_to(rxn_x + arrow_width - 8, TOP_HARPOON_Y) + rxn_cr.move_to(rxn_x + arrow_width - 14, TOP_HARPOON_Y - 3.0) + rxn_cr.line_to(rxn_x + arrow_width - 8, TOP_HARPOON_Y) + + # Draw bottom harpoon + rxn_cr.move_to(rxn_x + arrow_width - 8, BOTTOM_HARPOON_Y) + rxn_cr.line_to(rxn_x + 8, BOTTOM_HARPOON_Y) + rxn_cr.move_to(rxn_x + 14, BOTTOM_HARPOON_Y + 3.0) + rxn_cr.line_to(rxn_x + 8, BOTTOM_HARPOON_Y) + + + else: # draw forward arrow + rxn_cr.move_to(rxn_x + 8, rxn_top + 0.5 * rxn_height) + rxn_cr.line_to(rxn_x + arrow_width - 8, rxn_top + 0.5 * rxn_height) + rxn_cr.move_to(rxn_x + arrow_width - 14, rxn_top + 0.5 * rxn_height - 3.0) + rxn_cr.line_to(rxn_x + arrow_width - 8, rxn_top + 0.5 * rxn_height) + rxn_cr.line_to(rxn_x + arrow_width - 14, rxn_top + 0.5 * rxn_height + 3.0) + rxn_cr.stroke() rxn_cr.restore() rxn_x += arrow_width diff --git a/rmgpy/molecule/resonance.pxd b/rmgpy/molecule/resonance.pxd index 64662a797c3..a28c0aa1adc 100644 --- a/rmgpy/molecule/resonance.pxd +++ b/rmgpy/molecule/resonance.pxd @@ -25,6 +25,8 @@ # # ############################################################################### +cimport numpy as cnp + from rmgpy.molecule.graph cimport Vertex, Edge, Graph from rmgpy.molecule.molecule cimport Atom, Bond, Molecule @@ -64,6 +66,4 @@ cpdef list generate_kekule_structure(Graph mol) cpdef list generate_clar_structures(Graph mol, bint save_order=?) -cpdef list _clar_optimization(Graph mol, list constraints=?, max_num=?, save_order=?) - -cpdef list _clar_transformation(Graph mol, list aromatic_ring) +cpdef list _clar_optimization(Graph mol, list recursion_constraints=?, int clar_number=?, bint save_order=?) diff --git a/rmgpy/molecule/resonance.py b/rmgpy/molecule/resonance.py index 24eccaa636e..3cb0bad9b88 100644 --- a/rmgpy/molecule/resonance.py +++ b/rmgpy/molecule/resonance.py @@ -55,15 +55,17 @@ from operator import attrgetter import cython +import numpy as np +from scipy.optimize import Bounds, LinearConstraint, milp import rmgpy.molecule.filtration as filtration import rmgpy.molecule.pathfinder as pathfinder -from rmgpy.exceptions import ILPSolutionError, KekulizationError, AtomTypeError, ResonanceError +from rmgpy.exceptions import AtomTypeError, KekulizationError, ResonanceError from rmgpy.molecule.adjlist import Saturator +from rmgpy.molecule.fragment import CuttingLabel from rmgpy.molecule.graph import Vertex from rmgpy.molecule.kekulize import kekulize from rmgpy.molecule.molecule import Atom, Bond, Molecule -from rmgpy.molecule.fragment import CuttingLabel def populate_resonance_algorithms(features=None): @@ -906,8 +908,7 @@ def generate_clar_structures(mol, save_order=False): Returns a list of :class:`Molecule` objects corresponding to the Clar structures. """ - cython.declare(output=list, mol_list=list, new_mol=Graph, aromatic_rings=list, bonds=list, solution=list, - y=list, x=list, index=cython.int, bond=Edge, ring=list) + cython.declare(output=list, mol_list=list, new_mol=Graph, aromatic_rings=list, bonds=list, index=cython.int, bond=Edge, ring=list) if not mol.is_cyclic(): return [] @@ -917,19 +918,18 @@ def generate_clar_structures(mol, save_order=False): mol.assign_atom_ids() try: - output = _clar_optimization(mol, save_order=save_order) - except ILPSolutionError: + solutions = _clar_optimization(mol, save_order=save_order) + except (RuntimeError, ValueError): # either a crash during optimization or the result was an empty tuple # The optimization algorithm did not work on the first iteration return [] mol_list = [] - for new_mol, aromatic_rings, bonds, solution in output: - + for new_mol, aromatic_rings, bonds, solution in solutions: # The solution includes a part corresponding to rings, y, and a part corresponding to bonds, x, using # nomenclature from the paper. In y, 1 means the ring as a sextet, 0 means it does not. # In x, 1 corresponds to a double bond, 0 either means a single bond or the bond is part of a sextet. - y = solution[0:len(aromatic_rings)] + y = solution[:len(aromatic_rings)] x = solution[len(aromatic_rings):] # Apply results to molecule - double bond locations first @@ -939,18 +939,21 @@ def generate_clar_structures(mol, save_order=False): elif x[index] == 1: bond.order = 2 # double else: - raise ValueError('Unaccepted bond value {0} obtained from optimization.'.format(x[index])) + raise ValueError(f'Unaccepted bond value {x[index]} obtained from optimization.') # Then apply locations of aromatic sextets by converting to benzene bonds for index, ring in enumerate(aromatic_rings): if y[index] == 1: - _clar_transformation(new_mol, ring) + for i, atom_1 in enumerate(ring): + for j, atom_2 in enumerate(ring): + if new_mol.has_bond(atom_1, atom_2): + new_mol.get_bond(atom_1, atom_2).order = 1.5 try: new_mol.update_atomtypes() except AtomTypeError: pass - else: + finally: mol_list.append(new_mol) return mol_list @@ -963,165 +966,175 @@ def _tuplize_bond(bond): return (bond.atom1.id, bond.atom2.id) -def _clar_optimization(mol, constraints=None, max_num=None, save_order=False): +def _clar_optimization(mol, recursion_constraints=None, clar_number=-1, save_order=False): """ - Implements linear programming algorithm for finding Clar structures. This algorithm maximizes the number - of Clar sextets within the constraints of molecular geometry and atom valency. + Implements Mixed Integer Linear Programming for finding Clar structures. + + First finds the Clar number (and an arbitrary structure with that number), then recursively + calls itself to enumerate more structures with that Clar number. No guarantees about + which structures will be found, or how many (we stop solving once solution would require + branching, which typically happens after at least a couple have already been found). Returns a list of valid Clar solutions in the form of a tuple, with the following entries: - [0] Molecule object + [0] Copy of mol [1] List of aromatic rings [2] List of bonds - [3] Optimization solution + [3] Solution vector - The optimization solution is a list of boolean values with sextet assignments followed by double bond assignments, - with indices corresponding to the list of aromatic rings and list of bonds, respectively. + The solution vector is a binary integer list of length (# aromatic rings + # bonds) indicating + if each ring is aromatic and if each bond is double. - Method adapted from: - Hansen, P.; Zheng, M. The Clar Number of a Benzenoid Hydrocarbon and Linear Programming. - J. Math. Chem. 1994, 15 (1), 93–107. + This implementation follows the original implementation very closely (Hansen, P.; Zheng, M. The + Clar Number of a Benzenoid Hydrocarbon and Linear Programming. J. Math. Chem. 1994, 15 (1), 93–107.) + with only slight modifications to prevent changing bond orders for exocyclic atoms. """ - cython.declare(molecule=Graph, aromatic_rings=list, exo=list, l=cython.int, m=cython.int, n=cython.int, - a=list, objective=list, status=cython.int, solution=list, innerSolutions=list) - - from lpsolve55 import lpsolve - - # Make a copy of the molecule so we don't destroy the original + cython.declare( + molecule=Graph, + aromatic_rings=list, + exo_info=list, + n_ring=cython.int, + n_atom=cython.int, + n_bond=cython.int, + A=list, + solutions=list, + ) + + # after we find all the Clar structures we will modify the molecule bond orders to be aromatic, + # so we make explicit copies to avoid overwrites. molecule = mol.copy(deep=True) aromatic_rings = molecule.get_aromatic_rings(save_order=save_order)[0] - aromatic_rings.sort(key=_sum_atom_ids) - if not aromatic_rings: + # doesn't work for system without aromatic rings, just exit + if len(aromatic_rings) == 0: return [] - # Get list of atoms that are in rings + # stability between multiple runs + aromatic_rings.sort(key=_sum_atom_ids) + + # Cython doesn't allow mutable defaults, so we just set it here instead + if recursion_constraints is None: + recursion_constraints = [] + + # Get set of atoms that are in rings atoms = set() for ring in aromatic_rings: atoms.update(ring) atoms = sorted(atoms, key=attrgetter('id')) - # Get list of bonds involving the ring atoms, ignoring bonds to hydrogen + # Get set of bonds involving the ring atoms, ignoring bonds to hydrogen bonds = set() for atom in atoms: bonds.update([atom.bonds[key] for key in atom.bonds.keys() if key.is_non_hydrogen()]) bonds = sorted(bonds, key=_tuplize_bond) - # Identify exocyclic bonds, and save their bond orders - exo = [] + # identify exocyclic bonds, and save their order if exocyclic + exo_info = [] for bond in bonds: if bond.atom1 not in atoms or bond.atom2 not in atoms: + # save order for exocyclic if bond.is_double(): - exo.append(1) + exo_info.append(1) else: - exo.append(0) - else: - exo.append(None) + exo_info.append(0) + else: # not exocyclic + exo_info.append(None) # Dimensions - l = len(aromatic_rings) - m = len(atoms) - n = l + len(bonds) + n_ring = len(aromatic_rings) + n_atom = len(atoms) + n_bond = len(bonds) + + # The aromaticity assignment problem is formulated as an MILP problem + # minimize: + # c @ x + # such that + # b_l <= A @ x <= b_u + # l <= x <= u + # x are integers # Connectivity matrix which indicates which rings and bonds each atom is in - # Part of equality constraint Ax=b - a = [] + # Part of equality constraint A_eq @ x = b_eq + A = [] for atom in atoms: in_ring = [1 if atom in ring else 0 for ring in aromatic_rings] in_bond = [1 if atom in [bond.atom1, bond.atom2] else 0 for bond in bonds] - a.append(in_ring + in_bond) + A.append(in_ring + in_bond) + initial_constraint = [LinearConstraint( # i.e. an equality constraint + A=np.array(A, dtype=int), lb=np.ones(n_atom, dtype=int), ub=np.ones(n_atom, dtype=int) + )] + + # on recursive calls we already know the Clar number, so we can additionally constrain the system + # to only find structures with this number. Allows detecting when all formulae have been found and, + # in theory, should solve faster + if clar_number != -1: + initial_constraint += [ + LinearConstraint( + A=np.array([1] * n_ring + [0] * n_bond), + ub=clar_number, + lb=clar_number, + ), + ] # Objective vector for optimization: sextets have a weight of 1, double bonds have a weight of 0 - objective = [1] * l + [0] * len(bonds) - - # Solve LP problem using lpsolve - lp = lpsolve('make_lp', m, n) # initialize lp with constraint matrix with m rows and n columns - lpsolve('set_verbose', lp, 2) # reduce messages from lpsolve - lpsolve('set_obj_fn', lp, objective) # set objective function - lpsolve('set_maxim', lp) # set solver to maximize objective - lpsolve('set_mat', lp, a) # set left hand side to constraint matrix - lpsolve('set_rh_vec', lp, [1] * m) # set right hand side to 1 for all constraints - for i in range(m): # set all constraints as equality constraints - lpsolve('set_constr_type', lp, i + 1, '=') - lpsolve('set_binary', lp, [True] * n) # set all variables to be binary - - # Constrain values of exocyclic bonds, since we don't want to modify them - for i in range(l, n): - if exo[i - l] is not None: - # NOTE: lpsolve indexes from 1, so the variable we're changing should be i + 1 - lpsolve('set_bounds', lp, i + 1, exo[i - l], exo[i - l]) - - # Add constraints to problem if provided - if constraints is not None: - for constraint in constraints: - try: - lpsolve('add_constraint', lp, constraint[0], '<=', constraint[1]) - except Exception as e: - logging.debug('Unable to add constraint: {0} <= {1}'.format(constraint[0], constraint[1])) - logging.debug(mol.to_adjacency_list()) - if str(e) == 'invalid vector.': - raise ILPSolutionError('Unable to add constraint, likely due to ' - 'inconsistent aromatic ring perception.') - else: - raise - - status = lpsolve('solve', lp) - obj_val, solution = lpsolve('get_solution', lp)[0:2] - lpsolve('delete_lp', lp) # Delete the LP problem to clear up memory + # we negate because the original problem is formulated as maximization, whereas SciPy only does min + c = -np.array([1] * n_ring + [0] * n_bond, dtype=int) + + # variable bounds + # - binary problem, so everything must be either zero or one + # - rings are simply 0 or 1 + # - bonds are also 0 or 1, except exocyclic bonds which must remain unchanged + bounds = Bounds( + lb=np.array([0] * n_ring + [0 if b is None else b for b in exo_info], dtype=int), + ub=np.array([1] * n_ring + [1 if b is None else b for b in exo_info], dtype=int), + ) + + result = milp( + c=c, + integrality=1, + bounds=bounds, + constraints=initial_constraint + recursion_constraints, + options={'time_limit': 10}, + ) + + if (status := result.status) != 0: + if status == 2: # infeasible + raise RuntimeError("All valid Clar formulae have been enumerated!") + else: + raise RuntimeError(f"Optimization failed (Exit Code {status}) for an unexpected reason '{result.message}'") - # Check that optimization was successful - if status != 0: - raise ILPSolutionError('Optimization could not find a valid solution.') + _clar_number, solution = -result.fun, result.x - # Check that we the result contains at least one aromatic sextet - if obj_val == 0: + # optimization may have reached a bad local minimum - this case is rare + if _clar_number == 0: return [] - # Check that the solution contains the maximum number of sextets possible - if max_num is None: - max_num = obj_val # This is the first solution, so the result should be an upper limit - elif obj_val < max_num: - raise ILPSolutionError('Optimization obtained a sub-optimal solution.') - + # on later runs, non-integer solutions occur - branching might be able to find actual solutions, + # but we just call it good enough here if any([x != 1 and x != 0 for x in solution]): - raise ILPSolutionError('Optimization obtained a non-integer solution.') - - # Generate constraints based on the solution obtained - y = solution[0:l] - new_a = y + [0] * len(bonds) - new_b = sum(y) - 1 - if constraints is not None: - constraints.append((new_a, new_b)) - else: - constraints = [(new_a, new_b)] + raise RuntimeError("Optimization obtained a non-integer solution - no more formulae will be enumerated.") + + # first solution, so the result should be an upper limit + if clar_number == -1: + clar_number = _clar_number + + selected_sextets = list(solution[:n_ring]) + # restrict those rings which were selected for the current clar structure + # from forming it again by requiring that their sum is 1 less than it used + # to be + recursion_constraints += [ + LinearConstraint( + A=np.array(selected_sextets + [0] * n_bond), + ub=clar_number - 1, + lb=0, + ), + ] # Run optimization with additional constraints try: - inner_solutions = _clar_optimization(mol, constraints=constraints, max_num=max_num, save_order=save_order) - except ILPSolutionError: + inner_solutions = _clar_optimization(mol, recursion_constraints=recursion_constraints, clar_number=clar_number, save_order=save_order) + except RuntimeError as e: + logging.debug(f"Clar Optimization stopped: {e}") inner_solutions = [] return inner_solutions + [(molecule, aromatic_rings, bonds, solution)] - - -def _clar_transformation(mol, aromatic_ring): - """ - Performs Clar transformation for given ring in a molecule, ie. conversion to aromatic sextet. - - Args: - mol a :class:`Molecule` object - aromaticRing a list of :class:`Atom` objects corresponding to an aromatic ring in mol - - This function directly modifies the input molecule and does not return anything. - """ - cython.declare(bondList=list, i=cython.int, atom1=Vertex, atom2=Vertex, bond=Edge) - - bond_list = [] - - for i, atom1 in enumerate(aromatic_ring): - for atom2 in aromatic_ring[i + 1:]: - if mol.has_bond(atom1, atom2): - bond_list.append(mol.get_bond(atom1, atom2)) - - for bond in bond_list: - bond.order = 1.5 diff --git a/rmgpy/molecule/util.py b/rmgpy/molecule/util.py index d1093f45821..46f8c3504c6 100644 --- a/rmgpy/molecule/util.py +++ b/rmgpy/molecule/util.py @@ -30,8 +30,8 @@ import itertools import re -from rmgpy.molecule.molecule import Molecule from rmgpy.molecule.fragment import Fragment +from rmgpy.molecule.molecule import Molecule def get_element_count(obj): @@ -47,7 +47,7 @@ def get_element_count(obj): match = re.match(r"([a-z]+)([0-9]*)", piece, re.I) if match: element, count = match.groups() - if count is '': + if count == '': count = 1 if element in element_count: element_count[element] += int(count) diff --git a/rmgpy/pdep/sls.py b/rmgpy/pdep/sls.py index a3474ac6e51..f8f66d5b836 100644 --- a/rmgpy/pdep/sls.py +++ b/rmgpy/pdep/sls.py @@ -31,17 +31,29 @@ Contains functionality for directly simulating the master equation and implementing the SLS master equation reduction method """ +import logging +import sys -from juliacall import Main -Main.seval("using ReactionMechanismSimulator.SciMLBase") -Main.seval("using ReactionMechanismSimulator.Sundials") import numpy as np import scipy.linalg import scipy.optimize as opt +import scipy.sparse as sparse import rmgpy.constants as constants from rmgpy.pdep.me import generate_full_me_matrix, states_to_configurations -from rmgpy.rmg.reactors import to_julia +from rmgpy.rmg.reactionmechanismsimulator_reactors import to_julia + +NO_JULIA = False +try: + from juliacall import Main + Main.seval("using ReactionMechanismSimulator.SciMLBase") + Main.seval("using ReactionMechanismSimulator.Sundials") +except Exception as e: + logging.info( + f"Unable to import Julia dependencies, original error: {str(e)}" + ". Master equation method 'ode' will not be available on this execution." + ) + NO_JULIA = True def get_initial_condition(network, x0, indices): @@ -153,6 +165,11 @@ def get_rate_coefficients_SLS(network, T, P, method="mexp", neglect_high_energy_ tau = np.abs(1.0 / fastest_reaction) if method == "ode": + if NO_JULIA: + raise RuntimeError( + "Required Julia dependencies for method 'ode' are not installed.\n" + "Please follow the steps to install Julia dependencies at https://reactionmechanismgenerator.github.io/RMG-Py/users/rmg/installation/anacondaDeveloper.html." + ) f = Main.seval( """ function f(du, u, M, t) diff --git a/rmgpy/qm/molecule.py b/rmgpy/qm/molecule.py index 81d771781d0..45cf457b672 100644 --- a/rmgpy/qm/molecule.py +++ b/rmgpy/qm/molecule.py @@ -408,19 +408,28 @@ def parse(self): 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 = ( + ) # cf. https://cclib.github.io/index.html#how-to-use-cclib + parser.molmass = None # give it an attribute and it won't delete it, leaving it on the parser object + parser.rotcons = ( # for cclib < 1.8.0 [] - ) # 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() + ) + parser.rotconsts = ( # for cclib >= 1.8.0 + [] + ) + cclib_data = parser.parse() # fills in either parser.rotcons or parser.rotconsts but not both + assert bool(parser.rotconsts) != bool(parser.rotcons) + if parser.rotcons: # for cclib < 1.8.0 + cclib_data.rotcons = ( + parser.rotcons + ) + else: # for cclib >= 1.8.0 + cclib_data.rotconsts = ( + parser.rotconsts + ) 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 + ) # this hack required because molmass is 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. diff --git a/rmgpy/qm/qmdata.py b/rmgpy/qm/qmdata.py index 976649ddf7a..6bf0b6a44e8 100644 --- a/rmgpy/qm/qmdata.py +++ b/rmgpy/qm/qmdata.py @@ -98,7 +98,11 @@ def parse_cclib_data(cclib_data, ground_state_degeneracy): molecular_mass = None energy = (cclib_data.scfenergies[-1], 'eV/molecule') atomic_numbers = cclib_data.atomnos - rotational_constants = (cclib_data.rotcons[-1], 'cm^-1') + if hasattr(cclib_data, 'rotconsts'): + rotational_constants = (cclib_data.rotconsts[-1], 'cm^-1') + else: + rotational_constants = (cclib_data.rotcons[-1], 'cm^-1') + atom_coords = (cclib_data.atomcoords[-1], 'angstrom') frequencies = (cclib_data.vibfreqs, 'cm^-1') diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index d902b844b9d..510f618c889 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -1542,8 +1542,8 @@ def generate_high_p_limit_kinetics(self): """ raise NotImplementedError("generate_high_p_limit_kinetics is not implemented for all Reaction subclasses.") -def same_object(object1, object2, _check_identical, _only_check_label, - _generate_initial_map, _strict, _save_order): +def _same_object(object1, object2, _check_identical=False, _only_check_label=False, + _generate_initial_map=False, _strict=True, _save_order=False): if _only_check_label: return str(object1) == str(object2) elif _check_identical: @@ -1573,10 +1573,10 @@ def same_species_lists(list1, list2, check_identical=False, only_check_label=Fal """ same_object_passthrough = partial( - same_object, + _same_object, _check_identical=check_identical, _only_check_label=only_check_label, - _generate_intial_map=generate_initial_map, + _generate_initial_map=generate_initial_map, _strict=strict, _save_order=save_order, ) diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index e04cd82e2ed..ba9b86f87a8 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -48,9 +48,10 @@ from rmgpy.solver.surface import SurfaceReactor from rmgpy.util import as_list from rmgpy.data.surface import MetalDatabase -from rmgpy.rmg.reactors import Reactor, ConstantVIdealGasReactor, ConstantTLiquidSurfaceReactor, ConstantTVLiquidReactor, ConstantTPIdealGasReactor from rmgpy.data.vaporLiquidMassTransfer import liquidVolumetricMassTransferCoefficientPowerLaw from rmgpy.molecule.fragment import Fragment +from rmgpy.rmg.reactionmechanismsimulator_reactors import Reactor, ConstantVIdealGasReactor, ConstantTLiquidSurfaceReactor, ConstantTVLiquidReactor, ConstantTPIdealGasReactor, NO_JULIA + ################################################################################ @@ -1558,6 +1559,11 @@ def read_input_file(path, rmg0): exec(f.read(), global_context, local_context) except (NameError, TypeError, SyntaxError) as e: logging.error('The input file "{0}" was invalid:'.format(full_path)) + if NO_JULIA: + logging.error( + "During runtime, import of Julia dependencies failed. To use phase systems and RMS reactors, install RMG-Py with RMS." + " (https://reactionmechanismgenerator.github.io/RMG-Py/users/rmg/installation/anacondaDeveloper.html)" + ) logging.exception(e) raise finally: diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index b5509012834..d9b3a1f3669 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -52,8 +52,6 @@ from scipy.optimize import brute import rmgpy.util as util -from rmgpy.rmg.model import Species, CoreEdgeReactionModel -from rmgpy.rmg.pdep import PDepNetwork from rmgpy import settings from rmgpy.chemkin import ChemkinWriter from rmgpy.constraints import fails_species_constraints @@ -61,26 +59,31 @@ from rmgpy.data.kinetics.family import TemplateReaction from rmgpy.data.kinetics.library import KineticsLibrary, LibraryReaction from rmgpy.data.rmg import RMGDatabase -from rmgpy.exceptions import ForbiddenStructureException, DatabaseError, CoreError, InputError -from rmgpy.kinetics.diffusionLimited import diffusion_limiter from rmgpy.data.vaporLiquidMassTransfer import vapor_liquid_mass_transfer -from rmgpy.kinetics import ThirdBody -from rmgpy.kinetics import Troe +from rmgpy.exceptions import ( + CoreError, + DatabaseError, + ForbiddenStructureException, + InputError, +) +from rmgpy.kinetics import ThirdBody, Troe +from rmgpy.kinetics.diffusionLimited import diffusion_limiter from rmgpy.molecule import Molecule from rmgpy.qm.main import QMDatabaseWriter from rmgpy.reaction import Reaction -from rmgpy.rmg.listener import SimulationProfileWriter, SimulationProfilePlotter +from rmgpy.rmg.listener import SimulationProfilePlotter, SimulationProfileWriter +from rmgpy.rmg.model import CoreEdgeReactionModel, Species from rmgpy.rmg.output import OutputHTMLWriter -from rmgpy.rmg.pdep import PDepReaction +from rmgpy.rmg.pdep import PDepNetwork, PDepReaction +from rmgpy.rmg.reactionmechanismsimulator_reactors import Reactor as RMSReactor from rmgpy.rmg.settings import ModelSettings -from rmgpy.solver.base import TerminationTime, TerminationConversion +from rmgpy.solver.base import TerminationConversion, TerminationTime from rmgpy.solver.simple import SimpleReactor from rmgpy.stats import ExecutionStatsWriter from rmgpy.thermo.thermoengine import submit from rmgpy.tools.plot import plot_sensitivity from rmgpy.tools.uncertainty import Uncertainty, process_local_results from rmgpy.yml import RMSWriter -from rmgpy.rmg.reactors import Reactor ################################################################################ @@ -269,12 +272,6 @@ def load_input(self, path=None): if self.solvent: self.reaction_model.solvent_name = self.solvent - if self.surface_site_density: - self.reaction_model.surface_site_density = self.surface_site_density - self.reaction_model.core.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si - self.reaction_model.edge.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si - self.reaction_model.coverage_dependence = self.coverage_dependence - self.reaction_model.verbose_comments = self.verbose_comments self.reaction_model.save_edge_species = self.save_edge_species @@ -507,6 +504,19 @@ def initialize(self, **kwargs): # Read input file self.load_input(self.input_file) + + # Check if ReactionMechanismSimulator reactors are being used + # if RMS is not installed but the user attempted to use it, the load_input_file would have failed + # if RMS is installed but they did not use it, we can avoid extra work + # if RMS is not installed and they did not use it, we avoid calling certain functions that would raise an error + requires_rms = any(isinstance(reactor_system, RMSReactor) for reactor_system in self.reaction_systems) + + if self.surface_site_density: + self.reaction_model.surface_site_density = self.surface_site_density + if requires_rms: + self.reaction_model.core.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si + self.reaction_model.edge.phase_system.phases["Surface"].site_density = self.surface_site_density.value_si + self.reaction_model.coverage_dependence = self.coverage_dependence if kwargs.get("restart", ""): import rmgpy.rmg.input @@ -550,10 +560,10 @@ def initialize(self, **kwargs): self.load_database() for spec in self.initial_species: - self.reaction_model.add_species_to_edge(spec) + self.reaction_model.add_species_to_edge(spec, requires_rms=requires_rms) for reaction_system in self.reaction_systems: - if isinstance(reaction_system, Reactor): + if isinstance(reaction_system, RMSReactor): reaction_system.finish_termination_criteria() # Load restart seed mechanism (if specified) @@ -583,8 +593,9 @@ def initialize(self, **kwargs): # Do all liquid-phase startup things: if self.solvent: solvent_data = self.database.solvation.get_solvent_data(self.solvent) - self.reaction_model.core.phase_system.phases["Default"].set_solvent(solvent_data) - self.reaction_model.edge.phase_system.phases["Default"].set_solvent(solvent_data) + if requires_rms: + self.reaction_model.core.phase_system.phases["Default"].set_solvent(solvent_data) + self.reaction_model.edge.phase_system.phases["Default"].set_solvent(solvent_data) diffusion_limiter.enable(solvent_data, self.database.solvation) logging.info("Setting solvent data for {0}".format(self.solvent)) @@ -618,12 +629,12 @@ def initialize(self, **kwargs): # Seed mechanisms: add species and reactions from seed mechanism # DON'T generate any more reactions for the seed species at this time for seed_mechanism in self.seed_mechanisms: - self.reaction_model.add_seed_mechanism_to_core(seed_mechanism, react=False) + self.reaction_model.add_seed_mechanism_to_core(seed_mechanism, react=False, requires_rms=requires_rms) # Reaction libraries: add species and reactions from reaction library to the edge so # that RMG can find them if their rates are large enough for library, option in self.reaction_libraries: - self.reaction_model.add_reaction_library_to_edge(library) + self.reaction_model.add_reaction_library_to_edge(library, requires_rms=requires_rms) # Also always add in a few bath gases (since RMG-Java does) for label, smiles in [("Ar", "[Ar]"), ("He", "[He]"), ("Ne", "[Ne]"), ("N2", "N#N")]: @@ -695,35 +706,37 @@ def initialize(self, **kwargs): # This is necessary so that the PDep algorithm can identify the bath gas for spec in self.initial_species: if not spec.reactive: - self.reaction_model.enlarge(spec) + self.reaction_model.enlarge(spec, requires_rms=requires_rms) for spec in self.initial_species: if spec.reactive: - self.reaction_model.enlarge(spec) + self.reaction_model.enlarge(spec, requires_rms=requires_rms) # chatelak: store constant SPC indices in the reactor attributes if any constant SPC provided in the input file # advantages to write it here: this is run only once (as species indexes does not change over the generation) if self.solvent is not None: for index, reaction_system in enumerate(self.reaction_systems): if ( - not isinstance(reaction_system, Reactor) and reaction_system.const_spc_names is not None - ): # if no constant species provided do nothing + not isinstance(reaction_system, RMSReactor) + ) and reaction_system.const_spc_names is not None: # if no constant species provided do nothing reaction_system.get_const_spc_indices(self.reaction_model.core.species) # call the function to identify indices in the solver self.initialize_reaction_threshold_and_react_flags() if self.filter_reactions and self.init_react_tuples: - self.react_init_tuples() + self.react_init_tuples(requires_rms=requires_rms) self.reaction_model.initialize_index_species_dict() self.initialize_seed_mech() + return requires_rms - def register_listeners(self): + def register_listeners(self, requires_rms=False): """ Attaches listener classes depending on the options found in the RMG input file. """ self.attach(ChemkinWriter(self.output_directory)) - self.attach(RMSWriter(self.output_directory)) + if requires_rms: + self.attach(RMSWriter(self.output_directory)) if self.generate_output_html: self.attach(OutputHTMLWriter(self.output_directory)) @@ -735,7 +748,7 @@ def register_listeners(self): if self.save_simulation_profiles: for index, reaction_system in enumerate(self.reaction_systems): - if isinstance(reaction_system, Reactor): + if requires_rms and isinstance(reaction_system, RMSReactor): typ = type(reaction_system) raise InputError(f"save_simulation_profiles=True not compatible with reactor of type {typ}") reaction_system.attach(SimulationProfileWriter(self.output_directory, index, self.reaction_model.core.species)) @@ -749,10 +762,10 @@ def execute(self, initialize=True, **kwargs): """ if initialize: - self.initialize(**kwargs) + requires_rms = self.initialize(**kwargs) # register listeners - self.register_listeners() + self.register_listeners(requires_rms=requires_rms) self.done = False @@ -779,7 +792,7 @@ def execute(self, initialize=True, **kwargs): # Update react flags if self.filter_reactions: # Run the reaction system to update threshold and react flags - if isinstance(reaction_system, Reactor): + if requires_rms and isinstance(reaction_system, RMSReactor): self.update_reaction_threshold_and_react_flags( rxn_sys_unimol_threshold=np.zeros((len(self.reaction_model.core.species),), bool), rxn_sys_bimol_threshold=np.zeros((len(self.reaction_model.core.species), len(self.reaction_model.core.species)), bool), @@ -822,6 +835,7 @@ def execute(self, initialize=True, **kwargs): unimolecular_react=self.unimolecular_react, bimolecular_react=self.bimolecular_react, trimolecular_react=self.trimolecular_react, + requires_rms=requires_rms, ) if not np.isinf(self.model_settings_list[0].thermo_tol_keep_spc_in_edge): @@ -834,7 +848,7 @@ def execute(self, initialize=True, **kwargs): ) if not np.isinf(self.model_settings_list[0].thermo_tol_keep_spc_in_edge): - self.reaction_model.thermo_filter_down(maximum_edge_species=self.model_settings_list[0].maximum_edge_species) + self.reaction_model.thermo_filter_down(maximum_edge_species=self.model_settings_list[0].maximum_edge_species, requires_rms=requires_rms) logging.info("Completed initial enlarge edge step.\n") @@ -900,7 +914,7 @@ def execute(self, initialize=True, **kwargs): prune = False try: - if isinstance(reaction_system, Reactor): + if requires_rms and isinstance(reaction_system, RMSReactor): ( terminated, resurrected, @@ -993,7 +1007,7 @@ def execute(self, initialize=True, **kwargs): # Add objects to enlarge to the core first for objectToEnlarge in objects_to_enlarge: - self.reaction_model.enlarge(objectToEnlarge) + self.reaction_model.enlarge(objectToEnlarge, requires_rms=requires_rms) if model_settings.filter_reactions: # Run a raw simulation to get updated reaction system threshold values @@ -1002,7 +1016,7 @@ def execute(self, initialize=True, **kwargs): temp_model_settings.tol_keep_in_edge = 0 if not resurrected: try: - if isinstance(reaction_system, Reactor): + if requires_rms and isinstance(reaction_system, RMSReactor): ( terminated, resurrected, @@ -1071,7 +1085,7 @@ def execute(self, initialize=True, **kwargs): skip_update=True, ) logging.warning( - "Reaction thresholds/flags for Reaction System {0} was not updated due " "to resurrection".format(index + 1) + "Reaction thresholds/flags for Reaction System {0} was not updated due to resurrection".format(index + 1) ) logging.info("") @@ -1094,13 +1108,14 @@ def execute(self, initialize=True, **kwargs): unimolecular_react=self.unimolecular_react, bimolecular_react=self.bimolecular_react, trimolecular_react=self.trimolecular_react, + requires_rms=requires_rms, ) if old_edge_size != len(self.reaction_model.edge.reactions) or old_core_size != len(self.reaction_model.core.reactions): reactor_done = False if not np.isinf(self.model_settings_list[0].thermo_tol_keep_spc_in_edge): - self.reaction_model.thermo_filter_down(maximum_edge_species=model_settings.maximum_edge_species) + self.reaction_model.thermo_filter_down(maximum_edge_species=model_settings.maximum_edge_species, requires_rms=requires_rms) max_num_spcs_hit = len(self.reaction_model.core.species) >= model_settings.max_num_species @@ -1127,6 +1142,7 @@ def execute(self, initialize=True, **kwargs): model_settings.tol_move_to_core, model_settings.maximum_edge_species, model_settings.min_species_exist_iterations_for_prune, + requires_rms=requires_rms, ) # Perform garbage collection after pruning collected = gc.collect() @@ -1248,8 +1264,9 @@ def run_uncertainty_analysis(self): ) self.uncertainty["global"] = False else: - import re import random + import re + from rmgpy.tools.canteramodel import Cantera from rmgpy.tools.globaluncertainty import ReactorPCEFactory @@ -1894,7 +1911,7 @@ def initialize_reaction_threshold_and_react_flags(self): if self.trimolecular: self.trimolecular_react[:num_restart_spcs, :num_restart_spcs, :num_restart_spcs] = False - def react_init_tuples(self): + def react_init_tuples(self, requires_rms=False): """ Reacts tuples given in the react block """ @@ -1927,6 +1944,7 @@ def react_init_tuples(self): unimolecular_react=self.unimolecular_react, bimolecular_react=self.bimolecular_react, trimolecular_react=self.trimolecular_react, + requires_rms=requires_rms, ) def update_reaction_threshold_and_react_flags( @@ -2221,7 +2239,7 @@ def __init__(self, reaction_system, bspc): if isinstance(value, list): self.Ranges[key] = [v.value_si for v in value] - if isinstance(reaction_system, Reactor): + if isinstance(reaction_system, RMSReactor): self.tmax = reaction_system.tf else: for term in reaction_system.termination: @@ -2439,7 +2457,16 @@ def make_profile_graph(stats_file, force_graph_generation=False): if display_found or force_graph_generation: try: - from gprof2dot import PstatsParser, DotWriter, SAMPLES, themes, TIME, TIME_RATIO, TOTAL_TIME, TOTAL_TIME_RATIO + from gprof2dot import ( + SAMPLES, + TIME, + TIME_RATIO, + TOTAL_TIME, + TOTAL_TIME_RATIO, + DotWriter, + PstatsParser, + themes, + ) except ImportError: logging.warning("Trouble importing from package gprof2dot. Unable to create a graph of the profile statistics.") logging.warning("Try getting the latest version with something like `pip install --upgrade gprof2dot`.") diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 1ab45dc67ca..ca67effe89c 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -57,7 +57,8 @@ from rmgpy.species import Species from rmgpy.thermo.thermoengine import submit from rmgpy.rmg.decay import decay_species -from rmgpy.rmg.reactors import PhaseSystem, Phase, Interface, Reactor +from rmgpy.rmg.reactionmechanismsimulator_reactors import PhaseSystem, Phase, Interface, NO_JULIA +from rmgpy.rmg.reactionmechanismsimulator_reactors import Reactor as RMSReactor from rmgpy.molecule.fragment import Fragment ################################################################################ @@ -75,7 +76,7 @@ def __init__(self, species=None, reactions=None, phases=None, interfaces={}): if phases is None: phases = {"Default": Phase(), "Surface": Phase()} interfaces = {frozenset({"Default", "Surface"}): Interface(list(phases.values()))} - self.phase_system = PhaseSystem(phases, interfaces) + self.phase_system = None if NO_JULIA else PhaseSystem(phases, interfaces) def __reduce__(self): """ @@ -349,9 +350,10 @@ def make_new_species(self, object, label="", reactive=True, check_existing=True, orilabel = spec.label label = orilabel i = 2 - while any([label in phase.names for phase in self.edge.phase_system.phases.values()]): - label = orilabel + "-" + str(i) - i += 1 + if self.edge.phase_system: # None when RMS not installed + while any([label in phase.names for phase in self.edge.phase_system.phases.values()]): + label = orilabel + "-" + str(i) + i += 1 spec.label = label logging.debug("Creating new species %s", spec.label) @@ -589,7 +591,7 @@ def make_new_pdep_reaction(self, forward): return forward - def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bimolecular_react=None, trimolecular_react=None): + def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bimolecular_react=None, trimolecular_react=None, requires_rms=False): """ Enlarge a reaction model by processing the objects in the list `new_object`. If `new_object` is a @@ -636,7 +638,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi # Add new species if new_species not in self.core.species: - reactions_moved_from_edge = self.add_species_to_core(new_species) + reactions_moved_from_edge = self.add_species_to_core(new_species, requires_rms=requires_rms) else: reactions_moved_from_edge = [] @@ -644,7 +646,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi pdep_network, new_species = new_object new_reactions.extend(pdep_network.explore_isomer(new_species)) - self.process_new_reactions(new_reactions, new_species, pdep_network) + self.process_new_reactions(new_reactions, new_species, pdep_network, requires_rms=requires_rms) else: raise TypeError( @@ -668,7 +670,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi if len(products) == 1 and products[0] == species: new_reactions = network.explore_isomer(species) - self.process_new_reactions(new_reactions, species, network) + self.process_new_reactions(new_reactions, species, network, requires_rms=requires_rms) network.update_configurations(self) index = 0 break @@ -693,7 +695,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi # Identify a core species which was used to generate the reaction # This is only used to determine the reaction direction for processing spc = spcTuple[0] - self.process_new_reactions(rxnList, spc) + self.process_new_reactions(rxnList, spc, requires_rms=requires_rms) ################################################################ # Begin processing the new species and reactions @@ -705,7 +707,7 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi # Do thermodynamic filtering if not np.isinf(self.thermo_tol_keep_spc_in_edge) and self.new_species_list != []: - self.thermo_filter_species(self.new_species_list) + self.thermo_filter_species(self.new_species_list, requires_rms=requires_rms) # Update unimolecular (pressure dependent) reaction networks if self.pressure_dependence: @@ -802,7 +804,7 @@ def clear_surface_adjustments(self): self.new_surface_spcs_loss = set() self.new_surface_rxns_loss = set() - def process_new_reactions(self, new_reactions, new_species, pdep_network=None, generate_thermo=True, generate_kinetics=True): + def process_new_reactions(self, new_reactions, new_species, pdep_network=None, generate_thermo=True, generate_kinetics=True, requires_rms=False): """ Process a list of newly-generated reactions involving the new core species or explored isomer `new_species` in network `pdep_network`. @@ -824,12 +826,12 @@ def process_new_reactions(self, new_reactions, new_species, pdep_network=None, g if spec not in self.core.species: all_species_in_core = False if spec not in self.edge.species: - self.add_species_to_edge(spec) + self.add_species_to_edge(spec, requires_rms=requires_rms) for spec in rxn.products: if spec not in self.core.species: all_species_in_core = False if spec not in self.edge.species: - self.add_species_to_edge(spec) + self.add_species_to_edge(spec, requires_rms=requires_rms) isomer_atoms = sum([len(spec.molecule[0].atoms) for spec in rxn.reactants]) @@ -857,9 +859,9 @@ def process_new_reactions(self, new_reactions, new_species, pdep_network=None, g # The reaction is not new, so it should already be in the core or edge continue if all_species_in_core: - self.add_reaction_to_core(rxn) + self.add_reaction_to_core(rxn, requires_rms=requires_rms) else: - self.add_reaction_to_edge(rxn) + self.add_reaction_to_edge(rxn, requires_rms=requires_rms) else: # Add the reaction to the appropriate unimolecular reaction network # If pdep_network is not None then that will be the network the @@ -1102,7 +1104,7 @@ def log_enlarge_summary( logging.info(" The model edge has {0:d} species and {1:d} reactions".format(edge_species_count, edge_reaction_count)) logging.info("") - def add_species_to_core(self, spec): + def add_species_to_core(self, spec, requires_rms=False): """ Add a species `spec` to the reaction model core (and remove from edge if necessary). This function also moves any reactions in the edge that gain @@ -1120,7 +1122,7 @@ def add_species_to_core(self, spec): if spec in self.edge.species: # remove forbidden species from edge logging.info("Species {0} was Forbidden and not added to Core...Removing from Edge.".format(spec)) - self.remove_species_from_edge(self.reaction_systems, spec) + self.remove_species_from_edge(self.reaction_systems, spec, requires_rms=requires_rms) # remove any empty pdep networks as a result of species removal if self.pressure_dependence: self.remove_empty_pdep_networks() @@ -1132,7 +1134,8 @@ def add_species_to_core(self, spec): rxn_list = [] if spec in self.edge.species: - self.edge.phase_system.pass_species(spec.label, self.core.phase_system) + if requires_rms: + self.edge.phase_system.pass_species(spec.label, self.core.phase_system) # If species was in edge, remove it logging.debug("Removing species %s from edge.", spec) self.edge.species.remove(spec) @@ -1152,31 +1155,26 @@ def add_species_to_core(self, spec): # Move any identified reactions to the core for rxn in rxn_list: - self.add_reaction_to_core(rxn) + self.add_reaction_to_core(rxn, requires_rms=requires_rms) logging.debug("Moving reaction from edge to core: %s", rxn) - else: - if spec.molecule[0].contains_surface_site(): - self.core.phase_system.phases["Surface"].add_species(spec, edge_phase=self.edge.phase_system.phases["Surface"]) - self.edge.phase_system.species_dict[spec.label] = spec - self.core.phase_system.species_dict[spec.label] = spec - else: - self.core.phase_system.phases["Default"].add_species(spec, edge_phase=self.edge.phase_system.phases["Default"]) - self.edge.phase_system.species_dict[spec.label] = spec - self.core.phase_system.species_dict[spec.label] = spec + elif requires_rms: + destination_phase = "Surface" if spec.molecule[0].contains_surface_site() else "Default" + self.core.phase_system.phases[destination_phase].add_species(spec, edge_phase=self.edge.phase_system.phases[destination_phase]) + self.edge.phase_system.species_dict[spec.label] = spec + self.core.phase_system.species_dict[spec.label] = spec return rxn_list - def add_species_to_edge(self, spec): + def add_species_to_edge(self, spec, requires_rms=False): """ - Add a species `spec` to the reaction model edge. + Add a species `spec` to the reaction model edge and optionally the RMS phase. """ self.edge.species.append(spec) - if spec.molecule[0].contains_surface_site(): - self.edge.phase_system.phases["Surface"].add_species(spec) - self.edge.phase_system.species_dict[spec.label] = spec - else: - self.edge.phase_system.phases["Default"].add_species(spec) - self.edge.phase_system.species_dict[spec.label] = spec + if not requires_rms: + return + destination_phase = "Surface" if spec.molecule[0].contains_surface_site() else "Default" + self.edge.phase_system.phases[destination_phase].add_species(spec) + self.edge.phase_system.species_dict[spec.label] = spec def set_thermodynamic_filtering_parameters( self, Tmax, thermo_tol_keep_spc_in_edge, min_core_size_for_prune, maximum_edge_species, reaction_systems @@ -1200,7 +1198,7 @@ def set_thermodynamic_filtering_parameters( self.reaction_systems = reaction_systems self.maximum_edge_species = maximum_edge_species - def thermo_filter_species(self, spcs): + def thermo_filter_species(self, spcs, requires_rms=False): """ checks Gibbs energy of the species in species against the maximum allowed Gibbs energy @@ -1215,13 +1213,13 @@ def thermo_filter_species(self, spcs): "greater than the thermo_tol_keep_spc_in_edge of " "{3} ".format(spc, G, Gn, self.thermo_tol_keep_spc_in_edge) ) - self.remove_species_from_edge(self.reaction_systems, spc) + self.remove_species_from_edge(self.reaction_systems, spc, requires_rms=requires_rms) # Delete any networks that became empty as a result of pruning if self.pressure_dependence: self.remove_empty_pdep_networks() - def thermo_filter_down(self, maximum_edge_species, min_species_exist_iterations_for_prune=0): + def thermo_filter_down(self, maximum_edge_species, min_species_exist_iterations_for_prune=0, requires_rms=False): """ removes species from the edge based on their Gibbs energy until maximum_edge_species is reached under the constraint that all removed species are older than @@ -1263,7 +1261,7 @@ def thermo_filter_down(self, maximum_edge_species, min_species_exist_iterations_ logging.info( "Removing species {0} from edge to meet maximum number of edge species, Gibbs " "number is {1}".format(spc, Gns[rInds[i]]) ) - self.remove_species_from_edge(self.reaction_systems, spc) + self.remove_species_from_edge(self.reaction_systems, spc, requires_rms=requires_rms) # Delete any networks that became empty as a result of pruning if self.pressure_dependence: @@ -1293,7 +1291,7 @@ def remove_empty_pdep_networks(self): del self.network_dict[source] self.network_list.remove(network) - def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_edge_species, min_species_exist_iterations_for_prune): + def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_edge_species, min_species_exist_iterations_for_prune, requires_rms=False): """ Remove species from the model edge based on the simulation results from the list of `reaction_systems`. @@ -1373,7 +1371,7 @@ def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_ed for index, spec in species_to_prune[0:prune_due_to_rate_counter]: logging.info("Pruning species %s", spec) logging.debug(" %-56s %10.4e", spec, max_edge_species_rate_ratios[index]) - self.remove_species_from_edge(reaction_systems, spec) + self.remove_species_from_edge(reaction_systems, spec, requires_rms=requires_rms) if len(species_to_prune) - prune_due_to_rate_counter > 0: logging.info( "Pruning %d species to obtain an edge size of %d species", len(species_to_prune) - prune_due_to_rate_counter, maximum_edge_species @@ -1381,7 +1379,7 @@ def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_ed for index, spec in species_to_prune[prune_due_to_rate_counter:]: logging.info("Pruning species %s", spec) logging.debug(" %-56s %10.4e", spec, max_edge_species_rate_ratios[index]) - self.remove_species_from_edge(reaction_systems, spec) + self.remove_species_from_edge(reaction_systems, spec, requires_rms=requires_rms) # Delete any networks that became empty as a result of pruning if self.pressure_dependence: @@ -1389,7 +1387,7 @@ def prune(self, reaction_systems, tol_keep_in_edge, tol_move_to_core, maximum_ed logging.info("") - def remove_species_from_edge(self, reaction_systems, spec): + def remove_species_from_edge(self, reaction_systems, spec, requires_rms=False): """ Remove species `spec` from the reaction model edge. """ @@ -1397,11 +1395,12 @@ def remove_species_from_edge(self, reaction_systems, spec): # remove the species self.edge.species.remove(spec) self.index_species_dict.pop(spec.index) - self.edge.phase_system.remove_species(spec) + if requires_rms: + self.edge.phase_system.remove_species(spec) # clean up species references in reaction_systems for reaction_system in reaction_systems: - if not isinstance(reaction_system, Reactor): + if not requires_rms or not isinstance(reaction_system, RMSReactor): try: reaction_system.species_index.pop(spec) except KeyError: @@ -1473,7 +1472,7 @@ def remove_species_from_edge(self, reaction_systems, spec): self.species_cache.remove(spec) self.species_cache.append(None) - def add_reaction_to_core(self, rxn): + def add_reaction_to_core(self, rxn, requires_rms=False): """ Add a reaction `rxn` to the reaction model core (and remove from edge if necessary). This function assumes `rxn` has already been checked to @@ -1482,7 +1481,8 @@ def add_reaction_to_core(self, rxn): """ if rxn not in self.core.reactions: self.core.reactions.append(rxn) - if rxn not in self.edge.reactions: + + if requires_rms and rxn not in self.edge.reactions: # If a reaction is not in edge but is going to add to core, it is either a seed mechanism or a newly generated reaction where all reactants and products are already in core # If the reaction is in edge, then the corresponding rms_rxn was moved from edge phase to core phase in pass_species already. rms_species_list = self.core.phase_system.get_rms_species_list() @@ -1499,7 +1499,7 @@ def add_reaction_to_core(self, rxn): if rxn in self.edge.reactions: self.edge.reactions.remove(rxn) - def add_reaction_to_edge(self, rxn): + def add_reaction_to_edge(self, rxn, requires_rms=False): """ Add a reaction `rxn` to the reaction model edge. This function assumes `rxn` has already been checked to ensure it is supposed to be an edge @@ -1508,6 +1508,8 @@ def add_reaction_to_edge(self, rxn): edge). """ self.edge.reactions.append(rxn) + if not requires_rms: + return rms_species_list = self.edge.phase_system.get_rms_species_list() species_names = self.edge.phase_system.get_species_names() bits = np.array([spc.molecule[0].contains_surface_site() for spc in rxn.reactants + rxn.products]) @@ -1564,7 +1566,7 @@ def get_stoichiometry_matrix(self): stoichiometry[i, j] = nu return stoichiometry.tocsr() - def add_seed_mechanism_to_core(self, seed_mechanism, react=False): + def add_seed_mechanism_to_core(self, seed_mechanism, react=False, requires_rms=False): """ Add all species and reactions from `seed_mechanism`, a :class:`KineticsPrimaryDatabase` object, to the model core. If `react` @@ -1636,9 +1638,9 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False): # This unimolecular library reaction is flagged as `elementary_high_p` and has Arrhenius type kinetics. # We should calculate a pressure-dependent rate for it if len(rxn.reactants) == 1: - self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0]) + self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0], requires_rms=requires_rms) else: - self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0]) + self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0], requires_rms=requires_rms) # Perform species constraints and forbidden species checks @@ -1675,7 +1677,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False): spec.get_liquid_volumetric_mass_transfer_coefficient_data() spec.get_henry_law_constant_data() - self.add_species_to_core(spec) + self.add_species_to_core(spec, requires_rms=requires_rms) for rxn in self.new_reaction_list: if self.pressure_dependence and rxn.is_unimolecular(): @@ -1686,7 +1688,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False): submit(spec, self.solvent_name) rxn.fix_barrier_height(force_positive=True) - self.add_reaction_to_core(rxn) + self.add_reaction_to_core(rxn, requires_rms=requires_rms) # Check we didn't introduce unmarked duplicates self.mark_chemkin_duplicates() @@ -1698,7 +1700,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False): new_edge_reactions=[], ) - def add_reaction_library_to_edge(self, reaction_library): + def add_reaction_library_to_edge(self, reaction_library, requires_rms=False): """ Add all species and reactions from `reaction_library`, a :class:`KineticsPrimaryDatabase` object, to the model edge. @@ -1763,9 +1765,9 @@ def add_reaction_library_to_edge(self, reaction_library): # This unimolecular library reaction is flagged as `elementary_high_p` and has Arrhenius type kinetics. # We should calculate a pressure-dependent rate for it if len(rxn.reactants) == 1: - self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0]) + self.process_new_reactions(new_reactions=[rxn], new_species=rxn.reactants[0], requires_rms=requires_rms) else: - self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0]) + self.process_new_reactions(new_reactions=[rxn], new_species=rxn.products[0], requires_rms=requires_rms) # Perform species constraints and forbidden species checks for spec in self.new_species_list: @@ -1802,7 +1804,7 @@ def add_reaction_library_to_edge(self, reaction_library): spec.get_liquid_volumetric_mass_transfer_coefficient_data() spec.get_henry_law_constant_data() - self.add_species_to_edge(spec) + self.add_species_to_edge(spec, requires_rms=requires_rms) for rxn in self.new_reaction_list: if not ( @@ -1817,7 +1819,7 @@ def add_reaction_library_to_edge(self, reaction_library): ) ): # Don't add to the edge library reactions that were already processed - self.add_reaction_to_edge(rxn) + self.add_reaction_to_edge(rxn, requires_rms=requires_rms) if self.save_edge_species: from rmgpy.chemkin import mark_duplicate_reaction diff --git a/rmgpy/rmg/pdep.py b/rmgpy/rmg/pdep.py index 0e1398e51a7..834f0f72432 100644 --- a/rmgpy/rmg/pdep.py +++ b/rmgpy/rmg/pdep.py @@ -736,7 +736,7 @@ def remove_products_from_reactants(self): self.reactants.remove(prod) self.products = self.products_cache - def update(self, reaction_model, pdep_settings): + def update(self, reaction_model, pdep_settings, requires_rms=False): """ Regenerate the :math:`k(T,P)` values for this partial network if the network is marked as invalid. @@ -917,7 +917,7 @@ def update(self, reaction_model, pdep_settings): f'from the {rxn.library} library, and was not added to the model') break else: - reaction_model.add_reaction_to_core(net_reaction) + reaction_model.add_reaction_to_core(net_reaction, requires_rms=requires_rms) else: # Check whether netReaction already exists in the edge as a LibraryReaction for rxn in reaction_model.edge.reactions: @@ -929,7 +929,7 @@ def update(self, reaction_model, pdep_settings): f'from the {rxn.library} library, and was not added to the model') break else: - reaction_model.add_reaction_to_edge(net_reaction) + reaction_model.add_reaction_to_edge(net_reaction, requires_rms=requires_rms) # Set/update the net reaction kinetics using interpolation model kdata = K[:, :, i, j].copy() diff --git a/rmgpy/rmg/reactors.py b/rmgpy/rmg/reactionmechanismsimulator_reactors.py similarity index 98% rename from rmgpy/rmg/reactors.py rename to rmgpy/rmg/reactionmechanismsimulator_reactors.py index 29c17095b09..252fa90523d 100644 --- a/rmgpy/rmg/reactors.py +++ b/rmgpy/rmg/reactionmechanismsimulator_reactors.py @@ -35,12 +35,6 @@ import logging import itertools -import juliacall -from juliacall import Main -Main.seval("using PythonCall") -Main.seval("using ReactionMechanismSimulator") -Main.seval("using ReactionMechanismSimulator.Sundials") - from rmgpy.species import Species from rmgpy.molecule.fragment import Fragment from rmgpy.reaction import Reaction @@ -57,6 +51,17 @@ from rmgpy.data.kinetics.family import TemplateReaction from rmgpy.data.kinetics.depository import DepositoryReaction +NO_JULIA = False +try: + import juliacall + from juliacall import Main + Main.seval("using PythonCall") + Main.seval("using ReactionMechanismSimulator") + Main.seval("using ReactionMechanismSimulator.Sundials") +except Exception as e: + logging.info("Unable to import Julia dependencies, original error: " + str(e) + ". RMS features will not be available on this execution.") + NO_JULIA = True + def to_julia(obj): """ @@ -76,7 +81,7 @@ def to_julia(obj): if isinstance(obj, dict): return Main.PythonCall.pyconvert(Main.Dict, obj) elif isinstance(obj, (list, np.ndarray)): - if obj.getattr("shape", False) and len(obj.shape) > 1: + if getattr(obj, "shape", False) and len(obj.shape) > 1: return Main.PythonCall.pyconvert(Main.Matrix, obj) return Main.PythonCall.pyconvert(Main.Vector, obj) else: # Other native Python project does not need special conversion. diff --git a/rmgpy/tools/fluxdiagram.py b/rmgpy/tools/fluxdiagram.py index 8a472046e06..81ecf44f89a 100644 --- a/rmgpy/tools/fluxdiagram.py +++ b/rmgpy/tools/fluxdiagram.py @@ -42,7 +42,7 @@ from rmgpy.kinetics.diffusionLimited import diffusion_limiter from rmgpy.rmg.settings import SimulatorSettings -from rmgpy.solver.base import TerminationTime, TerminationConversion +from rmgpy.solver.base import TerminationConversion, TerminationTime from rmgpy.solver.liquid import LiquidReactor from rmgpy.tools.loader import load_rmg_job @@ -308,7 +308,10 @@ def generate_flux_diagram(reaction_model, times, concentrations, reaction_rates, '-pix_fmt', 'yuv420p', # Pixel format 'flux_diagram.avi'] # Output filename - subprocess.check_call(command, cwd=output_directory) + try: + subprocess.run(command, check=True, capture_output=True, cwd=output_directory) + except subprocess.CalledProcessError as err: + raise RuntimeError(f"{err} {err.stderr.decode('utf8')} {err.stdout.decode('utf8')}") from err ################################################################################ diff --git a/test/arkane/encorr/isodesmicTest.py b/test/arkane/encorr/isodesmicTest.py index 83e36e06de1..02e5684f440 100644 --- a/test/arkane/encorr/isodesmicTest.py +++ b/test/arkane/encorr/isodesmicTest.py @@ -31,7 +31,6 @@ This script contains unit tests of the :mod:`arkane.isodesmic` module. """ - import numpy as np from rmgpy.molecule import Molecule @@ -137,12 +136,11 @@ def test_initializing_constraint_map(self): """ Test that the constraint map is properly initialized when a SpeciesConstraints object is initialized """ - caffeine_consts = SpeciesConstraints(self.caffeine, [self.butane, self.benzene]) - assert set(caffeine_consts.constraint_map.keys()) == { - "H", - "C", - "O", - "N", + consts = SpeciesConstraints(self.caffeine, [self.butane, self.benzene]) + caffeine_features = consts._get_all_constraints(self.caffeine) + caffeine_constraint_list = [feat.__repr__() for feat in caffeine_features] + + assert set(caffeine_constraint_list) == { "C=O", "C-N", "C-H", @@ -154,55 +152,69 @@ def test_initializing_constraint_map(self): } no_rings = SpeciesConstraints(self.caffeine, [self.butane, self.benzene], conserve_ring_size=False) - assert set(no_rings.constraint_map.keys()) == {"H", "C", "O", "N", "C=O", "C-N", "C-H", "C=C", "C=N", "C-C"} - - atoms_only = SpeciesConstraints(self.caffeine, [self.butane], conserve_ring_size=False, conserve_bonds=False) - assert set(atoms_only.constraint_map.keys()) == {"H", "C", "O", "N"} + caffeine_features = no_rings._get_all_constraints(self.caffeine) + caffeine_constraint_list = [feat.__repr__() for feat in caffeine_features] + assert set(caffeine_constraint_list) == {"C=O", "C-N", "C-H", "C=C", "C=N", "C-C"} def test_enumerating_constraints(self): """ Test that a SpeciesConstraints object can properly enumerate the constraints of a given ErrorCancelingSpecies """ spcs_consts = SpeciesConstraints(self.benzene, []) - assert set(spcs_consts.constraint_map.keys()) == {"C", "H", "C=C", "C-C", "C-H", "6_ring"} - - # Now that we have confirmed that the correct keys are present, overwrite the constraint map to set the order - spcs_consts.constraint_map = { - "H": 0, - "C": 1, - "C=C": 2, - "C-C": 3, - "C-H": 4, - "6_ring": 5, - } + benzene_features = spcs_consts._get_all_constraints(self.benzene) + benzene_constraint_list = [feat.__repr__() for feat in benzene_features] + assert set(benzene_constraint_list) == {"C=C", "C-C", "C-H", "6_ring"} + + target_constraints, _ = spcs_consts._enumerate_constraints([benzene_features]) + benzene_constraints = target_constraints assert np.array_equal( - spcs_consts._enumerate_constraints(self.propene), - np.array([6, 3, 1, 1, 6, 0]), + benzene_constraints, + np.array([1, 3, 6, 3]), # ['6_ring', 'C-C', 'C-H', 'C=C'] ) + + spcs_consts.all_reference_species = [self.propene] + propene_features = spcs_consts._get_all_constraints(self.propene) + _, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, propene_features]) + propene_constraints = reference_constraints[0] assert np.array_equal( - spcs_consts._enumerate_constraints(self.butane), - np.array([10, 4, 0, 3, 10, 0]), + propene_constraints, + np.array([0, 1, 6, 1]), # ['6_ring', 'C-C', 'C-H', 'C=C'] ) + + spcs_consts.all_reference_species = [self.butane] + butane_features = spcs_consts._get_all_constraints(self.butane) + _, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, butane_features]) + butane_constraints = reference_constraints[0] assert np.array_equal( - spcs_consts._enumerate_constraints(self.benzene), - np.array([6, 6, 3, 3, 6, 1]), + butane_constraints, + np.array([0, 3, 10, 0]), # ['6_ring', 'C-C', 'C-H', 'C=C'] ) - # Caffeine and ethyne should return None since they have features not found in benzene - assert spcs_consts._enumerate_constraints(self.caffeine) is None - assert spcs_consts._enumerate_constraints(self.ethyne) is None + # Caffeine and ethyne should return empty list since they have features not found in benzene + spcs_consts.all_reference_species = [self.caffeine] + caffeine_features = spcs_consts._get_all_constraints(self.caffeine) + _, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, caffeine_features]) + assert len(reference_constraints) == 0 + + spcs_consts.all_reference_species = [self.ethyne] + ethyne_features = spcs_consts._get_all_constraints(self.ethyne) + _, reference_constraints = spcs_consts._enumerate_constraints([benzene_features, ethyne_features]) + assert len(reference_constraints) == 0 def test_calculating_constraints(self): """ Test that a SpeciesConstraints object can properly return the target constraint vector and the constraint matrix """ spcs_consts = SpeciesConstraints(self.caffeine, [self.propene, self.butane, self.benzene, self.ethyne]) - assert set(spcs_consts.constraint_map.keys()) == { - "H", - "C", - "O", - "N", + caffeine_features = spcs_consts._get_all_constraints(self.caffeine) + propene_features = spcs_consts._get_all_constraints(self.propene) + butane_features = spcs_consts._get_all_constraints(self.butane) + benzene_features = spcs_consts._get_all_constraints(self.benzene) + ethyne_features = spcs_consts._get_all_constraints(self.ethyne) + + caffeine_feature_list = [feat.__repr__() for feat in caffeine_features] + assert set(caffeine_feature_list) == { "C=O", "C-N", "C-H", @@ -213,36 +225,20 @@ def test_calculating_constraints(self): "6_ring", } - # Now that we have confirmed that the correct keys are present, overwrite the constraint map to set the order - spcs_consts.constraint_map = { - "H": 0, - "C": 1, - "O": 2, - "N": 3, - "C=O": 4, - "C-N": 5, - "C-H": 6, - "C=C": 7, - "C=N": 8, - "C-C": 9, - "5_ring": 10, - "6_ring": 11, - } - target_consts, consts_matrix = spcs_consts.calculate_constraints() # First, test that ethyne is not included in the reference set assert spcs_consts.reference_species == [self.propene, self.butane, self.benzene] # Then, test the output of the calculation - assert np.array_equal(target_consts, np.array([10, 8, 2, 4, 2, 10, 10, 1, 1, 1, 1, 1])) + assert np.array_equal(target_consts, np.array([ 1, 1, 1, 10, 10, 1, 1, 2, 0, 8, 10, 4, 2])) # ['5_ring', '6_ring', 'C-C', 'C-H', 'C-N', 'C=C', 'C=N', 'C=O'] assert np.array_equal( consts_matrix, np.array( [ - [6, 3, 0, 0, 0, 0, 6, 1, 0, 1, 0, 0], - [10, 4, 0, 0, 0, 0, 10, 0, 0, 3, 0, 0], - [6, 6, 0, 0, 0, 0, 6, 3, 0, 3, 0, 1], + [ 0, 0, 1, 6, 0, 1, 0, 0, 0, 3, 6, 0, 0], # ['5_ring', '6_ring', 'C-C', 'C-H', 'C-N', 'C=C', 'C=N', 'C=O'] + [ 0, 0, 3, 10, 0, 0, 0, 0, 0, 4, 10, 0, 0], # ['5_ring', '6_ring', 'C-C', 'C-H', 'C-N', 'C=C', 'C=N', 'C=O'] + [ 0, 1, 3, 6, 0, 3, 0, 0, 0, 6, 6, 0, 0], # ['5_ring', '6_ring', 'C-C', 'C-H', 'C-N', 'C=C', 'C=N', 'C=O'] ] ), ) @@ -255,12 +251,6 @@ class TestErrorCancelingScheme: @classmethod def setup_class(cls): - try: - import pyomo as pyo - except ImportError: - pyo = None - cls.pyo = pyo - lot = LevelOfTheory("test") cls.propene = ErrorCancelingSpecies(Molecule(smiles="CC=C"), (100, "kJ/mol"), lot, (105, "kJ/mol")) cls.propane = ErrorCancelingSpecies(Molecule(smiles="CCC"), (75, "kJ/mol"), lot, (80, "kJ/mol")) @@ -276,10 +266,12 @@ def setup_class(cls): def test_creating_error_canceling_schemes(self): scheme = ErrorCancelingScheme( - self.propene, - [self.butane, self.benzene, self.caffeine, self.ethyne], - True, - True, + target=self.propene, + reference_set=[self.butane, self.benzene, self.caffeine, self.ethyne], + isodesmic_class="rc2", + conserve_ring_size=True, + limit_charges=True, + limit_scope=True, ) assert scheme.reference_species == [self.butane] @@ -288,26 +280,20 @@ def test_creating_error_canceling_schemes(self): assert isodesmic_scheme.reference_species == [self.butane, self.benzene] - # def test_find_error_canceling_reaction(self): - # """ - # Test that the MILP problem can be solved to find a single isodesmic reaction - # """ - # scheme = IsodesmicScheme( - # self.propene, - # [self.propane, self.butane, self.butene, self.caffeine, self.ethyne], - # ) - - # # Note that caffeine and ethyne will not be allowed, so for the full set the indices are [0, 1, 2] - # rxn, _ = scheme._find_error_canceling_reaction([0, 1, 2], milp_software=["lpsolve"]) - # assert rxn.species[self.butane] == -1 - # assert rxn.species[self.propane] == 1 - # assert rxn.species[self.butene] == 1 - - # if self.pyo is not None: - # rxn, _ = scheme._find_error_canceling_reaction([0, 1, 2], milp_software=["pyomo"]) - # assert rxn.species[self.butane] == -1 - # assert rxn.species[self.propane] == 1 - # assert rxn.species[self.butene] == 1 + def test_find_error_canceling_reaction(self): + """ + Test that the MILP problem can be solved to find a single isodesmic reaction + """ + scheme = IsodesmicScheme( + self.propene, + [self.propane, self.butane, self.butene, self.caffeine, self.ethyne], + ) + + # Note that caffeine and ethyne will not be allowed, so for the full set the indices are [0, 1, 2] + rxn, _ = scheme._find_error_canceling_reaction([0, 1, 2]) + assert rxn.species[self.butane] == -1 + assert rxn.species[self.propane] == 1 + assert rxn.species[self.butene] == 1 def test_multiple_error_canceling_reactions(self): """ @@ -328,20 +314,13 @@ def test_multiple_error_canceling_reactions(self): ) reaction_list = scheme.multiple_error_canceling_reaction_search(n_reactions_max=20) - assert len(reaction_list) == 20 + assert len(reaction_list) == 6 reaction_string = reaction_list.__repr__() # Consider both permutations of the products in the reaction string rxn_str1 = " 1*CCC + 1*C=CCC >" rxn_str2 = " 1*C=CCC + 1*CCC >" assert any(rxn_string in reaction_string for rxn_string in [rxn_str1, rxn_str2]) - if self.pyo is not None: - # pyomo is slower, so don't test as many - reaction_list = scheme.multiple_error_canceling_reaction_search(n_reactions_max=5, milp_software=["pyomo"]) - assert len(reaction_list) == 5 - reaction_string = reaction_list.__repr__() - assert any(rxn_string in reaction_string for rxn_string in [rxn_str1, rxn_str2]) - def test_calculate_target_enthalpy(self): """ Test that ErrorCancelingScheme is able to calculate thermochemistry for the target species @@ -360,10 +339,6 @@ def test_calculate_target_enthalpy(self): ], ) - target_thermo, rxn_list = scheme.calculate_target_enthalpy(n_reactions_max=3, milp_software=["lpsolve"]) - assert target_thermo.value_si == 115000.0 + target_thermo, rxn_list = scheme.calculate_target_enthalpy(n_reactions_max=3) + assert target_thermo.value_si == 110000.0 assert isinstance(rxn_list[0], ErrorCancelingReaction) - - if self.pyo is not None: - target_thermo, _ = scheme.calculate_target_enthalpy(n_reactions_max=3, milp_software=["pyomo"]) - assert target_thermo.value_si == 115000.0 diff --git a/test/arkane/ess/molproTest.py b/test/arkane/ess/molproTest.py index fcc1c27ea86..d0b1b3c5591 100644 --- a/test/arkane/ess/molproTest.py +++ b/test/arkane/ess/molproTest.py @@ -144,6 +144,14 @@ def test_load_non_f12_e0(self): e0 = molpro_log.load_energy() assert round(abs(e0 - -301585968.58196217), 7) == 0 + def test_load_rs2c_e0(self): + """ + Load E0 for RS2C from a molpro output file + """ + molpro_log = MolproLog(os.path.join(self.data_path, "HO2_RS2C.out")) + e0 = molpro_log.load_energy() + assert round(abs(e0 - -395663227.2323160), 6) == 0 + def test_load_mrci_e0(self): """ Load the MRCI and MRCI+Davidson energies from a molpro output file diff --git a/test/arkane/ess/orcaTest.py b/test/arkane/ess/orcaTest.py index bf9a0e878f9..616fa800571 100644 --- a/test/arkane/ess/orcaTest.py +++ b/test/arkane/ess/orcaTest.py @@ -65,6 +65,8 @@ def test_number_of_atoms_from_orca_log(self): Uses Orca log files to test that number of atoms can be properly read. """ + log = OrcaLog(os.path.join(self.data_path, "N_MRCI.log")) + assert log.get_number_of_atoms() == 1 log = OrcaLog(os.path.join(self.data_path, "Orca_opt_freq_test.log")) assert log.get_number_of_atoms() == 3 log = OrcaLog(os.path.join(self.data_path, "Orca_dlpno_test.log")) @@ -144,6 +146,13 @@ def test_load_modes_from_orca_log(self): assert len([mode for mode in conformer.modes if isinstance(mode, HarmonicOscillator)]) == 1 assert len(unscaled_frequencies) == 3 + log = OrcaLog(os.path.join(self.data_path, "N_MRCI.log")) + conformer, unscaled_frequencies = log.load_conformer() + assert len([mode for mode in conformer.modes if isinstance(mode, IdealGasTranslation)]) == 1 + assert len([mode for mode in conformer.modes if isinstance(mode, NonlinearRotor)]) == 0 + assert len([mode for mode in conformer.modes if isinstance(mode, HarmonicOscillator)]) == 0 + assert len(unscaled_frequencies) == 0 + def test_spin_multiplicity_from_orca_log(self): """ Test that molecular degrees of freedom can be properly read. diff --git a/test/regression/minimal_surface/input.py b/test/regression/minimal_surface/input.py index ddfea2dd872..9054a32827d 100644 --- a/test/regression/minimal_surface/input.py +++ b/test/regression/minimal_surface/input.py @@ -104,7 +104,8 @@ toleranceKeepInEdge=0.0, toleranceMoveToCore=1e-1, toleranceInterruptSimulation=0.1, - maximumEdgeSpecies=100000, + maximumEdgeSpecies=100, + maxNumSpecies=10, ) options( units='si', diff --git a/test/rmgpy/data/solvationTest.py b/test/rmgpy/data/solvationTest.py index 584eb0b7e41..eb936522241 100644 --- a/test/rmgpy/data/solvationTest.py +++ b/test/rmgpy/data/solvationTest.py @@ -29,6 +29,8 @@ import os +import pytest + from rmgpy import settings from rmgpy.data.solvation import ( DatabaseError, @@ -36,14 +38,12 @@ SolvationDatabase, SolventLibrary, get_critical_temperature, - get_liquid_saturation_density, get_gas_saturation_density, + get_liquid_saturation_density, ) -from rmgpy.molecule import Molecule -from rmgpy.rmg.main import RMG -from rmgpy.rmg.main import Species from rmgpy.exceptions import InputError -import pytest +from rmgpy.molecule import Molecule +from rmgpy.rmg.main import RMG, Species class TestSoluteDatabase: @@ -141,8 +141,8 @@ def test_saturation_density(self): """ compound_name = "Hexane" temp = 400 # in K - assert round(abs(get_liquid_saturation_density(compound_name, temp) - 6383.22), 2) == 0 - assert round(abs(get_gas_saturation_density(compound_name, temp) - 162.99), 2) == 0 + assert round(abs(get_liquid_saturation_density(compound_name, temp) - 6385.15), 2) == 0 + assert round(abs(get_gas_saturation_density(compound_name, temp) - 163.02), 2) == 0 # Unsupported compound name with pytest.raises(DatabaseError): get_gas_saturation_density("Hexadecane", temp) diff --git a/test/rmgpy/molecule/drawTest.py b/test/rmgpy/molecule/drawTest.py index 03e68a0c212..6c416f171d8 100644 --- a/test/rmgpy/molecule/drawTest.py +++ b/test/rmgpy/molecule/drawTest.py @@ -33,9 +33,9 @@ import os import os.path +import itertools - -from rmgpy.molecule import Molecule +from rmgpy.molecule import Molecule, Atom, Bond from rmgpy.molecule.draw import MoleculeDrawer from rmgpy.species import Species @@ -144,3 +144,157 @@ def test_draw_hydrogen_bond_adsorbate(self): from cairo import PDFSurface surface, _cr, (_xoff, _yoff, _width, _height) = self.drawer.draw(molecule, file_format="pdf") assert isinstance(surface, PDFSurface) + + def test_draw_bidentate_adsorbate(self): + try: + from cairocffi import ImageSurface + except ImportError: + from cairo import ImageSurface + + test_molecules = [ + Molecule().from_adjacency_list( + """ +1 C u0 p0 c0 {2,S} {4,S} {5,S} {6,S} +2 C u0 p0 c0 {1,S} {3,S} {7,S} {8,S} +3 C u0 p0 c0 {2,S} {9,S} {10,S} {11,S} +4 X u0 p0 c0 {1,S} +5 H u0 p0 c0 {1,S} +6 H u0 p0 c0 {1,S} +7 H u0 p0 c0 {2,S} +8 X u0 p0 c0 {2,S} +9 H u0 p0 c0 {3,S} +10 H u0 p0 c0 {3,S} +11 H u0 p0 c0 {3,S} + """), + Molecule().from_adjacency_list( +""" +1 C u0 p0 c0 {2,S} {4,S} {5,S} {6,S} +2 C u0 p0 c0 {1,S} {3,S} {7,S} {8,S} +3 C u0 p0 c0 {2,S} {9,S} {10,S} {11,S} +4 X u0 p0 c0 {1,S} +5 H u0 p0 c0 {1,S} +6 H u0 p0 c0 {1,S} +7 H u0 p0 c0 {2,S} +8 H u0 p0 c0 {2,S} +9 H u0 p0 c0 {3,S} +10 H u0 p0 c0 {3,S} +11 X u0 p0 c0 {3,S} +"""), + ] + for molecule in [Molecule(smiles="CC(=O)CCO"), Molecule(smiles="C1CC=CC1CC"), Molecule(smiles="C=CCC(O)=O")]: + bondable = [a for a in molecule.atoms if a.is_non_hydrogen() and any(b.is_hydrogen() for b in a.bonds)] + for b1, b2 in itertools.combinations(bondable, 2): + # find a hydrogen atom bonded to each of the two atoms + for h1 in b1.bonds: + if h1.is_hydrogen(): + break + for h2 in b2.bonds: + if h2.is_hydrogen(): + break + molecule.remove_atom(h1) + molecule.remove_atom(h2) + x1 = Atom(element='X', radical_electrons=0, charge=0, label='', lone_pairs=0) + x2 = Atom(element='X', radical_electrons=0, charge=0, label='', lone_pairs=0) + molecule.add_atom(x1) + molecule.add_atom(x2) + molecule.add_bond(Bond(b1, x1, order=1)) + molecule.add_bond(Bond(b2, x2, order=1)) + test_molecules.append(molecule.copy(deep=True)) + molecule.remove_atom(x1) + molecule.remove_atom(x2) + molecule.add_atom(h1) + molecule.add_atom(h2) + molecule.add_bond(Bond(b1, h1, order=1)) + molecule.add_bond(Bond(b2, h2, order=1)) + + for b1, b2, b3 in itertools.combinations(bondable, 3): + # find a hydrogen atom bonded to each of the two atoms + for h1 in b1.bonds: + if h1.is_hydrogen(): + break + for h2 in b2.bonds: + if h2.is_hydrogen(): + break + for h3 in b3.bonds: + if h3.is_hydrogen(): + break + molecule.remove_atom(h1) + molecule.remove_atom(h2) + molecule.remove_atom(h3) + x1 = Atom(element='X', radical_electrons=0, charge=0, label='', lone_pairs=0) + x2 = Atom(element='X', radical_electrons=0, charge=0, label='', lone_pairs=0) + x3 = Atom(element='X', radical_electrons=0, charge=0, label='', lone_pairs=0) + molecule.add_atom(x1) + molecule.add_atom(x2) + molecule.add_atom(x3) + molecule.add_bond(Bond(b1, x1, order=1)) + molecule.add_bond(Bond(b2, x2, order=1)) + molecule.add_bond(Bond(b3, x3, order=1)) + test_molecules.append(molecule.copy(deep=True)) + molecule.remove_atom(x1) + molecule.remove_atom(x2) + molecule.remove_atom(x3) + molecule.add_atom(h1) + molecule.add_atom(h2) + molecule.add_atom(h3) + molecule.add_bond(Bond(b1, h1, order=1)) + molecule.add_bond(Bond(b2, h2, order=1)) + molecule.add_bond(Bond(b3, h3, order=1)) + + test_molecules.append(Molecule().from_adjacency_list( +""" +1 C u0 p0 c0 {2,S} {3,S} {9,S} {10,S} +2 X u0 p0 c0 {1,S} +3 C u0 p0 c0 {1,S} {4,S} {5,S} {11,S} +4 X u0 p0 c0 {3,S} +5 C u0 p0 c0 {3,S} {6,S} {7,S} {12,S} +6 X u0 p0 c0 {5,S} +7 C u0 p0 c0 {5,S} {8,S} {13,S} {14,S} +8 X u0 p0 c0 {7,S} +9 H u0 p0 c0 {1,S} +10 H u0 p0 c0 {1,S} +11 H u0 p0 c0 {3,S} +12 H u0 p0 c0 {5,S} +13 H u0 p0 c0 {7,S} +14 H u0 p0 c0 {7,S} +""")) + test_molecules.append(Molecule().from_adjacency_list( +""" +1 O u0 p2 c0 {4,S} {5,S} +2 O u0 p2 c0 {5,S} {6,S} +3 O u0 p2 c0 {6,S} {26,S} +4 C u0 p0 c0 {1,S} {7,S} {8,S} {19,S} +5 C u0 p0 c0 {1,S} {2,S} {10,S} {21,S} +6 C u0 p0 c0 {2,S} {3,S} {9,S} {20,S} +7 C u0 p0 c0 {4,S} {15,S} {16,S} {24,S} +8 C u0 p0 c0 {4,S} {17,S} {18,S} {25,S} +9 C u0 p0 c0 {6,S} {11,S} {12,S} {22,S} +10 C u0 p0 c0 {5,S} {13,S} {14,S} {23,S} +11 H u0 p0 c0 {9,S} +12 H u0 p0 c0 {9,S} +13 H u0 p0 c0 {10,S} +14 H u0 p0 c0 {10,S} +15 H u0 p0 c0 {7,S} +16 H u0 p0 c0 {7,S} +17 H u0 p0 c0 {8,S} +18 H u0 p0 c0 {8,S} +19 X u0 p0 c0 {4,S} +20 X u0 p0 c0 {6,S} +21 X u0 p0 c0 {5,S} +22 X u0 p0 c0 {9,S} +23 X u0 p0 c0 {10,S} +24 X u0 p0 c0 {7,S} +25 X u0 p0 c0 {8,S} +26 X u0 p0 c0 {3,S} +""")) + test_molecules.append(Molecule(smiles="*CC(*)(C*)OCC#*")) + test_molecules.append(Molecule(smiles="*CC(*)(C*)C*")) + for number, molecule in enumerate(test_molecules, 1): + path = f"test_polydentate_{number}.png" + if os.path.exists(path): + os.unlink(path) + self.drawer.clear() + surface, _cr, (_xoff, _yoff, width, height) = self.drawer.draw(molecule, file_format="png", target=path) + assert os.path.exists(path), "File doesn't exist" + os.unlink(path) + assert isinstance(surface, ImageSurface) \ No newline at end of file diff --git a/test/rmgpy/molecule/resonanceTest.py b/test/rmgpy/molecule/resonanceTest.py index 2130307e25c..8c620ff8a35 100644 --- a/test/rmgpy/molecule/resonanceTest.py +++ b/test/rmgpy/molecule/resonanceTest.py @@ -31,7 +31,6 @@ from rmgpy.molecule.molecule import Molecule from rmgpy.molecule.resonance import ( _clar_optimization, - _clar_transformation, generate_clar_structures, generate_kekule_structure, generate_optimal_aromatic_resonance_structures, @@ -1369,15 +1368,6 @@ class ClarTest: Contains unit tests for Clar structure methods. """ - def test_clar_transformation(self): - """Test that clarTransformation generates an aromatic ring.""" - mol = Molecule().from_smiles("c1ccccc1") - sssr = mol.get_smallest_set_of_smallest_rings() - _clar_transformation(mol, sssr[0]) - mol.update_atomtypes() - - assert mol.is_aromatic() - def test_clar_optimization(self): """Test to ensure pi electrons are conserved during optimization""" mol = Molecule().from_smiles("C1=CC=C2C=CC=CC2=C1") # Naphthalene diff --git a/test/rmgpy/rmg/mainTest.py b/test/rmgpy/rmg/mainTest.py index cea42c859cb..2a68efc73d7 100644 --- a/test/rmgpy/rmg/mainTest.py +++ b/test/rmgpy/rmg/mainTest.py @@ -32,15 +32,12 @@ import shutil from unittest.mock import patch -import pytest - import pandas as pd +import pytest -from rmgpy.rmg.main import RMG, initialize_log, make_profile_graph -from rmgpy.rmg.main import RMG_Memory -from rmgpy import get_path -from rmgpy import settings +from rmgpy import get_path, settings from rmgpy.data.rmg import RMGDatabase +from rmgpy.rmg.main import RMG, RMG_Memory, initialize_log, make_profile_graph from rmgpy.rmg.model import CoreEdgeReactionModel originalPath = get_path() @@ -195,7 +192,10 @@ def setup_class(cls): cls.outputDir = os.path.join(cls.testDir, "output_w_filters") cls.databaseDirectory = settings["database.directory"] - os.mkdir(cls.outputDir) + try: + os.mkdir(cls.outputDir) + except FileExistsError: + pass # output directory already exists initialize_log(logging.INFO, os.path.join(cls.outputDir, "RMG.log")) cls.rmg = RMG( diff --git a/test/rmgpy/solver/liquidTest.py b/test/rmgpy/solver/liquidTest.py index dd03ce323b9..3559a595b36 100644 --- a/test/rmgpy/solver/liquidTest.py +++ b/test/rmgpy/solver/liquidTest.py @@ -29,7 +29,6 @@ import os - import numpy as np from rmgpy.kinetics import Arrhenius @@ -635,5 +634,7 @@ def teardown_class(cls): import rmgpy.data.rmg rmgpy.data.rmg.database = None - - os.remove(os.path.join(cls.file_dir, "restart_from_seed.py")) + try: + os.remove(os.path.join(cls.file_dir, "restart_from_seed.py")) + except FileNotFoundError: + pass # file will not be present if any tests failed diff --git a/test/rmgpy/statmech/ndTorsionsTest.py b/test/rmgpy/statmech/ndTorsionsTest.py index 80186b1f3e3..5df907f074b 100644 --- a/test/rmgpy/statmech/ndTorsionsTest.py +++ b/test/rmgpy/statmech/ndTorsionsTest.py @@ -109,4 +109,4 @@ def test_hindered_rotor_nd(self): hdnd.read_scan() assert round(abs(hdnd.Es[0]) - 8.58538448, 4) == 0 hdnd.fit() - assert round(abs(hdnd.calc_partition_function(300.0)) - 2.899287634962152, 5) == 0 + assert round(abs(hdnd.calc_partition_function(300.0)) - 2.899287634962152, 4) == 0 diff --git a/utilities.py b/utilities.py index f84963a6324..87d39e73845 100644 --- a/utilities.py +++ b/utilities.py @@ -48,7 +48,6 @@ def check_dependencies(): print('{0:<15}{1:<15}{2}'.format('Package', 'Version', 'Location')) missing = { - 'lpsolve': _check_lpsolve(), 'openbabel': _check_openbabel(), 'pydqed': _check_pydqed(), 'pyrdl': _check_pyrdl(), @@ -80,23 +79,6 @@ def check_dependencies(): """) -def _check_lpsolve(): - """Check for lpsolve""" - missing = False - - try: - import lpsolve55 - except ImportError: - print('{0:<30}{1}'.format('lpsolve55', - 'Not found. Necessary for generating Clar structures for aromatic species.')) - missing = True - else: - location = lpsolve55.__file__ - print('{0:<30}{1}'.format('lpsolve55', location)) - - return missing - - def _check_openbabel(): """Check for OpenBabel""" missing = False