Using Monte Carlo to Simulate Complex Polymer Systems: Recent Progress and Outlook

Metropolis Monte Carlo has been employed with remarkable success over the years to simulate the dense phases of polymer systems. Owing, in particular, to the freedom it provides to accelerate sampling in phase space through the clever design and proper implementation of even unphysical moves that take the system completely away from its natural trajectory, and despite that it cannot provide any direct information about dynamics, it has turned to a powerful simulation tool today, often viewed as an excellent alternative to the other, most popular method of Molecular Dynamics. In the last years, Monte Carlo has advanced considerably thanks to the design of new moves or to the efficient implementation of existing ones to considerably more complex systems than those for which these were originally proposed. In this short review, we highlight recent progress in the field (with a clear emphasis in the last 10 years or so) by presenting examples from applications of the method to several systems in Soft Matter, such as polymer nanocomposites, soft nanostructured materials, confined polymers, polymer rings and knots, hydrogels and networks, crystalline polymers, and many others. We highlight, in particular, extensions of the method to non-equilibrium systems (e.g., polymers under steady shear flow) guided by non-equilibrium thermodynamics and emphasize the importance of hybrid modeling schemes (e.g., coupled Monte Carlo simulations with field theoretic calculations). We also include a short section discussing some key remaining challenges plus interesting future opportunities.


INTRODUCTION
Metropolis Monte Carlo is a powerful simulation technique for equilibrating the dense phases of complex systems and predicting their key physicochemical properties because of the freedom it provides to sample new points in phase space thanks to the design of artificial (even fictitious) trial moves that can take the system completely away from its natural trajectory. It is particularly suited for simulating the bulk phases of chain-like (polymer or macromolecular) systems, because one can think of several such moves (both simple and complex) and combinations thereof that can design in order to generate trial states. Combined, in particular, with methods such as replica exchange, these moves can accelerate system equilibration by several orders of magnitude compared to dynamic methods such as Molecular Dynamics (MD). This is especially important at low temperatures or for systems characterized by highly dense structures, since dynamics becomes too slow to be followed reliably and ergodically by a detailed dynamic method.
In the last years, the method has expanded considerably through numerous applications in a variety of systems, often in the form of hybrid schemes with theoretic methods, thus providing invaluable insight into their thermodynamic, structural and conformational properties (and indirectly the dynamic ones). Our goal in this rather short report is to highlight important advances over the last 10 years or so, with some more emphasis on new applications to complex systems. We start by presenting the basic concepts underlying the method, then we review progress in the development of new Monte Carlo moves or interesting implementations of existing moves to new systems, then we discuss the development of new software for the more friendly execution of Monte Carlo simulations, and finally we discuss findings from several applications to polymer systems. These include the use of Monte Carlo in addressing polymer self-assembly, structure and conformation in polymer nanocomposites, polymers under flow, confined polymers, grafted polymers, polymer rings and knots, semiflexible polymers, polymer networks, polymer crystallization, and polymerization reactions. We discuss recent advances in all these areas and outline interesting future directions. Given the brief character of the review, we would like to deeply apologize in advance if we failed to include some significant contributions over the period (last decade or so) covered. For exactly the same reason, we have restricted our presentation to applications related with system equilibration and relaxation for soft materials made up of oligomers or polymers. Thus, extremely interesting extensions such as implementations to calculate rate constants or to track physicochemical processes in the form of kinetic Monte Carlo are not covered.

BASIC CONCEPTS
Let ρ eq be the probability density in the equilibrium ensemble wherein the Monte Carlo simulation is carried out, and let us assume that in a given Monte Carlo step a move is attempted from an old state (o) to a new or trial state (n). Let us also denote by a (o → n) the corresponding stochastic matrix a of attempt probabilities, which is usually symmetric, i.e., it satisfies Then, in the Metropolis Monte Carlo method, the new state is accepted with probability [1]: with ρ eq (o) and ρ eq (n) denoting the probabilities of the system to be in state o or in state n, respectively. The corresponding Markov chain of states thus generated samples asymptotically the probability distribution ρ eq . In the canonical (NVT) ensemble, Equation (2) reduces to [2][3][4][5]: where U (n) pot is the difference U pot in the potential energies between new and old states, k B denotes the Boltzmann constant, and T is the absolute temperature.
The power of the Monte Carlo method in simulating complex physical systems lies in the fact that one has considerable freedom in choosing the matrix a as long as the requirement a (o → n) = a (n → o) is satisfied. For example, one can even think of applying totally unphysical moves for the system at hand, as far as the internal geometry and the molecular architecture of the constituent molecules are not destroyed. This can substantially speed up the rate with which the system moves through configuration space, which is particularly advantageous in the case of the dense phases of chain-like molecules (e.g., synthetic polymers and biopolymers) that are known to be characterized by a very broad distribution of relaxation times that renders their direct MD simulation a formidable task.
An interesting point in the Metropolis Monte Carlo scheme is that one can replace Equation (1) with the following more general condition: also known as detailed balance or condition of microscopic reversibility [3][4][5]. Then, the Metropolis criterion, Equation (2), becomes If, in addition, the elementary move is designed in a set of coordinates that differs from the configuration-space coordinates wherein the equilibrium probability density ρ eq was defined, then Equation (5) must be modified to account for the Jacobian of transformation J from one coordinate system to the other: From the elements of the underlying matrix a, the transition probability matrix π is defined next according to and, by definition, is stochastic [2]. Then, as the Monte Carlo iterations progress, the row vector .., ρ t (n) , ...) containing the a priori probabilities of all states after step t of the simulation converges to the desired equilibrium distribution lim t→∞ ρ t = ρ eq satisfying [4,5]: Equation (8) implies that the equilibrium distribution ρ eq is an eigenvector of the transition probability matrix π , with the corresponding eigenvalue equal to one [4,5].

