# CP2K Open Source Molecular Dynamics

### Site Tools

exercises:2015_uzh_molsim:h2o_ff

# Differences

This shows you the differences between two versions of the page.

 exercises:2015_uzh_molsim:h2o_ff [2015/05/04 11:25]yakutovich exercises:2015_uzh_molsim:h2o_ff [2020/08/21 10:15] Line 1: Line 1: - ====== 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. - - **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) - - - 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 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 perform a geometry optimization of a water molecule. - - - cp2k.sopt -i gopt.in -o gopt.out - - - 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. - **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? - - - 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. - 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. - **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)? - - - 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. - - - Before writing the results of the vibrational analysis, CP2K may take some time without writing any output. Don't worry, this is normal. - - **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) - - - 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. - This simulation will take some time. You may continue with the next exercise, while it is running. - **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) - - - 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''. - **TASK 6** - - - 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. - - **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? -