[gmx-users] potential energy NaN and strange dependence on cut-offs
Mark Abraham
Mark.Abraham at anu.edu.au
Tue Sep 8 18:22:00 CEST 2009
Jennifer Williams wrote:
> Hi users,
>
> 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.
>
> I start with an energy minimisation-however this runs to completion
> almost instantly and I keep get NaN for my potential energy:
>
> Steepest Descents converged to machine precision in 18 steps,
> but did not reach the requested Fmax < 10.
> Potential Energy = nan
> Maximum force = 6.5738518e+01 on atom 2133
> Norm of force = 1.5461593e+00
This nan suggests some kind of severe atomic overlap. Reconsider your
coordinates and the box size implied by your coordinate file.
Mark
> 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:
>
> Energies (kJ/mol)
> LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energy
> nan 0.00000e+00 nan 3.36749e+01 nan
> Conserved En. Temperature Pressure (bar)
> nan 3.00010e+02 nan
>
> and calculating the Diffusion coefficient gives:
> D[ CH4] 613.6682 (+/- 97.0563) 1e-5 cm^2/s
>
> If I do the same calculation but reduce the cut-offs to 0.9. I get
>
> Energies (kJ/mol)
> LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energy
> nan 0.00000e+00 nan 3.36750e+01 nan
> Conserved En. Temperature Pressure (bar)
> nan 3.00011e+02 nan
>
> D[ CH4] 237.8712 (+/- 53.5975) 1e-5 cm^2/s
>
> And for a cut-off of 1.3nm I get
>
> Energies (kJ/mol)
> LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energ
> y
> nan 0.00000e+00 nan 3.36737e+01 na
> n
> Conserved En. Temperature Pressure (bar)
> nan 2.99999e+02 nan
>
>
> D[ CH4] 19.7953 (+/- 154.0168) 1e-5 cm^2/s
>
>
> 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.
>
> 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:
>
> I?m using gromacs 4.0.5 ?at the moment running in serial.
>
> Thanks for any advice,
>
> Top file
>
> [ defaults ]
> ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> 1 2 yes 1.0 1.0
> ;
> ;
> [ atomtypes ]
> ; type mass charge ptype c6 c12
> OSM 15.9994 0.00 A 0.2708 1.538176
>
> ;
> ; Include forcefield parameters
> #include "CH4.itp"
> ;
> ;
> [ moleculetype ]
> ; Name nrexcl
> MCM 3
> [ atoms ]
> ; nr type resnr residue atom cgnr charge
> mass
> 1 OSM 1 MCM OSM 1 0 15.9994
> 2 OSM 1 MCM OSM 2 0 15.9994
> ..etc
> 2127 OSM 1 MCM OSM 2127 0 15.9994
> 2128 OSM 1 MCM OSM 2128 0 15.9994
>
>
> [ system ]
> ; Name
> CH4 in MCM
>
> [ molecules ]
> ; Compound #mols
> MCM 1
> CH4 10
>
> CH4.itp file
>
> [ atomtypes ]
> ; type mass charge ptype c6 c12
> CH4 16.043 0.00 A 0.3732 1.24650457
> ;
> [ moleculetype ]
> ; name nrexcl
> CH4 2
>
> [ atoms ]
> ; nr type resnr residu atom cgnr charge mass
> 1 CH4 1 CH4 CH4 1 0.00 16.043
>
>
>
> .mdp file
>
> ;
> ; File 'mdout.mdp' was generated
> ; By user: jwillia4 (353773)
> ; On host: vlxhead2
> ; At date: Fri Jun 26 15:47:37 2009
> ;
> ; VARIOUS PREPROCESSING OPTIONS
> ; Preprocessor information: use cpp syntax.
> ; e.g.: -I/home/joe/doe -I/home/mary/hoe
> include = -I../top
> ; e.g.: -DI_Want_Cookies -DMe_Too
> define =
>
> ; RUN CONTROL PARAMETERS
> integrator = steep
> ; Start time and timestep in ps
> tinit = 0
> dt = 0.0001
> nsteps = 100000
> ; For exact run continuation or redoing part of a run
> ; Part index is updated automatically on checkpointing (keeps files
> separate)
> simulation_part = 1
> init_step = 0
> ; mode for center of mass motion removal
> comm-mode = linear
> ; number of steps for center of mass motion removal
> nstcomm = 1
> ; group(s) for center of mass motion removal
> comm-grps =
>
> ; LANGEVIN DYNAMICS OPTIONS
> ; Friction coefficient (amu/ps) and random seed
> bd-fric = 0
> ld-seed = 1993
>
> ; ENERGY MINIMIZATION OPTIONS
> ; Force tolerance and initial step-size
> emtol =
> emstep = 0.001
> ; Max number of iterations in relax_shells
> niter =
> ; Step size (ps^2) for minimization of flexible constraints
> fcstep =
> ; Frequency of steepest descents steps when doing CG
> nstcgsteep =
> nbfgscorr =
>
>
> ; OUTPUT CONTROL OPTIONS
> ; Output frequency for coords (x), velocities (v) and forces (f)
> nstxout = 100
> nstvout = 100
> nstfout = 0
> ; Output frequency for energies to log file and energy file
> nstlog = 100
> nstenergy = 100
> ; Output frequency and precision for xtc file
> nstxtcout = 100
> xtc-precision = 100
> ; This selects the subset of atoms for the xtc file. You can
> ; select multiple groups. By default all atoms will be written.
> xtc-grps =
> ; Selection of energy groups
> energygrps =
>
> ; NEIGHBORSEARCHING PARAMETERS
> ; nblist update frequency
> nstlist =
> ; ns algorithm (simple or grid)
> ns_type = grid
> ; Periodic boundary conditions: xyz, no, xy
> pbc = xyz
> periodic_molecules = yes
> ; nblist cut-off
> rlist = 1.7
>
> ; OPTIONS FOR ELECTROSTATICS AND VDW
> ; Method for doing electrostatics
> coulombtype = Cut-off
> rcoulomb-switch = 0
> rcoulomb = 1.7
> ; Relative dielectric constant for the medium and the reaction field
> epsilon_r =
> epsilon_rf =
>
> ; Method for doing Van der Waals
> vdw-type = Cut-off
> ; cut-off lengths
> rvdw-switch = 0
> rvdw = 1.7
> ; Apply long range dispersion corrections for Energy and Pressure
> DispCorr = No
> ; Extension of the potential lookup tables beyond the cut-off
> table-extension =
> ; Seperate tables between energy group pairs
> energygrp_table =
>
>
> ; Spacing for the PME/PPPM FFT grid
> fourierspacing = 0.12
> ; FFT grid size, when a value is 0 fourierspacing will be used
> fourier_nx = 0
> fourier_ny = 0
> fourier_nz = 0
> ; EWALD/PME/PPPM parameters
> pme_order =
> ewald_rtol = 1e-05
> ewald_geometry = 3d
> epsilon_surface = 0
> optimize_fft = yes
>
> ; IMPLICIT SOLVENT ALGORITHM
> implicit_solvent = No
>
>
> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
> ; Temperature coupling
> tcoupl = no
> ; Groups to couple separately
> tc-grps =
> ; Time constant (ps) and reference temperature (K)
> tau_t =
> ref_t =
>
> ; Pressure coupling
> Pcoupl = No
> Pcoupltype =
> ; Time constant (ps), compressibility (1/bar) and reference P (bar)
> tau-p =
> compressibility =
> ref-p =
> ; Scaling of reference coordinates, No, All or COM
> refcoord_scaling = no
> ; Random seed for Andersen thermostat
> andersen_seed =
>
>
> ; GENERATE VELOCITIES FOR STARTUP RUN
> gen_vel = no
> gen_temp = 300
> gen_seed = 173529
>
> ; OPTIONS FOR BONDS
> constraints = none
> ; Type of constraint algorithm
> constraint-algorithm = Lincs
> ; Do not constrain the start configuration
> continuation = no
> ; Use successive overrelaxation to reduce the number of shake iterations
> Shake-SOR = no
> ; Relative tolerance of shake
> shake-tol = 0.0001
> ; Highest order in the expansion of the constraint coupling matrix
> lincs-order = 4
> ; Number of iterations in the final step of LINCS. 1 is fine for
> ; normal simulations, but use 2 to conserve energy in NVE runs.
> ; For energy minimization with constraints it should be 4 to 8.
> lincs-iter = 1
> ; Lincs will write a warning to the stderr if in one step a bond
> ; rotates over more degrees than
> lincs-warnangle = 30
> ; Convert harmonic bonds to morse potentials
> morse = no
>
> ; ENERGY GROUP EXCLUSIONS
> ; Pairs of energy groups for which all non-bonded interactions are excluded
> energygrp_excl =
>
>
> ; Non-equilibrium MD stuff
> acc-grps =
> accelerate =
> freezegrps = MCM
> freezedim = Y Y Y
> cos-acceleration = 0
> deform =
>
>
>
>
>
>
More information about the gromacs.org_gmx-users
mailing list