Bayesian Optimization for High-Dimensional Coarse-Grained Model Parameterization: A Case Study on Pebax Polymer

TL;DR

This paper shows that Bayesian optimization can effectively tune complex polymer models, enabling accurate simulations of material properties.

Contribution

The study demonstrates Bayesian optimization's viability for high-dimensional polymer model parameterization.

Findings

Bayesian optimization with TPE successfully parametrized a 41-parameter Pebax model.

The optimized model accurately reproduces structural and thermodynamic properties of the atomistic system.

BO-TPE outperforms traditional methods in convergence speed and accuracy.

Abstract

Coarse-grained (CG) force field models are extensively
utilized
in material simulations because of their scalability. Ordinarily,
these models are parametrized using hybrid strategies that sequentially
integrate top-down and bottom-up approaches. However, this combination
restricts the capacity to jointly optimize all parameters. Although
Bayesian optimization (BO) has been explored as an alternative search
strategy to identify well-optimized CG parameters, its application
has conventionally been limited to low-dimensional scenarios. This
has contributed to the assumption that BO is unsuitable for more complex
CG models, which often involve a large number of parameters. In this
study, we challenge this assumption by successfully extending BO,
using the tree-structured Parzen estimator (TPE) model, to optimize
a high-dimensional CG model. Specifically, we show that a 41-parameter
CG model of Pebax-1657, a copolymer composed of alternating polyamide
and polyether segments, can be effectively parametrized using BO,
resulting in a model that accurately reproduces the key physical properties
of its parent atomistic representation. Our optimization framework
simultaneously targets structural and thermodynamic properties, namely,
density, radius of gyration, and glass transition temperature. Compared
to traditional search algorithms, BO-TPE not only converges faster
but also delivers consistent improvements over more standard parametrization
approaches.

Figures

Chemical structure of Pebax and its polyamide (PA) and polyether (PE) segments.

Panel (a), chemical structure of Pebax-1657 composed of PA6 and PEO segments. Panel (b), AA model of the Pebax-1657 chain with a composition of 40% PA6 and 60% PEO, corresponding to n = 1, x = 2, and y = 3. The gray spheres represent carbon atoms, while the red, blue, and white spheres represent oxygen, nitrogen, and hydrogen atoms, respectively.

CG model of Pebax-1657 chain considering five distinct beads. The numbers 1 to 10 clearly visualize the structural organization and the total number of beads. The gray spheres represent carbon atoms, while the red, blue, and white spheres represent oxygen, nitrogen, and hydrogen atoms, respectively.

Hybrid strategy combining top-down and bottom-up methodologies to develop CG models.

(Panel a) Convergence of BO as a function of the iterations. The solid blue line and shaded blue region indicate the logarithm of L(θ) (eq ) and the standard deviation computed with three independent optimization trajectories. The insight is the histogram of the sampled values of log⁡⁡L(θ) of one of the BO trajectories. (Panel b) t-SNE visualization of the sampled parameters of a single BO run. Each point represents a sampled parameter set, color-coded by its total objective function value. ▲ and ★ represent the parameters of the initial and best-found CG model. We used a logarithmic scale to improve the visualization of the ranges for each property’s error function.

Panel (a) Cumulative explained variance as a function of the number of principal components for the CG parameters sampled by BO. The variance curves indicate that most of the variability in the optimized parameter space can be captured by a reduced number of effective dimensions; 20 to 25 components account for 90% of the variance. Panel (b) Distribution of sampled parameters obtained by BO for each component of the total objective function (eq ). The ▲ and ★ symbols denote the initial parameter set and the optimal configuration found by BO, respectively.

Comparison of the density (Panel a) and the radius of gyration (Panel b) for a single Pebax chain using the atomistic model as a reference, alongside the CG-hybrid strategy and CG-BO models.

Panels (a) and (b) show the predicted density and radius of gyration, respectively, for the three optimal CG models, CG-BO­(w i ), obtained using different weight vectors w i in the total loss function. Panels (c) and (d) display the sampled values of log⁡⁡Lρ and log⁡⁡LRg corresponding to each optimization. ★ symbols mark the final CG-BO model obtained for each choice of w. Each optimization ran for 600 iterations to ensure convergence. See the main text for further details.

Tables

1. Predicted Glass Transition Temperature and Its Relative Error (Ltg(θ)) with the CG-BO and CG-Hybrid Strategy Models

