A New Transferable Forcefield for Simulating the Mechanics of CaCO3

ARTICLE pubs.acs.org/JPCC

A New Transferable Forcefield for Simulating the Mechanics of CaCO3 Crystals Shijun Xiao,† Scott A. Edwards,† and Frauke Gr€ater*,†,‡ †

CAS-MPG Partner Institute and Key Laboratory for Computational Biology, Chinese Academy of Sciences, 320 Yueyang Road, Shanghai, China ‡ Heidelberg Institute for Theoretical Studies, Schloss-Wolfsbrunnenweg 35, Heidelberg, Germany

bS Supporting Information ABSTRACT:

Many sets of forcefield parameters for calcium carbonate (CaCO3) and CaCO3water interactions have been developed for thermodynamic calculations, but growing interest in modeling the molecular-scale mechanics of biomineral nanocomposite materials such as nacre has led to a need for interaction parameters that accurately model the anisotropic mechanical properties of CaCO3. A novel forcefield for aragonite, one polymorph of CaCO3, has been fitted to the structure and elastic constants of the mineral, and the validation of these interaction parameters demonstrates that the forcefield can well capture the shear and elastic moduli of aragonite and also performs well when transferred to other CaCO3 polymorphs. The corresponding aragonitewater and aragoniteprotein parameters are also obtained and utilized in force probe molecular dynamics (FPMD) simulations of the forced desorption of an acidic polypeptide from an aragonite crystal surface, resulting in a rupture force of roughly 60 pN per amino acid residue at pulling speeds characteristic of Atomic Force Microscope experiments. Our forcefield for CaCO3 and CaCO3protein interactions can be applied to study the physical and mechanical properties of organicinorganic composite systems, especially for the next generation of bionanocomposite materials.

’ INTRODUCTION Spiders and silkworms secrete tough silk fibers to form cocoons and webs for protection and hunting,1 invertebrates produce stiff shells to resist pressure and harm from the environment,2 and many species, including humans, can grow hard tissues, such as bones and teeth,3 for supporting the body and breaking up food. All of these mechanically enhanced biomaterials are composites with unique and complex hierarchical structures. In the past half-century, biomaterials have come under ever-growing attention from researchers in several fields, not only because of the inherent interest of explaining the molecular mechanisms of these biosynthesis processes but also because of potential industrial applications. Among these biomaterials, nacre from the mollusk shell is perhaps one of the most thoroughly studied.48 However, a detailed understanding of the mechanisms for nacre’s fracture resistance and stress dissipation still remains elusive. Ninety five percent by weight of nacre is calcium carbonate (CaCO3), which is hard but brittle, and the r 2011 American Chemical Society

remaining components are organic, mainly silk-like proteins and chitin. Nacre overcomes the apparent disadvantages of its constituents by combining the strengths of both CaCO3 and soft matter to become stiff and tough. The elastic modulus of nacre is about 70 (dry) to 90 (wet) GPa, and the fracture work is more than 3000 times that of pure CaCO3.9 Nacre’s excellent mechanical performance is attributed to its lamellar “mortar and brick” structure. Unearthing the dominant factors underlying its outstanding mechanical properties will greatly enhance our capacity for designing industrial analogues. The molecular dynamics (MD) simulation technique was first developed to deal with problems in thermodynamics and physical chemistry but is becoming more and more a central tool in computational mechanics, as it can illuminate processes at the Received: March 24, 2011 Revised: August 29, 2011 Published: September 02, 2011 20067

dx.doi.org/10.1021/jp202743v | J. Phys. Chem. C 2011, 115, 20067–20075

The Journal of Physical Chemistry C

ARTICLE

developed mainly to reproduce structural and thermodynamic properties, this force field also has limited capabilities in capturing the main elastic properties of CaCO3 (as discussed in Results). An accurate representation of the anisotropic elastic properties is needed when manipulating 3D models in MD simulations. Hence, there is an outstanding need to develop a set of forcefield parameters for CaCO3 bulk crystal, including cross-term parameters between CaCO3, water, and protein, that is suitable for accurate simulation of the response of carbonate biomaterials to mechanical perturbation in three dimensions. While CaCO3 has three polymorphs, calcite, aragonite, and vaterite, we focused on aragonite because it is the CaCO3 polymorph found in nacre. We made use of aragonite’s crystal structure and elastic properties during the optimization of forcefield parameters for intracrystal interactions. The cross-term parameters between aragonite and water were generated by fitting to the unit-cell structure parameters of monohydrocalcite (CaCO 3 3 H 2O). Atom-wise parameters, deduced from our aragonitewater forcefield, were then combined with peptide parameters from the OPLS-AA forcefield24 to generate aragonite protein van der Waals interactions. After the forcefield parametrization, additional simulations were performed to verify the fitted forcefield. Finally, the new forcefield was used to simulate the forced unbinding of a polypeptide from an aragonite crystal surface in aqueous solution. Our set of parameters represents the first forcefield developed specifically for accurately modeling the mechanical response of aragonite in an aqueous and/or protein environment.

’ METHODS Figure 1. (A) Overview of the simulation systems of aragonite and aragonite tablet in water for FPMD simulations, and (B) elastic modulus for aragonite and calcite with the forcefield in Table 1 from GULP and Gromacs, compared with previous developed forcefield and experiments.

