[lammps-users] rdf related question

Dear all,

I consider a box of 10000 water molecules, and randomly distribute equal number of Na and Cl ions in it. The system is equilibrated at 300K and 1 atm and then I compute rdf for the different species in the box. While the rdf for Na-O, Cl-O, Na-H, Cl-H, Na-Cl and O-H seem to give results in agreement with the literature, that for Na-Na and Cl-Cl is confusing. I tested with two different cases considering 2 and 10 pairs of Na, Cl ions in the system. Both cases depict an extremely dilute solution, and so hydration of the individual ions should be observed (and that is well represented through the rdf of Na-O and Cl-O), and Na-Na, Na-Cl and Cl-Cl rdf values should remain zero typically. But for the 10 pair case, I find that each of the Na-Na and Cl-Cl rdf shows a huge (unphysical) spike at a single location and then remains zero throughout. I have tested the cases with different pressure/temperature combinations, and also by switching the long range interactions on/off. I was wondering if anyone has faced a similar issue and if there are possible suggestions of removing this spurious rdf profile.

Thank you
Ganesh

Have you looked at the dump file coords of your 2 or 10 Na and Cl ions?
Is the computed RDF spike consistent with those coords? Note that
you cannot get the RDF of atom pairs outside the cutoff distance on-the-fly
in LAMMPS. This is discussed on the compute rdf doc page.

Steve

Dear all,

I consider a box of 10000 water molecules, and randomly distribute equal
number of Na and Cl ions in it. The system is equilibrated at 300K and 1 atm
and then I compute rdf for the different species in the box. While the rdf
for Na-O, Cl-O, Na-H, Cl-H, Na-Cl and O-H seem to give results in agreement
with the literature, that for Na-Na and Cl-Cl is confusing. I tested with
two different cases considering 2 and 10 pairs of Na, Cl ions in the system.
Both cases depict an extremely dilute solution, and so hydration of the
individual ions should be observed (and that is well represented through the
rdf of Na-O and Cl-O), and Na-Na, Na-Cl and Cl-Cl rdf values should remain
zero typically. But for the 10 pair case, I find that each of the Na-Na and
Cl-Cl rdf shows a huge (unphysical) spike at a single location and then
remains zero throughout. I have tested the cases with different
pressure/temperature combinations, and also by switching the long range
interactions on/off. I was wondering if anyone has faced a similar issue and
if there are possible suggestions of removing this spurious rdf profile.

i have two questions/comments.

- how long did/do you equilibrate before you start accumulating the g(r)?
  if you place ions randomly, they may initially be _very_ close.
  if you have your trajectory, you can compute the rdfs for parts of the
  trajectory efficiently (and up to arbitrary r) with VMD.

- if you look at very dilute systems, the peak height at low r _will_
be very high.
  in this case it is always advisable to also calculate the running volume
  integral to get an impression as to how many average atoms this peak
  actually represents. i have students in my computational chemistry classes
  regularly compute all possible g(r) for salt solutions the unexpected peak
  heights for Na-Cl but also Na-Na and Cl-Cl g(r)s is often an item that
  requires discussion and clarifications.

cheers,
    axel.

The initial configuration is defined at the appropriate density and I equilibrate the system using NVT at 300K for 100 ps and then start collecting the results for the next 50 ps (with the NVT intact). To be sure that the g® results are not caused by initial closeness of the ions in the system, i also computed g® from the very start. I find that initially g® is zero for Na-Na and Cl-Cl when using 2 or 10 pairs; but with time, peaks appear for the 10 pair case, and it seems that the peak is progressing from a larger value of r towards r=0 (which indicated that though the like ions were apart when the simulation was started, they started moving towards each other with time). However, the peak seems to stop progressing with time after reaching a fixed r value (close to the sigma of the ion). i tried using VMD to plot the g® but apparently the XYZ dump from lammps needs some atom names to be set appropriately during the selection process of the vmd-rdf-plugin.

For the second point, could you suggest a suitable technique to compute the running volume integral, or should that be done by dumping the coordinates at different time intervals. Further, I checked the interaction energy values between Cl-Cl, Na-Na and Na-Cl, and find that while Na-Na and Cl-Cl have finite values, the Na-Cl is zero. This is also confusing since for dilute systems one would expect oppositely charged ions to come close and interact (if at all they are not completely hydrated), and not the like charged ions.

The initial configuration is defined at the appropriate density and I
equilibrate the system using NVT at 300K for 100 ps and then start
collecting the results for the next 50 ps (with the NVT intact). To be sure

