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.
For our computer experiments, we are going to use the CP2K software. CP2K is under active development and some exercises require a recent CP2K version.
cp2k.sopt -h # get version and revision number of your cp2k executable
CP2K does not have a graphical user interface, instead we use input files to tell CP2K what to do.
Our input files have the file extension .in and below you see a commented example of such a file.
Many parameters are quite technical and not important at the moment, still it is very helpful to get an overview of the basic structure.
Can you guess what it does?
- gopt.in
@SET SYSTEM run-1H2O-GOPT @SET LIBDIR /Network/Servers/130.60.15.139/Users/c329s00/cp2k-files &GLOBAL 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 &END GLOBAL &FORCE_EVAL METHOD QuickStep &DFT # Use density functional theory BASIS_SET_FILE_NAME ${LIBDIR}/GTH_BASIS_SETS POTENTIAL_FILE_NAME ${LIBDIR}/GTH_POTENTIALS &MGRID CUTOFF 250 # plane-wave cutoff for the charge density [Rydbergs] &END MGRID &SCF EPS_SCF 1.0E-6 # convergence threshold for total energy &END SCF &XC # parameters for exchange-correlation functional &XC_FUNCTIONAL BLYP &END XC_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 &END DFT &SUBSYS &CELL # box containing the molecule: 15x15x15 Angstroms ABC [angstrom] 15 15 15 &END CELL &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 &END COORD &KIND H # basis sets and pseudo-potentials for atomic species BASIS_SET TZVP-GTH POTENTIAL GTH-BLYP-q1 &END KIND &KIND O BASIS_SET TZVP-GTH POTENTIAL GTH-BLYP-q6 &END KIND &END SUBSYS &END FORCE_EVAL
The input file gopt.in 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 find the optimal geometry of the water molecule.
cp2k.sopt -i gopt.in # run cp2k, writing output to screen cp2k.sopt -i gopt.in -o gopt.out # run cp2k, writing output to file gopt.out cp2k.sopt -i gopt.in -o gopt.out & # as before, but run in background
CP2K creates several new files in the folder, most of which contain information on previous steps of the calculation (allowing e.g. to restart the calculation after a crash) and are not relevant at the moment.
We are interested in the file run-1H2O-GOPT-pos-1.xyz, which contains the atomic positions during the steps of the geometry optimization. Visualize the optimized geometry with VMD.
- 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.
- 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 resp.in with the optimized geometry and and start the fit.
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.
- Besides fitting the electrostatic potential, one often defines additional constraints (to be fulfilled exactly) or restraints (to be fulfilled approximately). Analyze
resp.into figure out which constraints and/or restraints are applied in our fit. - Compare the fitted charges with the ones used by Praprotnik et al.
- 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 vibr.in with the optimized geometry and start the calculation of the normal modes.
- How many normal modes does the water molecule have and why?
- 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)
- 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 md.in input file, replace the initial coordinates for the AIMD by the ones determined from the geometry optimization
and start the simulation.
- What is the physical time duration you have simulated? How much computing time was required per MD step?
- What type of ensemble are you simulating? What are the initial and final temperatures of the simulation and why are they different? (2P)
- Why is the final average temperature approximately 1/2 of the initial temperature? Hint: Use the virial theorem. (2P)
- 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 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. 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.
Note: Since we are considering an isolated water molecule, the Lennard-Jones interactions are explicitly set to zero in ff_match.in.
- 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.
- Compare the output (geometry, vibrational frequencies, average bond length in MD) with the ab initio results. (3P)
- How much faster is the classical MD compared to the ab initio one?