Project

General

Profile

Bug #98

error in calculate the potential between a pair of particles by using tabulated potentials

Added by Dongsheng Zhang over 13 years ago. Updated over 13 years ago.

Status:
Closed
Priority:
High
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

I have written a short program to calculate the non-bonded potential energy. I
found my result is different from gromacs. After I did several tests, I think
gromacs has made a mistake somewhere. Please refer my email in the gmx-developer
email list to get the detailed information of my tests. The email title is
"likely an error in gromacs".

The submitted table is table 5 in my email. You can vary the entries to see how
gromacs to respond.

rerun_test.tpr (3.11 MB) rerun_test.tpr tpr file Dongsheng Zhang, 08/20/2006 06:52 AM
table_PEO_EOP.xvg (121 KB) table_PEO_EOP.xvg tabulated potentials. When you check the energy, please choos LJ(SR)_PEO_EOP. Dongsheng Zhang, 08/20/2006 06:55 AM
table.xvg (121 KB) table.xvg tabulated potentials Dongsheng Zhang, 08/20/2006 06:59 AM
rerun_1.gro (395 KB) rerun_1.gro rerun input Dongsheng Zhang, 08/20/2006 03:44 PM
rerun.tpr (3.11 MB) rerun.tpr run input tpr file Dongsheng Zhang, 08/20/2006 11:02 PM
rerun_test_PEO_EOP.ndx (270 KB) rerun_test_PEO_EOP.ndx ndx file for making tpr file Dongsheng Zhang, 08/21/2006 06:06 PM

History

#1 Updated by Dongsheng Zhang over 13 years ago

Created an attachment (id=52)
tpr file

you need to use table***.xvg files to run it.

#2 Updated by Dongsheng Zhang over 13 years ago

Created an attachment (id=53)
tabulated potentials

#3 Updated by Dongsheng Zhang over 13 years ago

Created an attachment (id=54)
tabulated potentials

Please copy table.xvg to:
tablep.xvg
table_CCT_CCT.xvg
table_CNT_CCT.xvg
table_CNT_CNT.xvg
table_CNT_EOP.xvg
table_CNT_POP.xvg
table_EOP_CCT.xvg
table_EOP_EOP.xvg
table_EOP_POP.xvg
table_PEO_CCT.xvg
table_PEO_CNT.xvg
table_PEO_PEO.xvg
table_PEO_POP.xvg
table_PEO_PPO.xvg
table_POP_CCT.xvg
table_POP_POP.xvg
table_POP_PPO.xvg
table_PPO_CCT.xvg
table_PPO_CNT.xvg
table_PPO_EOP.xvg
table_PPO_POP.xvg
table_PPO_PPO.xvg

#4 Updated by David van der Spoel over 13 years ago

please provide a simpler example, or at least clear instructions what to do.
- command line for mdrun
- the energies you get from mdrun
- the energies you think one should get
- a tpr file with only one step and preferably for a much smalle system and
without position restraints
- the source code to your table energy program
- by the way, your system has a non-integral charge

#5 Updated by Dongsheng Zhang over 13 years ago

(In reply to comment #4)

please provide a simpler example, or at least clear instructions what to do.
- command line for mdrun
- the energies you get from mdrun
- the energies you think one should get
- a tpr file with only one step and preferably for a much smalle system and
without position restraints
- the source code to your table energy program
- by the way, your system has a non-integral charge

The command is :
mdrun -s rerun_test.tpr -rerun rerun_1.gro -s e rerun -table -tablep

gromacs gives energy of LJPEO_EOP as 7.212971, it should be 7.86254128.
The two particles are #3483 (coordinates:42.04700 37.95200 42.47800
) and #8985 (corrdinates:42.99200 37.88300 42.77700). The distance between
them is 0.993572.

I am sorry I don't know what you mean about "the source code to your table
energy program".

It is true that my system has a non-integral charge. It is a trick to calculate
the 1-4 interaction (read tablep.xvg). My force field is course-grained model.
There is no COUL term.

If you need more information or something is not clear to you, please let me
know. Thank you!

#6 Updated by Dongsheng Zhang over 13 years ago

Created an attachment (id=55)
rerun input

#7 Updated by Dongsheng Zhang over 13 years ago

