Benchmarking Reaction Energies across AMS Engines¶

Run the same set of reference reactions with several AMS engines and compare how the predicted reaction energies vary across methods.

Downloads: Notebook | Script ?

Additional files required by this example: molecules_ReactionEnergyBenchmark.zip

Requires: AMS2026 or later

Initial Imports¶

from scm.plams import read_molecules, Settings, AMSJob, init, config, JobRunner, view

Summary Text File¶

Set up a results file, summary.txt, which will hold a table with reaction energies.

summary_fname = "summary.txt"
with open(summary_fname, "w", buffering=1) as summary_file:
    summary_file.write(
        "#Method HC7_1 HC7_2 HC7_3 HC7_4 HC7_5 HC7_6 HC7_7 ISOL6_1 ISOL6_2 ISOL6_3 ISOL6_4 ISOL6_5 (kcal/mol)\n"
    )
    summary_file.write(
        "Truhlar_ref 14.34 25.02 1.90 9.81 14.84 193.99 127.22 9.77 21.76 6.82 33.52 5.30\n"
    )
    summary_file.write(
        "Smith_wB97X_ref 28.77 41.09 1.75 6.26 9.30 238.83 157.65 9.32 20.80 1.03 26.43 0.40\n"
    )

Benchmark Calculation Setup¶

The molecules used in the benchmark calculations are loaded, and the settings for each engine created.

from pathlib import Path
import zipfile

# dependency: {} molecules_ReactionEnergyBenchmark.zip

archive = Path("molecules_ReactionEnergyBenchmark.zip")
if archive.exists() and not Path("HC7-11").exists():
    with zipfile.ZipFile(archive) as zf:
        zf.extractall(Path.cwd())

mol_dict = read_molecules("HC7-11")
mol_dict.update(read_molecules("ISOL6"))

s_reaxff = Settings()
s_reaxff.input.ReaxFF.ForceField = "CHON-2019.ff"

s_ani2x = Settings()
s_ani2x.input.MLPotential.Backend = "TorchANI"
s_ani2x.input.MLPotential.Model = "ANI-2x"

s_ani1ccx = Settings()
s_ani1ccx.input.MLPotential.Backend = "TorchANI"
s_ani1ccx.input.MLPotential.Model = "ANI-1ccx"

s_dftb = Settings()
s_dftb.input.DFTB.Model = "SCC-DFTB"
s_dftb.input.DFTB.ResourcesDir = "DFTB.org/mio-1-1"

s_band = Settings()
s_band.input.BAND.XC.libxc = "wb97x"
s_band.input.BAND.Basis.Type = "TZP"
s_band.input.BAND.Basis.Core = "None"

s_adf = Settings()
s_adf.input.ADF.Basis.Type = "TZP"
s_adf.input.ADF.Basis.Core = "None"
s_adf.input.ADF.XC.libxc = "WB97X"

s_mp2 = Settings()
s_mp2.input.ADF.Basis.Type = "TZ2P"
s_mp2.input.ADF.Basis.Core = "None"
s_mp2.input.ADF.XC.MP2 = ""

The engines in s_engine_dict below will be used for the calculation - you can remove the engines that you would not like to run (for example, if you do not have the necessary license).

s_engine_dict = {
    "ANI-1ccx": s_ani1ccx,
    "ANI-2x": s_ani2x,
    "DFTB": s_dftb,
    "ReaxFF": s_reaxff,
    "ADF": s_adf,
    "MP2": s_mp2,
    "BAND": s_band,
}

s_ams = Settings()
s_ams.input.ams.Task = "SinglePoint"
# s_ams.input.ams.Task = 'GeometryOptimization'
# s_ams.input.ams.GeometryOptimization.CoordinateType = 'Cartesian'

jobs = dict()

Running Benchmark Calculations¶

Benchmark calculations are configured to run in parallel with one core per job.

import multiprocessing

