Systematic Bottom-Up Molecular Coarse-Graining via Force Matching using Anisotropic Particles

S3 Force Basis Functions and Least-Squares Solver S9 S3.1 Scalar Variables of Basis Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S9 S3.2 Multivariate Potential and Force Basis Functions . . . . . . . . . . . . . . . . . . . . . . . S9 S3.3 Basis Functions for Benzene and Perylene Parametrizations . . . . . . . . . . . . . . . . . S11 S3.4 Least-Squares Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . S11


S1 Derivation of Force-Matching Equation
In the main paper, the coarse-grained (CG) potential energy function U was derived from matching the fine-grained (FG) and CG equilibrium configurational probability distributions. An expression for the force on CG particle I was subsequently derived using the expression for the CG potential and the definitions of the CG configurational mapping from the FG coordinates. Here, a full derivation of the CG force equation is provided.
Enforcing consistency between FG and CG equilibrium configurational distributions requires that which can be written in terms of the FG and CG potential energy functions, u(r n ) and U (R N , Ω N ), respectively, using eqs (3) and (5) of the main paper as where and the second term in eq (S3) is independent of FG or CG coordinates. Using eqs (S3) and (S4), the force on CG particle I is Using the definition of the linear mapping from FG to CG positions of the n I FG particles in the set ζ I that are mapped onto CG particle I, and applying the chain rule, the delta function derivative in eq (S5) can be recast in terms of FG coordinates as 1 Inserting eq (S7) into eq (S5) and integrating by parts, gives The second equality in eq (S8) follows from the fact that derivatives with respect to positions of FG particles mapped on to CG particle I vanish for terms involving CG particles J = I, because no FG particle contributes to more than one CG particle. The factor operated on by the derivatives in eq (S8) includes the delta function for the mapping to CG orientation coordinates for particle I. As the orientational coordinates are defined by the inertia tensor with respect to the center-of-mass of the FG particles mapped onto CG particle I, by the chain rule the derivatives of this delta function are proportional to the derivatives applied to the inertia tensor. Using the formula for the inertia tensor of the FG particles mapped onto CG particle I, where E is the identity matrix, the sum of derivatives of the inertia tensor with respect to positions of FG particles mapped on to CG particle I is All terms in eq (S10) contain a sum of partial derivatives with respect to the FG coordinates relative to the center-of-mass of the FG particles. It is straightforward to show this sum of derivatives is zero by inserting the expression for ∆r i ≡ r i − R I using eq (S6) as follows: The derivation in eq (S11) uses the fact that the derivative ∂ri ∂rj is an identity matrix when i = j and a zero matrix otherwise, and likewise for ∂r k ∂rj . From eq (S11), the right-hand side of eq (S10) is zero, which implies that the sum of derivatives with respect to FG positions of the delta function for the mapping to CG orientational coordinates in eq (S8) is zero. This means that the only factor affected by the derivatives with respect to the FG positions that map on CG particle I in eq (S8) is the Boltzmann factor of the FG potential energy function. Thus, eq (S8) becomes ∂ri is the force on particle i in the FG system. Inserting Z(R N , Ω N ) from eq (S4) gives where · · · R N ,Ω N denotes an equilibrium ensemble average over FG configurations that are mapped to CG configuration (R N , Ω N ) and F FG,I (r n ) ≡ i∈ζ I f i (r n ) is the total force acting on the FG particles mapped onto CG particle I. Equation (S13) generalizes the force-matching condition of the MS-CG method for isotropic CG particles 1 to anisotropic particles.

S2 Derivation of Expressions for CG Mass and Inertia Tensor
The mapping from FG to CG linear and angular momenta are derived from Hamilton's equations and the definitions of the CG configuration as described in the main paper. The expressions for the CG linear and angular momenta for CG particle I in terms of FG variables are where I I I I is the inertia tensor with respect to the center-of-mass of CG particle I, I I I FG,I is the inertia tensor relative to the center-of-mass of the set of FG particles mapped onto CG particle I, and ∆r i is the position vector of the FG particles relative to their center-of-mass. The mapping M M M LI (p n , r n ) for the angular momentum depends not only on the momenta but also on the positions of the FG particles due to the dependence of the angular momentum on FG particle positions with respect to the center-of-mass. In the last step in eq (S15), the cross product of the FG position relative to the center-of-mass and the FG linear momentum has been written as a matrix product using the skew-symmetric matrix [∆r i ] defined in terms of the Cartesian components of ∆r i as The FG inertia tensor can be expressed in terms of this matrix as Enforcing consistency of the FG and CG equilibrium phase-space distributions requires that P CG,m (P N , L N )P CG,c (R N , Ω N ) = dr n dp n P FG,m (p n )P FG, where the normalized CG and FG momentum distributions are respectively. To derive eqs (S19) and (S20), we have used the identity 2 where A A A is an m × m positive-definite matrix, B and x are column vectors of length m, and det(A A A) denotes the determinant of A A A, with B = 0 in the cases of eqs (S19) and (S20).

