Structure and Dynamics of Solvated Polymers near a Silica Surface: On the Different Roles Played by Solvent

TL;DR

This study uses molecular dynamics simulations to explore how solvent and polymer interactions influence adhesion and structure of solvated polymers near silica surfaces, revealing the importance of explicit solvent modeling.

Contribution

It demonstrates the critical role of solvent in polymer-silica interactions and shows limitations of implicit solvent models in capturing these effects.

Findings

Polymer-surface adhesion varies with polymer and solvent interactions.

Explicit solvent models are necessary to accurately describe polymer adsorption.

Polymer chains exhibit significant dynamics and structural modifications near silica surfaces.

Abstract

Whereas it is experimentally known that the inclusion of nanoparticles in hydrogels can lead to a mechanical reinforcement, a detailed molecular understanding of the adhesion mechanism is still lacking. Here we use coarse-grained molecular dynamics simulations to investigate the nature of the interface between silica surfaces and solvated polymers. We show how differences in the nature of the polymer and the polymer--solvent interactions can lead to drastically different behavior of the polymer--surface adhesion. Comparing explicit and implicit solvent models, we conclude that this effect cannot be fully described in an implicit solvent. We highlight the crucial role of polymer solvation for the adsorption of the polymer chain on the silica surface, the significant dynamics of polymer chains on the surface, and details of the modifications in the structure solvated polymer close to the interface.

Tables

Table 1. Table 1: Interaction Parameter ϵ i ​ j subscript italic-ϵ 𝑖 𝑗 \epsilon_{ij} (kJ.mol -1 ) Matrix with Explicit and Implicit Solvent.

Explicit solvent	C	APAAm	APDMA	S	P4
C	2.625	2.042	2.71	2.042	2.3
APAAm	2.042	3.375	—	3.375	4.5
APDMA	2.71	—	4.0	4.5	4.5
S	2.042	3.375	4.5	3.0	4.5
P4	2.3	4.5	4.5	4.5	4.5
Implicit solvent	C	APAAm	APDMA	S
C	3.375	2.042	2.71	2.042
APAAm	2.042	2.042	—	2.042
APDMA	2.71	—	2.3	2.71
S	2.042	2.042	2.71	1.725

Table 2. Table 2: Bonded Interaction Parameters for Bonds and Angles, for PAAm and PDMA.

Bonds	Kbond (kJ.mol-1.nm2)	r0 (nm)
C–C	566.7	0.249
APAAm–C	666.7	0.237
APDMA–C	666.7	0.271
Angles	Kangle (kJ.mol-1)	${\theta }_0$ (degrees)
C–C–C	116.7	127.5
APAAm–C–C	233.3	85.5
APDMA–C–C	233.3	85.5

Table 3. Table 3: Dynamical quantities of PAAm and PDMA’s A Beads Interacting through an Implicit Solvent.

Dynamic	A beads of PAAm	A beads of PDMA
Number of active beads	38.8($\pm$12.3)	21.4($\pm$13.1)
Averaged number of events per active bead	1.9($\pm$0.4)	3.8($\pm$2.0)
Frequency of the event "bead leave" (Ghz)	1.0	0.8

Table 4. Table 4: Interaction Energy Between Beads and the Surface (in kcal.mol -1 ).

System	Zone I	Zone II	Zone III
A beads of PAAm	-2.8	-0.1	0
A beads of PDMA	-4.2	-0.3	0
Solvent beads	-4.6	-0.3	0

Table 5. Table 5: Dynamic of PAAm and PDMA’s A beads.

Dynamic	A beads of PAAm	A beads of PDMA
Number of active beads	4.1($\pm$1.7)	13($\pm$5.8)
Averaged number of events per active bead	1.5($\pm$0.6)	2.5($\pm$0.7)
Frequency of the event "bead leave"	0.1	0.4

Table 6. Table 6: Ratio of First Solvent Neighbors of PDMA and PAAm.

Region	Ratio PDMA/PAAm
Full box	1.4
Region III	2.1
Region II	0.7
Region I	0.8

Table 7. Table 7: Part of the Interaction Energy Between A and Solvent Beads among the Interaction Energy Between A Beads and every other Beads.

System	Part of the interaction energy in region I (%)
PAAm	31.1
PDMA	20.3

Equations

$$V_{bond}(r)=\frac{1}{2}{K}_{bond}(r-r_0{)}^2,$$

$$V_{angle}(\theta )=\frac{1}{2}{K}_{angle}[cos(\theta )-cos({\theta }_0){]}^2,$$

$$E$$

$$\displaystyle=4\epsilon\bigg{[}\Big{(}\frac{\sigma}{r}\Big{)}^{12}-\Big{(}\frac{\sigma}{r}\Big{)}^{6}\bigg{]}\frac{[r^{2}_{out}-r^{2}]^{2}[r^{2}_{out}+2r^{2}-3r^{2}_{in}]}{[r^{2}_{out}-r^{2}_{in}]^{3}}\qquad\qquad{}$$

