Skip to main content

Improved volume variable cluster model method for crystal-lattice optimization: effect on isotope fractionation factor


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 (cadmium-containing calcite) and their aqueous solutions under superficial conditions. We summarized a detailed procedure and used it to re-evaluate published theoretical results for cadmium-containing 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–Al-layered double hydroxide was 2.10 Å, with a relative volume change of 1.96%. Cadmium-containing 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–Al-layered double hydroxide, cadmium-containing calcite, and their corresponding aqueous solutions were \(\Delta ^{30/28} {\text{Si}}_{{{\text{Qtz-H}}_{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 LDH-Zn}} ( {\text{H}}_{2} {\text{O}} )_{\text{n}}^{2+}} = 1.12{\permil}\) and \(\Delta^{114/110} {\text{Cd}}_{ {\text{(Cd--Cal)-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.


In nature, the isotopic composition of an element varies with rock type, source, and age, especially those determined by high-resolution 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 δ7Li and δ114Cd 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 laboratory-derived 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, first-principles 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 second-order 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 strength-conserving quasi-atoms” along the clipped bonds, which make the core charge of each quasi-atom equal to the Pauling bond valence, that is, the charge of the removed cation divided by its coordination number. The distance between the quasi-atom 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 quasi-atoms 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 well-established. 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 four-parameter Vinet equation of state [35, 36] can be expressed as follows:

$${\text{E}}({\text{V}}) = {\text{a}} - \frac{{4{\text{B}}_{0} {\text{V}}_{0} }}{{\left( {{\text{B}}_{0}^{\prime } - 1}\right)^{2} }}\left\{ {1 - \frac{3}{2}\left( {{\text{B}}_{0}^{\prime } - 1} \right)\left[ {1 - \left( {1 + \frac{{\Delta {\text{V}}}}{{{\text{V}}_{0} }}} \right)^{{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}}} } \right]} \right\} \times \exp \left\{ {\frac{3}{2}\left( {{\text{B}}_{0}^{\prime } - 1} \right)\left[ {1 - \left( {1 + \frac{{\Delta {\text{V}}}}{{{\text{V}}_{0} }}} \right)^{{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}}} } \right]} \right\}$$

where the fitting parameter a = E0 + 4B0V0/(B0-1)2, V0, B0, and B0 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 cadmium-containing 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.

Fig. 1
figure 1

Flow chart of the improved VVCM method

Constructing molecular cluster and adding background point charges

We reconstructed the three-dimensional crystal lattices of the minerals from X-ray single-crystal 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 electron-neutral environment of the crystal. Existence of free interlayer anions makes LDH different from others. The general formula of LDH is [M2+1−x M3+x (OH)2]x+(An−)x/n·mH2O, where M2+ and M3+ represent divalent and trivalent metal cations, respectively; An− is an interlayer anion including CO32−, Cl, SO42−, NO3; x represents the molar ratio M3+/(M2+  + M3+); and m is the number of water molecules [44, 45]. In nature, the M2+/M3+ molar ratio is usually two [46]. Therefore, we focused on a typical Zn-containing 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 non-shared electron pair of the oxygen atom. Meanwhile, oxygen atom can contribute its two non-shared 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/V0) 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].

Fig. 2
figure 2

Schematic diagram of PCA

