Structure and Mobility of Nanoconfined Polyamide-6, 6 Oligomers

Mar 31, 2009 - Structure and Mobility of Nanoconfined Polyamide-6,6 Oligomers: Application of a Molecular Dynamics Technique with Constant Temperature...
0 downloads 0 Views 553KB Size
5568

J. Phys. Chem. B 2009, 113, 5568–5581

Structure and Mobility of Nanoconfined Polyamide-6,6 Oligomers: Application of a Molecular Dynamics Technique with Constant Temperature, Surface Area, and Parallel Pressure Hossein Eslami*,†,‡ and Florian Mu¨ller-Plathe† Eduard-Zintl Institut fu¨r Anorganische and Physikalische Chemie, Technische UniVersita¨t Darmstadt, Petersenstraβe 20, D-64287, Germany, and Department of Chemistry, College of Sciences, Persian Gulf UniVersity, Boushehr 75168, Iran ReceiVed: December 20, 2008; ReVised Manuscript ReceiVed: February 19, 2009

A molecular dynamics simulation method with coupling to an external bath is used to simulate polyamide6,6 trimers confined between graphite surfaces. In this simulation method, the temperature and the parallel component of pressure are kept fixed, and the distance between the confining graphite surfaces is changed to achieve equilibrium. The simulation results on the oscillatory behavior of solvation force, the number density of confined oligomers, and stepwise variation of the oligomer numbers as a function of distance between the confining graphite surfaces are reported and discussed. The hydrogen bonding in confined oligomers has also been studied in detail, and it is shown that hydrogen-bond formation depends on the layering effect and on the geometrical restrictions and reveals an oscillatory behavior like the solvation force oscillations. The autocorrelation functions for NH, CO, and end-to-end vectors are also studied, and it is concluded that the confined fluid has a considerably lower relaxation time than that of the bulk fluid, and the relaxation times for confined fluid show an oscillatory behavior with maxima corresponding to well-formed structures parallel to the surfaces. The study of local dynamics via calculating the autocorrelation functions for bins parallel to the surfaces reveals that the fluid layers close to the surfaces have higher relaxation times than the fluid in the central region. Our calculated center-of-mass diffusion coefficients also show oscillatory behavior with outof-phase oscillations with respect to the solvation force oscillations. Introduction An understanding of the atomic processes occurring at the interfaces plays a central role in many technologically important phenomena such as adhesion, lubrication, coating, friction, wear, wetting, spreading, chromatography, and membrane separation.1,2 Polyamide (PA) thin films have been studied as one of the suitable membrane materials because of their high thermal stability, excellent mechanical strength, and high resistance to organic solvents.3 PA thin film composite membranes show high selectivity in dehydration of alcohols, at a wide range of water concentrations, and are of practical importance in such phenomena as reverse osmosis4 and nanofiltration separation processes.5 The applicability of PA-imide polymers to thin film optical application has also been reported.6 These films are known as highly robust and stable films, which are excellent candidate materials for integrated optics applications.6 Advances in experimental techniques, such as the use of the surface forces apparatus, atomic force microscopy, and friction force microscopy to study the properties of confined fluids, have considerably increased the amount of experimental information on these systems.1,2 Although experimental studies have revealed many aspects and properties of confined liquid films, there is still a lot to be known about these systems. Moreover, the experimental works on the confined fluids are very difficult to perform, due to the existence of only a small amount of material in the confined region. * Corresponding author. E-mail: [email protected]. † Technische Universita¨t Darmstadt. ‡ Persian Gulf University.

Confined fluids, on the other hand, have been studied using both analytical7-11 and computer simulation methods.12-29 Monte Carlo (MC) and molecular dynamics (MD) simulation studies allow direct investigation of such complex systems at molecular level. Early MD simulation methods in this field on model monatomic liquids reveal interfacial local liquid density oscillations and variations in transport and dynamic properties of the liquid as a function of distance from the interface.12-19 In more recent studies, MC and MD simulations are applied to study more complex fluids.20-29 Although experimental as well as simulation methods have revealed a wealth of information regarding the energetic, structural, dynamical, and rheological properties of interfacial and confined molecular fluids under static and shear conditions, still less is known about the molecular level properties of confined oligomers and polymers. Experimentally, examples of polymer orientation, confinement, and frustration at small scales are reported as shifts in measured glass transition temperatures,30,31 changes in polymer thermophysical properties,32 or changes in mechanical viscoelastic properties.33,34 On the simulation front, some of the pioneering studies focusing on the behavior of abstract Lennard-Jones oligomers in nanoscopic confinements consist of studies of the relaxation times of different modes, of the transport properties like the self-diffusion coefficients and the mean square displacements,35-38 and of the rheological properties and the viscosity.38 Yethiraj and Hall39,40 performed a MC simulation to study the density and conformational behavior of polymers, modeled as tangenthard-spheres, confined between hard surfaces. In a similar study, Aoyagi et al.41 studied by MD simulation a coarse-grained

10.1021/jp8112655 CCC: $40.75  2009 American Chemical Society Published on Web 03/31/2009

Structure of Nanoconfined Polyamide-6,6 Oligomers bead-spring polymer model confined between two solid surfaces, and the surface effect was studied by adjusting distance between the surfaces and the surface-polymer interactions. They also studied the position-dependent dynamics by analyzing the trajectory of the center of mass of polymer chains. In most of the simulations performed on confined polymer systems, the investigators used a generic bead-spring model to study, mostly, the equilibrium structure of the confined polymer melts, or in some cases, as addressed above, they studied the diffusion coefficient, or the position-dependent diffusion coefficient. However, recently there have been some reports of detailed MD simulations of polymers confined between solid surfaces. Borodin et al.42 have studied the interface between poly(ethylene oxide) and TiO2 surface and concluded that the polymer density increases near the surface. They also reported that the relaxation of poly(ethylene oxide) at TiO2 surface is affected by the density increasing near the surface and by surface topography. More recently, Daoulas et al.43 and Harmandaris et al.44 have performed united-atom MD simulations of thin polyethylene films supported by crystalline graphite. They concluded that the local mobility of the polymer melt near the surface is highly anisotropic. Confined complex molecules under shear have also been simulated to study the viscosity coefficients. In a series of papers, Jabbarzadeh et al.45-49 have done MD simulation of a united-atom hexadecane confined between three layers of atoms of a body-centered-cubic lattice, as the confining surface, to study the effect of surface roughness, hydrocarbon branching, and shear rate on the density and the viscosity coefficient. Khare et al.50 and Jeng et al.51 performed MD simulation of bead-spring polymer models to study the behavior of polymer melts under Couette flow. In most of the afore-cited simulation studies, the authors have tried to simulate the surface force apparatus experiments. In these experiments, a confined liquid is in thermodynamic equilibrium with a surrounding bulk fluid. Consequently, the measured force (or pressure) is more appropriately described as a disjoining force (or pressure), which is the difference between the force (pressure) in the confinement and that of the bulk fluid with which the confined film is in equilibrium.52 The statistical-mechanical ensemble corresponding to the surface force apparatus experiments is the grand canonical ensemble (GCE), in which the chemical potential, volume, and temperature are constant. In the GCE, the solvation force, the force exerted on the confining surfaces by the confined fluid, as a function of distance between the confining surfaces can be calculated directly at a fixed chemical potential. Although there are some reports in the literature on the MD simulation in the GCE,53-55 essentially most of the simulations in this ensemble have been performed using MC simulation.56-59 This is because of the fact that particle insertion or deletion in the GCE simulation fits with the stochastic moves in MC, but they cannot easily be incorporated into the equations of motion in MD simulation. The probabilities for insertion and deletion of entire particles depend on the presence of fluctuations, which generate a cavity in a configuration for particle insertion, or a high energy configuration, from which a particle can be deleted. In the case of big molecules, the probability of insertion becomes very low; as a result, grand canonical ensemble MC or MD simulation methods become very inefficient for simulating polymeric systems. In almost all of the simulation studies addressed above, the simulation is not done in the GCE, because of these problems. In the studies by Yethiraj and Hall,39,40 the authors have tried to calculate the solvation force at constant chemical potential by invoking the superposition approximation of van

