Project

General

Profile

Bug #294

Erroneus Buckingham pairs

Added by Luca De over 11 years ago. Updated over 11 years ago.

Status:
Closed
Priority:
Normal
Assignee:
Erik Lindahl
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

I want to perform a Buckingham simulation using penthane moelcule as an example.
I've selected nbfunc = 2 in topol.top and comb-rule=2

Now I want to add Lennard-Jones pairs. First, I've explicitly excluded all the pairs in [exclusions]. Then I've added the following strings:

[pairs ]
1 5 1 0.33 0.69

to activate 1-5 intramolecular LJ, with epsilon and sigma values, respectively.
The bug I've noted is this: instead of a normal Lennard Jones, the potential calculated by Gromacs has this structure:

4 epsi * [ (2/3 1/sqrt(3))^12 * (sigma/r)^12 - (sigma/r)^6]

whereas the standard Lennard Jones is

4 epsi * [ (sigma/r)^12 - (sigma/r)^6]

so the only one modification is in the prefactor (2/3 1/sqrt(3))^12; it must be noted that this is my approximation, and this value can vary a little if sigma (and only sigma) is changed. I also tested a lot of sigma / epsilon values and it seems to confirm my approximation.

The same prefactor (whatever it is) is present again with "comb-rule = 1", i.e. with C12, C6 potential. In this case it is the C12 term to be modified by the same prefactor.

Infact if I add that [pairs] section to a Buckingham simulation, the simulation quickly become instable and explodes.
If I turn off the whole [pairs] section no problem occurs (but obviously I need 1-4 terms). If I exploit the "genpairs" key, Gromacs calculates an exact intramolecular potential, but I do need to add my proper terms.
If I use the same conf.gro but "normal" LJ topology, no error occurs.

I attach a 10 molecule system with:
1. the initial configuration (conf.gro)
2. The buckingham topology (topol.buck.top)
3. The LJ topology (topol.LJ.top)
4. grompp.mdp

It is probably better to run testes on 1 molecules systems (as I done).

buckingham.tgz (12.3 KB) buckingham.tgz test case Luca De, 02/11/2009 04:23 PM

History

#1 Updated by Luca De over 11 years ago

Created an attachment (id=352)
test case

#2 Updated by Berk Hess over 11 years ago

For you pairs you have: sigma=0.33, eps=0.69
I ran grompp (of 4.0.3) and then gmxdump -s, which gives:
functype8=LJ14, c6A= 3.56445252e-03, c12A= 4.60337742e-06, c6B= 3
.56445252e-03, c12B= 4.60337742e-06
These are correct c6 and c12 values.

How did you conclude that the LJ pairs are incorrect?

Berk

#3 Updated by Luca De over 11 years ago

(In reply to comment #2)

For you pairs you have: sigma=0.33, eps=0.69
I ran grompp (of 4.0.3) and then gmxdump -s, which gives:
functype8=LJ14, c6A= 3.56445252e-03, c12A= 4.60337742e-06, c6B= 3
.56445252e-03, c12B= 4.60337742e-06
These are correct c6 and c12 values.

How did you conclude that the LJ pairs are incorrect?

Berk

I agree with you with the gmxdump output, but I still have problem with the simulation.
So, I'll send you a 1 molecole pentane simulation. The distance between sites 1-5 is 3.70Angstrom (I attached an xyz molden file).
There is only one 1-5 interaction.
For simplicity, I set sigma = 1 and epsilon = 0.69; comb rule = 2
From gmxdump -s C6B = C12B = 2.76 and it is OK.
But the md.log file output is: LJ-14 = -3.50895e+02

With correct sigma/epsilon it should be:
4 * 0.69 *( (1/0.37)12 - (1/0.37)6) = 418189 (even the sign is opposite!!)
or altrnatively
2.76 / (0.3712) - 2.76 / (0.376) = circa 419000
A problem seems to be.
Am I wrong?
Luca De Gaetani

#4 Updated by Berk Hess over 11 years ago

This was a tricky bug.
LJ14 used the exponential repulsion table with buckingham
instead of r^-12.
I fixed it for 4.0.4.

To fix it, line 789 of src/mdlib/tables.c should be changed from
if (fr->bBHAM) {
to
if (fr->bBHAM && !b14only) {

Thanks for reporting this.

Berk

Also available in: Atom PDF