Anisotropy in Stable Conformations of Hydroxylate Ions between the {001} and {110} Planes of TiO2 Rutile Crystals for Glycolate, Lactate, and 2-Hydroxybutyrate Ions Studied by Metadynamics Method.

Control over TiO2 rutile crystal growth and morphology using additives is essential for the development of functional materials. Computer simulation studies on the thermodynamically stable conformations of additives at the surfaces of rutile crystals contribute to understanding the mechanisms underlying this control. In this study, a metadynamics method was combined with molecular dynamics simulations to investigate the thermodynamically stable conformations of glycolate, lactate, and 2-hydroxybutyrate ions at the {001} and {110} planes of rutile crystals. Two simple atom-atom distances were selected as collective variables for the metadynamics method. At the {001} plane, a conformation in which the COO- group was oriented toward the surface was found to be the most stable for the lactate and 2-hydroxybutyrate ions, whereas a conformation in which the COO- group was oriented toward water was the most stable for the glycolate ion. At the {110} plane, a conformation in which the COO- group was oriented toward the surface was the most stable for all three hydroxylate ions, and a second most stable conformation was also observed for the lactate ion at positions close to the {110} plane. For all three hydroxylate ions (α-hydroxycarboxylate ions), the stability of the most stable conformation was higher for the {110} plane than for the {001} plane. At both planes, the stability of the most stable conformation was highest for the 2-hydroxybutyrate ion and lowest for the glycolate ion. Supposing that all three hydroxylate ions serve to decrease the surface free energy at the rutile surface and that a more stable conformation at the rutile surface leads to a greater decrease in the surface free energy, the present results partially explain experimentally observed differences in the changes in growth rate and morphology of rutile crystals in the presence of glycolic, lactic, and 2-hydroxybutyric acids.

Introduction

Control over crystal growth and morphology is essential for the development of functional materials because the functionality of crystals depends on their size and shape. The use of additives is a promising method for achieving such control, and a diverse range of additives have been applied, such as peptides, amino acids, alcohols, and metal ions. In general, the binding of additives to a crystal surface changes the growth rate at that surface, and anisotropy in the binding of additives between different crystal planes changes the morphology.[1] Therefore, when considering the use of a particular additive to control crystal growth or morphology, it is important to investigate the thermodynamically stable conformations of the additives at various crystal surfaces. Although certain experimental techniques, such as infrared spectroscopy[2] and solid-state NMR,[3] can reveal information about the conformations of additives at wide crystal surfaces, the detailed conformations and dynamics of additives at specific planes of crystals grown in solution are difficult to determine experimentally.

Computer simulations, such as molecular dynamics (MD), are helpful tools for investigating the stable conformations of additives at crystal surfaces. To date, MD simulations have been conducted to examine the stable conformations of additives at the surfaces of several inorganic crystals, such as CaCO3[4−8] and TiO2.[1,9−25] Several previous MD simulation studies of these inorganic crystals have discussed the thermodynamic stabilities of a few selected conformations by calculating the potential of mean force.[16,17,20]

The run of an MD simulation is typically on the order of a microsecond or less. Such a short run time is not always sufficient to reproduce all of the possible conformations of an additive at a crystal surface, especially for large additives possessing a complex structure. Ideally, the thermodynamically stable conformations of an additive at a crystal surface should be determined definitively by comparing the free energies of all possible conformations. In addition, to elucidate the effects of additives on crystal growth and morphology, their conformations on all possible exposed facets should be considered. It is impossible to estimate the free energies of all possible conformations on exposed facets using an MD simulation with a short run. Calculating the potential of mean force for all possible conformations on each facet is extremely time-consuming, if the number of possible conformations is large.

Metadynamics (MTD)[26] is a method for enhancing the transition of a system between different states in an MD simulation by increasing the probability of reaching high-energy states. The MTD method affords a free-energy landscape that represents the thermodynamic stabilities of all possible states of a system. Thus, if the MTD method is combined with an MD simulation of an additive at a crystal surface, the thermodynamic stabilities of all possible conformations can be determined definitively, even if the run of the simulation is not very long. The applicability of the MTD method to investigating the stable conformation of an additive at a crystal plane has been confirmed in recent simulation studies.[25,27] However, to the best of our knowledge, the MTD method has not yet been applied to investigate the anisotropy in the stable conformation of an additive between different crystal planes to infer the mechanism of crystal growth or morphology control.