Model	Tg (K)	LTg(θ)
Atomistic	378.10	-
CG-BO	336.69	0.012
CG-Hybrid Strategy	247.18	0.119

Full text

Introduction

1

Molecular dynamics (MD) is a well-known computational technique widely used to study and simulate physical and chemical phenomena. ?,? By employing molecular models such as all-atom or united-atom representations, MD has allowed prediction of the behavior of complex systems under a wide range of thermodynamic conditions. ?−? ? All-atom models offer the highest level of atomic detail by individually representing each atom, making them ideal for capturing phenomena governed by precise atomic-scale interactions. ?,? In contrast, united-atom models consolidate specific groups of atoms into a single interaction site, effectively simplifying the system and reducing computational demands. In both approaches, accuracy relies on the choice of force fields, which describe interatomic interactions and enable the calculation of structural, thermodynamic, mechanical, or transport properties. However, some phenomena extend beyond the typical size and time treated at the molecular scale, requiring a mesoscale perspective to capture their dynamics effectively. This necessity underscores the importance of coarse-grained (CG) models, as emphasized in several recent reviews. ?−? ?

CG models simplify molecular representations by grouping atoms or functional groups into beads. This simplification inevitably results in a loss of chemical resolution, but when carefully designed, CG models can reproduce key properties of the system.? Two main approaches guide the development of CG models: (i) bottom-up and (ii) top-down. Bottom-up frameworks derive CG parameters by mapping detailed quantum-level or atomistic information to a coarser scale, often using techniques such as iterative Boltzmann inversion,? force matching? or relative entropy.? This family of methodologies emphasizes accuracy in reproducing molecular-level properties but lacks transferability across diverse thermodynamic conditions.? In contrast, top-down approaches leverage experimental macroscopic data to optimize CG parameters directly, aiming for broader applicability. Advanced top-down methods include the MARTINI force field,? widely used for biochemical applications, and the Statistical Associating Fluid Theory (SAFT-γ Mie) CG force field, which has found applications in fluid phase equilibria. ?,?

Hybrid strategies combine both top-down and bottom-up methodologies in a sequential optimization process to develop CG models that are computationally efficient and capable of accurately representing molecular interactions. For example, Rahman et al.? integrated the SAFT-γ Mie with atomistic simulations to model linear alkanes. The SAFT-γ Mie EoS represents molecules as chains of tangentially bonded spherical segments interacting via the Mie potentiala generalization of the Lennard-Jones potential that allows for independent tuning of repulsive and attractive forces. In the top-down optimization process, SAFT-γ Mie is employed to derive nonbonded interaction parameters by fitting experimental thermophysical data, such as densities, vapor pressures, and phase equilibria, ensuring that CG models align with macroscopic observations. The bottom-up stage complements this by refining bonded interactions through atomistic simulations, maintaining consistency with structural features at the molecular level.

As emphasized in previous works,? the importance of accurately capturing both intra- and intermolecular interactions is critical for obtaining reliable predictions of thermodynamic properties such as vapor–liquid equilibrium and densities. Similarly, Fayaz-Torshizi and Müller? applied these methodologies to polymers, demonstrating that incorporating detailed interactions significantly enhances the model’s ability to predict structural and dynamic behavior. By mapping atomistic configurations to CG beads, the models ensure that the bond lengths and angles closely reflect the equilibrium states observed experimentally. The harmonic potentials governing bond stretching and angular bending were calibrated using distributions of bond lengths and angles derived from atomistic simulations. These distributions were typically fitted with weighted Gaussian functions, ensuring consistency with key structural features such as end-to-end distances and radii of gyration. Although effective, this sequential optimization process may occur in isolated stages, providing no assurance of achieving a global solution, i.e., the optimal set of CG parameters. The decoupled nature of the optimization may lead to suboptimal CG models. It should be noted that, although decoupling the optimization problem in the development of CG models can lead to suboptimal parametrizations, there are successful examples of general-purpose CG force fields that employ this approach. A notable example is the SPICA family of force fields, based on the Shinoda-DeVane-Klein (SDK) multiproperty fitting framework, ?−? ? which provides one prominent example of such a hybrid strategy, demonstrating good transferability for surfactants, lipids, and biomolecular systems.

