AA, CG, TMA, MMT, TPS, $T_{g}$, $R_{g}$, MSD

# MARTINI Coarse-grained Force Field for Thermoplastic Starch Nanocomposites

Ankit Patidar  Gaurav Goel  Department of Chemical Engineering, Indian
Institute of Technology Delhi, New Delhi, India<EMAIL_ADDRESS>
###### Abstract

Thermoplastic starch (TPS) is an excellent film-forming material, and the
addition of fillers such as tetramethylammonium-montmorillonite (TMA-MMT) clay
has significantly expanded its use in packaging applications. We first used
all-atom (AA) simulation to predict several macroscopic (Young’s modulus,
glass transition temperature, density) and microscopic (conformation along 1–4
and 1–6 glycosidic linkages, composite morphology) properties of TPS melt and
TPS–TMA-MMT composite. The interplay of polymer-surface (weakly repulsive),
plasticizer-surface (attractive), and polymer-plasticizer (weakly attractive)
interactions lead to conformational and dynamics properties distinct from
those in systems with either attractive or repulsive polymer-surface
interactions. A subset of AA properties was used to parameterize the MARTINI
coarse-grained (CG) force field (FF) for the melt and composite systems.
Specifically, we determined the missing bonded parameters of amylose and
amylopectin and rationalized the bead types for 1–4 and 1–6 linked $\alpha$-D
glucose using two-body excess entropy, density, and bond and angle
distributions in AA TPS melt. The MARTINI CG model for TPS was combined with
an existing parameter set for TMA-MMT. The liquid-liquid partitioning-based
MARTINI-2 FF shows freezing and compaction of polymer chains near the sheet
surface, further accentuated by lowering of dispersive interactions between
pairs of high covalent-coordination ring units of TPS polymers and MMT sheet.
A rescaling of the effective dispersive component of TPS-MMT cross-
interactions was used to optimize the FF for the composite system, with
structural (chain size distribution), thermodynamic (chain conformational
entropy, density), and dynamic (self-diffusion coefficient) properties
obtained from long AA simulations forming the constraints for optimization.
The obtained CG FF parameters provided excellent estimates for several other
properties of the melt and composite systems not used in parameter estimation,
thus establishing the robustness of the developed model.

###### keywords:

American Chemical Society, LaTeX

##  1 Introduction

Thermoplastic starch (TPS) is obtained by gelatinization of native starch, a
polysaccharide consisting of amylose and amylopectin, by addition of
plasticizers, such as glycerol, sorbitol, water, etc. 1. Amylose is a linear
polymer of 1-4 linked $\alpha$-D glucose and amylopectin has additional 1-6
branching on the 1-4 backbone (Figure S4). Low cost, excellent film-forming
properties (low glass transition temperature, easy melt processability,
acceptable elongation at break) 2, 3, low oxygen permeability, and
biodegradability 4 make TPS an attractive option for food packaging, drug
delivery, paper additive, etc. 5, 6, 7, 8, 9. However, its poor mechanical
properties (low Young’s modulus and tensile strength, limited elongation at
break) 10, 11, 12 and low moisture resistance 13 necessitate the addition of
components with complementary properties, such as addition of nanofillers 14,
15, 16, blending with other polymers 17, 18. Several studies have shown
improvement in tensile strength and moisture barrier of TPS films on adding
surface polarity-matched phyllosilicate clays, such as organically-modified
montmorillonite (MMT) 19, 15, 20, 10, 21, 22, 23. In general, wide
availability, low cost, and surface tunability make MMT a particularly
suitable nanofiller for modulating properties of polymeric materials 14, 24,
25. For example, blending with polylactic acid (PLA) was shown to nominally
reduced water sorption of a TPS film, but microphase separation of two polymer
phases led to poorer mechanical properties 26. Here, addition of
$1\text{\,}\mathrm{\char 37\relax}$ Na-MMT) resulted in a lower number-average
domain diameter of PLA (from $21.3\text{\,}\mathrm{\SIUnitSymbolMicro m}$ to
$13.45\text{\,}\mathrm{\SIUnitSymbolMicro m}$) and an improved tensile
strength (from $5.63\text{\,}\mathrm{MPa}$ to $7.95\text{\,}\mathrm{MPa}$).
The properties of composite materials depend on various factors, including
polymer chemistry, polymer-filler interactions, and filler aspect ratio, which
consequently determine filler particle dispersion and composite’s morphology
27, 23, 14, 16. The identification of specific nanoparticles and other
topological parameters for improving properties of interest and optimization
in the composition-property space requires a time- and resource-intensive
exploration 27, 23, 14, 16, made more challenging for starch-based composites
due to a large diversity in the base components (TPS is a multi-component
material) and chemical heterogeneity present in each component (e.g., primary
structure diversity of amylopectin). Molecular dynamics (MD) simulations
provide a viable first step to characterize structure-property relationships
of multi-component systems and obtain narrow bounds on the composition space
for a target property, thus complementing efforts to develop new polymer
composites 28, 29, 30, 31. MD simulations using all-atom (AA) representation
of polysaccharides 32, 33, 34, 35, 36 have been shown to accurately predict a
variety of properties of some simple systems, for example, glass transition
temperature and Young’s modulus of TPS melt 37, 38, chain size distribution of
amylose in water 39, and crystalline packing of cellulose fibers 40. However,
slow relaxation dynamics and disparate length scales of the components
preclude the applicability of AA simulations to several systems of fundamental
and technological interest. For example, the characteristic relaxation time of
a small 20-unit amylose
($3260.86\text{\,}\mathrm{g}\text{\,}{\mathrm{mol}}^{-1}$) in its melt was
reported to be $404.7\text{\,}\mathrm{ns}$ at $700\text{\,}\mathrm{K}$ (much
above typical processing temperature) 41, eight orders of magnitude higher
than the AA time step even for this limiting case. The presence of high
molecular weight polymer chains (radius of gyration
$\sim$50\text{\,}\mathrm{nm}$$) and highly anisotropic fillers such as
$\sim$200\text{\,}\mathrm{nm}$$ diameter bentonite clay) bentonite clay in TPS
film make property determination using AA simulations computationally
intractable. The coarse-grained (CG) force fields, wherein some of the AA
degrees of freedom are integrated out, can be $\sim$ 100-fold faster than AA
simulations even at the lower-end of coarse-graining. CG approaches employed
to study polymer composites at some level of relevant length- and timescales
can be broadly divided into two categories based on the treatment of
interaction between various components. In one case, the enthalpic
contribution is modeled by a minimal, single-parameter effective interaction
between large units (a Kuhn monomer of $\sim$ $1\text{\,}\mathrm{nm}$, a
nanoparticle, etc.). The effective interaction parameters are then used in
mesoscale dissipative particle dynamics (DPD) simulations 42, 43, 44, 45, 46,
47, 48, 49, 50 or field-based / integral-equation theories 51, 52, 53, 54, 55,
56, 57, 58, 59, 60, 61, 62, 63, to characterize structure and thermodynamics
of diverse polymeric systems. In a multiscale DPD scheme, integral equation
theory was used to obtain an effective softcore interaction between polymer
chains, enabling use of a very large time step ($\text{\,}\mathrm{ps}$ \-
$\text{\,}\mathrm{ns}$) 64, 65, and its application to polyethylene melt 66,
67 and blends 68 showed a fairly good comparison with structural (radius of
gyration and radial distribution function) and dynamic properties (diffusion
coefficient of chain center of mass) obtained in united atom (UA) simulations.
These methods, based on a simplified model for enthalpy, have provided
qualitative and, in several cases, quantitatively accurate predictions for the
phase diagram and the equilibrium properties of large-scale polymer systems
49, 69, 70, 71, 72. The interaction parameter is dependent on the
thermodynamic state variables (temperature, pressure, composition) and needs
to be obtained semi-empirically from molecular simulations, liquid-state
theories, or experiments at a sufficiently large number of state points,
wherein an increase in system heterogeneity significantly expands the phase
space for interaction parameter determination.

Chemical-specific CG models provide an alternative with possible state-point
transferability and good accuracy across highly diverse systems. Herein, 4-10
atoms are grouped in a single bead, with effective pair interaction obtained
from a combination of AA simulations and experimental data. This mapping from
AA to CG can be obtained by using structure distribution functions as
constraints, as done in iterative Boltzmann inversion (IBI) 73, 74 and force
matching 75 or by using thermodynamic properties as constraints, such as
equilibrium partitioning between polar and apolar solvents used for MARTINI CG
force field 76. In structure-based methods, parameters need to be optimized
using data from multiple state points 77, and can still exhibit environmental
dependence. The MARTINI force field has been demonstrated to mitigate these
limitations. For instance, the MARTINI CG parameters for polypropylene were
shown to accurately reproduce (values obtained in AA simulation) the radius of
gyration in various chemical environments (good, bad, and theta solvent
conditions) as well as the partitioning of polyethylene-polypropylene blends
near phospholipid membranes 78. However, a reduction in self- (saccharides
molecules) 79, 80 or cross-interaction (graphene-small organic molecules) 81,
82, 83 was required in some cases, highlighting a need for further parameter
optimization.

Given the wide interest in phyllosilicate-based composites, we have developed
MARTINI CG parameters for a TPS–[tetra methyl ammonium(TMA)]-MMT system in
this study. The CG parameters for some individual TPS components (amylose 84,
sorbitol 85, water 86) and TMA-MMT 87 have been shown to reproduce several
properties determined from AA simulations. However, to the best of our
knowledge, MARTINI CG FF for a TPS melt or its composites is not available.
Here, we have re-parameterized and developed missing bonded parameters for the
CG model of two polymeric components (amylose and amylopectin) by using data
from AA simulation of TPS melt, rationalized the bead type assignment for 1-4
and 1-6 linked $\alpha$-D glucose, and determined the cross-interaction levels
for all TPS and TMA-MMT beads. The developed parameter set was extensively
validated by comparing structural (chain size distribution), thermodynamic
(chain conformational entropy, density, excess entropy), and dynamic (self-
diffusion coefficient) properties against sufficiently long AA simulations.
Finally, predictions for radial distribution function (RDF) and two-body
excess entropy in the composite system, two properties not used in CG
parameterization, were used to test the robustness of developed CG parameters.

##  2 Methods

###  2.1 System and Simulation Details

####  2.1.1 AA MD Simulations

The initial structures and CHARMM36 88, 89 parameters of amylopectin (thirty-
six monomers), amylose (eighteen monomers), and sorbitol were obtained using
CHARMM-GUI 90, 89, 88. We used a TIP3P water model and the INTERFACE force
field 91 for the TMA-MMT sheet. A low density
($\rho=$117.62\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$$) TPS
configuration (55327 atoms) consisting of (weight percent)
$50\text{\,}\mathrm{\char 37\relax}$ amylopectin, $20\text{\,}\mathrm{\char
37\relax}$ amylose, $29\text{\,}\mathrm{\char 37\relax}$ sorbitol, and
$1\text{\,}\mathrm{\char 37\relax}$ water was prepared in Packmol 92. A high-
density system was obtained by rescaling the box in multiple steps followed by
temperature-pressure ($T$-$P$) annealing cycles (details in Section 5.1).
Densities of $1227.75\text{\,}\mathrm{a}$nd
$1227.62\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ at
$613\text{\,}\mathrm{K}$, $1\text{\,}\mathrm{bar}$ were obtained after the
first two $T$-$P$ cycles, respectively, indicating the adequacy of only one
cycle. The obtained value compares well with a TPS melt density of
$1369.80\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ from AA simulation at
$550\text{\,}\mathrm{K}$, $1\text{\,}\mathrm{bar}$ 37 and
$1442\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ at
$300\text{\,}\mathrm{K}$, $1\text{\,}\mathrm{bar}$ from experiments 93. For
composite simulations, a
$$10.32\text{\,}\mathrm{nm}$\times$10.71\text{\,}\mathrm{nm}$\times$1.2\text{\,}\mathrm{nm}$$
TMA-MMT sheet, connected across $x$-$y$ periodic boundaries to model an
infinitely large sheet, was placed between two equilibrated TPS melt boxes
(Figure S11 (a)). An $82.5\text{\,}\mathrm{\char 37\relax}$ TPS and
$17.5\text{\,}\mathrm{\char 37\relax}$ TMA-MMT system with dimensions
$10.32\text{\,}\mathrm{nm}$$\times$$10.71\text{\,}\mathrm{nm}$$\times$$11.40\text{\,}\mathrm{nm}$
was obtained after one $T$-$P$ annealing cycle. All melt and composite
properties were calculated from a $1100\text{\,}\mathrm{ns}$ MD simulation
started from the configuration at the end of the first $T$-$P$ annealing
cycle. MD simulations were performed using GROMACS 2021.4 94, with a time step
of $1\text{\,}\mathrm{fs}$, V-rescale thermostat
($\tau$=$1\text{\,}\mathrm{ps}$) 95, and either the Berendsen barostat
($\tau$=$2\text{\,}\mathrm{ps}$, first $100\text{\,}\mathrm{ns}$) 96 or the
Parrinello-Rahman barostat ($\tau$=$5\text{\,}\mathrm{ps}$) 97.

####  2.1.2 CG MD Simulations

The initial CG structures for TPS melt and TPS–TMA-MMT composite were obtained
from the last frame of the corresponding $1100\text{\,}\mathrm{ns}$ AA
simulation by placing CG beads at the center of mass positions (Equation S3)
of constituent atoms (mapping scheme in Figure 1). Force field parameters at
the level of MARTINI v2.276 were either directly used, if available and found
suitable, or were obtained using reference AA simulations. For all CG
simulations, a time step of $5\text{\,}\mathrm{fs}$ was used, the system was
coupled to the V-rescale thermostat ($\tau$=$1\text{\,}\mathrm{ps}$) and
either the Berendsen barostat ($\tau$=$8\text{\,}\mathrm{ps}$, first
$100\text{\,}\mathrm{ns}$) or the Parrinello-Rahman barostat
($\tau$=$12\text{\,}\mathrm{ps}$). The CG potential is given as per Equation
1.

