Reparametrizing ReaxFF with the CMA-ES optimizer¶

Introduction¶

The tutorial is inspired by the work done by A. Vashisth and coworkers on re-fitting an existing ReaxFF force field for usage with the novel bond boost acceleration method:

A. Vashisth, C. Ashraf, W. Zhang, C. E. Bakis, and Adri C. T. van Duin, J. Phys. Chem. A 122 (2018) 6633

Like in the paper, the parameters related to C-O, C-N and N-H in the force field CHON-2017_weak.ff are going to be refitted. Both the initial force field and the refitted force field from the paper (CHON-2017_weak_bb.ff) are shipped with AMS2019, so you can compare both with the force field you are going to create in this tutorial.

The training set will contain the following types of training data:

• Geometry optimizations
• Bond scans
• Conformers
• Trajectory snapshots

Preparation¶

Download the tutorial package CMA-ES-FIT.zip and unpack its contents into a location where you have read and write permissions. On a windows computer, the ADF_DATA directory is a good location for the tutorial package.

Preparation: Set up a DFTB3 preset¶

Each of the training set entries is explained for one compound in detail within its own subsection. It’s your task to apply the same calculation for the other compounds to successively build up a training set.

To guide the process, each subsection concludes with a small checklist of training set entries that you should create before moving on to the next chapter. It might seem like a lot for just a tutorial, but keep in mind that not all entries are used for the fitting, some will be used for cross-validating the fitted force fields!

Throughout this tutorial the following level of theory will be used

• Task: Single Point
• DFTB3
• Dispersion correction: D3-BJ
• Parameterset: DFTB.org/3ob-3-1

This setup - referred to as DFTB3 from now on - marks a good compromise between computational speed and efficiency.

For fast access to these settings, save a new preset in AMSinput:

File → Preset → Save as preset…

Enter DFTB3 and press Save

Now a DFTB3 single point can be requested with just one click:

File → Preset → DFTB3

Optimized geometries¶

The optimized reference structure is added to the training set and a set of target values (bond lengths, angles, energies, etc…) has to be defined as a measure. During each(!) evaluation of the objective function, ReaxFF will perform a geometry optimization with the current force field guess. This optimization starts from the reference structure in the training set and once the geometry optimization has ended - converged or not - the target values are compared with the references in the training set. Since geometry optimizations are computationally demanding, a good approach is to limit the max. number of geometry optimization steps to a small number…- just large enough to be able to detect whether or not it diverges from the reference.

As a first step in building up a training set, add the optimized geometries of the following three epoxy resins and the DETDA hardener to the training data:

The structures can be drawn in AMSinput or just be imported from the folder xyz inside the tutorial package CMA-ES-FIT. After importing the coordinates into AMSinput, run a DFTB3 geometry optimization:

File → Import Coordinates

File → Preset → DFTB

In the DFTB panel of AMSinput:

Task → Geometry Optimization

Save and Run

Repeat the calculation for all four compounds (BisF, DETDA, TGAP and TGDDM). After the calculation has finished, open AMSjobs to create a new ReaxFF training set from this geometry optimizations:

Select the finished jobs in AMSjobs (hold down SHIFT to select multiple jobs)

Job → Add To Training Set

The four structures will be added to a new training set automatically. In the AMStrain window that pops up, the entries are marked bright red to indicate that they are selected. Geometries that have been added but are not yet used in the training set are considered an error and thus highlighted in light red.

To add all bonds and angles of the BisF molecule to the training set:

Select the BisF entry by clicking on the System field next to it

Add → Add Bonds

Add → Add Angles

Note how the BisF entry is not colored red anymore, now that we have defined a set of references for the training set from it. Repeat these steps for the remaining three compounds.

To reduce the computational cost of the four geometry optimizations, we set the maximum number of iterations of the optimizer to 10:

Double Click into one of the Details fields

Enter 10 as the Maximum iterations

Click OK

