Development of Coarse-Grained Models for Poly(4-vinylphenol) and Poly(2-vinylpyridine): Polymer Chemistries with Hydrogen Bonding

In this paper, we identify the modifications needed in a recently developed generic coarse-grained (CG) model that captured directional interactions in polymers to specifically represent two exemplary hydrogen bonding polymer chemistries—poly(4-vinylphenol) and poly(2-vinylpyridine). We use atomistically observed monomer-level structures (e.g., bond, angle and torsion distribution) and chain structures (e.g., end-to-end distance distribution and persistence length) of poly(4-vinylphenol) and poly(2-vinylpyridine) in an explicitly represented good solvent (tetrahydrofuran) to identify the appropriate modifications in the generic CG model in implicit solvent. For both chemistries, the modified CG model is developed based on atomistic simulations of a single 24-mer chain. This modified CG model is then used to simulate longer (36-mer) and shorter (18-mer and 12-mer) chain lengths and compared against the corresponding atomistic simulation results. We find that with one to two simple modifications (e.g., incorporating intra-chain attraction, torsional constraint) to the generic CG model, we are able to reproduce atomistically observed bond, angle and torsion distributions, persistence length, and end-to-end distance distribution for chain lengths ranging from 12 to 36 monomers. We also show that this modified CG model, meant to reproduce atomistic structure, does not reproduce atomistically observed chain relaxation and hydrogen bond dynamics, as expected. Simulations with the modified CG model have significantly faster chain relaxation than atomistic simulations and slower decorrelation of formed hydrogen bonds than in atomistic simulations, with no apparent dependence on chain length.


Introduction
Advances in modeling and simulation of polymers over the past few decades have enabled many valuable studies of macromolecular materials over a broad range of relevant length and time scales-from oscillations in bonds and angles at the monomer level, to relaxation and diffusion at the chain level, to the assembly of chains into ordered domains [1][2][3][4][5][6][7][8][9][10][11][12]. Polymer simulations with atomistic models provide chemically detailed representations of monomers but are limited due to the larger computational resources and longer run times needed to predict experimentally relevant phenomena (e.g., disorder to order transition of high molecular weight polymer chains). To probe experimentally relevant length and time scales with reasonable computational resources and run time, one can use coarse-grained (CG) models. CG polymer models reduce the degrees of freedom by representing a polymer as a string of CG beads, where each CG bead represents either groups of atoms within capture directional interactions in polymers are needed that can ultimately guide experiments in screening suitable candidates for material design.
In recent work, Kulshreshtha et al. presented a generic CG model that captured directional interactions in polymers in general and used this model to study structure in polymer nanocomposites (PNCs) containing generic homopolymer-grafted nanoparticles in a homopolymer matrix [109]. In their CG model, the hydrogen bonding "acceptor" and "donor" CG beads were embedded in the graft and matrix "monomer" CG beads, respectively. An effective directional interaction between the graft and matrix CG beads was achieved despite the use of isotropic attractive interaction potential between the acceptor and donor CG beads by careful tuning of the relative size, placement, and bonded interactions of acceptor and donor beads with respect to the graft and matrix polymer beads. Using this CG model in molecular dynamics (MD) simulations, Kulshreshtha et al. showed that attractive directional interactions between graft and matrix chains improved the penetration of the grafted layer by matrix chains (i.e., grafted layer "wetting") in PNCs, as seen with isotropic graft-matrix attractive interaction; however, the directional graft-matrix attraction led to a lesser number of matrix chains interacting with each graft chain and a lower free volume per graft chain at equivalent wetting compared to isotropic graft-matrix attraction [109]. The implications of these results are that the thermomechanical properties for PNCs with hydrogen bonding graft and matrix polymers would be different from those of PNCs with isotropically attractive graft-matrix interaction. This generic CG model of Kulshreshtha et al., capturing directional interactions in polymers, enabled simulation studies of universal structural behavior common to many hydrogen bonding polymers independent of specific polymer chemistry. However, this generic CG model of Kulshreshtha et al. [109] did not include specific bonded constraints (e.g., angle or dihedral potentials to mimic local orientational penalty) that could alter the ability to form a hydrogen bond between two monomers. Furthermore, it also did not account for intra-chain hydrogen bonds. In principle, this generic CG model of Kulshreshtha et al. [109] can be extended to include unique structural modifications to faithfully model specific polymer chemistries, which is the focus of this article.
In this article, we use atomistic MD simulations to guide the modifications needed in this generic CG model of Kulshreshtha et al. [109] to represent two specific polymer chemistries-namely poly (4-vinylphenol) and poly (2-vinylpyridine) in tetrahydrofuran (THF). THF is expected to be a good solvent for both polymers [110][111][112][113]. We choose poly(4-vinylphenol) as an example polymer chemistry since it is capable of forming both intra-and inter-chain hydrogen bonds and previous studies have shown the role of hydrogen bonds in promoting miscibility in blends of poly (4-vinylphenol) with other hydrogen bonding polymer chemistries [114][115][116]. The second polymer poly(2-vinylpyridine) does not exhibit intra-chain hydrogen bonds but is capable of forming inter-chain hydrogen bonds with donor polymers (e.g., poly(4-vinylphenol)) [117][118][119], thus making it another suitable candidate for model development. Rather than conduct a bottom-up development for a completely new CG model using atomistic to CG mapping approaches described earlier, we want to demonstrate in an incremental step-by-step manner what few modifications (e.g., intra-chain interactions and torsional constraint) are needed in the generic CG model of Kulshreshtha et al. [109] to reproduce the atomistic structure of these polymer chemistries. First, we achieve this specifically for the chain length of 24 monomers. We then test how well this modified CG model predicts the structure for polymer chain lengths of 12, 18, and 36 monomers. We also show that this modified CG model, tailored to reproduce atomistic structure, does not reproduce atomistic chain relaxation or hydrogen bond dynamics. Additionally, to motivate the need for our modified CG model over atomistic models, we quantify the computational speed up that we achieve by using simulations with the modified CG model as compared to atomistic simulations.
This article is organized as follows. In Section 2, we describe details pertaining to the atomistic and CG models, MD simulation details, and the data analyses. In Section 3, we first present the model development for poly (4-vinylphenol) and poly (2-vinylpyridine), followed by testing the transferability of the CG model for different chain lengths and then comparing the dynamic behavior of the CG model to atomistic simulations. In Section 4, we conclude with a summary of key results and an outlook for future work.