MONTE CARLO MOVES FOR POLYMERS
Due to difficulties associated with excluded volume interactions, chain connectivity, and conformational stiffness [6], early Monte Carlo simulations were performed on lattice models [7], but soon implementations in continuous space appeared. Today, we can categorize Monte Carlo moves developed for polymers roughly into three groups: simple, complex, and advanced.

Simple Monte Carlo Moves
Simple Monte Carlo moves include [4]: reptation [8], endmer rotation [9], libration or flip [9], dimer flip [10,11], configurational bias (CB) [12][13][14][15], extended configurational bias [16][17][18][19][20][21][22], concerted rotation (ConRot) [23], intra-chain and interchain concerted rotation [24,25], generalized reptation [9,26], parallel rotation [27], (g) pivot [28,29], atom identity exchange [30], and volume fluctuation [9]. These moves apply only to a single polymer chain; moreover, reptation, configurational bias, end-mer rotation and generalized reptation apply only at chain ends. With these simple moves, only chains up to about 70 units long can be simulated. Special mention should be made of the reptation move (Figure 1), perhaps the very first Monte Carlo move devised for polymer chains; it is realized by cutting out a monomer from one end of the chain and appending it to the other end. Thus, the move simulates the diffusive (slithering-snake) motion of the chain. In the configurational bias move (Figure 2), on the other hand, one cuts out an entire segment composed of many monomers at one end of the chain and regrows it monomer by monomer in a biased way so as to avoid overlaps with the monomers on the same or nearby chains. The bias introduced in the construction of the new configuration is removed at the final stage of the move when the acceptance criterion is applied, by appropriately modifying it. As far as the local conformation at the interior of a chain is concerned, this is typically equilibrated with the help of the ConRot move, which involves the concerted rotation of an internal segment of the chain, of size equal to 5 monomers (Figure 3).

Complex Monte Carlo Moves
Complex Monte Carlo moves induce drastic reconfigurations of large internal sections within one or two polymer chains simultaneously. Initially they appeared in the form of chain breaking moves [31] in Monte Carlo simulations of polymers on lattices; later they gave rise to the so-called family of chainconnectivity altering moves. Following some rigorous schemes, these moves effect large conformational changes at the level of the end-to-end distance of one or two chains at the same time, which can dramatically accelerate the rate with which the long-range conformational features of the polymer are sampled, often however at the expense of introducing polydispersity in the polymer. The family of complex Monte Carlo moves includes [9,26,32,33]: end-bridging, directed internal bridging, directed endbridging, self end-bridging, and fusion-scission [34]. Introduction of these moves almost revolutionized the field of the molecular simulation of polymers [4,5]. Figure 4 presents a design of the end-bridging move. From a technical point of view, and in order not to violate the very important condition of microscopic reversibility, in a variable connectivity Monte Carlo move proper care should be taken to: (a) evaluate all possible geometric solutions to the underlying bridging problem, (b) incorporate appropriate Jacobians in the acceptance criterion (because the solution of the geometric problem is typically carried out in the space of generalized coordinates), and (c) address both the forward and the reverse problem pertinent to the move.
As already mentioned, chain connectivity altering algorithms typically induce polydispersity in the sample because the molecular lengths of the chains involved in these moves are altered. Because of this, the corresponding Monte Carlo simulations are carried out in a semigrand canonical ensemble, wherein the following quantities are kept constant: the total number of chains N ch , the total number of mers n, the pressure P, the temperature T, and the spectrum of relative chemical potentials µ * of all chain species within the system except two that are taken as reference species [26,35]. Such an ensemble is typically denoted as [N ch nPTµ * ] and the average molecular length is estimated from the distribution of chain lengths eventually sampled, which in turn is dictated by the applied spectrum of relative chemical potentials µ * .

Advanced Monte Carlo Moves: Simulation of Non-linear Chain Systems
Advanced Monte Carlo moves include Double Bridging (DB) and Intramolecular Double Rebridging (IDR) [36,37]. These are generalizations of the end-bridging move in the sense that they involve the sequential or simultaneous construction of two trimer bridges (instead of one). The latter family of moves has found tremendous applications in simulations of polymers with a non-linear chain architecture such as long-chain branched, triarm star, and H-shaped polymers [10,11,37,38]; they have helped to understand how chain dimensions and other important conformational quantities depend on the precise molecular architecture of the chains [39][40][41] (frequency and length of branches). A typical design of the double bridging move for linear polymers is shown in Figure 5.

