[gmx-users] Doubt about energies in a very simple system
IÑIGO SAENZ
inigo.saenz01 at estudiant.upf.edu
Tue Feb 24 17:42:46 CET 2015
Hi,
I designed a very simple system that is composed of only one Glutamine with
Tleap. I've transformed the corresponding .prmtop and .inpcrd into .top and
.gro files, using a conversor that I'm developing.
I attach .top .gro and .mdp files
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 yes 0.5 0.83333
[ atomtypes ]
;name bond_type mass charge ptype sigma epsilon
Amb
AA AA 0.00000 0.00000 A 3.25000e-01
7.11280e-01 ; 1.82 0.1700
AB AB 0.00000 0.00000 A 1.06908e-01
6.56888e-02 ; 0.60 0.0157
AC AC 0.00000 0.00000 A 3.39967e-01
4.57730e-01 ; 1.91 0.1094
AD AD 0.00000 0.00000 A 2.47135e-01
6.56888e-02 ; 1.39 0.0157
AE AE 0.00000 0.00000 A 2.64953e-01
6.56888e-02 ; 1.49 0.0157
AF AF 0.00000 0.00000 A 3.39967e-01
3.59824e-01 ; 1.91 0.0860
AG AG 0.00000 0.00000 A 2.95992e-01
8.78640e-01 ; 1.66 0.2100
[ moleculetype ]
;name nrexcl
sys 0
[ atoms ]
; nr type resi res atom cgnr charge mass ; qtot
1 AA 1 xxx AA 1 -0.516300 14.01000 ; -0.516
2 AB 1 xxx AB 2 0.293600 1.00800 ; -0.223
3 AC 1 xxx AC 3 0.039700 12.01000 ; -0.183
4 AD 1 xxx AD 4 0.110500 1.00800 ; -0.073
5 AC 1 xxx AC 5 0.056000 12.01000 ; -0.017
6 AE 1 xxx AE 6 -0.017300 1.00800 ; -0.034
7 AE 1 xxx AE 7 -0.017300 1.00800 ; -0.051
8 AC 1 xxx AC 8 0.013600 12.01000 ; -0.038
9 AE 1 xxx AE 9 -0.042500 1.00800 ; -0.080
10 AE 1 xxx AE 10 -0.042500 1.00800 ; -0.123
11 AF 1 xxx AF 11 0.805400 12.01000 ; 0.683
12 AG 1 xxx AG 12 -0.818800 16.00000 ; -0.136
13 AG 1 xxx AG 13 -0.818800 16.00000 ; -0.955
14 AF 1 xxx AF 14 0.536600 12.01000 ; -0.418
15 AG 1 xxx AG 15 -0.581900 16.00000 ; -1.000
[ pairs ]
1 6 1
1 7 1
1 8 1
1 15 1
2 4 1
2 5 1
2 14 1
3 9 1
3 10 1
3 11 1
4 6 1
4 7 1
4 8 1
4 15 1
5 12 1
5 13 1
5 15 1
6 9 1
6 10 1
6 11 1
6 14 1
7 9 1
7 10 1
7 11 1
7 14 1
8 14 1
9 12 1
9 13 1
10 12 1
10 13 1
[ exclusions ]
1 2 3 4 5 14
2 1 3
3 1 2 4 5 6 7 8 14 15
4 1 3 5 14
5 1 3 4 6 7 8 9 10 11 14
6 3 5 7 8
7 3 5 6 8
8 3 5 6 7 9 10 11 12 13
9 5 8 10 11
10 5 8 9 11
11 5 8 9 10 12 13
12 8 11 13
13 8 11 12
14 1 3 4 5 15
15 3 14
[ system ]
sys
[ molecules ]
; Compound nmols
sys 1
I have ommited [Bond] [angle] and [dihedral] section because they aren't
neccesary for my question. Now the .gro
15
1 xxx AA 1 0.3326 0.1548 -0.0000
1 xxx AB 2 0.3909 0.0724 -0.0000
1 xxx AC 3 0.3970 0.2846 -0.0000
1 xxx AD 4 0.3672 0.3400 -0.0890
1 xxx AC 5 0.3577 0.3654 0.1232
1 xxx AE 6 0.2497 0.3801 0.1241
1 xxx AE 7 0.3877 0.3116 0.2131
1 xxx AC 8 0.4267 0.4996 0.1195
1 xxx AE 9 0.5347 0.4850 0.1186
1 xxx AE 10 0.3967 0.5535 0.0296
1 xxx AF 11 0.3874 0.5805 0.2429
1 xxx AG 12 0.4595 0.5679 0.3454
1 xxx AG 13 0.2856 0.6542 0.2334
1 xxx AF 14 0.5486 0.2705 -0.0000
1 xxx AG 15 0.6009 0.1593 -0.0000
20.000000 20.000000 20.000000
and the SPE.mdp
dt = 0.001000
gen-vel = no
gen-temp = 0.000000
pbc = xyz
integrator = md
nsteps = 0
constraints = none
constraint-algorithm = SHAKE
rlist = 1.200000
rcoulomb = 1.200000
rvdw = 1.200000
coulombtype = Cut-off
vdwtype = Cut-off
vdw-modifier = None
coulomb-modifier = None
nstlog = 1
nstenergy = 1
nstcalcenergy = 1
cutoff-scheme = Verlet
comm-mode = Linear
nstcomm = 1
continuation = yes
My question is that when I execute a SPE of the original system in ACEMD, I
obtain the following VDW and coulomb energies
LJ 15.5754 KJ/Mol
COULOMB 123.8486 KJ/Mol
and in gromacs;
LJ-14 14.3788 (kJ/mol)
Coulomb-14 207.034 (kJ/mol)
LJ (SR) 29.9432 (kJ/mol)
Coulomb (SR) 165.181 (kJ/mol)
Question/Observation number 1:
Gromacs LJ (SR) includes LJ-14 energies in its result, no? I mean, LJ (SR)
= "the energy of the LJ interactions of those pairs of atoms that are
within the cutoff radius, including those that are enlisted in the [ pairs
] section without applying the fudgeLJ to them"
LJ-14 = "the sum of the energies of the interactions between the pairs of
atoms enlisted in the [ pairs ] section, applying to each of these
interactions the fudgeLJ scaling factor"
ACEMD LJ energy is equal to Gromacs LJ(SR) - Gromacs(LJ-14)
15.5754 is approximately equal to 29.9432 - 14.3788 = 15.5644
Does my "hypothesis" make sense?
Question/Observation number 2:
Continuing with the logic that I have exposed and supposing that I'm right,
ACEMD Coulomb energy should be equal to Gromacs Coulomb(SR) -
Gromacs(Coulomb-14). But this is false...
ACEMD Coulomb energy = 123.8486
Gromacs Coulomb energy ? 165.181 - 207.034 = -41.853
This doesn't make any sense for me... Gromacs and ACEMD coulomb energies
should be approximately the same. Coulomb energy only depends on the
distance between atoms and their charges, it is even simpler than the
calculation of VDW energy, which depends on more parameters ( sigma,
epsilon ) and it gives very distinct results.
I'm missing something and I don't know what, I desperately need help
I also attach the input file of ACEMD
amber on
parmfile system.prmtop
coordinates system.pdb
celldimension 200.0 200.0 200.0
outputunits kj
cutoff 12
pairlistdist 12
timestep 1
rigidbonds none
pairlistdist 12.0
exclude scaled1-4
1-4scaling 0.8333
pme off
hydrogenscale 1
switching off
Please, let me know if you need to know something about how I generate the
.top or the .gro, or if you need to see the .prmtop
I've tried everything and this is driving me crazy.
Thank you very much for your help.
More information about the gromacs.org_gmx-users
mailing list