The general settings of the geometry optimizer, e.g. which optimizer will be used or whether or not the cell will be optimized, are defined at a later stage in a ReaxFF input file.

Add the following entries to your training set:

Checklist (Geo. Opt.)
---------------------
[ ] DETDA
[ ] BisF
[ ] TGAP
[ ] TGDDM

The weights assigned to the bond lengths and angles, 0.1 and 5.0, mark a good starting point and can be kept for now. Before continuing to the next chapter, save your training set

File → Save

PES/Bond Scans¶

Scanning internal coordinates is very easily done with the help of the graphical user interface of AMS2020. The AMS driver provides a potential energy surface (PES) scan, so the setup and results will always look the same, regardless of the engine that was used to calculate them. As the name suggests, multiple coordinates can be scanned at the same time to construct a multi-dimensional surface. However, for this tutorial the use of one dimensional bond scans is sufficient. If needed, PES scans can be added to your training data in exactly the same way as the bond scans.

For the generation of ReaxFF training data, scanning a bond distance is probably the most common task. For our purpose we start by scanning the bonds that we know will be broken and formed during the polymerization reaction. Rather than use whole molecules we opt to use fragments to reduce computation cost, and improve reproducibility versus large flexible molecules. At the refinement stage - when comparing the first fitted force fields - it might show that more scans are need, e.g. different bond-, angle- or dihedral scans.

Start by importing the coordinates of fragment1 from the xyz directory of the tutorial package, which is a fragment of the DETDA hardener. Then, from within the DFTB input panel of AMSinput:

Select the DFTB3 preset from File → Preset → DFTB3

Change the task to PES Scan

Click on

Hold down the SHIFT key of the keyboard to select an N-H pair

Press the “+” button next to N(11)H(21) distance

Scan from 0.8 to 1.8 A in 15 steps

Save & Run (name: fragment1-scan)

Click yes if prompted to adjust the initial structure

When setting up the bond scan calculations, we try to choose the scan parameters such that we capture enough of the anharmonicity of the curve. The gradient displayed in AMSmovie is a good indicator for this criterion. Do not fully dissociate the molecule to avoid the bond breaking problems associated with HF and KS-DFT.

Once the calculation has finished, you can add the scan points to your training data. Open AMSjobs:

Select the calculation, e.g. fragment1-scan

Job → Add to AMStrain

If an instance of AMStrain is already opened, the new entries will be appended to training set. If not, a new instance of AMStrain is opened. In the latter case, just edit your training data as described and merge the two training sets using the Add command in the File menu. You can either save the current set and add it later to the training set that contains the geometry optimizations or you load the geometry optimizations into this training set.

Per default the energies are added relative to the point of lowest energy in the scan curve or on the PES and assigned a weight of 5.0. The default run type is the same as the one of the PES scan: A restricted geometry optimization, meaning that all coordinates except the scan coordinate are held fixed.

In this tutorial we shall use a single point calculation only for the bond scans. To change the run type for all entries at the same time:

Select all new System entries (Hold down CTRL to select multiple entries)

Click into the Detail field of the first entry

Select Single Point as Task

Click OK

Continue by adding the bond scans in the below checklist to your training set. Try to come up with reasonable scan ranges, don’t compress too far and don’t scan too far:

Checklist (Bond scans)
----------------------
[ ] N-H bond in fragment 1
----------------------
[ ] C-O bond of epoxy in fragment 2
--------------------------------------------------
[ ] N-C bond between epoxy and amino group in fragment 3
---------------------------------------------------
[ ] O-H bond in fragment 4

Before continuing to the next chapter, save your training set

File → Save

Conformers¶

We can make use of the conformer functionality of the GUI to create conformers. The conformers are stored in an sdf file format that can directly be used from within AMSinput.

Begin by importing the structure of BisF from into AMSinput and generating 10 of its conformers. Switch to the conformer panel in AMSinput and request these settings

Enter 10 as the Number of conformations

Enter 0.5 as the RMS filter

Enter 0.5 as the Final RMS filter

Save as bisf_conformers and Run

