Interactions between calcium ion and functional groups of organic scale inhibitors in aqueous solutions: an ab initio study

Abstract

The anti-scaling performance of organic inhibitors is closely related to the Ca-trapping capacity of functional groups in polymer frameworks. In this study, we performed static density functional theory calculations and ab initio metadynamics simulations to systematically investigate the association behaviors of Ca²⁺ with carbonate species (CO₃²⁻ and HCO₃⁻) and representative functional groups of organic inhibitors (carboxylate, phosphonate, and sulfonate) in aqueous solutions. By calculating the binding free energies and activation barriers for separating contact ion pairs, we propose a new strategy for predicting the Ca-trapping capacity of organic additives. Simulation results revealed that carboxylate and phosphonate groups possess potential calcium sequestration abilities, while the formation of the calcium–sulfonate pair is almost impossible in aqueous solutions. For the carboxylate group, which possesses a planar structure akin to that of the carbonate species, its presence has less impact on the solvation structure of Ca²⁺. The monodentate and bidentate configurations have similar stabilities for the contact ion pair. However, for the phosphonate and sulfonate groups that possess a trigonal pyramid structure, their association with Ca²⁺ disrupts the symmetry of the coordination structure to a greater extent. The monodentate state is clearly favored over the more tightly coordinated bidentate state. Moreover, since CO₃²⁻ exhibits the strongest affinity for Ca²⁺, a rising CO₃²⁻ concentration at a high pH diminishes the anti-scaling efficiency of all inhibitors. Overall, we offer new insights into the inhibiting effects of organic additives on the initial stage of CaCO₃ formation by analyzing calcium-organic compound interactions. These findings provide theoretical support for the development of more effective scale inhibitors in the industry.

---

1. Introduction

Calcium carbonate (CaCO₃) scale deposition has been recognized as one of the major issues in industrial production.¹ The development of scale layers in circulating cooling water systems blocks pipes and boilers and lowers the efficiency of heat exchangers, thus resulting in severe machine damage.² In the oil and gas industries, water-flooding operations are commonly required to enhance secondary recovery from mature fields.³ However, excessive precipitation of scale impedes fluid flow in the wellbore, reservoir, tubing, casing, and even the oil-producing formation matrix, reducing well productivity and ultimately leading to considerable treatment costs and deterioration of the oil extraction process.⁴ Additionally, membrane distillation has emerged as one of the most promising technologies in desalination for water separation and purification.⁵ Nevertheless, undesired scale deposits progressively clog the pores and increase the wettability of the membrane, resulting in a decline in both the permeate flux and separation efficiency.^6,7 Therefore, in recent years, increasing attention has been directed toward exploring the nucleation and growth mechanisms of CaCO₃ and how to regulate such crystallization events in a particular context.^8–10

A typical method for hindering CaCO₃ precipitation is the use of environmentally friendly chemicals, which act as anti-scaling reagents, that is, inhibitors.¹¹ The commonly used organic inhibitors include, but are not limited to, carboxylic acids (e.g., polyacrylic acid (PAA),¹² polyepoxysuccinic acid (PESA),¹³ and polyaspartic acid (PASP)¹⁴), phosphonic acids (e.g., amino trimethylene phosphonic acid (ATMP),¹⁵ 2-phosphonobutane-1,2,4-tricarboxylic acid (PBTCA),¹⁶ and diethylenetriamine pentamethylene phosphonic acid (DTPMP)¹⁷), and sulfonates (e.g., polystyrene sulfonate (PSS)¹⁸).

As discussed in detail in our previous review article,¹⁹ inhibitors can restrain scale deposition by influencing one or multiple stages of CaCO₃ crystallization, ranging from the aggregation of prenucleation clusters (PNCs), formation of amorphous CaCO₃ (ACC) particles, and crystallization of ACC to the growth of crystals. Notably, the anti-scaling performance of an inhibitor is largely determined by the type, concentration, and arrangement of the functional groups within the polymer structure.²⁰ In titration experiments conducted by Verch et al.,²¹ carboxylate groups were found to possess a stronger affinity toward Ca²⁺ than sulfonate groups in aqueous solutions, so the carboxylate-containing additive PASP can reduce the concentration of free Ca²⁺ ions involved in CaCO₃ crystallization to a greater extent and therefore manifests a stronger inhibiting effect on scale deposition than the sulfonate-containing additive PSS. Li et al.⁸ employed scanning electron microscopy (SEM) to observe the crystal structure of scale deposits in both the presence and absence of inhibitors. They discovered that the carboxylate-containing additive PAA exhibits a rather insignificant effect on the surface structure of CaCO₃ crystals, while the phosphonate-containing additives PBTCA and DTPMP fully open the crystal layers and effectively retard the growth of CaCO₃ in a three-dimensional direction. Cui et al.²² synthesized a terpolymer scale inhibitor incorporating carboxylate, sulfonate, and amide groups and proposed that the synergistic effects of multiple functional groups in the polymer structure are beneficial in improving inhibition efficiency.

Although exceptional progress has been made in experimental techniques for evaluating the anti-scaling performance of inhibitors, the cardinal fact that one has to acknowledge is that the discovery of new inhibitors still relies on the labor-intensive “trial-and-error” paradigm, which consumes too much time and resources. In the past two decades, with the rapid development of computing facilities, there have been major advances in computer simulation, which plays a momentous role in providing atomic-level insights to uncover the physical mechanisms underlying the convoluted phenomena and thus exhibits remarkable potential to offer suggestions for the continuous development of more effective scale inhibitors in the industry.

To the best of our knowledge, the earliest attempt that employed atomistic simulations to investigate the inhibiting effects of organic additives on CaCO₃ formation was the classical molecular dynamics (MD) studies conducted by de Leeuw and Cooper.^23–25 They used a combination of three empirical potential models to describe the interatomic interactions in the system: the model by Pavese et al.²⁶ for CaCO₃ crystals, the consistent valence force field (CVFF)²⁷ for organic molecules, and the water potential model by de Leeuw and Parker.²⁸ Based on an analysis of the adsorption properties of a selection of organic molecules on calcite surfaces, they suggested that carboxylic acids can block the active growth sites of CaCO₃ crystals, reduce the direct contact between surfaces and solutions, and thus retard crystal growth. Duchstein et al.²⁹ utilized the assisted model building and energy refinement (AMBER) force field³⁰ for the organic part and performed MD simulations to examine the dynamic interaction modes of carboxylate- and phosphonate-containing additives with PNCs. Their findings revealed the significance of scale inhibitors as “surfactants” in the initial stage of CaCO₃ formation, reducing interactions among PNCs, restraining their aggregation, and delaying the formation of ACC. Recently, Schuitemaker et al.³¹ developed and characterized a new polarizable model for both calcium carbonate and organic functional groups (carboxylate and ammonium) in water based on sophisticated electrostatic treatment within the atomic multipole optimized energetics for biomolecular application (AMOEBA) force field.³² Compared with conventional rigid-ion (RI) models, not only the thermodynamics but also the solvation structures around the species are improved. Although numerous force fields have been proposed in MD studies to evaluate the effects of organic inhibitors on CaCO₃ precipitation, the applicability of these models is severely restricted to chemical environments intended for parametrization. With the ever-growing complexity of the investigated problems, some degree of compromise between accuracy and transferability is inevitable.³³

Compared with classical MD, quantum-mechanics-based density functional theory (DFT) describes interatomic interactions by directly solving the time-independent Schrödinger equation, thereby providing accurate and consistent predictions for chemical systems without relying on empirical data. Saharay et al.³⁴ examined the interactions of Ca²⁺ and CO₃²⁻ with representative amino acid functional groups, including guanidinium, acetate, acetic acid, and ethanol, in the gas phase. Delgado et al.³⁵ and Dooley et al.³⁶ incorporated implicit solvation models into their DFT studies and successfully reproduced the solvation structures of Ca²⁺ and CO₃²⁻. Lopez-Berganza et al.³⁷ and Di Tommaso et al.³⁸ explored the structural and energetic properties of CaCO₃(H₂O)_n and CaHCO₃(H₂O)_n⁺ clusters, respectively, advancing the understanding of ion association behaviors in the initial stage of CaCO₃ nucleation. Unfortunately, no DFT study has evaluated the interaction strength between Ca²⁺ and anionic groups in aqueous environments. Despite the broad application of DFT, it is worth noting that this method is fundamentally static because it computes the energies and forces for fixed configurations and cannot capture the real-time evolution of microstructures under specific conditions.

