A Computational Model for Molecular Interactions Between Curved Slender Fibers Undergoing Large 3D Deformations With a Focus on Electrostatic, van der Waals and Repulsive Steric Forces

This contribution proposes the first 3D beam-to-beam interaction model for molecular interactions between curved slender fibers undergoing large deformations. While the general model is not restricted to a specific beam formulation, in the present work it is combined with the geometrically exact beam theory and discretized via the finite element method. A direct evaluation of the total interaction potential for general 3D bodies requires the integration of contributions from molecule or charge distributions over the volumes of the interaction partners, leading to a 6D integral (two nested 3D integrals) that has to be solved numerically. Here, we propose a novel strategy to formulate reduced section-to-section interaction laws for the resultant interaction potential between a pair of cross-sections of two slender fibers such that only two 1D integrals along the fibers' length directions have to be solved numerically. This section-to-section interaction potential (SSIP) approach yields a significant gain in efficiency, which is essential to enable the simulation of relevant time and length scales for many practical applications. In a first step, the generic structure of SSIP laws, which is suitable for the most general interaction scenario (e.g. fibers with arbitrary cross-section shape and inhomogeneous atomic/charge density within the cross-section) is presented. Assuming circular, homogeneous cross-sections, in a next step, specific analytical expressions for SSIP laws describing short-range volume interactions (e.g. van der Waals or steric interactions) and long-range surface interactions (e.g. Coulomb interactions) are proposed. The validity of the SSIP laws as well as the accuracy and robustness of the general SSIP approach to beam-to-beam interactions is thoroughly verified by means of a set of numerical examples considering steric repulsion, electrostatic or van der Waals adhesion.

also for the interaction of a carbon nanotube with a Lennard-Jones (LJ) wall in 3D. 29 In both cases, the influence of the rigid half space can be evaluated analytically and formulated as a distributed load on the beam.
To the best of the authors' knowledge, no approach for describing molecular interactions between curved 3D beams for arbitrary configurations and large deformations has been proposed yet. Notable previous approaches to similar problems have made simplifying assumptions. Ahmadi and Menon 30 study the clumping of fibers due to vdW adhesion by means of an analytical 2D beam method, yet only include vdW interaction between the hemispherical tips based on an analytical expression for the interaction of two spheres. A numerical study of the influence of interfiber adhesion on the mechanical behavior of 2D fiber networks assumes an effective adhesion energy per unit length of perfectly parallel fiber segments and solves for the unknown contact length in a second, nested minimization algorithm. 31 In this article, we propose the first model specifically for molecular interactions between arbitrarily curved and oriented slender fibers undergoing large deformations in 3D. While the general model is not restricted to a specific beam formulation, in the present work, it is combined with the geometrically exact beam theory and discretized via the finite element method. This novel approach is based on reduced SSIP laws that describe the resulting interaction potential between a pair of cross-sections as a closed-form analytical expression. Thus, the two-body interaction potential follows from two nested 1D integrals over this SSIP law along both fibers' axes, which are evaluated numerically. In this way, the proposed, so-called SSIP approach significantly reduces the dimensionality of the required numerical integration from six to two, and hence the associated computational cost. As compared to methods for 3D solid bodies or even to all-atom methods, this gain in efficiency opens up new fields of applications, for example, the complex biological systems mentioned above. Regarding the practicability of the SSIP approach, it is also important to emphasize that it can be seamlessly integrated into an existing finite element framework for solid mechanics. In particular, it does neither depend on any specific beam formulation nor the applied spatial discretization scheme, and in the context of the present work, it has exemplarily been used with geometrically exact Kirchhoff-Love as well as Simo-Reissner type beam finite elements. Likewise, it is independent of the temporal discretization and we have used it along with static and (Lie group) Generalized-alpha time stepping schemes as well as inside a Brownian dynamics framework.
For the proposed SSIP laws, which can either be derived analytically or postulated and fitted to experimental data, we first present the most general form describing the interaction between arbitrarily shaped cross-sections with inhomogeneous distribution of the elementary interacting points (eg, atoms or charges). Subsequently, we focus on homogeneous, circular cross-sections and propose specific, ready-to-use SSIP laws for vdW adhesion, steric repulsion, and electrostatic interaction. Based on the fundamental distinction into either short-range or long-range interactions, we present the required steps and theoretical considerations underlying the analytical derivation of the SSIP laws in a general manner, starting from first principles in form of a point pair interaction potential that is described by a power law with general exponent. Besides the expressions for the total interaction potential, also the corresponding virtual work contributions, its finite element discretization and the consistent linearization are presented. Due to the characteristic singularity of molecular interactions in the limit of zero separation, also, a suitable numerical regularization of the SSIP laws will be proposed.
The remainder of this article is structured as follows. Section 2 briefly summarizes the fundamental concepts and theory of molecular interactions. Along with the fundamentals of the geometrically exact beam theory to be introduced in Section 3, this forms the basis for the novel SSIP approach to be proposed in Section 4. This general approach will then be applied to specific types of physical interactions, namely, vdW, steric, and electrostatic interactions in Section 5. In Section 6, we turn to the finite element discretization of the newly developed numerical methods and discuss some important algorithmic aspects such as the regularization of the reduced interaction laws and the algorithmic complexity of the reduced approach. Finally, the accuracy of the proposed SSIP laws, as well as the general SSIP approach to beam-beam interaction, is validated by means of analytical as well as numerical reference solutions for academic test cases in Section 7.1. In a series of numerical examples including steric repulsion, electrostatic, or vdW adhesion, the effectiveness and robustness of the novel approach are verified in the remainder of Section 7. We conclude the article in Section 8 and provide an outlook to promising future enhancements of the novel approach.

FUNDAMENTALS OF INTERMOLECULAR FORCES AND POTENTIALS
Interactions between molecules may result from various physical origins and are a complex and highly active field of research within the community of theoretical as well as experimental physics. The methods to be derived in this work make use of the most essential and well established findings as summarized, for example, in the textbooks given as References 2 and 3. This section briefly presents a selection of aspects relevant for this work.

Characterization, terminology, and disambiguation
To begin with, a number of universal aspects characterizing molecular interactions, especially with regard to the numerical methods to be developed in this work, shall be presented. A few simple facts about the examples mostly considered throughout this work, namely, electrostatic and vdW interaction, are presented straight away, whereas the details on these and further types of molecular interactions are to be discussed in Section 2.2.

A collection of characteristics of molecular interactions with high relevance for this work
• Type of elementary interaction partners: Interaction may originate from unit charges as in the case of electrostatics.
Other popular examples are vdW effects that are caused by fluctuating dipole interactions occurring in every molecule and hence are related to the molecular density of the material.
• Spatial distribution of elementary interaction partners: Thinking of the resulting interaction between two bodies as accumulation of all molecular interactions, the question for the locations of all elementary interaction partners arises. Charges can often be found on the bodies' surfaces, whereas molecules relevant for vdW interactions spread over the entire volume of the bodies. This work focuses on solid bodies (ie, condensed matter) that are nonconducting such that interaction partners will not redistribute, that is, change their position within a body.
• Distance dependency of the fundamental potential law: Generally, the strength of molecular interactions decays with increasing distance. Most commonly, inverse power laws with different exponents or exponential decay can be identified.
• Range of interactions: As a result of the previous aspects, a range of significant strength of an interaction can be defined. Rather than an inherent property, the classification of long-versus short-range interactions is a theoretical concept to judge the perceptible impact in specific scenarios. Moreover, it is a decisive factor in the derivation of the well-suited numerical methods.
• Additivity and higher order contributions: Many approaches including the one presented in this work make use of superposition, that is, accumulating all the individual contributions from elementary interaction partners to obtain the total effect of interaction. This assumes that the interactions behave additively, that is, the sum of all pairwise interactions describes the overall interaction sufficiently well. More specifically, the presence of other elementary interaction partners in the surrounding must not have a pronounced effect as compared to an isolated system of an interacting pair. Otherwise, the sum of all pairwise interactions would need to be extended by contributions from sets of three, four, and more elementary interaction partners.
As a matter of course, this list is not exhaustive but represents a selection of the most relevant aspects considered in the development of our methods throughout Sections 4 and 5.

Interaction potential and corresponding force
An interaction potential Φ(r), also known as (Gibbs) free energy of the interaction, is defined as the amount of energy required to approach the interaction partners starting from a reference configuration with zero energy at infinite separation. Hence, the following relations between the interaction potential Φ(r) and the magnitude of the force f(r) acting upon each of the partners, each in terms of the distance between both interacting partners r, hold true: Although the final quantities of interest often are the resulting, vectorial forces on slender bodies, it is yet convenient and sensible to consider the scalar interaction potential throughout large parts of this work. This is underlined by the fact that nonlinear finite element methods in the context of structural dynamics can be formulated on the basis of energy and work expressions. Equation (1) expresses the direct and inherent relation between force and potential. Note that the forces emanating from such interaction potentials are conservative and the integral value in Equation (1) is path-independent. Furthermore, the interaction is symmetric in a sense that the force acting upon the first interacting partner f 1 (r) has the same magnitude yet opposite direction as compared to the force acting on the second partner f 2 (r). Using the partners' position vectors x 1 , x 2 ∈ R 3 , we can formulate the vectorial equivalent of the formula above: and f 2 (r) = dΦ(r) dr

Disambiguation
In order to particularize the very general term molecular interactions, we may note that we solely consider interactions between distinct, solid (macro)molecules, that is, no covalent or other chemical bonds, but rather what is sometimes referred to as physical bonds. Thus, we restrict ourselves to intermolecular forces as opposed to intramolecular ones.

Interactions between pairs of atoms, small molecules, or point charges
First principles describing molecular interactions are formulated for a pair of atoms, molecules, or point charges. In the following, all types of interactions to be considered in this work are thus first presented for a minimal system consisting of one pair of these elementary interaction partners.

Electrostatics
Coulomb's law is one of the most fundamental laws in physics and describes the interaction of a pair of point charges under static conditions by Φ elstat (r) = Q 1 Q 2 4 0 1 r , ||f elstat (r)|| = Q 1 Q 2 4 0 1 r 2 , f elstat,1 (r) = where 0 and are the vacuum and dielectric permittivity, respectively. For the sake of brevity in any later usage, let us define the abbreviation C elstat ∶= (4 0 ) −1 . Depending on the signs of Q 1 and Q 2 , electrostatic forces may either be repulsive or attractive. Besides a pair of point charges, Coulomb's law likewise holds for a pair of spherically symmetric charge distributions with resulting charges Q 1 and Q 2 , respectively. This is an important insight since ultimately we are interested in interactions between two bodies with finite extension rather than points. Furthermore, interactions between rigid spheres and rigid bodies are of interest for applications such as particle diffusion in hydrogels. Throughout the entire work, no electrodynamic effects shall be considered. This is a valid assumption as long as bodies are nonconductive and the motion of bodies carrying the attached charges happens on much larger time scales than relevant eigenfrequencies in electrodynamics. Due to the inverse-first power law, the electrostatic potential has quite a long range, meaning that two point charges at a large distance still experience a considerable interaction force as compared to small distances. This behavior is even more pronounced for the interaction of two extended bodies, where the whole lot of all distant point pairs dominates the total interaction energy as compared to the few closest point pairs. This property is crucially different as compared, for example, to vdW interactions considered in the next section. We will account for and indeed make use of this important property in the development of the methods to be presented in this work.

VdW interactions
VdW forces originate from charge fluctuations, thus being an electrodynamic effect caused by quantum-mechanical uncertainties in positions and orientations of charges. Depending on the interaction partners, three subclasses can be distinguished as Keesom (two permanent dipoles), Debye (one permanent dipole, one induced dipole), and London dispersion interactions (two transient dipoles). The ubiquitous nature of vdW interactions is due to the fact that the latter contribution even arises in neutral, nonpolar, yet polarizable matter that means basically every atom or molecule. All three kinds of dipole interactions can be unified in that their interaction free energy follows an inverse-sixth power law in the separation 3 : This is a pleasantly simple expression, yet intricate when it comes to transferring it to two-body interactions, as we will discuss in Section 2.3.2. In general, vdW forces are always attractive for two identical or similar molecules, yet may be repulsive for other material combinations.

Steric exclusion
Two approaching atoms or molecules will at some very small separation suddenly experience a seemingly infinite repulsive force. This effect is attributed to the overlap of electron clouds and referred to as steric repulsion, steric exclusion, or hard core repulsion. Without thorough theoretical foundation, several (almost) infinitely steep repulsive potential laws are empirically used to model this phenomenon. The first option is a hard wall/core/sphere potential, which has a singularity at zero separation lim r→0 Φ c,hs (r) = ∞ and Φ c,hs (r) ≡ 0 for r > 0.
Other common choices include a power-law potential with a large integer exponent n c,pow and finally an exponential potential Φ c,exp (r) = C c,exp e −r∕r c,exp .
Note that the former two coincide in the limit n c,pow → ∞ and increase indefinitely for r → 0 while the exponential one does not. Generally, this behavior of steric exclusion is in good agreement with our intuition based on macroscopic solid bodies coming into contact.

Total molecular pair potentials and force fields
In many systems of interest, any two or more of the aforementioned effects may be relevant at the same time such that a combination of the pair potentials is required. This is typically done by summation of the individual potential contributions and leads to a total intermolecular pair potential. Among the large number of possible combinations, * the LJ potential is probably the most commonly used variant (see Figure 1).
It is a special case of the Mie potential Φ Mie (r) = C Mie,m r −m − C Mie,n r −n with exponents being chosen to model the inverse-sixth vdW attraction on the one hand and a strong repulsion on the other hand. The parameters can be identified as the minimal value Φ LJ,eq < 0 that the LJ potential takes at equilibrium separation r LJ,eq > 0, that is, at the separation where the resulting force is zero. *See p. 138 of Reference 2 for a comprehensive list of combinations used as total pair potential laws.
F I G U R E 1 Lennard-Jones interaction potential for a pair of points, that is, atoms Other important quantities characterizing the LJ force law are the minimal force value f LJ,min and corresponding distance r LJ,f min f LJ,min ≈ 2.6899 Φ LJ,eq r LJ,eq and r LJ,f min = ( 13 7 ) 1∕6 r LJ,eq ≈ 1.1087 r LJ,eq .
The minimal force, that is, the maximal adhesive force, is commonly referred to as pull-off force. Israelachvili 2 also points out the chance of a fortunate cancellation of errors in total pair potentials, especially close to the limit r → 0. In this regime, attractive forces tend to be underestimated by the simplified inverse-sixth term but likewise the steric repulsion is probably stronger than estimated from the power law. Both errors cancel rather than accumulate and increase the model accuracy.