NEW MONTE CARLO SOFTWARE
In their majority, the Monte Carlo moves discussed in the preceding section have been implemented for specific polymer chemistries and architectures using home-made algorithms. However, there exist a few software packages offering these Monte Carlo moves for generic polymer structures. Available Monte Carlo software includes: (a) the Enhanced Monte Carlo code [42] which was used in simulations of semicrystalline polymers [43][44][45][46]; (b) Towhee [47]; (c) Cassandra [48]; and (d) RASPA [49].
Recently, Alexiadis et al. [50] presented Chameleon, a new simulation software in C++ wherein several chain connectivity altering Monte Carlo moves were implemented. The new software can handle several polymer structures such as polyethylene (PE), polystyrene (PS) and polyvinyl chloride (PVC), and many different polymer architectures (e.g., linear and branched), described either through an allatom, a united-atom, or a coarse-grained model. Direct comparison of several structural and volumetric properties of the above systems obtained with Chameleon against available experimental data and previous computational works demonstrated excellent agreement.
For solute-solvent systems in liquid or gas phase, Cezar et al. [51] have introduced a new version of DICE (https://portal.if.usp. br/dice), a Monte Carlo algorithm for the molecular simulation of long flexible molecules, based on an improved implementation of the configurational bias move.

MONTE CARLO SIMULATION OF SELF-ORGANIZED POLYMER PHASES
A family of polymer systems for which molecular modeling could provide significant insight is that of soft nanostructured materials (e.g., organic semiconducting polymers, polypeptides, polymerosomes, micellar surfactants and many others) because of the capabilities it offers to follow chain self-organization into a variety of structures and morphologies depending on operating (e.g., temperature) and physicochemical (e.g., concentration, type of solvent present, pH) conditions. Because self-organization can occur over long time scales, many of these systems cannot be simulated by a brute force application of the MD method. This explains why people resort to coarse-graining [52], an approach, however, that suffers from many drawbacks [53,54]. In this case, atomistic Monte Carlo offers an excellent alternative, as one can still maintain the detailed molecular representation while overcoming the problem of long relaxation times through the    application of the artificial moves presented in section Monte Carlo Moves for Polymers. Directly simulating ordered phases of polymers or oligomers with Monte Carlo is best exemplified by the work of Alexiadis et al. [30] for an alkanethiol selfassembled monolayer on Au. These authors proposed a Monte Carlo algorithm involving a mix of moves such as reptation, flip, concerted rotation, atom identity exchange, configurational bias, and intramolecular double rebridging which led to the formation of a self-assembled monolayer (SAM) on the gold substrate starting from a completely random configuration of a certain number of alkanethiols above the substrate. The conformational characteristics of the formed monolayer (especially the tilt angle) came out to be in very favorable agreement with measured data already known in the literature. This was an excellent demonstration of the limitless capabilities of Monte Carlo, as the formation of the corresponding monolayer is practically impossible to achieve with MD except if one starts from a preassembled structure very close to the final (equilibrium) one. In a later study, a similar Monte Carlo algorithm was developed for the bulk phase self-assembly of semi-fluorinated alkanes [55].
Recently, Tsourtou et al. [56,57] redesigned several of the Monte Carlo moves discussed in section Monte Carlo Moves for Polymers for bulk models of oligo and polythiophenes by properly accounting for the regular presence of thiophene rings along their backbones. Monte Carlo moves implemented included bias reptation of an end thiophene ring, flip of an internal thiophene ring, rotation of an end thiophene ring, concerted rotation of three consecutive thiophene rings, rigid translation of an entire molecule, rotation of an entire molecule, and volume fluctuation. A schematic representation of the ring ConRot move is shown in Figure 6. According to this, the C-S-C trimer of a randomly selected thiophene ring (e.g., segment C12-S13-C14 in Figure 6) is re-positioned by having the two atoms C9 and C17 neighboring the trimer being displaced following small rotations of the rings on the left and right of the central ring around their respective bonds; this alters also the positions of atoms C6, C10, C16, and C20. At the end of a successful ring ConRot move, eleven inner atoms (C6, C9, C10, C11, C12, S13, C14, C15, C16, C17, and C20) have been displaced to new positions.
In the work of Tsourtou et al. [56], thiophene ring atoms in all moves implemented were assumed to remain rigid and strictly coplanar while inter-ring torsion and bond bending angles were assumed to be fully flexible governed by suitable potential energy functions.
With the new algorithm, the authors studied ( Figure 7) the different phases formed when α-sexithiophene (α-6T), an important thiophene oligomer, is cooled down to lower temperatures isobarically, starting from the isotropic (Iso) phase at a relatively high temperature (above 700 K). To avoid system size effects, a rather large simulation cell was used in the simulations. The Monte Carlo simulations were performed with a new united-atom model specifically developed for the purpose of the study. For comparison, the authors carried out similar MD simulations with a detailed, well-validated all-atom model, which showed four phase transitions: an isotropic-to-nematic (Iso-to-Nem) at 640 K, a nematic-to-smectic A (Nem-to-SmA) at 630 K, a smectic A-to-smectic C (SmA-to-SmC) at 620 K (demonstrating smectic polymorphism), and a SmC-to-crystallike (SmC-to-Cry) at 600 K. Similar transitions were observed in the Monte Carlo simulations except that no Nem phase was seen, which was attributed to the stiffer nature of the corresponding forcefield. Indirectly, this provides evidence that predictions of morphology from coarse-grained models should not be fully trusted because of the highly approximate character of the effective potentials utilized in such models.
Monte Carlo moves similar to those discussed above have been adapted [58][59][60][61] over the last years for simulations of jamming and crystallization of polymer molecules modeled as freelyjointed chains of tangent hard spheres of uniform size, both in the bulk and under confinement. Simulations carried out over a wide range of concentrations, from the very dilute up to the maximally random jammed state, addressed how factors like chain length, chain flexibility and volume fraction affect the structure and packing of polymer chains, with emphasis on phase transitions from disordered to more ordered structures. Recently, the work was extended to polymer chains interacting with the square wellpotential [62]. The authors also tested highly confined systems of such chains by letting the inter-wall distance approach the size of the polymer bead [61]. Due to attraction, distinct morphologies were developed, spanning the entire spectrum of structures from purely amorphous to well-ordered ones depending on the specific values of the model parameters assigned [62].
For similar systems (dense hard sphere polymer melts), Kampmann et al. [63] proposed a Monte Carlo algorithm for their off-lattice simulation using both cluster and swap moves. The algorithm was validated by comparing against other Monte Carlo and MD simulations for the same system. At short time scales, the event chain Monte Carlo algorithm was shown to exhibit Rouse dynamics. At intermediate or long time scales, on the other hand, and in the absence of swap moves, it followed reptation dynamics.
Monte Carlo moves have also been implemented by Reith and Virnau [64] for single flexible globular homopolymer chains. For chains larger than a few hundreds of monomers, analysis of correlation times between unknotted globular states showed that bridging moves become more efficient than slithering-snake ones. However, the performance of the moves should depend on long-range correlations (which were not considered in the analysis) as well as on the specific molecular model considered. Also, if chain stiffness is included, efficiency is expected to drop.

MONTE CARLO AS A TOOL FOR UNDERSTANDING THE STRUCTURE OF POLYMER NANOCOMPOSITES
Connectivity-altering Monte Carlo algorithms have been coupled with preferential sampling techniques by Pandey and Doxastakis [65] and Pandey et al. [66,67], together with an extended reptate algorithm to facilitate polymer mass transfer from the nanoparticle surface to the bulk polymer, to explore structural and conformational features of polymers in a polymer nanocomposite melt containing highly curved nanoparticles. Using a rather detailed molecular model, the authors found that when the size of nanoparticle becomes comparable to the Kuhn segment length of the polymer, long train segments are disfavored.
Vogiatzis et al. [68] and Vogiatzis and Theodorou [69] combined Monte Carlo simulations with a field theoretic  approach to address structural effects of long polymer chains next to nanoparticles. For example, the authors investigated the dependence of several physical properties (such as local polymer density around the nanoparticle, thickness of the depletion layer, bond orientation and scattering patterns) on chain length and grafting density. This was one of the very first times that a real experimental nanocomposite system was studied, since the coarse-grained nature of the model employed could address length scales comparable to those accessed by small angle neutron scattering experiments, thereby allowing for a direct comparison between simulation and experiments. For example, the scattering of the whole corona was analyzed (Figure 8). An important conclusion of the work was that the brush thickness increases with increasing molecular weight of the grafted chains, thus causing a shift of scattering peaks to lower q-values.
Dodd and Jayaraman [72] employed Monte Carlo to study polydispersity effects in the structure of polymers grafted on spherical surfaces. They examined the conformation of grafted polymer chains, the thickness of the polymer layer formed, and the distribution of free-end monomers within the grafted layer. At brush-like grafting densities, with increasing polydispersity index, the scaling exponent describing the variation of the radius of gyration of the grafted chains with their molecular weight was observed to approach that of a single chain grafted on the same nanoparticle; this happens because polydispersity decongests the brush from monomer crowding. At high polydispersity indices, chains shorter than the number average chain length had more compressed conformations while chains longer than the number average chain length stretched less than in the monodisperse case. Monte Carlo simulations were also employed to obtain intramolecular structure factors needed in theoretical approaches to polymer nanocomposites through the self-consistent polymer reference interaction site model (PRISM) [73]. Such a combined approach allows mapping out equilibrium structure and phase behavior of polymer nanocomposites involving polymer-grafted nanoparticles in the matrix. Heterogeneous (e.g., copolymergrafted nanoparticles) as well as polydisperse (e.g., bidisperse graft chain lengths) systems were studied [73].
At low enough temperatures, nanoparticles can act as effective heterogeneous nucleation sites, thus promoting polymer crystallization [74][75][76]. Here, Monte Carlo can provide information concerning the role of factors such as nanoparticle shape and size on polymer crystallization. According to the work of Gu et al. [74], one-dimensional nanoparticles are the most effective in inducing crystallization, leading to uniformlyoriented crystals. Additional Monte Carlo simulations by Nie et al. [75] addressed the competition for crystallization of mixed polymers grafted on a substrate. It was found that stereo-complex crystallites have higher thermal stability compared to homocrystallites, which appears to be supported by experimental observations. Such simulations are of relevance (e.g.) in studies involving stereo-complex formation in poly(L-lactide)-poly(Dlactide) (PLLA/PDLA) blends. Ming et al. [76] used Monte Carlo simulations to study how polymer crystallization is affected by chain grafting onto the surface of the filler. It was found that the induction period for nucleation in a grafted system is shorter than in the system where chains are not grafted, but eventually the same degree of crystallinity is obtained. With increasing molecular weight, the induction period for nucleation decreases followed by an increase in the degree of crystallinity finally observed.

MONTE CARLO SIMULATION OF POLYMERS UNDER FLOW
Metropolis Monte Carlo was originally designed to sample equilibrium states. However, as Kikuchi et al. [77] showed, the method can also be used as a numerical technique for solving the Fokker-Planck equation, thus providing information for the diffusive motion of a Brownian particle. Later, Sanz and Marenduzzo [78] presented a comparative study between Monte Carlo and Brownian Dynamics simulations for colloidal suspensions and reported that for the results from the two methods to be in agreement, the Monte Carlo time must be rescaled by the acceptance probability. A subsequent study [79] demonstrated that one can achieve quantitative agreement between Brownian Monte Carlo and Brownian Dynamics simulations also in the case of systems with orientational degrees of freedom.
From a strict statistical mechanics point of view, using Metropolis Monte Carlo to sample states in systems beyond equilibrium requires the introduction of expanded (generalized) ensembles wherein, in addition to the typical macroscopic variables routinely employed, a coarse-grained variable should be included accounting for the average polymer conformation that develops in response to the external field. For melts of unentangled polymer chains, a good such structural variable is the dimensionless conformation tensorc [80,81]. Then, guided from non-equilibrium thermodynamics [80][81][82][83][84][85][86], one can postulate the following expression for the internal energy function U of the polymer melt under the external field a: dU (S, V, N ch ,c) = TdS − PdV + µdN ch + N ch k B Ta : dc (9) In Equation (9), T is the absolute temperature, k B the Boltzmann constant, S the entropy, P the pressure, V the volume, µ the chain chemical potential, N ch the number of chains in the system,c the dimensionless conformation tensor, and a a tensorial variable (the synthetic field) carrying information about the underlying (true) flow field that cannot be used directly in the Monte Carlo simulation. Equation (9) denotes the corona structure factor, n g the number of grafted chains on the nanoparticle, and N g the length of grafted chains in Kuhn segments) of a system consisting of an 8 nm-radius silica nanoparticle on the surface of which atactic polystyrene chains of molecular weight equal to 20 kg/mol have been grafted at surface density equal to 0.5 nm −2 . The grafted corona is inside a 100 kg/mol atactic polystyrene matrix. A piecewise cubic Hermite (PcH) spline interpolation scheme has been used. Results from two theoretical models [the form factor of a spherical shell of uniform density and thickness equal to the estimated brush thickness, and the model proposed by Pedersen and Gerstenberg [70] and Pedersen [71] for block copolymer micelles] are also included. Reprinted (adapted) with permission from Vogiatzis and Theodorou [69]. Copyright (2013) American Chemical Society.
typical variables N ch , T, and P) the synthetic field a. If we further allow for changes in the connectivity of chains implying a polydisperse system (this is important in cases where drastic Monte Carlo moves are incorporated in the algorithm such as the generalized reptation and the end-bridging ones), the simulation must be carried out in the expanded semi-grand canonical ensemble [N ch , n, P, T, µ * , a] where the spectrum of chemical potentials µ * also appears controlling the chain length distribution together with the total number of atomistic segments (or monomers) n which remains constant in the course of the simulation. The corresponding probability density function is [82]: implying that system configurations are sampled according to the following generalized Metropolis criterion: or, for the case of a system simulated under conditions of constant volume V (i.e., constant density ρ) according to: where β ≡ 1 k B T . In the above equations, {r } = {r 1 , r 2 , ..., r n } denotes the space of atomic position vectors, U pot is the potential energy of the system, µ * k N ch k=1 the set of chain relative chemical potentials, and {c k } N ch k=1 the set of chain conformation tensors. For the method to be applied, one needs to provide input data for the tensor a. Defining, however, a for a given flow field is not an easy task. Here, we can get help by looking at the corresponding evolution equation for the variablec as derived from the generalized bracket [80] and GENERIC [81] formalisms of non-equilibrium thermodynamics. Through this, we find that a conveys information related to the underlying dissipative or relaxation matrix of the viscoelastic model describing the polymeric fluid under study. Any proposition, therefore, for a from theory would be model-dependent, and thus also approximate. For example, in the case of steady shear flow described by the following velocity gradient tensor and for the simple upper-convected Maxwell model, we find that where λ 0 denotes the longest relaxation time of the polymer.
To compute model-independent values of a one should run GENERIC Monte Carlo simulations in parallel with direct nonequilibrium molecular dynamics (NEMD) simulations for the specified value of shear rateγ using (e.g.) the following more general form of a a =   a xx a xy 0 a xy a yy 0 0 0 0   (16) and then match the resulting chain conformation tensors from the two methods [84]. An interesting feature of such an approach is that by comparing the values of a computed from the GENERIC Monte Carlo simulations with those suggested by the macroscopic model, one can validate the model or deduce information how to modify it so that model predictions and simulation data for the same flow superimpose. It suggests therefore a solid framework for improving available macroscopic models for the viscoelasticity of polymer melts [84,85]. It can also be used to test fundamental thermodynamic equations far away from equilibrium. For example, the form of a provided by Equation (16) for a yy = 0 implies the following generalized Maxwell relation [86]: Test GENERIC Monte Carlo simulations provide direct evidence for the validity of this equation over a wide range of field strengths (Figure 9). The methodology was implemented initially for unentangled polymer melts although recently Roh and Baig [87] proposed a simple extension to an entangled polyethylene melt by treating each chain as a sequence of entanglement segments (each one consisting of approximately 68 carbon atoms), with the synthetic field a coupled simultaneously and independently with the conformation tensors of all these entanglement segments. Although all components of the tensor a had to be recomputed (compared to those estimated for an unentangled PE melt at the same temperature and flow conditions) in order for the GENERIC Monte Carlo predictions to match the direct NEMD ones for the overall conformational properties of the polymer (Figure 10), gratifying agreement was observed between the two methodologies in the regime of weak flows, but some systematic deviations were noted in the regime of weak-tostrong flows. These deviations were attributed to the fact that the synthetic field a was coupled independently to all entanglement strands along a chain, which fails to account for the interdependence of deformation between entanglement strands either intra-molecularly or inter-molecularly.

