Free and constrained molecular dynamics

We have prepared an input file free.in to perform an unconstrained molecular dynamics of $\text{NaCl}$ dissolved in water.

TASK 1

Ideally, we would like to study isolated $\text{NaCl}$ ion pairs. Determine the concentration of $\text{NaCl}$ in the simulation in units of mol per liter. Do you consider the solution to be dilute?

When collecting statistics by successive MD runs, one often wants initial conditions to be random in order for the runs to be statistically independent. At the same time, it is often desirable for MD simulations to be exactly reproducible.

This is, where pseudorandom number generators (PRNGs) come in. A PRNG produces a sequence of numbers that seems random (sampling e.g. the normal distribution), but is in fact completely determined by the so-called seed given to the PRNG. Given the same seed, the PRNG produces exactly the same sequence of numbers.

Here, we use a pseudorandom number generator to initialize the velocities of the atoms. The velocities are drawn from the Boltzmann statistics, ensuring that the corresponding kinetic energy corresponds to the desired temperature.

TASK 2

Run MD simulations of $\text{NaCl}$ in water using three different SEED values for the PRNG. Note: You may want to run the simulations in different directories to avoid overwriting files. We have already provided results for the seeds 2,4,6,7,8,9,12,13,15,17, each spanning 100 ps.
Use VMD to compute the distance between $\text{Na}$ and $\text{Cl}$ in periodic boundary conditions (see tip below). Can you identify different regimes? (2P)
Compute the average time needed for the dissociation of $\text{NaCl}$. Do you observe re- association of the ions? Note: In order to set a standard, let's count a dissociation, when the ion separation increases beyond 5 $\unicode{x212B}$.

In order to compute the Na-Cl distance within vmd, use the Tcl interface

pbc set { lx ly lz } -all # define simulation box lx*ly*lz
# center the box around Na and wrap all atoms into the box
pbc wrap -centersel "element Na" -center com -all

You may use the dynamic bonds representation to recompute the bonds after the wrapping procedure. Use a different representation for NaCl.

Hint: In order to save time with the next trajectory, save the visualization state and reuse it. The pbc commands are not save into the visualization state file, but you can edit the file with a text editor and append them at the end.

Next, we are going to constrain the Na-Cl distance in the MD simulation to a value of X.Y $\unicode{x212B}$, where X.Y ranges from 2.5 to 7.0. We have already performed the MD simulations and computed the radial distribution functions.

TASK 3

Use gnuplot to analyze the rdfs as a function of the constraint. Note: You may use either the individual gofr_A_B_X.Y or the joined files gofr_A_B (here, use splot).
Describe the effect of the Na-Cl distance on the different rdfs. Note: The radial distribution function computed here only considers atom pairs A and B from different molecules. (2P)
How is the coordination number $n_1$ defined? Give an analytical formula to compute it from the rdf.
Calculate the first coordination number of Na from the Na-O rdf using the script ./integrate.py --N=<Natoms> --L=<BoxLength> < gofr_A_B-X.Y > nc_A_B-X.y where <Natoms> is the number of B atoms per cell and <BoxLength> the length of the simulation box in $\unicode{x212B}$.
How does the coordination number of Na with respect to O vary with the constrained Na-Cl distance? What does this mean for the coordination of Na? How to explain the linear increase of Na-O coordination number after 4 angstroms (should be cubic, right?)? (2P)