The Martini Model in Materials Science

The Martini model, a coarse-grained force field initially developed with biomolecular simulations in mind, has found an increasing number of applications in the field of soft materials science. The model's underlying building block principle does not pose restrictions on its application beyond biomolecular systems. Here we highlight the main applications to date of the Martini model in materials science, and we give a perspective for the future developments in this field, in particular in light of recent developments such as the new version of the model, Martini 3.


I. INTRODUCTION
Coarse-grained (CG) force fields have gained a lot of popularity in the field of molecular dynamics (MD) simulations [1][2][3]. By averaging out some of the atomistic degrees of freedom, they allow for a substantial alleviation of both the spatial and temporal limitations of allatom models. The Martini model [4][5][6] is an example of a popular CG force field that has been incorporated by the worldwide user community to study a large variety of (bio)molecular processes [6,7].
In the Martini model, typically four heavy atoms with their associated hydrogen atoms are grouped into one functional group, denoted as a CG bead. This effectively reduces the number of particles to be simulated in a system, increasing the simulation speed. In addition, the smoother CG energy landscape leads to faster overall dynamics and allows the use of larger time steps compared to all-atom simulations. Together, this results in a significant increase in accessible length and time scales of a few orders of magnitude, albeit at somewhat reduced level of accuracy.
The CG beads represent small chemical fragments, and are used as building blocks for larger molecules. Parametrization of the non-bonded interactions between the CG beads is based on reproducing thermodynamic data such asfree energies of transfer of organic compounds. In addition, reference all-atom simulation data are used to derive effective bonded terms. Such a combination of top-down and bottom-up approaches enables the Martini model to distinguish different chemical species and form a useful bridge between atomistic and macroscopic scales.
From the first applications, purely concerned with lipids [8,9], the Martini model has been applied to a vast amount of biomolecular systems. The compatibility with a wide library of existing molecules, which includes all major biomolecules such as proteins [10], sugars [11], DNA [12], or RNA [13], as well as an increasing amount of synthetic molecules [6], is one of the key strengths of the Martini model. It enables researchers to easily simulate * s.j.marrink@rug.nl complex many-component systems and focus on more advanced simulation methodologies.
In recent years, the Martini model has found more and more applications in the field of materials science. In light of the underlying building block principle of the Martini model, there is no reason to restrict its applications to biomolecular systems. Indeed already more than a decade ago the Martini model has been applied to simulate polymeric systems [14]. In principle, any molecule can be represented by Martini CG beads, as illustrated in Figure 1. Based on this versatility, the Martini model is nowadays used to simulate a wide range of materials, including (block co)polymers, nanoparticle-polymer composites, organic electronic materials, ion-conducting materials, self-assembled supramolecular materials, ionic liquids, and others.
The use of CG models in material science is of course not new. In fact, some of the very first applications of CG models, such as the freely jointed chain models of Binder and coworkers [15], already targeted polymer dynamics. Since those early days, a large plethora of CG models have been developed to model an equally large variety of materials [16][17][18][19][20][21][22][23]. Two major assets of the Martini model, that set it apart from most other CG models, are (i) the retaining of near atomic resolution, as opposed to more generic models that are frequently used to capture global system properties but are unable to offer a direct connection to chemical specificity, and (ii) the broad range of compatible parameters available for different classes of molecules, enabling simulations of the ever expanding group of complex and hybrid materials as well as the interaction of materials with biological systems.
In the following sections, after summarizing Martini parametrization strategies tailored to material systems, we discuss the main application areas to date of the Martini model in materials science, illustrated with selected examples. While there are numerous Martini applications which involve the interaction of synthetic materials with biomolecules, we will not cover those here. The interested reader is referred to recent reviews and references therein [6,24,25]. We conclude with a perspective for the future developments in this field, in particular in light of the new version of the model, Martini 3, as well as recently developed tools to generate starting structures for polymeric systems and to allow constant pH simulations.

II. PARAMETERIZATION STRATEGIES
A. General guidelines Martini typically gathers groups of four non-hydrogen atoms in CG beads-see Figure 1 for some representative mappings. The interactions between beads-described by a 12-6 Lennard-Jones (LJ) potential-represent the nature of the underlying chemical groups and have been systematically parametrized to reproduce free energies of transfer of solutes between polar and non-polar solvents. There are four main particle types: polar (P), intermediately polar (N), apolar (C) and charged (Q). These types are in turn divided in subtypes based on their hydrogenbonding capabilities (with a letter denoting: d = donor, a = acceptor, da = both, 0 = none) or their degree of polarity (with a number from 1 = low polarity to 5 = high polarity). This gives a total of 18 particle types: the Martini building blocks. Such a building block approach-a discrete number of particles which interact using a limited number of interaction levels-provides compatibility of different Martini models and facilitates parameteriza-tion of new molecules, albeit limiting the quantitative accuracy of the force field. In addition to the regular Martini beads, smaller bead sizes (small and tiny beads) are used for groups that are represented at higher resolution such as ring-like fragments. [4,6].
In general, a Martini model for an arbitrary molecule can be generated as follows. (1) The atomistic structure is partitioned into a number of beads, maximizing the four-to-one mapping while preserving the symmetry of the molecule; smaller beads may be used to better represent the geometry of small ring-like fragments. (2) Bead types describing the nonbonded interactions of the models are determined by comparing to already existing fragments or by computing free energies of transfer and selecting the best matching bead. (3) Bonded interactions, defined by a standard set of potential energy functions typical of classical force fields, are parametrized by comparing to atomistic simulations or experimental data. The basic assumption underlying the Martini approach is that the thoroughly parameterized properties of the individual Martini bead types are transferable to the molecule as a whole when linked together to reproduce the overall topology of the desired molecule. This basic assumption entails some pitfalls [26], and hence requires validation, which commonly comes from either comparison to higher resolution atomistic simulations or to experimental data [6].
Below we describe strategies based on the outlined general procedure but which have been found helpful for specific classes of molecules relevant to material systems such as polymers, nanoparticles, and surfaces. These may include specific reference data coming from atomistic simulations or experimental measurements which have been found especially helpful to validate Martini models in this area.