OTHER APPLICATIONS
In addition to addressing soft nanostructured materials, polymer nanocomposites, and steady-state polymer flows, Metropolis Monte Carlo has been used in the past few years to provide a wealth of information for other polymer systems and/or properties. We briefly review many of them in the next paragraphs.

Polymer Chain Stiffness
On the basis of single chain Metropolis Monte Carlo simulations, Tzounis et al. [88] proposed a methodology for computing the unperturbed chain dimensions of any polymer chain, irrespective of its chemical or architectural complexity. The authors tested many Monte Carlo moves and the most efficient ones in terms of their capability to relax both linear and branched polymer chains were found to be the rotate strand (pivot) and rotate branch ones, as they could induce drastic changes to the conformation of long strands and branches along the polymer chain. In addition, an iterative scheme was proposed to define local interactions. The method was used to predict the characteristic ratio of a series of polymers (Figure 11) and satisfactory agreement with experimental data was obtained.

Confined Polymers
Lattice Monte Carlo simulations have been combined with lattice self-consistent field (SCF) theory [89] to study athermal homopolymer solutions confined between parallel, nonabsorbing surfaces at equilibrium with a bulk solution. Calculations of the effective interaction between the two surfaces provided support for a fluctuation-induced repulsion between the confining surfaces at intermediate separation as had been suggested by Obukhov and Semenov [90] and Semenov and Obukhov [91].
Despite being inherently a non-dynamic method, Monte Carlo has been exploited to capture also the main features of polymer dynamics. It has been used [92] to simulate polymer diffusion in narrow periodic channels with alternating attractive and repulsive parts, and to study [93] the dynamics of a semiflexible polymer chain in the presence of an array of periodically distributed nanoparticles. In the former work, the diffusion coefficient was found to change periodically with the polymer length. In the latter work, the authors investigated polymer dynamics for repulsive, weak attractive, and strong attractive nanoparticles. An interesting finding in the case of strongly attractive nanoparticles was that a stiff polymer may move faster than a flexible one, since chain stiffness effectively weakens nanoparticle attraction. Monte Carlo has also been employed [94] to study the rupture of ultra-thin polymer film melts under strong confinement.