atomic scale that are inaccessible to direct study by experiments or continuum methods such as finite element analysis (FEA).10 In the past decade, force probe molecular dynamics (FPMD) and force clamp molecular dynamics (FCMD) have been developed to study the direct mechanical perturbation of molecular systems, in parallel with complementary single-molecule experimental techniques such as atomic force microscopy (AFM).1113 Recently, the force distribution analysis (FDA) method, which reveals internal strain distributions in complex molecular systems, was introduced by our research group and successfully applied to several biological protein molecules.1416 However, few MD studies of the combined organic and mineral components of nacre have been reported so far, partly because of the lack of reliable CaCO3water and CaCO3protein forcefield parameters. Although bulk CaCO3 has been well studied via both experiments and theoretical calculations,1723 most existing forcefields do not provide a satisfactory description of the elastic properties of bulk CaCO3 crystal. One of the most widely used is the forcefield of Pavese et al.,17 fitted to elastic and vibrational data, which leads to a stable unit cell structure and elastic constant matrix values that agree well with experiment in Table 3; however, agreement is less satisfactory for the elastic moduli along the three unit-cell axes (Figure 1B), as deduced by extended Hook’s law from the elastic constant matrix. The latest forcefield of CaCO3 published by Raiteri et al.21 has facilitated successful atomic-resolution simulations of the growth of amorphous calcium carbonate; however, because the forcefield was

Forcefield Development. Mineral Forcefield. Calcium carbonate shares the characteristics of both ionic and molecular crystals, due to the strong covalent bonds within the carbonate ion and the high positive charge of the calcium ion. Existing forcefields for CaCO3 treat interatomic interactions as a combination of a purely attractive Coulomb potential and a van der Waals potential that is attractive at large separations and repulsive when atoms overlap.17,20 For calculations in lattice systems, such as minerals and metals, the van der Waals interaction is typically treated with a Buckingham potential25 (eq 1), in contrast to protein forcefields which mostly rely on LennardJones (LJ) potentials26 (eq 2).

r C VBH ¼ A expð Þ 6 F r

ð1Þ

C12 C6 r 12 r 6

ð2Þ

VLJ ¼

In eqs 1 and 2, r represents the interatomic separation, and A, F, and C (C6 and C12) are the forcefield parameters for a Buckingham (LJ) potential. In developing our forcefield, we decided to use LJ potentials rather than Buckingham potentials to describe all van der Waals interactions. This approach eases integration of the mineral forcefield with existing protein and water forcefields, as well as allowing simulations of larger systems due to the faster calculation of LJ energies compared with exponential Buckingham potentials. To improve the description of the structural and elastic properties of aragonite, we included an LJ potential for the CC interaction, which was neglected in Pavese’s work.17,20 For the covalent interactions between CO in the carbonate ion, 20068

dx.doi.org/10.1021/jp202743v |J. Phys. Chem. C 2011, 115, 20067–20075

The Journal of Physical Chemistry C harmonic angular potentials and a dihedral potential were employed to restrain the OCO angle and carbonate plane, respectively. The total set of potential functions is: Vtotal ¼ Vnonbond þ Vangle þ Vdihedral ( !) N1 N e 2 zi zj c12 c6 ¼ þ rij rij12 rij6 j¼1 i¼j þ 1

∑ ∑

þ

∑

1

angle 2

Kθ ðθ θo Þ þ

∑ Kψð1 cos 2ψÞ

plane

where zi, zj are the atomic charges, θ is the OCO angle, θ0 is the equilibrium angle value of 120°, and Ψ is the torsional angle of the O1CO2 and O2CO3 planes in carbonate. The General Utility Lattice Program 3.4.9 (GULP)27 was used to perform the forcefield development and optimization, which was carried out in the NpT ensemble with a distance cutoff of 1 nm for the LJ potential. In previous reports of CaCO3 forcefield development, the charge of calcium was fixed at its full ionic value of +2e; however, we found that sightly lowering the charge of calcium leads to a better fitting to the experimental structural and elastic values of aragonite, which was also noted in Pavese’s work.17 In our present work, we allowed the charge of all atoms, calcium, carbon, and oxygen, to vary during the fitting procedure, implicitly taking charge-transfer effects into account. The experimental static unit-cell parameters of a bulk aragonite crystal,28 listed in Table 2, and the fractional coordinates of the atoms were used in an initial round of conventionalmode fitting. The elastic constant matrix was then added as an additional optimization target in a subsequent round of relaxationmode fitting. Previous forcefield developments17,19,23 relied on elastic constant matrix data published a century ago by Voigt29 and later by Hearmon;30 we instead chose the more recent values reported by Liu et al.,31 because they are measured for natural single-crystal aragonite by Brillouin spectroscopy at ambient conditions, rather than with the previously used ultrasonic technique, which is believed to be very sensitive to the impurities of the natural crystal. MineralWater Interactions. As for pure CaCO3, several sets of forcefield parameters for CaCO3water interactions have been derived in the past for thermodynamic studies, especially for studying the process of crystallization from solution.21,22,32,33 These earlier works have demonstrated that the cross-term parameters between mineral and water are highly dependent on the forcefield of the mineral. Previously, calcite has always been chosen as the mineral phase in computational studies of CaCO3 and water interactions. We derived cross-term parameters for aragonitewater interactions based on our aragonite forcefield. The derivation of cross-term interaction parameters between mineral and organic molecules is still a big challenge for twophase simulations, especially for mineralpeptide interactions. The widely used method proposed by Schr€oder et al.34 and developed by Freeman et al.35 allows the inference of cross-term forcefield parameters from existing single-phase forcefields by fitting to the known structures of crystals containing atoms from both phases (e.g., a hydrated crystal incorporating both mineral and water). A number of forcefields have been successfully derived using this procedure.3537 We utilized a modified form of this method, taking advantage of the TIP3P water model38 and the structure of monohydrocalcite (CaCO3 3 H2O)39 (Table 2) to obtain CaCO3water interaction parameters.

ARTICLE

