Coarse-grained conformational surface hopping: Methodology and transferability

Coarse-grained (CG) conformational surface hopping (SH) adapts the concept of multisurface dynamics, initially developed to describe electronic transitions in chemical reactions, to accurately describe classical molecular dynamics at a reduced level. The SH scheme couples distinct conformational basins (states), each described by its own force field (surface), resulting in a significant improvement of the approximation to the many-body potential of mean force [Phys. Rev. Lett. 121, 256002 (2018)]. The present study first describes CG SH in more detail, through both a toy model and a three-bead model of hexane. We further extend the methodology to non-bonded interactions and report its impact on liquid properties. Finally, we investigate the transferability of the surfaces to distinct systems and thermodynamic state points, through a simple tuning of the state probabilities. In particular, applications to variations in temperature and chemical composition show excellent agreement with reference atomistic calculations, introducing a promising"weak-transferability regime,"where CG force fields can be shared across thermodynamic and chemical neighborhoods.


I. INTRODUCTION
In the realm of multiscale models for soft matter and biomolecular systems, particle-based coarse-grained (CG) resolutions have offered tremendous insight. [1][2][3][4][5][6][7][8][9] CG models average over the faster degrees of freedom by lumping several atoms into super-particles or beads. When adequately built and parametrized, these models can strike an excellent balance between accuracy and computational efficiency. Their success stems largely from a mapping commensurate with the system's scale separation and an adequate use of physics-based modeling. The latter aspect is the main topic of this study.
Coarse-graining replaces the coveted potential-energy surface (PES) by a configuration-dependent free-energy function known as the many-body potential of mean force (MB-PMF). 10,11 Over the last two decades the community has been developing and improving a number of systematic methods aimed at targeting the MB-PMF. [11][12][13][14][15][16][17] While early efforts established a strong theoretical and practical foundation for these methods, a number of fundamental challenges have arisen, which largely prevent a more widespread utilization of systematic (i.e., bottom-up) CG models. 9,18 Transferability-the capability of a given model to be accurately applied to systems and thermodynamic state points distinct from those used for parametrization-is an intrinsic problem for coarsegraining, since the MB-PMF is inherently state-point dependent. 9,19 As a consequence, there has been a continued effort to systematically investigate the temperature, density, and solvent-mixture transferability properties of CG models. [20][21][22][23][24][25][26] In limited cases, it has been demona) Electronic mail: rudzinski@mpip-mainz.mpg.de b) Electronic mail: t.bereau@uva.nl strated that CG interactions can reproduce temperaturedependence of liquid structure through an ad-hoc linear interpolation, 27,28 although a systematic approach has been lacking. Recent work has begun to fill this gap either through Bayesian techniques 29 or approaches that approximate the entropic contributions to the effective potentials, allowing explicit predictions of state-point dependence. [30][31][32] These studies have focused on the CG representations without significant intramolecular flexibility. Beyond thermodynamic-state-point dependence, few studies have reported detailed characterizations of the chemical transferability of bottom-up CG models. [33][34][35] Even for a single system or thermodynamic state point, persistent efforts have not led to steady improvements in the quality of the force fields-the accuracy being limited less by the performance of the methods, and more by the molecular-mechanics terms used to approximate the MB-PMF. Because these terms only offer an incomplete projection of the full MB-PMF, a CG model's accuracy critically depends on two aspects: (i) an optimized mapping that most effectively simplifies the form of the MB-PMF, [36][37][38][39][40] and (ii) interaction-potential forms that are flexible enough to describe complex physical phenomena, such as interfaces or environment-dependent conformational changes. 41,42 Going beyond the typical interaction terms-especially non-bonded pairwise interactions-can have significant impact, as seen by recent investigations that considered physics beyond pairs, such as threebody interactions 43,44 or local density-dependent potentials. [45][46][47][48][49] However, these approaches are also limited by the functional forms applied, and it may not be obvious how to generalize them for optimal improvement in modeling accuracy. Recent applications using machine learning can provide a more accurate reproduction of the MB-PMF either through a multi-body decomposition or by a direct interpolation of the many-body forces. [50][51][52][53][54] This improved accuracy typically comes with added computa-tional cost-a significantly larger evaluation time needed for the force prediction, 55 unless projected onto tabulated potentials. 54 We recently introduced a complementary approach to improve the description of cross-correlations between interaction terms in a force field. 56 This approach was inspired by the modeling of chemical reactions, where distinct electronic configurations are decomposed onto separate surfaces in order to overcome limitations of the force field by coupling distinct PESs-notable examples include (multisurface) empirical valence bond and surfacehopping schemes. [57][58][59] Instead of describing transitions between electronic states, our method considers switching between conformational basins: Distinct force fields describe interactions for a subset of conformational space. While several approaches have been devised to couple internal states in various ways in the context of classical molecular simulations, 60-63 explicit hopping schemes have been avoided through simple approaches that either (i) linearly interpolate between two force fields (e.g., multi-state Gō models) or (ii) describe the forcefield change as an analytic function of a continuous order parameter (e.g., local density-dependent potentials). In contrast to these approaches, our surface-hopping method allows for a more general representation of the distinct force fields describing interactions within particular conformational basins. Voth and coworkers have proposed a so-called "ultra coarse-graining" methodology to switch stochastically and discretely between internal states. 64, 65 Sharp et al. extended the surface-hopping idea (within the ultra coarse-graining framework) to describe transitions between conformational basins defined along a set of collective variables. 66 By identifying conformational basins according to the internal CG degrees of freedom, our surface-hopping method effectively couples local features along the MB-PMF and describes transitions between conformation states that occur continuously on time scales comparable to the molecular dynamics. 56 Coupling interaction terms of the Hamiltonian offers the possibility to rescue cross correlations beyond the typical global separation of variables. By focusing on the coupling of intramolecular interactions, we reported significant improvements in the accuracy of the MB-PMF for a three-bead model of hexane, as compared to the baseline force-matching-based multiscale-coarsegraining (MS-CG) method. Furthermore, the surfacehopping model for a tetra-alanine peptide in water not only showed a virtually exact reproduction of a twodimensional projection of the MB-PMF but also demonstrated excellent agreement in the ratio of mean-first passage times between helical and extended states. The latter is significant: it shows that a faithful representation of the MB-PMF can offer an accurate reproduction of the barrier-crossing dynamics up to a speedup factor. While equilibrium properties depend exponentially on the freeenergy minima, an accurate reproduction of the barriercrossing dynamics critically depend on the free-energy barriers. [67][68][69] The present report extends our previous work in several ways. We first provide a more detailed account of the methodology starting from a toy example-a single particle in a double-well potential. Next, we extend the methodology to non-bonded pairwise interactions and report results on liquid properties. Finally, we investigate "weak-transferability" properties, corresponding to the transfer of surfaces while specifically tuning the state probabilities. Applications to both temperature and chemical composition show excellent agreement for a series of small molecules.

