[gmx-users] Energy miscalculation

Thanasis Koukoulas koukoulas_th at yahoo.gr
Mon Jun 21 13:18:57 CEST 2010


hi gmx-users
In order to discover the miscalculation that gromacs computes in a system of 450 molecules in vacuum and after a gmx-user guidance i ran a simulation of one molecule in vacuum just to calculate by hand all the energies. the mdp file that i use is 
; VARIOUS PREPROCESSING OPTIONStitle                    = Yocpp                      = /usr/bin/cppinclude                  = define                   = 
; RUN CONTROL PARAMETERSintegrator               = md; Start time and timestep in pstinit                    = 0dt                       = 0.001nsteps                   = 500000; For exact run continuation or redoing part of a runinit_step                = 0; mode for center of mass motion removalcomm-mode                = Linear; number of steps for center of mass motion removalnstcomm                  = 1; group(s) for center of mass motion removalcomm-grps                = 
; ENERGY MINIMIZATION OPTIONS; Force tolerance and initial step-sizeemtol                    = 10emstep                   = 0.01; Max number of iterations in relax_shellsniter                    = 20; Step size (1/ps^2) for minimization of flexible constraintsfcstep                   = 0; Frequency of steepest descents steps when doing CGnstcgsteep               = 1000nbfgscorr                = 10
; NEIGHBORSEARCHING PARAMETERS; nblist update frequencynstlist                  = 5; ns algorithm (simple or grid)ns_type                  = grid; Periodic boundary conditions: xyz (default), no (vacuum); or full (infinite systems only)pbc                      = xyz; nblist cut-off        rlist                    = 1.5domain-decomposition     = no
; OPTIONS FOR ELECTROSTATICS AND VDW; Method for doing electrostaticscoulombtype              = PMErcoulomb-switch          = 0rcoulomb                 = 1.5; Dielectric constant (DC) for cut-off or DC of reaction fieldepsilon-r                = 1; Method for doing Van der Waalsvdw-type                 = Cut-off; cut-off lengths       rvdw-switch              = 0rvdw                     = 1.5; Apply long range dispersion corrections for Energy and PressureDispCorr                 = EnerPres; Extension of the potential lookup tables beyond the cut-offtable-extension          = 1; Spacing for the PME/PPPM FFT gridfourierspacing           = 0.12; FFT grid size, when a value is 0 fourierspacing will be usedfourier_nx               = 0fourier_ny               = 0fourier_nz               = 0; EWALD/PME/PPPM parameterspme_order                = 4ewald_rtol  
             = 1e-05ewald_geometry           = 3depsilon_surface          = 0optimize_fft             = no
; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)implicit_solvent         = No
; OPTIONS FOR WEAK COUPLING ALGORITHMS; Temperature coupling  Tcoupl                   = nose-hoover; Groups to couple separatelytc-grps                  = System; Time constant (ps) and reference temperature (K)tau_t                    = 0.1ref_t                    = 298; Pressure coupling     Pcoupl                   = BerendsenPcoupltype               = isotropic; Time constant (ps), compressibility (1/bar) and reference P (bar)tau_p                    = 1compressibility          = 4.5e-5ref_p                    = 1.0; Random seed for Andersen thermostatandersen_seed            = 815131
; SIMULATED ANNEALING  ; Type of annealing for each temperature group (no/single/periodic)annealing                = no; Number of time points to use for specifying annealing in each groupannealing_npoints        = ; List of times at the annealing points for each groupannealing_time           = ; Temp. at each annealing point, for each group.annealing_temp           = 
; GENERATE VELOCITIES FOR STARTUP RUNgen_vel                  = yesgen_temp                 = 300gen_seed                 = 1993
; OPTIONS FOR BONDS    constraints              = all-bonds; Type of constraint algorithmconstraint-algorithm     = Lincs; Do not constrain the start configurationunconstrained-start      = no; Use successive overrelaxation to reduce the number of shake iterationsShake-SOR                = no; Relative tolerance of shakeshake-tol                = 1e-04; Highest order in the expansion of the constraint coupling matrixlincs-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 thanlincs-warnangle          = 30; Convert harmonic bonds to morse potentialsmorse                    = no
; ENERGY GROUP EXCLUSIONS; Pairs of energy groups for which all non-bonded interactions are excludedenergygrp_excl           = 


the top file is
[ defaults ]; nbfunc	comb-rule	gen-pairs	fudgeLJ	fudgeQQ  1		1		no              0      0.5

[ atomtypes ];   type      mass    charge    ptype       c6            c12    CH3     15.035     0.000    A       9.0638e-3      2.5206e-5    CH2     14.027     0.000    A       5.8108e-3      2.2071e-5    OS      15.999     0.000    A       8.8147e-4      4.2477e-7     OA      15.999     0.000    A       2.3465e-3      1.7802e-6    HO       1.008     0.000    A       0.0000e00      0.0000e00 

