Anisotropic molecular coarse-graining by force and torque matching with
  neural networks
We develop a machine-learning method for coarse-graining condensed-phase
molecular systems using anisotropic particles. The method extends currently
available high-dimensional neural network potentials by addressing molecular
anisotropy. We demonstrate the flexibility of the method by parametrizing
single-site coarse-grained models of a rigid small molecule (benzene) and a
semi-flexible organic semiconductor (sexithiophene), attaining structural
accuracy close to the all-atom models for both molecules at considerably lower
computational expense. The machine-learning method of constructing the
coarse-grained potential is shown to be straightforward and sufficiently robust
to capture anisotropic interactions and many-body effects. The method is
validated through its ability to reproduce the structural properties of the
small molecule's liquid phase and the phase transitions of the semi-flexible
molecule over a wide temperature range.
Marltan O. Wilson, David M. Huang
# Anisotropic molecular coarse-graining by force and torque matching with neural networks

###### Abstract

We develop a machine-learning method for coarse-graining condensed-phase molecular systems using anisotropic particles. The method extends currently available high-dimensional neural network potentials by addressing molecular anisotropy. We demonstrate the flexibility of the method by parametrizing single-site coarse-grained models of a rigid small molecule (benzene) and a semi-flexible organic semiconductor (sexitiophene), attaining structural accuracy close to the all-atom models for both molecules at considerably lower computational expense. The machine-learning method of constructing the coarse-grained potential is shown to be straightforward and sufficiently robust to capture anisotropic interactions and many-body effects. The method is validated through its ability to reproduce the structural properties of the small molecule's liquid phase and the phase transitions of the semi-flexible molecule over a wide temperature range.

## I Introduction

Machine learning is quickly becoming an invaluable tool in the search, analysis, and development of new materials.[1; 2] Neural networks, in particular, have had major recent success in areas ranging from predicting the folded conformation of biological macromolecules such as proteins[3] to developing highly accurate temperature-transferable interatomic potentials.[4; 5] The latter is an important advance in the field of molecular dynamics (MD) simulations. Improvements in these machine-learning models aim to expand the length and time scale of simulations without sacrificing accuracy.[6; 7] Currently used ab initio molecular dynamics simulation models are generally accurate but are computationally expensive, limiting their ability to probe long time scales.[8; 9] However, machine-learning potentials can produce ab initio accuracy at the computational cost of classical atomistic models.[10; 11]

Even though simulations at the classical MD level are faster than ab initio MD, the speedup is still insufficient to model the large length and time scales needed to fully understand certain phenomena and processes such as supramolecular assembly. It is well known that explicit modeling of high-frequency motion is not critical for describing many phenomena in molecular systems. These simplifications have led to the development of molecular coarse-grained models to study large, complex materials and biological systems.[12] Parameterization of coarse-grained interaction potentials commonly uses one of two strategies: the top-down approach in which parameters are tuned to match macroscopic observables, as exemplified by the Martini model,[13] and the bottom-up approach in which interactions are derived from the properties of a fine-grained model with more degrees of freedom.[12] By following a similar bottom-up process used to apply machine learning to ab initio MD data, neural-network approaches have been extended to coarse-grained molecular models, further extending the length and time scale of simulations with atomistic accuracy.[14; 15]

Neural-network potentials using isotropic coarse-grained particles have several advantages over their pair-wise additive analytical counterparts since they are constructed as many-body potentials. This many-body potential can become costly when multiple coarse-grained particles are needed to preserve the shape anisotropy of a molecular fragment. It is sometimes more accurate and computationally efficient to represent these groups of atoms as a single anisotropic coarse-grained particle such as an ellipsoid, such as in the case of large, rigid, anisotropic molecular fragments. Analytical anisotropic coarse-grained potentials such as the Gay-Berne potential[16; 17] were developed to address the poor performance of spherically symmetric potentials in replicating intrinsic anisotropic interactions such as \(\pi\)-stacking. By modeling rigid anisotropic groups of atoms as ellipsoids, the anisotropic properties of the group are preserved in a single-site model. Shape and interaction anisotropy is especially important for the study of organic semiconductor molecules, which typically consist of highly anisotropic and rigid \(\pi\)-conjugated units and often form liquid-crystal phases whose morphology strongly affects their performance in devices such as solar cells, transistors, and light-emitting diodes.[18]

Unlike analytical pair-wise additive potentials such as the Gay-Berne potential, high-dimensional neural-network potentials are constructed based on the immediate neighborhood of a molecule and thus account for many-body effects as well as local density variations. Notable machine-learning implementations of interatomic and intermolecular potentials include the neural-network potentials developed by Behler et al.[19] The Behler neural-network potentials are constructed from a set of symmetry functions used to represent the invariant properties of the atomic environment of each atom taken from ab initio simulations. DeepMD[10] and DeepCG[14] are two other neural-network codes constructed for atomistic and coarse-grained simulations, respectively. All these neural-network potentials rely on an invariant representation of the atomic/molecular environment. The CGnets deep-learning approach[15] employs a prior potential to account for areas in a coarse-grained data set that may not be properly sampled due to high repulsive energies. These interactions are especially important to reproduce the local structure of the simu
lated material.

Machine learning has previously been applied to the parameterization of coarse-grained models with anisotropic particles,[20] but no such implementation has used a nonlinear neural-network optimization method to construct the coarse-grained potential. In this work, we address this gap in knowledge by using a neural network to construct a high-dimensional anisotropic coarse-grained potential. We parameterize the neural-network potential using a recently derived systematic and general bottom-up coarse-graining method called anisotropic force-matching coarse-graining (AFM-CG),[21] which generalizes the multi-scale coarse-graining (MS-CG) method[22] for isotropic coarse-grained particles to anisotropic particles. The method rigorously accounts for finite-temperature, many-body effects without assuming a specific functional form of the anisotropic coarse-grained potential. It yields general equations relating the forces, torques, masses, and moments of inertia of the coarse-grained particles to properties of a fine-grained (e.g. all-atom) MD simulation based on a mapping between fine-grained and coarse-grained coordinates and momenta, and by matching the equilibrium coarse-grained phase-space distribution with the mapped distribution of the fine-grained system. The previous implementations of the AFM-CG method approximated the coarse-grained potential as a sum of pair interactions between particles.[21] Here, we extend this approach to more general many-body anisotropic interactions described by a neural-network potential. We also extend the approach, which was derived for constant-volume systems in the canonical ensemble to constant-pressure systems by applying a virial-matching condition previously derived for the MS-CG method.[23]

A general coarse-grained potential should capture any temperature-dependent phase transition associated with either melting, annealing, or glass transition temperatures as well as the local structure and density of the material. The focus is on the development of a model for which trained parameters can be easily obtained and one capable of reproducing interaction anisotropy, temperature transferability, and many-body effects. The flexibility of the new model is demonstrated through the matching of structural and thermodynamic properties of condensed-phase systems of a small anisotropic molecule, benzene, and of a larger, more flexible organic semiconductor molecule, sexithiophene. These two molecules were chosen to determine the conditions under which coarse-grained structural inaccuracy outweighs the computational efficiency of a single-anisotropic-site model.

## II Theory

The key aspects of the theory that underpins the AFM-CG method and its extension to constant pressure via virial matching are summarized below. The reader is referred to Ref. [21] for a more detailed description and the full derivation of the method.

The positions \(\mathbf{r}^{n}=\{\mathbf{r}_{1},\mathbf{r}_{2},\ldots,\mathbf{r}_{n}\}\) of the \(n\) fine-grained particles are mapped onto the positions \(\mathbf{R}^{N}=\{\mathbf{R}_{1},\mathbf{R}_{2},\ldots,\mathbf{R}_{N}\}\) and orientations \(\mathbf{\Omega}^{N}=\{\mathbf{\Omega}_{1},\mathbf{\Omega}_{2},\ldots,\mathbf{\Omega}_{N}\}\) of the \(N\) anisotropic coarse-grained particles. Each fine-grained particle \(i\) is mapped to a single coarse-grained particle by defining \(N\) non-intersecting subsets, \(\zeta_{1},\zeta_{2},\ldots,\zeta_{N}\), of the fine-grained particle indices such that \(\zeta_{I}\) contains the indices of the fine-grained particles mapped onto coarse-grained particle \(I\). The position \(\mathbf{R}_{I}\) of coarse-grained particle \(I\) is defined to be equal to the center-of-mass of the group of fine-grained particles that are mapped onto it, i.e.

\[\mathbf{R}_{I}=\frac{\sum_{i\in\zeta_{I}}m_{i}\mathbf{r}_{i}}{\sum_{i\in\zeta_{I}}m_{ i}}, \tag{1}\]

where \(m_{i}\) is the mass of fine-grained particle \(i\). The orientation

