[gmx-users] Re: Choice of value of rlist, rvdw and rcoulomb

[Lum Nforbi]

Dear Mark and Amit,

    I appreciate your replies. I have read this portion of the gromacs
manual which talks about cut-offs, but I don't clearly understand the
physical significance of the choice of values. I am using for my experiments
vdwtype = cutoff, which is twin-range cut-off  with rlist = neighbour list
cut off and rvdw = vdw cut-off where rvdw >= rlist.
     I have run three different experiments (NPT simulation of 200 particles
of oxygen at 80K ie liquid oxygen) but with different rlist and rvdw values
to evaluate the influence of these values on my Lennard-Jones plots. The
literature Lennard-Jones parameters for oxygen which I have been using are
epsilon = 0.939482 KJ/mol and sigma = 0.3433 nm. The book "Computer
Simulation of Liquids" by Allen and Tildesley recommends the use of rvdw =
2.5xsigma. In my case, if I use this recommendation, my rvdw value would be
0.86.

    The values I have used and the results are as follows:

RLIST = 0.3, RVDW = 0.7
RESULT: RDF Plot converges to 1 (WORKS).

RLIST = 0.3, RVDW = 1.0
RESULT: RDF Plot converges to 0 after third peak (DOES NOT WORK).

RLIST = 0.7, RVDW = 0.7
RESULT: RDF Plot converges to 0 after third peak (DOES NOT WORK).

I know that the problem is not about equilibration, because the first
experiment which worked did not even equilibrate fully but gave the right
plot.
    My question is: what is the physical significance of choosing these
values? I know that rlist is used as cut-off for neighbour list generation
and to specify which atom pairs are interacting, so that these interactions
can be calculated at every step of the simulation. rvdw specifies the radius
of atoms which will interact in a vdw manner. Why would a rvdw = 0.7 work
and rvdw = 1.0 not work looking at a system in a physical sense? 0.7 and 1.0
are both close to the recommended 0.86 value for rvdw, so why does one work
and not other? What could be happening in the system? If I don't understand
physically what these values really mean, I may never fully grasp what is
going on in the system. Please help me out.

The .mdp file I have used below is the same in all three cases, the only
modifications being the rlist and rvdw values.

integrator               = md         ; a leap-frog algorithm for
integrating Newton's equations of motion
dt                       = 0.002      ; time-step in ps
nsteps                   = 500000     ; total number of steps; total time (1
ns)

nstcomm                  = 1          ; frequency for com removal
nstxout                  = 1000       ; freq. x_out
nstvout                  = 1000       ; freq. v_out
nstfout                  = 0          ; freq. f_out
nstlog                   = 100         ; energies to log file
nstenergy                = 100         ; energies to energy file

nstlist                  = 10         ; frequency to update neighbour list
ns_type                  = grid       ; neighbour searching type
rlist                    = 0.3        ; cut-off distance for the short range
neighbour list
pbc                      = xyz        ; Periodic boundary conditions:xyz,
use periodic boundary conditions in all directions
periodic_molecules       = no         ; molecules are finite, fast molecular
pbc can be used

vdw-type                 = Cut-off    ; van der Waals interactions
rvdw                     = 0.7        ; distance for the LJ or Buckingham
cut-off

fourierspacing           = 0.135       ; max. grid spacing for the FFT grid
for PME
fourier_nx               = 0          ; highest magnitude in reciprocal
space when using Ewald
fourier_ny               = 0          ; highest magnitude in reciprocal
space when using Ewald
fourier_nz               = 0          ; highest magnitude in reciprocal
space when using Ewald
pme_order                = 4          ; cubic interpolation order for PME
ewald_rtol               = 1e-5       ; relative strength of the
Ewald-shifted direct potential
optimize_fft             = yes        ; calculate optimal FFT plan for the
grid at start up.
DispCorr                 = no         ;

Tcoupl                   = v-rescale  ; temp. coupling with vel. rescaling
with a stochastic term.
tau_t                    = 0.1        ; time constant for coupling
tc-grps                  = OXY        ; groups to couple separately to temp.
bath
ref_t                    = 80         ; ref. temp. for coupling

Pcoupl                   = berendsen  ; exponential relaxation pressure
coupling (box is scaled every timestep)
Pcoupltype               = isotropic  ; box expands or contracts evenly in
all directions (xyz) to maintain proper pressure
tau_p                    = 0.5        ; time constant for coupling (ps)
compressibility          = 4.5e-5     ; compressibility of solvent used in
simulation
ref_p                    = 1.0        ; ref. pressure for coupling (bar)

gen_vel                  = yes        ; generate velocities according to a
Maxwell distr. at gen_temp
gen_temp                 = 80         ; temperature for Maxwell distribution
gen_seed                 = 173529     ; used to initialize random generator
for random velocities


I appreciate your answers.

Lum



> Hi all,
>
>   Please, can someone let me know if the choice of the value of rlist,
rvdw and rcoulomb is related to or depends in someway to the distance
between atoms on a lennard-jones potential plot?

Yes, but it sounds like you should do some background reading to understand
this better. The early chapters of the GROMACS manual are an excellent
starting point.

In short, these distance cut-offs define the sizes of the sphere centred on
each atom over which various potentials get evaluated, as determined by
other .mdp options. The treatment of the boundary region depends on yet
other .mdp options. Various constraints between the rxxx values exist in
various cases to permit efficient implementation.

Mark

Hi Lum,

These values have to be chosen very icarefully.
Many artifacts can show up due to poorly chosen cutoff. I would suggest to
get the parameters from a reproducible publication or something.

amit
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100301/831fc026/attachment.html>