**1. Introduction**

The development of new materials with enhanced properties is one of the most interesting and important topics in research in materials science and engineering. To achieve this ambitious goal, one has to relate the behavior of atoms and molecules to the macroscopic properties of the end material. In this perspective, molecular simulation is of paramount importance, since it allows the study of materials at the atomistic/molecular level without needing an experimental process, which, in specific cases, can become expensive, time consuming, and environmentally hazardous. Over the years, different molecular simulation techniques and methodologies have risen to answer relevant questions of general atomic and particulate systems [1–5].

A system composed of macromolecules is a very challenging case from the perspective of molecular simulation. This stems from the fact that polymers are characterized by a wide spectrum of characteristic time and length scales. Their simulation can become prohibitively difficult when very long and well-entangled chains are involved due to the very slow dynamics. Added to this is the fact that atomistic simulations have to take into full account the chemical constitution of the repeat units and the corresponding bonded and non-bonded interactions. To address this problem, a large amount of different molecular simulation methods has been developed and constantly improved over the last

**Citation:** Herranz, M.; Martínez-Fernández, D.; Ramos, P.M.; Foteinopoulou, K.; Karayiannis, N.C.; Laso, M. Simu-D: A Simulator-Descriptor Suite for Polymer-Based Systems under Extreme Conditions. *Int. J. Mol. Sci.* **2021**, *22*, 12464. https://doi.org/ 10.3390/ijms222212464

Academic Editor: Małgorzata Borówko

Received: 5 October 2021 Accepted: 12 November 2021 Published: 18 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

decades. The choice of simulation approach/scheme depends on the system/phenomenon, its physical-chemical details, size, and properties of interest. For example, Molecular Dynamics (MD) provides dynamical information at the local level of segments and global one of chains. However, as it follows the evolution of the equations of motion in time, it can be too slow to be effective when very long chains are involved. Monte Carlo (MC), by resorting to different stochastic algorithms ("moves"), can offer rapid equilibration at all length scales. However, MC cannot provide any information about the real dynamics. Accordingly, it is not uncommon for different techniques to be combined together into powerful hierarchical modeling approaches. The MD approach is widely used when dynamical or temporal evolutions are of interest. One of the most widely used software packages for the simulation of synthetic polymers is LAMMPS [6], which has been further used in other tools such as Polymatic for the polymerization of amorphous polymers [7]. Regarding the modeling of biomolecular systems, NAMD [8] and GROMACS [9] are two of the most popular simulation software, both placing special emphasis on parallelization in order to enhance performance. Other relevant open MD software suites are ms2 [10], to extract thermodynamical properties of homogeneous fluids using hybrid parallelization on MPI and OpenMP [11]; MOLDY [12] for solids and liquids under periodic boundary conditions; or GULP [13] for solids. Commercial suites include, among others, CHARMM [14], AMBER [15], and HyperChem [16].

With respect to Monte Carlo simulations, homemade software programs usually targe<sup>t</sup> a specific type of polymer structure, either of its chemistry or the architecture of the chain, but most of them follow rather similar approaches. Monte Carlo simulations are applied when equilibrium structural properties, including phase transitions, constitute the main research focus. The Enhanced Monte Carlo code [17,18] is a multi-purpose modular environment for particle simulations using force fields such as COMPASS, CHARMM, or Born. This open tool has been used to study the effect of semicrystalline interphase polyethylene under different conditions of tensile deformation [19,20] or chain branching [21]. MCCCS Towhee [22] was initially developed as a simulator suitable for computing phase equilibria in the Gibbs ensemble, but later extended to different force fields and ensembles. As an example, this open tool has been used to study gas-mixture separations on clathrate hydrates [23] among many other studies. DL\_MONTE [24] is a very recent MC-based open tool that can be applied to general atomistic systems under different force fields and ensembles, as well as introducing transition pathways of umbrella sampling and Wang–Landau [25]. Furthermore, it is compatible with the molecular dynamic tool DL\_POLY or chain branching [21].

