Reaction energies with many different engines

#!/usr/bin/env plams
"""
Script for evaluating the HC7-11 and ISOL6 benchmark sets with different computational methods.

Modify the s_engine_dict variable below to specify which settings will be used.

A table with reaction energies will be printed to summary.txt
"""

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")

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 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()

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)


# summary.txt will contain a table similar to the below
##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)
#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
#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
#ANI-1ccx         15.8    30.7   -0.2    7.7   11.7   196.5   127.9     9.2    21.5     3.8    35.4     6.9
#ANI-2x           42.4    48.1   -1.9    5.4    8.0   238.4   156.9     7.9    20.8    -1.5    26.7     0.5
#DFTB              7.1    19.0    0.2    3.6    5.4   223.5   150.8    11.7    22.6     7.9    35.2     8.6
#ReaxFF           34.3    23.0   13.4    2.6    5.5   289.2   168.3    -9.2    24.2    70.7    30.9    89.3
#ADF              20.8    31.4   -1.9    6.7    9.9   214.5   139.7     9.5    20.5     4.7    31.0     4.8
#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
#BAND             19.2    30.2   -2.7    7.3   11.0   218.5   141.3     9.9    18.5     5.3    31.7     4.7

Note

The Amsterdam Modeling Suite requires the installation of additional Python packages to run the machine learning potential backends. You can use the following command to install all available machine learning potentials:

$ "$AMSBIN"/amspackages install mlpotentials

Note

To execute this PLAMS script:

This PLAMS script might take a while to run – around 4 hours on a modern laptop. Feel free to grab a coffee or take a break while it executes!