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[triple bond, length as m-dash]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. This journal is © The Royal Society of Chemistry.

Introduction

The construction of the carbon framework is one of the main goals of organic synthesis, and can be achieved using simple building blocks such as HCCH or, better, HCC− and even CC2−. To obtain the latter two, one needs to use an acetylide source. Many metal acetylides are known, including acetylides of practically all classes of metals,[1-10] and bi-metallic acetylides.[11-13]

Among metal acetylides, widely produced CaC2 now appears to be the most versatile choice for the synthesis of organic substances, including those that are biologically active,[14-24] as well as monomers.[25-27] Moreover, CaC2 is envisioned to become the feedstock for the sustainable, carbon-neutral chemical industry.[16,21,28] It is also considered valuable or promising in the synthesis of nanostructured materials,[29,30] agriculture,[31-34] and metallurgy (alloy making, see Section 2.3.8 in ref. 28). However, CaC2 is insoluble in organic solvents, which hampers its reactivity in the liquid phase.[19,21,23,35,36]

Preformed acetylides and acetylide intermediates play a key role in organic synthesis. Copper-catalyzed azide–alkyne cycloaddition (CuAAC)[37,38] is a widely used reaction, in which the main intermediate is the unstable Cu acetylide having a Cu–CC–R moiety. The use of a preformed acetylide makes the reaction significantly more facile.[39] Other organometallic acetylides such as those of Au,[40] Bi,[39] and Pt[41] also undergo dipolar cycloaddition to azides; the corresponding AAC reaction products are potent precursors to a wide range of substituted heterocyclic compounds.[39]

Activation of HCCH or RCCH via acetylide formation is necessary in CuAAC[42-45] and other[46] reactions. It was proposed that “any s-acetylide that can effectively recruit a p-bound copper atom will undergo annulation with a compatible dipolar partner.”[43] Ca acetylides undergo dipolar cycloadditions as well.[14,47-49]

Acetylide species, like HCC–Ca–OH, are often assumed to be intermediates in solution-phase organic reactions with pan class="Chemical">CaC2 that is insoluble by itself;[19,23,50] however, it is hard to detect these species in the liquid phase. Detection of soluble alkaline acetylides was reported under extremely basic conditions.[35,51] Ca acetylide was experimentally detected with Fourier transform infrared spectroscopy in solid CaC2 in a KBr matrix when subjected to trace amounts of H2O.[52] Acetylene chemistry in dimethyl sulfoxide (DMSO) under basic and super-basic conditions is a valuable and indispensable tool of modern organic chemistry.[46,53-56] Greater potential of practically valuable synthesis with CaC2 can be realized through understanding the unique performance of DMSO solutions.

Quantum chemical modeling of Ca acetylides in pan class="Chemical">DMSO, reported below, required innovative consideration of ionic pairs in solution that have strong solute–solvent interactions. To obtain consistent models, we combined conformational sampling by molecular dynamics (MD) with the density-functional tight-binding (DFTB) Hamiltonian followed by DFT post-treatment of the conformations and free energy computations. Conformational sampling with fast semi-empirical methods has seen tremendous development recently.[57-60] Combining them with DFT post-treatment allows, e.g., estimation of realistic IR spectra in solution[61] and reliable exhaustive conformational sampling of organic macrocycles.[62]

Given the importance of CaC2 as a sustainable pan class="Chemical">carbon source for organic synthesis and Ca acetylides as potent intermediates, we performed this computational study. Obtaining active acetylide intermediates is key to new solution-phase organic reactions with solid CaC2. As the main result, we propose a strategy for the development of new sustainable solution-phase transformations based on the utilization of CaC2.

Results and discussion

Thermodynamic model

The suggested strategy for modeling the dissolution of ionic solids, possibly including partial solvolysis and (or) solvent coordination, combines DFTB molecular dynamics and static DFT computations. Fig. 1 depicts a schematic description of the methodology. All parameters are listed in the Computational details section below and in the ESI.† By using the suggested methodology, it is possible to obtain dissolution free energies (following all three stages in Fig. 1), as well as to model chemical transformations of ionic pairs in polar coordinating solutions (performing computations in the last two stages).

Fig. 1

The core methodology and software used to perform the calculations outlined in Fig. 2.