Recently, the ab initio molecular dynamics (AIMD) method has emerged as an appealing alternative for overcoming the intrinsic limitations of both classical MD and static DFT calculations.³⁹ On the one hand, interatomic interactions are determined using DFT, which goes beyond the feasibility of empirical models. On the other hand, atomic trajectories are generated by integrating Newton's laws of motion, which enables the analysis of dynamical properties. In particular, in recent studies,⁴⁰ AIMD simulations were combined with some enhanced sampling methods, e.g., umbrella sampling (US)⁴¹ and metadynamics (MetaD),⁴² to study the thermodynamics and kinetics of chemical reactions and material processes. Henzler et al.⁴³ employed US in their AIMD simulations to construct potentials of mean force between Ca²⁺ and carbonate species (CO₃²⁻ and HCO₃⁻). Li et al.⁴⁴ utilized a two-dimensional set of collective variables characterized by the coordination number in their ab initio MetaD simulations to monitor the association processes of calcium silicate aqua complexes. In studies conducted by Gale and co-workers,^45,46 ab initio MetaD simulations were carried out to calculate the ion-pairing free energy profiles between Ca²⁺ and either carbonate species (CO₃²⁻ and HCO₃⁻) or SO₄²⁻. Notably, there are also some issues associated with AIMD, including the charge delocalization error induced by exchange–correlation functionals⁴⁶ and the limited time and length scales due to high computational costs. However, integrating the AIMD results with those from MD, DFT, and experimental studies can facilitate a more comprehensive understanding of the interactions between Ca²⁺ ions and organic-scale inhibitors in aqueous solutions.

In this study, we first carried out static DFT calculations to examine the pairwise interactions between Ca²⁺ and organic scale inhibitors in implicit solvents. Given that the majority of inhibitors, including carboxylic acids,⁴⁷ phosphonic acids,⁴⁸ and sulfonates,⁴⁹ exist in deprotonated forms in slightly alkaline environments, where the formation efficiency of CaCO₃ reaches its maximum,⁵⁰ we opted to use acetate (CH₃CO₂⁻, ACE⁻), mesylate (CH₃SO₃⁻, MES⁻), and two methylphosphonates (CH₃PO₃²⁻, MP²⁻ and CH₃PO₃H⁻, MP⁻) as the basic units representing the functional groups within the polymer frameworks of organic inhibitors. Hereafter, using AIMD simulations and the MetaD method, we obtained the Helmholtz free energy surfaces for ion pairing, calculated the activation barriers for separating the contact ion pairs, and analyzed the dynamic structures of the ion complexes in real aqueous environments. The discussion of the association behaviors of calcium–organic ion pairs has provided new insights into the functional mechanisms of organic inhibitors on CaCO₃ formation.

2. Computational methods

Pairwise interactions were investigated using the quantum chemistry package Gaussian 16.⁵¹ All structures were fully optimized at the level of B3LYP/def2-TZVP^52,53 with the DFT-D3 dispersion correlation.⁵⁴ Vibration frequency analyses were performed to ensure that the obtained stable minima had no imaginary frequencies. Based on the optimized geometries, single-point energies were calculated within the framework of the expensive ωB97X-D functional⁵⁵ and the def2-QZVPP basis set^53,56 of the quadruple zeta quality. At such a high level of theory, the influence of the basis set superposition error (BSSE) can be neglected, as recommended in previous studies.^56–58 In addition, the implicit solvation model based on solute electron density (SMD)⁵⁹ with the water parameter (ε = 78.4) was used to depict the solvation status. As utilized in the parametrization of SMD, the solvation free energies of all species were calculated at the level of M05-2X/6-31G(d).^60,61 The binding free energies between Ca²⁺ and various anions were calculated following Scheme S1 and eqn (S1)–(S6). Please refer to Section S1 of the SI for more details.

Hereafter, the optimized stable complexes, i.e., aqueous Ca²⁺–CO₃²⁻, Ca²⁺–HCO₃⁻, Ca²⁺–ACE⁻, Ca²⁺–MP²⁻, Ca²⁺–MP⁻, and Ca²⁺–MES⁻ ion pairs, were selected as candidates for the explicit solvation model in further AIMD simulations to explore their association behaviors in real aqueous environments. All the AIMD simulations reported in this work were carried out using the CP2K/quickstep package^62,63 with a mixed Gaussian and plane wave (GPW) approach.⁶⁴ The generalized gradient approximation (GGA) proposed by Perdew, Burke, and Ernzerhof (PBE)⁶⁵ was applied to describe the exchange–correlation energy functional. Grimme's DFT-D3⁵⁴ dispersion correlation was used to depict the van der Waals interactions. Core electrons were described using Goedecker–Teter–Hutter (GTH) pseudopotentials.^66–68 For valence electrons, a DZVP–MOLOPT–SR–GTH basis set⁶⁴ was chosen for real space representation, and plane waves up to a cutoff energy of 500 Ry were employed to expand the electron density in reciprocal space. Brillouin zone sampling was performed only at the Γ-point.

In CP2K, periodic boundary conditions (PBCs) are employed to approximate an infinite system using a finite unit cell. However, if the investigated system carries a net charge, the long-range Coulomb interactions between the cell and its periodic replicas become divergent.⁶⁹ To address this issue, a uniform background charge, often referred to as a “jellium” background, was introduced to neutralize the net charge,⁷⁰ enabling charged systems, such as the aqueous Ca²⁺–HCO₃⁻, Ca²⁺–ACE⁻, Ca²⁺–MP⁻, and Ca²⁺–MES⁻ pairs in this work, to be studied within the periodic framework. In contrast, Gaussian 16 treats the system as an isolated finite molecule without applying PBCs, thereby avoiding divergence of the electrostatic energy. The change in the total charge is implemented by adjusting the number of electrons in the system and constructing the corresponding spin state. In Section S2, we discuss the influence of the system charge on the calculated results.

In AIMD simulations, the initial structures of bulk aqueous systems were constructed by randomly solvating an optimized ion pair in a 14.41 × 14.41 × 14.41 ų periodic cubic box with 100 water molecules (for a density of 1 g cm⁻³) using the PACKMOL code.⁷¹ Energy minimization was first performed to optimize the initial structures. Hereafter, the AIMD simulations were carried out within the isothermal–isochoric (NVT) ensemble at 330 K for 270 ps, where the last 240 ps run was used for the statistical analysis. As suggested by previous studies,^72–74 elevating the temperature by ∼30 K is an economical and effective approach that mitigates the over-structured description of liquid water by GGA functionals and mimics the real properties of aqueous solutions at room temperature. During AIMD, the time step for the integration of the equations of motion was set to 0.5 fs. A three-chain Nosé–Hoover thermostat^75–77 was used to regulate the temperature with a coupling constant of 100 fs.

The well-tempered MetaD method⁷⁸ was combined with AIMD simulations to accelerate the sampling of the configuration space. In MetaD, a chemical reaction or material process can be represented by a set of carefully selected geometric parameters, typically referred to as collective variables (CVs). During the simulation, the bias potential is constructed by periodically depositing Gaussian hills along the CVs, which forces the system away from the traps and drives it toward unexplored regions of the energy landscape. Once the CV space is sufficiently sampled, the accumulated bias potential with its sign inverted provides an estimate of the free energy surface (FES). In this study, we used the distance d between the Ca²⁺ ion and the central atom of the anionic group (C for CO₃²⁻, HCO₃⁻, and ACE⁻; P for MP²⁻ and MP⁻; and S for MES⁻) as the CV to investigate the thermodynamics of ion pairing. As suggested in previous studies,^43,79–81 the interionic distance d provides a physically unambiguous descriptor for distinguishing the bidentate, monodentate, solvent-shared, and dissociated states of an ion pair. In the resulting FES, as a function of d, the depths of the local minima quantify the stability of each state, while the swells between neighboring minima indicate the reaction barriers along the corresponding pathways. The Gaussian hills were deposited every 60 time steps (30 fs) with an initial height of 2.5 kJ mol⁻¹, a width of 0.1 Å, and a bias factor of 6. A quadratic upper wall with a force constant of 500 kJ mol⁻¹ nm⁻² was introduced at d = 6.2 Å in all systems to prevent the exploration of distances beyond half the box length. The time evolutions of CV during the whole simulation time and the convergence tests for FESs are shown in Section S3 of the SI.

