[gmx-users] Doubt about energies in a very simple system
Justin Lemkul
jalemkul at vt.edu
Tue Feb 24 19:31:12 CET 2015
On 2/24/15 11:42 AM, IÑIGO SAENZ wrote:
> 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?
>
Well, you're double-counting interactions, that's why the difference lines up.
A [pairs] directive is effectively a bonded potential. But you have no bonds,
so the atomic configuration is really just all nonbonded. But then you're
adding more bonded interactions on top...
>
> 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.
Everything is still dependent on the coordinates, regardless of the number of
parameters.
> 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.
>
The first thing to do is make sure you're properly approaching the SPE. You
showed the .mdp file, but *how* are you doing the calculation? If the answer
isn't mdrun -rerun on a configuration of equivalent precision to that of the
other programs, then you're not doing the SPE right.
http://www.gromacs.org/Documentation/How-tos/Single-Point_Energy
-Justin
--
==================================================
Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow
Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201
jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul
==================================================
More information about the gromacs.org_gmx-users
mailing list