Inconsistent results when using pair table interactions

[CreditDefaultSwap]

Hello,

I am applying a protocol on a simple fluid system whose pair interactions use some customized function that I input using a pair.Table method. In particular, I need to implement a protocol that continously changes the cut off radius, r_cut, of this interaction. However I noticed that if I redefined the pair.Table object with a new cut-off radius every time, some wierd 'flickering' of the energy would occur, where the values would jump periodically after a certain number of steps, with longer simulation times (meaning slower changes in r_cut) leading to longer periods in between the jumps.

To get around that I am using a table.Pair object that instead of changing r_cut explictly, now uses a fixed r_fix at some large distance and pads the table with 0s past the r_cut distance. Problem is that pair-based quantities differ strongly quantitatively between when using the padded table interaction with the one using a pair.Table object defined with the same r_cut value, where this difference depends on the value of r_fix. I understand that the pair.Table object applies linear interpolation between points, but even when holding constant for the number of points in the table, and the corresponding distance between points, the results are nonetheless noticeably different.

Could anyone help me understand why this is the case? What could I be doing wrong with the padded table pair potential that causes these random issues, or is this approach with the padded potential simply wrong because it breaks the internal (black box) interpolation scheme somehow?

Thank you for any tips or information

[joaander]

Yes, I would expect that changing the pair potential and/or cutoff radius would introduce a discontinuity in the potential energy of the system. That discontinuity would occur each time you make the change. Why do you expect the change to be continuous?

Both your workaround and original method should produce nearly identical results. The only difference is that the first method will result in a changing step size $\Delta r$ which the second keeps fixed.

I confirm the expected behavior with this example:

# For demonstration purposes only. U and F are not consistent!

import hoomd
simulation = hoomd.util.make_example_simulation()

nlist = hoomd.md.nlist.Cell(buffer=0.4)
table1 = hoomd.md.pair.Table(nlist)
table1.params[('A', 'A')] = dict(r_min=0, U=[1.0], F=[1.0])
table1.r_cut[('A', 'A')] = 2.0

table2 = hoomd.md.pair.Table(nlist)
table2.params[('A', 'A')] = dict(r_min=0, U=[1.0, 0.0], F=[1.0, 0.0])
table2.r_cut[('A', 'A')] = 4.0

simulation.operations.computes.extend([table1, table2])
simulation.run(0)

print("Energies with r_cut=2.0:", table1.energy, table2.energy)

simulation.run(1)
table1.r_cut[('A', 'A')] = 4.0
table2.params[('A', 'A')] = dict(r_min=0, U=[1.0, 0.5], F=[1.0, 0.5])

print("Energies with r_cut=4.0:", table1.energy, table2.energy)

As expected, this script prints:

Energies with r_cut=2.0: 0.0 0.0
Energies with r_cut=4.0: 0.5 0.5

You don't provide a similar minimal working example, so I cannot even guess at what you might be doing wrong in your implementation of the 2nd approach.

Also, the interpolation scheme is not a black box. The documentation clearly states:

Provide F_table(r) and U_table(r) on evenly spaced grid points points between r_min​ and r_cut​. Table linearly interpolates values when r lies between grid points and between the last grid point and r=r_cut.

If you don't believe the documentation, you can read the code itself:

hoomd-blue/hoomd/md/EvaluatorPairTable.h

Lines 162 to 185 in 26810d6

	const Scalar rcut = fast::sqrt(rcutsq);
	const Scalar delta_r = (rcut - rmin) / static_cast<Scalar>(width);
	// precomputed term
	const Scalar value_f = (r - rmin) / delta_r;
	// compute index into the table and read in values
	unsigned int value_i = static_cast<unsigned int>(slow::floor(value_f));
	// unpack the data
	const Scalar V0 = V_table[value_i];
	const Scalar F0 = F_table[value_i];
	Scalar V1 = 0;
	Scalar F1 = 0;
	if (value_i + 1 < width)
	{
	V1 = V_table[value_i + 1];
	F1 = F_table[value_i + 1];
	}
	// compute the linear interpolation coefficient
	const Scalar f = value_f - Scalar(value_i);
	// interpolate to get V and F;
	const Scalar V = V0 + f * (V1 - V0);
	const Scalar F = F0 + f * (F1 - F0);