According to previous studies, Ca acetylide can undergo organic transformations in pan class="Chemical">DMSO[19,20,23,48,50] and dimethylformamide[63] solutions upon the addition of water. That is why the consideration of partial hydrolysis is essential. We compare direct CaC2 solvation and solubilization involving partial hydrolysis by considering the elementary steps depicted in Fig. 2: the consideration starts from solid CaC2 and proceeds in the clockwise direction to solvated species. Note that the states of the intermediates (solid, gas, solvated) are explicitly defined in Fig. 2. We use an analog of the Born–Haber cycle, and model the solvation and hydrolysis as the sequence of hypothetical sublimation (ΔGsub), reactions in the gas phase (ΔGbind, ΔGprot, ΔGprotbind), and the subsequent solvation (ΔGsolv, ΔGprotsolv).

Fig. 2

Thermodynamic model: direct CaC2 solvation in DMSO vs. partial hydrolysis in the water/DMSO system (hydrolytic solubilization). The heterogeneity of CaC2 was accounted for by considering two polymorphs with fractions equal to their Boltzmann weights (Q).

We calculated Boltzmann weights for the stable CaC2 polymorphs, isomers and conformers:where ΔG is the relative Gibbs energy of the i-th isomer, conformer, or polymorph. Q values obtained in this way were used to compute the average free energies of species in solution.

Solid Ca carbide and its sublimation

CaC2 is reported to be a mixture of three polymorphs pan class="Chemical">CaC2-I, CaC2-II, and CaC2-III under ambient conditions.[64-66] Because of the inconsistencies in previous studies,[64,65] we present a more detailed discussion of our findings on the relative stability of CaC2 polymorphs in Section S2.†

We compared the stability (relative ΔG) of these three phases at 200, 300, and 400 K and found that the well-known tetragonal CaC2-I form is the most stable. The equilibrium distribution of stable CaC2 phases under standard conditions was estimated by computing Boltzmann weights according to the calculated ΔG values (see Fig. 2 and the ESI† for details).

According to the harmonic vibrational mode analysis at the PBE0-D3(BJ)/pob-TZVP-gCP level (see Section S2† for details), CaC2-III has an imaginary frequency at Γ point, so we excluded it from the set of allowed thermodynamic states for the sake of model consistency. Excluding CaC2-III from the calculation of ΔGsub in Fig. 2 resulted in a negligible correction of less than 0.1 kcal mol−1 due to its relatively high free energy. In contrast, we did not observe any imaginary modes at the chosen level in the cases of CaC2-I and CaC2-II. The relative and absolute stability of CaC2 polymorphs remains unclear under theoretical considerations with computational methods (see the discussion of the relevant literature in Section S2.1†).

It was hypothesized that anharmonic effects may affect the stability of pan class="Chemical">CaC2 phases.[65] We believe that further investigation of potential energy surfaces of CaC2 polymorphs may be worthwhile, ideally, with Born–Oppenheimer MD (BOMD), to elucidate possible anharmonicity of atomic vibrations. As long as the proposed methodology (Fig. 1) is modular, any refinements of ΔG values can easily be incorporated.

The first elementary reaction to consider is the sublimation of CaC2 (ΔGsub). Computed at the pan class="Chemical">PBE0-D3(BJ)/pob-TZVP level, the free energy of sublimation only slightly varied for the stable polymorphs: from 184.5 (CaC2-II) to 185.8 kcal mol−1 (CaC2-I). After Boltzmann averaging over stable CaC2 polymorphs we obtained 185.6 kcal mol−1 for the two-phase acetylide.

Ionic pairs in realistic solvent

DMSO solvates cations very efficiently, even more strongly than water.[67,68] Exergonicity of cation solvation in DMSO can be attributed to the formation of strong cation–oxygen bonds in the coordination shell. Using implicit solvent models and neglecting direct Ca–O-bonding when modeling the solvation of CaC2 and HCC–Ca–OH ionic pairs in DMSO leads to inconsistent results, as demonstrated in Section S3.†

We used Born–Oppenheimer molecular dynamics with the dispersion-corrected DFTB3-D3(BJ) Hamiltonian to determine the solvation shell of ionic pairs HCC–Ca–OH and [pan class="Chemical">Ca2+][C22−] in DMSO. First, we performed 10 ps-long isobaric-isothermal MD runs with the Berendsen thermostat and barostat to equilibrate the systems. Plots depicting the relaxation of the thermodynamic parameters V, P, and T, as well as of the sum U + PV + TSelec, are given in Section S5.†