Machine learning (ML)-based CG models have also emerged as an alternative parametrization of CG models. ?−? ? However, despite significant progress, challenges remain, such as the need for extensive data sets and ensuring model transferability across diverse systems.? Other alternative schemes based on ML models use Bayesian optimization (BO), a search algorithm, to identify the well-optimized parameters for CG models. BO is a powerful framework for design optimization, with successful applications spanning diverse fields such as robotics, environmental monitoring, and experimental design.? In materials science, BO has proven highly effective across a range of tasks. For example, in the context of physical models, it has been employed to screen chemical compounds, ?−? ? minimize the energy of the Ising model,? and optimize laser pulses for molecular control.? In computational chemistry, BO has been applied to calibrate density functional models,? refine physical models for the cis-trans photoisomerization of retinal in rhodopsin, ?,? and design potential energy surfaces for reactive molecular systems.?

BO provides a systematic and efficient approach for exploring high-dimensional parameter spaces, ?−? ? ? ? ? ? ? making it particularly well-suited for CG models of molecular and nanostructured materials, where traditional methods often struggle. Its application to CG models with low-dimensional parameter space has been explored in refs ?−? ? , where the optimization process involved extracting physical properties from MD simulations and incorporating them into an objective function, typically defined as the sum of relative errors compared to an atomistic model or, by following a multiobjective optimization framework. The resulting CG models demonstrated a good agreement with their atomistic counterparts, highlighting BO’s potential for improving CG model development. However, these studies were restricted to CG models with relatively few parameters ?,?,? and did not fully account for all inter- and intramolecular interactions during optimization.

Due to recent studies demonstrating the successful application of BO in high-dimensional optimization problems involving up to several thousand parameters, ?−? ? ? ? ? ? ? we extend its use to the development of a system-specific CG model for Pebax. Unlike previous works limited to low-dimensional parameter spaces, our approach leverages BO’s scalability to address the complex, high-dimensional optimization landscape associated with this copolymer or similar systems. Here, we demonstrate that BO can efficiently optimize CG models with a substantially larger number of parameters, capturing both intra- and intermolecular interactions. In this system, the parameter search space exceeds 40 dimensions, underscoring BO’s scalability and its ability to navigate and refine intricate molecular interaction landscapes, ultimately enhancing the accuracy and versatility of CG models.

The remainder of this paper is organized as follows: Section introduces the CG modeling framework for copolymers. Section outlines the BO algorithm applied for CG model optimization, with emphasis on the physical properties used for calibration. Finally, Section presents and discusses the results of our optimization strategy.

Coarse-Grained Models for Copolymers

2

Pebax Physical Model

2.1

Pebax copolymers,? composed of alternating polyamide and polyether segments, are widely used in membrane technologies due to their tunable mechanical and transport properties. Pebax consists of polyamide (PA) and polyether (PE) segments, as shown in Figure. The PA segment, typically polyamide 6 (PA6) or polyamide 12 (PA12), provides mechanical strength and is covalently linked to the PE segment via ester groups.? The PE component, which can be poly­(ethylene oxide) (PEO) or poly­(tetramethylene oxide) (PTMO),? enhances gas transport and separation. As a case study, we focus on Pebax-1657, a specific grade composed of PEO and PA segments, as illustrated in Figurea. The repeat unit of the Pebax-1657 chain consists of 40% PA6 and 60% PEO, ?−? ? with each chain composed of one repeat unit only, as illustrated in Figureb. The PA-to-PEO ratio plays a crucial role in balancing mechanical flexibility and gas permeability, making Pebax-1657 particularly suitable for membrane applications such as CO_2_ capture and gas separation. ?,? Additionally, the incorporation of nanostructured materials, including MoS_2_ nanosheets,? nonionic surfactants,? zeolitic imidazolate frameworks-8 (ZIF-8) particles,? and Azo@MOF-199?–has been shown to enhance its adsorption and transport properties further. ?,? Given its technological potential and the challenges associated with fully atomistic simulations, Pebax-1657 is the polymer of choice in this work. An accurate coarse-grained (CG) model is essential for studying its structural and transport properties at larger scales.

Chemical structure of Pebax and its polyamide (PA) and polyether (PE) segments.