S4
Substituting eqs (S19) and (S20) into eq (S18) and re-arranging gives where in the last step we have used the fact that no FG particle is mapped to more than one CG particle to factor the right-hand side of eq (S22) into a product of factors that each involves only the FG particles that are mapped to a specific CG particle. Integration of the factor that involves momenta of the n I FG particles mapped onto CG particle I can be carried out analytically using the mapping operators in eqs (S14) and (S15) and the Fourier representation of the Dirac delta function, which in m dimensions is 3 where x, a, and k are vectors of length m. Using the auxiliary variables Ψ and Φ in the Fourier representations of the two delta functions gives where and C C C L ≡ I I I I I I I −1 FG,I .
Evaluating the integrals over the FG particle momenta p i in eq (S24) using the Gaussian identity in S5 eq (S21) gives where in the second step we have used [∆r i ] T = −[∆r i ], in the third step we have used i∈ζ I m i ∆r i = 0 (from definition of the position of the CG particle in terms of the center-of-mass of the constituent FG particles) to eliminate the last two terms in the exponent, and in the last step we have used the definitions of I I I FG,I , c p , and C C C L in eqs (S17), (S25), and (S26), respectively. Inserting eq (S27) into eq (S24) and rearranging gives where we have applied the Gaussian identity in eq (S21). Inserting eq (S28) into eq (S22) and using det I I I I I I I −1 FG,I I I I I canceling common terms on either side of eq (S29), and pulling the factor involving the CG linear momentum P I out of the integral as it does not depend on the FG coordinates, gives For this equation to hold for all values of the CG linear momentum for all CG particles, each factor involving the linear momentum P I of a CG particle I must be equal on the two sides of the equation, i.e.
which requires the masses of the CG particles to satisfy i.e. the mass of each CG particle must equal the total mass of the FG particles that are mapped to it. This consistency condition for the CG linear momenta had previously been derived, in a more general form for more general linear mappings of FG to CG positions, in the development of the MS-CG method for isotropic CG particles. 1 Inserting eq (S32) into eq (S31) and rearranging gives where we have used the consistency condition between the FG and CG configurational probability distributions, P FG,c (r n ) and P CG,c (R N , Ω N ), respectively, given in eq (S1) to write the integral as an equilibrium average · · · R N ,Ω N over FG configurations that are mapped to CG configuration (R N , Ω N ).
Assuming that the inertia tensor I I I FG,I of the group of FG particles mapped onto a particular CG particle does not depend on the configuration of the other particles in the FG system, the average on the right-hand side of eq (S34) can be factored into independent averages for each CG particle, each of which must equal the corresponding factor on the left-hand side, which gives for I = 1, 2, . . . , N , where · · · R I ,Ω I denotes an equilibrium average over FG configurations consistent with the coordinate mapping of CG particle I. Writing eq (S35) in terms of moments of inertia, I FG,I,q and I I,q , and angular momenta, L I,q , about the principal axes, q = a, b, c, in the body-fixed frame, assuming that the factors for each principal axis are independent of one another, gives S7 for I = 1, 2, . . . , N and q = a, b, c, where I FG,I,q , I I,q , and L I,q are the FG moment of inertia, CG moment of inertia, and angular momentum about the q axis for CG particle I. This consistency condition for the CG angular momenta can be simplified to eliminate the CG moment of inertia I I,q from the right-hand side by using eq (S15) to write it in terms of CG angular velocities as for I = 1, 2, . . . , N and q = a, b, c, whereΩ I,q is the angular velocity about the q axis.
Defining δI = I FG,I,q − I FG,I,q R I ,Ω I , writing for simplicity I FG,I,q ≡ I FG,I,q R I ,Ω I , and writing the square root and exponential on the right-hand side as Taylor series, eq (S37) can be written as where (δI) m ≡ (δI) m R I ,Ω I , where m is an integer, and we have used the fact that δI R I ,Ω I = 0. For the left-and right-hand sides of eq (S38) to be consistent, they must have the same functional dependence on the angular velocityΩ I,q . This is the case if the terms on the right-hand side inside the curly braces that depend onΩ I,q are negligible for values ofΩ I,q that occur with significant probability. An order-of-magnitude estimate of such values can be obtained by assuming a Gaussian distribution oḟ Ω I,q for a rigid body with principal moment of inertia I FG,I,q about the q axis, for which the largest values ofΩ 2 I,q that occur with significant probability are on the order of k B T / I FG,I,q . Inserting this value into eq (S38), it can be seen that consistency between the left-and right-hand sides of eq (S38) is achieved if (δI) m I FG,I,q m for integers m ≥ 2. Using the definition of δI, this condition is i.e. the m-th order cumulant of the distribution of I FG,I,q is much smaller than the m-th power of its mean value. For a Gaussian distribution of I FG,I,q , for which all cumulants can be written in terms of the second moment (variance), this condition reduces to the standard deviation of I FG,I,q being smaller than its mean, i.e. (I FG,I,q − I FG,I,q R I ,Ω I ) satisfied, eq (S38) reduces to for I = 1, 2, . . . , N and q = a, b, c, from which it can be identified that for I = 1, 2, . . . , N and q = a, b, c. This means that the angular-momentum distributions of the FG and CG systems will be consistent if the fluctuations of the principal moments of inertia of the group of FG particles that is mapped onto a CG particle are small compared with their mean values, in which case the principal moments of inertia of the CG particle are approximately equal to the corresponding averages of the FG principal moments of inertia.

