The Gay-Berne documentation and implementation do not match

[RainierBarrett]

Description

In the docs for the Gay-Berne anisotropic potential, there is a description of the model used to calculate the pairwise energy between two anisotropic particles that includes a $\zeta_{cut}$ term:

Two points of confusion arise here.

First, I think this $r$ is meant to be $r_{cut}$, not just $r$ -- I don't think the cutoff in $\zeta$ units depends on the center-to-center distance, right?

Second, this equation is not what is used in the code:

Scalar zetacut = rcut / sigma_max;

I don't know which of these two is the desired behavior, but from my understanding of the paper linked in the docs, it seems like the physics there expect $\zeta_{cut}$ to be calculated the way the docs show, rather than as rcut/sigma_max. Regardless, either the docs should be changed to reflect what math the code is doing, or the code should be adjusted to do the math the docs show.

Script

# not sure what to put here but this is the docstring
import hoomd
help(hoomd.md.pair.aniso.GayBerne)

Input files

No response

Output

Help on class GayBerne in module hoomd.md.pair.aniso:

class GayBerne(AnisotropicPair)
 |  GayBerne(nlist, default_r_cut=None, mode='none')
 |
 |  Gay-Berne anisotropic pair force.
 |
 |  Args:
 |      nlist (hoomd.md.nlist.NeighborList): Neighbor list
 |      default_r_cut (float): Default cutoff radius :math:`[\mathrm{length}]`.
 |      mode (str): energy shifting/smoothing mode.
 |
 |  `GayBerne` computes the Gay-Berne anisotropic pair force on every particle
 |  in the simulation state. This version of the Gay-Berne force supports
 |  identical pairs of uniaxial ellipsoids, with orientation-independent
 |  energy-well depth. The potential comes from the following paper `Allen et.
 |  al. 2006`_.
 |
 |  .. _Allen et. al. 2006: http://dx.doi.org/10.1080/00268970601075238
 |
 |  .. math::
 |      U(\vec r, \vec e_i, \vec e_j) =
 |      \begin{cases}
 |      4 \varepsilon \left[ \zeta^{-12} - \zeta^{-6} \right]
 |      & \zeta < \zeta_{\mathrm{cut}} \\
 |      0 & \zeta \ge \zeta_{\mathrm{cut}} \\
 |      \end{cases}
 |
 |  where
 |
 |  .. math::
 |
 |      \zeta &= \left(\frac{r-\sigma+\sigma_{\mathrm{min}}}
 |                         {\sigma_{\mathrm{min}}}\right),
 |
 |      \sigma^{-2} &= \frac{1}{2} \hat{\vec{r}}
 |          \cdot \vec{H^{-1}} \cdot \hat{\vec{r}},
 |
 |      \vec{H} &= 2 \ell_\perp^2 \vec{1}
 |          + (\ell_\parallel^2 - \ell_\perp^2)
 |            (\vec{e_i} \otimes \vec{e_i} + \vec{e_j} \otimes \vec{e_j}),
 |
 |  and :math:`\sigma_{\mathrm{min}} = 2 \min(\ell_\perp, \ell_\parallel)`.
 |
 |  The cut-off parameter :math:`r_{\mathrm{cut}}` is defined for two particles
 |  oriented parallel along the **long** axis, i.e.
 |  :math:`\zeta_{\mathrm{cut}} = \left(\frac{r-\sigma_{\mathrm{max}}
 |  + \sigma_{\mathrm{min}}}{\sigma_{\mathrm{min}}}\right)`
 |  where :math:`\sigma_{\mathrm{max}} = 2 \max(\ell_\perp, \ell_\parallel)` .
 |
 |  The quantities :math:`\ell_\parallel` and :math:`\ell_\perp` denote the
 |  semi-axis lengths parallel and perpendicular to particle orientation.
 |
 |  Example::
 |
 |      nl = nlist.Cell()
 |      gay_berne = md.pair.aniso.GayBerne(nlist=nl, default_r_cut=2.5)
 |      gay_berne.params[('A', 'A')] = dict(epsilon=1.0, lperp=0.45, lpar=0.5)
 |      gay_berne.r_cut[('A', 'B')] = 2 ** (1.0 / 6.0)
 |
 |  .. py:attribute:: params
 |
 |      The Gay-Berne potential parameters. The dictionary has the following
 |      keys:
 |
 |      * ``epsilon`` (`float`, **required**) - :math:`\varepsilon`
 |        :math:`[\mathrm{energy}]`
 |      * ``lperp`` (`float`, **required**) - :math:`\ell_\perp`
 |        :math:`[\mathrm{length}]`
 |      * ``lpar`` (`float`, **required**) -  :math:`\ell_\parallel`
 |        :math:`[\mathrm{length}]`
 |
 |      Type: `TypeParameter` [`tuple` [``particle_type``, ``particle_type``],
 |      `dict`]
 |
 |  Method resolution order:
 |      GayBerne
 |      AnisotropicPair
 |      hoomd.md.pair.pair.Pair
 |      hoomd.md.force.Force
 |      hoomd.operation.Compute
 |      hoomd.operation.Operation
 |      hoomd.operation.AutotunedObject
 |      hoomd.operation._HOOMDBaseObject
 |      hoomd.operation._HOOMDGetSetAttrBase
 |      hoomd.operation._DependencyRelation
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  __init__(self, nlist, default_r_cut=None, mode='none')
 |      Initialize self.  See help(type(self)) for accurate signature.
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties defined here:
 |
 |  type_shapes
 |      Get all the types of shapes in the current simulation.
 |
 |      Example:
 |
 |          >>> gay_berne.type_shapes
 |          [{'type': 'Ellipsoid', 'a': 1.0, 'b': 1.0, 'c': 1.5}]
 |
 |      Returns:
 |          A list of dictionaries, one for each particle type in the system.
 |
 |
 |      (`Loggable <hoomd.logging.Logger>`: category="object")
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from hoomd.md.pair.pair.Pair:
 |
 |  compute_energy(self, tags1, tags2)
 |      Compute the energy between two sets of particles.
 |
 |      Args:
 |          tags1 (``ndarray<int32>``): a numpy array of particle tags in the
 |              first group.
 |          tags2 (``ndarray<int32>``): a numpy array of particle tags in the
 |              second group.
 |
 |      .. math::
 |
 |          U = \sum_{i \in \mathrm{tags1}, j \in \mathrm{tags2}} V_{ij}(r)
 |
 |      where :math:`V_{ij}(r)` is the pairwise energy between two particles
 |      :math:`i` and :math:`j`.
 |
 |      Assumed properties of the sets *tags1* and *tags2* are:
 |
 |      - *tags1* and *tags2* are disjoint
 |      - all elements in *tags1* and *tags2* are unique
 |      - *tags1* and *tags2* are contiguous numpy arrays of dtype int32
 |
 |      None of these properties are validated.
 |
 |      Examples::
 |
 |          tags=numpy.linspace(0,N-1,1, dtype=numpy.int32)
 |          # computes the energy between even and odd particles
 |          U = mypair.compute_energy(tags1=numpy.array(tags[0:N:2]),
 |                                    tags2=numpy.array(tags[1:N:2]))
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from hoomd.md.force.Force:
 |
 |  additional_energy
 |      float: Additional energy term :math:`U_\mathrm{additional}`         :math:`[\mathrm{energy}]`.
 |
 |      (`Loggable <hoomd.logging.Logger>`: category="scalar")
 |
 |  additional_virial
 |      (1, 6) `numpy.ndarray` of ``float``: Additional virial tensor         term :math:`W_\mathrm{additional}` :math:`[\mathrm{energy}]`.
 |
 |      (`Loggable <hoomd.logging.Logger>`: category="sequence")
 |
 |  cpu_local_force_arrays
 |      hoomd.md.data.ForceLocalAccess: Expose force arrays on the CPU.
 |
 |      Provides direct access to the force, potential energy, torque, and
 |      virial data of the particles in the system on the cpu through a context
 |      manager. All data is MPI rank-local.
 |
 |      The `hoomd.md.data.ForceLocalAccess` object returned by this property
 |      has four arrays through which one can modify the force data:
 |
 |      Note:
 |          The local arrays are read only for built-in forces. Use `Custom` to
 |          implement custom forces.
 |
 |      Examples::
 |
 |          with self.cpu_local_force_arrays as arrays:
 |              arrays.force[:] = ...
 |              arrays.potential_energy[:] = ...
 |              arrays.torque[:] = ...
 |              arrays.virial[:] = ...
 |
 |  energies
 |      (*N_particles*, ) `numpy.ndarray` of ``float``: Energy         contribution :math:`U_i` from each particle :math:`[\mathrm{energy}]`.
 |
 |      Attention:
 |          In MPI parallel execution, the array is available on rank 0 only.
 |          `energies` is `None` on ranks >= 1.
 |
 |
 |      (`Loggable <hoomd.logging.Logger>`: category="particle")
 |
 |  energy
 |      float: The potential energy :math:`U` of the system from this force         :math:`[\mathrm{energy}]`.
 |
 |      (`Loggable <hoomd.logging.Logger>`: category="scalar")
 |
 |  forces
 |      (*N_particles*, 3) `numpy.ndarray` of ``float``: The         force :math:`\vec{F}_i` applied to each particle         :math:`[\mathrm{force}]`.
 |
 |      Attention:
 |          In MPI parallel execution, the array is available on rank 0 only.
 |          `forces` is `None` on ranks >= 1.
 |
 |
 |      (`Loggable <hoomd.logging.Logger>`: category="particle")
 |
 |  gpu_local_force_arrays
 |      hoomd.md.data.ForceLocalAccessGPU: Expose force arrays on the GPU.
 |
 |      Provides direct access to the force, potential energy, torque, and
 |      virial data of the particles in the system on the gpu through a context
 |      manager. All data is MPI rank-local.
 |
 |      The `hoomd.md.data.ForceLocalAccessGPU` object returned by this property
 |      has four arrays through which one can modify the force data:
 |
 |      Note:
 |          The local arrays are read only for built-in forces. Use `Custom` to
 |          implement custom forces.
 |
 |      Examples::
 |
 |          with self.gpu_local_force_arrays as arrays:
 |              arrays.force[:] = ...
 |              arrays.potential_energy[:] = ...
 |              arrays.torque[:] = ...
 |              arrays.virial[:] = ...
 |
 |      Note:
 |          GPU local force data is not available if the chosen device for the
 |          simulation is `hoomd.device.CPU`.
 |
 |  torques
 |      (*N_particles*, 3) `numpy.ndarray` of ``float``: The torque         :math:`\vec{\tau}_i` applied to each particle         :math:`[\mathrm{force} \cdot \mathrm{length}]`.
 |
 |      Attention:
 |          In MPI parallel execution, the array is available on rank 0 only.
 |          `torques` is `None` on ranks >= 1.
 |
 |
 |      (`Loggable <hoomd.logging.Logger>`: category="particle")
 |
 |  virials
 |      (*N_particles*, 6) `numpy.ndarray` of ``float``: Virial tensor         contribution :math:`W_i` from each particle :math:`[\mathrm{energy}]`.
 |
 |      Attention:
 |          To improve performance `Force` objects only compute virials when
 |          needed. When not computed, `virials` is `None`. Virials are computed
 |          on every step when using a `md.methods.NPT` or `md.methods.NPH`
 |          integrator, on steps where a writer is triggered (such as
 |          `write.GSD` which may log pressure or virials), or when
 |          `Simulation.always_compute_pressure` is `True`.
 |
 |      Attention:
 |          In MPI parallel execution, the array is available on rank 0 only.
 |          `virials` is `None` on ranks >= 1.
 |
 |
 |      (`Loggable <hoomd.logging.Logger>`: category="particle")
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from hoomd.operation.AutotunedObject:
 |
 |  tune_kernel_parameters(self)
 |      Start tuning kernel parameters.
 |
 |      After calling `tune_kernel_parameters`, `AutotunedObject` will vary the
 |      kernel parameters in subsequent time steps, check the run time of each,
 |      and lock to the fastest performing parameters after the scan is
 |      complete.
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from hoomd.operation.AutotunedObject:
 |
 |  is_tuning_complete
 |      bool: Check if kernel parameter tuning is complete.
 |
 |      ``True`` when tuning is complete and `kernel_parameters` has locked
 |      optimal parameters for all active kernels used by this object.
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from hoomd.operation.AutotunedObject:
 |
 |  kernel_parameters
 |      dict[str, tuple[float]]: Kernel parameters.
 |
 |      The dictionary maps GPU kernel names to tuples of integers that control
 |      how the kernel executes on the GPU. These values will change during the
 |      tuning process and remain static after tuning completes. Set the kernel
 |      parameters for one or more kernels to force specific values and stop
 |      tuning.
 |
 |      Warning:
 |          The keys and valid values in this dictionary depend on the hardware
 |          device, the HOOMD-blue version, and the value of class attributes.
 |          Keys and/or valid values may change even with new patch releases of
 |          HOOMD-blue.
 |
 |          Provided that you use the same HOOMD-blue binary on the same
 |          hardware and execute a script with the same parameters, you may save
 |          the tuned values from one run and load them in the next.
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from hoomd.operation._HOOMDBaseObject:
 |
 |  __getstate__(self)
 |      Helper for pickle.
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from hoomd.operation._HOOMDBaseObject:
 |
 |  loggables
 |      dict[str, str]: Name, category mapping of loggable quantities.
 |
 |  ----------------------------------------------------------------------
 |  Methods inherited from hoomd.operation._HOOMDGetSetAttrBase:
 |
 |  __dir__(self)
 |      Expose all attributes for dynamic querying in notebooks and IDEs.
 |
 |  __getattr__(self, attr)
 |
 |  __setattr__(self, attr, value)
 |      Implement setattr(self, name, value).
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from hoomd.operation._HOOMDGetSetAttrBase:
 |
 |  __dict__
 |      dictionary for instance variables (if defined)
 |
 |  __weakref__
 |      list of weak references to the object (if defined)

Expected output

See above.

Platform

CPU, Linux

Installation method

Conda package

HOOMD-blue version

3.11.0, but this discrepancy is true of 4.x as well

Python version

3.11.4

[joaander]

Feel free to submit a pull request that fixes this issue. Users of this potential, let us know which of the two proposed fixes is correct.

[github-actions[bot]]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs.

[github-actions[bot]]

This issue has been automatically closed because it has not had recent activity.