(1) Overview
Introduction
The computational study of polymeric systems has become increasingly important, driven by advances in computational hardware that enable simulation of increasingly large and architecturally complex polymeric systems. Such simulations can reproduce experimental data and show promise for the architectural design of polymer networks, potentially accelerating research and development in the polymer industry.
When simulating a polymer system, a typical workflow involves:
setting up an initial structure through a Monte Carlo (MC) procedure [1, 2], a hierarchical strategy [3, 4], or a molecular dynamics (MD) diffusion-collision-based generator [5, 6, 7],
running simulations using MC, dissipative particle dynamics (DPD), or MD methods,
analyzing and visualizing the results.
The Kremer-Grest (KG) bead-spring model [8, 9, 10] provides an optimal balance between structural detail and computational efficiency for studying the mechanical properties of polymers. The computer codes for KG simulations are usually highly optimized to support parallel computations on both CPUs and GPUs in order to maximize the productivity of the researchers while allowing for the study of larger systems and longer trajectories, as is required in particular for polymeric systems.
Usually, general MD simulators or programs for atomistic simulations are used for such simulations of polymeric systems. Examples include Large-scale Atomic/Molecular Massively Parallel Simulator [11] (LAMMPS), GROningen MAchine for Chemical Simulations [12] (GROMACS), COarse-Grained molecular dynamics program by NAgoya Cooperation [13] (COGNAC) (part of the Open, flexible and expandable system OCTA [14] (OCTA) project), or CP (Car-Parrinello = ab initio MD) code for the new Millennium; note that CP2K [15] didn’t implement the CP method (CP2K). As a consequence of the programs being general, even though these programs provide many possibilities to output various properties of the simulated system, they are not specific for polymeric systems. Therefore, various network properties such as p (extent of reaction), pgel (gelation point), or r (stoichiometric imbalance) have to be computed externally, through clever use of the commands provided by the program, or by providing a custom addition to the program. This becomes even more pronounced for analyzing the network structure, such as determining ΦD (mass fraction of dangling chains), Φel (mass fraction of network backbone), or the distribution and length of loops in the structure, as these are properties that are key only for polymeric simulations and therefore have no immediate applications in other MD projects, such as simulating metals, crystals, glasses, or galaxies.
Some specialized tools exist for specific aspects of polymer simulation. For example, the Pizza.Py toolkit offers a Python [16] programming interface for LAMMPS, which can be used to create and manipulate structures. The MDAnalysis Python package [17, 18] also provides a toolkit to read and process MD data from various simulation programs. For generating worm-like three-dimensional chains, the PolymerCpp Python package [19] is available. A hierarchy of polymer models can be generated using the Polymer Structure Predictor [20] (PSP). And the toolbox by Barrett [21] provides polymer generation functionality based on molecular templates, and methods to compute the static structure factor and the radial distribution function.
However, these tools often lack interoperability and comprehensive coverage of polymer-specific analysis methods. To address these limitations, we have developed pylimer-tools, a comprehensive Python package that consolidates polymer-specific computational tools into a unified, efficient framework designed specifically for bead-spring polymer network (see Figure 1) studies.
Figure 1
2D illustration of a bead-spring network with characteristic defects highlighted: in cyan dangling structures and chains, in blue free structures and chains, in pink a secondary loop and in purple a primary loop. The cross-links are shown enlarged with a black border. The distribution of functionalities, end-to-end distances and defects are for illustration purposes only and do not necessarily represent a realistic polymer network [22].
Implementation and architecture
This package, pylimer-tools, is implemented primarily in C++ for optimal performance of computationally intensive calculations and input/output operations. Python bindings are provided through pybind11 [23], ensuring seamless integration with Python workflows while maintaining computational efficiency. The package is cross-platform compatible, tested on Linux, macOS, and Windows with Python 3.9 and later versions. Comprehensive unit tests and continuous integration ensure code quality and correctness across different platforms and use cases.
Core Architecture
The software leverages the igraph library [24] to represent molecular structures as graphs, where atoms are vertices and bonds are edges. This graph-based representation enables efficient analysis of network topology, including clustering, connectivity traversal, angle and dihedral detection, and manipulation of periodic boundary conditions. The graph structure facilitates rapid determination of atom coordination numbers and the decomposition of systems into disconnected networks or individual polymer chains.
Data Structures
The central Universe class encapsulates complete polymer network structures. Users can build these programmatically by adding atoms and bonds, or import existing structures from LAMMPS data and dump files. A trajectory from a dump file or from a sequence of data files is represented by the UniverseSequence class, which allows memory-efficient on-demand reading of the respective Universe at a specific time-step. The Universe can be decomposed into Molecules, which represent the strands and chains in the structure. This is done by either splitting the graph at junction points as determined by their functionality, or by manually specifying an atom type to use as junctions and splitting at those. After manipulating the structure, this library also provides the functionality to output the structure back into a data file compatible with LAMMPS. The package provides additional comprehensive file I/O capabilities, including reading various simulation output files of LAMMPS (correlation data, vector outputs, thermodynamic logs).
Network Generation
The MCUniverseGenerator class implements flexible Monte Carlo network generation algorithms [1, 25] supporting end-linked and vulcanized networks, as well as complex architectures like bottle-brush or comb-like polymers [26]. The beads between the cross-links are placed by default using a Brownian Bridge process [27, 28, 29, 30, 31], though additional MC steps can be taken to progress with the equilibration of the Gaussian springs between the beads.
Force Balance and Energy Minimization
A key feature is the implementation of the Maximum Entropy Homogenization Procedure (MEHP) [1, 2] and the related Force Balance procedure. There are three different implementations available. In the class MEHPForceRelaxation, phantom network shear moduli can be determined efficiently with alternative, nonlinear spring potentials. In the class MEHPForceBalance, the package extends the phantom network with slip-links (see Figure 2) to model entanglements, providing estimates of real network moduli. The third implementation is given in MEHPForceBalance2 through the Force Balance procedure with fixed entanglement links, that is, entanglements modeled as tetrafunctional cross-links, as published recently [32]. This is our fastest implementation of the MEHP, but is consequently less flexible than the other two implementations. Apart from making accurate predictions of the modulus of real networks, it also enables a systematic investigation of different entanglement detection methods and their effects on mechanical properties, for example.
Figure 2
Illustration of a slip-link in orange. The other colors serve as a reference of the surrounding polymer beads.
Simulation Capabilities
An integrated DPD simulator with slip-spring (see Figure 3) entanglement modeling [33, 34] supports multi-core parallel execution via OpenMP. Normal mode analysis for dynamic mechanical properties (stress autocorrelation function, loss and storage moduli [35]) is implemented using LAPACK [36] for efficient eigenvalue computations [35].
Figure 3
Illustration of a slip-spring in orange. The other colors serve as a reference of the surrounding polymer beads.
Advanced Algorithm Implementations
An efficient loop detection algorithm based on Johnson [37] was implemented and can be used to detect loops of arbitrary order in the network. The package includes a neighbor list implementation using spatial partitioning that optimizes performance for entanglement sampling and force calculations.
Periodic Boundary Conditions
All the simulation methods mentioned so far allow switching between two modes of periodic boundary conditions (PBC). The classical mode implements the periodic boundaries simply by using the shortest distance, as is also implemented, for example, in LAMMPS. However, there are some scenarios where this is not sufficient. For example in the MEHP, to reduce the number of degrees of freedom, multiple bonds are combined into one long spring (bifunctional beads are not represented in the degrees of freedom). Another example is the DPD simulations, which may lead to extraordinarily long bonds due to random fluctuations and the soft bond potential. For these scenarios, one has to prevent atoms or beads from crossing into another periodic image (see Figure 4). This is achieved by deriving a periodic offset at the beginning of the procedure. At that point, the springs have not yet replaced the bonds, and the bonds are assumed to still be shorter than half the box length. Consequently, one can follow the graph structure to find the bonds that must have an offset to account for the periodic boundary condition. Throughout the remainder of the procedure, this offset will remain associated with the corresponding bonds or springs. Apart from maintaining this constant offset, no additional steps are necessary to achieve a functional implementation of periodic boundary conditions.
Figure 4
An illustration of the dangers of the common periodic boundary conditions for optimization procedures implemented in pylimer-tools: when strands are converted into single springs, if the strands are longer than half the box size, the springs might collapse (step 3) even if they should not have. To prevent this, the initial state is used to derive multiples of the box to be used as offsets when computing the spring lengths. In this example, the vector r1 from A to B would have an offset of negative one box length.
Unit Management
Integration with the pint library [38] provides automatic unit conversion between LAMMPS unit systems [39] and système international (d’unités) (SI) units, using conversion factors from Everaers et al. [40], preventing unit-related errors in calculations.
Theoretical Framework Implementation
Readily available methods are provided to use the Miller-Macosko Theory (MMT) [41, 42] for relating the quantities p, r, equilibrium shear modulus Geq and f of end-linked polymer networks, which were shown to be highly accurate [6]. We have implemented the methods for predicting the wsol, wdang, pgel, and both the phantom and the entangled equilibrium shear modulus. For all of them, we go further than the previously published literature: for f > 4, we solve the polynomial equation to predict the probability numerically. We include the possibility of having a fraction of monofunctional chains for arbitrary f to predict any of the relevant properties. These derivations were validated using those published previously for the determination of wsol for f = 3 [43] and f = 4 [44].
Quality control
Software quality is ensured through a comprehensive testing framework using Catch2 for C++ components and pytest for Python code. The test suite includes unit tests, integration tests, and performance benchmarks. Continuous integration through GitHub Actions automatically runs a representative set of tests on every code change, ensuring that new contributions do not introduce regressions or break existing functionality. The complete test suite is run on multiple platforms (Linux, macOS, Windows) and Python versions (3.9–3.12) released before a new version of pylimer-tools is published on PyPI.
The testing pipeline includes the following.
Unit tests: Individual function and class testing with edge case validation
Integration tests: End-to-end workflow testing including file I/O operations
Numerical validation: Comparison with analytical solutions and published results
Performance benchmarks: Monitoring computational efficiency for critical algorithms
Memory safety: Leak detection and bounds check using sanitizers
Version releases follow semantic versioning with automated deployment to PyPI after passing all tests. The package includes extensive documentation with API references, tutorials, and example workflows. Code quality is maintained through automated formatting (clang-format, ruff) and static analysis.
For end users, the software provides clear error messages and warnings for invalid inputs or edge cases. Example datasets and benchmark problems are included to verify the correct installation and functionality. The test suite can be run locally using the build scripts provided, enabling users to validate their installation environment.
(2) Availability
Operating system
This package is at least compatible with any MacOS, Windows and Linux version supporting Python 3.9 and newer.
Programming language
The core of the package is implemented in C++20 and exposed to Python 3.9+ through pybind11 bindings.
Additional system requirements
The package installation requires approximately 10 MB of disk space. Memory and CPU requirements scale with system size: typical polymer networks with around 5 × 105 atoms require 1 GB to 8 GB RAM for a workflow that includes network generation using the MC procedure, entanglement sampling and the Force Balance procedure to determine the equilibrium shear modulus. Multi-core processors are recommended for DPD simulations and normal mode analysis.
Dependencies
Python packages that are required for basic functionality:
pandas[45, 46]: reading simulation output files will return Pandas data frames for comfortable data handling
numpy[47]: for efficient vectorized calculations
pint[38]: for unit conversions
Build dependencies (required only for source compilation):
CMake: Cross-platform build system
C++20-compatible compiler
Python development headers (python3-dev on Linux)
External libraries (automatically downloaded during build if not found on the system):
List of contributors
Tim Bernhard, ETH Zürich: main developer. These codes are the means used for all results published for and during his PhD
Andrei A. Gusev, ETH Zürich: initial Mathematica and C versions of the MC network generator and the MEHP force relaxation codes that had a large influence on the current C++ codes
Fabian Schwarz, Uppsala University: initial Matlab versions of the MC network generator and the MEHP force relaxation codes that had a large influence on the current C++ codes
Jorge Ramírez, Universidad Politécnica de Madrid: multiple-tau autocorrelation [51] C codes used as a basis of the C++ implementation
Martin Kröger, ETH Zürich: discussions regarding implementations of periodic boundary conditions for small systems
Software location
Archive
Name: Zenodo
Persistent identifier: https://doi.org/10.5281/zenodo.17582031
Licence: GNU General Public License v3.0 or later
Publisher: Tim Bernhard
Version published: 0.3.13
Date published: 11/11/2025
Code repository
Name: GitHub
Persistent identifier: https://github.com/GenieTim/pylimer-tools
Licence: GNU General Public License v3.0 or later
Date published: 31/07/2025
Language
English
(3) Reuse potential
Other researchers can use this package already if their sole intent is to read LAMMPS output files, since it provides fast and ergonomic functions to do so, or if they use the MMT thanks to the respective formulas being implemented and easily available. The real benefit becomes apparent for polymer researchers who do KG MD simulations or other studies using bead-spring polymer models. For these researchers, this package provides a large number of tools that can be used, from generating structures using MC procedures, writing these files for usage with LAMMPS, using the structures with the built-in DPD simulator or analyzing the loop distributions, end-to-end distances, bond lengths, or reducing it to its minimum energy state using the Force Balance procedure to compute the equilibrium shear modulus.
As the software is licensed under an open source license, further adaptations are also possible. Feature requests and problems or issues can be submitted by raising an issue in the GitHub repository. The distribution includes a detailed manual that covers the installation and use of the software.
Notation
| NOTATION | DESCRIPTION | SYMBOL |
|---|---|---|
| ΦD | mass fraction of dangling chains | ΦD |
| Φel | mass fraction of network backbone | Φel |
| f | cross-link functionality | f |
| pgel | gelation point | pgel |
| p | extent of reaction | p |
| r | stoichiometric imbalance | r |
| wdang | dangling chain fraction | wdang |
| wsol | soluble chain fraction | wsol |
| COGNAC | COarse-Grained molecular dynamics program by NAgoya Cooperation [13] | |
| CP2K | CP (Car-Parrinello = ab initio MD) code for the new Millennium; note that CP2K [15] didn’t implement the CP method | |
| DPD | dissipative particle dynamics | |
| equilibrium shear modulus | stress relaxation modulus for a viscoelastic solid [52] | Geq |
| GROMACS | GROningen MAchine for Chemical Simulations [12] | |
| KG | Kremer-Grest | |
| LAMMPS | Large-scale Atomic/Molecular Massively Parallel Simulator [11] | |
| MC | Monte Carlo | |
| MD | molecular dynamics | |
| MEHP | Maximum Entropy Homogenization Procedure | |
| MMT | Miller-Macosko Theory | |
| OCTA | Open, flexible and expandable system OCTA [14] | |
| PBC | periodic boundary conditions | |
| PSP | Polymer Structure Predictor [20] | |
| SI | système international (d’unités) |
Acknowledgements
The authors thank Martin Kröger for the discussion on alternative implementations of PBC.TB also wishes to thank Szabolcs Horvát for his positive interactions on the use and extension of igraph, and Jorge Ramírez for providing access to his implementation of the multiple-tau autocorrelation and answering questions about its license.
Competing Interests
The authors have no competing interests to declare.