Lawrence Livermore National Laboratory



Nir Goldman (16-LW-020)

Abstract

This research project focused on the development of new molecular-dynamics models to capture chemical reactivity in condensed phases under extreme conditions, including astrobiological synthesis in prebiotic mixtures. We created a completely new molecular-dynamics force-field development scheme for atomistic simulations of reactive materials. Our models, which include many-body interactions, are generated by fitting linear combinations of Chebyshev polynomials through force matching to trajectories from Kohn-Sham Density Functional Theory. This approach enabled an increase of five orders of magnitude in computational efficiency over standard density functional theory while retaining most of its accuracy. We applied our method to a number of different materials and conditions, including liquid water under dissociative conditions, metallic liquid carbon, and amino acid oligomerization in aqueous solution. The simplicity of our approach enables the extension of quantum simulations to experimental time and length scales for any type of reactive systems of interest.

Background and Research Objectives

The origins of simple prebiotic hydrocarbons was first studied over 60 years ago, when electric discharges in a model of an early terrestrial system of methane, water, ammonia, and hydrogen led to milligram quantities of amino acids including glycine, alanine, and alpha amino-N-butyric acid (Miller 1953; Miller et al 1959). Reducing mixtures subjected to ultraviolet radiation (Sagan and Kahare 1971; Bernstein et al 2002; Muñoz Caro et al 2002), shock compression (Bar Nun et al 1970; McKay and Borucki 1997), and geothermal conditions (Harada and Fox 1964; LaRowe and Regnier 2008) produced prebiotic organic material such as amino acids and precursors to DNA and RNA nucleobases. However, early Earth's atmosphere was probably a more oxidizing environment (Brack 2007), with CO2 as its major constituent (Engel 2011). An inventory of possible life-building chemical reactivity on early Earth has yet to be established.

Planetary prebiotic chemistry may have been triggered by cometary impact (Oro 1961; Chiba et al 1990) (Figure 1).


Figure 1.
Figure 1. Proposed mechanism for impact synthesis in a prototypical cometary ice mixture (red = oxygen, white = hydrogen, blue = nitrogen, grey = carbon). The images show three sets (regions) of thermodynamic conditions. Shock compression, or impact, (Region I) caused the formation of large, carbon–nitrogen bonded nanoclusters (right). Upon adiabatic expansion (Region II) and cooling to ambient conditions (Region III), these clusters broke apart to form exotic amino acids (upper left corner) and long-chain amines (right side).

Comets and astrophysical ices are heterogeneous materials predominantly composed of water, with significant concentrations of both oxidized and reduced forms of carbon (Mumma and Charnley 2011). These compounds can comprise as much as 22 percent by concentration of a comet (Fomenkova et al 1994; Herd et al 2011) and include simple alkanes such as methane (CH4) and ethane (C2H6), and small-chain alcohols such as methanol (CH3OH). During the period of Late Heavy Bombardment, approximately four billion years ago, as much as 1,013 kg/yr of extraterrestrial organic material could have been delivered to early Earth, yielding several orders-of-magnitude greater mass of organic materials than previously existed on the planet (Blank et al 2001). However, to the best of our knowledge, a detailed survey of possible shock synthesis of complex prebiotic carbohydrates, from both exogenous and terrestrial starting material, has yet to be conducted.

Understanding shock synthesis in the laboratory may require the investigation of a large number of permutations of different starting materials and thermodynamic conditions. Frequently, insufficient data exists for the equation of state (EOS) and chemical reactivity of many materials under extreme conditions (Barrios et al 2010), and experiments require theoretical studies to elucidate measured physical and chemical properties. For example, studies on carbon-rich materials under pressures up to 30 GPa (1GPa = 10 kbar) indicate that slow chemical kinetics can extend beyond the time-scales of nanosecond laser-driven compression experiments (Manaa et al 2009). In addition, reported experimental temperatures can contain large uncertainties, making it difficult to adequately constrain the EOS on the basis of experiment alone (Goldman et al 2009).

