3.1. ReaxFF (basic): H₂O bond scan¶

def run_reference():
    # choose a reference engine. If you do not have a license for ADF, use DFTB, ForceField, or another engine for which you have the license.
    # from_smiles('O') creates an H2O molecule
    job = AMSJob.from_input(
        name="H2O_pesscan",
        molecule=from_smiles("O"),
        text_input="""
    Task PESScan
    PESScan
        # Scan O-H bond length between 0.85 and 1.35 angstrom (11 points)
        ScanCoordinate
            nPoints 11
            Distance 1 2 0.85 1.35
        End
    End
    Engine ADF
        XC
            GGA PBE
        End
    EndEngine
    #Engine ForceField
    #EndEngine
    #Engine DFTB
    #EndEngine
    """,
    )
    job.run()
    return job

def import_results():
    # If you ran the PESScan job via the graphical user interface, set
    # pesscan_job = '/path/to/ams.rkf'
    # Otherwise, the below line runs the reference job
    pesscan_job = run_reference()
    ri = ResultsImporter(settings={"units": {"energy": "kcal/mol"}})
    ri.add_singlejob(
        pesscan_job,
        task="PESScan",
        properties={
            "pes": {
                # apply a weights scheme (optional).
                # This sets the weight for the larger energies (very short or very long bond distance) to smaller numbers.
                # stdev is the standard deviation of the Gaussian in energy units (here kcal/mol)
                "weights_scheme": WeightsSchemeGaussian(normalization=11.0, stdev=20.0)
            }
        },
    )
    ri.store(backup=False, store_settings=False)

def create_parameter_interface():
    # start from the Water2017.ff force field
    parameter_interface = ReaxFFParameters("Water2017.ff", bounds_scale=1.2)
    parameter_interface.header["head"] = "Reparametrization of Water2017.ff"
    # only activate the "Standard" bond parameters for O-H
    for p in parameter_interface:
        p.is_active = (
            (p.metadata["block"] == "BND")
            and (p.metadata["category"] == "Standard")
            and (p.metadata["atoms"] == ["H", "O"] or p.metadata["atoms"] == ["O", "H"])
        )
    parameter_interface.yaml_store("parameter_interface.yaml")

Task Optimization

DataSet
	Name training_set
	Path training_set.yaml
End

LoggingInterval
   General 5
   Parameters 5
   History 5
   Flush 5
End

ExitCondition
	Type MaxTotalFunctionCalls
	MaxTotalFunctionCalls 100
End

Optimizer
	Type CMAES
	CMAES
		Sigma0 0.1
		Popsize 10
		MinSigma 0.0001
		Sampler full
		Verbose True
	End
End

#!/usr/bin/env amspython
# Usage: "$AMSBIN/amspython" optimization.py
from scm.plams import *
from scm.params import *
import os
import matplotlib.pyplot as plt


def main():
    init()  # initialize PLAMS workdir
    examples_dir = os.path.expandvars("$AMSHOME/scripting/scm/params/examples/ReaxFF_water")

    # load the optimization job from examples directory  and run
    optimization_job = ParAMSJob.from_inputfile(examples_dir + "/params.in", name="optimization")
    optimization_job.run()
    # or load a finished previously run job:
    # optimization_job = ParAMSJob.load_external('plams_workdir/optimization/results')

    print("Optimized force field stored in " + optimization_job.results.get_ffield())

    # loop over all PES scans and show plots
    # in this case there is only 1 PES Scan
    all_pesscans = optimization_job.results.get_all_pes_prediction_names()
    for scan in all_pesscans:
        x_headers, x_values, ref_values, pred_values, unit = optimization_job.results.get_pes_prediction(scan)
        plt.plot(x_values[0], ref_values, "-")  # ref. vs distance
        plt.plot(x_values[0], pred_values, "-")  # pred. vs distance
        plt.legend(["Reference", "Prediction"])
        plt.xlabel(x_headers[0])
        plt.ylabel(f"Relative energy ({unit})")
        plt.title(scan)
        plt.savefig(f"plot_{scan}_from_python.png")
        plt.show()

    finish()


if __name__ == "__main__":
    main()

#!/usr/bin/env amspython
# Usage: "$AMSBIN/amspython" optimization.py
from scm.plams import *
from scm.params import *
import os
import matplotlib.pyplot as plt


def main():
    init()  # initialize PLAMS workdir
    examples_dir = os.path.expandvars("$AMSHOME/scripting/scm/params/examples/ReaxFF_water")

    # load the optimization job from examples directory  and run
    optimization_job = ParAMSJob.from_inputfile(examples_dir + "/params.in", name="optimization")
    optimization_job.run()
    # or load a finished previously run job:
    # optimization_job = ParAMSJob.load_external('plams_workdir/optimization/results')

    print("Optimized force field stored in " + optimization_job.results.get_ffield())

    # loop over all PES scans and show plots
    # in this case there is only 1 PES Scan
    all_pesscans = optimization_job.results.get_all_pes_prediction_names()
    for scan in all_pesscans:
        x_headers, x_values, ref_values, pred_values, unit = optimization_job.results.get_pes_prediction(scan)
        plt.plot(x_values[0], ref_values, "-")  # ref. vs distance
        plt.plot(x_values[0], pred_values, "-")  # pred. vs distance
        plt.legend(["Reference", "Prediction"])
        plt.xlabel(x_headers[0])
        plt.ylabel(f"Relative energy ({unit})")
        plt.title(scan)
        plt.savefig(f"plot_{scan}_from_python.png")
        plt.show()

    finish()


if __name__ == "__main__":
    main()