In this study, we used the MTD method to examine the thermodynamically stable conformations of hydroxylate ions at the {001} and {110} planes of a TiO2 rutile crystal. TiO2 is a highly functional material,[28−33] and the functionality of TiO2 crystals depends on their size and shape.[34−46] Rutile crystals are the most thermodynamically stable polymorph of TiO2 and have been widely employed in practical applications.[47−52] Thus, the use of additives to control the growth and morphology of rutile crystals is an important area of research for the development of functional materials.

Kobayashi et al. reported that the growth of rutile crystals in the presence of hydroxy acids (α-hydroxycarboxy acids), such as glycolic acid,[47−49] lactic acid,[52] and 2-hydroxybutyric acid,[52] resulted in a higher growth rate at the rutile crystal surface and an increase in the {110} area ratio compared with that for the equilibrium shape, which is a tetragonal prism bounded by {110} planes and terminated by a pair of tetragonal pyramids bounded by {101} planes.[53−55] The authors speculated that the increased growth rate originated from the binding of hydroxylate ions (α-hydroxycarboxylate ions) at the rutile crystal surface and the increased {110} area ratio originated from the anisotropic binding of the ions between different crystal planes.[1,47,55] Since the increase in the {110} area ratio was caused by an increase in the growth rate at the {001} planes of rutile crystals,[50] we previously performed MD simulations of a glycolate ion at the {001} and {110} planes of a rutile crystal.[24] The results suggested a difference in the stable conformation between the two planes, which was considered to be related to the experimental results reported by Kobayashi et al. Detailed knowledge regarding the difference in the stable conformations between different planes of a rutile crystal for various hydroxylate ions obtained using the MTD method would be expected to significantly advance our understanding of the mechanisms of crystal growth and morphology control of actual rutile crystals in the presence of hydroxy acids.

The purpose of this study was to evaluate the applicability of the MTD method to studying the mechanisms of crystal growth and morphology control using additives. In this paper, we demonstrate that the difference in the stable conformations of glycolate, lactate, and 2-hydroxybutyrate ions between the {001} and {110} planes, which were analyzed using the MTD method, partially explains the experimentally observed differences in the growth rate and morphology of rutile crystals grown in the presence of hydroxy acids, namely, glycolic, lactic, and 2-hydroxybutyric acids.

It has been debated whether water molecules on the surfaces of rutile crystals are dissociated.[56−65] Recently, Wang et al. indicated that an isolated water molecule on the rutile {110} plane was dissociated.[58] Kang et al. reported that the dissociation of water molecules at the rutile crystal surface can change the adsorption conformation of a protein molecule.[59] In contrast, González et al. demonstrated that if water molecules on the rutile surface formed a molecular layer, the layer structure in which no water molecules are dissociated was energetically stable for the {110} and {011} planes.[57] In this study, we used a rigid water molecule model by assuming that water molecules at the rutile surface were not dissociated. Simulations in which the dissociation of water molecules is allowed, such as first-principles MD simulations, are needed to confirm this assumption. However, the purpose of the present simulations was qualitative elucidation of the differences in the stable conformations of hydroxylate ions between different rutile planes. In this paper, we demonstrate that, as with earlier simulation studies that used a rigid water molecule model,[12,14,18,24] the present simulation results are qualitatively in good agreement with previous experimental investigations.[47−49,52]

Simulation Method

Potential Models

The potential models used in this study were the same as those used in our previous study.[24] The interactions acting on Ti and O in the rutile crystal were estimated using the potential model proposed by Matsui and Akaogi.[66] This model reproduces Ti–O bond lengths in the bulk and on the surface of the rutile crystal and is therefore suitable for the present study. The potential parameters for the lactate and 2-hydroxybutyrate ions were determined using a general AMBER force field (GAFF).[67] The charge on each atom of the hydroxylate ions was determined using ANTECHAMBER,[68] which is an extensively used module for determining the parameters in the input files needed to execute an MD simulation using a GAFF. Since the potential parameters for the GAFF were optimized using the TIP3P model,[69] the interaction between a pair of water molecules was also estimated using the TIP3P model.

The potential parameters of the glycolate ion were already determined using the GAFF in the previous study.[24] However, in that study, the charge on each atom of the glycolate ion was determined such that the minimum-energy structure generated using the GAFF corresponded to the optimum structure obtained via first-principles calculations without using ANTECHAMBER.[24] In the present study, to impartially evaluate the relative stabilities of the conformations of the three hydroxylate ions at the rutile surface, the simulation of the glycolate ion was conducted using the charge on each atom, which was redetermined using ANTECHAMBER in the present study, and the parameters for bond potentials, valence angle potentials, dihedral angle potentials, and Lennard-Jones potentials, which were determined in our previous study.[24] The minimum-energy structures of the hydroxylate ions, which were obtained via first-principles calculations using the MP2 method with the 6-311+G** basis function, are presented in Figure , and the atomic charge values determined for each hydroxylate ion using ANTECHAMBER are listed in Table S1.