We should also mention other relevant open-source MC-packages, such as Cassandra [26], that can be applied to obtain thermodynamic properties of fluids and solids; RASPA [27] for simulating adsorption and diffusion phenomena; GOMP [28,29] for GPU optimized phase equilibria simulations; or DICE [30] that uses a configurational bias scheme to study flexible molecules in solute-solvent systems. Most relevant MC software packages are benchmarked in terms of computational efficiency using adsorption simulations [31]. Regarding realistic polymeric systems, Chameleon [32] is one of the latest available pieces of software. This tool employs different chain connectivity altering moves to simulate atomistically detailed polyethylene (PE), polystyrene (PS), and polyvinyl chloride (PVC) for different polymer architectures.

Usually, the development of a commercial or open code, especially when built around Monte Carlo algorithms (moves), requires a major effort and programming in order to make it user-friendly, efficient, and of general applicability. Besides, it is very common that clever MC-based or general structure-optimization algorithms have and are being developed for specific applications or general classes of physical problems in continuous or lattice cells and in systems of varied chemical detail, in the bulk and under confinement [33–52].

Equally important to the simulation itself is the post-simulation analysis. This step can include visualization, including 3-D representation and animation, of the computergenerated system configurations and calculation of relevant quantities through proper

interpretation of the raw simulation data. Corresponding suites also exist for interactive visualization, description, and analysis including, among others, ParaView [53], VMD [54], disLocate [55], UCSF chimera [56], OVITO [57], and i-Rheo GT [58].

In the present manuscript, we analyze the main features of Simu-D, an MC-based simulator and structural descriptor suite for the molecular modeling of polymer-based systems under extreme conditions. The simulator, which is the central component of the present software, is effectively the accumulation of successive expansions, modifications, and improvements implemented on the MC code [59], originally built for the simulation of dense and jammed athermal polymer-based systems in the bulk. The structural descriptor is the latest version of the Characteristic Crystallographic Element (CCE) norm [60,61], a metric used to gauge the similarity of local structure with respect to reference crystals in general atomic and particulate systems. Over the last years, the MC suite has been extended to simulate athermal polymers under confinement [62] and more recently macromolecules whose monomers interact with the square well (SW) or square shoulder (SS) potential [63]. In the corresponding research studies, emphasis was placed on how the employed conditions affect the ability of chains to pack at the local and global level [64,65], the topological network of entanglements [66–68], and the entropy- or energy-driven phase behavior (crystallization) in the bulk and under extreme confinement [63,69–73]. Here, the suite is further extended to include additional factors: chain stiffness, blends of chains and monomers, spherical or cylindrical confinement, the varied potential for bonded and nonbonded interactions, nanofillers in the form of cylinders and spheres, and combinations of the above. The ongoing effort is to create a general-purpose simulator-descriptor suite that will be as efficient, general, and user-friendly as possible given the variety of simulation conditions to be considered and the stochastic nature of the underlying MC method.

The manuscript is organized as follows: Section 2 presents the molecular model, the interspecies interactions, and the systems under study. Section 3 presents the moves behind the MC simulator and briefly discusses the features of the CCE-based structural descriptor. Section 4 discusses results from representative applications of Simu-D. Finally, Section 5 summarizes the main conclusions and lists current efforts and plans.

#### **2. Molecular Model/Systems Studied**

The current version of Simu-D allows the simulation of atomistic systems composed of *N*at spherical monomers. These monomers can be part of macromolecules and/or exist as individual particles. In the general case, the system contains *N*ch chains with the average length of *N* and *N*s individual particles with *N*ch × *N* + *N*s = *N*at. Obviously, the two limiting cases correspond to the pure polymer matrix ( *N*s = 0, *N*at = *N*ch × *N*) and a system composed entirely of monomers ( *N*ch = 0, *N*at = *N*s).

Non-bonded atoms interact with a pair-wise potential, which can be discontinuous such as the hard sphere (HS) or the square well/shoulder (SW/ SS) ones or continuous such as Lennard–Jones (LJ) with the corresponding formulas being displayed in Equation (1).