A. Intramolecular interactions
We recall the example presented in our previous publication: 56 We consider a two-dimensional potential U = U (x, y) with corresponding canonical equilibrium distribution p = p(x, y) ∝ exp (−βU (x, y)), where β = (k B T ) −1 is the inverse temperature. Standard molecularmechanics force fields apply a global separation of variables on the potential, such that U (x, y) ≈ U x (x)+U y (y), also impacting the equilibrium distribution: p(x, y) ≈ p x (x)p y (y). As a result, we cannot ensure an accurate reproduction of cross correlations between x and y. For instance, the intramolecular interactions of a 3-particle, linear molecule made of two bonds, b 1 and b 2 , and one bending angle θ will typically be modeled by a potential of the form U (b 1 , b 2 , θ) = U b1 (b 1 )+U b2 (b 2 )+U θ (θ). While significantly advantageous from a computational standpoint, the separation of variables can drastically hamper the accuracy of the (free-)energy landscape. Fig. 1a illustrates the potential issues of such an approach. In particular, if there exists two local minima along each degree of freedom, a model which employs the global separation of variables will likely sample all four combinations of these minima, regardless of the true underlying distribution.
The conformational surface-hopping (SH) scheme retains the same form of the Hamiltonian, as well as the separation of variables, but ascribes a local force field for a subset of conformational space-a conformational basin, say. In the case of two surfaces, the SH equilibrium distribution takes the form p(x, y) ≈ p (1) x (x)p (1) y (y) + p (2) x (x)p (2) y (y), allowing the description of a wider range of cross correlations between the degrees of freedom (see Fig. 1b). This prescription trivially generalizes to n surfaces. An SH force-field parametrization thereby consists of the following steps: 1. A clustering of (intramolecular) conformational space is performed (here with respect to variables x and y) to identify homogeneous regions, ideally leading to unimodal one-dimensional distribution functions along each (intramolecular) degree of freedom. Each cluster is assigned a center, , corresponding to the local maximum of probability density, and a spatial extent, , related to the standard deviation of configurations belonging to the cluster.

2.
A linear transformation is applied to the conformational space in order to enhance the isotropy of the clusters: 3. n − 1 surfaces are defined according to the clustering, while an additional surface is introduced which covers the remaining configurations. This surface will be referred to as the "fallback" surface. 4. A structure-based parametrization of n force fields is performed (e.g., via force matching), one for each surface.
Each SH force field, f i (R) = −∇U i (R), is related to a typical molecular mechanics potential, U i (R), which employs a global separation of variables. In the SH method, the net force for any configuration of the system can be written as a linear combination of the individual force fields: where the coefficients or weights are restricted to 0 < w i < 1. Force field i will contribute to the net force according to the proximity of the system's instantaneous configuration to cluster i. Practically w i is computed as a Euclidean distance of the system's CG interaction variables (x, y) to the cluster center (2) d i is then compared to the spatial extent of the cluster |σ (i) |. When d i < |σ (i) |, the system is completely within cluster i and its force field receives the full weight, w i = 1, while all other force fields are neglected. In the case that the system's configuration is not inside one of the clusters, the SH approach will connect surfaces together to ensure a smooth hopping. To this end, the force field weight is exponentially suppressed with respect to the distance from the boundary of the cluster, d i − |σ (i) |: The sharpness of this suppression is determined via the scaling parameter α, which can assist in avoiding numerical instabilities in the simulations. On the other hand, we stress the importance of keeping α small, as it blurs the force-field boundaries.
Mixing different force fields can lead to unphysical behavior, for instance if the aggregate contributions yield large net forces. This is especially relevant at the boundaries between conformational basins, where a localized force field will have large restoring forces at the boundaries (see, for instance, panels a and b of Fig. 2). We hinder this behavior by restricting mixing to occur between only two force fields: one corresponding to the closest cluster and one corresponding to the fallback surface. More specifically, we first compute the initial w i for each of the first n − 1 surfaces according to Eq. 3. The largest weight, w l = max i<n w i , is kept, while the remaining weights are set to zero. Then, the final weight is assigned to the fallback surface: w n = 1 − w l . This approach assumes that the fallback surface is well connected to all of the surfaces, and thereby be well-defined broadly across the conformational space of the system.
The algorithm so far leads to surface hopping but does not ensure the correct probability of sampling each conformational basin. To this end, we enforce that the time average of the probabilities to be within each state roughly matches a set of target reference probabilities, available upon partitioning conformational space. This approach, both simple to implement and effective, is described in more detail in our previous work, 56 as well as below.