Figure 1

Minimum-energy structures of the glycolate, lactate, and 2-hydroxybutyrate ions.

As in our previous study,[24] the interactions between rutile and the hydroxylate ions were determined via the same method as that used by Carravetta and Monti.[11] Briefly, the interaction energy between a TiO5H6 cluster with the same structure as a rutile crystal, which was used instead of a rutile crystal, and the hydroxylate ion was evaluated for different intermolecular distances via first-principles calculations using the MP2 method with the 6-311+G** basis function. These calculations were performed for two different orientations of each hydroxylate ion. Then, the calculated interaction energies were fitted to a potential function, which was assumed to be a Coulomb potential plus a Buckingham potential. The values of the Buckingham potential parameters for the lactate and 2-hydroxybutyrate ions are listed in Table S2, and those for the glycolate ion were provided in our previous paper.[24]

Simulation Systems

Two simulation systems, namely, a system for the simulation of the {001} plane (the {001} system) and a system for the simulation of the {110} plane (the {110} system), were established as in our previous study.[24] Each system was a rectangular parallelepiped consisting of a rutile crystal containing 2016 TiO2 units, a water phase containing 3000 H2O molecules, and a vapor phase. The rutile crystal was placed at the center of the system such that each of the two crystal planes was perpendicular to the z-axis of the system. A single hydroxylate ion was placed in the water phase, which was juxtaposed to one of the two crystal planes, and the remainder of the system was the vapor phase.

The system dimensions were 38.8 × 45.0 × 140.0 Å3 for the {001} system and 35.4 × 38.8 × 140.0 Å3 for the {110} system. The dimensions in the x and y directions were determined from a separate MD simulation of a bulk rutile crystal at 298 K and 1 atm. The thickness of the rutile crystal in the z direction was 35.4 Å for the {001} system and 45.0 Å for the {110} system and that of the water phase was approximately 50 Å for the {001} system and 64 Å for the {110} system. A schematic illustration of the simulation system and the structures of the rutile surface in each system are depicted in Figure S1. As each hydroxylate ion had a charge of −e, a Na+ ion was added to the vacuum phase to maintain the system in an electrically neutral state. The position of the Na+ ion was fixed at the center of the x–y plane at a distance of 3 Å in the z direction from the crystal plane that was not adjacent to the water phase. In this study, the position of the crystal plane in each system was defined as the position of the outermost O atom layer of the rutile crystal in the z direction.

The surface energy, γ, of each plane was estimated using the following equationHere, Usurf and Ubulk are the potential energies of the system including the rutile crystal only and a bulk rutile crystal comprising the same number of TiO2 units, respectively, and S is the surface area of one side of the rutile crystal in the system. The estimated γ values were 0.148 eV/Å2 for the {001} system and 0.121 eV/Å2 for the {110} system, indicating that the structure of the {110} plane was more stable than that of the {001} plane. The absolute values of γ in the present potential model were greater than the values obtained via density functional theory calculations.[57,65,70] However, we confirmed that the relative stabilities of the structures of the {001} and {110} planes were satisfactorily reproduced by the present potential model.[65,70]

MTD Method

The MTD method is an enhanced sampling method that increases the probability of visiting high free-energy states by adding a bias potential, which is a function of several collective variables (CVs), to the minima of the potential energy surface of the system.[26] Selecting appropriate CVs is crucial for obtaining a desirable free-energy landscape for a system of interest. In this study, two CVs were chosen, namely, the distance between the C1 atom and the surface, d_C1, and that between the C2 atom and the surface, d_C2. These CVs were selected because they permitted the convenient determination of the thermodynamically stable orientation and position of the hydroxylate ion on each plane. In a similar manner to the present study, YazdanYar et al. investigated the binding conformations of amino acids at a rutile {110} plane using two CVs, that is, the distances between pairs of centers of mass for different atomic groups.[25] In this paper, we demonstrate that the present CVs of simple atom–atom distances are sufficient for determining the stable binding conformation of an additive at a crystal surface.

With reference to earlier studies,[25,27,71,72] the well-tempered MTD method[73] was used to obtain the free-energy landscape with sufficiently high accuracy.[73] The bias factor parameter in the well-tempered MTD method was set to 5.

