Ions perturb water structure

FIg. 1.  Nearby water molecules orient their negatively charged oxygens towards a Na+ ion.

Ions in solution perturb the arrangements of nearby water molecules, both sterically and electrostatically.  Ions take up space, which in itself leads to a nonuniform average concentration of nearby water molecules, with an exclusion radius, a shell of near neighbors, and decaying oscillations as the local concentration relaxes to its average value.  And the ionic charge orients nearby water molecules, with positive ions like Na+ tending to attract the negatively charged oxygen atoms and repel positively charged hydrogen atoms on water. 

The electrostatic potential is likewise perturbed by the presence of an ion.  In vacuum, the ion would have a Coulomb field, decaying as 1/r^2.  In water, the oriented dipoles of nearby water molecules tends to screen the field of the ion, as the oppositely charged atoms collect in a shell around the ion and reduce its net charge. 

We can use MD simulations to explore this screening behavior of water in the vicinity of ions,  by performing a simulation in which a single Na+ ion is fixed in the center of a box of water.  From the simulation trajectory, we would like to compute the average electric field, dipole orientation, charge density, and distribution of molecules near the fixed ion. 

Custom measurements and analysis

Gromacs includes many useful analysis utilities.  However, not every measurement we dream up can be made with the in-built utilities.  Here, we demonstrate techniques for a) measuring quantities not provided by Gromacs analysis routines, b) exporting molecular trajectories, and c) writing analysis programs in Mathematica.  

We use Mathematica here because it is a high-level language, in which complicated tasks can be performed in a few lines of code.  Writing in a high-level language like Mathematica is the fastest way to get a working program.  For all but the most computationally intensive tasks, Mathematica code is fast enough for our purposes, compared to more computationally efficient languages like C++ that take much longer to write.

Electric field from “test charges”

We want to measure the average electric field near the ion.  By symmetry, we expect the field to point in the radial direction.  Gromacs provides no utility for this purpose.  We obtain the electric field from the simulation by means of a trick:  we add charge to the oxygen on a given water molecule; rerun the simulation trajectory using -rerun; output the molecular position and force versus time; and take the difference of force with and without the added charge.  In effect, we have added a “test charge” to the oxygen that does not perturb the trajectory,  because we rerun the trajectory obtained without the added charge.  The extra force F equals qE, so we have measured the field E wherever the oxygen goes.  To measure the field throughout the system, we repeat this for many randomly chosen water molecules, saving files of position and force for each.     

Export and analyze molecular trajectories

From the molecular positions, we can also compute the dipole orientation, charge density, and distribution of molecules near the fixed ion.  We want these quantities as a function of distance from the ion; however, there are no Gromacs utilities for this.  So we must export atomic coordinates from the trajectory, and analyze them outside Gromacs.  Since we already write trajectories for the water molecules with “test charges”, we use these to obtain the dipole orientation, charge density, and molecular concentration near the ion.

[Page 2:  “how-to” help.]