====== 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. * 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''. ===== 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: Shake Lagrangian Multipliers: -0.054769270 Rattle Lagrangian Multipliers: -0.020937479 Shake Lagrangian Multipliers: -0.020937479 Rattle Lagrangian Multipliers: -0.020937479 ... * From these files you can calculate the average Lagrange multiplier of the Shake-algorithm like this: grep Shake NACL-XXX.LagrangeMultLog | awk '{c++ ; s=s+$4}END{print s/c}' * 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} Make sure that you get the units right. The Largange multipliers are written in atomic units (Hartree/bohr), while the distances are in Angstrom. 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. &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 ==== Motion section TO ADD for constrained MD ==== This section has to be added to the above input file for ** Task 2 ** and ** Task 3 ** &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 ==== Average Largange multiplier for NaCl in water at 350K (incomplete) ==== This is the Lagrange Multipliers table to be completed for ** Task 3 ** # 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