Panel (a), chemical structure of Pebax-1657 composed of PA6 and PEO segments. Panel (b), AA model of the Pebax-1657 chain with a composition of 40% PA6 and 60% PEO, corresponding to n = 1, x = 2, and y = 3. The gray spheres represent carbon atoms, while the red, blue, and white spheres represent oxygen, nitrogen, and hydrogen atoms, respectively.

The all-atom (AA) representation of the Pebax-1657 polymer was constructed using the Material Exploration and Design Analysis (MedeA) simulation software (version 3.5).? The system consists of 50 polymer chains arranged in an amorphous configuration, generated using the Amorphous Builder module. The dimensions of the initial simulation box were set to 32.6 Å × 32.6 Å × 32.6 Å (length × width × height), yielding an experimental density of 1.136 g/cm^3^.?

Pebax CG Model

2.2

Given the technological relevance of Pebax-1657, an accurate coarse-grained (CG) model is essential for studying its structural and transport properties at larger scales. The CG model, shown in Figure, was developed based on available data from small oligomers, assuming the transferability of parameters. The atom-to-bead mapping followed a protocol to determine the intermolecular parameters of the CG model using the SAFT-γ Mie group contribution method. The model classifies the beads into five types: T1, T2, T3, T4, and T5, reflecting their chemical structure and functional roles within the polymer. The T1, T2, and T3 beads represent the polyamide segment, while T5 corresponds to the polyether segment. The T4 bead serves as a linker between the PA and the PEO regions.

CG model of Pebax-1657 chain considering five distinct beads. The numbers 1 to 10 clearly visualize the structural organization and the total number of beads. The gray spheres represent carbon atoms, while the red, blue, and white spheres represent oxygen, nitrogen, and hydrogen atoms, respectively.

For this standard CG model, a hybrid strategy combining top-down and bottom-up methodologies was employed to derive both intermolecular and intramolecular FF parameters. Figure summarizes the usual procedure to obtain the nonbonded intermolecular parameters [ϵ, σ, λ] using the SAFT-γ Mie group contribution method, where ϵ, σ, and λ represent the depth of the potential well, bead diameter, and Mie potential parameter that controls repulsion and attraction contributions, respectively. All intramolecular interactions were modeled using harmonic potentials for bond stretching and angle bending,? and these parameters [K _ b _, K ϕ, ϕ] were derived through fully atomistic simulations using MD simulations based on the PCFF+ potential.?

Hybrid strategy combining top-down and bottom-up methodologies to develop CG models.

In the SAFT formalism, the SAFT-γ Mie framework determines nonbonded intermolecular parameters by minimizing an objective function $[eqn]$ that compares experimental data with predictions from the equation of state (EoS) for target properties such as saturated liquid densities and vapor pressures.? This approach frames the determination of top-down interaction parameters as an optimization over macroscopic bulk properties.? A Python-based implementation provides a practical framework for such optimizations.?

For simplicity, we employed Bottled SAFT.? Bottled SAFT is a web application that provides SAFT-γ Mie parameters for a wide range of molecular fluids, allowing researchers to search for a specific molecule or upload a structure file to generate a set of parameters. Bottled SAFT is freely available at http://www.bottledsaft.org and is widely recognized for its accessibility and utility. For beads T1, T3, and T4, the parameters were obtained, considering heptanoic acid (CAS 111-14-8), N-methylacetamide (CAS 79-16-3), and pentyl pentanoate (CAS 2173-56-0), respectively. While Bottled-SAFT parameters derived from n segments are typically used to construct CG models for homonuclear molecules, several studies have adopted an alternative approach to build CG models. Notably, in the CG modeling of polystyrene,? a single representative segment for toluene, from the n segments, was selected to construct the CG model. In line with this precedent, we believe that selecting a representative segment from each relevant molecule, along with its corresponding Bottled-SAFT parameters, offers a more accurate representation for our CG model. Meanwhile, data for bead T2 and the PEO segment (T5) were derived from previous studies of similar chain groups.? For intramolecular force fields, all the steps illustrated in Figure were performed to derive the bond stretching parameters (K _ b _) and angle bending parameters ([K ϕ, ϕ]).

Bayesian Optimization

2.3

The development of CG models commonly depends on hybrid strategies that combine top-down and bottom-up methodologies, as shown in Figure. In these hybrid strategy approaches, the optimization is divided into two subproblems: (i) the intermolecular and (ii) the intramolecular parameters. First, the well-optimized parameters of the intermolecular potential [ϵ, σ, λ] are determined by minimizing the error between the reference data and the EoS calculations. Once the intermolecular parameters are optimized, a bottom-up approach follows to tune the intramolecular potential parameters, [K _ b _, K _ ϕ _, ϕ], for harmonic force fields.

