Differences

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

--- exercise:nacl_free_energy [2014/05/28 12:25] – exercise:13_2 renamed to exercise:nacl_free_energy oschuett
+++ exercise:2014_ethz_mmm:nacl_free_energy [2014/10/15 13:36] – oschuett
@@ Line 1: / Line 1: @@
-====== Dissociation Profile of NaCl ======
+====== Free Energy Profile of NaCl Dissociation======
-Plot a gas phase dissociation profile of NaCl based on potential energy using the input file ''NaCl_gasphase.inp''.
+<note tip>
+  * You'll have to run many similar simulations. Try to automatize as much as possible.
+  * The first two task can be run directly on the login node, i.e. without using bsub.
+  * The third task should be run on 4 cores with ''bsub -n 4''.
+</note>
-Plot a gas phase dissociation profile of NaCl based on thermodynamic integration
+===== 1. Task: Potential energy curve (gas phase) =====
+Plot the gas phase dissociation profile of NaCl of the potential energy. For this you have to run the input file ''NaCl_gasphase.inp'' at a range of Na-Cl distances.
- using the input file ''NaCl_gasphase.inp''.
-The constrained MD runs will produce output ''.LagrangeMultLog''-files, which look like this:
+===== 2. Task: Free energy curve at 1K (gas phase) =====
+Plot the gas phase dissociation profile of NaCl of the free energy at 1K.
+For this you have to run constrained MD simulations at 1K for a range of Na-Cl distances. You have to add the ''MOTION''-section provided below to the ''NaCl_gasphase.inp'' and change the ''RUN_TYPE''.
+Each constrained MD will produce a ''.LagrangeMultLog''-files, which look like this:
 <code>
 Shake  Lagrangian Multipliers:            -0.054769270
@@ Line 16: / Line 25: @@
 </code>
-For the thermodynamic integration the average values of the Shake Lagrange multiplier has to be calculate.
+From these files you can calculate the average Lagrange multiplier of the Shake-algorithm like this:
 <code>
-grep Shake NACL-DIMER-XXX.LagrangeMultLog | awk '{c++ ; s=s+$4}END{print s/c}'
+grep Shake NACL-XXX.LagrangeMultLog | awk '{c++ ; s=s+$4}END{print s/c}'
 </code>
-<note tip>
+The average Lagrange multiplier is the average force $F(x)$ required to constrain the atoms at the distance $x$.
-Running many similar simulations is tedious. Try to automatize as much as possible.
+From these forces the free energy difference can be obtained via integration:
+\begin{equation}
+\Delta A = -\int_a^b F(x)\, dx
+\end{equation}
+A dissociation profile can be obtained by choosing the closest distance $d_{min}$ as lower integration-bound:
+\begin{equation}
+A(d) = -\int_{d_{min}}^d F(x)\, dx
+\end{equation}
+<note warning>
+Make sure that you get the units right. The Largange multipliers are written in atomic units (Hartree/bohr), while the distances are in Angstrom.
 </note>
+Compare the free-energy dissociation curve at 1K with the potential energy curve. What do you expect? What do you observe?
+===== 3. Task: Free energy curve of NaCl in water at 350K =====
+Take the solvated system from the [[nacl_md | first exercise]] and add the constraint for a distance of 2.9 Å. Then run 100.000 MD steps MD at 350K. From the MD output calculate the average Largange multiplier. As a check for convergence you can divide the trajectory into two parts and calculate the average for each part separately. Once you are convinced of the result you can use it to complete the table given below. From the complete table calculate the free energy dissociation profile via numerical integration.
 ===== Required Files =====
-<code - motion.inp>
-&MOTION
- &CONSTRAINT
-    &COLLECTIVE
+==== Input file for NaCl in gasphase ====
-      COLVAR 1
-      INTERMOLECULAR
-      TARGET [angstrom] MYDIST
-    &END COLLECTIVE
-    &LAGRANGE_MULTIPLIERS
-      COMMON_ITERATION_LEVELS 1
-    &END
- &END CONSTRAINT
- &MD
-   ENSEMBLE  NVE
-   TIMESTEP  0.5
-   STEPS      100
-   TEMPERATURE 1
-   &THERMOSTAT
-     &NOSE
-       LENGTH 3
-       YOSHIDA 3
-       TIMECON 1000
-       MTS 2
-     &END NOSE
-   &END THERMOSTAT
-   &PRINT
-     &TRAJECTORY OFF
-     &END TRAJECTORY
-     &VELOCITIES OFF
-     &END VELOCITIES
-     &FORCES OFF
-     &END FORCES
-     &RESTART_HISTORY OFF
-     &END RESTART_HISTORY
-     &RESTART OFF
-     &END RESTART
-   &END PRINT
-&END MOTION
-</code>
 <code - NaCl_gasphase.inp>