maxjobs = multiprocessing.cpu_count()
config.default_jobrunner = JobRunner(parallel=True, maxjobs=maxjobs)
config.job.runscript.nproc = 1

print(f"Running up to {maxjobs} jobs in parallel")

Running up to 12 jobs in parallel

Note: calculations will take some time to run, around 1-2 hours on a modern laptop.

with open(summary_fname, "a", buffering=1) as summary_file:
    for engine_name, s_engine in s_engine_dict.items():
        s = s_ams.copy() + s_engine.copy()
        jobs[engine_name] = dict()

        # call .run() for *all* jobs *before* accessing job.results.get_energy() for *any* job
        for mol_name, mol in mol_dict.items():
            jobs[engine_name][mol_name] = AMSJob(
                settings=s, molecule=mol, name=engine_name + "_" + mol_name
            )
            jobs[engine_name][mol_name].run()

    for engine_name in s_engine_dict:
        # for each engine, calculate reaction energies
        E = dict()
        for mol_name, job in jobs[engine_name].items():
            E[mol_name] = job.results.get_energy(unit="kcal/mol")
        deltaE_list = [
            ##### HC7/11 ######
            E["22"] - E["1"],
            E["31"] - E["1"],
            E["octane"] - E["2233tetramethylbutane"],
            5 * E["ethane"] - E["hexane"] - 4 * E["methane"],
            7 * E["ethane"] - E["octane"] - 6 * E["methane"],
            3 * E["ethylene"] + 2 * E["ethyne"] - E["adamantane"],
            3 * E["ethylene"] + 1 * E["ethyne"] - E["bicyclo222octane"],
            ###### start ISOL6 ######
            E["p_3"] - E["e_3"],
            E["p_9"] - E["e_9"],
            E["p_10"] - E["e_10"],
            E["p_13"] - E["e_13"],
            E["p_14"] - E["e_14"],
        ]

        out_str = engine_name
        for deltaE in deltaE_list:
            out_str += " {:.1f}".format(deltaE)
        out_str += "\n"

        print(out_str)
        summary_file.write(out_str)

[23.03|12:29:26] JOB ANI-1ccx_ethylene STARTED
[23.03|12:29:26] JOB ANI-1ccx_hexane STARTED
[23.03|12:29:26] JOB ANI-1ccx_methane STARTED
[23.03|12:29:26] JOB ANI-1ccx_2233tetramethylbutane STARTED
[23.03|12:29:26] JOB ANI-1ccx_bicyclo222octane STARTED
[23.03|12:29:26] JOB ANI-1ccx_31 STARTED
[23.03|12:29:26] JOB ANI-1ccx_adamantane STARTED
[23.03|12:29:26] JOB ANI-1ccx_ethyne STARTED
[23.03|12:29:26] JOB ANI-1ccx_22 STARTED
[23.03|12:29:26] JOB ANI-1ccx_ethylene RUNNING
... output trimmed ....
[23.03|13:17:08] JOB BAND_2233tetramethylbutane SUCCESSFUL
[23.03|13:17:18] JOB MP2_e_10 FINISHED
[23.03|13:17:18] JOB MP2_e_10 SUCCESSFUL
MP2 20.5 27.1 5.4 9.9 15.0 214.8 141.3 11.0 24.6 8.4 36.8 6.2

[23.03|13:17:18] Waiting for job BAND_p_9 to finish
[23.03|13:21:17] JOB BAND_p_9 FINISHED
[23.03|13:21:17] JOB BAND_p_9 SUCCESSFUL
BAND 19.3 30.4 -2.6 7.2 10.9 218.8 141.5 9.9 18.5 5.3 31.7 4.7

Results¶

Results can be loaded from the summary file and displayed using pandas. The pandas package is installable with amspackages.

with open(summary_fname) as summary_file:
    import pandas as pd

    df = pd.read_csv(summary_file, delim_whitespace=True).iloc[:, :-1]