In addition, we utilized a two-dimensional set of CVs characterized by the coordination number (CN) to monitor the ion association processes. Notably, the CVs defined here were used only for analysis and not for constructing the bias potential. CN(Ca–O_a) is the coordination number of the Ca²⁺ ion with O atoms from anionic groups, while CN(Ca–O_w) is the coordination number of the Ca²⁺ ion with O atoms from water molecules. As defined in the PLUMED plugin,⁸² CN has the following expression:

where r_ij represents the distance between atoms i and j; d₀ and r₀ denote the central value and acceptance distance of the switching function s_ij, respectively; and n and m are two tunable exponents that control the curvature of the function. Here, we used the same parameter settings from previous literature,^10,44,83 in which d₀ and r₀ are 2.42 and 0.4 Å, and n and m are 6 and 12, respectively. Information for CNs was generated every 30 fs, yielding 8000 sets of data for each system during the 240 ps production run. The heat maps of the CNs are shown in Section 3.3.

3. Results and discussion

3.1. Pairwise interactions

We first investigated the pairwise interactions between the Ca²⁺ ion and either the carbonate species (CO₃²⁻ and HCO₃⁻) or the functional groups of organic scale inhibitors (ACE⁻, MP²⁻, MP⁻, and MES⁻) in implicit solvents. The structures and energetics of the ion complexes were explored by initiating simulations with a number of starting geometries in which the anionic groups were placed around Ca²⁺ via a monodentate or bidentate binding mode. Hereafter, twelve explicit water molecules were introduced into the systems using the PACKMOL code⁷¹ to generate solvation shells around ion pairs. As reported in previous studies,^84–86 such a mixed explicit/implicit solvation method captures both local and bulk solvent effects and thus enhances the accuracy of thermodynamic properties. In Section S4, we discussed the influence of the number of explicit water molecules on the binding free energy. The most stable configurations of all aqueous calcium–anion complexes after optimization are illustrated in Fig. 1(a)–(f), while the corresponding atomic coordinates are included in Table S1. The structural and energetic information of several representative, less stable calcium–anion complexes is displayed in Table S2.

Fig. 1 The most stable configurations of aqueous (a) Ca²⁺–CO₃²⁻, (b) Ca²⁺–HCO₃⁻, (c) Ca²⁺–ACE⁻, (d) Ca²⁺–MP²⁻, (e) Ca²⁺–MP⁻, and (f) Ca²⁺–MES⁻ ion pairs.

Table 1 summarizes the binding free energies of the contact ion pairs. The calculated results agree well with the existing experimental and theoretical data from previous studies (see Table 1), indicating the accuracy achieved in our calculations. Generally, a more negative binding free energy is associated with a stronger interaction. Hence, among four organic functional groups, Ca²⁺ ion prefers to bind MP²⁻ with the lowest binding free energy of −14.0 kJ mol⁻¹, followed by ACE⁻, MP⁻, and MES⁻ ions of −6.2, −5.3, and +4.4 kJ mol⁻¹, respectively. For the Ca²⁺–ACE⁻ pair, Ca²⁺ is coplanar with the COO⁻ group and coordinates with the carboxylate O in a bidentate mode. For MP²⁻, MP⁻, and MES⁻ ions that contain a trigonal pyramid –XO₃(H) (X = P, S) structure, even if the phosphonate or sulfonate group is placed in a bidentate coordination with Ca²⁺ in the initial setup, one O atom from anions is detached from Ca²⁺ during the energy minimization, forming the monodentate binding mode in the end.

Table 1 Binding free energies (kJ mol⁻¹) for contact ion pairs compared with previous experimental and theoretical data

Ion pair	This work	Exp.	Theor.
a Ref. 87.b Ref. 88.c Ref. 89.d Ref. 46.e Ref. 81.f Ref. 90.g Ref. 91.h Ref. 80.i Ref. 92.j Ref. 31.
Ca^2+–CO3^2−	−17.0	−19.0^a, −18.0^b	−20.3^c, −17.3^d, −12.5^e
Ca^2+–HCO3^−	−8.1	−6.2^a, −7.5^f	−11.3^c, −8.0^d
Ca^2+–ACE^−	−6.2	−4.4 to −7.1^g, −3.2^h	−9.5^i, −8.0^j
Ca^2+–MP^2−	−14.0
Ca^2+–MP^−	−5.3
Ca^2+–MES^−	+4.4

---

To ascertain the competition among carbonate species and organic inhibitors for interactions with Ca²⁺, we examined the binding of CO₃²⁻ and HCO₃⁻ ions to Ca²⁺, as shown in Fig. 1(a) and (b), respectively. Relative to the HCO₃⁻ ion, three organic functional groups, including MP²⁻, ACE⁻, and MP⁻, exhibit a stronger or comparable affinity to Ca²⁺, indicating the potential capacity of these groups as Ca-trapping agents in near-neutral environments where HCO₃⁻ dominates the carbonate species.⁹³ With the binding of Ca²⁺ to organic inhibitors, fewer ions are involved in PNC formation, and the crystallization process is therefore retarded. In contrast, the positive binding free energy for the Ca²⁺–MES⁻ pair demonstrates that the association of the sulfonate group with Ca²⁺ is energetically unfavorable in aqueous solutions. This is one of the reasons that the sulfonate-containing additive PSS cannot practically reduce the concentration of free Ca²⁺ ions in carbonate buffer, as observed in the titration experiments by Verch et al.²¹

With increasing solution pH, CO₃²⁻ ion gradually dominates the carbonate species.⁹³ In this case, the binding energy of the Ca²⁺–CO₃²⁻ pair reaches −17.0 kJ mol⁻¹, and all interactions between Ca²⁺ and inhibitors lose their advantages. This finding can explain the diminishment in the anti-scaling performance of organic inhibitors in alkaline environments.⁹⁴ However, since a substantial portion of the phosphonic acids exist as the fully deprotonated form (–PO₃²⁻) at high pH⁴⁸ and Ca²⁺ displays a stronger affinity toward MP²⁻ than MP⁻ with a lower binding free energy, the Ca-trapping capacity of the phosphonate-containing inhibitors can be maintained to a greater extent. As evidenced by the experiments by Yao et al.,⁹⁵ the inhibition rate of the carboxylate-containing additive PAA decreases by 21%, while the phosphonate-containing additive PBTCA merely decreases by 6% in the pH ranging from 5 to 10.

Moreover, in Section S6, we performed an analysis of the charge density difference to estimate the bonding strength between Ca²⁺ and the anionic groups. Combining the values of binding free energy in Table 1 with the results of the charge density difference in Fig. S4–S6, one can notice that stronger calcium–anion interactions typically result in enhanced stability. The only exception is the Ca²⁺–MP⁻ pair, which exhibits a stronger Ca–O bonding strength but less negative binding free energy in comparison with the Ca²⁺–HCO₃⁻ pair. We attributed this phenomenon to differences in the geometries of the functional groups. Compared to the planar (bi)carbonate or carboxylate, the attachment of the trigonal pyramidal phosphonate group with Ca²⁺ squeezes the living space of water molecules in the first solvation shell of Ca²⁺ and thus weakens the stability of the contact ion pair. The bulky structure of the phosphonate group also restricts its coordination with Ca²⁺. For the Ca²⁺–MP²⁻ pair that exhibits stability second only to Ca²⁺–CO₃²⁻, the monodentate state is favored over the more tightly coordinated bidentate state.

3.2. Free energy profiles

After analyzing the pairwise interactions from the static DFT calculations, we performed ab initio MetaD simulations with explicit water molecules to investigate the thermodynamics of ion pairing in real aqueous solutions. Fig. 2 shows the free energy profiles as a function of ion separation for all systems. The free energy minima, or shoulders, at calcium–anion distances of approximately 3 and 3.5 Å correspond to the bidentate and monodentate binding modes of ion complexes, respectively. The activation barriers for the decomposition of the contact ion pairs are shown in Table 2.

