Fortran Code to Calculate the Hydration Shell of an Ion in Water

(Citation: Instantaneous, Parameter-Free Methods to Define a Solute's Hydration Shell, A. Chatterjee, J. Higham and Richard H Henchman, J. Chem. Phys., 2015, 143, 234501.)

This works for a single-atom solute (cation, anion or noble gas solute) from a GROMACS simulation. Note that not all methods work for all kinds of solute. See text for details.

Fortran code files:

(detailed descriptions of what all codes do are given as comments in the fortran files):

  • LJO/SAHB.f: Lennard-Jones Overlap/Strongest-Acceptor Hydrogen Bond (LJO/SAHB) definition for cations and noble gas solutes only.
  • NN_anion.f: Nearest Neighbors (NN) definition for anions.
  • NN_cation_ng.f: Nearest Neighbors (NN) definition for cations and noble gases.
  • NNNW_anion.f: Nearest Neighbor’s Nearest Water (NNNW) definition for anions.
  • NNNW_cation_ng.f: Nearest Neighbor’s Nearest Water (NNNW) definition for cations and noble gas solutes.
  • SAHBanion.f: Strongest-Acceptor Hydrogen Bond (SAHB) definition for anions only.
  • WD-SND.f: Weakest Donor (WD), Strongest Non-Donor (SND) and Weakest Donor/Strongest Non-Donor (WD/SND) definitions for cations and noble gas solutes only.


How to use the codes starting from a GROMACS simuation:

These programs will give hydration numbers for a system of one solute in any number of SPC/E water molecules with variable parameters which can be defined at the beginning of each fortran file. In fact, any 3-site water model cal be used instead of SPC/E. 4 site water models can also be used by defining the C.O.Mass of the Lennard-Jones site as the coordinate of water's oxygen atom, but the codes have only been explicitly tested for the SPC/E water model. The following steps should be followed to get hydration numbers:

  1. For each solute-water system, run a GROMACS simulation in the NVT ensemble with cubic Periodic Boundary Conditions and generate a trajectory file (.trr) with a specified number of total frames.

  2. Using the gromacs tool trjconv (refer latest GROMACS manual for details), and the full trajectory file (traj.trr) as input, print 4 formatted trajectory files (in the .gro format) for the 4 different subsystems of the box: 'ionout.gro' for the positions of only the solute molecule, 'hw1out.gro' for one hydrogen atom of each water molecule, 'hw2out.gro' for the second hydrogen atom of each water molecule, and 'owout.gro' for the oxygen atom of each water molecule. As we are printing out all different subsystems in separate .gro files, the relative order of molecules (whether solute first and waters later or waters first and solute later) in the original trajectory file
    (traj.trr) will not matter in the end.

  3. Using the 4 .gro files as input, compile the fortran code corresponding to the definition required, using the gfortran compiler, and run the executable a.out.

  4. Record/plot the output given in output files described inside each fortran file.