Differences

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

--- exercises:common:ensemble [2022/10/31 16:07] – created jglan
+++ exercises:common:ensemble [2022/10/31 17:28] (current) – [NVT Ensemble] jglan
@@ Line 1: / Line 1: @@
-====== Lennard-Jones liquids ======
+====== Ensembles (Lennard-Jones liquids) ======
 In this exercise, you will simulate a fluid of monoatomic particles that interact with a Lennard-Jones potential. The method to be used is molecular dynamics (MD) with periodic boundary conditions using CP2K. The aim is to explore the method, calculate the radial distribution function  $g(r)$  and investigate a variety of ensembles.
-You are expected to hand in the short report via OLAT, ONLY in PDF format.
- ====== Run your first simulation using CP2K ======
-When you check CP2K's [[::features]] and the outline of the lecture you will notice that there are many levels of theory, methods, and possibilities to combine them. This results in a large number of possible options and coefficients to setup, control and tune a specific simulation.
-Together with the parameters for the numerical solvers, this means that an average CP2K configuration file will contain quite a number of options (even though for many others the default value will be applied) and not all of them will be discussed in the lecture or the exercises.
-The [[https://manual.cp2k.org/|CP2K Manual]] is the complete reference for all configuration options. Where appropriate you will find a reference to the respective paper when looking up a specific keyword/option.
-To get you started, we will do a simple exercise using Molecular Mechanics (that is: a classical approach). The point is to get familiar with the options, learn how to organize and edit the input file and analyze the output.
 ====== Background ======
@@ Line 24: / Line 16: @@
 The radial distribution function, (or pair correlation function)  $g(r)$ , in a system of particles (atoms, molecules, colloids, etc.), describes how the density varies as a function of distance from a reference particle.
-===== Part I:  Set up MD simulation  =====
+===== NVE Ensemble  =====
 In this section, we provide you with an example CP2K input for an MD calculation.
-Extensive comments have been added to the file, which start with a has symbol '#'.
+Extensive comments have been added to the file, which starts with the symbol '#'.
-=== 1. Step ===
 Load the CP2K module as explained in Exercise 0, create a directory ''ex1'' and change to it:
@@ Line 48: / Line 39: @@
 &GLOBAL
   PROJECT ar108           #Project Name
-  RUN_TYPE md             #Calculation Type : MD (molecular dynamics), GEO_OPT (Geometry Optimization), Energy (Energy Calculation)
+  RUN_TYPE md             #Calculation Type: MD (molecular dynamics), GEO_OPT (Geometry Optimization), Energy (Energy Calculation)
 &END GLOBAL
@@ Line 66: / Line 57: @@
       &CHARGE             #charge of the MM atoms
         ATOM Ar           #Defines the atomic kind of the charge
-        CHARGE 0.0        #Defines the charge of the MM atom in electron charge unit
+        CHARGE 0.0        #Defines the charge of the MM atom in an electron charge unit
       &END
       &NONBONDED
@@ Line 222: / Line 213: @@
   *Run calculations with different timesteps (0.5, 2, 5fs), different temperatures(84, 300, 400K), different densities( $\rho$ =0.25, 0.5, 1), and a different number of steps (100, 1000, 5000), and inspect the geometry in each case. Plot the total energy, temperature, and potential energy. Try to comment and explain what you observe.
 </note>
-===== Part II:  Force Field Parameter  =====
-In this section we investigate the dependence of the L-J potential on the its different parameters. To do se, we investigate a small system of only two Ar atoms. The following code snippet shows the changes that you need to make to the input file in order to run such a calculation.
-      &GLOBAL                  ! section to select the kind of calculation
-         RUN_TYPE ENERGY       ! select type of calculation. In this case: ENERGY (=Single point calculation)
-      &END GLOBAL
-      ...
-      &NONBONDED
-        &LENNARD-JONES    #LENNARD-JONES potential type.Functional form: V(r) = 4.0 * EPSILON * [(SIGMA/r)^12-(SIGMA/r)^6]
-          atoms Ar Ar     #Defines the atomic kind involved in the nonbonded potential
-          EPSILON 119.8   #Defines the EPSILON parameter of the LJ potential (K_e)
-          SIGMA 3.405     #Defines the SIGMA parameter of the LJ potential (Angstrom)
-          RCUT 8.4        #Defines the cutoff parameter of the LJ potential
-        &END LENNARD-JONES
-      &END NONBONDED
-      ...
-      &SUBSYS                 ! system description
-       &CELL
-         ABC [angstrom] 10 10 10
-         PERIODIC NONE
-       &END CELL
-      &COORD
-        UNIT angstrom
-        Ar  0 0 0
-        Ar  4 0 0
-      &END COORD
-  &END SUBSYS
-In order to investigate the effect of the force-field parameters on the L-J potential you need to vary multiple parameters of your calculation:
+===== NVT Ensemble  =====
-  - you need to actually vary the L-J force-field parameters,  $\epsilon$  and  $\sigma$
-  - you need to run each calculation at different distances between the two Ar atoms
-A simple way to generate the different input files is using shell scripting in combination with ''sed'' (the stream editor):
+In the previous sections, you have already run NVE ensemble molecular dynamics calculations for liquid Ar. In this section, we will focus on the NVT ensembles.
-<code>
+Although the most popular Nose-Hoover thermostat is commonly-used in other MD codes, the original Nose-Hoover thermostat has an ergodic issue. This has been solved by Mark Tuckerman et al. See [[https://doi.org/10.1063/1.463940 |Nose-Hoover Chain Thermostat]]
-for d in $(seq 2 0.1 4); do
+In CP2K, the default length of the Nose-Hoover chain is set to 3. (See Manual [[https://manual.cp2k.org/cp2k-2022_1-branch/CP2K_INPUT/MOTION/MD/THERMOSTAT/NOSE.html| %MOTION%MD%THERMOSTAT%NOSE]])
-  sed -e "s|4 0 0|${d} 0 0|" energy.inp > energy_${d}A.inp
-  cp2k.popt -i energy_${d}A.inp -o energy_${d}A.out
-  awk '/Total FORCE_EVAL/ { print  $9; }' energy_$ {d}A.out
+To set up an NVT calculation, change the settings in the &MD section as shown below:
-done
-</code>
+  &MD
+    ENSEMBLE NVT
+    STEPS 3000
+    TIMESTEP 5
+    TEMPERATURE 298
+    &THERMOSTAT
+      REGION MASSIVE
+      &NOSE                    #Uses the Nose-Hoover thermostat
+        TIMECON 100           #timeconstant of the thermostat chain, how often does the thermostat adjust your system
+      &END NOSE
+    &END
+  &END MD
+Alternatively, one can also use [[https://manual.cp2k.org/cp2k-2022_1-branch/CP2K_INPUT/MOTION/MD/THERMOSTAT/CSVR.html| Canonical sampling through velocity rescaling
+(CSVR)]] as developed by Giovanni [[https://doi.org/10.1063/1.2408420 |Bussi et al. ]]
+  &MD
+    ENSEMBLE NVT
+    STEPS 3000
+    TIMESTEP 5
+    TEMPERATURE 298
+    &THERMOSTAT
+      &CSVR
+        TIMECON 100           #timeconstant of the CSVR, how often does the thermostat adjust your system
+      &END CSVR
+    &END
+  &END MD
+===== NPT Ensemble  =====
+To set up an NPT calculation, change the settings in the &MD section as shown below:
+  &FORCE_EVAL
+  ...
+  STRESS_TENSOR ANALYTICAL
+  ...
+  &END FORCE_EVAL
+  ...
+  ...
+  &MD
+    ENSEMBLE NPT_I                #constant temperature and pressure using an isotropic cell
+    STEPS 3000
+    TIMESTEP 5.
+    TEMPERATURE 85.0
+    &BAROSTAT
+      PRESSURE 0.                 # PRESSURE, unit[bar]
+      TIMECON 1000
+    &END BAROSTAT
+    &THERMOSTAT
+      &NOSE
+        TIMECON 1000
+      &END NOSE
+    &END THERMOSTAT
+   &END MD
-  * The command ''seq 2 0.1 4'' generates the numbers ''2.0'', ''2.1'', ''2.2'', ... , ''4.0'' (try it out!)
-  * With ''for d in  $(seq 2 0.1 4); do'' we use the shell to run all commands which follow once for every number (stored in ''$ d'')
-  * ''sed -e "s|4 0 0| $d 0 0|" energy.inp'' looks for ''4 0 0'' in the file ''energy.inp'' (the original file from above) and replaces ''4 0 0'' by ''$ d 0 0'' (that is: ''2.0'', ''2.1'', ''2.2'', ...)
-  * ... and using ''> energy_${d}A.out'' we redirect the output of the ''sed'' command to new files ''energy_2.0A.out'', ''energy_2.1A.out'', etc.
-  * Then we run ''cp2k.popt'' as shown in Exercise 0 on those new input files and write the output to new output files as well
-  * Using ''awk'' we extract the energy from the output file and print it
 <note>
 **TASK**
-  *Plot the Lennard-Jones potential against the Ar-Ar distance  $r$  (2-4 Å) using different values for  $\epsilon$  and  $\sigma$ .
+   *Run a calculation using the NVT ensemble at 300K. Check the temperature and energy of the whole system, and compare the result to an NVE ensemble (300K). Rationalize and discuss the difference.
-  *Repeat the L-J MD calculation with different  $\epsilon$  and  $\sigma$ , compare the potential energy and temperature evolution, and explain the relation between the calculated properties and force field parameters.
+   *Run a calculation using the NVT ensemble (300K) until the system is equilibrated, then run an NVE ensemble calculation. Check the temperature and energy of the whole system, and compare to the previous NVE simulation. Discuss your observations.
+   *Remove some atoms and run an NPT ensemble simulation. Then, check the size of the simulation box. Discuss your observations.
+</note>
+<note tip>
+You have multiple options on how to restart a CP2K calculation off of a previous one. What all approaches have in common, is that you need to make use of the RESTART-files which are automatically written by CP2K (unless you explicitly disable them).
+For the purposes of this example, you should see a file called ''ar108-1.restart'' (at least for Part 1 of this exercise). //Note:// there will also be multiple backup files (''.bak-#'' suffixes) which you do not need to care about.
+These files are nothing but another input file. However, their parameters are set such that they continue a CP2K calculation from the last step of the simulation which generated the RESTART file.
+Here are two options for how you can use these RESTART-files:
+. Directly using the RESTART as an input.
+   - You can copy the RESTART file to a new input file: <code>cp ar108-1.restart argon_follow_up.inp</code>
+   - Now you can change the input to your liking (e.g. change the ensemble, etc.)
+   - And finally simply run CP2K with the new input file: <code>cp2k -i argon_followup.inp -o argon_followup.out</code>
+. You can also tell CP2K to load a specific RESTART-file.
+   - Write a new input file as usual: <code>argon_followup.inp</code>
+   - Add an [[https://manual.cp2k.org/trunk/CP2K_INPUT/EXT_RESTART.html|EXT_RESTART]] section:
+      <code>
+&EXT_RESTART
+  RESTART_FILE_NAME ar108-1.restart
+&END EXT_RESTART</code>
+   - And now, again, simply run CP2K:  <code>cp2k -i argon_followup.inp -o argon_followup.out</code>
 </note>
-===== Part III:  Radial distribution function  =====
+===== Radial distribution function  =====
 In this section we analyze the dependence of the radial distribution function (rdf),  $g(r)$ , on the temperature of the system. To do so, you should plot  $g(r)$  against various temperatures and examine the effects.
 You can use VMD (as explained below) or write your own program (Fortran, C, C++, Python etc.) to calculate the rdf.
@@ Line 455: / Line 491: @@
 </code>
-===== Part IV: Other Ensembles  =====
-In the previous sections, you have already run NVE ensemble molecular dynamics calculations for liquid Ar. In this section, we will focus on the NVT and NPT ensembles.
-To set up an NVT calculation, change the settings in the &MD section as shown below:
-  &MD
-    ENSEMBLE NVT
-    STEPS 3000
-    TIMESTEP 5
-    TEMPERATURE 298
-    &THERMOSTAT
-      REGION MASSIVE
-      &NOSE                    #Uses the Nose-Hoover thermostat
-        TIMECON 1000           #timeconstant of the thermostat chain, how often does thermostat adjust your system
-      &END NOSE
-    &END
-  &END MD
-To set up an NPT calculation, change the settings in the &MD section as shown below:
-  &FORCE_EVAL
-  ...
-  STRESS_TENSOR ANALYTICAL
-  ...
-  &END FORCE_EVAL
-  ...
-  ...
-  &MD
-    ENSEMBLE NPT_I                #constant temperature and pressure using an isotropic cell
-    STEPS 3000
-    TIMESTEP 5.
-    TEMPERATURE 85.0
-    &BAROSTAT
-      PRESSURE 0.                 # PRESSURE, unit[bar]
-      TIMECON 1000
-    &END BAROSTAT
-    &THERMOSTAT
-      &NOSE
-        TIMECON 1000
-      &END NOSE
-    &END THERMOSTAT
-   &END MD
-<note>
-**TASK**
-   *Run a calculation using the NVT ensemble at 300K. Check the temperature and energy of the whole system, and compare the result to an NVE ensemble (300K). Rationalize and discuss the difference.
-   *Run a calculation using the NVT ensemble (300K) until the system is equilibrated, then run an NVE ensemble calculation. Check the temperature and energy of the whole system, and compare to the previous NVE simulation. Discuss your observations.
-   *Remove some atoms and run an NPT ensemble simulation. Then, check the size of the simulation box. Discuss your observations.
-</note>
-<note tip>
-You have multiple options on how to restart a CP2K calculation off of a previous one. What all approaches have in common, is that you need to make use of the RESTART-files which are automatically written by CP2K (unless you explicitly disable them).
-For the purposes of this example, you should see a file called ''ar108-1.restart'' (at least for Part 1 of this exercise). //Note:// there will also be multiple backup files (''.bak-#'' suffixes) which you do not need to care about.
-These files are nothing but another input file. However, their parameters are set such that they continue a CP2K calculation from the last step of the simulation which generated the RESTART file.
-Here are two options for how you can use these RESTART-files:
-. Directly using the RESTART as an input.
-   - You can copy the RESTART file to a new input file: <code>cp ar108-1.restart argon_follow_up.inp</code>
-   - Now you can change the input to your liking (e.g. change the ensemble, etc.)
-   - And finally simply run CP2K with the new input file: <code>cp2k -i argon_followup.inp -o argon_followup.out</code>
-. You can also tell CP2K to load a specific RESTART-file.
-   - Write a new input file as usual: <code>argon_followup.inp</code>
-   - Add an [[https://manual.cp2k.org/trunk/CP2K_INPUT/EXT_RESTART.html|EXT_RESTART]] section:
-      <code>
-&EXT_RESTART
-  RESTART_FILE_NAME ar108-1.restart
-&END EXT_RESTART</code>
-   - And now, again, simply run CP2K:  <code>cp2k -i argon_followup.inp -o argon_followup.out</code>
-</note>