\[\mathbf{\Omega}_{I}=\begin{bmatrix}\mathbf{\Omega}_{I,1}\\ \mathbf{\Omega}_{I,2}\\ \mathbf{\Omega}_{I,3}\end{bmatrix} \tag{2}\]

of coarse-grained particle \(I\) is specified by the rotation matrix whose components are the particle's three normalized principal axes of inertia, \(\mathbf{\Omega}_{I,q}\) for \(q=1,2,3\). These axes are defined to be equal to the corresponding principal axes relative to the center-of-mass of the group of fine-grained particles that are mapped onto the coarse-grained particle. Thus, these axes are the normalized eigenvectors of the inertia tensor

\[\mathbf{\mathrm{I}}_{\mathrm{FG},I}=\sum_{i\in\zeta_{I}}m_{i}(||\mathbf{\alpha}\mathbf{r}_ {i}||^{2}\mathbf{E}-\mathbf{\alpha}\mathbf{r}_{i}\mathbf{\Delta}\mathbf{r}_{i}^{\mathrm{T}}), \tag{3}\]

where \(\mathbf{\Delta}\mathbf{r}_{i}=\mathbf{r}_{i}-\mathbf{R}_{I}\) is the position of fine-grained particle \(i\) relative to the center-of-mass (coarse-grained particle position) and \(\mathbf{E}\) is the \(3\times 3\) identity matrix. From these coordinate mappings and the relationship between generalized coordinates and momenta from Hamilton's equations,[24] mappings from the linear momenta \(\mathbf{p}^{n}=\{\mathbf{p}_{1},\mathbf{p}_{2},\ldots,\mathbf{p}_{n}\}\) of the fine-grained particles to the linear momenta \(\mathbf{P}^{N}=\{\mathbf{P}_{1},\mathbf{P}_{2},\ldots,\mathbf{P}_{N}\}\) and angular momenta \(\mathbf{L}^{N}=\{\mathbf{L}_{1},\mathbf{L}_{2},\ldots,\mathbf{L}_{N}\}\) of the anisotropic coarse-grained particles can also be defined.[21] The mappings for coarse-grained particle \(I\) are

\[\mathbf{P}_{I}=\frac{M_{I}}{\sum_{i\in\zeta_{I}}m_{i}}\sum_{i\in\zeta_{I}}\mathbf{p}_{i} \tag{4}\]

and

\[\mathbf{L}_{I}=\mathbf{\mathrm{I}}_{I}\mathbf{\mathrm{F}}_{\mathrm{FG},I}^{-1}\sum_{i\in \zeta_{I}}\mathbf{\Delta}\mathbf{r}_{i}\mathbf{\times}\mathbf{p}_{i}, \tag{5}\]

respectively, where \(\mathbf{\mathrm{I}}_{I}\) is the inertia tensor of coarse-grained particle \(I\).

Given these mappings, several conditions can be derived that the coarse-grained model must satisfy for its equilibrium coarse-grained phase-space distribution to match the corresponding mapped distribution of the fine-grained system. Consistency between the configuration-space distributions gives the following matching conditions between the forces \(\mathbf{F}_{I}\) and torques \(\mathbf{\tau}_{I}\) on coarse-grained particle \(I\) and the forces on the fine-grained particles mapped onto it:[21]

\[\mathbf{F}_{I}(\mathbf{R}^{N},\mathbf{\Omega}^{N})=-\frac{\partial U}{\partial R_{I}}= \left\langle\sum_{i\in\zeta_{I}}\mathbf{f}_{i}\right\rangle_{\mathbf{R}^{N},\mathbf{\Omega} ^{N}} \tag{6}\]
and

\[\mathbf{\tau}_{I}(\mathbf{R}^{N},\mathbf{\Omega}^{N})=-\sum_{q}\Omega_{I,q}\times\frac{\partial U }{\partial\Omega_{I,q}}=\left\langle\sum_{i\in\mathcal{U}_{I}}\Delta\mathbf{r}_{i} \times\mathbf{f}_{i}\right\rangle_{\mathbf{R}^{N},\mathbf{\Omega}^{N}} \tag{7}\]

where \(U(R^{N},\mathbf{\Omega}^{N})\) is the coarse-grained potential, \(\mathbf{f}_{i}(\mathbf{r}^{n})=-\frac{\partial\mathbf{n}}{\partial\mathbf{r}_{i}}\) is the force on fine-grained particle \(i\), with \(\mathbf{u}(\mathbf{r}^{n})\) the fine-grained potential and \(\langle\cdots\rangle_{\mathbf{R}^{N},\mathbf{\Omega}^{N}}\) denoting an average over fine-grained configurations mapped to coarse-grained configuration \((R^{N},\mathbf{\Omega}^{N})\).

Consistency between the momentum-space distributions requires the mass \(M_{I}\) of coarse-grained particle \(I\) to be the sum of the masses of its constituent fine-grained particles, i.e. [21]

\[M_{I}=\sum_{i\in\mathcal{U}_{I}}m_{i}. \tag{8}\]

In addition, provided that the inertia tensor \(\mathbf{\mathrm{I}}_{\mathrm{FG},J}\) of the group of fine-grained particles mapped to this coarse-grained particle does not depend on the configuration of the other particles [21],

\[I_{I,q}^{1/2}\exp\left(-\frac{I_{I,q}\omega_{I,q}^{2}}{2k_{\mathrm{B}}T} \right)\approx\left\langle I_{\mathrm{FG},I,q}^{1/2}\exp\left(-\frac{I_{ \mathrm{FG},I,q}\omega_{I,q}^{2}}{2k_{\mathrm{B}}T}\right)\right\rangle_{\mathbf{R }_{I},\mathbf{\Omega}_{I}}, \tag{9}\]

where \(I_{I,q}\), \(I_{\mathrm{FG},I,q}\), and \(\omega_{I,q}\) are the components of the coarse-grained moment of inertia, fine-grained moment of inertia, and angular velocity about the \(q\) axis, and \(\langle\cdots\rangle_{\mathbf{R}_{I},\mathbf{\Omega}_{I}}\) denotes an equilibrium average of fine-grained configurations consistent with the coordinate mapping of coarse-grained particle \(I\). Furthermore, if the fluctuations in \(I_{\mathrm{FG},I,q}\) are small compared to its mean, it can be shown that [21]

\[I_{I,q}\approx\left\langle I_{\mathrm{FG},I,q}\right\rangle_{\mathbf{R}_{I},\mathbf{ \Omega}_{I}}, \tag{10}\]

i.e. the principal moment of inertia of a coarse-grained particle about each principal axis \(q\) is approximately equal to the equilibrium average of the corresponding principal moment of the fine-grained particles mapped onto it.

The AFM-CG method was derived only for the constant-volume conditions of the canonical ensemble, but is straightforwardly generalized to constant-pressure conditions by analogy with the MS-CG method for spherical coarse-grained particles in the isothermal-isobaric ensemble [23]. Thus, the force- and torque-matching conditions at constant pressure are the same as those in Eqs. (6) and (7), except that the coarse-grained forces, torques, and potential are in general functions of the coarse-grained system volume \(V\) and the equilibrium average is constrained to configurations in which the fine-grained system volume \(v=V\). The coarse-grained potential must also satisfy a virial-matching condition [23],

\[W(\mathbf{R}^{N},\mathbf{\Omega}^{N},V) =-\frac{\partial U}{\partial V}\] \[=\left\langle\frac{(n-N)k_{\mathrm{B}}T}{v}+\frac{1}{3v}\sum_{l=1 }^{n}\mathbf{f}_{i}\cdot\mathbf{r}_{i}\right\rangle_{\mathbf{R}^{N},\mathbf{\Omega}^{N},V} \tag{11}\]

In summary, for the equilibrium phase-space distribution of the coarse-grained model to match that of the fine-grained model in the isothermal-isobaric ensemble, the coarse-grained potential should satisfy Eqs. (6), (7), and (11), while the coarse-grained masses and principal moments of inertia should satisfy Eqs. (8) and (9), respectively. As shown below, using the more approximate Eq. (10) to parameterize the moments of inertia gives almost the same results as Eq. (9), even for a flexible molecule, so we have used this simpler equation for parameterization later on.

## III Methods

### Force-, torque-, and virial-matching algorithm

The analytical expression for the coarse-grain potential \(U\) is not usually known. However, an approximation to the functional form can be obtained using a neural-network optimization algorithm with Eqs. (6), (7), and (11) acting as necessary constraints. In general, \(U(\mathbf{R}^{N},\mathbf{\Omega}^{N},V)\) is a function of the particle configuration and system volume. In this work, we have assumed that \(U\) does not depend explicitly on \(V\), in which case [23]

\[\frac{\partial U}{\partial V}=\frac{1}{3V}\sum_{l=1}^{N}\frac{\partial U}{ \partial\mathbf{R}_{l}}\cdot\mathbf{R}_{l}. \tag{12}\]