Computer simulations such as molecular dynamics (MD) hold promise as an independent route to determining the EOS and chemical reactivity during reactive conditions. Such studies can provide simple chemical pictures of ionized intermediates and reaction mechanisms, and can help identify atomic-scale properties that determine observed macroscopic kinetics. These types of results can make experiments more tractable by aiding in their interpretation and helping to narrow the number of different materials and reactive conditions to investigate. Accurate modeling of the breaking and forming of bonds in condensed phases frequently requires the use of quantum theories such as density functional theory (DFT). DFT remains one of the most popular and widely used modeling methods for condensed-phase chemistry. It has been shown to accurately reproduce the phase boundaries and thermal decomposition of many materials (Goldman et al 2009), particularly at higher thermodynamic conditions like impact events (Goldman 2015). Our initial quantum simulations on a carbon-rich astrophysical mixture have shown the possibility of formation of amino acids due to shock synthesis (Goldman 2010; Goldman and 2013; Koziol and Goldman 2015). A subset of our results has since been confirmed by shock-recovery experiments on a similar mixture, where several linear and methyl-amino acids (including glycine) have been produced (Martins et al 2013). However, due to the extreme computational cost of these methods, many of our simulations were limited to tens-of-picoseconds time-scales and nanometer system sizes (Goldman et al 2010), far less than the nanosecond time scales and micrometer length scales needed to accurately model the chemical reactivity (Goldman and Bastea 2014). In addition, our simulations were limited to two specific icy mixtures only, and experiments so far have been designed to detect synthesized amino acids exclusively therein. Further experimentation would greatly benefit from detailed knowledge of the types of prebiotic molecules in addition to amino acids that can occur as a function of different starting mixtures as well as peak thermodynamic conditions.

Difficulty arises in determining models for chemical bonding that are both accurate and computationally efficient. Classical MD force fields treat interatomic interactions by analytic functional forms that are parameterized to specific physical properties over narrow ranges of thermodynamic conditions (van Duin et al 2001). These approaches usually exhibit a substantial increase in computational efficiency and scalability over standard quantum-simulation methods, enabling simulation of real experimental time and length scales (up to nanoseconds and micrometers, respectively). However, classical models are generally fit to the energetics of small molecular clusters or very specific bulk data; whereas the potential energy surface in these systems is rapidly evolving as new chemical products are formed, including large C–N bonded oligomers can coagulate and dissociate on less than picosecond timescales (Goldman et al 2010; Goldman and Tamblyn 2013; Koziol and Goldman 2015). This high degree of variability poses one of the most challenging aspects of creating an MD model for an astrobiological system under impact or other reactive conditions.

We were able to achieve a significant improvement in classical force-field development through force matching. Force matching (Figure 2) maximizes the data set that can be obtained from DFT by fitting parameters of a potential function (linear in this case, as discussed below) to each individual atomic force in an MD trajectory (Izvekov et al 2004), yielding a large quantity of data points for determining atomic level models (i.e., 3N data points per configuration, where N is the number of atoms in the system). Our approach is more general and simpler than previous efforts, enabling the efficient determination of accurate MD models for reactive systems through comparison with high-level quantum calculations.


Figure 2.
Figure 2. Flow chart for determination of our self-adaptive force field. Classical potential parameters are updated at fixed intervals (through force matching) to results from density-functional-theory (DFT) calculations.

Scientific Approach and Accomplishments

The MD approach we developed is referred to as the Chebyshev Interaction Model for Efficient Simulation (“ChIMES”) method (Figure 3). In this model, interatomic forces are represented by linear combinations of two-body and three-body Chebyshev polynomials. Our model can be self-adaptive, where potential-energy parameters can be modified during a simulation at fixed intervals through systematic comparison to high-level quantum calculations. This reduced reliance on quantum calculations results in an increase in computational efficiency of 5 orders of magnitude while retaining most of their accuracy. Through creation of the ChIMES formalism, we have established a general, high-throughput capability for the simulation of reactive chemistry that is accurate for a wide range of conditions related to shock-compressed organic materials.


Figure 3.
Figure 3. Vibrational spectra (left) and self-diffusion coefficients (right) for our models (ChIMES) for liquid metallic carbon at 15 GPa and 5,000 K. Our results showed markedly improved density functional theory (DFT), compared to previous molecular dynamics models. (ChIMES = Chebyshev Interaction Model for Efficient Simulation; LCBOP = long-range carbon bond-order potential; REBO = reactive empirical bond-order.)

