Profiles of potential energy and free energy
We are going to start with the simplified example of isolated $\text{Na}^+$ and $\text{Cl}^-$ in the gas phase, where we can directly compare the results of our computer simulation against the analytical formulae used to describe the interaction potential.
We have provided an input file NaCl_pot.in and a script potential_energy.sh that uses this input file to calculate the potential energy as a function of Na-Cl distance.
- Look into
NaCl_pot.inand write down the formula used for the potential energy of the interaction between $\text{Na}^+$ and $\text{Cl}^-$ in Hartree atomic units. (2P) - Use
./potential_energy.shto calculate the potential energy as a function of Na-Cl distance. Create a plot of the resulting potential energy profile inpot_profileand the mathematical formula. - What do you observe, when the distance approaches 1/2 of the simulation box? How might the finite size of the simulation box have impacted the MD simulation in the previous exercise? (2P)
For the next task, we remain with our simple system, but now perform molecular dynamics at $T=1\,\text{K}$.
We have prepared a script free_energy.sh, which runs MD simulations with constrained Na-Cl distance at $1\,\text{K}$.
It then integrates the average value of the Shake Lagrange multiplier to calculate the (low-temperature) free energy profile.
- What is a Lagrange multiplier? How can we obtain the free energy profile as a function of the Na-Cl distance using the associated Lagrange multiplier? (2P)
- Run the simulation. What kind of motion does the NaCl dimer perform?
- Compare the low-temperature free energy profile in
fe_profilewith the potential energy profile. Do the two profiles agree? Note: The profiles are shifted with respect to each other. What would be a reasonable reference point for both profiles? (2P) - What effects would you expect at higher temperature? Hint: If you like, you can adapt the temperature in the input file and give it a go.
Now, we are ready to move to a more realistic system – NaCl in water. We have performed constrained MD of NaCl in water and saved the trajectory of the corresponding Lagrange multipliers (ask your teaching assistant).
The script ./integrate.sh computes the average values of the Shake Lagrange multipliers and uses them to perform the free energy integration.
- Perform the free energy integration and plot the free energy profile.
- In the previous exercise, you determined the average time required for dissociation of Na-Cl. Is the free energy barrier consistent with the time scale determined before? Hint: Use the Arrhenius equation. You can obtain an estimate for the attempt frequency from the high-frequency oscillations in the Na-Cl distance in the previous exercise. (2P)
Another way to gain access to the free energy is through the radial distribution function (rdf) of the unconstrained system. The rdf $g(r)$ is related to the free energy $F(r)$ through the following set of equations $$\begin{eqnarray} g(r)4\pi r^2 &\propto& \int \delta(r-r') \exp(-\beta H(r'))\,dr \\ P(r) &\propto& \int \delta(r-r') \exp(-\beta H(r'))\,dr \\ F(r) &=& -k_BT \ln\,P(r) \end{eqnarray}$$
We have performed a trajectory spanning 50 ns of unconstrained molecular dynamics of NaCl in water (ask your teaching assistant). The individual frames are spaced by 1 ps in order to reduce correlation between subsequent frames.
- In the previous exercise, we computed the O-O radial distribution function for water with acceptable statistics using just 20 ps of simulated time. Give two reasons, why collecting enough statistics for the Na-Cl radial distribution function requires much longer simulation times (with our setup).
- Compute the radial distribution function for the provided trajectory and plot it as a function of Na-Cl distance.
- Use the equations above to compute the free energy profile. Does it agree with the one constructed from the Shake Lagrange multipliers?