[gmx-users] Re: Lennard-Jones Rdf plot

Lum Nforbi lumngwegia at gmail.com
Tue Feb 9 00:23:17 CET 2010


  This is in response to your mail on this subject. My objective was to do
an NPT simulation of a 200 particle fluid using Lennard-Jones parameters
only. I also used the atomic number, atomic mass, and Lennard-Jones
parameters(sigma in nm and epsilin in KJ/mol) for the oxygen atom, and the
force field files which I generated. I have included my procedure and  the
parameter and topology files I used below.

(1) I chose one set of 3 random coordinates(x,y,z) and multiplied them using
packmol to generate 200x3 random coordinates.
data1.pdb file is below:
ATOM      1  O   OXY     1       5.120  10.240  25.600  1.00  0.00
Command used to generate 200 coordinates: packmol < data1.inp. The cubic box
size is 27A which was later converted to 2.7nm in a .gro file as shown in
the command :   editconf -f data1_packmol.pdb -o data1_packmol.gro

(2) The energy of the system was minimized using the following 2 commands:
grompp and mdrun

 (a) grompp -f oxy_min.mdp -c data200_packmol.gro -p oxy.top -o oxymin.tpr
  Result from grompp:
Generated 0 of the 1 non-bonded parameter combinations
Excluding 3 bonded neighbours molecule type 'OXY'
processing coordinates...
double-checking input for internal consistency...
renumbering atomtypes...
converting bonded parameters...
initialising group options...
processing index file...
Analysing residue names:
Opening library file /usr/local/gromacs/share/
There are:   200      OTHER residues
There are:     0    PROTEIN residues
There are:     0        DNA residues

(b) mdrun -nice 0 -v -s oxymin.tpr -o oxymin.trr -c oxymin.gro -e
    Result of mdrun shows system was well minimized:
 Step=   24, Dmax= 1.2e-01 nm, Epot= -5.46417e+02 Fmax= 1.16183e+03, atom=
 Step=   26, Dmax= 6.9e-02 nm, Epot= -5.72183e+02 Fmax= 2.29392e+02, atom=
 Step=   27, Dmax= 8.3e-02 nm, Epot= -5.82409e+02 Fmax= 1.04963e+02, atom=
 Step=   30, Dmax= 2.5e-02 nm, Epot= -5.87203e+02 Fmax= 8.81632e+01, atom=

writing lowest energy coordinates.

Steepest Descents converged to Fmax < 100 in 31 steps
Potential Energy  = -5.8720325e+02
Maximum force     =  8.8163155e+01 on atom 93
Norm of force     =  1.4542491e+01

(3) My topology and force field files  are below. As you can see, I used
only nonbonding parameters.

; The force field files to be included
#include "ffgmxO.itp"
[ moleculetype ]
; molname        nrexcl
OXY            3
[ atoms ]
; id    at type res nr residu name      at name     cg nr charge
1     O      1     OXY            O          1  0.000
[ system ]
Pure Oxygen
[ molecules ]
;molecule name number
OXY 200

#define _FF_GROMACS
#define _FF_GROMACS1
[ defaults ]
; nbfunc        comb-rule       gen-pairs  FudgeLJ
  1             2               no
#include "ffgmxO_nb.itp"

[ atomtypes ]
;name at.num mass                 charge ptype sigma     epsilon
    O 8         15.999400    0.000     A 0.3433   0.939482
[ nonbond_params ]
  ; i j func     sigma    epsilon
    O O 1 0.3433          0.939482
[ pairtypes ]
  ;i     j func     sigma      epsilon
    O O 1 0.3433          0.939482

(4) My .mdp file for mdrun is below:

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                   = 5000000    ; total number of steps; total time
(10 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

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

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.
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
ref_p                    = 1.0        ; ref. pressure for coupling (bar)

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

(5) Results at the end of the simulation

Statistics over 5000001 steps [ 0.0000 thru 10000.0000 ps ], 7 data sets
All averages are exact over 5000001 steps

Energy                      Average       RMSD     Fluct.      Drift
Potential                  -1044.97    22.9468    22.7766 -0.000966612
Kinetic En.                 198.354    11.4282    11.4279 -3.0823e-05
Total Energy                -846.62    25.8194    25.6583 -0.000997435
Temperature                 79.9208    4.60467    4.60453 -1.24192e-05
Pressure (bar)             0.278412     89.842    89.8327 0.000447827
Volume                      10.0719    0.42992   0.426807 -1.78902e-05
Density (SI)                528.134    13.6698    13.5568 0.000607606
Heat Capacity Cv:      12.5342 J/mol K (factor = 0.00331954)
Isothermal Compressibility:  0.0016631 /bar
Adiabatic bulk modulus:        601.288  bar

(6) I then plotted a Potential Energy graph and a radial distribution
function for the system at the end of the simulation. I have attached the
graphs. The commands for g_rdf are below:

make_ndx -f oxymdrun.gro -o oxy.ndx
trjconv -f oxymdrun_traj.trr -o oxymdrun_traj.xtc
g_rdf -f oxymdrun_traj.xtc -n oxy.ndx -o oxy_rdf.xvg

Select a reference group and 1 group
Group     0 (      System) has   200 elements
Group     1 (         OXY) has   200 elements
Select a group: 1
Selected 1: 'OXY'
Select a group: 1
Selected 1: 'OXY'
Last frame      10000 time 10000.000

The first two peaks on the radial distribution function plot are ok, but the
third one starts slanting down ie it does not converge to 1.0. I am not sure
what I am not setting right and I don't know what could be the cause of the
graph slanting downwards.
In section 2.3 of page 9 of the manual, it is mentioned that reduced units
are good for simulating LJ systems, but Gromacs already has a fixed set
units that it uses and I am not sure how to set the reduced units such that
Gromacs can recognize them. Could it be a reduced units problem?

I appreciate your answers.

- Show quoted text -

Lum Nforbi wrote:
> Dear all,
>    Please, could someone tell me what could be the problem with the
> attached LJ plot?
> Could it be due to wrong assignment of periodic boundary conditions?

First, do not attach files without file extensions.  For those of us who are
security-conscious (some say "paranoid"), we won't open it.  It is much more
useful to post an image file to some public site, like Photobucket or Google

Second, there is likely no way anyone can diagnose your problem unless you
a very detailed description of what you're doing.  You mention an "LJ plot"
the file name of your attachment would suggest a radial distribution
What is it that you are trying to analyze?  How did you produce the plot?
 Is it
indeed an RDF?  If so, what was simulated, and for how long?  What were your
.mdp parameters?  It is very doubtful that PBC alone caused issues, but
your incorrect assignment of, e.g., cutoffs or some other feature.  But we
know any of this yet.

Remember, we can't get inside your head to know what you're doing.  You've
to make it easy for us.


> Thank you,
> Lum
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100208/8349b6b7/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: RDF_oxyoxy.agr
Type: application/x-grace
Size: 20170 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100208/8349b6b7/attachment.bin>

More information about the gromacs.org_gmx-users mailing list