All-Atom Model of Atactic 2-Vinyl Pyridine Polymer: Structural Properties Investigation by Molecular Dynamics Simulations

Atactic poly(2-vinyl pyridine) (P2VP) is widely used in several 3D-Printing applications, in blends or, as one of the blocks, in copolymers. Moreover, several applications have been designed by exploiting the stimuli response to the pH shown by P2VP. In this paper we propose an all atom model of P2VP, based on the well-known OPLS-AA force field, which ensures wide compatibility to model complex mixtures and/or interfaces involving P2VP in composite materials. The proposed all-atom model was checked in the reproduction of structural properties and compared with experimental data. Good reproductions of mass density and X-ray scattering pattern confirm the accuracy of the proposed model. © The Author(s) 2019. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0471909jes]

Poly(2-vinyl pyridine) (P2VP) is a polymer characterized by multiple 2-vinylpyridine units. The presence of the nitrogen atom on the aromatic ring provides peculiar physico-chemical features compared to similar polymers, such as poly(styrene) (PS), in terms of morphology and polar character (of the pyridine ring) that allows, for example, to interact with other polar molecules or to act as metal ligand. 1 This peculiar characteristic has opened many fields of application for this polymer. As example, P2VP-ABS (acrylonitrile-butadienestyrene) blend (with about 12% w/w of ABS) has been successfully used in the fabrication by 3D printing of pH responsive 3D hydrogels for microfluidic applications. In particular, P2VP has been used in the production of filaments that exhibit a coil-globular transition occurring at pH lower than 4. As further application, polyethylene glycol (PEO)-P2VP block-copolymers have been used as stimuli responsive material, induced by pH, to form hydrophilic vesicles in solution. 2 In the PS-b-P2VP block co-polymer, the capability of the pyridine group to form hydrogen bonds was used to control the supramolecular assembly of hydroxylated gold nanoparticle. 3 More in general, the stimuli-response behavior exhibited by P2VP material, mainly when it is in co-polymer, make it a good candidate for these new "smart" materials for the so-called "4D printing". [4][5][6][7] The P2VP-PS block co-polymers, characterized by the coexistence of polar and non-polar phases, opened new possibilities for ink-free lithography transfer applications. 8 In particular, the polar nature of pyridyl group allows the P2VP units to be functionalized while the PS units serve as rigid glassy scaffolds. P2VP is also used to build binary polymer brushes (usually with PS), to control and tune adhesion and friction behavior of solid surfaces. 9 Metal nanoparticle arrays 10,11 are efficiently obtained from block copolymer micelle lithography, in which P2VP is employed. 12 Morphological behavior of thin films of P2VP block co-polymers were studied in the context of fabrication of nanoporous membranes, lithography and nanophotonics. [13][14][15] Experimental studies and characterization of P2VP in terms of its synthesis, reactivity, formation of complexes, P2VP-based catalysis and photothermal reactivity, are reported in literature. 1,16 Despite several applications of P2VP in the technological field, computational studies, that could allow to understand in detail the P2VP features and to improve its performance, are not huge in number. To the best of our knowledge, only two atomistic models of P2VP have been reported in literature, a united-atom reported by He, 17 and an all-atom model reported by Soldera which employed Polymer Consistent Force Field (PCFF) 18 to investigate the reproduction of glass transition temperature. 19 However, none of these studies is specifically designed for P2VP, indeed an exhaustive investigation on the main P2VP structural properties and related comparison with experimental data is missing. Also for coarse-grained models, only few studies have been found. 20 In this study, we propose an all-atom model of P2VP, based on the OPLS-AA. 21 We have chosen to base our model on the well-known OPLS-AA force field because of its high adaptability to include and combine other organic molecules, 21-23 polymer models, 24-29 nanoparticles and other possible interfaces 30,31 which represents the common employment of the P2VP technological applications. In particular, equilibrated melt configurations were successfully obtained by applying a procedure, based on the hybrid particle-field molecular dynamics (PF-MD) approach which is described in our previous works. 32 A short introduction of the PF-MD approach is reported in the Computational Method and Model section. Then, structural properties compared with available experiments, are reported and discussed for proposed model. The Conclusion Section summarizes the main results.