Atomistic Molecular Dynamics Simulation
We conduct atomistic molecular dynamics (MD) simulations of a single chain of either poly (4-vinylphenol) or poly (2-vinylpyridine), abbreviated as pvpH and pvpY, respectively, in explicitly represented tetrahydrofuran (THF) molecules in an isothermal-isobaric (NPT) ensemble at a constant pressure and temperature of 1 bar and 298 K, respectively, using GROMACS 5.1.2 package [120][121][122]. For both pvpH and pvpY, THF is expected to act as a good solvent [110][111][112][113]. We consider pvpH and pvpY, comprised of 12, 18, 24, and 36 monomers, denoted as 12-mer, 18-mer, 24-mer, and 36-mer. The intra-and inter-chain interactions of pvpH, pvpY, and THF are modeled using OPLS-AA force field [123,124]. We choose OPLS-AA force field because (a) it is applicable for a wide range of organic molecules, such as organic liquids and ligands, aromatic biaryls, proteins, or nucleic acids [123][124][125][126] and (b) it has been used in many atomistic simulation studies aimed at development of CG polymer models [57,127,128]. We note however, that to the best of our knowledge, the OPLS-AA force field has not been optimized to reproduce correctly the polymer-solvent (pvpH-THF and pvpY-THF) interactions that are relevant to this paper.
We generate all the OPLS-AA parameters using LigParGen web server with 1.14*CM1A-LBCC model for assigning partial charges [125,126,129,130]. As proposed for the OPLS-AA model [123,124], the 1-2 and 1-3 non-bonded interactions are excluded while the 1-4 non-bonded interactions are reduced by a factor of 2, and geometric-mean combining rule is used for computing both energy and size, Lennard-Jones (LJ) [131] interaction parameters of unlike pairs. Analytical long-range tail corrections [132] accounting for dispersion are applied for the non-bonded LJ interactions while electrostatic interactions are handled using particle mesh Ewald (PME) method [133], using a fourth-order cubic interpolation, each with a potential cutoff of 12 Å.
We generate initial configurations using PACKMOL [134], by randomly placing a single pvpH or pvpY polymer chain of a particular chain length and 5000 THF molecules in a cubic simulation box of size 10 nm with periodic boundary conditions in x, y, and z directions. These initial configurations are subjected to the steepest descent energy minimization to remove overlaps. Then, the configuration is simulated in canonical (NVT) ensemble for a duration of 2 ns, followed by NPT ensemble MD equilibration for 10 ns, which allows the system to reach appropriate equilibrium density and potential energy, and a subsequent NPT ensemble production run for 100 ns. During the production run, the temperature and pressure are controlled using Nosé-Hoover [135,136] thermostat and Parrinello-Rahman [137] barostat with a coupling time constant of 0.4 ps and 2.0 ps, respectively. A time step of 0.001 ps is used for integrating the equations of motion using leap-frog integrator. The higher frequency bonds containing hydrogen atoms are constrained using LINCS [138] method.
For data analyses, we use the configurations obtained from the 100 ns production trajectory with coordinates saved every 10 ps. We perform, for each system, 5 independent trials with distinct initial configurations. When we report a single ensemble structural analysis value, we show the average and the standard deviations from the 50,000 total configurations from 5 trials. When we report probability distributions, we calculate the distribution from the 10,000 configurations in each trial and report the average distribution and standard deviation from 5 trials.

Coarse-Grained (CG) Model
As this paper is focused on showing the modifications that need to be made to the previously published generic CG polymer model of Kulshreshtha et al. [109] to specifically model pvpH and pvpY, we first describe the features that we inherit from the generic CG model of Kulshreshtha et al. [109] and then describe the modifications. In the generic CG model of Kulshreshtha et al. [109], each monomer Polymers 2020, 12, 2764 5 of 29 along the polymer (pvpH or pvpY) chain is represented with one CG backbone (B) bead and a CG hydrogen bonding (H) bead ( Figure 1). The B bead diameter is set to 1 d, with d serving as the reduced unit of distance; as we separately simulate single chain(s) of pvpH or pvpY in implicit solvent, this value of 1 d is equivalent to 6 Å in simulations of pvpH chain and 5.29 Å in simulations of pvpY. If one simulated a system with both these chains, 1 d would be equal to 5.29 Å and the sizes of the CG beads for the two chemistries would be scaled accordingly; the reader interested in this scenario is directed to Section I of the Supplementary Information (SI) and Figure S1. The H beads on pvpH and pvpY are included to model the hydrogen bonding interactions that one could observe between hydroxy (-OH) groups in systems with pvpH chains or between an acceptor nitrogen atom of pvpY and donor -OH group in systems involving a blend of pvpH and pvpY chains. The H bead diameter is set to 0.3 d and is placed at 0.37 d from the center of the B bead. This selection allows the H bead to be partially embedded within the B bead, exposing only a small volume of H bead to allow for effectively directional interactions, as shown in the work of Kulshreshtha et al. [109]. . Each CG backbone (B) bead represents a monomer and is shown in gray color overlaid on the atomistic monomer representation. Each CG hydrogen bonding (H) bead is shown in yellow color. Following the work of Kulshreshtha et al. [109], the position and size of H bead with respect to the B bead are set to be the values shown in (c), where d represents the reduced unit of distance and is equivalent to 6 Å and 5.29 Å for pvpH and pvpY, respectively. As done by Kulshreshtha et al. [109], the polymer chain is modeled as a bead-spring [139] chain. The bond between the monomers is represented by a harmonic potential between bonded B and B' beads (apostrophe denotes an adjacent bonded monomer), as shown in Equation (1).
Similarly, the bond between H bead and its parent B bead is modeled via a harmonic bonded B-H potential. The equilibrium bond length, r 0 , is set to 1 d for B-B and 0.37 d for B-H and the force constant, k bond , is equal to 50 and 1000 kT/d 2 for B-B and B-H, respectively. The angle potential between a H bead, its parent B bead, and the adjacent bonded B bead, denoted as H-B-B angle, with the form of Equation (2), is defined to constrain the rotation of H bead with respect to the B bead.
In Equation (2), k angle and θ 0 are set to 50 kT/radian 2 and 90 • respectively. The B-B -H angles are unrestricted, except at the last monomer bead in the chain. Moreover, it is important to mention that these B-B and B-H bonded potentials and H-B-B angle potential parameters are chosen to maintain directionality of the hydrogen bonding interactions [109].
In the early stages of the model development, as done in Kulshreshtha et al. [109], (a) we do not have 3-and 4-body bonded potentials along the backbone of the chains (i.e., no B-B -B angle potential or B-B -B -B dihedral angle potential) to mimic a flexible polymer chain, and (b) we do not have a H-B-B -H dihedral angle potential to allow free rotation of H beads along the polymer chain. In later stages of the modified CG model development, to better match the CG model chain conformations with those from atomistic results, we modify the above two choices. When needed, the dihedral constraints are incorporated into the model using the following steps. We first obtain the energy distribution by direct Boltzmann inversion of the target probability distribution functions: where C ∅ is the constant that sets the minima of the potential to zero. Based on the profile of U(∅) in Equation (3), we choose to fit a 4-term Fourier-type dihedral potential of the form and obtain dihedral coefficients, k dihedral, i , and equilibrium dihedral angles, d i .
In early stages of the modified CG model development, to mimic the THF (good) solvent implicitly in the CG model, the non-bonded B-B , B-H, and H-H interactions are modeled as isotropic and purely repulsive using Weeks-Chandler-Andersen (WCA) [140] potential described as: The pairwise non-bonded interaction parameters are set as ε BB = ε BH = ε HH = 0.1 kT and σ BB = 1 d, σ HH = 0.3 d and σ BH = 0.65 d, where σ ij is set according to the arithmetic mean diameter of the interacting bead pair. In later stages of the modified CG model development, to better match the CG model chain conformations with those from atomistic results, we introduce an attractive interaction between H-H Polymers 2020, 12, 2764 7 of 29 beads for pvpH polymer and between B-B beads for pvpY polymer. These attractive non-bonded interactions are modeled using cut and shift LJ [131] interaction, which takes the following form: The values of ε ij of the interacting bead pair are varied from 0.1 kT to higher values till we obtain a good match between the CG simulations and atomistic simulations for the target conformational property distribution.