$$=0$$

Full text

\alsoaffiliation

Fakultät für Mathematik und Naturwissenschaften, Stranski-Laboratorium für Physikalische und Theoretische Chemie, Sekr. C7, Technische Universität Berlin, Straße des 17. Juni 135, 10623 Berlin, Germany

\alsoaffiliationDepartment of Chemical and Biomolecular Engineering, Engineering Building I, Box 7905, North Carolina State University, 911 Partners Way, Raleigh, North Carolina 27695, United States

{tocentry}

1 Introduction

The design of nanocomposite materials is a fast-growing field, both for fundamental research and industrial applications, with many novel and exciting materials such as nanocomposite polymer hydrogels — reticulated polymer matrices, filled with water —, or nanoporous materials embedded in an amorphous matrix. Such materials can display an amazing array of mechanical 1, 2, 3, thermal, optical, electrical 4, responsiveness 5, frictional 6 and chemical properties.

Recently the use of nanoparticle solutions as adhesives has been demonstrated as a promising method for gluing soft matter, such as hydrogels or biological tissues 7, 8. Such a method has promising practical applications in the field of surgical adhesives 9, 10, 11, 12. At its core, this method relies on adsorption of the polymer onto the surface of the nanoparticles, with the nanoparticles acting as connectors between polymer gels. Local polymer chain rearrangements allow efficient dissipation of energy under stress and retard the fracture of the hydrogel. Although there have been early demonstrations of the feasibility of the method, a clear understanding of the reinforcement mechanisms involving both physical and chemical properties is still lacking and will be important for long-term progress in this area. To better understand the system, Rose and coworkers used hydrogels 13 as a simplified model of biological tissues 7. They compared the behavior of two kinds of hydrogels with regard to the surface: polydimethylacrylamide (PDMA) that glues to silica nanoparticles 14 and polyacrylamide (PAAm) which does not adhere to the surface of the silica nanoparticles 15.

Looking at the chemical structure of the two polymers (Figure 1(a)), this strong difference in behavior can appear surprising when thinking in terms of hydrogen bonds between silanol terminations (Si–O–H) that are present on the silica surface and in the polymers. PAAm is less sterically hindered than PDMA and one could expect stronger hydrogen bonds with the silica nanoparticles. However PAAm does not adsorb on the silica surface 15. One explanation could be that isolated silanol groups, which are the only binding groups for PAAm 16 (with available hydrogen-bond donnor), are not present on the surface of the particles used by Rose and coworkers. Another explanation could rely on the fact that the probability of forming hydrogen bonds between a silica surface and the polymer is lower than the probability of forming hydrogen bonds between silica and water (due to the high water content of such a hydrogel). This could explain why PAAm does not glue, however it does not explain why PDMA glues to the silica nanoparticles. Some authors have instead proposed looking at the hydrogel-silica hydrophobic interactions 17, 18. According to them PDMA adsorbs on the silica surface by hydrogen bonds between oxygen of the dimethylacrylamide group and hydrogen of the silanol group of the silica surface. The ability of PDMA to adsorb on the silica depends on the strength of the hydrogen bond, but also on the hydrophobic interactions of the two alkyl groups from PDMA. The latter reinforces the hydrogen bond through the $p$–$\pi$ conjugation effect between a nitrogen atom and the silanol group and removes water molecules from the silica surface 18, 17, 19. This effect does not exist for PAAm because there are no hydrophobic groups. As these hypotheses are not supported by direct experiments, a better understanding of the gluing of a hydrogel to a silica surface can only be acquired by simulation.

Therefore we investigated this phenomenon with coarse grained (CG) simulations in order to compute rather large systems at long timescales (50-100 ns). We decided to use a semi-coarse-grained approach in order to keep some of the chemical details. This is a compromise between a detailed atomistic approach and the efficiency of a CG model. We designed a simple and transferable model between the two polymers in order to focus on the origin of the different behavior between PAAm and PDMA (by limiting the number of parameters).

Previous CG simulations have concentrated on mechanical properties but were not focused on their chemical origin 20, 21. The coupling between the local dynamic of the polymer chains on the surface and the resulting mechanical properties is therefore only poorly understood. The solvation of polymer chains, which prevents them from adsorbing on the surface, and the competition between adsorption of solvent or polymer on the surface have, to the best of our knowledge, not been studied to date. We focus on a rather simple system where polymer chains are free of crosslinks and we neglect the curvature of the nanoparticle and consider a flat silica surface.

The paper is organized as follows: first, we present the CG model employed herein. Then, we use a simple model of a flat silica surface and a polymer solution in implicit water to investigate the local structure and the dynamic of PDMA and PAAm close to the surface. We then introduce an explicit water model to look at the role of the polymer solvation with regard to their adsorption on the silica surface.

2 Methods

2.1 Coarse-grained model