Grafted Polymers
The Monte Carlo moves described in section Monte Carlo Moves for Polymers were extended quite early by Daoulas et al. [95] to address melts of polymer chains grafted by one of their ends onto a solid substrate. The simulations were used to test theoretical scaling laws for the conformation and orientational order of grafted chains as a function of chain length and grafting density [96]. To account for the effect of grafting on the polydispersity of the melt (as induced by the chain connectivity moves included in the Monte Carlo algorithm), an iterative scheme was presented [35] capable of controlling the spectrum of chain relative chemical potentials discussed in section Complex Monte Carlo Moves so that the chain length distribution sampled in the course of the simulation is the desired one. This was achieved by accounting explicitly for interatomic interactions and is very important in reproducing and maintaining in the course of the simulation the desired polydispersity in applications of the method to inhomogeneous or anisotropic polymers.
Mendonça et al. [97] developed a configurational bias Monte Carlo method, as an alternative to dissipative particle dynamics (DPD) method, to study friction between grafted polymers in good solvent. The friction was examined as a function of intramolecular flexibility (controlled through the bondstretching and bond-angle bending potentials imposed in the simulation) based on predictions for the tangential component of the pressure induced by imposing a certain degree of mismatch in the registry of the two grafting surfaces. The main conclusion was that, the more flexible the polymer layer, the much lower the value of the shear force at which slip occurs.
Monte Carlo simulations with a single-site bond fluctuation model have also been used [98] to study protein adsorption on end-grafted polymers and investigate the role of polymer  hydrophilicity and grafting density on the conformation and height of the brush.
We also mention the Monte Carlo method introduced by Pakula [99] and Polanowski and Pakula [100] based on the dynamic lattice liquid model which was employed by Polanowski et al. [101] to simulate the growth of polymer brushes by the grafting-from method using atom transfer radical polymerization (ATRP). It allowed the authors to study the influence of the overall polymerization rate and the relationship between monomer attachment probability and activation/deactivation probabilities on polymer dispersity, polymer concentration profile and distribution of active chain ends in space. The work was extended to include also the dynamics of such realistic, polydisperse brushes starting from the "as obtained" (and usually non-equilibrium) state. This allowed investigating chain relaxation independently of chain growth, which in real ATRP synthesis is a much slower process, and results for various grafting densities and polymerization rates were presented.