Fig. 2 Free energy profiles as a function of distances between Ca²⁺ and the central atom of anionic groups for the (a) Ca²⁺–CO₃²⁻, (b) Ca²⁺–HCO₃⁻, (c) Ca²⁺–ACE⁻, (d) Ca²⁺–MP²⁻, (e) Ca²⁺–MP⁻, and (f) Ca²⁺–MES⁻ systems.

Table 2 Activation barriers (kJ mol⁻¹) for the decomposition of contact ion pairs in aqueous solutions

Ion pair	This work	AIMD	Force field
			RI	AMOEBA	ML
a Ref. 43.b Ref. 46.c Ref. 81.d Ref. 97.e Ref. 98.f Results for the Ca^2+–HCOO^− ion pair from ref. 79.g Ref. 31.
Ca^2+–CO3^2−	20.9	16.3^a, 10–16^b, 16^c	18^d	25^b	20^c
Ca^2+–HCO3^−	12.9	13.4^a	16^d	10^b
Ca^2+–ACE^−	16.8	14^e, 19^f	15^g	21^g
Ca^2+–MP^2−	17.2
Ca^2+–MP^−	11.5
Ca^2+–MES^−	5.7

---

For the Ca²⁺–CO₃²⁻ system, the bidentate configuration is 4.5 kJ mol⁻¹, which is preferred over the monodentate counterpart. For the Ca²⁺–HCO₃⁻ system, the Ca²⁺ ion favors coordinating with only one O atom from HCO₃⁻, which is consistent with the findings displayed in our static DFT calculations. According to the simulation results presented in Table 2, the total activation barrier for separating the contact Ca²⁺–CO₃²⁻ pair reaches 20.9 kJ mol⁻¹, which is 8.0 kJ mol⁻¹ higher than that of the Ca²⁺–HCO₃⁻ pair. Therefore, we can predict that CO₃²⁻ exhibits a stronger affinity to Ca²⁺ compared with HCO₃⁻, leading to the formation of a more stable precursor phase in the prenucleation stage, which benefits subsequent CaCO₃ nucleation. This finding is consistent with the experiment of Andritsos and Karabelas,⁹⁶ who reported that an increased CO₃²⁻ concentration in the pH range of 8.8–10 promotes CaCO₃ scale deposition from 2 to 12 mg cm⁻² after a 2 h run.

However, when comparing the results between the present work and previous AIMD studies, one notices that there are some clear differences, even though all calculations were performed using a similar-sized box. Henzler et al.⁴³ proposed that both the CO₃²⁻ and HCO₃⁻ ions favor coordinating with Ca²⁺ in a bidentate manner. This viewpoint is not completely consistent with our findings because the monodentate state of the Ca²⁺–HCO₃⁻ pair is effectively a point of inflection in our calculations. Although the lowest-energy configuration is debated, their study also reported a higher activation barrier for separating the Ca²⁺–CO₃²⁻ pair than for the Ca²⁺–HCO₃⁻ pair (16.3 kJ mol⁻¹ vs. 13.4 kJ mol⁻¹). This suggests that we are consistent with each other on the most essential conclusion; for example, CO₃²⁻ ion exhibits a stronger affinity to Ca²⁺ compared with HCO₃⁻. Scrutinizing the computational methodology, it is evident that the primary difference lies in the system temperature. Henzler et al.⁴³ used a temperature of 300 K, while in this study, we elevated the temperature to 330 K to mitigate the over-structured description of liquid water by GGA functionals. Here, the “over-structured” behavior indicates that the hydrogen-bond network described within the GGA framework is more persistent and less flexible than in real liquid water. Therefore, if Henzler et al. initialized the Ca²⁺–HCO₃⁻ pair in a bidentate configuration, the overly strong hydrogen-bond network may restrict the structural evolution of the contact ion pair and result in an underestimation of the free energy at this state. Compared with implementing expensive meta-GGA or hybrid functionals, elevating the temperature by ∼30 K is an economical and effective approach that mimics the real properties of aqueous solutions at room temperature.

Notably, system temperature is merely one of the factors that can potentially affect AIMD results. In the ab initio MetaD simulations conducted by Raiteri et al.,⁴⁶ the free energy landscape at 330 K suggests that the monodentate state of the Ca²⁺–CO₃²⁻ pair is ∼1 kJ mol⁻¹ more stable than the bidentate state, which is different from both our results and those of Henzler et al.⁴³ This discrepancy is partly due to the selection of CVs. Raiteri et al.⁴⁶ employed a two-dimensional set of CVs characterized by the calcium–anion distance and the coordination number of calcium by water to monitor ion pairing. Although this setup considers various solvation structures of ion pairs at the same distance, offering deeper insights into the ion association processes in theory, it poses additional challenges for the convergence of MetaD calculations. Notably, uneven sampling of different solvation structures within the same binding mode may influence the determination of the lowest-energy configuration. In the present study, we utilized the calcium–anion distance as the sole CV to provide an overall assessment of the stability of the monodentate and bidentate states. In a more recent study, Piaggi et al.⁸¹ utilized the same one-dimensional CV and obtained similar results to ours; that is, the bidentate configuration is the most stable state for the contact Ca²⁺–CO₃²⁻ pair. Apart from the parameter set in MetaD, some intrinsic weaknesses of the DFT method, such as the charge delocalization error induced by GGA functionals, have also been reported to affect ion association behaviors.⁴⁶

Given the systematic errors associated with AIMD simulations, it is reasonable for researchers to question whether the classical MD method, which utilizes a carefully parametrized force field, is a better choice for investigating the association behaviors of ion pairs. The MD method exhibits a striking advantage in computational speed, allowing the configuration space to be fully sampled within a sufficient simulation time. In addition, AIMD simulations can be performed only in small periodic cells. Even at the farthest calcium–anion distance considered in our simulations (∼6 Å), the ion pairs are not completely dissociated but exist in a solvent-shared state. Hence, AIMD simulations are incapable of mimicking the entire ion association process, as is carried out in MD simulations. However, it is worth noting that the quality of MD simulations is largely determined by the accuracy of the applied force field. In the case of the Ca²⁺–CO₃²⁻ pair, the polarizable AMOEBA force field⁴⁶ suggests that bidentate coordination is more stable than its monodentate counterpart. In the conventional RI force field,⁹⁷ the opposite order of stability is offered with a weaker minimum for the bidentate binding arrangement. Overall, our findings are more consistent with those obtained from the AMOEBA force field although the activation barrier for separating the contact Ca²⁺–CO₃²⁻ pair provided by this model is ∼4 kJ mol⁻¹ higher than our result. In 2025, Piaggi et al.⁸¹ developed a new force field based on advanced machine learning (ML) technology. In this model, the activation barrier for separating the Ca²⁺–CO₃²⁻ pair was corrected to 20 kJ mol⁻¹, closely aligning with our result of 20.9 kJ mol⁻¹.

Based on the aforementioned discussion, it can be concluded that there is no perfect approach for studying ion-pairing thermodynamics in aqueous solutions. The selection of the simulation method largely depends on the problem being investigated. In this study, we aim to devise a simple strategy for screening organic functional groups with excellent Ca-trapping capacity to offer suggestions for the development of effective scale inhibitors in the industry. Therefore, constructing a reliable and transferable force field that encompasses all organic functional groups in CaCO₃ solutions is unrealistic and uneconomical for us. In this case, we must use some general-purpose force fields to describe the calcium–organic interactions in MD simulations. A study by Kahlen et al.⁹⁸ has certified that the general-purpose force fields are far from accurate in describing the ion-pairing thermodynamics, resulting in a deviation of ∼30 kJ mol⁻¹ on binding free energy for the Ca²⁺–ACE⁻ pair. In view of this, the AIMD method is a more practical choice for us. Integrating the ab initio results with available MD and experimental studies can facilitate a more comprehensive understanding of the calcium–anion interactions. Actually, the main concerns regarding the ab initio MetaD simulations mainly focus on areas with calcium–anion distances above 5 Å,⁴⁶ where the rapidly expanding configuration space cannot be fully sampled within the limited time scale of AIMD simulations. However, as shown in Fig. 2, both the lowest energy states of the contact ion pairs and the transition states during the ion separation process are located in the area with calcium–anion distances below 4.5 Å. Therefore, the derived activation barriers should be reliable. Notably, the small peaks or minima observed at larger calcium–anion distances on free energy profiles should be considered numerical artifacts arising from insufficient sampling rather than specific configuration states. In Section S7 of the SI, we conducted three independent ab initio MetaD simulations for the Ca²⁺–CO₃²⁻ system, each starting from distinct initial structures, to confirm the consistency of our calculations within the given methodological framework.

