Print version ISSN 0103-9733
Braz. J. Phys. vol.34 no.1 São Paulo Mar. 2004
Molecular Dynamics simulation of the sodium octanoate micelle in aqueous solution: comparison of force field parameters and molecular topology effects on the micellar structure
André Farias de Moura; Luiz Carlos Gomide Freitas
Departamento de Química, Centro de Ciências Exatas e de Tecnologia, Universidade Federal de São Carlos, Rodovia Washington Luiz, km 235, C.P. 676, 13565-905, São Carlos, SP, Brazil
We have performed a series of 10 ns Molecular Dynamics simulations of the sodium octanoate micelle in aqueous solution in the constant NpT ensemble, at p = 1 bar and T = 300 K. Two molecular topologies were studied, one with all internal degrees of freedom and the other constraining bond stretching and angle bending degrees of freedom. Two Lennard-Jones parameters for sodium ions, namely the OPLS and ° Aqvist parameters, were used. The results show an artificial enhancement of stable sodium bridges between octanoate anions when the OPLS parameters for sodium are used. The ° Aqvist parameters give a micellar structure in good agreement with experimental and thermodynamical evidences. It is also observed that the aggregation of monomers is strongly dependent on the molecular topology. When the ° Aqvist parameters were employed, the model system without constraining geometry had one dissociated monomer after 10 ns, while the model system with bond length and bond angle constraining had five dissociated monomers after a 10 ns trajectory.
Molecular Dynamics simulation may be regarded as the method of choice to study complex liquids in atomic detail, yielding both structural and dynamical information of systems ranging from micelles to large biomolecules and polymers in solution. Model systems as large as 104-105 particles may be tackled with single CPU machines and it is possible to treat larger system if parallel processing is available.
As a general trend the interaction between particles is represented by parameterized classical potential energy surfaces (PES), whose parameters are developed to reproduce thermochemical and other experimental and theoretical data available. The set of parameters needed to reproduce the PES is called force field. Several force fields have been developed up to now and it is a common practice to optimize parameters to reproduce quantum chemical calculations of dimers in vacuum and thermochemical measurements in dilute solution. Nonetheless, complex liquids are generally concentrated solutions and it is not clear whether or not the usual force field parameters are adequate to simulate such systems.
Micellar systems have been studied by means of computer simulations in the last 30 years. There has been some concern about methodological issues, but there has not been a thorough evaluation of the interaction parameters and molecular topologies used so far.
Binary systems formed by sodium octanoate and water have been extensively studied by Molecular Dynamics simulations in the last 15 years. In these studies the most common model system used consisted of a single pre-assembled micelle comprising 15 octanoate anions [1-9], although some simulations evaluated the liquid crystalline phase at higher surfactant contents [10-11].
Jönsson et al. performed the first Molecular Dynamics simulation of the sodium octanoate micelle in aqueous solution . The interaction parameters were taken from several sources and the molecules were modeled including all atoms and all internal degrees of freedom. Two simulations were performed to compare the effect of an ad hoc screening introduced in the interaction between ionic pairs. A screening factor was achieved by means of cutting down the interactions due to atomic charges to half of their values. This procedure was justified by assuming that the SPC water model  is not able to account for the dielectric constant of pure water. Jönsson et al. reported a very interesting and detailed set of information regarding both structural and dynamical aspects of the micelle, but the molecular dynamics trajectories were too short and the precision was not good enough (the trajectory without charge screening lasted only 30 ps, while the run with screening lasted 72 ps). Another drawback of these simulations was the absence of any scheme to treat long-range interactions, which were truncated at 1.0 nm.
Klein and co-workers have performed a series of Molecular Dynamics simulations using different model systems for the aqueous solution of sodium octanoate micelle [2-5]. All their simulations have been carried out in the constant NVT ensemble, employing OPLS pseudoatoms parameters for octanoate anion  and sodium cation  with bond lengths constraints.
Shelley et al.  compared two similar water models, namely SPC  and SPC/E . The later has a slightly larger dipole moment, but no significant structural difference was reported. Two previous simulations were performed with SPC water model and minor differences have also been reported, but the model systems differed in size and composition. Therefore, it is difficult to make a detailed comparison of the results obtained in these calculations.
Shelley et al. also evaluated a polarizable model to represent the surfactant molecule . Charge polarization was achieved with two variable point charges with the same magnitude but opposite sign assigned to each heavy atom. These variable charges changed their magnitude in response to the local electrical field, while moving around the heavy atoms as rigid rotors governed by an extended Lagrangian scheme. Despite the fact that polarization was treated classically, it is an interesting work. The most striking feature reported regards the counter-ion distribution around the octanoate headgroups. The first intense peak previously reported [2,4] in the radial distribution function of sodium around the headgroup carbon atom was suppressed in the polarizable model simulation.
The polarizable model also increases the probability of finding water molecules near the aggregate center of mass, although a dry core region was preserved. The values obtained for the hydration numbers of hydrophobic carbon atoms were also larger, meaning that polarization increases the hydrocarbon-water affinity when compared with the nonpolarizable models [2-4].
Other Molecular Dynamics simulations have been also reported in the literature. Kuhn and Rehage [6-7] and Kuhn et al.  performed two simulations in the constant NpT ensemble, using the AMBER force field [16-18] without geometry constraints. Laaksonen and Rosenholm performed a Molecular Dynamics simulation in the constant NVT ensemble, using the CHARMM forcefield  with semi-empirical charges and the TIP3P water model , both without geometry constraints.
Altogether, these simulations of the sodium octanoate micelle in aqueous solution gave some insight into the micellar structure. The main feature was the existence of a small dry core, even though apolar tail were in contact with surrounding water molecules. Most studies reported the aggregate size in good agreement with experimental data available and found out a shape resembling a prolate ellipsoid. It is worth to note that the different studies do not agree about one very important issue: the distribution of counter-ions around the micellar aggregate. Further effort will be necessary to enlighten the structure of sodium cations near octanoate micelles.
Other aqueous micellar systems have also been studied by Molecular Dynamics simulations, such as octyl glucoside , dodecylphosphocholine [22-24], lysophosphatidylethanolamine , sodium pentadecafluorooctanoate , sodium dodecyl sulfate [27-29] and 2,3,6,7,10,11-hexakis(1,4,7-trioxaoctyl)triphenylene . Except for sodium dodecyl sulfate [27-29] all micellar systems studied are non-ionic, what makes computer simulations comparatively less time consuming, since long range interactions may be safely ignored when ionic species are absent.
In the present work, we report some structural results of Molecular Dynamics simulations on the sodium octanoate micelle in aqueous solution using different Lennard-Jones parameters for the sodium counter-ion: the OPLS  and the Åqvist parameters . The main purpose is to investigate structures containing sodium ions in solution, trying to reconcile theoretical and experimental results reported so far. We also report the differences arising when the molecular geometries are constrained, using these two parameters sets.
We have performed a set of Molecular Dynamics using the GROMACS software [32-33]. The simulations were performed on the constant NpT ensemble with the Berendsen temperature bath and pressure bath . The OPLS-AA parameters for carboxylate anions  were employed to describe the octanoate interaction, while the SPC model  was chosen for water molecules. The sodium cations were described by two parameters sets: the OPLS  and the Åqvist . Two different topologies were used: one with all internal degrees of freedom and the other with constrained bond lengths and angles. We shall refer to the former as the flexible model and to the latter as the rigid model. Octanoate geometry was constrained using LINCS  while water molecules were constrained with SETTLE .
To save computer time a cut off at 1.5 nm was used for the interaction potential and a neighbor list updated every tenth step. The particle mesh Ewald method [38-39] was used to treat long range electrostatic interactions beyond cut off distance.
In the flexible model simulations an integration time step of 1.0 fs, a temperature coupling of 0.1 ps and a pressure coupling of 1.0 ps were used. For the rigid model, a time step of 2.0 fs, a temperature coupling of 0.02 ps and a pressure coupling of 2.0 ps were used. For each model system a molecular dynamic trajectory with a length of 10 ns was generated.
The model system consisted of 15 octanoate anions, previously assembled into a spherical aggregate, surrounded by 15 sodium counter-ions and 2190 water molecules. This composition is around the critical micelle concentration, cmc . It is not a common practice to work near transition regions, but we have chosen to study the micellar system near its cmc to increase the number of water molecules, avoiding as much as possible undesirable effects due to the system size.
The simulation box was built as follows. The aggregate was assembled into a nearly spherical shape. The first monomer was aligned along the x-axis, the methyl carbon atom being at the origin. This monomer was rotated around the z-axis, forming an equatorial plane of six equally spaced octanoate anions. The first monomer was then aligned in the +/- 45o direction in the xz-plane and was rotated around the z-axis, thus forming two intermediate layers of octanoate anions, one with five molecules and the other with four. At last, the first monomer was tilted +/- 90o, adding two axial monomers to the aggregate. As a final preparation step, the monomers were displaced 0.5 nm away from the origin to avoid any strong repulsive contact between them.
The system density was chosen to be 1.0 g cm3. First, the aggregate was placed in the middle of the simulation box. Then, water molecules and counter-ions were randomly placed into the box, avoiding any repulsive contact. The resulting model system was made up of 6960 atoms.
Both steepest descent and conjugated gradient algorithms were employed to minimize the energy gradient of the system below 25 kJ mol1 nm1. After energy minimization, a series of simulated annealing runs were performed as follows: (i) the system was heated to T = 50 K; (ii) the temperature was linearly decreased to T = 0 K within a time interval of 1 ps. The procedure was repeated until no significant energy drift was observed at the final temperature.
The lowest energy structure was taken as the initial simulation box (Fig. 1). The flexible and the rigid systems were heated to the working temperature by performing three controlled short range Molecular Dynamics simulations as follows: (i) in the first simulation the position of the micelle was constrained and the constant NVT ensemble was used; (ii) in the second simulation, the position restraint was kept but the constant NpT ensemble was used; (iii) the third simulation was performed in the same conditions used in the production run, that is, without any position constraint.
Quantum chemical calculations of dimers in vacuum were carried out with Gaussian 98 . In contrast with OPLS force field, which is parameterized to reproduce Hartree-Fock potential energy surfaces in vacuum, we performed ab initio calculations using second order Møller-Plesset perturbation theory with 6-311+G** basis set, applying counter-poise correction to account for the basis set superposition error.
3 Results and discussion
In this work, we report a series of Molecular Dynamics Simulation of the sodium octanoate micelle in aqueous solution. We aim to investigate the effects of the molecular topology and the force field parameters on the micellar structure. Two molecular topologies were assigned to the octanoate anion, (i) a fully flexible model were all internal degrees of freedom are considered, and (ii) partially rigid model, were bond stretching and angle bending degrees of freedom are constrained to their equilibrium values. We shall refer to the former as the flexible model and to the later as the rigid model. Besides, two sets of force field parameters were assigned to the sodium cation, the OPLS parameters and the Åqvist parameters. Altogether, four Molecular Dynamics simulations were accomplished, and in each one averaging trajectories with 10 ns were generated.
The results of density and radius of gyration presented below are running averages of the original data. To remove the high frequency noise and to enhance the trends in the micelle behavior the averaged values were calculated using a 50 ps step between consecutive points in the trajectory.
4 OPLS parameters for sodium cation
Molecular Dynamics simulations have been carried out in the constant NpT ensemble. As expected from literature reports the pressure coupling scheme should lead the systems to their equilibrium densities in a short time period. Nonetheless, the systems densities were found out to converge slowly when the OPLS parameters for sodium ion were employed (Fig. 2). The relaxation was accomplished after ca. 1000 ps for the flexible model and after ca. 6000 ps for the rigid model. Such poor convergence rate indicates that slow structural changes are taking place in the model systems. This unexpected behavior was investigated and the nature of the slow structural changes will be discussed.
The micelle radius of gyration for the flexible model (Fig. 3) does not show any trend that might be correlated with the density relaxation shown in Fig. 2. On the other hand, the rigid model presented an increase in the aggregate radius of gyration after ca. 6000 ps (Fig. 3). Comparing with the data in Fig. 2, this increase in the aggregate radius of gyration is observed simultaneously with relaxation of system density.
It is also observed that the distribution of sodium counter-ions around the aggregates center of mass depends on the model topology: The rigid model presented an increase of the sodium density near the micellar center of mass as compared to the flexible model (Fig. 4).
The number of sodium atoms in the first coordination shell of the octanoate oxygen atoms increased in both simulations (Fig. 5). The increase in the coordination numbers is due to the counter-ion condensation on the micellar surface, which is characterized by the formation of sodium bridges between octanoate anions (Figs. 6 and 7). The coordination numbers suggest that the slow structural relaxation mentioned earlier can be attributed to a complex condensation phenomenon, possibly taking place in two steps: (i) the counter-ion reaches the micellar surface and is trapped by an octanoate headgroup; then, (ii) these ion pairs slowly rearrange themselves to form clusters resembling ionic lattices on the micellar surface (Figs. 6 and 7).
The behavior depicted from the average data reported above does not agree with the accepted view of micelle structure supported by theoretical and experimental evidences . In the theoretical and experimental data presented so far there is strong support to models in which the counter-ions are distributed around the micelle surface. It is also expected that the counter-ion concentration measured from the micelle center of mass decrease with the distance down to a non-zero bulk value. Nevertheless, due to size dependence, any model system that might be used in a Molecular Dynamics simulation would not be totally adequate, no matter the parameters or the topologies one might choose. The small and finite size of the system implies that bulk concentrations of molecular species cannot be fully attained. Nevertheless, a model to be acceptable should give at least a qualitatively correct picture of the micelles and that means that such ion condensation should not be observed. Instead, sodium ions should be only transiently bound to the micellar surface.
The flexible model also presented an ionic bridging between octanoate headgroups, but three counter-ions remained in solution after a 10 ns trajectory (Fig. 7). These free cations contributed to decrease the sodium density near the micelle center of mass (Fig. 4). The coordination number of counter-ions around the octanoate oxygen is stable after a ca. 1000 ps trajectory (Fig. 5), the same time lag in which density relaxed (Fig. 2).
The solvent structure around sodium ions relaxed along both trajectories. The number of water oxygen atoms in the first coordination shell around the sodium atoms decreases continuously for both models (Fig. 8). The decrease in coordination number is more pronounced for the rigid model, but clearly the values for the flexible model did not converge.
Altogether, these structural and relaxation data lead us to conclude that the OPLS parameters for sodium ion are not adequate to simulate the sodium octanoate micellar system. The main reason behind such anomalous behavior can be attributed to the fact that the force field favors ionic pairs instead of solvated ions in solution.
5 Åqvist parameters for sodium cation
In this section, we shall present the structural data of the sodium octanoate micelle in aqueous solution obtained using the Åqvist parameters for sodium ions. Åqvist parameters for alkaline and alkaline earth cations have been optimized to reproduce free energies of hydration for these ions in dilute solution. This approach somehow includes contributions of the bulk phase into the parameters, what does not happen if ab initio quantum chemistry calculations of dimers in vacuum are the main reference data, as it is the case for the OPLS parameterization process. The OPLS parameters also are optimized to reproduce some thermochemical solution data, but the results we have reported in the previous section indicate that the sodium ion parameters are not adequate to reproduce the microheterogeneous structure of an anionic micellar system.
The simulations with Åqvist parameters for sodium ion attained stable density values within a few picoseconds (Fig. 9), as opposed to the simulations conduced using the OPLS parameters (Fig. 2). The stable density suggests that no major structural change is taking place in the simulation box. This is indeed the case, as we shall demonstrate in the following discussion.
The final structure of the flexible model simulation presented one free monomer in solution (Fig. 10) while the rigid model simulation presented five (Fig. 11). Although the monomers dissociation might be regarded as major structural change, it does not affect the qualitative features characterizing the micellar structure as most of the hydrocarbon tails remained in a small, dry region, bearing the hydrophilic headgroups exposed to the aqueous medium.
The average density values are (1.023 ± 0.003) g cm3 for the flexible model and (1.002 ± 0.002) g cm3 for the rigid model. The difference reflects the fact that a more flexible molecule may assume a more compact configuration. It is interesting to note that both values are quite near the experimental density at this concentration, 1.011 g cm3 .
The radius of gyration also varied very little in the course of both simulations (Fig. 12), with an average value of (0.81 ± 0.03) nm for the flexible model and (0.94 ± 0.11) nm for the rigid model. Hayter and Zemb  reported a micellar radius of 1.17 nm from a neutron diffraction experiment. Zemb et al.  have reported a radius of 1.08 nm in a light scattering experiment. Friman and Rosenholm  obtained a larger micellar radius using an X-ray scattering experiment, 1.8 nm. The agreement is quantitative considering that parameters derived from scattering experiments are obtained assuming particular features of the micellar system structure. For instance, it is common to disregard any intermicellar interaction and to consider micelle size and shape to be uniform.
The counter-ion distribution around the micellar center of mass differed somewhat in the two simulations. The density is larger near the center of mass for the rigid model simulation, as compared to the flexible model simulation (Fig. 13).
The number of sodium atoms in the first coordination shell around the octanoate oxygen atoms confirmed the increase in the probability of finding a cation ion near the micelle in the rigid model simulation (Fig. 14). Nonetheless, the average number of sodium ions in the first coordination shell around octanoate oxygen atoms was rather low in both simulations, (0.08 ± 0.04) for the flexible model and (0.14 ± 0.06) for the rigid model.
The low coordination values contrast with the counter-ion condensation observed in the simulations with the OPLS parameters for sodium ion. The final structure of both simulations with Åqvist parameters for sodium ions presented no lattice-like structure, on the contrary, most counter-ions remained in solution (Figs. 10 and 11).
The number of water molecules in the first coordination shell around sodium atoms was very stable in both simulations (Fig. 15). The flexible model yielded an average value of (5.31 ± 0.09) water oxygen atoms and the rigid model an average of (5.25 ± 0.12) water oxygen atoms. This finding contrasts with both simulations with the OPLS parameters for sodium ions, in which the coordination numbers dropped to a value below 4.0 for the flexible model and below 2.0 for the rigid model (Fig. 8).
The molecular picture emerging from these data resembles the accepted view of anionic micelles in aqueous solution , suggesting that Åqvist parameters for sodium ions are adequate to simulate such systems. Of course, a quantitative agreement will depend on further improvements in the force field parameters as well as in the model systems.
6 Comparison with ab initio results
To further investigate the large differences between these two parameters sets ab initio calculations using second order Møller-Plesset perturbation theory with 6-311+G** basis set were performed for the sodium-water and sodium-propanoate dimers.
The sodium-water dimer was built as follows: (i) all atoms were in the same plane; (ii) the sodium atom approached the water oxygen along a line bisecting the HOH angle; (iii) the water molecule geometry was constrained to the SPC equilibrium values.
The sodium-propanoate dimer was build as follows: (i) all the heavy atoms were in the same plane; (ii) the sodium atom approached the carboxylic carbon atom along a line bisecting the OCO angle; (iii) the propanoate molecule geometry was constrained to the OPLS-AA equilibrium values.
These calculations were performed using the Gaussian 98  program. The ab initio potential energy profile obtained were compared with the force field energy profile for the same molecular geometries and the results are presented in Figs. 16 and 17. One observes that at large distances the energy profile calculated using both parameter sets agree with the ab initio energy profile for both dimers. However, at short distances, the force fields yielded steeper energy increases. One also observes that the energy obtained with the Åqvist parameters is in better agreement with ab initio results in the minima regions while the values calculated with OPLS parameters shows a deeper minimum for both dimers.
The energy difference between the two force field calculations at the minimum position was found to be ca. 1.8 kcal mol1 for the sodium-water dimer and ca. 5.6 kcal mol1 for the sodium-propanoate dimer. These energy differences explain the sodium cation condensation on the micellar surface in the simulations carried out with the OPLS force field. The average kinetic energy is ca. 0.9 kcal mol1 at 300 K, so the energy differences above are high enough to account for the degree of ionic association we found out in the Molecular Dynamics simulations we have performed with the OPLS parameters.
7 Concluding remarks
Molecular Dynamics is a powerful technique to unravel structural and dynamical features of complex liquids. Nonetheless, great care must be taken when choosing simulation conditions and interaction parameters.
Our simulation results point to the dependence of micellar structure on both the molecular topology and the interaction parameters employed. Both sodium ion parameters were developed to model dilute aqueous solutions, but Åqvist parameterization reproduced hydration free energies instead of hydration enthalpies, as is the case for the OPLS parameters.
The micellar structures obtained with Åqvist parameters for sodium cation are in good agreement with experimental and theoretical evidences. On the other hand, OPLS parameters resulted in the sodium condensation on the micellar surface, forming stable ionic cluster, in disagreement with the accepted micellar picture arising from experimental data.
Besides the effect of sodium ion parameters, we showed that molecular topology also affects the average micellar structure, but differences due to topology are less significant than those due to the force field parameterization.
We have evaluated the adequacy of two methodological issues to study complex liquids by means of computer simulation. There are many other features that should be evaluated to validate the simulation protocols. It is important to obtain accurate molecular models to describe complex liquids as the experimental techniques are also model dependent. The structural and dynamical information derived from reliable simulation models are complementary to the experimental data, since they provide a detailed description at time and distance intervals that cannot be reached by experimental means.
Gratitude is expressed to FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo) for financial support. One of us (AFM) also thanks FAPESP for the award of a scholarship.
 B. Jönsson, O. Edholm, and O. Teleman, J. Chem. Phys. 85, 2259 (1986). [ Links ]
 K. Watanabe, M. Ferrario, and M. L. Klein, J. Phys. Chem. 92, 819 (1988). [ Links ]
 K. Watanabe and M. L. Klein, J. Phys. Chem. 93, 6897 (1989). [ Links ]
 J. Shelley, K. Watanabe, and M. L. Klein, Electrochimica Acta 36, 1729 (1991). [ Links ]
 J. Shelley, M. Sprik, and M. L. Klein, Langmuir 9, 916 (1993). [ Links ]
 H. Kuhn and H. Rehage, Ber. Bunsenges. Phys. Chem. 101, 1485 (1997). [ Links ]
 H. Kuhn and H. Rehage, Ber. Bunsenges. Phys. Chem. 101, 1493 (1997). [ Links ]
 H. Khun, B. Breitzke, and H. Rehage, Colloid Polym. Sci. 276, 824 (1998). [ Links ]
 L. Laaksonen and J. B. Rosenholm, Chem. Phys. Lett. 216, 429 (1993). [ Links ]
 K. Watanabe and M. L. Klein, J. Phys. Chem. 95, 4158 (1991). [ Links ]
 J. C. Shelley, M. Sprik, and M. L. Klein, Progr. Colloid Polym. Sci. 103, 146 (1997). [ Links ]
 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gusteren, and J. Hermans, in Intermolecular Forces, B. Pullman, Ed., Reidel, Dordrecht, 1981. [ Links ]
 W. L. Jorgensen, J. Am. Chem. Soc. 103, 335 (1981). [ Links ]
 J. Chandrasekhar, D. C. Spellmeyer, and W. L. Jorgensen, J . Am. Chem. Soc. 106, 903 (1984). [ Links ]
 H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987). [ Links ]
 P. K. Weiner and P. A. Kollmann, J. Comp. Chem. 2, 287 (1981). [ Links ]
 S. J. Weiner, P. A. Kollmann, D. A. Case, U. C. Singh, C. Ghio, G. Alagona, S. Profeta, and P. Weiner, J. Am. Chem. Soc. 106, 765 (1984). [ Links ]
 S. J. Weiner, P. A. Kollmann, D. T. Nguyen, and D. A. Case, J. Comp. Chem. 7, 230 (1986). [ Links ]
 B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus, J. Comp. Chem. 4, 187 (1983). [ Links ]
 W. L. Jorgensen, J. Chandrasekhar, J. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983). [ Links ]
 S. Bogusz, R. M. Venable, and R. W. Pastor, J. Phys. Chem. B 104, 5462 (2000). [ Links ]
 T. Wymore, X. F. Gao, and T. C. Wong, Journal of Molecular Structure 485-486, 195 (1999). [ Links ]
 D. P. Tieleman, D. van der Spoel, and H. J. C. Berendsen, J. Phys. Chem. B 104, 6380 (2000). [ Links ]
 S. J. Marrink, D. P. Tieleman, and A. E. Mark, J. Phys. Chem. B 104, 12165 (2000). [ Links ]
 J. J. Wendoloski, S. J. Kimatian, C. E. Schutt, and F. R. Salemme, Science 243, 636 (1989). [ Links ]
 S. Balasubramanian and B. Bagchi, J. Phys. Chem. B 106, 3668 (2002). [ Links ]
 J. Shelley, K. Watanabe, and M. L. Klein, Int. J. Quantum Chem., Quantum Biol. Symp. 17, 103 (1990). [ Links ]
 A. D. MacKerell Jr., J. Phys. Chem. 99, 1846 (1995). [ Links ]
 M. Tarek, S. Bandyopadhyay, and M. L. Klein, Journal of Molecular Liquids 78, 1 (1998). [ Links ]
 T. Bast and R. Henttschke, J. Phys. Chem. B 100, 12162 (1996). [ Links ]
 J. Åqvist, J . Phys. Chem. 94, 8021 (1990). [ Links ]
 H. J. C. Berendsen, D. van der Spoel, and R. van Drunen, Comp. Phys. Comm. 91, 43 (1995). [ Links ]
 E. Lindahl, B. Hess, and D. van der Spoel, J. Mol. Mod. 7, 306 (2001). [ Links ]
 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gusteren, A. Dinola, and J. R. Haak, J. Phys. Chem. 81, 3684 (1984). [ Links ]
 W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, J. Am. Chem. Soc. 118, 11225 (1996). [ Links ]
 B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije, J. Comput. Chem. 18, 1463 (1997). [ Links ]
 S. Miyamoto and P. A. Kollmann, J. Comput. Chem. 13, 952 (1992). [ Links ]
 T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993). [ Links ]
 U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995). [ Links ]
 P. Ekwall, H. Eikrem, and L. Mandell, Acta. Chem. Scand. 17, 111 (1963). [ Links ]
 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, V. G. Zakrzewski, J. A. Montgomery, Jr., R. E. Stratmann, J. C. Burant, S. Dapprich, J. M. Millam, A. D. Daniels, K. N. Kudin, M. C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G. A. Petersson, P. Y. Ayala, Q. Cui, K. Morokuma, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. Cioslowski, J. V. Ortiz, A. G. Baboul, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, J. L. Andres, C. Gonzalez, M. Head-Gordon, E. S. Replogle, and J. A. Pople, Gaussian 98, Revision A.7, Gaussian, Inc., Pittsburgh PA, 1998. [ Links ]
 J. N. Israelachvili, Intermolecular & Surface Forces, 2nd Ed., Academic Press Limited, London, 1995. [ Links ]
 J. B. Hayter and T. Zemb, Chem. Phys. Lett. 93, 91 (1982). [ Links ]
 T. Zemb, M. Drifford, M. Hayoun, and A. Jehanno, J. Phys. Chem. 87, 4524 (1983). [ Links ]
 R. Friman and J. B. Rosenholm, Colloid & Polymer Sci. 260, 545 (1982). [ Links ]
Received on 01, September, 2003