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 >
RDF_oxyoxy.agr
Description: application/grace
-- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php