df

	#Method	HC7_1	HC7_2	HC7_3	HC7_4	HC7_5	HC7_6	HC7_7	ISOL6_1	ISOL6_2	ISOL6_3	ISOL6_4	ISOL6_5
0	Truhlar_ref	14.34	25.02	1.90	9.81	14.84	193.99	127.22	9.77	21.76	6.82	33.52	5.3
1	Smith_wB97X_ref	28.77	41.09	1.75	6.26	9.30	238.83	157.65	9.32	20.80	1.03	26.43	0.4
2	ANI-1ccx	15.80	30.70	-0.20	7.70	11.70	196.50	127.90	9.20	21.50	3.80	35.40	6.9
3	ANI-2x	42.40	48.10	-1.90	5.40	8.00	238.40	156.90	7.90	20.80	-1.50	26.70	0.5
4	DFTB	7.10	19.00	0.20	3.60	5.40	223.50	150.80	11.70	22.60	7.90	35.20	8.6
5	ReaxFF	34.30	23.00	13.40	2.60	5.50	289.20	168.30	-9.20	24.20	70.70	30.90	89.3
6	ADF	20.30	30.90	-2.00	5.50	8.20	211.60	139.70	9.50	20.50	4.70	31.20	4.8
7	MP2	20.50	27.10	5.40	9.90	15.00	214.80	141.30	11.00	24.60	8.40	36.80	6.2
8	BAND	19.30	30.40	-2.60	7.20	10.90	218.80	141.50	9.90	18.50	5.30	31.70	4.7

Python Script¶

#!/usr/bin/env python
# coding: utf-8

# ## Initial Imports

from scm.plams import read_molecules, Settings, AMSJob, init, config, JobRunner, view


# ## Summary Text File

# Set up a results file, `summary.txt`, which will hold a table with reaction energies.

summary_fname = "summary.txt"
with open(summary_fname, "w", buffering=1) as summary_file:
    summary_file.write(
        "#Method HC7_1 HC7_2 HC7_3 HC7_4 HC7_5 HC7_6 HC7_7 ISOL6_1 ISOL6_2 ISOL6_3 ISOL6_4 ISOL6_5 (kcal/mol)\n"
    )
    summary_file.write("Truhlar_ref 14.34 25.02 1.90 9.81 14.84 193.99 127.22 9.77 21.76 6.82 33.52 5.30\n")
    summary_file.write("Smith_wB97X_ref 28.77 41.09 1.75 6.26 9.30 238.83 157.65 9.32 20.80 1.03 26.43 0.40\n")


# ## Benchmark Calculation Setup

# The molecules used in the benchmark calculations are loaded, and the settings for each engine created.

from pathlib import Path
import zipfile

# dependency: {} molecules_ReactionEnergyBenchmark.zip

archive = Path("molecules_ReactionEnergyBenchmark.zip")
if archive.exists() and not Path("HC7-11").exists():
    with zipfile.ZipFile(archive) as zf:
        zf.extractall(Path.cwd())

mol_dict = read_molecules("HC7-11")
mol_dict.update(read_molecules("ISOL6"))


s_reaxff = Settings()
s_reaxff.input.ReaxFF.ForceField = "CHON-2019.ff"

s_ani2x = Settings()
s_ani2x.input.MLPotential.Backend = "TorchANI"
s_ani2x.input.MLPotential.Model = "ANI-2x"

s_ani1ccx = Settings()
s_ani1ccx.input.MLPotential.Backend = "TorchANI"
s_ani1ccx.input.MLPotential.Model = "ANI-1ccx"

s_dftb = Settings()
s_dftb.input.DFTB.Model = "SCC-DFTB"
s_dftb.input.DFTB.ResourcesDir = "DFTB.org/mio-1-1"

s_band = Settings()
s_band.input.BAND.XC.libxc = "wb97x"
s_band.input.BAND.Basis.Type = "TZP"
s_band.input.BAND.Basis.Core = "None"