Once the calculation has finished the results can be inspected in AMSmovie:

The energies displayed in AMSmovie stem from a UFF calculation and need to be recalculated using DFTB3. This can easily be done by using the sdf file from within AMSinput:

In a new AMSinput, switch to the DFTB panel

Model → Coordinates

Use Selected File from the dropdown menu

Click on the folder icon Molecule from: and select bisf_conformers.sdf

The DFTB3 calculation we are about to set up now, will be applied individually to all entries in the sdf file. To re-load the DFTB3 preset:

File → Preset → DFTB3

Save as bisf_conf_dftb3 and Run

Once the calculation has finished you can import the results into the training set much as you did with the geometry optimizations earlier:

Select the finished jobs in AMSjobs (hold down SHIFT to select multiple jobs)

Job → Add To Training Set

Before setting up the energy entries in the training data we shorten the text of the System IDs as ReaxFF doesn’t allow IDs longer than 20 characters. To remove the string _mol from all IDs:

Select all Systems

Training Set → Change SystemIDs

Type _mol and hit ENTER

Leave the next field empty and hit ENTER

Again the System entries are marked red since the ReaxFF training data has not yet been defined. For this objective function the relative energies between the conformers and the according optimized geometry are used. To automatically add an energy curve to the training set:

Select all of the conformer entries

Hold down the CTRL key and additionally select the BisF Geo Opt system entry

Add → Energy Curve

Per default the lowest energy is chosen as the reference structures, in this case this will be the optimized geometry.

Add the following entries to your training set:

Checklist (Conformers)
----------------------
[ ] 10 conformers of BisF
[ ] 10 conformers of TGAP
[ ] 10 conformers of TGDDM
[ ]  5 conformers of DETDA

Before continuing to the next chapter, save your training set

File → Save

Transition states / Trajectory snapshots¶

It is possible to optimize a transition state and include it in the training data, given that expected reactions are known beforehand. A good way to get started with transition state searches is by following the DFTB transition state and PES tutorial.

All three entries, reactants, transition state and products can then be added to the training data by exporting the according geometries from the GUI and adding the energies, charges, etc. to the training set. For the sake of simplicity, this current training set will not contain actual transition states. Though it might be worthwhile to try and optimize a transition state for the uncatalyzed reaction between one amine and an epoxide.

Instead an approximate energy profile of the reaction is used that we obtain by simply extracting some snapshots from two bond boost trajectories. To generate two suitable bond boost trajectories, just run the two input files located in the folder xyz/bond-boost-trajectories

Simply open the input file with AMSinput:

File → Open (choose .ams input file)

File → Run (or CTRL+R)

The calculation will take only a couple of seconds. Inspect the trajectories in AMSmovie, if you don’t see a reaction, re-run the calculation. Open the trajectory in AMSmovie and inspect the number of molecules to locate the moment the new cross linking bond is formed

Properties → Number of Molecules

Save approx. 10-15 frames (only saving every 5th frame or so) that capture the reactive event. Use the arrow of your keyboard to move forward or backward a single frame. For example, these snapshots could look like the following

The structure shown in AMSmovie can be added to the training set directly:

File → Add coordinates → To AMStrain (or directly via keyboard short CTRL+T)

Continue adding frames until you see the final H-transfer in these trajectory snapshots (see image above). Once all desired frames have been added to AMStrain, a DFTB3 job needs to be created and run on these structures to generate the reference data. From within AMStrain

Select new systems (hold down CTRL or SHIFT to select multiple lines)

Training Set → Generate Ref jobs…

Select an .ams input file from a previous DFTB3 calculation

Choose a directory for the new jobs

In the pop-up dialog you will be asked to select an example job. This job will be used as a pre-set. In this case this should be a DFTB3 single point calculation, so you can just select the .ams file you generated for the conformer calculation.

Switch to AMSjobs to run the reference jobs:

Select all reference jobs

Job → Run (or CTRL+R)

