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

[Justin A. Lemkul]

Lum Nforbi wrote:
> 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.

Just because you got the desired result does not mean what you did was right. 
Haphazard combinations of parameters lead to unpredictable and unreliable 
results.  I think your rlist is absurdly short - what motivated this choice? 
Why have you yet to use the value calculated based on the Allen and Tildesley 
rule of thumb (0.9 nm)?

>     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. 

To quote from the manual (section 4.6.3): "In the neighbor list all
interaction pairs that fall within rlist are stored. Furthermore, the 
interactions between pairs that do not fall within rlist but do fall within 
max(rcoulomb,rvdw) are computed during NS, and the forces and energy are stored 
separately, and added to short-range forces at every time step between 
successive NS."

The neighbor list is updated every step for particles within rlist.  The 
neighbor list between rlist and rvdw is updated every nstlist steps; between 
these neighbor-searching steps, the neighbor list is considered to be constant. 
  Beyond rvdw, the van der Waals interaction is assumed to be zero (hence "cut 
off").  I see that you're setting "DispCorr = no," so you may be getting cutoff 
artifacts.  All long-range interactions should be attractive, so the problem 
with plain cutoffs is system expansion which is not compensated for when 
dispersion correction is off.

> 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 

Close enough does not mean good enough.  For a system of only LJ particles, your 
choice of rvdw will likely have a huge impact on your results.  You may just be 
getting some random success based on some fortuitous combination.

> 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.

I would suggest to you that you read some basic literature on force field 
derivation.  While not necessarily directly related to what you're doing, you 
could benefit from understanding how reliable parameter sets are determined.  If 
you don't comprehend what you're doing yet, you shouldn't be diving in head first.

-Justin

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


-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================