CG MD Simulation Details
Using the CG model described above, we perform Langevin dynamics simulations using LAMMPS (August 2018 version) package [141] in NVT ensemble. The choice of simulation package, i.e., LAMMPS for CG MD simulations and GROMACS for atomistic simulations, is purely based on the ease of implementation of the chosen models in the respective packages and our results should be independent of the software package used for running MD simulations.
At the start of the simulation, we randomly place a single chain of pvpH or pvpY (12-mer, 18-mer, 24-mer, or 36-mer) in an extended rod-like configuration, with the B-B and B-H distance set to 1 d and 0.37 d, respectively, in a cubic simulation box of size 100 d, with periodic boundary conditions in x, y, and z directions. To relax the chain away from this unphysical initial configuration, we run the simulation for 10 7 time steps at temperature T* = 1 (in reduced units) using a Nosé-Hoover [135,136] thermostat, with all the non-bonded interactions, including 1-3 and 1-4 interactions that prevent intra-chain bead overlap, modeled as purely repulsive using WCA interaction potential. We note that one simulation time step is set to ∆t = 0.0001τ (in reduced units), where τ is equivalent to 4.18 ps for simulations with pvpH and 3.45 ps for simulations with pvpY (see Section I in the SI for conversion from reduced time to real time units). After the initialization stage, the system is equilibrated for another 10 7 time steps, where the non-bonded interactions are set to those in the CG model specifications and the temperature is maintained at T* = 1 using a Langevin thermostat with the damping parameter (i.e., "damp" in LAMMPS package) of 10 time steps to model the solvent effect implicitly. Our choice of this damping parameter should not impact the values of the equilibrated ensemble structural properties presented in the Results section. We tested a range of damping parameters between 10 and 100 time steps and found that the chosen value of 10 time steps allowed frictional forces due to the implicit solvent to be commensurate with conservative forces, allowing sampling of configurations in an implicit solvent environment. However, given the inverse relationship of this damping parameter and simulated viscous effects of the solvent, we expect this value to impact the CG model dynamics versus atomistic model dynamics.
The equilibration stage is followed by a production stage of 5 × 10 8 time steps, which is equivalent to 209 ns for pvpH and 172.5 ns for pvpY, during which we sample configurations every 10 5 time steps. We repeat, for each system, 10 independent trials with different initial configurations and random number seeds (used for initial velocities and damping forces in Langevin equations). When we report a single ensemble structural analysis value, we show the average and the standard deviations from the total configurations from 10 trials, and when we report probability distributions, we calculate the distribution from the 5000 configurations in each trial and report the average distribution and standard deviation from 10 trials. the center of the B bead corresponds to the center of mass of the (pvpH or pvpY) monomer and the center of the H bead corresponds to the relative positions of hydroxy (-OH) group in pvpH and nitrogen atom in pvpY monomer from the center of mass of the monomer. For CG simulations, the center of the CG B bead is mapped to the center of mass of the corresponding monomer in the atomistic simulations, whereas the center of the CG H bead is pre-defined relative to the CG B monomer, similar to the generic CG model of Kulshreshtha et al. [109].

Analyses
We quantify the polymer chain conformation and chain backbone stiffness in both atomistic and CG MD simulations. The chain conformations sampled are plotted as probability distributions of end-to-end distance (R ee ), as shown in Equation (7): where r 0 and r l are the positions of the first and last monomer beads of the chain. For comparison between atomistic and CG simulations' results, the B-B , B-H, and R ee probability distributions obtained from atomistic simulations are scaled by the average B bead diameter (6 Å for pvpH and 5.29 Å for pvpY), to convert the distributions from real units (Å) to reduced distance units (d).
We quantify the chain backbone stiffness with the persistence length (L P ) calculated using the autocorrelation function of bond vectors along the polymer chain [142]: with → b 1 being the bond vector for the first bond (from bead 0 to 1), → b i being the i th bond vector (from bead i-1 to i) of the chain, and b being the average bond length, where . . . denotes an ensemble average. L P is solved by fitting an exponential function to the autocorrelation function C(i) in Equation (8) and is also reported in reduced distance units (d).
We analyze the chain relaxation dynamics by calculating the autocorrelation function of the end-to-end vector [ACF(R ee (t))], described as: where R ee (0) is the end-to-end vector at any initial time t = 0, R ee (t) is the end-to-end vector at any time t, and . . . denotes an ensemble average. To quantify the dynamic behavior of the intra-chain hydrogen bonds in pvpH, we calculate time autocorrelation function shown in Equation (10).
In Equation (10), the variable h ij takes on the value 1 when the pair of i and j H beads are hydrogen bonded, and 0 otherwise. The subscript x in C x (t) refers to the "continuous" definition of the hydrogen bonds, i.e., a hydrogen bond once broken is considered broken for the remainder of the time, thus providing information on short-time scale behavior of hydrogen bonds. For atomistic simulations, hydrogen bonds are considered to be formed when the distance between donor and acceptor atoms is less than or equal to 3.0 Å and donor-hydrogen-acceptor angle is less than or equal to 30 • . For CG simulations, we consider a pair of CG H beads to be hydrogen bonded when they are within 1.50 σ HH (0.45 d) of each other, which is large enough to ensure that all the hydrogen bonding pairs within the first coordination shell are taken into account.