Our approach differed from that of Freeman et al. in two ways. First, we included both attractive and repulsive terms for the van der Waals interaction between the calcium ion and the water oxygen Ow, in contrast to earlier works that included only a repulsive term for CaOw. This neglect of an attractive van der Waals interaction is generally acceptable because of the high CaOw Coulombic attraction, but in our case the calcium charge was decreased from its fully ionic charge (+2e), and thus we allowed an attractive term in the CaOw van der Waals potential to compensate for the electrostatic attraction. Second, in Freeman’s original work,35 the structure of ikaite (CaCO3 3 6H2O) was chosen as the model to fit the cross-term forcefield parameters for CaCO3 and water; however, we instead used monohydrocalcite (CaCO3 3 H2O) as the target hydrated CaCO3. The CaCO3water cross-term parameters thus obtained were then tested by applying them to model the structure of ikaite40 (Table 2). The good agreement obtained between the model and the experimental structure stands as a validation of our new mineralwater interaction forcefield. MineralProtein Interactions. Freeman et al.35 have used a similar method to that described above for mineralwater interactions to derive a forcefield for mineral and organic molecules, with which they have successfully performed simulations on the protein-induced polymorphic transformation of CaCO3.41 But their forcefield was fitted to the structure of calcite and not designed primarily for mechanics; thus, it is necessary for us to derive a new set of mineralprotein interaction parameters based on our own aragonite forcefield. The simplest means of doing this is to use a combination rule to generate pairwise van der Waals interaction parameters from atomwise parameters taken from existing forcefields. One problem with this approach is that atomic charges on atoms shared between the two systems (such as oxygen) are typically not consistent between the two single-phase forcefields, which can invalidate the simple combination-rule approach for obtaining LJ parameters.35 However, in the case of our new CaCO3 forcefield, the carbonate oxygen atomic charge is relatively low at 0.889e, which is very close to the oxygen charge used in protein forcefields. Compared with earlier CaCO3 forcefields that assign a higher charge to oxygen, our forcefield suffers little from such an inconsistency of the oxygen charge among the different phases. Therefore, we adopted the simple approach of deducing pairwise mineralprotein LJ parameters from atomwise mineral parameters (taken from our forcefield) and atomwise protein parameters (from the OPLS_AA forcefield24), using the following combination rule:42 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C12ij ¼ C12i C12j pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð3Þ C6ij ¼ C6i C6j where C12 and C6 are the parameters LJ potential in eq 2. Molecular Dynamics Simulations. All-atom molecular dynamics (MD) simulations of pure aragonite, aragonitewater and aragonitewaterpeptide systems were carried out with Gromacs 4.0.443 The OPLS_AA24 forcefield was used to describe all proteinprotein and proteinwater interactions, and water was treated using the TIP3P model.38 The following simulation environment was used unless specified otherwise. Periodic boundary conditions were employed. A distance cutoff of 1.0 nm was used for LJ interactions, and the time step was set to 0.002 ps. Simulations were performed in the NpT ensemble with a temperature of 300 K and a pressure of 1 bar. We used NoseHoover44,45 temperature coupling with a time constant of 20069


The Journal of Physical Chemistry C 0.1 ps, and Berendsen pressure coupling46 with a time constant of 1 ps. After a steepest-descent energy minimization to remove atomic overlaps, systems incorporating water were subjected to a 500 ps solvent-relaxation run, in which the mineral and peptide atoms were position-restrained by a harmonic potential. Then a 5 ns equilibration with no position restraints was carried out before proceeding with production runs. To calculate elastic moduli for the three crystal axes, FPMD simulations11 were carried out in which the crystal was subjected to a uniaxial stretching force in vacuum. To prepare the system, a nonpulling simulation was first carried out of a cubic box containing 125 aragonite unit-cells. Given the 3D periodic boundary conditions, this system was equivalent to an infinite bulk aragonite crystal. The root-mean-square deviation (rmsd) and the total energy of the crystal were both found to be stable over 50 ps. To facilitate pulling the mineral crystal, we then rebuilt the mineral system by lengthening the box along the pulling direction, leaving a vacant gap on either side, resulting in a semi-infinite mineral crystal slab extending indefinitely (via periodic boundary conditions) along the two nonpulled axes but with finite extent along the pulled axis (Figure 1A, top). Three such slabs were constructed, with the finite axis along the x, y, and z directions, and each contained a total of 20003600 atoms. A pair of harmonic potentials with spring constants of 500 kJ mol1 nm2 were attached to the center of mass of the first layer of atoms in each of the pulled crystal surfaces. The harmonic potentials were then separated with a constant velocity of 0.2 nm/ns. The pulling force and elongation were monitored during simulations, and the elastic modulus was then calculated from the ratio of stress and engineering strain according to: E¼

stress F=A ¼ strain ΔL=Lo

where F is the real-time force applied to the atoms, A is the crosssection area, ΔL is the elongation of the separation between the two ends of the mineral crystal, and L0 is the initial separation. To assess the ability of our forcefield to describe the mechanical properties of a crystal in an aqueous environment, we also performed FPMD simulations for finite aragonite nanocrystals surrounded by water. Three different tablets were built for studying mechanical perturbation along each of the three axial directions, with a space of about 2.5 nm between each mineral surface and the nearest box boundary to ensure that the crystal did not interact with its periodic image. The space around the crystal was filled with TIP3P water,38 resulting in a total of roughly 30 000 to 60 000 atoms (depending on the system), as illustrated in Figure 1A (bottom). The first three layers of one surface along the long direction were position-restrained to prevent tablets from swaying during simulations. The elastic moduli of the mineral slabs in water were calculated using the same technique described above for pure aragonite in vacuum. Ab Initio Calculations. Ab inito optimization and interaction energy calculations of a single hydrated calcium ion were performed to compare with MD results as a validation of the lower calcium charge used in our forcefield. These calculations were performed with the Gaussian 03 package47 using the B3LYP method, which is based on an empirically parametrized hybrid functional including HartreeFock exchange and gradient corrections to the density. We used the 6-311+G** basis set, which has previously been demonstrated to well describe the CaH2O interaction.48 The ab initio optimization was followed

ARTICLE

Table 1. Forcefield for Pure Aragonite, WaterAragonite, and the Extracted C6, C12 Parameters of Calcium, Carbon, and Oxygen in CaCO3 with a LJ Cutoff of 1 nma Aragonite ZCa(e) = +1.668 ZC(e) = +0.999

Kθ (kJ/mol/rad2) = 1852

ZO(e) = 0.889

KΨ (kJ/mol) = 28.9 C12 (kJ/mol/nm12)