With this approximation, the virial-matching condition in Eq. (11) can be written, using \(v=V\), as

\[-\sum_{l=1}^{N}\frac{\partial U}{\partial\mathbf{R}_{l}}\cdot\mathbf{R}_{l}=\left\langle 3 (n-N)k_{\mathrm{B}}T+\sum_{i=1}^{n}\mathbf{f}_{i}\cdot\mathbf{r}_{i}\right\rangle_{\mathbf{R }^{N},\mathbf{\Omega}^{N},V}. \tag{13}\]

Despite this approximation, we show later on that the coarse-grained models parameterized accurately match the average density of the corresponding all-atom fine-grained system at constant pressure.

To ensure that all equivalent configurations are assigned the same position in coordinate space, a transformation was made from the set of Cartesian coordinates to a vector \(\mathbf{D}_{IJ}\) that was invariant under translation, rotation, and permutation of any pair of coarse-grained particles \(I\) and \(J\)[25, 26, 27, 10], which was defined in terms of the positions, \(\mathbf{R}_{I}\) and \(\mathbf{R}_{J}\), and orientations, \(\mathbf{\Omega}_{I}\) and \(\mathbf{\Omega}_{J}\), of the two particles by

\[\mathbf{D}_{IJ}=\{R_{IJ},\mathbf{R}_{IJ}\cdot\mathbf{\Omega}_{I,1},\mathbf{R}_{IJ }\cdot\mathbf{\Omega}_{I,2},\mathbf{R}_{IJ}\cdot\mathbf{\Omega}_{I,3},\] \[\mathbf{R}_{IJ}\cdot\mathbf{\Omega}_{J,1},\mathbf{R}_{IJ}\cdot\mathbf{\Omega}_{J, 2},\mathbf{R}_{IJ}\cdot\mathbf{\Omega}_{J,3},\] \[\mathbf{\Omega}_{I,1}\cdot\mathbf{\Omega}_{J,1},\mathbf{\Omega}_{I,1}\cdot\bm {\Omega}_{J,2},\mathbf{\Omega}_{I,1}\cdot\mathbf{\Omega}_{J,3},\] \[\mathbf{\Omega}_{I,2}\cdot\mathbf{\Omega}_{I,1},\mathbf{\Omega}_{I,2}\cdot\bm {\Omega}_{J,2},\mathbf{\Omega}_{I,2}\cdot\mathbf{\Omega}_{J,3},\] \[\mathbf{\Omega}_{I,3}\cdot\mathbf{\Omega}_{J,1},\mathbf{\Omega}_{I,3}\cdot\bm {\Omega}_{J,2},\mathbf{\Omega}_{I,3}\cdot\mathbf{\Omega}_{J,3}\}, \tag{14}\]

where \(R_{IJ}\equiv\|\mathbf{R}_{IJ}\|\), \(\mathbf{R}_{IJ}\equiv\mathbf{R}_{I}-\mathbf{R}_{J}\) and \(\mathbf{\Omega}_{I}\) and \(\mathbf{\Omega}_{J}\) are specified by rotation matrices of the form of Eq. (2). The coordinates of each neighbor within the cut-off distance of particle \(I\) were transformed to a \(\mathbf{D}_{IJ}\) vector. All the \(\mathbf{D}_{IJ}\) vectors for a given neighborhood were concatenated into a 2D matrix
of size \(N\times\dim(\mathbf{D}_{IJ})\) representing a unique configurational fingerprint for coarse-grained particle \(I\).

The potential function could then be written in terms of a set of neural network trainable parameters and activation functions transforming \(\mathbf{\mathsf{D}}_{I}\) to a potential energy value. While \(\mathbf{\mathsf{D}}_{I}\) is a sufficient specification of the coarse-grained coordinates to enforce relevant invariant properties of the molecular environment, it does not possess all the symmetries of the potential energy surface that it aims to fit.[25; 28] For each molecular environment, it was assumed that the interactions were predominantly short-ranged such that neighbors beyond a certain cut-off distance, \(R_{\mathrm{c}}\), did not contribute to the potential.[19] This condition can be enforced by a cut-off function of the form

\[g_{\mathrm{c}}(R_{IJ})=\begin{cases}\frac{1}{2}\left[\cos\left(\frac{\pi R_{ IJ}}{R_{\mathrm{c}}}\right)+1\right],&R_{IJ}\leq R_{\mathrm{c}},\\ 0,&R_{IJ}>R_{\mathrm{c}}.\end{cases} \tag{15}\]

A set of these cut-off functions can enforce the radial symmetry conditions of the underlying potential energy surface by storing information about the radial distribution of neighbors according to[19]

\[G_{I}^{\mathrm{j}}=\sum_{J\neq I}g_{\mathrm{c}}(R_{IJ}). \tag{16}\]

Continuity of the potential along angular dimensions was ensured by using a compression layer to learn a set of collective variables from vector \(\mathbf{D}_{IJ}\) which were constrained by the well-behaved modified \(G^{5}\) symmetry function[19] given by

\[G_{I}^{5}=\sum_{J\neq I}\prod_{\mu=1}^{M}2^{1-\nu}\left(1+\lambda\cos\theta_{ IJ,\mu}\right)^{\nu}e^{-\eta(R_{IJ}-R_{\mathrm{c}})^{2}}g_{\mathrm{c}}(R_{IJ}), \tag{17}\]

where \(\lambda\in\{-1,1\}\) and \(R_{\mathrm{s}}\), \(\nu\), and \(\eta\) are tunable hyperparameters, while \(\{\cos\theta_{IJ,\mu}\}\), is the set of machine-learned collective variables with the same properties as the angular component of the underlying potential and \(M\) is the total number of machine-learned angular variables. These angular symmetry functions store information about the angular-radial distribution of neighbors in the local environment of coarse-grained particle \(I\). Unlike the case of spherically symmetric particles, in a local reference frame, a neighboring anisotropic particle requires a minimum of seven independent scalar variables to fully describe its position and orientation. However, previous implementations of analytical potentials, including the Gay-Berne potential,[16; 17] have used fewer coordinates for the calculation of the potential and forces. Similarly, for the neural-network potential, an additional compression layer was included to remove the redundant angles from the \(\mathbf{D}_{IJ}\) vectors, since the combination of translation and rotation in 3D is parameterized by at most 7 unique coordinates. The Behler symmetry functions were enforced on the output of the compression layer, ensuring that the learned compression had the same symmetry and continuity of the underlying potential. The reduction in the dimension of \(\mathbf{D}_{IJ}\) also decreases the amount of data that is needed to train a sufficiently accurate potential. By removing the redundant angles in \(\mathbf{D}_{IJ}\) there is a reduced possibility of over-fitting on a small data set.

A set of these symmetry functions with tuned hyperparameters \((\lambda,\nu,\eta,R_{\mathrm{s}},R_{\mathrm{c}})\) can be used to uniquely represent the structural fingerprint of the molecular environment. Symmetry functions used to represent the local environment were constructed using all possible permutations of values from a specified set of hyperparameters. Training of the neural network started with 8 symmetry functions and hyperparameters tuned to minimize the loss function, which is defined below. New symmetry functions were added to the set if they resulted in a significant reduction in the neural-network loss compared with the preceding iteration. The set of hyperparameters in the symmetry functions used in the anisotropic coarse-grained models parameterized in this work can be found in the Supplementary Material.

To further reduce the amount of data needed to train the neural network, a prior repulsive potential was defined with pairwise additive properties. This potential was used to ensure physical behavior in regions of the potential where the forces are large and thus are rarely sampled in an equilibrium MD simulation. This prior potential only needs to satisfy two conditions: firstly, it must be repulsive at short radial separations, and, secondly, the position of its repulsive barrier must be orientationally dependent. A simple equation satisfying these conditions is

\[U_{\mathrm{prior},I}=\sum_{J\neq I}B_{1}\sigma_{\mathrm{c}}\left(\mathbf{\mathsf{D }}_{I}\right)^{-B_{2}}, \tag{18}\]

where \(\sigma_{\mathrm{c}}\) is a neural-network compression layer function and \(B_{1}\) and \(B_{2}\) are strictly positive trainable parameters. It is also possible to achieve a similar large repulsive barrier through a more advanced nonlinear sampling of the MD simulation data. The prior potential fits the purely repulsive part of the angular-dependent potential to the molecular environment, while the neural-network potential fits the attractive and oscillatory corrections to the environment. The final prediction for the potential energy of the environment of coarse-grained particle \(I\) is therefore the sum of the neural-network potential \(U_{\mathrm{NN},I}\) and the prior repulsive potential \(U_{\mathrm{prior},I}\),[15]

\[U_{I}=U_{\mathrm{NN},I}+U_{\mathrm{prior},I}, \tag{19}\]

and, thus, the total coarse-grained potential is

\[U=\sum_{I=1}^{N}U_{I} \tag{20}\]