Table 1 Comparison of theoretical results for three optimal PCAs

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 Si4+, 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 Zn2+, Al3+, and H atoms. Specifically, the total amount of charge was the summation of the charge contributed by one hexa-coordinated Zn2+ (+ 2/6), one hexa-coordinated Al3+ (+ 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 hexa-coordinated Al3+ (+ 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 hexa-coordinated Ca2+ (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.

Table 2 PCA on the anion or anion group in the outermost layer of a molecular cluster

The configuration of Cd-containing 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 nine-coordinated Ca2+ and one seven-coordinated Ca2+ (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 seven-coordinated Ca2+ and one nine-coordinated Ca2+, 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 nine-coordinated Ca2+ (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 seven-coordinated Ca2+ (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 seven-coordinated Ca2+ and a nine-coordinated Ca2+ (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 seven-coordinated Ca2+ is removed; and the point charge on μ3-O obtains a charge of + 2/9 when the nine-coordinated Ca2+ 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.


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) is approximately 100–120 pm [22]. Thus, we adjusted d 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.

Fig. 3
figure 3

Electronic energy of the Si–Qtz system varies with d

Fig. 4
figure 4

Electronic energy of the Zn–Al LDH system varies with d and d

Fig. 5
figure 5

Electronic energy of the Cd–Cal system varies with d and d

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) and/or the distance between the 2 × point charge and μ2-O (d). Subsequently, we successively adjusted d until we determined the optimal value. For the Zn–Al LDH and Cd–Cal molecular clusters, we determined the optimal d or d after fixing d5×. However, confirming the most suitable distance d (d) in the Cd–HAp molecular cluster requires d (d) and d 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 all-electronic 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/V0) 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/Bohr3, 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/Bohr3. Finally, we counted the number of points (n) with electron densities larger than x e/Bohr3, 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/Bohr3 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/Bohr3 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/Bohr3 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 × 2i 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/Bohr3 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/Bohr3. We used the following formula for these calculations:

$${{\Delta {\text{V}}} \mathord{\left/ {\vphantom {{\Delta {\text{V}}} {{\text{V}}_{0} }}} \right. \kern-\nulldelimiterspace} {{\text{V}}_{0} }} = \frac{{{\text{V}}_{{{\text{opt}}}} - {\text{V}}_{0} }}{{{\text{V}}_{0} }} \times 100\%$$

where Vopt denotes the optimized volume and V0 is the unoptimized volume. A positive ΔV/V0 indicates expansion of the molecular cluster, whereas a negative ΔV/V0 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:

$${\text{AX}}^{\prime} + {\text{X}} = {\text{AX}} + {\text{X}}^{\prime}$$

where AX and AXʹ represent compounds with heavy and light isotopes, respectively; X and Xʹ denote the ideal single-atom 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:

$${\text{K}} = \frac{{{\text{Q}}_{trans}^{{{\text{AX}}}} {\text{Q}}_{rot}^{{{\text{AX}}}} Q_{vib}^{{{\text{AX}}}} Q_{elec}^{{{\text{AX}}}} }}{{{\text{Q}}_{{{\text{trans}}}}^{{{\text{AX}}}^{\prime}} {\text{Q}}_{{{\text{rot}}}}^{{AX}}}^{\prime} Q_{{{\text{vib}}}}^{{{\text{AX}}}^{\prime}} Q_{{{\text{elec}}}}^{{AX}^{\prime}}}/\frac{{{\text{Q}}_{{{\text{trans}}}}^{{\text{X}}} {\text{Q}}_{{{\text{elec}}}}^{{\text{X}}}}}{{Q_{{{\text{trans}}}}^{{X}^{\prime}}}Q_{{{\text{elec}}}}^{{X}^{\prime}}}$$

where Qtrans is the transitional partition function, Qrot is the rotational partition function, Qvib is the vibrational partition function, and Qelec is the electronic partition function [61]. For light elements, the difference in the ground-state electronic energy of the compound before and after isotope substitution is indistinguishable, so Qelec can be ignored. The translational, rotational, and vibrational partition functions can be respectively written as follows:

$${\text{Q}}_{{{\text{trans}}}} = {\text{V}}\left( {\frac{{2\uppi {\text{MkT}}}}{{{\text{h}}^{2} }}} \right)^{3/2}$$
$${\text{Q}}_{{{\text{rot}}}} = \frac{{\uppi ^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}}} \left( {8\uppi ^{2} {\text{kT}}} \right)^{{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-\nulldelimiterspace} 2}}} \left( {{\text{I}}_{{\text{A}}} {\text{I}}_{{\text{B}}} {\text{I}}_{{\text{C}}} } \right)^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}}} }}{{\sigma {\text{h}}^{3} }}$$
$${\text{Q}}_{{{\text{vib}}}} = \mathop \prod \limits_{{\text{i}}}\frac{{{\text{e}}^{{ - {\text{u}}_{{{{_{{\text{i}}}} \mathord{\left/ {\vphantom {{\text{i}} 2}} \right. \kern-\nulldelimiterspace} 2}}} }} }}{{1-{\text{e}}^{{ - {\text{u}}_{_{\text{i}}} }} }}$$