Development of CG Model Using 24-mer Chains
As our CG model is extended from the generic CG model of Kulshreshtha et al. [109], we first compare the structures generated by the CG model of Kulshreshtha et al. [109], denoted as the "original" CG model without any modification, against the atomistic simulation results of pvpH and pvpY. Figure [109], such orientational preference is lacking, as no torsional constraints are imposed. Despite these differences in local monomer-level structure between the atomistic and the original CG model, the correlation of bonds (calculated as in Equation (8)) from atomistic and CG simulations ( Figure 3d) and the values of persistence lengths (LP) of 1.69 ± 0.54 d, 1.41 ± 0.12 d, and 2.00 ± 0.02 d for atomistic simulations of 24-mer pvpH and 24-mer pvpY chains, and the 24-mer chain modeled with the original generic CG model, respectively, show reasonable agreement. We also find that the original CG of Kulshreshtha et al. [109] can correctly capture the polymer scaling exponent ~0.6 for a polymer chain in good solvent (see Section II and Figure S2 in SI) as well as the expected site has an orientational preference for both pvpH and pvpY chains. However, in the original CG model of Kulshreshtha et al. [109], such orientational preference is lacking, as no torsional constraints are imposed. Despite these differences in local monomer-level structure between the atomistic and the original CG model, the correlation of bonds (calculated as in Equation (8)) from atomistic and CG simulations (Figure 3d) and the values of persistence lengths (L P ) of 1.69 ± 0.54 d, 1.41 ± 0.12 d, and 2.00 ± 0.02 d for atomistic simulations of 24-mer pvpH and 24-mer pvpY chains, and the 24-mer chain modeled with the original generic CG model, respectively, show reasonable agreement. We also find that the original CG of Kulshreshtha et al. [109] can correctly capture the polymer scaling exponent 0.6 for a polymer chain in good solvent (see Section II and Figure S2 in SI) as well as the expected distribution of mean-squared internal distances (see Section II and Figures S3 and S4 of the SI). These put together demonstrate the power of simple generic CG polymer models to correctly capture universal polymer physics for a broad range of chain lengths [12,14,[143][144][145][146].  Interestingly, the Ree distances ( Figure 3e) sampled in atomistic simulations of pvpH and pvpY chains in THF are smaller than those sampled with the original CG model and exhibit more fluctuations with high standard deviations. These Ree distributions sampled in the atomistic simulations suggest that the selected force field parameters (that are not optimized for pvpH and pvpY with THF) are likely modeling a solvent quality that is poorer than the expected good solvent quality and that there could be some kinetic trapping of configurations in the atomistic simulations. The original CG model, based on the expected good solvent scaling behavior (see Section II and Figure S2 in SI) and the smooth Ree distribution with small standard deviations, demonstrates  Interestingly, the R ee distances (Figure 3e) sampled in atomistic simulations of pvpH and pvpY chains in THF are smaller than those sampled with the original CG model and exhibit more fluctuations with high standard deviations. These R ee distributions sampled in the atomistic simulations suggest that the selected force field parameters (that are not optimized for pvpH and pvpY with THF) are likely modeling a solvent quality that is poorer than the expected good solvent quality and that there could be some kinetic trapping of configurations in the atomistic simulations. The original CG model, based on the expected good solvent scaling behavior (see Section II and Figure S2 in SI) and the smooth R ee distribution with small standard deviations, demonstrates sufficient sampling of the equilibrated states. Assuming the atomistic simulation results to be correct, to match the atomistically observed R ee distances for 24-mer pvpY and pvpH chains, we need to modify the original CG model to include attractive interactions within the chain either due to intra-chain H-bonds or due to solvent induced B-B interactions. We note that despite similarities in the atomistic R ee distributions of pvpH and pvpY chains, for the pvpH chemistry, we see the formation of intra-chain hydrogen bonds between the hydroxy (-OH) groups (see snapshots included in the figure). We see in the atomistic simulations that as the number of intra-chain hydrogen bonds increases, the R ee distance decreases, as shown in Figure S5 of SI. These intra-chain hydrogen bonds are absent in pvpY chemistry as pvpY only has hydrogen bond accepting nitrogen atoms. Thus, for pvpY, any chain collapse is likely driven by van der Waals interactions, interactions between the pvpY and THF, and potential intra-chain π-π stacking interactions between the pyridine aromatic rings.
In summary, the atomistic simulation results at the monomer-level and chain-level for chains of pvpH and pvpY show differences between the two chemistries and some monomer-level features, unique to these chemistries. Next, we describe in a step-by-step manner the modifications that one would need to incorporate into the original generic CG model of Kulshreshtha et al. [109] to reproduce the atomistically observed structures for each polymer chemistry.

24-mer pvpH Chain
For pvpH, to capture the intra-chain hydrogen bonds seen in the atomistic simulations, in the CG model, we introduce an attractive non-bonded interaction between H-H beads and systematically vary the ε HH from 6 to 10 kT, keeping all the other parameters the same as the original CG model of Kulshreshtha et al. [109]. Anticipating this modification to alter the R ee sampled, we plot the probability distribution of R ee distance as a function of changing hydrogen bonding strength along with the reference atomistic results ( Figure S6 of SI). Use of ε HH = 6 kT barely changes the R ee from that obtained via the generic CG model of Kulshreshtha et al. [109]. As the ε HH is systematically increased, the range of R ee sampled shifts to smaller values with the ε HH = 7 kT, producing the same range of sampled R ee in the CG simulations as that in the atomistic simulations. In Figure S7 of SI, we plot the structural features of the CG pvpH chains for this case of ε HH = 7 kT versus that of the reference atomistic results. In Figure S7a, we see overall agreement between the R ee sampled by CG and atomistic simulations, but the CG simulation distribution exhibits the expected mono-peaked shape in the R ee distribution, while the atomistic simulations exhibit large error due to lesser sampling relative to CG simulations. The introduction of H-H attraction also affects the overall stiffness of the chain.  Figure S7e) also exhibits a peak at~60 • as in atomistic simulations (not seen in the original CG model in Figure 3a); however, the overall agreement with the atomistic B-B -B angle distribution is still not good. This is because the pvpH chain simulated with this modified CG model still samples higher angles than in atomistic simulations due to the higher flexibility and excluded volume effects in the CG model. Overall, we see significant improvement over the original CG model of Kulshreshtha et al. [109] in reproducing pvpH specific atomistic simulation results for 24-mer chains. We emphasize that rather than conducting a complete CG model parameterization of every bonded and non-bonded potential using one of the many methods described in the Introduction, we want to show what few modifications are needed to ensure that the generic "original CG model" better represents a specific chemistry. By incorporating a single modification-in this case, an appropriate attractive interaction between hydrogen bonding (H) beads to enable intra-chain hydrogen bonds-it is possible to reproduce most of the structural features of pvpH chains, as seen in the atomistic simulations, including chain conformations and stiffness.
To test if it is possible to further improve the performance of the modified CG model, in addition to the attractive interaction between H-H beads, we also impose H-B-B -H dihedral as a second modification to maintain the orientation of the H beads. This H-B-B -H dihedral potential is parameterized using direct Boltzmann inversion of the corresponding probability distribution function obtained from the atomistic model, as shown in Equation (3). Figure S8a of SI presents the quality of fit obtained, using Fourier-style dihedral with four terms (Equation (4)). We, again, systematically vary the attractive strength, ε HH , from 6 to 10 kT. Clearly, explicit incorporation of H-B-B -H dihedral in the CG model does not alter the R ee trends and the best match results are again obtained for ε HH = 7 kT (see Figure S9 of SI). The results of all the structural properties of this modified CG model are plotted in Figure 4 and we only highlight a few key observations next.
The snapshots of the CG model chain conformation (Figure 4a) at the average R ee demonstrate that the attractive interaction between H beads lead to the formation of intra-chain contacts and a compact conformation. The relatively small size of the H bead and its placement with respect to the B bead center forces the H-H beads to interact directionally. We note, however, that unlike Kulshreshtha et al. [109], where the authors also captured the specificity of hydrogen bonding interactions (i.e., allowing an H-bead to interact with at most one other H-bead) by judiciously choosing repulsive interactions between like acceptor H-acceptor H and donor H-donor H bead pairs, in this modified CG model of pvpH, the specificity of the hydrogen bonding interaction is lost and the formation of H-H-H interactions is now a possibility. However, the formation of such H-H-H interactions is also seen in the atomistic model, as observed in the snapshot shown in Figure 3 (see black circles around such contacts). Furthermore, the calculated value of L P from the CG model, 1.47 ± 0.14 d, and the probability distributions for B-B -B -B dihedral and B-B -B angle remain unaffected. As expected, the agreement of H-B-B -H dihedral between the CG and atomistic model is significantly enhanced (Figure 4c versus Figure S7 of SI). Therefore, depending upon the choice of structural property that needs to be reproduced correctly, one can judiciously introduce additional modifications in our CG model to facilitate future studies.