From the matching conditions in Eqs. (6), (7), and (13), optimization of the neural-network weights and biases requires a loss function of the form

\[L = \left\langle\sum_{I=1}^{N}\left(\alpha\left|\mathbf{F}_{\mathrm{FG},I} +\frac{\partial U}{\partial R_{I}}\right|^{2}+\beta\left|\tau_{\mathrm{FG},J}+ \sum_{q}\Omega_{I,q}\times\frac{\partial U}{\partial\Omega_{I,q}}\right|^{2} \right)\right.\] \[+ \left.\gamma\left|3\left(n-N\right)k_{\mathrm{B}}T+\sum_{I=1}^{N} \left(\Psi_{\mathrm{FG},J}+\frac{\partial U}{\partial\mathbf{R}_{I}}\cdot\mathbf{R}_{I }\right)\right|^{2}\right\rangle_{R^{N},\Omega^{N},V}, \tag{21}\]
where

\[\mathbf{F}_{\text{FG},I}\equiv\sum_{i\in\zeta_{I}}\mathbf{f}_{i},\,\tau_{\text{FG},I}\equiv \sum_{i\in\zeta_{I}}\Delta\mathbf{r}_{i}\times\mathbf{f}_{i},\,W_{\text{FG},I}\equiv \sum_{i\in\zeta_{I}}\mathbf{f}_{i}\cdot\mathbf{r}_{i}, \tag{22}\]

and \(\alpha,\beta\), and \(\gamma\) are weights. These weights specify the fraction of each loss that is used for backpropagation and were free to change with the learning rate during optimization.[14] Even though there have been significant efforts to develop methods to fit the constrained averaged coarse-grained forces directly,[29; 30] the average total fine-grained forces subject to the constraint of matching fine-grained and coarse-grained configurations are not easily obtained. An indirect means of minimizing the loss function in Eq. (21) above is possible by replacing the constrained ensemble average with an average over instantaneous unconstrained simulation configurations,[14]

\[L_{\text{inst}} =\sum_{t=1}^{N_{\text{t}}}\left[\sum_{I=1}^{N}\left(\alpha\left| \mathbf{F}_{\text{FG},I}(\mathbf{r}_{t}^{n})+\frac{\partial U(\mathbf{\xi}_{I})}{\partial R _{I}}\right|^{2}\right.\right.\] \[+\left.\left.\beta\left|\mathbf{\tau}_{\text{FG},I}(\mathbf{r}_{t}^{n}) +\sum_{q}\Omega_{I,q}(\mathbf{\xi}_{I})\right.\right)\times\frac{\partial U(\mathbf{ \xi}_{I}))}{\partial\Omega_{I,q}}\right|^{2}\right)\] \[+\left.\left.\gamma\left|3\left(n-N\right)k_{\text{B}}T+\sum_{I= 1}^{N}\left(\tilde{W}_{\text{FG},I}(\mathbf{r}_{t}^{n})+\frac{\partial U(\mathbf{ \xi}_{I})}{\partial\mathbf{R}_{I}}\cdot\mathbf{R}_{I}(\mathbf{\xi}_{I})\right)\right|^{2 }\right], \tag{23}\]

since it can be shown, for a sufficiently large dataset that comprehensively samples the equilibrium ensemble of the fine-grained system, that \(L\) and \(L_{\text{inst}}\) have the same global minimum. Here, \(N_{\text{t}}\) is the number of simulation configurations in the dataset, \(\mathbf{r}_{t}^{n}\) and \(v_{t}\) are the fine-grained coordinates and system volume for configuration \(t\), and \(\mathbf{\xi}_{I}=(\mathbf{R}^{N}(\mathbf{r}_{t}^{n}),\mathbf{\Omega}^{N}(\mathbf{r}_{t}^{n}),V(v_{t }))\) is the mapped coarse-grained configuration for this fine-grained configuration. The loss function was optimized using the minibatch gradient descent as implemented in TensorFlow.

The feedforward neural network shown in Fig. 1 was then trained, where the forward propagation used matrix \(\mathbf{\mathbbm{D}}_{I}\) as an input to predict the coarse-grained potential \(U\), after which TensorFlow's computational derivative was used to calculate the outputs, namely the predicted forces, torques, and virial. In the backpropagation stage, the loss function was used to calculate the error between the true and predicted values, which was then used to update the network weights and biases. The errors between the true and predicted parameters were calculated using TensorFlow's mean squared error, and gradient descent was implemented using TensorFlow's Adam optimizer.[31] Once the error of the neural network was minimized, the neural-network model was used to predict the forces, torques, and virial. However, removing the output and derivative layers gives access to the predicted potential of mean force. By optimizing the partial derivatives of the potential instead of the potential itself, by the nature of the derivative, there will be less oscillation in the potential at the edges of the data set close to the cut-off distances.

### LAMMPS modification and neural network implementation

The neural network was constructed in Tensorflow (version 2.3.0)[32] using the Keras (version 2.4.3) functional API[33] and saved using the Tensorflow SavedModel format. The trained neural network was implemented in LAMMPS using the Tensorflow C API and cppflow wrapper. All simulations were carried out using the LAMMPS MD software package (version 20Nov19).[34; 35; 36] The Optimized Potentials for Liquid Simulations-All Atom (OPLS-AA) force field[37; 38; 39; 40] was used for all all-atom simulations with a cut-off distance of 10 A for short-ranged non-bonded interactions; long-ranged electrostatic interactions were calculated with the particle-particle particle-mesh (PPPM) method[36; 41] The bonds that include hydrogen were constrained using the SHAKE algorithm.[42] Simulations were carried out in the isothermal-isobaric (NPT) ensemble at a pressure of 1 atm, with the temperature and pressure controlled by a Nose-Hoover thermostat and barostat.[43; 44]

Neural-network training was carried out using data from a 25 ns all-atom simulation in which simulation configurations and forces and velocities were saved at 2 ps intervals. The simulation snapshots from the last 20 ns were shuffled and then divided into 4 groups of equal size, \(\{g_{0},g_{1},g_{2},g_{3}\}\). The neural network was initially trained on \(g_{0}\) and validated on \(g_{3}\). The validation set \(g_{3}\) was further divided into an 8:2 ratio where the lesser was reserved as the test set. New snapshots were added from \(g_{1}\) and \(g_{2}\) if the mean errors of their predicted forces and torques were larger than that of the test set. The accuracy of the trained neural network was then compared to the expected accuracy determined from k-fold cross-validation.[45; 46] During k-fold cross-validation, the last 20 ns of simulation data was shuffled and divided into 10 folds, \(\{\psi_{0},...,\psi_{b}\}\). The model was validated on \(\psi_{i}\) and trained on \(\bigcup_{j\neq i}\psi_{j}\) for all \(i,j\in\{0,1,\ldots,9\}\). The loss of the iterative training method was found to be identical to the k-fold cross-validation loss.

The coarse-grained simulations were done using a modified version of LAMMPS in which the trained neural network was

Figure 1: Schematic of anisotropic force-matching neural network architecture.
introduced to calculate the forces and energies. The dimensions of the coarse-grained sites used in the simulations were derived from the inertia tensor of the all-atom model. To test the ability of the coarse-grained model to capture the properties of the all-atom model under a variety of conditions in addition to the single temperature at which the neural network was trained, the equilibrium structural properties of equivalent coarse-grained and all-atom systems were compared in simulations at several different temperatures. In all cases, the total length of the coarse-grained simulation was 25 ns, with the last 20 ns used to calculate structural and thermodynamic properties.

## IV Results and Discussion

To demonstrate the flexibility of the method, we have used our neural-network model to construct coarse-grained interaction potentials for benzene, an archetypal anisotropic small molecule, and \(\alpha\)-sexithiophene, an organic semiconductor with significant applications in organic electronic devices [47, 48, 49] (Fig. 2). These molecules were selected to demonstrate the neural network's ability to handle anisotropic molecules of varying complexity, flexibility, and aspect ratio while still reproducing the structural and phase behavior.

The shape of a coarse-grained particle obtained from the anisotropic coarse-graining method is determined by the "average" shape of the fine-grained molecule or molecular fragment that is mapped to it under the parameterization conditions. Thus, the variation of the aspect ratio of the molecule or molecular fragment with temperature in the all-atom simulations can potentially be a qualitative indicator of the temperature transferability of the coarse-grained model. Here, the aspect ratio of the molecule was calculated as the ratio of the length to the breadth of the molecule, where the length was defined as the longest principal axis and the breadth was defined as the sum of the remaining two semi-axes. Unlike benzene, the thiophene-thiophene dihedral angles also have a temperature-dependent effect on the aspect ratio of sexithiophene.