pairwise LJ

C6 (kJ/mol/nm6)

Ca

Om

9.49 107

0

Cm

Cm

4.61 106

1.43 102

Cm

Om

Om

Om

10

3.08 104

7

5.21 105

9.04 10 5.94 10

AragoniteWater pairwise LJ

C12 (kJ/mol/nm12)

C6 (kJ/mol/nm6)

Ca

Ow

8.85 107

2.00 103

Cm

Ow

5.91 106

3.39 103

Om

Ow

2.07 106

2.25 103

Om

Hw

7.98 109

9.65 107

atomwise LJ

C12 (kJ/mol/nm12)

C6 (kJ/mol/nm6)

Ca

2.52 107

1.42 103

1.44 10

5

4.61 103

1.77 10

6

2.03 103

Cm Om a

C6 and C12 parameters for calcium, carbon, and oxygen in CaCO3 are used to describe the VDW interaction between aragonite and peptides from OPLS_AA.

by vibration analysis to ensure the absence of negative frequencies and that the true minimum had been reached. From the resulting minimized structure, the CaOw distance was measured as the reference for our forcefield performance in classical MD simulations, and the single point energy of the whole hydrated system was calculated for further energy comparison with our MD simulations.

’ RESULTS Our new CaCO3 forcefield parameter set is shown in Table 1. Using the atomwise LJ potential parameters of Ca, C, and O in CaCO3, along with known atomwise parameters for protein atoms, pairwise mineralprotein LJ interactions can be obtained from the combination rule given in eq 3. The structure of pure CaCO3 and hydrated CaCO3, elastic properties of CaCO3, CaCO3water interaction and CaCO3peptide adsorption were simulated and compared with ab initio DFT calculations and experimental results to validate our forcefield. New Forcefield Accurately Reproduces Anisotropic Elastic Properties and Surface Energy of Aragonite. To verify our

newly developed forcefield parameters, a number of properties such as the unit-cell structural parameter, axial elastic modulus, and shear modulus of bulk aragonite were calculated by static calculations in GULP27 and are compared with those from previous simulations and experiments in Table 2 and Figure 1B. For the shear modulus, GULP provides three different calculation methods; we used the average of the three values as the final 20070



ARTICLE

Table 2. Unit-Cell Parameters of Aragonite, Calcite, Vaterite, Monohydrocalcite, and Ikaite, Determined Using GULP and the Forcefield in Table 1 mineral

unit-cell parameter

aragonite

vaterite

C33

C12 C13 C23 C44 C55 C66

Gb

aragonite 174.8 112.9 104.7 67.9 41.1 56.1 40.1 26.6 45.6 35.9c

0.00

98.4 60.3 27.8 41.9 39.3 24.2 40.2 35.8d31 59.2 65.3 39.0 48.2 40.5 33.9 49.0 33.8e17

0.00

155.3 104.2

89.9 55.9 54.7 48.0 36.7 12.4 23.3 25.0e23

162.8

87.5 71.9 75.1

32.5

45.5 29.0c

149.4

85.2 57.9 53.5

34.1

45.8 32.7d61

152.7

77.8 56.8 53.9

30.4

48.0 33.6e17

86.6

170.6 40.6 61.2

9.9

0.00

b (Å)

7.97

7.97

c (Å)

5.74

5.74

3

226.90

226.92

0.01

a (Å) b (Å)

5.03 5.03

4.99 4.99

0.80 0.80

c (Å)

17.53

17.07

2.69

3

V(Å )

383.94

368.20

4.27

a (Å)

7.30

7.29

0.14 0.14

b (Å)

7.30

7.29

c (Å)

25.46

25.30

0.63

1173.71

1164.50

0.79

a (Å) b (Å)

10.58 10.58

10.55 10.55

0.28 0.28

c (Å)

7.48

7.56

calcite

vaterite

23.0 18.1c

C22 = C11, C23 = C13, C55 = C44, and C66 = (C11 C12)/2 for calcite and vaterite crystals. b Shear modulus. c This work. d Experimental results. e Other theoretical results. a

Table 4. Lattice Energy (kJ/mol) of the Minerals from GULP this

Pavese’s

Raiteri’s

1.06

mineral

forcefield

forcefield17a

forcefield21a

expt62 2820

725.37

729.79

0.61

aragonite

3391.76

3645.62

2887.55

a (Å)

8.93

8.87

0.68

calcite

3381.99

4046.23

2887.82

2814

b (Å)

8.41

8.23

2.19

vaterite

3362.30

2876.58

2777

c (Å)

11.14

11.02

1.09

β (deg)

112.17

110.2

1.79

V (Å3)

774.98

754.98

2.64

result. Force probe molecular dynamics (FPMD) simulations on semi-infinite aragonite slabs in vacuum were also performed, as described in Methods. The elastic moduli obtained from GULP and from MD simulations are both reported in Figure 1B. The comparison of elastic moduli from MD, GULP, and experiments demonstrates not only that FPMD simulations of mineral slabs are able to reproduce the mechanical properties of a bulk material but also that our forcefield captures the elastic properties of aragonite in all-atom MD simulations rather better than the earlier forcefields of Pavese and Raiteri, especially along the z [001] axis. We also simulated stretching of finite aragonite prisms surrounded by water and found that the resulting elastic moduli agree with those for the semi-infinite slabs in vacuum. The properties of the surfaces of materials are crucial for determining the interaction between the substance and the surrounding environment. The most fundamental surface property, the surface energy γ, is defined as the thermodynamic energy penalty for creating a new surface from a bulk crystal: γ¼

C22

171.1 110.1 164.4 112.0

4.96

V (Å3) ikaite

C11

Δ (%)

4.96

V (Å3) monohydrocalcite

expt

a (Å)

V (Å ) calcite

theor

Table 3. Elastic Constant Matrix for Aragonite, Calcite, and Vaterite Obtained from GULPa

Usurf Ubulk A