Computational Method and Model
Hybrid PF-MD approach.-The hybrid PF-MD approach employed in this work has been widely validated and successfully employed in many previous works 25,[33][34][35][36] to relax polymer melt systems of both atomistic and coarse-grained models, and to set-up models of polymer composite materials. Only a brief description of such method will be given here, while the full description can be found in Ref. 33. The core feature of PF-MD method is that the non-bonded forces calculation (most computationally expensive term in a MD simulation), is obtained by the evaluation of an external potential depending on the local density at position r. 32 In the spirit of self-consistent field theory (SCFT), a many-body problem can be reduced to a problem of deriving the partition function of a single particle in an external potential V(r). 32 The non-bonded forces acting on particles can be obtained from a suitable expression of V(r) and its derivatives. In particular, in the framework of the SCF theory, the interaction of a particle with its surrounding is not calculated by pair-wise interactions but through a mean field modeled by an external potential. Assuming this, the density dependent interaction  potential W can be written as: where φ K (r) is the coarse-grained density of a particle type K at position r, and χ KK are the mean field parameters describing the interaction between a particle of type K and the density field due to the particles of type K'. It can be shown, by using the saddle point approximation, that the external potential V K (r) can be written as: To connect the particles and the field for the hybrid PF-MD scheme, it is central to obtain a smooth coarse-grained density function φ K (r) straightforwardly from the position of the particles. For this purpose, a mesh-based approach was used to obtain both the coarse-grained density function and its density derivatives (needed to calculate the forces acting on the particles). The derivation of the Eq. 2 and the details of the implementation of the PF-MD approach are reported in Refs. 32,37. Atomistic Model of Atactic P2VP. The bonded and non-bonded interaction parameters for the atomistic model of the atactic P2VP are listed in Tables II-VI. A schematic picture of the P2VP chemical structure along with the employed atom types and charges scheme is reported in Figure 1. All the parameters used for pyridine rings were taken from the OPLS-AA force field, 21 which is optimized to fit experimental properties of liquids (such as density and heat of vaporization). Instead, for the aliphatic carbons (backbone) and hydrogens, the set of parameters reported for atactic polystyrene by Müller-Plathe 38 was used.
Our choice is based on the wide range of applicability of OPLS-AA and the possibility to easily combine its force field parameters to study complex polymer systems. In fact, OPLS-AA was successfully adopted to study different polymer families, such as: polyacrylamides, 27 polyaniline, 28 polyglutamine, 39 macrocyclic polyketidies, 40 polystyrene, 38 polymer/carbon nanotube interfaces, 41 and long hydrocarbons. 23,29 According to the OPLS-AA force-field, the bond stretching between two atoms and the angle bending, are both represented by a harmonic potential (Equations 3,4): where K bond and K angle are the bond and angle force field constants, respectively, while r 0 and θ 0 are the bond length and angle equilibrium values, respectively. The full list of parameters is reported in Table II  and Table III. All the dihedral angles ψ involving the heavy atoms of the pyridine ring were represented with the Ryckaert-Bellemans (RB) potential function (Eq. 5), where ψ = (θ -180). The RB torsion parameters are reported in Table IV. V rb θ ijkl = 5 n=0 C n (cos (ψ)) n [5] All the other proper and improper dihedral angles ϕ, were represented by Eq. 6, and the full list of dihedral parameters is reported in Table V.
According to the PF-MD approach described above, the nonbonded interactions were represented by the interaction between a particle and the external field depending on the density (Eq. 2). In the procedure to relax the polymer melt, described in the following Results and discussion section, the final stage consists in the reintroduction of the non-bonded pair interactions, that in the OPLS-AA force field are described by the Lennard-Jones (LJ) and Coulombic potentials. Because the proposed model, based on the OPLS-AA, is independent from the procedure used to obtain a suitable and well relaxed initial configuration, in the following the LJ and Coulombic are reported (Eq. 7). [7] where i and j are the indexes for two different atoms, while σ ij = (σ i + σ j )/2 and ε ij = (ε ii ε jj ) 1 2 are the LJ evaluated by considering the arithmetic and geometric average, respectively. In Table VI  Simulation details.-All the PF-MD simulations were performed with the OCCAM software package. 42 The simulations were performed in the NVT ensemble at 293 K, by keeping the temperature constant through the Andersen thermostat. 43 A time step of 1 fs was used all the simulations. According with the PF-MD approach, the smooth coarse-grained density function φ K (r), that is crucial to evaluate the external potential V K (r) (see Eq. 2), is calculated using a mesh-based approach. In particular, as function of the mesh size l, different levels of coarsening of the density were obtained. In this work we used three different grid sizes: 0.8, 0.4 and 0.2 nm, in the relaxation procedure of P2VP10 melt, while, for P2VP48 system an additional grid size of 2.5 nm was employed. A scheme illustrating the mesh approach to calculate the density function is reported in the Figure 2. The MD simulations were performed using the GROMACS package 44 in the NVT and then in the NPT ensemble for 5 and 200 ns, respectively. The temperature was held constant at 293 K by using a velocity rescale algorithm 45 with a coupling constant τ T = 0.02 ps while the pressure was held constant at 1.013 bar using the Berendsen algorithm 46 with a coupling constant τ P = 0.2 ps. A time step of 2 fs was used for all the simulations. A cutoff distance of 1 nm was employed for both van der Waals and Coulomb interactions. The nonbonded interactions were excluded between first and third neighbors. In addition, non-bonded interactions among atoms within an aromatic ring were excluded. The LINCS constraint algorithm 47 was employed to fix all the distances involving hydrogen atoms and between the carbons of the backbone. The constraint between backbone carbons was removed in the simulations with the original OPLS-AA force field.

