Theoretical Studies of a Silica Functionalized Acrylamide for Calcium Scale Inhibition

The calcium carbonate (CaCO3) scale is one of the most common oilfield scales and oil and gas production bane. CaCO3 scale can lead to a sudden halt in production or, worst-case scenario, accidents; therefore, CaCO3 scale formation prevention is essential for the oil and gas industry. Scale inhibitors are chemicals that can mitigate this problem. We used two popular theoretical techniques in this study: Density Functional Theory (DFT) and Ab Initio Molecular Dynamics (AIMD). The objective was to investigate the inhibitory abilities of mixed oligomers, specifically acrylamide functionalized silica (AM-Silica). DFT studies indicate that Ca2+ does not bind readily to acryl acid and acrylamide; however, it has a good binding affinity with PAM and Silica functionalized PAM. The highest binding affinity occurs in the silica region and not the –CONH functional groups. AIMD calculations corroborate the DFT studies, as observed from the MD trajectory that Ca2+ binds to PAM-Silica by forming bonds with silicon; however, Ca2+ initially forms a bond with silicon in the presence of water molecules. This bonding does not last long, and it subsequently bonds with the oxygen atoms present in the water molecule. PAM-Silica is a suitable calcium scale inhibitor because of its high binding affinity with Ca2+. Theoretical studies (DFT and AIMD) have provided atomic insights on how AM-Silica could be used as an efficient scale inhibitor.


Introduction
Scale deposition is the bane of the oil and gas industry. These oilfield scales often contain calcium carbonates and sulfates, including barium sulfate, iron sulfides, and many others [1,2]. Calcium carbonate (CaCO 3 ) is one of the most frequently occurring scales [3]. This scaling can occur during oil production, stimulation, and even transportation. CaCO 3 scale formation occurs when Ca 2+ is supersaturated with CO 3 2− or HCO 3 − carbonates, as illustrated in Equations (1) and (2): Different factors cause ion supersaturation. Seawater can react with the formation water in offshore fields such as the North Sea, leading to CaCO 3 deposits that hinder oil production, which is a significant problem [4]. The High Pressure and High Temperature (HPHT) operating conditions during the production process also contribute to scaling deposition [5] The consequences of calcium scale deposition can bring production capacity to a halt in a few hours, with a very high treatment cost [6]. A notable example is the North Sea DFT calculations to better understand how they interact. Combining these two techniques gives an accurate and detailed understanding of Ca 2+ adsorption by the studied polymers, including silica functionalized PAM.

Quantum Chemical Calculations
The Gaussian 09 code [30] was used for all DFT quantum calculations. All calculations were completed using the B3LYP level of theory with the TZVP basis set. The level of theory B3LYP is robust and well-known in theoretical scale inhibitor studies since it gives reliable results with reasonable computational costs [31,32]. We chose TZVP (triple-zeta valence with polarization) to minimize basis set superposition errors, despite it being computationally expensive compared to other reliable basis sets such as 631g-(d) [33,34]. All calculations were performed using water as a solvent with the PCM-SCRF (Polarizable Continuum Model-Self-Consistent Reaction Field). Vibrational frequencies were calculated, and the absence of negative values implied true minima and no imaginary frequency present. Gauss View 5.0 [35,36] was used to visualize molecular electrostatic and frontier molecular orbital maps. The Quantum Theory Atom In Molecules (QTAIM) developed by Bader [37] was done on the optimized structures from Gaussian 09 to validate the binding interactions with the aid of Multiwfn program [38]. QTAIM analysis has proven to be a useful tool in analyzing bond interactions [39]. The VMD program [40] was employed for the visualization of the QTAIM Bond Critical Points (BCP).