24-mer pvpY Chain
For pvpY chains, unlike pvpH, we do not need to introduce attraction between H beads because we do not have intra-chain hydrogen bonds; the nitrogen in the pvpY monomer serves only as a hydrogen bonding acceptor atom. However, in the atomistic simulations, the chain samples more collapsed conformations than the original generic CG model (Figure 3e). To produce these collapsed conformations, we introduce into the original generic CG model the first modification-a weak non-bonded attractive interaction between backbone B beads. We systematically vary the ε BB from 0.6 to 1.0 kT, while keeping all the other parameters same as the original generic CG model of Kulshreshtha et al. [109]. One would expect that upon increasing the attraction between B-B beads, pvpY should exhibit smaller chain conformations. Figure S10 of SI illustrates changes in probability distribution of R ee as a function of changing B-B attractive strength, plotted against the reference atomistic results. The range of R ee sampled by this modified CG model best matches the atomistic results for values of ε BB within 0.7 and 0.8 kT. We choose ε BB of 0.7 kT as the best case for further analysis.
To test if it is possible to further improve the performance of the modified CG model, in addition to the attractive interaction between H-H beads, we also impose H-B-B´-H´ dihedral as a second modification to maintain the orientation of the H beads. This H-B-B´-H´ dihedral potential is parameterized using direct Boltzmann inversion of the corresponding probability distribution function obtained from the atomistic model, as shown in Equation (3). Figure S8a of SI presents the quality of fit obtained, using Fourier-style dihedral with four terms (Equation (4)). We, again, systematically vary the attractive strength, , from 6 to 10 kT. Clearly, explicit incorporation of H-B-B´-H´ dihedral in the CG model does not alter the Ree trends and the best match results are again obtained for = 7 kT (see Figure S9 of SI). The results of all the structural properties of this modified CG model are plotted in Figure 4 and we only highlight a few key observations next.  In Figure S11 of SI, we plot and compare structural features such as R ee distance, stiffness of the chain , H-B-B -H and B-B -B -B dihedrals, and B-B -B angle obtained for the simulations with the modified CG model (ε BB = 0.7 kT) and the reference atomistic simulations. We see good agreement between the R ee sampled by the modified CG model and the atomistic simulations ( Figure S11a). The calculated value of L P = 1.74 ± 0.02 d with the modified CG model better agrees with the atomistic model L P of 1.41 ± 0.12 d than the original CG model (2.0 ± 0.02 d). The H-B-B -H dihedral distribution with the modified CG model ( Figure S11c) is effectively the same as the original CG model (Figure 3c) as we have not introduced any attractive interaction between H-H beads or imposed a H-B-B -H dihedral. Moreover, weak attraction between B beads (ε BB = 0.7 kT) does not improve agreement between the probability distribution of B-B -B -B dihedral angles ( Figure S11d) and the corresponding atomistic distribution; however, we note the emergence of slight bimodality in Figure S11d that is absent in the original CG model results (Figure 3b). As discussed in the section on pvpH model development, the B-B -B angle distribution does not agree with the atomistic distribution well (Figure S11e) due to chain flexibility and excluded volume effects of the chain in the CG model. These results suggest that this modified CG model needs one or more additional modification(s).
As a second modification, in addition to attractive interaction between B-B beads, we introduce the B-B -B -B dihedral angle constraint. The dihedral potential is parameterized using direct Boltzmann inversion and Figure S8b of SI presents the quality of fit obtained. We, again, systematically vary the attractive strength, ε BB , from 0.6 to 1.0 kT and find that R ee trends remain unaltered and the best match between the atomistic and the modified CG model results is again obtained for ε BB = 0.7 kT. We present all the structural properties in Figure S12 of SI. As expected, the agreement in B-B -B -B dihedral distributions from CG and atomistic is significantly improved; however, the probability distributions/average values of R ee distance, B-B -B angle, H-B-B -H dihedral, and chain stiffness (L P = 1.72 ± 0.01 d) remain mostly unaffected as compared to Figure S11, suggesting that the overall performance of the CG model is only marginally improved.
To further improve the orientation of the H beads in the CG model, in addition to attractive interaction between B-B beads and B-B -B -B dihedral angle constraint, we also impose H-B-B -H dihedral as a third modification and systematically re-evaluate the R ee trends. Interestingly, the incorporation of H-B-B -H dihedral into the CG model does not alter the R ee trends and the best match results are again obtained for ε BB = 0.7 kT. The results of all the other structural properties are plotted in Figure 5. As expected, the agreement in H-B-B -H dihedral distributions from CG and atomistic improves; however, the probability distributions/average values of R ee distance, B-B -B angle, B-B -B -B dihedral, and chain stiffness (L P = 1.72 ± 0.01 d) do not change at all, suggesting that the overall performance of the CG model is slightly improved over the second modification ( Figure S12 of SI). We further note that a higher value of ε BB , e.g., ε BB = 0.9 kT (1.0 kT), for which the corresponding R ee distribution does not agree with that from atomistic simulations ( Figure S10 of SI), can lead to the L P value of 1.66 ± 0.01 d (1.63 ± 0.01 d), which is closer to the atomistic value. Nonetheless, this modified CG model of pvpY produces atomistically observed chain conformations and stiffness (R ee and L P ) within 10% deviation.

