Dear Mark and Amit,<br><br> 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. <br>
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. <br>
<br> The values I have used and the results are as follows:<br><br>RLIST = 0.3, RVDW = 0.7<br>RESULT: RDF Plot converges to 1 (WORKS).<br><br>RLIST = 0.3, RVDW = 1.0<br>RESULT: RDF Plot converges to 0 after third peak (DOES NOT WORK).<br>
<br>RLIST = 0.7, RVDW = 0.7<br>RESULT: RDF Plot converges to 0 after third peak (DOES NOT WORK).<br><br>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. <br>
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.<br>
<br>The .mdp file I have used below is the same in all three cases, the only modifications being the rlist and rvdw values.<br><br>integrator = md ; a leap-frog algorithm for integrating Newton's equations of motion<br>
dt = 0.002 ; time-step in ps<br>nsteps = 500000 ; total number of steps; total time (1 ns)<br><br>nstcomm = 1 ; frequency for com removal<br>nstxout = 1000 ; freq. x_out<br>
nstvout = 1000 ; freq. v_out<br>nstfout = 0 ; freq. f_out<br>nstlog = 100 ; energies to log file<br>nstenergy = 100 ; energies to energy file<br>
<br>nstlist = 10 ; frequency to update neighbour list<br>ns_type = grid ; neighbour searching type<br>rlist = 0.3 ; cut-off distance for the short range neighbour list<br>
pbc = xyz ; Periodic boundary conditions:xyz, use periodic boundary conditions in all directions<br>periodic_molecules = no ; molecules are finite, fast molecular pbc can be used<br>
<br>vdw-type = Cut-off ; van der Waals interactions<br>rvdw = 0.7 ; distance for the LJ or Buckingham cut-off<br><br>fourierspacing = 0.135 ; max. grid spacing for the FFT grid for PME<br>
fourier_nx = 0 ; highest magnitude in reciprocal space when using Ewald<br>fourier_ny = 0 ; highest magnitude in reciprocal space when using Ewald<br>fourier_nz = 0 ; highest magnitude in reciprocal space when using Ewald<br>
pme_order = 4 ; cubic interpolation order for PME<br>ewald_rtol = 1e-5 ; relative strength of the Ewald-shifted direct potential<br>optimize_fft = yes ; calculate optimal FFT plan for the grid at start up.<br>
DispCorr = no ; <br><br>Tcoupl = v-rescale ; temp. coupling with vel. rescaling with a stochastic term.<br>tau_t = 0.1 ; time constant for coupling<br>tc-grps = OXY ; groups to couple separately to temp. bath<br>
ref_t = 80 ; ref. temp. for coupling<br><br>Pcoupl = berendsen ; exponential relaxation pressure coupling (box is scaled every timestep)<br>Pcoupltype = isotropic ; box expands or contracts evenly in all directions (xyz) to maintain proper pressure<br>
tau_p = 0.5 ; time constant for coupling (ps)<br>compressibility = 4.5e-5 ; compressibility of solvent used in simulation<br>ref_p = 1.0 ; ref. pressure for coupling (bar)<br>
<br>gen_vel = yes ; generate velocities according to a Maxwell distr. at gen_temp<br>gen_temp = 80 ; temperature for Maxwell distribution<br>gen_seed = 173529 ; used to initialize random generator for random velocities<br>
<br><br>I appreciate your answers.<br><br>Lum<br><br><br><br>> Hi all,<br>
><br>
> 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?<br>
<br>
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.<br>
<br>
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.<br>
<br>
Mark<br><br>Hi Lum,<br>
<br>
These values have to be chosen very icarefully.<br>
Many artifacts can show up due to poorly chosen cutoff. I would suggest to<br>
get the parameters from a reproducible publication or something.<br>
<br>
amit<br>
<br>