Although hybrid approaches are successful, ?,?,?,?,? the division of the global optimization problem into separate subproblems can lead to inefficient exploration of the search space. This decomposition often results in nonoptimal solutions, as it fails to capture the interdependencies between parameter groups. To address this limitation, an alternative optimization approach is needed: one that explicitly accounts for these interdependencies while remaining computationally feasible, even for high-dimensional parameter spaces and complex physical systems. BO is often regarded as unsuitable for high-dimensional problems, primarily due to the challenge of constructing accurate surrogate models in large parameter spaces with limited data. This sense is reflected in prior applications to CG model optimization, where the number of parameters is typically limited to ∼10 or fewer. ?,?,? Despite these dimensional limitations, BO has been successfully applied to developing CG models by optimizing force fields based on physical properties extracted from molecular dynamics simulations. ?,? While existing applications of BO to CG models have been constrained to relatively low-dimensional problems, recent advances have demonstrated that BO can be effective in high-dimensional settings. ?−? ? ? ? ? ? ? For example, BO was successfully applied to a set of high-dimensional optimization problems involving up to 1,000 variables,? demonstrating its scalability and adaptability to large, complex search spaces. These findings suggest that BO could be an effective tool for optimizing high-dimensional CG models and overcoming prior limitations of existing optimization algorithms, which do not fully account for all inter- and intramolecular interactions, thereby limiting the robustness of the resulting models.

For the Pebax polymer, we speculate that the specification of three key physical properties, namely (i) density (ρ), (ii) radius of gyration (Rg), and (iii) glass transition temperature (Tg) is sufficient to bootstrap the CG representation. The proposed optimization framework integrates these three properties into a single objective function,

$$[eqn]$$

where w = [w ρ, w _ Rg _, w _ Tg _] are weight coefficients that balance the relative importance of each property. These weights prevent any single property from disproportionately influencing the optimization, which is particularly important here since density has more fitting points than the other properties. We used w ρ = 1, w _ Rg _ = 1,875, and w _ Tg _ = 1.5 × 10^4^, see Supporting Information (SI) for more discussion regarding these values. The individual loss terms in eq are computed as relative errors between the CG model predictions and the reference data,

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where T _ i _ denotes the temperatures at which MD simulations were performed, while N ρ and N _ Rg _ represent the number of data points used for density and radius of gyration calculations, respectively. These properties play distinct roles in the validation of CG models. The density provides direct experimental comparability and is directly obtained in MD. The radius of gyration reflects the spatial conformation and compactness of polymer chains, and the glass transition temperature integrates thermodynamic and structural information, offering insight into thermal behavior.

For the copolymer Pebax, $[eqn]$ is influenced by the glass temperature. To ensure representative sampling, we selected eight discrete temperatures $[eqn]$ , four below and four above the glass transition temperature observed in the atomistic model; see Figure S1 in the SI. The glass transition temperature was then determined as the intersection of these two fitted curves, following the methodology outlined by Patrone et al.? Further numerical details are provided in the SI. Temperatures near the glass transition were excluded, as the density–temperature relationship is discontinuous at the glass temperature.? We used density measurements at each of the eight temperatures to fit two linear curves. This approach provided a denser set of target points on either side of the glass transition, refining the density–temperature relationship for optimization. Initially, we included points near the glass transition in the fits. However, we found that excluding them led to improved results. The glass transition temperature was then determined as the intersection of these two fitted curves, following the methodology outlined by Patrone et al.?

For $[eqn]$ , the radius of gyration component, we also ensured reliable sampling by simulating a single Pebax chain instead of a fully populated polymer box. This approach mitigated the risk of the sampling of the radius of gyration becoming confined to a limited region of phase space, which could otherwise lead to inaccurate property estimates. By isolating a single chain, we enhanced sampling across phase space. The radius of gyration was then computed at 15 different temperature points. Additional details on the MD simulations are provided in the SI. All target properties used as reference data in eq were determined using an atomistic model based on the PCFF+ force field. The data supporting the findings of this study are openly available in the following https://github.com/camjjr/bo_cgff. The parameters of the atomistic model used in this work can be found in the following link.

