CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
exercises:2021_uzh_acpc2:ex03

Differences

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

Link to this comparison view

Next revision		Previous revision
exercises:2021_uzh_acpc2:ex03 [2021/04/13 11:02] – external edit 127.0.0.1exercises:2021_uzh_acpc2:ex03 [2021/05/21 10:13] (current) – [Free energy surface] Add tip about RuntimeWarning mrossmannek
Line 1: Line 1:
 ====== Nudged elastic band and free energy calculations ====== ====== Nudged elastic band and free energy calculations ======
  
-In this exercise we will study the $S_N2$ nucleophilic substitution reaction +In this exercise we will study the $S_N2$ nucleophilic substitution reaction: 
  
 $$ \text{Cl}^- + \text{CH}_3\text{Cl} \longleftrightarrow \text{Cl}\text{CH}_3 + \text{Cl}^-. $$ $$ \text{Cl}^- + \text{CH}_3\text{Cl} \longleftrightarrow \text{Cl}\text{CH}_3 + \text{Cl}^-. $$
  
-For energy and force evaluation, we will use the Parameterization Method 6 (PM6), which is a relatively inexpensive semiempirical model for electronic structure evaluation. For accurate reaction characterizations, ab initio methods, such as DFT or higher level, should be used.+For the energy and force evaluation, we will use the Parameterization Method 6 (PM6), which is a relatively inexpensive semiempirical model for electronic structure evaluation. For accurate reaction characterizations, ab initio methods, such as DFT or higher level, should be used instead.
  
 ===== NEB activation barrier ===== ===== NEB activation barrier =====
  
-In order to find the activation barrier for the reaction we will use the NEB method. The NEB method provides a way to find the minimum energy path (MEP) between two different conformations of the system corresponding local potential energy minima. The MEP is a one-dimensional path on the $3N$-dimensional potential energy landscape and every point on the MEP is a potential energy minimum in all directions perpendicular to the path. MEP passes at least one saddle point and the energy of the highest saddle point is the peak of the activation barrier of the reaction.+In order to find the activation barrier for the reaction we will use the NEB method. The NEB method provides a way to find the minimum energy path (MEP) between two different conformations of the system which correspond to local potential energy minima. The MEP is a one-dimensional path on the $3N$-dimensional potential energy landscape for which every point along its path is a potential energy minimum in all directions **perpendicular** to the path. Any MEP passes through at least one saddle point and the energy of the highest saddle point is the peak of the activation barrier of the reaction.
  
 === Geometry optimizations === === Geometry optimizations ===
Line 76: Line 76:
 <note>**TASK 1** <note>**TASK 1**
   *Run the geometry optimization calculation   *Run the geometry optimization calculation
-  *Modify the initial geometry guess in the input and run a second optimization to obtain the other local minimum geometry (The other Cl needs to make a covalent bond with the C). +  *Modify the initial geometry guess in the input and run a second optimization to obtain the other local minimum geometry (The other ''Cl'' needs to make a covalent bond with the ''C'') 
 +</note> 
 + 
 +<note tip> 
 +When you run the second geometry optimization you must make sure that you do not overwrite the results of your previous calculation! The easiest way to achieve this is to run the second simulation in another directory. 
 + 
 +Another possibility is to change the ''PROJECT'' name in the input file which will result in different names for the output files.
 </note> </note>
  
 === NEB === === NEB ===
  
-Now we are ready to start the NEB. Obtain the two optimized geometries from the GEO_OPT trajectories (last step in the ''ch3cl-pos-1.xyz'' file) to a separate folder and name them ''init.xyz'' and ''final.xyz''. +Now we are ready to start the NEB calculation. Obtain the two optimized geometries from the ''GEO_OPT'' trajectories (the last step in the corresponding ''*-pos-1.xyz'' file) and place them in a separate folder in files named ''init.xyz'' and ''final.xyz'', respectively. 
-Now the NEB calculation can be run by the following input script:+Now, the NEB calculation can be run by the following input script:
  
 <code - neb.inp> <code - neb.inp>
Line 106: Line 112:
     &END     &END
     &REPLICA     &REPLICA
-      COORD_FILE_NAME init.xyz+      COORD_FILE_NAME init.xyz       # Filename of the initial coordinate file
     &END     &END
     &REPLICA     &REPLICA
-      COORD_FILE_NAME final.xyz+      COORD_FILE_NAME final.xyz      # Filename of the final coordinate file
     &END     &END
     &PROGRAM_RUN_INFO     &PROGRAM_RUN_INFO
Line 128: Line 134:
     &END CELL     &END CELL
     &TOPOLOGY     &TOPOLOGY
-      COORD_FILE_NAME init.xyz+      COORD_FILE_NAME init.xyz       # Filename of the initial coordinate file
       COORDINATE xyz       COORDINATE xyz
     &END TOPOLOGY     &END TOPOLOGY
Line 151: Line 157:
  *******************************************************************************  *******************************************************************************
 </code> </code>
-These sections show for every replica geometry along the NEB trajectory, the distance to its neighbors and its energy. The final section corresponds to the converged NEB trajectory.+These sections show, for every replica geometry along the NEB trajectory, the distance to its neighbors and its energy. The last of these sections corresponds to the converged NEB trajectory.
  
 <note>**TASK 2** <note>**TASK 2**
   * Run the NEB calculation   * Run the NEB calculation