AIMD Calculations
The Ab Initio Molecular Dynamics (AIMD) calculations were performed with VASP code (v5.4.4) [41][42][43] using periodic boundary conditions (PBC). The Projected-Augmented-Waves Perdew-Burke-Ernzerhof (PAW-PBE) pseudopotentials [44] under the Generalized Gradient Approximation (GGA) exchange form of correlation were used for all elements in the system studied. The temperature was set to begin at 298 K and increased to 373 K at a 1.25 × 10 −4 K/ps rate since many operating conditions in the oil industry are approximately 373 K. The number of ionic steps and time steps were 20,000 and 0.5 fs, respectively, with a total of 10 ps simulation time. However, in following best practices, the first 2 ps was not considered in the trajectory to enable the system get equilibrated. The inclusion of Grimme's DFT+D3 [45][46][47] was included in the calculation has it has shown to improve the accuracy of the results by describing the dispersion forces in the studied system.
The implicit solvent was considered by implementing the VASPsol models [27,48]. This was done to simulate a realistic model. This method has proved to be reliable as it is better than the Finite-Difference Poisson Boltzmann (FDPB) approach and on par with the Solvation Model Density (SMD) approach [49]. Moreover, this method has also been applied to polymer systems [50] Though it is important to note that the perfect simulation of reality based on quantum-mechanical description remains a tricky situation [51]. A 1 × 1 × 1 Gamma (γ) grid was employed for the k-point and the planewave energy cut-off was set to the default which is the largest ENMAX in the POTCAR. The Quantum ATK Virtual NanoLab builder and MD analyzer [52] were used to build and visualize the model.

Binding Affinities
The Ca 2+ ion was complexed with the DFT optimized structures of six oligomers ( Figure 1). The binding affinities of these oligomers were calculated using Equation (3). The first three sets of oligomers were a mix of acryl acid (AA) and acrylamide (AM) monomers. The first oligomer contained 70% acryl acid. Three of the four monomer units were acryl acid, and the last was acrylamide (AA-AA-AA-AM). The second oligomer had 50% each acryl acid and acrylamide (AA-AA-AM-AM), while the third oligomer contained three

Binding Affinities
The Ca 2+ ion was complexed with the DFT optimized structures of six oligomers (Figure 1). The binding affinities of these oligomers were calculated using Equation (3). The first three sets of oligomers were a mix of acryl acid (AA) and acrylamide (AM) monomers. The first oligomer contained 70% acryl acid. Three of the four monomer units were acryl acid, and the last was acrylamide (AA-AA-AA-AM). The second oligomer had 50% each acryl acid and acrylamide (AA-AA-AM-AM), while the third oligomer contained three units of acrylamide and one unit of acryl acid (AA-AM-AM-AM). The mixed oligomers all had positive values, denoting that the binding affinity with the Ca 2+ ion was weak. The remaining three oligomers ( Figure 2) had negative values, denoting good binding affinity with the Ca 2+ ion. The fourth oligomer comprised four acrylamide (PAM) units The fifth oligomer of the hydrolyzed_PAM had a binding affinity of −289.03 kcal/mol, which is more than double the binding affinity for PAM alone; however, this double-fold increase is possibly due to the bond type formed. Two bonds are formed in hydrolyzed_PAM; however, the two bonds formed are Ca-O bonds with bond lengths of 2.220 Å and 2.205 Å, unlike PAM, where one of the bonds formed was Ca-N. These results imply that Ca-O bonds are shorter and stronger than Ca-N bonds. A silica-functionalized AM was the last oligomer studied; however, the binding affinities were studied at two positions. The first position studied was the binding affinity of PAM_Silica with Ca 2+ ions using the oxygen ing two Ca-O bonds with 2.185 Å and 2.420 Å bond lengths, one Ca-N bond (2.868 Å), and a Ca-Si bond (3.039 Å). The presence of these bonds is responsible for its high binding affinity, implying that functionalizing PAM with silica will significantly increase its binding affinity to Ca 2+ ions. These results agree with what has been observed, where the adsorption of SiO2-PAM on Ca 2+ is a spontaneous process, primarily due to chemical interaction and chelation [53].  Nevertheless, its relatively high binding affinity compared to ordinary PAM (−106.05 kcal/mol) is possible due to intra-hydrogen bonding, which contributes positively to its binding affinity. The second position studied the binding of Ca 2+ to the silica part of the oligomer (Figure 2d). The binding affinity calculated was the highest amongst all studied oligomers (−313.43 kcal/mol). Ca 2+ was within the bonding distance of four atoms, including two Ca-O bonds with 2.185 Å and 2.420 Å bond lengths, one Ca-N bond (2.868 Å), and a Ca-Si bond (3.039 Å). The presence of these bonds is responsible for its high binding affinity, implying that functionalizing PAM with silica will significantly increase its binding affinity to Ca 2+ ions. These results agree with what has been observed, where the adsorption of SiO 2 -PAM on Ca 2+ is a spontaneous process, primarily due to chemical interaction and chelation [53].

