Understanding the solubilization of Ca acetylide with a new computational model for ionic pairs

The unique reactivity of the acetylenic unit in DMSO gives rise to ubiquitous synthetic methods. We theoretically consider CaC2 solubility and protolysis in DMSO and formulate a strategy for CaC2 activation in solution-phase chemical transformations. For this, we use a new strategy for the modeling of ionic compounds in strongly coordinating solvents combining Born–Oppenheimer molecular dynamics with the DFTB3-D3(BJ) Hamiltonian and static DFT computations at the PBE0-D3(BJ)/pob-TZVP-gCP level. We modeled the thermodynamics of CaC2 protolysis under ambient conditions, taking into account its known heterogeneity and considering three polymorphs of CaC2. We give a theoretical basis for the existence of the elusive intermediate HC 
<svg xmlns="http://www.w3.org/2000/svg" version="1.0" width="23.636364pt" height="16.000000pt" viewBox="0 0 23.636364 16.000000" preserveAspectRatio="xMidYMid meet"><metadata>
Created by potrace 1.16, written by Peter Selinger 2001-2019
</metadata><g transform="translate(1.000000,15.000000) scale(0.015909,-0.015909)" fill="currentColor" stroke="none"><path d="M80 600 l0 -40 600 0 600 0 0 40 0 40 -600 0 -600 0 0 -40z M80 440 l0 -40 600 0 600 0 0 40 0 40 -600 0 -600 0 0 -40z M80 280 l0 -40 600 0 600 0 0 40 0 40 -600 0 -600 0 0 -40z"/></g></svg>
 C–Ca–OH and show that CaC2 insolubility in DMSO is of thermodynamic nature. We confirm the unique role of water and specific properties of DMSO in CaC2 activation and explain how the activation is realized. The proposed strategy for the utilization of CaC2 in sustainable organic synthesis is outlined.


Computational Details
A notice should be given on calculations ∆G of reactions involving charged species in the gas phase. In this case, one needs to calculate molar Gibbs energies of naked cations or anions, which is only a hypothetical state of matter. It was suggested to use the formulae for the ideal gas as, for example, the Sackur-Tetrode equation. 1,2 Indeed, ∆G computations treating charged gas-phase species as the ideal gas are ubiquitous and offer well acceptable accuracy, compared to experiment. [3][4][5][6] Therefore, we used the IG approximation for gas-phase cations and anions in this work, while being inherently dissatisfied with such an approach.
Except for the computations described in Section S1.1, we run all calculations on a simple workstation PC with an Intel Core i7-9700K CPU overclocked to 5 GHz (gaming-grade CPU), 48 Gb of RAM, and an entry-level MSI GeForce GTX 1050 Ti GPU (latter was used for the computations in the DFTB+ program). Notably, two runs were performed most of the time simultaneously (on CPU and on GPU).

Solid State and Gas-phase CaC2
All solid-state density functional theory (DFT) calculations were performed using the CRYSTAL17 software. 7 In CRYSTAL17, crystalline orbitals (CO) are approximated as a linear combination of Bloch sums of Gaussian type orbitals (GTOs). We used a triple-zeta valence basis set with polarization quality (pob-TZVP) optimized for solid-state calculations 8 in combination with the geometrical counterpoise correction (gCP). 9, 10 We used the hybrid PBE0 functional 11 with Grimme's dispersion D3 correction, including three-body term and Becke-Johnson dumping function. 12,13 S4 Overlap thresholds for Coulomb and exchange integrals were tightened to 10 -8 , 10 -8 , 10 -8 , 10 -8 , and 10 -16 . Integration over the Brillouin zone was performed on the Monkhorst-Pack grid with the shrinking factor 6 corresponding to the primitive tetragonal cell with inversely scaling to the cell constants. We performed full optimization of cell parameters and atomic displacements. For the self-consistent calculations (SCF) the tolerance on total energy change was set to 10 -10 Hartree.
Gas-phase CaC2 calculations were performed on the same level of theory as for the solid-state (PBE0-D3/pob-TZVP-gCP).

∆Gsolv Computations with CSMs
In all computations with CSMs applied, ∆Gsolv is the difference between single point (total) energies of gas-phase geometries with a CSM applied, and without.