Remarks.
1. Many of the presented point-point interaction potentials decay rapidly with the distance as shown exemplarily for a law Φ(r) ∝ r −12 in Figure 2. In anticipation of the numerical methods to be proposed in this work, we can already state at this point that these extreme gradients are very challenging for numerical quadrature schemes that will, therefore, be discussed in a dedicated Section 6.4. 2. In molecular dynamics, a force field is typically used instead of the potential law to model the total interaction of a pair of atoms. Specific forms are being proposed for (coarse-grained) force fields modeling the interaction of macromolecules such as DNA instead of atoms. Since these all-atom approaches follow an entirely different approach as compared to the continuum model proposed here, we will not discuss force fields any further at this point.

Two-body interaction: Surface vs volume interaction
In this section, we take the important step from interacting point pairs to interactions between two bodies with defined spatial extension containing many of the fundamental pointlike interaction partners considered throughout the preceding Section 2.2. The obvious question for the spatial distribution of the interaction partners leads to the important distinction F I G U R E 2 Example of a point-point interaction potential Φ(r) ∝ r −12 plotted over a circular domain (blue circle) with (A) small and (B) large distance to the pointlike interaction partner (red dot). Note the huge difference in scales between surface and volume interactions. As the name suggests, in the first case, the elementary interaction partners are distributed over the surface of the bodies but not in the interior. The most important examples from this category are electrostatic interactions between bodies where the charges sit on the surfaces and are not free to move around. This applies to a large number of charged, nonconductive biopolymer fibers such as actin or DNA. In the second case of volume interactions, the elementary interacting partners are distributed over the entire volume of the bodies. The most important examples here are vdW interactions and steric exclusion. As compared to surface interactions, this further increases the dimensionality of the problem making it more challenging to tackle, both by analytical as well as numerical means. In terms of notation, one may also find the expressions body forces or bulk interaction referring to this category of interactions. Let us briefly look at volume and surface interactions as an abstract concept, leaving aside the specifics of the underlying physical effects that are to be discussed in Sections 2.3.1 and 2.3.2. Likewise, we assume additivity here and discuss the applicability later with the physical type of interaction. Since volume interactions are the more general and challenging case, we will discuss most aspects and approaches first for volume interactions and later only point out the differences considering surface interactions throughout this article. Figure 3 schematically visualizes the distribution of elementary interaction partners within two macromolecular or macroscopic bodies.
Assuming additivity, we apply pairwise summation to arrive at the two-body interaction potential Further assuming a continuous atomic density i , i = 1, 2, the total interaction potential can alternatively be rewritten as integral over the volumes V 1 and V 2 of both bodies  1 and  2 : It can be shown that this continuum approach is the result of coarse-graining, that is, smearing-out the discrete positions of atoms in a system into a smooth atomic density function (x). 23 In the case of surface interactions, the dimensionality of the problem reduces and summation or integration is carried out over both bodies' surfaces  1 ,  2 : Accordingly, surface densities i (x i ), i = 1, 2, replace the volume densities in this case.
The range of two-body interaction forces originating from point pair potentials: Let us assume a general inverse power law Φ(r) = kr −m for the point pair interaction potential. It is obvious that the potential becomes infinitely large if the separation of the two individual points r approaches zero and, on the other hand, that the potential rapidly decays with increasing distance. Turning to two bodies of finite size, that is, two clouds of points, things are more involved as the following theoretical considerations demonstrate. In short, it can be shown that there is a fundamental difference between potentials with an exponent m ≤ 3 on the one hand and m > 3 on the other hand. Starting with the case m > 3, for example, vdW interactions, the two-body interaction potential goes to infinity if the bodies approach until their surfaces touch each other. This can be illustrated by the simple example of two spheres of radius R, where the vdW interaction potential scales with Π vdW ∝ g −1 (cf. Reference 2 [p. 255]) with surface-to-surface F I G U R E 3 Two arbitrarily shaped, deformable bodies  1 and  2 with volumes V 1 and V 2 and continuous particle densities 1 and 2 separation or gap g = d − 2R and the scalar distance between the spheres' centers d. This singularity of the two-body interaction potential in the limit of zero separation g → 0 is due to the fact that potentials with m > 3 decay so rapidly that the few point pairs with smallest separation outweigh the potentially very large number of all other, distant point pairs in terms of their potential contributions. Therefore, we can conclude that potentials with m > 3 have no significant large distance contribution and the two-body interaction potential is governed by the separation of any two closest points (and their immediate surrounding). Considering the example of two cylinders later on, we will also see that the vdW interaction potential of two perpendicular cylinders does not change perceptibly if the length of the cylinders is increased (cf. Equation (16)), which can again be attributed to the short range of vdW interactions.
The situation is substantially different for potentials with m ≤ 3, for example, Coulombic interactions. Here, the total contribution of all distant point pairs dominates over the few closest point pairs and the total interaction potential remains finite even if both bodies are in contact. Looking once again at the simple example of two spheres, Coulomb law (cf. Equation (3)) directly shows Π elstat ∝ d −1 , and thus, no singularity occurs for (nearly) contacting surfaces g → 0, that is, d → 2R. Also, in contrast to the case of vdW mentioned above, the Coulomb interaction potential of two perpendicular cylinders would increase if their length is increased. The underlying theoretical derivations revealing also the turning point, that is, the exponent m = 3 was first noted by Newton and can be found in p. 11 of Reference 2. Due to this crucial difference, potentials with m > 3 will be denoted as short-range interactions (eg, repulsive as well as attractive part of LJ) and potentials with m ≤ 3 as long-range interactions (eg, Coulomb) throughout this work.

Electrostatics of nonconductive bodies: An example for long-range surface interactions
The Coulomb interaction is additive such that the net force acting on an individual point charge in a system of point charges can be calculated from superposition of all pair-wise computed force contributions. 2 Equivalently, the net interaction potential results from summation of all pair potentials. A large body of literature deals with the problem of electrostatic multibody interaction. One concept of high relevance for the present work is a well-known strategy called multipole expansion, which aims to express the resultant electrostatic potential of a (continuous) charge distribution as an (infinite) series (see Reference 32 for details). The individual terms of the series expansion generally are inverse power laws in the distance with increasing exponent and referred to as mono-, di-, quadru-, up to n-pole moments. At points far from the location of the charge cloud, the series converges quickly and can thus be truncated in good approximation. Regarding the total interaction potential of two charged bodies as formulated in Equation (12) or (13), this already outlines the way how to determine Π elstat for trivial geometries of the interacting bodies, where the integrals can potentially be solved analytically. We will return to this concept in the context of deformable slender fibers when proposing the general SSIP approach in the beginning of Section 4 and make use of a (truncated) multipole expansion of the charged cross-sections for the (simplified) SSIP law for long-range surface interactions to be proposed in Section 5.3.

VdW interaction: An example for short-range volume interactions
Here, we want to discuss vdW interactions as one example of physically relevant volume interactions, which is based on the inverse-sixth power law given in Equation (4). However, very similar considerations and formulas apply to steric interactions as well as LJ interactions. Today, we know that vdW interactions are generally nonadditive. The latest and most accurate models for two-body vdW interactions are based on the Lifshitz theory and, among other effects, include retardation, anisotropy, and differences in polarizability. Nevertheless, a "happy convergence" of old, that is, Hamaker's pairwise summation, and new, that is, Lifshitz theory allows determining the distance dependency from pairwise summation and then estimate the prefactor, that is, Hamaker "constant" A Ham , from more advanced modern theory. This approach yielding a so-called Hamaker-Lifshitz hybrid form (cf. Reference 3, 2 (p. 257)) is what motivates us to use pairwise summation in the derivation of the numerical methods to be proposed in the present work. Also, there are some special scenarios (negligible retardation, negligible difference in optical properties of the bodies, interaction in vacuum, etc), where additivity can be assumed as a good approximation even without adaption of the Hamaker constant.
Generally, even the simple approach of pairwise summation requires two nested 3D integrals over both bodies' volumes, that is, 6D integration. Mainly due to this high dimensionality of the problem, unfortunately, (closed-form) analytical solutions can only be obtained for some simple special cases. Still, careful considerations and selection allow us to exploit some of these analytical expressions in order to develop efficient, reduced methods in Section 4.2. To get a concise overview of all expressions relevant for the remainder of this work, we will provide a collection of closed-form analytical solutions in the following.
First, we want to look at two cylinders representing the simplest model for two interacting straight, rigid fibers with circular cross section. A number of publications consider this scenario, and due to the simplicity of the geometry, they were able to derive analytical solutions for some special cases. The resulting expressions are summarized in Table 1 and will be used for verification purposes in Section 7.1. A second, highly relevant scenario is the one considering two disks. These analytical expressions, summarized in Table 2, will be beneficial, and in fact provide the main ingredient, for the SSIP approach to describe molecular interactions between deformable fibers modeled as 1D Cosserat continua.

Two cylinders.
To begin with, we consider the cases of parallel and perpendicular cylinders. Generally, the cylinders are assumed to be infinitely long, such that the boundary effects at its ends may be neglected. As the interaction potential for parallel cylinders would be infinite, one typically considers a length-specific interaction potential̃v dW,cyl||cyl with dimensions of energy per unit length. This quantity thus describes the interaction of one infinitely long cylinder with a section of unit length of the other infinitely long cylinders. For perpendicular orientation (and all other mutual angles apart from = 0), on the other hand, the total interaction potential Π vdW,cyl⊥cyl remains finite.
Even for this simple case of two cylinders, no closed-form analytical solution for the vdW interaction energy can be found for all mutual angles and all separations. One thus resorts to the consideration of the limits of small and large surface-to-surface separations for which the general solution, an infinite series, converges to the expressions presented in Table 1.
Despite the different dimensions of the quantities for parallel and perpendicular cylinders, we can still compare these expressions as becomes clear by the following thought experiment. Considering two "sufficiently long" cylinders TA B L E 1 A collection of analytical solutions for the cylinder-cylinder interaction potential derived via pairwise summation

Limit of Small Separations
Limit of Large Separations see p. 255 of Reference 2 and p. 172 of Reference 3̃v dW,cyl||cyl,ls = − see p. 255 of Reference 2 see p. 16 of Reference 3 Note: R i denotes the cylinder radii, d denotes the closest distance between the cylinder axes, g denotes the surface-to-surface separation, that is, gap, and A Ham ∶= 2 1 2 C vdW is the commonly used abbreviation known as Hamaker constant, where i denotes the particle densities and C vdW denotes the constant prefactor of the point-pair potential law (see Equation (4)).

Limit of Small Separations
Limit of Large Separations see Reference 33̃v see Reference 33 Note: R i denotes the disk radii, d denotes the closest distance between the disk midpoints, g denotes the surface-to-surface separation, that is, gap, and A Ham ∶= 2 1 2 C vdW is the commonly used abbreviation known as Hamaker constant, where i denotes the particle densities and C vdW denotes the constant prefactor of the point-pair potential law (see Equation (4)). of length L in parallel orientation, the total interaction potential is well described by Π vdW,cyl||cyl =̃v dW,cyl||cyl ⋅ L and thus shows the same scaling behavior in the separation as Equations (14) and (15). In addition, Equations (16) and (17) are also a good approximation for the perpendicular orientation of these cylinders of finite length L since the difference in the distant point pairs is negligible.
We would like to point out just a few interesting aspects of these equations. First, it is remarkable how the expressions differ in the exponent of the power law describing the distance dependency of the potential. This relates to a diverse and highly nonlinear behavior already for this simplest model system of fiber-fiber interactions composed of two cylinders. Second, the parallel orientation is a very special orientation that gives rise to the strongest possible adhesive forces between two cylinders and at the same time is the only stable equilibrium configuration. Third, the distance scaling behavior of two parallel cylinders at small separations̃v dW,cyl||cyl,ss ∝ g − 3 2 lies between the fundamental solutions known for two infinite half spaces̃∝ g −2 (double tilde indicates a potential per unit area) and two spheres Π ∝ g −1 . Note that again multiplication of these laws by a length or area does not alter the scaling law in the distance. Looking at the equations for large separations, we see similar relations, once again with a stronger distance scaling behavior in the parallel case. Generally, the solutions for large separations are expressed more naturally in the interaxis separation d rather than the surface-to-surface separation g.

Two disks.
This problem has been studied in the literature on vdW interaction of straight, rigid cylinders of infinite length. 33 In analogy to the cylinder-cylinder scenario, it turns out that not even in the simplest case of parallel oriented disks, that is, two disks with parallel normal vectors, a closed analytical solution can be found for all separations. Instead, two expressions for the limit of small and large separations of the disks g as compared to their radii R 1 and R 2 are presented in Table 2, respectively.
To summarize, a closed-form expression for the two-body vdW interaction potential is only known for some rare special cases and the ones relevant for fiber-fiber interactions have been identified in the voluminous literature on this topic and presented here in a brief and concise manner.
We would like to conclude this section on two-body vdW interactions with a note on the analogy to steric exclusion, that is, contact interactions, as already discussed for point pairs in Section 2.2.3. This class of physical interactions shares the two central properties of being extremely short in range and being volume interactions. Starting from an inverse-12 power law as in the repulsive part of the LJ interaction law, one may apply very similar solution strategies and finally obtain very similar expressions as the ones presented in this section. For the sake of brevity, we refer to the derivations in Appendix A and the analysis of the resulting total LJ interaction that will also be used for the regularization of the reduced potential laws in Section 6.3.

FUNDAMENTALS OF GEOMETRICALLY EXACT 3D BEAM THEORY
This section aims to provide a brief and concise introduction to well-known concepts of beam theory to be used in the remainder of this article. As commonly used in engineering mechanics, we refer to the term beam as a mathematical model for a 3D, slender, deformable body for which the following assumption can be made. The much larger extent of the body in its axial direction as compared to all transverse directions often justifies the Bernoulli hypothesis of rigid and therefore undeformable cross sections. This, in turn, allows for a reduced dimensional description as a 1D Cosserat continuum embedded in the 3D Euclidian space. The so-called Simo-Reissner beam theory dates back to the works of Reissner, 34 Simo, 35 and Simo and Vu-Quoc, 36 who generalized the linear Timoshenko beam theory 37 to the geometrically nonlinear regime. Since the Simo-Reissner model, which accounts for the deformation modes of axial tension, bending, torsion, and shear deformation, is the most general representative of geometrically exact beam theories, we choose it as the one to be used exemplarily throughout this work. Nevertheless, the novel approach to be proposed is not restricted to a specific beam formulation. We have likewise applied it to formulations of Kirchhoff-Love type, which are known to be advantageous in the regime of high slenderness ratios where the underlying assumption of negligible shear deformation is met. 6,38 Refer to Section 6.1 for more details.

