Crossroads at the Origin of Prebiotic Chemical Complexity: Hydrogen Cyanide Product Diversification

Products of hydrogen cyanide (HCN) reactivity are suspected to play important roles in astrochemistry and, possibly, the origin of life. The composition, chemical structure, and mechanistic details for formation of products from HCN’s self-reactions have, however, proven elusive for decades. Here, we elucidate base-catalyzed reaction mechanisms for the formation of diaminomaleonitrile and polyimine in liquid HCN using ab initio molecular dynamics simulations. Both materials are proposed as key intermediates for driving further chemical evolution. The formation of these materials is predicted to proceed at similar rates, thereby offering an explanation of how HCN’s self-reactions can diversify quickly under kinetic control. Knowledge of these reaction routes provides a basis for rationalizing subsequent reactivity in astrochemical environments such as on Saturn’s moon Titan, in the subsurface of comets, in exoplanet atmospheres, and on the early Earth.


Convergence of Barrier to Polymerization into Polyimine, as a Function of Chain Length
. Gibbs free energy barriers (278 K) for the early stages of polyimine formation as a function of chain length. Grey and black circles correspond to PCM-PBE-D3/6-31G(d,p) and PCM-B3LYP-D3/6-31+G(d,p) levels of theory, respectively. The implicit PCM solvent model is that of water but with an adjusted dielectric constant of 144.8.

Construction of path collective variables for steered simulations.
Umbrella sampling simulations were performed along reaction paths defined by path collective variables. 9 The path collective variables s, and z are defined as: where N k refers to the number of reference structures and is a parameter determining the distance between the reference structures in s,z space.
is a distance function. For this study we made use of the coordination-based distance function by Pietrucci and Saitta, 10 ( ) = ∑ [ ( ) − ] 2 , where ( ) is the coordination between atom I and atoms of type S and is the corresponding coordination in reference structure k. The coordination patters are calculated as, where N and M are set to 8 and 16 following the recommendations of Pérez-Villa et al., 11 and ( ) corresponds to the Euclidian distance between atoms I and J. ′ 0 is a reference distance which was set to 1.8 Å for distances between non-hydrogen atoms and 1.5 Å for distances between hydrogen atoms and non-hydrogen atoms. In this work, we used two separate reference structures for each studied reaction, corresponding to reactants and products. The -parameter was adjusted so that = 2.3 ( ) . 11 We initially ran a reference simulation of 2 performed with a cutoff of 200 Ry and a Nose-Hoover thermostat of length 4 and a time constant of 500 fs. 12 In all subsequent reference states, we used the canonical sampling through velocity rescaling thermostat 13 for both equilibration and production runs (with a time constant of 50 fs for production runs and 1 fs for equilibrations) since it worked better with our steered simulations. We used a more accurate cutoff of 280 Ry for all other simulations besides those of 2. We nonetheless still utilized the reference state simulation of 2 for analysis. We do not expect the cutoff and thermostat to change the atomic coordination of the atoms in the reference simulations in a meaningful way. However, we note that convergence with a cutoff of 280 Ry is better (< 1 meV/atom) than with 200 Ry ( < 2 meV/atom). 280 Ry has also been shown to produce accurate forces and energies for NVT ensembles of water. 14 Tables S1-S5 shows the coordination patterns used to construct the s and z variables for all reaction steps in pathways 1-3. For each step, we defined coordination patterns of the atoms of the cyanide anion. We also included the coordination of the carbon subjected to the nucleophilic attack by the cyanide anion as well as that of the nitrogen bonded to it. The coordination in the reference states was calculated based on the average distance between atoms during 20 ps production runs preceded with a 5 ps equilibration.
In the reference state simulations of 2 and 3 there occurred a proton transfer to the cyanide anion. Analysis of the C-H coordination based on the whole production run of the cyanide anion gave of 0.69 and 0.73 for the simulation of 2 and 3 respectively. This coordination pattern was used during the transition state equilibration. For subsequent committor analysis and umbrella sampling simulations a coordination pattern based on the section of the production run prior to the proton transfer was used. Prior to the proton transfer was equal to 0.21 in the simulations of 2 and 3 both.