$U=4\epsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}-\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right]+\frac{q_{i}q_{j}}{4{\pi}\epsilon_{0}\epsilon_{rel}{r_{ij}}}\\\
+\frac{1}{2}K_{l}(l_{ij}-l_{0})^{2}+\frac{1}{2}K_{\theta}(cos\theta_{ijk}-cos\theta_{0})^{2}\\\
+K_{\phi}[1+(cos\phi_{ijkl}-cos\phi_{0})],$ (1)

where $\sigma_{ij}$, $\epsilon_{ij}$, $q$, and $\epsilon_{rel}$ are bead
diameter, interaction strength, bead charge, and relative dielectric constant,
respectively. The force constants for bond ($l_{ij}$), angle ($\theta_{ijk}$),
and torsion ($\phi_{ijkl}$) potentials are $K_{l}$, $K_{\theta}$, and
$K_{\phi}$, respectively. $\sigma_{ij}=$0.47\text{\,}\mathrm{nm}$$ for the
regular beads (four heavy atoms) and $0.43\text{\,}\mathrm{nm}$ for the small
beads (s-type, 2-3 heavy atoms). The choice of $\epsilon_{ij}$ is linked to
the functional group(s) in the CG bead, with primary categorization defined as
polar (P), neutral (N), apolar (C), or charged (Q), and level within a
category assigned by an integer in range 0–6 (e.g., P4 is more polar than P1,
C4 is more apolar than C1). We used a P4c bead for three diol units of
sorbitol (see Figure 1), a CG model shown to be in qualitative agreement with
experiments on the kinetics of aggregation of amyloid- and elastin-like
peptides in presence of sorbitol 85. The refined polarizable model was used
for water since the standard MARTINI water model (one P4 bead) exhibits
unphysical freezing near surfaces 98, 86. The parameters for the TMA-MMT sheet
were taken from Khan and Goel 87 who found several properties of hydrophobic
polymers–TMA-MMT composites to be in a good agreement with AA simulations
(Supplementary Information Figure S10). The bonded parameters for 1–4 linked
$\alpha$-D glucose were either taken from Lopez et al.’s study 84 on short
amylose chains in polar (water) and non-polar (nonane) solvents (labeled CG1
model) or obtained as best-fit models to the corresponding pair, triplet, and
quadruplet distributions in AA simulations of TPS melt (labeled CG2 model).
The bonded parameters for 1-6 linked $\alpha$-D glucose were similarly
obtained from the AA simulations of TPS melt (for both CG1 and CG2).

Figure 1: Mapping scheme for TPS components. (a) A trimeric 1–4 bonded
$\alpha$-D glucose oligomer (represents terminal and middle monomer of
amylose), (b) A trimeric 1–4 or 1–6 bonded $\alpha$-D glucose oligomer
(represents the branch point in amylopectin), and (c) sorbitol. Red, white,
and cyan represent oxygen, hydrogen, and carbon atoms, respectively, with
atoms in a single CG bead grouped within shaded circles. The bead number and
type are indicated by black and red color labels, respectively.

###  2.2 Property Calculation

In the composite system, the near- and far-sheet regions were classified based
on density profiles of TPS components along the normal to the TMA-MMT sheet
($z$-axis). For polymer chains, the $z$-component of the distance of the
center of mass (COM) of a given trimeric segment ($\sim$ persistence length
($0.86\text{\,}\mathrm{nm}$) of amylose) from the sheet was used to classify
it as a near or a far segment, with a chain labeled as a near/far-region chain
if $>$70\text{\,}\mathrm{\char 37\relax}$$ segments were in the near/far
region. Unless otherwise indicated, all properties are calculated at
$613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$ using the last
$900\text{\,}\mathrm{ns}$ trajectory (a $200\text{\,}\mathrm{ns}$ convergence
time for some properties in total $1100\text{\,}\mathrm{ns}$ trajectory).

####  2.2.1 Structural and Dynamic Properties

Density, radial distribution function ($RDF$), persistence length (${l_{p}}$),
and radius of gyration (${R_{g}}$) of TPS melt and composite systems were
estimated using GROMACS utilities. COMs of atoms making the CG beads were used
for $RDF$ calculation in the AA simulations to allow comparison with the CG
simulations. An exclusion of up to four CG-bonded neighbors or up to 16 AA-
bonded neighbors, equivalent to an average COM separation of
$1.29\text{\,}\mathrm{nm}$, was used for $RDF$ calculation. For the composite
system, a 2D $RDF$ was calculated in $0.37\text{\,}\mathrm{nm}$ slices
parallel (in the $x-y$ plane) to the TMA-MMT sheet. The data is reported as
average over four slices in the near region and five slices in the far region.

The diffusion coefficients of sorbitol and trimeric polymer chain segments
(amylose and amylopectin) were calculated from the corresponding temporal
profile of mean-square-displacement (${MSD}(t)$) using the Einstein’s
diffusion equation (Equation S2). The AA and CG data were compared using the
normalized diffusion coefficient, $D_{\mathrm{n}}$, defined as the ratio of
lateral diffusion coefficient in the composite (in the $x-y$ plane) to the
melt diffusion coefficient.

####  2.2.2 Thermodynamic Properties

The TPS melt was cooled from $613\text{\,}\mathrm{K}$ to
$150\text{\,}\mathrm{K}$ in steps of $25\text{\,}\mathrm{K}$, wherein a
$4\text{\,}\mathrm{ns}$ NPT simulation was performed at each step, with last
$500\text{\,}\mathrm{ps}$ trajectory at each step used for density
calculation. The intersection of linear fits to the temperature-density data
in the rubbery (high temperature) and glassy (low temperature) regions was
used as an estimate for $T_{\mathrm{g,sim}}$. The Williams-Landel-Ferry (WLF)
equation 99 (Equation 2) was used to obtain an estimate for the glass
transition temperature at a typical experimental cooling rate of
$10\text{\,}\mathrm{K}\text{\,}{\mathrm{min}}^{-1}$ ($T_{\mathrm{g,exp}}$)
from the simulation estimate $T_{\mathrm{g,sim}}$ obtained at
$3.75\text{\times}{10}^{11}\text{\,}\mathrm{K}\text{\,}{\mathrm{min}}^{-1}$.

$\mathrm{log_{10}}\frac{\tau_{g,\mathrm{sim}}}{\tau_{g,\mathrm{exp}}}=\frac{\mathrm{-C_{1}*(T}_{g,\mathrm{sim}}-\mathrm{T}_{g,\mathrm{exp}})}{\mathrm{C_{2}+T}_{g,\mathrm{sim}}-\mathrm{T}_{g,\mathrm{exp}}},$
(2)

where $C_{1}=17.44$ and $C_{2}=$51.6\text{\,}\mathrm{K}$$ (constant for
several polymers), $\tau_{\mathrm{g,sim}}$ and $\tau_{\mathrm{g,exp}}$ are the
inverse of cooling rates in simulation and experiment, respectively.
The pairwise ($S_{2,\alpha\beta}$) and total ($S_{\mathrm{2}}$) two-body
excess entropy was calculated using the corresponding $RDFs$
($g_{\alpha\beta}$) as per Equations 3 and 4 (integral calculated up to
$2.5\text{\,}\mathrm{nm}$), respectively 100.

$\mathrm{S_{2,\alpha\beta}}=-2\pi\rho\int_{0}^{\infty}g_{\alpha\beta}(r)\mathrm{ln}\,g_{\alpha\beta}(r)-\\\
[g_{\alpha\beta}(r)-1]r^{2}dr$ (3)
$\mathrm{S_{2}}=\sum_{\alpha,\beta}{x_{\alpha}x_{\beta}}\mathrm{S_{\alpha\beta}},$
(4)

where $x_{\alpha}$ and $x_{\beta}$ are the mole fractions of $\alpha$ and
$\beta$ components, respectively. Equation 3 was suitably modified to
integrate 2D RDFs for the composite system over cylindrical disks (Equation
S1).

The configurational entropy ($S_{\mathrm{c}}$) of amylose and amylopectin
chains was calculated from the covariance matrix of each chain (after removing
rigid-body translational and rotational motion) using the Schlitter’s equation
101, implemented with GROMACS utility anaeig.

####  2.2.3 Mechanical properties

Young’s modulus of TPS melt and TMA-MMT sheet was calculated from AA
simulations at $300\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$. A
uniaxial strain along $x$-axis was added in steps of
$0.1\text{\,}\mathrm{\char 37\relax}$ of the initial box length, with a
$2\text{\,}\mathrm{ns}$ $NPT$ simulation at each step (anisotropic pressure
coupling for $y$ and $z$). The stress at every step was calculated from the
${Pxx}$ Virial component averaged over the last $300\text{\,}\mathrm{ps}$ of
the trajectory, and the slope of the stress-strain curve in the linear regime
was used to estimate the Young’s modulus. Maximum strain value
($\sigma_{\mathrm{max}}$) marks the onset of the plateau region in the stress-
strain profile.

##  3 Results and discussion

###  3.1 Validation of TPS AA parameters

Figure 2: Polymer chain conformation at inter-ring links in a TPS melt.
$\phi$-$\psi$ torsion angle distribution in TPS melt for (a) 1–4 and (b) 1–6
glycosidic linkages, shown as the potential of mean force (PMF) profile (color
scale in $\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$) obtained from the
probability distributions averaged over $900\text{\,}\mathrm{ns}$ AA
simulation at $613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$.

The local chain conformation in TPS melt was characterized using the torsion
angles for 1–4 and 1–6 glycosidic linkages, viz.
$\phi=\angle{\mathrm{O_{5}-C_{4}-O_{1}-C_{1}}}$ and
$\psi=\angle{\mathrm{C_{1}-O_{1}-C_{4}-C_{3}}}$ for 1–4, and
$\phi=\angle{\mathrm{O_{5}-C_{1}-O_{1}-C_{6}}}$ and
$\psi=\angle{\mathrm{C_{1}-O_{1}-C_{6}-C_{5}}}$ for 1–6 (Figure S5). Figure 2
(a) shows syn- and anti-states as two dominant conformers for the 1–4 linkage,
with peak at ($\phi$,$\psi$) equal to
($90\text{\,}\mathrm{\SIUnitSymbolDegree}$,$100\text{\,}\mathrm{\SIUnitSymbolDegree}$)
and
($100\text{\,}\mathrm{\SIUnitSymbolDegree}$,$-50\text{\,}\mathrm{\SIUnitSymbolDegree}$),
respectively. This agrees well with the AA simulations of an amylose oligomer
in water using the same force field ($syn$:
($90\text{\,}\mathrm{\SIUnitSymbolDegree}$,$100\text{\,}\mathrm{\SIUnitSymbolDegree}$)
and $anti$:
($80\text{\,}\mathrm{\SIUnitSymbolDegree}$,$-50\text{\,}\mathrm{\SIUnitSymbolDegree}$))
39 and XRD spectroscopy of a hydrated cycloamylose ($syn$:
($103.6\text{\,}\mathrm{\SIUnitSymbolDegree}$,$115.1\text{\,}\mathrm{\SIUnitSymbolDegree}$)
and $anti$:
($88.1\text{\,}\mathrm{\SIUnitSymbolDegree}$,$-48.4\text{\,}\mathrm{\SIUnitSymbolDegree}$))
102. The conformational space of 1–6 linkage has a single peak at
($70\text{\,}\mathrm{\SIUnitSymbolDegree}$,$180\text{\,}\mathrm{\SIUnitSymbolDegree}$)
(Figure 2(b)), in close agreement with replica exchange simulation and NMR
spectroscopy on $\alpha$-1–6 glucopyranoside dimer
(($\phi$,$\psi=$)($70\text{\,}\mathrm{\SIUnitSymbolDegree}$,$180\text{\,}\mathrm{\SIUnitSymbolDegree}$))
103.

Figure 3: AA and CG representation of TPS melt. A 18-mer amylose chain in (a)
AA (oxygen, hydrogen, and carbon atoms in red, white, and cyan, respectively)
and (b) CG (P1, SP1, P4c, and Na beads in blue, pink, yellow, and ice blue,
respectively) representations. (c) and (d) TPS melt simulation box obtained at
end of $1100\text{\,}\mathrm{ns}$ AA simulation at $613\text{\,}\mathrm{K}$
and $1\text{\,}\mathrm{bar}$ in AA and CG, respectively. Polymer, sorbitol,
and water beads in red, blue, and yellow, respectively.

Table 1: A comparison of various properties of TPS melt obtained from AA
simulations with previous simulations and experimental measurements.