s_adf = Settings()
s_adf.input.ADF.Basis.Type = "TZP"
s_adf.input.ADF.Basis.Core = "None"
s_adf.input.ADF.XC.libxc = "WB97X"

s_mp2 = Settings()
s_mp2.input.ADF.Basis.Type = "TZ2P"
s_mp2.input.ADF.Basis.Core = "None"
s_mp2.input.ADF.XC.MP2 = ""


# The engines in `s_engine_dict` below will be used for the calculation - you can remove the engines that you would not like to run (for example, if you do not have the necessary license).

s_engine_dict = {
    "ANI-1ccx": s_ani1ccx,
    "ANI-2x": s_ani2x,
    "DFTB": s_dftb,
    "ReaxFF": s_reaxff,
    "ADF": s_adf,
    "MP2": s_mp2,
    "BAND": s_band,
}


s_ams = Settings()
s_ams.input.ams.Task = "SinglePoint"
# s_ams.input.ams.Task = 'GeometryOptimization'
# s_ams.input.ams.GeometryOptimization.CoordinateType = 'Cartesian'

jobs = dict()


# ## Running Benchmark Calculations

# Benchmark calculations are configured to run in parallel with one core per job.

import multiprocessing

maxjobs = multiprocessing.cpu_count()
config.default_jobrunner = JobRunner(parallel=True, maxjobs=maxjobs)
config.job.runscript.nproc = 1

print(f"Running up to {maxjobs} jobs in parallel")


# Note: calculations will take some time to run, around 1-2 hours on a modern laptop.

with open(summary_fname, "a", buffering=1) as summary_file:
    for engine_name, s_engine in s_engine_dict.items():
        s = s_ams.copy() + s_engine.copy()
        jobs[engine_name] = dict()

        # call .run() for *all* jobs *before* accessing job.results.get_energy() for *any* job
        for mol_name, mol in mol_dict.items():
            jobs[engine_name][mol_name] = AMSJob(settings=s, molecule=mol, name=engine_name + "_" + mol_name)
            jobs[engine_name][mol_name].run()

    for engine_name in s_engine_dict:
        # for each engine, calculate reaction energies
        E = dict()
        for mol_name, job in jobs[engine_name].items():
            E[mol_name] = job.results.get_energy(unit="kcal/mol")
        deltaE_list = [
            ##### HC7/11 ######
            E["22"] - E["1"],
            E["31"] - E["1"],
            E["octane"] - E["2233tetramethylbutane"],
            5 * E["ethane"] - E["hexane"] - 4 * E["methane"],
            7 * E["ethane"] - E["octane"] - 6 * E["methane"],
            3 * E["ethylene"] + 2 * E["ethyne"] - E["adamantane"],
            3 * E["ethylene"] + 1 * E["ethyne"] - E["bicyclo222octane"],
            ###### start ISOL6 ######
            E["p_3"] - E["e_3"],
            E["p_9"] - E["e_9"],
            E["p_10"] - E["e_10"],
            E["p_13"] - E["e_13"],
            E["p_14"] - E["e_14"],
        ]

        out_str = engine_name
        for deltaE in deltaE_list:
            out_str += " {:.1f}".format(deltaE)
        out_str += "\n"

        print(out_str)
        summary_file.write(out_str)


# ## Results

# Results can be loaded from the summary file and displayed using pandas. The pandas package is installable with `amspackages`.

with open(summary_fname) as summary_file:
    import pandas as pd

    df = pd.read_csv(summary_file, delim_whitespace=True).iloc[:, :-1]

print(df)

Benchmarking Reaction Energies across AMS Engines¶

Initial Imports¶

Summary Text File¶

Benchmark Calculation Setup¶

Running Benchmark Calculations¶

Results¶

See also¶

Related examples¶

Related tutorials¶

Related documentation¶

Python Script¶