collecting results for only 50ps seems mightily short to me.
you may get half-way decent O-O, O-H, and H-H g(r)s,
and perhaps even Na-O, Na-H, Cl-O, and Cl-H, but for Na-Na,
Cl-Cl, and Na-Cl this would result in very, very bad statistical sampling.

that the g(r) results are not caused by initial closeness of the ions in the
system, i also computed g(r) from the very start. I find that initially g(r)
is zero for Na-Na and Cl-Cl when using 2 or 10 pairs; but with time, peaks
appear for the 10 pair case, and it seems that the peak is progressing from
a larger value of r towards r=0 (which indicated that though the like ions
were apart when the simulation was started, they started moving towards each
other with time). However, the peak seems to stop progressing with time

how many different systems with 2 and 10 randomly placed ion pairs
did you generate? from what you write it almost sounds as if you created
only one of each.

after reaching a fixed r value (close to the sigma of the ion). i tried
using VMD to plot the g(r) but apparently the XYZ dump from lammps needs
some atom names to be set appropriately during the selection process of the
vmd-rdf-plugin.

using .xyz dump output is the worst possible out of all options
supported by LAMMPS as it contains the least amount of information
(and still consumes a lot of space). especially the system cell size
information is missing which makes proper use of PBC for calculation
of correct bulk system g(r)s difficult. outside of that, VMD can only
use information that you hand it. since LAMMPS doesn't us or store
any kind of atom type labels but only numbered atom types, you either
have to use that information in VMD or find other way to reconstitute
this information. the topotools plugin in VMD has a reader module
for data files, that will under certain circumstances allow to pass
atom (and residue) label information to VMD from suitably formatted
comments in the data file. there are also ways to infer atom names
via their mass or use manual VMD scripting to reconstitute the required
information and store it in a some suitable format (e.g. .psf) for later use.

For the second point, could you suggest a suitable technique to compute the
running volume integral, or should that be done by dumping the coordinates
at different time intervals. Further, I checked the interaction energy

the volume integral is best computed from the raw data of the
g(r) calculation before the normalization step. the gofr support
in VMD has this implemented.

values between Cl-Cl, Na-Na and Na-Cl, and find that while Na-Na and Cl-Cl
have finite values, the Na-Cl is zero. This is also confusing since for
dilute systems one would expect oppositely charged ions to come close and
interact (if at all they are not completely hydrated), and not the like
charged ions.

the screening of charges by the water molecules is pretty good,
and thus it is not only a question of whether unlike charged ions
attract each other, but also how strongly like-charged repel each
other. since you have good screening, there is always a good
chance that a single ion pair will come fairly close.

that all being said, there may, of course, also be an oversight
in your input and yet again, without sufficient statistical sampling
you may just be seeing something that is less likely to happen.
the fact that it is less likely to happen, doesn't mean that it never
happens. if that were the case, people would not play the lottery.

cheers,
     axel.

Dear Axel and lammps-users,

In continuation of the calculation of rdf, I would like to share some observations and seek suggestions on correcting them. I considered a system of 10 Na and 10 Cl ions distributed (with at least 10 Ang separation) in a box of 9980 water molecules, and observed the rdf. The system was equilibrated for 500ps followed by 500ps of production run, both using NVE+temp/rescale. I find that the Na-Na or Na-Cl or Cl-Cl peaks do not appear, which meant that the results obtained earlier were an artefact of the simulation geometry considered earlier. Now, looking at the mean squared displacement, I find the diffusion coefficient for oxygen is 2 orders of magnitude less than the experimentally reported values, while those for Na and Cl were of the order of 1 Ang^2. This meant that the ions were not moving freely in the system.
Also, increase the rescaling timestep increases the msd but still 2 orders of magnitude lower than the reported values. Using Langevin thermostat also does not help. All simulations were at 300K with the correct density.

I was wondering if there was a particular thermostat that needs to be applied to such systems? And suggestions about possible errors that may be impeding the free motion of the ions and molecules would really help.

Thanks
Ganesh

Thermostats, in general, will screw up diffusion. Some more than others - all depending on the parameters chosen. I don't use thermostats when I can get away with it. Equilibrate to the energy where your system will fluctuate about the temperature you want (some trial and error needed), then run NVE with a timestep that reasonably conserves energy over the trajectory length you want. That would be my approach. But you could probably get a good result with Nose-Hoover and a very weak damping. I'm guessing you haven't equilibrated long enough and/or used a too aggressive thermostat.
Matt