Next, we studied the ion-pairing thermodynamics of Ca²⁺ with organic functional groups in aqueous solutions. For the Ca²⁺–ACE⁻ system, the bidentate configuration is found to be 1.5 kJ mol⁻¹ preferred over the monodentate counterpart. The nearly equal stability between bidentate and monodentate configurations aligns well with the findings presented in the AIMD study of Mendes de Oliveira⁸⁰ although they reported that the monodentate configuration is ∼1 kJ mol⁻¹ favored. Moreover, the activation barrier for separating the Ca²⁺–ACE⁻ pair is predicted to be 16.8 kJ mol⁻¹, which falls within the range of results obtained from previous studies (14–21 kJ mol⁻¹).^31,79,98 Notably, this value is 3.9 kJ mol⁻¹ higher than that of the Ca²⁺–HCO₃⁻ pair, indicating that the carboxylate group exhibits a potential calcium sequestration ability in near-neutral environments. However, as the CO₃²⁻ ion exhibits a stronger affinity to Ca²⁺ than ACE⁻, with a higher separation barrier of 20.9 kJ mol⁻¹, an increase in the concentration of CO₃²⁻ at alkaline pH diminishes the inhibition efficiency of carboxylate-containing additives.

Additionally, for phosphonate and sulfonate groups that contain a trigonal pyramid structure, their association with Ca²⁺ favors the monodentate binding mode. The bidentate state in Ca²⁺–MP²⁻ and Ca²⁺–MP⁻ pairs are still local free energy minima, which are separated by small activation barriers of 3.6 and 0.5 kJ mol⁻¹, respectively. However, for the Ca²⁺–MES⁻ system, the bidentate state essentially becomes a point of inflection rather than a distinct minimum. Regarding the decomposition processes, the activation barrier for separating the contact Ca²⁺–MP⁻ pair is only 11.5 kJ mol⁻¹, which is 1.4 kJ mol⁻¹ lower than that in the Ca²⁺–HCO₃⁻ system. Combining the less negative binding free energy between Ca²⁺ and MP⁻ shown in Section 3.1, we can confirm that the phosphonate group is actually a substandard Ca-trapping agent at near-neutral pH, whose affinity to Ca²⁺ has no obvious advantage compared with HCO₃⁻. This finding provides crucial evidence in support of the viewpoint of Duchstein et al.,²⁹ where they proposed that the inhibiting effect of the phosphonate-containing additives is mainly achieved by stabilizing the prenucleation phases in the initial stage of CaCO₃ formation, which is independent of a strong calcium sequestration ability. Although the strong interaction between Ca²⁺ and MP²⁻ enhances the decomposition barrier to 17.2 kJ mol⁻¹, it is still 3.7 kJ mol⁻¹ lower than that of the Ca²⁺–CO₃²⁻ pair.

For the Ca²⁺–MES⁻ system, the free energy profile illustrates that the Ca²⁺ ion favors keeping a distance from MES⁻. Any contact ion pairs, whether following the monodentate or bidentate binding mode, are energetically unfavorable in aqueous solutions. Although ab initio MetaD simulations exhibit limited accuracy in calculating the free energy of states where the calcium–anion distances exceed 5 Å, such a special phenomenon is sufficient to demonstrate the lability of the contact Ca²⁺–MES⁻ pair. Even if minor amounts of ion pairs can be formed occasionally, it is effortless to separate them with a decomposition barrier of 5.7 kJ mol⁻¹. As shown in Table 2, this value is the lowest among all ion pairs considered, indicating that the Ca²⁺–MES⁻ pair is almost impossible to form in the presence of carbonate species in aqueous solutions. This finding is supported by the crystallization experiments by Verch et al.,²¹ who estimated the amount of free and bound Ca²⁺ ions in solutions by means of a Ca²⁺-ion selective electrode and discovered that only 0.04 Ca²⁺ ions are bound per sulfonate group in the presence of a carbonate buffer. In contrast, one carboxylate group can adsorb 0.21 Ca²⁺ ions in the same environment, demonstrating a stronger calcium sequestration ability.

3.3. Structural analysis

During ab initio MetaD simulations, we utilized a two-dimensional set of CVs characterized by CN(Ca–O_a) and CN(Ca–O_w) to monitor the solvation structures of all ion pairs in real aqueous environments. The probability distributions of CNs and the configurations of the representative states on the heat maps are illustrated in Fig. 3. The coordinates of each state are presented in the form of (CN(Ca–O_a), CN(Ca–O_W)). As depicted in Fig. S8, the spatial distributions of Ca²⁺ around the anionic groups illustrate the effect of the anion type on the binding mode of the ion pair.

Fig. 3 Probability distributions of CN(Ca–O_a) and CN(Ca–O_w) for (a) Ca²⁺–CO₃²⁻, (b) Ca²⁺–HCO₃⁻, (c) Ca²⁺–ACE⁻, (d) Ca²⁺–MP²⁻, (e) Ca²⁺–MP⁻, and (f) Ca²⁺–MES⁻ systems. For simplicity, the atoms of the solution molecules are colored white.

At CN(Ca–O_a) equals 0, the Ca²⁺ ion tends to be surrounded by several water molecules to form the hydrated Ca(H₂O)_n²⁺ (n = 5–7) clusters, wherein the six-coordinated (0, 6) states with a nearly octahedral symmetry are commonly regarded as the starting points for the association processes of all ion pairs. For carbonate species and the carboxylate group that possess a planar structure, as shown in Fig. 3(a)–(c), their association with Ca²⁺ is initiated by introducing one O atom from anions into the first solvation shell of Ca²⁺. After going through the seven-coordinated (1, 6) states, the systems reach (1, 5) with one water molecule detached from Ca²⁺, and the total coordination number returns to six. Moreover, since the occurrence probability of (1, 5) is higher than that of (1, 6), the seven-coordinated (1, 6) states with distorted pentagonal bipyramidal structures serve as intermediates of the association processes. Therefore, the formation of the contact Ca²⁺–CO₃²⁻, Ca²⁺–HCO₃⁻, and Ca²⁺–ACE⁻ pairs with the CN(Ca–O_a) value increasing from 0 to 1 can be simplified as a water–anion exchange process, which follows an associative ligand substitution mechanism. The derived formation process of the Ca²⁺–CO₃²⁻ pair agrees well with the findings presented in the AIMD study carried out by Raiteri et al.⁴⁶ Hereafter, the Ca²⁺ ion can connect with the second O atom from the anionic groups to form the bidentate binding arrangement at CN(Ca–O_a) equals 2. In particular, for Ca²⁺–CO₃²⁻ and Ca²⁺–ACE⁻ systems, the bidentate configurations are actually ion complexes that encompass 4–5 water molecules. Notably, for contact ion pairs, the total coordination number of Ca²⁺ generally falls in the range of 6–7. This result is consistent with previous AIMD simulations of Ca²⁺ in water,³⁸ demonstrating that the presence of carbonate species and the carboxylate group does not change the coordination preference of Ca²⁺ in its first solvation shell.

However, for the phosphonate and sulfonate groups, their association with Ca²⁺ is much more complicated. As shown in Fig. 3(d)–(f), the intermediate states (1, 6) almost disappear. Instead, a large number of transition states emerge in the area where CN(Ca–O_a) and CN(Ca–O_w) vary from 0.5 to 1 and from 5.5 to 6, respectively. Based on the analysis of the dynamic structures of the calcium–organic ion complexes, we primarily ascribed this difference to the varying geometries of functional groups. Compared to the planar (bi)carbonate/carboxylate, the attachment of the trigonal pyramidal phosphonate or sulfonate group with Ca²⁺ undoubtedly disrupts the symmetry of the coordination structure to a greater extent. To accommodate the bulky phosphonate or sulfonate, the Ca²⁺ ion at (0, 6) has to simultaneously liberate one water ligand to prepare ample room for ion association. More interestingly, the formation of the contact ion pairs can potentially follow the dissociative ligand substitution mechanism. In other words, the Ca²⁺ ion discards one water ligand at first and coordinates one anion O afterward to keep the total coordination number no more than six. Actually, it has been reported that the release of water ligands from the metal ion's solvation shell can increase the translational and rotational degrees of freedom of water molecules and, therefore, enhance the entropy of the system.^88,99 As proposed by Byrne et al.,⁴⁵ such entropic terms drive the association of ions.