Geometry representation
A certain configuration of the 1D Cosserat continuum is uniquely defined by the centroid position and the orientation of the cross-section at every point of the continuum. The set of all centroid positions is referred to as centerline or neutral axis and expressed by the curve in space and time t ∈ R. Each material point along the centerline is represented by a corresponding value of the arc-length parameter s ∈ [ 0, l 0 ] =∶ Ω l ⊂ R. Note that this arc-length parameter s is defined in the stress-free, initial configuration of the centerline curve r 0 (s) = r(s, t = 0). Thus, the norm of the initial centerline tangent vector yields but generally, in the presence of axial tension, ||r ′ (s, t)|| = || r(s, t)∕ s|| ≠ 1. Furthermore, the cross-section orientation at each of these material points is expressed by a right-handed orthonormal frame often denoted as material triad: The second and third base vectors follow those material fibers representing the principal axes of inertia of area. Such a triad can equivalently be interpreted as rotation tensor transforming the base vectors of a global Cartesian frame E i ∈ R 3 , i ∈ 1, 2, 3 into the base vectors of the material triad g i ∈ R 3 , i ∈ 1, 2, 3 via In summary, a beam's configuration may be uniquely described by a field of centroid positions r(s, t) and a field of associated material triads (s, t), altogether constituting a 1D Cosserat continuum (see Figure 4). According to this concept of geometry representation, the position x of an arbitrary material point P of the slender body is obtained from Here, the additional convective coordinates s 2 and s 3 specify the location of P within the cross section, i. e., as linear combination of the unit direction vectors g 2 and g 3 . For a minimal parameterization of the triad, for example, the three-component rotation pseudovector may be used such that we end up with six independent degrees of freedom (r, ) to define the position of each material point in the body by means of Equation (24).

F I G U R E 4 Geometry description and kinematics of the
Cosserat continuum formulation of a beam: initial, that is, stress-free (blue) and deformed (black) configuration. Straight configuration in initial state is chosen exemplarily here without loss of generality Remark on notation. Unless otherwise specified, all vector and matrix quantities are expressed in the global Cartesian basis E i . Differing bases, for example, the material frame, are indicated by a subscript [.] g i . Quantities evaluated at time t = 0, that is, the initial stress-free configuration, are indicated by a subscript 0, for example, in r 0 (s). Differentiation with respect to the arc-length coordinate s is indicated by a prime, for example, for the centerline tangent vector r ′ (s, t) = dr(s, t)∕ds. Differentiation with respect to time t is indicated by a dot, for example, for the centerline velocity vectoṙr(s, t) = dr(s, t)∕dt. For the sake of brevity, the arguments s, t will often be omitted in the following.
Remark on finite 3D rotations. To a large extent, the challenges and complexity in the numerical treatment of the geometrically exact beam theory can be traced back to the presence of large rotations. In contrast to common vector spaces, the rotation group SO(3) is a nonlinear manifold (with Lie group structure) and lacks essential properties such as additivity and commutativity, which makes standard procedures such as interpolation or update of configurations much more involved. While Section 4 introduces the concept of section-section interaction laws in the most general manner, in Section 5, some additional (practically relevant) assumptions are made, which allow formulating the interaction laws as pure function of the beam centerline configuration. In turn, this strategy will allow to avoid the handling of finite rotations and to achieve simpler and more compact numerical formulations. Figure 4 summarizes the kinematics of geometrically exact beam theory. Based on these kinematic quantities, deformation measures as well as constitutive laws can be defined. Finally, the potential of the internal (elastic) forces and moments Π int is expressed uniquely by means of the set of six degrees of freedom (r, ) at each point of the 1D Cosserat continuum. See References 4-6 for a detailed presentation of these steps.

THE SECTION-SECTION INTERACTION POTENTIAL APPROACH
Based on the fundamentals of molecular interactions (Section 2) as well as geometrically exact beam theory (Section 3), this section will propose the novel SSIP approach to model various types of molecular interactions between deformable fibers undergoing large deflections in 3D.

Problem statement
For a classical conservative system, the total potential energy of the system can be stated taking into account the internal and external energy Π int and Π ext . The additional contribution from molecular interaction potentials Π ia is simply added to the total potential energy as follows: Note that the existing parts remain unchanged from the additional contribution. One noteworthy difference is that internal and external energy are summed over all bodies in the system, whereas the total interaction free energy is summed over all pairs of interacting bodies.
According to the principle of minimum of total potential energy, the weak form of the equilibrium equations is derived by means of variational calculus. The very same equation may alternatively be derived by means of the principle of virtual work which also holds for nonconservative systems: Clearly, the evaluation of the interaction potential Π ia , or rather its variation Π ia , is the crucial step here. Recall Equation (12) to realize that it generally requires the evaluation of two nested 3D integrals. † The direct approach using 6D numerical quadrature turns out to be extremely costly and in fact inhibits any application to (biologically) relevant multibody systems. See Section 6.4 for more details on the complexity and the cost of this naive, direct approach as well as the novel SSIP approach to be proposed in the following.

4.2
The key to dimensional reduction from 6D to 2D We propose a split of the integral in the length dimensions l 1 , l 2 on the one hand and the cross-sectional dimensions A 1 , A 2 on the other hand: Exploiting the characteristic slenderness of beams, the 4D integration over both undeformable cross-sections shall be tackled analytically and only the remaining two nested 1D integrals along the centerline curves shall be evaluated numerically to allow for arbitrarily deformed configurations. Generally, we follow the key idea of reduced dimensionality from beam theory and thus aim to express the relevant information about the cross-sectional dimensions by the pointwise six degrees of freedom (r i , i ) of the 1D Cosserat continua without loss of significant information. To this end, we need to consider the resulting interaction between all the elementary interaction partners within two cross-sections expressed by an interaction potential̃(r 1−2 , 1−2 ) that depends on the separation r 1−2 of the centroid positions and the relative rotation 1−2 between both material frames attached to the cross-sections. For this reason, the novel approach is referred to as the SSIP approach. The SSIP̃is a double length-specific quantity in the way that it measures an energy per unit length of beam 1 per unit length of beam 2, which is indicated here by the double tilde. The sought-after total interaction potential Π ia of two slender deformable bodies thus results from double numerical integration of the double length-specific SSIP along both centerline curves: This relation suggests another, alternative interpretation of the SSIP̃. In analogy to the term intersurface potential, introduced by Reference 22,̃can be understood as an interaxis potential, that is, describing the interaction of two spatial curves (with attached material frames).
To further illustrate this novel concept, a simple, demonstrative example is shown in Figure 5. In this scenario of two beams with circular cross-section, the SSIP̃(r 1−2 , 1−2 ) describes the interaction of two circular disks at arbitrary mutual distance and orientation. To evaluate the two nested 1D integrals along the beam axes numerically, the SSIP needs to be evaluated for all combinations of integration points (denoted here as Gaussian quadrature points (GP), without loss of generality). For one of these pairs ( 1,GP , 2,GP ), the geometrical quantities are shown exemplarily.
While analytical integration of the inner 4D integral of Equation (27) has already been suggested above as one way to find a closed-form expression for the SSIP̃, we would like to stress the generality of the SSIP approach at this point. The question of how to find̃is independent of the strategy to determine the interaction energy Π ia of two slender bodies via numerical double integration as proposed in this section. This is important to understand because the SSIP̃will obviously depend on the type of interaction as well as the cross-section shape and a number of other factors and there might also be cases where no analytical solution can be obtained and one wants to resort to relations fitted to experimental data. In the scope of this work, several specific expressions of̃, for example, for vdW as well as electrostatic interactions, will be derived analytically in Section 5. † It is important to mention that, assuming additivity of the involved potentials, systems with more than two bodies can be handled by superposition of all pairwise two-body interaction potentials. It is thus sufficient to consider one pair of beams in the following. The same reasoning applies to more than one type of physical interaction, that is, potential contribution. F I G U R E 5 Illustration of the novel section-section interaction potential (SSIP) approach: two cross-sections at integration points 1/2,GP of beams 1 and 2, respectively, their separation r 1−2 and relative rotation 1−2 In its most general form,̃will be a function of the relative displacement r 1−2 and the relative rotation 1−2 between both cross-sections, that is, three translational and three rotational degrees of freedom. This becomes clear if one recalls that the position x P of every material point in a slender body can be uniquely described by the six degrees of freedom of a 1D Cosserat continuum (cf. Equation (24)). Thus, keeping one cross-section fixed, the position x P1−P2 of every material point in the second cross-section relative to the (centroid position and material frame of the) first cross-section is again uniquely described by six degrees of freedom (r 1−2 , 1−2 ). This insight naturally leads to the interesting question under which conditions the SSIP̃can be described by a smaller set of degrees of freedom, thus simplifying the expressions. Rotational symmetry of the interacting cross-sections is one common example where the SSIP would be invariant under rotations around the cross-section's normal axis. We will return to this topic in Section 5.1 as a preparation for the following derivation of specific expressions for the SSIP.
Remark on the included special case of surface interactions. It is very convenient that the practically highly relevant case of surface potentials is already included as a simpler, special case in the proposed SSIP approach to model molecular interactions between the entire volume of flexible fibers. In simple words, it is sufficient to omit one spatial dimension of analytical integration on each interacting body in the analytical derivation of the required SSIP̃. More specifically, this means that̃may be obtained from solving analytically two nested 1D integrals along both, for example, ring-shaped, contour lines of the fiber cross-sections.
Remark on self-interactions. Note that self-interaction, that is, the interaction of distinct parts of the same fiber, can be accounted for if we evaluate also interactions between those cross-sections that belong to one and the same beam. Leaving everything else of the SSIP approach unchanged, this is naturally achieved by including those fiber pairs where the first fiber is identical to the second fiber. Including self-interaction in this manner is considered to be important whenever the fibers likely exhibit self-contact or generally whenever distinct (nonneighboring) parts of the same fiber come close to each other. If this is not the case, an alternative modeling strategy is to include these contributions from self-interaction in the constitutive model of the beam theory and thus consider them as part of an effective elastic stiffness as has been suggested in References 23 and 39. The latter approach has been applied in the numerical examples of Section 7.

APPLICATION OF THE GENERAL SSIP APPROACH TO SPECIFIC TYPES OF INTERACTIONS
At this point, we would like to return to the fact that the SSIP approach proposed in Section 4 is general in the sense that it does not depend on the specific type of physical interaction. This section provides the necessary information and formulae to apply the newly proposed, generally valid approach from Section 4.2 to certain types of real-world, physical interactions such as electrostatics or vdW. As mentioned above, the approach requires a closed-form expression for the SSIP̃. We basically see two alternative promising ways to arrive at such a reduced interaction law̃: 1. analytical integration, for example, as presented in Section 5.2; 2. postulatẽas a general function of separation r 1−2 and mutual orientation 1−2 and determine the free parameters via fitting to: (a) experimental data for specific section-section configurations, that is, discrete values of separation r 1−2 and mutual orientation 1−2 ; (b) data from (one-time) numerical 4D integration for specific section-section configurations, that is, discrete values of r 1−2 and 1−2 ; (c) experimental data for the global system response, for example, of the entire fiber pair or a fiber network.
As a starting point, we will restrict ourselves to the first option based on analytical integration throughout the remainder of this work. See Section 5.2 for an example of the further steps required to derive the final, ready-to-use expressions in case of vdW interactions. To give but one example for a recent experimental work, which could serve as the basis for the second option listed above, we refer to Reference 40 measuring cohesive interactions between a single pair of microtubules. Postulating an SSIP and studying the global system response in numerical simulations could also be used as a verification of theoretical predictions for the system behavior. To give an example, we refer to the work of theoretical biophysicists studying the structural polymorphism of the cytoskeleton resulting from molecular rod-rod interactions, 41 which is based on a postulated model potential "that captures the main features of any realistic potential". In summary, we see a large number of promising future use cases for the proposed SSIP approach.