Property | Present | Previous | Previous
---|---|---|---
| Simulations | Experiments | simulations
1–4 $\psi$$\text{\,}\mathrm{\SIUnitSymbolDegree}$ |  $100\text{\,}\mathrm{\SIUnitSymbolDegree}$,$50\text{\,}\mathrm{\SIUnitSymbolDegree}$111The $\phi,\psi$ torsion angle distribution from the AA simulation of TPS melt at $613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$ |  $115.1\text{\,}\mathrm{,}$$48.4\text{\,}\mathrm{\SIUnitSymbolDegree}$222XRD spectroscopy of hydrated cycloamylose 102 |  $90\text{\,}\mathrm{\SIUnitSymbolDegree}$, $50\text{\,}\mathrm{\SIUnitSymbolDegree}$333AA simulation of 9-mer Amylose in water 39
1–4 $\phi$$\text{\,}\mathrm{\SIUnitSymbolDegree}$ |  $90\text{\,}\mathrm{\SIUnitSymbolDegree}$,$100\text{\,}\mathrm{\SIUnitSymbolDegree}$111The $\phi,\psi$ torsion angle distribution from the AA simulation of TPS melt at $613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$ |  $103.6\text{\,}\mathrm{\SIUnitSymbolDegree}$,$88.1\text{\,}\mathrm{\SIUnitSymbolDegree}$222XRD spectroscopy of hydrated cycloamylose 102 |  $90\text{\,}\mathrm{\SIUnitSymbolDegree}$,$80\text{\,}\mathrm{\SIUnitSymbolDegree}$333AA simulation of 9-mer Amylose in water 39
1–6 $\psi$$\text{\,}\mathrm{\SIUnitSymbolDegree}$ |  $90\text{\,}\mathrm{\SIUnitSymbolDegree}$,$180\text{\,}\mathrm{\SIUnitSymbolDegree}$111The $\phi,\psi$ torsion angle distribution from the AA simulation of TPS melt at $613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$ |  $89.8\text{\,}\mathrm{\SIUnitSymbolDegree}$,$180\text{\,}\mathrm{\SIUnitSymbolDegree}$444Replica exchange simulation and NMR spectroscopy on $\alpha$-1–6 glucopyranoside dimer 103 | -
1–6 $\phi$$\text{\,}\mathrm{\SIUnitSymbolDegree}$ |  $70\text{\,}\mathrm{\SIUnitSymbolDegree}$111The $\phi,\psi$ torsion angle distribution from the AA simulation of TPS melt at $613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$ |  $70\text{\,}\mathrm{\SIUnitSymbolDegree}$444Replica exchange simulation and NMR spectroscopy on $\alpha$-1–6 glucopyranoside dimer 103 | -
$\rho$,$\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$ | 1.43555TPS melt density ($\rho$) calculated at $300\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$ | 1.4766670:30 ($wt$$\text{\,}\mathrm{\char 37\relax}$) amylopectin-sorbitol melt 37, 38 | 1.4766670:30 ($wt$$\text{\,}\mathrm{\char 37\relax}$) amylopectin-sorbitol melt 37, 38
$T_{\mathrm{g}}$,$\text{\,}\mathrm{K}$ | 325777Glass transition temperature ($T_{\mathrm{g}}$) at the experimental cooling rate was obtained from the combination of temperature-density variation plot and WLF equation | 300-340888Native barley starch having $29\text{\,}\mathrm{\char 37\relax}$ glycerol 104 | 32166670:30 ($wt$$\text{\,}\mathrm{\char 37\relax}$) amylopectin-sorbitol melt 37, 38
${E_{\mathrm{{TPS}}}}$,$\text{\,}\mathrm{GPa}$ | 5.2999The Young’s moduli ($E_{\mathrm{MMT}}$, $E_{\mathrm{TPS}}$) was calculated from the stress-strain curve at $300\text{\,}\mathrm{K}$. | 1.2566670:30 ($wt$$\text{\,}\mathrm{\char 37\relax}$) amylopectin-sorbitol melt 37, 38 | 4.366670:30 ($wt$$\text{\,}\mathrm{\char 37\relax}$) amylopectin-sorbitol melt 37, 38
${E_{\mathrm{{MMT}}}}$,$\text{\,}\mathrm{GPa}$ | 171.95999The Young’s moduli ($E_{\mathrm{MMT}}$, $E_{\mathrm{TPS}}$) was calculated from the stress-strain curve at $300\text{\,}\mathrm{K}$ |  | 171101010Density functional theory (DFT) calculation on MMT 105

A TPS melt density of $1.43\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$ at
$300\text{\,}\mathrm{K}$ compares well with a value of
$1.47\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$ obtained from both
experimental measurements and simulations for a 70:30
($wt$$\text{\,}\mathrm{\char 37\relax}$) amylopectin-sorbitol blend at
$300\text{\,}\mathrm{K}$ 37, 38. The available experimental data on the glass
transition temperature ($T_{\mathrm{g}}$) and Young’s modulus ($E_{i}$), both
directly linked to the nature of intermolecular interactions and atomic-scale
packing, allow an additional test for force field parameters and simulation
convergence. We obtained $T_{\mathrm{g}}=$405\text{\,}\mathrm{K}$$ from
temperature-density variation obtained at a cooling rate of
$3.75\text{\times}{10}^{11}\text{\,}\mathrm{K}\text{\,}{\mathrm{min}}^{-1}$ in
AA simulation of TPS melt (Figure S3), which yields a
$T_{\mathrm{g}}=$325\text{\,}\mathrm{K}$$ on scaling to the experimental
cooling rate of $10\text{\,}\mathrm{K}\text{\,}{\mathrm{min}}^{-1}$ using the
WLF equation (Equation 2). This is in excellent agreement with a value of
$321\text{\,}\mathrm{K}$ obtained in a simulation of 70:30 amylopectin-
sorbitol blend 38 and experimental estimates in range
$300\text{\,}\mathrm{K}$–$340\text{\,}\mathrm{K}$ for barley starch with
$29\text{\,}\mathrm{\char 37\relax}$ glycerol (a small molecule plasticizer
like sorbitol) 104. $E_{\mathrm{TPS}}=$5.2\text{\,}\mathrm{GPa}$$ and
$\sigma_{\mathrm{max}}=$270\text{\,}\mathrm{MPa}$$ obtained from the stress-
strain variation for TPS melt at $300\text{\,}\mathrm{K}$ (Figure S2) compare
well with $E=$4.3\text{\,}\mathrm{GPa}$$ and
$\sigma_{\mathrm{max}}=$285\text{\,}\mathrm{MPa}$$ from AA simulation of 70:30
amylopectin-sorbitol blend at $300\text{\,}\mathrm{K}$, but is only within the
same order of magnitude of the experimental estimate of
$E=$1.25\text{\,}\mathrm{GPa}$$ 38. The MMT sheet modulus of
$171.95\text{\,}\mathrm{GPa}$ from AA simulation matches the density
functional theory (DFT) estimate of $178.4\text{\,}\mathrm{GPa}$ 105 (both
calculated at $300\text{\,}\mathrm{K}$). Table 1 gives an itemized comparison
of various properties for TPS melt calculated in this work and those reported
in other studies using experimental measurements and simulations. These
favorable comparisons for the local conformation of the polymer chains as well
as several macroscopic properties establish the high accuracy of TPS AA
parameters used in this study.

###  3.2 Development of CG Parameters

The developed parameter set is intended for integration with the MARTINI-2 CG
forcefield as available for several polymers and nanoparticles of interest,
enabling direct application to their composites with TPS. We have not used the
more recent MARTINI-3 CG parameterization 106, since it does not yet have a
polarizable model for water.

####  3.2.1 CG Non-bonded Parameters for TPS melt

The three diol groups of sorbitol (Figure 1 (c)) were modeled by a P4 bead
with its self-interaction lowered from level I to level II
($\epsilon$=$4.5\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$), since an
unphysical aggregation was observed in a low-concentration aqueous solution
85. The refined polarizable model was used for water 98, 86, where one neutral
bead with self-interaction level V and cross-interaction level
($60\text{\,}\mathrm{\char 37\relax}$ to $95\text{\,}\mathrm{\char 37\relax}$
of standard MARTINI water P4 bead) is connected to two virtual beads
($q$=$\pm$0.46, no LJ interaction). This was shown to provide a better
representation of dielectric behaviour of water. The initial bead choice in
polymer (amylose and amylopectin) was taken from an existing set of MARTINI-2
CG parameters for polysaccharides (referred to as CGLITa), obtained by
parameterization on AA simulations of mono- and disaccharides in polar and
apolar solvents 84. Accordingly, hemiacetal (Figure 1(a): beads 3, 4, and 9),
acetal (beads 1, 6, and 7), and diol groups (beads 2, 5, and 8) were modeled
by P1, P2, and P4 beads, respectively. This is referred to as CGLITa parameter
set. In this study, we used the melt density and the two-body excess entropy
(higher values indicate stronger pair correlations) for polymer-polymer (given
higher weight in evaluation) and polymer-sorbitol bead pairs to further
optimize the bead-type assignments for amylose and amylopectin.

The CGLITa melt density of
$1136\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ underestimates the melt
density of $1214.69\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ in AA
simulations. The $S_{2}$ (Table S1) for polymer-polymer beads is in excellent
agreement with AA values ($3\text{\,}\mathrm{\char 37\relax}$ higher), but
there is a significant deviation for the polymer-sorbitol beads
($108\text{\,}\mathrm{\char 37\relax}$ lower). An improved density estimate of
$1253\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ is obtained on using
s-type for all polymer beads (CGLITa-S), a recommended choice for several ring
molecules 76. All $S_{2}$ values involving a polymer bead became higher
(resulting from a weakening of interaction potential), leading to a poorer
estimate for polymer-polymer beads ($29\text{\,}\mathrm{\char 37\relax}$
higher than AA) but an improved estimate for polymer-sorbitol beads
($11\text{\,}\mathrm{\char 37\relax}$ lower than AA). A similar improvement in
density, but an incorrect prediction of the crystallization behavior of a
highly concentrated sugar solution was observed on replacing regular beads
with small s-type beads 84. Alternatively, on using the s-type bead only for
the hemiacetal group (referred to as CGLITb), as per the bead assignment used
for lipopolysaccharides 107, a density of
$1159.36\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ is obtained. This
represented a small improvement over CGLITa, retained almost equivalent
accuracy for polymer-polymer $S_{2}$, but polymer-sorbitol $S_{2}$ still
showed a high deviation ($86\text{\,}\mathrm{\char 37\relax}$ lower than AA).
We made two additional changes in the bead assignments: acetal group as Na
bead (recommended for obtaining accurate crystal behavior of cellulose 108)
and one level reduction (from level I to level II) in the polymer-polymer diol
bead (P4) self-interaction (shown to provide an accurate Virial coefficient
and prevent unphysical aggregation of polysaccharides 79, 80). The parameter
set obtained after these two changes in CGLITb is referred to as the CG1
parameter set (Table S3). The CG1 melt density of
$1141.75\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ represented only a
small improvement over CGLITa, a slightly increased polymer-polymer $S_{2}$
deviation ($9\text{\,}\mathrm{\char 37\relax}$ lower than AA), and a somewhat
improved polymer-sorbitol $S_{2}$ ($61\text{\,}\mathrm{\char 37\relax}$
lower). Since CG1 provided the best balance in accuracy of $S_{2}$ for
polymer-polymer and polymer-sorbitol pairs, we took it forward for further
evaluation and optimization.

####  3.2.2 CG Bonded Parameters for TPS melt

Figure 4:  Bonded distribution for amylopectin AA ( ), CG1 ( ), and CG2 ( )
probability distributions, $P(\phi)$ for torsions ($\phi$), and $P(\theta)$
for three-bead angles ($\theta$) were estimated from TPS melt simulations at
$613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$. I(a) 2–1–4–5, I(b)
2–1–4–6, and I(c) 3–1–4–5 represent torsion distribution around the 1–4
linkage, and I(d) 3–1–6–4 around the 1–6 linkage. II(a) 1–4–7 and II(b) 4–6–1
represent angle distribution along 1–4 and 1–6 linkages, respectively. II(c)
1–4–5 and II(d) 1–4–6 represent intra-ring angle distributions. Bead numbers
are as per Figure 1. The beads involved in the calculated quantity are shown
in red and rest in cyan.

The average 1–4 glycosidic bond length in AA simulation was
$0.48\text{\,}\mathrm{nm}$, considerably smaller than
$0.56\text{\,}\mathrm{nm}$ in CG1, while the average 1-4-7 angle of
$125\text{\,}\mathrm{\SIUnitSymbolDegree}$ is close to the CG1 value for an
amylose oligomer in an apolar solvent
($120\text{\,}\mathrm{\SIUnitSymbolDegree}$). Accordingly, we set
$l_{0}=$0.48\text{\,}\mathrm{nm}$$ and
$\theta_{0}=$125\text{\,}\mathrm{\SIUnitSymbolDegree}$$ for the corresponding
bonded potential parameters in Equation 1. Further, the force constant
($k_{\theta}$) for 1-4-5 and 1-4-6 angles needed to be doubled (w.r.t. to the
CG1 value) to reproduce the AA angle distribution in TPS melt (Figure 4 (II)).
This re-parameterized set, referred to as CG2 (Table S2), provided an improved
representation of the three-bead angle distributions for amylose and
amylopectin in TPS melt (Figure 4 (II)) and S6). The torsion distribution for
the 1-4 and 1-6 linkages are well-represented by a proper dihedral potential
with multiplicity 1 (Figure 4 (I)) and are taken without change from CGLITa.
Secondary peaks in the AA distribution, as seen here also, are not present in
the CG distributions since the dihedral potential of multiplicity 1 is used in
MARTINI FF to prevent overfitting. The CG2 parameter set also provided a
slightly improved TPS melt density of
$1153.74\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ (comparable to CGLITb,
better than CGLITa and CG1) and is the only set that provided a good estimate
for $S_{2}$ of both polymer-polymer ($2\text{\,}\mathrm{\char 37\relax}$
higher) and polymer-sorbitol ($22\text{\,}\mathrm{\char 37\relax}$ lower).
Further, amylose chain persistence length, $l_{p}$, of
$0.86\text{\,}\mathrm{nm}$ in AA simulations compares more favorably with CG2
($0.90\text{\,}\mathrm{nm}$) than CG1 ($0.96\text{\,}\mathrm{nm}$). For
comparison, simulations of amorphous amylose melt at $713\text{\,}\mathrm{K}$
gave a value of $0.9\text{\,}\mathrm{nm}$ with AA and $1\text{\,}\mathrm{nm}$
with CG MARTINI 41. The chain radius of gyration, ${\langle
R_{\mathrm{g}}\rangle}$, also follows the same trend for both amylose and
amylopectin: within $5\text{\,}\mathrm{\char 37\relax}$ of AA in CG2 and
$15\text{\,}\mathrm{\char 37\relax}$ of AA in CG1 (Table 2). The 1–4
glycosidic bond length mismatch for the CG1 model (Figure 2(II)) is the likely
reason for these small but consistent differences. A lower angle stiffness and
larger bond length for 1–4 bonds also lead to $\sim$ $10\text{\,}\mathrm{\char
37\relax}$-$15\text{\,}\mathrm{\char 37\relax}$ higher $S_{\mathrm{c}}$ for
amylose and amylopectin chains with CG1 parameters compared to AA (Table 2). A
correction for both bonded parameters in CG2 set gave $S_{\mathrm{c}}$ values
within $5\text{\,}\mathrm{\char 37\relax}$ of AA. Simultaneous improvements in
several structural and thermodynamic properties make the CG2 parameter set a
distinct improvement over existing parameters for TPS.