-  * Find the activation barrier of the reaction in eV.+  * Find the activation barrier of the reaction **in eV**.
 </note> </note>
  
Line 163: Line 169:
 Sampling the free energy surface (FES) of a chemical system is a convenient method to explore various stable conformations and possible reaction pathways. To calculate the FES for complicated systems, advanced sampling methods (such as umbrella sampling, metadynamics, parallel tempering, ...) have to be used. However, for our simple $S_N2$ reaction, we will use unbiased Molecular Dynamics (MD).  Sampling the free energy surface (FES) of a chemical system is a convenient method to explore various stable conformations and possible reaction pathways. To calculate the FES for complicated systems, advanced sampling methods (such as umbrella sampling, metadynamics, parallel tempering, ...) have to be used. However, for our simple $S_N2$ reaction, we will use unbiased Molecular Dynamics (MD). 
  
-The FES is a projection of the high-dimensional free energy landscape into a small number, usually two, dimensions. These two dimensions are called collective variables (CV) and they must be chosen such that various stable conformations can be distinguished and the reaction pathways can be adequately described by the FES. For complex systems, the choice of CVs is a non-trivial task. Fortunately for our simple system, the choice is simple: we take the distances of the two Cl anions from the central C as the CVs.+The FES is a projection of the high-dimensional free energy landscape onto a small number, usually two, dimensions. These two dimensions are called collective variables (CV) and they must be chosen such that various stable conformations can be distinguished and the reaction pathways can be adequately described by the FES. For complex systems, the choice of CVs is a non-trivial task. Fortunately, for our simple system, the choice is simple: we take the distances of the two ''Cl'' anions from the central ''C'' as the CVs.
  
-To help the calculation sample the parts of FES we care about, we will include restraints in the MD run, which prevent the Cl anions from going too far from the molecule.+To help the simulation sample the parts of the FES which we care about, we will include restraints in the MD run, which prevent the ''Cl'' anions from going too far away from the molecule.
  
 The following CP2K input script runs our MD calculation and prints out the CV values for every step: The following CP2K input script runs our MD calculation and prints out the CV values for every step:
Line 172: Line 178:
 &GLOBAL &GLOBAL
   PRINT_LEVEL LOW   PRINT_LEVEL LOW
-  PROJECT ch3f+  PROJECT ch3cl
   RUN_TYPE MD                 # Molecular Dynamics   RUN_TYPE MD                 # Molecular Dynamics
 &END GLOBAL &END GLOBAL
Line 268: Line 274:
 $$ F(s) = -k T \log(P(s)),$$ $$ F(s) = -k T \log(P(s)),$$
  
-where $s$ is the set CVs and $P(s)$ is the probability that the system has the set of CV values $s$.+where $s$ is the set of CVs and $P(s)$ is the probability that the system has the set of CV values $s$.
  
-The following Python script can be used to calculate the FES from the ''ch3f-COLVAR.metadynLog'' file produced by the MD run. (Don't forget to change the temperature in Python script!)+The following Python script can be used to calculate the FES from the ''ch3f-COLVAR.metadynLog'' file produced by the MD run. (**Don't forget to change the temperature in the Python script!**)
  
 <code> <code>
Line 279: Line 285:
 kb = 8.6173303e-5 # eV * K^-1 kb = 8.6173303e-5 # eV * K^-1
  
-temperature = 1000.0                          #Change temperature according to your MD simulations!+temperature = 1000.0                          # Change the temperature according to your MD simulations!
 colvar_path = "./ch3f-COLVAR.metadynLog" colvar_path = "./ch3f-COLVAR.metadynLog"
  
Line 310: Line 316:
 </code> </code>
  
-Here is an example output for temperature 1000K. We clearly see the two local minima corresponding to the one of the chlorine atoms being covalently bonded (distance 1.8 Å) while the other one is around distance 2.5 Å.+Here is an example output for a temperature of 1000K. We clearly see the two local minima corresponding to one of the chlorine atoms being covalently bonded (distance 1.8 Å) while the other one is around a distance of 2.5 Å.
  
 {{  :exercises:2019_uzh_acpc2:fes1000.png ?direct&500 |}} {{  :exercises:2019_uzh_acpc2:fes1000.png ?direct&500 |}}
  
 <note>**TASK 3** <note>**TASK 3**
-  * Run the MD calculation for 400K, 800K, 1200K and 1600K. (The calculations can a take a while.)+  * Run the MD calculation for 400K, 800K, 1200K and 1600K. (These calculations can a take a while.)
   * Create the corresponding FES plots and discuss the temperature dependence.   * Create the corresponding FES plots and discuss the temperature dependence.
   * In general, how does potential energy differ from free energy? For our reaction, what are the activation barriers from the different free energy surfaces? How and why do they differ from the NEB barrier?   * In general, how does potential energy differ from free energy? For our reaction, what are the activation barriers from the different free energy surfaces? How and why do they differ from the NEB barrier?
 +</note>
 +
 +<note tip>
 +You may encounter a warning like the following:
 +<code>
 +plot.py:24: RuntimeWarning: divide by zero encountered in log
 +  fes = -kb * temperature * np.log(prob)
 +</code>
 +Don't worry about this, the script still works as expected and produces the plot in the file ''fes.png''.
 </note> </note>
exercises/2021_uzh_acpc2/ex03.1618311776.txt.gz · Last modified: 2021/04/13 11:02 by 127.0.0.1