Time evolution of the radial distribution functions (RDFs) demonstrates the equilibration of Ca2+ coordination number (CN, see Fig. 3 and Section S5†). Four DMSO molecules rapidly coordinate Ca2+ in the system with [Ca2+][C22−]. In the system with HCC–Ca–OH, in contrast, the fourth DMSO molecule bonded to Ca2+ only in the last picosecond of the equilibration run (Fig. 3d).

Fig. 3

Radial distribution functions in the [Ca2+][C22−]/DMSO (a–c) and HCC–Ca–OH/DMSO (d–f) systems. Time evolution of the interatomic distances is shown in (a), (b), (d) and (e); (a) and (d) depict the evolution in the equilibration runs, (b) corresponds to the simulated annealing run, and (d) corresponds to the sampling NPT run. Equilibrated structures of the Ca2+ solvation sphere are shown in (c) and (f). In (c), the RDFs in the last 5 ps of the run (b) are shown (cooled to 300 K after the annealing); in (f), the RDFs in the whole run (e) are depicted. Integral functions I(r) show the time-averaged CN of Ca2+. For clarity, the structural formula is shown only for one solvent molecule. Other solvent molecules are abbreviated as DMSO; they are also coordinated via O atoms.

Next, we subjected the HCC–Ca–OH system to another 10 ps NPT run, now using the Nosé–Hoover chain thermostat and Berendsen barostat, to sample the configuration space (Fig. 3e). The model of pan class="Chemical">CaC2 in DMSO was subjected to simulated annealing (NVT ensemble, Nosé–Hoover chain thermostat) by gradually heating the system to 600 K for 3 ps, preserving the temperature for 5 ps, gradually cooling the system for 3 ps, and then keeping the temperature at 300 K for another 5 ps.

In the sampling run, we observed no additional binding of pan class="Chemical">DMSO molecules in both systems (Fig. 3e); analogously, no additional DMSO molecules were bound as a result of the annealing (Fig. 3b). The resulting CNs are obtained from the integrals of the RDFs (Fig. 3c and f). Evidently, Ca2+ is six-coordinated in the HCC–Ca–OH system, which agrees with the experimentally observed CN of six for this cation in DMSO solutions.[69] One may consider C22− as a κ2- or, equally, η2-ligand. In dynamics at 300 K, however, C22− mostly resides in the singly coordinated mode, which is why the second peak is present on the corresponding RDF at ∼340 pm (Fig. 3c and Section S6†). Since such behavior of C22− was unexpected, we performed simulated annealing of [Ca2+][C22−] in DMSO instead of an NPT run to check if the solvent shell would equilibrate to the same CN after the annealing and no more DMSO molecules would bind to Ca2+. We suppose that C22− strongly electrostatically repels O-centers in DMSO, so only 4 DMSO molecules could bind to Ca2+ under the selected computational protocol.

We performed Boltzmann averaging over the ensembles of solvated [Ca2+][pan class="Chemical">C22−] and HCC–Ca–OH to obtain a conformationally sampled structure of Ca2+ solvation shell. For each system, we took 5 snapshots at distant trajectory points and cut Ca2+ with its first solvation shell representing a new model system for step 3 in Fig. 1, right (see Section S1.6† for details). Also, for both systems, we manually constructed conformations of the solvation shell by symmetrically placing 4 DMSO molecules in the equatorial plane of [(DMSO)4Ca(CCH)(OH)] and in the base of the tetragonal pyramidal [(DMSO)4Ca(CC)]. The latter artificial conformations were included as a stress test of the presented methodology. As shown below, these artificial conformations are negligible contributors to the pool of conformers. Geometries of all snapshot conformations obtained in this way were optimized at the PBE0-D3(BJ)/pob-TZVP-gCP level.

