Fitting a water potential

In this example we will use a brief molecular dynamics simulation of a single water molecule to fit a (in MoleCuilder terminology) compound potential, consisting of constant, Morse, and angle potentials.

Molecular dynamics simulation of water

First of all, we need to perform a molecular dynamics simulation of water.

Below you found the necessary command to perform such an MD simulation:

./molecuilder \
        --load water.pdb \
        --input md.pdb \
        --select-all-atoms \
        --molecular-dynamics \
                --steps 100 \
                --keep-bondgraph 1 \
                --output-every-step 1 \
                --deltat 0.1 \
                --order 1 \
                --DoLongrange 0 \
                --fragment-resultfile BondFragmentResults.dat \
        --save-homologies homologies.dat

The command consists of loading the file water.pdb (see attached files), preparing to store all subsequent steps in the output file md.pdb, then selecting all atoms for the dynamics simulation and performing it. There, we use 100 steps, keep the bond graph (in case the initial configuration is not equilibrated and oscillations are too strong) and do not include long range forces. Furthermore, a bond order of 1 is sufficient and a small time steps (deltat) ensures proper time integration of the dynamics.

Note especially the last command that stores the so-called homologies in a separate file. Homologies are the set of all equivalent configurations and fragments, i.e. when performing a simulation with many water molecules, then each water molecule is homologous (having the same structure, the same chemical formula) to any other water molecule. Hence, each molecule is stored under the same formula, namely H2O, in an internal container of MoleCuilder, including the resulting ab-initio energy and all interatomic distances. At the end of the dynamics simulation, many different configurations of a single water molecules are gathered under this key H2O. And each of them will have a slightly different energy; for example because the hydrogen atoms are slightly further away from the oxygen. In essence, we have many snapshots (with energies) of the oscillating water molecule.

For the curious ones, the fragment configurations and their energies are stored at the end of the analyse-fragment-results Action. This Action is contained in molecular-dynamics which is simply a so-called MakroAction​, combining various other Actions into a single one. This analysis action is responsible for placing results from the fragment calculations into various places: first everything goes into two result containers (short and long range ones), then forces (aka gradients) go to their respective atoms, energies finally go to the homology container.

This homology file homologies.dat serves then as the data set to the subsequent potential fitting.

Fitting a potential

While you can also fit a potential directly, for example as follows,

./molecuilder \
        --parse-homologies homologies.dat \
        --fit-potential \
                --fragment-charges 8 1 1 \
                --potential-charges 8 1 \
                --potential-type morse \
                --take-best-of 3 \
                --training-file training.dat

which would fit a Morse-type potential (potential-type) to the oxygen-hydrogen bond (set by the potential-charges "8 1") of the water molecule (specified by the fragment--charges "8 1 1"). The threshold for terminating the iterative fit procedure after three rounds and we write all input training data to a file training.dat.

However, it is usually better to fit a so-called compound potential. A compound potential consists of a combination of potentials -- for example a constant potential, a Morse potential, and an harmonic angle potential in our case -- and fits the whole lot together.

Fitting single potentials is sensible if you create snapshots that only oscillate in the degrees of freedom of this single potential, e.g. if only an elongation of the oxygen-hydrogen bond has been excited and the angle between the two O-H bonds does not change. However, as we have created a molecular dynamics simulation of an unknown initial configuration, we cannot ascertain that such a precondition holds. This is also the reason why take-best-of has been specified and not set-threshold as you will probably have noticed that the remaining L2 error is still quite large.

A compound potential is defined by the file water.potentials (we use the .potentials format use by the ​tremolo package), see the attached file. You can either create such a file from scratch (the format is quite simple) or you can use the above fit-potential calls followed by a save-potentials action to store the fitted potential to a potentials file. Afterwards you can simply combine multiple potential files into a single one.

Finally, we are able to fit the compound potential.

./molecuilder \
        --parse-homologies homologies.dat \
        --parse-potentials water.potentials \
        --fit-compound-potential \
                --fragment-charges 8 1 1 \
                --take-best-of 3 \
                --training-file training.dat \
        --save-potentials water_final.potentials

Note that the fitting of a compound potential will take longer than for a single potential as there are more degrees of freedom to consider: Each parameter of each potential represents one degree of freedom. Also note that the fit may fail if no sensible default value for each parameter is specified. If it fails, it is recommended to play around with the parameters a bit. Note that increasing verbosity (--verbose 2) will give you more information on how the iterative fitting procedure progresses.

After a minute or so, the iteration will have stopped and in the file water_final.potentials the resulting potential parameters are given. Note that the contents of the file is also printed to the screen on default verbosity including information on remaining L2 error and so on.