J. Phys. Chem. B, Vol. 113, No. 16, 2009 5569

Figure 1. Chemical structure of PA-6,6 oligomer studied in this work.

Megen and Snook,58 in which the density profile of a fluid in a small pore at a fixed chemical potential is estimated by knowing the density profile of the fluid at a single wall at the same chemical potential. However, this method is claimed to be just in qualitative agreement with solvation force calculations using GCE simulations.12 Confined fluids can also be simulated by explicitly including the bulk fluid in the simulation box, with which the confined fluid can establish its equilibrium.25,28,60,61 While this method necessitates simulating many bulk molecules, outside the confined region, the slow relaxation dynamics associated with the motion of fluid molecules into or out of the confined region makes the equilibration rate too slow. Recently, we have reported a simulation method in which only the confined region is simulated at a constant number of confined particles, N, constant surface area, A, constant temperature, T, and constant parallel component of the pressure, P|, which we call NAPT ensemble simulation hereafter.62 In this method, the MD simulation method of Berendsen et al.63 with coupling to an external bath is extended to keep the parallel component of pressure fixed by changing the distance between the confining surfaces. In this Article, we employ our new MD simulation method62 in the NAPT ensemble to study the details of structural, conformational, and dynamical properties of PA-6,6 oligomers confined between graphite surfaces. The motivations for this study stem from: (1) in almost all of the simulation studies of equilibrium thermodynamic and transport properties of confined polymers, as addressed above, a bead-spring model is used to describe the polymer matrix, and simulation studies of confined realistic polymers are scarce; (2) even in some cases where the investigators have tried to simulate realistic polymer chains, they have dealt with flexible polymer chains of relatively simple structure, such as polyethylene43,44 or poly(ethylene oxide),42 but the simulation results on the polymers consisting of stiff chains and/or hydrogen-bonded chains are rare; and (3) in most of the existing simulations in the literature, the fluid in the confined region is not kept in equilibrium with the bulk phase, which necessitates either performing simulations in GCE or explicitly including the bulk fluid in the simulation box. Our new NAPT simulation method62 is shown to reproduce the results of GCE simulation method,12,64-67 which guarantees the existence of equilibrium between the fluid in the confined region and the bulk fluid. Model and Simulation Details Molecular dynamics simulations were performed for ethyland butyl-terminated PA-6,6 trimers, containing 6 amide groups, for which the chemical structure is shown in Figure 1. The oligomers simulated in this work may seem too short to be a representative of a real PA-6,6 chain, but the sluggish dynamics of the polymer in the confined region prevents us from simulating bigger chains. Meanwhile, it is worth considering that each PA-6,6 trimer, studied in this work, consists of 116 atoms, which looks big enough in size to reveal many

5570 J. Phys. Chem. B, Vol. 113, No. 16, 2009

Eslami and Mu¨ller-Plathe

TABLE 1: Force Field Parameters for Polyamide-6,6 and Graphitea Ubond ) (1)/(2)kbond(rij - rij0 )2

bond vibrations i-j

r (nm)

10-3kbond (kJ mol-1 nm-2)

Ca-Cb Cb-Cb Cb-H Ca-N Cb-N N-H Ca-O Cg-Cg

0.152 0.1522 0.109 0.1335 0.1449 0.101 0.1229 0.1420

265.44 259.58 284.70 410.3 282.2 363.42 477.3 258.70

0

0 2 Uangle ) (1)/(2)kangle(θijk - θijk )

angle bendings i-j-k

θ (deg)

kangle (kJ mol-1 rad-2)

Ca/b-Cb-Cb H-Cb-H Ca-Cb-H Cb-Cb-H Cb-Ca-O Cb-Ca-N N-Ca-O Ca-N-Cb Ca-N-Hn Cb-N-Hn Cg-Cg-Cg

111.0 109.5 109.5 109.5 120.0 116.6 120.0 120.0 120.0 120.0 120.0

250 145 210 150 335 293 335 125 125 209 418.7

0

0 Utor ) ∑p(ktor)/(2)[1 - cos(p(φijkl - φijkl ))]

torsions i-j-k-l

periodicity

φ0 (deg)

ktor (kJ mol-1)

Ca-Cb-Cb-Cb

1 2 1 3 1 2 1 2 1 1 2 3

180 90 180 180 180 90 180 90 180 180 90 60

8.1 3.3 2.0 3.6 25.0 40.0 9.0 3.42 35.0 8.2 3.2 3.5

Cb-Cb-Cb-Cb Cb-N-Ca-Cb Ca-N-Cb-Cb O-Ca-N-Hn N-Cb-Cb-Cb N-Cb-Cb-H

0 2 Ud ) (1)/(2)kd(φijkl - φijkl )

harmonic dihedrals φ (deg)

i-j-k-l Cg-Cg-Cg-Cg nonbonded interactions

a

kd (kJ mol-1)

0

0.0

334.95

Unb ) 4εij[((σij)/(rij)) - ((σij)/(rij)) ] + (qiqj)/(4πε0)((1)/(rij) + (εrf - 1)/(2εrf + 1)(rij2 )/(r3c )) 12

6

atom

σ (nm)

ε (kJ mol-1)

q (esu)

H Hn Ca Cb N O Cg

0.25 0.10 0.35 0.35 0.33 0.32 0.3469

0.0657 0.0657 0.360 0.458 0.7118 0.8792 0.276

0.08007 0.28190 0.59730 -0.1506 -0.4057 -0.5479 0.0

Ca, Cb, and Cg represent an amide carbon, an aliphatic carbon, and graphite carbon, respectively. Hn represents an amide hydrogen.

characteristics of the confined polymers. We utilized atomistic MD simulation for better understanding of the structure and dynamics of polymer chains in the confinement. Recently, MD simulations of bulk melts of PA-6,6 20-mers have been performed in this group.68,69 The atomistic force-field parameters for oligomers studied in this work are mainly those published earlier.68 The only difference is that the bond lengths have not been constrained in this simulation, and the equilibrium torsional angles, which were incorrectly reported in Table 2 of ref 68, are corrected in this work. Leaving the bond lengths unconstrained has the effect of reducing the time step to 1.0 fs in the

present work, as compared to 2.0 fs in former simulations.68,69 The force constants for bond vibrations are taken from Cornell et al.70 The force field parameters for PA-6,6 oligomers as well as the corresponding potential energy formula are tabulated in Table 1. The confining surfaces are graphite monolayers. The simulation box is periodic in the x and y directions, but not in the z direction. A snapshot of the simulation cell is shown in Figure 2. The force-field parameters for graphite, tabulated in Table 1, were taken from Bedrov and Smith.71 It is shown71,72 that simulation of graphite with these parameters provides a spacing of 0.34 nm between the graphene sheets, which is in a

Structure of Nanoconfined Polyamide-6,6 Oligomers