B. Intermolecular interactions
Having described an SH model that switches between force fields according to the order parameters governing intramolecular interactions, we now turn to the treatment of intermolecular interactions. In this work, the intermolecular interactions rely on the SH state definition, determined by the intramolecular order parameters, which effectively couples the two types of interactions. However, the local (non-bonded) environment of each molecule does not play a role in defining the SH state. For instance, consider two particles of type A and B belonging to distinct molecules. For each particle we compute their most contributing surface-a function of, for example, the bond distances and bending angles of each molecule. Let these surfaces be denoted j and k for particle types A and B, respectively. Then, the resulting pairwise non-bonded interaction between these particles will not only depend on the pair of particle types A-B, but also on the combination j-k. The parametrization of the intermolecular interactions consists of appropriate filtering of the reference trajectory: we gather statistics between particles A and B that also have internal state j and k, respectively. Computationally, the non-bonded interaction switches nearly instantaneously according to the pair of internal states, as defined by the bonded interactions. The relatively small difference between nonbonded potentials helps avoid numerical instabilities.

III. COMPUTATIONAL METHODS
The protocol applied here largely follows our previous study. 56

A. All-atom simulations
In this work, we consider four small molecules: hexane, octane, hexanediamine, and hexanediol. For each molecule, we performed simulations of (i) a single molecule in vacuum and (ii) 267 molecules in the liquid phase at various temperatures. These simulations employed the OPLS-AA force field 70 to model interactions and were performed with the Gromacs 4.5.3 simulation suite 71 according to standard procedures, described in more detail in the Supporting Information.

B. Coarse-grained representation and interactions
For hexane, we considered a 3-site representation, which represents subsequent pairs of carbon atoms with a CG site. The CG potential included two identical bonded interactions between subsequent pairs of sites along the chain and an angle bending interaction between the three CG sites. This representation and set of interactions has been applied in several previous studies. 36,56,72 Octane, hexanediamine, and hexanediol can be considered "extensions" to the hexane molecule, through the addition of a functional group on each end of the molecule. To assess to what extent the SH state definitions can be transfered between molecules, we employed the hexane mapping and interaction set to these other three molecules.
That is, each pair of carbon atoms were represented by a CG site, while the terminal functional groups were not explicitly represented (see Fig. 7(d)). For each of these 3-site models, the terminal CG sites were represented by identical types, denoted as CT. The center CG site was represented with a distinct type, denoted CM. Distinct pairwise non-bonded interactions were then employed between each unique pair of bead types.

C. Partitioning of conformational space
To obtain the SH state definitions, we performed a density-based clustering analysis 73 to the atomistic trajectories of single molecules in vacuum, after mapping each configuration to the CG representation. This clustering analysis was performed along the order parameters governing the intramolecular CG interactions, i.e., the two bond distances and bending angle. Before clustering, these intramolecular order parameters were transformed to mean-centered and normalized values for regularity. The clustering used a search radius R = 0.1. The initial clusters were grouped into coarser states manually via visualization of the cluster distributions along each order parameter, although an automated dynamics-based algorithm 74 yielded similar results.
For the three-bead representations of both hexane and octane, the clustering resulted in a set of 7 clusters, representing different combinations of bond and bending-angle values. This result already indicates some consistency between the intramolecular states sampled by these distinct molecules. In the following, we will consider a 3-state SH model, where the two most populated clusters determined the states denoted 3S-1 and 3S-2, while the rest of the configurations were lumped into the fallback surface (3S-3). In the Results section, we first assess the properties of the SH model for hexane, both in vacuum and in the liquid phase, using the state definitions determined from the AA simulations of hexane. Subsequently, the transferability of the SH state definitions across chemistry is assessed. For this investigation, a single set of state definitions, determined from the AA simulations of octane, were applied for each molecule.

