CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
exercises:2021_uzh_acpc2:ex01

Differences

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

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
exercises:2021_uzh_acpc2:ex01 [2021/04/26 06:38] – [Lennard-Jones liquids] mrossmannekexercises:2021_uzh_acpc2:ex01 [2021/05/17 11:35] (current) – [Part III: Radial distribution function] Fix type mrossmannek
Line 10: Line 10:
 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. 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, organizing and editing the input file and analyze the output.+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 ====== ====== Background ======
 You are expected to carry out an MD simulation of Lennard-Jones (L-J) fluid containing mono-atomic You are expected to carry out an MD simulation of Lennard-Jones (L-J) fluid containing mono-atomic
-particles. The [[https://en.wikipedia.org/wiki/Lennard-Jones_potential | L-J potential ]] energy expression used is+particles. The [[https://en.wikipedia.org/wiki/Lennard-Jones_potential | L-J potential ]] energy expression in use is
  
  
 $U(x)=4\epsilon \left [\left ( \frac{\sigma }{r_{ij}} \right )^{12}- \left ( \frac{\sigma }{r_{ij}} \right )^{6}  \right ]$ $U(x)=4\epsilon \left [\left ( \frac{\sigma }{r_{ij}} \right )^{12}- \left ( \frac{\sigma }{r_{ij}} \right )^{6}  \right ]$
  
-where is $\epsilon$ the well depth, σ is related to the minimum of Lennard-Jones potential, rij is the +where $\epsilon$ is the well depth, $\sigma$ is related to the minimum of the Lennard-Jones potential, and $r_{ij}$ is the distance between atoms i and j.
-distance between atoms i and j.+
  
-[[https://en.wikipedia.org/wiki/Periodic_boundary_conditions|Periodic boundary conditions (PBC)]] should be used in this simulation. The atom near the ”edge” of the simulation box interacts with atoms contained in the neighboring image of the box. In computer simulations, one of these is the original simulation box, and others are copies called images. During the simulation, only the properties of the original simulation box need to be recorded and propagated. The minimum-image convention is a common form of PBC particle bookkeeping in which each individual particle in the simulation interacts with the closest image of the remaining particles in the system. To prevent artifacts, it requires a cut-off value of rij for the L-J potential. For realistic results, the cut-off should be less than the half of the simulation box size and over σ Radial distribution +[[https://en.wikipedia.org/wiki/Periodic_boundary_conditions|Periodic boundary conditions (PBC)]] should be used in this simulation. The atom near the ”edge” of the simulation box interacts with atoms contained in the neighboring image of the box. In computer simulations, one of these is the original simulation box, and others are copies called images. During the simulation, only the properties of the original simulation box need to be recorded and propagated. The minimum-image convention is a common form of PBC particle bookkeeping in which each individual particle in the simulation interacts with the closest image of the remaining particles in the system. To prevent artifacts, it requires a cut-off value for $r_{ij}$ of the L-J potential. For realistic results, the cut-off should be less than half of the simulation box size and larger than $\sigma$. 
-function, (or pair correlation function) $g(r)$ in a system of particles (atoms, molecules, colloids, +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.
-etc.), describes how density varies as a function of distance from a reference particle.+
  
 ===== Part I:  Set up MD simulation  ===== ===== Part I:  Set up MD simulation  =====
  
-In this section, a commented CP2K input example for an MD calculation is provided. +In this section, we provide you with an example CP2K input for an MD calculation. 
-Comments are added, they start with a hash mark '#'.+Extensive comments have been added to the file, which start with a has symbol '#'.
  
 === 1. Step === === 1. Step ===
  
-Load the CP2K module as shown before, create a directory ''ex1'' and change to it:+Load the CP2K module as explained in Exercise 0, create a directory ''ex1'' and change to it:
  
 <code> <code>
Line 43: Line 41:
            
 <code - argon.inp> <code - argon.inp>
-##  It's highly recommended to go +##  It's highly recommended to go to 
 ##  https://manual.cp2k.org/  ##  https://manual.cp2k.org/ 
-##  and learn how to set up CP2K  +##  and learn how to set up a CP2K  
-##  calculation correctly using manual.+##  calculation correctly, using the manual.
  
 &GLOBAL &GLOBAL
Line 222: Line 220:
 <note>**TASK** <note>**TASK**
   *Run the calculation and visualize the trajectories using VMD   *Run the calculation and visualize the trajectories using VMD
-  *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 timesteps (100,1000,5000), and inspect geometry in each case, plot the total energies, temperature, potential energies. Try to comment and explain what you observe.+  *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> </note>
 ===== Part II:  Force Field Parameter  ===== ===== Part II:  Force Field Parameter  =====
-You need to modify the L-J force field parameter. Run calculations with different $\sigma$ and $\epsilon$ and inspect the insight of L-J potential (Energies, velocities, temperatures etc ). Then use the force field parameter provided to plot the L-J potential curve.+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       &GLOBAL                  ! section to select the kind of calculation
          RUN_TYPE ENERGY       ! select type of calculation. In this case: ENERGY (=Single point calculation)          RUN_TYPE ENERGY       ! select type of calculation. In this case: ENERGY (=Single point calculation)
Line 251: Line 250:
   &END SUBSYS   &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: 
-any times you will have to run the same simulation with different parameters (here the distance).+  - 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): A simple way to generate the different input files is using shell scripting in combination with ''sed'' (the stream editor):
Line 258: Line 258:
 <code> <code>
 for d in $(seq 2 0.1 4); do for d in $(seq 2 0.1 4); do
-  sed -e "s|4 0 0|${d} 0 0|" argon.inp > energy_${d}A.inp+  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   cp2k.popt -i energy_${d}A.inp -o energy_${d}A.out
   awk '/Total FORCE_EVAL/ { print $9; }' energy_${d}A.out   awk '/Total FORCE_EVAL/ { print $9; }' energy_${d}A.out
Line 268: Line 268:
   * ''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'', ...)   * ''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.   * ... 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 before on those new input files and write the output to new output files as well +  * 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+  * Using ''awk'' we extract the energy from the output file and print it
              
 <note> <note>
 **TASK** **TASK**
-  *Plot the Lennard-Jones potential against Ar-Ar distance $r$ (2-4 Å) using different $\epsilon$ and $\sigma$. +  *Plot the Lennard-Jones potential against the Ar-Ar distance $r$ (2-4 Å) using different values for $\epsilon$ and $\sigma$. 
-  *Repeat the L-J MD calculation with different $\epsilon$ and $\sigma$, and compare the potential energy, temperature evolution, and explain the relation between calculated properties and force field parameters.+  *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.
  
  
 </note> </note>
-===== Part III:  Radial distribution functions  =====+===== Part III:  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.
  
-Use VMD or write your own program (Fortran, C, C++, Python etc.) to calculate radial distribution $g(r)$. Plot $g(r)$, and against various the temperatures to examine the effects. +VMD comes with an extension for exactly this purpose: In the VMD Main window open “Extensions → Analysis” click on “Radial Pair Distribution function $g(r)$“. In the appearing window use “Utilities → Set unit cell dimensions” to tell VMD the size of the simulation box you used. After that use Selection 1 and 2 to define the atomic types that you want to calculate the rdf for, for example “element Ar”. In the plot window, use the "File" menu to save the plot data.
-VMD comes with an extension for exactly this purpose: In the VMD Main window open “Extensions → Analysis” click on “Radial Pair Distribution function $g(r)$. In the appearing window use “Utilities → Set unit cell dimensions” to let VMD know the simulation box you used. After that use Selection 1 and 2 to define the atomic types that you want to calculate the rdf for, for example “element Ar”. In the plot window, use "File", you can save the plot data.+
  
 <note>**TASK** <note>**TASK**
  
   * Plot $g(r)$ at 84, 300 and 400 K into the same graph.   * Plot $g(r)$ at 84, 300 and 400 K into the same graph.
-  * What are the differences in the height of the first peak, and why does temperature contribute to the differences? +  * What are the differences in the height of the first peak, and why/how does the temperature contribute to the differences? 
-  * Compared to experimental data ''exp_gr.dat'' taken at 84 K, what does this say about the structure of the liquid and is this expected? . +  * Compare your results to the experimental data taken at 84 K (given in ''exp_gr.dat'', below). What does this say about the structure of the liquid and did you expect this?
 </note> </note>
 <code - exp_gr.dat> <code - exp_gr.dat>
Line 441: Line 442:
 10.2150    1.1041                                               10.2150    1.1041                                              
 10.2831    1.0977                                               10.2831    1.0977                                              
-10.3512    1~0904                                              +10.3512    1.0904                                              
 10.4193    1.0835                                               10.4193    1.0835                                              
 10.4874    1.0774                                               10.4874    1.0774                                              
Line 454: Line 455:
 </code>                                            </code>                                           
  
-===== Part IV:  Ensembles  =====+===== Part IV: Other Ensembles  =====
  
-In previous section, you have already run NVE ensemble molecular dynamics for Ar liquid. In this section, we will focus on the NVT, NPT 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.
  
-Step up NVT calculation, change the setting in &MD section.  +To set up an NVT calculation, change the settings in the &MD section as shown below:
  
  
Line 475: Line 476:
      
      
-Step up NPT calculation, change the setting in &MD section.  +To set up an NPT calculation, change the settings in the &MD section as shown below:
  
   &FORCE_EVAL   &FORCE_EVAL
Line 502: Line 503:
 <note> <note>
 **TASK** **TASK**
-   *Run calculation using NVT at 300K, check the temperature, and energy of the whole system, and compare the result to NVE (300K), and rationalize the difference.  +   *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 calculation using NVT (300K) until the system is equilibrated then run NVE, check the temperature, and energy of the whole system, and  compare to previous NVE simulation. +   *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 NPT simulation then check the size of the simulation box.+   *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:
  
 +1. 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>
  
 +2. 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> </note>
exercises/2021_uzh_acpc2/ex01.1619419107.txt.gz · Last modified: 2021/04/26 06:38 by mrossmannek