where V is the volume, M is the mass, IA is the moment of inertia around the axis A of rotation, σ is the symmetry number of the molecule, and ui = 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 Teller-Redlich [62] approximation is done as follows:

$$\left( {\frac{{\text{M}}}{{\text{M}}^{\prime}}} \right)_{{{\text{AX}}}}^{3/2} \left( {\frac{{{\text{I}}_{{\text{A}}} {\text{I}}_{{\text{B}}} {\text{I}}_{{\text{C}}} }}{{{{\text{I}}^{\prime}}_{{\text{A}}} {{\text{I}}^{\prime}}_{{\text{B}}} {{\text{I}}^{\prime}}_{{\text{C}}}}}} \right)_{{{\text{AX}}}}^{1/2} \left( {\frac{{{\text{m}}^{\prime}}}{{\text{m}}}} \right)_{{\text{X}}}^{{3{\text{n}}/2}} = \mathop \prod \limits_{{\text{i}}} \frac{{{\text{u}}_{_{{\text{i}}}} }}{{{{\text{u}}^{\prime}}_{_{{\text{i}}}} }}$$

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:

$${\text{RPFR}}({\text{AX}}/{\text{AX}}^{\prime}) = \frac{\sigma }{{\sigma}^{\prime}}{\text{K}} = \mathop \prod \limits_{}^{{3{\text{n}} - 6}} \frac{{\text{u}}}{{{\text{u}}}^{\prime}}\left( {\frac{{{\text{e}}^{{ - \frac{{\text{u}}}{2}}} }}{{{\text{e}}^{{ - \frac{{{\text{u}}}^{\prime}}{2}}} }}} \right)\left( {\frac{{1 - {\text{e}}^{{ - {\text{u}}^{\prime}}} }}{{1 - {\text{e}}^{{ - {\text{u}}}}}}} \right)$$

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: α = βAB, where A and B denote different compounds or minerals. The isotopic equilibrium fractionation between the two phases can then be obtained as follows:

$$\Delta_{{_{\text{A}} - _{\text{B}}}} \approx 10^{3} \ln {\upalpha }_{{{\text{A}} - {\text{B}}}} = 10^{3} ({\text{ln}}\beta_{{_{\text{A}}}} - {\text{ln}}\beta _{_{\text{B}}} )$$


Influence of d, d, and d 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, d, and d, and the R2 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.

Fig. 6
figure 6

Electronic energy of the Cd–HAp system varies with d, d and d

Effect of d, d, and d on relative volume change

During the geometric optimization process, the volume of the Si–Qtz molecular cluster decreased, with ΔV/V0 ranging from − 2.05 to  − 0.11% (Fig. 7). We observed an extremely strong linear correlation between ΔV/V0 and d (R2 = 0.9895). ΔV/V0 gradually increased with an increase in d. We recorded the ΔV/V0 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/V0 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/V0 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].

Fig. 7
figure 7

ΔV/V0 variations in the Si–Qtz molecular cluster

Fig. 8
figure 8

ΔV/V0 variations in the Zn–Al LDH molecular cluster. The red points were the relative volume changes for molecular cluster without interlayer anions, while the black point was that for molecular cluster with interlayer anions

Fig. 9
figure 9

ΔV/V0 variations in the Cd–Cal molecular cluster

Fig. 10
figure 10

ΔV/V0 variations in the Cd–HAp molecular cluster


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 X-ray diffraction. Our results are consistent with those of previous theoretical calculations and experimental studies.

Table 3 Isotopic equilibrium fractionation between molecular clusters and their corresponding aqueous solutions

The central Zn2+ 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 X-ray diffraction and EXAFS experiments, obtained an average Zn–O bond length of 2.08 Å (coordination number = 6.5) in the Zn–Al layer-like 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/V0 values of Ba-containing species by using the periodic boundary method. Except for Ba(OH)2(H2O)8, the relative volume change of which was -0.67%, the ΔV/V0 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/V0 of various Zn-containing 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/V0 values calculated using the molecular cluster method are much smaller (Table 4), approximately one-tenth of those obtained by the periodic boundary method.