Using the gas-phase optimized geometries, we calculated ΔGsolv for every conformer structure within the pan class="Disease">SMD approach (Solvation Model based on Density). We listed relative ΔG of the conformers and the corresponding Q values in the ESI .xlsx table.† The most populated states (those with the highest Q) in DMSO and vacuum mostly do not coincide; in all cases except iso1 of [(DMSO)4Ca(CCH) (OH)] (shown in Fig. 4), the highest Q-conformers in DMSO are minor in a vacuum. Such a discrepancy can be expected because polar DMSO stabilizes polar conformations of the solute.

Fig. 4

Optimized structures of conformers: [(DMSO)4Ca(CC)] (top) and [(DMSO)4Ca(CCH)(OH)] (bottom). Relative Gibbs energies and Boltzmann weights at 300 K are given below the structures. The most abundant conformers iso1 and iso3 are depicted with marked close noncovalent C(sp)–H and O–H contacts. Note that the sum of the van der Waals radii for the C(sp)–H and O–H contacts is 2.88 and 2.62 Å, according to Bondi.[70]

The computed ΔG of the conformers in DMSO, the corresponding Boltzmann weights, and the optimized structures are given in Fig. 4. The relative free energies of conformers vary within 8.1 kcal mol−1 for [(DMSO)4Ca(CC)], and 11 kcal mol−1 for [(DMSO)4Ca(CCH)(OH)].

The DMSO molecules of the solvation shells can form hydrogen bonds with C22−, HCC−, and OH− ligands, thereby giving this considerable spread in relative ΔG in solution with the selected model systems and at the chosen level of theory. Close C–H and O–H contacts, as well as the reference sum of the van der Waals radii, are given in Fig. 4.

In contrast to the case of [(DMSO)4Ca(CC)] that is predominantly represented by iso1, the model conformer space of [(DMSO)4Ca(CCH)(OH)] has two significant structures iso1 and iso3, and the somewhat minor iso4. All this emphasizes the importance of conformational sampling for cluster-continuum modeling of species in solutions.

Boltzmann averaging over the conformers negligibly shifts the Gibbs energy of the ensemble of [(DMSO)4Ca(CC)] by 0.02 kcal mol−1, relative to the lowest energy conformer iso1. Similarly, the ensemble-averaged Gibbs energy of [(DMSO)4Ca(CCH)(OH)] is 0.19 kcal mol−1 higher than that of iso3. Even though the averaging correction at 300 K is minor, we still suggest using the presented two-step conformational sampling (BOMD plus static DFT). Therefore, in the absence of the sampling, if one considers only a minor conformer with low Q, ΔG of elementary reaction steps can be inaccurate by several kcal mol−1.

Now we can estimate ΔG of the following reactions using the averaged free energies of the solvated species:

The reactions in eqn (1) and (2) are among the model steps in Fig. 2. We attribute the extremely exergonic effect of reactions (1) and (2) to the formation of strong Ca–O bonds, and—equally importantly—to the formation of many hydrogen bonds in the solvation shell. Even anionic centers of HCC− and CC2− ligands are hydrogen bond acceptors, as can be seen from the abundance of close contacts in the structures in Fig. 4.

The last step in Fig. 2 is to compute solvation energies of [(DMSO)4Ca(CC)] and [(DMSO)4Ca(CCH)(OH)] (see eqn (3) and (4) below). The process of the immersion of electro-neutral species [(DMSO)4Ca(CC)] and [(DMSO)4Ca(CCH)(OH)] into DMSO is moderately exergonic, in contrast to the gas-phase formation of the coordination shell, as in eqn (1) and (2).

In this work, we selected M06-2X/6-31+G** as the underlying method for SMD computations of ΔGsolv since 6-31+G**[71] was included in the original SMD parameterization,[72] and since this CSM is often used in conjunction[73-75] with the M06-2X functional.[76] In Section S4,† we demonstrate that predictions of ΔGsolv with SMD at the M06-2X/6-31+G** level deviate by only 0.7 kcal mol−1 from those obtained at the M05-2X/6-31+G** level that was used in the original parameterization of SMD.[72]