Testing the Transferability of the CG Model for Describing Structural Properties at Other Chain Lengths
To test the chain length transferability of the modified CG model developed in the previous section for 24-mers, in describing the structural properties, we conduct atomistic and CG model simulations for single pvpH and pvpY chains with chain lengths that are longer (36-mer) and shorter (18-mer and 12-mer) than the 24-mer chain used for determining the CG model parameters. We remind the reader of the few modifications made to the original generic CG model of Kulshreshtha et al. that the overall performance of the CG model is slightly improved over the second modification ( Figure S12 of SI). We further note that a higher value of , e.g., = 0.9 kT (1.0 kT), for which the corresponding Ree distribution does not agree with that from atomistic simulations ( Figure S10 of SI), can lead to the LP value of 1.66 ± 0.01 d (1.63 ± 0.01 d), which is closer to the atomistic value. Nonetheless, this modified CG model of pvpY produces atomistically observed chain conformations and stiffness (Ree and LP) within 10% deviation. In Figure 6, we plot and compare the R ee distance of 36-mer, 18-mer, and 12-mer pvpH chains obtained from the simulations with the modified CG model and the reference atomistic simulations. We see excellent agreement between the R ee sampled by the modified CG model and the atomistic simulations. We note, however, that smaller chain lengths have the tendency to sample collapsed conformations. For 36-mer pvpH ( Figure S13 of SI) and 18-mer pvpH ( Figure S14 of SI), there is excellent agreement between other structural properties obtained from best modified CG model and atomistic model simulations. The calculated value of L P from the CG model, 1.54 ± 0.12 d (36-mer) and 1.63 ± 0.12 d (18-mer), agrees with the atomistic model L P of 1.47 ± 0.23 d (36-mer) and 2.18 ± 0.54 d (18-mer). For 12-mer pvpH ( Figure S15 of SI), although the chain conformation and stiffness are in good agreement, including the calculated value of L P (1.53 ± 0.13 d) from the CG model that agrees with the atomistic model L P of 1.48 ± 0.39 d, the local monomer-level agreement is slightly reduced from that seen for 36-mer and 18-mer pvpH. To analyze if the performance of our CG model worsens with a reduction in the chain length, we also perform atomistic and CG model simulations for single 10-mer pvpH chains. For the 10-mer pvpH chain ( Figure S16 of SI), although the range of R ee sampled by CG simulations is similar to that in the atomistic simulations, the chains sample more collapsed conformations in the CG simulations than the atomistic simulations. Furthermore, the calculated value of L P , 1.24 ± 0.08 d for the 10-mer, also does not agree well with the corresponding atomistic L P of 2.17 ± 0.57 d. This means that either the chosen value of ε HH that worked well for 24-mer, 36-mer, 18-mer, and 12-mer chains needs to be refined in the CG model for 10-mer or that the atomistic simulations of 10-mer chains are not as prone to the kinetically trapped collapsed configurations as the longer chains are. Similarly, for pvpY chains, Figure 7 shows that the modified CG model can reproduce the atomistically observed Ree distance for 36-mer, 18-mer, and 12-mer chains. We also note that, unlike 36-mer, the small differences in the range of Ree sampled by 18-mer and 12-mer pvpY chains are similar to those observed in 24-mer pvpY chains (see Figure 5a). In Figures S17-S19 (18-mer), and 1.58 ± 0.33 d (12-mer). As discussed for pvpH chains, further decrease in the chain length leads to disagreement between results from CG and atomistic simulations. For 10-mer pvpY, as shown in Figure S20   Similarly, for pvpY chains, Figure 7 shows that the modified CG model can reproduce the atomistically observed R ee distance for 36-mer, 18-mer, and 12-mer chains. We also note that, unlike 36-mer, the small differences in the range of R ee sampled by 18-mer and 12-mer pvpY chains are similar to those observed in 24-mer pvpY chains (see Figure 5a). In Figures S17-S19 of SI, we plot and compare all the other structural features for 36-mer, 18-mer, and 12-mer pvpY chains obtained from atomistic and CG simulations, respectively. The figures show that, similar to 24-mer, there is good agreement for monomer-level structure, chain conformations, and chain stiffness, including the calculated values of L P 1.73 ± 0.02 d (36-mer), 1.73 ± 0.01 d (18-mer), and 1.72 ± 0.01 d (12-mer) from the CG model, which statistically agree with the atomistic model L P of 1.91 ± 0.16 d (36-mer), 1.91 ± 0.42 d (18-mer), and 1.58 ± 0.33 d (12-mer). As discussed for pvpH chains, further decrease in the chain length leads to disagreement between results from CG and atomistic simulations. For 10-mer pvpY, as shown in Figure S20 of SI, the calculated L P , 1.72 ± 0.01 d, is within the standard deviation of the atomistic model L P of 1.59 ± 0.26 d, but the distribution of the R ee sampled by CG simulations is much broader than that of the atomistic simulations. For significantly shorter chain lengths, one could either conduct atomistic simulations that are computationally less expensive than longer chains or consider developing and refining the parameters of the modified CG model for the particular oligomer of interest. We also conduct simulations with the modified CG model for pvpH and pvpY for a broader range of chain lengths (12-mer to 500-mer) to establish the polymer scaling exponent and the implicit solvent quality. Our calculation of the scaling exponent shows that the modified CG model for pvpH samples nearly good solvent conformations (see Figure S21a in the SI). The modified CG model for the pvpY chain suggests a poor solvent scaling exponent ( Figures S21b in the SI). We refrain from placing too much emphasis on the corresponding scaling exponents and mean-squared internal distances obtained from atomistic simulations over a small range of chain lengths (12-mer to 36-mer). However, representative atomistic simulation snapshots ( Figure S22) show that THF molecules form explicit hydrogen bonds with pvpH but do not form such hydrogen bonds with pvpY; the latter leads to more collapsed pvpY conformations and poor solvent behavior in atomistic simulations. As the modified CG model for pvpY chains is developed to reproduce atomistically observed chain conformations, the modified CG model, not surprisingly, exhibits poor solvent scaling.

Chain Relaxation and Hydrogen Bonding Dynamics: CG Model versus Atomistic Model
We expect that the modified CG model developed to reproduce atomistically observed (a) 36  We also conduct simulations with the modified CG model for pvpH and pvpY for a broader range of chain lengths (12-mer to 500-mer) to establish the polymer scaling exponent and the implicit solvent quality. Our calculation of the scaling exponent shows that the modified CG model for pvpH samples nearly good solvent conformations (see Figure S21a in the SI). The modified CG model for the pvpY chain suggests a poor solvent scaling exponent ( Figure S21b in the SI). We refrain from placing too much emphasis on the corresponding scaling exponents and mean-squared internal distances obtained from atomistic simulations over a small range of chain lengths (12-mer to 36-mer). However, representative atomistic simulation snapshots ( Figure S22) show that THF molecules form explicit hydrogen bonds with pvpH but do not form such hydrogen bonds with pvpY; the latter leads to more collapsed pvpY conformations and poor solvent behavior in atomistic simulations. As the modified CG model for pvpY chains is developed to reproduce atomistically observed chain conformations, the modified CG model, not surprisingly, exhibits poor solvent scaling.