Neural networks in general are very good at interpolation but struggle with extrapolation.[50, 51, 52, 53] The accuracy of the model is therefore expected to decrease as the aspect ratio of the molecule deviates from that at the parameterization temperature, as well as when the density distribution is sufficiently different from the parameterization temperature. By parameterizing the systems in the liquid phase, the model can capture a wider variety of fluctuations in the density of the system and the dimensions of the molecules. The average size of a flexible molecule in the isotropic phase will be different from the size of the molecule when locked in a rigid crystal structure.[54, 55] However, this temperature-dependent size difference should decrease with increased rigidity of the molecule.

### Benzene

Simulations consisting of 500 benzene molecules were carried out at 280, 300, 320, 330, and 350 K, and the coarse-grained neural-network model was parameterized at 300 K. The time step was 2 fs in the all-atom simulations and 12 fs in the coarse-grained simulations. The cut-off distance hyperparameter \(R_{\mathrm{c}}\) was 10 A. The root mean squared validation error for the forces was \(2.55\,\mathrm{kcal}\,\mathrm{mol}^{-1}\,\mathrm{\AA}^{-1}\) and that of the torque was \(4.35\,\mathrm{kcal}\,\mathrm{mol}^{-1}\). The average post-training error in the pressure was 0.0092 atm.

The average principal moments of inertia in the all-atom simulation at 300 K were used to determine the principal moments of the coarse-grained model using Eqn. (10) (values given in the Supplementary Material) since fluctuations in the moments at the parameterization temperature were small.[21] The variation of the molecular aspect ratio of the all-atom benzene model with temperature is shown in Fig. 3. The distribution of possible dimensions observed for benzene is narrow and remains fairly constant with temperature, making benzene an ideal case where molecular flexibility does not contribute significantly to the overall error of the model.[56]

Fig. 4 shows that the coarse-grained neural-network model accurately captures the liquid density of the all-atom model over a wide range of temperatures from just above the freezing point to just below the boiling point, with only slight deviations for the temperature furthest from the parameterization temperature. As shown in Fig. 5, the coarse-grained model also accurately predicts the radial distribution function (RDF) of the all-atom model over the same temperature range.

To further elucidate the accuracy of the neural-network coarse-grained model, the angular-radial distribution function (ARDF) was analyzed. The ARDF is defined by

\[g(r,\theta)=\frac{\langle n(r,\theta)\rangle}{\frac{4}{3}\pi\rho[(r+\Delta r) ^{3}-r^{3}]\sin\theta\Delta\theta}, \tag{24}\]

where \(\langle n(r,\theta)\rangle\) is the average number of molecules in the spherical shell within the bounds \(r\) to \(r+\Delta r\) of the center-of-mass of a chosen molecule and having an out-of-plane axis rotation of \(\theta\) with respect to the out-of-plane axis of the chosen molecule,[57] and \(\rho\) is the bulk number density. Fig. 6

Figure 2: Chemical structures of (a) benzene and (b) \(\alpha\)-sexithiophene with a coarse-grained ellipsoid superimposed on a possible configuration of each molecule.
shows the 2D heatmap of the ARDF along with 1D slices of this function at specific angles at 300 K (the parameterization temperature) for the all-atom and coarse-grained models. The ARDFs at the other simulated temperatures are compared in the Supplementary Material. At all simulated temperatures between 280 and 350 K, the coarse-grained model captures all the major features of the fine-grain structure of the fluid. The only difference is a slight underestimation of the peak heights by the coarse-grained model. The neural-network model is, however, able to more faithfully capture the angular radial distribution of benzene at all temperatures compared with the coarse-grained benzene model previously parameterized with the AFM-CG method using a pair potential to describe the interparticle interactions.[21] This improvement can be attributed to the greater flexibility of the neural-network potential in describing the intermolecular interactions. The neural-network model can demonstrate temperature transferability through careful selection of the neural network hyperparameters to prevent overfitting of the local number density variations.

The coarse-grained simulation of anisotropic molecules using a neural-network potential is more suited for large, preferably rigid, molecules, for which a high degree of coarse-graining can be achieved with reasonable accuracy. However, the model was still able to achieve a modest 20\(\times\) speedup compared with the atomistic simulations, through a combination of reduced computation time per timestep and a larger timestep. This poor performance for a small molecule such as benzene is due to the small reduction in the number of degrees of freedom from the all-atom model to the coarse-grained model, coupled with a neural-network potential that is more computationally expensive than an analytical potential. Nevertheless, computational savings are obtained even in this suboptimal case. Simulations were carried out on a 4-core Intel i7-4790K CPU, but, further speedups could be achieved by taking advantage of the GPU-enabled version of TensorFlow.

Figure 4: Density versus temperature of the all-atom (AA) and coarse-grained (CG) benzene models at 1 atm. Error bars are smaller than the symbols.

Figure 5: Radial distribution function (RDF) of the all-atom (solid lines) and coarse-grained (dashed lines) benzene models at 1 atm and various temperatures. The RDFs have been shifted vertically for clarity.

Figure 3: Length-to-breadth ratio of the all-atom benzene model at 1 atm and various temperatures.
### Sexithiophene

Simulations of 512 sexithiophene molecules were carried out at 570, 590, 640, and 680 K temperatures, corresponding to temperatures previously identified in all-atom MD simulations to correspond to crystalline (K), smectic-A (Sm-A), nematic (N), and isotropic (I) phases respectively.[58] The time step was 1 fs in the all-atom simulations and 12 fs in the CG simulations. Although we have used the OPLS-AA force field for our all-atom simulations, whereas these previous MD simulations[58] used the related AMBER force field[59; 60; 61] the structural properties of systems simulated with these two force fields (in particular the density, orientational order parameter, and radial distribution function discussed below) are very similar for the temperature range studied. The cut-off distance hyperparameter \(R_{\mathrm{c}}\) was set to 21 A. The neural-network model was parameterized using simulation snapshots from the isotropic phase at 680 K, where the molecular mobility was highest. The conditions of the isotropic bulk phase are advantageous in efficiently sampling the configuration space, especially rare high-energy configurations necessary for the accurate reproduction of the repulsive part of the coarse-grained potential.