D. Generation of the coarse-grained potentials
In this work, all CG force fields are derived using the framework of the force-matching-based multiscale coarsegraining (MS-CG) method. The MS-CG method approximates the MB-PMF via a mathematical projection of the many-body mean force, i.e., the negative gradient of the MB-PMF, into the space of force fields spanned by the chosen basis set representation for the CG force field. 11,12 This corresponds to matching the average force on each CG particle sampled in the simulation of the underlying, higher-resolution model. Practically, the projection is expressed as a linear least squares problem in the basis function coefficients, i.e., the CG force-field parameters, φ, and can be written in the normal equation representation as In Eq. 4, b AA is a vector of ensemble averages that can be expressed as a set of either force 11,75 or structural 16,76,77 correlation functions. The latter is possible through a generalized Yvon-Born-Green (g-YBG) framework, which connects the MS-CG method to traditional liquid state theory. 16,76 For a non-bonded, pairwise interaction represented by a set of spline basis functions, b AA is directly related to the corresponding radial distribution function generated by the reference model but mapped to the CG representation. 77 G AA is a matrix that quantifies the cross correlations between pairs of CG degrees of freedom generated by the reference model. If the model derived from the MS-CG method fails to reproduce the target vector of these equations, i.e., b AA , it implies that the cross-correlation matrix generated by the higher resolution model does not accurately represent the correlations that would be generated by the resulting CG model. This indicates a fundamental limitation of the model representation and interaction set. Nevertheless, the system of equations can be solved self-consistently to determine the force field φ * that reproduces the target vector, albeit at the expense of the representation of the cross correlations of the underlying model 36 : This approach has been previously denoted as an iterative g-YBG (iter-gYBG) method. 36,78,79 In the following, we consider both φ MS-CG and φ * (denoted as the iter-gYBG model), simulated according to standard techniques, for comparison with the SH simulation method described above.
Using the partitioning of configuration space described above, we also determine sub-ensemble-specific CG potentials by solving Eq. 4, but employing trajectories containing only configurations from a specific sub-ensemble to calculate each of the correlation functions. Although structure-based methods are often applied to conformational ensembles at equilibrium, several studies have demonstrated the benefit of performing parametrizations over sub-ensembles or biased ensembles. 17,80,81 The formal theory for such calculations in the context of the MS-CG method has been detailed by Voth and coworkers. 64,65 All force-field calculations in this work were performed using the BOCS package. 82 Further numerical details for these calculations are provided in the Supporting Information.

E. Coarse-grained simulations
We performed CG simulations of the SH models using a modified version of ESPResSo++. 83 Simulations in the canonical (N V T ) ensemble were performed using a Langevin thermostat at various temperatures (more details below and in the SI), where a friction constant Γ = 10 τ −1 , was applied. Here, τ corresponds to the intrinsic unit of time of the CG model. We integrated the equations of motion with a time-step δt = 0.001 τ . All cluster sizes {σ i } were scaled by a factor of 0.4 to significantly localize each surface. The smoothness scaling parameter was set to a small value, α = 0.05, to ensure numerical stability of the dynamics while minimally distorting the individual force fields. An ESPResSo++ implementation of the CG surface-hopping scheme, including support for non-bonded interaction, is available online. 84 We performed CG simulations of the MS-CG and iter-gYBG models using version 4.5.3 of the GROMACS package, 71 according to standard procedures (see Supporting Information).

A. Toy model
We first illustrate the method using a toy model: a single particle dynamically evolving in a one-dimensional double-well potential. 85 The potential U (x), which we will refer to as the global surface, is shown in Fig. 2a, while panel (b) displays the corresponding force −∇U (x). Throughout this section we express the results in the natural units of the model: energy ( ), length (σ), mass (M), and time (T = σ M/ ). The system is simulated using Brownian dynamics according to the stochastic differential equation where R(t) is a white-noise process, T = /k B is the temperature of the system, and D = 10 2 σ 2 /T is the diffusion constant. We use an integration time-step δt = 10 −7 T . Integration of Eq. 6 leads to a time trajectory of the coordinate, x(t) in Fig. 2c, and an equilibrium distribution, P eq (x) in Fig. 2e, featuring the expected two peaks. We now turn to a surface-hopping model of the system. We split the global surface into two components: a surface corresponding to the global minimum S 1 (violet curves in panels (a) and (b) of Fig. 2) and a higher-energy surface S 2 (orange curves). Two distinct potentials are fitted to best reproduce the local basins of U (x) within a harmonic approximation. The resulting energy functions, and also the corresponding force curves, show high fidelity to all parts of the global surface except around the barrier (x ≈ 0.25 σ). We connect the two surfaces S 1 and S 2 by means of an instantaneous switching at x = 0.25, leading to a discontinuity in the force (red dashed line in Fig. 2b). This generates a cusp in the  potential energy, leading to inaccuracies in the shape of the potential energy around the barrier. A straightforward integration of the equations of motion of this surface-hopping model (denoted "SH-noprob") qualitatively samples the two surfaces by regularly switching between them (Fig. 2d), but leads to noticeable discrepancies in the equilibrium distribution. Fig. 2e demonstrates that the SH-noprob simulation slightly overpopulates S 1 . This overrepresented sampling of S 1 is clearly displayed in the time evolution of the probability of that surface (Fig. 2f), which converges to around P S1 (t → ∞) = 0.70 instead of 0.65.
A correction to the inaccurate representation of the barrier can be obtained by enforcing the probability of sampling S 1 . To this end, we restrict the hopping between surfaces by adjusting the force interpolation scheme based on the instantaneous time average of the probability of sampling each surface in the simulation. More specifically, once the system completely enters a cluster, the weight given to the corresponding force field is fixed to be 1 until the probability of sampling the cluster exceeds a given target probability. 56 The surfacehopping simulations with this restriction, denoted simply "SH," converges by construction to the target probability P S1 (t → ∞) = 0.65, leading to an improved description of equilibrium distribution (Fig. 2e). Thus, enforcing the target probabilities mitigates potential issues due to an inaccurate modeling of the boundaries between surfaces. We emphasize the need for a small interpolation regime between surfaces: too large of a region would lead to the inclusion of unreasonably large forces from the less favored surface, resulting in artifacts at the interface. We also note that an alternative approach could consist of interpolating between potential energies, although this would require to shift each surface by an appropriate amount.
To further probe the dynamics, Fig. 3 presents the probability distribution of escape times between basins S 1 and S 2 . Assuming single-exponential kinetics, we focus on the characteristic time scales of the forward and backward processes, k S1→S2 and k S2→S1 , respectively. While the integration of the global surface and the surface-hopping surfaces (panels a and b, respectively) lead to different characteristic time scales, their ratio are in good agreement with one another: k S1→S2 /k S2→S1 ≈ 1.47 and 1.55, respectively. This is on par with our previous conclusions about the method's capability to conserve the barrier-crossing dynamics, as illustrated on a tetraalanine peptide. 56