Ca 2+ Hydration
Since we were unable to model systems with Ca 2+ cations in the periodic CRYSTAL17 program, we had to use the ORCA program at the level of theory closely resembling that one used for the computations in CRYSTAL17. PBE0 functional was selected in conjunction with the def2-TZVP basis set. 14 The RI approximation was used for both Coulomb and exchange terms; 33,34 accordingly, the def2/JK basis set 35 was employed. The gCP empirical correction for the BSSE was used. 9 Also, we used D3-correction for dispersion interactions with the Becke-Johnson damping function, and three-body terms included. 12,13 The TightSCF, GRID6, and S6 NOFINALGRID options were set. Vibrational frequencies were computed numerically. Fock matrices were recalculated each 8-th KS-SCF iteration. See the "Ca2+ and HCC-Ca-OH" tab of the supporting .xslx table for the formulae used in the calculation of the Ca 2+ Gibbs free energy.

Mechanisms of C≡C 2− and − C≡CH Protonation by DMSO
This Subsection is relevant to the "Barriers" tab of the supporting .xslx table. For a rough estimation of free energy barriers of the protonation within the Born-Oppenheimer approximation, we selected the B97-3c method. 36 It is a combination of a re-parameterized B97 37 functional and empirical corrections. 12 In B97-3c, a specially optimized triple-ζ basis set is used. We performed the B97-3c computations within the RI approximation and used the default auxiliary basis set (def2-mTZVP/J).
To find the transition states of the proton transfer, first, relaxed potential surface scans were performed in which H + migrated from the donor to acceptor atoms. Second, the highest energy structures from the scans were used for geometry optimization to transition states. All transition state structures exhibited a single imaginary mode corresponding to the proton vibration along the donor-acceptor line. The structures are included in the supporting PDF file in the XYZ format.
Each .xyz file contains the value of the imaginary mode within it as a commentary line.
The bulk solvent effects we accounted for with SMD. The ∆ !"#$ values were computed for the stationary points and transition states, as described in Section S1.2. S7

DFTB Modeling of DMSO Solutions
The integration of the equations of motion was done with the velocity Verlet algorithm and timestep of 0.5 fs. The Berendsen thermostat and barostat 38 were used in the equilibration MD runs. The Nosé-Hoover chain thermostat (chain length equal to 3) 39 was used for sapling runs and in the simulated annealing (SA). Except in the SA run, the target temperature of the thermostat was set to 300 K, while the external target pressure was set to 10 5 Pa. No constraints on cell vectors and angles were applied during the NPT runs. Initial velocities in the equilibration runs were set according to the Boltzmann distribution at 300 K. In the sampling runs, velocities from the last iteration of the equilibration runs were used as initial.
During the equilibration runs, the coupling strength of the Berendsen thermostat, as well as that of the barostat, was controlled by setting the timescale parameter as equal to 100 fs. Simulated annealing was performed at the cell volume relaxed in the equilibration run (19.996 ✕ 19.996 ✕ 19.996 Å, α = β = γ = 90°). The temperature profile for the SA is described in the main text; the coupling strength of the Nosé-Hoover chain thermostat was set to 10 THz. In the sampling run, the coupling strength of the Nosé-Hoover chain thermostat was 1 THz, and the coupling strength of the barostat was 1000 fs.
The SCC convergence tolerance was equal to 10 -6 Hartree; modified Broyden's mixing scheme was used. 40 The parameter α in the Ewald electrostatic summation was determined automatically by the DFTB+. The tolerance for the Ewald summation was 10 -8 Hartree. The MAGMA library Hydrogen atom interactions were damped with the parameter ζ equal to 4.05. 46 Dispersion interactions were accounted for by the inclusion of the empirical corrections. 12,47 The models of Ca(C≡C) and HO-Ca-C≡CH in DMSO were constructed in the following way.
As a model structure of liquid DMSO, we took a box obtained via a classic MD simulation with the OPLS forcefield, 48 as provided at the "virtualchemistry.org" website. 49 We put the ionic pairs

