Accurate coarse-graining of small organic molecules in melts and thin films using density-dependent potentials

TL;DR

This work addresses the challenge that standard coarse-grained (CG) models often fail to capture interfacial inhomogeneity in melts and thin films of small organic molecules. It introduces local-density-dependent potentials (LDPs) within a bottom-up force-matching CG framework, decomposing the non-bonded interactions as $U_ ext{nb}(oldsymbol{R}) = U_2(oldsymbol{R}) + U_ ho(oldsymbol{R})$ and using local densities $ ho_{ au|I}$ with weights $w_{ au|t_I}$ and lengths $L_{ au| au'}$. The results show that these CG models reproduce bulk density, compressibility, radial distribution functions, and interfacial orientational order, with good transferability across temperatures; for thin films, they capture density profiles, molecular orientations, and, in some cases, require cross-pair LD terms to accurately describe cohesion at exposed sites. The approach yields speed-ups of roughly two orders of magnitude relative to all-atom simulations, enabling efficient exploration of processing conditions and non-equilibrium phenomena like vapor deposition while maintaining atomically resolved interfacial detail. Overall, density-dependent CG models offer a reliable, scalable route to design and predict interfacial behavior in conjugated organic melts and thin films.

Abstract

Conjugated organic molecules play a central role in a wide range of optoelectronic devices, including organic light-emitting diodes, organic field-effect transistors, and organic solar cells. A major bottleneck in the computational design of these materials is the discrepancy between simulation and experimental time and length scales. Coarse-graining (CG) offers a promising solution to bridge this gap by reducing redundant degrees of freedom and smoothing the potential energy landscape, thereby significantly accelerating molecular dynamics simulations. However, standard CG models are typically parameterized from homogeneous bulk simulations and assume density-independent effective interactions. As a consequence, they often fail to replicate inhomogeneous systems, such as (free-standing) thin films, due to an incorrect representation of liquid-vacuum interfacial properties. In this work, we develop a CG parametrization strategy that incorporates local-density-dependent potentials to capture material heterogeneities. We evaluate the methodology by simulating free-standing films and comparing interfacial orientational order parameters between all-atom and CG simulations. The resulting CG models accurately reproduce bulk densities and radial distribution functions as well as molecular orientations at the liquid-vacuum interface. This work paves the way for reliable, computation-driven predictions of atomically resolved interfacial ordering in organic molecular systems.

Figures (8)

• Figure 1: Chemical structures of molecules with CG mapping schemes. (a-c) Low resolution two-site models, (d,e) Low resolution four-site models, and (f,g) high resolution multi-site models.
• Figure 2: Comparison of probability distribution $P$ of bonded particles from AA mapped (solid red lines) and CG simulations (dashed blue lines) for (a) mCBP and (b-d) Tm3PyPB. The black lines show the bonded interactions used in the CG simulations. All data computed at $T=550\,\text{K}$.
• Figure 3: (a,b) Radial distribution functions $g(r)$ (right $y$-axes) and weight functions $w(r)$ (left $y$-axes) from the AA simulations of mCBP, using the (a) two-site and (b) eight-site mapping for various cutoff radii $r_\text{c}$, as indicated. (c,d) Local-density distributions $P(\rho)$ corresponding to the choices of $w(r)$ shown in (a,b). The vertical dashed lines indicate the mean number densities $\bar{\rho} \approx 2.45\,\text{nm}^{-3}$ and $\approx 7.88\,\text{nm}^{-3}$ of the two- and eight-site model, respectively. All simulations performed in bulk at $T=550\,\text{K}$.
• Figure 4: (a) Average mass density $\bar{\rho}_\text{m}$ at different temperatures $T$ as functions of cutoff radius $r_\text{c}$. The horizontal dotted line represents the reference AA mass density, $\bar{\rho}^\text{AA}_\text{m}$. The arrows indicate the chosen optimum cutoff radius $r_\text{c}^*$. (b) Pair potential $U_2(r)$ and (c) LDP $U_\rho$ at different $T$, optimized at $r_\text{c}^*$. The dashed lines indicate the extrapolated $\rho$-region of $U_\rho$. The two vertical lines indicate $\bar{\rho}^\text{CG}_\text{m}$ from CG simulations at $T=480\,\text{K}$ and $680\,\text{K}$, respectively. (d) Local-density distributions $P(\rho)$ from mapped AA simulations (solid lines) and from CG simulations (dashed lines with symbols). (e) Intermolecular radial distribution function $g(r)$ from mapped AA simulations (solid lines) and from CG simulations (symbols). All data for two-site mCBP model.
• Figure 5: Comparison of radial distribution functions $g(r)$ from CG and mapped AA simulation for (a-c) TPBi and (d-e) Tm3PyPB at $T=550\,\text{K}$.