Exercise Week 1: MD

The first week’s exercise will be to construct a force field for the H₂O molecule by fitting the parameters to ab initio calculations. We will first carry out a geometry optimization with CP2K, then we will perform a charge and normal mode analysis, and then we will run an ab initio MD simulation for a single water molecule. After this, we will do force matching to obtain refined parameters for our water model. The second week’s exercise will be to run a classical MD simulation for liquid water system using the force field we parametrized, and then computing the radial distribution function (RDF) for it.

The required iput files can be obtained from here.

An example input file for the geometry optimization is given below. This input tells CP2K to perform a geometry optimization for H₂O molecule using density functional theory with the BLYP exchange-correlation functional.

First go to the directory 1.1-1H2O-GOPT, the input file (gopt.in) for the geometry optimization is given there.

The job can be submitted by the command:

After the calculation is completed, you can check the output file to see if the SCF steps are converged or not.

The file we are mostly interested in is run-1H2O-GOPT-pos-1.xyz. This file contains the atomic positions for each geometry optimization step. You can visualize how the geometry changes using VMD.

Now we can also compute atomic point charges. There are many different models for associating point charges to atoms. The one we will be using is RESP, which fits the electrostatic potential directly in real space. The starting point is the electrostatic potential, which is provided by a quantum mechanical calculation. One then tries to place charges on atoms in such a way that the potential generated by them matches the original potential as closely as possible. Since the charges are intended to be used for interactions with other molecules, the potential is fitted only in the spatial region outside the molecule. The excluded volume is typically defined by the union of the van der Waals spheres around each atom of the molecule.

First go to the directory 1.2-1H2O-RESP and edit the input file resp.in. Change the dummy coordinates with the optimized ones which you can directly take from the last optimization step of run-1H2O-GOPT-pos-1.xyz. Then run the calculation again.

Once the calulation is finished, check the atomic charges from the output file resp.out.

There are several options to determine the force constants for our water model. We calculate it from the frequencies of the vibrational normal modes of our molecule.

Go to the directory 1.3-1H2O-VIBR and replace the dummy coordinates in vibr.in with the optimized geometry and start the calculation of the normal modes.

In the case of more complex model potentials (e.g. such as the non-analytic embedded atom model) an 'analytic fit' can become impossible. Alternatively, one may follow the more generally applicable force matching approach: missing parameters of the model potential are chosen such as to best reproduce the forces computed along a MD trajectory (or any other given set of atomic configurations).

Therefore, we will now perform ab initio molecular dynamics (AIMD) of the single water molecule in order to sample its potential energy landscape and compute the forces along the trajectory.

Go to the directory 1.4-1H2O-AIMD and in the md.in input file, replace the initial coordinates for the AIMD by the ones determined from the geometry optimization and start the simulation.

Once the AIMD simulation is done, we can match the force constants. In order to be applicable to a wide range of potentials, the implementation of the parameter search in CP2K does not rely on analytic derivatives of the potentials with respect to the parameters, but uses the trial- and error-like Powell search algorithm.

Go to the directory 1.5-1H2O-MATCH. In the input file ff_match.in, replace the dummy values for the preferred O-H bond length, H-O-H angle and the point charges by the parameters determined from the geometry optimization and the RESP fit. Also, replace the dummy variables for the force constants by the parameters given by Praprotnik et al.. Alternatively, if you did the Task 1.3, you can also use the force constants you calculated, or the ones calculated by CP2K. However, pay attention to the units!

Copy the trajectory run-1H2O-AIMD.xyz and the forces run-1H2O-AIMD.force to the new directory and replace the dummy filenames in ff_match.in. Run ff_match.in to fit the stretching and bending force constants via force matching.

Finally, you are ready to use your own force field! The next week’s exercise will be to perform a classical MD simulation for liquid water using the fitted parameters.

For more MD exercises with CP2K, please look at here.