$$\mathrm{LI}\_{\mathrm{HS}}(r\_{\vec{\eta}}) = \left\{ \begin{array}{c} 0, \ r\_{\vec{\eta}} \ge \sigma \\ \infty, \ r\_{\vec{\eta}} < \sigma \end{array}, \mathrm{LI}\_{\mathrm{SW/SS}} = \left\{ \begin{array}{c} 0, \ r\_{\vec{\eta}} \ge \sigma\_{2} \\ -\varepsilon\_{\mathrm{SW}}, \ \sigma \le \tau\_{\vec{\eta}} < \sigma\_{2} \\ \infty, \ r\_{\vec{\eta}} < \sigma \end{array}, \mathrm{LI}\_{\mathrm{L}} = \varepsilon\_{\mathrm{L}} \right\} \left[ \begin{array}{c} \sigma\_{\mathrm{L}} \\ \tau\_{\vec{\eta}} \end{array} \right)^{12} - \left( \frac{\sigma\_{\mathrm{L}}}{r\_{\vec{\eta}}} \right)^{6} \right\} \tag{1}$$

where *rij* is the distance of the centers of atoms *i* and *j* and *σ* is the collision diameter, which is further considered as the characteristic length of the system. *σ*2 and *εSW* correspond, respectively, to the range and intensity of the repulsive (SS) and attractive (SW) potentials. *<sup>ε</sup>LJ* and *<sup>σ</sup>LJ* are the depth and zero-energy point of the LJ potential. As in any traditional molecular simulation, depending on the type of the applied non-bonded potential, the original simulation cell is split automatically into overlap cells (HS), or into overlap and cut-off cells (SW/SS, LJ) to expedite the calculation of interactions.

Polymers are modeled as linear sequences of monomers of varying chain stiffness. Bond lengths can be longer (bond gaps), equal (bond tangency), or shorter (fused spheres) than the collision distance, *σ*. Chain stiffness is introduced through a potential governing

bending angle (supplement of bending angle, *θ*) formed by triplets of successive atoms along the chain backbone. The formula for the energetic calculations can be general. Configurations of semi-flexible chains have been generated in the present work with the following bending angle potential:

$$\mathcal{U}\_{bmd}(\theta) = k\_{\theta}(\theta - \theta\_0)^2 \tag{2}$$

where *kθ* is the bending constant and *θ*0 is the equilibrium bending angle supplement (i.e., a fully extended bending angle corresponds to *θ*0 = 0◦). For fixed bond lengths, setting *kθ* = 0 allows the simulation of freely jointed chains while *kθ* → ∞ corresponds to the freely rotating model. In the current version of the suite, torsion angles, *φ*, can also be controlled through the implementation of a torsional potential, *Utor*(*φ*). However, in all results presented below, torsion angles are allowed to fluctuate freely and thus chain stiffness is governed solely by the bending potential.

The presence and activation of specific MC moves, as will be described in the continuation, enforces dispersity in chain lengths. Such polydispersity is controlled by casting the simulation in the *Nat NchVTμ*\* ensemble where *V* is the total volume of the simulation cell, *T* is the temperature, and *μ*\* is the spectrum of relative chemical potentials of all chain species, as explained in detail in References [59,74]. The uniform and Flory (most probable) distributions of molecular lengths can be selected in the simulation of polydisperse systems. In the case that strictly monodisperse samples are required, then all moves that vary the chain length (sEB, x-reptation, and IdEx3, see below) are deactivated from the mix.

Depending on the system under study, initial configurations are generated under very dilute conditions and the system is brought to the desired density through compressions or simulations in the isothermal-isobaric (*NPT*) ensemble. For the latter, conventional volume fluctuation moves are attempted at regular intervals. For the former, cell compaction is achieved by a combination of volume fluctuation moves, and in the case of confined systems, the wall wrapping "MRoB" algorithm as explained in [75].