####  3.2.3 CG Non-bonded Parameters for TPS-MMT composite

Table 2: Comparison of properties determined from AA and CG simulations of TPS
melt and TPS-MMT-TMA composite. The radius of gyration (${\langle
R_{\mathrm{g}}\rangle}$), conformational entropy per CG bead
($S_{\mathrm{c}}$), and normalized diffusion coefficient of the polymer
trimeric segments ($D_{\mathrm{n}}$) were calculated at
$613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$, with average and
standard deviation calculated from $900\text{\,}\mathrm{ns}$ trajectories. The
green, cyan, red, and blue values indicate a relative error (in percent) of
$\leq$$5\text{\,}\mathrm{\char 37\relax}$, 5–$10\text{\,}\mathrm{\char
37\relax}$, 10–$15\text{\,}\mathrm{\char 37\relax}$, and
$\geq$$15\text{\,}\mathrm{\char 37\relax}$, in CG values (w.r.t. AA),
respectively.
| Amylose | Amylopectin
---|---|---
| near | far | melt |  | near | far | melt
${\langle R_{\mathrm{g}}\rangle}$($\text{\,}\mathrm{nm}$) | AA | 1.35(0.09) | 1.37(0.08) | 1.38(0.02) |  | 1.74(0.15) | 1.77(0.16) | 1.83(0.02)
| CG1 | 1.49(0.18) | 1.57(0.10) | 1.58(0.05) |  | 1.95(0.25) | 2.01(0.19) | 2.07(0.03)
| CG1red | 1.53(0.13) | 1.56(0.15) | 1.58(0.05) |  | 1.99(0.23) | 1.99(0.23) | 2.07(0.03)
| CG2 | 1.34 (0.19) | 1.45(0.12) | 1.43(0.03) |  | 1.78(0.19) | 1.84(0.22) | 1.83(0.06)
| CG2red | 1.38(0.16) | 1.40(0.16) | 1.43(0.03) |  | 1.78(0.22) | 1.84(0.21) | 1.83(0.06)
${\langle S_{c}\rangle}$ | AA | 46.11(1.39) | 46.73(1.64) | 46.26(0.14) |  | 45.82(2.70) | 46.63(3.10) | 46.23(0.11)
$($\text{\,}\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}\text{\,}\mathrm{K}$)$ | CG1 | 44.44(4.67) | 48.68(2.98) | 52.16(0.11) |  | 47.09(3.12) | 51.41(3.54) | 51.90(0.09)
| CG1red | 48.23(1.69) | 48.31(2.52) | 52.16(0.11) |  | 50.30(3.61) | 51.58(2.99) | 51.90(0.09)
| CG2 | 41.58(2.92) | 45.71(2.51) | 47.00(0.29) |  | 43.25(2.21) | 48.07(2.12) | 45.59(0.25)
| CG2red | 44.53(2.17) | 45.18(3.18) | 47.00(0.29) |  | 46.21(2.77) | 48.21(2.44) | 45.59(0.250)
|  | near | overall | $D_{\mathrm{melt}}$ |  | near | overall | $D_{\mathrm{melt}}$
${\langle D_{n}\rangle}$ | AA | 0.96 | 0.83 | 0.063(0.003) |  | 1.16 | 0.79 | 0.033(0.003)
| CG1 | - | 0.47 | 0.033(0.001) |  | - | 0.54 | 0.020(0.001)
| CG1red | 0.85 | 0.78 | 0.033(0.001) |  | 1.11 | 0.98 | 0.020(0.001)
| CG2 | - | 0.41 | 0.017(0.001) |  | - | 0.41 | 0.012(0.001)
| CG2red | 0.88 | 0.92 | 0.017(0.0010) |  | 0.97 | 0.79 | 0.012(0.001)

Three polymer-surface interaction dependent properties, viz. polymer
normalized diffusion coefficient, $D_{\mathrm{n}}$, radius of gyration,
$\langle R_{g}\rangle$, and conformational entropy, $S_{c}$, were
simultaneously optimized to parameterize the cross-interaction for TPS-MMT-TMA
beads. A lower $\langle R_{g}\rangle$ and $S_{c}$ near the MMT sheet in CG
simulations compared to the bulk (not observed in AA) suggest a disc-like
flattened configuration near the surface to maximize polymer-sheet contacts
(Table 2), which also gave a much lower $D_{n}$ for polymer trimeric segments
and sorbitol molecules in CG simulations (Figure S9 and Table 2). This change
in polymer chain conformation also affects the RDF and $S_{2}$ for all self-
and cross-pairs of polymer and sorbitol beads, as discussed later in section
3.3.3. All these observations indicate that MARTINI-2 forcefield parameters
overestimate solid-fluid interactions for the TPS-MMT system. A reduction in
solid-fluid interaction parameter ($\epsilon_{ij}$) was shown to be essential
for accurate prediction of adsorption enthalpy of small organic molecules
(butane, benzene, hexadecane), order-disorder transitions of hexadecane, and
preferential adsorption of long-chain alkanes on graphite surfaces 82. Here,
the effective dispersive component of $\epsilon_{ij}$ was re-scaled for all
TPS-MMT pairs as per Equation 5, maintaining consistency with the
parameterization strategy used for the TMA-MMT sheet 87.

${\epsilon_{\mathrm{i,j}}}={\gamma(\epsilon_{\mathrm{N0,N0}})+(\epsilon_{\mathrm{i,j}}-\epsilon_{\mathrm{N0,N0}})},$
(5)

where, $i$ and $j$ represent TPS and MMT beads, respectively, and $\gamma$ is
the scaling factor. A previous study has shown that taking the N0 bead self-
interaction, $\epsilon_{\mathrm{N0,N0}}$, as purely dispersive and scaling the
dispersive component for all interaction pairs as per Equation 5 gave an
accurate distribution of polar and dispersive contributions to the cleavage
free energy of two TMA-MMT sheets 87. $D_{\mathrm{n}}$, $\langle
R_{g}\rangle$, and $S_{c}$ were calculated for the TPS-MMT system at six
$\gamma$ values in range 0.15-1 (Table S5 and Figure S9). The CG parameters
were evaluated by using the property value in the composite normalized by the
corresponding value in the melt (for the self-diffusion coefficient) or the
property value in the near-sheet region normalized by the value in the far
region (for the radius of gyration and the conformational entropy). This
separated the effect of differences between AA and CG properties in the bulk
systems (melt, far region) from the differences in the composite system,
allowing specific optimization of TPS and MMT cross-interactions. For the CG1
set, $\langle R_{g}\rangle$ and $S_{c}$ ratios at $\gamma=0.21$ and
$D_{\mathrm{n}}$ at $\gamma=0.25$ were closest to those in AA simulations,
with $\langle R_{g}\rangle$ and $S_{c}$ at $\gamma=0.25$ and $D_{\mathrm{n}}$
at $\gamma=0.21$ also in very good agreement with AA. A re-scaling factor of
$\gamma=0.21$ was found most suitable for self-interaction of MMT beads 87,
and therefore, the same factor was chosen here for the TPS-MMT cross-
interactions for the CG1 parameter set (referred to as CG1red). For the CG2
parameter set, all three properties are in excellent agreement with AA at
$\gamma=0.15$ and represent a distinct improvement over all other choices, and
therefore, selected as CG2 parameter set for the composite system (referred to
as CG2red).

The selected parameter set that concurrently optimize structural,
thermodynamic, and dynamic properties are expected to provide an accurate
estimate for other properties not used in optimization and also have better
state-point transferability 77, 109, 110. For example, the addition of
miscibility data in MARTINI-3 FF development 106 allowed for a significant
improvement in some of the key limitations of MARTINI-2 FF concerning large
self-interaction of some species (such as polysaccharides) in solution 79, 80,
poor reproduction of liquid properties at interfaces 82.

###  3.3 TPS-MMT Composite Properties from the CG Model

Several properties of the composite system obtained using the different CG FFs
were compared with the corresponding values obtained from AA simulations, and
the inter-relationship between these properties was investigated. A single
clay sheet sandwiched between two large TPS regions forms the model system
S11, which provides the simplest setup for separate determination of
properties in the near-sheet and far-sheet (”bulk”) regions while also
allowing for establishing a chemical equilibrium between these regions.

####  3.3.1 TPS Components’ Density Profiles Above the TMA-MMT Sheet

Figure 5:  Density profile in TPS-MMT composite. Normalized density
perpendicular to the TMA-MMT sheet (placed in $xy$ plane),
$\rho_{\mathrm{n}}(z)$, for (a,d) polymer, (b,e) sorbitol, and (c,f) water. CG
parameter set 1 is in (a-c) and set 2 in (d-f), with AA as  , CG as  , and
CGred as  .

The two plasticizers, water and sorbitol, show a sharp density peak near the
TMA-MMT surface (Figure 5), driven by the interaction between their highly
polar hydroxyl groups and TMA-MMT silicate layer with mobile ions. This
preferential adsorption of water and sorbitol leads to depletion of polymer
chains at the sheet surface, but polymer density shows a sharp increase within
the first solvation layer, increasing to $\sim$50\text{\,}\mathrm{\char
37\relax}$$ of the bulk value. This mixing of polymer and low molecular weight
plasticizers in the interfacial region is possibly driven by the presence of
polar beads, viz. diol (P4), hemiacetal (P1), and acetal (Na or P2), groups in
the polymer chain. This provides one reason for an observed increase in TPS
intercalation with increased water content in a TPS-cloisite (bentonite
nanoclay) composite 19. Figure 5(a, d) shows that the scale of interactions
(regular ($\gamma=1$) or reduced) have essentially no effect on the normalized
density of polymer beads. This directly follows from high polymer depletion at
the sheet surface, which reduces the direct effect of TPS-TMA-MMT bead
interactions. The choice of CG model (CG1 versus CG2) also has a negligible
effect, partly because of the density scaling by corresponding values in the
melt system. Overall, all four CG models (CG1/2, regular/reduced) show an
excellent agreement with CG and AA density profiles, $\rho_{\mathrm{n}}(z)$,
for polymer and sorbitol. However, there is a significant increase in water
density at the clay surface (Figure 5(c, f)) for the CGred models. This
implies that TPS-TMA-MMT cross-interactions are not fully optimized, with very
low water content ($1\text{\,}\mathrm{\char 37\relax}$) possibly accentuating
the effect of sub-optimal interactions. All component densities reach a
plateau (close to the bulk value) at an average $z$-distance of
$2.5\text{\,}\mathrm{nm}$ from the TMA-MMT sheet, indicating minimal influence
of the clay sheet beyond this distance. Therefore, the plane at
$z=$2.5\text{\,}\mathrm{nm}$$ was used to distinguish between the near-clay
and far regions for subsequent calculation of all reported properties.

####  3.3.2 Polymer Chain Conformation and Dynamics

Average chain conformation was characterized using radius of gyration,
topological constraints and correlations using conformational entropy, and
overall polymer dynamics using segment diffusion coefficients, all expected to
show a significant dependence on polymer-surface interactions. In the AA
simulations, the near-to-far ${\langle R_{\mathrm{g}}\rangle}$ and
$S_{\mathrm{c}}$ ratios, and amylose and amylopectin trimeric segment $D_{n}$
were close to unity. It is expected that favorable polymer-surface
interactions lead to chain expansion, lower $S_{\mathrm{c}}$, and lower
$D_{n}$ in the near-surface regions (w.r.t. bulk) while unfavorable
interactions lead to the opposite effect. For example, AA simulations showed
that strongly attractive polymer-surface interactions for polyethylene oxide-
Na-MMT 111 lead to doubling of ${\langle R_{\mathrm{g}}\rangle}$ and a
decrease in diffusion coefficient in this system as well as in polyethylene-
graphite composite 112. Weakly attractive interactions for polyethylene,
polypropylene, and polystyrene near the TMA-MMT sheet lead to a
$\sim$10\text{\,}\mathrm{\char 37\relax}$$ decrease in both $S_{\mathrm{c}}$
and diffusion coefficients 87. Conversely, the addition of repulsive
nanoparticles to a LJ polymer accelerated chain dynamics 113, 114 but had a
negligible effect on ${\langle R_{\mathrm{g}}\rangle}$ 115, 113. For the TPS-
MMT system, the interplay of polymer-surface (weakly repulsive), plasticizer-
surface (attractive), and polymer-plasticizer (weakly attractive) interactions
lead to polymer depletion at the clay surface and a sharp increase in the
density in the first solvation shell of plasticizers at the surface (Figure
5). This tempering of direct polymer-clay interactions and polymer-plasticizer
mixing in the near clay region gave comparable ${\langle
R_{\mathrm{g}}\rangle}$, $S_{\mathrm{c}}$, and $D_{n}$ values in the near-clay
and the far regions. In a similar ternary composite system of cellulose
nanocrystals (CNC)-polyvinylpyrrolidone (PVP)-polycaprolactone (PCL), PVP
${\langle R_{\mathrm{g}}\rangle}$ increased from $\sim$
$0.90\text{\,}\mathrm{nm}$ (unfavorable PVP-CNC interactions lowered ${\langle
R_{\mathrm{g}}\rangle}$ w.r.t. bulk) to $1.65\text{\,}\mathrm{nm}$ on addition
of dichloromethane (DCM), a good solvent for PVP that pushed the PVP chains
away from CNC 116. The use of standard TPS-MMT cross-interactions in CG1 and
CG2 gave smaller near to far ${\langle R_{\mathrm{g}}\rangle}$ and
$S_{\mathrm{c}}$ ratios than AA: in the range 0.92-0.97 for ${\langle
R_{\mathrm{g}}\rangle}$ and 0.9-0.92 for $S_{\mathrm{c}}$ in CG compared to
$\sim$ 0.98 for both ratios in AA. The $D_{\mathrm{n}}$ in the composite was
significantly lower ($\sim 0.4-0.5$) than that in AA ($\sim 0.8-0.9$) (Figure
S9). As discussed in Section 3.2.3, the solid-fluid interactions are typically
overestimated in MARTINI FF 98, 82, leading to chain flattening and lower
$S_{\mathrm{c}}$ and $D_{\mathrm{n}}$. Therefore, TPS-MMT cross-interactions
were rescaled as per Equation 5, which has built-in asymmetry leading to a
smaller reduction in the pair interaction strength involving more polar beads.
An excellent agreement with AA values was obtained for ${\langle
R_{\mathrm{g}}\rangle}$, $S_{\mathrm{c}}$, and $D_{n}$ at $\gamma=0.21$ for
CG1 (referred to as CG1red) and $\gamma=0.15$ for CG2 (referred to as CG2red).
The high level of rescaling is directly linked to the high covalent
coordination state of the polymer chain and the MMT beads 111A different value
of re-scale parameter $\gamma$ for cross-interactions involving non-ring
sorbitol and water beads would have been more appropriate but potentially
affected transferability and therefore, was not pursued, shown to
significantly lower the dispersive interactions in earlier studies 87, 91.
This use of multiple properties in CG parameter optimization is expected to
provide improved state-point transferability as shown in previous studies, for
example, use of multiple state point optimization 77 or addition of density
correction 109 in IBI-derived potentials, use of polymer hydration
thermodynamics and structural properties in a CG FF combining MARTINI FF with
IBI 110.