Results and Discussion
Relaxation of atactic P2VP melt with PF-MD method.-It is well known that the equilibration of polymer melts is a challenging task, even for low molecular weights. [48][49][50] In addition to such difficulties, a further technical problem of MD simulations is connected to setup and obtain a suitable initial set of coordinates of polymer melt system. This is mainly due to the large number of atom overlaps, at typical melt density, which causes divergence of force calculations in a MD simulation. In order to avoid such technical problem, a well-tested and validated procedure based on PF-MD 33 was applied in this study. In particular, the procedure adopted for two different molecular weights (10 and 48) for PV2P is the following: I) The initial configurations where prepared by randomly distributing P2VP chains in the simulation box, with two molecular weights (10 and 48 repeating units, P2VP10 and P2VP48, respectively). The initial density of the systems was set equal to the experimental density at 293 K; II) For the P2VP10 case, three sequential relaxations runs were performed at 293 K. In particular, a first relaxation, 5 ns long, was performed using a grid size l = 0.8 nm, a value that is close to the radius of gyration (Rg) of the polymer chain. Then, two subsequent relaxation runs, by adopting smaller grid sizes of l = 0.4 and 0.2 nm, were performed. III) Starting from the last configurations of the PF-MD simulations with the smallest grid size, a short minimization to remove residual atom superpositions followed by a short MD simulation (NVT at 293 K) of 5 ns were performed before running the NPT production runs (200 ns). In case of P2VP48, an additional grid size of l = 2.5 nm (close to Rg), was used first. Then, three subsequent runs were performed by using the same grid sizes employed in the P2VP10 case. All these computational details are summarized in Table I. The adopted scheme to get the density function from particle's positions is reported in Figures 2A-2B. In the panel C of the same figure, representations of the density isosurfaces, as function of the grid size l, are reported. As can be seen from the Figure 2C, asperities due to the bulkiness of the pyridine backbone substituents are found in the density isosurface shape (l = 0.2 nm).
In order to correctly apply the relaxation procedure to obtain wellrelaxed polymer melt, 33 the simulation time of each step of the procedure must be checked. To this aim we analyzed different quantities. In particular, since the equilibrium condition is reached when the diffu- Table IV. Ryckaert-Bellemans dihedral potential (for atoms in pyridine ring).  Dihedrals k dihedral (kJ mol −1 ) Multiplicity (n)  sion of the chains center of mass covers a length equal or greater than its Rg (0.8 nm for P2VP10 and 2.5 nm for P2VP48), we calculated the root-mean square displacement (RSMD) of the center of mass, as reported in Figure 3, from which it is clear that this condition is achieved for both chain lengths. The RMSD is calculated according to the definition reported in the Eq. 8: Where the angle brackets indicate an ensemble average over all the molecules in the simulations and all the time origins. By time origins, we denote that any timestep can be considered the time t 0 in the Eq. 8. As shown in the figure, this quantity is calculated at different time origins to make evident that the dynamics is time translational invariant.
Different grid sizes provide relaxation at different scales of the polymer melt, as described in the Ref. 32. Indeed, by analyzing the autocorrelation function (ACF) of the end-to-end vector, we can verify the relaxation time in which independent chain conformations are explored. In Figure 4, the ACF calculated for PF-MD simulations  are reported. As for the RMSD case, also the ACF is calculated at two different time origins. It is worth noting that the behavior we found is in agreement with the behavior reported for polymer melt 51 and polymer nanocomposite systems. [52][53][54] A fast decay to zero of the ACF is observed for all the grid sizes l. In particular, slower decays are observed as the grid size is decreased. This fact, also observed for the RMSD of the chains center of mass, is related to a decreased smoothness of the potential as the grid becomes finer, since increased detail on the polymer structure is added. The relaxation time τ of the end-to-end vector was evaluated by fitting its autocorrelation function with a stretched exponential: Table VII reports the results. In case of PF-MD simulations, the relaxation of a P2VP10 chain occurs in a time interval going approximately from 0.16 to 0.86 ns, as function of the grid size. A similar behavior is found for the P2VP48 ( Figure 4B). With the finest grid, the dependence of the relaxation time on chain length agrees well with the prediction of the Rouse model where τ∼N 3/2 : (48/10) 3/2 = 10.5.
The application of the procedure described above allows to obtain a suitable initial set of coordinates for a polymer melt in which the chains are directly packed at a density equal to the experimental value. The short-range correlations between atoms, that cannot be reproduced by the PF-MD method, can be easily reintroduced by performing a short MD simulation in which the LJ and Coulombic interactions are treated explicitly. In particular, in less than 1 ns the radial distribution functions (RDF) become indistinguishable with respect to RDF calculated after 200 ns ( Figures 5B-5C). This full recovery of short-range correlations was observed also for other polymers that were modelled with the same procedure. 25,33,34 Atomistic Model of Atactic P2VP.-In order to check if the model we proposed was able to reproduce experimental data, such as the mass density and X-Ray diffraction pattern, an MD simulation in the NPT ensemble was performed to this aim. In particular, starting from the last configuration obtained from the relaxation procedure (described above), a production run of 200 ns was performed for both P2VP10 and P2VP48, after a short NVT simulation of 5 ns.
The calculated mass densities are compared with the experimental data 55 in Table VIII. In particular, an error of about +3.5% for P2VP10 and less than -1% for P2VP48 is found for the proposed model. In particular, the largest error we found (3.5%) is in the range   of error (from 1.8 to 4.2%) usually obtained from the OPLS-AA. 22 In addition to the mass density, also the X-Ray scattering pattern was calculated and compared with the experimental data. 16 In Figure 6, a comparison of both X-Ray scattering patterns is reported for the P2VP48. In particular, the experimental X-ray data were measured on a P2VP sample with M w = 9 kg/mol 16 , which corresponds to almost two times the M w we simulated (P2VP48). As it can be seen from the comparison in Figure 6, the main peaks, at a q ∼ 8 nm −1 and at q ∼14 nm −1 representing the inter-chains distances and the Van der Waals (VDW) contact between atoms, 16 are centered at similar q values with respect to experimental ones. In particular, the difference in peak positions correspond to a very small distance difference in real space of about 0.2 Å, suggesting that the X-Ray scattering pattern was well-reproduced in our simulation. The equilibrium Rg values obtained from MD simulations calculated for P2VP10 and P2VP48 are: 0.62 and 1.51 nm for, respectively. Comparing our results with the ones calculated from the characteristic ratio of the mean-square of Rg (equal to 1.11), 56 we find a good agreement. In particular, from the characteristic ratio a value of Rg equal to 0.632 and 1.414 nm has been calculated for 10-48-P2VP chains, which are in the error range of 2-3%. In addition, considering the oligomer of PS with 10 monomers, which has similar Mw with respect to P2VP10, the calculated Rg (0.68) from MD simulations of the polymer melt, 38 is comparable with the value (0.62 nm) obtained from our simulations.

Conclusions
An all-atom model of atactic P2VP has been presented in this work. A procedure for an efficient relaxation of the P2VP melt, which is based on hybrid particle-field PF-MD method, has been applied with satisfactory results. The proposed model, based on the OPLS-AA force field, has been tested to reproduce a number of structural properties such as the mass density and X-Ray pattern. In particular, the mass density for the short oligomers (P2VP10) was reproduced within an error of +3.5%, which is in the error range expected from the OPLS-AA force field. For longer polymer chains (P2VP48) the density was reproduced within an error less than -1%, lower than the average error found for the OPLS-AA. Comparing the experimental and calculated X-Ray scattering data, we found a good agreement, which indicates, together with the density mass result, a good description of the structural correlations within the P2VP melt.