QTAIM Analysis
The QTAIM analysis confirmed the presence of 3, 2, 2, 3 BCPs (3, −1) for Ca-PAM, Ca-PAM-hydro, Ca-PAM-Silica I and Ca-PAM Silica_2 respectively ( Table 1). The DFT calculated energy values and binding energies for all the studied complexes can be seen in the supplemental information's Table S1. The BCPs are illustrated in red circles in Figure 3. Usually, low electron density (ρ(r)) −0.001 to −0.059 a.u, positive values of the Laplacian (∇ 2 ρ(r)), and zero or near zero energy density (H b )-0.000 to 0.005 a.u. at BCPs denote typical non-covalent interactions [54,55]. The interaction energies which are defined based on the works of Espinosa et al. [56] and Vener et al. [57] for E int a and E int b respectively show

Quantum Chemical Analysis
The electronic structure calculations from DFT give rise to other interesting studies, which provide atomistic insights into how molecules are formed. These insights include the molecular electrostatic potential (MEP) and frontier molecular orbitals (FMO) maps. The MEP denotes the charge distribution within the molecule, designated by yellow, orange, and red, with red being the most electronegative. Regions of electropositivity are represented as green, blue, and violet, with violet corresponding to the most electropositive region. The MEP maps of the silica functionalized PAM before and after binding with Ca 2+ are compared in Figure 4. The MEP before the binding is primarily electronegative, with the oxygen atoms appearing in the yellow region except for where silicon is located, which seems very light yellow. This change in electronegativity is possibly a neutralization effect by the silicon atom since it is less electronegative than the surrounding oxygen atoms, which are highly electronegative.
A small blue region denotes electropositivity on the hydrogen atom of the -CONH functional group preceding silica. This electropositivity is possibly due to the close proximity of the silicon atom to the -CONH functional group. This group is also coupled with the many carbon and hydrogen atoms that are less electronegative; therefore, the high electronegativity of the nitrogen and oxygen atoms of the -CONH group is neutralized. The second position has the highest binding affinity, involving a Ca-Si bond. The MEP map for Ca-PAM-Silica indicates a solid electropositive (deep blue) region at the location of Ca 2+ and its surrounding atoms (silicon and oxygen), indicating that Ca 2+ is electropositive and forms strong bonds with the surrounding oxygen atoms, dominating the oxygen atoms' Polymers 2022, 14, 2333 7 of 14 electronegativity. A yellow region is located where the only oxygen atom (SiO 4 − ) is not binding to Ca 2+ , indicating that this oxygen atom still retains its high electronegativity. The electronegativity is somewhat neutralized; therefore, the color is yellow instead of red due to its nearness to Ca 2+ . The rest of the molecule is primarily light blue, denoting an overall electropositive charge.
are illustrated in red circles in Figure 3. Usually, low electron density (ρ(r)) −0.001 to −0.059 a.u, positive values of the Laplacian (∇ 2 ρ(r)), and zero or near zero energy density (Hb)-0.000 to 0.005 a.u. at BCPs denote typical non-covalent interactions [54,55]. The interaction energies which are defined based on the works of Espinosa et al. [56] and Vener et al. [57] for Eint a and Eint b respectively show that non-covalent interactions exist between Ca-H in Ca-PAM and Ca-N in Ca-PAM_Slica1 and Ca-PAM_Silica 2. However, the high interaction energy observed in Ca-PAM_Silica 2 implies a degree of covalency and this corroborates the high binding energy observed from the DFT calculations. The trend in both the QTAIM interactions enegies correlates with the same trend observed in the binding energy calculation, with Ca-PAM_Silica2 having the highest in both cases. There is no significant difference between the Eint a and Eint b except in Ca-PAM_Silica 2, where Eint a is about 3 kcal/mol higher than Eint b .    The Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO) make up the Frontier Molecular (FMO) Orbitals. These orbitals give insight into the electron delocalization of the electrons in the molecule. The PAM HOMO map indicates that the electrons are delocalized at the -CONH end of the molecule (Figure 5a). Electron delocalization upon binding with Ca 2+ occurs where Ca 2+ bonds with the silicon and oxygen atoms, denoting electron exchange (Figure 5b). The LUMO of PAM-Silica (Figure 6a) reveals that the electrons are mostly delocalized at the silica region, similar to the HOMO of PAM-Silica when bound to Ca 2+ (Figure 5b). This phenomenon indicates that the HOMO and LUMO of PAM-Silica occur at both molecule ends. The electrons are delocalized at the atoms close to the calcium bonded region when the LUMO of PAM-Silica is bound to Ca 2+ (Figure 6b). The number of lobes is higher than ordinary PAM-silica alone in the HOMO and LUMO maps of the Ca 2+ bonded PAM-Silica, implying that electrons are more delocalized when PAM-Silica bounds to Ca 2+ .  Figure  5a). Electron delocalization upon binding with Ca 2+ occurs where Ca 2+ bonds with the silicon and oxygen atoms, denoting electron exchange (Figure 5b). The LUMO of PAM-Silica (Figure 6a) reveals that the electrons are mostly delocalized at the silica region, similar to the HOMO of PAM-Silica when bound to Ca 2+ (Figure 5b). This phenomenon indicates that the HOMO and LUMO of PAM-Silica occur at both molecule ends. The electrons are delocalized at the atoms close to the calcium bonded region when the LUMO of PAM-Silica is bound to Ca 2+ (Figure 6b). The number of lobes is higher than ordinary PAM-silica alone in the HOMO and LUMO maps of the Ca 2+ bonded PAM-Silica, implying that electrons are more delocalized when PAM-Silica bounds to Ca 2+ .