Tree-Structured Parzen Estimator

2.4

Grounded in Bayesian inference, BO systematically balances exploration and exploitation by iteratively updating a surrogate model that navigates the search space efficiently. By leveraging information from previous evaluations, BO aims to identify the most favorable regions and direct the search toward the global optimum. Two key components constitute the core of BO: the surrogate model, which approximates the objective function, and the acquisition function, which guides the search. We used as a surrogate model the Tree-structured Parzen Estimator (TPE).? Unlike traditional Gaussian process-based approaches, TPE constructs a nonparametric density estimator that partitions the search space into two regions: (i) a top quantile of observations with favorable objective function values and (ii) a lower quantile containing the remaining observations.? This structure allows TPE to model the search space adaptively, prioritizing promising regions while maintaining diversity in exploration. Formally, TPE models? parametrize the conditional probability distribution as,

$$[eqn]$$

where p(θ|y, D) represents the probability density function of the parameters, and D is the set of observed objective function values {y 1, y 2, y 3, ···, y _ n _} obtained from evaluations of eq. The subsets D ^ b ^ and D ^ w ^ correspond to the best (b) observations and worst (w), respectively. These groups are dynamically updated at each iteration based on the hyperparameter γ, which controls the balance between exploration and exploitation.

The probability density function within each group is estimated as follows,

$$[eqn]$$

where ω^ i ^ represents the weights assigned to observations, i denotes either the best (b) or worst (w) group, and N _ i _ is the number of observations within the group. The term p 0 defines a prior distribution that influences the level of exploration, while k _ n _(θ, θ _ n _|b ^ i ^) is the kernel function used for density estimation, employing a Gaussian kernel for numerical variables and an Aitchison-Aitken kernel for categorical variables. b ^ i ^ is simply the bandwidth for each group. For this study, we used the following acquisition function,

$$[eqn]$$

Using eq, the acquisition function is,?

$$[eqn]$$

which is equivalent to the well-established PI acquisition function. ?,? The TPE model and the BO algorithm were both implemented using the Optuna library.?

Results

3

CG Optimization

3.1

Optimization of CG models, even with BO, is a dynamic process, as at each iteration BO infers more about $[eqn]$ . Figurea illustrates the lowest value of $[eqn]$ found by BO as a function of the iterations. To demonstrate BO’s robustness, we considered three independent runs, each with an initial random θ. All runs successfully converged to similar optimal solutions with similar $[eqn]$ , in less than 600 iterations (≈2000 min), as shown in Figurea. Furthermore, the optimal solutions identified in each run exhibit objective function values several orders of magnitude lower than the initial guesses, and by the hybrid strategy. These results suggest that BO effectively avoids regions where $[eqn]$ is high (red points in Figureb) while concentrating trials in the low areas (blue points in Figureb). This is further supported by the histogram of the sampled values of $[eqn]$ of a single BO run; insight panel in Figurea.

(Panel a) Convergence of BO as a function of the iterations. The solid blue line and shaded blue region indicate the logarithm of L(θ) (eq ) and the standard deviation computed with three independent optimization trajectories. The insight is the histogram of the sampled values of log⁡⁡L(θ) of one of the BO trajectories. (Panel b) t-SNE visualization of the sampled parameters of a single BO run. Each point represents a sampled parameter set, color-coded by its total objective function value. ▲ and ★ represent the parameters of the initial and best-found CG model. We used a logarithmic scale to improve the visualization of the ranges for each property’s error function.

It is important to emphasize the capability of this CG model, which consists of 41 parameters. Conventional search methods, such as grid-based approaches, typically require more than 10 trials or grid points per parameter of the CG model. Given the complexity of the Pebax CG model, a grid search approach would need more than 10^40^ iterations, an infeasible number considering the computational cost of MD simulations to evaluate the three target physical properties. Recent studies have also demonstrated the efficiency of BO in optimizing CG models. ?,? These previous studies focused on simpler CG models, resulting in lower-dimensional optimization problems. To our knowledge, this study is the first to optimize all parameters of a complex CG model, resulting in a high-dimensional optimization problem, due to the number of parameters in the CG model. The results demonstrate that BO-TPE remains effective even for optimization problems with more than a dozen parameters, reinforcing its suitability for complex CG model development.

