[gmx-users] Doubt about energies in a very simple system
Mark Abraham
mark.j.abraham at gmail.com
Tue Feb 24 20:54:18 CET 2015
Hi,
In addition to Justin's points about comparing apples with apples on the
same coordinates, by default in the Verlet scheme GROMACS also shifts the
potential to be zero at the cutoff, so that the potential is actually the
integral of the force (mentioned in manual 3.4.2, along with a sea of other
things, sorry). You can find several threads in the archive where people
have tried to compare with other MD programs, chiefly ACEMD, IIRC. So that
effect is more or less noteworthy depending on the system and settings.
Mark
On Tue, Feb 24, 2015 at 7:30 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>
> 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
>
> ==================================================
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>
More information about the gromacs.org_gmx-users
mailing list