Chain Relaxation and Hydrogen Bonding Dynamics: CG Model versus Atomistic Model
We expect that the modified CG model developed to reproduce atomistically observed structures need not replicate the atomistically observed chain relaxation or hydrogen bond decorrelation dynamics. Nonetheless, to understand how the chain relaxation for both pvpH and pvpY with the modified CG polymer model (in implicit solvent) compares to that of the atomistic polymer model (in explicit solvent), we compare end-to-end vector autocorrelation functions (ACFs) for both atomistic and CG simulations for 12-mer, 18-mer, 24-mer, and 36-mer (Figure 8). For atomistic simulations, considering the average correlations and their standard deviations, evidently, all the end-to-end vector ACFs decay to zero, which indicates chain relaxation. However, slower decay (on an average) and rising standard deviations as the chain length increases show how difficult it is to relax chain conformations using unbiased atomistic simulations. An overlap of the standard deviations for all chain lengths (except 12-mer) suggests no apparent chain length dependent trend ( Figure S23 of SI displays the zoomed-in correlation decay within 0-20 ns). Unlike the atomistic model, the (structurally best performing) modified CG model in implicit solvent exhibits the expected qualitative trend of relaxation time increasing with increasing chain length. Furthermore, as expected, the CG chain relaxation times for the chosen Langevin dynamics damping parameter are significantly faster as compared to atomistic simulations of the polymer chain in explicit solvent. This is expected of CG models and has been discussed in prior review articles on CG modeling of polymer dynamics [9,39,40,147]. Figure 8 simply confirms that the best performing CG model, structurally, is not optimized to reproduce the chain relaxation dynamics seen with atomistic For atomistic simulations, considering the average correlations and their standard deviations, evidently, all the end-to-end vector ACFs decay to zero, which indicates chain relaxation. However, slower decay (on an average) and rising standard deviations as the chain length increases show how difficult it is to relax chain conformations using unbiased atomistic simulations. An overlap of the standard deviations for all chain lengths (except 12-mer) suggests no apparent chain length dependent trend ( Figure S23 of SI displays the zoomed-in correlation decay within 0-20 ns). Unlike the atomistic model, the (structurally best performing) modified CG model in implicit solvent exhibits the expected qualitative trend of relaxation time increasing with increasing chain length. Furthermore, as expected, the CG chain relaxation times for the chosen Langevin dynamics damping parameter are significantly faster as compared to atomistic simulations of the polymer chain in explicit solvent. This is expected of CG models and has been discussed in prior review articles on CG modeling of polymer dynamics [9,39,40]. Figure 8 simply confirms that the best performing CG model, structurally, is not optimized to reproduce the chain relaxation dynamics seen with atomistic models. We note, however, that the Langevin dynamics damping parameter, used in CG simulations to model the solvent effect implicitly, can be adjusted to tune this dynamic behavior.
For pvpH chains, as intra-chain hydrogen bonds are present, we also analyze the hydrogen bond decorrelation dynamics. The behavior of C cont (t) (as described in Section 2.4 in Equation (10)) is strongly affected by short time fluctuations due to fast motion of molecules, as it requires the continuous presence of hydrogen bonds. Thus, to calculate C cont (t), we run three independent short simulations, for both atomistic and CG, generate a 100 ps trajectory, and sample configurations at a resolution of 0.01 ps for each system. For CG simulations, real-time units are converted to time steps (in reduced units) using the calculations shown in SI. We plot the C cont (t) for 12-mer, 18-mer, 24-mer, and 36-mer chains obtained from atomistic and CG simulations in Figure 9. Although there is no apparent chain length dependent trend, markedly, all the autocorrelation functions in CG simulations decay slower than atomistic simulations, indicating that short-time scale hydrogen bond dynamics in CG simulations are slower than atomistic simulations. At the classical limit, in atomistic simulations with non-polarizable force fields, a hydrogen bond is a non-bonded electrostatic interaction between fixed atomic charges (represented as partial charges on donor and acceptor atoms) with directionality achieved via atomistic bond vectors around the donor and acceptor atoms. In contrast, our CG model captures this directional hydrogen bond attraction effectively using isotropic, attractive, van der Waals-type interaction potentials and through the choice of the small-sized CG H bead with limited exposure on the CG B bead. To achieve the same structural effects as in atomistic simulations, the choice of the H-H attraction in the CG model is optimized (as discussed in Section 3.1.1). We do not expect to achieve the same hydrogen bond dynamics. However, the choice of the distance cutoff in the CG model's H-H attraction potential to qualify an hydrogen bond to be formed or broken in the CG model does impact C cont (t); Figure S24 of SI shows that as the cutoff distance in the CG model's H-H attraction potential decreases, the C cont (t) from the CG simulations decay faster and become more comparable to atomistic simulations.
Polymers 2020, 12, x FOR PEER REVIEW 19 of 30 (in reduced units) using the calculations shown in SI. We plot the Ccont(t) for 12-mer, 18-mer, 24-mer, and 36-mer chains obtained from atomistic and CG simulations in Figure 9. Although there is no apparent chain length dependent trend, markedly, all the autocorrelation functions in CG simulations decay slower than atomistic simulations, indicating that short-time scale hydrogen bond dynamics in CG simulations are slower than atomistic simulations. At the classical limit, in atomistic simulations with non-polarizable force fields, a hydrogen bond is a non-bonded electrostatic interaction between fixed atomic charges (represented as partial charges on donor and acceptor atoms) with directionality achieved via atomistic bond vectors around the donor and acceptor atoms.
In contrast, our CG model captures this directional hydrogen bond attraction effectively using isotropic, attractive, van der Waals-type interaction potentials and through the choice of the smallsized CG H bead with limited exposure on the CG B bead. To achieve the same structural effects as in atomistic simulations, the choice of the H-H attraction in the CG model is optimized (as discussed in Section 3.1.1). We do not expect to achieve the same hydrogen bond dynamics. However, the choice of the distance cutoff in the CG model's H-H attraction potential to qualify an hydrogen bond to be formed or broken in the CG model does impact Ccont(t); Figure S24 of SI shows that as the cutoff distance in the CG model's H-H attraction potential decreases, the Ccont(t) from the CG simulations decay faster and become more comparable to atomistic simulations. Overall, we conclude that the CG model is not intended to and not optimized to reproduce the inherent local (hydrogen bonding) and chain relaxation dynamics of the atomistic models.

Computational Performance
CG simulations have an inherent speedup advantage over atomistic simulations, as it allows sampling of large conformational space more efficiently within the same simulation time. We Overall, we conclude that the CG model is not intended to and not optimized to reproduce the inherent local (hydrogen bonding) and chain relaxation dynamics of the atomistic models.

Computational Performance
CG simulations have an inherent speedup advantage over atomistic simulations, as it allows sampling of large conformational space more efficiently within the same simulation time. We highlight the speedup in Table 1, where we present the average wall time (simulation run time) for pvpH chain simulations with the atomistic and structurally best performing CG model. We run short independent simulations for 10,000 time steps, with a time step of 0.001 ps in real units, for each system. For CG simulations, we convert the real time units to time steps (in reduced units), using the calculations shown in SI. Each simulation is run on 1 CPU core of 2 ×18C Intel E5-2695 v4 cluster-the standard architectural node on a local cluster [147]. From these results, we find that simulations with the CG model with implicit solvent are 30 to 100 times faster than the atomistic simulations with explicit solvent, with speedup smaller for larger polymer chain lengths. Even though these results are only shown for pvpH, as the total numbers of pairwise interactions evaluated are similar for both pvpH and pvpY chains in our study, the results presented here hold for pvpY chains as well.

