ARTICLE pubs.acs.org/JPCB
Comparative Study of Commonly Used Molecular Dynamics Force Fields for Modeling Organic Monolayers on Water Michael B. Plazzer, David J. Henry,† George Yiapanis, and Irene Yarovsky* Applied Physics, School of Applied Sciences, RMIT University, GPO Box 2476 V, Victoria, 3001, Australia ABSTRACT: This study compares the performance of the allatom molecular dynamics force fields OPLS-AA and COMPASS, and the united-atom GROMOS96 ff53a6 force field, for organic monolayers at aqueous interfaces, as a function of surface density, temperature, and system size. Where possible, comparison with experimental data was undertaken and used to scrutinize the performance of each force field. We find close agreement between the all-atom force fields (OPLS and COMPASS) and experiment for the description of organic monolayers on water. However, the united-atom force field 53a6 tends to exhibit poorer agreement than the all-atom force fields.
’ INTRODUCTION Naturally occurring and synthetic monolayer/bilayer systems have been an active area of research for decades. Mono- and bilayers of phospholipids form biological membranes that are used by nature to define a chemical or biochemical environment that differs from its surroundings. Considerable research has been directed to understand how these interfaces function within biological systems and how this can be modified for targeted drug delivery and treatment of various medical conditions. For example, the biocompatible and biodegradable composition of phospholipid-based liposomes has led to their use in research and indeed clinical practice as safe, versatile, and effective drug delivery agents for an increasing number of illnesses.1 At the same time, synthetic monolayers are becoming increasingly important materials for industry in the areas of sensor arrays2 and corrosion mitigation3 and as potential evaporation suppressants.4 There has also been increasing interest in the use of synthetic monolayers for industrial applications such as boundary layer lubricants5,6 for magnetic storage devices7,8 and microelectrochemical systems.9 Molecular dynamic simulation has been a useful tool in the study of biological monolayers/membranes systems,10,11 allowing microscopic insight into how these materials aggregate at interfaces in biological systems. Lindahl and Edholm performed molecular dynamic simulations of dipalmitoylphosphatidylcholine (DPPC) bilayers over a range of temporal and spatial scales and observed mesoscopic phenomena such as undulatory and thickness fluctuation (peristaltic) modes.12 This allowed them to calculate mesoscopic continuum properties that compared well with experiment and that are important properties of biological cells, e.g., the ability of cells to change shape and their interaction with other membranes. Self-assembled monolayers of bistable2 rotaxanes on Au (111) surfaces have been modeled by Jang et al.13 using molecular dynamics for their application as molecular machines.14 They determined optimal surface densities of these monolayers in both metastable and ground states which corresponded, respectively, to ON and OFF states of an r 2011 American Chemical Society
electronic memory device. These monolayer systems and others15,16 have benefited from molecular dynamics modeling; however, to the best of our knowledge this is the first comparative, fully atomistic modeling study of organic monolayers at the air/water interface for industrial use. Previously we have applied the COMPASS17 force field to model organic monolayers18 on water, with the purpose of relating monolayer structure to water evaporation resistance. The COMPASS force field has been developed specifically for the study of polymers in the condensed phase.19 Our results indicated that COMPASS is capable of reproducing the basic structural properties of water-based organic monolayers. However, due to high computational demand associated with the use of COMPASS, our simulations were limited in size and time. It is therefore of interest to explore how accurately the more computationally efficient MD algorithms using the GROMOS and OPLS force field sets describe nonbiological organic monolayers at the airwater interface. Comparative studies of these types of force fields can guide force field selection for subsequent studies. Todorova et al.20 investigated OPLS and GROMOS force fields, among others, for the most accurate representation of the experimentally observed characteristics of some peptides. Rog et al.21 assessed a number of different force field configurations with the carbohydrate-optimized OPLS22 force field and found that it reproduced the phase and surface density properties of glycolipids. In this study, we compare the performance of the all-atom OPLS and united-atom GROMOS96 53a6 force fields against COMPASS for the well-characterized organic monolayer composed of octadecanol.
’ COMPUTATIONAL DETAILS Two systems were constructed by placing 20 and 80 octadecanol molecules aligned parallel to the z-direction in periodic unit Received: December 8, 2010 Revised: February 14, 2011 Published: March 21, 2011 3964
dx.doi.org/10.1021/jp1116867 | J. Phys. Chem. B 2011, 115, 3964–3971
The Journal of Physical Chemistry B cells with dimensions 30 30 140 Å and 60 60 140 Å, respectively. The cell size in the z-direction was chosen to accommodated a 20 Å water layer situated directly beneath the octadecanol’s hydrophilic headgroups and a ∼100 Å vacuum space adjacent to the hydrophobic tails to eliminate interaction of the chains with water in the z-direction due to 3D periodic boundary conditions. A quasi 2-D interface was therefore constructed. Three force fields were employed to model the interatomic potentials of the system: the all-atom force fields COMPASS17 and OPLS,23 and the united-atom GROMOS96 ff53a6.24 The single point charge (SPC) model was used to represent the water molecules25 in the OPLS and 53a6 systems. A modified SPC water model was used with the COMPASS systems.26 We also performed additional simulations with the OPLS force field using the TIP4P27 model. The particle mesh Ewald method28 was employed to calculate the electrostatic interactions, with the real space part of the Ewald sum and the Lennard-Jonesvan der Waals interaction cutoff at 9 Å. To eliminate any potential strain inherent in the initial system, the steepest descent method was used to minimize the potential energy with the convergence criterion set to 1 kJ mol1 Å1. Following this, molecular dynamics was performed in the NVT ensemble, annealed for 0.5 ns to overcome local energy minima before an equilibration of a further 0.5 ns. The data acquisition was performed over 10 ns with an integration time step of 1 fs. The temperature was maintained using the Berendsen thermostat29 for simulations at 273, 298, and 368 K. The annealing occurred at 323 and 348 K prior to the target 273 and 298 K simulations, respectively, while the 368 K systems were not annealed due to the already high target temperature. In order to assess the performance of the force fields over a range of monolayer surface densities, multiple systems were generated in a stepwise process as follows: from an initial system 30 30 140 Å the cell size was decreased by 4 Å each in the xand y-directions, the system annealed and equilibrated (0.5 ns), and then data acquired at the new surface density for 10 ns. These steps were repeated until the system size reached 22 22 140 Å at which point the step size was reduced to 2 Å for one step, then to 0.5 Å steps. Data was acquired until 19 19 140 Å at which point signs of monolayer collapse were observed. The following properties were collected and examined across different systems and conditions; specific interactions were characterized by averaging the radial distribution functions g(r) over 10 ns of the trajectory for different interaction pairs. Values of g(r) were normalized to account for the vacuum spacer introduced in the cell to mimic a 2D interface of the monolayer and water. The g(r) for carbon to carbon excludes contributions between bonded carbons within the same chain. Key interactions between chain headgroups were characterized by determining changes in the average number of H-bonds per headgroup as a function of surface density and through changes in the oxygenoxygen RDFs for headgroupheadgroup and headgroupwater interactions. We also elucidated the system’s structural properties; these include the tilt angle which was calculated as the radial angle between the aliphatic chains and the normal to the interface, and the length of the aliphatic carbon chain (C1C18) was measured as a function of the monolayer surface density. The octadecanol headgroup conformation was characterized by the dihedral angle O1C1C2C3, and how this changed with temperature and system size. The analysis tools used for OPLS and 53a6 were from the GROMACS software package.30 The COMPASS systems were analyzed using Accelrys Materials Studio.
ARTICLE
The results are compared with experimental work, and a brief discussion of factors associated with monolayer stability follows. The analysis is composed of three sections; the three force fields are compared against each other and against experiment where it applies. Following this, the behavior of the systems at different temperatures was compared, and then the effect of system size was analyzed by comparing 20 and 80 chain systems. Unless otherwise stated, the surface density of the monolayers presented in the Results and Discussion section is 20 Å2 molecule1. This density was found to be the highest packing density at which none of the three force field systems showed signs of monolayer collapse. We found the time associated with running the simulations to be approximately 2 orders of magnitude longer with COMPASS systems compared to the equivalent OPLS and ff53a6 systems. The functional forms used in the COMPASS force field include several cross-coupling terms that are absent from OPLS and ff53a6. In addition, the equivalent bond and angle functional forms are of a higher order. The modified SPC water model used by COMPASS includes flexible bonds, which were absent in TIP4P and standard SPC. We believe these factors led to the increased computational efficiency in OPLS and ff53a6 force fields.
’ RESULTS AND DISCUSSION We present results obtained by OPLS and GROMOS force fields for this study and compare these with the COMPASS results we published previously for the same system.18 The OPLS with SPC water model was found to give very similar results to the OPLS with TIP4P water, consequently we present the OPLS/TIP4P results unless otherwise specified. We begin our assessment of the performance of the three force fields for the simulation of octadecanol monolayers at constant temperature (298 K) by comparing the monolayer/water interface profiles as a function of surface density. Octadecanol molecules were initially orientated perpendicular to the water interface at a surface coverage of 45.0 Å2 molecule1, with the hydroxyl head groups partially immersed into the water layer. Figure 1 illustrates typical equilibrium structures of watermonolayer interfaces obtained by MD simulations at different surface densities. Visual observation of the monolayer structures suggests that at lower surface densities (Figure 1ac), tilted and more disordered thin monolayers are formed. At higher surface
Figure 1. Equilibrium structures of octadecanol monolayer/water interface as modeled by OPLS, at the surface densities (a) 45.00 Å2 molecule1, (b) 33.80 Å2 molecule1, (c) 24.20 Å2 molecule1, (d) 20.00 Å2 molecule1, (e) 19.01 Å2 molecule1, and (f) 18.05 Å2 molecule1. 3965
dx.doi.org/10.1021/jp1116867 |J. Phys. Chem. B 2011, 115, 3964–3971
The Journal of Physical Chemistry B
ARTICLE
Figure 2. Monolayer properties at 298 K: (a) OPLS concentration profiles along z-axis of water and monolayer at 20 Å2 molecule1; (b) average tilt angle as a function of surface density for the three force fields with standard deviations included; and (c) carboncarbon radial distribution functions (RDF) for 20 Å2 molecule1.
densities (Figure 1e,f), we found that more ordered, perpendicularly arranged and thick monolayers are formed. Figure 2a displays the concentration profiles at different surface densities, attained with the OPLS force field, though the GROMOS force field concentration profile is nearly indistinguishable. The profiles demonstrate the relationship between surface density of the monolayer, the layer thickness in the z-direction, and packing order. Water being an incompressible fluid maintains a constant density, as is demonstrated in Figure 2a by the water profile plateauing at ∼0.12 atoms Å3. The monolayers at relatively low surface density have a lower thickness. This can be seen in the concentration profiles by a comparatively narrower band for decreased surface densities. As the surface density increases, the chain molecules become more ordered and the thickness of the monolayer increases. At the high surface density of ∼19 Å2 molecule1 a well-ordered uniform layer is achieved. This is shown in Figure 2a where the concentration profile exhibits a uniform density across the extent of the monolayer region. Generally, the concentration profiles obtained with OPLS, GROMOS, and COMPASS18 force fields display similar behavior as a function of surface density. However, there are slight differences in the optimum packing densities predicted by the three force fields. For the all-atom force fields, the optimum packing densities are defined as those densities that immediately precede the solidlike phase formation, we defined this solidlike phase as that density where characteristic monolayer properties (tilt, RDFs, H-bonding, chain length) did not change with a further increase in surface density. The united-atom GROMOS force field displayed this behavior, though the determining factor for optimum surface density was the occurrence of a chain ejection below the surface density of ∼18.05 Å2 molecule1. Figure 2b displays the average of the distribution of the chains’ radial tilt angle as a function of surface density obtained via the different force fields. The standard deviation (SD) of the radial angle distribution is given; in the lower surface density region the radial tilt of the chains is stable with little or no fluctuation, and this is reflected in the SD. At ∼20.00 Å2 molecule1 the SD increases for all force fields; this is a consequence of the fluctuating tilt around the axis, and the SD reduces again for the lowest surface densities as the fluctuations decrease. The general trend is the same for all three force fields with a gradual decrease in tilt as the surface density is initially increased from lower surface densities (>24.2 Å2 molecule1). In the higher surface density region the tilt angle decreases more rapidly until 19.01 Å2 molecule1 after which the change in angle is significantly reduced, suggesting a solidlike phase formation. At 18.05
Å2 molecule1 all three force fields agree at a value of ∼2°, i.e., essentially perpendicular chains with respect to the water surface. At low surface densities OPLS and 53a6 predict a similar degree of tilt of the alkyl chains of the monolayer molecules with COMPASS giving slightly lower values. However, as the surface density increases, the variation between the three force fields generally decreases. Our observations are in qualitative agreement with experiment; Lautz et al.31 measured the chain tilt angle of an octadecanol monolayer as a function of the surface pressure at 9 and 20 °C. In the latter case, a weak first-order or secondorder phase transition was seen as the surface density was increased. This is consistent with our results where a rapid decrease in tilt angle at higher surface densities was observed, with OPLS presenting the greatest decrease. However, it should be noted that our simulations were run in the NVT ensemble and so have a fixed cell geometry, and therefore are not able to accommodate a change in unit cell symmetry as a result of the phase transition. Another important measure of the degree of order within the monolayer is given by the carboncarbon radial distribution function (RDF). The RDF is a measure of the probability of finding one group of atoms between a distance r and r þ dr from another group of atoms, in this case the aliphatic carbon. As the surface density increases, the carbon chains form a more ordered state which is reflected in the carboncarbon intermolecular RDF shown in Figure 2c. RDFs generated by all three force fields share the same general shape with an initial broad peak around 45 Å and a broader peak around 9 Å. However, the all-atom force fields show greater detail around ∼6.8 Å, particularly OPLS. The differences between the united and all-atom based force fields are most evident around ∼4.25 and ∼5.3 Å. The distance corresponding to the largest peak (∼5.3 Å) is significantly larger than the nonbonded equilibrium separation between alkyl-type carbons (∼3.86 Å). It demonstrates that the intermolecular interactions between the alkyl chains are dominated by next nearest neighbors (NNN), in agreement with experimental studies of tilted phases.32 The differences between the united-atom and all-atom force fields may be associated with the presence of the explicitly represented hydrogen atoms in the all-atom force fields. In the united atom force field, the potential due to the H atoms is incorporated into the potential of carbon and therefore is uniform in all directions from the carbon. In Figure 3a, the average number of hydrogen bonds per time frame is shown as a function of surface density. We considered average number of hydrogen bonds between headgroups, displayed on the left axis, as well as hydrogen bonds between 3966
dx.doi.org/10.1021/jp1116867 |J. Phys. Chem. B 2011, 115, 3964–3971
The Journal of Physical Chemistry B
ARTICLE
Figure 3. Octadecanol monolayer headgroup properties at 298 K. (a) Average hydrogen bond number per time frame as a function of surface density with standard deviations included, and the radial distribution function (RDF) for octadecanol monolayer at 20.00 Å2 molecule1 for (b) headgroup oxygenheadgroup oxygen and (c) headgroup oxygenwater oxygen.
headgroups and water, which are displayed on the right axis. Hydrogen bonds were determined based on the cutoffs for the hydrogen donoracceptor angle (e30°) and the donoracceptor distance (e3.5 Å).33 For headgroupheadgroup bonding, the 53a6 and OPLS force fields present a similar trend at lower densities but begin to diverge at higher surface densities, while the converse occurs for headgroupwater hydrogen bonding. Overall, the 53a6 systems exhibit greater hydrogen bonding than the all-atom force fields. At low density, the octadecanol monolayer headgroups are fully accessible to water, enabling significant hydrogen bonding. However, as the surface density is increased, water is excluded from between the headgroups and the average number of hydrogen bonds with water decreases. Consequently, as the system is compressed, the headgroups on different chains come into closer proximity and the average number of H-bonds between headgroups increases. The hydration of the monolayer head groups can also be assessed by the radial distribution function between water oxygen and monolayer headgroup oxygen atoms (Figure 3c). This can be seen in Figure 3a where the hydrogen bonding between these molecules increases with decreasing surface density. The converse occurs when measuring the affinity between headgroups (Figure 3b): these interactions increase with increasing surface density, as illustrated in Figure 3a. Focusing on data obtained at 20 Å2 molecule1, we find that the all-atom force fields agree well with each other in terms of relative peak intensity for both RDFs, while the united-atom force field presents a distinctly higher peak for the headgroup oxygen pairs. However, there are slight variations in the location of the first peak, OPLS (∼2.8 Å), COMPASS (∼2.65 Å), and 53a6 (∼2.60 Å) for the headgroupheadgroup RDF. This is most likely a reflection of the particular parametrization of the interaction term between oxygen atoms in these force fields. This shift is not present in the carboncarbon RDFs, suggesting similar parametrization of the force field for aliphatic carbon. Overall, the RDFs for headgroup oxygenwater oxygen show much less variation, with 53a6 closely agreeing with the all-atom force fields. Analysis from torsions within the alkyl chains provides an additional measure of order and packing. The two all-atom based force fields present a similar headgroup conformation distribution (O1C1C2C3) shown in Figure 4a for molecular density 20 Å2 molecule1. In particular, we observe a major peak at 180° reflecting trans conformation and two smaller gauche conformation peaks at ∼60° and ∼300°. The 53a6 force field presents a distinctly different distribution, where less order is reflected by the broader distribution which slightly favors a trans conformation. This difference in behavior is a direct result
Figure 4. (a) O1C1C2C3 dihedral conformation distribution at 20.00 Å2 molecule1, and (b) the average length of the carbon chains as a function of surface density for OPLS and 53a6 force fields at 298 K with standard deviations included.
of the periodic potential used to describe this dihedral, whereas a RyckartBellemans type potential is used for OPLS which has a stronger potential corresponding to the trans conformation. The 53a6 force field also has a lower force constant which leads to less restrictive conformational orientations which are more subject to change compared to the OPLS force field, thus resulting in a broader dihedral distribution. As a result of the extra freedom the headgroup can more easily adopt conformations that enable it to form hydrogen bonds with water. Indeed the same reasoning could be applied to the interchain oxygen RDF which has a higher peak for the 53a6 system. Figure 4b displays the average lengths for the alkyl chains over the range of surface densities investigated. At the lower surface densities, the monolayer molecules are generally not fully extended due to the presence of gauche dihedrals in the chains. As the surface density increases, there is an associated decrease in the aliphatic gauche conformations as the chains are straightened. This trend is consistent with Gericke et al.34 who found the antisymmetric stretching vibration of the similar 1-hexadecanol to change only “slightly” upon compression (3219.5 Å2 molecule1). The actual differences between the force fields are insignificant over the length of 18 aliphatic carbons. Effect of Temperature. Studies of monolayers in the literature will often look at how a particular property changes with temperature,35 which can be used to explore phase behavior of a monolayer. Three temperatures were chosen that would represent the limiting cases of temperature effects on water (freezing, ambient, and boiling conditions) and span a large enough range of values such that results were appreciably different and trends, if any, would be evident. In this section we compare the performance of the three force fields for describing monolayers at the chosen temperatures 3967
dx.doi.org/10.1021/jp1116867 |J. Phys. Chem. B 2011, 115, 3964–3971
The Journal of Physical Chemistry B
ARTICLE
Figure 5. Octadecanol monolayer average tilt angle as a function surface area for the three force fields (a) OPLS, (b) 53a6, and (c) COMPASS. For each force field, the monolayer system at the three different temperatures 273, 298, and 368 K is shown. Standard deviations are included. (0° angle corresponds to monolayer perpendicular to the interface).
Figure 6. Intramolecular carboncarbon radial distribution function of octadecanol monolayer at 19.01 Å2 molecule1 for the three force fields (a) OPLS, (b) 53a6, and (c) COMPASS. For each force field, the RDFs at the three different temperatures 273, 298, and 368 K are shown.
Figure 7. RDFs of headgroup oxygenwater oxygen of octadecanol monolayer at 20.00 Å2 molecule1 for the three force fields (a) OPLS, (b) 53a6, and (c) COMPASS. For each force field, the monolayer system at the three different temperatures 273, 298, and 368 K is shown.
(273, 298, 368 K). Figure 5ac demonstrates how the tilt of the monolayer responds to different temperatures across the range of surface densities as simulated by the three force fields. Work by Lautz et al.31 established a first-order change for tilt angle in the higher density octadecanol monolayers; this change was seen to be stronger at 9 °C compared to 20 °C, where the change was described as weak first order or second order. Similar behavior was observed here for all three force fields; however, it is better reproduced by the two all-atom force fields. At 368 K all three force fields show a more shallow increase in tilt angle at the higher surface densities (>24 Å2 molecule1) compared with 273 and 298 K. For the lower surface densities at this temperature, the tilt angle trend for each force field was different. The deviation in the 368 K trend for OPLS (Figure 5a) was a result of the chains no longer forming an ordered monolayer at these temperatures and surface densities. The carbon chains were not straightened and a valid axis of tilt was
not present. This phenomenon also occurred in the 53a6 system at 45 Å2 molecule1. The low-density COMPASS systems did not exhibit this behavior; instead, the trend continued to increase.18 The response of the aliphatic carbons to changes in temperature can be analyzed with the intermolecular carboncarbon RDF as seen in Figure 6 with the surface density held constant at 19.01 Å2 molecule1. It is clear from Figure 6 that the three force fields show negligible differences in the CC RDFs as a function of temperature. This could be a consequence of the fixed geometry of the cell in the NVT ensemble which is limiting the lateral chain expansion and therefore maintains the average distance between the molecules. Figure 7 displays the RDFs between monolayer headgroup oxygens and water oxygens at varying temperatures. All three force fields yield a similar trend—an increase in temperature results in a decrease in the intensity of the H-bond peak at ∼2.7 Å. 3968
dx.doi.org/10.1021/jp1116867 |J. Phys. Chem. B 2011, 115, 3964–3971
The Journal of Physical Chemistry B
ARTICLE
Figure 8. Octadecanol O1C1C2C3 dihedral angle distribution at 20 Å2 molecule1 for (a) OPLS, (b) 53a6, and (c) COMPASS. For each force field, the monolayer system at the three different temperatures 273, 298, and 368 K is shown.
Figure 9. Octadecanol monolayer size comparison for (a) carboncarbon RDF at 20 Å2 molecule1 and 298 K, (b) headgroupheadroup oxygen RDF for OPLS and COMPASS force fields for 20 Å2 molecule1 at 298 K, and (c) dihedral distribution for OPLS and COMPASS force fields for 20 Å2 molecule1 at 298 K.
The temperature effects on the conformational behavior of the headgroup at constant surface density (20 Å2/molecule) can be seen in the O1C1C2C3 dihedral distribution shown in Figure 8. At 273 and 298 K there is very little change in the predominantly trans-oriented monolayer for the OPLS and COMPASS system. At 368 K there is a small reduction in the trans conformation which can be attributed to the higher kinetic energy of the system enabling it to overcome the barrier for rotation to the gauche conformation. This was not the case for the 53a6 force field where the dihedral distribution shows no discernible difference across the three temperatures. Effect of System Size. The systems discussed so far consisting of 20 octadecanol molecules have generally reproduced the behavior and properties observed experimentally. However, the size of these systems may be insufficient to describe variations in the phase space across the interface and to capture localized effects. Therefore, systems were considered in which the number of octadecanol molecules was increased to 80 which equates to doubling of the x- and y-dimensions of the cell. By increasing the scale of the system and comparing it with the above results, the differences arising from modeling different system sizes can be determined and the necessary system size identified for properties modeled. The subsequent analysis of system size will focus on the two allatom force fields OPLS and COMPASS. These two force fields have shown a stronger consistency with each other and experiment compared with the united-atom based 53a6 force field. The COMPASS system was run at a single surface density value (20 Å2 molecule1) with 80 molecules, due to the prohibitive computational cost of running these systems at this scale. The carbon RDFs in Figure 9a show there is no significant difference in the behavior of the carbon chains for the different size systems, particularly so for the COMPASS force field systems.
The finer detail in the smaller OPLS system was smoothed out in the larger system to more closely resemble the COMPASS distribution. The radial distribution functions for the oxygen groups did not vastly change with an increasing system size. Figure 9b demonstrates that the RDF obtained with the COMPASS force field is indistinguishable over the two sizes while the OPLS force field shows an increase in interchain oxygen affinity in the large system, with a corresponding decrease in the magnitude of the broader next-nearest-neighbor peak. Figure 9c displays the dihedral distribution of the headgroup region (O1C1C2C3) for the larger and smaller system with the two all-atom force fields. The similar distribution suggests the conformational configuration of the headgroup is not dependent on system size for the OPLS force field. Figure 10a displays the tilt angle as a function of surface density for the 20 and 80 molecule systems, respectively, obtained with OPLS. The trends in the two systems agree well; however, a notable difference is observed at the surface density of 45 Å2 molecule1. The 80-molecule system exhibits a significant reduction in tilt angle (∼13°) compared to its 20-molecule counterpart. At this particular surface density, the monolayer was initially in a tilted but wellordered configuration (Figure 10b) with all molecules tilted at approximately the same angle and direction. As the system evolved, the molecules began to aggregate, leaving a region of uncovered water (Figure 10c) in the cell. This process occurred within ∼1 ns and remained stable for the full duration of the simulation. To further confirm the validity of our results, we repeated the simulations with different random seeds to generate varying initial velocities. All simulations resulted in the formation of stable aggregates, irrespective of the initial velocities generated. 3969
dx.doi.org/10.1021/jp1116867 |J. Phys. Chem. B 2011, 115, 3964–3971
The Journal of Physical Chemistry B
ARTICLE
Figure 10. (a) Tilt angle as a function of system size for OPLS force field at 298 K with standard deviations included, and image of 45.00 Å2 molecule1 system’s initial configuration (b) and final configuration (c).
Figure 11. Top view of octadecanol molecules with backbone orientation overlaid, at 19.01 Å2 molecule1 at 298 K for OPLS system with SPC water.
We began the characterization of the aggregate by first calculating its surface density, which was found to be ∼31 Å2 molecule1. The monolayer tilt angle was 33.4 ( 2.9° which was less than the tilt angle of the equivalent surface density in the smaller system (∼46°) as determined from Figure 10a. The number of hydrogen bonds per frame as defined previously was 2.6 ( 0.8 which was again higher than the corresponding number obtained for the equivalent surface density of the smaller system. The OPLS system with SPC water did not exhibit this phenomenon in the large system; the trend continued to follow that of the smaller system. Indeed, this was the only discernible difference in the results presented thus far. We also observed a particular orientational ordering of the molecules which can be seen in Figure 11a, and which occurred solely in the OPLS/SPC system. The molecules aligned themselves into columns with a particular backbone orientation common to all molecules in that column. The behavior is repeated in the two adjacent columns, though the backbone orientation angle was different. This behavior occurred in the 19.01 system at 298 K, and in the 273 K systems at 19.01 and 20.00 and 24.20 Å2 molecule1, where the molecules aligned approximately parallel.
’ CONCLUSIONS We have compared the three force fields, OPLS, 53a6, and COMPASS, for molecular dynamics simulations of 1-octadecanol monolayers at an aqueous interface. Monolayer properties as a function of surface density were compared at three different temperatures (273, 298, and 368 K) and across two different system sizes (20 and 80 chain molecules per unit cell).
Properties investigated include average radial tilt angle, dihedral distribution, average lengths, and intermolecular interactions including hydrogen bond and nonspecific carbon chain interactions. Where possible, we have compared the results to experiment to assess the efficacy of the force fields for these types of systems. Hydrogen bonding was characterized by headgroup oxygenheadgroup oxygen RDFs, and headgroup oxygenwater oxygens RDFs with the 53a6 systems showing a much more significant hydrogen-bonding peak at ∼2.7 Å compared to the two all-atom-based force fields. A possible explanation for this can be found in the 53a6 dihedral distribution which shows little preference to any one discrete conformation. The carboncarbon RDFs provided insight into the difference in chain ordering predicted by the united-atom-based force field 53a6 and the all-atom OPLS and COMPASS force fields. The all-atom force fields predict greater ordering of the chains compared to the united-atom force field. The tilt angle was measured as a function of surface density, and all three force fields yielded similar trends: an initially slowly decreasing tilt angle with increasing surface density, then a rapid decrease in tilt angle at higher surface density values. This is in qualitative agreement with experiment which shows that octadecanol monolayers undergo a first-order phase change at high surface densities. This phase change is accompanied by a sudden decrease in tilt angle, as we have observed. The OPLS force field agreed best with this observation. The three monolayer systems were examined at 273, 298, and 368 K to determine how they respond to changes in temperature. The carboncarbon RDFs remained essentially unchanged over the three temperatures due to spatial limitations of a fixed cell geometry. However, appreciable changes in tilt angle and headgroup conformity were observed. The two all-atom force fields (OPLS and COMPASS) displayed good agreement with experiment. For a larger system model, the OPLS results deviated only slightly from those obtained with the smaller models while no obvious changes were evident in the COMPASS force field simulation. The tilt angle in the larger OPLS system was shown to be a few degrees smaller compared to the smaller system. This may be important when studying phase behavior of monolayers, where the tilted/untilted transition point indicates a 2-dimensional phase change.31,36,37 The interchain headgroup oxygen RDF in the larger OPLS system suggested a stronger interaction between monolayer headgroups. Whether or not there is a relationship between tilt angle and headgroupheadgroup/water affinity is yet to be determined, though it is not unreasonable that an untilted monolayer with a hydroxyl headgroup in trans conformation will have a greater affinity to an water surface than a tilted one. There were no changes observed in the larger COMPASS system. 3970
dx.doi.org/10.1021/jp1116867 |J. Phys. Chem. B 2011, 115, 3964–3971
The Journal of Physical Chemistry B We have found that large systems (g80 molecules) are required to observe monolayer aggregation at low surface densities, though it should be noted there are small changes in interchain headgroup behavior and monolayer tilt angle compared to the smaller systems (20 molecules) with equivalent surface densities. The all-atom-based force fields can be recommended for MD monolayer simulations that are aiming to investigate the effects of temperature. There was little difference in results obtained with the all-atom-based OPLS and COMPASS force fields, though OPLS better represented the monolayer tilt angle behavior which may have implications for MD simulations involving phase behavior. Also, the greater computational efficiency of OPLS and 53a6 force fields compared to COMPASS, as presented in the context of organic monolayer/ water systems, should not be discounted as an important factor when selecting force fields.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected]. Present Addresses †
School of Chemical and Mathematical Sciences, Murdoch University, 90 South St., Perth, Western Australia, 6150 Australia.
’ ACKNOWLEDGMENT We thank the National Computational Infrastructure (NCI) and Victorian Partnership for Advanced Computing (VPAC) for allocation of computational resources. ’ REFERENCES
ARTICLE
(19) Prathab, B.; Subramanian, V.; Aminabhavi, T. Polymer 2007, 48, 409–416. (20) Todorova, N.; Legge, F.; Treutlein, H.; Yarovsky, I. J. Phys. Chem. B 2008, 112, 11137–11146. (21) Rog., T.; Vattulainen., I.; Karttunen., M. Cell. Mol. Biol. Lett. 2005, 10, 625–630. (22) Damm, W.; Frontera, A.; Tirado-Rives, J.; Jorgensen, W. J. Comput. Chem. 1997, 18, 1955–1970. (23) Friesner, R.; Tirado-Rives, J.; Jorgensen, W.; Kaminski, G. A. J. Phys. Chem. B 2001, 105, 6474–6487. (24) Oostenbrink, C.; Villa, A.; Mark, A.; Van Gunsteren, W. J. Comput. Chem. 2004, 25, 1656–1676. (25) Ferguson, D. J. Comput. Chem. 1995, 16, 501–511. (26) Berendsen, H.; Postma, J.; van Gunsteren, W.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981; p 331. (27) Jorgensen, W.; Chandrasekhar, J.; Madura, J.; Impey, R.; Klein, M. J. Chem. Phys. 1983, 79, 926. (28) Essmann, U.; Perera, L.; Berkowitz, M.; Darden, T.; Lee, H.; Pedersen, L. J. Chem. Phys. 1995, 103, 8577–8593. (29) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (30) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput 2008, 4, 435–447. (31) Lautz, C.; Fischer, T. M.; Weygand, M.; Losche, M.; Howes, P. B.; Kjaer, K. J. Chem. Phys. 1998, 108, 4640–4646. (32) Steitz, R.; Peng, J.; Peterson, I.; Gentle, I.; Kenn, R.; Goldmann, M.; Barnes, G. Langmuir 1998, 14, 7245–7249. (33) Luzar, A. J. Chem. Phys. 2000, 113, 10663. (34) Gericke, A.; Simon-Kutscher, J.; Huehnerfuss, H. Langmuir 1993, 9, 3115–3121. (35) Kaganer, V.; Brezesinski, G.; Mohwald, H.; Howes, P.; Kjaer, K. Phys. Rev. E 1999, 59, 2141–2152. (36) Lautz, C.; Fischer, T. J. Phys. Chem. B 1997, 101, 8790–8793. (37) Riviere-Cantin, S.; Henon, S.; Meunier, J. Phys. Rev. E 1996, 54, 1683–1686.
(1) Torchilin, V. Nat. Rev. Drug Discovery 2005, 4, 145–160. (2) Hirsch, T.; Kettenberger, H.; Wolfbeis, O.; Mirsky, V. Chem. Commun. 2003, 432–433. (3) Feng, Y.; Teo, W.-K.; Siow, K.-S.; Gao, Z.; Tan, K.-L.; Hsieh, A.K. J. Electrochem. Soc. 1997, 144, 55–64. (4) Barnes, G. Adv. Colloid Interface Sci. 1986, 25, 89–200. (5) Chandross, M.; Grest, G.; Stevens, M. Langmuir 2002, 18, 8392–8399. (6) Tutein, A.; Stuart, S.; Harrison, J. Langmuir 2000, 16, 291–296. (7) Himpsel, F.; Ortega, J.; Mankey, G.; Willis, R. Adv. Phys. 1998, 47, 511–597. (8) Henelius, P.; Fr€obrich, P.; Kuntz, P.; Timm, C.; Jensen, P. Phys. Rev. B 2002, 66, 94407. (9) Collet, J.; Vuillaume, D. Appl. Phys. Lett. 1998, 73, 2681. (10) Baoukina, S.; Monticelli, L.; Amrein, M.; Tieleman, D. Biophys. J. 2007, 93, 3775–3782. (11) de Vries, A.; Mark, A.; Marrink, S. J. Am. Chem. Soc. 2004, 126, 4488–4489. (12) Lindahl, E.; Edholm, O. Biophys. J. 2000, 79, 426–433. (13) Jang, S.; Jang, Y.; Kim, Y.; Goddard, W., III; Flood, A.; Laursen, B.; Tseng, H.; Stoddart, J.; Jeppesen, J.; Choi, J. J. Am. Chem. Soc. 2005, 127, 1563–1575. (14) Hawthorne, M.; Zink, J.; Skelton, J.; Bayer, M.; Liu, C.; Livshits, E.; Baer, R.; Neuhauser, D. Science 2004, 303, 1849. (15) Marrink, S.; Mark, A. J. Am. Chem. Soc. 2003, 125, 11144–11145. (16) Srivastava, P.; Chapman, W.; Laibinis, P. Langmuir 2009, 25, 2689–2695. (17) Sun, H. J. Phys. Chem. B 1998, 102, 7338–7364. (18) Henry, D. J.; Dewan, V. I.; Prime, E. L.; Qiao, G. G.; Solomon, D. H.; Yarovsky, I. J. Phys. Chem. B 2010, 114, 3869–3878. 3971
dx.doi.org/10.1021/jp1116867 |J. Phys. Chem. B 2011, 115, 3964–3971