Simulations can be conducted in two or three dimensions under periodic boundary conditions or on flat surfaces. Confinement is realized through the presence of such impenetrable surfaces. The current implementation allows confinement in the form of (i) flat, parallel walls in at least one dimension, (ii) a cylinder with closed or open ends (subjected to periodic boundary conditions), and (iii) a sphere (full confinement). The intensity of confinement is controlled by the distance between the confining surfaces, i.e., the cylinder or sphere diameter or the inter-wall distance. The latter can, in general, be different in each confined dimension *i*, *dwall*(*i*). Simulation cells are always orthogonal but can be anisotropic, and the number of confined dimensions, *dconf*, ranges from 0 (bulk cell with periodic boundary conditions) to 3 (full confinement). The cell aspect ratio, *ζ*, is defined as the ratio of the maximum inter-wall distance divided by the minimum one [75].

Nanocomposites can be simulated with the fillers taking the form of spherical or cylindrical particles of varied sizes and populations. In the current implementation of the suite, each nanocylinder spans the whole simulation cell and its direction is held fixed throughout the simulation. Nanospheres can, in principle, move in space, but in all computer-generated polymer nanocomposite configurations to be presented in the continuation, they are treated as immobile inclusions.

For a bulk system of pure polymer, the matrix number density, *ρ*, is trivially defined as *ρ* = *Nat*/*V*, while for non-overlapping entities (such as hard spheres), packing density, *ϕ*, is given by:

$$
\varphi = \frac{V\_{\text{mou}}}{V} = \frac{\pi}{6} \frac{N\_{\text{at}}}{V} \sigma^3 = \frac{\pi}{6} \rho \sigma^3 \tag{3}
$$

where *V* is the volume of the simulation cell and *Vmon* is the volume occupied by the monomers, either as individual entities ("single monomers") or by being part of polymer chains ("chain monomers").

For interfacial/confined/composite systems, the above definition provides little information on the free or accessible volume given that for very large nanofillers, the volume

occupied by the nanofiller can be up to four orders of magnitude higher than the one of the monomers. Thus, we can further define an effective packing density, *ϕeff*, considering the reduction of the accessible volume due to the presence of the nanofillers as:

$$
\varphi\_{eff} = \frac{V\_{\text{mou}}}{V\_{\text{acc}}} = \frac{V\_{\text{mou}}}{\left(V - V\_{fill}\right)}\tag{4}
$$

where *V*acc is the volume accessible to the spherical monomers, *Vfill* (= *Vcyl* + *Vsph*) is the volume occupied by the nanofillers, being the summation of the volume occupied by *Ncyl* cylinders (*Vcyl*) and of *Nsph* spheres (*Vsph*). Additionally, in the calculations above, one could further incorporate a depletion layer as monomer centers cannot lie closer than σ/2 from the surface of nanofillers or walls. In the general case of a system under confinement and being composed of nanofillers, if *dconf* is the number of confined dimensions, the depleted effective packing density, *ϕdep*, including the effect of all nano-entities, can be calculated as:

$$\mathcal{P}\_{\rm dep} = \frac{V\_{\rm mon}}{V\_{\rm dep}} = \frac{V\_{\rm mon}}{\left(\prod\_{i=1}^{d\_{\rm conf}} (d\_{\rm wall} - \sigma) \prod\_{j=d\_{\rm config}+1}^{3} l\_j - \frac{\pi}{6} \left(d\_{\rm sph} + \sigma\right)^3 \mathcal{N}\_{\rm spl} - \frac{\pi}{4} \left(d\_{\rm cyl} + \sigma\right)^2 L\_{\rm cyl} N\_{\rm cyl}\right)} \tag{5}$$

where *dsph* and *dcyl* are the diameter of the nanospheres and nanocylinders, respectively, *Lcyl* is the nanocylinder length, index *i* runs over all confined dimensions, index *j* over all unrestricted ones, and *lj* is the length of the simulation cell in dimension *j*.

#### **3. Simulator-Descriptor Suite**
