[gmx-users] potential energy NaN and strange dependence on cut-offs

David van der Spoel spoel at xray.bmc.uu.se
Tue Sep 8 19:57:13 CEST 2009


Manik Mayur wrote:
> 
> 2009/9/8 Mark Abraham <Mark.Abraham at anu.edu.au 
> <mailto: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
> 

Check your tpr file: are the C6 and C12 reasonable?
If you get NaN something is wrong.

>         ;
>         ; 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
>     <mailto: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
>     <mailto: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
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> 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


-- 
David van der Spoel, Ph.D., Professor of Biology
Molec. Biophys. group, Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205. Fax: +4618511755.
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se



More information about the gromacs.org_gmx-users mailing list