We present here the CG model used for the polymer chain, for the silica surface, and for the water. Coarse-graining a system consists in finding a representation that maps a group of atoms into one CG bead. We designed the CG model within the framework of the Martini force field. The first version of this force field was developed by the Marrink group in 2004 22 and improved in the 2007 version 23, which is the version we use in this work. The overall aim of the Martini CG approach is to provide a simple model that is computationally fast and easy to use, yet flexible enough to be applicable to a large range of biomolecular systems. The Martini force field has been successfully used for a variety of biomolecules, including lipids 22, sugars 24, proteins 25, and polymers such as: polystyrene 26, polyethylene oxide and polyethylene glycol 27, polycaprolactam 28, and poly(methyl methacrylate) 29. The Martini force field is based on a four-to-one mapping, meaning that one CG bead represents on average four heavy atoms and connected hydrogens. This choice is a compromise between computational efficiency and chemical representativeness. The Martini model has four main types of CG bead: polar (P), non-polar (N), apolar (C) and charged (Q). Each particle type has several subtypes. This allows for a more detailed description of the underlying atomistic structure. There is a total of eighteen subtypes. For example, polar beads are distinguished by a number indicating the degree of polarity (from 1 for low polarity to 5 for high polarity). Non-polar beads are distinguished by a letter indicating the hydrogen-bonding capability ($d$ for donor and $a$ for an acceptor of a hydrogen bond, $da$ for both and 0 for non). The use of eighteen subtypes of beads allows one to set up a rather complicated molecule, with a limited number of building blocks. We used them to design the CG model of the polymers PAAm and PDMA, of the silica surface and of the water.

CG model of PAAm and PDMA. There are two steps to design a CG model of an all-atoms structure using the Martini force field. The first step is to map chemical groups onto CG beads. There are several ways to coarse-grain a system. One of the most widely used is the Kremer-Grest bead-spring model 30 where one CG bead corresponds to one or more monomer of the polymer chain. In this description, the CG representation of a polymer is a chain of beads connected with a spring potential. In this work we use a more detailed CG model where one monomer is mapped onto two beads: one bead for the backbone chain (grey circle in figure 1(a)) and one for the side chemical function (red circle in figure 1(a)). The "C" bead (C for Chain) accounts for two carbon atoms of the backbone chain of PAAm and PDMA. The "A" bead (A for Amide), represents the side groups of PAAm and PDMA. For PAAm, the side group is an amide function. For PDMA, it is a dimethylacrylamide function. Therefore a C bead is similar for PAAm and PDMA, whereas an A bead differs between PAAm and PDMA, with an APDMA bead being bigger than an APAAm bead. The second step is to associate the corresponding type of bead, using the 18 subtypes we just mentioned. According to Martini’s type of particles, a C bead is a bead of type C1: it is an apolar bead with a low degree of polarity. APAAm is a bead of type Nda: a non-polar bead which can donate hydrogen bonds (through NH2) or accept hydrogen bonds (through oxygen). APDMA is a bead of type Na which can only accept hydrogen bonds. Once the type of bead is determined, its size must also be evaluated as the size of the bead plays an important role in the interaction parameters of the system that will be presented later. Martini is based on a four-to-one mapping: beads that map four or more atoms are considered as "normal" beads. Beads that map less than four atoms are "small" beads. C bead is a small bead because it maps only two carbon atoms and connected hydrogens. An APAAm bead maps three heavy atoms: it is also a small bead. APDMA is a normal bead because it maps four heavy atoms.

CG model of the silica surface. The model we use for the silica surface has a simple geometry and CG model. First, the geometry of the surface is flat to get rid of the effect of the nanoparticle’s curvature at this stage. This choice is motivated by the small curvature of the silica surface in the considered system. The silica nanoparticle we perform matches with silica nanoparticles that are experimentally used 7, Ludox TM-50. It is a water solution with concentration of 52 wt% of silica nanoparticles and radius of 15 nm. Our silica surface has dimensions of $l_{x}$ = $l_{y}$ = 63 Å. Considering a silica nanoparticle with a radius of 15 nm would lead to a height difference of 3 Å between the edge of our surface and its center. Such curvature, considering the properties we are interested in (local dynamic quantities, interaction energy, solvation of the polymer chain) can be neglected and a flat silica surface is considered in this work. Silica surfaces exhibit several types of silanol terminations as one can see from figure 1(b)) 31, 32. There are isolated silanols that do not interact with other silanols of the silica surface, H-bonded silanols that interact via hydrogen bonds, vicinal silanols which are two silanols forming a siloxane Si–O–Si bridge and geminal silanols wherein one silicon atom bears two hydroxyls. Our CG model is simple in the sense that we choose to map one silanol group into one CG bead called "S" (for silanol). This is shown in the bottom part of the figure 1(a). We designed our flat surface in order to reproduce the averaged silanol density of an amorphous silica surface, corresponding to 5 OH/nm2 33, 34, 35. Within the Martini framework, silanol beads are non-polar beads (N) which are hydrogen-bond donors, Nd. One last thing to note is that the beads of the polymer and of the silica surface are comparable in terms of size. We use a simple surface made of one layer. We checked that using a slab containing several layers (up to eleven layers) does not change the considered properties. Figure LABEL:fig:1-11layers is a comparison of the potential of mean force of PAAm when one or eleven layers are used for the silica surface. They are similar, showing that the use of one silica layer is reasonable.