Thermostats, in general, will screw up diffusion. Some more than others - all depending on the parameters chosen. I don't use thermostats when I can get away with it. Equilibrate to the energy where your system will fluctuate about the temperature you want (some trial and error needed), then run NVE with a timestep that reasonably conserves energy over the trajectory length you want. That would be my approach. But you could probably get a good result with Nose-Hoover and a very weak damping. I'm guessing you haven't equilibrated long enough and/or used a too aggressive thermostat.

indeed. temp/rescale is useless in this context, particularly, if the
temperature delta is too small.
every time it rescales, it will give the system a kick and that will
mess up self-diffusion massively.

what matt suggests is the only acceptable solution. also, you have to
keep in mind that
self-diffusion has significant finite size effects, and is _very_
dependent on the chosen
water potential and thus may be far from the experimental value. to
get a well converged
SDC you have to simulate for a _very_ long time (tens to hundreds of
nanoseconds).
results in the literature for simulations of multi nanosecond time can
easily vary by a
factor of two when using the same(!) water potential.

cheers,
   axel.

If anything, a Langevin thermostat, applied to water, might
increate the MSD and diffusivity. So if you are seeing 100x too
small a diffusivity, I would guess something is fundamentally
wrong with your model. E.g. the pressure and density is too
high so the ions are not moving. Aside from MSD what
other evidence do you have that you have setup your
model correctly?

Steve

Thank for the input everyone. I am not sure what to check apart from MSD. I checked pressure and density and that corresponds to values for water 300K and 1 atm. I however like to point that no matter how long I thermostat the system, once the thermostat is released (ie simple NVE) temperature starts increasing (much beyond boiling temperature of water), though system dimensions remain the same. I am not sure how this can happen and what phase is the system at, because if it is a gas at that high temperature, it must also expand in volume and lower in pressure, but that does not happen.

Please find attached (and below) a copy of the MD code and the corresponding water file I have used. The long range interactions were neglected, and potential params were used as in the lammps manual. The waterbox was generated using the solvate-plugin in vmd. Please suggest if corrections need to be made and how can we overcome this strange behavior?

Thanks
Ganesh

in.hydration (1.27 KB)

waterfromvmd.lammps (1.87 MB)

Thank for the input everyone. I am not sure what to check apart from MSD. I checked pressure and density and that corresponds to values for water 300K and 1 atm. I however like to point that no matter how long I thermostat the system, once the thermostat is released (ie simple NVE) temperature starts increasing (much beyond boiling temperature of water), though system dimensions remain the same. I am not sure how this can happen and what phase is the system at, because if it is a gas at that high temperature, it must also expand in volume and lower in pressure, but that does not happen.

Below is a copy of the MD code. The long range interactions and SHAKE algorithm for bonds and angles were not used, and potential params were used as in the lammps manual. The waterbox was generated using the solvate-plugin in vmd. Please suggest if corrections need to be made and how can I overcome this strange behavior?

Thanks
Ganesh

Bala:

I’m not sure what you’re trying to do, but your parameters aren’t quite right for any of the variants of TIP3P.

In particular, the standard TIP3P doesn’t include a hydrogen LJ term. If you’re using the CHARMM version of TIP3P, then you need to use the specific coefficients given on the LAMMPS homepage (the O-H term does not follow either the arithmetic or the geometric mixing rules). You probably also need to use the CHARMM-modified version of the LJ potentials if you’re looking to use the TIP3P-CHARMM model.

I don’t remember off the top of my head if the CHARMM version is meant to be flexible or not; certainly the original TIP3P was meant to have constraints.

So some or all of these could be contributing to your problem.

–AEI

I accept the error on my part. However, I get the same problems when I employ the SPC model on a corresponding data file. The parameter section reads:

## SPC model for water
pair_coeff 1 1 0.00 0.00
pair_coeff 2 2 0.1553 3.166

pair_modify mix arithmetic

bond_style harmonic
bond_coeff 1 450 1.0 ## O-H bonds

angle_style harmonic
angle_coeff 1 55 109.47 ## H-O-H

special_bonds charmm
fix 11 all shake 0.0001 200 100000 b 1 a 1

I am a little confused about applying the angle and bond interactions. It seems from the Berendsen paper cited in the Lammps manual, that angle and bond length are constrained by shake algorithm, and use of harmonic bond or angle interactions are not specifically mentioned. Possibly I am making an error understanding that. It seems that to apply shake in lammps, angle and bond coeff need to be specified. Any suggestions would really help.

Thank you
Ganesh