Polymer Rings and Knots
Polymer rings is an intriguing class of polymers whose properties cannot be described by the reptation theory due to lack of chain ends. They have been the subject of intense research work in the past few years, especially after the pioneering work of Kapnistos et al. [102] was published according to which contamination by linear chains even at very small levels can dramatically affect the stress relaxation properties. Given the capabilities offered by simulations to precisely control sample composition, molecular simulation techniques such as Molecular Dynamics and Monte Carlo have emerged as excellent tools for studying the conformational and dynamic properties of polymer rings and ring-linear blends including their dependence on the relative concentration of the two species (ring and linear).
Lee and Jung [103] used lattice Monte Carlo simulations to study the connection between slow diffusional processes in melts of polymer rings and topological constraints associated with threading events between different ring molecules. Monte Carlo simulations on a face-centered-cubic lattice were employed by Suzuki and collaborators in a series of papers [104][105][106][107][108] to study: (a) the size and conformation of non-concatenated polymer rings in the melt state [104], (b) molecular conformations of trivial, 3 1knot, and 5 1 -knot ring polymers at their theta points [105,106], (c) the temperature dependence of the second virial coefficient of ring polymers [107], and (d) the size and conformation of catenated ring polymers in dilute solution, over a wide range of chain lengths [108].
Lattice Monte Carlo simulations have been used by Reigh and Yoon [109] to study concentration effects on the conformation of ring polymers. Shanbhag [110] and Henke and Shanbhag [111] have made use of Monte Carlo simulations with the bond fluctuation model to study dynamics (chain center-of-mass diffusion) and conformational properties (radius-of-gyration and end-to-end distance) in symmetric and asymmetric ringlinear polymer blends.
In a very recent study, Monte Carlo simulations based on the reptation (slithering-snake) and flip (crankshaft) moves were employed to sample configurational space and thus parameterize a generic soft repulsive potential of a mesoscopic, worm-like chain model for polymer knots (i.e., closed loops typically categorized according to the minimum number of crossings in a projection onto a plane) [112]. A similar method had been followed in the past to determine knotting probabilities and typical sizes of knots in double-stranded DNA containing up to half a million base pairs [113].

