 Methodology
 Open Access
 Published:
Improved volume variable cluster model method for crystallattice optimization: effect on isotope fractionation factor
Geochemical Transactions volume 23, Article number: 1 (2022)
Abstract
The isotopic fractionation factor and element partition coefficient can be calculated only after the geometric optimization of the molecular clusters is completed. Optimization directly affects the accuracy of some parameters, such as the average bond length, molecular volume, harmonic vibrational frequency, and other thermodynamic parameters. Here, we used the improved volume variable cluster model (VVCM) method to optimize the molecular clusters of a typical oxide, quartz. We documented the average bond length and relative volume change. Finally, we extracted the harmonic vibrational frequencies and calculated the equilibrium fractionation factor of the silicon and oxygen isotopes. Given its performance in geometrical optimization and isotope fractionation factor calculation, we further applied the improved VVCM method to calculate isotope equilibrium fractionation factors of Cd and Zn between the hydroxide (Zn–Al layered double hydroxide), carbonate (cadmiumcontaining calcite) and their aqueous solutions under superficial conditions. We summarized a detailed procedure and used it to reevaluate published theoretical results for cadmiumcontaining hydroxyapatite, emphasizing the relative volume change for all clusters and confirming the optimal point charge arrangement (PCA). The results showed that the average bond length and isotope fractionation factor are consistent with those published in previous studies, and the relative volume changes are considerably lower than the results calculated using the periodic boundary method. Specifically, the average Si–O bond length of quartz was 1.63 Å, and the relative volume change of quartz centered on silicon atoms was − 0.39%. The average Zn–O bond length in the Zn–Allayered double hydroxide was 2.10 Å, with a relative volume change of 1.96%. Cadmiumcontaining calcite had an average Cd–O bond length of 2.28 Å, with a relative volume change of 0.45%. At 298 K, the equilibrium fractionation factors between quartz, Zn–Allayered double hydroxide, cadmiumcontaining calcite, and their corresponding aqueous solutions were \(\Delta ^{30/28} {\text{Si}}_{{{\text{QtzH}}_{4} {\text{SiO}}_{4} }} = 2.20{\permil} \), \(\Delta^{18/16} {\text{O}}_{ {\text{Qtz}}{} ( {\text{H}}_{2} {\text{O}} )_{\text{n}}} = 36.05{\permil}\), \(\Delta^{66/64} {\text{Zn}}_{ {\text{Zn}} {} {\text{Al LDHZn}} ( {\text{H}}_{2} {\text{O}} )_{\text{n}}^{2+}} = 1.12{\permil}\) and \(\Delta^{114/110} {\text{Cd}}_{ {\text{(CdCal)Cd}} ( {\text{H}}_{2} {\text{O}} )_ {\text{n}}^{2 +} } =  0.26{\permil}\) respectively. These results strongly support the reliability of the improved VVCM method for geometric optimization of molecular clusters.
Introduction
In nature, the isotopic composition of an element varies with rock type, source, and age, especially those determined by highresolution mass spectroscopy since the 2000s. This has considerably expanded the usefulness of isotopic tools, which have been widely used to study processes occurring in terrestrial, atmospheric, and aquatic environments [1,2,3,4,5,6,7]. Specifically, the composition of isotopes can be used for geological temperature measurement and geological dating [8,9,10,11], as well as to determine the genesis of ore deposits [12, 13]. Moreover, several important geological events are reflected by the sudden change in isotope composition, such as the largest mass extinction at the Permian–Triassic boundary, revealed by seawater δ^{7}Li and δ^{114}Cd content [14, 15]. Therefore, isotopes also play an important role in tracing paleoclimatic variations [16]. Although isotopes have been widely studied in various fields of geochemistry, the magnitude of isotope fractionation is unknown. This has substantially hindered the development of stable isotope geochemistry. However, the magnitude of isotope fractionation is a key parameter for establishing an isotope evolution model. At present, three methods are used to obtain the isotope equilibrium fractionation factor: experimental determination, empirical estimation, and theoretical calculation. However, laboratory experiments are often conducted under conditions that deviate from natural conditions, such as elevated background concentrations, high temperatures, and short durations. Hence, extrapolation from the laboratoryderived results to geological settings requires further evaluation. The accuracy of an empirical estimation is dependent on the parameters adopted in the calculations, which may result in a large error in the obtained fractionation factor. Therefore, theoretical calculations have become a useful tool for obtaining precise isotope equilibrium fractionation factors. With the rapid increase in storage and computing efficiency in recent decades, supercomputers are now being used to model increasingly large systems of up to tens of thousands of atoms. Computational simulations, particularly quantum chemical simulations, consider various chemical interactions, mimicking key reaction steps. The accuracy of the calculated results is equivalent to that of the experimental analyses. Simulations have proven particularly useful, especially when sampling and/or experimental constraints are difficult to obtain. Among the various simulation methods, firstprinciples calculations based on density functional theory (DFT) have been widely used to calculate isotopic fractionation factors, element partition coefficients, and vibrational spectra. However, obtaining stable geometrical configurations is a prerequisite for calculating these parameters. The harmonic vibrational frequency required to calculate the isotopic fractionation factor is a secondorder partial derivative of the electronic energy of the system with respect to its atomic coordinates.
Research on solid earth science involves the geometric optimization of solid structures. The periodic boundary method is the most commonly used method for simulating solids. Mineral cells are used to construct periodic systems. In general, supercells are required to ensure that the constructed system is comparable to the real mineral lattice, but the total number of unit cells is limited by computational power. In addition, improper handling of atoms at the boundary of the supercells may considerably affect the accuracy of the calculated results. Given the shortcomings of the periodic boundary method, Rustad and Dixon [17], Rustad and Yin [18], and Rustad et al. [19, 20] introduced an embedded cluster model with a Pauling bond strength (PBS)conserving termination of the atoms in their rigid shell. In this model, a large molecular cluster with three shells is employed to represent the minerals. A small core and a second shell are optimized at high and low levels of theory and basis sets, respectively. The third shell is fixed at the measured lattice positions. Chemical bonds entering the third shell are clipped and replaced by “Pauling bond strengthconserving quasiatoms” along the clipped bonds, which make the core charge of each quasiatom equal to the Pauling bond valence, that is, the charge of the removed cation divided by its coordination number. The distance between the quasiatom and the third shell is maintained at 1 Å. Finally, the frequency calculations are performed only for the small core. As an alternative, the volume variable cluster model (VVCM) method was proposed by Liu [21], which was developed during the optimization of silicates, carbonates, oxides, hydroxides, and sulfides [22,23,24,25,26,27,28]. In contrast to the embedded cluster model, in the VVCM model all atoms of the cluster are geometrically optimized, and all or partial atoms are involved in the calculations of the vibrational frequency according to the computational cost. The total number of background point charges amounts to several hundreds, which is much higher than the number of quasiatoms in the embedded cluster model. This ensures the coverage of point charges on the entire molecular cluster. Additionally, the positions of the background point charges are adjusted during the construction of the initial molecular clusters. More importantly, geometric optimization and harmonic vibrational frequency calculations can be performed for both periodic systems (minerals) and nonperiodic systems (solutions) with identical exchange–correlation functionals and basis sets. Although the main two limitations of the molecular cluster method are the same as those of the periodic boundary method (the size of the system and the treatment of atoms on the boundary) and its representation of the periodic mineral lattice is not as accurate as that of the periodic boundary method, it is able to precisely treat the local configuration around the atom of interest, which is the dominant factor controlling isotope fractionation. In recent years, static quantum chemistry calculations based on molecular clusters have accurately calculated the stable isotope equilibrium fractionation factors [22,23,24,25,26,27,28,29,30,31,32]. The calculated results are in agreement with previous experimental and theoretical calculations. However, the critical question is to treat the local configuration around the atom of interest. In previous studies, a detailed procedure to obtain an optimized local configuration was absent. In addition, evaluation proxies for the optimized local configuration are not wellestablished. The energy of a solid (E) is a function of the molecular volume (V) [33], total number of atoms, and interactions among the atoms. Additionally, the relationship between energy and volume is linear and/or nonlinear [34]. For example, the nonlinear fourparameter Vinet equation of state [35, 36] can be expressed as follows:
where the fitting parameter a = E_{0} + 4B_{0}V_{0}/(B_{0}^{’}1)^{2}, V_{0}, B_{0}, and B^{’}_{0} represent the equilibrium volume, bulk modulus, and its first derivative with respect to pressure, respectively; and ΔV is the deviation from the equilibrium volume.
Given the role of molecular volume in the energy of the solid, the expansion or contraction of molecular clusters will influence the accuracy of harmonic vibrational frequencies and the isotope fractionation factor. However, the degree of expansion or contraction of molecular clusters was not considered in the previous molecular cluster modeling method.
In this study, we chose the typical oxide quartz (Qtz) present on the Earth’s surface [37, 38] as an example. By geometrically optimizing the Qtz molecular cluster, extracting its parameters (converged electronic energy, average bond length, and relative volume change), and calculating the equilibrium oxygen and silicon isotopic fractionation factors between molecular clusters and the aqueous solution, we aimed to verify the usefulness of the improved VVCM method after comparison of methods developed in previous studies. Based on the performance of the improved VVCM method for Qtz, we extended the method to the calculation of other nontraditional stable isotopes (Zn and/or Cd) between layered double hydroxide (Zn–Al LDH), carbonate (Cd–containing calcite, Cd–Cal), and their corresponding solutions. After rechecking previous theoretical results for cadmiumcontaining hydroxyapatite (Cd–HAp), emphasizing the relative volume change for all clusters, and confirming the optimal point charge arrangement (PCA), we constructed a detailed procedure to implement the improved VVCM method.
Theory and methods
To obtain an optimal geometric structure, we systematically used the optimization process based on the improved VVCM method, as shown in Fig. 1.
Constructing molecular cluster and adding background point charges
We reconstructed the threedimensional crystal lattices of the minerals from Xray singlecrystal diffraction or neutron diffraction data obtained from previous studies [39,40,41,42,43]. The N–N–N principle must be followed when cutting molecular clusters. This requires the atoms (nearest neighbor) to form chemical bonds with the central atom (atom of interest), and the atoms (the second neighbor atom) forming chemical bonds with the nearest neighbor must be retained. In specific cases, the size of the molecular clusters depends on the computational cost and accuracy. Generally, the outermost layer of a cluster is composed of anions or anionic groups, while the metal ions and/or protons connected to these anions or anionic groups are deleted. The hanging bonds are typically saturated with positive point charges to fix the outermost atoms and to simulate the electronneutral environment of the crystal. Existence of free interlayer anions makes LDH different from others. The general formula of LDH is [M^{2+}_{1−x} M^{3+}_{x} (OH)_{2}]^{x+}(A^{n−})_{x/n}·mH_{2}O, where M^{2+} and M^{3+} represent divalent and trivalent metal cations, respectively; A^{n−} is an interlayer anion including CO_{3}^{2−}, Cl^{−}, SO_{4}^{2−}, NO_{3}^{−}; x represents the molar ratio M^{3+}/(M^{2+} + M^{3+}); and m is the number of water molecules [44, 45]. In nature, the M^{2+}/M^{3+} molar ratio is usually two [46]. Therefore, we focused on a typical Zncontaining LDH with a Zn:Al ratio of two. We constructed two kinds of molecular clusters for LDH without and with interlayer anions.
Oxygen atom can strongly attract the electron of hydrogen atom due to its strong electronegativity. It means that each hydrogen atom of a water molecule can accept a nonshared electron pair of the oxygen atom. Meanwhile, oxygen atom can contribute its two nonshared electron pairs to two hydrogen atoms. Therefore, water molecules can form 4 H bonds with each other in solution [47], which constructs a tetrahedron structure found in ice crystals (Additional file 1: Figure S1). Inspired by this observation, we used point charges and connected chemical bonds to encompass the terminal anions. Different terminal anions have various arrangements of positive point charges. In a previous study, three PCAs were proposed [25], and for the terminated oxygen (η^{1}O) capped by five positive point charges (5 ×), two positive point charges were used as two wings (2 ×) distributed on the bridge oxygen (µ_{2}O) or the terminated hydroxyl (η^{1}OH), and one positive point charge was set to the apex (1 ×) of the bridge oxygen or hydroxyl (µ_{3}O or µ_{2}OH) (Fig. 2). Considering η^{1}O as the predominant anion on the surface of the cluster, we added point charges in four different ways to evaluate their effect on the theoretical results. The one was to add one point charge in the direction of the deleted atoms [17,18,19,20]. In other three ways, three, four, and five point charges were added to η^{1}O, respectively. However, point charges were not fixed at the directions of chipped bonds. Alternatively, these point charges were averagely distributed in the space surrounding the η^{1}O. Obviously, η^{1}O is more tightly capped by five point charges than by three or four point charges (Additional file 1: Figure S1). Furthermore, we performed geometrical optimizations for the Qtz cluster considering 1 ×, 3 × , and 5 × PCAs (Table 1). Although the average Si–O bond length optimized in three PCAs was slightly larger than that in the experimental data (1.61 Å [39] and 1.62 Å [48]), the more accurate relative volume change (ΔV/V_{0}) and \(\Delta^{{{{30} \mathord{\left/ {\vphantom {{30} {28}}} \right. \kern\nulldelimiterspace} {28}}}} {\text{Si}}_{{{\text{Qtz  H}}_{4} {\text{SiO}}_{4} }}\) indicated that 5 × was the best PCA for η^{1}O. Similarly, we safely capped µ_{2}O and µ_{3}O with two and one point charges, respectively, according to a previous study [25].
Before optimizing the configuration, we set the charge quantity and distance of the point charges of the outermost oxygen atoms. The quantity of charge was equal to the amount of charge provided by the outermost anions or anionic groups of the cluster. Specifically, the total amount of point charge on each terminal anion was equal to the sum of the contributions of all “cut” metal ions and/or protons. The distance between the point charge and the outermost oxygen atom was set empirically. Upon adjustment of the distance between the point charge and outermost oxygen atoms, the electronic energy of the system also changed, and each atom of the inner part of the cluster relaxed under its three freedoms (X, Y, and Z dimensions). When the inner atoms amounted to N, 3 N freedoms were observed. The outermost layers of the Qtz clusters centered on silicon or oxygen were all η^{1}O (Fig. 2). In the Qtz cluster, the total electrical charge of the point charges on each η^{1}O was equal to 1/4 of that of the removed Si^{4+}, that is, + 1 (Table 2). Therefore, the partial electrical charge of each point charge was set to + 1/5. In the Zn–Al LDH cluster, the two types of PCAs were 5 × and 1 × (Fig. 2). The total amount of charge on η^{1}O was equal to the total amount of charge contributed by the removed Zn^{2+}, Al^{3+}, and H atoms. Specifically, the total amount of charge was the summation of the charge contributed by one hexacoordinated Zn^{2+} (+ 2/6), one hexacoordinated Al^{3+} (+ 3/6), and one H atom (+ 1). Thus, the total amount of charge on η^{1}O was + 11/6. We set the electric quantity of each point charge to + 11/30. The total charge on µ_{2}OH was equal to the removed charge of one hexacoordinated Al^{3+} (+ 3/6). The electric quantity of each point charge was + 1/2. Furthermore, for Cd–Cal, the two PCA modes were 5 × and 2 × (Fig. 2). The total charge on η^{1}O was equal to the total amount of charge of the two removed hexacoordinated Ca^{2+} (2/6 + 2/6), so the amount of charge divided by each point charge was + (2/6 + 2/6)/5. The total amount of charge on μ_{2}O was + 2/6, and the electric quantity of each point charge was + (2/6)/2.
The configuration of Cdcontaining hydroxyapatite is further complicated by the simultaneous existence of η^{1}O, µ_{2}O, and µ_{3}O in its outermost layer (Fig. 2). Two ways can be used to set the electric quantity for each point charge in the 5 × mode. One fraction of η^{1}O obtains the total charge of two ninecoordinated Ca^{2+} and one sevencoordinated Ca^{2+} (2 × 2/9 + 2/7), setting the electric quantity of each point charge to + (2 × 2/9 + 2/7)/5. The other fraction of η^{1}O obtains the total charge (2 × 2/7 + 2/9) of two sevencoordinated Ca^{2+} and one ninecoordinated Ca^{2+}, with the electric quantity of each point charge set to + (2 × 2/7 + 2/9)/5. The three ways to set the electric quantity for each point charge when adding two points of charge are as follows: (1) the µ_{2}O obtains the total charge of two ninecoordinated Ca^{2+} (2/9 + 2/9) and the electric quantity of each point is set to + (2/9 + 2/9)/2; (2) the µ_{2}O obtains the total charge of two sevencoordinated Ca^{2+} (2/7 + 2/7) and the electric quantity of each point is set to + (2/7 + 2/7)/2; and (3) the µ_{2}O obtains a total charge of a sevencoordinated Ca^{2+} and a ninecoordinated Ca^{2+} (2/7 + 2/9) and the electric quantity of each point is set to + (2/7 + 2/9)/2. The two types of electric quantities for each point charge in the 1 × modes are as follows: the electric quantity of each point charge on μ_{3}O is set to + 2/7 when the sevencoordinated Ca^{2+} is removed; and the point charge on μ_{3}O obtains a charge of + 2/9 when the ninecoordinated Ca^{2+} is removed.
Molecular cluster optimization
We systematically adjusted the distance between the point charge and outermost oxygen atom to change the distribution of the background charge and to change the molecular cluster volume and geometry.
Qtz
Only one PCA exists in the Qtz molecular cluster, that is, 5 × PCA. Here, we considered the Qtz molecular cluster centered on a silicon atom as an example to provide a detailed description. According to previous studies, the range of distances between the 5 × point charge and η^{1}O (d_{5×}) is approximately 100–120 pm [22]. Thus, we adjusted d_{5×} to 100 pm. Each time, we placed the point charge with two sets of distances, with a difference of only 1 pm, for instance, 100 and 101 pm. We determined the direction of further distance adjustment by determining the lower electronic energy (corresponding to 101 pm) in the two groups. Therefore, increasing the distance between the 5 × point charge and η^{1}O decreased the electronic energy of the cluster (Fig. 3). Subsequently, we set the point charge at two other distances (119 and 120 pm) and compared their electronic energies. Assuming that the electronic energy value of each configuration obeys a parabolic (concave upward) distribution, we determined that the point charge distance of the lowest electronic energy configuration should be between 101 and 119 pm. By further comparing the electronic energies of the system at medium distances of 111 and 112 pm, we further constrained the position of the parabola vertex within a narrower range (112–119 pm). Finally, we geometrically optimized all configurations within a small range to determine the optimal point charge distance corresponding to the lowest electronic energy value.
Zn–Al LDH, Cd–Cal, and Cd–HAp
We used the optimal point charge distance search method, similar to that used for Qtz, for evaluating the Zn–Al LDH, Cd–Cal, and Cd–HAp molecular clusters. In the case of the coexistence of multiple PCA modes, we determined the most suitable distance for each PCA mode. Specifically, we arbitrarily fixed the distance between the 1 × point charge and μ_{2}OH or μ_{3}O (d_{1×}) and/or the distance between the 2 × point charge and μ_{2}O (d_{2×}). Subsequently, we successively adjusted d_{5×} until we determined the optimal value. For the Zn–Al LDH and Cd–Cal molecular clusters, we determined the optimal d_{1×} or d_{2×} after fixing d_{5×.} However, confirming the most suitable distance d_{2×} (d_{1×}) in the Cd–HAp molecular cluster requires d_{1×} (d_{2×}) and d_{5×} to be both fixed.
In this study, we used the exchange–correlation functional B3LYP for Qtz and Zn–Al LDH, and BP86 for Cd–Cal and Cd–HAp. We used the allelectronic basis set 6–311G (d) [49, 50] to describe the wavefunction of the Qtz molecular cluster and a mixed basis set to describe the wavefunction of other systems: LanL2DZ for Zn [51,52,53]; LanL2MB for Cd [52]; 6–311 + G (d, p) for Al, O, and H in Zn–Al LDH; and 6–311G (d) for C, Ca, P, O, and H in Cd–Cal and Cd–HAp. We also performed harmonic vibrational frequency analyses to check whether there was an imaginary frequency and to ensure that the most stable configuration corresponded to at least the local minima on the potential electronic energy surface. We calculated electronic energy and performed geometrical optimization and harmonic vibrational frequency analyses using Gaussian09 code [54].
Relative volume change
We calculated the relative volume change (ΔV/V_{0}) to quantify the expansion or contraction of the molecular clusters during geometric optimization. By comparing the relative volume change (relative to the volume of the unoptimized molecular cluster) and the system electronic energy, we screened the configuration with little volume change with respect to the natural mineral lattice and the lowest electronic energy. Only this method could ensure that the simulation matched the real scenario to the largest extent. Because a single molecular cluster is not a unit cell or supercell of a mineral, accurately calculating the molecular cluster volume using cell parameters is impossible. The wavefunction or electronic density of a molecule is a function of space points. Moreover, the Monte Carlo method can be used to calculate the molecular cluster volume (molecular van der Waals volume), which is the space enclosed by an isosurface with an electron density of x e/Bohr^{3}, where x represents the electron density value of the isosurface. We set up a rectangular box with volume V for molecular volume integration. Subsequently, we placed many random points (m) in the box for sampling and tested the electron densities at these space points to determine if they were larger than x e/Bohr^{3}. Finally, we counted the number of points (n) with electron densities larger than x e/Bohr^{3}, for which the van der Waals volume of the molecule was n/m*V [55, 56]. Bader et al. [57] proposed a more accurate definition of the van der Waals volume, stating that if a molecule is in the gas phase, an isosurface with an electron density of 0.001 e/Bohr^{3} is considered a van der Waals surface. Such surfaces typically contain more than 98% the electron density of a molecule. For molecules in a condensed state, intermolecular extrusion and various forms of interaction cause the van der Waals surface to penetrate. Using an isosurface with an electron density of 0.002 e/Bohr^{3} as the van der Waals surface is generally recommended. Among the molecular clusters constructed in this study, Cd–HAp had the largest number of atoms (no more than 169 atoms). Hence, we treated the molecular clusters as small molecules in the gas phase and used an isosurface with an electron density of 0.001 e/Bohr^{3} as the van der Waals surface of the molecular clusters. Cd–HAp was studied by He et al. [25], but they did not select the optimal PCA, by strictly adhering to our procedure. We rechecked the optimal PCA for this mineral, emphasizing the relative volume changes for all clusters.
Calculating the van der Waals volume requires the molecular wavefunctions. We obtained the wavefunction files using Gaussian09 code [54]. By inputting i, x, and k values, we calculated the volume using the Multiwfn program [58], where i is a proxy for the number of random points distributed in the box and k is the space reserved by the box around the molecule. Alternatively, the shortest distance between the outermost atom and the boundary of the box is k times the van der Waals radius of the atom. For specific calculations, we distributed 100 × 2^{i} random points in a box and set the value of k to 1.7. We used an isosurface with an electron density of 0.001 e/Bohr^{3} to envelop the molecule. Normally, to obtain a more accurate molecular van der Waals volume, the density of random points in the box must not be less than 2000 points/Bohr^{3}. We used the following formula for these calculations:
where V_{opt} denotes the optimized volume and V_{0} is the unoptimized volume. A positive ΔV/V_{0} indicates expansion of the molecular cluster, whereas a negative ΔV/V_{0} suggests contraction of the molecular cluster.
Isotopic equilibrium fractionation factor
The Bigeleisen–Mayer equation and Urey model [59, 60] provide the equilibrium constants for the isotopic exchange reaction, laying the foundation for the theory and calculation of stable isotope geochemistry. Consider the following isotope exchange reaction as an example:
where AX and AXʹ represent compounds with heavy and light isotopes, respectively; X and Xʹ denote the ideal singleatom gas with heavy and light isotopes, respectively. The reaction equilibrium constant K can be obtained as follows from the ratio of the partition functions of the products and reactants:
where Q_{trans} is the transitional partition function, Q_{rot} is the rotational partition function, Q_{vib} is the vibrational partition function, and Q_{elec} is the electronic partition function [61]. For light elements, the difference in the groundstate electronic energy of the compound before and after isotope substitution is indistinguishable, so Q_{elec} can be ignored. The translational, rotational, and vibrational partition functions can be respectively written as follows:
where V is the volume, M is the mass, I_{A} is the moment of inertia around the axis A of rotation, σ is the symmetry number of the molecule, and u_{i} = hcω_{i}/kT, where h represents the Planck constant, k is the Boltzmann constant, T is the temperature in degrees Kelvin, c is the speed of light, and ω is the harmonic frequency in cm^{−1}. The TellerRedlich [62] approximation is done as follows:
Therefore, the ratio of the translational partition function and the ratio of the rotational partition function can be expressed by the vibrational frequency. Through algorism transformation, the reduced partition function ratio (RPFR(AX/AX^{’})) is defined as follows:
When only one atom is to be exchanged in the molecular cluster, the RPFR is equivalent to the β factor [61], which is the isotopic fractionation factor between a compound and the ideal atomic gas [63]. Given the β factors of a pair of compounds or minerals in equilibrium, the isotope fractionation factor α can be derived from the ratio of their β factors: α = β_{A}/β_{B}, where A and B denote different compounds or minerals. The isotopic equilibrium fractionation between the two phases can then be obtained as follows:
Results
Influence of d_{5×}, d_{2×}, and d_{1×} on system electronic energy
We collected the electronic energy values of the four systems and performed polynomial fitting. The results showed that the system electronic energy had a good parabolic relationship with d_{5×}, d_{2×}, and d_{1×}, and the R^{2} value ranged from 0.60 to 0.99 (Figs. 3, 4, 5, 6). The PCA modes corresponding to the most stable configurations of the molecular clusters were 5 × 116, 5 × 116, 5 × 111/1 × 76, 5 × 123/2 × 122, and 5 × 123/2 × 119/1 × 115 pm for Si–Qtz, O–Qtz, Zn–Al LDH, Cd–Cal, and Cd–HAp, respectively.
Effect of d_{5×}, d_{2×}, and d_{1×} on relative volume change
During the geometric optimization process, the volume of the Si–Qtz molecular cluster decreased, with ΔV/V_{0} ranging from − 2.05 to − 0.11% (Fig. 7). We observed an extremely strong linear correlation between ΔV/V_{0} and d_{5×} (R^{2} = 0.9895). ΔV/V_{0} gradually increased with an increase in d_{5×}. We recorded the ΔV/V_{0} value as − 0.39% when the Si–Qtz molecular cluster reached the optimal structure. However, the volume of the Zn–Al LDH molecular cluster without interlayer anions expanded during the optimization process, with ΔV/V_{0} ranging from 3.17 to 3.54% (Fig. 8). The relative volume change was 3.54% when the system reached its lowest electronic energy configuration (5 × 111/1 × 76 pm). But when interlayer anions were introduced, the relative volume change reached 1.96% for the lowest electronic energy configuration (5 × 111/1 × 76 pm). As shown in Fig. 9, the volume of Cd–Cal expanded during the optimization process ranging from 0.24 to 0.50%. When the structure reached the maximum stability (5 × 123/2 × 122 pm), the volume expansion was 0.45% (Table 4). The volume of the Cd–HAp molecular clusters also expanded during the optimization process, with a volume expansion of 0.64–0.92%. We found that ΔV/V_{0} was 0.74% when the system reached the most stable configuration (5 × 123/2 × 119/1 × 115 pm) (Fig. 10), which is identical to the 0.7% reported by He et al. [25].
Discussion
Geometrically optimized configuration
In this study, we obtained an average Si–O bond length in the Qtz cluster with the lowest electronic energy of 1.63 Å (Table 3). Purton et al. [64] used ab initio Local Density Functional (LDF) theory to study the structure of α–quartz, and obtained an average Si–O bond length of 1.62 Å. He and Liu [22] calculated an average Si–O bond length of 1.61 Å on the 6–311G (2df) level. Hazen et al. [39] and Kihara [48] obtained the average Si–O bond lengths of 1.61 Å and 1.62 Å, respectively, using Xray diffraction. Our results are consistent with those of previous theoretical calculations and experimental studies.
The central Zn^{2+} in the Zn–Al LDH cluster had a coordination number of six, with an average Zn–O bond length of 2.10 Å. Gou et al. [3], by Xray diffraction and EXAFS experiments, obtained an average Zn–O bond length of 2.08 Å (coordination number = 6.5) in the Zn–Al layerlike double hydroxides. Thus, the findings of these two studies are in good agreement with each other. For Cd–Cal, the average Cd–O (2.28 Å) bond length calculated in this study is consistent with the results obtained using EXAFS data (2.30 Å), as measured by Bailey et al. [65]. According to our procedure, the average Cd–O bond length in the Cd–HAp cluster was confirmed to be 2.39 Å, which is exactly the same as the average Cd–O (2.39 Å) bond length measured by Hata et al. [66] and the average Cd–O (2.38 Å) bond length calculated by He et al. [25].
Relative volume changes
The relative volume change of mineral lattices during geometrical optimization has only been studied by a few researchers. Wang et al. [4] calculated the ΔV/V_{0} values of Bacontaining species by using the periodic boundary method. Except for Ba(OH)_{2}(H_{2}O)_{8}, the relative volume change of which was 0.67%, the ΔV/V_{0} of the other species ranged from 1.26% to 6.30%, with an average of 4.46%. Ducher et al. [67] also used periodic boundary method to calculate the ΔV/V_{0} of various Zncontaining minerals (sulfide, carbonate, oxide, silicate, sulfate, and arsenate), which ranged between 1.40% and 5.80%, with an average of 3.56%. Thus, we concluded that for the majority of the studied systems, the ΔV/V_{0} values calculated using the molecular cluster method are much smaller (Table 4), approximately onetenth of those obtained by the periodic boundary method.
We found that the ΔV/V_{0} of quartz using the periodic boundary method calculated by the PBE and LDA functionals, primarily used by Méheut et al. [68], were 7.17% and 2.66%, respectively. Hamann [69] obtained the ΔV/V_{0} of quartz using GGA (4.61%, PW91) and LDA (2.75%). The ΔV/V_{0} obtained in this study was 0.39%, which is substantially less than that obtained in previous studies. The Zn–Al LDH has a slightly larger expansion (1.96%) during geometric optimization, compared to that of other minerals in our study. But it is still smaller than those obtained by the periodic boundary method. The volume of the calculated Cdcontaining calcite also increased (0.45%). However, the relative volume change of the Mncontaining calcite (3.54%) [70] and tetrahedral AO_{4}^{2−} group (A = S, Cr, Se)doped calcite (3.47%) [71] calculated by the periodic boundary method were larger than that obtained during this study. For HAp without impurities, Sailuam et al. [72] calculated a ΔV/V_{0} of 3.43% using GGA–PBE as the exchange–correlation functional. Bhat et al. [73] obtained ΔV/V_{0} (4.37%) using the exchange–correlation functional GGA–PW91. In contrast, the ΔV/V_{0} of Cd–HAp (Cd/(Cd + Ca) = 1/13) calculated in this study was 0.74%, which is much smaller than the previously reported values.
\(\Delta^{{{{30} \mathord{\left/ {\vphantom {{30} {28}}} \right. \kern\nulldelimiterspace} {28}}}} {\text{Si}}_{{{\text{Qtz  H}}_{4} {\text{SiO}}_{4} }}\) and \(\Delta^{{{{18} \mathord{\left/ {\vphantom {{18} {16}}} \right. \kern\nulldelimiterspace} {16}}}} {\text{O}}_{{{\text{Qtz  }}\left( {{\text{H}}_{2} {\text{O}}} \right)_{{\text{n}}} }}\)
With different exchange–correlation functionals and basis sets, 1000lnβ may be different. At 298 K, Méheut et al. [74] calculated the silicon βfactor of quartz as 69‰ and the oxygen βfactor as 102‰. At the same temperature, our silicon β–factor (73.65‰) and oxygen βfactors (110.31‰) were slightly larger than those reported by Méheut et al. [74]. However, when calculating the fractionation factor α, the difference in the βfactor for the same mineral can be offset. No substantial difference was observed between the \(\Delta^{{{{18} \mathord{\left/ {\vphantom {{18} {16}}} \right. \kern\nulldelimiterspace} {16}}}} {\text{O}}_{{{\text{quartz  }}\left( {{\text{H}}_{2} {\text{O}}} \right)}}\) value reported by Méheut et al. [68] and that obtained in this study (31.37‰ vs. 36.05‰). We compared the silicon and oxygen fractionation factors calculated in this study with previous experimental and theoretical data.
He and Liu [22] calculated the Si isotopic equilibrium fractionation between quartz and H_{4}SiO_{4} solution at 273–673 K using B3LYP/6–311G(2df). They thought that a large Si isotopic equilibrium fractionation factor of approximately 3.3‰ exists between quartz and H_{4}SiO_{4} solution at 298 K. Dupuis et al. [75] calculated the fractionation factors of Si isotopes between quartz and aqueous H_{4}SiO_{4} at 273–323 K using firstprinciples methods, and found that at 300 K, the fractionation between quartz and H_{4}SiO_{4} was + 2.1 ± 0.2‰. Figure 11 shows that our calculation results are in good agreement with those of Dupuis et al. [75]. Under neutral and acidic conditions, Si often exists in the form of aqueous H_{4}SiO_{4} solutions. We considered the experimentally determined ∆^{30/28}Si_{quartzaqueous solution} to be \( \Delta^{30/28}{\text{Si}}_{{\text {quartz}}\text{H}_{4}{\text{SiO}}_{4}}\). From this perspective, the available experimental data [5, 7, 76, 77] match the existing theoretical data.
To verify the reliability of the VVCM method, we further calculated the equilibrium oxygen isotope fractionation factor between quartz and the aqueous solution and then compared our result with previous experimental and theoretical results (Fig. 12). In the lowtemperature range, compared to the theoretical results of the linear fitting reported by Méheut et al. [68], our calculated results are more consistent with the experimental values [6, 83] and the fitting equation based on experimental data [84]. However, at high temperatures, both theoretical results are less than the experimental value [85,86,87,88,89,90,91]. For example, at 673 K, our calculated oxygen isotope equilibrium fractionation factor and the experimental values reported by Méheut et al. [68] and Sharp et al. [84] were 1.79‰, 2.77‰, and 4.25‰, respectively.
In conclusion, the quartz geometric optimization and calculation of traditional/nontraditional stable isotope fractionation factors was accurate using the improved VVCM method. Finally, we can reliably apply this method to other systems such as Zn–Al LDH and Cd–Cal.
\(\Delta^{{{{66} \mathord{\left/ {\vphantom {{66} {64}}} \right. \kern\nulldelimiterspace} {64}}}} {\text{Zn}}_{{{\text{Zn}}  {\text{Al LDH  Zn}}\left( {{\text{H}}_{2} {\text{O}}} \right)_{{\text{n}}}^{2 + } }}\) and \(\Delta^{{{{114} \mathord{\left/ {\vphantom {{114} {110}}} \right. \kern\nulldelimiterspace} {110}}}} {\text{Cd}}_{{\left( {{\text{Cd}}  {\text{Cal}}} \right)  {\text{Cd}}\left( {{\text{H}}_{2} {\text{O}}} \right)_{\text{n}}^{2 +}}}\)
Because no experimental or theoretical data on Zn isotope fractionation between Zn–Al LDH and the solution have been published, we chose the Zn isotope fractionation produced during Zn^{2+} adsorption on the aluminum oxide/hydroxide surface for comparison. As reported by Komárek et al. [92] and Gou et al. [93], at higher surface coverage, Zn–Al LDH was detected by XRD on the surface of γAl_{2}O_{3} with a fractionation factor of ∆^{66/64}Zn_{solidaq} = 0.02 ± 0.07‰. Pokrovsky et al. [94] also reported lower fractionation on corundum (∆^{66/64}Zn_{solidaq} = 0.19 ± 0.05‰) and gibbsite (∆^{66/64}Zn_{solidaq} = 0.10 ± 0.05‰) surfaces, suggesting that the potential formation of Zn–Al LDH on the surface of gibbsite and corundum is similar to that of γAl_{2}O_{3.} Although the size of calculated isotope fractionation (1.12‰) was larger than that obtained experimentally, the fractionation directions were consistent with each other. In fact, experimentally determined Zn isotope fractionation factor for adsorption on aluminum oxides/hydroxides may be complicated by mixed adsorption and precipitation. In future, we hope that Zn isotope fractionation can be experimentally determined for the individual Zn–Al LDH precipitation.
The calculated isotopic equilibrium fractionation between Cd–Cal and Cdcontaining aqueous solution was 0.26‰. Previous studies on precipitated calcite from artificial seawater solutions reported the entry of Cd^{2+} into the calcite from seawater, with the Cd equilibrium isotope fractionation, ∆^{114/110}Cd_{calciteseawater} = 0.45 ± 0.12‰ [95]. Furthermore, Xie et al. [96] revealed apparent Cd isotope fractionation in the coprecipitation of Cd–calcite, with ∆^{114/110}Cd _{calciteaqueous} = 0.38 ± 0.18‰ to 0.20 ± 0.12‰. Our theoretical calculation results are consistent with previous experimental results.
Conclusions
We used the firstprinciples calculation method to geometrically optimize the molecular clusters of Qtz, Zn–Al LDH, and Cd–Cal and obtained relative volume changes of 0.39%, 1.96%, and 0.45%, respectively. Compared to the periodic boundary method, the improved VVCM method slightly changes the molecular cluster volume during optimization, which increased the accuracy of our simulations. We also calculated the isotopic equilibrium fractionation between Qtz, Zn–Al LDH, Cd–Cal, and their corresponding aqueous solutions at 298 K, which were as follows: \(\Delta^{{{{30} \mathord{\left/ {\vphantom {{30} {28}}} \right. \kern\nulldelimiterspace} {28}}}} {\text{Si}}_{{{\text{Qtz  H}}_{4} {\text{SiO}}_{4} }} = 2.20{\permil}\), \(\Delta^{{{{18} \mathord{\left/ {\vphantom {{18} {16}}} \right. \kern\nulldelimiterspace} {16}}}} {\text{O}}_{{{\text{Qtz  }}\left( {{\text{H}}_{2} {\text{O}}} \right)_{{\text{n}}} }} = 36.05{\permil}\), \(\Delta^{{{{66} \mathord{\left/ {\vphantom {{66} {64}}} \right. \kern\nulldelimiterspace} {64}}}} {\text{Zn}}_{{{\text{Zn}}  {\text{Al LDH  Zn}}\left( {{\text{H}}_{2} {\text{O}}} \right)_{{\text{n}}}^{2 + } }} = 1.12{\permil}\)and \(\Delta^{{{{114} \mathord{\left/ {\vphantom {{114} {110}}} \right. \kern\nulldelimiterspace} {110}}}} {\text{Cd}}_{{\left( {{\text{Cd}}  {\text{Cal}}} \right)  {\text{Cd}}\left( {{\text{H}}_{2} {\text{O}}} \right)_{{\text{n}}}^{2 + } }} =  0.26{\permil}\). These results are consistent with those of previous studies, confirming the reliability of the improved VVCM method for crystallattice optimization.
Availability of data and materials
The data sets that support the conclusions of this study are included in the study.
Abbreviations
 VVCM:

Volume variable cluster model
 DFT:

Density Functional Theory
 Si–Qtz:

Quartz centered on silicon atoms
 O–Qtz:

Quartz centered on oxygen atoms
 LDH:

Layered double hydroxide
 Cd–Cal:

Cd–containing calcite
 Cd–HAp:

Cd–containing hydroxyapatite
 PCA:

Point charge arrangement
References
Zhu XC, Di D–R, Ma Ming–G, Shi Wei–Y, (2019) Stable isotopes in greenhouse gases from soil: a review of theory and application. Atmosphere 10:377
Erkaev NV, Scherf M, Thaller SE, Lammer H, Mezentsev AV, Ivanov VA, Mandt KE (2020) Escape and evolution of Titan’s N_{2} atmosphere constrained by ^{14}N/^{15}N isotope ratios. Mon Not R Astron Soc 500:2020–2035
Gou WX, Siebecker MG, Wang ZM, Li W (2018) Competitive sorption of Ni and Zn at the aluminum oxide/water interface: an XAFS study. Geochem Trans 19:9
Wang WZ, Wu ZQ, Huang F (2021) Equilibrium barium isotope fractionation between minerals and aqueous solution from first–principles calculations. Geochim Cosmochim Acta 292:64–77
Wang W, Geilert S, Wei HZ, Jiang SY (2021) Competition of equilibrium and kinetic silicon isotope fractionation during silica precipitation from acidic to alkaline pH solutions in geothermal systems. Geochim Cosmochim Acta 306:44–62
Labeyrie LJ (1974) New approach to surface seawater paleotemperatures using^{18}O/^{16}O ratios in silica of diatom frustules. Nature 248:40–42
Dai S, Li D, Ren D, Tang Y, Shao L, Song H (2004) Geochemistry of the late Permian No. 30 coal seam, Zhijin Coalfield of Southwest China: influence of a siliceous lowtemperature hydrothermal fluid. Appl Geochem 19:1315–1330
Clayton RN, Epstein S (1961) The use of oxygen isotopes in high–temperature geological thermometry. J Geol 69:447–452
Klein RT, Lohmann KC, Thayer CW (1996) Bivalve skeletons record seasurface temperature and δ^{18}O via Mg/Ca and ^{18}O/^{16}O ratios. Geology 24:415
Mao JW, Zhang ZC, Zhang ZH, Du A (1999) Re–Os isotopic dating of molybdenites in the Xiaoliugou W (Mo) deposit in the northern Qilian mountains and its geological significance. Geochim Cosmochim Acta 63:1815–1818
Wan Y, Liu D, Wilde SA, Cao J, Chen B, Dong C, Song B, Du L (2010) Evolution of the Yunkai terrane, south China: evidence from shrimp zircon UPb dating, geochemistry and Nd isotope. J Asian Earth SCI 37:140–153
Zhang Z, Zuo R (2014) Sr–Nd–Pb isotope systematics of magnetite: implications for the genesis of Makeng Fe deposit, southern China. Ore Geol Rev 57:53–60
Ma YB, Bagas L, Xing SW, Zhang ST, Wang RJ, Li N, Zhang ZJ, Zou YF, Yang XQ, Wang Y (2016) Genesis of the stratiform Zhenzigou Pb–Zn deposit in the North China Craton: Rb–Sr and CO–S–Pb isotope constraints. Ore Geol Rev 79:88–104
Zhang Y, Wen H, Zhu C, Fan H, Cloquet C (2018) Cadmium isotopic evidence for the evolution of marine primary productivity and the biological extinction event during the PermianTriassic crisis from the Meishan section, South China. Chem Geol 481:110–118
Sun H, Xiao Y, Gao Y, Zhang G, Casey JF, Shen Y (2018) Rapid enhancement of chemical weathering recorded by extremely light seawater lithium isotopes at the PermianTriassic boundary. Proc Natl Acad of Sci 115:3782–3787
Eiler JM (2011) Paleoclimate reconstruction using carbonate clumped isotope thermometry. Quat Sci Rev 30:3575–3588
Rustad JR, Dixon DA (2009) Prediction of ironisotope fractionation between hematite (α–Fe_{2}O_{3}) and ferric and ferrous iron in aqueous solution from density functional theory. J Phys Chem A 113:12249–12255
Rustad JR, Yin QZ (2009) Iron isotope fractionation in the Earth’s lower mantle. Nat Geosci 2:514–518
Rustad JR, Casey WH, Yin QY, Bylaska E, Felmy AR, Bogatko SA, Jackson VE, Dixon DA (2010) Isotopic fractionation of Mg^{2+}(aq), Ca^{2+}(aq), and Fe^{2+}(aq) with carbonate minerals. Geochim Cosmochim Acta 74:6301–6323
Rustad JR, Nelmes SL, Jackson VE, Dixon DA (2008) Quantumchemical calculations of carbon–isotope fractionation in CO_{2}(g), aqueous carbonate species, and carbonate minerals. J Phys Chem A 112:542–555
Liu Y (2013) On the test of a new volume variable cluster model method for stable isotopic fractionation of solids: equilibrium Mg isotopic fractionations between minerals and solutions. In: Abstracts of the Goldschmidt conference (Vol. 1632)
He HT, Liu Y (2015) Si isotope fractionations during the precipitation of quartz and the adsorption of H_{4}SiO_{4}(aq) on Fe(III)–oxyhydroxide surfaces. Chin J Geochem 34:459–468
He HT, Zhang ST, Zhu C, Liu Y (2016) Equilibrium and kinetic Si isotope fractionation factors and their implications for Si isotope distributions in the Earth’s surface environments. Acta Geochim 35:15–24
Gao CH, Cao XB, Liu Q, Yang YH, Zhang ST, He YY, Tang M, Liu Y (2018) Theoretical calculation of equilibrium Mg isotope fractionations between minerals and aqueous solutions. Chem Geol 488:62–75
He HT, Xing L–C, Qin S–J, Ji X–Y, Sun P–F, (2020) Equilibrium Cd isotopic fractionation between Cd(OH)_{2}(s), apatite, adsorbed Cd^{2+}, and Cd^{2+}(aq): Potential application of δ^{114}Cd in evaluating the effectiveness of Cd–contamination remediation. Geochem J 54:289–297
Zhang JX, Liu Y (2018) Zinc isotope fractionation under vaporization processes and in aqueous solutions. Acta Geochim 37:27–39
Li X, Liu Y (2015) A theoretical model of isotopic fractionation by thermal diffusion and its implementation on silicate melts. Geochim Cosmochim Acta 154:18–27
Zhang J (2021) Equilibrium sulfur isotope fractionations of several important sulfides. Geochem J 55:135–147
Fujii T, Moynier F, Albarède BT, F, (2014) Density functional theory estimation of isotope fractionation of Fe, Ni, Cu, and Zn among species relevant to geochemical and biological environments. Geochim Cosmochim Acta 140:553–576
Li XF, Zhao H, Tang M, Liu Y (2009) Theoretical prediction for several important equilibrium Ge isotope fractionation factors and geological implications. Earth Planet Sci Lett 287:1–11
Wu F, Qin T, Li XF, Liu Y, Huang JH, Wu ZQ, Huang F (2015) First–principles investigation of vanadium isotope fractionation in solution and during adsorption. Earth Planet Sci Lett 426:216–224
Ducher M, Blanchard M, Balan E (2018) Equilibrium isotopic fractionation between aqueous Zn and minerals from first–principles calculations. Chem Geol 483:342–350
Blanco MA, Francisco E, Luana V (2004) GIBBS: isothermalisobaric thermodynamics of solids from energy curves using a quasiharmonic Debye model. Comput Phys Commun 158:57–72
Shang SL, Wang Y, Kim D, Liu ZK (2010) First–principles thermodynamics from phonon and Debye model: application to Ni and Ni_{3}Al. Comput Mater Sci 47:1040–1048
Vinet P, Ferrante J, Smith JR, Rose JH (1986) A universal equation of state for solids. J Phys C 19:L467–L473
Vinet P, Smith JR, Ferrante J, Rose JH (1987) Temperature effects on the universal equation of state of solids. Phys Rev B 35:1945–1953
Liu WT, Shen YR (2008) Surface Vibrational Modes of α–Quartz (0001) probed by sum–frequency spectroscopy. Phys Rev Lett 101:016101
Tang C, Zhu J, Zhou Q, Wei J, Zhu R, He H (2014) Surface heterogeneity of SiO_{2} polymorphs: an XPS investigation of α–Quartz and α–Cristobalite. J Phys Chem C 118:26249–26257
Hazen RM, Finger LW, Hemley RJ, Mao HK (1989) High–pressure crystal chemistry and amorphization of α–quartz. Solid State Commun 72:507–511
Paquette J, Reeder RJ (1990) Single–crystal X–ray structure refinements of two biogenic magnesian calcite crystals. Am Mineral 75:1151–1158
Wilson RM, Elliot JC, Dowker SEP (1999) Rietveld refinement of the crystallographic structure of human dental enamel apatites. Am Mineral 84:1406–1414
Fleet ME, Liu XY, Pan YM (2000) Rare–earth elements in chlorapatite [Ca_{10}(PO_{4})_{6}Cl_{2}]: uptake, site preference, and degradation of monoclinic structure. Am Mineral 85:1437–1446
Zhu SD, Khan MA, Wang FY, Bano Z, Xia MZ (2020) Rapid removal of toxic metals Cu^{2+} and Pb^{2+} by amino trimethylene phosphonic acid intercalated layered double hydroxide: a combined experimental and DFT study. Chem Eng J 392:123711
Xiang YL, Xiang YX, Jiao YR, Wang LP (2019) Surfactant–modified magnetic CaFe–layered double hydroxide for improving enzymatic saccharification and ethanol production of Artemisia ordosica. Renew Energ 138:465–473
Yang F, Zhang SS, Sun YQ, Tsang DCW, Cheng K, Ok YS (2019) Assembling biochar with various layered double hydroxides for enhancement of phosphorus recovery. J Hazard Mater 365:665–673
Costa DG, Rocha AB, Diniz R, Souza WF, Chiaro SSX, Leitão AA (2010) Structural model proposition and thermodynamic and vibrational analysis of hydrotalcitelike compounds by DFT calculations. J Phys Chem C 114:14133–14140
Ignatov I, Mosin O (2014) Mathematical models describing water clusters as interaction among water molecules. distributions of energies of hydrogen bonds. J Med Physiol Biophys 3:48–70
Kihara K (1990) An X–ray study of the temperature dependence of the quartz structure. Eur J Mineral 2:63–77
McLean AD, Chandler GS (1980) Contracted Gaussian basis sets for molecular calculations. I. Second row atoms, Z=11–18. J Chem Phys 72:5639–5648
Blaudeau JP, McGrath MP, Curtiss LA, Radom L (1997) Extension of Gaussian–2 (G2) theory to molecules containing third–row atoms K and Ca. J Chem Phys 107:5016–5021
Hay PJ, Wadt WR (1985) Ab initio effective core potentials for molecular calculations. Potentials for the transition metal atoms Sc to Hg. J Chem Phys 82:270–283
Hay PJ, Wadt WR (1985) Ab initio effective core potentials for molecular calculations–potentials for K to Au including the outermost core orbitals. J Chem Phys 82:299–310
Wadt WR, Hay PJ (1985) Ab initio effective core potentials for molecular calculations. Potentials for main group elements Na to Bi. J Chem Phys 82:284–298
Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Scalmani G, Barone V, Mennucci B, Petersson GA, Nakatsuji H, Caricato M, Li X, Hratchian HP, Izmaylov AF, Bloino J, Zheng G, Sonnenberg JL, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Vreven T, Montgomery JA, Peralta JE Jr, Ogliaro F, Bearpark M, Heyd JJ, Brothers E, Kudin KN, Staroverov VN, Keith T, Kobayashi R, Normand J, Raghavachari K, Rendell A, Burant JC, Iyengar SS, Tomasi J, Cossi M, Rega N, Millam JM, Klene M, Knox JE, Cross JB, Bakken V, Adamo C, Jaramillo J, Comperts R, Stratmann RE, Yazyev O, Austin AJ, Cammi R, Pomelli C, Ochterski JW, Martin RL, Morokuma K, Zakrzewski VG, Voth GA, Salvador P, Dannenberg JJ, Dapprich S, Daniels AD, Farkas O, Foresman JB, Ortiz JV, Cioslowski J, Fox DJ (2010) Gaussian 09 (Revision C01). Gaussian Inc Wallingford.
Wong MW, Wiberg KB, Frisch MJ (1995) Ab initio calculation of molar volumes: comparison with experiment and use in solvation models. J Comput Chem 16:385–394
Till MS, Ullmann GM (2010) McVol–A program for calculating protein volumes and identifying cavities by a Monte Carlo algorithm. J Mol Model 16:419–429
Bader RFW, Carroll MT, Cheeseman JR, Chang C (1987) Properties of atoms in molecules: atomic volumes. J Am Chem Soc 109:7968–7979
Lu T, Chen F (2012) Multiwfn: a multifunctional wavefunction analyzer. J Comput Chem 33:580–592
Bigeleisen J, Mayer MG (1947) Calculation of equilibrium constants for isotopic exchange reactions. J Chem Phys 15:261–267
Urey HC (1947) The thermodynamic properties of isotopic substances. J Chem Soc 85:562–581
Liu Q, Tossell JA, Liu Y (2010) On the proper use of the Bigeleisen–Mayer equation and corrections to it in the calculation of isotopic fractionation equilibrium constants. Geochim Cosmochim Acta 74:6965–6983
Redlich O (1935) Eine allgemeine beziehung zwischen den schwingungsfrequenzen isotoper molekeln. Z Phys Chem B 28:371–382
Richet P, Bottinga Y, Javoy M (1977) Review of hydrogen, carbon, nitrogen, oxygen, sulfur, and chlorine stable isotope fractionation among gaseous molecules. Annu Rev Earth Planet Sci 5:65–110
Purton J, Jones R, Heggie M, Öberg S, Catlow CRA (1992) LDF pseudopotential calculations of the αquartz structure and hydrogarnet defect. Phys Chem Miner 18:389–392
Bailey EH, Mosselmans JFW, Young SD (2005) Time–dependent surface reactivity of Cd sorbed on calcite, hydroxylapatite and humic acid. Mineral Mag 69:563–575
Hata M, Okada K, Iwai S, Akao M, Aoki H (1978) Cadmium hydroxyapatite. Acta Crystallogr Sect B: Struct Sci 34:3062–3064
Ducher M, Blanchard M, Balan E (2016) Equilibrium zinc isotope fractionation in Zn–bearing minerals from first–principles calculations. Chem Geol 443:87–96
Méheut M, Lazzeri M, Balan E, Mauri F (2007) Equilibrium isotopic fractionation in the kaolinite, quartz, water system: prediction from first–principles densityfunctional theory. Geochim Cosmochim Acta 71:3170–3181
Hamann DR (1996) Generalized gradient theory for silica phase transitions. Phys Rev Lett 76:660–663
Son SB, Newton AG, Jo KN, Lee JY, Kwon KD (2019) Manganese speciation in Mn–rich CaCO_{3}: a density functional theory study. Geochim Cosmochim Acta 248:231–241
Dompablo ME, Fernández MA, Fernández L (2015) Computational investigation of the influence of tetrahedral oxoanions (sulphate, selenate and chromate) on the stability of calcium carbonate polymorphs. RSC Adv 5:59845–59852
Sailuam W, Phacheerak K, Bootchanont A, Fongkaew I, Limpijumnong S (2020) Elastic and mechanical properties of hydroxyapatite under pressure: a first–principles investigation. Compu Condens Matter 24:e00481
Bhat SS, Waghmare UV, Ramamurty U (2014) Firstprinciples study of structure, vibrational, and elastic properties of stoichiometric and calciumdeficient hydroxyapatite. Cryst Growth Des 14:3131–3141
Méheut M, Lazzeri M, Balan E, Mauri F (2009) Structural control over equilibrium silicon and oxygen isotopic fractionation: a first–principles densityfunctional theory study. Chem Geol 258:28–37
Dupuis R, Benoit M, Nardin E, Méheut M (2015) Fractionation of silicon isotopes in liquids: the importance of configurational disorder. Chem Geol 396:239–254
Dai S, Chou CL, Yue M, Luo K, Ren D (2005) Mineralogy and geochemistry of a Late Permian coal in the Dafang Coalfield, Guizhou, China: influence from siliceous and ironrich calcic hydrothermal fluids. Int J Coal Geol 61:241–258
Douthitt C (1982) The geochemistry of the stable isotopes of silicon. Geochim Cosmochim Acta 46:1449–1458
Ding T, Jiang S, Wan D, Li Y, Li J, Song H, Liu Z, Lao X (1996) Silicon Isotope Geochemistry. Geological Publishing House, Beijing
Opfergelt S, Burton KW, Pogge von Strandmann PAE, Gislason SR, Halliday AN (2013) Riverine silicon isotope variations in glaciated basaltic terrains: implications for the Si delivery to the ocean over glacial–interglacial intervals. Earth Planet Sci Lett 369–370:211–219
Geilert S, Vroon PZ, Keller NS, Gudbrandsson S, Stefansson A, van Bergen MJ (2015) Silicon isotope fractionation during silica precipitation from hot–spring waters: evidence from the Geysir geothermal field. Geochim Cosmochim Acta 164:403–427
De La Rocha CL, Brzezinski MA, DeNiro MJ (2000) A first look at the distribution of the stable isotopes of silicon in natural waters. Geochim Cosmochim Acta 64:2467–2477
O’Reilly C, Parnell J (1999) Fluid flow and thermal histories for CambrianOrdovician platform deposits, New York: evidence from fluid inclusion studies. Geol Soc Am Bull 111:1884–1896
Kita I, Taguchi S, Matsubaya O (1985) Oxygen isotope fractionation between amorphous silica and water at 34–93 °C. Nature 314:83–84
Sharp ZD, Gibbons JA, Maltsev O, Atudorei V, Pack A, Sengupta S, Shock EL, Knauth LP (2016) A calibration of the triple oxygen isotope fractionation in the SiO_{2}–H_{2}O system and applications to natural samples. Geochim Cosmochim Acta. https://doi.org/10.1016/j.gca.2016.04.047
Clayton RN (1972) Oxygen isotope exchange between quartz and water. J Geophys Res 77:3057
Matthews A, Beckinsale RD (1979) Oxygen isotope equilibration systematics between quartz and water. Am Mineral 64:232–240
Matsuhisa Y, Goldsmith JR, Clayton RN (1978) Mechanisms of hydrothermal crystallization of quartz at 250 °C and 15 kbar. Geochim Cosmochim Acta 42:173–182
Matsuhisa Y, Goldsmith JR, Clayton RN (1979) Oxygen isotopic fractionation in the system quartz–albite–anorthite–water. Geochim Cosmochim Acta 43:1131–1140
Matthews A, Goldsmith JR, Clayton RN (1983) On the mechanisms and kinetics of oxygen isotope exchange in quartz and feldspars at elevated temperatures and pressures. Geol Soc Am Bull 94:396–412
Hu G, Clayton RN (2003) Oxygen isotope salt effects at high pressure and high temperature and the calibration of oxygen isotope geothermometers. Geochim Cosmochim Acta 67:3227–3246
Pollington AD, Kozdon R, Anovitz LM, Bastian Georg R, Spicuzza MJ, Valley JW (2016) Experimental calibration of silicon and oxygen isotope fractionations between quartz and water at 250 °C by in situ microanalysis of experimental products and application to natural samples. Chem Geol 421:127–142
Komárek M, Ratié G, Vaňková Z, Šípková A, Chrastný V (2021) Metal isotope complexation with environmentally relevant surfaces: opening the isotope fractionation black box. Crit Rev Environ Sci Technol. https://doi.org/10.1080/10643389.2021.1955601
Gou WX, Li W, Ji JF, Li WQ (2018) Zinc isotope fractionation during sorption onto Al Oxide: atomic level understanding from EXAFS. Environ Sci Technol 52:9087–9096
Pokrovsky OS, Viers J, Freydier R (2005) Zinc stable isotope fractionation during its adsorption on oxides and hydroxides. J Colloid Interface Sci 291:192–200
Horner TJ, Rickaby REM, Henderson GM (2011) Isotopic fractionation of cadmium into calcite. Earth Planet Sci Lett 312:243–253
Xie XJ, Yan L, Li JX, Guan LR, Chi ZY (2020) Cadmium isotope fractionation during Cd–calcite coprecipitation: insight from batch experiment. Sci Total Environ 760:143330
Acknowledgements
The authors are grateful for the financial support from the Hebei Natural Science Foundation, Hebei Education Department Key Program and the State Natural Sciences Foundation.
Funding
This work was financially supported by the Hebei Natural Sciences Foundation (Nos. D2020402004 and D2021402020), Hebei Education Department Key Program (No. ZD2018086), and State Natural Sciences Foundation (No. 41603011).
Author information
Authors and Affiliations
Contributions
YFW carries out calculations and analyses, and prepares the first draft; PDW helps to collect data; XYJ, LCX, JL, TDZ and HNZ assist in analysis and discussion; YFW and HTH design the simulation frame, discuss and edit the manuscript. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Competing interests
There is no conflict of interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1: Figure S1
. Hydrogen bonds connected water molecules and scheme of point charge arrangements (PCA) for 1×, 3×, 4×, and 5×.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Wang, YF., Ji, XY., Xing, LC. et al. Improved volume variable cluster model method for crystallattice optimization: effect on isotope fractionation factor. Geochem Trans 23, 1 (2022). https://doi.org/10.1186/s12932022000786
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12932022000786
Keywords
 Molecular cluster
 VVCM
 Geometric optimization
 Relative volume change
 Isotopic equilibrium fractionation factor