RDF of NaCl in water

[horlust]

Hi everyone,

I have a simulation containing 884 water molecules, 1 Na+ atom and 1 Cl- atom. I am trying to compute the RDF of Na+ and Cl- in water and also the coordination number.

At first I did:

from MDAnalysis.analysis.rdf import InterRDF

# select atoms
na_selection = u.select_atoms("name NA")
water_selection = u.select_atoms("resname SOL")

# define the range for RDF calculation
r_min = 0.0  # minimum distance (in Angstroms)
r_max = 10.0  # maximum distance (in Angstroms)
nbins = 100  # number of bins for RDF

# calculate RDF
rdf_nacl = InterRDF(na_selection, 
                        water_selection, 
                        nbins=nbins, 
                        range=(r_min, r_max), 
                        periodic=True  # account for periodic boundary conditions
                        )
rdf_nacl.run()

And using a function I calculated the coordination number of Na+ in this case, which was close to 6. However, if you look at my code you will see that I've calculated RDF with respect to the water molecules (water_selection is "resname SOL"). Is the RDF function on MDAnalysis capable of performing this calculation for a whole molecule, instead of just for atoms? I am asking because pretty much all RDF calculations that I've seen consider the oxygen atoms of the water molecule instead of the whole water molecule. So I also did this:

# select atoms
na_selection = u.select_atoms("name NA")
water_selection = u.select_atoms("resname SOL and name OW")

# define the range for RDF calculation
r_min = 0.0  # minimum distance (in Angstroms)
r_max = 10.0  # maximum distance (in Angstroms)
nbins = 100  # number of bins for RDF

# calculate RDF
rdf_nacl = InterRDF(na_selection, 
                        water_selection, 
                        nbins=nbins, 
                        range=(r_min, r_max), 
                        periodic=True  # account for periodic boundary conditions
                        )
rdf_nacl.run()

In this case, I consider only the oxygen atoms of the water molecules around Na+.

As I'm having two very different results (images attached), I'd like to understand the handling of the RDF function. Which way of calculating the RDF and the coordination number makes more sense in the context of MDAnalysis?

Thank you in advance!
Best,
Maria

[hmacdope]

Hi!

Could you clarify which graph is which experiment?

I am assuming that the 5.2 coordination number is the one using water oxygens. From my knowledge use of the oxygens as the reference is standard and gives a number that makes a lot of sense, octahedral with ~ 6 * 1 coordinating atoms (likely has 0.8 fractional Cl- coordination).

Using all the solvent molecule as a whole gives a number that also makes sense (~6 * 3 atoms per molecule), but is not the standard way to talk about water coordination, as the oxygen is normally the coordinating atom.

You can see from the data tables in a paper parametrising ions such as this one: https://pubs.acs.org/doi/10.1021/jp505875v that coordination numbers are normally reported as directly chelating atoms only.

If you are looking at data of this kind, I would reccomend our solvation subpackage solvation-analysis: https://github.com/MDAnalysis/solvation-analysis

Hope this helps.

Hugo

[horlust]

Hi Hugo!

Thank you for your answer :-)

In fact, for coordination number = 5.2 I used the whole molecule of water. For coordination number = 17 I used the water oxygens. I also used Solvation Analysis and I got this for the same NaCl simulation:

I think the Solvation Analysis output is closer to the calculation considering the whole water molecule. That's why I'm a bit confused. Maybe my coordination number calculation from RDF data is wrong too? From what I've read it should be close to 6 for Na+, but with the MDAnalysis RDF calculation I get 17. I did something like:

`# numerical density calculation
def calc_density(u):
    pbc = u.atoms.bbox(pbc=True)
    x_l = (pbc[1][0] - pbc[0][0])
    y_l = (pbc[1][1] - pbc[0][1])
    z_l = (pbc[1][2] - pbc[0][2])
    volume = x_l*y_l*z_l
    n_part = len(u.atoms)
    density = n_part/volume
    return density

# coordination number calculation
def calc_coord_number(r_polyol, g_r_polyol, r_min, r_max, rho):
    # interpolate g(r) function
    g_interpolated = np.interp

    # define the function to integrate
    def integrand(r):
        return g_interpolated(r, r_polyol, g_r_polyol) * r**2

    # perform numerical integration
    result, _ = integrate.quad(integrand, r_min, r_max)

    # calculate CN(r)
    CN_r = 4 * np.pi * rho * result

    return CN_r

  # find the first local minimum in g(r)
  minima_index = np.where(np.r_[True, g_r_nacl[1:] < g_r_nacl[:-1]] & np.r_[g_r_nacl[:-1] < g_r_nacl[1:], True])[0][0]
  
  # define the integration limits
  r_min_integrate = 0
  r_max_integrate = r_nacl[minima_index]
  
  CN_r = calc_coord_number(r_nacl, g_r_nacl, r_min_integrate, r_max_integrate, calc_density(u))`

Thanks again, your answer have been very helpful.
Maria

[rodpollet]

Hi Maria,

I checked your code and found two mistakes.

First, I think you should not use bbox but u.trajectory.ts.dimensions, namely:
x_l = u.trajectory.ts.dimensions[0]
y_l = u.trajectory.ts.dimensions[1]
z_l = u.trajectory.ts.dimensions[2]
This should give a better estimate of the volume.

Second, I think you should not set n_part to the total number of atoms in your universe but only the number of atoms in the second species, namey water oxygen atoms only. Again, this will correct the density.

With these corrections, I think you should obtain the correct result. Do you confirm?

Hope this helps,
Rodolphe

[horlust]

Hi Rodolphe,

I modified the code according to your corrections and I got a coordination number of 5.65 for Na+ and 6.84 for Cl-, for the ion regarding the oxygen atoms. I think it is correct now. Thanks!

[hmacdope]

Amazing! Thanks @rodpollet!

[orbeckst]

Now that we know that it's the volume, I'd just use the built-in unit cell volume calculation, which will be correct for any box shape:

volume = u.trajectory.ts.volume

[horlust]

@rodpollet, I kept thinking about your answer and I don't understand why n_part should be the number of oxygen atoms and not the number of atoms "tout court". Can you explain? :)

[orbeckst]

For RDF calculations if ions in water I would choose ion-oxygen. The MDA InterRDF() class will compute RDFs for atom-atom distances; therefore in your first example you get a superposition of the Na-O and the NA-H RDFs. There is not option for atom-molecule or molecule-molecule evaluations at the moment.

[hmacdope]

@horlust is your question answered? If so can you mark one of the answers as such.