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: AMS2027 or later

Initial Imports¶

from scm.base import ChemicalSystem
from scm.inputs import ADF, AMS, BAND, DFTB, MLPotential, ReaxFF
from scm.plams import AMSJob, 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())

def read_chemical_systems(folder):
    return {
        path.stem: ChemicalSystem.from_xyz(str(path))
        for path in sorted(Path(folder).glob("*.xyz"))
    }


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

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

s_ani2x = MLPotential(Backend="TorchANI", Model="ANI-2x")

s_ani1ccx = MLPotential(Backend="TorchANI", Model="ANI-1ccx")

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

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

s_adf = ADF()
s_adf.Basis.Type = "TZP"
s_adf.Basis.Core = "None"
s_adf.XC.LibXC = "WB97X"

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

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 = AMS(Task="SinglePoint")
# s_ams.Task = 'GeometryOptimization'
# s_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 = 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)

[12.10|16:44:28] JOB ANI-1ccx_1 STARTED
[12.10|16:44:28] JOB ANI-1ccx_22 STARTED
[12.10|16:44:28] JOB ANI-1ccx_2233tetramethylbutane STARTED
[12.10|16:44:28] JOB ANI-1ccx_31 STARTED
[12.10|16:44:28] JOB ANI-1ccx_adamantane STARTED
[12.10|16:44:28] JOB ANI-1ccx_bicyclo222octane STARTED
[12.10|16:44:28] JOB ANI-1ccx_ethane STARTED
[12.10|16:44:28] JOB ANI-1ccx_ethylene STARTED
[12.10|16:44:28] JOB ANI-1ccx_ethyne STARTED
[12.10|16:44:28] JOB ANI-1ccx_hexane STARTED
... output trimmed ....
[12.10|17:46:02] Waiting for job BAND_e_9 to finish
[12.10|17:46:30] JOB BAND_p_3 FINISHED
[12.10|17:46:30] JOB BAND_p_3 SUCCESSFUL
[12.10|17:48:03] JOB BAND_e_9 FINISHED
[12.10|17:48:03] JOB BAND_e_9 SUCCESSFUL
[12.10|17:48:03] Waiting for job BAND_p_10 to finish
[12.10|17:49:08] JOB BAND_p_10 FINISHED
[12.10|17:49:08] JOB BAND_p_10 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.base import ChemicalSystem
from scm.inputs import ADF, AMS, BAND, DFTB, MLPotential, ReaxFF
from scm.plams import AMSJob, 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())

def read_chemical_systems(folder):
    return {
        path.stem: ChemicalSystem.from_xyz(str(path))
        for path in sorted(Path(folder).glob("*.xyz"))
    }


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


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

s_ani2x = MLPotential(Backend="TorchANI", Model="ANI-2x")

s_ani1ccx = MLPotential(Backend="TorchANI", Model="ANI-1ccx")

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

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

s_adf = ADF()
s_adf.Basis.Type = "TZP"
s_adf.Basis.Core = "None"
s_adf.XC.LibXC = "WB97X"

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


# 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 = AMS(Task="SinglePoint")
# s_ams.Task = 'GeometryOptimization'
# s_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 = 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)

Required additional files¶

molecules_ReactionEnergyBenchmark.zip

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¶

Required additional files¶