# Open SourceMolecular Dynamics

### Sidebar

#### For Developers

exercises:2018_ethz_mmm:h2o_md

# Molecular dynamics of water

In this exercise we will focus on the calculation of the self-diffusion coefficient for water

This exercise is mostly taken by a previous lab session by Marcella Iannuzzi, UZH, who should be credited and acknowledged here.

The simulation is run in classical molecular dynamics, using the force field used in 10.1063/1.1884609 to ab initio calculations of the $\text{H}_2\text{O}$ molecule.

The input file looks like that:

md300.inp
#
# MD of liquid water,
# using FF parameters from Praprotnik et al.
#

@SET SYSTEM T300

&GLOBAL
PROJECT ${SYSTEM} RUN_TYPE MD IOLEVEL LOW &PRINT PHYSCON &END PRINT &END GLOBAL &FORCE_EVAL METHOD FIST &MM &FORCEFIELD &BEND ATOMS H O H KIND HARMONIC K [rad^-2kcalmol] 110.0 THETA0 [deg] 104.52 &END BEND &BOND ATOMS O H KIND HARMONIC K [angstrom^-2kcalmol] 900.0 R0 [angstrom] 0.9572 &END BOND &CHARGE ATOM O CHARGE -0.834 &END CHARGE &CHARGE ATOM H CHARGE 0.417 &END CHARGE &NONBONDED &LENNARD-JONES atoms O O EPSILON [kcalmol] 0.152073 SIGMA [angstrom] 3.1507 RCUT [angstrom] 11.4 &END LENNARD-JONES &LENNARD-JONES atoms O H EPSILON [kcalmol] 0.0836 SIGMA [angstrom] 1.775 RCUT [angstrom] 11.4 &END LENNARD-JONES &LENNARD-JONES atoms H H EPSILON [kcalmol] 0.04598 SIGMA [angstrom] 0.400 RCUT [angstrom] 11.4 &END LENNARD-JONES &END NONBONDED &END FORCEFIELD &POISSON &EWALD EWALD_TYPE spme ALPHA .3 GMAX 12 O_SPLINE 6 &END EWALD &END POISSON &END MM &SUBSYS &CELL ABC 9.865 9.865 9.865 &END CELL &COORD O 4.0990771611739696E+00 2.8633390287706875E+00 -1.6649750693973509E+01 H2O H 3.6198825143472750E+00 2.9036335565407616E+00 -1.5786648392717291E+01 H2O H 3.3938346757245652E+00 2.6337360930791434E+00 -1.7258942999974757E+01 H2O O 2.0552971764987920E+01 -6.0317716888860673E-01 1.5803161825113377E+00 H2O H 2.1120490846876077E+01 1.8081944903699243E-01 1.5972671732140640E+00 H2O H 2.0893915175102396E+01 -1.1019809737906978E+00 2.2954442847330800E+00 H2O O -7.9877938931079253E+00 -4.1105886694798226E+00 5.1753650941730527E+00 H2O H -8.8790878338731307E+00 -4.3113400305239749E+00 4.7757679757154712E+00 H2O H -8.2891098658871574E+00 -3.8515146034780807E+00 6.0560547138619309E+00 H2O O -1.5699295004190399E+00 -1.1419760510561048E+01 3.6129811814629091E+00 H2O H -1.2515766648456859E+00 -1.1923591199013609E+01 4.3578597300341348E+00 H2O H -1.6014794034952906E+00 -1.0516592744819420E+01 3.8884643815900621E+00 H2O O 1.4346361535163705E+01 5.7379597570369754E+00 2.6136480784137825E+00 H2O H 1.3652893679436762E+01 6.3630082232359664E+00 2.2771706319341090E+00 H2O H 1.4142392462670498E+01 4.8901695938881256E+00 2.2500603744101371E+00 H2O O -1.2978502407406938E+00 -4.4197376154960653E+00 -1.0108922244726434E+01 H2O H -1.8209963244428125E+00 -4.0747449322863023E+00 -9.3928887615783125E+00 H2O H -1.8451970174557182E+00 -4.0605585415637524E+00 -1.0828278436385499E+01 H2O O 2.3505758212524612E+00 6.5988051487995394E+00 -9.0122820030981963E+00 H2O H 1.4916378847116123E+00 6.4099106638859480E+00 -8.5658667433350448E+00 H2O H 2.3259896427433664E+00 7.5392171863950361E+00 -9.0180780131236329E+00 H2O O -1.4312545906816201E+00 1.0119870177802614E+01 4.2355709830113553E-01 H2O H -1.9812309973608044E+00 9.3726251166317347E+00 7.5378044133890876E-01 H2O H -5.2677028459445807E-01 9.7434365832171981E+00 6.8082337693327510E-01 H2O O -7.2000125111789446E+00 -1.7297909920535137E+00 3.5519625527895298E+00 H2O H -7.5167017909363727E+00 -2.5625554479762518E+00 3.9663954244597615E+00 H2O H -7.4506425314478069E+00 -1.0888225237769447E+00 4.2397624267581158E+00 H2O O 1.0026382924586821E+01 6.4767710974957611E+00 2.3941661233133869E+00 H2O H 9.5657822147267115E+00 5.6711814281190902E+00 2.7635232090888051E+00 H2O H 9.6384789304574774E+00 7.1925052122997686E+00 2.8006359052237411E+00 H2O O 1.1844324850463268E+01 -2.0927058882474009E-01 -1.4065649668312794E+01 H2O H 1.2314749496538584E+01 4.4265839944891588E-01 -1.3587211076047117E+01 H2O H 1.1036049490850663E+01 -1.3482532936801822E-01 -1.3581077198077274E+01 H2O O -1.3112158228748623E+01 -6.2864535685630267E+00 2.0201791079602431E+00 H2O H -1.3342241444014668E+01 -6.7879338549113566E+00 1.2161643714725265E+00 H2O H -1.3917490102408042E+01 -6.4745955342905610E+00 2.5361963291835736E+00 H2O O 6.5504055434638353E+00 2.1722091632679059E+00 -5.2540920378685352E-01 H2O H 7.1181690032685019E+00 2.5160448913501052E+00 -1.2237604523620926E+00 H2O H 7.2032005514839970E+00 1.6132595514874122E+00 -8.4915330775898967E-02 H2O O 2.8728641994261084E+00 1.2463894784822655E+01 1.5766802157132808E+01 H2O H 3.5319661130868818E+00 1.2345097305742453E+01 1.6480886597851171E+01 H2O H 2.7260654590898588E+00 1.3425501065702791E+01 1.5689058285651134E+01 H2O O -4.7232088660258364E+00 6.1874970257059410E+00 -4.2010201320925882E+00 H2O H -5.7015580792394100E+00 6.3229099912716578E+00 -4.3313313722768019E+00 H2O H -4.5258103994562235E+00 6.2059968095552209E+00 -5.1454834292675624E+00 H2O O -1.1207108202024731E+01 1.2262911915342658E+01 -2.3602598361478542E+00 H2O H -1.1272681819601486E+01 1.2282332810869732E+01 -3.3705272042851795E+00 H2O H -1.0522707298067429E+01 1.1539271180629608E+01 -2.3426251777429385E+00 H2O O 1.0157679618370272E+01 -1.6491057241445276E+01 -3.8891647715560207E-01 H2O H 9.7050356527735868E+00 -1.6842658877059659E+01 -1.1648717278703189E+00 H2O H 9.5777475747466223E+00 -1.5735770228942069E+01 -1.4504538901147770E-01 H2O O -3.3867315013774526E+00 7.8050895970982230E+00 -8.3184507808782868E+00 H2O H -2.5883287464162126E+00 7.8752895157315681E+00 -7.7780603012568621E+00 H2O H -3.9367339883859653E+00 7.2104800149382875E+00 -7.8563335584074725E+00 H2O O 4.1388042326267591E+00 -1.5072806578102765E-01 1.4303824695133505E+00 H2O H 5.0472500573404178E+00 -4.8722470110749433E-01 1.4589770925283638E+00 H2O H 3.7467405608582829E+00 -6.3998028017083075E-01 2.1394965295434982E+00 H2O O 4.2595384151575244E+00 8.0363289325704235E-01 7.8672617022531943E+00 H2O H 4.8130766850341962E+00 1.2024067826714520E+00 7.1413687326139179E+00 H2O H 4.7953746600283598E+00 9.3924952891872560E-01 8.6335903258658764E+00 H2O O -4.5191738819706266E+00 4.8674040686154907E+00 -7.7424914554764745E-01 H2O H -3.9584835582169862E+00 4.1916762292422938E+00 -3.5632919560932813E-01 H2O H -5.3361240627326687E+00 4.6593170146920473E+00 -2.7748760122726546E-01 H2O O -4.6486190763083393E+00 1.0559684456132420E+01 2.4666697824632244E+01 H2O H -5.1718646071776799E+00 9.9061307184318412E+00 2.5191864943109294E+01 H2O H -5.2151318309590131E+00 1.0863848671746059E+01 2.4008154438222835E+01 H2O O -1.8920336095563112E-01 -1.9723082418168051E+01 -1.2331751859031851E+01 H2O H -6.9853839676937413E-01 -2.0510391457743804E+01 -1.2648106790669067E+01 H2O H -6.7178927685170767E-02 -1.9904341227137060E+01 -1.1401598461365145E+01 H2O O 8.6092522191222010E+00 -1.5300879172778661E+01 -6.2889540677378104E+00 H2O H 8.6164121625349033E+00 -1.5988040732713319E+01 -5.6163515969455950E+00 H2O H 7.8631907409852495E+00 -1.5565758302286353E+01 -6.8912232725417359E+00 H2O O 1.4533087231921096E+00 -1.3305185412448356E+01 -1.1916886354803323E+01 H2O H 1.5704975438853019E+00 -1.2389522343553942E+01 -1.1638778267823120E+01 H2O H 2.1592317355561943E+00 -1.3688032462538523E+01 -1.1445000291043685E+01 H2O O 2.9687705641225186E+00 -5.7847368452718690E+00 1.0198124958568776E+01 H2O H 2.1055313600020069E+00 -6.1369494659703072E+00 9.8643510074097875E+00 H2O H 2.7673476315489345E+00 -4.8175890711447638E+00 1.0261860604897310E+01 H2O O 2.1307179526981219E+00 9.0141614782957706E+00 8.8029803669936140E+00 H2O H 2.6906994240494826E+00 9.4300032870082422E+00 9.5306845691514326E+00 H2O H 2.5357829972259256E+00 9.3708237774629808E+00 7.9904382164776395E+00 H2O O -2.7590496648261693E+00 -3.1955215144759475E+00 -2.3261918496018064E+00 H2O H -3.3546076991918565E+00 -3.3511138435382928E+00 -3.1173095619864135E+00 H2O H -3.4116085527171047E+00 -3.3919867560013448E+00 -1.5671279742947930E+00 H2O O -7.9159449798759551E+00 1.8733590643400650E+00 -8.2162396984353716E+00 H2O H -8.7859157283547180E+00 2.2519209017166939E+00 -8.5544872930068880E+00 H2O H -7.2213059662322143E+00 2.0590146256161748E+00 -8.9557258250561773E+00 H2O O 7.4839506528220578E+00 1.2172260395768921E+01 5.0433745483595249E+00 H2O H 6.8305009483129693E+00 1.1483218223760431E+01 4.6766276844537620E+00 H2O H 6.9209089699397879E+00 1.2599276062084877E+01 5.7102955698020521E+00 H2O O -2.0365405094837961E+01 -1.2256879300334079E+01 6.2150050341365217E+00 H2O H -1.9797664340014357E+01 -1.2734530157707740E+01 6.7789740677763728E+00 H2O H -2.1192917559895324E+01 -1.2576379478142982E+01 6.7398537673836154E+00 H2O O 5.3894209467850960E+00 3.4139122557112409E+00 -3.3441897135960801E+00 H2O H 5.4941839731265079E+00 4.0025859745240160E+00 -2.5460134060473312E+00 H2O H 5.3456451321990039E+00 4.1287379135378721E+00 -3.9781343543418788E+00 H2O &END COORD &END SUBSYS &END FORCE_EVAL &MOTION &MD ENSEMBLE NVE TIMESTEP 0.5 STEPS 50000 TEMPERATURE 300 &PRINT &ENERGY &EACH MD 20 &END &END ENERGY &PROGRAM_RUN_INFO &EACH MD 20 &END &END PROGRAM_RUN_INFO &END PRINT &END MD &PRINT &TRAJECTORY SILENT FILENAME =${SYSTEM}.xyz
&EACH
MD 20
&END EACH
&END TRAJECTORY
&VELOCITIES ON
&EACH
MD 20
&END EACH
&END VELOCITIES
&FORCES OFF
&END FORCES
&RESTART_HISTORY OFF
&END RESTART_HISTORY
&RESTART OFF
&END RESTART
&END PRINT
&END MOTION
1. Read the paper through section III to see which parameters are used for which interaction.
2. Recognize the parameters in the input file
3. Run the md300.in: the command is
cp2k.ssmp -i md300.in > md300.out
4. Postprocess the trajectory to compute mean square displacement and velocity-velocity correlation function
5. Extract the diffusion coefficients and compare the two approaches (the instruction about the parameters: ./command.py -h)
./vel_autocorr.py T300-vel-1.xyz 3 10 100
./vel_autocorr_plot.py vel_autocorr.out
./mean_square_disp.py T300.xyz 3 10 100
./msd_fit_and_plot.py mean_square_disp.out

Remember: to display the png files, the command is

display file.png
Assignments
1. look at the file T300-1.ener at the end of the simulation. It contains several quantities. Check the header!
2. Using gnuplot, check the stability of the MD, you can also use
 ./simpleplot.py Col1 Col2
3. How do you interpret the oscillation of the potential energy? And of temperature?
4. Compare the values of the diffusion coefficient with the two methods. Are they compatible?
5. Look at the trajectories with vmd. Why are the particles “exiting” the box with time? What is the importance of this for the algorithm?
6. Copy md300.in into md200.in and change the system name and the initial temperature. Run the md. What is the final temperature? Why?
7. Copy md300.in into md400.in and change to 400 in the two places. Run the md. What is the final temperature? Why?
8. Check the temperature dependence of the diffusion coefficient. Plot the result.