B. Hexane
In the following we consider the coarse-graining of hexane to a three-bead representation. We first simulate a single molecule in vacuum, effectively focusing on the intramolecular interactions. Later we turn to inter molecular interactions by probing the liquid state.
3. Probability distribution of escape times of the toy model between basins S1 and S2 for (a) the global surface and (b) surface hopping. While the characteristic time scales are different between the models, the ratio of characteristic time scales is conserved.

Hexane in vacuum
The modeling of hexane in vacuum using a 3-site CG representation, though presumably straightforward at first sight, displays remarkably rich cross correlations between the bond and bending angle degrees of freedom. This offers a stringent test for molecular mechanics force fields. The system was first studied by Rühle et al. 72 using the force-matching based multiscale coarsegraining (MS-CG) method and later by Rudzinski and Noid, focusing on the cross correlations and presenting results based on the iterative generalized Yvon-Born-Green (iter-gYBG) scheme. 36 Some of the analysis presented here was described in previous work, 56 although the present work provides additional details and uses the previous analysis as a basis to dive further into various features of the method.
To build the SH model of hexane, we first partitioned the conformational space defined by the two order parameters governing CG interactions: bond, b, and bending angle, θ. The torsional degrees of freedom at the atomistic level give rise to a bimodal distributions of CG bond distances and an approximately trimodal distribution of angles (violet curves in Fig. 4). The angle distribution also displays a tail at short distances, which corresponds to a partially hidden fourth mode, described further below. By separating each order pa-rameter into distinct states based on these distributions, the intramolecular state of the molecule can be described as a discrete set of two bond states and an angle state. The AA model then samples approximately six unique intramolecular states, with varying equilibrium probabilities. 36 The surface-hopping model simplifies this description with a 3-state representation for the intramolecular configuration of the hexane molecule. This leads to the definition of three surfaces denoted 3S-1, 3S-2, and 3S-3, which we will characterize below. Notably, an analysis of the reference AA simulation provides the probability of sampling each surface: 0.45, 0.14, and 0.41, respectively.    Fig. 4 show that the MS-CG model is capable of reproducing the bond distribution, characterized by a short bond (b ≈ 0.24 nm) and a long bond (b ≈ 0.26 nm). The 3S model generates essentially the same distribution, interestingly using a rich combination of bond potentials (Fig. 4a). While 3S-1 is dedicated to describing the long bond, 3S-2 is shifted to values that are even larger-the small probability of sampling this surface leads to a virtually negligible impact on the bond probability distribution. Lastly, 3S-3 describes both the long bond-with a basin aligned with 3S-1-and a short bond that is energetically offset. This surface alone is responsible for the smaller peak in Fig. 4b.
Turning to the bending-angle potential and probability distribution (panels c and d of Fig. 4), the MS-CG model displays severe discrepancies: it significantly under-samples the two larger angles (θ ≈ 170 • and θ ≈ 155 • ) and over-stabilizes the two lower angles (θ ≈ 105 • and θ ≈ 125 • ). This discrepancy has been demonstrated to be due to complex AA cross correlations between the bond and angle degrees of freedom, which are used as a proxy for CG correlations within the MS-CG procedure. 36 Unlike the MS-CG model, the iter-gYBG model presented by Rudzinski and Noid is capable of reproducing the one-dimensional distribution function p θ (θ), but does not accurately reproduce the cross correlations p(b, θ). 56 The 3S model also matches the AA bendingangle distribution nearly quantitatively, but describes the sub-populations of the distribution with greater detail through the multi-surface representation. The 3S-1 surface focuses solely on the largest-angle state, while 3S-2 focuses on the two intermediate angles. We note that despite the predominance of the lower intermediate angle (θ ≈ 125 • ) within 3S-2, the higher intermediate angle displays a higher population due to 3S-3. The 3S-3 fallback surface does not target a particular conformational basin, but instead broadly covers the entire dynamic range of populated angles with various weights.
The major improvements of the SH model can be seen through the cross-correlations, namely the free-energy surface −k B T ln p(b, θ), displayed in Fig. 5. Compared to the reference distribution (panel c), the CG iter-gYBG model (panel a) yields exceedingly symmetric features, illustrative of the additivity of the interactions. On the other hand, the 3-state SH model displays a much more accurate free-energy surface. 56