CG model of water. There are two ways to simulate water within the Martini framework: either explicitly or implicitly. The Martini model of explicit water comes down to mapping four water molecules onto one bead (P4 beads, according to the type of beads) 23. The use of the explicit Martini water in a system containing polymer chains has already been done with polystyrene26, polyethylane oxide and polyethylene glycol27. The interaction parameter $\epsilon$ between two solvent beads is proposed to be 5.0 kJ.mol*-1*. However, is has been reported that water, modeled as P4 particles, has a freezing temperature that is too high compared to real water. This is linked to the use of a Lennard-Jones 12-6 potential for non-bonded interactions which over-emphasizes the formation of structures 36. This is particularly observed in systems where a nucleation site is already present (like a solid surface in our case). We note, for a system containing a surface and 2064 water beads, a freezing temperature around 360 K. To avoid this problem, we decided to decrease the $\epsilon$ interaction parameter of water beads to 4.5 kJ.mol*-1*. With this new interaction parameter, water close to the surface does not freeze at 300 K (see Supplementary Fig. LABEL:fig:water+surface online). Because we modify the interaction parameter given by Martini for water and we do not strictly use the water from Martini, we decide to use the term "solvent" instead of "water" in our case. However, using an explicit solvent is computationally expensive, especially in our case where experimental hydrogels are composed of 90 % of water in weight. Almost all of the considerable computational time is spent on calculating interactions involving solvent beads.

One solution is to use implicit solvent. There is a wide range of implicit solvents37 employed in polymer science such as Brownian Dynamics (BD)38, Dissipative Particles Dynamics (DPD)39 or Lattice-Boltzmann40. In this work, molecular dynamics simulations are performed. We do not account for the random force due to the solvent as in BD or in DPD. The implicit solvent we use consists of tuning the interaction parameters to take into account the effect of explicit solvent. Indeed, the interactions of a solvophilic bead with other beads would be screened by the presence of a solvation shell in explicit solvent. Therefore, the interaction parameter of solvophilic beads is reduced to take into account the screening due to solvent. On the contrary, there is a depletion of solvent molecules around a solvophobic bead: the interaction parameter of solvophobic beads is increased in an implicit solvent model. We use the implicit solvent CG model developed by Marrink and coworkers called Dry Martini 41. This implicit solvent developed within the Martini framework has extensively been used to study systems containing polymers with molecular dynamics 42, 43, 44, 45, 46. Such implicit solvent is easy to use and inexpensive. However we will show that the use of this implicit solvent fails to reproduce details of the microscopic interactions, making the use of explicit solvent necessary in some cases.

2.2 Nonbonded interactions

Once the CG model is designed, the interactions occurring between the beads must be considered. There are two kinds of interactions: bonded interactions between chemically connected sites and nonbonded interactions. In the Martini force field, non-bonded interactions are described by a Lennard-Jones (LJ) 12-6 potential and there are typically no Coulombic interactions explicitly considered for such noncharged polymers 26, 27, 28. The strength of the interaction, determined by the value of the LJ welldepth $\epsilon_{ij}$, depends on the interacting beads $i$ and $j$. The $\epsilon_{ij}$ for normal beads (beads that map four or more atoms) are found in the interaction matrix given in the paper from Marrink and coworkers 23. For small beads (beads that map two or three atoms), $\epsilon_{ij}$ is scaled to 75% of the standard value. The effective size of the particles is governed by the LJ parameter $\sigma_{ij}$. $\sigma_{ij}$, which is 0.47 nm for normal particles and 0.43 nm for small particles. The LJ potential is shifted to zero between $r_{shift}=0.9~{}nm$ and $r_{shift}=1.2$ nm. Both self and cross interaction parameters $\epsilon_{ij}$ are given in Table 1. The upper part of the table summarizes the interaction parameter $\epsilon_{ij}$ for a system containing the explicit solvent we described before. For instance, a C bead interacts with APAAm bead with an $\epsilon_{ij}=2.042$ kJ.mol*-1*. But APAAm and APDMA do not interact because our system contains eitheir pur APAAm or pur of APDMA.

2.3 Bonded interactions

Bonded interactions between two chemically connected beads are described by a harmonic potential $V_{bond}(r)$:

[TABLE]

where $r$ is the distance between two beads, $K_{bond}$ is the force constant of the harmonic bonding potential, and $r_{0}$ is the equilibrium distance. Angles between three neighboring beads are governed by the potential $V_{angle}(\theta)$:

[TABLE]