Computations

The Newton’s equations of motion for each atom of the TiO2 units, water molecules, and hydroxylate ion were integrated with a time step of 1 fs using a leapfrog algorithm.[74] The temperature was maintained at 298 K by means of a Berendsen thermostat with a coupling parameter of 0.1 ps.[75] The C–H and O–H distances of each hydroxylate ion were kept constant at their equilibrium values using the SHAKE algorithm.[76] An MD simulation combined with the MTD method was performed using DL_POLY 2.20,[77] in which PLUMED 1.3[78] was implemented.

The long-range Coulomb interactions were estimated using the Ewald method, and the short-range Lennard-Jones and Buckingham interactions were cut off at an interatomic distance of 10 Å. The accuracy of the potential energy and force calculated using the Ewald method was set to 10–6. The real space cutoff distance was set to 10 Å. The Ewald convergence parameter was 0.3208 Å–1. The maximum indices of the reciprocal lattice vector in the x, y, and z directions for the reciprocal space sum were 13, 15, and 42, respectively, for the {001} system and 12, 13, and 45, respectively, for the {110} system.

The bias potential for the MTD method was represented as a function of Gaussians. The height of the Gaussians was set to 0.03 kBT, and the sigma parameter of the Gaussians was set to 0.35 Å. The bias potential was added to the Hamiltonian of the system at every 100 MD steps. The total run of the MD simulation was 12 ns, which was sufficiently long to determine the free-energy landscape because the CVs frequently changed within their entire ranges during the run (Figure S2). This run of 12 ns was considerably shorter than the total run of the MD simulation in the previous study (30 ns).[24] As we focused only on the stable binding conformations at the rutile surface, the variable range of each CV was restricted to the region near the surface by setting the directive UWALL in PLUMED 1.3 to CV = 8 Å. It is known that the Buckingham potential function generates a spurious maximum at a short interatomic distance and sharply decreases to negative potential energies as the interatomic distance approaches zero.[79] Therefore, if a particular atom of the hydroxylate ion became too close to the rutile surface when the MTD method searched for a high free-energy state, the force acting on a particular atom of the hydroxylate ion may be extremely large, causing computation errors, such as the convergence error of the SHAKE algorithm.[77] To avoid such errors, we also set the directive LWALL to CV = 3 Å. These directives produced the following external wall potential, Vwall, acting on the systemHere, Δs was the distance between the positions s (d_C1 ≥ 8 Å and d_C2 ≥ 8 Å for UWALL and d_C1 ≤ 3 Å and d_C2 ≤ 3 Å for LWALL) and slimit on the d_C1–d_C2 plane. slimit was (8 Å, 8 Å) for UWALL and (3 Å, 3 Å) for LWALL. The energy constant κ and rescaling factor ε were set to 1 kJ/mol and 1 Å, respectively. p was set to 4, such that Vwall gradually increased with increasing |s – slimit|.

Notably, the free energy at each CV, which was computed in this study, was the difference from the free energy at the upper level of the deposited bias potentials in the space of the CVs, ΔF. Thus, the relative free energy between different CVs is meaningful rather than the absolute value of the free energy at each CV. We confirmed the reliability of the relative free energy between different CVs by comparing the free-energy landscapes, which were created using the simulation data for various periods in the simulation (Figure S3). In this study, we assumed that the upper level of the deposited bias potentials was equivalent for the simulations of both the {001} and {110} systems.

Simulation Results

Glycolate Ion

First, we examined the binding conformation of the glycolate ion to confirm that the present method properly provides the stable conformation of a hydroxylate ion at a rutile surface. Figure shows the free-energy landscapes along the d_C1 and d_C2 axes for the glycolate ion conformation. The free-energy minima appeared at (d_C1, d_C2) = (5.8 Å, 4.7 Å) for the {001} plane and (3.9 Å, 5.0 Å) for the {110} plane. These results indicate that the conformation in which the COO– group was oriented toward water corresponded to the free-energy minimum for the {001} plane, whereas the conformation in which the COO– group was oriented toward the rutile surface corresponded to the free-energy minimum for the {110} plane. These most stable conformations are the same as those indicated by our previous MD simulation study,[24] thus confirming the applicability of the present method to studying the stable conformation of a hydroxylate ion at a rutile surface. The lower value of the free-energy minimum for the {110} plane (−51.02 kJ/mol) than for the {001} plane (−50.11 kJ/mol) indicates that the stability of the most stable conformation was greater for the {110} plane than for the {001} plane. The lower free-energy minimum for the {110} plane than for the {001} plane is also consistent with the lower potential energy between the glycolate ions of the rutile crystal for the {110} plane than for the {001} plane, which was indicated by our previous MD simulation study.[24]