Table 4 Relative volume change of solid molecular clusters

We found that the ΔV/V0 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/V0 of quartz using GGA (4.61%, PW91) and LDA (-2.75%). The ΔV/V0 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 Cd-containing calcite also increased (0.45%). However, the relative volume change of the Mn-containing calcite (3.54%) [70] and tetrahedral AO42− 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/V0 of 3.43% using GGA–PBE as the exchange–correlation functional. Bhat et al. [73] obtained ΔV/V0 (4.37%) using the exchange–correlation functional GGA–PW91. In contrast, the ΔV/V0 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 H4SiO4 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 H4SiO4 solution at 298 K. Dupuis et al. [75] calculated the fractionation factors of Si isotopes between quartz and aqueous H4SiO4 at 273–323 K using first-principles methods, and found that at 300 K, the fractionation between quartz and H4SiO4 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 H4SiO4 solutions. We considered the experimentally determined ∆30/28Siquartz-aqueous 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.

Fig. 11
figure 11

Theoretical Si isotope equilibrium fractionation factors between quartz and aqueous H4SiO4 solution [5, 7, 22, 75,76,77,78,79,80,81,82]

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 low-temperature 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.

Fig. 12
figure 12

Theoretical O isotope equilibrium fractionation factors between quartz and liquid water [6, 68, 83,84,85,86,87,88,89,90,91]

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 Zn2+ 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 γ-Al2O3 with a fractionation factor of ∆66/64Znsolid-aq = 0.02 ± 0.07‰. Pokrovsky et al. [94] also reported lower fractionation on corundum (∆66/64Znsolid-aq = 0.19 ± 0.05‰) and gibbsite (∆66/64Znsolid-aq = 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 γ-Al2O3. 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 Cd-containing aqueous solution was -0.26‰. Previous studies on precipitated calcite from artificial seawater solutions reported the entry of Cd2+ into the calcite from seawater, with the Cd equilibrium isotope fractionation, ∆114/110Cdcalcite-seawater = -0.45 ± 0.12‰ [95]. Furthermore, Xie et al. [96] revealed apparent Cd isotope fractionation in the coprecipitation of Cd–calcite, with ∆114/110Cd calcite-aqueous = -0.38 ± 0.18‰ to -0.20 ± 0.12‰. Our theoretical calculation results are consistent with previous experimental results.


We used the first-principles 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 crystal-lattice optimization.

Availability of data and materials

The data sets that support the conclusions of this study are included in the study.



Volume variable cluster model


Density Functional Theory


Quartz centered on silicon atoms


Quartz centered on oxygen atoms


Layered double hydroxide


Cd–containing calcite


Cd–containing hydroxyapatite


Point charge arrangement


  1. Zhu X-C, 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

    Article  Google Scholar 

  2. Erkaev NV, Scherf M, Thaller SE, Lammer H, Mezentsev AV, Ivanov VA, Mandt KE (2020) Escape and evolution of Titan’s N2 atmosphere constrained by 14N/15N isotope ratios. Mon Not R Astron Soc 500:2020–2035

    Article  Google Scholar 

  3. 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

    Article  Google Scholar 

  4. 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

    Article  Google Scholar 

  5. 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

    Article  Google Scholar 

  6. Labeyrie LJ (1974) New approach to surface seawater paleotemperatures using18O/16O ratios in silica of diatom frustules. Nature 248:40–42

    Article  Google Scholar 

  7. 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 low-temperature hydrothermal fluid. Appl Geochem 19:1315–1330

    Article  Google Scholar 

  8. Clayton RN, Epstein S (1961) The use of oxygen isotopes in high–temperature geological thermometry. J Geol 69:447–452

    Article  Google Scholar 

  9. Klein RT, Lohmann KC, Thayer CW (1996) Bivalve skeletons record sea-surface temperature and δ18O via Mg/Ca and 18O/16O ratios. Geology 24:415

    Article  Google Scholar 

  10. 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

    Article  Google Scholar 

  11. 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 U-Pb dating, geochemistry and Nd isotope. J Asian Earth SCI 37:140–153

    Article  Google Scholar 

  12. 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

    Article  Google Scholar 

  13. 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 C-O–S–Pb isotope constraints. Ore Geol Rev 79:88–104

    Article  Google Scholar 

  14. 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 Permian-Triassic crisis from the Meishan section, South China. Chem Geol 481:110–118

    Article  Google Scholar 

  15. 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 Permian-Triassic boundary. Proc Natl Acad of Sci 115:3782–3787

    Article  Google Scholar 

  16. Eiler JM (2011) Paleoclimate reconstruction using carbonate clumped isotope thermometry. Quat Sci Rev 30:3575–3588

    Article  Google Scholar 

  17. Rustad JR, Dixon DA (2009) Prediction of iron-isotope fractionation between hematite (α–Fe2O3) and ferric and ferrous iron in aqueous solution from density functional theory. J Phys Chem A 113:12249–12255

    Article  Google Scholar 

  18. Rustad JR, Yin QZ (2009) Iron isotope fractionation in the Earth’s lower mantle. Nat Geosci 2:514–518

    Article  Google Scholar 

  19. Rustad JR, Casey WH, Yin QY, Bylaska E, Felmy AR, Bogatko SA, Jackson VE, Dixon DA (2010) Isotopic fractionation of Mg2+(aq), Ca2+(aq), and Fe2+(aq) with carbonate minerals. Geochim Cosmochim Acta 74:6301–6323

    Article  Google Scholar 

  20. Rustad JR, Nelmes SL, Jackson VE, Dixon DA (2008) Quantum-chemical calculations of carbon–isotope fractionation in CO2(g), aqueous carbonate species, and carbonate minerals. J Phys Chem A 112:542–555

    Article  Google Scholar 

  21. 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)

  22. He H-T, Liu Y (2015) Si isotope fractionations during the precipitation of quartz and the adsorption of H4SiO4(aq) on Fe(III)–oxyhydroxide surfaces. Chin J Geochem 34:459–468

    Article  Google Scholar 

  23. He H-T, 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

    Article  Google Scholar 

  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

    Article  Google Scholar 

  25. He H-T, Xing L–C, Qin S–J, Ji X–Y, Sun P–F, (2020) Equilibrium Cd isotopic fractionation between Cd(OH)2(s), apatite, adsorbed Cd2+, and Cd2+(aq): Potential application of δ114Cd in evaluating the effectiveness of Cd–contamination remediation. Geochem J 54:289–297

    Article  Google Scholar 

  26. Zhang JX, Liu Y (2018) Zinc isotope fractionation under vaporization processes and in aqueous solutions. Acta Geochim 37:27–39

    Google Scholar 

  27. 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

    Article  Google Scholar 

  28. Zhang J (2021) Equilibrium sulfur isotope fractionations of several important sulfides. Geochem J 55:135–147

    Article  Google Scholar 

  29. Fujii T, Moynier F, Albarède B-T, 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

    Article  Google Scholar 

  30. 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

    Article  Google Scholar 

  31. 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

    Article  Google Scholar 

  32. Ducher M, Blanchard M, Balan E (2018) Equilibrium isotopic fractionation between aqueous Zn and minerals from first–principles calculations. Chem Geol 483:342–350

    Article  Google Scholar 

  33. Blanco MA, Francisco E, Luana V (2004) GIBBS: isothermal-isobaric thermodynamics of solids from energy curves using a quasi-harmonic Debye model. Comput Phys Commun 158:57–72

    Article  Google Scholar 

  34. Shang SL, Wang Y, Kim D, Liu ZK (2010) First–principles thermodynamics from phonon and Debye model: application to Ni and Ni3Al. Comput Mater Sci 47:1040–1048

    Article  Google Scholar 

  35. Vinet P, Ferrante J, Smith JR, Rose JH (1986) A universal equation of state for solids. J Phys C 19:L467–L473

    Article  Google Scholar 

  36. 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

    Article  Google Scholar 

  37. Liu WT, Shen YR (2008) Surface Vibrational Modes of α–Quartz (0001) probed by sum–frequency spectroscopy. Phys Rev Lett 101:016101

    Article  Google Scholar 

  38. Tang C, Zhu J, Zhou Q, Wei J, Zhu R, He H (2014) Surface heterogeneity of SiO2 polymorphs: an XPS investigation of α–Quartz and α–Cristobalite. J Phys Chem C 118:26249–26257

    Article  Google Scholar 

  39. Hazen RM, Finger LW, Hemley RJ, Mao HK (1989) High–pressure crystal chemistry and amorphization of α–quartz. Solid State Commun 72:507–511

    Article  Google Scholar 

  40. Paquette J, Reeder RJ (1990) Single–crystal X–ray structure refinements of two biogenic magnesian calcite crystals. Am Mineral 75:1151–1158

    Google Scholar 

  41. Wilson RM, Elliot JC, Dowker SEP (1999) Rietveld refinement of the crystallographic structure of human dental enamel apatites. Am Mineral 84:1406–1414

    Article  Google Scholar 

  42. Fleet ME, Liu XY, Pan YM (2000) Rare–earth elements in chlorapatite [Ca10(PO4)6Cl2]: uptake, site preference, and degradation of monoclinic structure. Am Mineral 85:1437–1446

    Article  Google Scholar 

  43. Zhu SD, Khan MA, Wang FY, Bano Z, Xia MZ (2020) Rapid removal of toxic metals Cu2+ and Pb2+ by amino trimethylene phosphonic acid intercalated layered double hydroxide: a combined experimental and DFT study. Chem Eng J 392:123711

    Article  Google Scholar 

  44. 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

    Article  Google Scholar 

  45. 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

    Article  Google Scholar 

  46. Costa DG, Rocha AB, Diniz R, Souza WF, Chiaro SSX, Leitão AA (2010) Structural model proposition and thermodynamic and vibrational analysis of hydrotalcite-like compounds by DFT calculations. J Phys Chem C 114:14133–14140

    Article  Google Scholar 

  47. 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

    Google Scholar 

  48. Kihara K (1990) An X–ray study of the temperature dependence of the quartz structure. Eur J Mineral 2:63–77

    Article  Google Scholar 

  49. 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

    Article  Google Scholar 

  50. 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

    Article  Google Scholar 

  51. 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

    Article  Google Scholar 

  52. 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

    Article  Google Scholar 

  53. 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

    Article  Google Scholar 

  54. 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.

  55. 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

    Article  Google Scholar 

  56. 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

    Article  Google Scholar 

  57. Bader RFW, Carroll MT, Cheeseman JR, Chang C (1987) Properties of atoms in molecules: atomic volumes. J Am Chem Soc 109:7968–7979

    Article  Google Scholar 

  58. Lu T, Chen F (2012) Multiwfn: a multifunctional wavefunction analyzer. J Comput Chem 33:580–592

    Article  Google Scholar 

  59. Bigeleisen J, Mayer MG (1947) Calculation of equilibrium constants for isotopic exchange reactions. J Chem Phys 15:261–267

    Article  Google Scholar 

  60. Urey HC (1947) The thermodynamic properties of isotopic substances. J Chem Soc 85:562–581

    Article  Google Scholar 

  61. 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

    Article  Google Scholar 

  62. Redlich O (1935) Eine allgemeine beziehung zwischen den schwingungsfrequenzen isotoper molekeln. Z Phys Chem B 28:371–382

    Article  Google Scholar 

  63. 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

    Article  Google Scholar 

  64. 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

    Article  Google Scholar 

  65. 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

    Article  Google Scholar 

  66. Hata M, Okada K, Iwai S, Akao M, Aoki H (1978) Cadmium hydroxyapatite. Acta Crystallogr Sect B: Struct Sci 34:3062–3064

    Article  Google Scholar 

  67. Ducher M, Blanchard M, Balan E (2016) Equilibrium zinc isotope fractionation in Zn–bearing minerals from first–principles calculations. Chem Geol 443:87–96

    Article  Google Scholar 

  68. Méheut M, Lazzeri M, Balan E, Mauri F (2007) Equilibrium isotopic fractionation in the kaolinite, quartz, water system: prediction from first–principles density-functional theory. Geochim Cosmochim Acta 71:3170–3181

    Article  Google Scholar 

  69. Hamann DR (1996) Generalized gradient theory for silica phase transitions. Phys Rev Lett 76:660–663

    Article  Google Scholar 

  70. Son SB, Newton AG, Jo KN, Lee JY, Kwon KD (2019) Manganese speciation in Mn–rich CaCO3: a density functional theory study. Geochim Cosmochim Acta 248:231–241

    Article  Google Scholar 

  71. 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

    Article  Google Scholar 

  72. 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

    Article  Google Scholar 

  73. Bhat SS, Waghmare UV, Ramamurty U (2014) First-principles study of structure, vibrational, and elastic properties of stoichiometric and calcium-deficient hydroxyapatite. Cryst Growth Des 14:3131–3141

    Article  Google Scholar 

  74. Méheut M, Lazzeri M, Balan E, Mauri F (2009) Structural control over equilibrium silicon and oxygen isotopic fractionation: a first–principles density-functional theory study. Chem Geol 258:28–37

    Article  Google Scholar 

  75. 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

    Article  Google Scholar 

  76. Dai S, Chou C-L, 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 iron-rich calcic hydrothermal fluids. Int J Coal Geol 61:241–258

    Article  Google Scholar 

  77. Douthitt C (1982) The geochemistry of the stable isotopes of silicon. Geochim Cosmochim Acta 46:1449–1458

    Article  Google Scholar 

  78. Ding T, Jiang S, Wan D, Li Y, Li J, Song H, Liu Z, Lao X (1996) Silicon Isotope Geochemistry. Geological Publishing House, Beijing

    Google Scholar 

  79. 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

    Article  Google Scholar 

  80. 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

    Article  Google Scholar 

  81. 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

    Article  Google Scholar 

  82. O’Reilly C, Parnell J (1999) Fluid flow and thermal histories for Cambrian-Ordovician platform deposits, New York: evidence from fluid inclusion studies. Geol Soc Am Bull 111:1884–1896

    Article  Google Scholar 

  83. Kita I, Taguchi S, Matsubaya O (1985) Oxygen isotope fractionation between amorphous silica and water at 34–93 °C. Nature 314:83–84

    Article  Google Scholar 

  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 SiO2–H2O system and applications to natural samples. Geochim Cosmochim Acta.

    Article  Google Scholar 

  85. Clayton RN (1972) Oxygen isotope exchange between quartz and water. J Geophys Res 77:3057

    Article  Google Scholar 

  86. Matthews A, Beckinsale RD (1979) Oxygen isotope equilibration systematics between quartz and water. Am Mineral 64:232–240

    Google Scholar 

  87. 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

    Article  Google Scholar 

  88. Matsuhisa Y, Goldsmith JR, Clayton RN (1979) Oxygen isotopic fractionation in the system quartz–albite–anorthite–water. Geochim Cosmochim Acta 43:1131–1140

    Article  Google Scholar 

  89. 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

    Article  Google Scholar 

  90. 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

    Article  Google Scholar 

  91. 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

    Article  Google Scholar 

  92. 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.

    Article  Google Scholar 

  93. 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

    Article  Google Scholar 

  94. 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

    Article  Google Scholar 

  95. Horner TJ, Rickaby REM, Henderson GM (2011) Isotopic fractionation of cadmium into calcite. Earth Planet Sci Lett 312:243–253

    Article  Google Scholar 

  96. 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

    Article  Google Scholar 

Download references


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.


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



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

Correspondence to Hong-Tao He.

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 The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Wang, YF., Ji, XY., Xing, LC. et al. Improved volume variable cluster model method for crystal-lattice optimization: effect on isotope fractionation factor. Geochem Trans 23, 1 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Molecular cluster
  • VVCM
  • Geometric optimization
  • Relative volume change
  • Isotopic equilibrium fractionation factor