where $\theta$ is the angle between the three beads, $K_{angle}$ is the force constant and $\theta_{0}$ is the equilibrium angle. We decided not to use dihedral potential for the sake of simplicity. The Martini paper 23 provides general values for $K_{bond}$, $K_{angle}$, $r_{0}$ and $\theta_{0}$. We find those values not adapted to our system, which is a rather fine grained model, $i.e$ containing small beads. We have thus reparametrized our CG model, based on input from all-atoms simulations. Parameters for the bond and angle potentials were obtained by comparing distributions from all-atom simulations with distributions from CG simulation.

We ran all-atoms simulations, using the LAMMPS code 47, 48 and the CHARMM22 force field 49 (the timestep was set to 1 fs). The CHARMM force field is widely used for all atoms simulations and is adapted to the system we study. Along the all-atoms simulation we compute lengths and angles between the center of mass of the groups of atoms that are mapped onto one CG bead. The resulting histograms of bond length $r$ and of the angle $\theta$ are the red curves in figure 2. The CG histograms are the green curves in 2. Both histograms are normalized. They indicate the probability of finding a bond at a given distance or an angle at a given value.

All-atoms and coarse-grained simulations are done with a system containing five polymer chains of twenty monomers. The averaged value of the all-atoms histograms is used to optimize the equilibrium bond length $r_{0}$ and angle $\theta_{0}$ of the CG potentials $V_{bond}(r)$ (Eq. 1) and $V_{angle}(\theta)$ (Eq. 2). We tune $r_{0}$ and $\theta_{0}$ in such a way so that the averaged value of the CG histogram matches the averaged value of the all-atoms histogram. The coupling constants $K_{bond}$ and $K_{angle}$ are optimized by comparing the width of the CG histogram with the corresponding width of the all-atoms histogram. End-to-end distances (the distance between the two extremities of a polymer chain) and radius of gyration between all-atoms and CG simulation are finally compared to validate our CG model. We find, for the radius of gyration, an average value of 7.1 $\pm$ 0.3 Å for the all-atoms system and of 7.3 $\pm$ 1.0 Å for the coarse-grained system. The values of the radius of gyration are in very good agreement between all-atoms and coarse-grained systems. As for the end-to-end distance, the average value is 15.07 $\pm$ 2.53 Å in the all-atoms system and 13.21 $\pm$ 3.33 Å in the coarse-grained system. End-to-end distances also show good agreement between all-atoms and coarse-grained simulations. The final Kbond, Kangle, r0 and $\theta_{0}$ parameters are summarized in Table 2. For the sake of simplicity and to have a CG model that differs as little as possible between PAAm and PDMA, Kbond, Kangle and $\theta_{0}$ are the same for PAAm and PDMA. Only APDMA–C differs. This arises from the comparison between all-atoms and CG distributions.

2.4 Numerical details

CG molecular dynamic simulations were done using the LAMMPS simulation package 48, 47 and the Martini force field developed by Marrink and coworkers 22, 23, 36, 41. The LJ interactions were smoothly shifted to zero between 0.9 and 1.2 nm by using the following formula:

[TABLE]

where $r_{in}$ = 0.9 nm and $r_{out}$ = 1.2 nm and $r$ is the distance between two particles. The timestep was set to 10 fs and coordinates were saved every 5 ps for analysis. Periodic boundary conditions in $x,y,z$ directions were used. The systems were first equilibrated in the NPT ensemble and then in the NVT ensemble during 5 ns at a pressure of 1 atm and a temperature of 300 K. For both systems, simulations are done in the NVT ensemble. Production simulations were performed in the NVT ensemble during 50 ns to 100 ns. Nose-Hoover thermostat and barostat were employed and the pressure was set to one atmosphere. The relaxation time was set to 100 fs for the temperature and to 1000 fs for the pressure. Several independent simulations (5 or 10) are performed and are used to compute averaged quantities. Input files for system with an implicit or with an explicit solvent are provided in the Supplementary material. No coupling method was used for the interaction parameters.

3 Results and discussion

3.1 Solvated polymer interacting through an implicit solvent on a surface

To understand the influence of the nanoparticles on the polymer network, we first look at the local structure and dynamics at the interface between the surface of the nanoparticles and the polymers. As a starting point, we studied a polymer solution near a surface interacting through an implicit solvent. The system is composed of a flat surface (200 S beads) and 24 chains of 90 monomers (2160 A beads and 2160 C beads). A snapshot of the simulation box is displayed on Fig. LABEL:fig:snapshotPAM_nolength of the Supplementary material. The simulation box dimensions are : $l_{x}$, $l_{y}$ = 63 Å, $l_{z}$ = 51.4 Å for PAAm and $l_{z}$ = 55.8 Å for PDMA after equilibration at a pressure of 1 atm. We checked that the polymer slab is large enough so that polymer beads located in the middle of the slab do not feel any interaction from the surface. The system was first equilibrated during 5 ns in the NPT ensemble at a pressure of one atmosphere (pressure was applied in the direction normal to the surface) and a temperature of 300 K, then during 5 ns in the NVT ensemble. We follow the dynamic for 50 ns and save the coordinates of the beads along the simulation. The results we present here are averaged over 5 independent simulations of 50 ns each. As A beads interact more strongly with the surface than C beads and are different between PAAm and PDMA, we concentrate only on A beads for the rest of this work. Figure 3(a) is a side view of the polymer conformation on the silica surface where only three over twenty-four PDMA chains are represented. It illustrates how polymer chains are well entangled and form trains and loops close to the silica surface. Figure 3(b) shows that A beads in region I are well dispersed on the silica surface.