Figure 2

Free-energy landscapes for the conformation of the glycolate ion at the {001} and {110} planes. The left- and right-hand panels in the upper part of the figure show snapshots of a typical conformation corresponding to the free-energy minimum, (d_C1, d_C2) = (5.8 Å, 4.7 Å) for the {001} plane and (3.9 Å, 5.0 Å) for the {110} plane, with and without the hydrogen-bonded networks of water molecules (white rods), respectively. The yellow, red, dark gray, and white spheres represent the Ti, O, C, and H atoms, respectively. Water molecules are not shown in the snapshots.

The free-energy landscape for the {001} plane also indicated a second free-energy minimum at (4.3 Å, 4.6 Å) corresponding to the second most stable conformation, in which the COO– group was oriented toward the rutile surface. This conformation was not observed in our previous MD simulation.[24] The free energy at the ridge between the minimum and second minimum (−48.39 kJ/mol) was only 1.72 kJ/mol higher than that at the minimum and only 1.32 kJ/mol higher than that at the second minimum. This result implies that transitions between the most stable and second most stable conformations may occur frequently at the {001} plane.

Figure presents the number density (ρ) profiles of the O (O in TiO2), Ow (O of water molecule), Oc (carboxyl O of the hydroxylate ion), Oh (hydroxyl O of the hydroxylate ion), and Ho (hydroxyl H of the hydroxylate ion) atoms of the glycolate ion along the z direction. The distinct peaks appearing in the ρ profiles of the Oc, Oh, and Ho atoms originated from the stable conformation of the glycolate ion at the rutile plane. It can be seen that the glycolate ion was separated from the surface by a layered structure of water present on the surface. This result is consistent with our previous MD simulation, in which indirect binding of the glycolate ion to the rutile plane led to greater stability than direct binding.[24]

Figure 3

Number density (ρ) profiles of the O, Ow, Oc, Oh, and Ho atoms of the glycolate ion along the z direction. The positions of the first and second layers of water molecules for the {001} plane and the first, second, and third layers of water molecules for the {110} plane are indicated by red arrows. The origin of the z-axis is the outermost layer of O atoms in the rutile crystal.

The layered structures of water at each rutile plane, which were obtained via MD simulations of each system without including the hydroxylate ion and Na+ ion,[24] are depicted in Figure S4. The layered structure of water can also be observed in the snapshot shown in Figure . The Oc or Oh atom of the glycolate ion in the stable conformation was located at the positions of the Ow atoms in the second layer of water molecules at z = 4 Å for the {001} plane and in the third layer of water molecules at z = 3.4 Å for the {110} plane. The longer distance between the outermost O atoms at the surface and the position at which the free-energy minimum appeared for the {001} plane than for the {110} plane originated from the longer distance between the surface and the second layer of water molecules for the {001} plane than between the surface and the third layer of water molecules for the {110} plane.

Notably, the binding conformation in which the hydroxylate ion was directly bound to the surface could not occur owing to Vwall. Ideally, the simulation should be performed without using the directive UWALL. In this study, we assumed that the conformation in which the hydroxylate ion was directly bound to the surface was not stable, as was indicated by our previous MD simulation study of the glycolate ion.[24]

Lactate Ion

Figure shows the free-energy landscapes for the lactate ion. The free-energy minima appeared at (4.6 Å, 5.8 Å) for the {001} plane and (4.0 Å, 5.0 Å) for the {110} plane. As in the case of the glycolate ion, the lactate ion was separated from the surface by a layered structure of water molecules and the stability of the most stable conformation was greater for the {110} plane than for the {001} plane. At both planes, the values of the free-energy minima for the lactate ion were lower than those for the glycolate ion, suggesting that the stability of the most stable conformation was greater for the lactate ion than for the glycolate ion.

Figure 4

Free-energy landscapes for the conformation of the lactate ion at the {001} and {110} planes. Snapshots of a typical conformation corresponding to the free-energy minimum, (4.6 Å, 5.8 Å) for the {001} plane and (4.0 Å, 5.0 Å) for the {110} plane, are also shown.

Notably, in contrast to the case of the glycolate ion, the conformation in which the COO– group was oriented toward the surface was the most stable for both the {001} and {110} planes. Furthermore, the free-energy landscape for the {110} plane indicates the existence of two metastable conformations at (5.3 Å, 4.3 Å) and (2.7 Å, 2.7 Å). In particular, the thermodynamic stability of the latter, which corresponded to the second most stable conformation, was almost equivalent to that of the most stable conformation. Notably, the free-energy landscape for the {110} plane in Figure shows the free energies from which Vwall, which was estimated with eq , was subtracted.

