[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