## Bug #294

### Erroneus Buckingham pairs

**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).

### 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:

functype^{8}=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:

functype^{8}=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