Semiflexible Polymers
Monte Carlo simulations of self-avoiding walks on a simple cubic lattice have been employed to study semiflexible polymers (chains with variable flexibility) and test the applicability of the Kratky-Porod model [114] and the structure of bottle-brush polymers [115], as well as to get information for the force-vs.extension behavior of flexible chains and semiflexible bottlebrush polymers adsorbed from a good solvent on a planar substrate using the bond fluctuation model [116,117]. The same bond fluctuation model was used to study transitions associated with orientational ordering in thin films [118] or near walls [119], at nematic ordering. For example, Monte Carlo simulations have confirmed [120] surface-induced nematic ordering in semi-dilute solutions of chains with a moderate stiffness in films that are thick enough for their central region to exhibit bulk behavior.
Greco et al. [121] used coarse-grained Monte Carlo simulations to study conjugated polymers forming nematic mesophases. In the method, polymer chains were described in the context of the discrete worm-like chain model subject to a soft non-bonded potential containing an isotropic repulsive and an anisotropic attractive term (of the type of Maier-Saupe) capable of inducing nematic order.
Manca et al. [122] used Monte Carlo to study how flexible and semiflexible single polymer chains anchored by one of their ends stretch due to an applied external field, and the results were used to validate a statistico-mechanical analytical model on the basis of the freely-jointed and worm-like chain models.
In the past, Monte Carlo simulations of long, single linear polymer chains subjected to uniform stretching by Pierleoni et al. [123] provided confirmation of the tensile blob concept introduced earlier by Pincus [124] and de Gennes [125]. In a later study, Titantah et al. [126] used Monte Carlo to study the thermo-elastic properties of individual polyethylene chains at theta and good solvent conditions by employing a realistic model for polyethylene, in two different ensembles (fixedexternal stretching force on chain and fixed relative extension of chain).

Physical Gels, Hydrogels, and Networks
Monte Carlo simulations with the bond fluctuation model have been employed [127][128][129] for star-polymer networks. Lange et al. [128] reported Monte Carlo simulations of the network structure of a tetra-PEG hydrogel (formed by an A-B reaction of two symmetric four-arm polymers), which was experimentally characterized using proton multiplequantum NMR measurements at low field. Different connectivity modes (regular single links and double links) between the macromonomers were suggested for individual stars, together with some other network structures characterized by lower order parameters. In this case, the simulations confirmed the concentration-dependence of the network structure and the fraction of double links probed experimentally. Monte Carlo simulations with the bond-fluctuation model were also performed by Lang et al. [129] to test the predictions of a rate theory for the formation of short cyclic structures in networks formed by the homo-and co-polymerization of ffunctional molecules.
Bergsma et al. [130] combined Monte Carlo with the Scheutjens-Fleer self-consistent field theory to address gels of ABA triblock copolymers as a function of polymer volume fraction. Associative A blocks were confined to small volumes called nodes, the number of polymers per node being a parameter. B blocks, on the other hand, were allowed to move freely as long as they were connected to A blocks. With Monte Carlo, the authors could sample node configurations on a lattice (Figure 12) while with the Scheutjens-Fleer theory they could determine changes in free energy. The proposed method (a hybrid Monte Carlo-SCF scheme) has the advantage that the interaction potential between nodes needs not to have been defined in advance. FIGURE 12 | An example of the system simulated by Bergsma et al. [130]. Dark red cubes indicate nodes (whose cores are shown with a slightly lighter red color). Due to steric repulsion, polymer molecules push each other away from the node, this drags anchoring groups outside the core, and therefore the density within the core decreases. Figure reprinted with permission from Bergsma et al. [130]. Notice to readers: further permissions related to the material excerpted should be directed to the ACS.

Polymer Crystallization
Monte Carlo has also been used to study the kinetics and morphology of polymers undergoing crystallization under an applied 3-d flow [131], or strain-enhanced stereo-complex polymer crystallization [132]. In another study [133], the method was extended to binary blends of symmetric crystallizable polymers in which the authors enhanced (separately) the driving forces leading to polymer-uniform and polymer-staggered crystals. Under parallel enhancements, polymer-uniform crystals were observed to exhibit faster nucleation and growth, more chain folds and less lamellar thickening than polymer-staggered ones. Tahara et al. [134] used Monte Carlo to simulate 2d small-angle X-ray scattering (SAXS) patterns exhibited by oriented crystalline polymer samples, and thus indirectly deduce information concerning the structure of lamellar aggregates (Figure 13). The 2-d model corresponded to the projection of the 3-d lamellar structure or to the averaged structure viewed along the direction of the beam.