Additional assumptions and possible simplification of the most general form of SSIP laws
Recall that the most general form of the SSIP is uniquely described by a set of six degrees of freedom, three for the relative displacement and three for the relative orientation of the two interacting cross-sections, as presented in Section 4.2. The following assumptions turn out to significantly simplify this most general form of the SSIP law by reducing the number of relevant degrees of freedom from six to four, two, or even one. This, in turn, eases the desirable derivation of analytical closed-form solutions of the SSIP̃based on the point pair potentials Φ(r) presented in Section 2.2. Specifically, these assumptions are: 1. undeformable cross-sections; 2. circular cross-section shapes; 3. homogeneous (or, more generally, rotationally symmetric) particle densities 1 , 2 in the cross-sections or surface charge densities 1 , 2 over the circumference.
The first assumption is typical for geometrically exact beam theory and the second and third assumptions are reasonable regarding our applications to biopolymer fibers such as actin or DNA that can often be modeled as homogeneous fibers with circular cross-sections. Based on these three assumptions, we can conclude that the interaction between two cross-sections is geometrically equivalent to the interaction of two homogeneous, circular disks (or rings in case of surface interactions). The rotational symmetry of the circular disks then implies that the interaction potential is invariant under rotations around their own axes and thus reduces the number of degrees of freedom to four. The relative importance of the remaining degrees of freedom, that is, modes, will be the crucial point in the following discussion, where we turn to the interaction of two slender bodies, that is, consider the entirety of all cross-section pairs. At this point, recall the fundamental distinction between either short-range or long-range interactions as outlined in Section 2.3.2.
In the case of short-range interactions, the cross-section pairs in the immediate vicinity of the mutual closest points of the slender bodies dominate the total interaction. As is known from macroscopic beam contact formulations, 12,20,21 the criterion for the closest point is that the distance vector r 1−2 is perpendicular to both centerline tangent vectors r i ′ , that is, assuming small shear angles, the normal vectors of the disks (see Figure 6A). Since only cross-section pairs in the direct vicinity of the closest points are relevant, arbitrary relative configurations (ie, separations and relative rotations) between those cross-sections shall in the following be discussed on the basis of six alternative degrees of freedom as illustrated in Figure 6A. By considering the cross-sections A 1c and A 2c at the closest points as reference, the relative configuration between cross-sections in the direct vicinity of A 1c and A 2c can be described via (small) rotations of A 1c around the axis r ′ 1 (angle 1 ) and r 1−2 × r ′ 1 (angle 1 ), (small) rotations of A 2c around the axes r ′ 2 (angle 2 ) and r 1−2 × r ′ 2 (angle 2 ), (small) relative rotations between A 1c and A 2c around the axis r 1−2 (angle ), and (small) changes in the (scalar) distance d = ||r 1−2 ||. As a consequence of Assumptions 1-3 discussed above, the considered interaction potentials are invariant under rotations 1 and 2 . From the remaining four degrees of freedom, the scalar distance d clearly has the most significant influence on the interaction potential because changes in d directly affect the mutual distance r of all point pairs in the body and, most importantly, the smallest surface separation g between both bodies. The second most significant influence is expected for the scalar relative rotation angle between the cross-section normal vectors, that is, cos( ) = n 1 A change in does not alter the gap g but influences the distance of all next nearest point pairs in the immediate vicinity of the closest surface point pair. For the remaining two relative rotations 1 and 2 , arguments for both sides, either significant or rather irrelevant influence on the total interaction potential, can be found at this point. On the one hand, the orthogonality conditions r ′ i ⋅ r 1−2 = 0, i = 1, 2 at the closest points are fulfilled in good approximation also for cross-sections in the direct vicinity of the closest points, such that the influence of 1/2 could be considered negligible. On the other hand, even small rotations 1/2 change the smallest separation of any two point pairs in the immediate neighborhood of the closest point pair as soon as the centroid distance vector r 1−2 rotates out of the two cross-section planes. Therefore, it seems hard to draw a final conclusion with respect to the influence and thus importance of 1/2 based on the qualitative theoretical considerations of this section.
To summarize, the scalar distance d between the cross-section centroids, the scalar relative rotation angle between the cross-section normal vectors, and possibly also the relative rotation components 1/2 are supposed to have a perceptible influence on the short-range interaction between slender beams fulfilling Assumptions 1-3, with a relative importance that decreases in this order. In favor of the simplest possible model, we will, therefore, assume at this point that the effect of this scalar relative rotations , 1 , and 2 is negligible as compared to the effect of the scalar separation d. This allows us to directly use the analytical, closed-form expression for the disk-disk interaction potential as presented in Section 2.3.2. The error for arbitrary configurations associated with this model assumption will be thoroughly analyzed in Section 7.1.1. In this context, it is a noteworthy fact that the first published method for 2D beam-rigid half space LJ interaction 27 likewise neglects the effect of cross-section orientation. In the subsequent publication, 28 the effect of cross-section rotation, that is, interaction moments, has finally been included and a quantitative analysis considering a peeling experiment of a Gecko spatula revealed that the differences in the resulting maximum peeling force and bending moment are below 8% and 2%, respectively. However, it is unclear whether this assessment also holds for beam-beam interactions modeled via the proposed SSIP approach. Including the orientation of the cross-sections thus is a work in progress and will be addressed in a subsequent publication. Finally, it is emphasized that by the discussed assumptions the SSIP law̃as well as the total two-body interaction potential Π ia can be formulated as pure function of the beam centerlines r 1 and r 2 without the necessity to consider cross-section orientations via rotational degrees of freedom. This is considered a significant simplification of the most general case of the SSIP approach and thus facilitates both the remaining derivations in the present work as well as potential future applications.
Remark on configurations with nonunique closest points. It is well-known from the literature on macroscopic beam contact that the location of the closest points is nonunique for certain configurations of two interacting beams, for example, the trivial case of two straight beams, where an infinite number of closest point pairs exists (see Reference 20). Note, however, that the reasoning presented above also holds in these cases since the cross-section pairs in either one or several of these regions will dominate the total interaction potential. In the case of long-range interactions, the situation is fundamentally different. Recall that here the large number of cross-section pairs with large separation d ≫ R outweighs the contributions from those few pairs in the vicinity of the closest point and dominates the total interaction. Thus, the regime of large separations is decisive in this case and it has already been shown in the literature considering disk-disk interaction (see the brief summary in Appendix A2.1) that in this regime, the exact orientation of the disks can be neglected as compared to the centroid separation d. In simple terms, this holds because the distance x P1−P2 between any point in disk 1 and any point in disk 2 may be approximated by the centroid separation d if d is much larger than the disk radii R i , which again holds for the large majority of all possible cross-section pairs. The validity of this assumption will be thoroughly verified by means of numerical reference solutions in Section 7.1.2.
Remark. The following, similar reasoning from the perspective of slender continua comes to the same conclusion. As visualized in Figure 6B, even pure (rigid body) rotations of slender bodies always entail large displacements of the centerline ‡ in the region far away from the center of rotation. The displacement of any material point due to cross-section rotation will be in the order of ΩR, where Ω is the angle of rotation and R denotes the cross-section radius, whereas the displacement due to centerline displacement will be in the order of ΩL, where L is the distance from the center of rotation and thus in the order of the beam length l. Due to the high slenderness l∕R ≫ 1 of beams, the displacement from translation of the centroid will dominate in the region of large separations with L ≫ R, which is the decisive one here, because it includes the large majority of all possible cross-section pairs, as outlined above. The original, analogous reasoning has been applied to the relative importance of translational versus rotational contributions to the mass inertia of beams.
To conclude, we have discussed the possibility of defining and using SSIP laws̃as a function of the scalar separation of the centroids d instead of the six degrees of freedom in the most general form. This significantly simplifies the theory because the analytical solutions for the planar disk-disk interaction from literature can directly be used and the complex treatment of large rotations is avoided. Having considered the additional assumptions above in the context of short-range interactions, the relative importance of cross-section rotations still needs to be verified in the subsequent quantitative analysis of Section 7.1.1. In the case of long-range interactions between slender bodies, we have argued that the application of such simple SSIP laws̃(d) is expected to be a good approximation, which will be confirmed by the quantitative analysis of Section 7.1.2.

Short-range volume interactions such as vdW and steric repulsion
In the following, a generic short-range volume interaction described by the point-pair potential law will be considered because it includes vdW interaction for exponent m = 6 (cf. Equation (4)) as well as steric repulsion as modeled by LJ for exponent m = 12 (cf. Equation (8)). As outlined already in the preceding section, only the regime of small separations is practically relevant in this case of short-range interactions and we neglect the effect of cross-section rotations throughout this article. At this point, we can thus return to the results for the disk-disk scenario obtained in literature on vdW interactions 33 and summarized in Table 2. In particular, we make use of expression (18) or rather the more general form (A12). The latter is valid for all power-law point pair interaction potentials with a general exponent m > 7∕2, that is, all interactions where the strength decays "fast enough." First, let us introduce the following abbreviation containing all constants in the lengthy expression: Using Equation (A12) in combination with the general SSIP approach (28) from Section 4.2, we directly obtain an expression for the total interaction potential of two deformable fibers in the case of short-range interactions: ‡ Disregarding rotations around its own axis, which are irrelevant here due to rotational symmetry, as mentioned above.
Here, the so-called gap g is the (scalar) surface-to-surface separation, that is, the beams' centerline curves r 1 (s 1 ) and r 2 (s 2 ) minus the two radii R i , as visualized in Figure 5. In general, the particle densities 1/2 may depend on the curve parameters s 1/2 , that is, vary along the fiber, without introducing any additional complexity at this point. For the sake of brevity, these arguments s 1/2 will be omitted in the remainder of this section. The variation of the interaction potential required to solve Equation (26) finally reads Here, we used the variation of the gap g, which is a well-known expression from the literature on macroscopic beam contact 12 and is identical to the variation of the separation of the beams' centerlines d to be used in Equation (36) because the cross-sections are assumed to be undeformable: Solving Equation (26) generally requires two further steps of discretization and subsequent linearization of this additional contribution Π m,ss to the total virtual work. The resulting expressions will be presented in Section 6.1.1 and Appendix B2, respectively. As discussed along with the general SSIP approach in Section 4.2, the remaining two nested 1D integrals are evaluated numerically, for example, by means of Gaussian quadrature. See Section 6.4 for details on this algorithmic aspect.
Remark on the regularization of the integrand. The inverse power law in the integrand of Equation (33) has a singularity for the limit of zero surface-to-surface separation g → 0. Consequently, a so-called regularization of the potential law is needed to numerically handle (the integration of) this term robustly as well as sufficiently accurate. This approach is well known from (beam) contact mechanics (see References 20,25,and 42) and will be further discussed and elaborated in Section 6.3. At the end of this section, we can conclude that we found specific, ready-to-use expressions for the free energy of interaction as well as virtual work of generic short-range interactions described via the SSIP approach. Thus, vdW interaction or steric exclusion of slender, deformable continua can now be modeled in an efficient manner, reducing the numerical integral to be evaluated from six to two dimensions. A detailed quantitative study of the approximation quality with regard to the assumptions discussed in the preceding Section 5.1 is the content of Section 7.1.1.
Note that speaking about the proposed SSIP laws in the remainder of this work refers to the fact that to the best of the authors' knowledge, these expressions have not been used as an interaction potential for beam cross-sections before and we thus propose to use specifically these expressions in this new context. In contrast to the general SSIP approach from Section 4, the specific expressions exemplarily chosen for the SSIP laws considered throughout this work are not original content of the present work and instead have been carefully selected from the literature and referenced accordingly.

Long-range surface interactions such as electrostatics
Having discussed short-range volume interactions, we now want to consider one example of long-range surface potentials.
Since electrostatic interaction is the prime example of surface potential interaction and at the same time of high interest for the application to biopolymers we have in mind, we will focus on this case throughout the following section and mostly speak of point charges as the elementary interaction partners. However, the required steps and formulas will be presented as general as possible in order to allow for a smooth future transfer to other applications.
Especially in this context, it is important to stress again that within this model the elementary interaction partners, that is, charges, must not redistribute within the bodies. Hence, only nonconducting materials can be modeled with the SSIP approach. This, however, covers our main purpose to model electrostatic interactions between biomacromolecules such as protein filaments and DNA because charges are not free to move therein.
According to the SSIP approach proposed in Section 4.2, we aim to use analytical expressions for the two inner integrals over the cross-section circumferences, while the integration along the two beam centerlines will be evaluated numerically (cf. Equation (27) in combination with the remark on surface interactions at the end of the corresponding section). As discussed in Section 5.1, the regime of large separations is the decisive one for beam-beam interactions in this case of long-range interactions and the SSIP law can be simplified in good approximation to depend only on the centroid separation d, which will be confirmed numerically in Section 7.1.2.
At this point, we again return to the expressions for the disk-disk interaction based on a generic point-pair potential Φ m (r) = k r −m , as derived in the literature on vdW interactions 33 and summarized in Appendix A2.1. In particular, the relation (A6) will be used, which is the same approximation used to derive Equation (19) that describes the practically rather irrelevant scenario of short-range vdW interactions in the regime of large separations. Note that in the context of electrostatics, this result is well known as the first term, that is, zeroth pole or monopole of the multipole expansion of the ring-shaped charge distribution on each of the disks' circumference, which represents the effect of the net charge of a (continuous) charge distribution and has no angular dependence (see Section 2.3.1). In simple terms, this monopole-monopole interaction means that the point pair interaction potential Φ(r) is evaluated only once for the distance between the centers of the distributions r = d = ||r 1 − r 2 || and weighted with the number of all point charges on the two circumferences of the circular cross-sections. The expression for the SSIP law to be used throughout this work would thus be exact for the scenario of the net charge of each cross-section concentrated at the centroid position (or distributed spherically symmetric around the centroid position). If the accuracy of the SSIP approach needs to be improved beyond the level resulting from this simplified SSIP law (see Section 7.1.2 for the analysis), one could simply include more terms from the multipole expansion of the (ring-shaped) charge distributions to the SSIP law, which would take the relative rotation of the cross-sections into account. Throughout this work and for the applications we have in mind, the simplified SSIP law, which is based on the monopole-monopole interaction of cross-sections, turns out to be an excellent approximation for the true electrostatic interaction law and we thus restrict ourselves to this variant.
Two nested 1D integrals over the beams' length dimensions then yield the two-body interaction potential for two fibers with arbitrary centerline shapes The surface (charge) densities j , j = 1, 2 have already been introduced in Equation (13). Particularly, for the case of electrostatics, the surface charge per unit length can be identified as j = 2 R j j , j = 1, 2, and is commonly referred to as linear charge density. Note, however, that Equation (35) holds for all long-range point pair potential laws Φ(r), for example, all power laws Φ m (r) = k r −m with m ≤ 3. In order to obtain the weak form of the continuous problem, the variation of this total interaction energy needs to be derived. This variational form can immediately be stated as as the consistent variation of the separation of the beams' centerlines d, which is well known from macroscopic beam contact formulations. 12 By inserting the generic (long-range) power law into Equation (36), we obtain the final expression for the variation of the two-body interaction energy of two deformable slender bodies The specific case of Coulombic surface interactions follows directly for m = 1 and k = C elstat (cf. Equation (3)). At this point, we have once again arrived at the sought-after contribution to the weak form (26) of the space-continuous problem. The steps of finite element discretization and linearization will again be presented later, in Section 6.1.2 and Appendix B3, respectively.
Remark on volume interactions. Note that there is no conceptual difference if long-range volume interactions were considered instead of the long-range surface interactions presented exemplarily in this section. The only difference lies in the constant prefactor c m,ls , which would read c m,ls = kmA 1 A 2 1 2 instead. Rather than the spatial distribution of the elementary interaction points in the volume or on the surface, it is the long-ranged nature of the interactions, which is important for the derivations in this section and allows the use of approximations for large separations (refer to the extensive discussion in Section 5.1).

FINITE ELEMENT DISCRETIZATION AND SELECTED ALGORITHMIC ASPECTS
Having discussed the space-continuous theory in Sections 4 and 5, we now turn to the step of spatial discretization by means of finite elements. Subsequently, the most important aspects of the required algorithmic framework will be presented briefly and discussed specifically in the light of the novel SSIP approach. This includes the applied regularization technique, multidimensional numerical integration, an analysis of the algorithmic complexity, as well as the topics of search for interaction partners and parallel computing.

