Justin,<br><br> 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.<br>
<br>(1) I chose one set of 3 random coordinates(x,y,z) and multiplied them using packmol to generate 200x3 random coordinates. <br>data1.pdb file is below:<br>ATOM 1 O OXY 1 5.120 10.240 25.600 1.00 0.00<br>
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<br>
<br>(2) The energy of the system was minimized using the following 2 commands: grompp and mdrun<br><br> (a) grompp -f oxy_min.mdp -c data200_packmol.gro -p oxy.top -o oxymin.tpr<br> Result from grompp:<br>Generated 0 of the 1 non-bonded parameter combinations<br>
Excluding 3 bonded neighbours molecule type 'OXY'<br>processing coordinates...<br>double-checking input for internal consistency...<br>renumbering atomtypes...<br>converting bonded parameters...<br>initialising group options...<br>
processing index file...<br>Analysing residue names:<br>Opening library file /usr/local/gromacs/share/<div id=":yl" class="ii gt">gromacs/top/aminoacids.dat<br>There are: 200 OTHER residues<br>There are: 0 PROTEIN residues<br>
There are: 0 DNA residues<br>
<br>(b) mdrun -nice 0 -v -s oxymin.tpr -o oxymin.trr -c oxymin.gro -e oxymin_ener.edr<br> Result of mdrun shows system was well minimized:<br> Step= 24, Dmax= 1.2e-01 nm, Epot= -5.46417e+02 Fmax= 1.16183e+03, atom= 39<br>
Step= 26, Dmax= 6.9e-02 nm, Epot= -5.72183e+02 Fmax= 2.29392e+02, atom= 157<br> Step= 27, Dmax= 8.3e-02 nm, Epot= -5.82409e+02 Fmax= 1.04963e+02, atom= 65<br> Step= 30, Dmax= 2.5e-02 nm, Epot= -5.87203e+02 Fmax= 8.81632e+01, atom= 93<br>
<br>writing lowest energy coordinates.<br><br>Steepest Descents converged to Fmax < 100 in 31 steps<br>Potential Energy = -5.8720325e+02<br>Maximum force = 8.8163155e+01 on atom 93<br>Norm of force = 1.4542491e+01<br>
<br><br>(3) My topology and force field files are below. As you can see, I used only nonbonding parameters.<br><br>OXY.TOP<br>; The force field files to be included<br>#include "ffgmxO.itp"<br>[ moleculetype ]<br>
; molname nrexcl<br>OXY 3<br>[ atoms ]<br>; id at type res nr residu name at name cg nr charge<br>1 O 1 OXY O 1 0.000<br>[ system ]<br>Pure Oxygen<br>[ molecules ]<br>
;molecule name number<br>OXY 200<br><br>FFGMXO.ITP<br>#define _FF_GROMACS<br>#define _FF_GROMACS1<br>[ defaults ]<br>; nbfunc comb-rule gen-pairs FudgeLJ <br> 1 2 no<br>#include "ffgmxO_nb.itp"<br>
<br>FFGMXO_NB.ITP<br>[ atomtypes ]<br>;name at.num mass charge ptype sigma epsilon<br> O 8 15.999400 0.000 A 0.3433 0.939482<br>[ nonbond_params ]<br> ; i j func sigma epsilon<br>
O O 1 0.3433 0.939482<br>[ pairtypes ]<br> ;i j func sigma epsilon<br> O O 1 0.3433 0.939482<br><br><br>(4) My .mdp file for mdrun is below:<br><br>title = NPT simulation of a LJ FLUID<br>
cpp = /lib/cpp<br>include = -I../top<br>define = <br>integrator = md ; a leap-frog algorithm for integrating Newton's equations of motion<br>
dt = 0.002 ; time-step in ps<br>nsteps = 5000000 ; total number of steps; total time (10 ns)<br><br>nstcomm = 1 ; frequency for com removal<br>nstxout = 500 ; freq. x_out<br>
nstvout = 500 ; freq. v_out<br>nstfout = 0 ; freq. f_out<br>nstlog = 50 ; energies to log file<br>nstenergy = 50 ; energies to energy file<br>
<br>nstlist = 10 ; frequency to update neighbour list<br>ns_type = grid ; neighbour searching type<br>rlist = 1.0 ; cut-off distance for the short range neighbour list<br>
<br>coulombtype = PME ; particle-mesh-ewald electrostatics<br>rcoulomb = 1.0 ; distance for the coulomb cut-off<br>vdw-type = Cut-off ; van der Waals interactions<br>
rvdw = 1.0 ; distance for the LJ or Buckingham cut-off<br><br>fourierspacing = 0.12 ; 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 = 300.0 ; temperature for Maxwell distribution<br>gen_seed = 173529 ; used to initialize random generator for random velocities<br>
<br><br>(5) Results at the end of the simulation<br><br>Statistics over 5000001 steps [ 0.0000 thru 10000.0000 ps ], 7 data sets<br>All averages are exact over 5000001 steps<br><br>Energy Average RMSD Fluct. Drift Tot-Drift<br>
-------------------------------------------------------------------------------<br>Potential -1044.97 22.9468 22.7766 -0.000966612 -9.66612<br>Kinetic En. 198.354 11.4282 11.4279 -3.0823e-05 -0.30823<br>
Total Energy -846.62 25.8194 25.6583 -0.000997435 -9.97435<br>Temperature 79.9208 4.60467 4.60453 -1.24192e-05 -0.124192<br>Pressure (bar) 0.278412 89.842 89.8327 0.000447827 4.47828<br>
Volume 10.0719 0.42992 0.426807 -1.78902e-05 -0.178902<br>Density (SI) 528.134 13.6698 13.5568 0.000607606 6.07606<br>Heat Capacity Cv: 12.5342 J/mol K (factor = 0.00331954)<br>
Isothermal Compressibility: 0.0016631 /bar<br>Adiabatic bulk modulus: 601.288 bar<br><br>(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:<br>
<br>make_ndx -f oxymdrun.gro -o oxy.ndx<br>trjconv -f oxymdrun_traj.trr -o oxymdrun_traj.xtc<br>g_rdf -f oxymdrun_traj.xtc -n oxy.ndx -o oxy_rdf.xvg<br><br>Select a reference group and 1 group<br>Group 0 ( System) has 200 elements<br>
Group 1 ( OXY) has 200 elements<br>Select a group: 1<br>Selected 1: 'OXY'<br>Select a group: 1<br>Selected 1: 'OXY'<br>Last frame 10000 time 10000.000 <br>
<br><br>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.<br>
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? <br>
<br>I appreciate your answers.<br><br>Lum<div><div><span id="q_126afe8c802b5f4e_1" class="h4">- Show quoted text -</span></div><div class="h5"><br><br>
<br><br>
Lum Nforbi wrote:<br>
> Dear all,<br>
><br>
> Please, could someone tell me what could be the problem with the<br>
> attached LJ plot?<br>
> Could it be due to wrong assignment of periodic boundary conditions?<br>
><br>
<br></div></div>
First, do not attach files without file extensions. For those of us who are<br>
security-conscious (some say "paranoid"), we won't open it. It is much more<br>
useful to post an image file to some public site, like Photobucket or Google Docs.<br>
<br>
Second, there is likely no way anyone can diagnose your problem unless you post<br>
a very detailed description of what you're doing. You mention an "LJ plot" but<br>
the file name of your attachment would suggest a radial distribution function.<br>
What is it that you are trying to analyze? How did you produce the plot? Is it<br>
indeed an RDF? If so, what was simulated, and for how long? What were your<br>
.mdp parameters? It is very doubtful that PBC alone caused issues, but rather<br>
your incorrect assignment of, e.g., cutoffs or some other feature. But we don't<br>
know any of this yet.<br>
<br>
Remember, we can't get inside your head to know what you're doing. You've got<br>
to make it easy for us.<br>
<br>
-Justin<br>
<br>
> Thank you,<br>
> Lum<br>
></div>