B. Polymers
Martini polymers are applied in a wide range of studies on biomolecular and materials science systems. As first suggested by Rossi et al., parameters for Martini polymers are ideally generated by matching: (1) the free energy of transfer of the monomeric repeat unit; (2) bond and angle distributions of atomistic reference simulations; and (3) long-range structural properties. [27] Overall many carefully parametrized Martini polymers, such as polystyrene (PS) [28], PEO [29], polyethylene (PE) [30], and polypropylene (PP) [30], have been validated by showing that they are able to reproduce a number of single-chain properties. For example, the PS, PEO, as well as PEO-PPO models reproduce structural properties such as radii of gyration in different solvents [28,29,31]. Indeed, validation of these properties in multiple solvents-so as to probe good, bad, and theta solvent conditions-is desirable. Furthermore, persistence lengths, structure factors, and polymer melt density consistute other properties which may be used as validation targets.
An important aspect to keep in mind when parametrizing (Martini) CG polymer models is that the use of torsion angle potentials borrowed from atomistic MD simulations is often unsuitable [32]. Given the softer nature of angle potentials in CG models, conformations which lead to numerical instabilities can be encountered much more often, leading to impractical simulations. Bulacu et al. have devised strategies, like the use of special bonded potentials such as the Restricted Bending Potential (ReB)-implemented in GROMACS [33]-or virtual site-aided definition of bonded terms, to combat this problem [32]. These strategies allow to avoid such instabilities and have been successfully applied to PEO [29,32], PE [30] and other models. Using these strategies, simulations with chains of 500 repeat units over several 10s of microseconds can easily be realized [29].

C. Nanoparticles
Martini models for fullerene [58][59][60], carbon nanotubes (CNTs) [61][62][63][64], graphene [65][66][67] or MXene [68] flakes, clay nanoparticles [69], and functionalized nanoparticles [70][71][72] have also been developed to study their interaction with both other synthetic or biomolecular systems. The C 60 fullerene model developed by Monticelli represents a good parametrization strategy for nanoparticles. The model has been developed by matching experimental free energies of transfer between a wide range of solvents of different polarity and potentials of mean force (PMF) of dimerization in water and octane [59]. This thorough parametrization allowed the Martini fullerene model to reasonably reproduce solid-state properties and lead to translocation PMF across a lipid membrane in good agreement with atomistic reference data. We note that the final model did not use a standard Martini CG bead but instead required some refinements. The thorough refinement of the parameters across a wide range of solvents, however, resulted in a model which could be used in other solvents, such as chlorobenzene, where a comparison to atomistic reference data showed excellent agreement [49] even though the C60 model had not been tested explicitly in chlorobenzene at the time of development. Hence, validation of nanoparticle models in different solvents, by means of comparison to experimental transfer free energies when available and PMFs of dimerization in different solvents is desirable when developing Martini models for these systems.

D. Surfaces
Simulations of Martini systems in materials science have also lead to the development of models for graphite [73,74], graphene [75], and silica [76]. As an example, the parametrization of the graphite model by Gobbo et al. used as reference experimental data enthalpies due to lack of experimental free energy data [74]. Namely, the authors used (1) enthalpies of adsorption of individual molecules from the gas phase on graphite, (2) wetting enthalpies of pure liquids, and (3) enthalpies of displacement of solutes (long-chain organic molecules) from different solvents (heptane and phenyloctane) to graphite. The model required the development of a custom bead representing graphite with a 2-to-1 mapping scheme, which eventually could achieve semi-quantitative reproduction of the experimental enthalpies [74].
Besides the parametrization strategies and validation targets outlined above, several other application-specific critical tests can be carried out to validate a particular Martini model. As we describe the applications in materials science to date in the following sections, we invite the reader to check the reference of interest to find out about further validation targets for the application of interest.

A. Polymeric hydrogels
Polymeric hydrogels are networks composed of hydrophilic polymers that are covalently or physically crosslinked. These polymer networks can swell taking up a multiple of their dry weight in water [77]. They are frequently employed in drug delivery either as nanogel or macroscopic material [77,78]. Martini simulations were employed to understand: (1) interactions of hydrogels with the cargo molecules at molecular detail [79] (2) the gel response to environmental effects such as pH [80]; (3) transport properties of the cargo inside a gel [81], and (4) effect of the salt concentration on the gel [82,83].
For example, Xu and Matysiak have developed a Martini model for chitosan and self-assembled a chitosan hydrogel. Their findings indicate that physical crosslinking patterns impact significantly the hydrogel's mechanical properties. In particular, increasing the polymer concentration or the pH translates into an increase of the elastic modulus of the system, as a consequence of changes in the crosslinking patterns ( Figure 2). In another ex- ample, using a multi-step protocol, protein imprinting of hydrogels was simulated using Martini. Protein imprinting proceeds by polymerizing monomers in the presence of a template protein to which monomers are reversibly coordinated. Following polymerization the gel is washed leaving -so the idea -specific coordination sides for the template protein. Subsequently, by swelling in solution the template protein or alike proteins adsorb more preferentially over random proteins [84]. To mimic this process Zadok and Srebnik first simulated coordination of acrylic Martini monomers with lysozyme and cytochrome c. Subsequently, using a reaction protocol within LAMMPS, monomers were cross-linked to form a gel. After removal of the unreacted monomers a hybrid MD-NVT grand-canonical Monte Carlo simulation was used to swell the hydrogel by allowing the water content to change. Characterizing the interactions of both proteins with different hydrogels, they found that protein binding and selectivity is largely dependent on the nature of the polymer. However, it appeared lysozyme overall has the tendency of forming stronger interactions [79].