Relative Stability of CaC2 Polymorphs
According to previous research, [53][54][55][56][57][58] there are three co-existing phases of calcium carbide at room temperature: CaC2-I (I4/mmm); CaC2-II (C2/c); CaC2-III (C2/m). The proportion of three polymorphic forms depends on the reaction conditions. 55 Tetragonal phase CaC2-I is prevalent at room temperature, but its thermodynamic stability remains uncertain as the previous theoretical works provide inconsistent results. CaC2-I is stable according to earlier reports. 57 54 The authors considered two monoclinic phases (CaC2-II, CaC2-III) to be stable, with CaC2-II being a low-temperature phase. In contrast, CaC2-III phase is metastable according to the experimental work by Knapp and Ruschewitz. 55 Doll, Jansen, et al. performed an extensive computational search for CaC2 polymorphs. 58 Several new structures were predicted. CaC2-I was found as particularly stable at standard pressure, as well as two of newly predicted polymorphs. Here we considered CaC2-I, CaC2-II, and CaC2-III since these three phases are often observed in experimental studies.
We performed phonon calculations for the three CaC2 modifications using appropriate supercells to consider multiple points of the Brillouin zone. Our computations indicate the dynamical stability of CaC2-I (no imaginary frequencies at the X point was found) and CaC2-II phases, while CaC2-III has an imaginary frequency at Г point. Relative electronic and Gibbs free energies of the phases are listed in Table S1. According to our calculations, CaC2-I is the most stable polymorph. CIF files enclosing the primitive cells are provided in the SI.

S10
We found no imaginary frequencies in the optimized supercell of CaС2-I. Briefly, we used the Gaussian basis set pob-TZVP 8 instead of the plane-wave BS, the hybrid PBE0 11 functional instead of its pure GGA counterpart, and included corrections for BSSE 9,10 and dispersion interactions. 12,13 All computational parameters are given below.
Authors of Ref. 54 Table S1). CaC2-I, as the most stable phase, is the most abundant in the equilibrium distribution. We performed calculations of a set of (CaC2)n clusters with the number of formula units (n) ranging from 1 to 8 (see figures in the supporting .xslx table). All obtained clusters are dynamically S11 stable (no imaginary frequencies were observed). XYZ files with the cluster structures are provided in the supporting PDF file.
The Gibbs free energy of the cluster formation from crystal are presented in Table S2. We averaged the Gibbs free energy of the solid CaC2 according to the Boltzmann distribution, as described above. Despite the energy of the cluster formation per formula unit decrease, the full energy increases significantly. Considering the fact that implicit solvation of molecules usually results in an exergonic effect of only 20-30 kcal/mol, we suggest crystal decomposition to the molecule CaC2 (n=1) is by far the most favorable process.

CaC2 Ionic Pair in Implicit (Continuum) Solvent
An overly simplified way is to model the solvation or hydrolysis of CaC2(s.) with a continuum solvation model (CSM). We calculated ∆Gsolv of the acetylide species HC≡C-Ca-OH in DMSO with three CSMs, namely, COSMO, 15 C-PCM, 61 and SMD 62 (see Table S3

Comparing SMD results at different levels of theory
We checked, if ∆Gsolv computed with SMD at the M06-2X/6-31+G** and PBE0/6-31+G** levels functionals are close to the results obtained at the M05-2X/6-31+G**, which was used in the original parameterization. ∆Gsolv were calculated as the difference in total energies with and without SMD applied; gas-phase geometries optimized at the PBE0-D3(BJ)/pob-TZVP-gCP level were used for single point energy evaluations. Evidently, using SMD at the M06-2X/6-31+G** level leads to a small deviation of 0.7 kcal/mol on a small test set of related species. On the contrary, using PBE0/6-31+G** ∆Gsolv estimations with SMD may lead to small but pronounced deviation of 1.7 kcal/mol, as shown in Table S4. a The 6-31+G** basis set was used, see computational details in Section S1.2; b Root-mean-square deviation, relative to the SMD at the original M05-2X/6-31+G** level; All values are in kcal/mol. Corresponding full electronic energies are given in the supporting .xslx table. S14

Details of the MD runs
In this Section, plots of thermodynamic parameters (G, T, P, V) and radial distribution functions are given. Note, that for consistency with the definitions of the thermodynamic parameters in DFTB+, the following equation for G was used: I.e., only electronic entropy was included.