Liquid hexane
We now turn to assessing the capabilities of the SH models to describe liquid properties. As a test system we focus on a homogeneous bulk liquid of hexane, comprised of 267 molecules in a cubic box of size L = 3.89 nm simulated at T = 300 K. In principle, the surface definitions could be extended to depend on additional order parameters, e.g., as a function of local density. However, since the benefit of local density-dependent potentials has already been characterized by others, [45][46][47][48][49] here we focus on the extent to which a more accurate treatment of intramolecular structure impacts the resulting properties of the liquid. To this aim, we employ the surface definitions derived from the vacuum case, described above. While the surface definition does not depend on the intermolec- ular environment, the intermolecular interactions do depend on the intramolecular state of the molecule. That is, we calculate distinct pairwise interactions as a function not only of the set of bead types but also as a function of the surface of each molecule. As an illustrative example, we focus on the interactions between terminal beads, CT, but additional results for the other interactions can be found in the Supporting Information. Fig. 6 (a-c) shows a comparison of the pairwise interaction potential U (r CT-CT ) for the MS-CG, iter-gYBG, and SH models. The MS-CG and iter-gYBG potentials are purely repulsive and only very weakly attractive, re- spectively, consistent with a variety of work which has demonstrated that structure-based methods tend to underestimate the cohesive energy of liquids. 18,[86][87][88] Interestingly, while the SH potentials do include some more repulsive interactions on par with the MS-CG and iter-gYBG potentials, there are also some significantly more attractive interactions. The 3-3 SH potential between fallback surfaces roughly resembles the MS-CG potential, although it displays an additional small distant attractive basin around r CT-CT ≈ 0.8 nm. On the other hand, the 1-3 and 2-3 SH potentials show a dip that is akin to the iter-gYBG model, with a depth of about 0.5 kJ/mol, albeit without a large barrier around r ≈ 0.7 nm. This is quite striking, since such barriers and secondary potential minima have been associated to a type of over-fitting that occurs in structure-based models. 89 Further, the 1-1, 1-2, and 2-2 SH potentials show a deeper minimum (1.3 kJ/mol). Critically, we emphasize that there is a clear clustering of the set of SH potentials into three families: (i) the interaction between two molecules both in the fallback surface (3-3), (ii) interactions of a molecule in the fallback surface with a molecule in one of the other two surfaces (1-3 and 2-3), and (iii) interactions between two molecules not in the fallback surface (1-1, 1-2, and 2-2). Natural groupings such as these provide a clear strategy for addressing the combinatorial explosion of the SH framework as the number of surfaces and bead types increase. Fig. 6d presents the CT-CT radial distribution functions (RDFs) generated by the various models, demonstrating that calculating the pairwise interactions as a function of intramolecular state within the forcematching framework is robust (i.e., does not result in errors in the structural properties). In fact, the SH model actually demonstrates an improvement with respect to the MS-CG model, which shows small deviations after the first solvation shell. These deviations are, of course, at least partially associated with the inaccurate determination of the MS-CG angle potential. Fig. 6 (e-h) further characterizes the structural properties of the models via the cross correlations between the bending angle and pairwise distance r CT-CT . These cross correlations correspond to sub-blocks of the correlation matrix employed in Eqs. 4 and 5, and are described in detail elsewhere. 90 Compared to the AA cross correlations (panel h), the cross correlations generated by the MS-CG model (panel e) demonstrate significant discrepancies, largely due to the inaccurate description of the bending-angle distribution (Fig. 4d). In contrast, the cross correlations generated by the iter-gYBG and SH models demonstrate good agreement with the AA model, with the notable exception of the feature at very short distances (r ≈ 0.4 nm). While traditional molecular mechanics potentials fail to describe the intramolecular cross correlations (Fig. 5), distance-dependent pairwise interaction potentials are capable of reasonably describing the intermolecular cross correlations of hexane. This is consistent with the good performance of the MS-CG model in terms of accurately describing the RDFs. Still, our results demonstrate that the description of intermolecular interactions as a function of intramolecular state may assist in alleviating some of the standard problems experienced with structure-based coarse-graining (e.g., overly repulsive and over-fitted potentials) while providing a straightforward approach for characterizing the environment dependence of CG interactions.

C. Surface transferability
While our previous study, 56 as well as the results so far, highlighted the improved accuracy of conformational surface hopping over traditional CG structurebased schemes for a single system or thermodynamic state point, this section explores prospects of transferability. Our definition of transferability deviates somewhat from previous CG studies, as we do not aim to recycle the entire force field for a different state point. Instead, our methodology allows us to explore a weaker transferability regime: carrying over identical conformational surfaces, while reparametrizing their state probabilities. We focus on two aspects: temperature and chemical composition.

Temperature and compositional variation
We consider a set of 3 molecules that are chemically similar to hexane: octane, hexanediamine, and hexanediol. They correspond to the same alkane backbone with different terminal substitutions of a methyl hydrogen on each end: carbon, nitrogen, and oxygen, respectively (with appropriate saturation), as shown in Fig. 7d. Fig. 7 shows the variation of the state probabilities as a function of both chemical composition and temperature. The former is described via the electronegativity parameter χ of the substitution atom H, C, N, and O corresponding to hexane, octane, hexanediamine, and hexanediol, respectively. While we do not provide a formal justification for the use of χ, it offers a convenient proxy to describe the change in chemical composition through a continuous variable. Further, we have observed monotonic changes of our results with respect to χ. We will see below that χ offers a convenient parameter for scaling the non-bonded interactions.
The three panels of Fig. 7 show a two-dimensional projection of the state probabilities for each conformational surface: P S1 (χ, T ), P S2 (χ, T ), and P S3 (χ, T ). The most significant difference between the three surfaces is their range of state probabilities: larger for S 1 and S 3 , while smaller for S 2 . Surface S 1 varies significantly with respect to both parameters, S 2 is more sensitive to composition, and S 3 varies mostly against temperature. Their unique behavior sheds light on the conformational basin they represent: for instance, the population of S 2 is sensitive to the chemistry, but its low temperature dependence suggests an enthalpic stabilization. On the other hand, S 3 is rather insensitive to the chemistry, but significantly varies with temperature. While this could mean that S 3 is stabilized by entropy, we also note that as the fallback surface it amounts to a collection of different conformational basins. P S3 (χ, T ) varies remarkably little with respect to chemical composition, given its heterogeneous nature.

