analysis.rdf Not Normalizing to 1 at Large Distances

[ahrazz99]

Hi,

I am trying to calculate the Radial Distribution Function (RDF) of salt ions to water molecules in a 2 phase system of ice and water. I allowed the system to progress such that some ions have been trapped in the ice and the remainder are in the solution. I have already located the phase boundary and use that to select the ions on in each phase. I also use it to select the oxygens of each water molecule in each phase. I then input these into necessary functions, run, and plot the results. However, the plots have been showing a normalization at 2 or higher rather than at 1 and I am not sure what to change about it. Is there something I am doing wrong/missing that I need to fix to get the correct normalization?

My code is shown below:

#Get the atom groups
#   Need groups for oxygens left of the ice, oxygens right of the ice, each ion left of the ice and the ions right
#       of the ice
interface = str(interfaceLocations[count])
selPhrase = "prop x < " + interface + " and type 1"
oxyInIce = system.select_atoms(selPhrase)
selPhrase = "prop x < " + interface + " and type 3"
NaInIce = system.select_atoms(selPhrase)
selPhrase = "prop x < " + interface + " and type 4"
ClInIce = system.select_atoms(selPhrase)
selPhrase = "prop x > " + interface + " and type 1"
oxyInWater = system.select_atoms(selPhrase)
selPhrase = "prop x > " + interface + " and type 3"
NaInWater = system.select_atoms(selPhrase)
selPhrase = "prop x > " + interface + " and type 4"
ClInWater = system.select_atoms(selPhrase)


#Set up list of atom groups for RDF calculation and other parameters
if len(NaInIce) > 0:
    rdfNaIce = rdf.InterRDF(NaInIce, oxyInIce, norm="rdf")
    rdfNaIce.run()
    rdf_results_naice.append(rdfNaIce.results.rdf)
if len(ClInIce) > 0:
    rdfClIce = rdf.InterRDF(ClInIce, oxyInIce, norm="rdf")
    rdfClIce.run()
    rdf_results_clice.append(rdfClIce.results.rdf)

#Calculate the rdf for each atom pair
rdfNaWater = rdf.InterRDF(NaInWater, oxyInWater, norm="rdf")
rdfClWater = rdf.InterRDF(ClInWater, oxyInWater, norm="rdf")

rdfNaWater.run()
rdfClWater.run()

#Add results to the set of results
bins = rdfNaWater.results.bins
#print(bins)
rdf_results_nawater.append(rdfNaWater.results.rdf)
rdf_results_clwater.append(rdfClWater.results.rdf)

#Graph the individual calculation results
fig, (ax1, ax2) = plt.subplots(2, sharex=True)
fig.set_figwidth(10)
fig.set_figheight(15)
ax1.plot(bins, rdfNaWater.results.rdf, "r-", label="Na in Solution")
ax1.plot(bins, rdfNaIce.results.rdf, "b--", label="Na in Ice")
ax2.plot(bins, rdfClWater.results.rdf, "m-", label="Cl in Solution")
ax2.plot(bins, rdfClIce.results.rdf, "g--", label="Cl in Ice")
label = "r (" + r'$\AA$' + ")"
ax2.set_xlabel(label)
ax1.set_xlabel(label)
ax1.set_ylabel("g(r)")
ax2.set_ylabel("g(r)")
ax1.set_title("Sodium Ion RDF")
ax2.set_title("Chloride Ion RDF")
ax1.set_xlim(0, 15.0)
ax2.set_xlim(0, 15.0)
ax1.legend(loc="upper right")
ax2.legend(loc="upper right")
figName = "output/testing_rdf_2_" + str(count) + ".png"
fig.savefig(figName)

The results are plotted in the following figure.

Thanks!

[hmacdope]

@ahrazz99 would you be able to examine a longer RDF cutoff (if needed in a larger system) for the ice system? Ice has a periodic structure with a lot more coherence and lower entropy than water and may have long range structure at > 14 A. Perhaps some snapshots of your system may also help.

[ahrazz99]

@hmacdope, thank you for the assistance.
I increased the cutoff to 50 A and it seems to drop from ~2 to fluctutuating around 0. I also increased the bin count to 500. I've included another image below. I am attempting to recreate the findings of Luo et al. in their article found here (see Fig. 2): https://pubs.rsc.org/en/content/articlelanding/2021/CP/D1CP01733K

With increased cutoff and bin count:

As for my system, I have included an image below. It has 20270 particles, 6720 H2O molecules, and 55 NaCl ion pairs. I allowed it to evolve in a NPT ensemble with pressure fluctuations allowed in the x direction. The models I used were the TIP4P/2005 for the water molecules and the Madrid-2019 Ion potentials. The system snapshot is below.

Thanks again for the help.

[hmacdope]

Hi, sorry for the slow response.

So it looks like the reason it is decaying to 0 is that there is empty space in your box! When the averaging is done for the RDF, even if the empty box is eventually filled, the time average will still see the empty box, so it will drop below one. An RDF that decays neatly to one is a function of the periodicity of the system.

[ahrazz99]

! I removed the gaps and recalculated the RDFs and they are now normalizing. Thank you for the help!