To determine the relative stabilities of these two conformations, we increased the simulation run for the lactate ion at the {110} plane to 16 ns and created the free-energy landscape using the data from this longer simulation. The free-energy landscape indicated the same relative stabilities of the two conformations as observed in Figure (Figure S5). However, the free-energy difference remained very small, even in the free-energy landscape obtained from the longer simulation. Although we consider that the relative stabilities of the different conformations shown in Figure are correct, more comprehensive studies are needed to confirm this.

Figure shows the ρ profiles for the lactate ion. As in the case of the glycolate ion, the Oc or Oh atom of the most stable conformation of the lactate ion at the {001} plane was located at the positions of the Ow atoms in the second layer of water molecules. For the {110} plane, the positions of the peaks that appeared in the ρ profiles of the Oc atoms were almost the same as those observed for the glycolate ion. However, in contrast to the case of the glycolate ion, the ρ profile of the Oh atom contained a sharp peak at z = 2.3 Å, where the Ow atom in the second layer of water molecules was located. The appearance of this sharp peak originated from the second most stable conformation. This result suggests that the lactate ion had a tendency to approach the {110} plane more closely than the glycolate ion.

Figure 5

ρ profiles of the O, Ow, Oc, Oh, and Ho atoms of the lactate ion along the z direction.

2-Hydroxybutyrate Ion

Figure shows the free-energy landscapes for the 2-hydroxybutyrate ion. The free-energy minima appeared at (4.6 Å, 5.8 Å) for the {001} plane and (4.0 Å, 5.0 Å) for the {110} plane, indicating that the most stable conformations were the same as those of the lactate ion. As in the case of the glycolate and lactate ions, the stability of the most stable conformation was greater for the {110} plane than for the {001} plane. For both planes, the values of the free-energy minima for the 2-hydroxybutyrate ion were lower than those for the glycolate and lactate ions, suggesting that the most stable conformation of the 2-hydroxybutyrate ion had the greatest stability among the most stable conformations of the three hydroxylate ions.

Figure 6

Free-energy landscapes for the conformation of the 2-hydroxybutyrate ion at the {001} and {110} planes. Snapshots of a typical conformation corresponding to the free-energy minimum, (4.6 Å, 5.8 Å) for the {001} plane and (4.0 Å, 5.0 Å) for the {110} plane, are also shown.

In contrast to the case of the lactate ion, the 2-hydroxybutyrate ion did not show a tendency to approach the {110} plane more closely than the glycolate ion. Figure shows the ρ profiles for the 2-hydroxybutyrate ion. It can be seen that the conformation of the 2-hydroxybutyrate ion was strongly influenced by the layered structure of water molecules, as with the stable conformations of the glycolate and lactate ions.

Figure 7

ρ profiles of the O, Ow, Oc, Oh, and Ho atoms of the 2-hydroxybutyrate ion along the z direction.

Discussion

Differences in Stable Binding Conformation between Hydroxylate Ions

The present results indicate that the most stable conformation at the {001} plane for the lactate and 2-hydroxybutyrate ions is quite different from that for the glycolate ion. Here, we discuss the reason for this difference. In our previous paper,[24] we rationalized why the COO– group of the stable conformation of the glycolate ion at the {001} plane was oriented toward water as follows. The COO– group interacts strongly with both the rutile surface and water. However, as the glycolate ion was located at a position relatively far from the surface, the COO– group could not interact strongly with the surface. Therefore, the conformation in which the COO– group was oriented toward water to form hydrogen bonds with water molecules was preferred. The present results suggest that this explanation does not apply to the lactate and 2-hydroxybutyrate ions.

If the lactate and 2-hydroxybutyrate ions adopted a conformation in which the COO– group was oriented toward water at the {001} plane, the methyl group of the lactate ion and the ethyl group of the 2-hydroxybutyrate ion would necessarily be oriented toward the surface. Consequently, the water molecules surrounding the methyl or ethyl group would form a hydration structure, resulting in a loss of entropy and disruption of the layered structure of water formed at the {001} plane. Thus, the conformation in which the COO– group was oriented toward water would lead to an increase in the free energy of the system. This hypothesis qualitatively explains why the conformation in which the COO– group was oriented toward the surface was rather preferred for the lactate and 2-hydroxybutyrate ions.