Structural analysis of ion complexes can also be used to explain why the phosphonate group fails to serve as a qualified calcium sequestration agent. In real aqueous environments, as the MP²⁻ or MP⁻ ion approaches Ca²⁺ and connects with it via bridging oxygen, the remaining –O or –OH moiety of the phosphonate group strongly squeezes the living space of water molecules in the first solvation shell of Ca²⁺ and even compels the Ca²⁺ ion to release more water ligands to balance the structure. As shown in Fig. 3(d) and (e), the probability of the monodentate configurations encompassing 4–5 water molecules increases significantly in comparison with the Ca²⁺–CO₃²⁻, Ca²⁺–HCO₃⁻, and Ca²⁺–ACE⁻ systems, indicating a decrease in the total coordination number of Ca²⁺. In this case, the Ca²⁺ ion has a tendency to restore coordination symmetry, and the stability of the Ca²⁺–MP²⁻ and Ca²⁺–MP⁻ pairs is, therefore, degraded. For the Ca²⁺–MES⁻ system, considering the weak interaction between Ca²⁺ and MES⁻ and the bulky structure of the sulfonate group, it is not surprising that the probability of contact ion pairs is lower than in the other systems.

4. Conclusions

In the current study, we performed static DFT calculations and ab initio MetaD simulations to systematically explore the association behaviors of Ca²⁺ ions with either carbonate species or organic inhibitors in aqueous solutions. Here, we employed acetate (ACE⁻), mesylate (MES⁻), and two methylphosphonates (MP⁻ and MP²⁻) as the basic units representing carboxylate, sulfonate, and phosphonate groups, respectively, within the polymer frameworks of organic inhibitors.

By analyzing the thermodynamic and structural properties of the calcium–anion ion pairs, we discovered that the Ca-trapping capacity of an inhibitor is determined by two factors. The first one depends on the strength of the interaction. Apparently, a strong interionic interaction is a crucial prerequisite for the stable existence of an ion complex. The second one is derived from the geometry of the functional group. For the carboxylate group, which possesses a planar structure akin to the carbonate species, its presence has less impact on the solvation structure of Ca²⁺. The monodentate and bidentate binding modes have similar stability for the contact ion pair. However, for phosphonate and sulfonate groups that contain a trigonal pyramid structure, their association with Ca²⁺ disrupts the symmetry of the coordination structure and decreases the capacity of Ca²⁺ for accommodating O. The monodentate binding arrangement is preferred for contact ion pairs.

By considering the synergistic effects of these two factors on ion pairing, we confirmed that carboxylate and phosphonate groups possess potential calcium sequestration ability, while the formation of the calcium–sulfonate pair is almost impossible in near-neutral environments, where HCO₃⁻ dominates the carbonate species. Moreover, since CO₃²⁻ exhibits the strongest affinity to Ca²⁺, an increase in the concentration of CO₃²⁻ from near-neutral to alkaline pH diminishes the Ca-trapping capacity of all organic functional groups. Because a substantial portion of the phosphonic acids exist as the fully deprotonated form (–PO₃²⁻) at high pH and Ca²⁺ displays a stronger affinity toward MP²⁻ than MP⁻ with a more negative binding free energy and a higher separation barrier, the calcium sequestration ability of the phosphonate group can be maintained to a greater extent compared with the carboxylate in alkaline environments. Notably, we discovered that all functional groups of organic inhibitors considered in this study have no obvious advantage for interactions with Ca²⁺ compared to carbonate species, indicating that Ca-trapping capacity is perhaps not the dominant factor that directly determines anti-scaling efficiency. A more comprehensive strategy is needed to evaluate the overall effects of organic inhibitors on CaCO₃ crystallization.

Overall, we provided new insights into the understanding of the inhibiting effects of organic additives on the initial stage of CaCO₃ formation by analyzing the calcium–organic interactions in aqueous environments. However, we must acknowledge the limitations of this ab initio study. First, we neglected the possible multibody interactions between the dissolved species in solutions, e.g., the synergistic effects of multiple functional groups and the interactions between organic inhibitors and larger molecular clusters (PNCs or ACC particles). Second, due to the limitations of AIMD studies on time and length scales, the impact of some macroscopic factors on anti-scaling efficiency (e.g., the concentrations of organic inhibitors) was not considered. MD simulations based on force field methods are more adequate for solving these complex problems. Our calculation results can be utilized as reference data in future research to accelerate the construction of more accurate force fields for describing the role of organic additives in CaCO₃ crystallization.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: pairwise interaction calculation, influence of the system charge on the calculated results, MetaD convergence tests, influence of the number of explicit water molecules on the calculated results, optimized structures and atomic coordinates of ion complexes, analysis of charge density difference, effects of initial structures on simulation results and spatial distributions of Ca²⁺ around the anionic groups. See DOI: https://doi.org/10.1039/d5cp03876f.

Acknowledgements

Y. L., H. Z., and H. B. Z. gratefully acknowledge the support of the Natural Sciences and Engineering Research Council of Canada under the Alliance program (ALLRP 557113 - 20).

References