We computed the normalized histogram $P(d)$ of the distance $d$ between the surface and the polymer A beads. The free energy profile of the A beads is then $F(d)=-k_{B}Tln[P(d)]$, plotted in figure 4.The surface introduces a symmetry breaking that forces the polymer beads to organize in layers. At small $d$, there is a well corresponding to a high stability of the polymer beads: they form a first layer in the direct vicinity of the surface, then a second layer. Beyond that second layer, there is no more influence of the surface on the polymer beads. Both PAAm and PDMA beads are similarly structured close to the surface.

From this structure we define different regions of stability that are represented in figure 4. Region I is defined as the first layer near the surface: it starts close to the surface and ends at the first maximum as one can see from figure 4. A beads are well dispersed on the silica surface, as can be seen in Fig. 3(b). Region II starts just after region I and ends at the second peak of the probability distribution. This second layer is less structured than the first layer. Region III corresponds to a bulk-like region where there is no interaction between the polymer and the surface and is defined by a slab of 10 Å in the middle of the box.

From an energetic point of view, in region I, the averaged interaction energy between the surface and A beads is 0.9 kcal.mol*-1* per bead for PAAm and 1.6 kcal.mol*-1* per bead for PDMA. It is consistent with the interaction parameters that are given in Table 1: APDMA beads interact more strongly with S beads than APAAm. This is interesting to compare those values with the free energy barrier $\Delta F_{B}$, obtained by subtracting the value of the PMF at the first peak from the PMF’s value at the first well. It characterizes the local exchange of the A beads, or the ability to detach/attach an A bead from the surface. $\Delta F_{B}$ turns out to be slightly higher for PDMA (1.7 kcal.mol*-1*) than for PAAm (1.4 kcal.mol*-1*). It is thus a little harder for an A bead of PDMA to leave the first adsorption layer to go to the second layer or to the bulk. The results obtained for the free-energy barrier and for the interaction energy between A beads and the surface in region I are then consistent. We can conclude that there is a slight difference between PAAm and PDMA when the behavior in the vicinity of the surface is considered. A beads of PDMA interact more strongly with the surface, making it more difficult for them to leave the first adsorption layer.

Our aim is now to go beyond energetic and thermodynamic information and to study whether they have consequences on the ability of the polymer beads to move in the vicinity of the surface or not. Indeed, the capacity of polymer beads to attach and detach and to reorganize near the surface is linked to the mechanical properties of the resulting nanocomposite system. If the system has the ability to rearrange in an efficient way, it will dissipate energy under stress and will lead to a stronger system 50, 51. In order to quantify the ability of the polymer beads to move close to the surface, we consider as "active beads" A beads that cross regions over the course of the simulation (regions being defined in Fig. 4). Polymer beads undergo two kinds of events: whether they leave region I (close to the surface) or whether they reach region I. Beads that leave region I go from region I to beyond region II. Beads that reach region I arrive in region I from outside region II. The averaged number of events is in the second row of Table 3. The number of occurrences of the second event divided by the simulation time is the frequency of this event and is in the third row of the table. We compare the numbers between A beads of PAAm and of PDMA.

First of all, it is worth noting that PAAm and PDMA have rather few active beads over the course of the simulations. In order to be sure that active beads are not only end-of-the-chain beads but are rather well dispersed along the polymeric chain, we checked where the active beads are located (see Supplementary Fig. LABEL:fig:active_beads online). Figure LABEL:fig:active_beads shows the position of active beads along doubled up polymer chains, averaged over the 24 chains and over the course of the simulations: position 1 is the end of the chains and position 45 is right in the middle of the chain. The figure shows an excess of active beads at the end of the chains, but also active beads all along the chains. Therefore we are confident that the dynamical behavior of polymer chains is due to the motion of the whole polymer chain rather than small motions of end beads. Moreover, one can wonder whether "nonactive" beads, which are beads that stay close to the surface during the entire course of the simulation, play a role with regards to the adsorption of polymer chains on the silica surface. To shed light on this point, we compute the number of nonactive beads, for both PAAm and PDMA, and averaged it over the five independent simulations. PAAm has 446 nonactive beads and PDMA has 427 beads. This difference is rather small and PAAm, which experimentally does not adsorbs on the silica surface, has even more nonactive beads than PDMA. Thus, nonactive beads do not explain the different behavior between PAAm and PDMA. One can see from Table 3 that PAAm has more active beads than PDMA, but that PDMA’s active beads undergo more events than PAAm’s active beads. The consequence is that the dynamical properties of A beads of PAAm and A beads of PDMA is the same: the resulting frequency of the event "bead leaves" is the same for PAAm and PDMA for instance. At this stage, the study of the dynamical properties of A beads does not allow us to draw a conclusion about a difference in the behavior between PAAm and PDMA that would lead to different mechanical properties.