We performed a DFT quantum simulation on a reactive mixture to determine initial configurations for our reactive force field (5–10 ps in length). We relied on the Vienna Ab initio Simulation Package (VASP; https://www.vasp.at) simulation suite for calculations of hundreds of atoms. A quantum code with linear computational scaling with system size (OpenMX; http://www.openmx-square.org) can be used to enable calculation of large system sizes (tens of thousands of atoms). ChIMES reactive force fields can then be determined through an iterative approach (Figure 2). Here, parameters for our interaction potentials can be computed through periodic comparison to the forces from our quantum calculations. We found that Chebyshev polynomials provide a simple and versatile functional form for our models. Chebyshev polynomials have standard power laws where the polynomials are orthonormal and are restricted to values over the interval [-1,1], minimizing the size and correlations of the linear coefficients. Chebyshev polynomials possess recursive relations, enabling easy generation of arbitrarily high orders. Linear combinations enable the use of singular value decomposition (SVD) to determine exact potential parameters in a single step, removing the need for nonlinear optimization methods that can be trapped in a local minimum or can be too computationally intensive to be of use (e.g., simulated annealing). MD simulations are then conducted with this new potential for a fixed interval. Depending on the desired accuracy, additional configurations from these simulations can be fed into quantum calculations, and the force-field parameters can be recomputed. This is in contrast to standard classical force-field development, where function parameters are predetermined by comparison with a fixed data set, and are held constant throughout the course of a simulation. Our method can thus achieve improved accuracy by adapting to the varying chemical reactivity that is known to occur in these prebiotic systems.

Our efforts have spanned several application areas related to the chemistry of planetary interiors, organic materials under extreme conditions, and astrobiological synthesis. Our pertinent results include the following.

(1) We created a simulation approach for water under conditions from ambient conditions up to the extreme pressures and temperatures of planetary interiors, as well as carbon near the graphite–liquid–diamond triple point (Figure 3). Our simulations retain the accuracy of DFT while yielding significant improvements in computational efficiency and scalability. This method shows significant improvement over previous force field approaches.

(2) We determined molecular dynamics models for the combustion chemistry of a phenolic resin (Figure 4). Our models predict accurate small molecule synthesis, both in terms of reaction rates and long-time scale concentrations.


Figure 4.
Figure 4. Numbers of specific chemical species from the phenolic polymer at 3,500 K. The solid black line corresponds to results from density functional theory (DFT) molecular dynamics; the dashed line corresponds to a previous semi-empirical quantum parameterization; and the solid blue line corresponds to our interaction potentials from force matching to the first 9 ps of the DFT trajectory (red line). Results from our potentials extrapolate to longer time scales.

(3) We used our methods to determine the free energies of formation of diglcyine (glycine-glycine peptide) in aqueous solution at conditions near hydrothermal vents in oceans on early Earth (Figure 5). We also determined the formation of complex nitrogen-containing polycyclic aromatic hydrocarbons in impacting glycine aqueous mixtures. These structures can form the basis for synthesis of vitamins, enzymes, nucleases, and other potentially life-building materials.


Figure 5.
Figure 5. Free-energy (F) calculation for diglycine formation in aqueous solution. Our results indicate that glycine can readily polymerize under moderate conditions.

Impact on Mission

Creating computationally efficient tools that are broadly applicable to shocked systems and materials directly supports Lawrence Livermore National Laboratory's stockpile stewardship science. Our project's focus on the simulation of the behavior and performance of complex systems, such as the development of biomaterials under extreme conditions and development of new simulation tools, also aligns with the Laboratory's core competency of high-performance computing, simulation, and data science.

The simulation methods we developed support the Laboratory's leadership position in modeling materials under extreme conditions and large-scale quantum mechanical simulations. The techniques we developed are very general and could be applied to fusion target materials such as diamond and beryllium, as well as planetary fluids such as hydrogen and water. Other materials, such as iron, are planned for study in dynamic compression experiments conducted at the Laboratory's National Ignition Facility. Our simulation capability will aid in the design and interpretation of future experiments. The simulation capability also could be applied to the chemistry of strongly shocked energetic materials. These are essential inputs for advanced explosive models (such as the Cheetah thermochemical code developed at the Laboratory) that are unattainable with conventional simulation techniques. Studies of the synthesis of prebiotic hydrocarbons under extreme pressures directly address several areas outlined in the NASA Exobiology Program Strategic Road Map (http://astrobiology.nasa.gov/roadmap). Shockless compression of materials under less extreme conditions occurs in a wide range of energy production environments, and will be of interest to DOE’s Basic Energy Science Program.

Conclusions

We have developed a high-throughput, highly scalable approach for simulation of chemical reactivity in condensed phases that will circumvent the critical time and length scale limitations of standard quantum codes while retaining most of their accuracy. Our efforts will help determine the self-assembly of carbon into life-building materials that could have been formed during impact events on early Earth. Computational chemical exobiology is an area of growing scientific and popular interest, and our group has established the Laboratory as a world leader in this area. The self-adaptive force-matching method could be successfully transitioned to a number of programs at Livermore, including studies of shocked explosives, polymers, and composite materials. The proposed capability may lead to collaborative partnerships with the NASA Exobiology program. We also anticipate possible collaborations with the Department of Defense to study the reaction mechanisms of novel energetic formulations.

References

Bar-Nun, A., et al. 1970. "Shock Synthesis of Amino Acids in Simulated Primitive Environments." Science 168:470.

Barrios, M. A., et al. 2010. "High-Precision Measurements of the Equation of State of Hydrocarbons at 1–10 Mbar Using Laser-Driven Shock Waves." Physics of Plasmas 17:056307.

Bernstein, M. P., et al. 2002. "Racemic Amino Acids from the Ultraviolet Photolysis of Interstellar Ice Analogues." Nature 416:401.

Blank, J. G., et al. 2001. "Experimental Shock Chemistry of Aqueous Amino Acid Solutions and the Cometary Delivery of Prebiotic Compounds." Origins of Life and Evolution of Biospheres B 31(15).

Brack, A. 2007. "From Interstellar Amino Acids to Prebiotic Catalytic Peptides: A Review." Chemistry and Biodiversity 4:665.

Chyba, C. F., et al. 1990. "Cometary Delivery of Organic Molecules to Early Earth." Science 249:366.

Engel, M. H. 2011. "Pre-Biotic Organic Synthesis: Laboratory Simulation Experiments and Their Significance for the Origin of Life in the Solar System." Proceedings of SPIE 8152:815208.

Fomenkova, M. N., et al. 1994. "Carbonaceous Components in the Comet Halley Dust." Geochimica et Cosmochimica Acta 58:4503–4512.

Goldman, N., et al. 2009. "Ab Initio Simulation of the Equation of State and Kinetics of Shocked Water." Journal of Chemical Physics 130:124517.

——— 2010. "Synthesis of Glycine-Containing Complexes in Impacts of Comets on Early Earth." Nature Chemistry 2:949.

Goldman, N. and I. Tamblyn. 2013. "Prebiotic Chemistry within a Simple Impacting Icy Mixture." Journal of Physical Chemistry A 117:5124–5131.28.

Goldman, N., and S. Bastea. 2014. "Nitrogen Oxides As a Chemistry Trap in Detonating Oxygen-Rich Materials." Journal of Physical Chemistry A 118, 2897—2903.

Goldman, N. 2015. "Multi-Center Semi-Empirical Quantum Models for Carbon Under Extreme Thermodynamic Conditions." Chemical Physics Letters 622:128.

Harada, K., and S. W. Fox. 1964. "Thermal Synthesis of Natural Amino-Acids from Postulated Primitive Terrestrial Atmosphere." Nature 201:335–336.

Herd, C. D. K., et al. 2011. "Origin and Evolution of Prebiotic Organic Matter as Inferred from the Tagish Lake Meteorite." Science 332:1304–1307.

Izvekov, S., et al. 2004. "Effective Force Fields for Condensed Phase Systems from Ab Initio Molecular Dynamics Simulation: A New Method for Force-Matching." Journal of Chemical Physics 120:10896.

Koziol, L., and N. Goldman. 2015. "Prebiotic Hydrocarbon Synthesis in Impacting Reduced Astrophysical Icy Mixtures." Astrophysical Journal 803:91.

LaRowe, D. E., and P. Regnier. 2008. "Thermodynamic Potential for the Abiotic Synthesis of Adenine, Cytosine, Guanine, Thymine, Uracil, Ribose and Deoxyribose in Hydrothermal Systems." Origins of Life and Evolution of the Biosphere 38:383–397.

Manaa, M. R., et al. 2009. "Nitrogen-Rich Heterocycles as Reactivity Retardants in Shocked Insensitive Explosives." Journal of the American Chemical Society 131:5493.

Martins, Z., et al. 2013. "Shock Synthesis of Organics from Simple Ice Mixtures." Nature Geoscience 6:1045.

McKay, C. P., and W. J. Borucki. 1997. "Organic Synthesis in Experimental Impact Shocks." Science 276:390.

Miller, S. L. 1953. "A Production of Amino Acids Under Possible Primitive Earth Conditions." Science 117:528.

Miller, S. L., and H. C. Urey. 1959. "Organic Compound Synthesis on the Primitive Earth." Science 130, 245.

Muñoz Caro, G. M., et al. 2002. Amino Acids from Ultraviolet Irradiation of Interstellar Ice Analogues." Nature 416:403.

Mumma, M. J., and S. B. Charnley. 2011. "The Chemical Composition of Comets - Emerging Taxonomies and Natal Heritage." Annual Review of Astronomy and Astrophysics 49:471–524.

Oro, P. J. 1961. "Comets and the Formation of Biochemical Compounds on the Primitive Earth." Nature 190:389–390.

Sagan, C., and B. N. Khare. 1971. "Long-Wavelength Ultraviolet Photoproduction of Amino Acids on the Primitive Earth." Science 173:417–420.

van Duin, A. C. T., et al. 2001. "ReaxFF: A Reactive Force Field for Hydrocarbons." Journal of Physical Chemistry A 105: 9396.

Publications and Presentations

Goldman, N., et al. 2015. “Using Force-Matched Potentials to Improve the Accuracy of Density Functional Tight Binding for Reactive Conditions.” Journal of Chemical Theory and Computation 11(10): 150909123846008. doi:10.1021/acs.jctc.5b00742. LLNL-JRNL-675158.

Koziol, L., et al. 2017. “Using Force-Matching to Determine Reactive Force Fields for Bulk Water Under Extreme Thermodynamic Conditions.” Journal of Chemical Theory and Computation 13(1):135–146. doi:10.1021/acs.jctc.6b00707. LLNL-JRNL-696539.