928
Langmuir 2002, 18, 928-931
Molecular Dynamics Simulation of (Octadecylamino)dihydroxysalicylaldehyde at Air/Water Interface A. Dhathathreyan*,† and S. J. Collins‡ Chemical Laboratory, Central Leather Research Institute, Adyar, Chennai 600 020, India, and Indira Gandhi Centre for Atomic Energy Research, Kalpakkam 603 102, India Received July 12, 2001. In Final Form: October 24, 2001 Surface pressure-molecular area (π-A) isotherms for monolayers of a new amphiphilic Schiff base, (octadecylamino)dihydroxysalicylaldehyde (ODADS), at an air/water interface were simulated by the molecular dynamics (MD) method for the first time. Potential energy, torsional angles, and hydration of the polar headgroups under periodic boundary conditions (PBC) have been used for this study. Hydration was done by a stacked array of water molecules relaxing under these boundary conditions. The simulations were carried out in a hexagonal lattice with the spatial PBC constraints in order to avoid clustering of the water molecules. The resulting isotherm with the transitions seen in the π-A isotherm are explained from the potential energy plots for various surface densities.
Introduction In recent times molecular dynamics simulation (MDS) has been widely attempted with Langmuir monolayers that are amphiphiles spread at an air/water interface.1-7 The structure and molecular orientation of the monolayers can be predicted theoretically by this method. These results agree to a great deal of accuracy with a number of experimental observations.8,9 Simulation of ordered domains often observed in monolayer organization and their dependence on length of the molecules and layer density have been carried out on simple long-chain carboxylic acids.10 One of the most important influences on the packing of the monolayer is certainly the polar headgroup, as a result of steric hindrance as well as its specific interactions with the monolayer environment. Freyer et al.11 have shown that aromatic hydrocarbons with hydrophilic substituents form stable layers provided there is a balance between hydrophobic and hydrophilic ends of the molecule. Langmuir-Blodgett (LB) films of such functionalized molecules with cross-sectional mismatches between head and tail are well characterized by synchrotron, electron microscopy (EM), X-ray diffraction (XRD), and atomic force microscopy * Author to whom correspondence should be addressed: fax +91(44)-4911589; e-mail
[email protected]. † Chemical Lab., CLRI. ‡ IGCAR. (1) Harris, J.; Rice, S. A. J. Chem. Phys. 1988, 89, 5898. (2) Karaborni, S.; Toxvaerd, S. J. Chem. Phys. 1992, 96, 5505. (3) Karaborni, S.; Toxvaerd, S. J. Chem. Phys. 1992, 97, 5876. (4) Karaborni, S.; Toxvaerd, S.; Olsen, O. H. J. Phys. Chem. 1992, 96, 4965. (5) Alper, H. E.; Bassolino, D.; Stouch, T. R. J. Chem. Phys. 1993, 98, 1. (6) Van Bauern, A. R.; de Vlieg, J.; Berendson, H. J. C. Langmuir 1995, 11, 2957. (7) Okamura, E.; Fukushima, N.; Hayashi, S. Langmuir 1999, 15, 3589. (8) Collins, S. J.; Dhathathreyan, A.; Ramasami, T. J. Colloid Interface Sci. 1998, 203, 249. (9) Gallet, X.; Deleu, M.; Razafindralambo, H.; Jacques, P.; Thonart, P.; Paquot, M.; Brasseur, R. Langmuir 1999, 15, 2409. (10) Rovillard, S.; Perez, E.; Ionov, R.; Voue, M.; De Coninck, J. Langmuir 1999, 15, 2749. (11) Freyer, J. H.; McConnell, C. H.; Grant, G. H.; Hann, R. A.; Gupta, S. K.; Eyres, B. L. Philos. Mag. B 1991, 63, 1193.
Chart 1 Structure of ODADS
(AFM).12-15 Interactions of metal ions with Schiff basecontaining amphiphiles have been extensively reported in the literature.16-18 The present work describes the MDS study carried out on a novel Schiff base amphiphile, (octadecylamino)dihydroxysalicylaldehyde (ODADS) (Chart 1), synthesized in our laboratory to study phase transitions in amphiphiles having a cross-sectional mismatch between the head and tail groups.19 We deal with MD simulation of a Langmuir monolayer of ODADS, and the experimentally observed surface pressure-molecular area (π-A) isotherm of ODADS monolayers at an air/water interface were simulated by the molecular dynamics (MD) method for the first time. Potential energy, torsional angles, and hydration of the polar headgroups under periodic boundary conditions (PBC) have been used for this study. The intermolecular potential functions include short-range van der Waals potentials, long-range Coulombic potentials, and many body potentials. Short-range interaction potentials include 12-6 Lennard-Jones, n-m Buckingham and BornHuggin-Meyer H-bond potentials, and three body systems are represented by truncated harmonic, screened harmonic, and Dreiding potentials. A system of charged point ions mutually interacting via a Coulomb potential using the Ewald sum for calculating the electrostatic interaction is used. (12) Vaknin, D.; Kjaer, K.; Als-Nielsen, J.; Losche, J. Biophys. J. 1991, 59, 1325. (13) Als-Nielsen, J.; Mohwald, H. In Handbook of synchrotron radiation; Ebashi, S., Rubenstien, E., Koch, M. I., Eds.; Vol. 4; Elsevier: Amsterdam, 1981. (14) Meyer, E.; Howald, I.; Overney, R. M. Nature 1991, 349, 398. (15) Bohm, C.; Mohwald, H.; Leiserowitz, L.; Als-Nielsen, J.; Kjaer, K. Biophys. J. 1993, 64, 553. (16) Hemakanthi, G.; Dhathathreyan, A. Langmuir 1999, 15, 3317. (17) Shyamalasundari, S.; Dhathathreyan, A.; Kanthimathi, M.; Nair, B. U. Langmuir 1997, 13, 3422. (18) Hemakanthi, G.; Dhathathreyan, A.; Moebius, D. Colloids Surf. (in press). (19) Collins, S. J.; Dhathathreyan, A.; Ramasami, T.; Mohwald, H. Thin Solid Films 2000, 358, 229.
10.1021/la011073e CCC: $22.00 © 2002 American Chemical Society Published on Web 01/09/2002
MD Simulation of ODADS at Air/Water Interface
Figure 1. (a, top panels) Simulated plots of dynamic total energy, dynamic potential energy, dynamic kinetic energy, and dynamic pressure of ODADS monolayer. (b, bottom panel) π-A isotherm of ODADS obtained experimentally (subphase water, T ) 22 °C).
Hydration was done by a stacked array of water molecules relaxing under these boundary conditions. The simulations were carried out in a hexagonal lattice with the spatial PBC constraints in order to avoid clustering of the water molecules. It is known that with hexagonal symmetry the monolayer is composed of unit cells containing a single hydrocarbon chain. The resulting isotherm with the transitions seen is explained from the potential energy plots for various surface densities. Electron diffraction patterns of the films transferred onto transmission electron microscopy (TEM) grids were studied by highresolution imaging and showed a clear agreement with the packing behavior predicted by the MDS.
Langmuir, Vol. 18, No. 3, 2002 929
Figure 2. (a, top panel) Plot of average bond angles as a function of bond number in ODADS. (b, bottom panels) Plot of bond variation for various packing densities.
with X-ray reflectivity measurements.22 The potential energy is then partitioned explicitly for a single molecule orientation potential, dipolar couplings, rotationalrotational, translational-rotational, and translationaltranslational interactions, and the simulations were carried out using the Insight II DISCOVER program.23 In this the simulation of the molecule can be done after the molecule is built and charges and potentials are assigned to each of the atoms or groups. The potential energy of the molecules consists of diagonal terms, off-diagonal or cross terms, and nonbonded interactions, and the total potential energy of the molecule is given by
∑D {1 - e (b - b ) } + / ∑H (θ - θ ) + / ∑H {1 + s cos (nφ)] + / ∑H χ + ∑∑F (b - b )(b′ - b′ ) + ∑∑F (θ - θ )(θ′ - θ′ ) + ∑∑F (b - b )(θ - θ ) + ∑F cos (φ) (θ - θ ) × (θ′ - θ′ ) + ∑∑F χχ′ + ∑ {(r′/r) - 2(r′/r) } + ∑q q /r + ∑(C /r - D /r ) -R
E)
There have been a number of attempts to use MDS on monolayers, using the interactions of headgroups and tails between themselves and with water by working out the interaction potentials of Lennard-Jones type with strong or weak attractive or repulsive terms as the case may be.20 In this study, we begin by defining the Langmuir monolayer as a system of close-packed rigid rod-type molecules of hexagonal symmetry. Such models of grafted rods where molecules are approximated by rigid rodlike particles grafted onto planar surface were first used by Kaganer et al. in their simulation work on amphiphiles21 and also in their studies relating the hexatic mesophase (20) Kraer, M.; Kremes, K.; Binder, K. J. Chem. Phys. 1990, 92, 6195. (21) Kaganer, V. M.; Osipow, M. A.; Peterson, I. R. J. Chem. Phys. 1993, 98, 3512.
1
0
2
2
1
1
2
2
φ
0
0
θ
bθ
0
0
θθ′
b′
b
χ
χ
φ
bb′
b
0
θ
θ
2
Model
2
b
b
0
θ′
0
0
φθθ′0
θ
φ
12
0
6
χχ′
χ
χ′
i j
ij
ij
ij
12 ij
10
ij
ij
ij
The first four terms in the expression correspond to the energy of deformation of bond length, valence angles, torsional angles and out-of plane interactions and are referred to as the diagonal terms of the valence field. Db (22) Kaganer, V. M.; Loginov, E. B. Phys. Rev. 1995, E51, 2237. (23) DISCOVER 2.5 manual, Biosym Inc., 1990.
930
Langmuir, Vol. 18, No. 3, 2002
Dhathathreyan and Collins
Figure 3. Configuration of ODADS at (a, left panels) low packing density and (b, right panels) high packing density.
denotes the dissociation energy, R the anharmonicity constant, b0 the equilibrium bond length, and Hθ and Hφ the force constants according to the earlier reported model.24 Conformation of the molecule is obtained by minimizing the starting structure and the molecular dynamic technique by use of F ) -∇(φ) and numerical integrations of a ) F/m simulate the global minima. DISCOVER adjusts the velocities of the atoms to maintain a constant average kinetic energy by use of a leapfrog integration algorithm.25 To minimize the molecular system, the conjugate gradient and Newton-Raphson methods were used. Initially the strained molecular system was minimized by the conjugate gradient method for 1000 iterations. This procedure is slow but steadily approaches a minimum. The next 1000 iterations were done by the Newton-Raphson method, which minimizes at a faster rate. It was found that the rms derivative in the energy at the final iteration was less than 0.001 kcal/mol. Initial dynamics calculations for 1000 steps were done to equilibrate the molecular system. The pressure and dynamic energy parameters were calculated from another 10 000 steps of dynamics. Before the MD is run, the system should relax the highly strained structures or overlapping atoms, which may take additional time to dissipate the unrealistic thermal energy or even break the molecules apart. Dynamics is initialized by assigning random velocities consistent with MaxwellBoltzmann distribution for the set temperature. Trajectories of the dynamics are produced at a time step around 1 fs. In the simulation of ODADS, all atoms and their full set of interactions are taken into account for the molecular system of two-dimensional packing of ODADS on a packed (24) Taut, C.; Pertsin, A. J.; Grunze, M. Langmuir 1996, 12, 3481. (25) Toxvaerd, S. Mol. Phys. 1991, 72, 159.
layer of water molecules. The water layer is simulated with dense packing of water molecules in a hexagonal lattice with spatial PBC constraints in order to avoid clustering of water molecules. The water layer is energyminimized and used for the subphase of the Langmuir monolayer. The simulation box consists of nine ODADS molecules arranged in a hexagonal lattice under PBC. A water layer of the same size is associated with the headgroup region in the simulation box. The whole system is minimized for the full potential energy by DISCOVER.25 The minimum energy is found to occur around 24 Å2 of packing area. The bond angle energy and coulomb energy due to partial charges also go to minima at this area. Figure 1 shows the energy parameters of dynamics. The dynamic total energy starts decreasing around 30 Å2/ molecule and reaches a minimum between 22 and 24 Å2/ molecule, and the kinetic energy is approximately constant. The corresponding pressure reaches a low value at 24 Å2/molecule and this is seen also from the experimental isotherm (Figure 1b). To analyze the conformation of the molecules, bond angles between the various bonds of the long alkyl chains are calculated. The orientations of the 18 bonds in the alkyl chain are studied for the nine molecules in the simulation box. Figure 2a shows the average bond angles of the nine molecules plotted as a function of the various bond numbers (starting from the carbon closest to the aromatic ring). Analyzing the plots for all the nine simulated systems at various packing densities, Figure 2b shows the bond variation at various packing densities. It is found that maximum bending of the bonds progresses from the tail end toward the aromatic ring in the polar plane when the packing density increases. At the highly packed state maximum bending takes place at the bond between the chain and the aromatic ring. Such an effect
MD Simulation of ODADS at Air/Water Interface
has been seen by Okamura et al.7 in their simulation studies of vibrational spectra of stearic acid monolayers at the air/water interface. The molecular configurations of ODADS monolayer at high and low packing densities in Figure 3 show the packing looking through x, y, and z axes with the full stretched view. At the highest packing density, the hexagonal lattice of the monolayer is maintained at the end of the simulation. The experimental results with the TEM micrograph of a LB film of ODADS transferred onto a carbon-coated copper grid (one layer, transferred at π ) 30 mN/m) showed a clear hexagonal packing for the ODADS molecules at the air/water interface.19 In conclusion, MDS has been carried out to study the conformation and tilt of ODADS molecules with respect to the water surface. The cross-sectional mismatch between the headgroup and the long chain in this amphiphile is analyzed with molecular parameters and bond angles of the simulated molecules in the monolayer
Langmuir, Vol. 18, No. 3, 2002 931
at various packing densities. The molecules undergo a transition to lower energy state at 24 Å2 and with a change in the orientation of the aromatic headgroup about the subphase plane. Molecular tilt is greater at the carbon of the alkyl chain nearest to the headgroup at the highest packing density. The simulated pressure-area curve is in agreement with the experimentally observed π-A isotherm. The possible limitation of the present simulation is the smaller dimension size. Much larger time and length scales need to be used to solve several problems, especially with respect to long-range interactions. Acknowledgment. We thank the Department of Science and Technology, Government of India, for a project grant under which this work was carried out. LA011073E