[ moleculetype ]; Name nrexclDRG      3
[ atoms ];   nr      type  resnr resid  atom  cgnr   charge     mass     1       CH3     1  DRG     CAA     1    0.000  15.035     2       CH2     1  DRG     CAC     1    0.000  14.027     3       CH2     1  DRG     CAE     1    0.000  14.027     4       CH2     1  DRG     CAF     1    0.000  14.027     5       CH2     1  DRG     CAG     2    0.000  14.027     6       CH2     1  DRG     CAI     2    0.250  14.027     7        OS     1  DRG     OAM     2   -0.500  15.999     8       CH2     1  DRG     CAK     2    0.250  14.027     9       CH2     1  DRG     CAJ     3    0.250  14.027     10        OS     1  DRG     OAL     3   -0.500  15.999      11       CH2     1  DRG     CAH     3    0.250  14.027      12       CH2     1  DRG     CAD     4    0.265
  14.027    13        OA     1  DRG     OAB     4   -0.700  15.999      14        HO     1  DRG     HAA     4    0.435   1.008
[ bonds ];  ai  aj funct           c0           c11       2	1	   0.15400       517488.482       3	1	   0.15400       517488.483       4       1          0.15400       517488.484       5       1          0.15400       517488.485       6       1          0.15400       517488.486       7       1          0.14100       618809.047       8       1          0.14100       618809.048       9       1          0.15400       517488.489      10       1          0.14100       618809.0410     11       1          0.14100       618809.0411     12       1          0.15400       517488.4812     13       1          0.14300       618809.0413     14       1          0.09450       618809.04
[ angles ];  ai    aj    ak       funct   c0      c1    1     2     3       1       114.00  519.65    2     3     4       1       114.00  519.65    3     4     5       1       114.00  519.65    4     5     6       1       114.00  519.65    5     6     7       1       112.00  418.22    6     7     8       1       112.00  502.19    7     8     9       1       112.00  418.22    8     9    10       1       112.00  418.22    9    10    11       1       112.00  502.19   10    11    12       1       112.00  418.22   11    12    13       1       109.47  419.05   12    13    14       1       108.50  460.62
[ pairs ]; ai aj funct  1   4   1 0 0 0 0 2   5   1 0 0 0 0 3   6   1 0 0 0 0 4   7   1 0 0 0 0 5   8   1 0 0 0 0 6   9   1 0 0 0 0 7  10   1 0 0 0 0 8  11   1 0 0 0 0 9  12   1 0 0 0 010  13   1 0 0 0 011  14   1 0 0 0 0;10 14    2 0 0 0 0 3.325789e-7

[ dihedrals ]; ai aj ak al funct c0           c1           c2            c3           c4           c5   1  2  3  4  3    8.231078 16.952626 1.133928 -26.317632 0.000000e+00 0.000000e+00   2  3  4  5  3    8.231078 16.952626 1.133928 -26.317632 0.000000e+00 0.000000e+00   3  4  5  6  3    8.231078 16.952626 1.133928 -26.317632 0.000000e+00 0.000000e+00   4  5  6  7  3    6.983076 17.736182 0.886988 -25.606246 0.000000e+00 0.000000e+00   5  6  7  8  3    7.949051  7.892513 2.722990 -18.564553 0.000000e+00 0.000000e+00   6  7  8  9  3    7.949051  7.892513 2.722990 -18.564553 0.000000e+00 0.000000e+00   7  8  9 10  3    8.368267 25.104800 4.184175 -33.473067 0.000000e+00 0.000000e+00   8  9 10 11  3    7.949051  7.892513 2.722990 -18.564553 0.000000e+00 0.000000e+00   9 10 11 12  3    7.949051  7.892513 2.722990 -18.564553 0.000000e+00 0.000000e+00  10 11 12 13
  3    8.368267 25.104800 4.184175 -33.473067 0.000000e+00 0.000000e+00  11 12 13 14  3    2.822015  2.943074 0.485066  -6.250155 0.000000e+00 0.000000e+00
[ system ]ciej 
[ molecules ]DRG	 1
and the gro file is
PRODRG COORDS   14    1DRG  CAA      1   0.313   0.879   0.001    1DRG  CAC      2   0.381   0.742   0.002    1DRG  CAE      3   0.534   0.757   0.000    1DRG  CAF      4   0.601   0.620   0.001    1DRG  CAG      5   0.754   0.634  -0.001    1DRG  CAI      6   0.821   0.497   0.001    1DRG  OAM      7   0.964   0.513  -0.001    1DRG  CAK      8   1.026   0.383   0.001    1DRG  CAJ      9   1.177   0.399  -0.001    1DRG  OAL     10   1.238   0.268   0.001    1DRG  CAH     11   1.382   0.282   0.000    1DRG  CAD     12   1.446   0.143   0.002    1DRG  OAB     13   1.588   0.156   0.001    1DRG  HAA     14   1.630   0.066   0.002   4.00000   4.00000   4.00000
as you can see at top file there is a line in pairs which syntaxis is
10 14    2 0 0 0 0 3.325789e-7
it is an extra interaction that must be applied to the molecule according to the formula 3.325789e-7/r^12 the only way i could think importing  this interaction as i have mentioned several times was the above. 
i did three simulations
1) leaving the top file as it is above2) deleting the  [pairs] and all the lines below3) deleting only the last line of the [pairs]
the results with usage of the command gmxdump -s topol.tpr | lessfor step 0 are below

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100621/0b338746/attachment.html>


More information about the gromacs.org_gmx-users mailing list