A closely related two-step model process is the hydration of Ca2+ (Table 1). The details of the performed modeling of pan class="Chemical">Ca2+ solvation in water are described in the ESI table.† As in the previous case with DMSO, most of the solvation exergonicity stems from the formation of the coordination sphere. The experimental value for the hydration of Ca2+ in water varies from −359.7 (ref. 77) to −386.2 (ref. 78) kcal mol−1 (the divergence is equal to 26.5 kcal mol−1), so the comparison with the experiment is possible, but cannot be performed reliably. Depending on the experimental reference, our computational estimation of ΔGsolv deviates from −1.3 to −27.8 kcal mol−1. The continuum models, used directly, i.e., without the explicit inclusion of a solvation shell, yield minimal deviations of +66.0, +64.7, and +75.9 kcal mol−1 for the COnductor-like Screening MOdel (COSMO), conductor-like polarizable continuum model (C-PCM), and SMD, respectively. Cluster-continuum computations, with our two-step calculation of ΔGsolv in H2O being one of this class, are a well-established approach to the modeling of ionic species in solution.[79-82]

Hydration of Ca2+a

Transformation	ΔGrxn, kcal mol−1
Ca(g.)2+ + 7H2O(g.) ⇌ [Ca(H2O)7](g.)2+	−205.6
[Ca(H2O)7](g.)2+ ⇌ [Ca(H2O)7](aq.)2+	−181.9
Ca(g.)2+ + 7H2O(g.) ⇌ [Ca(H2O)7](aq.)2+	−387.5
Experimental reference	−386.2 (ref. 78) to −359.7 (ref. 77)
Classical (non-quantum) electrostatic models	−377.3,[83] −403.2 (ref. 78)

The binding of H2O to Ca2+ was modeled at the RIJK-PBE0-D3(BJ)/def2-TZVP-gCP level; the hydration was modeled using SMD (M06-2X/6-31+G**).

Understanding CaC2 solubilization: direct solvation vs. hydrolytic solubilization

Table 2 summarizes the cumulative thermodynamic effect of the CaC2 dissolution in pure pan class="Chemical">DMSO (+25.1 kcal mol−1) and its favorable hydrolysis in the DMSO/water solvent system (−9.1 kcal mol−1). The former is markedly endergonic, in accordance with the experimental observations of CaC2 inactivity in pure DMSO.[19,21,23,35,36] Therefore, CaC2 dissolution in DMSO is thermodynamically forbidden. The solubilization can be achieved via the steady protonation of CaC2 at the solid–liquid interface, and the concomitantly formed HCC− can participate in subsequent transformations.

Summary: the unfavorable CaC2 solvation vs. protolysis-assisted solubilization of CaC2a

Transformation	ΔGrxn, kcal mol−1
Direct solvation
CaC2(s.) ⇌ [Ca2+][C22−](g.)	ΔGsub = 185.6
[Ca2+][C22−](g.) + DMSO(g.) ⇌ [(DMSO)4Ca(CC)](g.)	ΔGbind = −138.0
[(DMSO)4Ca(CC)](g.) ⇌ [(DMSO)4Ca(CC)](solv.)	ΔGsolv = −22.5
CaC2(s.) ⇌ CaC2(solv.) (same as [(DMSO)4Ca(CC)](solv.))	ΔGsub + ΔGbind + ΔGsolv = 25.1
Protolysis-assisted solubilization
[Ca2+][C22−](g.) ⇌ HCC–Ca–OH(g.)	ΔGprot = −104.5
HCC–Ca–OH(g.) + 4 DMSO(g.) ⇌ [(DMSO)4Ca(CCH)(OH)](g.)	ΔGprotbind = −72.9
[(DMSO)4Ca(CCH) (OH)](g.) ⇌ [(DMSO)4Ca(CCH)(OH)](solv.)	ΔGprotsolv = −17.4
CaC2(s.) ⇌ HCC–Ca–OH(solv.) (same as [(DMSO)4Ca(CCH) (OH)](solv.))	ΔGsub + ΔGprot + ΔGbindH+ + ΔGsolvH+ = −9.1

Gas-phase thermochemistry was modeled at the PBE0-D3(BJ)/pob-TZVP-gCP level; bulk solvent effects were modeled using SMD (M06-2X/6-31+G**).

H2O can easily protonate C22−, as it is a much stronger acid. At the same time, H2O is less acidic than HCCH in DMSO/water solutions, as seen from Table 3. We used a coarse quantum chemical approach to calculate free energies of H2O, DMSO, HCCH, and PhC≡CH deprotonation in DMSO (see also the ESI .xlsx table†). The reference pKa values show that H2O is ∼102 times less acidic than HCCH, and ∼1010 times less acidic according to our calculations. The deprotonation of HCC− yielding C22− should be as unfavorable as DMSO autoprotolysis.