J. Phys. Chem. B, Vol. 113, No. 16, 2009 5571 where ∆t is the time step, T0 is the target temperature, T is the actual instantaneous temperature, and τT is the time constant for temperature coupling, which determines the strength of coupling to the bath. Coupling to a constant pressure bath is accomplished according to our new method of MD simulation in the NAPT ensemble. The details of the method are explained elsewhere.62 Here, it is sufficient to say that in this method the parallel component of the pressure, P|, which is defined below (eq 2), is kept fixed and is assumed to be equal to the corresponding bulk pressure.

Figure 2. Snapshot of a simulation box with PA-6,6 oligomers confined between graphene surfaces.

TABLE 2: Description of Systems Simulated in This Worka system

number of oligomers

Ng

S1 S2 S3 S4 S5 S6 S7 (bulk)

75 75 100 100 125 150 250

4300 2232 1972 1500 1500 1144

A (nm ) 〈h〉 (nm) 2

N

17 300 112.633 13 164 58.465 15 544 51.654 14 600 39.291 17 500 39.291 19 688 29.966

0.9931 1.7932 2.3888 3.2027 3.8417 5.9354

a Ng, N, A, and 〈h〉 are the number of C atoms per each graphene sheet, the total number of atoms, surface area (of one surface), and the average surface separation, respectively.

good agreement with experiment.73 The parameters for unlike interactions were determined using Lorentz-Berthelot mixing rules.74 MD simulations were carried out using our simulation package, YASP.75,76 Because we did not constrain the C-C bonds in the graphene sheets, the bonds in PA oligomers were also left unconstrained. Nonbonded interactions were applied fully between atoms interacting through torsion terms (1-4 interactions). All nonbonded interactions were truncated at 0.95 nm with a reaction field correction for the Coulombic interactions.74 The effective dielectric constant was taken to be 5.5.68,69 An atomic Verlet neighbor list was used, which was updated every 15 time steps, and the neighbors were included if they were closer than 1.0 nm. The most severely confined system simulated in this work is composed of 75 PA-6,6 trimers confined between two graphene sheets, each composed of 4300 C atoms with a surface area of 112.6 nm2 and an average separation distance of less than 1.0 nm. To be able to calculate properties of confined oligomers as a function of surface separation, we reduced the surface area of the graphene surfaces and increased the number of confined oligomers. A total of six systems, with varying ratios of polymer/surface atoms, as well as the oligomer bulk, were simulated in this work. The systems studied, as well as the number of polymer molecules, number of surface atoms, and the dimensions of surfaces, are listed in Table 2. In this simulation, the polymer melt in the confined region as well as the surface atoms are simulated in NAPT ensemble.62 For this purpose, the system is coupled to a heat bath at a fixed temperature, T0, according to the method of Berendsen et al.,63 in which a scaling of velocities at each time step from V to λV is done, where the scaling constant, λ, is defined as:

λ)1+

(

∆t T0 -1 2τT T

)

(1)

Pxx + Pyy 2 1 1 ) m V2 + [ 3V i i i 2V

P|| )



∑ ∑ (Xij · Fx,ij + Yij · Fy,ij)+ i

j

∑ ∑ (Xis · Fx,is + Yis · Fy,is) i

s

+

∑ ∑ (Xst · Fx,st + Yst · Fy,st)] s

(2)

t

where Pxx and Pyy are the x- and y-components of pressure tensor, respectively, m is the atomic mass, V is the volume, subscripts i and j show the atoms in the confined region, subscripts s and t stand for the surface atoms, X and Y are the relative distances between particles in the x and y directions, respectively, and Fx and Fy are their corresponding forces. To keep P| fixed, a proportional scaling of the z-coordinates of all particles per time step from z to µz is done. The scaling constant µ is defined as:62

µ ) 1 - β∆t

(P0,| - P|) τP

(3)

where β is the isothermal compressibility, P| is the target value of the parallel component of pressure, and τP is the time constant for pressure coupling, determining the strength of coupling to the barostat. The distance between the confining surfaces varied accordingly during the simulation.62 From the surface force apparatus measurements, it is known that the external field created by the configuration of confined surfaces varies extremely slowly from the center of confined region to the bulk.1,2 Because variation of P| from the bulk fluid to the confined region depends on the gradient of fluid-surface interactions parallel to the surface, as it is already reported,77 P| is essentially the same as the bulk pressure. This point has been validated practically in the NAPT ensemble MD simulations of Wang and Fichthorn,77,78 who extended the constraint method of Evans79 and Evans and Moriss80 to develop a NAPT ensemble MD simulation. It is worth mentioning that in the constraint method of Evans,79,80 the Hamiltonian does not represent a physical system, and the extent of the error introduced is not easily evaluated. However, the results of Wang and Fichthorn77,78 on the confined Lennard-Jones fluid, n-decane, and 2,2-dimethyloctane well reproduce the GCE simulation results.12,64-67 Our new NAPT MD simulation method is the extension of the Berendsen et al.63 method with coupling to an external bath. In this method, the principle of least local perturbation consistent with the required global coupling is used. This approximates the perturbation that would occur in a physical nonequilibrium system.63 Moreover, the effect of coupling can be evaluated and controlled by varying the

5572 J. Phys. Chem. B, Vol. 113, No. 16, 2009

Eslami and Mu¨ller-Plathe

Figure 4. Calculated solvation force as a function of surface separation at T ) 400 K and P| ) 101.3 kPa. The points represent the calculated average solvation force as a function of average surface separation, 〈h〉, for systems S1-S6, and the dashed curve indicates the calculated solvation force based on the superposition approximation for system S6. Figure 3. Local density profiles for systems S1-S6 at T ) 400 K and P| ) 101.3 kPa. The number of oligomers, surface area, and average surface separations are given in Table 1.

coupling strength.62,63 The results of our NAPT ensemble MD simulation have also been validated by comparing the results on the simulation of Lennard-Jones fluid confined between fcc surfaces with GCE simulation results12,64-67 and also with the results generated by Wang and Fichthorn.77,78 The method has also been checked by simulating water confined between graphite surfaces and comparing the results with recent MD simulation results of Choudhury and Pettitt,28 who included the bulk water molecules in the simulation cell. The equality of parallel pressure to bulk pressure in the confined region has also been reported in a recent study by Delhommelle and Cummings.81 The principle of our NAPT ensemble MD simulation is that by keeping fixed the temperature and the parallel pressure, by coupling to a constant temperature heat bath and a constant pressure bath, the equilibrium between the confined region and a bulk reservoir is achieved. Therefore, the simulation method employed in this work simulates the surface force apparatus experiments.1,2 The simulation is easy to perform, especially in the case of solvation force calculations. In this simulation, the temperature was fixed at 400 K using a Berendsen thermostat63 with a temperature coupling time of 0.2 ps. The time step for the leapfrog integration scheme74 was 1.0 fs, and P| was kept constant at 101.3 kPa, with the time constant for pressure coupling of 5.0 ps. The simulations were performed for 5.0 ns to get equilibration and 10.0 ns for data collections. During the production runs, the trajectory frames were recorded every 0.5 ps. Results and Discussion Local Density Profile. The local density profile of PA-6,6 trimers studied in this work can be shown on the basis of the center of mass of the molecules or the position of all atoms in the simulation box. The density profile based on the position of all atoms probably shows more structure. Figure 3 shows the distribution of the mass density of atoms with respect to the distance from the graphene surfaces, z. The position of