To conclude on the implicit solvent model, even if it gives a slight difference of $\Delta F_{B}$: it is more difficult for APDMA beads to move away from the surface (which is confirmed by the interaction energy between A beads and surface), it does not lead to significant difference of the dynamical properties of PAAm and PDMA. Thus, we need a more complex system, including explicit solvent, to understand the differing behavior between PAAm and PDMA. This indicates that experimentally, water not only plays a role by screening the interactions between polymer chains and surface. Water molecules also play a steric role by competing with the adsorption of the polymers on the silica surface; they will prevent polymer chains, or not, to adsorb on the silica surface. Therefore, water molecules have to be taken into account to explain the different behavior of PAAm and PDMA. The resulting strengthening of polymer network by the addition of nanoparticles is not only explained by energetic considerations. In the resulting strengthening of polymer network by the addition of nanoparticles, explicit solvent plays an important role.

3.2 Solvated polymer interacting through an explicit solvent on a surface

We use systems containing 50% of solvent by weight. This is a compromise between the wish to be as close to experimental systems as possible (which contain around 90 % of water 7) and the computational cost that the simulation of an explicit solvent implies. Using 50% of solvent in weight insures having the correct number of solvent beads in the simulation box. We use the explicit solvent described above. We still consider 24 chains containing 90 monomers and 200 beads for the surface. As we work with a percentage in weight, the PAAm system contains 960 solvent beads and PDMA system contains 1340 solvent beads (Fig. 5). A snapshot of the simulation box is displayed on Fig. LABEL:fig:snapshotPAM+solvent_nolength of the Supplementary material. The simulation box dimensions are : $l_{x}$, $l_{y}$ = 63 Å, $l_{z}$ = 80.9 Å for PAAm and $l_{z}$ = 103.4 Å for PDMA. The system was first equilibrated during 5 ns in the NPT ensemble at a pressure of one atmosphere and a temperature of 300 K, then during 5 ns in the NVT ensemble. Figure 6 shows APDMA beads in region I. There are fewer A beads in the vicinity of the surface than for the polymer solution interacting through an implicit solvent (figure 3(b)). The rest of the surface is covered by solvent beads. The same behavior is observed for PAAm. The results we present are averaged over 10 independent simulations of 50 ns. The initial configuration is the same for PAAm and PDMA. First, the behavior of the polymer beads near the surface is studied, as have already been done for the polymer solution interacting with an implicit solvent, to see if it is modified by the presence of explicit solvent beads. Then, the way polymer beads are solvated by the explicit solvent beads is probed.

Polymer behavior with regard to the surface. The free energy barrier $\Delta F_{B}$ of A beads is higher for PDMA (2.1($\pm$0.3) kcal.mol-1) than for PAAm (1.6($\pm$0.6) kcal.mol*-1*). It is more difficult for APDMA beads to leave the region I. Moreover, the first two rows of Table 4 show that in region I, APDMA interact more strongly with the surface than APAAm. This is also what we observed with an implicit solvent. Using an explicit or an implicit solvent does not change the ranking of interactions between the polymers and the surface. The strength of interaction energy with an explicit solvent is stronger than with an implicit solvent because the implicit solvent, by screening the interactions between A and S beads which are both solvophilic beads, decreases their interaction parameters. However, what are the consequences of the use of an explicit solvent on the residence time of PAAm and PDMA near the surface?

One can see from Table 5 that PDMA has more active beads than PAAm and that these active beads are more active than the active beads of PAAm. Not only they undergo more events, but their frequency is also higher. Therefore A beads of PDMA are slightly faster than A beads of PAAm. However, we see from table 4 that A beads of PDMA interact more strongly with the surface than PAAm. The fact that PDMA has a stronger interaction with the surface but moves faster in the vicinity of the surface is surprising. This is certainly due to the fact that the interaction energy between APAAm and the surface is low below the interaction energy of the solvent beads and the surface: PAAm is then replaced by solvent beads near the surface. This is not the case for APDMA, which has an interaction energy with the surface in the same range as the interaction energy between the solvent beads and the surface. Explicit solvent apparently plays a role in this system and modifies the dynamical properties of the polymer beads, even if the interaction energy ranges are not changed.

