Draft: Modernize tabulated potentials

Gregory Lazarski (Golui)requested to merge
golui-tabulated-potential/calculate into main

Refs #1347

This merge request re-introduces non-bonded tabualted potentials to GROMACS (currently CPU only, see "Caveats" for details).

Changes presented in !5194 !5439 are included in this MR.

Highlights

  • Table definitions, for both non-bonded and bonded interactions, in topology, as either lists of values or muParser formulae
  • Tables stored in the TPR, including preserving formulas
  • Includes tests

User interface

The format follows concepts introduced in !5194. Since then, the syntax has evolved slightly.

The most significant differences are:

  • Two new entries in the MDP: user-tables and user-table-tolerances. Both are lists, with one being for table names, and the other for the tolerances.
    • A special * identifier specified in user-tables allows overriding default values (e.g. tolerance) for all. Values for tables specified explicitly still take precedence, regardless of the order the names appear in the list. See example below.
  • The use of different variable names, depending on which table type is being defined
  • Introduction of other variables and constants into the parser context (e. g. VdW cutoff as r_cut and switchign radius as r_sw for non-bonded tables)
  • introduction of two special keywords that can precede table data for both the numerical and analytic description: limits and accuracy-range
  • Note that for non-bonded tables, the second (dispersive) term is specified without absorbing the preceding negative sign. This convention is persists until the kernel combines the potential/force results using the fms operation. In addition, derivatives are specified directly as well, without absorbing any constants.

The concept of non-bonded groups for atoms and the non-bonded table to group assignment remain the same as in the previous MRs.

Example of a table tolerance specification:

; Specify tolerances for tables "Bonded1" and "Default_NB" as 1e-4 and 1e-6, respectively, and 5e-2 for all other tables.
user-tables = Bonded1 * Default_NB
user-table-tolerances = 1e-4 5e-2 1e-6

limits is only required for atom-atom bond interactions (not all bonded interactions) defined with a formula. This is used to set the interval over which to generate the table data. Example:

limits 0 10

accuracy-range can be used for all tables, but finds most use for tables with singularities at 0, such as non-bonded tables. This limits the region that will be checked by GROMACS for accuracy and consistency. In most cases, the values are automatically inferred. For non-bonded tables, GROMACS will look for the smallest sigma defined in the force field, and set the lower limit to minSigma / 10, and the upper limit to vdwCutoff, unless explicitly overriden by the user. Example:

accuracy-range 0.005 1.2

Full example of tables defined with a formula:

[ table ]
LJ		1

[ tableexpressions ]
; Define the range where validation should be performed (optional)
accuracy-range 0.005 1.2
; Repulsive 
r_ij^-12 - r_cut^-12
; Dispersive
r_ij^-6 - r_cut^-6
; Repulsive derivative (optional)
-12 * r_ij^-13
; Dispersive derivative (optional)
-6 * r_ij^-7
[ table ]
; Example table using more complex pre-defined functions to easily implement a CHARMM-like force-switched potential.
LJForceSwitch		1

[ tableexpressions ]
; Define the range where validation should be performed (optional)
accuracy-range 0.005 1.2
; Repulsive 
force_switch_potential(r_ij, 12, r_sw, r_cut)
; Dispersive
force_switch_potential(r_ij, 6, r_sw, r_cut)
; Repulsive derivative (optional)
force_switch_force(r_ij, 12, r_sw, r_cut)
; Dispersive derivative (optional)
force_switch_force(r_ij, 6, r_sw, r_cut)

And examples of tables defined as a list of values:

[ table ]
LJForceSwitchList		1

[ tablevalues ]
; Define the range where validation should be performed (optional)
accuracy-range 0.01 1.2
; Data can be included directly, for example:
; #include "lj_explicit.tsv"
; Below are the first few lines of the file:
;   r_ij	           V_rep	          V_disp	           F_rep	          F_disp
 0.00000	 0.000000000e+00	 0.000000000e+00	 0.000000000e+00	 0.000000000e+00
 0.00050	 4.096000000e+39	 6.400000000e+19	-9.830400000e+43	-7.680000000e+23
 0.00100	 1.000000000e+36	 1.000000000e+18	-1.200000000e+40	-6.000000000e+21
 0.00150	 7.707346629e+33	 8.779149520e+16	-6.165877303e+37	-3.511659808e+20
 0.00200	 2.441406250e+32	 1.562500000e+16	-1.464843750e+36	-4.687500000e+19
