User Tools

Site Tools


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

All files of this exercise be downloaded directly from this link

Download the 2.1 exercise into your EXERCISES folder and unzip it.

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.
exercises/2018_ethz_mmm/h2o_md.txt · Last modified: 2020/08/21 10:15 by 127.0.0.1