Solvation of the polymer. To understand how the presence of solvent beads modifies the dynamics of polymer beads, it is interesting to take a close look at the way polymers are solvated by the explicit solvent. We thus analyse the solvation of PAAm and PDMA to get some insight on a possible competition between solvent beads or polymer beads with regard to the adsorption on the surface. We quantify the solvation of PAAm and PDMA by computing the number of first solvent neighbors of PDMA and PAAm A beads in regions I, II and III. To do so, we count the number of solvent beads that stand at a certain distance $d_{shell}$ from A beads, where $d_{shell}$ corresponds to the radius of the first solvation shell of A beads. This procedure is done for A beads of PAAm and of PDMA and in regions I, II and III. NPAAm is the number of solvent first neighbors of A beads of PAAm; NPDMA is the number of solvent first neighbors of A beads of PDMA and we present the ratio NPAAm/NPDMA in Table 6.

One can first note that the ratio of first neighbors of PDMA and of PAAm in the full simulation box is about 1.4, which gives the number of solvent beads for PDMA divided by the number of solvent beads for PAAm. This "full box" ratio, that we use as a reference, is compared with the ratio in regions III, II and I. In region III, the ratio goes up to 2.1, meaning that APDMA has more first solvent neighbors than APAAm. The value of the ratio is above 1.4: there is an "excess" of solvent beads for PDMA, or a default of solvent beads for PAAm in the bulk-like region (region III), compared to the full box ratio. The tendency is well understood when we move toward the surface, in regions I and II, where the ratio lowers to 0.7 (region II) or 0.8 (region I). There is clearly, as one can see from figure 5, a large excess of solvent beads close to APAAm beads close to the surface (left side of the figure 5), compared with the box containing PDMA (right side of the figure 5). It is worth noting that we started the simulation of the systems containing PAAm and PDMA from the same initial configuration and led to different configurations of the solvent around the polymer. Therefore, in the box containing PAAm, solvent beads leave the bulk-like region to reach the interface with the surface resulting in a default of solvent beads in region III and an excess of solvent beads in regions I and II. In the PAAm system, solvent beads go in between the polymer chains and the surface, preventing PAAm chains from adsorbing to the surface. For PAAm, there is a competition between solvent beads and PAAm beads with regard to the adsorption on the surface. However there is also an attraction between PAAm beads and solvent beads. The attraction between PAAm and solvent is confirmed by the study of the part of the interaction energy between A beads and solvent among the interaction energy between A beads and all the other beads, which is summarized in table 7. Indeed, as PAAm and PDMA do not interact with the same strength with the surrounding beads, we can not directly compare A/solvent interactions energy between PAAm and PDMA. The first row of table 7 shows that, in region I, the part due to A/solvent interactions among the A/everything interactions is higher for PAAm than for PDMA. Therefore APAAm interacts more forcefully with solvent beads than PDMA close to the wall.

To conclude on this part, we show that explicit solvent plays an important role. Even if, at first glance, the presence of explicit solvent beads does not seem to perturb the interaction between A beads and surface beads, it makes the adsorption of APAAm more difficult by solvating the polymer’s beads and by moving toward the interface and preventing the polymer from adsorbing on the silica surface. It is indeed easier for solvent beads to adsorb on the surface because they are smaller and move faster than APAAm beads.

4 Conclusions and perspectives

In this work, we performed coarse-grained molecular dynamic simulations of a solvated polymer near a planar model silica surface, comparing implicit and explicit solvent models. We conclude that the implicit solvent model gives reasonable results regarding the interaction energy between polymer and surface, as well as the ability for the polymer beads to move away from the surface. However, it fails to describe the dynamical properties of the polymer beads near the silica surface. The explicit solvent model, which has much higher computational cost, is necessary to describe the competition between solvent and polymer beads near the interface. In the PAAm-based system, interactions between solvent and surface are stronger than those between PAAm and surface. As PAAm is well solvated by water and prefers to be surrounded by solvent beads than close to silica beads, this prevents adhesion of the polymer on the silica surface. This is in stark contrast with the PDMA chain, where polymer–surface interactions dominate and lead to adhesion on the surface. Moreover, we highlighted the crucial role of polymer solvation for the adsorption of the polymer on the silica surface, the significant dynamical properties of fragments of polymer on the surface, and detail the modifications in the structure of the polymer close to the interface.

{suppinfo}

The following files are available free of charge on the online Supplementary material.

• •

Figure LABEL:fig:1-11layers is a comparison of the potential of mean force of PAAm with one and eleven layers of silica surface.

• •

Figure LABEL:fig:water+surface is a snapshot of a simulation box containing solvent beads and a silica surface.

• •

Figure LABEL:fig:snapshotPAM_nolength is a snapshot of a simulation box containing PAAm and silica surface interacting through an implicit solvent.

• •

Figure LABEL:fig:active_beads is the distribution of active beads along polymer chains.

• •

Figure LABEL:fig:snapshotPAM+solvent_nolength is a snapshot of a simulation box containing PAAm, the silica surface and an explicit solvent.

{acknowledgement}

This work was supported by a PSL-Chimie grant and founded by the French Ministry of Higher Education and Research. The authors thank Alba Marcellan for insightful discussions. The computational resources of the Technical University of Berlin are also acknowledged. This work benefitted from access to HPC platforms provided by a GENCI grant (A0030910299).