====== 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 [[doi>10.1063/1.1884609]] to //ab initio// calculations of the $\text{H}_2\text{O}$ molecule.

<note>**TASK 1**

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) 
</note>

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?

<code bash 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
</code>

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 perform a geometry optimization of a water molecule.

<code bash>
cp2k.sopt -i gopt.in -o gopt.out 
</code>

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.
<note>**TASK 2** 

  - 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?
</note>

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 ([[doi>10.1063/1.1740588|Mulliken]], [[http://theory.cm.utexas.edu/henkelman/research/bader/|Bader]]), some perform a multipole expansion of the charge density and determine the atomic charges that best reproduce it ([[doi>10.1063/1.470314|Blöchl]], [[doi>10.1007/BF00549096|Hirshfeld]]), while others fit the electrostatic potential directly in real space (ESP, [[doi>10.1002/jcc.540110311|CHELPG]], [[doi>10.1021/j100142a004|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.
<note tip>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.</note>
<note>**TASK 3** 
  - Besides fitting the electrostatic potential, one often defines additional //constraints// (to be fulfilled exactly) or //restraints// (to be fulfilled approximately). Analyze ''resp.in'' to 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 [[http://en.wikipedia.org/wiki/Electric_dipole_moment|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)?
</note>

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.

<note tip>
Before writing the results of the vibrational analysis, CP2K may take some time without writing any output. Don't worry, this is normal.
</note>
<note>**TASK 4** 
  - 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:// [[doi>10.1021/jp046158d|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. [[http://books.google.ch/books?id=bE-9tUH2J2wC&hl=de|Landau-Lifshitz]], chapter V – //Vibrations of molecules//. (3P)
</note>

In the case of more complex model potentials (e.g. such as the non-analytic [[http://en.wikipedia.org/wiki/Embedded_atom_model|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.
<note tip>This simulation will take some time. You may continue with the next exercise, while it is running.</note>
<note>**TASK 5** 

  - 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 [[http://en.wikipedia.org/wiki/Virial_theorem|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)
</note>

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 [[doi>10.1093/comjnl/7.2.155|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''.
<note>**TASK 6** 

  - Provide a short description of all fitted force field parameters in the report.
</note>


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.

<note>**TASK 7** 

  - 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?
</note>