Temperature transferability
We first explore surface transferability across temperature. Starting from the three conformational surfaces obtained from reference AA simulations at T = 300 K (see Sec. IV B 1), we retain these surfaces and only tune the state probabilities to transfer to the other temperatures T = {250, 350, 400, 450} K. A comparison of the bond and bending-angle distributions generated by the AA and SH models are shown in Fig. 8 (a-d). The distributions show similar features as found in Fig. 4, monotonically evolving as a function of temperature. In particular, we find a strong temperature dependence of the long bond (b ≈ 0.26 nm) and the longest angle (θ ≈ 170 • ), while the other features show virtually no temperature dependence. Fig. 8 (e-f) also presents the CT-CT RDFs, which show reasonable agreement, although the SH distributions are somewhat too temperature dependent. In comparison to standard transferability properties, we note that the iter-gYBG model parametrized at T = 300 K and extrapolated to the other state points leads to similar performance for the one-dimensional distributions (Figs. S7 and S8). This is consistent with previous studies that have demonstrated that temperature-dependent, often linearly-scaled, interactions are necessary for accounting for the entropic contributions to the effective potentials. 27,28,[30][31][32] However, the SH model really shines when considering the description of cross correlations involving the bending angle (Figs. 5, S10 and S11), which standard parametrizations cannot reproduce.
Our weak transferability scheme offers an accurate reproduction of the distribution functions for all temperatures, despite the use of a single set of conformational surfaces. The results strongly suggest a large overlap in conformational space between the temperatures, adequately captured by retaining the conformational surfaces and simply adapting the state probabilities to each state point. We defer a deeper analysis of the temperature dependence of the state probabilities to Sec. IV C 1.

Compositional transferability
Beyond the transfer of force fields across temperatures, we now turn to the more challenging case of compositional transferability-across chemistry. We first assessed the transferability of surfaces in the gas phase, by employing the surface definitions obtained from octane to each of the other molecules. In this case, all molecules were simulated at T = 300 K. We note that hexane stands as an outlier in the set of compounds, due to its absence of heavy atoms beyond the six carbons. The impact of this difference will be illustrated below.
Panels a and b of Fig. 9 show a comparison of the bond distributions generated by the AA and SH models for the four molecules. All curves display overall similar behavior. Most strikingly, we observe a shift in the reference AA distributions: hexane shows its largest peak at larger g AA values of b, while the others are shifted to lower values, by up to 0.5Å. The reason for this shift is due to our choice of mapping: for consistency reasons we have kept the terminal CG bead defined as the center of mass of the two same carbons on the chain. Because of sterics, the presence of heavy atoms in octane, hexanediamine, and hexanediol have pushed these carbons slightly inward, resulting in the shifts observed in Fig. 9a. The interesting aspect here is the impact this has on our CG models: the use of a single set of surfaces will necessarily collapse all CG curves, only allowing for vertical shifts (by varying state probabilities). Interestingly we see little to no such vertical shifts, unlike what we had observed for temper- Hexane Octane Hexanediamine Hexanediol g AA ature variations (Fig. 8).
The bending-angle distributions show variations more in line with the above-mentioned temperature-variation study: all curves show similar behavior, i.e., peaks at the same places, with only variations in their heights. The results are different than the temperature variation with respect to the relative height differences between peaks: while varying T led to strong variations in the largest peak, it had virtually no effect on the second. In contrast, here we observe variations of similar magnitude between these two peaks. This strengthens the idea that a local change in chemical composition can be associated to a perturbation of the conformational space, akin to changing temperature. However, the local changes between peaks indicate that alterations occur at a more local level than an overall temperature rescaling. As a result, it would seem unlikely to reach compositional transferability of CG force fields by merely scaling it by a global prefactor. The CG distributions generated by the SH model reproduce the reference AA behavior remarkably accurately. They offer a useful compromise between a limited prefactor rescaling and state-point dependent potentials, and highlight the overlap in conformational space of the different molecules.
We also assessed chemical transferability in the liquid state. We first directly transfered the non-bonded force field for octane, while adjusting the state probabilities as described above. Each SH force field was probed at a distinct temperature T ref . T ref -corresponding to 300, 350, 435, and 470 K for hexane, octane, hexanediamine and hexanediol, respectively-was chosen to lie in approximately the same location with respect to the liquid phase existence for each molecule. This simple transfer of non-bonded interactions resulted in an underestimate of the changes in the CT-CT and CT-CM RDFs and an overestimate of the changes in the CM-CM RDF, while providing a good description of the intramolecular distributions (see Figs. S14 and S15). The discrepancies in the RDFs are not surprising, as we expect that the non-bonded interaction strengths associated with the CT bead should change as a function of chemistry. To test the impact of such changes, we performed a simple scaling of the octane non-bonded interactions. In particular, we applied a scaling factor to each of the CT-CT potentials equal to the ratio of electronegativity values of the corresponding substituted atoms: U M;CT-CT = (χ M /χ octane ) U octane;CT-CT , where M = hexane, hexanediamine, or hexanediol. Similarly, the CT-CM potentials were scaled by the square root of the same ratio (assuming an effective geometric mean combination rule): U M;CT-CM = ( χ M /χ octane ) U octane;CT-CM . The original octane CM-CM interactions were employed without adjustment. The full set of scaled potentials are presented in Figs. S12 and S13.
Remarkably, as shown in panels e and f of Fig. 9, this heuristic scaling of potentials along with the adjusted state probabilities results in an accurate description of the local packing as a function of changes in chemistry, despite employing a single set of surfaces for the molecules. The accuracy of the CT-CM RDFs are also improved (relative to the non-scaled SH model) while retaining a good description of the intramolecular distributions, although the discrepancy in the CM-CM RDF is somewhat exacerbated (Fig. S15). We note that the absolute accuracy of all CG CM-CM RDFS (i.e., also for the MS-CG and iter-gYBG models) is slightly degraded due to the challenging representation applied to the non-hexane molecules, as demonstrated in detail in the Supporting Information.

