User Tools

Site Tools


exercises:2017_ethz_mmm:nacl_free_energy

Differences

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


exercises:2017_ethz_mmm:nacl_free_energy [2020/08/21 10:15] (current) – created - external edit 127.0.0.1
Line 1: Line 1:
 +====== Free Energy Profile of NaCl Dissociation======
 +
 +In this exercise, you will run different simulations to compute the NaCl dissociation curve in both gas and solution environments.
 +
 +<note tip>
 +  * You'll have to run many similar simulations. Try to automatize as much as possible (we can help you).
 +  * To avoid confusion, try to perfrom every task in a new directory 
 +  * 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>
 +
 +===== 1. Task: Potential energy curve (gas phase) =====
 +This case is very similar to the computation of the Lennard Jones curve (See: [[exercises:2015_ethz_mmm:single_point_calculation|Computation of the Lennard Jones curve]] ). \\ 
 +  * For this you have to run the input file ''NaCl_gasphase.inp'' at a range of Na-Cl distances. This can be automathized, so we provide with an template and you have to vary the MYDIST parameter in the input. 
 +  * At the end, plot the potential energy dissociation profile of NaCl.
 +
 +===== 2. Task: Free energy curve at 1K (gas phase) =====
 +
 +For this you have to run constrained MD simulations at 1K for a range of Na-Cl distances. 
 +  * You have to modify the input file in the following way:
 +    - Copy the ''NaCl_gasphase.inp'' file to a new directory and rename it to something like: ''NaCl_MD.inp''.
 +    - Change the ''RUN_TYPE''  in the new input file, from "ENERGY" to "MD".
 +    - Add the ''MOTION''-section provided (end of this page) to the new ''NaCl_MD.inp'' file. 
 +
 +  * Then, as usual, run the simulation for a range of NaCl distances. This is a constrained MD simulation, meaning that you have to vary the MYDIST parameter at three points in the file:
 +     - In the COORD section of the new ''NaCl_MD.inp'' file
 +     - In the CONSTRAINT section of the new ''NaCl_MD.inp'' file
 +     - Where the PROJECT_NAME keyword is
 +
 +⇒ Each constrained MD will produce a ''.LagrangeMultLog''-files, which look like this:
 +<code>
 +Shake  Lagrangian Multipliers:            -0.054769270
 +Rattle Lagrangian Multipliers:            -0.020937479
 +Shake  Lagrangian Multipliers:            -0.020937479
 +Rattle Lagrangian Multipliers:            -0.020937479
 +...
 +</code>
 +
 +  * From these files you can calculate the average Lagrange multiplier of the Shake-algorithm like this:
 +<code>
 +grep Shake NACL-XXX.LagrangeMultLog | awk '{c++ ; s=s+$4}END{print s/c}'
 +</code>
 +
 +  * The average Lagrange multiplier is the average force $F(x)$ required to constrain the atoms at the distance $x$.
 +  * From these forces the free energy difference can be obtained via integration:
 +\begin{equation}
 +\Delta A = -\int_a^b F(x)\, dx
 +\end{equation}
 +
 +  * The 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 =====
 +In this section, we provide an incomplete list of average Lagrange multipliers. You will have to run a single constrained MD simulation, get the average Lagrange Multiplier. In this way you can complete the list and compute the free energy profile in water. 
 +
 +  * Use the same input as ** Task 2**.
 +  * BUT take the forcefield and the solvated system coordinates from the previous exercise (See: [[exercises:2015_ethz_mmm:nacl_md|Observe NaCl dissociation]]). In practice, you have to substitue the whole '' FORCE EVAL '' section in the ** Task 2** input with the '' FORCE EVAL '' section of [[exercises:2015_ethz_mmm:nacl_md|Observe NaCl dissociation]].
 +  * Other slight modifications to your input: 
 +         - In the MOTION-CONSTRAINT section set TARGET to 2.9.  
 +         - In the MOTION-MD section set STEPS 100.000 MD and T 350. 
 +  * Run the simulation in the same way you did for ** Task 2**.
 +  * From the MD output calculate the average Largange multiplier,in the same way you did for ** Task 2**. 
 +  * Complete the Lagrange Multiplier list we've provided (end of this page)
 +  * From the complete table calculate the free energy dissociation profile via numerical integration. 
 +
 +===== Required Files =====
 +
 +
 +
 +==== Input file for NaCl in gasphase ====
 +This is the basic input.
 +Note that for **Task 2** and ** Task 3 ** it should be modified.
 +<code - NaCl_gasphase.inp>
 +&FORCE_EVAL
 +  METHOD FIST
 +  &MM
 +    &FORCEFIELD
 +      &SPLINE
 +        EPS_SPLINE 1.0E-8
 +        EMAX_SPLINE 300000.0
 +      &END
 +      &CHARGE
 +        ATOM Na
 +        CHARGE 1.0
 +      &END CHARGE
 +      &CHARGE
 +        ATOM Cl
 +        CHARGE -1.0
 +      &END CHARGE
 +      &NONBONDED
 +        &LENNARD-JONES
 +          atoms Na Cl
 +          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
 +      &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 12.4138 12.4138 12.4138
 +    &END CELL
 +     &COORD
 +Na    0.0    0.0 0.0 NAP
 +Cl    MYDIST 0.0 0.0 CLM
 +     &END COORD
 +     &COLVAR
 +      &DISTANCE
 +        ATOMS 1 2
 +      &END DISTANCE
 +      &PRINT
 +      &END
 +     &END COLVAR
 +    &TOPOLOGY
 +      CONNECTIVITY GENERATE
 +      &GENERATE
 +        BONDLENGTH_MAX 7
 +      &END
 +    &END
 +  &END SUBSYS
 +&END FORCE_EVAL
 +
 +&GLOBAL
 +  PROJECT NACL-MYDIST
 +  RUN_TYPE ENERGY
 +&END GLOBAL
 +</code>
 +
 +==== Motion section TO ADD for constrained MD ====
 +This section has to be added to the above input file for ** Task 2 ** and ** Task 3 **
 +<code - motion section>
 +&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) ====
 +This is the Lagrange Multipliers table to be completed for ** Task 3 **
 +<code>
 +# dist     avg. Shake Lagrange multiplier
 +  2.5         0.0896372
 +  2.6         0.0469698
 +  2.7         0.0231717
 +  2.8         0.0100625
 +  2.9                          <--- Take missing value from your trajectory
 +  3.0        -0.000996937
 +  3.1        -0.00271078
 +  3.2        -0.00335324
 +  3.3        -0.00348111
 +  3.4        -0.00303697
 +  3.5        -0.00259636
 +  3.6        -0.00201541
 +  3.7        -0.00119027
 +  3.8        -0.000408723
 +  3.9        -8.19056e-05
 +  4.0         0.000972204
 +  4.1         0.00136578
 +  4.2         0.0016246
 +  4.3         0.00212447
 +  4.4         0.00199128
 +  4.5         0.00183284
 +  4.6         0.00188221
 +  4.7         0.00166909
 +  4.8         0.00137179
 +  4.9         0.00114308
 +  5.0         0.000671159
 +  5.1         0.000780625
 +  5.2         0.000556307
 +  5.3         0.000397211
 +  5.4         0.000237853
 +  5.5         0.000119549
 +  5.6        -0.000220194
 +  5.7        -0.000332539
 +  5.8        -0.000674227
 +  5.9        -0.00075852
 +  6.0        -0.00043128
 +</code>