; ...

Example of a bonded table:

[ table ]
; Defines a bonded interaction that acts as the repulsive component of a Mie 9-6 potential, for sigma = 0.5nm
; The function is shifted such that it reaches 0 at the minimum, and is 0 beyond the minimum.
RepulsiveBondedSigma0.5     3

[ tableexpressions ]
limits 0 5
accuracy-range 0.06 0.66
((0.5 / rb_ij)^9 - (0.5 / rb_ij)^6 + 4 / 27) * (rb_ij < 0.5 * (3/2)^(1/3))
(-9 * (0.5 / rb_ij)^9 / rb_ij + 6 * (0.5 / rb_ij)^6 / rb_ij) * (rb_ij < 0.5 * (3/2)^(1/3))

Implementation details

Non-bonded kernels are adapted from the code by @berkhess !5171.

Bonded kernels use the same code for interaction computation as previously. I have introduced a new code path that reads bonded tables specified in the topology to generate bondedtable_t objects. This code was made in a rush and duplicats a lot of the logic from *SplineTable classes, and therefore will need a cleanup. Old table reading code is retained - meaning that tables for bonded interactions can still be specified at mdrun invocation time as previously, with no changes.

Spline calculation uses the existing QuadraticSplineTable and CubicSplineTable code, with a few small tweaks:

  • Construction has been streamlined by using std::variant
  • An additional logical "layer" has been added. Previously, a single *SplineTable could hold only a single table that contained a fixed number of terms (with optimisations up to 3 terms). Changes introduced in 2ca0ed41 add support for storing multiple functions in a single object, which allows for the computation of multiple tables indexed by an SimdInt32. This is at a cost of an integer multiplication. Another drawback of this approach is that every function inside a *SplineTable needs to have the same spacing. This means that the minimum spacing for of all the defined functions will be used to ensure a uniform size of the internal array. This can lead to performance degradation and wasted memory. This can be observed when running the tabulated kernel tests in double precision. This necessitated the implementation of a new SIMD primitive for gathered load from arbitrary locations in memory.

Caveats

  • Consistency checks between the function and its derivative, and the automatic spacing calculations, seem to be overly sensitive. This is especially apparent in functions defined numerically.
    • There is a slight asymmetry between analytical and numeric definitions. Since analytical definitions do not define any spacing a priori, computation of the minimum quotient of the function and one of its derivative (for spacing inference) happens with an arbitrarily defined granularity of 500 points. This is in contrast with the same procedure for numerical tables, where every step is checked. This leads denser spacing requirements for such tables.
  • Currently, only CPU kernels are supported.
  • These changes are known to not compile with AVX_128_FMA and AVX_256, due to missing Simd*Int32 types or other functions. In addition, currently the IBM VSX SIMD and ARM SVE architectures is not supported, as I don have hardware that supports the ISA. Other architectures have been tested, with AVX_512, AVX2_256, and ARM NEON having been tested the most.
    • This is because of the need for the new gatherLoad (and its Hsimd variant) primitive. This is used for indexing tables in the NB kernels.
  • Due to how bonded interaction definitions are parsed, tables are still specified with numerical indices. This means that the order of the table definitions in the itp matters. This also means that older topologies that use tabulated bonded potentials are compatible with these changes.
  • While all tests pass, I feel that tabulated angles and dihedrals have not been tested extensivel enough yet.

Example usage

Full example systems for comparison between analytical runs are included here: CHARMMWater.tar.gz.

Remaining work

  • Documentation (code)
  • Documentation (manual, partially complete)
  • Support for AVX_128_FMA and AVX_256 OR refactoring code to not require the gatherLoad primitive
  • GPU implementation

... and likely more after feedback.

Supersedes: !5194 !5439

Merge request reports