V. CONCLUSIONS
This work extends our previous presentation of the coarse-grained (CG) conformational surface-hopping (SH) methodology: analogous to switching between different electronic states, we define one force field per conformational basin and hop between them. 56 Each force field is parametrized by applying force matching (i.e., the MS-CG method) while using only configurations from the corresponding basin. Our illustration of the method using a toy example highlights the benefits of enforcing a set of target state probabilities, which avoids possible errors due to an inaccurate description of the global surface at the barrier between two conformational states. While the SH models employ standard molecular mechanics interaction functions in the Hamiltonian, the focus on reproducing local properties of each surface results in increased accuracy relative to standard models. The results are particularly striking for the gas-phase properties of a three-bead hexane representation: the correlations between bond and bending angle, notoriously problematic for the MS-CG method, are accurately represented with the SH approach while employing only three surfaces.
We have also presented an extension of the SH method to intermolecular interactions: conformational surfaces are defined based on the intramolecular state of the molecule, while the intermolecular interactions depend on the pair of surfaces involved. For instance, the 3surface model for hexane consists of two distinct bead types, which corresponds to a total of 18 unique interactions (i.e., 6 interactions for each pair of bead types). The resulting SH models retain an accurate description of the local packing, while also demonstrating slight improvements in the RDFs compared with the MS-CG model. Perhaps more interestingly, an assessment of the SH potentials demonstrated promising properties with respect to the other structure-based potentials. In particular, the SH potentials tended to be more attractive with a single local minimum, counteracting two common problems with structure-based models: (i) the underestimation of the cohesive energy in liquids and (ii) an over-fitting of the features at the state point of parametrization.
We further investigated the capabilities of the SH models by examining their transferability properties. We focused on a so-called "weak-transferability regime," in which one state point determines the surface definitions; these surfaces are then transferred to other state points, while adjusting their state probabilities. In particular, we considered the transfer of state definitions across both temperature and chemical composition. In the latter case, where the strength of the interactions are expected to change as a function of chemistry, the use of the electronegativity parameter, χ, provided a useful proxy to scale the non-bonded interactions. Our results demonstrate that the SH models not only accurately describe the trends in the intramolecular distributions, which are largely reproduced with traditional models, but also better represent intramolecular cross correlations through-out the liquid state. The SH approach demonstrates similar results with respect to the description of local packing as a function of temperature for the molecules considered, but slightly overestimates the temperature dependence of the RDFs. It would be interesting in this context to explore the entropic contributions to the SH potentials. [30][31][32] The investigation of chemical transferability focuses on terminal substitutions via the comparison between hexane, octane, hexanediamine, and hexanediol. Notably, we find limitations in modeling the bond distributions, as the substitution of hydrogens to heavy atoms (i.e., moving from hexane to one of the other three molecules) shifts the distribution. Aside from this limitation, the tuning of individual state probabilities appears to be a promising framework for considering the construction of CG models that are not restricted to one state point, but rather applicable to a neighborhood of thermodynamic parameters and chemical compositions. We foresee the weak-transferability regime brought forward here to be of use when parametrizing not just one reference simulation, but collections of state points or compounds. This will be of use in the context of parametrizing CG models across subsets of chemical space.
Finally, we stress the conceptual and practical advantage of parametrizing the SH models using the MS-CG technique. The combined approach offers an enhanced capability to describe complex cross correlations between degrees of freedom that arise additively in the Hamiltonian, while using a direct inverse parametrization scheme. Since the MS-CG method results in errors whenever the AA cross correlations represent an inappropriate proxy for the cross correlations of the resulting CG model, the approach provides an automatic validation of the surface definitions. In other words, if there remain cross correlations within a single surface that cannot be reproduced by a molecular mechanics force field, errors will likely appear in the description of the modes along each distribution function corresponding to the inadequate surface. Moreover, the potentially large number of force fieldsup to one per conformational basin-can be derived independently, an aspect that would not be straightforward using iterative methods.