As shown in Fig. 7a, the distributions of the principal moments of inertia of sexithiophene in the all-atom simulation at the parameterization temperature are broad, indicating that Eqn. (10) may not be adequate for parameterizing the moments of inertia of the coarse-grained model. However, we found that using the more general Eqn. (9) to parameterize the coarse-grained moments of inertia (by fitting the distributions in Fig. 7b-d gave values within <1%. So we used the values from Eqn. (10) in the coarse-grained model.

The root mean squared validation error for the sexithiophene forces was \(3.95\,\mathrm{kcal\,mol^{-1}\,\AA^{-1}}\) and that of the torque was \(9.8\,\mathrm{kcal\,mol^{-1}}\). The sexithiophene final force and torque losses were larger than those of benzene because the model was not complex enough to account for the bending of the polymer and the rotation of the individual thiophene rings. The loss is also skewed to larger values when compared with benzene because sexithiophene is a larger molecule and so the interactions between molecules are stronger overall.

The structural properties of the coarse-grained model were compared with those of its all-atom counterpart at each of the simulated temperatures. The nonlinear change in density with respect to temperature is associated with the phase changes that occur at the simulated temperatures (Fig. 8).[58] The density of the coarse-grained system agrees well with that of the all-atom system, with minimal deviations from the fine-grained system with increasing distance from the parameterization temperature. Compared with benzene, sexithiophene has a much larger change in density between the crystalline and the isotropic phase. This difference results in less overlap between the local density variations in the crystalline phase at the lowest temperature and the training data set in the isotropic phase at the highest temperature. The sexithiophene molecule is also much more flexible than benzene, as seen in the wide distribution of the aspect ratio in the all-atom model at all the simulated temperatures shown in Fig. 9, and its dimensions change significantly with temperature over the range studied. Another limitation of representing sexithiophene as a single-site ellipsoid is the loss of thiophene-thiophene torsional information. That is, for any given position and orientation of the coarse-grained ellipsoid there are multiple different relative orientations between the thiophene groups.[62] This loss of information is significant because the anisotropic interactions of the thiophene subunits are lost, which reduces the neural network's ability to isolate which of the two short axes corresponds to the \(\pi\)-stacking direction.

To further confirm that the density changes were associated with transitions from the crystalline phase through the nematic

Figure 6: Angular–radial distribution function (ARDF) of the all-atom (AA) (top) and coarse-grained (CG) (bottom) benzene models at 300 K and 1 atm depicted as a heat map (left) and 1D slices at constant angle (right). Face-to-face, edge-to-edge, and parallel displaced configurations occur when the angle is 0\({}^{\circ}\), while edge-to-face configurations occur at 90\({}^{\circ}\).
and smectic phases to the isotropic phase, the scalar orientational order parameter \(P_{2}\) was introduced. For a given simulation snapshot at time \(t\), \(P_{2}\) can be found by diagonalizing the ordering matrix

\[\mathbf{Q}=\frac{1}{2N}\sum_{I=1}^{N}{(3\mathbf{u}_{I}\otimes\mathbf{u}_{I}-\mathbf{E})}, \tag{25}\]

Figure 8: Density versus temperature of the all-atom (AA) and coarse-grained (CG) sexithiophene models at 1 atm. Error bars are smaller than the symbols.

Figure 7: (a) Principal moment of inertia distributions for all-atom (AA) sexithiophene model at 680 K and 1 atm. The corresponding angular velocity distributions of each principal axis along with the coarse-grained (CG) fit to the distribution given by Eq. (9) is shown in (b)–(d).

Figure 9: Length-to-breadth ratio of all-atom sexithiophene model at 1 atm and various temperatures. The simulated phase is given in parentheses after each temperature in the legend (I = isotropic, N = nematic, SmA = smectic A, K = crystal).
where \(\mathbf{u}_{I}\) is the unit vector along the molecular axis and \(\mathbf{E}\) is the identity matrix. \(\langle P_{2}\rangle\) is the average over the largest eigenvalue of this matrix for all snapshots of equilibrium configurations.[58] Larger values of the scalar orientational order parameter close to one indicate an ordered crystalline structure while values close to zero correspond to an isotropic disordered phase. The coarse-grained model reproduces the orientational order parameter of the all-atom model reasonably well over the temperature range simulated, as shown in Fig. 10. The coarse-grained model underestimates the degree of orientational ordering observed in the all-atom model away from the parameterization temperature, likely because it does not capture the increasing molecular shape anisotropy that is observed in the all-atom model as the temperature decreases (Fig. 9). As expected, the largest difference occurs in the predicted crystalline phase.

The same trend is seen in the radial distribution functions shown in Fig. 11, in which the agreement between the coarse-grained and all-atom models at most temperatures is excellent, with the largest deviations for the crystalline phase. The underestimation and broadening of the peaks in the crystalline radial distribution function explain the discrepancy between the order parameter of the all-atom and coarse-grained models. The observed differences are most likely due to the effect on molecular packing of the aforementioned discrepancy in molecular shape between the two models as temperature decreases.[63] Nevertheless, even in the crystalline phase, the coarse-grained model captures the peak positions of the radial distribution function very well.

The coarse-grained model also accurately describes orientational correlations in condensed-phase sexithiophene, as illustrated by a comparison with the angular-radial distribution function of the all-atom model. At the parameterization temperature, the coarse-grained model is able to capture all major features when compared to the all-atom model (Fig. 12). The neural-network model is also able to capture the relevant features in the structure of sexithiophene's smectic liquid-crystal phase at 590 K, as shown in Fig. 13. The discrepancies in the width and height of the peaks are likely due to the differences in molecular shape away from the parameterization temperature that was mentioned earlier. The ARDFs of the two models in the nematic phase at 640 K are compared in the Supplementary Material and show similarly good agreement.

Despite sexithiophene not strictly meeting the conditions to be coarse-grained to a single anisotropic particle due to its significant flexibility, the coarse-grained neural-network model is still able to reproduce its condensed-phase structural properties and phase behavior with remarkable accuracy. The limitation of the single-site model is only evident under conditions where the conformation of the molecule is highly temperature-dependent. One way to construct a neural network model that is independent of temperature would be to extract the training data from multiple temperatures and define the molecular dimensions as the average over the crystalline and isotropic phases. While the results for sexithiophene are substantially better than expected given its flexibility, improvements can be made to the model by considering a coarse-grained mapping consisting of more than one site.[64]

The coarse-grained simulation of sexithiophene demon

Figure 11: Radial distribution function (RDF) of the all-atom (solid lines) and coarse-grained (dashed lines) sexithiophene models at 1 atm and various temperatures. The RDFs have been shifted vertically for clarity. The simulated phase is given in parentheses after each temperature in the legend (I = isotropic, N = nematic, SmA = smectic A, K = crystal).

Figure 10: Orientational order parameter versus temperature for the all-atom (AA) and coarse-grained (CG) sexithiophene models at 1 atm. Typical simulation configurations are shown at each temperature for each system (AA model above the data points and CG model below), in which the molecules have been colored according to their orientation with respect to the phase director (blue = parallel, red = perpendicular). Error bars are smaller than the symbols.
strated a speed-up of 132\(\times\) compared with the all-atom simulation using the same hardware employed for the benzene simulations. This speedup is primarily due to the large reduction in the number of degrees of freedom in coarse-graining this molecule.

## Conclusions

We have applied machine learning and a recently derived systematic coarse-graining method for anisotropic particles to develop a single-site anisotropic coarse-grained potential of a molecular system. The iterative training of the neural-network potential is able to reproduce the forces, torques, and pressure of the fine-grained all-atom system. The final loss of the it

Figure 12: Angular-radial distribution function (ARDF) of the all-atom (AA) (top) and coarse-grained (CG) (bottom) sexithiophene models at 680 K and 1 atm (isotropic phase) depicted as a heat map (left) and 1D slices at constant angle (right). Face-to-face, edge-to-edge, and parallel displaced configurations occur when the angle is 0\({}^{\circ}\), while edge-to-face configurations occur at 90\({}^{\circ}\).

Figure 13: Angular-radial distribution function (ARDF) of the all-atom (AA) (top) and coarse-grained (CG) (bottom) sexithiophene models at 590 K and 1 atm (smectic phase) depicted as a heat map (left) and 1D slices at constant angle (right). Face-to-face, edge-to-edge, and parallel displaced configurations occur when the angle is 0\({}^{\circ}\), while edge-to-face configurations occur at 90\({}^{\circ}\).
erative training model was identical to the loss obtained from k-fold cross-validation. The CG model performs well for a rigid molecule like benzene but, remarkably, it also describes the phase behavior and molecular-scale structural correlations of a flexible molecule like sexithiophene with comparable accuracy, even though the aspect ratio of the molecule changes significantly over the simulated temperature range. We have demonstrated the versatility of the coarse-graining method by parameterizing models of benzene and sexithiophene at a single temperature and then studying their accuracy in capturing the structural properties of the corresponding all-atom model at different temperatures. The sexithiophene model was also used to show the ability of the model to reproduce the phase behavior of the all-atom model, with the lowest fidelity coming from the crystalline phase, where the aspect ratio of the molecule had the largest deviation from the parameterization data set. A natural extension to this work would be to generalize the method to a multi-site anisotropic coarse-grained model for flexible molecules and polymers.

## Author's contributions

MOW: Investigation, methodology, formal analysis, visualization, writing - original draft. DMH: Conceptualization, methodology, supervision, writing - review and editing.

## Acknowledgments

This work was supported by the Australian Research Council under the Discovery Projects funding scheme (DP190102100). This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government, and from the University of Adelaide's Phoenix High-Performance Computing service.

## Author declarations

### Conflict of Interest

The authors have no conflicts to disclose.

## Data availability

The Supplementary Material provides more extensive structural comparisons between the all-atom and coarse-grained simulations, more details on the methods used for the calculations, and a list of the software needed to implement and train the neural-network potential, along with links to their GitHub repositories. Simulation files and neural network scripts can be found at [https://doi.org/10.25909/21218057.v1](https://doi.org/10.25909/21218057.v1)

## References

* Butler _et al._ (2020)K. T. Butler, D. W. Davies, H. Cartwright, O. Isayev, and A. Walsh, "Machine learning for molecular and materials science," Nature **559**, 547-555 (2018).
* Moosavi _et al._ (2020)S. M. Moosavi, K. M. Jablonka, and B. Smit, "The role of machine learning in the understanding and design of materials," J. Am. Chem. Soc. **142**, 20273-20287 (2020).
* Tungasuvunakool _et al._ (2021)K. Tungasuvunakool, J. Adler, Z. Wu, T. Green, M. Zielinski, A. Zidek, A. Bridgland, A. Cowie, C. Meyer, A. Laydon, S. Velankar, G. J. Kleywegt, A. Bateman, R. Evans, A. Pritzel, M. Figurnov, O. Ronneberger, R. Bates, S. A. A. Kohl, A. Potapenko, A. J. Ballard, B. Romera-Paredes, S. Nikolov, R. Jain, E. Clancy, D. Reiman, S. Petersen, A. W. Senior, K. Kavukcuoglu, E. Birney, P. Kohli, J. Jumper, and D. Hassabis, "Highly accurate protein structure prediction for the human proteome," Nature **596**, 590-596 (2021).
* Rowe _et al._ (2020)P. Rowe, V. L. Deringer, P. Gasparotto, G. Csanyi, and A. Michaelides, "An accurate and transferable learning potential for carbon," J. Chem. Phys. **153**, 034702 (2020).
* Stocker _et al._ (2002)S. Stocker, J. Gasteiger, F. Becker, S. Gunnemann, and J. T. Margraf, "How robust are modern graph neural network potentials in long and hot molecular dynamics simulations?" Mach. Learn. Sci. Tech. **3**, 045010 (2002).
* Friederich _et al._ (2021)P. Friederich, F. Hise, J. Proppe, and A. Aspuru-Guzik, "Machine-learned potentials for next-generation matter simulations," Nat. Mater. **20**, 750-761 (2021).
* Noe _et al._ (2020)F. Noe, A. Tkatchenko, K.-R. Muller, and C. Clementi, "Machine learning for molecular simulation," Annu. Rev. Phys. Chem. **71**, 361-390 (2020).
* Guo _et al._ (2022)Z. Guo, D. Lu, Y. Yan, S. Hu, R. Liu, G. Tan, N. Sun, W. Jiang, L. Liu, Y. Chen, L. Zhang, M. Chen, H. Wang, and W. Jia, "Extending the limit of molecular dynamics with ab initio accuracy to 10 billion atoms," in _Proc. 27th ACM SIGPLAN Symp. Princ. Pract. Parallel Program._, PPoPP '22 (Association for Computing Machinery, New York, NY, USA, 2022) pp. 205-218.
* Marx and Hutter (2009)D. Marx and J. Hutter, _Ab initio molecular dynamics: basic theory and advanced methods_ (Cambridge University Press, 2009).
* Wang _et al._ (2018)H. Wang, L. Zhang, J. Han, and W. E, "DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics," Comput. Phys. Commun. **228**, 178-184 (2018).
* Unke _et al._ (2021)O. T. Unke, S. Chmiela, H. E. Sanceda, M. Gastegger, I. Poltavsky, K. T. Schutt, A. Tkatchenko, and K.-R. Muller, "Machine learning force fields," Chem. Rev. **121**, 10142-10186 (2021).
* Jin _et al._ (2022)J. Jin, A. J. Pak, A. E. Durumeric, T. D. Loose, and G. A. Voth, "Bottom-up coarse-graining: Principles and perspectives," J. Chem. Theory Comput. **18**, 5759-5791 (2022).
* Marrink and Tielleman (2013)S. J. Marrink and D. P. Tielleman, "Perspective on the Martini model," Chem. Soc. Rev. **42**, 6801-6822 (2013).
* Zhang _et al._ (2018)L. Zhang, J. Han, H. Wang, R. Car, and W. E. Weinan, "DeePCG: Constructing coarse-grained models via deep neural networks," J. Chem. Phys. **149**, 034101 (2018).
* Wang _et al._ (2019)J. Wang, S. Olsson, C. Wehmeyer, A. Perez, N. E. Charron, G. D. Fabritiis, F. Noe, and C. Clementi, "Machine learning of coarse-grained molecular dynamics force fields," ACS Cent. Sci. **5**, 755-767 (2019).
* Gay and Berne (1981)J. G. Gay and B. J. Berne, "Modification of the overlap potential to mimic a linear site-site potential," J. Chem. Phys. **74**, 3316-3319 (1981).
* Berardi _et al._ (1995)R. Berardi, C. Fava, and C. Zannoni, "A generalized Gay-Berne intermolecular potential for biaxial particles," Chem. Phys. Lett. **236**, 462-468 (1995).
* Boehm _et al._ (2019)B. J. Boehm, H. T. Nguyen, and D. M. Huang, "The interplay of interfaces, supramolecular assembly, and electronics in organic semiconductors," J. Phys.: Condens. Matter **31**, 423001 (2019).
* Behler (2011)J. Behler, "Atom-centered symmetry functions for constructing high-dimensional neural network potentials," J. Chem. Phys. **134**, 074106 (2011).
* Campos-Villalobos _et al._ (2022)G. Campos-Villalobos, G. Giunta, S. Marin-Aguilar, and M. Dijkstra, "Machine-learning effective many-body potentials for anisotropic particles using orientation-dependent symmetry functions," J. Chem. Phys. **157**, 024902 (2022).
* Nguyen and Huang (2022)H. T. L. Nguyen and D. M. Huang, "Systematic bottom-up molecular coarse-graining via force and torque matching using anisotropic particles," J. Chem. Phys. **156**, 184118 (2022).
* Noid _et al._ (2008)W. G. Noid, J.-W. Chu, G. S. Ayton, V. Krishna, S. Izvekov, G. A. Voth, A. Das, and H. C. Andersen, "The multiscale coarse-graining method. I. A rigorous bridge between atomistic and coarse-grained models," J. Chem. Phys. **128**, 244114 (2008).
* Das and Andersen (2010)A. Das and H. C. Andersen, "The multiscale coarse-graining method. V. Isothermal-isobaric ensemble," J. Chem. Phys. **132**, 164106 (2010).
* Goldstein (2002)H. Goldstein, _Classical Mechanics_ (Addison-Wesley San Francisco, 2002).
* Bartok _et al._ (2013)A. P. Bartok, R. Kondor, and G. Csanyi, "On representing chemical environments," Phys. Rev. B **87**, 184115 (2013).
* Han _et al._ (2018)J. Han, L. Zhang, R. Car, and W. E., "Deep Potential: A general representation of a many-body potential energy surface," Commun. Comput. Phys. **23**, 629-639 (2018).
* Zhou _et al._ (2019)Y. Zhou, C. Barnes, J. Lu, J. Yang, and H. Li, "On the continuity of rotation representations in neural networks," in _2019 IEEECVF Conf. Comput. Vis. Pattern Recognit. CVPR_ (IEEE, 2019) pp. 5738-5746.
* Gastegger _et al._ (2018)M. Gastegger, L. Schwiedrzik, M. Bittermann, F. Berzsenyi, and P. Marquetad, "wACSF--Weighted atom-centered symmetry functions as descriptors in machine learning potentials," J. Chem. Phys. **148**, 241709 (2018).
* Ciccotti _et al._ (2005)G. Ciccotti, R. Kapral, and E. Vanden-Eijnden, "Blue moon sampling, vectorial reaction coordinates, and unbiased constrained dynamics," ChemPhysChem **6**, 1809-1814 (2005).
* Abrams and Tuckerman (2008)J. B. Abrams and M. Tuckerman, "Efficient and direct generation of multidimensional free energy surfaces via adiabatic dynamics without coordinate transformations," J. Phys. Chem. B **112**, 15742-15757 (2008).
* Kingma and Ba (2014)D. P. Kingma and J. Ba, "Adam: A method for stochastic optimization," arXiv, 1412.6980 (2014).
* Abadi _et al._ (2016)M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, _et al._, "TensorFlow: A system for large-scale machine learning," in _12th USENIX Symp. Oper. Syst. Des. Implement. OSDI 16_ (2016) pp. 265-283.
* Chollet _et al._ (2015)Y. Chollet _et al._, "Keras," [https://github.com/fchollet/keras](https://github.com/fchollet/keras) (2015).
* Plimpton (1995)S. Plimpton, "Fast parallel algorithms for short-range molecular dynamics," J. Comput. Phys. **117**, 1-19 (1995).
* short range forces," Comput. Phys. Commun. **182**, 898-911 (2011).
* particle-particle particle-mesh," Comput. Phys. Commun. **183**, 449-459 (2012).
* Jorgensen, Maxwell, and Tirado-Rives (1996)W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, "Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids," J. Am. Chem. Soc. **118**, 11225-11236 (1996).
* Jorgensen and McDonald (1998)W. L. Jorgensen and N. A. McDonald, "Development of an all-atom force field for heterocycles. Properties of liquid pyridine and diazenes," J. Mol. Struct. THEOCHEM **424**, 145-155 (1998).
* Rizzo and Jorgensen (1999)R. C. Rizzo and W. L. Jorgensen, "OPLS all-atom model for amines: Resolution of the amine hydration problem," J. Am. Chem. Soc. **121**, 4827-4836 (1999).
* Price, Ostrovsky, and Jorgensen (2001)M. L. Price, D. Ostrovsky, and W. L. Jorgensen, "Gas-phase and liquid-state properties of esters, nitriles, and nitro compounds with the OPLS-AA force field," J. Comput. Chem. **22**, 1340-1352 (2001).
* Hockney and Eastwood (1998)R. Hockney and J. Eastwood, _Computer Simulation Using Particles_ (CRC Press, 1998).
* Ryckaert, Ciccotti, and Berendsen (1977)J.-P. Ryckaert, G. Ciccotti, and H. J. Berendsen, "Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes," J. Comput. Phys. **23**, 327-341 (1977).
* Hoover (1985)W. G. Hoover, "Canonical dynamics: Equilibrium phase-space distributions," Phys. Rev. A **31**, 1695 (1985).
* Nose (1984)S. Nose, "A molecular dynamics method for simulations in the canonical ensemble," Mol. Phys. **52**, 255-268 (1984).
* Marcot and Hanea (2021)B. G. Marcot and A. M. Hanea, "What is an optimal value of k in k-fold cross-validation in discrete Bayesian network analysis?" Comput. Statist. **36**, 2009-2031 (2021).
* Bengio and Grandvalet (2003)Y. Bengio and Y. Grandvalet, "No unbiased estimator of the variance of k-fold cross-validation," Adv. Neural Inf. Process. Syst. **16** (2003).
* Katz (1997)H. Katz, "Organic molecular solids as thin film transistor semiconductors," J. Mater. Chem. **7**, 369-376 (1997).
* Fichou (2000)D. Fichou, "Structural order in conjugated oligothiophenes and its implications on opto-electronic devices," J. Mater. Chem. **10**, 571-588 (2000).
* Dong _et al._ (2020)Y. Dong, V. C. Nikolis, F. Talnacky, Y.-C. Chin, J. Bendudin, G. Loudi, J. Kublitski, X. Zheng, S. C. B. Mannstfeld, D. Spolotre, L. Muccioli, J. Li, X. Blase, D. Beljonne, J.-S. Kim, A. A. Bakulin, G. D'Avino, J. R. Durrant, and K. Vandewal, "Orientation dependent molecular electrostatics drives efficient charge generation in homojunction organic solar cells," Nat. Commun. **11**, 4617 (2020).
* Haley and Soloway (1992)P. J. Haley and D. Soloway, "Extrapolation limitations of multilayer feed-forward neural networks," in _Proc. 1992 IJCNN Int. Jt. Int. Conf. Neural Netw._, Vol. 4 (IEEE, 1992) pp. 25-30.
* Na _et al._ (2022)G. S. Na, S. Jang, and H. Chang, "Nonlinearity encoding to improve extrapolation capabilities for unobserved physical states," Phys. Chem. Chem. Phys. **24**, 1300-1304 (2022).
* Ding _et al._ (2021)Y. Ding, A. Pervaiz, M. Carbin, and H. Hoffmann, "Generalizable and interpretable learning for configuration extrapolation," in _Proc. 29th ACM J. Met. Eur. Softw. Eng. Conf. Symp. Found. Softw. Eng._ (2021) pp. 728-740.
* Rigol, Jarvis, and Stuart (2001)J. P. Rigol, C. H. Jarvis, and N. Stuart, "Artificial neural networks as a tool for spatial interpolation," Int. J. Geogr. Inf. Sci. **15**, 323-343 (2001).
* Mueller _et al._ (2015)M. Mueller, J. Zierenberg, M. Marenz, P. Schierz, and W. Janke, "Probing the effect of density on the aggregation temperature of semi-flexible polymers in spherical confinement," Phys. Procedia **68**, 95-99 (2015).
* Seaton, Mitchell, and Landau (2006)D. Seaton, S. Mitchell, and D. Landau, "Monte Carlo simulations of a semi-flexible polymer chain: a first glance," Braz. J. Phys. **36**, 623-626 (2006).
* Sinitskiy, Saunders, and Voth (2012)A. V. Sinitskiy, M. G. Saunders, and G. A. Voth, "Optimal number of coarse-grained sites in different components of large biomolecular complexes," J. Phys. Chem. B **116**, 8363-8374 (2012).
* Falkowska _et al._ (2016)M. Falkowska, D. T. Bowron, H. G. Manyar, C. Hardacre, and T. G. Youngs, "Neutron scattering of aromatic and aliphatic liquids," ChemPhysChem **17**, 2043-2055 (2016).
* Pizzirusco _et al._ (2011)A. Pizzirusco, M. Savini, L. Muccioli, and C. Zannoni, "An atomistic simulation of the liquid-crystalline phases of sexithiophene," J. Mater. Chem. **21**, 125-133 (2011).
* Weiner _et al._ (1984)S. J. Weiner, P. A. Kollman, D. A. Case, U. C. Singh, C. Ghio, G. Alagona, S. Profeta, and P. Weiner, "A new force field for molecular mechanical simulation of nucleic acids and proteins," J. Am. Chem. Soc. **106**, 765-784 (1984).
* Weiner _et al._ (1986)S. J. Weiner, P. A. Kollman, D. T. Nguyen, and D. A. Case, "An all atom force field for simulations of proteins and nucleic acids," J. Comput. Chem. **7**, 230-252 (1986).
* Cornell _et al._ (1995)W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, and P. A. Kollman, "A second generation force field for the simulation of proteins, nucleic acids, and organic molecules," J. Am. Chem. Soc. **117**, 5179-5197 (1995).
* Tsoutrou _et al._ (2018)F. D. Tsoutrou, S. D. Peroukidis, L. D. Peristeras, and V. G. Mavrantzas, "Monte Carlo algorithm based on internal bridging moves for the atomistic simulation of thiophene oligomers and polymers," Macromolecules **51**, 8406-8423 (2018).
* Xia _et al._ (2019)W. Xia, N. K. Hansoge, W.-S. Xu, F. R. Phelan Jr,
**Supplementary Material:**

**Anisotropic molecular coarse-graining by force and torque matching with neural networks**

Marltan O. Wilson and David M. Huang

_Department of Chemistry, School of Physics, Chemistry and Earth Sciences,_

_The University of Adelaide, Adelaide, South Australia 5005, Australia_

###### Contents

* SI. Benzene network parameters
* S2. Additional benzene structural distributions
* S3. Sexithiophene network parameters
* SIV. Additional sexithiophene structural distributions
* SV. TensorFlow and LAMMPS implementation requirements
Benzene network parameters

The cut-off distance \(R_{\rm c}\) for the benzene neural network was set at 10 A for the \(G^{1}\) type symmetry function. For the \(G^{5}\) angular symmetry function, \(\lambda\) had values of -1 and 1 and \(\nu\) has values of \(2^{n}\), where \(n\in\mathbb{Z}\). Hyperparameters \(\alpha\), \(\beta\), and \(\gamma\) in the loss function were adjusted to improve the speed of convergence, but this did not usually affect the global minimum of the optimization when the number of training epochs was large.
## Sii Additional Benzene Structural Distributions
Figure S2: Angular–radial distribution function (ARDF) of the all-atom (AA) (top) and coarse-grained (CG) (bottom) benzene models at 280 K and 1 atm depicted as a heat map (left) and 1D slices at constant angle (right). Face-to-face, edge-to-edge, or parallel displaced configurations occur when the angle is 0\({}^{\circ}\), while edge-to-face configurations occur at 90\({}^{\circ}\).
Figure S4: Angular–radial distribution function (ARDF) of the all-atom (AA) (top) and coarse-grained (CG) (bottom) benzene models at 330 K and 1 atm depicted as a heat map (left) and 1D slices at constant angle (right). Face-to-face, edge-to-edge, and parallel displaced configurations occur when the angle is 0\({}^{\circ}\), while edge-to-face configurations occur at 90\({}^{\circ}\).
## Siii **Sexithiophene network parameters**

The cut-off distance \(R_{\rm c}\) for the sexithiophene neural network was set to 21 A for the \(G^{1}\) type symmetry function. For the \(G^{5}\) angular symmetry function \(\lambda\) had values of -1 and 1, \(\nu\) has values of \(2^{n}\) where \(n\in\mathbb{Z}\). Hyperparameters \(\alpha\), \(\beta\), and \(\gamma\) in the loss function were adjusted to improve the speed of convergence but did not usually affect the global minimum of the optimization when the number of training epochs was large.
## Sv Tensorflow and Lammps implementation requirements

The following list of software is needed to train and use the neural network model in coarse-grained simulations.

1. Tensorflow C API: [https://github.com/tensorflow/tensorflow/blob/master/tensorflow/c/c_api.h](https://github.com/tensorflow/tensorflow/blob/master/tensorflow/c/c_api.h)
2. Cpp Flow: [https://github.com/serizba/cppflow](https://github.com/serizba/cppflow)
3. Tensorflow Python: [https://github.com/tensorflow/tensorflow](https://github.com/tensorflow/tensorflow)
4. Keras: [https://github.com/keras-team/keras](https://github.com/keras-team/keras)
5. LAMMPS: [https://github.com/lammps/lammps](https://github.com/lammps/lammps)

The training and testing of the neural network potential was done with TensorFlow in Python using the Keras functional API. The force and torque calculations were obtained through TensorFlow's Gradient Tape feature, which provides computational derivatives with respect to the network parameters. The tanh activation function was used for all standard neural network layers
except the output layer since the tanh activation produced a smooth differentiable potential energy surface. The mean squared error was used when calculating the loss for the forces, torques, and virials. The Adam optimizer[1] was used as the gradient descent algorithm since it was able to reach the global minimum without manually updating the learning rate during training. The machine-learning potential was deployed with the TensorFlow C API and Cpp Flow wrapper. Cpp Flow allows the TensorFlow C model to be accessed directly as a force and torque calculator in a LAMMPS pair-style function.