Polymer Reaction Engineering
The Monte Carlo method has also been very popular in polymer reaction engineering because of its capability to describe polymerization reactions. This is an interesting subject by itself and we refer readers to the detailed review by Brandão et al. [135].

SUMMARY AND OUTLOOK
Out of our short review we hope it has been shown that, thanks to the development and efficient implementation of advanced moves that can induce large conformational changes in the interior of a polymer chain and, consequently, considerably accelerate the rate with which configurational space is sampled, Metropolis Monte Carlo has advanced today to a powerful computational tool that can be used to address several physicochemical properties of real polymers (structural, conformational, thermodynamic, and morphological) with remarkable accuracy. We should emphasize, in particular, the capability to treat polymer molecular organization and morphology at several scales. Direct observation of self-organization starting from a purely random configuration of an ensemble of polymer chains with the detailed Molecular Dynamics method is often impossible for polymers with high molecular weight due to the long time scales involved in the course of the phase transition, despite that molecular organization in condensed phases is determined predominantly by short-range repulsive or attractive forces between polymer atoms or groups of atoms. Resorting to a non-dynamic method such as atomistic Monte Carlo in which the path followed in going from the original completely disordered state to the final self-organized one departs considerably from the natural trajectory seems to be much more promising.
The benefits of using a non-dynamic method to simulate ordered phases without making big compromises in the degree of atomistic detail accounted for are far more reaching and important than what has been discussed here, since tiny details in the chemical structure (e.g., exact length of a side alkyl group, or exact degree of flexibility of the main chain) can dramatically affect the phase behavior of the system. For example, it can influence the precise morphological structure or even the number of intermediate ordered phases predicted as the value of a physical parameter (e.g., temperature, concentration or volume fraction) is varied. Failure to take proper account of the full complexity of interactions may lead to erroneous predictions concerning the full spectrum of morphologies manifested by the system under study. This renders atomistic Monte Carlo a much more reliable and attractive method to employ than techniques based on structural coarse-graining and the use of effective potentials that offer only an approximate treatment of the problem, the possible dependence of the parameters of the effective potentials on the conditions of the simulation (e.g., temperature and pressure or density) notwithstanding.
Despite its unlimited potential (the method is only limited by our imagination to devise new moves), there is still many issues and problems to be addressed for the method to reach the level of maturity of Molecular Dynamics. First of all, most Monte Carlo simulation algorithms today are highly system-specific, in the sense that the corresponding Monte Carlo moves are tailored for a particular polymer chemistry (e.g., polyethylene or polypropylene or polybutadiene) and/or molecular architecture (e.g., linear polymer, branched polymer, cyclic polymer). This strongly prohibits code transferability from the original polymer system for which the moves were developed to another. We are clearly in the need of further algorithmic developments before Monte Carlo can be regarded as a friendly, readyto-use simulation tool for different polymers. The problem becomes more acute considering the great complexity of certain polymers such as the simultaneous presence on the same chain of side groups, branches, rings, functional groups etc., or the emphasis today on multicomponent systems (e.g., polymer nanocomposites, polymer mixtures, block copolymers, etc.). Even if the corresponding moves are rigorously adapted for a new polymer ensuring an ergodic simulation, it is not at all clear that their acceptance rates will be sufficient to ensure complete relaxation at all length scales. This is often the case when the FIGURE 13 | Stacked lamellar models proposed [134] on the basis of information extracted from the Monte Carlo simulations of observed SAXS patterns of a biaxially oriented LLDPE (linear low density polyethylene) sample. Different sets of data are shown corresponding to SAXS patterns measured at different positions in the sample. Also shown is the comparison between experimentally observed and Monte Carlo simulated SAXS patterns. Reprinted (adapted) with permission from Tahara et al. [134]. Copyright (2020) American Chemical Society.
polymer contains relatively long side groups or stiff units or when specific structural features exclude the implementation of some key Monte Carlo moves like ConRot or end-bridging.
The obstacles discussed in the previous paragraph doubtless pose significant limitations to the applicability of Metropolis Monte Carlo to polymer systems of considerable industrial and technological interest presently. At the same time, however, they call for important developments in the next years along several directions. One such direction would be to develop new moves for the more efficient treatment of polymer segments around branch points or junctions; the latter would be of relevance to polymer networks. Although procedures for treating such points have already appeared, devising moves to further enhance equilibration around branch points or junctions would be an important step forward, since these points constitute the Achilles heel of the majority of Monte Carlo algorithms presently available. Even under the hypothesis that no new moves are introduced, significant progress could be made by extending existing ones (in particular, the most advanced chain-connectivity altering moves) to systems that depart structurally from the original polymer system for which the moves were developed (e.g., linear or branched polyethylene, or polybutadiene) but whose complexity remains simple enough to offer a relatively straightforward implementation. This is the case of block copolymers. To the best of our knowledge, casting Monte Carlo moves already developed for specific pairs of homopolymers (e.g., polystyrene and polyisoprene or polybutadiene) to their copolymers with the object of exploring their rich phase diagram has not been undertaken thus far. A direct treatment of the full phase diagram of block copolymers with Monte Carlo would not be prohibitively difficult and would constitute a significant step forward in the field indeed.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.