Spatial discretization based on beam finite elements
As presented in Section 3, the centerline position r and the triad arise as the two primary fields of unknowns. Within the Simo-Reissner beam theory, both fields are uncorrelated and their discretization can hence be considered independently as follows. The Simo-Reissner finite beam element used throughout this work originates from References 4 and 5, although we apply a different centerline interpolation scheme here. We employ cubic Hermite polynomials based on nodal position vectorsd 1 ,d 2 and tangent vectorst 1 ,t 2 as the primary variables. See Reference 43 for a detailed discussion of Hermite centerline interpolation in the context of geometrically exact (Kirchhoff) beams and Reference 38 for the details on the Hermitian Simo-Reissner beam element that is used within this article. Applying this interpolation scheme results in the following discretized centerline geometry and variation: Here, all the degrees of freedom of one element relevant for the centerline interpolation, that is, nodal positionsd i and tangentst i , i = 1, 2 are collected in one vectord and H is the accordingly assembled matrix of shape functions, that is, Hermite polynomials H i d and H i t . The newly introduced element-local parameter ∈ [−1; 1] is biuniquely related to the arc-length parameter s ∈ [s ele,min ; s ele,max ] describing the very same physical domain of the beam as follows and the scalar factor defining this mapping between both infinite length measures is called the element Jacobian J( ): Our motivation to use Hermite interpolation is that it ensures C 1 -continuity, that is, a smooth geometry representation even across element boundaries. This property turned out to be crucial for the robustness of simulations in the context of macroscopic beam contact methods, 38 and is just as important if we include molecular interactions as proposed in this article. See Reference 44 for a comprehensive discussion of (non)smooth geometries and adhesive, molecular interactions using 2D solid elements. Note, however, that neither the SSIP approach proposed in Section 4 nor the specific expressions for the interaction free energy and the virtual work are limited to this Hermite interpolation scheme. In fact, all of the following discrete expressions will be equally valid for a large number of other beam formulations, where the discrete centerline geometry is defined by polynomial interpolation, which can generally be expressed in terms of the generic shape function matrix H introduced above.
Recall also that the proposed SSIP laws from Section 5 solely depend on the centerline curve description, that is, the rotation field does not appear in the additional contributions, and hence, its discretization is not relevant in the context of this work. It is, therefore, sufficient to apply the discretization scheme for the centerline field stated in Equation (39) to the expressions for the virtual work contributions Π ia presented in Section 5 and finally end up with the discrete element residual vectors r ia,1/2 . The latter need to be assembled into the global residual vector R as it is standard in the (nonlinear) finite element method. Note that the linearization of all the expressions presented in Section 6.1 is provided in Appendix B.

Short-range volume interactions such as vdW and steric repulsion
Discretization of the centerline curves according to Equation (39), that is, r j ≈ r h,j = H jdj and r T j ≈ r T h,j =d T j H T j , for both beam elements j = 1, 2 turns the space-continuous form (33) of the two-body virtual work contribution from molecular interactions Π m,ss into its discrete counterpart Refer to Equation (30) for the definition of the constant c m,ss . Note that Equation (41) only contributes to those scalar residuals associated with the centerline, that is, translational degrees of freedomd. This is a logical consequence of the fact that the SSIP law solely depends on the beams' centerline curves, as discussed in detail in Section 5. For the sake of brevity, the index "h," indicating all discrete quantities, will be omitted from here on since all the following quantities are considered discrete. In Equation (41), the discrete element residual vectors of the two interacting elements j = 1, 2 can finally be identified as See Section 6.4 for details on the numerical quadrature required to evaluate these expressions.

Long-range surface interactions such as electrostatics
In analogy to the previous section, we discretize Equation (38) and obtain the discrete-element residual vectors

Objectivity and conservation properties
It can be shown that the proposed SSIP approach from Section 4 in combination with the SSIP laws from Section 5 fulfills the essential mechanical properties of objectivity, global conservation of linear and angular momentum, as well as global conservation of energy. Due to the equivalent structure of the resulting space-discrete contributions, for example, Equation 41, as compared to the terms obtained in macroscopic beam contact formulations, we refer to the proof and detailed discussion of these important aspects in Appendix B of Reference 21. The fulfillment of conservation properties will furthermore be verified by means of the numerical examples in Sections 7.2 and 7.4.

Regularization of SSIP laws in the limit of zero separation
The singularity of inverse power laws for zero separation is a well-known pitfall when dealing with this kind of interaction laws. See p. 137 of Reference 2 for a discussion of this topic in the context of point-point LJ interaction as compared to a hard-sphere model. In numerical methods, one, therefore, typically applies a regularization that cures the singularity and ensures the robustness of the method. Sauer 44 gives an example for a regularized LJ force law between two half-spaces, where the force is linearly extrapolated below a certain separation, which is chosen as 1.05 times the equilibrium spacing of the two half spaces. Also, existing macroscopic beam contact formulations rely on the regularization of the seemingly instantaneous and infinite jump in the contact force when two macroscopic beams come into contact (see e.g. References 21 and 45). However, the SSIP laws derived for disk-disk vdW or LJ interaction from Section 5. 2 have not yet been considered in the literature. Note that LJ is the most general and challenging case considered in this work since strong adhesive forces compete with even stronger repulsive forces whenever two fibers are about to come into contact. To be more precise, it is not only the strength of these competing forces but also the high gradients in the force-distance relation that lead to a very stiff behavior of the governing partial differential equations. This alone places high demands on the nonlinear solver, which in combination with the already mentioned singularity at zero separation g = 0, and the fact that LJ interaction laws are not defined for configurations g < 0 where both fibers penetrate each other, makes it extremely demanding to solve the problem numerically.
The results and conclusions discussed throughout this section are mainly based on the extensive numerical peeling and pull-off experiment with two adhesive fibers, which is presented in the authors' recent contribution. 46 In the absence of a regularization, only the pragmatic yet effective approach of applying a very restrictive upper bound of the displacement increment per nonlinear iteration (see Appendix C for details) proved successful to solve for the quasi-static equilibrium configurations without occurrence of any invalid configuration g ≤ 0 for any integration point pair in any nonlinear iteration. It must be emphasized that even a single occurrence of the latter is fatal and aborts the simulation, such that the mentioned approach is the only way to compute a solution for the full LJ interaction law, which can, in turn, serve as a reference solution during the validation of the regularization to be proposed and applied. However, the mentioned approach severely deteriorates the convergence behavior and leads to a large number of nonlinear iterations per time step. Thus, the regularization to be proposed in this section is superior in two respects: it guarantees the avoidance of singular/undefined values and saves a factor of five in the number of iterations of the nonlinear solver.
Specifically, we apply a linear extrapolation of the total LJ force law below a certain separation g reg,LJ in a manner very similar to Reference 44 with the only difference that it is applied to the length-specific disk-disk force law instead of the force law between two half spaces. Linear extrapolation means that the original expression (m − 7∕2) c m,ss g −m+ 5 2 in Equations (42) and (43) is replaced by a linear equation a g + b in the gap g for all g < g reg,LJ . The two constants a and b are determined from the requirements that the force value as well as the first derivative of the original and the linear expression are identical for the regularization separation g = g reg,LJ . Figure 7 shows both the original (blue) and the regularized (red) LJ disk-disk force law as a function of the smallest surface separation g.
The numerical experiment of adhesive fibers studied in Reference 46 reveals that this regularization yields the already mentioned great enhancement in terms of robustness as well as efficiency without any change in the system response. As shown in the comparison of the force-displacement curves therein, the results obtained with the full LJ and with the regularized LJ force law do, indeed, coincide down to machine precision. This is reasonable and expected because we chose a regularization parameter g reg,LJ ≤ g LJ,eq,cyl||cyl that is smaller than any separation value g occurring anywhere in the system in any converged equilibrium state. Thus, the solution never "sees" the modification to the vdW force law in the interval g < g reg,LJ and the results are identical. However, since during the nonlinear iterations also nonequilibrium F I G U R E 7 Comparison of regularized (red) and full (blue) LJ disk-disk force law. Here, g reg,LJ = g LJ,eq,disk||disk is shown exemplarily. § Note that the same idea has previously been applied in the context of LJ interaction between half spaces 44 instead of disks as considered here configurations with g < g reg,LJ occur, the nonlinear solution procedure is influenced in an extremely positive way, leading to an overall saving of a factor of five in the number of nonlinear iterations as compared to the full LJ interaction without any regularization. For more details on the comparison, including all parameter values, we kindly refer the reader to Reference 46.

Numerical evaluation of n-dimensional integrals of intermolecular potential laws
Generally, we use n nested loops of a 1D Gauss-Legendre quadrature scheme that is the well-established and de-facto standard method in nonlinear finite element frameworks and has been used also in the previous publications in the context of molecular interactions. 22,25 Due to the strong nonlinearity, that is, high gradients of the power laws, a large number of quadrature points are required in each dimension to achieve sufficient accuracy. This effect is most critical for high exponents of the potential law, that is, vdW and steric interactions, and small separations of the interacting bodies. We thus implemented the possibility to subdivide the domain of a finite element into n IS integration segments and apply an n GP -point Gauss rule on each of them in order to achieve sufficient density of quadrature points in every case.