It will typically take a couple of seconds to calculate the references. Once the calculation has finished, switch back to AMStrain and add the Energy curve like you did with the conformers:

With the new systems selected

Add → Energy Curve

Next, add the following entries to your training set:

Checklist (Conformers)
----------------------
[ ] 10-15 snapshots from BisF-DETDA_BB
[ ] 10-15 snapshots from BisF_2-DETDA_BB

Before you continue…¶

At this stage you have created a first guess of your training set and it’s time to prepare the training set for the optimizer. Before doing so, let’s double check if all training set entries are present.

Your training set should now contain:

Reactants:

Bond scans:

Snapshots of the reaction:

It is possible to continue the tutorial with less or different training data than the above but the results will be different in this case.

Assigning weights¶

Some remarks on weights

The force field optimization algorithm requires the evaluation of an objective function. For ReaxFF the objective function takes the form:

where the sum runs over all training set entries. Each difference between the reference property (x_i,TS) and the ReaxFF value - calculated with your new force field - (x_i,ReaxFF) is weighted individually via the weightings σ_i. The weightings, σ_i, can be set automatically but will most likely need to be refined manually during the fitting process.

If you look at the weights in the above function more closely, you might spot that they can be interpreted in terms of an accuracy or rather a desired accuracy as you don’t know beforehand which accuracy the fitting process will yield in the end. When coming up with weights this view can guide your initial guesses. However, it should also be kept in mind that the relative values of the weights are important too as discussed in the below example.

Let’s assume you have a medium sized organic molecule and you assign a weight of 0.01 to all C-H and all C-O bonds. This scheme will probably not result in a very high accuracy of the C-O bond energies. Why? Simply because there will be much, much more C-H bonds in the system than C-O bonds. Since all errors are going to be summed up, even small changes in the C-H bond energies will affect the value of the objective function more than a medium change in the C-O bond. The optimizer “sees” only a single number, so it’s important to make sure that the objective function is balanced (hence the usage of the terms weight) - or - if it’s biased, it should be biased towards the entries you consider important for your system.

If this left you puzzled, why not take a look at the weights that were used in the published Co trainingset to gain some further inspiration?

Assigning weights with AMStrain

Within AMStrain the weights are referred to as Accuracy and can directly be entered into the according field. Multiple entries can be edited at the same time as follows:

Select the entries (hold down SHIFT to select multiple entries)

Edit the Accuracy

As an initial guess, assign the following weights to the training set entries:

Splitting the training data¶

The objective of this approach is to test the transferability and robustness of the newly generated parameters.

As opposed to the more sophisticated cross-validation schemes, we will just use our chemical intuition to generate a validation set that should ideally exhibit a similar but not exactly the same chemistry as the set used for training. This criteria can be met by simply moving one of the epoxy compounds, TGDDM, from the training data into the validation set alongside some randomly select entries of the other epoxies:

Instead of actually removing entries from the training set it is easier to just disable entries that should not be used when constructing the error function. Before continuing with editing, save the training set in a location outside the folder CMA-ES-FIT:

File → Save

Let’s start with generating the training set first. To remove all training set entries related the geometry optimization of TGDDM, disable the according system entry. Thus all related entries that use that system will be undefined and easy to locate in the GUI:

Type system into the filter

Select the TGDDM geometry optimization

Training Set → Disable

Remove system from the filter

Scroll through complete training set and disable all entries highlighted in light red

Remove the other TGDDM entries as well as selected conformer entries (see Table above)

• All TGDMM Conformers
• 2 random BisF conformers
• 2 random TGAP conformers
• 1 random DETDA conformer

Your training set should now contain only only entries related to BisF, TGAP and DETDA.

Now save the training set into the folder labeled training_data inside the folder CMA-ES-FIT you downloaded and extracted at the beginning.

To generate the validation set, we now need to invert the training set, i.e. all entries that have been disabled need to be enabled and vice versa. Simply select all disabled entries, enable them again, invert the selection and disable:

Training Set → Select Disabled

Training Set → Activate