To test this hypothesis, we performed an MD simulation of the {001} system in which the conformation of the lactate ion was fixed such that its COO– group was oriented toward water and its methyl group was oriented toward the surface (Figure a). In this simulation, the C3 atom constituting the methyl group was located at a position close to the second layer of water molecules (Figure b). Figure c shows a plot of the pair distribution function, g, as a function of the distance between the C3 and Ow atoms, rC3–Ow. For comparison, a plot of the values of g obtained from a separate MD simulation in which the methyl group of the lactate ion was oriented toward water is also shown. The appearance of the large peak around rC3–Ow = 3.7 Å provides evidence for the formation of a hydration structure around the methyl group. Notably, the water molecules in parts of the hydration structure were those contained within the second layer of water molecules. A small peak also appeared around rC3–Ow = 3.3 Å because the C3 atom displaced some of the water molecules in the second layer of water molecules from their ideal positions, and these water molecules were also part of the hydration structure. These results suggest that the layered structure of water was partially disrupted by the formation of the hydration structure. The higher peak observed for the second layer of water molecules in Figure b compared with that in Figure originated from a loss in entropy of the water molecule arrangement due to the formation of the hydration structure.

Figure 8

MD simulation of the {001} system in which the conformation of the lactate ion was fixed such that its COO– group was oriented toward water and its methyl group was oriented toward the surface. (a) Snapshot of the lactate ion in the system. (b) ρ profiles of the O, Ow, Hw, Oc, and C3 atoms of the lactate ion along the z direction. (c) Pair distribution function, g, as a function of the distance between the C3 and Ow atoms, rC3–Ow (solid line). For comparison, the values of g obtained from a separate MD simulation in which the methyl group of the lactate ion was oriented toward water are also plotted (dashed line).

Notably, the stable conformations observed for the lactate and 2-hydroxybutyrate ions at both planes and the glycolate ion at the {110} plane were consistent with a well-known stable conformation for organic molecules containing a COO– group at ionic crystal surfaces, in which the COO– group is oriented toward the surface.[80−82] However, the present result for the glycolate ion at the {001} plane suggests that the conformation in which the COO– group is oriented toward water can also be the most stable, depending on the structure of the additive, the thickness and stability of the layered structure of water, and the interaction of the additive with the surface.

Relationship between the Present Results and Experimental Rutile Crystals Grown in the Presence of Hydroxy Acids

According to the experimental study performed by Kobayashi et al.,[52] glycolic, lactic, and 2-hydroxybutyric acids, all induce an increase in the growth rate at the rutile surface. The extent of the increase in the growth rate was reported to be considerably greater for lactic and 2-hydroxybutyric acids than for glycolic acid. Supposing that a more stable conformation of hydroxylate ions at the rutile surface leads to a greater increase in the growth rate, the present results indicating the greater stabilities of the most stable conformations of the lactate and 2-hydroxybutyrate ions compared with the glycolate ion explain the differences in the experimentally observed growth rates in the presence of the three hydroxy acids.

Typically, the growth rate at a crystal surface decreases if additives bind to lattice sites, such as kinks, and inhibit the incorporation of crystal atoms at these sites.[1] However, we do not believe that this applies to the present systems because the stable conformations of all three hydroxylate ions appeared at positions relatively far from the surface owing to the layered structure of water molecules. Conversely, the growth rate at a crystal surface may increase if the additives serve to decrease the surface free energy.[1] We therefore speculate that the increase in the growth rate at the rutile surface originates from a decrease in the surface free energy. It is natural to speculate that a more stable conformation at the rutile surface leads to a greater decrease in the surface free energy and, hence, a greater increase in the growth rate.

The experimental study also indicated that the extent of the increase in the growth rate at the {110} plane was slightly greater for lactic acid than for 2-hydroxybutyric acid.[52] This cannot be explained by the present results, in which the most stable conformation at the {110} plane was more stable for the 2-hydroxybutyrate ion than for the lactate ion. This experimental observation may be attributable to the existence of the second most stable conformation at positions close to the {110} plane only for the lactate ion. In general, the decrease in the surface free energy mediated by additives is enhanced by increasing the additive concentration.[1] We speculate that the existence of both the most stable and second most stable conformations at the {110} plane only for the lactate ion leads to a higher surface concentration for the lactate ion than for the 2-hydroxybutyrate ion.

