Alchemical FEP and NBFIX

User discussions
,
1

GROMACS version: 2021.2
GROMACS modification: No

GROMACS uses the [ nonbond_params ] section of the topology to implement custom mixing rules (also referred to as NBFIX in CHARMM). How do these custom mixing rules (which are defined based on atom type) get treated when doing an alchemical FEP transformation (which is also based on atom type, ie atom_i → atom_j)

I’m trying to define my own mixing rules for a few atometypes and then trying to look at relative energy between the two. I’m doing this via an alchemical transformation, but I’m not sure how the NBFIX terms get applied. I got the below results from BAR without the custom mixing rules:

Final results in kJ/mol:

point 0 - 1, DG 0.00 +/- 0.00
point 1 - 2, DG -0.00 +/- 0.00
point 2 - 3, DG 0.00 +/- 0.00
point 3 - 4, DG 0.00 +/- 0.00
point 4 - 5, DG -0.85 +/- 0.01
point 5 - 6, DG -0.95 +/- 0.01
point 6 - 7, DG -1.05 +/- 0.01
point 7 - 8, DG -1.14 +/- 0.01
point 8 - 9, DG -1.22 +/- 0.01
point 9 - 10, DG -1.32 +/- 0.01
point 10 - 11, DG -1.43 +/- 0.01
point 11 - 12, DG -1.56 +/- 0.01
point 12 - 13, DG -1.67 +/- 0.01
point 13 - 14, DG -1.78 +/- 0.01
point 14 - 15, DG -1.89 +/- 0.01
point 15 - 16, DG -1.99 +/- 0.01
point 16 - 17, DG -2.08 +/- 0.01
point 17 - 18, DG -2.18 +/- 0.01
point 18 - 19, DG -2.28 +/- 0.01
point 19 - 20, DG -2.39 +/- 0.01
point 20 - 21, DG -2.48 +/- 0.00
point 21 - 22, DG -2.57 +/- 0.00
point 22 - 23, DG -2.67 +/- 0.01
point 23 - 24, DG -2.78 +/- 0.00

total 0 - 24, DG -36.28 +/- 0.05

And this is the result when I included my custom mixing rules:

Final results in kJ/mol:

point 0 - 1, DG -0.00 +/- 0.00
point 1 - 2, DG -0.00 +/- 0.00
point 2 - 3, DG 0.00 +/- 0.00
point 3 - 4, DG 0.00 +/- 0.00
point 4 - 5, DG 1780.79 +/- 0.67
point 5 - 6, DG 0.77 +/- 0.02
point 6 - 7, DG 0.51 +/- 0.04
point 7 - 8, DG 0.25 +/- 0.03
point 8 - 9, DG -0.06 +/- 0.03
point 9 - 10, DG -0.21 +/- 0.02
point 10 - 11, DG -0.35 +/- 0.03
point 11 - 12, DG -0.44 +/- 0.02
point 12 - 13, DG -0.56 +/- 0.02
point 13 - 14, DG -0.70 +/- 0.02
point 14 - 15, DG -0.83 +/- 0.01
point 15 - 16, DG -1.06 +/- 0.01
point 16 - 17, DG -1.25 +/- 0.02
point 17 - 18, DG -1.41 +/- 0.03
point 18 - 19, DG -1.50 +/- 0.01
point 19 - 20, DG -1.61 +/- 0.03
point 20 - 21, DG -1.78 +/- 0.02
point 21 - 22, DG -1.92 +/- 0.02
point 22 - 23, DG -2.07 +/- 0.03
point 23 - 24, DG -2.18 +/- 0.02

total 0 - 24, DG 1764.39 +/- 0.78

The first 4 lambda steps in my setup are changing the mass and starting at lambda 5 I start varying the VDW parameters. I assume the large jump at “point 4-5” is something going on with changing the atomtype which changes the mixing rules all at once while the VDW parameters are changed gradually through the rest of the lambda steps but I’m not sure. Any insights would be appreciated!

Thank you!

2

I don’t understand your results. I assume that in the second case you are using your custom mixing rules for all lambda values, right?
I don’t understand what can cause such a large free energy difference.

3

Sorry for the poor explanation on my part. I’m testing a system where I’m alchemically changing a single atom into another single atom of a different type (typeA → typeB) The first case used the standard mixing rules (for both typeA and typeB) with the surrounding solvent and the second used my attempt at custom mixing rules (again for both typeA and typeB) with the surrounding solvent (via the nonbond_params section of the topology). Both typeA and typeB have the same atomic charge.

My lambda pathway looks like the following:

mass_lambdas = 0.00 0.25 0.50 0.75 1.00 1.00 …. 1.00
vdw_lambdas = 0.00 0.00 0.00 0.00 0.00 0.05 …. 1.00

The big jump on question happens on the first instance where the vdw parameters are being altered.

Thanks again
Kevin

4

Yes, that is what I understood. But this can’t explain the massive free-energy difference. What is the change in parameters between state A and B?