User Tools

Site Tools


exercises:2015_uzh_molsim:h2o_ff

Differences

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

Link to this comparison view

exercises:2015_uzh_molsim:h2o_ff [2015/05/04 11:20]
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. 
- 
-<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. 
- 
- 
-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 find the optimal geometry of the water molecule. 
-<code bash> 
-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 
-</code> 
- 
-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. 
-<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> 
exercises/2015_uzh_molsim/h2o_ff.txt · Last modified: 2020/08/21 10:15 (external edit)