User Tools

Site Tools


Constructing a force field for the $\text{H}_2\text{O}$ molecule

In this exercise, we are going to construct a force field for the $\text{H}_2\text{O}$ molecule by fitting the parameters of the water model used in 10.1063/1.1884609 to ab initio calculations of the $\text{H}_2\text{O}$ molecule.


Download the paper and read (at least) through section III, part A. Describe the water model in your own words. Which potentials are used for which interaction? What are the free parameters? (3P)

We start with ab initio calculations that will form the basis for our fitting procedure.

Have a look at cp2k input file. Can you guess what it does?
@SET LIBDIR /Network/Servers/
  PROJECT ${SYSTEM}          # Name of calculation: run-1H2O-GOPT 
  RUN_TYPE GEO_OPT           # Perform a geometry optimization
  PRINT_LEVEL LOW            # Do not output too much information
  METHOD QuickStep
  &DFT                       # Use density functional theory
      CUTOFF 250             # plane-wave cutoff for the charge density [Rydbergs]
    &END MGRID  
      EPS_SCF 1.0E-6         # convergence threshold for total energy
    &END SCF
    &XC                      # parameters for exchange-correlation functional
      &XC_GRID               # some tricks to speed up the calculation of the xc potential
        XC_DERIV       SPLINE2_smooth
        XC_SMOOTH_RHO  NN10
      &END XC_GRID
    &END XC
    &CELL                    # box containing the molecule: 15x15x15 Angstroms
      ABC [angstrom] 15 15 15
    &COORD                   # coordinates of the atoms in the box [Angstroms]
      O  0 0 0
      H  0.7 0.7 0
      H -0.7 0.7 0
    &KIND H                  # basis sets and pseudo-potentials for atomic species
    &KIND O

The input file tells CP2K to perform a geometry optimization of a $\text{H}_2\text{O}$ molecule using density functional theory with the BLYP exchange-correlation functional.

Run CP2K to perform a geometry optimization of a water molecule.

cp2k.sopt -i -o gopt.out 

We are interested in the file, which contains the atomic positions during the steps of the geometry optimization. Visualize the optimized geometry with VMD.


  1. Which parameters of our water model can be obtained directly from the optimized geometry? Measure them with VMD and include a picture in the report. Hint: You can adjust the position of labels through Graphics → Labels → Properties.
  2. How do the values compare to the ones used in the paper?

Another set of parameters that can be obtained from our calculation are the atomic point charges. There are many different models for associating point charges to atoms. Some determine, which part of the electron density belongs to a given atom (Mulliken, Bader), some perform a multipole expansion of the charge density and determine the atomic charges that best reproduce it (Blöchl, Hirshfeld), while others fit the electrostatic potential directly in real space (ESP, CHELPG, RESP).

Here, we use a restrained electrostatic potential (RESP) fit: The starting point is the electrostatic potential, which is provided by a quantum mechanical calculation. One then tries to place charges on atoms (here: Gaussian-shaped charge distributions centered on the atomic positions) in such a way that the potential generated by them matches the original potential as best 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.

Replace the dummy coordinates in with the optimized geometry and and start the fit.

In order to speed up your calculation, you can first copy the file run-1H2O-GOPT-RESTART.wfn to your working directory and rename it to run-1H2O-RESP-RESTART.wfn. This file already contains the converged wave function of the system at the optimized geometry, thus saving CP2K the effort to recompute it.


  1. Besides fitting the electrostatic potential, one often defines additional constraints (to be fulfilled exactly) or restraints (to be fulfilled approximately). Analyze to figure out which constraints and/or restraints are applied in our fit.
  2. Compare the fitted charges with the ones used by Praprotnik et al.
  3. Using your geometry and fitted point charges, calculate the dipole moment of $\text{H}_2\text{O}$. How does it compare to the experimental value of 1.85 Debye (measured for isolated water molecules)?

In order to determine the force constants for our water model, we have several options. The first option is to calculate the frequencies of the vibrational normal modes of our molecule.

Replace the dummy coordinates in with the optimized geometry and start the calculation of the normal modes.

Before writing the results of the vibrational analysis, CP2K may take some time without writing any output. Don't worry, this is normal.


  1. How many normal modes does the water molecule have and why?
  2. Praprotnik et al. provide vibrational frequencies of the TIP3P model as well as experimental results from gas phase infrared spectroscopy. Compare them to your results and discuss the differences. Hint: Hint. (2P)
  3. BONUS In order to complete our water model, we need the force constants and not the frequencies. Calculate the force constants directly from the frequencies. See e.g. Landau-Lifshitz, chapter V – Vibrations of molecules. (3P)

In the case of more complex model potentials (e.g. such as the non-analytic EAM) 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 water molecule in order to sample its potential energy landscape and compute the forces along the trajectory.

In the input file, replace the initial coordinates for the AIMD by the ones determined from the geometry optimization and start the simulation.

This simulation will take some time. You may continue with the next exercise, while it is running.


  1. What is the physical time duration you have simulated? How much computing time was required per MD step?
  2. What type of ensemble are you simulating? What are the initial and final temperatures of the simulation and why are they different? (2P)
  3. Why is the final average temperature approximately 1/2 of the initial temperature? Hint: Use the virial theorem. (2P)
  4. Calculate a histogram of the O-H bond length and H-O-H bending angle (e.g. using gnuplot) and determine their averages. Hint: Use the gnuplot script plot.gnu, adjusting filename and bin-width as needed. (2P)

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.

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. Copy the trajectory and the forces run-1H2O-AIMD.force to the new directory and replace the dummy filenames in Run to fit the stretching and bending force constants via force matching.

Note: Since we are considering an isolated water molecule, the Lennard-Jones interactions are explicitly set to zero in


  1. Provide a short description of all fitted force field parameters in the report.

Finally, you are ready to use your own force field! Using the provided input files to perform a geometry optimization, vibrational analysis and MD with your force field, replacing the DUMMY values in the input files with the fitted values ones you have determined.


  1. Compare the output (geometry, vibrational frequencies, average bond length in MD) with the ab initio results. (3P)
  2. How much faster is the classical MD compared to the ab initio one?

exercises/2015_uzh_molsim/h2o_ff.txt · Last modified: 2015/05/04 13:25 by yakutovich