####  3.3.3 Radial Distribution Function and Two-body Excess Entropy

Figure 6:  Radial distribution functions in TPS melt and TPS–MMT-TMA composite. RDFs for polymer-polymer (P-P), polymer-sorbitol (P-S), and sorbitol-sorbitol (S-S) beads in the melt and in near- and far-regions in the composite (2D RDF): (a-c) $g_{\alpha\beta}$ for AA ( ), CG1 ( ), and CG1red ( ), and (d-f) $g_{\alpha\beta}$ for AA ( ), CG2 ( ), and CG2red ( ). Four bonded atoms were excluded from all RDF calculations. Table 3:  Normalized two-body excess entropy in TPS–TMA-MMT composite system. ${S_{2,\alpha,\beta}}$ were calculated from the average RDFs shown in Figure 6. The composite ${S_{2,\alpha,\beta}}$ values are reported after dividing with the respective ${S_{2,\alpha,\beta}}$ of TPS melt (Table S1). P and S refer to polymer (amylose and amylopectin) and sorbitol, respectively. Pairs | near | far
---|---|---
| P-P | P-S | S-S | P-P | P-S | S-S
AA |  | 0.87 | 1.23 | 1.07 | 0.79 | 1.19 | 1.23
CG1 |  | 1.81 | 2.03 | 1.32 | 0.63 | 0.79 | 0.81
CG1red |  | 1.69 | 1.05 | 0.83 | 0.60 | 0.79 | 0.81
CG2 |  | 2.00 | 2.60 | 1.16 | 0.60 | 0.91 | 0.86
CG2red |  | 1.14 | 1.45 | 0.95 | 0.60 | 0.91 | 0.84

We have further evaluated CG1 and CG2 parameter sets by comparing estimates
for radial distribution function and two-body excess entropy in the composite
system to corresponding AA simulations, neither of which were used in
parameter optimization. In the TPS-MMT composite system, the polymer-polymer
(PP), polymer-sorbitol (PS), and sorbitol-sorbitol (SS) 2D $RDFs$ in the far
region were identical to the melt (Figure 6), indicating a negligible effect
of clay on the structural properties of TPS at a distance (along sheet normal)
of $\sim 2\langle R_{\mathrm{g}}\rangle$ from the sheet. The corresponding
RDFs for CG1 and CG2 parameters are in good agreement, with the size of the
first solvation shell and the order of first-peak heights (SS ¿ PS ¿ PP) the
same as AA simulations. However, at the full level of TPS-MMT interaction
(used in CG1 and CG2), the near-sheet $RDFs$ show long-range structuring for
all three pairs (Figure 6). A 1:1 comparison of AA and CG RDFs is not easy
since the peaks show a broadening in AA RDFs because of the use of COM of
atoms making the CG bead. Therefore, individual pairwise
($S_{\mathrm{2,\alpha,\beta}}$) and total two-body excess entropy
($S_{\mathrm{2}}$), integrals over RDF calculated as per Equations
LABEL:eqn:tbe,_eq:overalltbe,_eqn:tbecomp, were used for a quantitative
comparison between AA and CG simulations (Table S1). The near-region
$S_{\mathrm{2,\alpha\beta}}$ values for CG1 and CG2 are up to
$108\text{\,}\mathrm{\char 37\relax}$ and $130\text{\,}\mathrm{\char
37\relax}$ higher (in magnitude) than AA simulations, respectively.An increase
in peak height and long-range structuring were directly linked to highly
attractive liquid-solid surface interactions 113, a known drawback of MARTINI
2 forcefield as discussed for CG parameterization. The use of CG1red and
CG2red parameter sets lead to both a lowering of first-peak height and a
reduction in long-range structuring in near-region RDFs (Figure 6). The
relative errors (w.r.t AA) in near-region $S_{\mathrm{2,\alpha,\beta}}$
calculated for CG1red simulations were $94\text{\,}\mathrm{\char 37\relax}$,
$15\text{\,}\mathrm{\char 37\relax}$, and $22\text{\,}\mathrm{\char 37\relax}$
for polymer-polymer, polymer-sorbitol, and sorbitol-sorbitol respectively (an
improvement of $14\text{\,}\mathrm{\char 37\relax}$, $50\text{\,}\mathrm{\char
37\relax}$, and $1\text{\,}\mathrm{\char 37\relax}$ over CG1), while for
CG2red simulations were $31\text{\,}\mathrm{\char 37\relax}$,
$18\text{\,}\mathrm{\char 37\relax}$, and $11\text{\,}\mathrm{\char 37\relax}$
(an improvement of $99\text{\,}\mathrm{\char 37\relax}$,
$93\text{\,}\mathrm{\char 37\relax}$, and $-3\text{\,}\mathrm{\char 37\relax}$
(this became slightly less accurate) over CG2), respectively. The success of
the developed parameter set in accurately capturing diverse material
properties beyond those used for optimization suggests their broad
applicability for material characterization.

##  4 Conclusions

Enthalpy-based coarse-graining methods, such as DPD, require a state-dependent
effective interaction parameter, wherein the phase space for its determination
is significantly expanded for high system heterogeneity, for example, in the
case of a TPS-MMT-TMA composite investigated in this work. Chemical-specific
coarse-grained (CG) models, such as IBI, force-matching, or MARTINI, provide
more effective parameterization for such systems. Here, we used AA simulation
to predict several TPS properties, spanning from macroscopic (the glass
transition temperature, Young’s modulus) to microscopic (conformation along
1–4 and 1–6 glycosidic linkages). These properties were used as constraints in
developing an accurate MARTINI-2 FF parameter set for TPS melt. While existing
MARTINI FF parameters obtained using the data from AA simulations of mono- and
disaccharides in polar and apolar solvents 84 form a good starting point,
several important differences between AA and CG simulations were observed for
chain conformation and thermodynamic properties of a TPS melt. We used the
data from AA simulations in the present work and bead-type assignments for
similar functional groups used in literature 84, 107, 108, 79, 80 to obtain a
new parameter set for TPS melt, referred to as the CG2 parameter set. It
included a shortening of the CG bond involving the 1–4 glycosidic link, an
s-type bead for the hemiacetal group, Na bead type for the acetal group, one
level reduction in the self-interaction of polymer diol bead, a small increase
in the 1-4-7 angle from its value for amylose in apolar solvent, and a
doubling of the force constants for 1-4-5 and 1-4-6 angles. The CG2 model
provided simultaneous improvement in several structural and thermodynamic
properties of TPS melt, making it a distinct improvement over existing MARTINI
FF parameters for TPS.

The AA simulation of the TPS-MMT system revealed that the interplay of
polymer-surface (weakly repulsive), plasticizer-surface (attractive), and
polymer-plasticizer (weakly attractive) interactions lead to preferential
adsorption of the plasticizer and depletion of polymer at the clay surface,
and mixing of polymer and plasticizer within the first solvation shell, a
feature captured very well by the optimized CG models used in this study. This
unique density profile also resulted in comparable $R_{g}$, $S_{c}$, and
$D_{n}$ in near-clay and far regions, making it distinct from the expected
behavior for attractive or repulsive polymer-surface interactions in simpler
systems. However, the CG simulation of the TPS-MMT composite showed compaction
and flattening of polymer chains (lower $R_{g}$) and significantly reduced
dynamics (lower $S_{c}$ and $D_{n}$) in the near-sheet region. The use of
liquid-liquid partition coefficients in MARTINI-2 parameterization makes it
inherently less suitable for solid-liquid interfaces, with an unphysical
freezing of liquid at a solid interface reported in several previous studies.
In the present case, a high covalent coordination state for the polymer chains
and MMT sheet, shown to significantly lower the dispersive interactions,
further accentuated this effect. The downscaling of the effective dispersive
component of the TPS-MMT MARTINI-2 interaction parameter led to simultaneous
improvement in $R_{g}$, $S_{c}$, and $D_{n}$, with the selected CG2red
parameter set giving $D_{n}$ for amylose within $10\text{\,}\mathrm{\char
37\relax}$ of AA and all other values (for both amylose and amylopectin)
within $5\text{\,}\mathrm{\char 37\relax}$. CG1red and CG2red models also
significantly improved over the corresponding full interaction models (CG1 and
CG2) in predicting two-body excess entropy (calculated from the radial
distribution function). Thus, the developed MARTINI-2 CG parameters for TPS-
MMT composites, which offer close to a hundred-fold speedup over AA
simulations, can be used to accurately estimate properties not used in
parameterization.

We acknowledge that rescaling the sheet-TPS interaction is a sub-optimal
solution that potentially limits transferability to systems with very
different components. The optimal bead type assignments determined in this
study, combined with a significant expansion of bead types and availability of
scale factors for cross- and self-interaction in the recently released
MARTINI-3 FF can allow for more rational development of parameters for various
TPS-clay composites.

The authors thank the IIT Delhi High-Performance Computing (HPC) facility for
providing the computational resources.