Training Set → Invert Selection

Training Set → Disable

Now save the validation set into the folder labeled error_function inside the folder CMA-ES-FIT you downloaded and extracted at the beginning.

Errors and Cross-validation¶

Once the calculation has finished, take a look at the results file to find out which optimization run yielded the best fitness value. Keep in mind that the outcome of the optimization runs are non-deterministic and hence your results will differ from the below, but the overall trend should be similar.

For example, a run with 4 independent CMA-ES optimizations, yielded the following results

/home/ole/Workspace/CMA-ES-FIT/RUN1
xbestever found after 1314 evaluations, function value 1.13854e+08
/home/ole/Workspace/CMA-ES-FIT/RUN2
xbestever found after 1116 evaluations, function value 4.77308e+07
/home/ole/Workspace/CMA-ES-FIT/RUN3
xbestever found after 1314 evaluations, function value 5.85587e+07
/home/ole/Workspace/CMA-ES-FIT/RUN4
xbestever found after 1260 evaluations, function value 7.06017e+08

So in this case CMA-ES run #2 yielded the best overall fitness value. Let’s see how well this force field does for the validation test set. From inside the folder error_function run the following command to generate a fort.99 file:

./calculate_errors.run  ../RUN2/ffield_best

Next, inspect the resulting RMSD errors:

RMSD (Bond)  :  0.0673 (errors_bond.csv)
RMSD (angle) :  3.4729 (errors_angle.csv)
RMSD (RMSG)  : 33.6747 (errors_rmsg.csv)
RMSD (Energy): 10.9002 (errors_energy.csv)

These errors reflect the weighting in the training set. By assigning a weight of 0.01 to the bond scans in the training data we have put quite some emphasis on the energies and thus it is not surprising that these are reproduced and predicted fairly well. The fit results become even more obvious when comparing with the two force fields CHON2017_weak.ff (our starting guess!) and a refitted version of our starting guess, CHON2017_weak_bb.ff (refitted for usage with the bond boost method).

To calculate the objective function for these force fields, run:

./calculate_errors.run $AMSHOME/atomicdata/ForceFields/OldReaxFF/CHON2017_weak.ff

for the force field we used as a starting guess, and

./calculate_errors.run $AMSHOME/atomicdata/ForceFields/OldReaxFF/CHON2017_weak_bb.ff

for a force field that was refitted from our starting guess to proper DFT training set. If compare against the refitted force field:

When comparing the error files generated by analysis of the errors you see that the geometries seem to be much better described by the CHON2017_weak_bb.ff. This is due to the fact that we put so much emphasis on the PES scans in our training set.

Refine the training set¶

Let’s try to adjust the weights in our training data to get better geometries:

• Increase the number of steps in the geometry optimizations to 15
• Assign a weight of 0.01 to every bond distance in the training.in
• Assign a weight of 0.1 to all the energies in trainset.in

and run the CMA-ES optimizer again, followed by a calculation of the objective function for the best force field of this run. The new results should look similar to these finding

The resulting force field is already much closer to the published force field. Again, taking closer look at the error files with AMSgraphs reveals that this weighting scheme seems to result in a much more balanced fore field

Keep in mind that the original force field - of which we have re-used quite some parameters - was fitted against DFT data. In other words, the force field has been trained to reproduce DFT geometries and DFT energies. This means deviations between our DFTB force field and the ‘real’ force fields are to be expected since the comparison isn’t really fair. After all this approach serves for illustrating the concepts and workflow of parameter fitting and not the generation of a production ready force field. If you want to optimize your force field even further, try some of the following:

Training set / Optimization tweaks:

• More tests: Apply new force field to MD and geometry optimizations
• Try to locate problematic trainset entries and change (decrease) their weights
• Increase max. steps for the geometry optimizer
• Add more entries
• Run more (~10) CMA-ES runs
• Check params file in folder input, maybe more parameters need to be fitted?
• Check parameter ranges in params file: should a range be increased (compare to your fitted FF)?