1. C. Li, C. Liu and W. Xu, et al., Formation mechanisms and supervisory prediction of scaling in water supply pipelines: A review, Water Res., 2022, 222, 118922.
2. Y. Guo, Z. Xu and S. Guo, et al., Selection of anode materials and optimization of operating parameters for electrochemical water descaling, Sep. Purif. Technol., 2021, 261, 118304.
3. S. Qazvini, A. Golkari and A. Azdarpour, et al., Experimental and modelling approach to investigate the mechanisms of formation damage due to calcium carbonate precipitation in carbonate reservoirs, J. Pet. Sci. Eng., 2021, 205, 108801.
4. A. A. Olajire, A review of oilfield scale management technology for oil and gas production, J. Pet. Sci. Eng., 2015, 135, 723–737.
5. L. D. Tijing, Y. C. Woo and J.-S. Choi, et al., Fouling and its control in membrane distillation—A review, J. Membr. Sci., 2015, 475, 215–244.
6. M. Gryta, The influence of magnetic water treatment on CaCO₃ scale formation in membrane distillation process, Sep. Purif. Technol., 2011, 80, 293–299.
7. K. Ali and M. I. Hassan Ali, A numerical study of CaCO₃ deposition on the membrane surface of direct contact membrane distillation, Desalination, 2024, 576, 117364.
8. A. Li, H. Zhang and Q. Liu, et al., Effects of chemical inhibitors on the scaling behaviors of calcite and the associated surface interaction mechanisms, J. Colloid Interface Sci., 2022, 618, 507–517.
9. Y. Li, H. Zeng and H. Zhang, Influence of impurity metal doping on calcite growth: A first-principles study, Appl. Surf. Sci., 2023, 637, 157927.
10. Y. Li, J. Zhang and H. Zeng, et al., Ion association behaviors in the initial stage of calcium carbonate formation: An ab initio study, J. Chem. Phys., 2024, 161, 014503.
11. M. Chaussemier, E. Pourmohtasham and D. Gelus, et al., State of art of natural inhibitors of calcium carbonate scaling. A review article, Desalination, 2015, 356, 47–55.
12. A. L. Kavitha, T. Vasudevan and H. G. Prabu, Evaluation of synthesized antiscalants for cooling water system application, Desalination, 2011, 268, 38–45.
13. D. Liu, W. Dong and F. Li, et al., Comparative performance of polyepoxysuccinic acid and polyaspartic acid on scaling inhibition by static and rapid controlled precipitation methods, Desalination, 2012, 304, 1–10.
14. J. Yang, Z. Hu and Z. Wang, et al., Preparation and scale inhibition performance of modified polyaspartic acid (M-PASP), J. Mol. Liq., 2024, 401, 124712.
15. Y. Ji, Y. Chen and J. Le, et al., Highly effective scale inhibition performance of amino trimethylenephosphonic acid on calcium carbonate, Desalination, 2017, 422, 165–173.
16. M. F. Mady, S. Abdel-Azeim and M. A. Kelland, Antiscaling evaluation and quantum chemical studies of nitrogen-free organophosphorus compounds for oilfield scale management, Ind. Eng. Chem. Res., 2021, 60, 12175–12188.
17. Z. Kiaei and A. Haghtalab, Experimental study of using Ca-DTPMP nanoparticles in inhibition of CaCO₃ scaling in a bulk water process, Desalination, 2014, 338, 84–92.
18. X.-R. Xu, A.-H. Cai and R. Liu, et al., The roles of water and polyelectrolytes in the phase transformation of amorphous calcium carbonate, J. Cryst. Growth, 2008, 310, 3779–3787.
19. Y. Li, H. Zeng and H. Zhang, Atomistic simulations of nucleation and growth of CaCO₃ with the influence of inhibitors: A review, Mater. Genome Eng. Adv., 2023, 1, e4.
20. H. Du and E. Amstad, Water: How Does It Influence the CaCO₃ Formation?, Angew. Chem., Int. Ed., 2020, 59, 1798–1816.
21. A. Verch, D. Gebauer and M. Antonietti, et al., How to control the scaling of CaCO₃: a “fingerprinting technique” to classify additives, Phys. Chem. Chem. Phys., 2011, 13, 16811–16820.
22. K. Cui, C. Li and B. Yao, et al., Synthesis and evaluation of an environment-friendly terpolymer CaCO₃ scale inhibitor for oilfield produced water with better salt and temperature resistance, J. Appl. Polym. Sci., 2020, 137, 48460.
23. N. H. de Leeuw, S. C. Parker and K. H. Rao, Modeling the competitive adsorption of water and methanoic acid on calcite and fluorite surfaces, Langmuir, 1998, 14, 5900–5906.
24. T. G. Cooper and N. H. de Leeuw, Adsorption of methanoic acid onto the low-index surfaces of calcite and aragonite, Mol. Simul., 2002, 28, 539–556.
25. N. H. de Leeuw and T. G. Cooper, A computer modeling study of the inhibiting effect of organic adsorbates on calcite crystal growth, Cryst. Growth Des., 2004, 4, 123–133.
26. A. Pavese, M. Catti and S. C. Parker, et al., Modelling of the thermal dependence of structural and elastic properties of calcite, CaCO₃, Phys. Chem. Miner., 1996, 23, 89–93.
27. P. Dauber-Osguthorpe, V. A. Roberts and D. J. Osguthorpe, et al., Structure and energetics of ligand binding to proteins: Escherichia coli dihydrofolate reductase-trimethoprim, a drug-receptor system, Proteins: Struct., Funct., Bioinf., 1988, 4, 31–47.
28. N. H. de Leeuw and S. C. Parker, Molecular-dynamics simulation of MgO surfaces in liquid water using a shell-model potential for water, Phys. Rev. B:Condens. Matter Mater. Phys., 1998, 58, 13901–13908.
29. P. Duchstein, P. I. Schodder and S. Leupold, et al., Small-molecular-weight additives modulate calcification by interacting with prenucleation clusters on the molecular level, Angew. Chem., Int. Ed., 2022, 61, e202208475.
30. W. D. Cornell, P. Cieplak and C. I. Bayly, et al., A second generation force field for the simulation of proteins, nucleic acids, and organic molecules, J. Am. Chem. Soc., 1995, 117, 5179–5197.
31. A. Schuitemaker, J. Aufort and K. B. Koziara, et al., Simulating the binding of key organic functional groups to aqueous calcium carbonate species, Phys. Chem. Chem. Phys., 2021, 23, 27253–27265.
32. J. W. Ponder, C. Wu and P. Ren, et al., Current status of the AMOEBA polarizable force field, J. Phys. Chem. B, 2010, 114, 2549–2564.
33. N. Yao, X. Chen and Z.-H. Fu, et al., Applying classical, ab initio, and machine-learning molecular dynamics simulations to the liquid electrolyte for rechargeable batteries, Chem. Rev., 2022, 122, 10970–11021.
34. M. Saharay and R. J. Kirkpatrick, Ab initio and metadynamics studies on the role of essential functional groups in biomineralization of calcium carbonate and environmental situations, Phys. Chem. Chem. Phys., 2014, 16, 26843–26854.
35. A. A. A. Delgado, D. Sethio and I. Munar, et al., Local vibrational mode analysis of ion–solvent and solvent–solvent interactions for hydrated Ca²⁺ clusters, J. Chem. Phys., 2020, 153, 224303.
36. M. R. Dooley and S. Vyas, Role of explicit solvation and level of theory in predicting the aqueous reduction potential of carbonate radical anion by DFT, Phys. Chem. Chem. Phys., 2025, 27, 6867–6874.
37. J. A. Lopez-Berganza, Y. Diao and S. Pamidighantam, et al., Ab initio studies of calcium carbonate hydration, J. Phys. Chem. A, 2015, 119, 11591–11600.
38. D. D. Tommaso and N. H. de Leeuw, The onset of calcium carbonate nucleation: A density functional theory molecular dynamics and hybrid microsolvation/continuum study, J. Phys. Chem. B, 2008, 112, 6965–6975.
39. A. A. Hassanali, J. Cuny and V. Verdolino, et al., Aqueous solutions: State of the art in ab initio molecular dynamics, Philos. Trans. R. Soc., A, 2014, 372, 20120482.
40. D. Ray and M. Parrinello, Kinetics from metadynamics: principles, applications, and outlook, J. Chem. Theory Comput., 2023, 19, 5649–5670.
41. G. M. Torrie and J. P. Valleau, Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling, J. Comput. Phys., 1977, 23, 187–199.
42. A. Laio and M. Parrinello, Escaping free-energy minima, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 12562–12566.
43. K. Henzler, E. O. Fetisov and M. Galib, et al., Supersaturated calcium carbonate solutions are classical, Sci. Adv., 2018, 4, eaao6283.
44. Y. Li, H. Pan and Z. Li, Ab initio metadynamics simulations on the formation of calcium silicate aqua complexes prior to the nuleation of calcium silicate hydrate, Cem. Concr. Res., 2022, 156, 106767.
45. E. H. Byrne, P. Raiteri and J. D. Gale, Computational insight into calcium–sulfate ion pair formation, J. Phys. Chem. C, 2017, 121, 25956–25966.
46. P. Raiteri, A. Schuitemaker and J. D. Gale, Ion pairing and multiple ion binding in calcium carbonate solutions based on a polarizable AMOEBA force field and ab initio molecular dynamics, J. Phys. Chem. B, 2020, 124, 3568–3582.
47. M. Müller, L. Wirth and B. Urban, Determination of the carboxyl dissociation degree and pK_a value of mono and polyacid solutions by FTIR titration, Macromol. Chem. Phys., 2021, 222, 2000334.
48. N. Aliouane, S. Chafaa and T. Douadi, et al., Novel polydentate phosphonic acids: Protonation and stability constants of complexes with Fe(III) and Cu(II) in aqueous medium, Heteroat. Chem., 2010, 21, 51–62.
49. K. C. Kwong, M. M. Chim and E. H. Hoffmann, et al., Chemical transformation of methanesulfonic acid and sodium methanesulfonate through heterogeneous OH oxidation, ACS Earth Space Chem., 2018, 2, 895–903.
50. M. D. G. de Luna, A. S. Sioson and A. E. S. Choi, et al., Operating pH influences homogeneous calcium carbonate granulation in the frame of CO₂ capture, J. Cleaner Prod., 2020, 272, 122325.
51. M. J. Frisch, G. W. Trucks, H. B. Schlegel, et al., Gaussian 16 Rev. C.01, Wallingford, CT, 2016.
52. A. D. Becke, Density-functional thermochemistry. III. The role of exact exchange, J. Chem. Phys., 1993, 98, 5648–5652.
53. F. Weigend and R. Ahlrichs, Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
54. S. Grimme, J. Antony and S. Ehrlich, et al., A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys., 2010, 132, 154104.
55. J.-D. Chai and M. Head-Gordon, Long-range corrected hybrid density functionals with damped atom–atom dispersion corrections, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
56. F. Weigend, F. Furche and R. Ahlrichs, Gaussian basis sets of quadruple zeta valence quality for atoms H–Kr, J. Chem. Phys., 2003, 119, 12753–12762.
57. R. Knochenmuss, R. K. Sinha and A. Poblotzki, et al., Intermolecular dissociation energies of hydrogen-bonded 1-naphthol complexes, J. Chem. Phys., 2018, 149, 204311.
58. S. P. Veccham and M. Head-Gordon, Density functionals for hydrogen storage: Defining the H2Bind275 test set with ab initio benchmarks and assessment of 55 functionals, J. Chem. Theory Comput., 2020, 16, 4963–4982.
59. A. V. Marenich, C. J. Cramer and D. G. Truhlar, Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions, J. Phys. Chem. B, 2009, 113, 6378–6396.
60. M. M. Francl, W. J. Pietro and W. J. Hehre, et al., Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements, J. Chem. Phys., 1982, 77, 3654–3665.
61. Y. Zhao, N. E. Schultz and D. G. Truhlar, Design of density functionals by combining the method of constraint satisfaction with parametrization for thermochemistry, thermochemical kinetics, and noncovalent interactions, J. Chem. Theory Comput., 2006, 2, 364–382.
62. J. VandeVondele, M. Krack and F. Mohamed, et al., Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach, Comput. Phys. Commun., 2005, 167, 103–128.
63. T. D. Kühne, M. Iannuzzi and M. Del Ben, et al., CP2K: An electronic structure and molecular dynamics software package – Quickstep: Efficient and accurate electronic structure calculations, J. Chem. Phys., 2020, 152, 194103.
64. J. VandeVondele and J. Hutter, Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases, J. Chem. Phys., 2007, 127, 114105.
65. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 1996, 77, 3865–3868.
66. S. Goedecker, M. Teter and J. Hutter, Separable dual-space Gaussian pseudopotentials, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
67. C. Hartwigsen, S. Goedecker and J. Hutter, Relativistic separable dual-space Gaussian pseudopotentials from H to Rn, Phys. Rev. B:Condens. Matter Mater. Phys., 1998, 58, 3641–3662.
68. M. Krack, Pseudopotentials for H to Kr optimized for gradient-corrected exchange-correlation functionals, Theor. Chem. Acc., 2005, 114, 145–152.
69. P. A. Schultz, Local electrostatic moments and periodic boundary conditions, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 60, 1551–1554.
70. G. Makov and M. C. Payne, Periodic boundary conditions in ab initio calculations, Phys. Rev. B:Condens. Matter Mater. Phys., 1995, 51, 4014–4022.
71. L. Martínez, R. Andrade and E. G. Birgin, et al., PACKMOL: A package for building initial configurations for molecular dynamics simulations, J. Comput. Chem., 2009, 30, 2157–2164.
72. R. A. DiStasio, Jr., B. Santra and Z. Li, et al., The individual and collective effects of exact exchange and dispersion interactions on the ab initio structure of liquid water, J. Chem. Phys., 2014, 141, 084502.
73. Y.-L. S. Tse, C. Knight and G. A. Voth, An analysis of hydrated proton diffusion in ab initio molecular dynamics, J. Chem. Phys., 2015, 142, 014104.
74. L. Ruiz Pestana, O. Marsalek and T. E. Markland, et al., The quest for accurate liquid water properties from first principles, J. Phys. Chem. Lett., 2018, 9, 5009–5016.
75. S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys., 1984, 81, 511–519.
76. W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A:At., Mol., Opt. Phys., 1985, 31, 1695–1697.
77. G. J. Martyna, M. L. Klein and M. Tuckerman, Nosé–Hoover chains: The canonical ensemble via continuous dynamics, J. Chem. Phys., 1992, 97, 2635–2643.
78. A. Barducci, G. Bussi and M. Parrinello, Well-tempered metadynamics: A smoothly converging and tunable free-energy method, Phys. Rev. Lett., 2008, 100, 020603.
79. T. Martinek, E. Duboué-Dijon and Š. Timr, et al., Calcium ions in aqueous solutions: Accurate force field description aided by ab initio molecular dynamics and neutron scattering, J. Chem. Phys., 2018, 148, 222813.
80. D. Mendes de Oliveira, S. R. Zukowski and V. Palivec, et al., Binding of divalent cations to acetate: Molecular simulations guided by Raman spectroscopy, Phys. Chem. Chem. Phys., 2020, 22, 24014–24027.
81. P. M. Piaggi, J. D. Gale and P. Raiteri, Ab initio machine-learning simulation of calcium carbonate from aqueous solutions to the solid state, Proc. Natl. Acad. Sci. U. S. A., 2025, 122, e2415663122.
82. G. A. Tribello, M. Bonomi and D. Branduardi, et al., PLUMED 2: New feathers for an old bird, Comput. Phys. Commun., 2014, 185, 604–613.
83. Z. Liu, T. He and X. Zhi, et al., Critical role of alkali-metal ions in a β-dicalcium silicate hydration reaction: Well-tempered metadynamics, J. Phys. Chem. C, 2024, 128, 8206–8218.
84. C. P. Kelly, C. J. Cramer and D. G. Truhlar, SM6:  A density functional theory continuum solvation model for calculating aqueous solvation free energies of neutrals, ions, and solute–water clusters, J. Chem. Theory Comput., 2005, 1, 1133–1152.
85. C. P. Kelly, C. J. Cramer and D. G. Truhlar, Adding explicit solvent molecules to continuum solvent calculations for the calculation of aqueous acid dissociation constants, J. Phys. Chem. A, 2006, 110, 2493–2499.
86. V. S. Bryantsev, M. S. Diallo and W. A. Goddard Iii, Calculation of solvation free energies of charged solutes using mixed cluster/continuum models, J. Phys. Chem. B, 2008, 112, 9709–9719.
87. L. N. Plummer and E. Busenberg, The solubilities of calcite, aragonite and vaterite in CO₂-H₂O solutions between 0 and 90 °C, and an evaluation of the aqueous model for the system CaCO₃-CO₂-H₂O, Geochim. Cosmochim. Acta, 1982, 46, 1011–1040.
88. M. Kellermeier, P. Raiteri and J. K. Berg, et al., Entropy drives calcium carbonate ion association, ChemPhysChem, 2016, 17, 3535–3541.
89. R. Demichelis, P. Raiteri and J. D. Gale, et al., Stable prenucleation mineral clusters are liquid-like ionic polymers, Nat. Commun., 2011, 2, 590.
90. O. I. Martynova, L. G. Vasina and S. A. Pozdniakova, Dissociation constants of certain scale forming salts, Desalination, 1974, 15, 259–265.
91. E. L. Shock and C. M. Koretsky, Metal-organic complexes in geochemical processes: Estimation of standard partial molal thermodynamic properties of aqueous complexes between metal cations and monovalent organic acid ligands at high pressures and temperatures, Geochim. Cosmochim. Acta, 1995, 59, 1497–1532.
92. P. Raiteri, R. Demichelis and J. D. Gale, et al., Exploring the influence of organic species on pre- and post-nucleation calcium carbonate, Faraday Discuss., 2012, 159, 61–85.
93. D. Gebauer, A. Völkel and H. Cölfen, Stable prenucleation calcium carbonate clusters, Science, 2008, 322, 1819–1822.
94. A. Li, J. Chang and T. Shui, et al., Probing interaction forces associated with calcite scaling in aqueous solutions by atomic force microscopy, J. Colloid Interface Sci., 2023, 633, 764–774.
95. Q. Yao, R. Zhan and H. Ren, et al., Carboxyl or phosphate functionalization polyamidoamine dendrimer efficient scale inhibitor: Preparation and properties, J. Mol. Struct., 2022, 1252, 132130.
96. N. Andritsos and A. J. Karabelas, The influence of particulates on CaCO₃ scale formation, J. Heat Transfer, 1999, 121, 225–227.
97. P. Raiteri, R. Demichelis and J. D. Gale, Thermodynamically consistent force field for molecular dynamics simulations of alkaline-earth carbonates and their aqueous speciation, J. Phys. Chem. C, 2015, 119, 24447–24458.
98. J. Kahlen, L. Salimi and M. Sulpizi, et al., Interaction of charged amino-acid side chains with ions: An optimization strategy for classical force fields, J. Phys. Chem. B, 2014, 118, 3960–3972.
99. X. Yang, M. Ji and C. Zhang, et al., Physical insight into the entropy-driven ion association, J. Comput. Chem., 2022, 43, 1621–1632.

---