## References

  * Carvalho 2008 Carvalho, A. J. Starch: major sources, properties and applications as thermoplastic materials. _Monomers, polymers and composites from renewable resources_ 2008, 321–342
  * Jiménez et al. 2012 Jiménez, A.; Fabra, M.; Talens, P.; Chiralt, A. Edible and biodegradable starch films: a review. Food Bioprocess Technol 5 (6): 2058–2076. 2012
  * Attaran et al. 2017 Attaran, S. A.; Hassan, A.; Wahit, M. U. Materials for food packaging applications based on bio-based polymer nanocomposites: A review. _Journal of Thermoplastic Composite Materials_ 2017, _30_ , 143–173
  * Babu et al. 2013 Babu, R. P.; O’connor, K.; Seeram, R. Current progress on bio-based polymers and their future trends. _Progress in biomaterials_ 2013, _2_ , 1–16
  * Khan et al. 2017 Khan, B.; Bilal Khan Niazi, M.; Samin, G.; Jahan, Z. Thermoplastic starch: A possible biodegradable food packaging material—A review. _Journal of Food Process Engineering_ 2017, _40_ , e12447
  * Stepto 2003 Stepto, R. The processing of starch as a thermoplastic. Macromolecular Symposia. 2003; pp 203–212
  * de Freitas et al. 2022 de Freitas, A. d. S. M.; da Silva, A. P. B.; Montagna, L. S.; Nogueira, I. A.; Carvalho, N. K.; de Faria, V. S.; Dos Santos, N. B.; Lemes, A. P. Thermoplastic starch nanocomposites: sources, production and applications–a review. _Journal of Biomaterials Science, Polymer Edition_ 2022, _33_ , 900–945
  * Kaboorani et al. 2021 Kaboorani, A.; Gray, N.; Hamzeh, Y.; Abdulkhani, A.; Shirmohammadli, Y. Tailoring the low-density polyethylene-thermoplastic starch composites using cellulose nanocrystals and compatibilizer. _Polymer Testing_ 2021, _93_ , 107007
  * Corrêa et al. 2022 Corrêa, A. C.; de Campos, A.; Claro, P. I. C.; Guimarães, G. G. F.; Mattoso, L. H. C.; Marconcini, J. M. Biodegradability and nutrients release of thermoplastic starch and poly ($\varepsilon$-caprolactone) blends for agricultural uses. _Carbohydrate Polymers_ 2022, _282_ , 119058
  * Cyras and Manfredi 2008 Cyras, V. P.; Manfredi, L. B. Physical and mechanical properties of thermoplastic starch / montmorillonite nanocomposite films. 2008, _73_ , 55–63
  * Chocyk et al. 2015 Chocyk, D.; Gladyszewska, B.; Ciupak, A.; Oniszczuk, T.; Moscicki, L.; Rejak, A. Influence of water addition on mechanical properties of thermoplastic starch foils. _International Agrophysics_ 2015, _29_
  * Domene-López et al. 2019 Domene-López, D.; García-Quesada, J. C.; Martin-Gullon, I.; Montalbán, M. G. Influence of starch composition and molecular weight on physicochemical properties of biodegradable films. _Polymers_ 2019, _11_ , 1084
  * Salaberria et al. 2015 Salaberria, A. M.; Diaz, R. H.; Labidi, J.; Fernandes, S. C. Role of chitin nanocrystals and nanofibers on physical, mechanical and functional properties in thermoplastic starch films. _Food Hydrocolloids_ 2015, _46_ , 93–102
  * Rezaei et al. 2015 Rezaei, M.; Ismail, A. F.; Bakeri, G.; Hashemifard, S. A.; Matsuura, T. Effect of general montmorillonite and Cloisite 15A on structural parameters and performance of mixed matrix membranes contactor for CO2 absorption. _Chemical Engineering Journal_ 2015, _260_ , 875–885
  * Zhou and Xu 2015 Zhou, M.; Xu, D. Starch-MMT composite films: Effects of bio-inspired modification on MMT. _Starch-Stärke_ 2015, _67_ , 470–477
  * Suter et al. 2015 Suter, J. L.; Groen, D.; Coveney, P. V. Chemically specifi C multiscale modeling of clay-polymer nanocomposites reveals intercalation dynamics, tactoid self-assembly and emergent materials properties. _Advanced Materials_ 2015, _27_ , 966–984
  * Altayan and Al Darouich 2022 Altayan, M. M.; Al Darouich, T. Toward reducing the food packaging waste impact: a study on the effect of Starch type and PE type in thermoplastic starch-polyethylene blends. _Chemical Papers_ 2022, _76_ , 2447–2457
  * Baumberger et al. 1998 Baumberger, S.; Lapierre, C.; Monties, B.; Della Valle, G. Use of kraft lignin as filler for starch films. _Polymer Degradation and Stability_ 1998, _59_ , 273–277
  * Park et al. 2002 Park, H.-M.; Li, X.; Jin, C.-Z.; Park, C.-Y.; Cho, W.-J.; Ha, C.-S. Preparation and properties of biodegradable thermoplastic starch/clay hybrids. _Macromolecular Materials and Engineering_ 2002, _287_ , 553–558
  * Schlemmer et al. 2010 Schlemmer, D.; Angélica, R. S.; Sales, M. J. A. Morphological and thermomechanical characterization of thermoplastic starch/montmorillonite nanocomposites. _Composite Structures_ 2010, _92_ , 2066–2070
  * Majdzadeh-Ardakani et al. 2010 Majdzadeh-Ardakani, K.; Navarchian, A. H.; Sadeghi, F. Optimization of mechanical properties of thermoplastic starch/clay nanocomposites. _Carbohydrate Polymers_ 2010, _79_ , 547–554
  * Rivadeneira-Velasco et al. 2021 Rivadeneira-Velasco, K. E.; Utreras-Silva, C. A.; Díaz-Barrios, A.; Sommer-Márquez, A. E.; Tafur, J. P.; Michell, R. M. Green nanocomposites based on thermoplastic starch: A review. _Polymers_ 2021, _13_ , 3227
  * Huang et al. 2004 Huang, M. F.; Yu, J. G.; Ma, X. F. Studies on the properties of Montmorillonite-reinforced thermoplastic starch composites. _Polymer_ 2004, _45_ , 7017–7023
  * Krishnamoorti et al. 1996 Krishnamoorti, R.; Vaia, R. A.; Giannelis, E. P. Structure and dynamics of polymer-layered silicate nanocomposites. _Chemistry of Materials_ 1996, _8_ , 1728–1734
  * Alias et al. 2021 Alias, A. H.; Norizan, M. N.; Sabaruddin, F. A.; Asyraf, M. R. M.; Norrrahim, M. N. F.; Ilyas, A. R.; Kuzmin, A. M.; Rayung, M.; Shazleen, S. S.; Nazrin, A.; others Hybridization of MMT/lignocellulosic fiber reinforced polymer nanocomposites for structural applications: a review. _Coatings_ 2021, _11_ , 1355
  * Ayana et al. 2014 Ayana, B.; Suin, S.; Khatua, B. Highly exfoliated eco-friendly thermoplastic starch (TPS)/poly (lactic acid)(PLA)/clay nanocomposites using unmodified nanoclay. _Carbohydrate polymers_ 2014, _110_ , 430–439
  * Mao et al. 2012 Mao, C.; Zhu, Y.; Jiang, W. Design of electrical conductive composites: Tuning the morphology to improve the electrical properties of graphene filled immiscible polymer blends. _ACS Applied Materials and Interfaces_ 2012, _4_ , 5281–5286
  * Huang and Zhou 2022 Huang, F.; Zhou, S. Molecular Dynamics Simulation of Coiled Carbon Nanotube Pull-Out from Matrix. _International Journal of Molecular Sciences_ 2022, _23_ , 9254
  * Guo et al. 2023 Guo, L.; Xu, H.; Wu, N.; Yuan, S.; Zhou, L.; Wang, D.; Wang, L. Molecular dynamics simulation of the effect of the thermal and mechanical properties of addition liquid silicone rubber modified by carbon nanotubes with different radii. _e-Polymers_ 2023, _23_ , 20228105
  * Lin and Wang 2023 Lin, K.; Wang, Z. Multiscale mechanics and molecular dynamics simulations of the durability of fiber-reinforced polymer composites. _Communications Materials_ 2023, _4_ , 66
  * Gartner III and Jayaraman 2019 Gartner III, T. E.; Jayaraman, A. Modeling and simulations of polymers: a roadmap. _Macromolecules_ 2019, _52_ , 755–786
  * Guvench et al. 2008 Guvench, O.; Greene, S. N.; Kamath, G.; Brady, J. W.; Venable, R. M.; Pastor, R. W.; Mackerell Jr, A. D. Additive empirical force field for hexopyranose monosaccharides. _Journal of computational chemistry_ 2008, _29_ , 2543–2564
  * Hatcher et al. 2009 Hatcher, E.; Guvench, O.; MacKerell Jr, A. D. CHARMM additive all-atom force field for aldopentofuranoses, methyl-aldopentofuranosides, and fructofuranose. _The Journal of Physical Chemistry B_ 2009, _113_ , 12466–12476
  * Geronimo et al. 2019 Geronimo, I.; Ntarima, P.; Piens, K.; Gudmundsson, M.; Hansson, H.; Sandgren, M.; Payne, C. M. Kinetic and molecular dynamics study of inhibition and transglycosylation in Hypocrea jecorina family 3 $\beta$-glucosidases. _Journal of Biological Chemistry_ 2019, _294_ , 3169–3180
  * Kirschner et al. 2008 Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; González-Outeiriño, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J. GLYCAM06: a generalizable biomolecular force field. Carbohydrates. _Journal of computational chemistry_ 2008, _29_ , 622–655
  * Plazinski et al. 2016 Plazinski, W.; Lonardi, A.; Hünenberger, P. H. Revision of the GROMOS 56A6CARBO force field: improving the description of ring-conformational equilibria in hexopyranose-based carbohydrates chains. _Journal of computational chemistry_ 2016, _37_ , 354–365
  * Özeren et al. 2020 Özeren, H. D.; Olsson, R. T.; Nilsson, F.; Hedenqvist, M. S. Prediction of plasticization in a real biopolymer system (starch) using molecular dynamics simulations. _Materials & Design_ 2020, _187_ , 108387
  * Özeren et al. 2020 Özeren, H. D.; Guivier, M.; Olsson, R. T.; Nilsson, F.; Hedenqvist, M. S. Ranking Plasticizers for Polymers with Atomistic Simulations: PVT, Mechanical Properties, and the Role of Hydrogen Bonding in Thermoplastic Starch. _ACS Applied Polymer Materials_ 2020, _2_ , 2016–2026
  * Koneru et al. 2019 Koneru, J. K.; Zhu, X.; Mondal, J. Quantitative Assessment of the Conformational Heterogeneity in Amylose across Force Fields. _Journal of Chemical Theory and Computation_ 2019, _15_ , 6203–6212
  * Mazeau and Heux 2003 Mazeau, K.; Heux, L. Molecular dynamics simulations of bulk native crystalline and amorphous structures of cellulose. _The Journal of Physical Chemistry B_ 2003, _107_ , 2394–2403
  * Gatsiou et al. 2022 Gatsiou, C. A.; Bick, A.; Krokidis, X.; Economou, I. G. Atomistic and Coarse-Grained Simulations of Bulk Amorphous Amylose Above and Below the Glass Transition. _Macromolecules_ 2022,
  * Groot and Warren 1997 Groot, R. D.; Warren, P. B. Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. _Journal of Chemical Physics_ 1997, _107_ , 4423–4435
  * Singh et al. 2018 Singh, A.; Chakraborti, A.; Singh, A. Role of a polymeric component in the phase separation of ternary fluid mixtures: A dissipative particle dynamics study. _Soft Matter_ 2018, _14_ , 4317–4326
  * Fu et al. 2013 Fu, Y.; Liao, L.; Yang, L.; Lan, Y.; Mei, L.; Liu, Y.; Hu, S. Molecular dynamics and dissipative particle dynamics simulations for prediction of miscibility in polyethylene terephthalate/polylactide blends. _Molecular Simulation_ 2013, _39_ , 415–422
  * Scocchi et al. 2007 Scocchi, G.; Posocco, P.; Fermeglia, M.; Pricl, S. Polymer - Clay nanocomposites: A multiscale molecular modeling approach. _Journal of Physical Chemistry B_ 2007, _111_ , 2143–2151
  * Long and Sotta 2006 Long, D.; Sotta, P. Nonlinear and plastic behavior of soft thermoplastic and filled elastomers studied by dissipative particle dynamics. _Macromolecules_ 2006, _39_ , 6282–6297
  * Sliozberg et al. 2020 Sliozberg, Y. R.; Gair Jr, J. L.; Hsieh, A. J. Dissipative particle dynamics simulation of microphase separation in polyurethane urea nanocomposites. _Polymer_ 2020, _193_ , 122339
  * Khani et al. 2015 Khani, S.; Jamali, S.; Boromand, A.; Hore, M. J.; Maia, J. Polymer-mediated nanorod self-assembly predicted by dissipative particle dynamics simulations. _Soft Matter_ 2015, _11_ , 6881–6892
  * Wang et al. 2021 Wang, J.; Han, Y.; Xu, Z.; Yang, X.; Ramakrishna, S.; Liu, Y. Dissipative particle dynamics simulation: A review on investigating mesoscale properties of polymer systems. _Macromolecular Materials and Engineering_ 2021, _306_ , 2000724
  * Ju et al. 2013 Ju, S.-P.; Wang, Y.-C.; Huang, G.-J.; Chang, J.-W. Miscibility of graphene and poly (methyl methacrylate)(PMMA): molecular dynamics and dissipative particle dynamics simulations. _Rsc Advances_ 2013, _3_ , 8298–8307
  * Schweizer et al. 1992 Schweizer, K. S.; Honnell, K. G.; Curro, J. G. Reference interaction site model theory of polymeric liquids: Self-consistent formulation and nonideality effects in dense solutions and melts. _The Journal of chemical physics_ 1992, _96_ , 3211–3225
  * Schweizer et al. 1997 Schweizer, K. S.; Fuchs, M.; Szamel, G.; Guenza, M.; Tang, H. Polymer-mode-coupling theory of the slow dynamics of entangled macromolecular fluids. _Macromolecular theory and simulations_ 1997, _6_ , 1037–1117
  * Zirkel et al. 2002 Zirkel, A.; Gruner, S.; Urban, V.; Thiyagarajan, P. Small-angle neutron scattering investigation of the Q-dependence of the Flory- Huggins interaction parameter in a binary polymer blend. _Macromolecules_ 2002, _35_ , 7375–7386
  * Wu 2006 Wu, J. Density functional theory for chemical engineering: From capillarity to soft materials. _AIChE journal_ 2006, _52_ , 1169–1193
  * Oxtoby 2002 Oxtoby, D. W. Density functional methods in the statistical mechanics of materials. _Annual Review of Materials Research_ 2002, _32_ , 39
  * Helfand and Sapse 1975 Helfand, E.; Sapse, A. M. Theory of unsymmetric polymer–polymer interfaces. _The Journal of chemical physics_ 1975, _62_ , 1327–1331
  * Helfand and Wasserman 1976 Helfand, E.; Wasserman, Z. R. Block copolymer theory. 4. Narrow interphase approximation. _Macromolecules_ 1976, _9_ , 879–888
  * Lowden and Chandler 1973 Lowden, L. J.; Chandler, D. Solution of a new integral equation for pair correlation functions in molecular liquids. _The Journal of Chemical Physics_ 1973, _59_ , 6587–6595
  * Sung and Yethiraj 2005 Sung, B. J.; Yethiraj, A. Integral equation theory of randomly coupled multiblock copolymer melts: Effect of block size on the phase behavior. _The Journal of chemical physics_ 2005, _123_ , 214901
  * Hall and Schweizer 2011 Hall, L. M.; Schweizer, K. S. Impact of monomer sequence, composition and chemical heterogeneity on copolymer-mediated effective interactions between nanoparticles in melts. _Macromolecules_ 2011, _44_ , 3149–3160
  * Martin and Jayaraman 2016 Martin, T. B.; Jayaraman, A. Using theory and simulations to calculate effective interactions in polymer nanocomposites with polymer-grafted nanoparticles. _Macromolecules_ 2016, _49_ , 9684–9692
  * Hsu and Chandler 1979 Hsu, C.; Chandler, D. RISM calculation of the structure of liquid chloroform. _Molecular Physics_ 1979, _37_ , 299–301
  * Martin et al. 2018 Martin, T. B.; Gartner III, T. E.; Jones, R. L.; Snyder, C. R.; Jayaraman, A. pyPRISM: a computational tool for liquid-state theory calculations of macromolecular materials. _Macromolecules_ 2018, _51_ , 2906–2922
  * Clark and Guenza 2010 Clark, A. J.; Guenza, M. G. Mapping of polymer melts onto liquids of soft-colloidal chains. 2010, 1–13
  * Clark et al. 2012 Clark, A. J.; Mccarty, J.; Lyubimov, I. Y.; Guenza, M. G. Thermodynamic Consistency in Variable-Level Coarse Graining of Polymeric Liquids. 2012, _168301_ , 1–5
  * Lyubimov et al. 2010 Lyubimov, I. Y.; Mccarty, J.; Clark, A.; Guenza, M. G. Analytical rescaling of polymer dynamics from mesoscale simulations. 2010, 1–6
  * Lyubimov and Guenza 2013 Lyubimov, I. Y.; Guenza, M. G. Theoretical reconstruction of realistic dynamics of highly coarse-grained cis-1 , 4-polybutadiene melts Theoretical reconstruction of realistic dynamics of highly coarse-grained. 2013, _546_
  * McCarty and Guenza 2010 McCarty, J.; Guenza, M. G. Multiscale modeling of binary polymer mixtures: Scale bridging in the athermal and thermal regime. _Journal of Chemical Physics_ 2010, _133_
  * Karatrantos et al. 2016 Karatrantos, A.; Clarke, N.; Kröger, M. Modeling of polymer structure and conformations in polymer nanocomposites from atomistic to mesoscale: A Review. _Polymer reviews_ 2016, _56_ , 385–428
  * Zeng et al. 2008 Zeng, Q.; Yu, A.; Lu, G. Multiscale modeling and simulation of polymer nanocomposites. _Progress in polymer science_ 2008, _33_ , 191–269
  * Gooneie et al. 2017 Gooneie, A.; Schuschnigg, S.; Holzer, C. A review of multiscale computational methods in polymeric materials. _Polymers_ 2017, _9_ , 16
  * Glotzer and Paul 2002 Glotzer, S. C.; Paul, W. Molecular and mesoscale simulation methods for polymer materials. _Annual Review of materials research_ 2002, _32_ , 401–436
  * Reith et al. 2001 Reith, D.; Meyer, H.; Müller-Plathe, F. Mapping atomistic to coarse-grained polymer models using automatic simplex optimization to fit structural properties. _Macromolecules_ 2001, _34_ , 2335–2345
  * Reith et al. 2003 Reith, D.; Pütz, M.; Müller-Plathe, F. Deriving effective mesoscale potentials from atomistic simulations. _Journal of computational chemistry_ 2003, _24_ , 1624–1636
  * Izvekov and Voth 2005 Izvekov, S.; Voth, G. A. A multiscale coarse-graining method for biomolecular systems. _The Journal of Physical Chemistry B_ 2005, _109_ , 2469–2473
  * Marrink et al. 2007 Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; De Vries, A. H. The MARTINI force field: Coarse grained model for biomolecular simulations. _Journal of Physical Chemistry B_ 2007, _111_ , 7812–7824
  * Moore et al. 2014 Moore, T. C.; Iacovella, C. R.; McCabe, C. Derivation of coarse-grained potentials via multistate iterative Boltzmann inversion. _The Journal of chemical physics_ 2014, _140_
  * Panizon et al. 2015 Panizon, E.; Bochicchio, D.; Monticelli, L.; Rossi, G. MARTINI coarse-grained models of polyethylene and polypropylene. _The Journal of Physical Chemistry B_ 2015, _119_ , 8209–8216
  * Schmalhorst et al. 2017 Schmalhorst, P. S.; Deluweit, F.; Scherrers, R.; Heisenberg, C. P.; Sikora, M. Overcoming the Limitations of the MARTINI Force Field in Simulations of Polysaccharides. _Journal of Chemical Theory and Computation_ 2017, _13_ , 5039–5053
  * Shivgan et al. 2020 Shivgan, A. T.; Marzinek, J. K.; Huber, R. G.; Krah, A.; Henchman, R. H.; Matsudaira, P.; Verma, C. S.; Bond, P. J. Extending the Martini coarse-grained force field to N-Glycans. _Journal of Chemical Information and Modeling_ 2020, _60_ , 3864–3883
  * Wu and Yang 2012 Wu, D.; Yang, X. Coarse-grained molecular simulation of self-assembly for nonionic surfactants on graphene nanostructures. _The Journal of Physical Chemistry B_ 2012, _116_ , 12048–12056
  * Gobbo et al. 2013 Gobbo, C.; Beurroies, I.; de Ridder, D.; Eelkema, R.; Marrink, S. J.; De Feyter, S.; van Esch, J. H.; de Vries, A. H. MARTINI model for physisorption of organic molecules on graphite. _The Journal of Physical Chemistry C_ 2013, _117_ , 15623–15631
  * Piskorz et al. 2019 Piskorz, T. K.; Gobbo, C.; Marrink, S. J.; De Feyter, S.; De Vries, A. H.; Van Esch, J. H. Nucleation mechanisms of self-assembled physisorbed monolayers on graphite. _The Journal of Physical Chemistry C_ 2019, _123_ , 17510–17520
  * López et al. 2009 López, C. A.; Rzepiela, A. J.; de Vries, A. H.; Dijkhuizen, L.; Hünenberger, P. H.; Marrink, S. J. Martini coarse-grained force field: Extension to carbohydrates. _Journal of Chemical Theory and Computation_ 2009, _5_ , 3195–3210
  * Sukenik et al. 2015 Sukenik, S.; Sapir, L.; Harries, D. Osmolyte induced changes in peptide conformational ensemble correlate with slower amyloid aggregation: a coarse-grained simulation study. _Journal of Chemical Theory and Computation_ 2015, _11_ , 5918–5928
  * Michalowsky et al. 2017 Michalowsky, J.; Schäfer, L. V.; Holm, C.; Smiatek, J. A refined polarizable water model for the coarse-grained MARTINI force field with long-range electrostatic interactions. _The Journal of Chemical Physics_ 2017, _146_ , 054501
  * Khan and Goel 2019 Khan, P.; Goel, G. Martini Coarse-Grained Model for Clay-Polymer Nanocomposites. _Journal of Physical Chemistry B_ 2019, _123_ , 9011–9023
  * Guvench et al. 2009 Guvench, O.; Hatcher, E.; Venable, R. M.; Pastor, R. W.; MacKerell Jr, A. D. CHARMM additive all-atom force field for glycosidic linkages between hexopyranoses. _Journal of chemical theory and computation_ 2009, _5_ , 2353–2370
  * 89 Guvench, O.; Mallajosyula, S.; Raman, E.; Hatcher, E.; Vanommeslaeghe, K.; Foster, T.; Jamison, F. II and Mackerell AD Jr (2011) CHARMM additive all-atom force field for carbohydrate derivatives and its utility in polysaccharide and carbohydrate-protein modeling. _J Chem Theory Comput_ _7_ , 3162–3180
  * Jo et al. 2008 Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: a web-based graphical user interface for CHARMM. _Journal of computational chemistry_ 2008, _29_ , 1859–1865
  * Heinz et al. 2005 Heinz, H.; Koerner, H.; Anderson, K. L.; Vaia, R. A.; Farmer, B. Force field for mica-type silicates and dynamics of octadecylammonium chains grafted to montmorillonite. _Chemistry of materials_ 2005, _17_ , 5658–5669
  * Martínez et al. 2009 Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: a package for building initial configurations for molecular dynamics simulations. _Journal of computational chemistry_ 2009, _30_ , 2157–2164
  * Lyckfeldt and Ferreira 1998 Lyckfeldt, O.; Ferreira, J. M. Processing of Porous Ceramics by ’Starch Consolidation’. _Journal of the European Ceramic Society_ 1998, _18_ , 131–140
  * Abraham et al. 2015 Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. _SoftwareX_ 2015, _1_ , 19–25
  * Bussi et al. 2007 Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. _The Journal of chemical physics_ 2007, _126_
  * Berendsen et al. 1984 Berendsen, H. J.; Postma, J. v.; Van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. _The Journal of chemical physics_ 1984, _81_ , 3684–3690
  * Parrinello and Rahman 1981 Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. _Journal of Applied physics_ 1981, _52_ , 7182–7190
  * Yesylevskyy et al. 2010 Yesylevskyy, S. O.; Schäfer, L. V.; Sengupta, D.; Marrink, S. J. Polarizable water model for the coarse-grained MARTINI force field. _PLoS computational biology_ 2010, _6_ , e1000810
  * Williams et al. 1955 Williams, M. L.; Landel, R. F.; Ferry, J. D. The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids. _Journal of the American Chemical society_ 1955, _77_ , 3701–3707
  * Sharma et al. 2008 Sharma, R.; Agarwal, M.; Chakravarty, C. Estimating the entropy of liquids from atom–atom radial distribution functions: Silica, beryllium fluoride and water. _Molecular Physics_ 2008, _106_ , 1925–1938
  * Schlitter 1993 Schlitter, J. Estimation of absolute and relative entropies of macromolecules using the covariance matrix. _Chemical physics letters_ 1993, _215_ , 617–621
  * Gessler et al. 1999 Gessler, K.; Usón, I.; Takaha, T.; Krauss, N.; Smith, S. M.; Okada, S.; Sheldrick, G. M.; Saenger, W. V-amylose at atomic resolution: X-ray structure of a cycloamylose with 26 glucose residues (cyclomaltohexaicosaose). _Proceedings of the National Academy of Sciences of the United States of America_ 1999, _96_ , 4246–4251
  * Patel et al. 2014 Patel, D. S.; Pendrill, R.; Mallajosyula, S. S.; Widmalm, G.; MacKerell, A. D. Conformational properties of $\alpha$\- Or $\beta$-(1→6)-linked oligosaccharides: Hamiltonian replica exchange MD simulations and NMR experiments. _Journal of Physical Chemistry B_ 2014, _118_ , 2851–2871
  * Forssell et al. 1997 Forssell, P. M.; Mikkilä, J. M.; Moates, G. K.; Parker, R. Phase and glass transition behaviour of concentrated barley starch-glycerol-water mixtures, a model for thermoplastic starch. _Carbohydrate Polymers_ 1997, _34_ , 275–282
  * Zartman et al. 2010 Zartman, G. D.; Liu, H.; Akdim, B.; Pachter, R.; Heinz, H. Nanoscale tensile, shear, and failure properties of layered silicates as a function of cation density and stress. _The Journal of Physical Chemistry C_ 2010, _114_ , 1763–1772
  * Souza et al. 2021 Souza, P. C.; Alessandri, R.; Barnoud, J.; Thallmair, S.; Faustino, I.; Grünewald, F.; Patmanidis, I.; Abdizadeh, H.; Bruininks, B. M.; Wassenaar, T. A.; others Martini 3: a general purpose force field for coarse-grained molecular dynamics. _Nature methods_ 2021, _18_ , 382–388
  * Hsu et al. 2016 Hsu, P.-C.; Jefferies, D.; Khalid, S. Molecular dynamics simulations predict the pathways via which pristine fullerenes penetrate bacterial membranes. _The journal of physical chemistry B_ 2016, _120_ , 11170–11179
  * Wohlert and Berglund 2011 Wohlert, J.; Berglund, L. A. A coarse-grained model for molecular dynamics simulations of native cellulose. _Journal of Chemical Theory and Computation_ 2011, _7_ , 753–760
  * Hu et al. 2019 Hu, C.; Lu, T.; Guo, H. Developing a transferable coarse-grained model for the prediction of thermodynamic, structural, and mechanical properties of polyimides at different thermodynamic state points. _Journal of chemical information and modeling_ 2019, _59_ , 2009–2025
  * Banerjee et al. 2018 Banerjee, P.; Roy, S.; Nair, N. Coarse-grained molecular dynamics force-field for polyacrylamide in infinite dilution derived from iterative Boltzmann inversion and MARTINI force-field. _The Journal of Physical Chemistry B_ 2018, _122_ , 1516–1524
  * Suter and Coveney 2009 Suter, J. L.; Coveney, P. V. Computer simulation study of the materials properties of intercalated and exfoliated poly (ethylene) glycol clay nanocomposites. _Soft Matter_ 2009, _5_ , 2239–2251
  * Daoulas et al. 2005 Daoulas, K. C.; Harmandaris, V. A.; Mavrantzas, V. G. Detailed atomistic simulation of a polymer melt/solid interface: structure, density, and conformation of a thin film of polyethylene melt adsorbed on graphite. _Macromolecules_ 2005, _38_ , 5780–5795
  * Smith et al. 2002 Smith, G. D.; Bedrov, D.; Li, L.; Byutner, O. A molecular dynamics simulation study of the viscoelastic properties of polymer nanocomposites. _The Journal of chemical physics_ 2002, _117_ , 9478–9489
  * Varnik et al. 2002 Varnik, F.; Baschnagel, J.; Binder, K. Reduction of the glass transition temperature in polymer films: A molecular-dynamics study. _Physical Review E_ 2002, _65_ , 021507
  * Karatrantos et al. 2015 Karatrantos, A.; Clarke, N.; Composto, R. J.; Winey, K. I. Polymer conformations in polymer nanocomposites containing spherical nanoparticles. _Soft Matter_ 2015, _11_ , 382–388
  * Voronova et al. 2021 Voronova, M. I.; Gurina, D. L.; Surov, O. V.; Zakharov, A. G. Interactions in solvent–polycaprolactone–cellulose nanocrystals–polyvinyl pyrrolidone system: Experiment and molecular dynamics simulation. _Journal of Molecular Liquids_ 2021, _341_ , 117409