surfaces is calculated by averaging the z-coordinates of all surface atoms. The local density profiles were computed by dividing the distance between the confining graphene surfaces, h, into a number of thin slabs and calculating the time averaged density for each slab during the simulation. The results in Figure 3 show that at a constant bulk pressure, P| ) 101.3 kPa, at surface separations of around 6.0 nm the local density shows four peaks in the vicinity of each surface, and the intermediate region does not show a layering effect. In this region, the local density has a mean value corresponding to the bulk liquid density. As the surfaces get closer, the fluid forms well-layered structures, parallel to the surfaces. At separation distances of about 2.0 nm, two peaks are seen in the vicinity of each surface. Further reduction of the distance between the confining graphene surfaces squeezes out a layer of the fluid, and at distances of about 1.0 nm, just one layer is seen in the vicinity of each surface. The layering effect seen in Figure 3 is in accordance with previous reports on the layering of Lennard-Jones fluid,62,77 alkanes,77,78 and polymers.39,40,43,44 Solvation Force. The solvation force is defined as the force exerted by the confined film on the confining surfaces, which is the same as the force that would be required to hold the two surfaces at the corresponding separation. The solvation force is expressed in terms of the pressure components as:

〈fs〉 ) A(〈P⊥〉 - 〈P|〉)

(4)

where fs is the solvation force, A is the surface area, P⊥ is the vertical component of pressure, and brackets indicate the average. The average solvation force per unit area, versus average surface separation, 〈h〉, is shown in Figure 4. It is wellknown that better layered structures with sharper density peaks correspond to more repulsive solvation forces. With increasing distance between the confining surfaces, the ordered layers get disordered, and therefore the solvation force decreases on changing the distance between the confining surfaces and shows an oscillatory behavior. The origin of spatial oscillations in the solvation force lies in the ordering/disordering of the chains in layers. The period of oscillations in the solvation force is thus

Structure of Nanoconfined Polyamide-6,6 Oligomers

J. Phys. Chem. B, Vol. 113, No. 16, 2009 5573

associated with the periodic capability of forming layered structures. A better resolution of the period of the oscillations would require more points, that is, more simulations with the number of confined oligomers and the average slit width varying in smaller steps. This is computationally expensive. However, we can calculate the qualitative behavior of the solvation force oscillations from the density profile of a fluid at a single wall using a superposition approximation. For hard chains between hard walls, the superposition approximation results in a simple expression for the solvation force, given as:39,58

fs ) L(F(z) - Fb) AkT

(5)

where k is the Boltzmann constant, F(z) is the number density at a distance z from one of the walls, Fb is the bulk density, and L is the chain length. We chose the system S6 with maximum separation between the confining surfaces as a system showing single-surface characteristics. Calculating the solvation force, using eq 5, allows us to have an estimate of the qualitative behavior of solvation force oscillations. Meanwhile, we can check the validity of the superposition approximation for polymeric systems, which is shown to be accurate for the force between surfaces in flexible hard chain fluids.82 The calculated results based on eq 5 for fs/A as a function of z are shown in Figure 4 and are compared with our calculations, which are also reported in the same figure as 〈fs〉/A versus 〈h〉. To be more accurate, the constants in eq 5 are fitted with our calculated solvation force for system S6. The results show that the superposition approximation is a crude approximation at small separations between the confining surfaces, but it provides a qualitative estimate of oscillation pattern. Surface Density of Confined Oligomers versus Distance. The layering and the solvation force oscillations can be elucidated more by studying the surface density of confined atoms as a function of distance between the confining surfaces. The amount of material held in a slit pore per surface area of pore is just the integral of the density profile along the direction normal to the surface, z, that is:

n(z) ) A

∫0z F(z) dz

(6)

where n is the number of oligomers. The calculated number of oligomers as a function of distance with respect to the confining surfaces, calculated for systems S1-S3, are shown in Figure 5. Also, we have shown the calculated results for the number of adsorbed oligomers as a function of average surface separations in Figure 5 for the sake of comparison. This is simply the number of chains in the system divided by the area. The results in Figure 5 show that on changing the distance between the confining surfaces, the oligomers adsorb into or desorb out of the confining region. When the oligomers desorb from a layer in a confined region, a new layer will be formed and vice versa. Therefore, a step height in Figure 5 corresponds to one layer, and a step width corresponds to one-layer’s width. A confined fluid with a specified number of particles can gain or lose particles, having the same number of layers, as long as a jump/ step in Figure 5 is not occurred. Comparing the results in Figures 4 and 5 reveals that the solvation force goes up when a new layer is forming and falls when pre-existing layers gain or lose particles.

Figure 5. The surface density of confined PA oligomers versus surface separation. The points represent the surface density versus the average surface separation for systems S1-S6. The full and dashed curves represent the integrated results for systems S5 and S3, respectively.

Figure 6. Oscillations of the number density of confined PA oligomers as a function of surface separation. The points represent the average number density versus the average surface separation for systems S1-S6. The dotted, dashed, and full curves represent the integrated results for system S3, S5, and S6, respectively.

More revealing is the number density of oligomers in the confined region. Plots of n/Az, for systems S3, S5, and S6 in which n is calculated employing eq 6, are shown in Figure 6. Also, the calculated results of n/A〈h〉 versus 〈h〉 are shown in the same figure. According to the results, the number density shows oscillatory behavior, like the pattern observed in the solvation force. The results in Figure 6 show that below a certain distance, the films exhibit certain features characteristics of solids. Moreover, the calculated average values (points) in Figures 5 and 6 for systems with surface separations corresponding to well-formed structures stay higher as compared to predicted results calculated by integrating eqs 5 and 6 for a system showing single-surface characteristics, which indicates that the maxima in number density curve correspond to wellformed structures. The reverse is true for systems with intersurface distances at which the better-layered structures are destroying. It is known that the number density increases when the particles are added to pre-existing layers and decreases when a new layer forms.52,62,77 Orientational Distribution. To investigate in more detail the structure of the confined oligomers, we have computed the

5574 J. Phys. Chem. B, Vol. 113, No. 16, 2009

Eslami and Mu¨ller-Plathe

Figure 7. Orientation distribution function of NH bond vectors with respect to z direction as a function of their distance from the surfaces. The full curves represent the first Legendre polynomial, p1, and the dashed curves represent the second Legendre polynomial, p2.

orientational distribution functions (ODF) as a function of z. The mutual orientation is described by ODF in terms of the first and/or second Legendre polynomials, which are expressed as:

p1(z) ) 〈ui · uj〉

(7)

1 p2(z) ) 〈3(ui · uj)2 - 1〉 2

(8)

and

where ui and uj are unit vectors, and p1 and p2 are the first and second Legendre polynomials, respectively. To calculate the orientations with respect to the surfaces, one of the unit vectors is taken as the surface normal z. The preferential orientations of NH, CO, and end-to-end vectors are calculated. The first Legendre polynomial shows the orientations of these vectors along the z direction and varies between 1, in the case of parallel orientation, to -1, in the case of complete antiparallel orientation with respect to the surface normal. The second Legendre polynomial is an order parameter of the angle between unit vectors and the vector normal to the surfaces. In this case, a positive value indicates parallel alignment, a negative sign indicates perpendicular alignment, and a zero value indicates random alignment with respect to the surface normal.