(In reply to comment #5)

(In reply to comment #4)

please provide a simpler example, or at least clear instructions what to do.
- command line for mdrun
- the energies you get from mdrun
- the energies you think one should get
- a tpr file with only one step and preferably for a much smalle system and
without position restraints
- the source code to your table energy program
- by the way, your system has a non-integral charge

The command is :
mdrun -s rerun_test.tpr -rerun rerun_1.gro -s e rerun -table -tablep

gromacs gives energy of LJPEO_EOP as 7.212971, it should be 7.86254128.
The two particles are #3483 (coordinates:42.04700 37.95200 42.47800
) and #8985 (corrdinates:42.99200 37.88300 42.77700). The distance between
them is 0.993572.

I am sorry I don't know what you mean about "the source code to your table
energy program".

It is true that my system has a non-integral charge. It is a trick to calculate
the 1-4 interaction (read tablep.xvg). My force field is course-grained model.
There is no COUL term.

If you need more information or something is not clear to you, please let me
know. Thank you!

Sorry, the second particle should be #8986. I don't know if you can get the
information of energy group from tpr file. In my energy group, PEO group only
has particle 3483, EOP group only has particle 8986.

#8 Updated by David van der Spoel over 13 years ago

0.000000 78.622437
I get the above energy, which is correct according to the information you
provided. Question then is, why do you get something else?

What platform are you running on and how did you get your mdrun binary

#9 Updated by Dongsheng Zhang over 13 years ago

(In reply to comment #8)

0.000000 78.622437
I get the above energy, which is correct according to the information you
provided. Question then is, why do you get something else?

What platform are you running on and how did you get your mdrun binary

Thank you for your help. In fact, the energy is slightly different from my
calculation result (78.625412). I know the difference is very small, so I won't
argue about it. In fact I upload a wrong tpr file. The new file rerun.tpr will
give you a big difference (it will be 77.973167 fromacs gromacs calculation).
The difference between rerun.tpr and rerun_test.tpr is as following:

When I made rerun_test.tpr, the index file only has group "protein", "PEO" and
"EOP". The index file has other groups (such as CNT, CCT, PPO ...) when I made
rerun.tpr. That's the reason I asked you to copy a bunch of table***.xvg files.
Unfortunately I uploaded a wrong tpr file. Sorry about that. I don't expect the
energy group assignment will change my calculation. But it does. Please try to
run as:
mdrun -s rerun.tpr -rerun rerun_1.gro -e rerun_new
g_energy -f rerun_new
choose LJ-SR:PEO-EOP to get the energy.

#10 Updated by Dongsheng Zhang over 13 years ago

Created an attachment (id=56)
run input tpr file

#11 Updated by Berk Hess over 13 years ago

In the tpr of 2006-08-20 23:02 none of the EOP atoms
is within 1.4 nm of the PEO atom.

#12 Updated by Dongsheng Zhang over 13 years ago

(In reply to comment #11)

In the tpr of 2006-08-20 23:02 none of the EOP atoms
is within 1.4 nm of the PEO atom.

Your find is really surprises me. The following information is from
rerun_1.gro:

93PEO C3 3483 42.047 37.952 42.478
691EOP C2 8986 42.992 37.883 42.777

#3483 is in energy group "PEO", and #8986 is in energy group "EOP". Therefore
gromacs can use them to calculate the interaction between PEO and EOP. The
distance between them is 0.993572. Could you please tell me how you find there
is no EOP within 1.4 nm of the PEO atom? By the way, could you please tell me
how to understand grpnrs=[0 2 0 0 0 0 0 0 0 0]?

#13 Updated by Berk Hess over 13 years ago

with gmxdump -s rerun.tpr or editconf -f rerun.tpr I get:
93PEO C4 3484 -9.749 55.348 42.434 -0.0249 -0.0247 -0.3532
691EOP C2 8986 -9.035 56.970 43.742 0.0584 0.2507 0.2936

The grpnrs thing is not really documented.
The second column are the energy group numbers.
They index into the "grp" group names, search for PEO in:
gmxdump -s rerun.tpr | less

#14 Updated by Dongsheng Zhang over 13 years ago

(In reply to comment #13)

with gmxdump -s rerun.tpr or editconf -f rerun.tpr I get:
93PEO C4 3484 -9.749 55.348 42.434 -0.0249 -0.0247 -0.3532
691EOP C2 8986 -9.035 56.970 43.742 0.0584 0.2507 0.2936

The grpnrs thing is not really documented.
The second column are the energy group numbers.
They index into the "grp" group names, search for PEO in:
gmxdump -s rerun.tpr | less

Yes. You are right. If you look up the corrdinates in tpr file. You result is
correct. Becasue the coordinates in tpr is from the gro file when we use command
grompp -c *gro. But when we do calculations, mdrun should use rerun_1.gro file,
is it right?

#15 Updated by Berk Hess over 13 years ago

Ah, sorry I overlooked the fact that you used mdrun -rerun.

But I guess the problem is that you have charge-groups and therefore
a whole charge-group (4 atoms) is within the cut-off and the energy
of 4 atom pairs ends up in LJ-SR:PEO-EOP

Berk.

#16 Updated by Dongsheng Zhang over 13 years ago

(In reply to comment #15)

Ah, sorry I overlooked the fact that you used mdrun -rerun.

But I guess the problem is that you have charge-groups and therefore
a whole charge-group (4 atoms) is within the cut-off and the energy
of 4 atom pairs ends up in LJ-SR:PEO-EOP

Berk.

Sorry I don't understand why you say 4 atoms (Initially I make four atoms as one
residue and one charge group). Here we only consider two atoms. And each atom is
in one specific charge group (I made it on purpose in the itp file). Let me make
it clear. Only atom 3483 is in charge group 68 (the residue group is 11) in the
polymer.itp, only atom 8986 is in charge group 67 (the residue group is 65).
Only atom 3483 is in energy group PEO, and only atom 8986 is in energy group EOP.

I hope my explanation is clear to you. If not or you need more information,
please let me know.

By the way, could you please try to run it to see what you can get? If your
result is not quite right, you might add some output statement in the source
code to see what is going on.

Thank you for your help!

Dongsheng

#17 Updated by David van der Spoel over 13 years ago

atom 8986 is in an energy group with 800 atoms in rerun.tpr
in the first tpr the charge groups are larger than one atom but in rerun.tpr
they are not.

#18 Updated by Berk Hess over 13 years ago

I have can reproduce your energy when I change only the 0.992 and 0.994
lines of the table.

My previous explanation was incorrect.
The real explanation is much simpler.
I did not realize your cut-off was 1.4.
There are many atoms within a distance of 1.4.
When I make the whole table zero, except for r=0.992 and 0.994
I get 7.862244

#19 Updated by Dongsheng Zhang over 13 years ago

(In reply to comment #18)

I have can reproduce your energy when I change only the 0.992 and 0.994
lines of the table.

My previous explanation was incorrect.
The real explanation is much simpler.
I did not realize your cut-off was 1.4.
There are many atoms within a distance of 1.4.
When I make the whole table zero, except for r=0.992 and 0.994
I get 7.862244

Yes. it is true that you can get 7.862244 when you change the table. I don't
understand why other entries in the table affect the calculation. It should be
irrelevant. Could you please check the tpr file to see if PEO energy group only
have atom 3483 and and EOP energy group only has atom 8986 (Please pay attention
that the atom number in tpr starts from 0, not 1)? If you check rerun_test.tpr,
there should be only two energy groups PEO and EOP besides protein. If you check
rerun.tpr, there are sveral energy groups (energygrps = PEO PPO CNT EOP POP CCT)
besides protein. If my enrgy group is correct is tpr file. I really want to ask
you a favor to output more information during the calculation.
I have attached my ndx file which I used to make my tpr file

#20 Updated by Dongsheng Zhang over 13 years ago

Created an attachment (id=62)
ndx file for making tpr file

#21 Updated by Berk Hess over 13 years ago

Sorry, but when I use rerun_test.tpr instead of rerun.tpr
I always get 7.862244, both with all other entries zero
or at the original values.

#22 Updated by Dongsheng Zhang over 13 years ago

(In reply to comment #21)

Sorry, but when I use rerun_test.tpr instead of rerun.tpr
I always get 7.862244, both with all other entries zero
or at the original values.

Yes. David got 7.862244 as well when he used rerun_test.tpr. That is the reason
I upload rerun.tpr. Please try rerun.tpr with the original table, then you can
get a different result. The difference between rerun.tpr and rerun_test.tpr is
following:
rerun.tpr: energygrps = PEO PPO CNT EOP POP CCT in the mdp file
rerun_test.tpr: energygrps = PEO PPO in the mdp file

Thank you for your help!

Dongsheng

#23 Updated by David van der Spoel over 13 years ago

With rerun_test.tpr I get: 78.622437
With rerun.tpr I get: 77.973167

In rerun.tpr there are two groups with one atom: 0 and 5, corresponding to atoms
3483 and 6292

In rerun_test.tpr there also are two groups with one atom, 0 and 1 corresponding
to 3483 and 8986

I conclude therefore that we are comparing apples and oranges.

#24 Updated by Dongsheng Zhang over 13 years ago

(In reply to comment #23)

With rerun_test.tpr I get: 78.622437
With rerun.tpr I get: 77.973167

In rerun.tpr there are two groups with one atom: 0 and 5, corresponding to atoms
3483 and 6292

In rerun_test.tpr there also are two groups with one atom, 0 and 1 corresponding
to 3483 and 8986

I conclude therefore that we are comparing apples and oranges.

atom 6292 is in group CCT, not EOP. so the energy is not for 3483 and 6292. I
have the answer now. In short I made a mistake to make the ndx file. I
unconsciously have two blocks for [EOP] group. The first entry has a bunch of
atoms. The second blocks has only one atom which is atom 8986 (This is what I
expect). gromacs takes the first block as my EOP group even though it is in
front of the second block (I don't know why the second block does not overwrite
the first one). Here I would like to kindly give one suggestion, how about
giving a warning message if there are two blocks with the same group name in
ndx file?

David and Berk, Thank you very much for your great help! I am sorry to give you
such a mess. All the best!

Also available in: Atom PDF