S3.1 Scalar Variables of Basis Functions
In the applications of the anisotropic force-matching coarse-graining (AFM-CG) algorithm presented in this paper, the basis functions for the CG potential were assumed to be functions of a set of scalar variables that are in turn functions of the relative positions and orientations of particle pairs. For two anisotropic CG particles A and B with position vectors R A and R B , respectively, and orientation defined by the matrices A A A ≡ (â 1 ,â 2 ,â 3 ) and B B B ≡ b 1 ,b 2 ,b 3 , respectively, whose columns are unit vectors aligned with the particle's principal axes in the space-fixed frame, the set of scalar variables is defined as ξ = {R, α, β, γ, φ, θ}, where the variables are the • three Euler angles {α, β, γ} that describe the rotation between the coordinate systems of particles A and B defined by their principal axes; for a CG particle that has only one unique rotation axis, such as a uniaxial ellipsoid, only β is needed to describe the relative orientation due to symmetry; • angles φ and θ between the inter-particle vector R AB and a particular principal axis of each particle; in this work the axis was taken to be the shortest principal axis of each molecule,â 3 andb 3 , respectively.
The Euler angles 4 {α, β, γ} are calculated from the components of the rotation matrix R R R that transforms the principal axes of particle A from the space-fixed frame to the body-fixed frame of particle B, which is obtained by multiplying B B B T by A A A. The angles are defined by cos α cos γ − sin α cos β sin γ − cos α sin γ − sin α cos β cos γ sin α sin β sin α cos γ + cos α cos β sin γ − sin α sin γ + cos α cos β cos γ − cos α sin β sin β sin γ sin β cos γ cos β   . (S42)

