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
--
gmx-users mailing list gmx-users@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-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