Force Constants in transition state equilibrations and umbrella sampling simulations.
The force constants were selected to not exceed 2 , where is the Boltzmann constant, T refers to the temperature (278 k) and is the standard deviation of the s coordinate in a typical simulation of the reactant or product state.
The input transition state (TS) guess structures were DFT optimized as minimal molecular models (later solvated in HCN prior to simulation). TS1-3 were optimized in gas phase with a 6-311++G(d,p) basis set using the B3LYP functional and Grimme D3 corrections. Transition states TS4 and TS5 were optimized with PBE-D3 and a 6-31G(d,p) basis set. In all cases the TSs were allowed to equilibrate for 20-25 ps. During TS equilibration, the s value of all TSs was restrained with a force constant of 10 000 kJ/mol. We added an additional weak bias of 3000 kJ/mol to the z variable during the equilibration of TS5. The z bias was introduced to hinder an otherwise observed immediate proton transfer to 16 at the start of the equilibration. A protonated 16 computes as 20 kcal/mol above the corresponding neutral molecule at the PCM-PBE-D3/6-31G(d,p) level of theory. Therefore, the spontaneous formation of a protonated 16 was likely due to the coarse initial solvent configuration prior to equilibration.
Input structures for the umbrella sampling windows were taken from trajectories reaching between the equilibrated TS to the reactant and product basins. The umbrella sampling windows were constructed so that they were at most 0.033 s units apart, corresponding to approximately 4 % of the distance between the reactants and products for all studied reaction steps. The s value in each umbrella sampling window was constrained using a harmonic potential on the form ( ( ) − 0 ) 2 . In Tables S6-S10 we tabulate the force constant, , and s value associated with each umbrella sampling window, 0 . Table S6. 0 and parameters used to describe the harmonic potentials along the reaction coordinate between 2,3-diiminopropanenitrile (16) and iminoacetonitrile (2). 0 denotes the center of the potential and is the force constant of the potential.  Figure S5. Histograms of sampled reaction coordinates in the umbrella sampling windows for all the reaction steps in pathways 1, 2 and 3. A good overlap between adjacent histogram is necessary for a reliable reconstruction of the relative Gibbs free energy. Figure S6. The free energy profiles for all reaction steps in pathways 1-3, obtained using umbrella sampling at the PBE-D3/DZVP-GTH level of theory. Error bars correspond to standard deviations that were obtained using block averaging as implemented in the weighted histogram analysis method (WHAM) program. 15

Statistical Error of Energy Profiles
The uncertainty of the relative energy between two points on a free energy curve depends on the uncertainty of the energy in both points. Therefore, we can use the standard deviations of the umbrella sampling curves (Fig. S7) to estimate the uncertainty of the relative energy of the transition state or product compared to the reactants. Under the assumption that the energies are independent and follow a Gaussian distribution, the error (Δ ( )) of the relative energy of products and transition states can be computed according to (Δ ( )) = √ ( ( )) + ( ( )), where i= 3 or 4 (corresponding to a trimerization (3) or tetramerization (4), respectively) and j refers to the transition state or product.
The relative energy error of the trimerization profile, Δ 3 , will influence the overall error of the relative energy of the tetramer compared to the dimer. To construct a complete energy profile Δ that connects 2 to the tetramers, the trimerization and tetramerization energy profiles ( Δ 3 and Δ 4 respectively) can be added according to Δ = Δ 3 + Δ 4 . The standard deviation of final relative energy of the states along the tetramerization energy profile are then computed as = √ (Δ 3 ( )) + (Δ 4 ( )) , where i is either a tetramer (15 or 17) or a second transition state of the pathway (either TS2, TS4, TS5). The final standard deviations of relative energies along the pathways lie between 1.1 and 2.9 kcal/mol (Fig. S8). Figure S7. The Gibbs free energy for formation of 2-amino-3-imino butanedinitrile (15) and 2,3,4triiminobutanenitrile (17) at 278 K as presented in Figure 3 of the main text, here shown together with the associated error bars. The errors correspond to standard deviation of the relative energy and were computed as described in the text above.

Cyanide Addition on Other Nitriles
The first dimerization step in Pathways 1-3 is rate-limiting and consists of a cyanide addition onto HCN. We have explored whether other nitriles could act as initiators for circumventing this first step. To accelerate reactivity, cyanide additions onto these nitriles would need to have lower barriers than the addition onto HCN. Possible initiators were selected following two criteria; they should be 1) small and 2) contain different functional groups attached to the -CN group. Figure 8 shows barrier heights for CN-addition relative to HCN (ΔΔE ). Only the electronic (Born-Oppenheimer) energy contribution is considered here, and we assume enthalpic and entropic effects to be relatively minor. The reference states for these barrier heights are free reactants. Calculations were made using Gaussian 16's default PCM model for water at the B3LYP-D3/6-311++G(d,p) level of theory. Figure S8. Activation energy barrier height for cyanide addition with a selection of small nitriles, shown relative to the reaction with HCN (ΔΔE , kcal/mol).