Conclusions
In this article, we have described the few modifications one needs to make to a recently published generic coarse-grained (CG) model of Kulshreshtha et al. [109] that captures directionally interacting polymers, in order to represent specific polymer chemistries that exhibit hydrogen bonding. The generic CG model represents directionally interacting polymer chains using two types of beads-backbone (B) beads representing the center of mass of the monomer and hydrogen bonding (H) beads denoting the atom/groups of atoms that participate(s) in hydrogen bonding. In this article, we use explicit solvent atomistic molecular dynamics (MD) simulations as reference, at ambient temperature and pressure, to modify the generic CG model of Kulshreshtha et al. [109] to mimic the structural behavior of 24-mer poly(4-vinylphenol) (pvpH) and poly(2-vinylpyridine) (pvpY) homopolymers in tetrahydrofuran (THF). We then test the transferability of the modified CG model to other chain lengths ranging from 12-mer to 36-mer. For pvpH, an agreement between the atomistic and CG (monomer-level to chain conformations) structures is observed by introducing an attractive interaction between CG H beads and by maintaining the orientation of H beads via H-B-B -H torsional constraint in the generic CG model of Kulshreshtha et al. [109]. For pvpY, multiple modifications are required, including an attractive interaction between CG B beads, B-B -B -B , and H-B-B -H torsional constraints to capture atomistic chain conformations and monomer-level structural arrangements. For both pvpH and pvpY, our CG model is transferable to chain lengths that are slightly smaller or larger than 24-mer.
The modified CG model is expected to reproduce atomistic monomer-level structure, chain conformations, and chain stiffness, but not the atomistically observed dynamics. We find, not surprisingly, that the CG model of the polymer chain in implicit solvent exhibits significantly accelerated chain relaxation dynamics as compared to atomistic simulations of the polymer chain in explicit solvent. We also find that the short-time scale behavior of hydrogen bonds in our CG model is slower as compared to atomistic simulations, with no chain-length dependent trend. Decreasing the cutoff distance used to qualify two H-beads to be considered as hydrogen bonded improves the agreement between the CG model and atomistic simulation. Cumulatively, this suggests that our CG model is not optimized to reproduce the inherent dynamics of the atomistic model, and in order to better reconstruct the dynamics using the CG model, one has to judiciously transform the CG time scale to atomistic time scale. Similarly, we have not tested the capability of our modified CG model to predict interfacial properties (e.g., surface tension) and other thermodynamic properties (e.g., density versus pressure relationship), which may require further optimization of the CG model parameters.
In terms of computational speedup, the MD simulations using the modified CG models developed in this work reproduce the atomistic monomer-level and chain-level structures of 12-36-mer pvpH and pvpY chains 30 to 100 times faster than corresponding atomistic simulations. Such computational speedup strongly motivates the development of such simple CG models for other chemistries as well. We note, however, that the steps described in this paper ensure that the modified CG model is only as good as the atomistic model and simulation. We noted some of the potential drawbacks of the atomistic model and (unbiased) simulations, including the force field not being optimized for the polymer-solvent interactions and the possibility of sampling kinetically trapped configurations. These issues with atomistic simulations can manifest in the modified CG model. Nevertheless, if atomistic force fields for other hydrogen bonding polymer chemistries and solvents are readily available and correct, the extension of our CG model for these other polymer chemistries will be easy if one follows the steps described in this paper. This will open up the possibility of modeling large numbers of polymer chemistries that have hydrogen bonding or other directional interactions with relatively small effort. = 2kT/rad 2 ) chains in a good solvent. The average and standard deviations are computed from the means of 10 independent trials for chain lengths of 24-mer and 36-mer, and 3 independent trials for chain lengths of 100-mer and 500-mer. Note that standard deviations wherever not visible are smaller than the symbol size; Figure S5: (a) Average number of intra-molecular hydrogen bonds, and (b) end-to-end distance (R ee ), and as a function of simulation time for 24-mer pvpH chain obtained from atomistic simulation. Note that trajectories obtained from all 5 atomistic simulation trials are combined for this analysis; Figure S6: Comparison of distribution of R ee distance for 24-mer pvpH chains, obtained from both atomistic and CG simulations, as the attractive interaction between any two hydrogen bonding beads (H-H), ε HH , in the CG model is systematically increased from 6kT -10kT. The standard deviations are computed from the means of 5 independent trials for atomistic simulations and 10 independent trials for CG simulations and the lines joining the symbols are drawn to guide the eye; Figure  obtained from atomistic simulations and CG simulations (for the best performing case of attractive interaction between any two hydrogen bonding beads (H-H), ε HH , in the CG model equal to 7kT). The standard deviations are computed from the means of 5 independent trials for atomistic simulations and 10 independent trials for CG simulations and the lines joining the symbols in parts (a), (c), (d), (e) are drawn to guide the eye; Figure S8: Parameterization of (a) H-B-B -H dihedral potential for pvpH and (b) B-B -B -B dihedral potential for pvpY, for 24-mer chains, obtained by taking the direct Boltzmann inversion of the target probability distributions of the atomistic simulations. Symbols show atomistic data, while solid lines show Fourier style dihedral fit using 4 terms; Figure S9: Comparison of distribution of R ee distance for 24-mer pvpH chains, obtained from both atomistic and CG simulations, as the attractive interaction between any two hydrogen bonding beads (H-H), ε HH , in the CG model is systematically increased from 6kT -10kT and H-B-B -H torsional constraint is simultaneously imposed. The standard deviations are computed from the means of 5 independent trials for atomistic simulations and 10 independent trials for CG simulations and the lines joining the symbols are drawn to guide the eye; Figure S10: Comparison of distribution of R ee distance for 24-mer pvpY chains, obtained from both atomistic and CG simulations, as the attractive interaction between any two backbone beads (B-B), ε BB , in the CG model is systematically increased from 0.6kT -1.0kT. The standard deviations are computed from the means of 5 independent trials for atomistic simulations and 10 independent trials for CG simulations and the lines joining the symbols are drawn to guide the eye; Figure S11:  Figure S21: Scaling behavior between end-to-end distance and chain length for (a) both atomistic and CG simulations (with best performing parameters) for pvpH chains, and (b) both atomistic and CG simulations (with best performing parameters) for pvpY chains. The symbols denote actual values and lines denote the fits. The figure also includes R 2 values to access the quality of the fits. For atomistic simulations, the average and standard deviations are computed from the means of 5 independent trials. For CG simulations, the average and standard deviations are computed from the means of 10 independent trials for chain lengths of 12-mer, 18-mer, 24-mer, 36-mer and 48-mer and 3 independent trials for chain lengths of 100-mer and 500-mer; Figure S22: Mean-squared internal distances from atomistic simulations for (a) pvpH, and (b) pvpY. The average autocorrelation functions and standard deviations (shown as shaded region) are computed from the means of 5 independent trials. Also included in (c) and (d) are instantaneous atomistic simulation snapshots for 24-mer pvpH, and 24-mer pvpY in explicit THF solvent, respectively. The snapshots are color-coded to reflect the THF solvent in first solvation shell as silver, and hydrogen bonds as red lines. Part (c) and (d) show that THF molecules make explicit hydrogen bonds with pvpH but do not form such hydrogen bonds with pvpY; Figure S23: Comparison of end-to-end vector autocorrelation functions obtained from atomistic simulations for (a) pvpH and (b) pvpY chains, of chain length 12-mer, 18-mer, 24-mer and 36-mer. The average autocorrelation functions and standard deviations (shown as shaded region) are computed from the means of 5 independent trials; Figure S24: Comparison of continuous hydrogen bonding time autocorrelation function as a function of changing hydrogen bond distance criteria, obtained from CG simulations (for the best performing case of attractive interaction between any two hydrogen bonding beads (H-H), ε HH , in the CG model equal to 7kT and H-B-B -H torsional constraint simultaneously imposed) for pvpH chains, of chain length (a) 36-mer, (b) 24-mer, (c) 18-mer and (d) 12-mer. The average autocorrelation functions and standard deviations (shown as shaded region) are computed from the means of 3 independent trials. Also shown in (e) is H-H interaction potential, of the form 6-12 LJ used in the CG simulation, as a function of distance r. Note that σ HH in the CG model is equal to 0.3d.