[gmx-users] Re: Lennard-Jones Rdf plot
Lum Nforbi
lumngwegia at gmail.com
Tue Feb 9 00:23:17 CET 2010
Justin,
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/
gromacs/top/aminoacids.dat
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
oxymin_ener.edr
Result of mdrun shows system was well minimized:
Step= 24, Dmax= 1.2e-01 nm, Epot= -5.46417e+02 Fmax= 1.16183e+03, atom=
39
Step= 26, Dmax= 6.9e-02 nm, Epot= -5.72183e+02 Fmax= 2.29392e+02, atom=
157
Step= 27, Dmax= 8.3e-02 nm, Epot= -5.82409e+02 Fmax= 1.04963e+02, atom=
65
Step= 30, Dmax= 2.5e-02 nm, Epot= -5.87203e+02 Fmax= 8.81632e+01, atom=
93
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.
OXY.TOP
; 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
FFGMXO.ITP
#define _FF_GROMACS
#define _FF_GROMACS1
[ defaults ]
; nbfunc comb-rule gen-pairs FudgeLJ
1 2 no
#include "ffgmxO_nb.itp"
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
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 = 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
Tot-Drift
-------------------------------------------------------------------------------
Potential -1044.97 22.9468 22.7766 -0.000966612
-9.66612
Kinetic En. 198.354 11.4282 11.4279 -3.0823e-05
-0.30823
Total Energy -846.62 25.8194 25.6583 -0.000997435
-9.97435
Temperature 79.9208 4.60467 4.60453 -1.24192e-05
-0.124192
Pressure (bar) 0.278412 89.842 89.8327 0.000447827
4.47828
Volume 10.0719 0.42992 0.426807 -1.78902e-05
-0.178902
Density (SI) 528.134 13.6698 13.5568 0.000607606
6.07606
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.
Lum
- 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
Docs.
Second, there is likely no way anyone can diagnose your problem unless you
post
a very detailed description of what you're doing. You mention an "LJ plot"
but
the file name of your attachment would suggest a radial distribution
function.
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
rather
your incorrect assignment of, e.g., cutoffs or some other feature. But we
don't
know any of this yet.
Remember, we can't get inside your head to know what you're doing. You've
got
to make it easy for us.
-Justin
> 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