Acidity in DMSO

Transformation	ΔGrxn, kcal mol−1	Calculated pKa	Reference pKa	Deviationa
HCCH(solv.) ⇌ HCC(solv.)− + H(solv.)+b	34.2	25.1	29.7[84]	−4.6
HCC(solv.)− ⇌ −CC(solv.)− + H(solv.)+	50.6	37.1	—
PhC≡CH(solv.) ⇌ PhC≡C(solv.)− + H(solv.)+	34.6	25.4	28.7[85]	−3.6
H2O(solv.) ⇌ HO(solv.)− + H(solv.)+	48.8	35.8	31.4[86]	4.4
CH3S(O)CH3(solv.) ⇌ CH3S(O)CH2(solv.)− + H(solv.)+	50.5	37.0	35.1[86]	1.9

Between calculated and reference values.

The Gibbs free energy of a proton in DMSO is taken as the sum of G in the gas phase at 298.15 K and 1 atm (ref. 87) and ΔGsolv of H+ in DMSO.[88]

Acidities (pKa) of DMSO and HCC− (second stage) are nearly equal, according to the PBE0/ma-def2-TZVP + SMD calculation. Therefore, we may suppose DMSO as a possible protolytic agent for C22− in solution. Indeed, C22− anions can undergo rapid protonation by DMSO (see the ESI .xlsx table†), meaning that the formation of free acetylide dianions in such a solution system is hardly possible. The solvent is not aprotic enough, even if we find a way to effectively solvate C22− with anion-sequestering host molecules, e.g., cavitands.

Modeled at the PBE0/ma-def2-TZVP level; bulk solvent effects were accounted for by applying SMD (M06-2X/6-31+G**) (see Section S1.3† for details).

We also estimated the favorability of HCC− protonation by the DMSO molecules of the Ca2+ solvation shell (as in Scheme 1). The free energies of activation of the two evaluated pathways are 20.7 and 21.5 kcal mol−1. Moreover, the process is endergonic by 17.7–20.3 kcal mol−1. Thus, the protonation of the acetylide in [(DMSO)4Ca(CCH)(OH)] is somewhat kinetically unfavorable, also being clearly unfavorable thermodynamically.

Scheme 1

The unfavorable process of the HCC− protonation by a DMSO molecule from the first solvation shell. All DMSO molecules are coordinated via O atoms.

Other protic molecules such as inorganic acids HX, pan class="Chemical">HClO4, and CF3SO3H (X = Cl, Br, I) are an inappropriate choice for the protolytic activation of CaC2. These acids are reported to be strong in DMSO.[89] That is why their DMSO solutions can protonate not only CC2− but also HCC−, thereby decomposing the reactive acetylide intermediate. Moderate acidity is crucial in our case.

Summary, conclusions, and outlook

A protolytic agent plays a crucial role in the activation of a carbide in solution reactions. pan class="Chemical">Water is a unique agent since it is less acidic than HCCH in DMSO. The formation of anionic monoacetylide intermediates is the way to activate CaC2 in liquid-phase organic transformations. The dynamics of CaC2(s.) protolysis, proceeding at the solid–liquid interface, may thus be of paramount importance for further understanding of CaC2 activation. Research on the dynamics of this interface process is currently underway in our group.

We tested a new modeling strategy for solvated ionic pairs formed in the process of dissolution or protolysis of ionic crystals. It allowed us to obtain ensemble-averaged ΔG of reactions in solutions with species for which no CSM parameterization is available. We plan to further use and test this computational methodology, as well as encourage its use in other groups.

The methodology is modular, as it consists of three distinct steps depicted in Fig. 1. Therefore, evaluating alternative tight-binding parameterizations (e.g., eXtended Tight-pan class="Chemical">Binding methods and GFNn-xTB)[90,91] and CSMs (such as COSMO-RS)[92] is advised to determine an optimal level of theory. Accounting for anharmonic effects in calculations of free energies can be another option for incremental improvement of the methodology. Such effects can be incorporated in the solid state step (Fig. 1, left),[93,94] as well as in the MD[95,96] and molecular DFT steps[97,98] (Fig. 1, middle and right, respectively), although in the latter two cases this may be technically non-trivial. It was shown in recent studies of solid state and surface systems that anharmonic effects may be crucial.[99,100] However, we should also mention a critique of existing approximations for computation of anharmonic free energies.[101]

