<br><div class="gmail_quote">2009/9/8 Mark Abraham <span dir="ltr"><<a href="mailto:Mark.Abraham@anu.edu.au">Mark.Abraham@anu.edu.au</a>></span><br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<div class="im">Jennifer Williams wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Hi users,<br>
<br>
I am running a very simple simulation of methane inside a pore (v.much like a carbon nanotube but in my case the tube is supposed to represent silica.) I keep this tube frozen.<br>
<br>
I start with an energy minimisation-however this runs to completion almost instantly and I keep get NaN for my potential energy:<br>
<br>
Steepest Descents converged to machine precision in 18 steps,<br>
but did not reach the requested Fmax < 10.<br>
Potential Energy = nan<br>
Maximum force = 6.5738518e+01 on atom 2133<br>
Norm of force = 1.5461593e+00<br>
</blockquote>
<br></div>
This nan suggests some kind of severe atomic overlap. Reconsider your coordinates and the box size implied by your coordinate file.<br>
<br>
Mark<div><div></div><div class="h5"><br>
<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Otherwise the trajectory looks OK (methane moving around inside the cylinder). If I go on to use the conf.gro file for an mdrun, it runs to completion and generates what looks like a reasonable trajectory, however the output again contains NaN i.e:<br>
<br>
Energies (kJ/mol)<br>
LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energy<br>
nan 0.00000e+00 nan 3.36749e+01 nan<br>
Conserved En. Temperature Pressure (bar)<br>
nan 3.00010e+02 nan<br>
<br>
and calculating the Diffusion coefficient gives:<br>
D[ CH4] 613.6682 (+/- 97.0563) 1e-5 cm^2/s<br>
<br>
If I do the same calculation but reduce the cut-offs to 0.9. I get<br>
<br>
Energies (kJ/mol)<br>
LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energy<br>
nan 0.00000e+00 nan 3.36750e+01 nan<br>
Conserved En. Temperature Pressure (bar)<br>
nan 3.00011e+02 nan<br>
<br>
D[ CH4] 237.8712 (+/- 53.5975) 1e-5 cm^2/s<br>
<br>
And for a cut-off of 1.3nm I get<br>
<br>
Energies (kJ/mol)<br>
LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energ<br>
y<br>
nan 0.00000e+00 nan 3.36737e+01 na<br>
n<br>
Conserved En. Temperature Pressure (bar)<br>
nan 2.99999e+02 nan<br>
<br>
<br>
D[ CH4] 19.7953 (+/- 154.0168) 1e-5 cm^2/s<br>
<br>
<br>
For this system, the cut-off shouldn?t need to be larger than 0.8 (I have plotted graphs of calculated V vs r) so it is worrying that the diffusion coefficient is showing such dependence on the cut-offs when they should all give the same result.<br>
<br>
Can anyone offer any insight into this? I?ve tried changing the timestep making it both larger and smaller and many other things. I?ve pasted the relevant parts of my files below:<br>
<br>
I?m using gromacs 4.0.5 ?at the moment running in serial.<br>
<br>
Thanks for any advice,<br>
<br>
Top file<br>
<br>
[ defaults ]<br>
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ<br>
1 2 yes 1.0 1.0<br>
;<br>
;<br>
[ atomtypes ]<br>
; type mass charge ptype c6 c12<br>
OSM 15.9994 0.00 A 0.2708 1.538176<br>
<br>
;<br>
; Include forcefield parameters<br>
#include "CH4.itp"<br>
;<br>
;<br>
[ moleculetype ]<br>
; Name nrexcl<br>
MCM 3<br>
[ atoms ]<br>
; nr type resnr residue atom cgnr charge mass<br>
1 OSM 1 MCM OSM 1 0 15.9994<br>
2 OSM 1 MCM OSM 2 0 15.9994<br>
.etc<br>
2127 OSM 1 MCM OSM 2127 0 15.9994<br>
2128 OSM 1 MCM OSM 2128 0 15.9994<br>
<br>
<br>
[ system ]<br>
; Name<br>
CH4 in MCM<br>
<br>
[ molecules ]<br>
; Compound #mols<br>
MCM 1<br>
CH4 10<br>
<br>
CH4.itp file<br>
<br>
[ atomtypes ]<br>
; type mass charge ptype c6 c12<br>
CH4 16.043 0.00 A 0.3732 1.24650457<br>
;<br>
[ moleculetype ]<br>
; name nrexcl<br>
CH4 2<br>
<br>
[ atoms ]<br>
; nr type resnr residu atom cgnr charge mass<br>
1 CH4 1 CH4 CH4 1 0.00 16.043<br>
<br>
<br>
<br>
.mdp file<br>
<br>
;<br>
; File 'mdout.mdp' was generated<br>
; By user: jwillia4 (353773)<br>
; On host: vlxhead2<br>
; At date: Fri Jun 26 15:47:37 2009<br>
;<br>
; VARIOUS PREPROCESSING OPTIONS<br>
; Preprocessor information: use cpp syntax.<br>
; e.g.: -I/home/joe/doe -I/home/mary/hoe<br>
include = -I../top<br>
; e.g.: -DI_Want_Cookies -DMe_Too<br>
define =<br>
<br>
; RUN CONTROL PARAMETERS<br>
integrator = steep<br>
; Start time and timestep in ps<br>
tinit = 0<br>
dt = 0.0001<br></blockquote></div></div></blockquote><div><br>dt = 0.1 fs are you sure??<br> <br></div><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<div><div class="h5"><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
nsteps = 100000<br>
; For exact run continuation or redoing part of a run<br>
; Part index is updated automatically on checkpointing (keeps files separate)<br>
simulation_part = 1<br>
init_step = 0<br>
; mode for center of mass motion removal<br>
comm-mode = linear<br>
; number of steps for center of mass motion removal<br>
nstcomm = 1<br>
; group(s) for center of mass motion removal<br>
comm-grps =<br>
<br>
; LANGEVIN DYNAMICS OPTIONS<br>
; Friction coefficient (amu/ps) and random seed<br>
bd-fric = 0<br>
ld-seed = 1993<br>
<br>
; ENERGY MINIMIZATION OPTIONS<br>
; Force tolerance and initial step-size<br>
emtol =<br>
emstep = 0.001<br>
; Max number of iterations in relax_shells<br>
niter =<br>
; Step size (ps^2) for minimization of flexible constraints<br>
fcstep =<br>
; Frequency of steepest descents steps when doing CG<br>
nstcgsteep =<br>
nbfgscorr =<br>
<br>
<br>
; OUTPUT CONTROL OPTIONS<br>
; Output frequency for coords (x), velocities (v) and forces (f)<br>
nstxout = 100<br>
nstvout = 100<br>
nstfout = 0<br>
; Output frequency for energies to log file and energy file<br>
nstlog = 100<br>
nstenergy = 100<br>
; Output frequency and precision for xtc file<br>
nstxtcout = 100<br>
xtc-precision = 100<br>
; This selects the subset of atoms for the xtc file. You can<br>
; select multiple groups. By default all atoms will be written.<br>
xtc-grps =<br>
; Selection of energy groups<br>
energygrps =<br>
<br>
; NEIGHBORSEARCHING PARAMETERS<br>
; nblist update frequency<br>
nstlist =<br>
; ns algorithm (simple or grid)<br>
ns_type = grid<br>
; Periodic boundary conditions: xyz, no, xy<br>
pbc = xyz<br>
periodic_molecules = yes<br>
; nblist cut-off<br>
rlist = 1.7<br>
<br>
; OPTIONS FOR ELECTROSTATICS AND VDW<br>
; Method for doing electrostatics<br>
coulombtype = Cut-off<br>
rcoulomb-switch = 0<br>
rcoulomb = 1.7<br>
; Relative dielectric constant for the medium and the reaction field<br>
epsilon_r =<br>
epsilon_rf =<br>
<br>
; Method for doing Van der Waals<br>
vdw-type = Cut-off<br>
; cut-off lengths<br>
rvdw-switch = 0<br>
rvdw = 1.7<br>
; Apply long range dispersion corrections for Energy and Pressure<br>
DispCorr = No<br>
; Extension of the potential lookup tables beyond the cut-off<br>
table-extension =<br>
; Seperate tables between energy group pairs<br>
energygrp_table =<br>
<br>
<br>
; Spacing for the PME/PPPM FFT grid<br>
fourierspacing = 0.12<br>
; FFT grid size, when a value is 0 fourierspacing will be used<br>
fourier_nx = 0<br>
fourier_ny = 0<br>
fourier_nz = 0<br>
; EWALD/PME/PPPM parameters<br>
pme_order =<br>
ewald_rtol = 1e-05<br>
ewald_geometry = 3d<br>
epsilon_surface = 0<br>
optimize_fft = yes<br>
<br>
; IMPLICIT SOLVENT ALGORITHM<br>
implicit_solvent = No<br>
<br>
<br>
; OPTIONS FOR WEAK COUPLING ALGORITHMS<br>
; Temperature coupling<br>
tcoupl = no<br>
; Groups to couple separately<br>
tc-grps =<br>
; Time constant (ps) and reference temperature (K)<br>
tau_t =<br>
ref_t =<br>
<br>
; Pressure coupling<br>
Pcoupl = No<br>
Pcoupltype =<br>
; Time constant (ps), compressibility (1/bar) and reference P (bar)<br>
tau-p =<br>
compressibility =<br>
ref-p =<br>
; Scaling of reference coordinates, No, All or COM<br>
refcoord_scaling = no<br>
; Random seed for Andersen thermostat<br>
andersen_seed =<br>
<br>
<br>
; GENERATE VELOCITIES FOR STARTUP RUN<br>
gen_vel = no<br>
gen_temp = 300<br>
gen_seed = 173529<br>
<br>
; OPTIONS FOR BONDS<br>
constraints = none<br>
; Type of constraint algorithm<br>
constraint-algorithm = Lincs<br>
; Do not constrain the start configuration<br>
continuation = no<br>
; Use successive overrelaxation to reduce the number of shake iterations<br>
Shake-SOR = no<br>
; Relative tolerance of shake<br>
shake-tol = 0.0001<br>
; Highest order in the expansion of the constraint coupling matrix<br>
lincs-order = 4<br>
; Number of iterations in the final step of LINCS. 1 is fine for<br>
; normal simulations, but use 2 to conserve energy in NVE runs.<br>
; For energy minimization with constraints it should be 4 to 8.<br>
lincs-iter = 1<br>
; Lincs will write a warning to the stderr if in one step a bond<br>
; rotates over more degrees than<br>
lincs-warnangle = 30<br>
; Convert harmonic bonds to morse potentials<br>
morse = no<br>
<br>
; ENERGY GROUP EXCLUSIONS<br>
; Pairs of energy groups for which all non-bonded interactions are excluded<br>
energygrp_excl =<br>
<br>
<br>
; Non-equilibrium MD stuff<br>
acc-grps =<br>
accelerate =<br>
freezegrps = MCM<br>
freezedim = Y Y Y<br>
cos-acceleration = 0<br>
deform =<br>
<br>
<br>
<br>
<br>
<br>
<br>
</blockquote></div></div><div class="im">
_______________________________________________<br>
gmx-users mailing list <a href="mailto:gmx-users@gromacs.org" target="_blank">gmx-users@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-users" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-users</a><br>
Please search the archive at <a href="http://www.gromacs.org/search" target="_blank">http://www.gromacs.org/search</a> before posting!<br></div>
Please don't post (un)subscribe requests to the list. Use the www interface or send it to <a href="mailto:gmx-users-request@gromacs.org" target="_blank">gmx-users-request@gromacs.org</a>.<div><div></div><div class="h5">
<br>
Can't post? Read <a href="http://www.gromacs.org/mailing_lists/users.php" target="_blank">http://www.gromacs.org/mailing_lists/users.php</a><br></div></div></blockquote><div> </div></div>Manik Mayur<br>Graduate student<br>
Microfluidics Lab<br>Dept. of Mechanical Engg.<br>IIT Kharagpur<br>INDIA<br>