##  5 Supplementary Material

###  5.1 Temperature-pressure annealing cycles

The initial box, obtained using Packmol software, underwent a multi-step
density increase followed by multiple temperature-pressure (T-P) cycles to
obtain an equilibrated TPS structure free from the memory effect of the
initial configuration. The coordinates of the simulation box were rescaled in
steps of $29.17\text{\,}{\mathrm{nm}}^{3}$ followed by energy minimization at
each step using the steepest decent algorithm. The T-P annealing cycles were
performed on the rescaled TPS system (state A) to accelerate the slow dynamics
of the starch chains, leading to a more rapid attainment of the equilibrated
configuration. (Figure S1). The simulation time at the lower
($613\text{\,}\mathrm{K}$) and higher temperature ($613\text{\,}\mathrm{K}$)
of the annealing cycle was chosen to ensure local equilibration (potential
energy plateau) and Brownian relaxation of the chain, respectively (Figure S1
(A)). Following the completion of each temperature cycle, the pressure was
incrementally doubled. Once a pressure of $1000\text{\,}\mathrm{bar}$ was
reached, the system was depressurized back to $1\text{\,}\mathrm{bar}$ using
the same stepwise approach (Figure S1 (C)). The complete T-P cycles were
repeated until the density of the initial (state A) and final (state A’)
states matched (Figure S1 (C)).

Figure S1: Temperature-pressure annealing simulations A The $MSD$ (at
$713\text{\,}\mathrm{K}$) and potential energy plateau (at
$613\text{\,}\mathrm{K}$), used to estimate simulation time and temperature in
annealing cycle required in temperature annealing cycles (shown in (B)) to
attain Brownian relaxation time and local equilibration of TPS melt,
respectively. (B) In the temperature annealing cycle, states A and B represent
the rescaled input and output simulation boxes, respectively. (C) T-P
annealing cycle workflow along with the pressure values.