To understand why BO is particularly effective for the Pebax CG model, we investigated whether its 41 parameters could be represented in a lower-dimensional space. We first applied Principal Component Analysis (PCA) to the parameters sampled by BO (Figure), finding that 90% of the variance is captured by the first 28 principal components. While this only reduces the dimensionality by about 32%, it does indicate some redundancy. To probe for deeper structure, we also performed an ISOMAP analysis, but it did not uncover clear patterns linking low-loss regions, suggesting the absence of a simple nonlinear manifold. Repeating the PCA on samples with $[eqn]$ slightly improved compression: 23 components sufficed to explain 90% of the variance (Figurea). This modest gain suggests that low-loss regions exhibit somewhat more structure. Recent work supports the idea that BO can benefit from linear embeddings in high-dimensional spaces.? Our results align partially with this viewparameter correlations do emerge in regions of low objective, but retaining 65% of the original dimensions remains necessary. Overall, these findings suggest that BO’s success here stems less from a highly compressible geometry and more from its ability to adaptively explore a moderately structured, high-dimensional space with some redundancy but no sharply defined low-dimensional manifold.

Panel (a) Cumulative explained variance as a function of the number of principal components for the CG parameters sampled by BO. The variance curves indicate that most of the variability in the optimized parameter space can be captured by a reduced number of effective dimensions; 20 to 25 components account for 90% of the variance. Panel (b) Distribution of sampled parameters obtained by BO for each component of the total objective function (eq ). The ▲ and ★ symbols denote the initial parameter set and the optimal configuration found by BO, respectively.

While the objective function in eq aggregates multiple target properties into a single global function, the CG model’s optimization can alternatively be approached as a multiobjective optimization problem. In this case, the goal is to identify the Pareto front, which characterizes the trade-offs between competing objectives. ?−? ? Although this approach provides a richer understanding of the parameter space, it is less practical for high-dimensional CG models, as this optimization framework requires a greater number of evaluations. Given the linear combination nature of the total objective function, we also analyze the sampled values of each component of $[eqn]$ . Figureb presents the individual values of $[eqn]$ , $[eqn]$ , and $[eqn]$ sampled by BO. These results show that BO effectively samples low values $[eqn]$ for all three components, without a clear preference for any one. Furthermore, these results also indicate that density and glass transition temperature exhibit lower overall relative errors compared to the radius of gyration, suggesting that Rg is the most challenging physical property to reproduce. The relatively low errors observed for density and glass transition temperature in some of the different sampled parameters of the CG model (each point in Figureb) may be attributed to their intrinsic relationship, as Tg can be determined from the slopes of the density curve; for more details, see refs ?,? . Overall, our findings suggest that none of the physical properties overwhelmingly dominate the objective function and indicate that BO successfully balances its optimization. This balance enhances the robustness of the resulting coarse-grained models, ensuring a more reliable reproduction of the target physical properties.

Physical Properties

3.2

Figure compares the density and radius of gyration for a single Pebax chain, as predicted by three different models: the atomistic reference model, the CG-BO model, and the CG-Hybrid Strategy model, the latter was developed using a hybrid optimization approach described in the SI. The well-optimized parameter set θ obtained by BO is reported in the SI, referred to as the CG-BO model. To mitigate the effects of discontinuity in the density profile, only four temperature below the glass temperature (150 K < T < 225 K) and four high-temperature above it (425 K < T < 500 K) were considered in the optimization, as detailed in SI.

Comparison of the density (Panel a) and the radius of gyration (Panel b) for a single Pebax chain using the atomistic model as a reference, alongside the CG-hybrid strategy and CG-BO models.

Regarding the density (Figurea), the CG-BO model closely follows the atomistic reference, with only minor deviations at low temperatures (T < 250 K). In contrast, the CG-Hybrid Strategy model exhibits significant discrepancies across the entire temperature range, suggesting poor transferability of its parameters, which were originally derived from different polymer systems. It is worth mentioning that the density of the atomistic model deviates from the experimental measurements by 4.3%.? The predicted radius of gyration (Figureb) further highlights the limitations of the hybrid approach. The CG-Hybrid Strategy model significantly deviates from the atomistic reference, particularly at temperatures above 200 K. Conversely, the CG-BO model shows strong agreement with the reference data, effectively capturing the general trend of Rg. This demonstrates BO’s capability to systematically explore the parameter space and identify a coarse-grained model that retains key structural characteristics, even within a high-dimensional optimization setting.

