[gmx-users] Doubt about energies in a very simple system
IÑIGO SAENZ
inigo.saenz01 at estudiant.upf.edu
Tue Feb 24 21:23:34 CET 2015
Hi Justin,
I always do the SPE as follows:
grompp -f SPE.mdp -p sys.top -c sys.gro
and after that I simply execute mdrun, i didn't know about the mdrun -rerun
function.
Now I have done: mdrun -s topol.tpr -rerun sys.gro
but the energy results are the exactly the same.
Thank you.
2015-02-24 19:30 GMT+01:00 Justin Lemkul <jalemkul at vt.edu>:
>
>
> 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