Algorithm complexity
Multidimensional numerical integration of the intermolecular potential laws as discussed above turns out to be the crucial factor in terms of efficiency. For the following analysis of efficiency, we consider the associated algorithmic complexity. Generally, all possible pairs of elements need to be evaluated, which has  ( n 2 ele ) complexity. Let us assume we apply a total of n GP,tot,ele-length integration points along the element length and n GP,tot,transverse integration points in the transversal, that is, cross-sectional in-plane directions. Thus, the complexity of an approach based on full 6D numerical integration over the 3D volumes of the two interacting bodies (cf. Equation (12)  .
In contrast to that, the novel SSIP approach proposed in Section 4.2 reduces the dimensionality of numerical integration from six to two (cf. Equation (28)) and thus yields complexity. The resulting difference between both clearly depends on the problem size, type of interaction, and other factors. To get an impression, typical numbers for the total number of quadrature points in transverse dimensions based on the numerical examples of Section 7 are given as n GP,tot,transverse = 10 … 100. The gain in efficiency thus easily exceeds a factor of 10 4 and can be as large as a factor of 10 8 . In addition to this tremendous saving from the inherent algorithmic complexity, the power law integrand has a smaller exponent due to the preliminary analytical integration in the case of the SSIP approach. This, in turn, allows for a smaller number of integration points n ele ⋅ n GP,tot,ele-length for each of the two nested 1D integrations along the centerline, given the same level of accuracy. To give an example, the vdW interaction force scales with an exponent of −7 if formulated for two points (cf. Equation (4)) as compared to an exponent of −7∕2 for two circular cross-sections (cf. Equation (33) for m = 6). This makes another significant difference, especially if very small separations-as typically observed for contacting bodies-are considered. The combination of high dimensionality and strong nonlinearity of the integrand renders the direct approach of 6D numerical quadrature to evaluate Equation (12) unfeasible for basically any problem of practical relevance. In fact, even a single evaluation of the vdW potential of two straight cylinders to serve as a reference solution turned out to be too computationally costly below some critical, small separation. See Section 7.1.1 for details on this numerical example. Note that although there might be more elaborate numerical quadrature schemes for these challenging integrands consisting of rational functions (see Equation 47), the basic problem and the conclusions drawn from this comparison of algorithmic complexities remain the same. These cost estimates based on the theoretical algorithm complexity and the experience from rather small academic examples considered in Section 7 show that the SSIP approach, indeed, makes the difference between feasible and intractable computational problems. This directly translates to the applicability to complex biopolymer as well as synthetic fibrous systems that we have in mind and thus significantly extends the range of (research) questions that are accessible by means of numerical simulation.

Search algorithm and parallel computing
In order to find the relevant pairs of interaction partners, the same search algorithms as in the case of macroscopic contact (between beams or 3D solids) may be applied; however, the obvious difference lies in the search radius. For contact algorithms, a very small search radius covering the immediate surrounding of a considered body is sufficient, whereas for molecular interactions, the search radius depends on the type of interaction and must be at least as large as the so-called cut-off radius. Only at separations beyond the cut-off radius, the energy contributions from a particular interaction are assumed to be small enough to neglect them. Depending on the interaction potential and partners, the range and thus cut-off radius can be considerably large, which underlines the importance of an efficient search algorithm. In the scope of this work, a so-called bucket search strategy has been used, which divides the simulation domain uniformly into a number of cells or buckets and assigns all nodes and elements to these cells to later determine spatially proximate pairs of elements based on the content of neighboring cells. This leads to an algorithmic complexity of (n ele ) and the search thus turned out to be insignificant in terms of computational cost as compared to the evaluation of pair interactions as discussed in the preceding section. See Reference 48 for an overview of search algorithms in the context of computational contact mechanics.
To speed up simulations of large systems, parallel computing is a well-established strategy of ever-increasing importance. Key to this concept is the ability to partition the problem such that an independent and thus simultaneous computation on several processors is enabled. In our framework, this partitioning is based on the same bucket strategy that handles the search for interaction partners. Regarding the evaluation of interaction forces, a pair (or set) of interacting beam elements is assigned to the processor that owns and thus already evaluates the internal and external force contribution of the involved elements. At processor boundaries, that is, if the two interacting elements are owned by different processors, one processor is chosen to evaluate the interaction forces and the required information such as the element state vector of the element owned by the other processor is communicated beforehand. Upon successful evaluation of the element pair interaction, the resulting contribution to the element residual vector and stiffness matrix is again communicated for the element whose owning processor was not responsible for the pair evaluation.

NUMERICAL EXAMPLES
The set of numerical examples studied in this section aims to verify the effectiveness, accuracy, and robustness of the proposed SSIP approach and the corresponding SSIP laws as a computational model for either steric repulsion, electrostatic or vdW adhesion, and also a combination of those. Supplementary information on the code framework and the algorithms used for the simulations is provided in Appendix C.

Verification of the simplified SSIP laws using the examples of two disks and two cylinders
As a follow-up to the general discussion of using simplified SSIP laws in Section 5.1 and the proposal of specific closed-form analytic expressions in Sections 5.2 and 5.3, this section aims to analyze the accuracy in a quantitative manner. The minimal examples of two disks or two cylinders are considered in order to allow for a clear and sound analysis of either the isolated SSIP laws or its use within the general SSIP approach to modeling beam-beam interactions, respectively.

Verification for short-range volume interactions such as vdW and steric repulsion
Throughout this section, we consider the example of vdW interaction, but analogous results are expected for steric interaction or any other short-range volume interaction. Specifically, we will focus on the approximation quality of the proposed SSIP law in Section 5.2, which is based on the assumptions and resulting simplifications discussed in-depth in Section 5.1. Recall that, beside the obviously most important surface-to-surface separation g, the rotation of the cross-sections around the closest point (quantified by the angle enclosed by their tangent vectors) and potentially also the rotation components 1 , 2 (see Figure 6A) have been identified as relevant degrees of freedom, yet are neglected in the simplified SSIP law proposed in Section 5.2. The influence of these factors, separation, and rotation, on the approximation quality will thus be analyzed numerically in the following.
Recall also from the discussions in Sections 5.1 and 5.2 that only the regime of small separations will be of practical relevance in the case of short-range interactions considered here. However, we include the regime of large separations in the following analyses, mainly because it will be interesting to see the transition from small to large separations and confirm that potential values, indeed, drop by several orders of magnitude as compared to the regime of small separations. Moreover, it is a question of theoretical interest and has been considered in the literature on vdW interactions. 33 This regime of large separations can be treated without any additional effort as described for the case of long-range interactions in Section 5.3 (where this regime is the decisive one) if we take into account the corresponding remark on volume interactions.
As presented in Section 2.3.2, analytical solutions for the special cases of parallel and perpendicular cylinders, for the regime of small and large separations, respectively, can be found in the literature 2,3 and thus serve as reference solutions in this section. To the best of the authors' knowledge, no analytical reference solution has yet been reported for the intermediate regime in between the limits of large and small separations. Another source for reference solutions is the full numerical integration of the point pair potential over the volume of the interacting bodies; however, it is limited due to the tremendous computational cost. Only a combination of both analytical and numerical reference solutions thus allows for a sound verification of the novel SSIP approach and the proposed SSIP laws.
In the following analyses, either the SSIP, that is, the interaction potential per unit length squared̃of a pair of circular cross-sections, the interaction potential per unit length̃of a pair of parallel cylinders or the interaction potential Π of a pair of perpendicular cylinders will be plotted as a function of the dimensionless surface-to-surface separation g∕R, respectively. For simplicity, the radii of the beams are set to R 1 = R 2 =∶ R = 1. Figure 8A shows the SSIP̃of two disks in parallel orientation, that is, their normal vectors are parallel with mutual angle = 0, as a function of the normalized separation g∕R. This is the simplest geometrical configuration and forms the basis of the proposed SSIP laws from Section 5. We thus begin our analysis with the verification of the used analytical solutions in the limit of small (green line, cf. Equation (18)) and large (red line, cf. Equation (19)) separations by means of a numerical reference solution (black dashed line with diamonds) obtained from 4D numerical integration of the point pair potential law (4). Figure 8A confirms that both analytical expressions match the numerical reference solution perfectly well in the limit of large and small separations, respectively. As predicted, the interaction potential of two circular disks follows a power law with (negative) exponent 2.5 for small and 6 for large separations. ¶ Note that all plots in this section are normalized with respect to the length scale R and the energy scale 1 2 C vdW . It is remarkable that the obtained values span several ¶ Note that in the double logarithmic plot, a power law with exponent m is a linear function with slope m. Regarding the full range of separations, one may ask where either of the two expressions may be used given a maximal tolerable error threshold. As can be concluded from Figure 8A, the resulting error is small for separations g∕R < 0.1 with a relative error below 8% and g∕R > 10 with a relative error below 7%. In the region of intermediate separations, the analytical solution for small separations seems to yield an upper bound, whereas the one for large separations seems to yield a lower bound for the interaction potential.

Parallel disks and cylinders
Let us have a look at the efficiency gain from using the analytical solutions. The numerical reference solution requires the evaluation of a 4D integral over both cross-sectional areas for a given separation g and has been carried out in polar coordinates. Assuming Gaussian quadrature with the same number of Gauss points n GP,tot,transverse in radial and circumferential dimension and for both cross-sections, this requires a total of (n GP,tot,transverse ) 4 function evaluations. By contrast, the analytical expressions for the large and small separation limit, respectively, require only one function evaluation. This significant gain in efficiency is most pronounced for small separations, where the number of required Gauss points increases drastically due to the high gradient of the power law that needs to be resolved (see Section 6.4 for details). If the number of Gauss points is not sufficient, this leads to the so-called underintegration and we observed that the obtained curve of the numerical reference solution erroneously flattens (because the contribution of the closest-point pair is not captured) or becomes steeper (because the contribution of the closest-point pair is overrated).
For these reasons, the computation of an accurate numerical reference solution shown in Figure 8A requires quite some effort. The integration domains were subdivided into integration sectors (see Figure 8B) in order to further increase the Gauss point density. However, even in this planar disk-disk scenario requiring only 4D integration, we reached a minimal separation of g∕R ≈ 5 × 10 −3 , below which the affordable number of Gauss points was not sufficient to correctly evaluate the SSIP̃by means of full numerical integration. # For these very small separations, only the exact analytical dimensional reduction from 4D to 2D according to Langbein (cf. Reference 33 and Equation (A7)) allowed to compute an accurate numerical reference solution. The analytical solutions for the disk-disk interaction potential from Equation (18) (and Equation (19)), used as SSIP law in Equation (33) (and Equation (36)  2 nm to be a typical value for contacting solid bodies and states that accurate numerical integration thus is the most challenging and, in fact, limiting factor for the numerical methods based on intersurface potentials. As a reference value for the applications we have in mind, the fiber radius R varies from several nanometers for DNA to millimeters for synthetic polymer fibers, resulting in a potentially very small normalized separation g∕R. An example for the simulation of adhesive fibers in contact can be found in the authors' recent contribution, 46 which studies the peeling and pull-off behavior of two fibers attracting each other either via vdW or electrostatic forces.
As a next step, the interaction potential per unit length̃of two parallel straight beams is considered. The length of the beams is chosen sufficiently high such that it has no perceptible influence on the results and meets the assumption of infinitely long cylinders made to derive the analytical reference solution from Reference 33. Accordingly, a slenderness ratio = l∕R = 50 is used in the regime of small separations, whereas = l∕R = 1000 is used for large separations. Based on the experience from the disk-disk scenario, it is not surprising that the full 6D numerical integration in this case exceeds the affordable computational resources by orders of magnitude and thus cannot serve as a reliable reference solution. In fact, we were not able to reproduce the theoretically predicted power law scaling in the regime of small separations despite using a number of Gauss points that led to computation times of several days. However, instead of the numerical reference solution, the analytical solution for infinitely long cylinders in the limit of very small (black dashed line, cf. Equation (14)) and very large separations (blue dashed line, cf. Equation (15)) serves as a reference in Figure 9. Note that as compared to the case of two circular disks, the exponent of the power laws and thus the slope of the curves drops by one due to the integration over both cylinders' length dimension.
Interestingly, the SSIP approach using the simplified SSIP law from Section 5.2 (green line with crosses) does not yield the correct scaling behavior even in this case of parallel cylinders. This confirms the concerns from Section 5.1 that the simplified SSIP law neglecting any relative rotations of the cross-sections deteriorates the accuracy of the approach in the case of short-ranged interactions in the regime of small separations. Due to this specific scenario of parallel cylinders, this deterioration can be attributed solely to the rotation components 1/2 (see Figure 6A) since the included angle of the cross-section normal vectors is zero for every of the infinitely many pairs of cross-sections. Despite the correct trend of the resulting interaction potential as an inverse power law of the surface separation, it must thus be stated that the simplified SSIP law overestimates the strength of interaction and that the error increases with decreasing separation. ‖ In the regime of large separations, however, the results for the SSIP approach (red line with circles) perfectly match the analytical reference solution (blue dashed line). This confirms the hypothesis from Section 5.1 that the relative rotation of cross-sections is negligible in this regime and a high accuracy can be achieved with the simplified SSIP law. Although being of little practical importance here due to the negligible absolute values, this is a first numerical evidence for the validity of the SSIP approach in general and its high accuracy even in combination with simplified SSIP laws in the particular case of long-range interactions to be considered in Section 7.1.2.
‖ Note that the numerical integration error has been ruled out as cause for this behavior by choosing a high number of Gauss points n GP,tot,ele-length = 2 × 50 = 100 for each of the 64 elements used to discretize each cylinder. A further increase of n GP by a factor of five does not change the results using double precision.

Perpendicular disks and cylinders
Up until now, we have only discussed the situation of parallel orientation of disks and cylinders. In the following, the accuracy of the simplified SSIP laws as well as the SSIP approach for twisted configurations will be analyzed by considering the most extreme configuration of perpendicular disks and cylinders. Again, computing a reference solution by means of full numerical integration was only affordable for the 4D case of two disks. The results for perpendicular disks shown in Figure 10A confirm that there is no difference between perpendicular and parallel orientations for large separations and the scaling behavior of the numerical reference solution (black dashed line with diamonds) with exponent 6 is met by the simplified SSIP law (red line, cf. Equation (19)). On the other hand, there is a remarkable difference in the scaling behavior for small separations. While the interaction potential of two parallel disks, which is the underlying assumption of the proposed SSIP law (green line, cf. Equation (18)), follows an inverse 2.5 power law, the numerical reference solution (black dashed line with diamonds) suggests that this behavior changes for perpendicular disks to an inverse 2 power law. This time, the difference in results can be attributed to the relative rotation , that is, the angle included by the cross-section normal vectors and again the error of the proposed simplified SSIP law increases with decreasing separation.
Finally, the scenario of perpendicular cylinders is considered and Figure 10B shows the total interaction potential Π as a function of the normalized smallest surface separation g∕R. As discussed before, the computational cost of the full 6D numerical integration is too high to compute a reliable reference solution in the case of two 3D bodies and we resort to the analytical solutions for the limits of very small and very large separations, respectively. Note that in contrast to the case of infinitely long parallel cylinders (cf. Figure 9), the total interaction potential of infinitely long perpendicular cylinders is finite and the result thus has dimensions of energy instead of energy per length. Perpendicular cylinders are worth to consider because they trigger both of the sources of errors that have been analyzed individually so far, neglecting relative rotations as well as 1/2 in the simplified SSIP law. In short, the resulting accuracy is similar as for either perpendicular disks or parallel cylinders. In the decisive regime of small separations, the SSIP approach based on the simplified SSIP law (green line with crosses) fails to reproduce the correct scaling behavior of the analytical reference solution (black dashed line, cf. Equation (16)), whereas in the regime of large separations, the SSIP approach based on the simplified SSIP law (red line with circles) perfectly matches the analytical reference solution (blue dashed line, cf. Equation (17)).

Conclusions
First, this section reveals that full 6D numerical integration to compute the total interaction potential of slender continua is by orders of magnitude too expensive and can not reasonably be used as a numerical reference solution even in minimal examples of one pair of cylinders. At most, 4D numerical integration required for disk-disk interactions allows computing numerical reference solutions for the intermediate regime of separations where no analytical solutions are known. This underlines the importance of reducing the dimensionality of numerical integration to 2D as achieved by the proposed SSIP approach in order to enable the simulation of large systems as well as a large number of time steps.
Second, the thorough analysis of the accuracy resulting from using the proposed simplified SSIP law, neglecting the cross-section rotations, reveals that one has to distinguish between the regime of small and large separations. In the decisive regime of small separations, we find that the scaling behavior deviates from the analytical prediction for perpendicular disks and parallel as well as perpendicular cylinders and that the resulting error increases with decreasing separation. A remedy of this limitation could be a calibration, that is, a scaling of the prefactor k in the simple SSIP law, to fit a given reference solution within a small range of separations (eg, around the equilibrium distance of the LJ potential). In the authors' recent contribution, 46 this pragmatic procedure is shown to reproduce the global system response very well. Still, it would be valuable to include the relative rotations of the cross-sections in the applied SSIP law to obtain the correct asymptotic scaling behavior. To the best of the authors' knowledge, no analytical closed-form expression has been published yet and the like is far from trivial to derive. Therefore, we leave this as a promising enhancement of the novel approach, which we are currently working on and will address in a future publication. As mentioned before, the regime of large separations is of little practical relevance in the case of short-range volume interactions; however, it is of some theoretical interest and the corresponding findings and conclusions of this regime will hold true also for long-range interactions such as electrostatics to be considered in the following section. Here, the results are in excellent agreement with the theoretically predicted power laws for parallel as well as perpendicular disks and cylinders.

Verification for long-range surface interactions such as electrostatics
Turning to long-range interactions, again parallel and perpendicular disks and cylinders will be considered in order to analyze the accuracy of the simplified SSIP law from Section 5.3 both individually as well as applied within the general SSIP approach proposed in Section 4.2. As before, Coulombic surface interactions are chosen as specific example; however, the conclusions are expected to hold true also for other types of long-range interactions. As compared to the preceding section, the computation of a numerical reference solution simplifies mainly due to the reduction from volume to surface interactions but also due to the smaller gradient values, which need to be resolved in the regime of small separations, thus requiring less integration points. This allows for a verification by means of a numerical reference solution also in the case of cylinder-cylinder interaction. Figure 11 shows the results for the simplified SSIP law obtained from the monopole-monopole interaction of two disk-shaped cross-sections in Section 5.3 (red line) and a numerical reference solution (black dashed line with diamonds). As expected, the proposed SSIP law excellently matches the reference solution in the regime of large separations, both for the parallel as well as perpendicular configuration. In both cases, the relative error is below 7% already for g∕R = 1. The most important and remarkable result of this section, however, is the following. The inevitable error of the simplified SSIP law in the regime of small separations does not carry over to beam-beam interactions as shown in Figure 12. For both parallel as well as perpendicular cylinders, the results from the SSIP approach using this simplified SSIP law from Section 5.3 (red line with crosses) agree very well with the numerical reference solution (black dashed line with diamonds) over the entire range of separations. This confirms the theoretical considerations from Section 5.1 arguing that the beam-beam interaction will be dominated by the large number of section pairs with large separations, which outweigh the contributions of the few section pairs with smallest separations.
A closer look reveals that the relative error for the parallel cylinders is below 0.3% even for the smallest separation g∕R = 10 −3 considered here. For the presumably worst case of perpendicular cylinders, this deviation is even smaller with a relative error of 0.03%, which can be explained by the following two reasons. First, a comparison of Figure 11A with Figure 11B reveals that the accuracy of the simplified SSIP law in the regime of small separations is higher for perpendicular orientation, which can be regarded as a happy coincidence. Second, the large majority of all section pairs has a larger separation (which according to Figure 11 is the regime of higher accuracy) as compared to the case of parallel cylinders.
Note that unlike in the case of short-range interactions, here the total interaction potential is also considered for the parallel cylinders. Due to the long range of interactions, the interaction potential per length depends on the length of the cylinders and is thus no representative quantity. For Figure 12, a slenderness ratio of = L∕R = 50 is chosen exemplarily. The two nested 1D integrals along the cylinder length dimensions are evaluated using n GP,ele-length = 5 Gauss points for each of the 64 elements used to discretize each cylinder. In addition, n GP,circ = 8 × 32 = 256 Gauss points over the circumference of each disk are used to compute the numerical reference solution. In all cases, it has been verified that the numerical integration error does not influence the results noticeably.
At this point, recall from Section 5.3 that the accuracy of the applied SSIP law can still be increased whenever deemed necessary by including more terms of the multipole expansion of the cross-sections. However, because the results of this section show a high level of accuracy and the resulting simplification is significant, the simplified SSIP law seems to be the best compromise for our purposes. To conclude this section, it can thus be stated that the novel SSIP approach as proposed in Section 4.2 in combination with the simplified SSIP law from Section 5.3 is a simple, efficient, and accurate computational model for long-range interactions of slender fibers. In the following, it will be applied to first numerical examples of deformable slender fibers in Sections 7.3 and 7.4.

Repulsive steric interaction between two contacting beams
This numerical example aims to demonstrate the general ability of our proposed method to preclude penetration of two slender bodies that come into contact under arbitrary mutual orientation in 3D. No adhesive forces are considered in this example. The setup is inspired by Example 1 in Reference 21 where the macroscopic, so-called all-angle beam contact (ABC) formulation is used to account for the nonpenetrability constraint. Here, we model the contact interaction based on the repulsive part of the LJ interaction potential (8). More specifically, we apply the novel SSIP approach as proposed in Section 4.2 in combination with the SSIP law proposed in Section 5.2. The parameter specifying the strength of repulsion is set to k 1 2 = 10 −16 . To be consistent throughout this article, we apply Hermitian Simo-Reissner beam elements instead of the torsion-free Kirchhoff elements used in Reference 21. As compared to the original example, this requires us to replace the hinged support of the upper beam by clamped end Dirichlet boundary conditions in order to eliminate all rigid body modes in this quasi-static example. The same number of three finite elements for the upper, deformable beam, and one element for the lower, rigid beam is used. A sequence of the resulting simulation snapshots is shown in Figure 13. As expected, the two beams do not penetrate each other in any of the various mutual orientations throughout the simulation. Figure 14 visualizes the contact force distributions ** in the most interesting time span before the beams reach the parallel orientation at time t = 2.0. The force distribution quickly changes from a pointlike force for large mutual angles to a broad distributed load for parallel beam axes. Note also that the line load has a 3D shape where the out-of-plane component decreases with decreasing mutual angle until both beam axes and thus also the line loads lie in one plane at t = 2.0. Another remarkable result is the symmetry between the line loads on both fibers. It nicely confirms that the novel approach indeed fulfills the **More precisely, the vectorial line load with dimensions of force per unit length is visualized as an arrow at each integration point. The force resultant, therefore, equals the integral over the contour curve defined by the arrows' tips (ie, the area under this curve), and not the vector sum of all arrows shown. This is important to understand because the number of visible arrows per unit length depends on the discretization and is thus higher for the upper, deformable beam. expected local equilibrium of interaction forces in good approximation. In contrast to existing, macroscopic formulations for beam contact, this is not postulated a priori in our approach and hence is a valuable verification at this point. See Reference 25 for a comprehensive discussion of this important topic in the context of contact between 3D solids described by intersurface potentials. The global equilibrium of contact forces, on the other hand, is fulfilled exactly, as can be concluded from the global conservation of linear momentum that can be shown analytically as outlined in Section 6.2. In this numerical example, we found that the sum of all reaction forces in either of the spatial dimensions is, indeed, zero with a maximal residuum of 10 −10 throughout all simulations considered here, which confirms the statement numerically. Figure 15 shows the resulting vertical reaction force as well as the interaction potential over time. Due to the inverse-12 power law and the extremely small separations of the interacting bodies, the numerical integration of the disk-disk interaction forces is very challenging and we studied the influence of the number of Gauss points. For this purpose, the number of integration segments per element with five Gauss points each is set to 30, 50, or 64, respectively. Interestingly, the interaction potential shown in Figure 15B seems to be more sensitive with respect to the integration error than the vertical reaction force shown in Figure 15A despite the fact that the latter has a higher inverse power law exponent. Presumably, this is due to the fact that the reaction force is dominated by the bending deformation of the beams. For reference, the reaction force obtained by using the macroscopic ABC formulation is shown as well and is in excellent agreement with the one resulting from the repulsive part of the LJ interaction potential.
A more comprehensive comparison of this novel SSIP approach to model contact between beams based on (the repulsive part of) the molecular LJ interaction and existing, macroscopic formulations based on heuristic penalty force laws is a highly interesting subject that is worth to investigate in the future.

Two initially straight, deformable fibers carrying opposite surface charge
The following example consists of two initially straight and parallel, deformable fibers that attract each other due to their surface charge of opposite sign. Its setup is kept as simple as possible to allow for an isolated and clear analysis of the physical effects as well as the main characteristics of the proposed SSIP approach. In a first step presented here, the interplay of elasticity and electrostatic attraction in the regime of large separations is studied. In addition, the authors' recent contribution 46 considers the scenario of separating these adhesive fibers starting from initial contact and studies a variety of physical effects and influences in depth, which would go beyond the scope of this work. In this numerical example, we are interested in the static equilibrium configurations for varying attractive strength. As shown in Figure 16A, two straight beams of length l = 5 are aligned with the global y-axis at an interaxis separation d = 5. Both are simply supported and restricted to move only within the xy-plane and rotate only around the global z-axis. The beams have a circular cross-section with radius R = 0.02, which results in a slenderness ratio of = 250. Cross-section area, area moments of inertia, and shear correction factor are computed using the standard formula for a circle. A hyperelastic material law with Young's modulus E = 10 5 and Poisson's ratio = 0.3 is applied. In terms of spatial discretization, we use five Hermitian Simo-Reissner beam elements per fiber (see Reference 38 for details on this element formulation). Electrostatic interaction is modeled via the SSIP approach as presented in Section 4.2 and applied to long-range Coulomb interactions in Section 5.3. Both beams are nonconducting with a constant surface charge density of 1 = 1.0 and 2 = −1.0, respectively. For simplicity, we vary the prefactor k of the underlying Coulomb law Φ(r) = k r −1 to vary the strength of attraction. However, as becomes clear from Equation (38), this is equivalent to a variation of surface charge densities because in our case the product of these quantities is a constant prefactor in all relevant equations. In order to evaluate the electrostatic force and stiffness contributions, Gaussian quadrature with two integration segments per element and ten Gauss points per integration segment is applied. This turns out to be fine enough to not change the presented results perceptibly. More precisely, the difference in the displacement of the beam midpoint for n GP = (2 × 10) 2 as compared to (2 × 32) 2 is below 10 −8 . No cut-off radius is applied here, that is, the contributions of all Gauss point pairs are evaluated and included. Figure 16B finally shows the resulting static equilibrium configurations for different levels of attractive strength. As expected, the beams are increasingly deflected and pulled toward each other if the prefactor of the applied Coulomb law k, and thus the attractive strength is increased. Like the problem definition, also all the solutions are perfectly symmetric with respect to the vertical axis of symmetry located at x = d∕2. Moreover, the centerline curves of each individual solution show a horizontal axis of symmetry defined by the position of the two beam midpoints in the respective deformed state. Consequently, the vertical force components in the system cancel and the vertical reaction forces vanish. This also becomes clear when looking at the visualization of the resulting electrostatic forces as shown exemplarily for k = 1.0 in Figure 17A. In addition, the forces acting on the Gauss point of one beam caused by the interaction with one finite element of the other beam are visualized individually in Figure 17B. This representation illustrates the nature of the SSIP approach, which is based on two nested 1D numerical integrals that are evaluated element pairwise. Accordingly, we can identify five force contributions at each Gauss point, one for each of the five beam elements on the opposing fiber. As expected, the magnitude of these individual forces decays with the distance and the contributions of the closest element pair shown in an isolated manner in Figure 17C constitute the largest part of the total electrostatic load on the beams and are clearly larger than the contributions of the next-nearest element pair shown in Figure 17D. However, the comparatively long range of electrostatic forces yields a smooth force distribution along the centerlines and we can identify nonzero force contributions even at the most distant Gauss points right next to the supports in Figure 17A. As mentioned above, a quantitative analysis of the resulting horizontal reaction forces is presented in Reference 46.
To conclude this example of two charged, attractive beams, we briefly look at the nonlinear solver. Newton's method without any adaptations is used here to allow for a clear and meaningful analysis of nonlinear convergence behavior. The solutions for k ≤ 0.4 can be found within one load step, which is a remarkable result given the resulting large deflection of the beams shown in Figure 16B and the strong nonlinear nature of the system. For stronger attractive forces, the strength of electrostatic attraction was ramped up in up to 10 equal steps Δk = 0.1. As convergence criteria, we enforced both for the Euclidean norm of the residual vector ||R|| < 10 −10 and for the norm of the iterative displacement update vector ||ΔX|| < 10 −8 . In fact, this combination leads to ||R|| < 10 −12 in almost all equilibrium configurations shown here.

Two charged deformable fibers dynamically snap into contact
Due to the high gradients in the inverse power laws, molecular interactions give rise to highly dynamic systems. This is a first and simple example for a dynamic system consisting of two oppositely charged fibers with a hinged support at one end each, which will snap into contact. In the initial configuration shown in Figure 18A, the straight fibers include an angle of 45 • and their axes are separated by 5R in the out-of-plane direction z. With a cross-section radius R = 0.02 and length l set to l = 5, they have a high slenderness ratio of = 250 and 354. Each of the fibers is discretized by 10 Hermitian Simo-Reissner beam elements and the material parameters are chosen to be E = 10 5 , = 0.3, and = 10 −3 . The fibers carry a constant, opposite surface charge 1/2 = ±1.0 and interact via the Coulomb potential law stated in Equation (3) with the prefactor set to C elstat = 0.4. In order to start from a stress-free initial configuration, the charge of one of the fibers is ramped up linearly within the first 100 time steps. We apply the SSIP approach as proposed in Section 4.2 and applied to Coulomb interactions in Section 5.3. A total of five integration segments per element with 10 Gauss points each is used to evaluate these electrostatic contributions. The contact interaction between the fibers is modeled by the line contact formulation proposed in Reference 20, using a penalty parameter = 10 3 and 20 integration segments per element with 5 Gauss points each for numerical integration. An undetected crossing of the fiber axes is prevented by applying the modified Newton method limiting the maximal displacement increment per nonlinear iteration to R∕2 (see Appendix C for details).
In terms of temporal discretization, we apply the Generalized-alpha scheme for Lie groups as proposed in Reference 49 and set the spectral radius at infinite frequencies to ∞ = 0.9 for small numerical damping. A small time step size of Δt = 10 −5 is applied to account for the highly dynamic behavior of this system. Figure 19 shows a sequence of simulation snapshots where the electrostatic forces on both fibers are visualized as green arrows. We observe a large variety of mutual orientations of the two fibers and a strong coupling of adhesive, repulsive, and elastic forces that demonstrate the effectiveness and robustness of the proposed SSIP approach. Most importantly, we see that the total system energy is preserved with very little deviation of ±2%, as shown in Figure 18B. Note that the negative energy values result from defining the zero level of the interaction potential at infinite separation, as described in Section 2.1. Based on this

F I G U R E 19
Sequence of simulation snapshots. Electrostatic forces acting on both fibers shown in green numerical example, we can thus conclude that the novel SSIP approach proves to be effective as well as robust in a highly dynamic example with arbitrary mutual orientations of the fibers in three dimensions.

CONCLUSIONS AND OUTLOOK
This contribution proposes the first 3D beam-beam interaction model for molecular interactions such as electrostatic, vdW, or repulsive steric forces between curved slender fibers undergoing large deformations. While the general model is not restricted to a specific beam formulation, in the present work, it is combined with the geometrically exact beam theory and discretized via the finite element method. A direct evaluation of the total interaction potential for general 3D bodies requires the integration of contributions from molecule or charge distributions over the volumes of the interaction partners, leading to a 6D integral (two nested 3D integrals) that has to be solved numerically. The central idea of our novel approach is to formulate reduced interaction laws for the resultant interaction potential between a pair of cross-sections of two slender fibers such that only the two 1D integrals along the fibers' length directions have to be solved numerically. This SSIP approach, therefore, reduces the dimensionality of the required numerical integration from 6D to 2D and yields a significant gain in efficiency, which only enables the simulation of relevant time and length scales for many practical applications.
Being the key to this SSIP approach, the analytical derivation of the specific SSIP laws is based on careful consideration of the characteristics of the different types of molecular interactions, most importantly their point pair potential law and the range of the interaction. In a first step, the most generic form of the SSIP law, which is valid for arbitrary shapes of cross-sections and inhomogeneous distributions of interacting points (eg, atoms or charges) within the cross-sections has been presented in Section 4 before we turn to specific SSIP law expressions for dedicated applications in Section 5. Considering the practically relevant case of homogeneous, disk-shaped cross-sections, Section 5.1 discusses the advantages as well as limitations of possible simplifications of the SSIP law, finally leading to pleasantly simple SSIP laws that solely depend on the scalar separation of the cross-section centroids and in particular neglect the cross-sections' mutual orientation. In this manner, simple SSIP laws for short-range volume interactions such as vdW or steric interactions (Section 5.2) and for long-range surface interactions such as Coulomb interactions (Section 5.3) are obtained, which have been used and further analyzed in the remainder of the article. We would like to stress that postulating the general structure of the SSIP law and fitting the free parameters to experimental data is one of the promising alternatives to the strategy of analytical derivation of the SSIP law as applied in this article.
It is also important to emphasize that the general SSIP approach can be seamlessly integrated into an existing finite element framework for solid mechanics. In particular, it does neither depend on any specific beam formulation nor the applied spatial discretization scheme, and in the context of the present work, we have exemplarily used it with geometrically exact Kirchhoff-Love as well as Simo-Reissner type beam finite elements. Likewise, it is independent of the temporal discretization and we have used it along with static and (Lie group) generalized-alpha time stepping schemes as well as inside a Brownian dynamics framework.
The accuracy of the proposed SSIP laws and the general SSIP approach has been studied in a thorough quantitative analysis using analytical as well as numerical reference solutions for the case of vdW and electrostatic interactions. We find that a very high level of accuracy is achieved for long-range interactions such as electrostatics both for the entire range of separations as well as all mutual angles of the fibers from parallel to perpendicular. In the case of short-range interactions, however, the derived SSIP law without cross-section orientation information slightly overestimates the asymptotic power-law exponent of the interaction potential over separation. As a pragmatic solution, a calibration of the simple SSIP law has been proposed to fit a given reference solution in the small yet decisive range of separations around the equilibrium distance of the LJ potential. In the authors' recent contribution, 46 this strategy led to very good agreement in the force response on the system level. While this accuracy might already be sufficient for certain real-world applications, our future research work will focus on the derivation of enhanced interaction laws including information about the cross-section orientation with the aim to achieve higher accuracy and the exact asymptotic scaling behavior.
The presented set of numerical examples finally demonstrates the effectiveness and robustness of the SSIP approach to model steric repulsion, electrostatic, or vdW adhesion. Several important aspects such as the influence of the Gauss integration error and the spatial discretization error as well as local and global equilibrium of forces and conservation of energy are studied in these simulations, including quasi-static and dynamic scenarios as well as arbitrary mutual orientations and separations of the interacting fibers. In order to remedy the characteristic singularity of inverse power interaction laws in the limit of zero separation, we have proposed a numerical regularization of the LJ SSIP law, which leads to a significant increase in robustness and efficiency, saving a factor of five in the number of nonlinear iterations while yielding identical results.

APPENDIX A. EXAMPLES FOR THE DERIVATION AND ANALYSIS OF THE TWO-BODY INTERACTION POTENTIAL AND FORCE LAWS FOR PARALLEL DISKS AND CYLINDERS
The aim of this appendix is to present the mathematical background of analytical solutions for two-body interaction potential as well as force laws. Generally, the strategy of pairwise summation, that is, integration of a point-pair potential, is applied. See Section 2.3.2 for a discussion of the applicability of this approach. Exemplarily, we consider the interaction between two parallel disks and two parallel cylinders since these scenarios proved to be most important throughout the derivation of SSIP laws as well as their verification in Sections 5 and 7.1, respectively. In addition, we are interested in the total LJ interaction potential and force law in the limit of small separations because the regularization proposed in Section 6.3 is based on these theoretical considerations. Finally, also, the equilibrium spacing g LJ,eq,cyl||cyl of two infinitely long cylinders interacting via the LJ potential will be derived and has proven helpful in order to choose an almost stress-free initial configuration of two deformable, straight fibers in the authors' recent contribution 46 studying the peeling and pull-off behavior.

A1 A generic interaction potential described by an inverse power law
Instead of Φ vdW (r) from Equation (4) or any other particular interaction type, here, we rather use the more general power law Φ m (r) = k m r −m for the point-pair potential. As noted already in Reference 33, this does not introduce any additional complexity in the derivations and the solutions can directly be used for other exponents m. We will make use of this fact when considering LJ interaction between two disks and two cylinders analytically in Appendix A3.1 and Appendix A3.2, respectively. These findings are to be used in the context of deriving a proper regularization of the potential laws in Section 6.3.

A1.1 Disk-disk interaction
The following refers to the analytical solutions for the disk-disk vdW interaction potential from the literature that is summarized in Section 2.3.2. Let us first state the underlying mathematical problem. We would like to find an analytical solution for the required 4D integral C m over the circular area of each disk in order to arrive at the disk-disk interaction potential

A1.1.1 Details on 2(a) in Section 2.3.2: The regime of large separations
For the limit of large separations g ≫ R 1 , R 2 , the solution is quite straightforward and can be explained in simple words as follows. The distance of any point in a disk to its center is of order (R) and thus much smaller than the disks' surface-to-surface separation g:r Figure A1 illustrates the introduced geometrical quantities. The distance r between any two points x 1 in disk 1 and x 2 in disk 2 may, therefore, be approximated by the interaxis distance d = g + R 1 + R 2 : Double integration over both disks is hence equivalent to a multiplication with the disks' areas and finally, we end up with the sought-after expression for the general disk-disk interaction potential in the limit of large separations̃m Note that this approximation is valid for arbitrary pair interaction functions Φ(r). Moreover, this solution does not even depend on the parallel orientation of the disks. It is valid for all mutual angles of the disks, which is important because we will apply it to arbitrary configurations of deflected beams. For the special case of parallel disks, this result can alternatively be obtained by the sound mathematical derivation of Reference 33, eq. (10). The leading term of his hypergeometric series is identical to the right hand side of Equation (A6).

F I G U R E A1
Two circular cross-sections, that is, disks in parallel alignment [Colour figure can be viewed at wileyonlinelibrary.com]

Remarks.
1 Note that this solution is valid for exponents m > 7∕2 only. This is in contrast to the approximation for large separations (A6) which is valid for arbitrary forms of the pair interaction potential Φ(r). 2 Note, however, the conceptual similarity of this expression to the one valid for the limit of large separations (A6). Here, we also find a power law, however, in the surface-to-surface distance g instead of the inter-axis distance d and with a different exponent.

A1.2 Cylinder-cylinder interaction
Considering the case of two parallel cylinders, we are interested in the length-specific interaction potential m,cyl||cyl = lim The integral over s 1 = −l 1 ∕2 … l 1 ∕2 yields a factor of l 1 since the integrand is constant along s 1 and thus immediately cancels with the normalization factor 1∕l 1 . Exemplarily, we want to discuss the more interesting and challenging regime of small separations here. Following Reference 33 (p. 63), one can interchange the order of integration, solve the integral over the infinitely long cylinder length analytically in a first step, and then make use of the generic solution for C m,ss from Equation (A11) Plugging in m = 6 for vdW interaction directly yields the two-body interaction potential per unit length for two parallel cylinders in the regime of small separations̃v dW,cyl||cyl,ss as stated in Equation (14). This generic expression (A15) will be exploited when deriving the total LJ interaction law in Appendix A3.2.

A2 Lennard-Jones force laws in the regime of small separations
As compared to the preceding sections, we now want to turn to the LJ interaction consisting of two power law contributions, one adhesive and one repulsive, respectively. Our motivation is to study the characteristics of the resulting, superposed force laws for disk-disk as well as cylinder-cylinder interactions by means of theoretical analysis of the analytical expressions. These findings shall prove valuable when deriving an effective yet accurate regularization of the LJ potential law for the limit of zero separation in Section 6.3. We, therefore, focus on the regime of small separations throughout this section. Coming from the expressions for the two-body interaction potential̃m ,disk||disk,ss and̃m ,cyl||cyl,ss derived for a generic point pair potential Φ m in Appendix A2, we will now sum the adhesive contribution m = 6 and the repulsive contribution m = 12 and differentiate once to arrive at the desired LJ force laws.
For later use in the analysis of the force law, let us restate the conversion from one set of parameters k 6 , k 12 specifying the point pair LJ potential to the other commonly used set Φ LJ,eq , r LJ,eq according to Equation (8): k 6 = 2Φ LJ,eq r 6 LJ,eq and k 12 = −Φ LJ,eq r 12 LJ,eq . (A18) Differentiation with respect to the separation yields the disk-disk LJ force law f LJ,disk||disk,ss = − d̃L J,disk||disk,ss dg = 5 2k 6 g − 7 2 + 17 2k 12 g − 19 2 .
See Section 6.3 for a plot of the function. This expression allows us to determine some very interesting, characteristic quantities like the equilibrium spacing g LJ,eq,disk||disk , that is, the distance where the force vanishes: Due to the fact that repulsive contributions from proximate point pairs decay faster than the adhesive contributions, we obtain a smaller equilibrium spacing as compared to the scenario of a point pair. Another differentiation allows us to determine the value of the force minimum, that is, the maximal adhesive force, and the corresponding separationf ≈ 0.7718448 r LJ,eq ≈ 1.18107 g LJ,eq,disk||disk .
These quantities turn out to be decisive for the choice of a regularized, that is, altered force law that is to be used instead of the original one in order to cure the numerical problems that come with the singularity at zero separation g = 0.
In summary, we have found an analytical, closed-form expression for the disk-disk LJ force law (A19), valid in the regime of small separations and for parallel disks. By means of elementary algebra, we were thus able to determine analytical expressions for the characteristic equilibrium spacing as well as value and spacing of the force minimum.

A2.2 Cylinder-cylinder interaction
As before in Appendix 1.2, we want to restrict ourselves to parallel, infinite cylinders and consider the length-specific interaction potential as well as force law. Again, starting from the expression for a generic interaction potential (A15), superposition yields̃L J,cyl||cyl,ss =k cyl,6 g − 3 2 +k cyl,12 g − 15 2 , where the following abbreviations have been introduced: that shall be further analyzed in the following. To begin with, the equilibrium spacing for two parallel cylinders interacting via a LJ potential can be derived as This is an extremely interesting and important result since it leads the way to the nontrivial stress-free configuration of two flexible, initially straight fibers. We make use of this knowledge for example in Reference 46. Again, since the repulsive contribution of proximate point pairs decays faster than the adhesive contribution, this equilibrium spacing is smaller than g LJ,eq,disk||disk for the disks, which, in turn, is smaller than r LJ,eq in the fundamental case of a point pair. The very same value of 57% of the point pair equilibrium spacing has already been mentioned as a side note by Langbein [p. 62], 33 however, without presenting the detailed, comprehensive derivation. In addition to the equilibrium spacing, we can again determine and look at the value and location of the force minimum ≈ 0.70104 r LJ,eq ≈ 1.22625 g LJ,eq,cyl||cyl .
Here, we find that the force minimum, that is, the maximal adhesive force is slightly shifted toward a smaller separation as compared to the disk-disk interaction. However, expressed in terms of the respective equilibrium spacing g LJ,eq,cyl||cyl , the value is slightly larger as compared to 1.18 g LJ,eq,disk||disk from Equation (A22). With these results, we conclude the derivation and analysis of LJ force laws in the regime of small separations and summarize the most important results in the following table to serve as a quick access reference. Table A1 gives an overview of some important quantities characterizing the LJ force laws for point-point, parallel disk-disk, and parallel cylinder-cylinder interactions.

APPENDIX B. LINEARIZATION OF THE VIRTUAL WORK CONTRIBUTIONS FROM MOLECULAR INTERACTIONS
Generally, the discrete element residual vectors r ia,j from molecular interactions between two beam elements j = 1, 2 depend on the primary variablesx k of both beam elements k = 1, 2. Consistent linearization thus yields the following four sub-matrices k jk to be considered and assembled into the global stiffness matrix, that is, system Jacobian K: Note that the linearization with respect to the primary variablesx k of both interacting beam elements simplifies due to the fact that the residuals r ia,j do not depend on the cross-section rotations as discussed along with the derivation of the specific SSIP laws in Section 5. Thus, only the linearization with respect to the centerline degrees of freedomd k yields nonzero entries and are, therefore, presented in the remainder of this section. B1 Short-range volume interactions such as vdW and steric repulsion The linearization of the residual contributions with respect to the primary variablesx of both interacting beam elements is directly obtained from differentiation of Equations (42) and (43): See again Equation (39) for the definition of the shape function matrices H j . As mentioned before, the discrete element residual vectors in the specific case of Coulombic interactions directly follow for m = 1 and c m,ls = C elstat 1 2 . See Section 2.2.1 for the definition of C elstat and Section 5.3 for the definition of the linear charge densities i . Again, as mentioned already in Section 5.3, the case of long-range volume interactions only requires to adapt the constant prefactor via c m,ls = kmA 1 A 2 1 2 .

APPENDIX C. SUPPLEMENTARY INFORMATION ON ALGORITHMS AND CODE FRAMEWORK USED FOR THE SIMULATIONS
Implementation All novel methods have been implemented in C++ within the framework of the multipurpose and multiphysics in-house research code BACI. 50 Integration into existing code framework The novel SSIP approach can be integrated very well into an existing nonlinear finite element solver for solid mechanics. In particular, it does not depend on a specific beam (finite element) formulation and has been used with geometrically exact Kirchhoff-Love as well as Simo-Reissner beam elements. Also, it is independent of the temporal discretization and has been used along with statics, Lie group Generalized-alpha, as well as Brownian dynamics. Load/time stepping We either applied a fixed step size or an automatic step size adaption that is outlined in the following. Starting from a given initial step size, a step is repeated with half of the previous step size if and only if the nonlinear solver did not converge within a prescribed maximum number of iterations. This procedure may be repeated until convergence is achieved (or until a given finest step size is reached which aborts the algorithm). After four subsequent converging steps with a new, small step size, the step size is doubled. Again, this is repeated until the initial step size is reached. Nonlinear solver The Newton-Raphson algorithm used throughout this work is based on the package NOX, which is part of the Trilinos project. 51 Unless otherwise stated, the Euclidean norms of the displacement increment vector and of the residual vector are used as convergence criteria. Typically, the corresponding tolerances were chosen as 10 −10 and 10 −7 , respectively. In some of the numerical examples, an additional Newton step size control is applied. It restricts the step size such that a specified upper bound of the displacement increment per nonlinear iteration is not exceeded. In simple terms, it is meant to prevent any two points on two beams from moving too far and eventually crossing each other without being detected from one iteration to the other. For this reason, the value for this upper bound is typically chosen as half of the beam radius. Linear solver We use the algorithm UMFPACK, 52 which is a direct solver for sparse linear systems of equations based on LU-factorization and included in the package Amesos which is part of the Trilinos project. 51 Parallel computing The implementation of the novel methods supports parallel computing and is based on the package Epetra, which is part of the Trilinos project. 51 See Section 6.6 for details on the partitioning of the problem in the context of the search algorithm applied to identify spatially proximate interaction partners. Post-processing and visualization The computer program MATLAB 53 was used to postprocess and plot simulation data. All visualizations of the simulation results were generated using Paraview. 54