Shown in Figures 7-9 are the first and second Legendre polynomials for NH, CO, and end-to-end vectors, respectively. As it is seen from the results in Figure 7, NH bonds in close vicinity of the surfaces orient perpendicular to the surfaces. The first Legendre polynomial shows that the NH bonds close to the surfaces direct antiparallel with respect to the surface normal. This type of orientation for NH bonds reinforces hydrogen bonding between amide groups. With increasing distance from the surfaces, the NH bond orientation shows oscillatory behavior. This means that preferential orientation of NH bonds in the first layer (the closest layer to the surfaces) causes similar orientations in the vicinal layers. A close look at the pattern of oscillations in the NH bond orientations in Figure 7 shows that the amplitude of oscillations is related to the same pattern as in the solvation force (see Figure 4), indicating that preferential orientation of NH bond vectors depends on the degree of layering of PA oligomers around the surfaces. In other words, better-layered structures indicate more perpendicular NH bond orientations with respect to the surfaces. The oscillations, however, are damped out at a distance of about 2.0 nm away from the surfaces (see Figure 7, S6). The results in Figure 8 show that CO bond vector orientations show a trend similar to that of the NH bond vectors, but to a smaller degree. Again, in the vicinity of surfaces, the CO bond vectors orient more antiparallel with respect to the surface normal to form hydrogen bonds with NH hydrogens of their neighboring layer. On the other hand, the orientation of end-

Structure of Nanoconfined Polyamide-6,6 Oligomers

J. Phys. Chem. B, Vol. 113, No. 16, 2009 5575

Figure 8. The same as Figure 7 for CO bond vectors.

to-end vectors, shown in Figure 9, expectedly reveals characteristics different from those of NH or CO bond vectors. The results in Figure 9 for systems with small intersurface separations, S1 and S2, show that the end-to-end vectors are aligning perpendicular to the surface normal (parallel to the surfaces) at almost all distances with respect to the confining surfaces. Increasing the distance between surfaces causes more randomness in the end-to-end vectors. In systems S4 and S5, end-toend vectors orient randomly just in the central region. In system S6, the oligomer end-to-end vectors align parallel to the surfaces at short distances and show bulk-like behavior at a distance about 1.5-2.0 nm with respect to the surfaces. The results in this section on the orientation of end-to-end vectors are in agreement with the results of Daoulas et al.43 on the conformational analysis of polyethylene in contact with graphite surfaces. Hydrogen Bonding. The discussions in the preceding section on the orientation of NH and CO bond vectors show that the hydrogen donor and acceptor groups orient in such a way as to form hydrogen bonds. To have a more quantitative analysis of the degree of hydrogen bonding in the confined region and its comparison with a bulk sample, we have analyzed the hydrogen bonding among the amide groups. We have considered amide NH as the hydrogen-bond donor and the amide oxygen as the hydrogen-bond acceptor. The results of previous studies in this group68,69 showed that including amide nitrogen as an acceptor leads to very small increase in the hydrogen bonds; therefore, the amide nitrogen has not been taken into account as a hydrogen-bond acceptor in this work.

Recently, hydrogen bonding in bulk PA-6,6 has been studied in detail in this group.83 The main structural property distinguishing the hydrogen bonds from van der Waals interaction is their preference for linearity. It is shown that the H · · · O bond distance for linear hydrogen bonds, θN-H · · · O > 130°, in bulk PA-6,6 shows a distribution with a distinct maximum around 0.218 nm at different temperatures.83 Increasing the temperature over a range from 300 to 600 K shifts the location of the maximum in the H · · · O distance distribution linearly to distances between 0.212 and 0.224 nm. Therefore, several criteria can be taken into account as the necessary condition for hydrogenbond formation. In their study of hydrogen bonding in bulk PA6,6, Karimi-Varzaneh et al.83 reported a geometric criterion, in which the distance between the hydrogen of the donor group, N, and the acceptor O has to be less than 0.297 nm and the donor-hydrogen-acceptor angle has to be bigger than 130°.69,83,84 In this work, two cases are studied: the donor-hydrogen-acceptor angle has to be above 130°, and the distance between the hydrogen of the donor group and the acceptor has to be less than (a) 0.297 nm and (b) 0.2 nm. Of the two cases studied here, criterion (a) has formerly been adopted for study of hydrogen bonding in PA-6,6, over a wide range of temperatures,83 and criterion (b) is close to the hydrogen-bond length in water. The results of our analysis on hydrogen bonding in confined PA-6,6 oligomers are plotted in Figure 10 for systems S1-S6 as the ratio of the average number of hydrogen bonds per confined PA oligomer, found over the 10 ns of simulation, to the corresponding quantity for bulk PA sample, S7. In all cases,

5576 J. Phys. Chem. B, Vol. 113, No. 16, 2009

Eslami and Mu¨ller-Plathe

Figure 9. The same as Figure 7 for end-to-end bond vectors.

Figure 10. Variation of the ratio of average number of hydrogen bonds per PA-6,6 oligomer to the same quantity in bulk sample as a function of average surface separation. The angle between N-H · · · O is above 130°, and the O · · · H distances are less than (]) 0.297 nm and (4) 0.2 nm.

system S1 shows considerably a smaller number of hydrogen bonds than that of the bulk polymer. This is in complete agreement with previous reports on the hydrogen bonding of

water confined between graphite surfaces.85 The same conclusions have been reported in computer simulations of the static properties of water confined in Vycor.86 The results in Figure 10 show that the ratio of hydrogen bonds with respect to the bulk fluid shows an oscillatory behavior. In the strongly confined cases, S1 and S2, the geometrical restrictions strongly influence the degree of hydrogen bonding. With increasing distance between the confining surfaces, the number of hydrogen bonds increases and shows an oscillatory behavior, with the oscillations damping out at larger slit widths. Choosing a shorter distance as the criterion for hydrogen-bond formation, a similar pattern is observed, but two points are worth considering: (1) The ratio of hydrogen bonds in the confined fluid to that of the bulk polymer increases; and (2) the amount of increase is higher for better-layered structures (see Figure 4, the points at 〈h〉 ) 0.99, 2.4, and 3.8 nm, corresponding to systems S1, S3, and S5, respectively). Increasing the ratio of hydrogen bonds at shorter distances for confined oligomers reveals that in confined region the average H · · · O distance is shorter than that of the bulk fluid. For well-formed structures with more positive solvation forces (Figure 4), the H · · · O distance is shorter with respect to the corresponding confined fluid showing low/negative solvation forces. These interpretations describe the oscillatory behavior of hydrogen-bond ratios observed in Figure 10, as well. The conclusions in this part show that, although the number of

Structure of Nanoconfined Polyamide-6,6 Oligomers

J. Phys. Chem. B, Vol. 113, No. 16, 2009 5577 11 for S3 and S4 shows that better layered structures have the capability of forming more hydrogen bonds. In the case of S1 and S2, the effect of geometrical restriction is more important than the layering effect. The results in this section are in complete agreement with the results of the former section on the perpendicular orientation of NH bond vectors with respect to the surfaces, the oscillatory behavior of NH bond orientations, and the preferential orientation of NH bond vectors depending on the degree of layering of oligomers parallel to the surfaces. Reorientation Correlation Function. The reorientation dynamics of local oligomer segments can be studied by means of the autocorrelation function of the first and/or second Legendre polynomials, that is:

Figure 11. Variation of the ratio of hydrogen bonds per PA-6,6 oligomer to the same quantity in bulk sample as a function of distance with respect to the surfaces.

