752
J. Phys. Chem. B 2009, 113, 752–758
Force Field for Molecular Dynamics Studies of Glycine/Water Mixtures in Crystal/Solution Environments Sivashangari Gnanasambandam,† Zhongqiao Hu,† Jianwen Jiang,†,‡ and Raj Rajagopalan*,†,‡ Department of Chemical and Biomolecular Engineering, National UniVersity of Singapore, Singapore 117576, and Chemical and Pharmaceutical Engineering Program, The Singapore-MIT Alliance, Singapore 117576 ReceiVed: April 5, 2008; ReVised Manuscript ReceiVed: NoVember 10, 2008
Defining a force field for the theoretical predictions of the structure and growth of pharmaceutical crystals is very difficult because of the complex intermolecular bonding found in these crystals. Here, we investigate the accuracy of the AMBER ff03 force field for R-glycine molecules using molecular dynamics. The validation of the force field is carried out in both solution and crystal/solution environments. The molecule R-glycine is chosen because it has a simple molecular structure while displaying complex intermolecular interactions typical of pharmaceutical crystals. We estimate the atomic charge for glycine from density-functional theory and represent water by the extended simple point charge (SPC/E) model. In a pure solution environment, our simulation results agree well with the experimental data available for solution density, radial distribution functions, hydration number, and diffusion coefficient of glycine. In the crystal/solution environment, the lattice energy calculations are in agreement with the available literature results and also the force field is able to identify the hydrophilic nature of the (010) surface of glycine crystal. The good agreement between simulation results and experimental data indicates that the chosen force field can be used to investigate the growth of glycine crystals from aqueous solutions. 1. Introduction Understanding the molecular mechanisms of crystal growth is an essential step toward controlling crystal growth, morphology, and shape, which are of prime interest and importance in pharmaceutical industries. Pharmaceutical crystals interact with their surroundings predominantly through their surfaces and, as a result, the shape and morphology of a crystal influence its chemical and physical properties such as dissolution rate and stability of a potential drug. Hence, a complete understanding of the crystal growth mechanism is vital for selecting the most suitable drug product. Crystal growth kinetics can be studied in a number of ways depending on the purpose and the scope of the investigation. Our interest in the present paper is on the growth kinetics at crystal/liquid interfaces. In recent years, a number of modern experimental probes have been employed to specifically focus on the interface. For example, atomic force microscopy (AFM) has been widely used to monitor the crystal/liquid interface on the subnanometer scale.1 Very recently, grazing incidence X-ray diffraction (GID) has become a powerful tool for getting information on the atomic and molecular levels.2 A number of theoretical studies have also focused on phenomena at crystal interfaces and their effects on crystal growth.3 These include the work of Gadomski and Sio´dmiak4 and Gadomski et al.,5 who studied the kinetics of crystal growth using a mass convection approach, instead of bulk diffusion, and the dependence of the growth on the electrostatic double layer created by charged ions at the interface. However, it is very useful to complement experimental and analytical (theoretical) studies such as the above by using molecular level computations, since such an approach provides a powerful means to identify specific mechanisms as well as to * Corresponding author. † National University of Singapore. ‡ The Singapore-MIT Alliance.
uncover details that are not always accessible in physical experiments. In this sense, computer experiments at an atomistic level play a key complementary role to physical experiments by making predictions and offering guidance prior to experimental studies. Computational modeling of a crystal requires, as the first step, a description of charge distribution in the molecule and selection of appropriate intermolecular force fields.6 Even though the physical nature of inter- and intramolecular forces acting between the molecules in the crystals are generally well-known, the known force fields all have implicit or explicit approximations and their suitability for computational modeling often rests on the type of properties of bulk phases needed and the type of physical phenomena investigated. Therefore, it is essential that a systematic examination of the force field be undertaken first prior to investigating the bulk structures and phases using the chosen force field. For the kinetics of crystallization and morphology of resulting crystals we seek to examine subsequently, we choose glycine (+H3NCH2COO-) as a model substance because of its simple molecular structure and the importance of its prochiral property and since the hydrogen bonding in its crystal structure is similar to that found in proteins and polypeptides. Moreover, glycine can act as a vital neurotransmitter in the central nervous system.7 Hence, a thorough knowledge about the behavior of glycine in solution and in the crystalline state could provide a glimpse of its interactions with neural receptors. As a simple amino acid, glycine usually crystallizes as the R-polymorph form from aqueous solutions. It adopts a zwitterionic form in aqueous solutions and in the crystalline state,8 but exists in neutral form in the gas phase.9 Selection of an appropriate force field for glycine molecule is quite challenging. For thermodynamically accurate studies, the force field for glycine/water mixture should satisfy the following two conditions: (1) the same force field should be sufficiently accurate for both the crystals and liquids, and (2) the force field chosen
10.1021/jp802949u CCC: $40.75 2009 American Chemical Society Published on Web 12/30/2008
MD Studies of Glycine/Water Mixtures should reproduce key (known) properties (relevant for the phenomena under investigation) sufficiently accurately. To explore the internal structure of glycine in a solution environment, many ab initio calculations have been carried out in the presence of water and ions. For example, Leung and Rempe10 studied the glycine intramolecular proton transfer in water by ab initio molecular dynamics simulations. They computed the potential of mean force associated with the direct intramolecular proton transfer event in glycine and found that the average hydration number of glycine in zwitterionic form was 8. Campo11 explored the hydration and internal structure of glycine using classical molecular dynamics (MD) with the use of the GROMOS9612 force field for glycine and solvent (water). The partial charges were calculated with the atoms-in-molecules (AIM) method13 using GAUSSIAN.14 He reported radial distribution functions, diffusion coefficients, and hydration number of glycine and water in the presence of ions. These earlier studies revealed the structural features of glycine in a solution environment by applying a specific force field, but they did not examine the applicability of that force field in crystalline environments. Some studies do exist for the crystalline environment. For example, the pioneering work for the force field in a crystalline environment was done by Momany et al.15 They presented intraand intermolecular interactions for a class of amino acids in the gas phase, where the acids are neutral, and then extended them to the solid state, where the acids are zwitterions. The disadvantage of this method is that it leads to a low binding energy due to the underestimation of the electrostatic interactions between the amino acid molecules in the solid state. For amino acids in crystalline form, the first force field was developed by Scheraga et al.,16 and it is an improved version of the force field of Momany et al., as it (the former) uses the available experimental information. Scheraga et al. used partial charges calculated by the molecular orbital CNDO/2 method17 and used atom-atom approximation to sum the pairwise interatomic interactions to calculate the intermolecular potential energy. The major drawback of this force field is that it consists of harmonic terms to describe the bonded interactions which are initially parametrized from experimental data. The second force field in this area is the Dreiding force field, developed by Mayo et al.18 Here, the charges were calculated using the Gastriger-Marsili method.19 The Dreiding force field overcomes the limitation of the Scheraga force field, as it has a broader applicability for a large group of compounds with only limited experimental data or without any data. However, it has the drawback of being of only moderate precision in the energy calculations. Thus, the calculated properties of compounds differ noticeably even with subtle differences in the interatomic potential. The aim of the present study is to examine the suitability of the AMBER ff03 force field20 for water/zwitterionic glycine mixtures in solution as well as in crystalline environments for subsequent studies of crystal growth. To this end, we study, using molecular dynamics simulations, the chosen force field for glycine and water (with the latter represented by SPC/E model) for its ability to predict the necessary properties in both solution and crystal interface environments. [SPC/E stands for the extended simple point charge water model.] For the solution environment, simulated values for density, radial distribution function, hydration number, and glycine diffusivity are compared with experimental values as functions of glycine mole fraction. For the crystal interface environment, the enthalpy of sublimation and proton transfer energies are compared with experimental and other computed values available in the literature.
J. Phys. Chem. B, Vol. 113, No. 3, 2009 753
Figure 1. A zwitterionic glycine molecule.
TABLE 1: Mulliken Atomic Charges of Zwitterionic Glycine atom index
atom
charge
1 2 3 4 5 6 7 8 9 10
N H1 H2 H3 CA HA1 HA2 C OC1 OC2
-0.127 0.199 0.218 0.224 0.007 0.064 0.061 0.483 -0.578 -0.552
This paper is organized as follows. In section 2, the molecular models for glycine and water are explained together with the details about the simulation setup and protocols for both the solution and crystal interface environments. The results obtained from the simulation are discussed and compared with experimental values for the thermodynamic and dynamic properties of the systems in section 3. Finally, a summary of observations is presented in section 4, along with relevant concluding remarks. 2. Models and Methodology As noted above, we begin with a description of the models used for glycine and water and further elaborate on the simulation methodologies used for the solution and crystal interface environments. 2.1. Models for Glycine and Water. Figure 1 shows a zwitterionic glycine molecule, which consists of three groups: the hydrophilic -COO- and -NH3+ groups, and the hydrophobic -CH2 group. We describe the intramolecular and intermolecular dispersive interactions of glycine by the AMBER ff03 force field.20 AMBER ff03 is a newly developed, thirdgeneration point-charge all-atom force field for organics and biomolecules. [“Third generation” refers to the fact that the force field parameters include polarization effect in the condensed phase (solvent).] The intermolecular Coulombic interactions are estimated on the basis of the Mulliken partial charges, as listed in Table 1. The charges are calculated from density functional calculations using Dmol3 in Material Studio (Accelrys, Materials Studio Release Notes, Release 4.1; Accelrys Software, Inc.: San Diego, CA, 2006), in which the Becke exchange plus Lee-Yang-Parr (BLYP)21,22 correlation functional and the double-numerical plus d- and p-polarization basis set (DNP) are used. As noted earlier in section 1, the purpose of the present paper is to examine the Amber ff03 force field and the parameters needed for this force field for subsequent use in studying the kinetics of growth and dissolution of the crystalline phase. Therefore, it is important that the chosen parameters are capable of reproducing the solubility of glycine sufficiently
754 J. Phys. Chem. B, Vol. 113, No. 3, 2009
Gnanasambandam et al.
TABLE 2: Glycine/Water Mixtures Simulated in the Solution Environmenta no. of molecules system no.
nglycine
nwater
nsystem
xglycine
1 2 3 4 5 6
0 1 8 27 64 125
511 506 858 1300 1882 2453
511 507 864 1327 1946 2578
0 0.0020 0.0093 0.0203 0.0329 0.0485
a nglycine: number of glycine molecules. nwater: number of water molecules. nsystem: total number of molecules in the system. xglycine: mole fraction of glycine.
accurately. (As an additional check, we have also chosen the enthalpy of sublimation of glycine crystals as a reference property.) Tests of Amber ff03 with a number of charge models including ESP (DNP) charges, Mulliken (DNP) charges, and Mulliken (CNDO/2) charges, in addition to RESP (HF/6-31G**) charges, show that Mulliken (DNP) charges reproduce the above two properties most accurately. In particular, the solubility (i.e., the saturation concentration) at 25 °C based on Mulliken (DNP) charges differs from the experimental value (249.9 ( 2 g/1000 g water; Dalton and Schmidt23) by only 10%, whereas the other charge models lead to 20 to 30% errors, with the largest (31%) coming from the RESP charges. Similarly, the estimated value of the enthalpy of sublimation based on RESP (HF/6-31G**) model for the charges is about one-third of the experimental value (136 ( 4 kJ/mol; Svec and Clyde24), whereas the Mulliken (DNP) charges lead to only about +19% error. For these reasons, we have chosen to use the Mulliken (DNP) charges, although this choice is semiempirical in nature. The well-known models for water are SPC,25 SPC/E,26 TIP3P,27 and TIP4P28 [TIP stands for transferable intermolecular potential functions]. Here, we have chosen to use the SPC/E model, since, when compared with other water models, it reproduces the experimental diffusivity data29 very well. (Note that accurate values for diffusivity are important in crystal growth studies, as the rate of transfer of solute from solution to crystal and vice versa are determined by diffusion.) 2.2. Simulation Methodology. 2.2.1. Solution EnWironment. For the solution environment, we consider five different glycine/water mixtures with the glycine mole fraction ranging from 0 to 4.85%, as listed in Table 2. Pure water is also included for comparison. MD simulations were carried out in the isothermal-isobaric (NPT) ensemble using Gromacs.30 Periodic boundary conditions were applied in the three directions for all the systems in order to simulate bulk systems. Each system was coupled to a Berendsen thermostat at 300 K and a Berendsen barostat31 at 1 atm with a relaxation time of 1 ps. The time step for the simulation was 2 fs. The Particle-Mesh Ewald (PME) method32 was used for treating the long-range electrostatic interactions with a cutoff radius of 0.9 nm. For Lennard-Jones interactions, the cutoff was 1 nm. The system was initially subjected to energy minimization for a total of 150 steps. The first 50 steps were carried out with positional restraints on the glycine molecules, followed by 100 steps without restraints. Thereafter, the velocities of all molecules including glycine and water were assigned according to the Boltzmann distribution at 100 K. The system was heated up and equilibrated for 1 ns, and the potential energy of the system was monitored to ensure equilibration. For each system, an additional NPT simulation of 10 ns production run was performed. The trajectory was saved every 0.1 ps for subsequent analysis.
Figure 2. A typical snapshot of water molecules above the (001) surface of the R-glycine crystal obtained from the simulations. (a) Top view normal to the (010) plane. (b) and (c) Side views normal to (100) and (001) planes. Water molecules are indicated by red spheres (O atoms) connected by gray spheres (H atoms).
2.2.2. Crystal Interface EnWironment. The (010) surface of R-glycine crystal was selected for the present study, as shown in Figure 2a-c. The R-glycine crystal belongs to a monoclinic space group of P21/n, and the lattice parameters are a ) 5.1054 Å, b ) 11.9688 Å, c ) 5.4645 Å, and β ) 111.697°.33 We used a thin slab of crystal consisting of a 5 × 5 supercell in the x and z directions with a depth of 3 unit cells in the y direction (corresponding to a thickness of 3.59 nm). Water molecules were added above the slab in the y direction for a height of 8.0 nm (i.e., more than 8 times the cutoff length). The use of the periodic boundary conditions in all three directions then implies that both sides of the slab are exposed to the solvent, with each exposed crystal surface corresponding to an infinitely large (010) surface. The above arrangement corresponds to a simulation box consisting of 300 glycine molecules and 1710 water molecules. We followed the same simulation methodology as in the solution environment, but an anisotropic pressure coupling in the y direction was used in this case, because anisotropic coupling is more appropriate than isotropic coupling, for the study of both solid and liquid phases.34 After the energy minimization, the glycine molecules were positionally restrained35 during both the equilibrium and the production steps while we examined the behavior of water at the glycine/water interface at 300 K and 1 atm. Applying positional restraint implies that the mobility of glycine in the crystal (including at the surface) is restrained, and this, in principle, could affect the calculated properties. Nevertheless, Boek et al.36 have already shown, by examining the effect of positional restraint on a urea/ water system, that the positional restraint has very little influence on the density and diffusivity of water both in the vicinity of the interface and beyond. Therefore, we expect the positional restraint in our case to have a negligible impact on the results. 3. Results and Discussion We report here the results of the molecular dynamics simulations in the above two environments at 300 K and 1 bar. First, the force field is examined in the solution environment to study the structural details and molecular properties of glycine in solution. Then the force field is tested for the crystal interface environment to check its suitability in investigating crystal growth. 3.1. Solution Environment. In the case of the solution environment, we compare the simulated density, glycine diffusivity, radial distribution functions, and hydration number to experimental data as functions of glycine molar fraction as discussed in detail below. 3.1.1. Density. The density of the glycine/water mixture as a function of glycine concentration at 300 K is reported in Table 3 and plotted in Figure 3. Some experimental data23 available in the literature (at 298 K) for glycine/water mixture are included
MD Studies of Glycine/Water Mixtures
J. Phys. Chem. B, Vol. 113, No. 3, 2009 755
TABLE 3: Comparison of Densities from Simulations (This Work) and Experiments (Dalton and Schmidt23) density (g/cm3) xglycine
simulationa
experimentb
0 0.0020 0.0093 0.0203 0.0329 0.0485
1.013 1.015 1.026 1.042 1.059 1.079
0.998 1.000 1.013 1.031 1.051 1.084
a Standard deviation in the simulation is less than 1%. b Standard deviation in the experiment is less than 1.5%.
Figure 3. Comparison between experimental and simulated densities of glycine/water solutions.
for comparison. For pure water, Wu et al.29 report the density based on the SPC/E model from molecular dynamics simulations as 0.999 ((0.017) g/cm3. The density for pure water from our simulation is 1.013 g/cm3, and differs from experimental value by 0.014 g/cm3, which is within the statistical error (0.017 g/cm3) reported by Wu et al. With increasing mole fraction xglycine, the density of the mixture increases slightly. It can be seen that the force field used in our simulations predicts the density of glycine/water mixtures well. The agreements between the simulated and experimental densities are good particularly at high xglycine when the contribution from glycine-glycine interaction plays an increasingly important role in the density calculation. The accuracy of the predicted density is sufficient for further use of the selected force field in crystallization studies. 3.1.2. Radial Distribution Functions. The radial distribution functions describe how the atoms are distributed radially around each other. Here, they are used to test the ability of the chosen force field to capture the structural features at the atomic level. Figure 4 shows the radial distribution functions between the atoms OC, H, and N in glycine and the oxygen atom in water (OW). There are two OC atoms (OC1 and OC2) in glycine; therefore, the calculated radial distribution functions for OC were averaged over the two OC atoms. Similarly, for the H atom, the average was taken over the three H atoms H1, H2, and H3. The results show that the locations of the first peaks in the radial distribution functions, denoted here by rOC-OW, rH-OW, and rN-OW, are 2.71, 2.05, and 2.95 Å, respectively. Each of these represents the position of first hydration shell for the corresponding atom in glycine. These values of rOC-OW, rH-OW, and rN-OW are in good agreement with the average hydrogen bond distances rO · · · OW ) 2.72 Å, rH · · · OW ≈ 2.0 Å, and rN-H · · · OW ) 2.89 Å reported for various organic crystals.37 One may note that, as shown in Figure 4, the heights of the peaks decrease slightly with increasing number of glycine molecules, since this increase perturbs the number of water molecules in the hydration shell.
Figure 4. Radial distribution functions between (a) OC, (b) H, and (c) N of glycine and OW of water in glycine solution. The arrows show the direction of increase in xglycine.
TABLE 4: Hydration Numbers of Carboxylic and Amino Groups of the Glycine Molecule xglycine
-COO-
-NH3+
0.0020 0.0093 0.0203 0.0329 0.0485
5.38 5.12 4.86 4.60 4.28
2.97 2.82 2.70 2.55 2.40
3.1.3. Hydration Number. The hydration number of a given molecule is defined as the number of water molecules around that molecule. Table 4 gives the hydration numbers of the amino and carboxylic and groups in glycine. At infinite dilution, the hydration numbers are 2.97 and 5.38, respectively, for the amino and carboxylic groups. On average, there are 3 water molecules around the amino group and 4 to 5.5 water molecules around the carboxylic group. This observation represents the more hydrophilic nature of the carboxylic group than the amino group in glycine. Further, the hydration number decreases with the increase in the concentration of glycine. This is attributed to the enhanced solute-solute interactions, which perturb the hydration shell to some extent. The calculated hydration numbers at infinite dilution agree well with the reported literature values. Kameda et al.38 reported 3 water molecules around the amino group from a neutron
756 J. Phys. Chem. B, Vol. 113, No. 3, 2009
Gnanasambandam et al. TABLE 5: Self-Diffusion Coefficients of Glycine and Water from Simulations (This Work) and Experiments29,42,43 D (10-5 cm2 s-1) xglycine a
0 0.0020 0.0093 0.0203 0.0329 0.0485
Figure 5. Mean-squared displacements of glycine in a solution. The arrow shows the direction of increase in xglycine.
diffraction measurement with H/D isotopic substitution method and also predicted the number of water molecules around the carboxylic molecule to be about 4 to 4.5 as found in formate.39 In a Monte Carlo simulation study, Alagona et al.40 found 2.3 and 5.3 water molecules around the amino and carboxylic groups, respectively, using the charges obtained from the best fit to the electrostatic potential produced by a 6-31G* basis set.41 Using ab initio molecular dynamics calculations and the BLYP21,22 correlation functional, Leung and Rempe10 obtained 3 and 4.7 water molecules around the amino and the carboxylic groups, respectively. Campo11 found 2.8 and 5.4 water molecules around the amino and carboxylic groups from MD simulations, using the GROMOS96 force field (as noted previously). Thus, the predicted values of hydration numbers in our calculations are very close to the experimental values found in the literature and those computed by independent methods. Therefore, we conclude that the force field chosen here predicts the structural features around the glycine molecules in a solution very accurately. 3.1.4. Self-DiffusiWity of Glycine. One of the key transport properties of glycine needed in crystallization studies is the selfdiffusivity. The self-diffusivity D(c) of glycine at a concentration c can be evaluated using the Einstein expression
D(c) ) lim tf∞
1 〈|r(t) - r(0)| 2〉 6t
(1)
where r(t) is the position of the molecule at time t. The ensemble-averaged mean-squared displacement (MSD) is represented by 〈|r(t) - r(0)|2〉. Figure 5 shows the MSDs of pure water and of glycine in mixtures as functions of time. We have determined the diffusivity from the slopes of the MSDs in the time range of 100-300 ps. These are reported in Table 5. The magnitudes of self-diffusion coefficient of pure water we obtain are in good agreement with the independent simulation result of Wu et al.,29 namely, (2.41 ( 0.08) × 10-5 cm2/s, as well as with the experimental value reported by Krynicki et al.,42 i.e., 2.3 × 10-5 cm2/s. For glycine in water, Figure 6 shows the comparison between simulated and experimental diffusivities43 as a function of mole fraction of glycine. As seen, the diffusivity of glycine decreases with increasing concentration. In general, the agreement between simulation and experiment improves with the glycine concentration, which shows the accuracy of the force field in predicting the dynamic properties of glycine in solution. (Please note that the value at infinite dilution is difficult to predict precisely statistically as the number of glycine molecules is very small at very low concentrations.) 3.2. Crystal Interface Environment. The lattice energy of a crystal is an important property to assess crystal/solution
simulationb
experimentc
2.01 1.21 1.03 0.93 0.78
2.30d 1.39 0.97 0.96 0.98 0.67
a Diffusion coefficients given are those of water. b Standard deviation for simulation is less than 1%. c Standard deviation for experiment is 1% (Ma et al.43). d Standard deviation for experiment is 5% (Krynicki et al.42).
Figure 6. Comparison between experimental and simulated diffusion coefficients of glycine in solution. Lines are to guide the eye.
kinetics and equilibrium. We first estimate the lattice energy of R-glycine crystal, and subsequently the sublimation enthalpy, which can then be compared with the available experimental values to validate the force field in crystalline environments. In addition to the sublimation enthalpy, we also examine the water density profile close to the crystal/liquid interface in order to provide more information about the nature of the surface for crystal growth studies. Moreover, calculations of diffusivity are also presented. 3.2.1. Lattice Energy and Sublimation Enthalpy. The lattice energy Ecry is the configurational energy of a molecule in the crystal. An accurate estimation of the lattice energy is very important for crystal growth studies as it determines the transfer of a molecule from a crystal to the solution (dissolution of the crystal) or depositing on a crystal surface from the solution (crystal growth). As the lattice energy cannot be measured directly experimentally, the validity of the force field has to be tested by comparing the theoretical values of the sublimation enthalpy (obtained from the calculated lattice energy) with the experimental sublimation enthalpy.44 The relation between the lattice energy and the sublimation enthalpy is described below. First, we note that our calculations for lattice energy lead to a value of -306 kJ/mol, obtained through energy minimization of glycine crystal structure. This value matches well with the simulated value reported by Boek et al.45 (also -306 kJ/mol), who used the same Mulliken population analysis for charges, as we do, but with a double-zeta plus polarization (DZP) basis set and Buckingham potential for the nonbonded interactions. Our calculated value is also very close to the value reported by No et al.46 (-285 kJ/mol), who used the same Mulliken population analysis for charges but with different basis set (631G* basis set). We now compare sublimation enthalpies obtained theoretically with experimental values. The lattice energy and sublima-
MD Studies of Glycine/Water Mixtures
J. Phys. Chem. B, Vol. 113, No. 3, 2009 757
tion enthalpy are connected through the protonation of the glycine molecules, which exist in the zwitterionic state in crystals8 and in the nonzwitterionic in the gas phase.9 Therefore, during crystallization from the gas phase or during sublimation, a proton is transferred from -NH3+ to -COO- and the energy required for this process is called the proton transfer energy Upr, and this needs to be incorporated in the calculation of the enthalpy of sublimation Hsub, as shown in eq 2:
Hsub ) -Ecry + Upr - 2RT
(2)
The term 2RT represents a correction for the difference between the vibrational contribution to the crystal enthalpy of 6RT47 and the gas-phase enthalpy (pV + 3RT). Our calculations of Upr using the density functional theory with the BLYP functional and DNP basis set leads to a value of -138.22 kJ/mol and falls within the range of -133 to -178 kJ/mol reported by Voogd et al.48 for the aliphatic R-amino acids. Also, Voogd et al.48 report that the values for Upr of glycine may range from -136.8 to -167 kJ/mol depending on the basis set used for estimating the partial charges. Therefore, the force field we use here leads to a reasonable estimate for Upr. Using the calculated Upr and Ecry in eq 2, we arrive at a sublimation enthalpy of 162.8 kJ/mol, which agrees with the experimental value 136 kJ/mol reported in the literature24 with an experimental error of about (4 kJ/mol. Slightly higher experimental values have been reported by others (e.g., 145 kJ/mol cited by Boek et al.45 without attribution and 143 kJ/mol by Docherty44), but our calculated value does not differ significantly from the experimental values when one considers the experimental and theoretical uncertainties involved. The above observations support the applicability of the chosen force field in crystal growth studies since the calculated sublimation enthalpy is closer to the available experimental data. Moreover, some previous studies in glycine crystals ignored the significance of proton transfer energy in the sublimation enthalpy calculation,49 but our calculations emphasize the importance of including the proton transfer energy in calculating sublimation enthalpy of glycine, as also noted by No et al.46 and Boek et al.45 3.2.2. Density Profile of Water RelatiWe to the Crystal Surface. We calculated the density distribution of the water normal to the (010) surface in order to examine the solvent structure relative to the surface. The simulation box is symmetric in the y direction due to periodic boundary conditions and crystal symmetry. Crystalline glycine molecules are located in the range from 0 to 3.5 nm in the y direction. The density was calculated by dividing the simulation box into thin slices, each with a thickness of 0.3 nm in the y direction and plotted against the distance as shown in Figure 7. The density profile obtained reveals a strong peak near the interface, and this peak is then followed by a small peak away from the surface. Further, the local density reaches the bulk density quickly as one moves away from the interface. The heights and the positions of the density oscillations in the vicinity of the crystal interface are insensitive to the changes in the bulk density, and therefore the density peaks can be safely described as corresponding to the adsorption of water at the interface, thereby indicating that the (010) surface of the glycine crystal is hydrophilic. This result provides essential information for carrying out crystal growth studies on the (010) surface, and such information can be useful for finding the effect of solvent on this surface.
Figure 7. Water density profile in the y-direction, normal to the (010) surface of the R-glycine crystal.
Figure 8. Mean-squared displacements of water in the x, y, and z directions and along the xz plane in the vicinity of the (010) surface (in a layer of 0.3 nm in thickness).
TABLE 6: Diffusion Coefficients of Water (10-5 cm2/s) in the x, y, and z Directions over the (010) Surface of r-Glycine Crystala layer
Dx
Dy
Dz
Dxz
layer 1 layer 2 layer 3 bulk
0.341 0.991 1.650 2.211
0.018 0.028 0.065 1.878
0.298 0.982 1.606 2.212
0.320 0.989 1.629 2.211
a
Dxz denotes the diffusion coefficient in xz plane.
3.2.3. Water DiffusiWity. Because of the interaction of a water molecule with the crystal surface, it is expected that the mobility of water would depend on the distance from the crystal surface. In order to examine the influence of the interface on the diffusivity of water, we divide the “water bath” over the crystal surface into four layers, as indicated in Figure 7. Layer 1 is the closest to the (010) surface, followed by layers 2, 3, and bulk. Figure 8 shows the mean-squared displacements of the water molecules in layer 1. As seen, the displacement in the y direction reaches a plateau after approximately 5 ps, indicating the relatively marginal mobility of water normal to the surface in the vicinity of the crystal surface. In contrast, the displacements in x and z directions, which are parallel to the surface, are far greater. The diffusivities calculated from Figure 8 are listed in Table 6. As evident from Table 6 (or, equivalently, from Figure 8), we observe that in all the layers Dxz ) D| (parallel to the surface) is larger than Dy ) D⊥ (perpendicular to the surface). The magnitude of diffusivity increases when going from layer 1 to 3, and the difference is smaller in the parallel and normal directions. Finally, as one moves away from the interface, one recovers the diffusivity of bulk water. The diffusivity calculation
758 J. Phys. Chem. B, Vol. 113, No. 3, 2009 again supports the hydrophilic nature of the (010) surface of glycine crystal as observed through the density profile of water. 4. Concluding Remarks In summary, we have examined, using molecular dynamics simulations, the suitability of the AMBER ff03 force field and partial charges based on density-functional theory to investigate glycine/water mixtures in both the solution and crystal interface environments. In the solution environment, the chosen force field predicts the physical properties such as density and diffusivity and structural properties such as the radial distribution functions and hydration number in good agreement with data available in the literature. In the crystal interface environment, the calculated lattice energy and proton transfer energy agree well with previously reported computational results. The calculated sublimation enthalpy is also predicted within experimental errors. Most importantly, the force field is able to elucidate the nature of the crystal surface (namely, that it is hydrophilic), based on the calculation of the solvent density profile relative to the crystal interface and the solvent diffusivity in the vicinity of the interface. Hence, we conclude that the chosen force field is suitable to describe glycine molecules in the solution as well as in the crystalline environments. In addition, the results lend confidence to the practical use of the force field in further study of glycine crystal growth from aqueous solutions. Acknowledgment. The authors are grateful for the support from the National University of Singapore (R279-000-238-112) and from the Chemical and Pharmaceutical Engineering Program of the Singapore-MIT Alliance. References and Notes (1) Binning, C. F. Q.; Gerber, C. Phys. ReV. Lett. 1986, 56, 930. (2) Gidalevitz, D.; Feidenhans’l, R.; Matlis, S.; Smilgies, D. M.; Christensen, J. M.; Leiserowitz, L. Angew. Chem., Int. Ed. Engl. 1997, 36, 955. (3) Interfacial Dynamics; Kallay, N., Ed.; Surfactant Science Series, Arthur, T. H., Ed.; Marcel Dekker, Inc.: New York, 1999; Vol. 88. (4) Gadomski, A.; Sio´dmiak, J. Cryst. Res. Technol. 2002, 37, 281. (5) Gadomski, A.; Sio´dmiak, J.; Santamaria-Holek, I.; Rubi, J. M.; Ausloos, M. Acta Phys. Pol., B 2005, 36, 1537. (6) Peters, M. H. Real-time Biomolecular Simulations: the behaVior of biological macromolecules from a cellular systems perspectiVe; McGrawHill Professional: New York, 2007. (7) Wolff, M. In Burger’s Medicinal Chemistry, 4th ed.; Wiley: New York, 1979. (8) Takagi, S.; Chihara, H.; Seki, S. Bull. Chem. Soc. Jpn. 1959, 32, 84. (9) Bonaccorsi, R.; Palla, P.; Tomasi, J. J. Am. Chem. Soc. 1984, 106, 1945. (10) Leung, K.; Rempe, S. B. J. Chem. Phys. 2005, 122, 184506. (11) Campo, M. G. J. Chem. Phys. 2006, 125, 114511. (12) GROMOS96, Groningen Molecular Simulation Package; Biomos b.v: Zurich, 1996.
Gnanasambandam et al. (13) Bader, R. F. W. Atoms in Molecules:A Quantum Theory; Oxford University Press: Oxford, UK, 1990. (14) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B. GAUSSIAN 98, Revision A.9; Gaussian, Inc: Pittsburgh, 1998. (15) Momany, F. A.; Carruthers, L. M.; Scheraga, H. A. J. Phys. Chem. 1974, 78, 1621. (16) Nemethy, G.; Pottle, M. S.; Scheraga, H. A. J. Phys. Chem. 1983, 87, 1883. (17) Pople, J. A.; Segal, G. A. J. Chem. Phys. 1966, 44, 3289. (18) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem. 1990, 94, 8897. (19) Gasteiger, J.; Marsili, M. Tetrahedron 1980, 36, 3219. (20) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P. J. Comput. Chem. 2003, 24, 1999. (21) Becke, A. D. J. Chem. Phys. 1988, 88, 2547. (22) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 786. (23) Dalton, J. B.; Schmidt, C. L. A. J. Biol. Chem. 1933, 103, 549. (24) Svec, H. J.; Clyde, D. D. J. Chem. Eng. Data 1965, 10, 151. (25) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981; p 331. (26) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (27) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (28) Jorgensen, W. L.; Madura, J. D. Mol. Phys. 1985, 56, 1381. (29) Wu, Y.; Tepper, H. L.; Voth, G. A. J. Chem. Phys. 2006, 124, 024503. (30) Van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701. (31) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (32) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (33) Li, L.; Rodriguez-Hornedo, N. J. Cryst. Growth 1992, 121, 33. (34) Piana, S.; Gale, J. D. J. Am. Chem. Soc. 2005, 127, 1975. (35) Boek, E. S.; Briels, W. J.; van Eerden, J.; Feil, D. J. Chem. Phys. 1992, 96, 7010. (36) Boek, E. S.; Briels, W. J.; Feil, D. J. Phys. Chem. 1994, 98, 1674. (37) Kuleshova, L. N.; Zorkii, P. M. Acta Crystallogr. 1981, B37, 1363. (38) Kameda, Y.; Ebata, H.; Usuki, T.; Uemura, O.; Misawa, M. Bull. Chem. Soc. Jpn. 1994, 67, 3159. (39) Sasaki, M.; Kameda, Y.; Yaegashi, M.; Usuki, T. Bull. Chem. Soc. Jpn. 2003, 76, 2293. (40) Alagona, G.; Ghio, C.; Kollman, P. A. J. Mol. Struct.: THEOCHEM 1988, 166, 385. (41) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650. (42) Krynicki, K.; Green, C. D.; Sawyer, D. W. Faraday Discuss. Chem. Soc. 1978, 66, 199. (43) Ma, Y. G.; Zhu, C. Y.; Ma, P. S.; Yu, K. T. J. Chem. Eng. Data 2005, 50, 1192. (44) Docherty, R. Crystal Growth of Organic Materials; In Conference Proceedings Series; Myerson, A. S., Green, D. A., Eds.; American Chemical Society: Washington, DC, 1995. (45) Boek, E. S.; Feil, D.; Briels, W. J.; Bennema, P. J. Cryst. Growth 1991, 114, 389. (46) No, K. T.; Cho, K. H.; Kwon, O. Y.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1994, 98, 10742. (47) Hagler, A. T.; Huler, E.; Lifson, S. J. Am. Chem. Soc. 1974, 96, 5319. (48) Voogd, J.; Derissen, J. L.; Van Duijneveldt, F. B. J. Am. Chem. Soc. 1981, 103, 7701. (49) Bisker-Leib, V.; Doherty, M. F. Cryst. Growth Des. 2003, 3, 221.
JP802949U