B. Polymer coatings and glues
Another area where Martini polymers have been used extensively is to study polymer behavior at surfaces and interfaces. Studies have targeted for example polymer conformations at oil/organic solvent-water interfaces [85][86][87], surface water interfaces [76,88], or even water/air interfaces [31,51]. Polymer behavior at interfaces is interesting in many applications among others for coatings or glues. [76,88,89] For example, Perrin et al. used Martini to study conformations of PDMA and PAAm absorbed to silica surfaces. The aim of the study was to investigate why PDMA glues to silica whereas PAAm does not regardless of their very similar chemical structure. According to their findings, PAAm is better solvated and therefore does not adhere to silica. They further found that dynamical properties of the polymer close to the surface can only be described by an explicit solvent model ( Figure 3) [76]. C. Microphase-separated polymers Although Martini reproduces well the self-assembly of small surfactants and polymers in aqueous solution [34,[90][91][92][93][94][95][96][97], self-assembly of large scale polymer -especially block copolymer systems -remains challenging. These micro-phase separated assemblies of copolymers are hugely important in many technological applications. For example, they are a key component for development of next generation batteries [98]. Hence, simulations of micro-phase separated block copolymers have been performed with Martini [50,[99][100][101][102] For example, in a first attempt, Johnson et al. managed to self-assemble in an unbiased fashion lamellar, micellar and cylindrical phases of a PDMS-γ-benzyl-L-glutamate block-copolymer (Figure 4). Their simulations showed a strong correlation between the obtained morphology and the geometry and type of the side chain [50]. In a similar study, Slimani et al. built a lamellar phase of a polyester gradient coblock-copolymer [100]. These works demonstrate that studying these microphase separated polymer systems is in principle possible with Martini.
In another study concerning block copolymer selfasssembly, Campos-Villalobos et al. simulated PEO-bpoly(butylmethacrylate) (PBMA) copolymers which are being studied for applications as nanostructured materials. [112] Martini simulations of the self-assembly of these block copolymers in water and tetrahydrofuran (THF) mixtures revealed the occurrence of a wide spectrum of mesophases ( Figure 5). The corresponding morphological phase diagram of this ternary system includes dis-persed sheets or disk-like aggregates, and spherical and rod-link vesicles at low block copolymer concentrations, and bicontinuous and lamellar phases at high concentrations. Moreover, the THF/water relative content is found to play a crucial role on the self-assembly kinetics and resulting morphologies [112].
FIG. 5. Morphology obtained from the self-assembly of PEOb-PBMA block copolymers (BCP) in water and THF mixtures [112]. The morphology changes from dissolved chains or monomers in THF, over dispersed sheets or disk-like aggregates, to vesicles as the fraction of water increases. The morphologies are also affected by the BCP concentration [112]. Reproduced with permission from Ref. 112. Copyright (2020) Elsevier.
For example, the behaviour of PEG grafted covalently or physically onto carbon nanotubes has been studied in some detail [131][132][133]135]. In addition to carbon-based nanoparticles, inorganic nanoparticles with grafted polymers are an interesting nanomaterial with potential applications in sensoring, microfluidics and smart surfaces to name some examples. [71] They are especially interesting for their response to different solvents. Within the Martini framework, in particular gold nanoparticles have received a lot of attention [71,85,86,136,138]. For example, Dong and Zhou have studied the solvent behavior of differently composed PEO-b-PS block copolymers attached to gold nanoparticles. They found varying the composition of the block copolymer leads to a variety of different morphologies. Some of them were the expected sphere-shell like conformations where the polymers extend or collapse in a trivial fashion onto the nanoparticles. On the other hand, Dong and Zhou also identified some non-trivial conformations described as rings, buckles, and sectorially arranged chains ( nanoparticles is another interesting effort. Solids formed by such nanoparticles form a versatile class of hybrid materials in which both the nanoparticle core and the organic ligand shell can be tuned, leading to a variety of materials. For example, Martini as been used to simulate single-layer coated nanoparticle membranes [139], or nanoparticle superlattices for energy applications [141]. Additionally, the dispersion of CNTs and 2D materials, such as graphene and MXene, in surfactant aqueous solutions [67,68,75,142] has been investigated. understanding and optimizing such dispersion is of paramount importance for the processing of these materials; in the case of the 2D material MXene, Li et al. probed three typical surfactants with different structural characteristics, finding that the surfactant with long hydrocarbon chain and positively charged head group can form stable bilayers at the surface with MXene, which has implications fot the thermal energy dissipation of the 2D material [68]. Other applications in the field of nanoparticle composite materials include polymer/graphene [66], polymer/graphite [143], polymer/clay [69] composites, and polymer-CNT-protein matrices for applications in the field of tissue regeneration [144].

F. Organic electronics
The morphology of the organic material which constitutes the active layer of organic electronic devices is a critical parameter for the functioning of such devices. Computational modeling of the morphology represents a fundamental step towards an increased rational approach to the design of high-performance organic materials for electronic applications. [145,146] It is possible to model the morphology of organic electronic materials with Martini, in particular to obtain and characterize morphologies, which are often composed of more than one organic semiconductor [49,[147][148][149][150][151][152][153]; and to subsequently backmap [154] the obtained CG morphologies to atomistic resolution, a step often useful in order to perform fine-grained calculations aimed at evaluating the electronic properties of such materials [49,147,[155][156][157]. Martini models have been already developed for many prototypical organic semiconductors used in organic electronic devices, such as conjugated polymers [49,95,158], small conjugated molecules [147,148,159,160], and C 60 fullerene [58,59] and some of its derivatives [49,148,156]. Arguably one of the most popular subfields of organic electronics is organic photovoltaics. Systems such as P3HT:DiPBI [147], P3HT:PCBM [49,[149][150][151][152]161], PBDB-T:F-ITIC [153], and P3HT:PTEG-1 [156] have already been simulated with Martini. Simulations of neat P3HT [155] have also been performed, while more organic semiconductor mixtures have been tested in the context of organic thermoelectric devices [148,159] and organic mixed ion-electron conductors (which will be described as part of the next section) [157,158,[162][163][164].
When modeling the morphology of organic electronic materials, simulating fabrication processes, such as solution-processing and thermal annealing, is an important step to be taken into account and which can be studied to obtain in silico insights. The solvent evaporation process which takes place during the fabrication of organic thin films can be simulated by simulating "bulk" evaporation, as first shown by Lee and Pao using a supra CG model. [165] Similar solvent evaporation simulations have been applied to simulate the prototypical polymer:fullerene photovoltaic blend-P3HT:PCBM-at the Martini level by Alessandri et al. (Figure 7) [49]. Other Martini organic semiconductor thin films have been solution-processed in silico [149-152, 156-158, 161]. In particular, Alessandri et al. studied the evolution of the morphology of P3HT:PCBM blends as a function of the molecular weight of P3HT, the solvent evaporation rate, and thermal annealing. In agreement with experiments, thermal annealing and slower evaporation rates lead to larger phase separation and increased crystallinity of the P3HT phase. The crystallinity of P3HT could be probed by computing scattering signals, which were found to be in qualitative agreement with experimental data [49]. The too-large size of Sbeads, however, prevents a quantitative reproduction of the stacking distance between the polythiophene backbones [49]. Besides allowing to quantify the degree of crystallinity, computing X-ray scattering signals of simulated morphologies allows for comparison to experimental data. Other works made comparison between CG and X-ray scattering data [49,150,157,166], or to scanning electron microscopy [49] or atomic force microscopy (AFM) [148,166] images. Martini simulations allow to scan a range of parameters which are known to affect the morphology of organic thin films in the lab, such as weight ratio of the components [150], polymer polydispersity [151], molecular weight [49,150], post- [49,150] and pre- [150]-evaporation heating treatments. Moreover, once morphologies have been generated, their macroscopic properties, such as mechanical properties [152], or microscopic features, such as the molecular orientations at the interfaces between the two blended organic semiconductors [156] (Figure 7b), can be investigated, possibly as a function of the above parameters. Finally, other works have looked at the solubility of small molecules used as dopants in environments of different polarity [159], another application which is suited to the Martini model.
An important advantage of obtaining morphologies at the Martini level is the possibility of directly backmapping [154] the CG morphologies to atomistic resolution (Figure 7, inset), hence obtaining atom-resolved structures which take into account the self-organization process which occurs during the processing of an organic blend [49]. Indeed, large-scale morphologies have been backmapped to atomistic resolution in order to compute, by means of (semi-empirical) quantum chemical calculations: UV/Vis spectra [147,155], energy levels taking into account the local molecular environment [156], and charge carrier hopping rates [157] for charge transport calculations.
A Martini model exists for the workhorse system of this field: poly(3,4-ethylenedioxythiophene):polystyrene sulfonate (PEDOT:PSS) [162]. Modarresi, Zozoulenko, and co-workers spearheaded Martini-based studies in this field by developing a model for PEDOT and investigating ion diffusion in morphologies of PEDOT:Tos, a system made of PEDOT chains and tosylate (Tos), a negatively charged molecular counterion. [158] The same authors, combining PEDOT with the available [171] PSS model, went on to investigate in detail PEDOT:PSS morphologies [162,164]. The simulations allowed to study the effect of pH on the morphology of in silico solution-processed PEDOT:PSS thin films. Changes in pH were found to greatly affect the morphology (Figure 8), and in turn the distribution of the 5-15 weight % of water content in the polymer film, which is of critical importance for ion diffusion. Once again, the possibility offered by Martini of easily combining mod-els allowed Mehandzhiyski and Zozoulenko to simulate PEDOT:PSS/cellulose [172] composite paper [166]; such paper can be used in applications such as fuel cells, sensors, and batteries, among others [173]. The authors could pinpoint the most likely configuration of PEDOT and PSS/PSSH chains around cellulose by comparing the simulation results to AFM images [166]. They could identify the most likely morphology observed in the experiments, namely a bead-like structure caused by PEDOT aggregates on the fibril which are separated by regions with a lower density. [166] The works presented above on mixed ion-electron conductors partly build on earlier developments of Martini models for polyelectrolytes, which include models for polystyrenesulfonate (PSS) [171,174], poly(diallyldimethylammonium) (PDADMA) [171] and, more recently, partially hydrolyzed polyacrylamide (HPAM) [175]. In the case of such highly charged systems, polarizability introduced for water and other beads has been shown to be important [176][177][178][179][180]. Hence, when charged interactions are expected to be important, it is recommended to develop models in the context of the polarizable [176,177] water model. For example, Vögele et al. have shown that the Martini model of PSS used in polarizable water is able to accurately describe ion distributions around the polymer and reduction in dielectric screening. These important properties are not reproduced in regular Martini water. [171] Furthermore, Martini has been used to model ionconducting polymeric materials used as ion exchange membranes for fuel cell applications [181][182][183][184][185][186][187]. Proton exchange membrane fuel cells use acid polyelectrolytes, such as Nafion, as the membrane material. The membranes are formed by solution processing techniques. In this context, Mabuchi et al. investigated dilute solutions of Nafion ionomers, the initial stage of solution processing: they studied ionomer self-assembly in mixtures of 1-propanol and water and probed the effects of ionomer concentration, alcohol content, and inclusion of salt [184,187]. Goncalves et al. studied instead how cavities nucleate and grow in a hydrated Nafion membrane subject to mechanical deformation, obtaining a nanoscale view on the mechanical properties of such membranes [182]. Next to proton exchange membranes, also alkaline anion-exchange membranes, which transport instead alkaline anions (usually hydroxide), have been modeled with Martini [183,185,186]. Pan et al., for instance, screened for different structural designs of polymer electrolytes which would increase hydroxide mobility, a key performance parameter. The prediction was implemented experimentally and found to lead to increased hydroxide mobility, which reached efficiencies as high as the proton mobility in the more developed Nafion-based proton exchange membranes [183]. Finally, recently Yang et al. used the Martini model to microscopically investigate the impact of adding quaternary ammonium phthalocyanine (Pc) groups into anion exchange membranes based on poly(2,6-dimethyl-1,4-phenylene oxide) for alkaline fuel cells [186]. The self-assembly of Pc helps structuring the anion exchange membrane, increasing its hydroxide conductivity (Figure 9) [186]. FIG. 9. Martini models of anion exchange membranes for fuel cells [186]. Introduction of (a) 10% of quaternary ammonium (QA) phthalocyanine (Pc) groups (Pc-PPO-10) induces a more structured self-assembly than the (b) random morphology formed by Pc-PPO-0 (where QA-Pc groups are not present), leading to enhanced hydroxide (OH − ) conductivity. Here, PPO stands for poly(2,6-dimethyl-1,4-phenylene oxide). Reproduced with permission from Ref. 186. Copyright (2020) Wiley-VCH.

H. Self-assembled supramolecular materials
Self-assembling of molecular building blocks into supramolecular materials holds much promise for a range of potential applications in nanotechnology. [188,189] Molecular building blocks which are very popular are short peptides (2-10 residues) and peptide conjugates, which can give rise to a large variety of biocompatible nanostructures [190]. Many groups have explored their self-assembly process by leveraging the Martini model . Besides allowing to simulate selfassembly and growth, Martini is also particularly suited for high-throughput applications. An example of such a high-throughput application in the area of peptidebased supramolecular materials is the work of Frederix et al. who simulated all 8000 combinations of tripep-tides [200]. The prediction of self-assembling and nonassembling peptide sequences coming from the Martini simulations were verified by a full experimental characterization, showing the predictive power of Martini in this area. The approach allowed to extract guidelines for new peptide materials [200]. The modularity of the Martini model also allows to easily combine peptides with other molecular moieties. Accordingly, Mansbach and Ferguson built models for π-conjugated peptides, which are promising bioelectronic materials due to their optoelectronic properties, and studied their self-assembly in detail [213][214][215][216].
Besides peptide-based compounds, studies on other supramolecular systems have been reported in more recent years.
An important example are 1,3,5benzenetricarboxamide (BTA)-based supramolecular polymers, synthetic supramolecular materials which have been studied in detail by Bochicchio, Pavan, and co-workers [217][218][219][220][221][222][223]. Two Martini models for the BTA core were developed, both capturing the step-wise cooperative polymerization mechanism which leads to the formation of supramolecular fibers ( Figure 10) [217]: BTA monomers initially aggregate quickly in water due to hydrophobic interactions; the disordered aggregates formed then reorganize into directional oligomers on a slower time scale; such ordered oligomers then fuse to form and elongate the supramolecular fiber on an even slower time scale. In the more refined model, additional charged particles were introduced to improve the stacking interactions between the BTA cores [217], similarly to the ones introduced within amino acid side chains by de Jong et al. in the polarizable version of the Martini protein model. [180] The refined model allows for accurate monitoring of hydrogen-bonding between the BTA monomers, while the simpler model-where such charges are omitted-is recommended for studies of interactions of BTA-based assemblies with other (macro)molecules [217]. A Dry Martini [52] version of the BTA model has also been put forward by the same authors [219]. There are several applications of structural variants of the BTA supramolecular polymer [218,[220][221][222], some in combination with enhanced sampling and machine-learning techniques [218,222]. For example, Martini-based well-tempered metadynamics simulations were used to investigate monomer exchange in and out of the fibers, finding a central role of defects on the supramolecular structure in this process [218]. Being able to characterize these defects may thus be important to control the dynamic behavior and properties of such systems: Gasparotto et al. used machine-learning techniques to systematically identify and compare such defects in this class of supramolecular materials [222].

I. Green solvents
Historically, the Martini model has been developed for the simulation of biological membranes formed by lipids. With lipids being one special class of surfactants early on Martini was extended to simulate the assem-bly and interaction of other synthetic ionic and non-ionic surfactants. [34,73,94,[235][236][237][238][239][240] Recently, interest in surfactants has renewed as ionic liquids (ILs) have attracted much attention for their use as bio-compatible and green solvents and co-solvents. This has led several authors to use Martini to simulate the self-assembly of IL mesophases [241,242], the process of IL mediated extractions [242,243] as well as to guide the design of de novo molecules [243]. The use of ionic liquids in applications such as extractions is directly linked to the phase behaviour of the ionic liquid as well as to the emerging phase behaviour when combined with co-solvents. Therefore, understanding and predicting these phases is an important step towards efficient computational solvent design. Martini simulations have very recently been used to unravel the phase behaviour of pure ILs, ILs in water, and mixtures of ILs with other molecules. [119,241,242] For example, Pérez-Sánchez et al. used Martini simulations to understand the temperature dependent effects of adding surface active ILs to Pluronic block-copolymer water mixtures. They were able to elucidate how micelle formation and aggregation changes with Pluronic blockcopolymer composition and type of IL. Their results were well in line with experimental cloud point measurements. This example demonstrates that Martini simulations can be used to study temperature dependence of ionic liquid phase behaviour, at least at a qualitative level. Also other studies have successfully used Martini simulations to investigate temperature dependent effects in the context of ILs. [241,242,244] For example, Crespo et al. investigated the phase diagram of [C n mim][Cl] water mixtures. Their simulations show that these type of ILs display a rich phase behaviour as function of the water content but also temperature in agreement with experimental data where available (Figure 11). The authors also compare the performance of Martini to bottom-up derived CG models. They found that Martini outperforms those models when it comes to transferability from the neat state to mixtures with water. A similar conclusion was reached by Potter et al. who investigated the phase behaviour of a type of a non-ionic chromonic molecule. In their study a direct comparison is made between a bottom-up CG model following the approach of Lu et al. and the top-down Martini approach to CGing. Only with Martini they were able to simulate the complete phase diagram, whereas the structural CG model showed severe limitations at higher concentrations. [240] Whereas the previous studies focused on ILs in cosolvency with water, Vazquez-Salazar et al. used the Martini 3 model to simulate pure [C n mim [CL]]. They observed that Martini well reproduces the system density as function of alky chain length and temperature. In addition, the simulations showed the clear formation of dynamic local organization of the IL, so-called nanodomains, which are an important feature of these type of ILs. Furthermore, the authors demonstrated that the Martini model is able to capture a phase transition of [C 12 mim[CL]] [242]. This phase transition has also been determined experimentally and takes place at around 324.75K. The CG model produced a clear phase transition at 325K, in excellent agreement with experiment.
While the phase behaviour of ILs and their mixtures is an important feature for extractions and applications, it is only indirect evidence for extraction efficiency of a particular IL. Vazquez-Salazar et al. also demonstrated that it is feasible to simulate the extraction process directly by creating a biphasic system of the IL and the solvent phase from which the solute is to be extracted ( Figure 12). Initially all solute molecules are in the solvent phase, but after 6 microseconds of simulation the solutes distribute between the two phases. From the analysis of the solute and solvent density profiles extraction efficiency and selectivity could be computed. Using this protocol, extraction of benzene and poly unsaturated fatty acids from model oil phases was characterized. It was found that the simulations, based on the new Martini3 version, are well in line with the trends observed in experiment.
In a different study, Huet et al. designed de novo ILs, so called zwitter ionic liquids, which contain the anion and cation in the same molecule. These ILs were synthesized as less toxic and more sustainable variant of [C 2 mim][OAc]. The authors used Martini simulations to assess the toxicity of their de novo designed ILs. It was found that the effect of the new ILs on model yeast membranes was less perturbing than for the original IL. Thus it was concluded that the new solvents are less toxic to microorganisms. These conclusions were verified experimentally by computing the minimum inhibitory concentration. [243] This study illustrates how Martini can be used in a more complete design process to also assess toxicity. However, ILs are not the only class of green solvents which can be simulated with Martini. Vainikka and co-workers have used Martini to simulate extraction A B processes with deep eutectic solvents (DES). [246] This class of comparatively new molecules show similar properties to ILs, however, have advantages in terms of cost efficiency and physical properties. [247] J. Oils Several Martini-based applications involving oils can also be found in the literature [248][249][250][251][252][253][254][255], which are of fundamental interest for a variety of industrial processes. An example are studies of the self-assembly of asphaltenes [248][249][250][251], a heavy aromatic fraction of crude oil. Their stability in solution strongly depends on temperature, pressure, and composition, and is important for the petroleum industry. Wang and Ferguson found the introduction of partial charges to represent the radial dipole moment of asphaltenes' aromatic core to be critical to reproduce the T-shaped stacking behavior observed in atomistic simulations [248]. Other studies have investigated the self-assembly and self-organization of various long-chain (functionalized) alkanes on graphite (Figure 13). [74,256,257] The Martini simulations performed by Piskorz et al. provided a microscopic view on the ad-sorption and subsequent rearrangement of alkanes on the surface to form long-range ordered lamellar structures ( Figure 13). The assembly of porphyrin nanorings on graphite has also been explored [258].

A. Martini 3: new opportunities
Recent identification of some of the limits of the current Martini version [6,26], opened the way for the development of a new version, coined Martini 3 [259]. This new version's more general re-parametrization strategy, which did not exclusively include biomolecules, is expected to further boost the application of Martini in soft materials science. Areas in materials science which are particularly expected to benefit from the new reparametrization are applications involving: polymers, which constitute the backbone of soft materials science given their high tunability; conjugated molecules, which are ubiquitous in materials given the possibility of exploiting them as self-assembling systems with interesting (opto)electronic properties; and charged systems, which are important for applications ranging from ionic liquids for green solvents to polyelectrolytes for exchange membranes and next-generation energy storage devices. Moreover, the Martini 3 parametrization has taken into account not only infinite dilution properties such as free energies of transfer but also miscibility data on binary mixtures [259]. As a consequence, the re-calibrated Martini interaction matrix is expected to perform better in applications involving relative miscibility, self-assembly, and aggregation propensities. Additionally, molecular packing is more accurate, as demonstrated for example in a recent biomolecular study where Martini 3 small molecules were able to find and bind to protein pockets in a wide range of systems with very high accuracy [260]. Such results are promising also in view of materials stud-ies were molecular packing is critical: moreover, the improved molecular packing implies that stacking distances between aromatic systems, which were off due to the size of Martini small beads [49], are more accurate in Martini 3. Overall, we anticipate the new version of the model to show improved predictions of molecular packing and interactions in general.
In Martini 2, the need for model refinement sometimes led to the development of custom beads (e.g., see Refs. 28,29,59,184). This need often emerged when parametrizing polymers, where a slight mismatch in the properties of a monomer can build up into relatively large deviations of the macromolecular properties, and in some cases due to a suboptimal parametrization of the smaller bead sizes of the Martini model [26]. However, the development of custom beads is a time-consuming process and can limit the applicability of a Martini model if the bead is not validated properly. The upcoming new version of Martini 3 includes more generic interaction modifiers, generally dubbed "labels". Besides the hydrogenbonding labels already present in Martini 2, which have been expanded and can now [259] be applied to all the new N-and P-bead types, also electron polarizability labels, which mimic the electron-donor or -acceptor character of certain aromatic fragments, and self-interaction labels, which more generically decrease/increase the selfinteraction of a certain bead type without changing its free energy of transfer, were introduced [259]. Such labels expand the capabilities of Martini by giving the user a wider selection of pre-calibrated bead types. Martini users can now use such extra bead types to fine-tune a certain model. Accordingly, we expect the introduction of the more generic interaction modifiers in Martini 3 to greatly reduce the need for custom beads and allow for quick model refinement.
With the improved balance of interactions and possibility of model refinements, we expect Martini 3 to be suited to the description of an even wider range of systems. For example, efforts ongoing in our group are tackling the description of polyelectrolyte complex coacervates, which have material applications in adhesives, coatings, and pharmaceutical applications, aedamers, aromatic molecules which mimic biomolecules as they self-assemble and fold into ordered states, or metalorganic frameworks (MOFs), which are of extreme interest for applications such as hydrogen storage or highcapacity adsorbents for various separation necessities. Although Martini 3 shows numerous improvements, limitations inherent to this CGing approach remain. Of general relevance to both material and biomolecular applications is the limited structural detail due to the CGing process itself. For applications requiring very fine descriptions, atomistic or structure-based CG approaches are more suitable [261,262]. Another limitation of general relevance is the entropy-enthalpy compensation. As the entropy of the system necessarily reduces due to the loss of internal degrees of freedom upon CGing, it is compensated by enthalpy to reproduce free energies. Such entropy-enthalpy imbalance is, for example, known to affect the temperature dependence of several properties, and should therefore be kept in mind. More specifically to materials systems, the description of bare metals, such as the one which may be needed to describe a metallic surface, has not been part of the parametrization and, although not impossible, requires careful validation of the chosen bead types.

B. High-throughput materials design
High-throughput screening of soft matter is an area of immense promise for materials science. [263]. Ideally, a small subset of soft materials would be obtained out of a computational screening procedure so as to speed up and lower the cost of the experimental step. Given the versatility and compatibility of Martini and the efficiency gain with respect to atomistic simulations, Martini simulations are in the position to contribute to the computational design of soft materials. Some applications we envision include: (1) Design of molecular dopants with tailored miscibility: molecular doping is an important strategy used to tune organic semiconductor properties [264]. There are many factors which affect the efficiency of molecular dopants, one of which is miscibility with the host semiconductor [264]. Martini simulations can be used to screen molecular dopants of different polarities for insights in their miscibility with a given host semiconductor. Pushing forward studies such as Refs. 148 and 159, which investigated the miscibility of only few molecular dopants, many dopant designs could be inexpensively explored with Martini and miscibility design rules extracted from such simulations. Such or similar efforts will have to be coupled with a parallel screening of said dopants' electronic properties, which could be obtained by quantum chemical methods. (2) Design of green solvents: the proof-of-concept showcases of benzene and omega-3 fatty acid extractions with Martini 3 [242] show promise for the usefulness of Martini in the computational design of green solvents for selective extraction using ionic liquids. Along the lines of Huet et al. [243], moreover, Martini can be used to design green, bio-compatible solvents which are less toxic and more sustainable. Different co-solvents, structural motifs of the ionic liquid, or mixtures could be computationally screened in order to optimize extraction of certain molecules or other properties of the solvent such as the predicted toxicity. Again, we anticipate the recalibrated Martini 3 interaction matrix, validated by taking into account miscibility data, to be a great asset for such studies. (3) Determining molecular structure-morphology relationships in self-assembling materials: in order to rationally design self-assembling materials for specific applications, one needs to derive robust molecular structure-morphology relationships of the final aggregate or melt. However, the chemical space of organic materials is extremely vast, even if one puts some molecular constraints dictated by the specific ap-plication. Many molecular features are known to affect the final morphology of soft materials, but extracting robust rules is very time and resource intensive with experiments alone. Given the compatibility of Martini models and the ease with which one can vary features such as polarity of molecular features, side chain lengths, molecular topology, we anticipate that Martini simulations aimed at exploring such parameter space in a high-throughput fashion are possible. These efforts, especially if combined with experimental feedback loops and validation, will help reach the overarching goal of establishing robust structure-morphology relationships in materials systems ranging from self-assembling supra molecular materials, over thin films for organic electronics, to nanoparticlepolymer composites.
The above envisioned applications necessarily reflect the own biases of the authors, and of course many more applications can be conceived along these lines and beyond.
To fully harness the compatibility of Martini models and the growing access to computational power, and hence realize high-throughput workflows, the development of tools to automate the Martini workflow will be of key importance. Such tools need to include programs able to automatically create Martini models or Martini building blocks for large molecules. This means they need to design the mapping, assign bead types, and derive bonded interactions from atomistic simulations. Furthermore tools to generate force field files from such building blocks are required as well as tools to create initial coordinates for a variety of target systems. Such tools would not only accelerate the making of Martini models needed for truly high-throughput pipelines, but also make the process more robust and reliable. Endeavors in this direction are already ongoing with tools such as AutoMAR-TINI [265] or Cartographer [266], which are able to derive atomistic-to-Martini mappings, as well as perform a basic bead type assignment. While the proof of concepts are promising more robust ways for assigning bead types and mappings will need to be found. Especially with the release of Martini 3, which includes more bead types, bead sizes, and refined mapping rules. For parametrization of CG bonded interactions from atomistic simulations Py-CGTOOL [267] and Swarm-CG [268] have been designed. Whereas PyCGTOOL is able to efficiently derive bonded interactions for small molecules, it is often impractical for large molecules such as polymers. For example, it cannot recognize and optimize redundant bonded interactions. This deficit is overcome by Swarm-CG, which uses a machine learning-based optimization approach. However, currently generating interactions for medium-sized molecules is still comparatively slow and not all bonded interaction types are implemented. This includes some of those especially important in materials science such as the restricted bending potential or the combined bending and torsion potential . [32] Swarm-CG tool at the time of writing is still actively being optimized on this aspect. Recently, the development of Martinize 2 [269] aims at mapping an entire system from an atomistic reference structure to Martini, generating both target coordinates and input files based on already parametrized fragments. This will allow to more rigorously transform all-atom systems to Martini resolution in one go without having to rebuild the system from the individual components. Together with backwards [154], this will also allow resolution transformation in both directions introducing atomistic detail when needed but also capturing that detail when transforming back to Martini.
Besides model parametrization and resolution transformation, tools to setup starting structures will also be of increasing importance as the complexity of the simulated systems grows. In this context, recently Grünewald and co-workers developed Polyply, a software suite for facilitating atomistic and CG polymer simulations. The tool can generate topologies of polymer systems ranging from simple or complex homopolymers, over branched and hyper-branched polymers, to block copolymers [56]. Not only single chain starting structure but also melts and more complex pre-assembled morphologies can be generated. The latter strategy is useful, because whereas self-assembly might be the preferred strategy to assemble a polymer morphology, it remains a challenging task even at the Martini CG level due to the slow dynamics of long polymer chains. For example, Figure 14 shows a PS melt system with a molecular weight of 1000 residues. Highlighted in gold is a single chain, which winds and twists almost from one edge to the other. The total dimensions of the system is about (30 nm) 3 and comprises half a million CG particles. While properly mixing such a melt from initially disentangled chains is an almost impossible undertaking, Polyply generates such structures within minutes. Polyply's internal library already contains several Martini polymer models from the literature but more models can be contributed via GitHub [56]. This tool is expected to standardize and greatly simplify the generation and setup of Martini polymer systems.
Other tools with a similar purpose are also being developed. For example, CHARMM-GUI already supports generation of Martini membranes [270] and sugars [271], and is currently being extended to polymers. Whereas Polyply and CHARMM-GUI can generate both topologies and structures, PACKMOL [272] is a tool which can be used to generate starting structures from existing molecule coordinates.
Maintenance of a library of available Martini models, and their curation, will also be important on this front.
Many Martini models are available at http://cgmartini.nl: besides an extensive number of biomolecules, and in particular lipids, a growing "polymerdome", and upcoming extensive Martini 3 solvent [259] and small-molecule [273] databases will be important for materials science applications. Moreover, both the upcoming Martinize 2 [269] and Polyply [56,57] tools rely on the Vermouth library [269], which offers a consistent way of defining building blocks and keeping track of model versions, thereby contributing to better FIG. 14.
Melt of Martini P S1000 generated using Polyply [56]. Highlighted in gold is a single chain with other chains (blue spheres) removed within a cut-off of 2.3 nm. The total system size is about (30 nm) 3 . data curation which is an important aspect of growing data sets.

C. Advanced Martini simulations
The growing computational power available and the increase in the complexity of the simulated systems demands for smart and (semi-)automated ways to drive, explore, and analyze the simulated system. Moreover, combining Martini with new method developments can open the way for advanced simulations which go beyond what is feasible with standard CG MD.
Despite the growing computational power available, brute-force sampling of the conformational space often is still not sufficient. Enhanced sampling techniques are available to increase the effective simulation time of MD simulations in general. A host of enhanced sampling techniques already exist and it is in continuous expansion, and many techniques are implemented in packages such as PLUMED [274] and SSAGES [275], which are compatible with many MD softwares, including GRO-MACS [33] and NAMD [276], and hence can be readily applied to Martini systems. Given the insights such techniques have already provided when applied to Martini simulations [218,224] and the active developments in this area, Martini simulations will surely benefit from coupling to such techniques.
Changes in bonded interactions and atom types, which reflect chemical reactions, cannot be captured by regular MD simulations. However, chemical reactions especially in the field of material science are ubiquitous. Examples include cross-linking reactions in polymers, dynamic changes in the protonation states of polyelectrolytes, or formation of gels. Efforts are being taken towards capturing these effects. [79,88,277,278] Rossi et al. studied the cross-linking of a polyester resin. To model a cross-linking reaction they used an ad-hoc empirical potential, which displayed harmonic-like features at close distance, but was able to dissociate at larger distance. In contrast, Zadok and Srebnik used a LAMMPS built-in feature to simulate the effects of gel formation for a protein imprinted hydrogel. Ghermezcheshme et al. used a similar approach, however, combining GROMACS with an in-house code to study the step-growth reaction of polyurethane. Here a bond is introduced after finding all reactive neighbors within a cut-off. Upon changing the topology a long equilibrium simulation is carried out after which another step of cross-linking is performed. In contrast, the recently developed titratable Martini [278] allows the dynamic representation of protonation reactions in a Martini simulations. Here the protonation state of a titratable functional group can change back and forth between protonated and deprotonated during the course of a continuous simulation. Using this approach for example, the pH induced collapse of the hyper-branched polymer poly(propylene imine) was simulated as shown in Figure 15. Like the approach by Rossi an empirical nonbonded potential is used to allow protons to tightly bind to titratable functional groups. Further progress in the field of reaction simulations is expected with the inclusion of lambda dynamics into GROMACS [33]. Protonation of a poly(propylene imine) (PPI) dendrimer as a function of pH as modelled with the titratable version of Martini 278. The protonation state of the core beads of the dendrimer, which represent tertiary amines, clearly change as function of pH becoming progressively less protonated as indicated with the color scale from red (protonated most of the time) to blue (deprotonated most of the time). The radius of gyration (Rg) quantifies the degree of polymer collapse as the charge density decreases at higher pH. Reproduced with permission from Ref. 278. Copyright (2020) American Institute of Physics.
Machine learning (ML) approaches coupled to molecular modeling are emerging in soft materials research [279,280]. Such approaches have the potential to support and augment traditional physics-based models in computational research [281]. In the realm of soft materials, ML techniques have been shown to be able to predict electronic properties, such as energy levels and absorption spectra, of soft materials directly from CG structures [282,283]. This strategy is a potentially quicker and less laborious alternative to the currently necessary backmapping and subsequent quantum chemical calculation step, and it is especially relevant to systems where Martini structures are usually backmapped to obtain electronic properties, such as organic electronic systems [155][156][157]. Another interesting avenue is the one of coupling CG simulations and ML techniques to efficiently explore the desired chemical space. For example, Shmilovic and co-workers [216] used Martini simulations within an active learning strategy to efficiently cover the chemical space span by π-conjugated peptides with a πcore flanked by two tripeptide units. Instead of simulating all possible tripeptide sequences (20 3 = 8000), the active learning strategy allowed to minimize the number of simulation data needed for training the ML model. Accordingly, by direct simulations of only 2.3% of the tripeptide space, the authors showed how a Gaussian process regression model could capture the statistical information necessary to represent this chemical space [216]. The model could then be leveraged to identify the πconjugated peptides predicted to exhibit superior assem-bly properties to those reported in previous work. Moreover, in this way the authors were able to reveal design rules governing assembly of these molecules. The fact that the Martini bead types discretize the chemical space, means that similar molecules will often map to the same Martini CG model. This introduces a degeneracy in the CG representation, which translates into a reduction of the size of the chemical (compound) space. [284] Accordingly, this reduction of the size of the chemical space represents also a further speed-up which can help for screening studies. [284] Hence, we expect hybrid Martini/machine-learning schemes to be highly promising in order to efficiently explore the chemical space for different applications.
In conclusion, we foresee a bright perspective for applications of the Martini model in the field of materials science, especially given the possibilities offered by the upcoming new version of the force field, existing and forthcoming tools to streamline model building and system preparation, and combination with existing and future method developments, such as constant pH simulations and hybrid simulation/machine learning schemes.