hydrogen bonds in the confined region reduces due to geometrical restrictions, the hydrogen bonds formed in the confined region have a higher strength than those of corresponding bonds in bulk fluid. To characterize the hydrogen-bond network, we have calculated the ratio of hydrogen bonds per confined PA-6,6 oligomer to the corresponding quantity in bulk sample for bins parallel to the surfaces as a function of distance with respect to the surfaces. The z coordinate of the midpoint of H · · · O bond, subject to criterion (a), was placed into a number of thin bins in the z-direction, and the time-averaged hydrogen bonds were counted. The results are shown for systems S1-S4 in Figure 11. The results in Figure 11 for the variation of hydrogen-bond ratios as a function of surface distance show that the position of the first (side) peaks, characterizing hydrogen bonds in closest layer to each surface, shows a small shift as compared to the position of the side peaks in the density profile of PA-6,6 oligomers, shown in Figure 3. The position of the remaining peaks in Figure 11 matches with the position of minima in the density profile graphs, Figure 3. The results indicate that hydrogen bonds are forming between polymer layers, which mean that the layers orient in such a way to form hydrogen bonds with their neighboring layers. PA-6,6 oligomers very close to the surfaces form hydrogen bonds with each other, resulting in the side peaks in Figure 11, and with the inner layer, resulting in the second peaks in Figure 11. The inner layers form hydrogen bonds with their neighboring layers, which leads to the hydrogen-bond peaks corresponding to the position of minima in the density profile peaks. In the case of S1, the hydrogen-bond ratios show two peaks in the vicinity of each surface, and their positions are displaced about the position of maxima in density profile peaks in Figure 3. This supports the idea of hydrogen-bond formation between PA-6,6 oligomers for layers close to the surfaces. In the case of S4, the inner peaks oscillate about one, indicating that the oligomers in the middle region of the box show bulk-like behavior. Comparison of the intensities of side peaks with inner peaks in Figure 11 shows that the closest layer of fluid to the surfaces forms less hydrogen bonds as compared to its neighboring layer, due to the geometrical restrictions. The closest layer to the surfaces mostly forms hydrogen bond with its neighboring layer. The inner layers form hydrogen bonds with their vicinal layers, and the intensity of peaks varies with the degree of layering of oligomers parallel to the surfaces (Figure 3). Comparison of the results in Figure

〈p1(t)〉 ) 〈u(t) · u(0)〉

(9)

1 〈p2(t)〉 ) 〈3(u(t) · u(0))2 - 1〉 2

(10)

and

Here, the reorienting vectors, u, are the NH, CO, and end-toend vectors, which provide key information about dynamics of amide and end groups. The reorientation of amide groups, resulting in intermolecular hydrogen-bond creation or breaking, is usually associated in the literature with the γ-process observed in dielectric spectra.87,88 The results for 〈p1(t)〉 are shown in Figure 12. Comparison of the results in Figure 12 for NH and CO bond vector correlation functions shows that in highly confined systems, for example, S1 and S2, NH and CO bond vectors show different relaxation times, but the correlation curves get more alike with increasing surface separations. Therefore, the results for NH and CO bond vectors are averaged together in the case of systems S5 and S6. This is in agreement with the results of a previous study in this group on bulk PA66 sample, which showed that NH and CO bond vectors had very similar autocorrelation functions.69 Our results also show that the decay rate of the first Legendre polynomial depends on the surface separation and on the layering effect. We used arbitrarily the values of 〈p1(t)〉 at 8000 ps as a measure of the amount of decorrelation achieved for NH and CO bond vectors. Using another time does not change the conclusions made in this section. Shown in Figure 13 are the results for 〈p1 (8000 ps)〉 as a function of surface separation for systems S1-S6. The results show an oscillatory behavior confirming that the reorientation correlation times depend on the degree of layering and on the surface separation. The oscillatory behavior seen in Figure 13 is similar to the trend seen in Figure 4, for solvation force oscillations; that is, better ordered structures show higher relaxation times. To obtain the position-dependent (along the z-axis) correlation functions, we have divided the distance between the confining surfaces into a few relatively narrow bins and calculated the reorientation correlation function of CO bond vectors over the 10 ns span of simulation. Symmetrically equivalent bins are averaged because of the symmetry, and the resulting first Legendre polynomials are plotted in Figure 14, in which the distance between the confining surfaces is divided into 4 bins, in the case of S1, to 8 bins, in the case of S6. The relaxation times decrease with increasing distance from the surface. The results for systems S4-S6 show that the central slabs show relatively close relaxation times, meaning that the fluid in the central region of the box shows bulk-like behavior.

5578 J. Phys. Chem. B, Vol. 113, No. 16, 2009

Eslami and Mu¨ller-Plathe

Figure 12. Reorientation correlation functions for NH, CO, and end-to-end bond vectors: full curves, NH; dashed curves, CO; and dotted curves, end-to-end vectors.

We have again selected the value of reorientation correlation function at t ) 8000 ps as a measure of decorrelation and compared the ratios of p1 (8000 ps) values for bins with approximately the same distance to the surfaces as well as the ratios of 〈p1 (8000 ps)〉, averaged over the whole slit width, in Table 3. The results in Table 3 show that the bins with approximately the same distance to the surfaces have different relaxation times, and an oscillatory behavior is seen in the ratios of relaxation times, corresponding to the degree of layering. The reorientation correlation function ratio for bins at relatively the same distance to the surfaces in some cases, like the first bins of S3 and S2, both locating at a distance of about 0.2 nm away from the closest surface, is nearly the same as the corresponding average ratio, averaged over the whole intersurface separations. This indicates that, although relaxation times change with the distance to the surfaces, in this separation range the ratios of relaxation times in the central region are nearly the same as the corresponding ratio for layers close to the surfaces. On the other hand, the reorientation correlation function

Figure 13. The ratio of 〈p1 (8000 ps)〉 for systems S1-S6 as a function of average separation distance. The filled and open markers represent the results for CO and end-to-end vectors, respectively.

Structure of Nanoconfined Polyamide-6,6 Oligomers

J. Phys. Chem. B, Vol. 113, No. 16, 2009 5579

Figure 14. Local reorientation correlation functions for CO bond vectors, calculated in a few slabs parallel to the surfaces, as a function of distance from the surfaces. The slower decaying curves correspond to slabs with shorter distances to the surfaces.

TABLE 3: Comparison of the Ratios of p1 (8000 ps), for Bins Located at Approximately the Same Distance to the Surfaces, to the Ratios of the Same Quantity Averaged over the Whole Slit Width systems (ratio)

ratio of bin distance to the surface (nm)

p1 (8000 ps) ratio

〈p1 (8000 ps)〉 ratio

S3/S2 S4/S2 S5/S2 S3/S4 S3/S4 S1/S6

0.202/0.224 0.200/0.224 0.240/0.224 0.202/0.200 0.606/0.600 0.372/0.371

1.805 0.990 1.591 1.832 1.853 3.051

1.752 0.401 1.153 4.271 4.271 7.490

ratio for the second bin of S1 and the first bin of S6, both located at around 0.37 nm away from the closest surface, is quite different from the ratio of corresponding average value for the same systems, due to the faster dynamics of the polymer melt in the central region of box in S6. The correlation time ratio for the second bins of S3 and S4, located around 0.6 nm of the closest surface, is approximately the same as the corresponding ratio for their first bins, but the result is different from the corresponding average quantity. This indicates that in the case of S4 with larger separation between the surfaces the bins at

the same distance to the surfaces show shorter relaxation times as compared to S3, and the trend continues to distances of about 0.6 nm to the surfaces, but system S4 shows on average much shorter relaxation times, because of the effect of central region of the box. The conclusions in this section are in agreement with previous results by Borodin et al.42 on the simulation of poly(ethylene oxide) near TiO2 surfaces and by Daoulas et al.43 and Harmandaris et al.44 on the simulations of thin polyethylene films near graphite surfaces. They concluded that the relaxation time is affected by density increasing near the surfaces.