AIMD Analysis
Ab Initio Molecular Dynamics (AIMD) further establishes how PAM-Silica would interact with calcium ions by observing the time evolution in the trajectory. Two AIMD systems were studied: the first system studied the interaction between PAM-Silica and Ca 2+ ions. In contrast, the second system was similar to the first system but included explicit water molecules. The MD trajectory analysis confirmed earlier observations from the DFT studies that Ca 2+ would preferentially bind to the silica region instead of the -CONH functional groups (Figure 7b). Ca forms Ca-O and Ca-Si bonds, as observed earlier from DFT. The radial distribution function (RDF), which is the probability of finding an atom at a distance r from another tagged atom, was analyzed between Ca-O and Ca-Si atoms (Figure 8). The RDF indicated that both Ca-Si and Ca-O were within bonding distance, less than 3 Å; however, Ca-Si has a more prominent and higher peak, implying that the probability of occurrence for the Ca-Si bond is higher than the Ca-O bonds.

AIMD Analysis
Ab Initio Molecular Dynamics (AIMD) further establishes how PAM-Silica would interact with calcium ions by observing the time evolution in the trajectory. Two AIMD systems were studied: the first system studied the interaction between PAM-Silica and Ca 2+ ions. In contrast, the second system was similar to the first system but included explicit water molecules. The MD trajectory analysis confirmed earlier observations from the DFT studies that Ca 2+ would preferentially bind to the silica region instead of the -CONH functional groups (Figure 7b). Ca forms Ca-O and Ca-Si bonds, as observed earlier from DFT. The radial distribution function (RDF), which is the probability of finding an atom at a distance r from another tagged atom, was analyzed between Ca-O and Ca-Si atoms (Figure 8). The RDF indicated that both Ca-Si and Ca-O were within bonding distance, less than 3 Å; however, Ca-Si has a more prominent and higher peak, implying that the probability of occurrence for the Ca-Si bond is higher than the Ca-O bonds.
ions. In contrast, the second system was similar to the first system but included explicit water molecules. The MD trajectory analysis confirmed earlier observations from the DFT studies that Ca 2+ would preferentially bind to the silica region instead of the -CONH functional groups (Figure 7b). Ca forms Ca-O and Ca-Si bonds, as observed earlier from DFT. The radial distribution function (RDF), which is the probability of finding an atom at a distance r from another tagged atom, was analyzed between Ca-O and Ca-Si atoms (Figure 8). The RDF indicated that both Ca-Si and Ca-O were within bonding distance, less than 3 Å; however, Ca-Si has a more prominent and higher peak, implying that the probability of occurrence for the Ca-Si bond is higher than the Ca-O bonds.   Table 2 ions were put in closer proximity to PAM-Silica to enable easier binding since AIMD is computationally expensive and can only be observed for a few picoseconds. The initial stage depicts two water molecules and one Ca 2+ near PAM-Silica ( Figure 8a); however, Ca 2+ first bonds to the silica region, not the water molecules (Figure 9b). There is competition between the oxygen molecules in the water and the oxygen molecules in the silica region as the simulation evolves; however, there is no longer a Ca-Si bond ( Figure  9c). RDF was analyzed to clarify the Ca-O and Ca-Si interactions ( Figure 10). The Ca-O peak is much larger and higher than the Ca-Si peak, unlike the first system studied. The Ca-Si peak occurs primarily after 3 Å, implying that the two atoms are not within bonding distance most of the time; therefore, Ca 2+ would bond more to water molecules in the pres-  Table 2 ions were put in closer proximity to PAM-Silica to enable easier binding since AIMD is computationally expensive and can only be observed for a few picoseconds. The initial stage depicts two water molecules and one Ca 2+ near PAM-Silica ( Figure 8a); however, Ca 2+ first bonds to the silica region, not the water molecules (Figure 9b). There is competition between the oxygen molecules in the water and the oxygen molecules in the silica region as the simulation evolves; however, there is no longer a Ca-Si bond (Figure 9c).
RDF was analyzed to clarify the Ca-O and Ca-Si interactions ( Figure 10). The Ca-O peak is much larger and higher than the Ca-Si peak, unlike the first system studied. The Ca-Si peak occurs primarily after 3 Å, implying that the two atoms are not within bonding distance most of the time; therefore, Ca 2+ would bond more to water molecules in the presence of water molecules. Ca 2+ will still bind to the oxygen atoms in the silica region of PAM-Silica. Table 2 ions were put in closer proximity to PAM-Silica to enable easier binding since AIMD is computationally expensive and can only be observed for a few picoseconds. The initial stage depicts two water molecules and one Ca 2+ near PAM-Silica ( Figure 8a); however, Ca 2+ first bonds to the silica region, not the water molecules (Figure 9b). There is competition between the oxygen molecules in the water and the oxygen molecules in the silica region as the simulation evolves; however, there is no longer a Ca-Si bond ( Figure  9c). RDF was analyzed to clarify the Ca-O and Ca-Si interactions ( Figure 10). The Ca-O peak is much larger and higher than the Ca-Si peak, unlike the first system studied. The Ca-Si peak occurs primarily after 3 Å, implying that the two atoms are not within bonding distance most of the time; therefore, Ca 2+ would bond more to water molecules in the presence of water molecules. Ca 2+ will still bind to the oxygen atoms in the silica region of PAM-Silica.  Summarily, Table 2 shows a list of all model systems calculated to explain how the calculations begin from DFT optimization and progress to AIMD. DFT optimization was initially done for all the oligomers followed by binding energy calculations. The workflow confirms how computational chemistry techniques such as DFT and AIMD can be employed in studying and validating a system. In this case, different oligomers were investigated and only those with PAM had good binding affinities with Ca 2+ . QTAIM analysis confirmed the binding interactions using the BCPs. Ca-PAM silica-II which had the highest binding energy with Ca 2+ , was further analyzed using AIMD to confirm further the way by which Ca 2+ binds to the silicon atom in PAM-Silica-II. The findings of this work are in agreement with the recent work Sun et al. [58] though in their case the binding of Ca 2+ with PAM was done with anionic PAM. However, their studies confirm the strong interaction of Ca with Oxygen atoms as observed in this work.