@@ Line 87: / Line 80: @@
           EPSILON [kcalmol]  .0838
           SIGMA   [angstrom] 3.63
+          RCUT    [angstrom] 11.4
+        &END LENNARD-JONES
+        &LENNARD-JONES
+          atoms Na Na
+          EPSILON [kcalmol]  0.0469
+          SIGMA   [angstrom] 2.7275
+          RCUT    [angstrom] 11.4
+        &END LENNARD-JONES
+        &LENNARD-JONES
+          atoms Cl Cl
+          EPSILON [kcalmol]  0.150
+          SIGMA   [angstrom] 4.54
           RCUT    [angstrom] 11.4
         &END LENNARD-JONES
@@ Line 108: / Line 113: @@
 Cl    MYDIST 0.0 0.0 CLM
      &END COORD
+     &COLVAR
+      &DISTANCE
+        ATOMS 1 2
+      &END DISTANCE
+      &PRINT
+      &END
+     &END COLVAR
     &TOPOLOGY
       CONNECTIVITY GENERATE
@@ Line 121: / Line 133: @@
   RUN_TYPE ENERGY
 &END GLOBAL
+</code>
+==== Motion section for constrained MD ====
+<code - motion.inp>
+&MOTION
+ &CONSTRAINT
+    &COLLECTIVE
+      COLVAR 1
+      INTERMOLECULAR
+      TARGET [angstrom] MYDIST
+    &END COLLECTIVE
+    &LAGRANGE_MULTIPLIERS
+      COMMON_ITERATION_LEVELS 1
+    &END
+ &END CONSTRAINT
+ &MD
+   ENSEMBLE NVT
+   TIMESTEP 0.5
+   STEPS    100
+   TEMPERATURE 1
+   &THERMOSTAT
+     &NOSE
+       LENGTH 3
+       YOSHIDA 3
+       TIMECON 1000
+       MTS 2
+     &END NOSE
+   &END
+   &PRINT
+     &ENERGY OFF
+     &END ENERGY
+     &PROGRAM_RUN_INFO OFF
+     &END PROGRAM_RUN_INFO
+   &END PRINT
+ &END MD
+ &PRINT
+  &TRAJECTORY OFF
+  &END
+  &VELOCITIES OFF
+  &END VELOCITIES
+  &FORCES OFF
+  &END FORCES
+  &RESTART_HISTORY OFF
+  &END RESTART_HISTORY
+  &RESTART OFF
+  &END RESTART
+ &END PRINT
+&END MOTION
 </code>
+==== Average Largange multiplier for NaCl in water at 350K (incomplete) ====
+<code>
+# dist     avg. Shake Lagrange multiplier
+.5         0.0896372
+.6         0.0469698
+.7         0.0231717
+.8         0.0100625
+.9                          <--- Take missing value from your trajectory
+.0        -0.000996937
+.1        -0.00271078
+.2        -0.00335324
+.3        -0.00348111
+.4        -0.00303697
+.5        -0.00259636
+.6        -0.00201541
+.7        -0.00119027
+.8        -0.000408723
+.9        -8.19056e-05
+.0         0.000972204
+.1         0.00136578
+.2         0.0016246
+.3         0.00212447
+.4         0.00199128
+.5         0.00183284
+.6         0.00188221
+.7         0.00166909
+.8         0.00137179
+.9         0.00114308
+.0         0.000671159
+.1         0.000780625
+.2         0.000556307
+.3         0.000397211
+.4         0.000237853
+.5         0.000119549
+.6        -0.000220194
+.7        -0.000332539
+.8        -0.000674227
+.9        -0.00075852
+.0        -0.00043128
+</code>