5580 J. Phys. Chem. B, Vol. 113, No. 16, 2009

Figure 15. Center-of-mass mean-square displacement for system S1.

Diffusion Coefficients. The mobility of confined PA-6,6 oligomers in a plane parallel to the confining graphene surfaces provides information about the transport properties of the confined melt. The parallel component of the diffusion coefficient, defined as the slope of the center-of-mass mean-square displacement, reads as:

〈D|〉 ) 〈(xcm(t) - xcm(0))2〉 + 〈(ycm(t) - ycm(0))2〉 1 (11) lim 4 t f∞ t where 〈D|〉 is the parallel component of the center-of-mass diffusion coefficient. Because a fluid in direct contact with a surface shows very slow dynamics perpendicular to the surfaces, we did not calculate the vertical component of the diffusion coefficient. Our calculated center-of-mass mean-square displacement as a function of time for system S1, as a typical example, is shown in Figure 15. Averaging is performed over all oligomers and all possible time origins. This leads to an increase of statistical fluctuations toward the end of the graph. The meansquare displacement is linear in time, confirming Einstein diffusion. The linear portion of the mean-square displacement is least-squares fitted by a straight line to calculate the parallel diffusion coefficient (Figure 16). The correlation coefficients of the fits are close to unity. The ratio of the average parallel diffusion coefficient to its corresponding bulk value, D0 ) 3.75 × 10-8 cm2 s-1, gives additional evidence for a confinementinduced slowing of the film’s mobility. The ratios 〈D|〉/D0 show oscillatory behavior with the slit width, which follows the solvation force pattern (Figure 4). The minima of 〈D|〉/D0 ratios are at surface separations corresponding to well-formed layers, and the maxima are at surface separations corresponding to the local minima in the solvation force, that is, when well-formed layers cannot be formed. Therefore, disordered structures allow larger diffusivities, while ordering reduces them. In other words, the higher density of ordered structures leaves a limited free volume for fluid molecules and thus decreases translational diffusivity.61,89,90 This is in agreement with the results observed for confined Lennard-Jones52 fluid and alkanes,78,89,90 and also with our above-cited results on the oscillatory behavior of autocorrelation function ratios. Conclusions A new NAPT ensemble simulation for confined fluids is employed to simulate PA-6,6 oligomers confined between

Eslami and Mu¨ller-Plathe

Figure 16. Average parallel component of diffusion coefficient for systems S1-S6 as a function of average surface separation.

graphite surfaces. It is assumed that the parallel component of pressure is uniform throughout the confining region and is equal to the bulk pressure; the distance between the surfaces is varied to get the target value of the parallel component of pressure. Our results show that oligomers form layered structures parallel to the surfaces, and the number density of oligomers shows a stepwise variation as a function of separation between the confining surfaces. The method allows an easy calculation of the solvation force. Like for Lennard-Jones and alkane liquids,52,62,77,78 the solvation force shows oscillations with the slit width, with maxima at surface separations corresponding to well-formed layers. Analysis of end-to-end vector alignments shows that the oligomers adopt a flattened conformation near the surfaces, whereas at distances about 1.5-2 nm from the surface the end-to-end vectors show bulk-like behavior. Despite geometrical restriction in the confined region, the NH and CO bond vectors orient in a way to form a hydrogen-bond network. The degree of hydrogen-bond formation is governed by layering and geometrical restrictions. It reveals an oscillatory behavior, like the solvation force oscillations. In highly confined systems, the hydrogen bonds weaken due to geometrical restrictions, but in well-formed structures at some intermediate distances the hydrogen bonding can even strengthen in comparison with the bulk fluid. The orientational relaxation times of NH, CO, and end-to-end vectors also show oscillatory behavior with maxima corresponding to well-formed structures parallel to the surfaces. Calculation of local relaxation times shows that the fluid layers close to the surfaces have higher relaxation times than the fluid in the central region. Our results on the center-of-mass diffusion coefficient of confined oligomers show that the parallel component of diffusion coefficient shows oscillatory behavior with local minima at distances corresponding to well-formed structures and local maxima at distances corresponding to distances at which the well-formed layers cannot develop due to incommensurability with the slit width. Acknowledgment. The support of this work by the Excellence Cluster Center of Smart Interfaces at the Technische Universita¨t Darmstadt is gratefully acknowledged. References and Notes (1) Bhushan, B.; Israelachvili, J. N.; Landman, U. Nature (London) 1995, 374, 607. (2) Israelachvili, J. N.; McGuiggan, P. M.; Homola, A. M. Science 1988, 240, 189.

Structure of Nanoconfined Polyamide-6,6 Oligomers (3) Huang, S. H.; Li, C. L.; Hu, C. C.; Tsai, H. A.; Lee, K. R.; Lai, J. Y. Desalination 2006, 200, 387. (4) Roh, I. J.; Park, S. Y.; Kim, J. J.; Kim, C. K. J. Polym. Sci., Part B: Polym. Phys. 1998, 36, 1821. (5) Kim, I. C.; Jegal, J.; Lee, K. H. J. Polym. Sci., Part B: Polym. Phys. 2002, 40, 2151. (6) Brycer, M.; Nguyen, H. T.; Nakeeran, P.; Clement, T.; Hauhen, C. J.; Tykwinski, R. R.; DeCorby, R. G.; McMullin, J. N. Thin Solid Films 2004, 458, 233. (7) Davis, H. T. Statistical Mechanics of Phases, Interfaces and Thin Films; VCH Publishers: New York, 1996. (8) Babuena, P. B.; Berry, D.; Gubbins, K. E. J. Phys. Chem. 1993, 97, 937. (9) Walley, K.; Schweizer, K. S.; Peanasky, J.; Cai, L.; Granick, S. J. Phys. Chem. 1994, 100, 3361. (10) Urbakh, M.; Daikhin, L.; Klafter, J. J. Chem. Phys. 1995, 103, 10707. (11) Rozman, M. G.; Urbakh, M.; Klafter, J. Phys. ReV. Lett. 1996, 77, 683. (12) Snook, I. K.; Van Megan, W. J. Chem. Phys. 1980, 72, 2907. (13) Cleveland, C. L.; Landman, U.; Bamett, R. N. Phys. ReV. Lett. 1982, 49, 790. (14) Magda, J. J.; Tirrell, M.; Davis, H. T. J. Chem. Phys. 1985, 83, 1888. (15) Abraham, F. F. AdV. Phys. 1986, 35, 1. (16) Bitsanis, I.; Magda, J. J.; Tirrell, M.; Davis, H. T. J. Chem. Phys. 1987, 87, 1733. (17) Schoen, M.; Rhykerd, C. L.; Diestler, D. J.; Cushman, J. H. Science 1989, 245, 1223. (18) Thompson, P. A.; Robbins, M. O. Phys. ReV. A 1990, 41, 6830. (19) Bitsanis, J.; Somars, S. A.; Davis, H. T.; Tirrell, M. J. Chem. Phys. 1990, 93, 3427. (20) Thompson, P. A.; Robbins, M. O. Science 1990, 250, 792. (21) Robbins, M. O.; Thompson, P. A. Science 1991, 253, 916. (22) Ribarsky, M. W.; Landman, U. J. Chem. Phys. 1992, 97, 1937. (23) Gupta, S.; Koopman, D. C.; Westerman-Clark, G. B.; Bitsanis, I. A. J. Chem. Phys. 1994, 100, 8444. (24) Ouyang, J.; Luedtke, W. D.; Landman, U. Bull. Am. Phys. Soc. 1995, 40, 425. (25) Gao, J.; Luedtke, W. D.; Landman, U. J. Chem. Phys. 1997, 106, 4309. (26) Curry, J. E. J. Chem. Phys. 2000, 113, 2400. (27) Cui, S. T.; Cummings, P. T.; Cochran, H. D. J. Chem. Phys. 2001, 114, 7189. (28) Choudhury, N.; Pettitt, B. M. J. Am. Chem. Soc. 2005, 127, 3556. (29) Marti, J.; Nagy, G.; Gordillo, M. C.; Guardia, E. J. Chem. Phys. 2006, 124, 094703. (30) Dokmeci, M.; Najafi, K. IEEE/MEMS 2001, 10, 197. (31) Keddie, J. L.; Jones, R. A. L.; Corey, R. A. Europhys. Lett. 1994, 27, 59. (32) Kurabayashi, K.; Asheghi, M.; Touzelbaev, M.; Goodson, K. E. IEEE/MEMS 1999, 8, 180. (33) Hu, H.-W.; Granick, S. Science 1992, 258, 1339. (34) Wyart, F. B.; deGennes, P.-G. Eur. Phys. J. E. 2000, 1, 93. (35) Manias, E.; Subbotin, A.; Hadziioannou, G.; ten Brinke, G. Mol. Phys. 1995, 85, 1017. (36) Mansfield, K. F.; Theodorou, D. N. Macromolecules 1989, 22, 3143. (37) Bitsanis, I. A.; Pan, C. J. Chem. Phys. 1993, 99, 5520. (38) Manias, E.; Hadziioannou, G.; ten Brinke, G. Langmuir 1996, 12, 4587. (39) Yethiraj, A.; Hall, C. K. Macromolecules 1990, 23, 1865. (40) Yethiraj, A. J. Chem. Phys. 1994, 101, 2489. (41) Aoyagi, T.; Takimoto, J.; Doi, M. J. Chem. Phys. 2001, 115, 552. (42) Borodin, O.; Smith, G. D.; Bandyopadhyaya, R.; Byutner, O. Macromolecules 2003, 36, 7873. (43) Daoulas, K. Ch.; Harmandaris, V. A.; Mavrantzas, V. G. Macromolecules 2005, 38, 5780. (44) Harmandaris, V. A.; Daoulas, K. Ch.; Mavrantzas, V. G. Macromolecules 2005, 38, 5796. (45) Jabbarzadeh, A.; Harrowell, P.; Tanner, R. I. Phys. ReV. Lett. 2005, 94, 126103.