S3.2 Multivariate Potential and Force Basis Functions
A CG pair potential between CG particles A and B can be written as a linear combination of N b basis functions that are multivariate functions of the scalar variables in section S3.1, as described in eq (35) in the main paper. In this work, each basis function B i for i = 1, 2, . . . , N b was taken to have the form where µ n (R) is a cubic spline basis function of the inter-particle distance variable, and k α , k β , k γ , k φ , k θ are parameters defining the cosine basis functions of the angular variables. A cubic spline was chosen for the radial function as it ensures the smoothness and continuity of the basis force vectors. 5 For the univariate case, a cubic spline is a piece-wise continuous function of cubic polynomials with continuous derivatives up to order 2 that interpolates the data between sub-intervals of the variable space. The subintervals are connected through predetermined points called knots between which the cubic polynomials are constructed, and the boundary knots are the boundaries of the variable interval. In addition to the continuity conditions, a cubic spline can have additional constraints, namely that the function is linear beyond the boundary knots, for which the function is called a natural cubic spline. In this work, natural cubic splines were used to define N R basis functions of the inter-particle distance R in terms of N R predetermined knots (R 0 , R 1 , ..., R NR ) as 5 where with f + denoting the positive part of f (negative values are set to zero).
The force basis functions are derived as minus the gradient for the basis functions with respect to the CG position. A basis function to construct the force that particle B exerts on particle A is derived using S9 the chain rule as Here the terms involving partial derivatives of the Euler angles with respect to R A vanish because these angles do not depend on the position coordinates of the particles, and derivatives have been expressed with respect to cos φ and cos θ instead of φ and θ, respectively, to facilitate the calculations. The partial derivatives of the basis function B i (ξ) with respect to each variables in eq (S46) can be easily derived, using and in eq (S47). The remaining partial derivatives with respect to and

S3.3 Basis Functions for Benzene and Perylene Parametrizations
Because benzene and perylene are coarse-grained into uniaxial particles in this work, only terms with the variables {R, β, φ, θ} appear in the expression of the basis functions (eq (S43)) (β is the angle between the unique molecular axes of the two particles if these axes are defined asâ 3 andb 3 , respectively). Details of the basis functions used in the parametrization of the benzene and perylene models are as follows: • Benzene -N R = 23 cubic splines were used for the basis functions the inter-particle distance. The knots for these cubic splines were between 3-10 Å and chosen so that the knot density was higher at small distances and lower at large distances. Specifically, the knot R n was taken to be R n = R NR − (R NR − R 0 ) cos nπ 2(NR−1) for n = 0, . . . , N R . -The parameters k β , k φ and k θ in the cosine terms were elements of the set {0, 2, 4, 6, 8}. Only even integers were used because the CG potential has a period of π with respect to β, φ, and θ due to the symmetry of the benzene molecule.
• Perylene -N R = 39 cubic splines were used for the basis functions of the inter-particle distance. The knots for these cubic splines are between 2.5-19 Å and chosen so that the knot density is higher at small distances and lower at large distances. The same function to determine the knots was used as for benzene: Only even integers were used because the CG potential has a period of π with respect to β, φ, and θ due to the symmetry of the perylene molecule.

S3.4 Least-Squares Solver
After evaluating the force matrix as described in the main paper, the basis coefficients that optimize the CG force field can be found using a least-squares method. Various strategies for solving the force matrix have been addressed and implemented in the MS-CG method for isotropic particles. 6 We chose to apply a strategy that involves the sequential householder transformation 7 of the force matrix to an equivalent matrix of smaller size, and solving the resulting matrix equation using singular value decomposition. 8 The key idea and advantages of this strategy have been discussed in detail previously by Lu et al. 6 . The method is especially useful in the applications of the AFM-CG method in this work, in which the number of atomistic configurations (which correspond to the number of matrix rows) is large compared to the number of basis functions (number of matrix columns).

S4 Modified S-function Fit to Anisotropic CG Potential
The general expression for the modified S-function pair potential U mod S for uniaxial particles used in the CG MD simulations of benzene and perylene is described in the main paper. Here the expressions for the various terms in the potential and the expressions for the resulting forces and torques are provided. Including the pressure correction ∆U , the total pair potential is where R IJ = R I − R J is the vector that connects the particle centers,û I andû J are the unit vectors along the unique axis of each particle, and R IJ = R IJ . The form of the S-function potential used in this work is (S54) which differs from the more conventional Lennard-Jones-like form used in previous studies 9,10 in including an orientation-dependent prefactor A(R ij ,û i ,û j ) and an orientation-dependent exponent p(R IJ ,û I ,û J ), whereR IJ = R IJ /R IJ , in the repulsive term instead of constant values of 1 and 12, respectively. These functions are defined as The exponential term has the form The oscillatory term has the form (S58) The pressure-correction term for the potential cut-off distance R c has the form The various anisotropic functions in eqs (S54)-(S58) are defined as σ R IJ ,û I ,û J = σ 0 [σ 000 S 000 + σ cc2 (S 202 + S 022 ) + σ 220 S 220 + σ 222 S 222 + σ 224 S 224 + σ cc4 (S 404 + S 044 )] , (S60) R IJ ,û I ,û J = 0 [ 000 S 000 + cc2 (S 202 + S 022 ) + 220 S 220 + 222 S 222 + 224 S 224 + cc4 (S 404 + S 044 )] , S12 and S 000 = 1, S 202 = 1 2 In total, there are 36 parameters to be optimized: σ 0 , σ 000 , σ cc2 , σ 220 , σ 222 , σ 224 , σ cc4 , 0 , 000 , cc2 , , 000 exp , cc2 exp , 220 exp , 222 exp , 224 , σ exp , κ 1 , osc0 , σ osc , k 2 , κ 2 , C 2 , and b. The list of fit parameters for the benzene and perylene CG potentials are presented below in Tables S2 and S4, respectively. Note that not all of the terms in U mod S are used in both models, with unused parameters set to zero.