Conclusions
The DFT studies indicate that Ca 2+ does not bind readily to mixed oligomers of acryl acid and acrylamide; however, it has a good binding affinity with PAM and silica functionalized PAM, with the latter having a very high binding affinity. The highest binding affinity occurs in the silica region and not the -CONH functional groups. AIMD calculations corroborate the DFT studies since it was observed from the MD trajectory that Ca 2+ binds to PAM-Silica by forming bonds with silicon. Ca 2+ initially forms a bond with silicon in the presence of water molecules, but this does not last long; therefore, it subsequently bonds with the oxygen atoms present in the water molecule. PAM-Silica is a suitable calcium scale inhibitor since it has a high binding affinity with Ca 2+ ; however, it will compete to capture the Ca 2+ in the presence of water molecules. Theoretical studies (DFT and AIMD) have provided atomic insights on how PAM_Silica may be used as an efficient scale inhibitor; however, it may be further modified considering the competition from water molecules when used as a scale removal in an aqueous environment. Nevertheless, future studies may explore the use of the ReaxFF molecular dynamics technique which combines the accuracy of AIMD whilst enabling the analysis of larger number of atoms at longer time scales. Hence, building a more realistic model particularly with respect to the oligomer. However, the challenge is to ensure an existing Force Field that contains all the atoms to be studied exist in the near future this would be possible as ReaxFF continues to develop.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/polym14122333/s1, Table S1. The DFT calculated energy values and binding energies of the studied complexes.