J. CHEM. SOC. FARADAY TRANS., 1994, 90(4), 641-647 Computer Modelling of Phosphate Biominerals Transfer of Parameters for Interatomic Potentials for Different Polymorphs of Divalent Metal Diphosphates Marina G. Taylor and Kenneth Simkiss Department of Pure and Applied Zoology, University of Reading, P.O. Box 228,Whiteknights, Reading, UK RG6 2AJ Maurice Leslie SERC Daresbury Laboratory, Daresbury, Warrington, UK WA4 4AD Manganese, zinc and calcium diphosphates have been modelled using atomistic simulation techniques. The ionic model was used for the cation-anion interactions. Parameters for the short-range interatomic potentials previously calculated using the electron gas methods for the diphosphate anion of or-magnesium diphosphate, were used together with specific metal oxygen parameters.The structural data for each compound were repro- duced well, indicating that the anion parameters were transferable. These studies present new opportunities for modelling important aspects of biominerals which are not readily accessible to experimental studies. Calcium phosphates form the major inorganic component of mineralised tissue such as bone and teeth in vertebrates. Calcium phosphate appears as hydroxyapatite, frequently with some carbonate included in the lattice. Precursors to the final mineral have been proposed ranging from an amorp- hous calcium phosphate phase, found in matrix vesicles close to the mineralising front, through to brushite and octacal- cium phosphate. In addition, in many invertebrates tissues intracellular granules, which also have an amorphous struc- ture, are found which have an ionic composition which can be an orthophosphate or a condensed phosphate such as diphosphate (pyrophosphate).' The cations present are mostly calcium and magnesium but a feature of these bio- minerals is their ability to incorporate other cations such as manganese, cobalt, zinc and transuranics when present in the environment and so these materials can be seen as pro- viding a sink for pollutants and so detoxifying them.In the bone precursor, in which magnesium and zinc have been detected in addition to calcium,' the magnesium may have a role as a crystal inhibitor. The amorphous materials are pro- duced from aqueous solutions at environmental temperatures but have many of the properties associated with glasses. They are of considerable biological importance since they possess unusual solubility and solid-phase properties that have been exploited by living systems in a variety of mineralising, acid- base regulating and detoxifying systems.Analytical studies have established the composition of the materials and provided information on the nature of the phosphate species. Structural studies using X-ray absorption spectroscopy have been used to determine the local atomic structure around calcium and other cations in granules which have orthophosphate and diphosphate anions. 1*3-5 The dis- tribution of cations in these materials is not known but may be similar to those found in mixed alkali-metal-silicate Atomistic simulation methods of compounds with known crystal structures are being used to determine the sub- stitution energies of the dopant metals and subsequently the surface reactivity of these biominerals.These methods have been used successfully to model the static and dynamic properties of the perfect lattices of compounds such as zeo- litic silica polymorphsg and calcium carbonate." We have previously modelled the perfect lattice of cr-magnesium diphosphate" and we now present our results on the trans- ferability of parameters, determined for this magnesium com- pound, to other anhydrous diphosphate structures with calcium, manganese and zinc as the counterions. Structures Inorganic diphosphates exhibit two major classes of polymorphism" with differing structural features.Diphos- phates are condensed dimers of phosphate units with a formula of P,O$- which consists of two PO, units linked by an oxygen bridge in the P-0-P form. The conformations of the PO, units, staggered or eclipsed and the angle of the P-0-P bridge are the main features which characterise the two major types, thortveitite and dichromate structures. The thortveitite structure has the staggered conformations of the PO, groups and the P-0-P bridging angle greater than 140". Only in the dichromate-type structure does the bridging oxygen come into the first coordination sphere of the cation. When the structure of thortveitite, Sc,Si,O,, was first deter- mined, it was reported that the Si-0-Si bond angle was linear, i.e.180°.'2 This has been the subject of much dis- cussion, especially as other similar compounds have mostly non-linear bridges, but the linear bridge has subsequently been confirmed." Among the diphosphates several linear P-0-P bridges have been reported'3,'4 but it is generally considered that this might arise as a result of thermal dis- order. We report the modelling of various metal diphosphate compounds. The initial study was of a-magnesium diphos- phate, a compound crystallising with a thortveitite type struc- ture.' ' These studies have been extended to manganese and zinc diphosphates which are isomorphous with magnesium diphosphate and to /3-calcium diphosphate which is not and has the dichromate structure.The study was extended, in part to test the transferability of the parameters for interatomic potentials, but also because these compounds are relevant to our experimental studies of biominerals, particularly those related to the incorporation of foreign cations into the calcium magnesium diphosphate granules. Manganese diphosphate is unique in being the only diphosphate with the thortveitite structure at room tem-perature. In early X-ray studies of manganese diphosphate, a linear P-0-P bridging bond was indicated.', A more recent investigation of the structure' has suggested that there may be thermal disorder in the bridge and a split-atom model was proposed retaining the C,,, symmetry but with a bridging angle of 165.9'.Only an occupancy factor of 0.80 was found and even further disorder was suggested for the remaining fifth. A further structure with a bent model with the symmetry reduced to C, and a bridging angle of 164.5" could not be refined in the X-ray study. We have attempted to model the structures with a linear P-0-P bridge and the split atom model. The second set of diphosphates we have considered are the a-and p-zinc diphosphates which also crystallise in thortveitite type structures. Both forms are pre- pared from a melt obtained from the decomposition of zinc ammonium phosphate. For the higher temperature form, p-zinc diphosphate, a linear bridging P-0-P has also been proposed with the suggestion that this too may reflect thermal disorder though in this case there have been no further attempts to refine the crystal structure taking this into account.l4 The lower temperature phase, a-zinc diphosphate is described as a hexpartite structure in which the a-axis is tripled and the c-axis is doubled relative to the thortveitite structure.16 It is believed to have the most complex structure of the series of metal diphosphates. The structure consists of a sequence of layers, one of which is like that found in a- copper diphosphate followed by two layers like those found in a-magnesium diphosphate. In a-magnesium diphosphate both the a-and c-axes are doubled compared with the man- ganese structure.' 5,1 Calcium diphosphates have very complicated structures with coordination numbers ranging from seven to nine, calcium oxygen distances ranging from 2.250 to 2.807 8, in two eight-coordinate sites in the high-temperature ct-form'* and 2.318 to 2.927 in four different coordination sites with seven, nine, seven and eight nearest-neighbour oxygens in the lower temperature p-form.lg In the latter case in sites one and four, bridging oxygens come within 3 8, of the central calcium.The a-form has the thortveitite-type structure and the p-form the dichromate type structure. All the crystalline calcium diphosphates are high-temperature phases. The p-calcium diphosphate was model- led partly because it is the lower temperature form, but also becuse it crystallises in the dichromate structure, and is tetragonal whereas the magnesium, manganese and zinc crys- tals are of the thortveitite type and are monoclinic.In the calcium structure there are two types of diphosphate units each with an eclipsed configuration of the PO, units. The structure shows infinite chains of diphosphate units linked in all directions by calcium ions with calcium found in four dis- tinct sites. Methods Computer Modelling Computer modelling was carried out using the University of Reading's AMDAHL 5870 and a SUN SPARC station 1. The programs used were THBREL and THBPHON made available by the SERC CCP5 scheme. These programs were developed initially at AEA Harwell to study defects in solids20*2' using the Mott-Littleton approximation.22 Theoretical Methods The compounds were modelled as perfect ionic lattices with cations interacting with internally covalent bound diphos- phate anions. Two- and three-body terms were included to model the phosphorus-oxygen bonds in the anion.The lattice energy can be computed at constant volume and constant pressure as the sum of the long-range Coulom- bic interactions, the short-range non-bonded potentials and the force field terms for the covalent anion. = ECoulornbic + Eshort-range + Etwo +three body The lattice energy is minimised using classical optimisation techniques23 such as the Newton-Raphson method as a func- tion of structural parameters, lattice parameters and atomic J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 coordinates.At the minimum the derivatives of the energy with respect to the geometric factors will be zero. Bond lengths and bond angles are calculated and compared with the crystallographic data to give an indication of the good- ness of the potentials used. A further check on the inter- atomic potentials is obtained by generating the phonon spectrum. The eigenvalues calculated from the THBPHON program using a wave vector of 0.0o0, O.OO0 and 0.001 have been compared with the infrared spectra of metal diphos- phate corn pound^.^^ The starting parameters for the non-bonded interatomic potentials were mean values obtained by the electron gas method as previously de~cribed.~~-~' Values were calculated for each atomic pair in magnesium, manganese, zinc and calcium diphosphates.When required, more than one range of parameters was used for different cation-oxygen distances. Initially the charges and parameters which were determined earlier by the electron gas technique for the diphosphate anion in magnesium diphosphate have been used' and have been maintained as far as possible throughout. Some changes were made to the charge distribution within the anion to improve the fit in some models. Specific values have been cal- culated for the parameters for each cation-anion non-bonded repulsion contribution to the interatomic potential of the compounds of interest. For the electron gas potential calcu- lations the wavefunctions were calculated in the Madelung potential well to account for the local crystal environment.The method used for magnesium diphosphate' ' was to set up different charge models for phosphorus and the bridging and terminal oxygens in each crystal and calculate the Madelung potential. The electron densities were then determined from expansions of Slater functions modified to account for the Madelung potential experienced in the crystal. The charge on the cations was maintained as 2+. The charges on the atoms in the anions were initially those which were used to model magnesium diphosphate. These were the formal charges in the Langmuir sense, i.e. 1+ on phosphorus, 1-on the ter- minal oxygen and 0 on the bridging oxygen. These were adjusted within the anion to improve the fit obtained from the calculation when necessary.The parameters for the non- bonded repulsive potentials in the anion were maintained throughout, except where indicated in Table 1. Clementi free- ion wavefunctions were used for manganese and calcium as had been used previously for magnesium. There were no Table 1 Parameters for the Buckingham non-bonded interatomic potentials for static lattice simulations of diphosphate compounds species AIeV pIA C/eV A6 range/A Mg2+-0, 3551.55 0.243 32 6.00 0-20 Mn2+ -0, 2538.565 0.267 56 4.00 0-2.2 3028.565 0.267 56 4.00 2.2-20.0 a-Zn2+ -0, 1052.000 0.287 82 0.00 0-20 P-Zn2+ -0, 1107.000 0.287 82 0.00 0-2.15 1125.000 0.287 82 0.00 2.15-20 P-Ca2+ -0, 3000.750 0.272 98 0.00 0-2.40 3100.833 0.272 98 0.00 2.40-2.60 3500.750 0.272 98 0.00 2.60-20.0 fi-Ca2+ -0, 500.53 1 0.276 00 30.00 0-3.10 750.531 0.276 00 00.00 3.10-20.0 OT-0, OT-0, OFl--OB 1394.529 108 1.8 17 764.00 0.297 82 0.303 57 0.320 00 120.00 80.00 53.00 0-20.0 0-20.0 0-20.0 P-P 1030.43 0.356 62 0.00 0-20.0 P-0, 388.475 0.399 78 2.00" 0-20.0 P-0, 300.7 16' 0.428 76 2.00* 0-20.0 a.b Values of 10.00 and 0.00, respectively, were used for the calcium compound.A value of 370.716 was used in the manganese com- pound. J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 643 7.5 7.5 6.5 6.5 5.5 5.5 4.5 4.5 23 3.5 233.5 5: 5: 2.5 2.5 1.5. 1.5 0.5. 0.5 -0.5. I 1 I 2 I 3 1 4 -0.5 0 1 2 3 4 r/A r/A 7.5 7.5- 6.5 6.5 - 5.5 5.5 - 4.5 4.5 - % i.3.5r -..>,.. T 3.5-5 2.5- 2.5 - 1.5 1.5- 0.5 0.5 - -0.5. I 1 I 2 I 3 I 4 - 0 0 . 1 5 2 . 1 3 , 4 r/A r/A Fig. 1 Mg,Mn and Zn overlap except where indicated. V(r)for non-bonded repulsive terms us. interatomic distance for (a)0,-0,, (b)0,-0,, (c) P-0, and (d)P-0, . The curves for Ca, parameters for Zn2 + in Clementi-Roetti28 so the cation density function was calculated n~merically~~ with a rela-tivistic correction. In the electron gas method exchange energy terms were calculated using the Handler formula and the correlation energy terms were calculated using Wigner’s function. The electron gas data were then fitted to the analy- tic Born-Mayer function.The electron gas method produces a range of parameters depending on the local environment of each individual atomic species. The starting point for fitting the structures used a mean value of the calculated values of the A and p parameters. Any modifications to the A parameter during the fitting were, as far as possible, kept within the range of the calculated values. Plots of the Born-Mayer function using the parameters for the anion in each compound for OT-OT, OT-O,, P-0, and P-0, interactions are com- pared in Fig. 1. 0, represents a terminal oxygen and 0, a bridging oxygen. Although the numerical values of A and p were not the same for each compound, these parameters are correlated and it is seen that the potential energy is virtually identical for O,-OT and 0,-0, interactions and there are only small differences for the P-0, and P-0, interactions. Therefore these parameters obtained for magnesium were used in modelling the anion in the other compounds.Van der Waals terms for the short-range dispersion effects were calculated using the formula derived by Slater and Kirk- For the interaction between species i and j C,(ij) = 3/2ai aj/[(ai/Pi)”2 + (c~~/P~)”~] Table 2 Force constants used in the perfect lattice simulations species Mgz+ Mn2+ a-Zn2+ b-ZnZ+ B-CaZf two-body P-0, k/eV k2 25.0 25.0 25.0 25.0 25.0 bond length, rlA 1.590 1.568 1.589 1.569 1.615 P-0, kleV A -35.0 30.0 35.0 30.0 35.0 bond length, r/A 1.516 1.526 1.52 1 1.555 1.518 fraction of Coulombic term excluded 0.35 1.0 0.45 0.50 0.50 three body 0,-P-0, k/eV rad-2 10.0 15.0 12.5 12.0 10.0 mean equilibrium angle, eldegrees 112.33 112.50 1 12.09 111.75 1 12.60 0,-P-0, k/eV rad-2 15.0 16.0 12.5 16.0 15.0 mean equilibrium angle, eldegrees 106.33 106.2 106.65 106.50 105.4 P-0B-P k/eV rad-2 20.0 20.0 20.0 20.0 25.0 mean equilibrium angle, eldegrees 144.0 180.0 143.0 180.0 134.2 where ai is the static polarizability of species i and Pi the effective number of electrons contributing to the polarizabil- ity.A damping factor zDi,(r) was also included to reduce the dispersion energies when the wavefunction overlap is not neg- ligible.31 The short-range terms can be collected into the Buck- ingham potential The long-range electrostatic interactions between each ion pair are included with an additional term 4i4j/rij where 4i represents the charge on ion i and rij is the interionic dis- tance.An approximate calculation of the two body force con- stants was made by considering the energy of assigned vibra- tions in the infrared spectrum. The P-0 bonds were treated as harmonic oscillators Fj = ik(r -ro)2 Table 3 Charge models for diphosphates species Mg Mn or-Zn B-Zn /?-Ca cation 2+ 2+ 2+ 2+ 2+ phosphorus 1+ 1,8+ 1.05+ 1.80+ 1.2+ bridging" 0.0 -0.10 -0.10 -0.10 -0.40 oxygen, 0,terminal -1.0 -1.25 -1.00 -1.25 -1.00 oxygen, 0, " THBPHON program requires a non-zero charge on each species. J. CHEM. SOC. FARADAY TRANS., 1994, VOL.90 Table 4 Perfect lattice simulation of manganese diphosphate (Mn,P,O,); comparison of experimental (X-ray) and calculated values (a) Cell parameters and lattice energy cell parameters calculated experimental alA 6.42 6.63 blA 8.59 8.58 CIA 4.42 4.55 /?/degrees 103.2 102.7 lattice energy/eV -88 x 2 -87 x 2" (Z = 2) (b) Bond distanceslA and anglesldegrees species X-ray CONV CONP Mn-0, X 2 2.164 2.149 2.148 Mn-0, x 2 2.132 2.147 2.134 Mn-0, x 2 mean value 2.316 2.204 k0.080 2.313 2.203 f0.078 2.315 2.199 f0.082 P-0, 1.568 1.564 1.561 P-P 3.140 3.128 3.121 P-0, 1.519 k0.001 1.519 1.516 k0.005 <P-0,-P 180.0 180.0 180.0 <OB-P-o+ 106.2 & 1.3 107.2 f0.2 102.3 & 0.5 mean value <o,-P-0, 112.3 f0.3 11 1.0 f0.8 11 1.5 f0.1 mean value " From Born-Haber calculation.where (r -r,) represents the displacement from the equi- librium distance ro/A and k is the spring constant. The three- body force constants were then estimated by considering the assigned vibrations in relation to the two-body terms. where k, is the bond bending force constant and 8, is the equilibrium bond angle. Results The parameters for the non-bonded terms of interatomic potentials used for magnesium, calcium, manganese and zinc diphosphates are listed in Table 1. The force constants for the two- and three-body terms are also given (Table 2). The geometric fitting of manganese diphosphate with the linear bridge was good.The crystal structure shows man- ganese in octahedral coordination with two Mn-0 distances at 2.132 8, and two at 2.164 A. The other two bond distances were longer at 2.316 A and so two ranges of A parameters were required for fitting the Mn-0 bond lengths. These were selected on either side of the calculated mean value. The distribution of charge within the diphosphate anion was modified from the values used for the a-magnesium diphos- phate to improve the fit but the overall charge of the anion was maintained (Table 3). Calculated bond lengths in the structure relaxed at constant volume and pressure are given in Table 4 and compared with those determined by X-ray crystallography and neutron diffraction studies. There was a decrease in the cell dimensions of the relaxed structure com- pared with the X-ray determinations (Table 4).The calculated lattice energy, -88 eV compares with the value determined J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 Table 5 Perfect lattice simulation of a-zinc diphosphate (Q-Zn2P20,); comparison of experimental (X-rays) and calculated values (a) Cell parameters cell parameters calculated experiment a1 a/A 19.67 20.07 blA 7.86 8.26 CIA 8.63 9.09 flldegrees 105.8 106.4 lattice energy/eV -79.7 x 6 -71.0 x 6" (2= 6) (b) Bond distance/A and anglesldegrees ~~ species X-ray CONV ~~ CONP Zn(Cu)-0, 1.956 1.955 1.956 2.024 2.053 2.034 2.060 2.057 2.040 2.060 2.065 2.063 2.089 2.117 2.090 mean value 2.038 f0.046 2.049 f0.052 2.037 f0.045 Zn(Mgl)-oT 2.032 2.05 1 2.025 2.056 2.022 2.044 2.083 2.104 2.089 2.097 2.108 2.104 2.102 2.111 2.107 2.177 2.132 2.1 14 2.090 f0.046 2.089 f0.037 2.080 f0.035 1.928 1.951 1.945 2.025 2.002 1.998 2.032 2.087 2.042 2.056 2.115 2.105 2.093 2.115 2.117 2.027 f0.055 2.054 f0.066 2.041 _+ 0.065 1.603 1.578 1.595 1.599 1.586 1.595 1.566 1.580 1.596 3.Ooo 3.012 3.045 3.046 3.022 3.049 1.529 & 0.012 1.519 f0.004 1.511 f0.004 1.512 f0.016 1.519 f0.006 1.511 f0.005 1.523 f0.006 1.527 f0.003 1.516 f0.005 Plopl 139 145 145 p2,3°p2.3 148 145 146 111.8 f1.5 113.1 f3.5 109.0 f 1.0 106.9 & 3.8 109.1 f2.5 109.9 f1.6 111.2 f2.0 109.9 f1.3 109.2 f0.6 106.2 k 3.5 109.0 f2.9 109.7 f1.7 112.0 f0.4 109.6 k 1.9 108.7 f0.5 106.8 f2.3 109.3 f2.2 110.2 f2.2 a From Born-Haber calculation. from experimental data in a Born-Haber type calculation -87 eV. Although there are no literature values for elastic constants and relative permittivities the values obtained were judged reasonable.We were unable to fit the split-atom model containing the disordered bridging oxygen with THBREL. The structure relaxed but a new site for man-ganese binding appeared which was inconsistent with the crystallographic structure proposed by Stefandis and Nord.' Both a- and /.?-forms of zinc diphosphate were fitted. The charge distribution in the diphosphate anion in the a-form was maintained as in magnesium diphosphate.The param- eters for the non-bonded interatomic potentials in the anion were those used for modelling magnesium diphosphate. Some Table 6 Perfect lattice simulation of ,!&zinc diphosphate (p-Zn,P,O,); comparison of experimental (X-ray) and calculated values (a) Cell parameters and lattice energy cell parameters calculated experimental 4 6.57 6.6 1 biA 8.15 8.29 CI'A 4.39 4.5 1 Ptdegrees 105.1 105.4 lattice energy/eV -77 x 2 -71 x 2" (2= 2) (b) Bond distances/A and bond anglesldegrees species X-ray CONV CONP ~ Zn-0, x 2 2.00 1 2.001 1.998 Zn-0, x 2 2.06 1 2.067 2.080 Zn-0, x 2 2.275 2.245 2.27 1 mean values 2.112 & 0.1 18 2.104 f0.103 2.1 16 & 0.1 14 P-0, 1.569 1.546 1.556 P- P 3.137 3.093 3.1 13 mean value P-0, 1.555 f0.001 1.574 f0.010 1.569 f0.008 angles POBP 180.0 180.0 180.0 mean value %PO, 111.5 0.8 107.4 f3.7 108.7 & 1.9 110.2 f2.0 108.3 f1.2 110.6 f1.4 " From Born-Haber calculation.small modifications were made to the force constants for the two-and three-body terms (Table 2) to improve the fit obtained. The parameters used initially for the zinc-oxygen interatomic potentials were those calculated using the elec- tron gas method. A lower A value, however, improved the fit. In the a-form16 zinc is in three different sites with coordi- nation numbers of five, six and five. This feature was repro- duced in the relaxed structure.The zinc-oxygen and phosphorus-oxygen distances are given in Table 5 and are in reasonable agreement with the crystallographic data. The lattice energy and cell dimensions of the relaxed structure are compared with the experimental data in Table 5. Once again reasonable elastic constants and relative permittivities were obtained although there is no reference data available. The crystallographic struct~re'~ of the p-zinc diphosphate was solved with a linear P-0-P bond although the ques- tion of thermal disorder was discussed. As with the structure of manganese diphosphate, also with a linear bridging oxygen, the structure was successfully modelled only when the charges in the diphosphate anion were similarly modified. Two ranges of the parameters for the zinc-oxygen non-bonded terms were also required to reproduce the crystallo- graphic geometry (Tables 1 and 6).The calcium diphosphate structure did not model as well as magnesium, manganese and zinc diphosphates. The charge model was changed slightly so that the bridging oxygen carried a partial charge of -0.4 and the charge on the phos- phorus atoms was then increased to + 1.2 to maintain charge neutrality. The A parameters are consequently from the higher range of calculated values. The bond distances for the four calcium sites and the phosphorus-oxygen distances are given in Table 7. Most of the values compare well with X-ray crystallography, especially the nearest neighbours. There was a tendency for distant terminal oxygen atoms in sites one and three to move in closer to calcium and increase the coordi- nation number.Despite the proximity of the bridging Table 7 Perfect lattice simulation of p-calcium diphosphate (p-Ca,P,O,); comparison of experimental (X-ray) and calculated values (a) Cell parameters and lattice energy Cell parameters calculated experimental ~ a/A 6.590 6.684 6.590 6.684blA 24.144 24.144CIA (a = B = y)/degrees 90 90 lattice energy/eV (2 = 8) -73.8 x 8 -76.5 x 8" (b) Bond distances/A and bond anglesldegrees species X-ray CONV CONP Ca(1)-0, 2.340 2.364 2.332 2.360 2.315 2.334 2.369 2.365 2.336 2.409 2.380 2.362 2.416 2.455 2.449 2.457 2.556 2.512 mean value 2.392 & 0.039 2.442 f0.105 2.385 f0.021 2.814 2.89 1 Ca(l)-oB Ca(2)-0, 2.780 2.342 2.397 2.364 2.42 1 2.351 2.378 2.414 2.455 2.469 2.509 2.470 2.518 2.557 2.55 1 2.539 2.640 2.569 2.565 2.668 2.649 2.666 2.745 2.661 2.667 mean value 2.534 f0.135 2.519 f0.109 2.855 2.680 mean value 2.570 k0.162 2.535 f0.107 Ca(3)-0, 2.318 2.330 2.332 2.341 2.309 2.332 2.343 2.383 2.352 2.356 2.469 2.477 2.462 2.57 1 2.501 2.539 2.670 2.636 2.692 2.747 2.683 mean value 2.434 f0.129 2.502 k 0.153 2.470 f0.138 2.760 2.727 ca(4)-0T 2.370 2.371 2.342 2.372 2.333 2.345 2.397 2.386 2.383 2.435 2.43 1 2.405 2.466 2.478 2.45 1 2.510 2.480 2.477 2.794 2.48 1 2.493 mean value 2.477 k0.138 2.424 f0.054 2.412 f0.058 ca(4)-0B P1-0, P3-0, '2-OE p4-0BP,--Pi 2.927 1.637 1.617 1.590 1.616 2.955 1.607 1.61 1 1.612 1.610 2.98 1 1.603 1.610 1.611 1.608 2.973 p3-4mean values 2.991 2.992 2.986 coplanar bonds 4 P-0, 1.497 f0.011 1.495 f0.007 1.491 f0.002 8 P-0, 1.529 f0.016 1.500 f0.003 1.497 k 0.002 angles ploEp2 p3°Bp4mean values 130.5 137.8 135.8 136.5 135.4 136.1 OJ'iOT OTPZ0T 112.6 f3.2 111.8 k 1.3 110.2 f2.5 109.7 k0.6 110.2 f2.7 109.6 k0.2 OTP3OT 111.9 & 2.5 109.6 ,+ 0.9 109.5 & 0.4 OTP4OT OBP2OT OBP3OT OBPiOT OBP~OT 113.9 f 1.2 106.1 & 3.2 106.0 k 1.5 106.8 f 1.5 104.7 k2.2 110.3 k 3.5 108.7 1.2 109.3 & 0.4 109.5 f0.9 108.6 & 1.3 110.2 f 3.7 108.7 f0.9 109.4 k0.3 109.4 f0.4 108.6 f1.3 a From Born-Haber calculation.J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 Table 8 Highest frequencies calculated from the THBPHON program compared with the highest experimental frequency band in the infrared spectrum assigned to v(P-0,) compound calculated frequency/cm- experimental infrared frequencylcm - 1210 1200 1230 1230 1212 oxygens to the central calcium in two sites as determined by X-ray crystallography our model did not allow such a close approach. Most of the lattice parameters decreased slightly on relaxation. The modelled structures were all stable in that there were no imaginary eigenvalues in the THBPHON calculations with a wave vector of q = O.OO0, O.OO0 and 0.001. The frequencies (wavenumbers) derived from the eigenvalues were in the same range as found in the experimental infrared spectra of all the compounds and compared with the range of values calcu- lated in normal coordinate analyses of divalent metal diphos- phate~.~~Because of the limitations of computation at that time these analyses used a model with D,, symmetry despite the fact that the symmetry of the diphosphate anions with a bent bridge can be as low as C,.For these reasons a more detailed analysis does not seem to be justified. The highest frequencies calculated from THBPHON relating to the phos- phorus terminal oxygen stretching vibrations do, however, compare reasonably with experimental values in Table 8. Discussion We have presented our results of modelling different metal diphosphates using the parameters for the non-bonded inter- atomic potentials obtained using the electron gas method- ology and fitted to the Born-Mayer analytic function and dispersion terms calculated by the Slater-Kirkwood method.We have previously used the parameters for the diphosphate unit in our modelling of magnesium diphosphate and have now tested their transferability to model other similar com- pounds. They were used most successfully in modelling the compounds which are structurally isomorphous with magne- sium diphosphate, the manganese and zinc diphosphates. The linear P-0-P bond angle was modelled successfully in the manganese and /%zinc structures. The P-0-P angle has been discussed in the context of the dr-pn overlap and results in the shorter P-0, bonds32 which are found in the manganese and zinc structures with the linear P-0-P bond.The linear bridging oxygen atom did present a slight problem as it is believed that this is an average position hiding the thermal disorder. However, an attempt to model the disordered structure, the split-atom model of manganese diphosphate proposed in an X-ray study was not successful. A double cell was set up to accommodate the extremes of the disorder but this resulted in the appearance of an additional site for manganese in our relaxed structure. Only in the model with a linear bridge was one unique MnO, polyhedron found. The structures with a linear P-0,-P bond in the manganese and the p-zinc compounds needed a modification in the charge distribution in the anion and two ranges for the metal oxygen parameters to reproduce the structure satisfac- torily. The rather more complicated a-zinc diphosphate mod- elled extraordinarily well with the same charge distribution and very minor changes in the parameters used in the magne- sium study.Calcium has a bigger ionic radius than magne- J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 sium, manganese or zinc and so is able to accommodate a higher number of nearest neighbours. In the X-ray structure of /?-calcium diphosphate a judgement was made as to what constituted the nearest neighbours of calcium. On the basis of a distance of under 3 A it was found that the coordination numbers were seven, nine, seven and eight for the four inde- pendent calcium sites.Included in the coordination sphere of Ca (1) and Ca (4) was a bridging oxygen atom at 2.78 and 2.93 A, respectively. All the sites are described as being dis- torted. Sites one and three are described as distorted pentag- onal bipyramids or an octahedron distorted by the approach of a more distant oxygen. Site three is described as a pentag- onal bipyramid with the apical atom replaced by three atoms and the fourth site as an octahedron with two oxygen atoms jammed into edges. Despite this challenging structure the model was reasonably successful especially at the closest neighbours. Conclusions Several metal diphosphates have been modelled using the same parameters for the non-bonded interatomic potentials in the anion as previously used to model magnesium diphos- phate.Phosphorus-oxygen bond distances tended to be slightly short but most were within 2% of the X-ray determi- nations. In the calcium structure the range of phosphorus ter- minal oxygen bond lengths was 1.480 to 1.562 A. Only the unusually long bond distance was not reproduced accurately. Bond angles in the diphosphate units tended to move towards the standard tetrahedral value of 109". The metal- oxygen distances, in general, were in good agreement with the crystallographic data. The bond distances in the manganese and /%zinc diphosphate structures relaxed at constant pres- sure were within 1.35% of the crystallographic values while those in a-zinc diphosphates, which had a more complicated structure, had a difference of 2.89% in one bond.Even in the calcium diphosphate structure most bond distances were within 2% of the values determined by X-ray crystallography. The work has demonstrated that the interatomic potentials used for the static simulation of the perfect lattice of diphos- phates can be used with confidence in various metal com- pounds provided specific potentials are used for the metal- oxygen interactions. Changes in the cell dimensions show that there are decreases in cell volumes on relaxation at constant pressure compared with the experimental values. Most of the experi- mental lattice parameters were obtained at ambient tem-perature except for B-zinc diphosphate when measurements were made on a heated cry~tal'~ (temperature not reported).The potentials have been derived from the electron gas meth- odology using expansions of Slater functions, where there is no compensation for the effects of temperature as there would be for empirically fitted parameters using the CCP5 program THBFIT which maintains cell volumes. Improvements may be obtained by using the shell model33 which takes account of the polarizability of the ions, notably in this case the oxygen which has a more diffuse electron cloud. These studies indicate that these potentials can be used to investigate cation substitutions in different diphosphate compounds and so extend our studies to model biological systems.The potentials are also being evaluated for ortho- phosphate compounds so that we may then model the struc- ture and reactivity of surfaces of the many important phosphatic biominerals such as hydroxyapatite. We acknowledge the help of Dr. A. H. Harker, AEA Tech-nology, Harwell in the calculation of interatomic potentials using the electron gas models. We thank Drs. M. G. B. Drew and P. C. H. Mitchell from the University of Reading, Department of Chemistry for helpful discussions on X-ray crystallography and modelling. This work has been sup- ported by NERC. M.G.T. is grateful to the Leverhulme Trust for funding. References 1 K. Simkiss, M. G.Taylor and G. N. Greaves, J. Inorg. Biochem., 1990,39, 17. 2 G. R.Sauer and R. E. Wuthier, Bone Mineral, 1992,17, 284. 3 G. N. Greaves, K. Simkiss, M. Taylor and N. Binsted, Biochem. J., 1984, 221, 855. 4 M. G. Taylor, K. Simkiss, G. N. Greaves and J. Harries, Proc. R. SOC.London, B, 1988,234,463. 5 M. G. Taylor, G. N. Greaves and K. Simkiss, Eur. J. Biochem., 1990,192,783. 6 G. N. Greaves, S. J. Gurman, C. R. A. Catlow, A. V. Chadwick, S. Houde-Walter, C. M. B. Henderson and B. R. Dobson, Philos. Mag. A, 1991,64, 1059. 7 C. Huang and A. N. Cormack, J. Chem. Phys., 1991,95,3634. 8 B. Vessal, G. N. Greaves, P. T. Marten, A. V. Chadwick, R. Mole and S. Houde-Walter, Nature (London), 1992,356, 504. 9 A. J. de Man, B. W. H. van Beest, M. Leslie and R. A. van Santen, J. Phys. Chem., 1990,94,2524. 10 M. T. Dove, B.Winkler, M. Leslie, M. J. Harris and E. K. H. Salje, Am. Mineral., 1992,77, 244. 11 M. G. Taylor, K. Simkiss, M. G. B. Drew, P. C. H. Mitchell and M. Leslie, Mol. Simul., 1992, 9, 129. 12 G. M. Clark and R. Morley, Chem. SOC.Rev., 1976,5,269. 13 K. Lukaszewicz and R. Smajkiewicz, Rocz. Chem., 1961,35,741. 14 C. Calvo, Can. J. Chem., 1965,43, 1147. 15 T. Stefanidis and A. G. Nord, Acta Crystallogr. Sect. C, 1984, C40, 1995. 16 B. E. Robertson and C. Calvo, J. Solid State Chem., 1970, 1, 120. 17 C. Calvo, Acta Crystallogr., 1967,23, 289. 18 C. Calvo, Inorg. Chem., 1968, 7, 1345. 19 N. C. Webb, Acta Crysrallogr., 1966,21,942 20 A. B. Lidiard, J. Chem. SOC.,Faraday Trans. 2, 1989,85,341. 21 M. Leslie, Physica B (Amsterdam), 1985, 131B, 145. 22 N. F. Mott and M. J. Littleton, Trans. Faraday SOC., 1938, 34, 485. 23 P. R. Adby and M. A. H. Dampster, in Introduction to Opti-mization Methods, Chapman and Hall, London, 1974. 24 A. Hezel and S. D. Ross, Spectrochim. Acta, Part A, 1967, 23, 1583. 25 J. H. Harding and A. H. Harker, Harwell Report, AERE-R10425, 1982. 26 R. G. Gordon and Y. S. Kim, J. Chem. Phys., 1972,56,3122. 27 P.T. Wedepohl, Proc. Phys. SOC.,1967,92, 79. 28 E. Clementi and C. Roetti, At. Nucl. Data Tables, 1974, 14. 29. F. Herman and S. Skillman, in Atomic Structure Calculations, Prentice-Hall, Englewood Cliffs, 1963. 30 J. C. Slater and J. G. Kirkwood, Phys. Rev., 1931,37,682. 31 N. C. Pyper, Philos. Trans. R. SOC. London, Ser. A, 1986, 320, 107. 32 D. W. J. Cruickshank, J. Chem. SOC., 1961, 5486. 33 B. G. Dick Jr. and A. W. Overhauser, Phys. Reu. 1958,112,90. Paper 3/04862D; Received 1lth August, 1993