###  5.2 Key Results

Figure S2: Stress-strain behavior of the TPS melt and MMT sheet (a) The linear
part of TPS stress-strain data (represented by  ) was used to estimate Young’s
modulus. The onset of the plateau (maximum stress) is represented by  . (b)
The AA ( ) and CG ( ) stress data of the TMA-MMT sheet. The system was
stretched using a series of $0.1\text{\,}\mathrm{\char 37\relax}$
unidirectional (in the $x$-direction) strains (at $300\text{\,}\mathrm{K}$,
and $1\text{\,}\mathrm{bar}$) and the corresponding Virial ($P_{xx}$) values
at each step were used to estimate the stress. Figure S3: Density-temperature
plot for TPS melt The TPS melt system was cooled from $613\text{\,}\mathrm{K}$
to $150\text{\,}\mathrm{K}$ in steps of $25\text{\,}\mathrm{K}$, where
$4\text{\,}\mathrm{ns}$ NPT simulation was performed at each step, with last
$500\text{\,}\mathrm{ps}$ trajectory at each step used for density
calculation. The intersection point (highlighted in red) of rubbery and glassy
region lines is used as an estimate of glass-transition temperature. Figure
S4: Chemical structure of amylose and amylopectin . A chemical representation
of (a) amylose and (b) amylopectin chain consisting of four and seven
$\alpha$-D glucose monomers, respectively. Figure S5: AA representation of 1–4
and 1–6 linkage . The 1–4 $\alpha$-D glucose linkage is represented by
$\phi=\angle\mathrm{O_{5}}-\mathrm{C_{4}}-\mathrm{O_{1}}-\mathrm{C_{1}}$,
$\psi=\angle\mathrm{C_{1}}-\mathrm{O_{1}}-\mathrm{C_{4}}-\mathrm{C_{3}}$, and
1–6 $\alpha$-D glucose linkage is represented by
$\phi=\angle\mathrm{O_{5}}-\mathrm{C_{1}}-\mathrm{O_{1}}-\mathrm{C_{6}}$,
$\psi=\angle\mathrm{C_{1}}-\mathrm{O_{1}}-\mathrm{C_{6}}-\mathrm{C_{5}}$.

Figure S6: AA and CG angle distribution of amylopectin y-axis show the
probability distribution ($P(l)$) of 1–4 glycosidic bond ($l$) for AA( ), CG1(
), and CG2( ) models. The distribution obtained from TPS melt simulation at
$613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{b}\mathrm{a}\mathrm{r}$

.

Figure S7: MSD for TPS-MMT composite system  The MSD of amylose and
amylopectin trimeric segments were estimated for CG1 (a,b) and CG2 (c,d)
models. MSDs in Near-MMT regions for AA ( ), CG ( ), and CGred ( ) and far-MMT
regions for AA ( ), CG ( ), and CGred ( ) are estimated from TPS–TMA-MMT
simulation at $613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$

.

Figure S8: MSD of TPS components in melt and composite systems MSD for AA ( ),
CG1 ( ), CG1red ( ), CG2 ( ), and CG2red ( ) were estimated at
$613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$. In the composite
system, the MSD was calculated in the $x$-$y$ plane using an equilibrated
trajectory of $900\text{\,}\mathrm{ns}$

.

Table S1:  Two-body excess entropy of TPS melt system The two-body excess entropies for polymer-sorbitol pairs were calculated by integrating the RDF (Figure 6) of the corresponding pairs in the TPS melt, using Equation 3. System |  |  | TPS melt
---|---|---|---
pairs |  |  | Polymer-Polymer | Polymer-Sorbitol | Sorbitol-Sorbitol
AA |  |  |  | -1.33 | -0.36 | -0.52
CGLITa |  |  |  | -1.29 | -0.75 | -0.57
CGLITb |  |  |  | -1.39 | -0.67 | -0.61
CGSSa-S |  |  |  | -0.95 | -0.40 | -0.95
CG1 |  |  |  | -1.21 | -0.58 | -0.83
CG2 |  |  |  | -1.30 | -0.44 | -0.98

The $S_{2,\alpha\beta}$ of polymer-sorbitol pair were calculated using
Equation S1

$\mathrm{S_{2,\alpha\beta}}=-\pi\rho\int_{0}^{\infty}g_{\alpha\beta}(r)\mathrm{ln}\,g_{\alpha\beta}(r)-[g_{\alpha\beta}(r)-1]rdr$
(S1)

The diffusion coefficients were calculated using the Equation S2.

$\mathrm{MSD(t)}=\mathrm{6\cdot D\cdot t}$ (S2)

Figure S9:  Normalized diffusion coefficient of TPS components The diffusion
coefficients of TPS components in melt ($D_{0}$) and composite ($D_{c}$) were
obtained from the slope of ${MSD}$ (shown in Figure S8) in the linear regime.
The $MSDs$ were calculated from $900\text{\,}\mathrm{ns}$ long simulation
trajectories, obtained at a temperature of $613\text{\,}\mathrm{K}$ and a
pressure of $1\text{\,}\mathrm{bar}$ Figure S10: AA and CG representation of
TMA-MMT sheet with the Mapping scheme. (a) shows the AA and CG representation
of the TMA-MMT sheet. (b) shows the mapping scheme and MARTINI bead type of
the highlighted portion of the sheet in (a) Figure S11: AA and CG
representation of TPS-TMA–MMT composite system. (a) shows the AA
representation of the composite system. (b) shows the same system in CG
representation. The red, blue, yellow, and lime colors represent the polymer,
sorbitol, water, and MMT, respectively.

###  5.3 CG-MARTINI parameters

Table S2: CG MARTINI bonded parameter 1-4 and 1-6 glycosidic links  The bonded
parameter obtained from the AA simulation of TPS melt system at
$613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$. The columns show force
constants $k_{l}$, $k_{\theta}$, and $k_{\phi}$ of bonds ($l_{0}$), angles
($\theta$), and dihedrals ($\phi_{0}$), respectively. Figure 1 illustrates the
nomenclature used to represent beads.

.

type | topological pattern | parameters |  |
---|---|---|---|---
bonds |  | $l_{o}$ |  $k_{l}$ [$\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}\text{\,}{\mathrm{nm}}^{-4}$] | function |
| 1-2 | 0.238 | 30000 | 1 |
| 1-3 | 0.276 | 30000 | 1 |
| 4-5 | 0.224 | 30000 | 1 |
| 4-6 | 0.244 | 30000 | 1 |
1-4 link | 1-4 | 0.482 | 30000 | 1 |
1-6 link | 1-6 | 0.300 | 30000 | 1 |
angle |  | $\theta_{0}$ |  $k_{\theta}$ [$\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$] | function |
| 2-1-4 | 70 | 100 | 2 |
| 3-1-4 | 50 | 25 | 2 |
| 1-4-5 | 140 | 100 | 2 |
| 1-4-6 | 140 | 100 | 2 |
1-4 link | 1-4-7 | 125 | 180 | 2 |
1-6 link | 4-6-1 | 100 | 55 | 2 |
| 6-1-2 | 140 | 45 | 2 |
| 6-1-3 | 140 | 35 | 2 |
dihedrals |  | $\phi_{0}$ |  $k_{\phi}$ [$\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$] | function | multiplicity
| 2-1-4-5 | 120 | 8 | 1 | 1
| 2-1-4-6 | -20 | 5 | 1 | 1
| 3-1-4-5 | -40 | 5 | 1 | 1
1-6 link | 4-6-1-2 | -100 | 5 | 1 | 1

The center of mass positions of functional groups (corresponding to CG beads)
in AA trajectories were converted to pseudo-CG trajectories using Equation S3.

$\mathrm{r^{CG}_{i}}=\frac{\sum_{j=1}^{n}{r_{j}}{m_{j}}}{\sum_{j=1}^{n}{m_{j}}}$
(S3)

Where, $m_{j}$ and $\mathrm{r}_{i}$ are the masses and positions of atoms,
$\mathrm{r^{CG}_{i}}$ are the positions for pseudo-CG beads.

Table S3: MARTINI non-bonded parameter of TPS melt. The listed bead diameter ($\sigma_{ij}$) and interaction strength ($\epsilon_{ij}$) were used in stimulating the TPS melt system (for both the CG1 and CG2 models). The bead assignments are shown in Figure 1. The P4c bead is a variant of P4 that employs a weaker level II self-interaction compared to level I in P4 Bead type |  $\sigma$ ($\text{\,}\mathrm{nm}$) |  $\epsilon$ ($\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$)
---|---|---
P4c |  P4c |  0.47 |  4.500
P4c |  Na |  0.47 |  4.000
P4c |  SP1 |  0.47 |  4.500
POL |  SP1 |  0.47 |  3.800
POL |  Na |  0.47 |  3.325
POL |  P4c |  0.47 |  4.275
POL |  POL |  0.47 |  3.100
Table S4: MARTINI non-bonded parameter of TPS-TMA–MMT composite system. The Table illustrates the bead diameter ($\sigma_{ij}$) and the rescaled interaction strengths ($\epsilon_{ij}$) between TPS and MMT beads ($\gamma$=0.15). The SN0, SQ0, and SNa are MMT sheet beads representing $\mathrm{AlO_{2}H}$, $\mathrm{MgO_{2}H}$, and $\mathrm{SiO_{2}}$, respectively.  Bead type |  $\sigma$ ($\text{\,}\mathrm{nm}$) |  $\epsilon$ ($\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$)
---|---|---
SN0 |  P1 |  0.47 |  1.025
SN0 |  P2 |  0.47 |  1.025
SN0 |  P4c |  0.47 |  1.025
SN0 |  Na |  0.47 |  0.525
SN0 |  SP1 |  0.43 |  0.769
SNa |  P1 |  0.47 |  1.525
SNa |  P2 |  0.47 |  1.525
SNa |  P4c |  0.47 |  1.525
SNa |  Na |  0.47 |  1.025
SNa |  SP1 |  0.43 |  1.144
SQ0 |  P1 |  0.47 |  1.025
SQ0 |  P2 |  0.47 |  1.525
SQ0 |  P4c |  0.47 |  2.625
SQ0 |  Na |  0.47 |  1.025
SQ0 |  SP1 |  0.43 |  0.769
SQ0 |  SQ0 |  0.43 |  0.560
SNa |  SNa |  0.43 |  0.935
SN0 |  SN0 |  0.43 |  0.560
Table S5: Comparison of AA and CG properties of TPS-MMT composite system The radius of gyration (${\langle R_{\mathrm{g}}\rangle}$) and configuration entropy per CG bead ($S_{\mathrm{c}}$) were calculated for six different resealing factors ($\gamma$) in range of 0.15 to 1. All properties calculated at $613\text{\,}\mathrm{K}$ and $1\text{\,}\mathrm{bar}$. | Amylose | Amylopectin
---|---|---
| $\gamma$ | near | far | melt |  | near | far | melt
${\langle R_{\mathrm{g}}\rangle}$($\text{\,}\mathrm{nm}$) | AA |  | 1.35 | 1.37 | 1.38 |  | 1.74 | 1.77 | 1.83
| CG1 | 1 | 1.49 | 1.57 | 1.58 |  | 1.95 | 2.01 | 2.07
|  | 0.65 | 1.47 | 1.58 | 1.58 |  | 2.00 | 2.05 | 2.07
|  | 0.44 | 1.51 | 1.58 | 1.58 |  | 2.01 | 2.03 | 2.07
|  | 0.25 | 1.52 | 1.58 | 1.58 |  | 1.97 | 2.01 | 2.07
|  | 0.21 | 1.53 | 1.56 | 1.58 |  | 1.99 | 1.99 | 2.07
|  | 0.15 | 1.52 | 1.59 | 1.58 |  | 2.03 | 2.09 | 2.07
| CG2 | 1 | 1.34 | 1.45 | 1.43 |  | 1.78 | 1.84 | 1.83
|  | 0.65 | 1.39 | 1.47 | 1.43 |  | 1.76 | 1.82 | 1.83
|  | 0.44 | 1.37 | 1.43 | 1.43 |  | 1.77 | 1.83 | 1.83
|  | 0.25 | 1.37 | 1.41 | 1.43 |  | 1.79 | 1.85 | 1.83
|  | 0.21 | 1.40 | 1.43 | 1.43 |  | 1.78 | 1.83 | 1.83
|  | 0.15 | 1.38 | 1.40 | 1.43 |  | 1.78 | 1.84 | 1.83
${\langle S_{c}\rangle}$ | AA |  | 46.11 | 46.73 | 46.26 |  | 45.82 | 46.63 | 46.23
$($\text{\,}\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}\text{\,}\mathrm{K}$)$ | CG1 | 1 | 44.44 | 48.68 | 52.16 |  | 47.09 | 51.41 | 51.90
|  | 0.65 | 45.95 | 49.02 | 52.16 |  | 47.07 | 51.66 | 51.90
|  | 0.44 | 44.81 | 48.74 | 52.16 |  | 48.71 | 51.93 | 51.90
|  | 0.25 | 47.64 | 48.65 | 52.16 |  | 49.68 | 51.61 | 51.90
|  | 0.21 | 48.23 | 48.31 | 52.16 |  | 50.30 | 51.58 | 51.90
|  | 0.15 | 48.31 | 48.54 | 52.16 |  | 50.72 | 51.99 | 51.90
| CG2 | 1 | 41.58 | 45.71 | 47.00 |  | 43.25 | 48.07 | 45.59
|  | 0.65 | 41.97 | 46.23 | 47.00 |  | 41.76 | 48.45 | 45.59
|  | 0.44 | 41.27 | 47.03 | 47.00 |  | 43.35 | 49.23 | 45.59
|  | 0.25 | 43.38 | 45.96 | 47.00 |  | 45.09 | 48.12 | 45.59
|  | 0.21 | 43.19 | 45.96 | 47.00 |  | 45.79 | 48.22 | 45.59
|  | 0.15 | 44.53 | 45.19 | 47.00 |  | 46.21 | 48.21 | 45.59