[gmx-users] Energy calculation
Thanasis Koukoulas
koukoulas_th at yahoo.gr
Mon Jun 21 13:05:50 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/dc2751e0/attachment.html>
More information about the gromacs.org_gmx-users
mailing list