S4.1 Modified S-function Force and Torque Calculations
Implementation of the modified S-function potential for MD simulations require calculations of the forces and torques on the CG particles. The expressions for the forces and torques for a pair of uniaxial particles are 9 and The derivative with respect to the inter-particle separation R IJ is where and The derivatives with respect to the orientational variables f i , where i = 0, 1, or 2, are and (S87) S14 S5 Benzene Simulations

S5.1 Distributions of Principal Moments of Inertia
The distributions of the principal moments of inertia from the fine-grained (FG) atomistic MD simulation of benzene at 300 K and 1 atm are shown in Figure S1 and the mean values are given in Table S1. The fluctuations of the principal moments are small compared to the mean values, with the standard deviation less than 2% of the mean value of each principal moment. Fig. S1: Distributions of the principal moments of inertia of FG benzene at 300 K and 1 atm.

S5.2 AFM-CG Model Parameters
The parameters used in the AFM-CG MD simulations of benzene are presented in Table S2.

S16
To obtain the optimized parameters for the modified S-function potential (eq (S53)) in Table S2 used in the AFM-CG MD simulations, a least-squares fit of the S-function term (eq (S54)) in this potential to the basis-expansion potential was first carried out for 300 000 dimer configurations with a separation distance up to 9 Å randomly sampled from the equilibrated FG simulation used for parametrization. This gave a reasonable fit, but the depth of the potential well for the edge-to-edge configuration was overestimated, which led to an overestimate of the probability of this configuration in CG simulations. To overcome this issue, the parameters in the orientation-dependent prefactor A(R ij ,û i ,û j ) and orientation-dependent exponent p(R IJ ,û I ,û J ) were optimized again by a least-squares fit with all other parameters in eq (S54) fixed using an augmented data set that included an additional 1000 edge-to-edge configurations with separation distances uniformly distributed between 6-9 Å. Then, to fit the significant positive hump in the potential for configurations close to face-to-face, the oscillatory term (eq (S58)) in eq (S53) was added to the fit potential, the set of dimer configurations was augmented with 2000 additional face-to-face and parallel-displaced (PD) configurations (PD-10, PD-30, PD-50 and PD-70 as denoted in Figure S2) with separation distances uniformly distributed up to 9 Å, and a least-squares fit was used to optimize the parameters in the oscillatory term with the previously optimized S-function parameters kept fixed. The lower limit of the separation distance of the dimer configurations added in this step was chosen such that the potential was less than 5 kcal mol −1 . Basis expansion and fitted modified S-function potential curves for various benzene-benzene dimer relative orientations are shown in Figure S2.  S2: Benzene AFM-CG pair potential (without pressure correction) for various dimer relative orientations using the basis set expansion (solid lines) and fitted modified S-function potential (dashed lines). The scalar variables β, φ, θ (defined in section S3) that describe the relative orientation for theses configurations are: face-to-face (0°,0°,0°), edge-to-face (90°,0,90°), edge-to-edge (0°,90°,90°), PD-x (parallel-displaced) (0°,x°,x°), Cross-x (x°,90°,90°), SP-x (π-stacking) (x°,0°,x°).
The pressure-corrected potential does not differ substantially the modified S-function potential fit to the uncorrected AFM-CG potential, as shown in Figure S3.

S5.3 NVT Simulations
Properties of the FG benzene model from NPT simulations at 1 atm are compared in this section with those from NVT simulations of CG benzene models in which the density was set to the average density of the FG simulation at the corresponding temperature.