where Ubulk is the bulk energy, Usurf is the energy for the system with the surface, and A is the created surface area. Three quasitwo-dimensional slabs (Figure 1A) were built according to this definition, and the mineral atoms in the center of the crystal were position-restrained to mimic bulk mineral. The surface energy of an aragonite mineral unit-cell was found to be 1.203 ( 0.007, 0.783 ( 0.005, and 0.773 ( 0.003 J/m2 for the x, y, and z surfaces, respectively, which agrees well with previous experimental and theoretical results of about 1 J/m2.18,49 This shows that our forcefield is capable of reproducing some thermodynamic properties of the crystal, although it was developed primarily to capture the mechanical properties. The positive

a

Estimated from the published forcefield.

surface energy on all three surfaces implies that aragonite tends to grow into a bulk crystal. The surface energy on the x-axis ([100]) is significantly higher than that of the y ([010]) and z ([001]) surfaces, implying that the [100] surface is less stable than the other two surfaces. According to Gibbs50 and Leeuw et al.,18 surface energies can be assumed to determine the equilibrium morphology of the crystal. Thus, the surface energy calculated from our atomic model is consistent with the result of aragonite crystallography that the [100] surface is not exposed in the aragonite equilibrium morphology.51 Transferability to Other CaCO3 Polymorphs. The structural and elastic properties and surface energies calculated above are evidence for the quality of our forcefield when applied to aragonite. It would be desirable to obtain a general forcefield applicable to every polymorph of CaCO3. Most calcium carbonate forcefields in the literature have been developed specifically for calcite. Because of differences in spatial atomic arrangement, these earlier forcefields have had only limited transferability to all three polymorphs of CaCO3. To check the transferability of our new aragonite forcefield, we used it to model the other two CaCO3 polymorphs, calcite and vaterite, in GULP. The structural and elastic matrix results are listed in Tables 2 and 3. Note that the calculated unit-cell parameters, the elastic moduli and average shear moduli all show excellent agreement with experiments, and thus support the transferability of our aragonite forcefield to other CaCO3 polymorphs. An important thermodynamical property, the lattice energy, which is the energy required to separate the ions from the crystal lattice and into the gas phase, is also estimated for all three polymorphs by our forcefield, and the values are compared with former theoretical and experimental results in Table 4. In this case, Raiteri’s forcefield21 produces the most accurate absolute value. Although the lattice energy obtained from our forcefield does not exactly overlap with the experimental values, the correct ordering is maintained for the lattice energy values of calcite, aragonite and vaterite, stressing 20071



ARTICLE

the CaOw distance derived from our forcefield (Figure 2A) indeed agrees better with the ab initio results than does the OPLS_AA forcefield,24 which uses the full charge of +2e. The interaction energy of calcium and water of the optimized structures from ab initio were also calculated by both ab initio and MD methods (Figure 2B). The interaction energy is defined as Eint ¼ Etotal EfCag EfðH2 OÞn g

Figure 2. Validation of the reduced calcium charge by comparison to ab initio simulations. CaH2O distance (A) and interaction energy (B) as a function of the water coordination number. Ab initio CaH2O distances and interactions are obtained using Gaussian 03 at the B3LYP/6311+G** level, and other data are from molecular dynamics simulations by Gromacs with the forcefield in Table 1 and OPLS_AA.

the power of our transferable forcefield to semiquantitatively reproduce physical quantities which were not used as fitting targets during forcefield optimization. Our aragonitewater cross-term forcefield was validated using the structures of monohydrocalcite (CaCO3 3 H2O) and ikaite (CaCO3 3 6H2O), with results shown in Table 2. Note that the deviations of the unit-cell parameters are smaller than 3% throughout, suggesting our forcefield can well reproduce the structure of a hydrated CaCO3 mineral. The calculated lattice energies for ikaite and pure calcite (aragonite) are about 7451 and 3382 (3391) kJ/mol, respectively. Assuming 621 kJ/mol for a gaseous water molecule,18 we can calculate the dissociation energy per water in ikaite to be 5557 kJ/mol. This value is comparable to previous theoretical results of 30 kJ/mol32 and 47 kJ/mol18 and the experimental value of around 50 kJ/mol.52 Lower Charge on Calcium Ion Validated by ab Initio Calculations. For further investigation of the robustness of our cross-term parameters, we used the new forcefield to simulate both a solvated calcium ion and a solvated carbonate ion, and the results were compared with those obtained using the OPLS_AA forcefield24 as well as ab initio calculations (see Methods for details). The solvated calcium system comprised a single calcium ion and several water molecules in a 12 nm3 cubic box. The average CaOw distance is plotted against the coordination number in Figure 2A, for both our forcefield and OPLS_AA.24 The ab initio results are also shown in the figure for reference. Our forcefield leads to a larger deviation of the CaOw distance from the ab initio results at small coordination number, compared with OPLS_AA,24 because of the decreased calcium charge. However, in nature, the charge of calcium would be slightly reduced with respect to its full ionic charge of +2e with increasing water coordination due to the charge-transfer effect. The actual charge of calcium in its solvated state was reported using density functional theory by Pavlov et al.53 According to their results, for a water coordination number between 5 and 7, the effective calcium charge would decrease to about +1.55e, which is close to our fitted charge of +1.668e. Because the coordination number of aqueous calcium is normally between 6 and 8 in experiments,48 our fitted calcium charge should be more realistic than the full ionic charge used in previous forcefields for solvated calcium. At higher coordination number,

where Etotal is the total single point energy of the cluster, E{Ca} and E{(H2O)n} is the single point energy of the separated calcium ion and water cluster of the frozen optimal structure. Both OPLS_AA24 and our new forcefield were used to calculate single point energies in MD simulations. Just as for the CaOw distance in Figure 2A, the OPLS_AA24 energy profile agrees better with the ab initio results at low coordination number, while our forcefield gives the better result for coordinate numbers in the range of 68 (Figure 2B), which is the appropriate water coordination number range for a calcium ion in aqueous solution. Further simulations of the mean water distribution around individual solvated calcium and carbonate ions, as well as at aragonite crystal surfaces, also show good agreement with previous simulations and experiments; see the Supporting Information for details. Proteins Incorporating Calcium or Carbonate Ions. Several proteins that incorporate calcium ions or carbonate anions were obtained from the RCSB Protein Data Bank (access IDs: 2I4C, 3GP2, 3NNG, 3N1U, and 3E3F). MD simulations on these proteins allowed a direct evaluation of our cross-term parameters for the interactions of calcium and carbonate with organic molecules. Equilibrium simulations of 10 ns were performed after energy minimization and position restraint simulations. To estimate the local structural stability, the rmsd of residues within 3 Å of the calcium ion were calculated and found to be all below 0.15 nm; that in proteins incorporating a carbonate anion cluster was slightly higher but still below 0.2 nm, demonstrating that our forcefield can successfully model proteincalcium and protein carbonate interactions. Forced Desorption of an Acidic Polypeptide from an Aragonite Crystal Surface. We applied our new mineral waterprotein forcefield to the simulation of an acidic amino acid adsorbed on the z surface ([001]) of an aragonite slab. Many experimental results have been reported on the organicinorganic interface in nacre.54,55 Pif 80,56 an 80-kD protein with dozens of Asp-rich domains, is one of the purported aragonite-binding proteins. To mimic the adsorption of one Pif 80 protein domain to an aragonite surface, we built an atomic model of an [Asp]7 peptide bound to the z surface of an aragonite slab. Two possible binding configurations were modeled: the peptide plane parallel to the mineral surface, in which case the main-chain atoms are closest to the aragonite surface, or perpendicular with the mineral surface, in which case the side-chain atoms are closest to the surface. The interaction energy of the main chain and side chains of the peptide and mineral for the two cases are shown in Figure 3. Independent of the configuration, the interaction energy between the mineral surface and the charged side chains exceeds that of the surface with the main chain. When the mainchain oxygen atoms are closer to the mineral calcium ion, i.e., the peptide plane is parallel to the surface, the adsorption of the backbone to the mineral surface is energetically unfavorable. Even if the peptide is initially placed parallel to the mineral surface and the side chains are relatively far from the mineral surface, the peptide backbone is found to rotate so that the acidic 20072



ARTICLE

Figure 3. Interaction energy of poly-Asp main chain (A) and side chain (B) with aragonite for different initial positions. The higher interaction energy between side chains and aragonite indicated that side chainaragonite adhesion is the crucial interaction between aragonite and peptides. This is consistent with the rotation of the peptide when side chains are placed far away from the aragonite surface (black line in B).

side chains can coordinate with the calcium ions, as illustrated by the decrease in side chainmineral interaction energy in Figure 3B. This phenomenon, the preference of aragonite to coordinate Asp side chains, agrees well with the fact that mineral-binding proteins are rich in amino acids with long acidic chains.57 FPMD simulations were performed to evaluate the adhesion force between the [Asp]7 peptide and the mineral surface. After equilibration of the mineralwaterpeptide system, the Cα at one end of the polypeptide was attached to a harmonic pulling potential, with a force constant of 500 kJ/mol/nm, and pulled away from the surface at a constant velocity of 0.2 nm/ns. The corresponding force is plotted versus time in Figure 4B. The three peaks in the force trace correspond to the sequential desorption of three side chains from the aragonite surface. The average rupture force required to pull one side chain away from the aragonite surface is around 700 pN, indicating a strong binding interaction between the peptide and the mineral surface. Note that the dissociation force of one Pif 80 protein from the aragonite surface should be at least 1 order of magnitude higher than this force because each protein has dozens of acidic binding domains.56 By comparing with the reported unfolding force of 700 pN at a similar pulling velocity for protein I27,14 one of the most robust protein domains known, we can conclude that when tensile forces are applied to nacre, the proteins and other organic components between the aragonite layers should become strongly deformed, with any native structure unfolded, before disassociation from the mineral surface occurs. This hypothesis based on comparing force magnitudes can also be explained from an energetic point of view: the large deformation in the organic layer effectively absorbs energy without unbinding from the mineral surface and thus protects the integrity of the whole composite; however, when the external perturbation is released, the organic components could return to their former relaxed state, slowly liberating the energy stored during deformation. Such a mechanism has already been proposed and discussed by many previous theoretical studies,3,58,59 but we are here able to support it for the first time with quantitative atomicscale simulations. To investigate the unbinding force of the acidic side chains from the CaCO3 surface at different pulling speeds, FPMD simulations were performed again covering a range of speeds,

Figure 4. Forced desorption of poly-Asp from aragonite. (A) The representive structures are shown, as indicated in the force profile. The atom subjected to the pulling force in shown as a yellow sphere. (B) Force profile for pulling poly-Asp away from the aragonite surface. Each sawtooth in the force profile corresponds to the dissociation of a side chain from aragonite. The dissociation of the first residue gives rise to the biggest rupture force, being further stabilized by the remaining Asp aragonite interactions. (C) Dependence of the rupture force in the pulling velocity. Error bars show standard errors of the mean. The red line represents a linear fit to the data and allows an approximate extrapolation to experimental time scales.

ranging from 0.05 to 8 nm/ns. At least five independent simulations were performed for each velocity. The measured increase in average rupture force relative to pulling speed is illustrated in Figure 4C. Following the revised Bell’s model,60 we have fitted the following linear model between the rupture force and the logarithm of pulling speed: Fr ¼ 52:7 lnðvÞ þ 785 where Fr is the rupture force and v is the pulling velocity. Extending the linear model to the much slower pulling speed of ∼1000 nm/s typical of AFM experiments,54 we predict a rupture force under AFM conditions of approximate 60 pN per acidic side chain. There are no experimental results at the single-molecule level available for direct comparison with our simulation result; however, using our predicted rupture force and a simplified model, we can estimate how much force is needed to peel an organic layer away from an aragonite tablet. We assume that the interaction range of a Pif 80 protein on the aragonite surface is 5 nm (the length of Pif 80 is 460 amino acids), that every Pif 80 protein has 10 domains bound to the mineral,54 and that the area of the proteinmineral interface is roughly 2.6 108 nm2.9 20073


The Journal of Physical Chemistry C Then the adhesion force between the organic layer and the whole aragonite tablet is calculated to be about 2 106 nN. According to the literature, the rupture force of aragonite tablets in nacre is 0.515 106 nN,3 which is comparable with our estimated force required to separate the organic layer from the mineral. This agreement supports the hypothesis given in many former theoretical studies: that the mineral itself and the mineralorganic interface may break simultaneously.58 The direct interaction between protein and mineral has been measured using AFM by Mohanty et al.,54 and a large rupture force (>5 nN) was reported for the desorption of proteins from aragonite under the AFM tip. They estimated the interaction area between the AFM tip and aragonite to be 1256 nm2. On the basis of our approximations above about the protein size and the number of acidic domains in aragonite-binding proteins, if the aragonite-binding protein Pif 80 were to completely cover the contact area of 1256 nm2, the rupture force to pull all proteins away from the aragonite tablet should be about 10 nN, a value consistent with the experimental results. Our assumption that Pif 80 proteins would completely cover the surface of aragonite is unlikely, and therefore our theoretical rupture force estimate should actually be treated as an upper bound.

’ CONCLUSION Many forcefields for aragonite and waterCaCO3 interactions have been developed, mainly focusing on the thermodynamic properties of the mineral and mineralization process.1723 For mechanical applications, one rather requires an atomic CaCO3 forcefield that can well reproduce the elastic properties of the material. We have derived a new set of forcefield parameters which accurately captures the elastic and shear modulus for bulk aragonite crystal. Although the forcefield is derived for aragonite, it has shown good transferability to the other two important CaCO3 polymorphs, calcite and vaterite. Based on our aragonite forcefield, the TIP3P water model,38 and the OPLS_AA forcefield24 for proteins, we have also developed a set of aragonitewater and aragoniteprotein cross-term parameters. Despite our forcefield being fitted primarily to the mechanical properties of aragonite, the parameters also give good agreements with a numbers of known thermodynamic properties from experiments. A notable aspect of our forcefield is the reduced calcium charge of +1.668e, in contrast to the full ionic charge of +2e used almost exclusively in the previous forcefield, and we have validated this reduced charge using ab initio results and MD simulations. To evaluate the interaction between proteins and aragonite, FPMD simulations of the forced desorption of a poly-Asp peptide chain from the aragonite surface have been performed, and the resulting rupture force agrees well with previous theoretical and AFM results. Overall, the reported structural and adsorption force analysis show that our new forcefield allows the realistic modeling of mineral mechanics and protein mineral interactions. Thus, these forcefield parameters are suitable to further study the molecular-scale mechanics of aragonite mineral tablets and biocomposite materials containing CaCO3 and protein. ’ ASSOCIATED CONTENT

bS

Supporting Information. Detailed structure analysis of the calcium, carbonate ion, and mineral surface solvation from

ARTICLE

MD simulations. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Phone: +49 6221-533267. Fax: +49 6221-533298.

’ ACKNOWLEDGMENT We acknowledge financial support from the Max Planck Society and from the Klaus Tschira Foundation and thank Julian D. Gale for sharing the forcefield parameters as well as the valuable and friendly communications and Fei Xia for helpful discussions. S.A.E. is further supported by a CAS Young International Scientist Fellowship (O91GC11401) and an NSFC Research Fellowship for International Young Scientists (O93DC11401). ’ REFERENCES (1) Gosline, J. M.; Guerette, P. A.; Ortlepp, C. S.; Savage, K. N. J. Exp. Biol. 1999, 202, 3295. (2) Luz, G. M.; Mano, J. F. Philos. Trans. R. Soc., A 2009, 367, 1587. (3) Ji, B.; Gao, H. J. Mech. Phys. Solids 2004, 52, 1963–1990. (4) Stempfle, P.; Pantale, O.; Rousseau, M.; Lopez, E.; Bourrat, X. Mater. Sci. Eng., C 2010, 30, 715–721. (5) Tushtev, K.; Murck, M.; Grathwohl, G. Mater. Sci. Eng., C 2008, 28, 1164–1172. (6) Katti, K. S.; Katti, D. R. Mater. Sci. Eng., C 2006, 26, 1317–1324. (7) Nassif, N.; Pinna, N.; Gehrke, N.; Antonietti, M.; J€ager, C.; C€olfen, H. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 12653. (8) Tang, Z.; Kotov, N. A.; Magonov, S.; Ozturk, B. Nat. Mater. 2003, 2, 413–418. (9) Meyers, M. A.; Lin, A. Y.; Chen, P. Y.; Muyco, J. J. Mech. Behav. Biomed. 2008, 1, 76. (10) Yamakov, V.; Wolf, D.; Phillpot, S. R.; Mukherjee, A. K.; Gleiter, H.; et al. Nat. Mater. 2002, 1, 45–49. (11) Grubm€uller, H.; Heymann, B.; Tavan, P. Science 1996, 271, 997. (12) Heymann, B.; Grubm€uller, H. Biophys. J. 2001, 81, 1295–1313. (13) Gr€ater, F.; Shen, J.; Jiang, H.; Gautel, M.; Grubm€uller, H. Biophys. J. 2005, 88, 790–804. (14) Stacklies, W.; Vega, M. C.; Wilmanns, M.; Gr€ater, F. PLoS Comput. Biol. 2009, 5, e1000306. (15) Xiao, S.; Stacklies, W.; Cetinkaya, M.; Markert, B.; Gr€ater, F. Biophys. J. 2009, 96, 3997–4005. (16) Baldauf, C.; Schneppenheim, R.; Stacklies, W.; Obser, T.; Pieconka, A.; Schneppenheim, S.; Budde, U.; Zhou, J.; Gr€ater, F. J. Thromb. Haemostasis 2009, 7, 2096–2105. (17) Pavese, A.; Catti, M.; Price, G. D.; Jackson, R. A. Phys. Chem. Miner. 1992, 19, 80–87. (18) de Leeuw, N. H.; Parker, S. C. J. Phys. Chem. B 1998, 102, 2914–2922. (19) Archer, T. D.; Birse, S. E. A.; Dove, M. T.; Redfern, S. A. T.; Gale, J. D.; Cygan, R. T. Phys. Chem. Miner. 2003, 30, 416–424. (20) Pavese, A.; Catti, M.; Parker, S. C.; Wall, A. Phys. Chem. Miner. 1996, 23, 89–93. (21) Raiteri, P.; Gale, J. J. Am. Chem. Soc. 2010, 132, 17623–17634. (22) Raiteri, P.; Gale, J. D.; Quigley, D.; Rodger, P. M. J. Phys. Chem. C 2010, 114, 5997–6010. (23) Fisler, D. K.; Gale, J. D.; Cygan, R. T. Am. Mineral. 2000, 85, 217. (24) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225–11236. (25) Buckingham, R. A. Proc. R. Soc. A - Math. Phys. Eng. Sci. 1938, 168, 264. 20074



ARTICLE

(26) Jones, J. E. Proc. R. Soc. A - Math. Phys. Eng. Sci. 1924, 106, 463–477. (27) Gale, J. D. J. Chem. Soc., Faraday Trans. 1997, 93, 629–637. (28) Jarosch, D.; Heger, G. Miner. Petrol. 1986, 35, 127–131. (29) Voigt, W. Lehrbuch der kristallphysik; BG Teubner: Leipzig, 1910; Vol. 34. (30) Hearmon, R. F. S. Rev. Mod. Phys. 1946, 18, 409. (31) Liu, L.; Chen, C.; Lin, C. C.; Yang, Y. Phys. Chem. Miner. 2005, 32, 97–102. (32) de Leeuw, N. H.; Parker, S. C. J. Chem. Soc., Faraday Trans. 1997, 93, 467–475. (33) de Leeuw, N. H.; Parker, S. C. Phys. Chem. Chem. Phys. 2001, 3, 3217–3221. (34) Schr€oder, K. P.; Sauer, J.; Leslie, M.; Richard, C.; Catlow, A.; Thomas, J. M. Chem. Phys. Lett. 1992, 188, 320–325. (35) Freeman, C. L.; Harding, J. H.; Cooke, D. J.; Elliott, J. A.; Lardge, J. S.; Duffy, D. M. J. Phys. Chem. C 2007, 111, 11943–11951. (36) De Leeuw, N. H.; Parker, S. C.; Rao, K. H. Langmuir 1998, 14, 5900–5906. (37) Duffy, D. M.; Harding, J. H. Langmuir 2004, 20, 7630–7636. (38) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (39) Swainson, I. P. Am. Mineral. 2008, 93, 1014. (40) Dickens, B.; Brown, W. E. Inorg. Chem. 1970, 9, 480–486. (41) Freeman, C. L.; Harding, J. H.; Quigley, D.; Rodger, P. M. Angew. Chem. 2010, 122, 5261–5263. (42) David, S.; Erik, L.; Hess, B. Gromacs user manual. version 4.0; University of Groningen: Groningen, Netherlands; 2006. (43) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701. (44) Nose, S. Mol. Phys. 1984, 52, 255–268. (45) Hoover, W. G. Phys. Rev. A 1985, 31, 1695. (46) Berendsen, H. J. C. Comput. Simul. Mater. Sci. 1991, 139–155. (47) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery , J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; et al. Gaussian 03; Gaussian; Inc., Pittsburgh, PA, 2003. (48) Megyes, T.; Grosz, T.; Radnai, T.; Bako, I.; Palinkas, G. J. Phys. Chem. A 2004, 108, 7261–7271. (49) Okumura, K.; De Gennes, P. G. Eur. Phys. J. E 2001, 4, 121–127. (50) Gibbs, J. W. The collected works of JW Gibbs; 1928; Vol. 2; pp 14. (51) Ford, W. E.; Dana, E. S. A textbook of mineralogy; John Wiley: New York, 1932. (52) Bischoff, J. L.; Fitzpatrick, J. A.; Rosenbauer, R. J. J. Geol. 1993, 21–33. (53) Pavlov, M.; Siegbahn, P. E. M.; Sandstr€om, M. J. Phys. Chem. A 1998, 102, 219–228. (54) Mohanty, B.; Katti, K. S.; Katti, D. R. Mech. Res. Commun. 2008, 35, 17–23. (55) Gilbert, P.; Abrecht, M.; Frazer, B. H. Rev. Mineral. Geochem. 2005, 59, 157. (56) Suzuki, M.; Saruwatari, K.; Kogure, T.; Yamamoto, Y.; Nishimura, T.; Kato, T.; Nagasawa, H. Science 2009, 325, 1388. (57) Weiner, S.; Addadi, L. Trends Biochem. Sci. 1991, 16, 252–256. (58) Gao, H.; Ji, B.; J€ager, I. L.; Arzt, E.; Fratzl, P. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 5597. (59) Smith, B. L.; Sch€affer, T. E.; Viani, M.; Thompson, J. B.; Frederick, N. A.; Kindt, J.; Belcher, A.; Stucky, G. D.; Morse, D. E.; Hansma, P. K. Nature 1999, 399, 761–763. (60) Ackbarow, T.; Chen, X.; Keten, S.; Buehler, M. J. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 16410. (61) Chen, C. C.; Lin, C. C.; Liu, L. G.; Sinogeikin, S. V.; Bass, J. D. Am. Mineral. 2001, 86, 1525. (62) Glasser, L.; Jenkins, H. D. B. J. Am. Chem. Soc. 2000, 122, 632–638.

20075


A New Transferable Forcefield for Simulating the Mechanics of CaCO3

Recommend Documents