[gmx-users] potential energy NaN and strange dependence on cut-offs
Manik Mayur
manik.mayur at gmail.com
Tue Sep 8 19:37:27 CEST 2009
2009/9/8 Mark Abraham <Mark.Abraham at anu.edu.au>
> 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
>>
>
dt = 0.1 fs are you sure??
> 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 =
>>
>>
>>
>>
>>
>>
>> _______________________________________________
> gmx-users mailing list gmx-users at 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-request at gromacs.org.
>
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
Manik Mayur
Graduate student
Microfluidics Lab
Dept. of Mechanical Engg.
IIT Kharagpur
INDIA
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090908/edd366c2/attachment.html>
More information about the gromacs.org_gmx-users
mailing list