S5.3.1 Radial Distribution Function (RDF)
RDFs of the FG, AFM-CG, and Bowen CG 9 models are compared at two different temperatures in Figure S4. ARDFs of the FG, AFM-CG, and Bowen CG 9 benzene models at 300 K and 330 K are compared in Figure S5 and Figure S6, respectively.

S5.4.1 Radial Distribution Function (RDF)
RDFs from NPT simulations at 1 atm of the FG, AFM-CG, and Bowen CG 9 models are compared at two different simulation temperatures in Figure S7. The mean squared displacement (MSD) was calculated as described in the main paper and used to calculate the translational diffusion coefficients of the FG and AFM-CG benzene models. The orientational time correlation function (OTCF) was calculated as described in the main paper and used to calculate the rotational diffusion coefficients with respect to the out-of-plane and in-plane molecular axes of the FG and AFM-CG benzene models.

S6.1 Distributions of Principal Moments of Inertia
Distributions of the moments of inertia about three principal axes for FG atomistic perylene are shown in Figure S12. The fluctuations of the principal moments are small compared to the mean values, with the standard deviation less than 2% of the mean value of each principal moment.

S6.2 AFM-CG Model Parameters
The parameters used in the AFM-CG MD simulations of perylene are presented in Table S4. To obtain the optimized parameters for the modified S-function potential (eq (S53)) used in the AFM-CG MD simulations, a least-squares fit of the sum of the S-function term (eq (S54)) and exponential term (eq (S57)) to the basis-expansion potential was first carried out for 500 000 dimer configurations up to a separation distance of 14 Å randomly sampled from the equilibrated FG simulation used for parametrization. Then, similarly to benzene, to fit the significant positive hump in the potential for configurations close to face-to-face, the oscillatory term (eq (S58)) in eq (S53) was added to the fit potential, the set of dimer configurations was augmented with 4200 additional face-to-face, various parallel-displaced (PD) (PD-10, PD-30, PD-50 and PD-70 as denoted in Figure S13), and pi-stacking (SP) (SP-10 and SP-30 as denoted in Figure S13) configurations with separation distances uniformly distributed up to 19 Å, and a least-squares fit was used to optimize the parameters in the oscillatory term with the previously S26 optimized parameters kept fixed. The lower limit of the separation distance of the dimer configurations added in this step was chosen such that the potential was less than 5 kcal mol −1 . Basis expansion and fitted modified S-function potential curves for various perylene-perylene dimer relative orientations are shown in Figure S13. Fig. S13: Perylene AFM-CG pair potential (without pressure correction) for various dimer relative orientations using the basis expansion (solid lines) and fitted modified S-function (dashed lines). The scalar variables β, φ, θ (defined in section S3) that describe the relative orientation for these configurations are face-to-face (0°,0°,0°), edge-to-face (90°,0,90°), edge-to-edge (0°,90°,90°), PD-x (parallel-displaced) (0°,x°,x°), Cross-x (x°,90°,90°), SP-x (π-stacking) (x°,0°,x°), and ST-x (x°,x°,90°).

S27
The pressure-corrected potential does not differ substantially to the modified S-function potential fit to the uncorrected AFM-CG potential, as shown in Figure S14.

S6.3 NVT Simulations
Properties of the FG perylene model from NPT simulations at 1 atm are compared in this section with those from NVT simulations of CG perylene models in which the density was set to the average density of the FG simulation at the corresponding temperature.

S6.3.1 Radial Distribution Function (RDF)
RDFs of the FG, AFM-CG, Babadi CG, 11 and Berardi CG 12 perylyene models are compared at three different temperatures in Figure S15. ARDFs from the NVT simulations at 570 K of the FG and AFM-CG perylene models are plotted in Figure S16, while those for the Babadi CG 11 and Berardi CG 12 perylene models are plotted in Figure S17. Analogous plots at 670 K are shown in Figures S18 and S19, respectively.    ARDFs from NPT simulations at 570 K and 1 atm of the FG and AFM-CG models are plotted in Figure S20, while those of the Babadi CG 11 and Berardi CG 12 perylene models at the same conditions are plotted in Figure S21. Analogous plots at 670 K are shown in Figures S22 and S23, respectively.