Strongly coordinating solvents such as DMSO form a well-defined solvation shell that should be sampled with BOMD. A very economical choice is to use a tight-pan class="Chemical">binding Hamiltonian such as DFTB3 with empirical corrections for non-covalent interactions. In our case, running even relatively short equilibration trajectories of 10 ps yielded ensemble-sampled structures of solvation shells. Free energy computations with MD methods require rather elaborate techniques.[102,103] That is why Boltzmann averaging over an MD-obtained set of solvation shell conformers can be a convenient option. The proposed combination of semi-empirical BOMD and static DFT computations of ΔG values is cost-efficient since the most demanding step—the sampling of conformer space with MD—is feasible even on a personal workstation. We performed most of the MD simulations on an entry-level graphics processing unit (GPU) and a gaming central processor (CPU, see the ESI† for details).

As a fundamental result, we propose a strategy for CaC2 activation in organic media that can boost further development of green and sustainable synthetic methodologies based on the use of pan class="Chemical">calcium carbide. DMSO, as well as dimethylformamide which is widely used in reactions with CaC2, is not a particularly “green” solvent. Less toxic polar aprotic solvents that allow water pKa higher than acetylene pKa would be a better choice for future organic synthesis; no less important is the propensity to effectively solvate Ca2+ by forming strong Ca-solvent bonds, such as, e.g., Ca–O. There are few such solvents. Here we assessed H2O as a suitable green protolytic agent for a solid acetylide. However, we hypothesize that any molecule less protic then HCCH in a given solvent can play its role, thereby allowing new synthetic transformations. Computational methods, as described in this work, can help in the evaluation of known green solvents for sustainable organic synthesis with CaC2 or in the search for new ones, as well as in the discovery of new protolytic agents for the activation of CaC2.

Computational details

Solvation free energies of species in Fig. 2 were estimated using ORCA 4.1.2.[104] The solvation model based on density (SMD) was selected for this purpose.[72] All implicit solvent computations were performed on gas-phase geometries, as in the original studies.[72,105,106] ΔGsolv is the difference between the single point (total) energies of gas-phase geometries with pan class="Disease">SMD applied, and without. We chose the diffuse basis set (BS) 6-31+G**,[71,107,108] since we modeled anionic species in pKa estimations, and, at the same time, this basis set was included in the SMD parameterization. For Ca, we accepted the default ORCA[104] choice and used diffuse exponents from 6-311+G**, which was adopted from the EMSL basis set exchange.[109-111] In the original work, SMD was parameterized for use at the M05-2X/6-31+G** level of theory;[112] the corresponding functional is, however, unavailable in ORCA. Therefore, we selected its successor, M06-2X,[76] that is successfully employed in computations with SMD[73-75] (see also Sections S1.2 and S4†).

The CRYSTAL17 (ref. 113) program was used for evaluation of gas-phase energies and thermodynamic corrections for reactions in Fig. 2. The pob-TZVP basis set was used.[114] The pan class="Chemical">PBE0 functional was selected. Empirical corrections for dispersion interactions (D3, including the Becke–Johnson dumping function) and geometrical counterpoise corrections (gCP)[115,116] were included (see Section S1.1† for details).

The self-consistent charge density-functional tight-binding method pan class="Chemical">DFTB3 (ref. 117 and 118) was used for Born–Oppenheimer molecular dynamics of model DMSO solutions. The computations were performed in the DFTB+ program (ver. 19.1).[119] The Third-Order Parametrization for Organic and Biological Systems (3OB) of SCC-DFTB was used.[120-122] All parameters selected in SCC-DFTB3 computations are given in Section S1.6,† together with a description of how model systems with explicit DMSO solvent were constructed.

We modeled CC2− and –CCH protonation by DMSO using the B97-3c method[123] for gas-phase calculations and SMD for the evaluation of solvation free energies (as described above). These computations were performed with ORCA 4.1.2 (see Section S1.5† for details).

Travis (update Jan 01, 2019)[124] was used to plot radial distribution functions.

Conflicts of interest

There are no conflicts to declare.