The experimental study also indicated that glycolic, lactic, and 2-hydroxybutyric acids all induce an increase in the {110} area ratio of rutile crystals.[52] This increase in the {110} area ratio originates from an increase in the growth rate at the {001} plane only for glycolic acid and an increase in the growth rate at both the {001} and {110} planes for lactic and 2-hydroxybutyric acids.[52] Thus, it is speculated that the anisotropy in the stable conformation of hydroxylate ions is not identical for glycolate, lactate, and 2-hydroxybutyrate ions. This speculation is consistent with the present finding that the anisotropy in the stable conformation of the glycolate ion between the {001} and {110} planes was different from those for the lactate and 2-hydroxybutyrate ions.

The experimentally observed increase in the {110} area ratio suggests that the extent of the decrease in the surface free energy mediated by hydroxylate ions is greater for the {001} plane than for the {110} plane. The origin of this difference remains unclear. The present simulation was performed using a simple model system in which a single hydroxylate ion was placed at an ideally flat rutile surface, whereas the growth of actual rutile crystals occurs at surfaces possessing steps and kinks in the presence of numerous hydroxylate ions. Further simulations using a more realistic system are therefore required to elucidate the origin of the differences in the extent of the decrease in the surface free energy mediated by hydroxylate ions at the various planes.

Conclusions

The thermodynamically stable conformations of glycolate, lactate, and 2-hydroxybutyrate ions at the {001} and {110} planes of a rutile crystal were investigated using the MTD method. For both planes, the stability of the most stable conformation was highest for the 2-hydroxybutyrate ion and lowest for the glycolate ion. For all three hydroxylate ions, the stability of the most stable conformation was higher at the {110} plane than at the {001} plane.

In this study, we assumed that the conformation in which the hydroxylate ion was directly bound to the surface was not stable, in accordance with our previous MD simulation study of the glycolate ion.[24] However, more comprehensive studies should be conducted to confirm this in the future.

The anisotropy in the most stable conformation between the planes was not identical for the three hydroxylate ions. For the {110} plane, all three hydroxylate ions adopted the same most stable conformation in which the COO– group was oriented toward the surface. However, for the {001} plane, only the lactate and 2-hydroxybutyrate ions adopted a most stable conformation in which the COO– group was oriented toward the surface, whereas the glycolate ion adopted a most stable conformation in which the COO– group was oriented toward water. For the {110} plane, a second most stable conformation was also found for the lactate ion at positions close to the surface.

Supposing that all three hydroxylate ions serve to decrease the surface free energy at the rutile surface and that a more stable conformation leads to a greater decrease in the surface free energy, the present results qualitatively explain the experimentally observed differences in the extent of the increase in the growth rate at the surfaces of rutile crystals in the presence of glycolic, lactic, and 2-hydroxybutyric acids.

In this study, two simple atom–atom distances, namely, that between a surface Ti ion and the C1 atom of the hydroxylate ion and that between the Ti ion and the C2 atom, were used as the CVs. The results demonstrate that the MTD method with these simple CVs is sufficient for determining the most stable and metastable conformations of hydroxylate ions at rutile planes and evaluating their relative stabilities.

To date, the differences in the stabilities of the conformations of various additives at crystal surfaces have mainly been evaluated based on the differences in the strengths of the electrostatic interactions between the additive and crystal surface, which originate from the differences in the valences of the additives or the functional groups present.[83] In this study, the valence was identical for the three hydroxylate ions. The functional groups present in the hydroxylate ion were also the same, namely, a single COO– group and a single OH group. Therefore, the differences in the strengths of the electrostatic interactions between the three hydroxylate ion additives and the surface can be considered to be small. The differences in the anisotropic stable conformations of the three hydroxylate ions, which were successfully detected using the present MTD method, originated mainly from the very small structural differences between the hydroxylate ions, namely, H atom, methyl group, or ethyl group.

In conclusion, the present MTD method is applicable to investigating fine differences in the stable conformations of additives between different crystal surfaces, which is expected to improve our understanding of the mechanisms of crystal growth and morphology control using additives.

Elucidation of the stable conformations of large additives possessing complex structures, such as polypeptides, at inorganic crystal surfaces is becoming increasingly important for the development of hybrid materials[84] and understanding the mechanism of biomineralization.[4,85,86] The present MTD method is expected to serve as a helpful tool for studying the stable conformations of such large and complex additives. This method should also assist in studying ice growth inhibition by antifreeze proteins.[87−89] Kobayashi et al. reported that the changes in the growth rate at the rutile surface and the {110} area ratio in the presence of carboxylic acids or alcohols were different from those in the presence of hydroxy acids.[52] The application of the present MTD method to elucidate this difference will also be an interesting topic for future work.