Add a function to equilibrate a compacted hydrogel

[jngrad]

Background

Diamond lattices are created by creating fully stretched polymer chains between crosslinker particles. It is sometimes desirable to compact the hydrogel to reduce pore sizes. There is more than one way to compact a polymer network, but not all of them are necessary computationally efficient, or safe. For example, progressively reducing the box size using the rescale method can lead to frequent retuning of the P3M parameters or Verlet skin, both of which being computationally demanding. Regarding safety, when using MPI parallelization and the regular decomposition cell system, the number of MPI ranks has an impact on the minimum box size[^domdec-formula], which can lead to cryptic error messages during the rescaling procedure.

There error messages are difficult to reproduce due to the stochastic nature of the Verlet skin tuning algorithm. They are also influenced by the number of MPI ranks in a non-linear manner[^integer-factorization]. It is important to set the skin to a small value before each call to system.change_volume_and_rescale_particles(), so that the skin remains under acceptable limits for the new box length.

ESPResSo has a sample that features a compression schedule that slowly rescales the box and relaxes the polymer network: minimal-diamond.py lines 98-104. The significance of this code is not properly conveyed to the user: the comment is not helpful and the docstring at the top of the file (that gets deployed to the user guide) is way too short to be useful[^sample-scripts-docs]. In its current stand, this code also runs the risk of disappearing if an ESPResSo maintainer stumbles into it and decides to simplify it or remove it, in an effort to reduce the complexity of the code samples.

Moving the polymer compression code to a dedicated library like pyMBE would make sense. There it could be encapsulated in a class and properly documented. It could also be extended with capabilities that would otherwise be out of scope for the ESPResSo project, e.g. adding a bending potential to the polymer chains. The ESPResSo sample would keep the current code snippet as a simplified version of that pyMBE feature, and the sample docstring would have an explanation of the compression code with a citation to pyMBE.

Proposed work

Introduce a feature in pyMBE to run a compression schedule on a polymer network. The compression schedule would shrink the box volume and run the integrator for a few time steps in a loop until the desired box size is obtained. The integrator would thermalize the polymer to absorb the kinetic energy released when the compressed bonds return to their equilibrium distance. The system.time variable could be reset at the end, since the compression is a form of energy minimization. The function or class could take an ESPResSo system and a diamond lattice object as arguments, as well as parameters relevant to the compression schedule, such as a number of integration steps per compression step, the compression factor in units of metres per second, or a class representing an external potential to induce curvature of the polymer chains[^collagen].

Improving safety

The feature could guarantee the correctness of the simulation by working out the optimal tuning parameters. In particular, the regular decomposition must have 2 cells per direction, and the maximal cell width must be greater than the maximal cutoff plus the Verlet list skin $s$. In other words, the maximal Verlet list skin can be derived as:

$s^{\mathrm{max}} = \min\left(\dfrac{\vec{l}}{\vec{n} \cdot \vec{c}}\right) - r^{\mathrm{cut}}$

where:

• $\vec{l} =$ system.box_ l
• $\vec{n} =$ system.cell_system.node_grid
• $\vec{c} =$ system.cell_system.get_state()["cell_grid"]
• $r^{\mathrm{cut}} =\sigma^{\mathrm{WCA}} \cdot 2^{1/6}$ or max(max(system.cell_system.max_cut_bonded, system.cell_system.max_cut_nonbonded), 0)

The values of $\vec{n}$ and $\vec{c}$ are obtained by an iterative solver in ESPResSo, such that each element $x_i$ in the Hadamard product $\vec{n} \odot \vec{c}$ satisfies the conditions $2 \le x_i \le 32$ and $l_i \cdot x_i^{-1} \ge (s + r^{\mathrm{cut}})$. The value of $s^{\mathrm{max}}$ can be passed as argument max_skin of system.cell_system.tune_skin() with adjust_max_skin=False.

For illustration purposes, consider the following simulation script:

import espressomd
import pprint
system = espressomd.System(box_l=[10, 10, 10])
system.non_bonded_inter[0, 0].lennard_jones.set_params(epsilon=1., sigma=1., cutoff=3., shift="auto")
system.cell_system.skin = 2.
pprint.pprint(system.cell_system.get_state())

With 2 MPI ranks:

• node_grid = [2, 1, 1] such that MPI rank 0 has the slab from x=0 to x=5 and rank 1 has the slab from x=5 to x=10
• cell_grid = [1, 2, 2] such that there is one cell per slab along the x-axis and two cells along the y- and z-axes
• cell_size = [5, 5, 5] such that each cell is a cube (this is obtained by dividing the box length by the product of the node grid and cell grid)
• max_cut_nonbonded = 3 and skin = 2 such that the interaction range is 5

Setting a skin of 2.1 would fail, because the interaction range 5.1 would exceed the cell size of 5. Here is the raised exception:

Traceback (most recent call last):
  File "script.py", line 5, in <module>
    system.cell_system.skin = 2.1
  File "script_interface.pyx", line 464, in espressomd.script_interface.ScriptInterfaceHelper.__setattr__
  File "script_interface.pyx", line 196, in espressomd.script_interface.PScriptInterface.set_params
  File "utils.pyx", line 213, in espressomd.utils.handle_errors
Exception: while setting parameter 'skin': ERROR: interaction range 5.1 in direction 0 is larger than the local box size 5

One can visualize the cell system partitioning with mpiexec -n 2 ./pypresso ../samples/visualization_cellsystem.py. The visualizer shows the local cells of MPI rank 0 as a wireframe and colors particles according to the MPI rank they reside in.

Improving performance

Since the compression algorithm takes the form of an energy minimization, the accuracy of the Hamiltonian can be decreased to increase performance. For example, charged hydrogels can use the simple Debye-Hückel potential to capture the effects of the electrostatic interaction on short distances. The P3M method is more accurate but also very computationally demanding in simulations where the box length changes frequently, since each change automatically triggers a retuning of the parameters to meet the requested accuracy of the method. In addition, the Verlet list skin probably doesn't need to be retuned at every compression step, instead one could simply rescale it and clamp it to the new $s^{\mathrm{max}}$ value for the new box length, and do a proper retune every 20 compression steps.

[^domdec-formula]: ESPResSo 4.2.1 user guide, section 4.2.2.2. Regular decomposition [^integer-factorization]: The number of cells in regular decomposition is constrained by the MPI Cartesian topology, which depends on the integer factorization of the number of MPI ranks. Some values are susceptible to yielding large prime factors, e.g. 20 MPI ranks will be decomposed into [5, 2, 2] domains. With 20 MPI ranks and a cubic box, the local domains will be slabs which are thinner along the x-direction by a factor 2.5 compared to the y- and z-directions. [^sample-scripts-docs]: All sample scripts in ESPResSo start with a docstring that supports reSt syntax and gets deployed to the user guide section 1.5. Sample scripts as a definition list. The reSt syntax offers many advanced features, such as citations with :cite:`frenkel02b` or math mode with :math:`\nu`. [^collagen]: Biological tissues like collagen are structured as polymer networks with curved chains to confer a high degree of elasticity. This can only be modeled with a compressed hydrogel, otherwise the polymer chains would naturally stretch back into their original linear topology to minimize their energy.

[kosovan]

When relaxing the hydrogel from the initial stretched conformation to the desired box size, electrostatics should be turned off completely. It should be activated only after the desired box size is reached. This improves the performance at no cost.

[kosovan]

@jngrad I was not aware that such a sample code exists but I would definitely like to review the updated version.

[pm-blanco]

@kosovan to put you in context: since we have the plan to extend pyMBE to support the setup of hydrogels we were discussing with @jngrad that it could be interesting to have a helper function in pyMBE to aid the equilibration of hydrogels