Add CMAP support

[JAMelendezD]

Implement CMAP energy, force, and virial calculations using B-spline interpolation. Add unit tests and update the example to use ubiquitin. Uses charmm style cmaps based on atomtypes and a 24x24 grid.

[jeg7]

Required changes

1. Preserve the selected calculation precision in the CMAP kernel

• In src/CudaBondedForce.cu, replace the precision-narrowing operations at the anchors below:

  const int ix = (int)floorf((float)x);
  const int iy = (int)floorf((float)y);

  and:

  const CT rg1 = sqrtf(...);
  const CT rabr1 = sqrtf(...);
  const CT phi1 = atan2f(st1, ct1) * ...;

  The current implementation narrows the CT=double specialization to single precision.

  Please use type-preserving floor and atan2 operations for both float and double. There is an existing sqrt_template<CT> mechanism that can be reused for square roots. In cmap_pot the bicubic coefficients should be cast to type CT.

2. Replace the one-vector-per-coefficient CMAP representation with contiguous storage

• In src/CharmmParameters.cu, replace the current packing at:

  for (const double coefficient : value.coeff)
    paramsVal.push_back({static_cast<float>(coefficient)});

  One 24 x 24 map contains 9,216 coefficients, so this representation creates 9,216 separate one-element std::vector<float> objects per unique CMAP type. That nested structure is then copied and flattened again in CudaBondedForce::setup_coef.

  Would it be a huge rewrite to store CMAP coefficients in a contiguous std::vector<float> payload in BondedParamsAndLists or another appropriate packed-data structure, and copy that buffer directly to device memory?

  The updated layout should make the following quantities explicit:

  • Number of unique CMAP parameter maps.
  • Number of coefficients per map: 24 * 24 * 16.
  • Total number of coefficient scalars.
  • Mapping from cmaplist_t::itype to the corresponding contiguous coefficient range.

3. Resolve the CMAP key orientation

• include/CharmmParameters.h describes CMAP keys as containing canonical dihedrals, but src/CharmmParameters.cu preserves the exact parameter-file orientation at:

  const DihedralKey dihe1(tokens[0], tokens[1], tokens[2], tokens[3]);
  const DihedralKey dihe2(tokens[4], tokens[5], tokens[6], tokens[7]);

  The PSF packing path also performs exact-order lookup without canonicalization.

  Choose and document one:

  1. CMAP keys are ordered and orientation-sensitive; or
  2. Reversed or swapped torsions are canonicalized, with the corresponding grid transformation or transposition applied correctly.

  Do not independently canonicalize the torsion keys without transforming the two-dimensional map.

• Add tests using an asymmetric synthetic CMAP grid for:

  • Exact parameter-file orientation.
  • Reversal of the first torsion.
  • Reversal of the second torsion.
  • Swapping the two torsions.

  An asymmetric grid is important because a symmetric production map may hide an accidental axis transposition.

4. Improve CMAP parser diagnostics and validation

• In src/CharmmParameters.cu, the context created at:

  const std::string valueContext =
      GetPrmValueContext("CMAP", line, fileName, lineNumber);

  is based on the header and then reused for every grid value. A malformed value on a later grid line is therefore reported as occurring on the header line.

  Construct the diagnostic context from the current cleaned grid line and current lineNumber before parsing that line's tokens.

Performance validation

5. Remove repeated spline work during parameter initialization

• In src/CmapSpline.cpp, the temporary psi interpolation and the two phi-direction spline solves are recomputed inside the phi loop even though they depend only on the current psi value.

  Can you make the psi loop outermost, construct and spline the temporary arrays once per psi value, and then evaluate all phi grid points from those precomputed spline arrays?

Style and documentation cleanup

• Update the stale include/CharmmParameters.h documentation that still states:

  CMAP output is always empty.
  

• Correct the comment in src/CudaBondedForce.cu that says each CMAP map has 576 consecutive coefficient values. The kernel indexes 16 coefficients for each of 576 cells, so the actual coefficient count is 9,216 floats per map?

• Add the omitted CMAP: row to CudaBondedForce::print().

• Update the CudaBondedForce API documentation to include the calc_cmap argument and CMAP energy component.

• Run clang-format, which should remove trailing whitespace, and add final newlines to the new spline source and header.