[gmx-users] Cut-offs in gromacs

[Mark Abraham]

On 23/02/2010 1:53 PM, Lum Nforbi wrote:
> Dear all,
>
>          I did two md simulations of 200 particles each of a
> lennard-jones fluid. One of them gave me the correct pair distribution
> function for a lennard-jones fluid (converging to 1) and one did not
> (converging to zero). I have attached the .mdp files for both systems
> below. The barostats are different but I don't think this is the cause.
> I think that one worked because of the cut-off specifications (rlist,
> rvdw and rcoulomb), but I am not sure of the explanation of how the
> cut-off values can influence the shape of a pair distribution function.
> The fourier spacing in both parameter files are also different.
>          Please, if someone knows how these cut-off values and maybe
> fourier spacing could influence the shape of a pair distribution
> function, let me know the explanation.

Shouldn't you be adopting a protocol from a reasonable-looking 
publication, or the paper of your force field, rather than haphazardly 
varying things? rvdw = 0.3 or 1.0 is a massive change.

If your LJ-fluid has no charges, then you don't want PME (though it 
won't hurt much) and its parameters will be irrelevant.

Otherwise it's pretty routine to expect a discontinuity in pair 
distributions at and above the cut-off. The model physics is changing 
there... With Coulomb interactions you get marked artefacts with 
cut-offs, not sure about LJ.

If you want to compare .mdp files, use the diff utility. If you want us 
to compare them, then provide its output!

Mark

> .mdp file which gave the plot which converges to zero:
>
> title                    = NPT simulation of a LJ FLUID
> cpp                      = /lib/cpp
> include                  = -I../top
> define                   =
> 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                  = 500        ; freq. x_out
> nstvout                  = 500        ; freq. v_out
> nstfout                  = 0          ; freq. f_out
> nstlog                   = 50         ; energies to log file
> nstenergy                = 50         ; energies to energy file
> nstlist                  = 10         ; frequency to update neighbour list
> ns_type                  = grid       ; neighbour searching type
> rlist                    = 1.0        ; 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
> coulombtype              = PME        ; particle-mesh-ewald electrostatics
> rcoulomb                 = 1.0        ; distance for the coulomb cut-off
> vdw-type                 = Cut-off    ; van der Waals interactions
> rvdw                     = 1.0        ; distance for the LJ or
> Buckingham cut-off
> fourierspacing           = 0.12       ; 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
>
> .mdp file which gave the plot which converges to 1:
>
> title                    = NPT simulation of a LJ FLUID
> cpp                      = /lib/cpp
> include                  = -I../top
> define                   =
> 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                   = 500          ; energies to log file
> nstenergy                = 500          ; 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 p b c in all directions
> periodic_molecules       = no           ; molecules are finite, fast
> molecular pbc can be used
> coulombtype              = PME          ; particle-mesh-ewald electrostatics
> rcoulomb                 = 0.3          ; distance for the coulomb cut-off
> 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                   = nose-hoover; temp. coupling with vel.
> rescaling with a stochastic term.
> tau_t                    = 0.5        ; time constant for coupling
> tc-grps                  = OXY        ; groups to couple separately to
> temp. bath
> ref_t                    = 80         ; ref. temp. for coupling
> Pcoupl                   = parrinello-rahman  ; 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                    = 5.0        ; 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 reply.
>
> Lum
>
>
>
>
> No virus found in this incoming message.
> Checked by AVG - www.avg.com
> Version: 8.5.435 / Virus Database: 271.1.1/2704 - Release Date: 02/22/10 19:34:00
>