J. Phys. Chem. B, Vol. 113, No. 16, 2009 5581 (46) Jabbarzadeh, A.; Atkinson, J. D.; Tanner, R. I. Phys. ReV. E 2000, 61, 690. (47) Jabbarzadeh, A.; Atkinson, J. D.; Tanner, R. I. J. Non-Newtonian Fluid Mech. 1998, 77, 53. (48) Jabbarzadeh, A.; Atkinson, J. D.; Tanner, R. I. Comput. Phys. Commun. 1997, 107, 123. (49) Jabbarzadeh, A.; Atkinson, J. D.; Tanner, R. I. J. Chem. Phys. 1999, 110, 2612. (50) Khare, R.; de Pablo, J. J.; Yethiraj, A. Macromolecules 1996, 29, 7910. (51) Jeng, Y.; Chen, C.; Shyu, S. Tribol. Lett. 2003, 15, 293. (52) Gao, J.; Luedtke, W. D.; Landman, U. J. Phys. Chem. B 1997, 101, 4023. (53) Eslami, H.; Mu¨ller-Plathe, F. J. Comput. Chem. 2007, 28, 1763. (54) Eslami, H.; Mu¨ller-Plathe, F. Macromolecules 2007, 40, 6413. (55) Lynch, C. J.; Pettitt, B. M. J. Chem. Phys. 1997, 107, 8594. (56) Rowley, L. A.; Nicholson, D.; Parsonage, N. G. J. Comput. Phys. 1975, 17, 401. (57) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187. (58) Van Megen, W.; Snook, I. K. J. Chem. Soc., Faraday Trans. 2 1979, 75, 1095. (59) Yao, J.; Greenkorn, R. A.; Chao, K. C. Mol. Phys. 1982, 46, 587. (60) Wang, Y.; Hill, K.; Harris, J. G. J. Chem. Phys. 1993, 100, 3276. (61) Gao, J.; Luedtke, W. D.; Landman, U. Phys. ReV. Lett. 1997, 79, 705. (62) Eslami, H.; Mozaffari, F.; Moghadasi, J.; Mu¨ller-Plathe, F. J. Chem. Phys. 2008, 129, 194702. (63) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (64) Schoen, M.; Diestler, D. J.; Cushman, J. H. J. Chem. Phys. 1994, 100, 7707. (65) Schoen, M.; Diestler, D. J.; Cushman, J. H. Phys. ReV. B 1993, 47, 5603. (66) Schoen, M.; Diestler, D. J.; Cushman, J. H. J. Chem. Phys. 1987, 87, 5464. (67) van Megen, W. J.; Snook, I. K. J. Chem. Phys. 1981, 74, 1409. (68) Goudeau, S.; Charlot, M.; Vergelati, C.; Mu¨ller-Plathe, F. Macromolecules 2004, 37, 8072. (69) Goudeau, S.; Charlot, M.; Mu¨ller-Plathe, F. J. Phys. Chem. B 2004, 108, 18779. (70) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (71) Bedrov, D.; Smith, G. D. Langmuir 2006, 22, 6189. (72) Li, L.; Bedrov, D.; Smith, G. D. J. Phys. Chem. B 2006, 110, 10509. (73) Nicklow, R.; Wakabayashi, N.; Smith, H. G. Phys. ReV. B 1972, 5, 4951. (74) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (75) Mu¨ller-Plathe, F. Comput. Phys. Commun. 1993, 78, 77. (76) Tarmyshov, K.; Mu¨ller-Plathe, F. J. Chem. Inf. Mod. 2005, 45, 1943. (77) Wang, J. C.; Fichthorn, K. A. J. Chem. Phys. 2000, 112, 8252. (78) Wang, J. C.; Fichthorn, K. A. Colloids Surf., A 2002, 206, 267. (79) Evans, D. J. J. Chem. Phys. 1983, 78, 3297. (80) Evans, D. J.; Morriss, G. P. Comput. Phys. Rep. 1984, 1, 297. (81) Delhommelle, J.; Cummings, P. T. Phys. ReV. B 2005, 72, 172201. (82) Yethiraj; Hall, C. K. Mol. Phys. 1991, 73, 503. (83) Karimi-Varzaneh, H. A.; Carbone, P.; Mu¨ller-Plathe, F. Macromolecules, 2008, 41, 7211. (84) Garcia, D.; Starkweather, H. W. J. Polym. Sci., Polym. Phys. Ed. 1985, 23, 537. (85) Rovere, M.; Gallo, P. Eur. Phys. J. E 2003, 12, 77. (86) Rovere, M. J. Chem. Phys. 2002, 117, 369. (87) Le Huy, H. M.; Rault, J. Polymer 1994, 35, 136. (88) Laredo, E.; Grimau, M.; Sanchez, F.; Bello, A. Macromolecules 2003, 36, 9840. (89) Bachl, F.; Ludemann, H.-D. Z. Naturforsch. 1986, 41a, 963. (90) Padilla, P. J. Chem. Phys. 1995, 103, 2157.

JP8112655