For the glass transition temperature (Tg), Table summarizes the predicted values obtained from the three models. Using the atomistic model as a reference, once again, the CG-BO approach outperforms the CG-Hybrid strategy, producing a significantly lower relative error (10.95%), while the CG-Hybrid strategy model exhibits a relative error of 34.62%.

1: Predicted Glass Transition Temperature and Its Relative Error (Ltg(θ)) with the CG-BO and CG-Hybrid Strategy Models

Lastly, we investigated the sensitivity of the objective function (eq) to changes in the weight parameters w = [w ρ, w _ Rg _, w _ Tg _], given our choice of a single-objective optimization approach over a full multiobjective framework. Although multiobjective optimization, as proposed by Sestito et al.? is a viable alternative for CG model tuning; it would require significantly more evaluations of $[eqn]$ and MD simulations, which becomes impractical for high-dimensional CG models. Furthermore, such an approach would yield a set of Pareto-optimal solutions, making model selection less straightforward if no clearly dominant Pareto point exists.

Figure shows the objective space corresponding to two different sets of weight combinations, w 2 = [10, 1875, 1.5 × 10^3^] and w 3 = [5, 18750, 1.5 × 10^3^], where we deliberately shifted the balance between the three physical properties. Each resulting CG model is denoted CG-BO­(w _ i _). As seen in Figurea and b, all three weight configurations lead to similarly accurate predictions for ρ and Rg. Figurec and d show that each w _ i _ emphasizes a different region of the individual loss components, yet the joint performance remains consistent. The predicted glass transition temperature in all three optimization cases deviates by 82.94 K for w 2 and 37.01 for w 3 from the atomistic model. Overall, the CG-BO­(w 3) reproduces more accurately ρ and Tg, but predicted Rg is ∼30% less accurate than CG-BO­(w). These results indicate that the optimization landscape is relatively robust to moderate changes in w and that the BO-CG model is flexible enough to capture relevant features regardless of the precise weighting scheme. For the W 2 and W 3 models, BO-TPE also converged to the well-optimized parameters in fewer than 600 total iterations; see Figure S2 in the SI.

*Panels (a) and (b) show the predicted density and radius of gyration, respectively, for the three optimal CG models, CG-BO­(w

i ), obtained using different weight vectors w

i in the total loss function. Panels (c) and (d) display the sampled values of log⁡⁡Lρ and log⁡⁡LRg corresponding to each optimization. ★ symbols mark the final CG-BO model obtained for each choice of w. Each optimization ran for 600 iterations to ensure convergence. See the main text for further details.*

Despite the demonstrated scalability and efficiency of the proposed BO-TPE framework, two main limitations remain. First, as with most Bayesian optimization methods, the accuracy and generality of the approach depend on the design of the objective function itself, as it defines how model performance is quantified across the parameter space. Second, the intrinsic cost of evaluating the objective function, which in our case involves molecular simulations. The current implementation was tailored to polymer systems through their glass-transition behavior and would require adaptation to other material classes or force fields. Future work will focus on developing transferable objective functions that extend the applicability of this framework to broader families of CG models.

Conclusion

4

We demonstrated that Bayesian optimization enables the efficient optimization of coarse-grained models with a large number of parameters. Unlike conventional optimization methods that often rely on parameter-splitting heuristics to navigate complex search spaces, BO adopts a fundamentally different strategy, using acquisition functions to guide efficient sampling and globally explore the parameter space. Using BO, we optimized a CG model for the copolymer Pebax without decomposing the problem into smaller subtasks, as is typically done in CG model development. Remarkably, BO-TPE required fewer than 600 iterations to converge to an optimal set of 41 parameters. The resulting CG model accurately reproduces key physical properties, showing strong agreement with atomistic simulations. It is important to emphasize that, although we used atomistic simulations as the reference in this study, the proposed methodology can be readily adapted to optimize CG models against experimental data.

This work is motivated by recent advances in extending BO to high-dimensional settings, highlighting its potential beyond low-dimensional parameter tuning. As BO techniques continue to evolve, particularly with the development of trust-region strategies and multifidelity frameworks, we anticipate that these tools will further improve the scalability and efficiency of CG model optimization, paving the way for faster, data-driven materials discovery and molecular design.

Supplementary Material