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