Complex Potential Energy Surface Research Articles

Nondirect-Product Local Diabatic Representation with Smolyak Sparse Grids.

Modeling nonadiabatic conical intersection dynamics is critical for understanding a wide range of photophysical, photochemical, and biological phenomena. Here we develop a nonadiabatic conical intersection wave packet dynamic method in the local diabatic representation using Smolyak sparse grids. Employing sparse grids avoids the direct-product grids in configuration space and alleviates the exponential scaling of computation costs with the molecular size. Numerical demonstrations are first performed for a two-dimensional vibronic model of pyrazine, where the results using sparse grids are in excellent agreement with those using direct-product grids, with sparse grids being much faster. Moreover, we demonstrate that for a four-dimensional pyrazine model, where direct-product grids are computationally infeasible, sparse grids can provide almost exact results. The sparse grid local diabatic representation method is further applied to a realistic model system of phenol photodissociation with much more complex potential energy surfaces; the results using sparse grids still agree very well with the direct-product grids. Finally, by combining with electronic structure calculations, we apply our method to the Shin-Metiu model without quasi-diabatization. The sparse grid and direct-product grid results are in good agreement, with the sparse grid computational cost being half of the full grid.

Theoretical and experimental study on the pyrolysis of N-methylpyrrolidone

As a widely used organic solvent, N-methylpyrrolidone (NMP) may suffer thermal decomposition during large-scale industrial usage. The complex potential energy surface and corresponding pressure-dependent rate coefficients of NMP and NMP radicals were obtained at a high level of theory. The intermediates and products of NMP pyrolysis within 900–1200 K were identified and measured by synchrotron vacuum ultraviolet photoionization mass spectrometry. Theoretical calculations show that the Keto-enol tautomerization and direct ring-opening channel yielding C2H4, CO, and CH3NCH2 exhibit the highest rate coefficient among NMP unimolecular reactions at low and high temperatures, respectively. However, the fuel decomposition process cannot be exclusively dominated by single low-barrier channel. The ring-opening reaction forming imino aldehyde radical and following decarbonylation reactions are favored channels during the NMP radical decomposition, while the formation of the dihydropyrrolidone via H or CH3 loss competes with them at higher temperatures. Various smaller heterocycle radicals can be produced and act as intermediates for the multi-step isomerization of the ring-opening products with a relatively low energy barrier, changing the position of the functional group and radical site. A variety of nitrogenous and oxygenated species were observed in the experiment, such as ketene, methyl isocyanate, pyrrole, benzene, and different kinds of imines. Among them, C2H4, HCN, and CO are the major products for NMP decomposition. Potential reaction pathways from NMP to the major experimentally observed products were inferred based on the calculated results and existing knowledge on the pyrolysis of nitrogenous and oxygenated compounds. These results shed light on the complex C/H/O/N reaction system. The calculated rate coefficients and extensive mole fractions data could provide good support for the modeling of a detailed reaction mechanism in the future.

Minimizing Redundancy and Data Requirements of Machine Learning Potential: A Case Study in Interface Combustion.

The machine learning potential has emerged as a promising approach for addressing the accuracy-versus-efficiency dilemma in molecular modeling. Efficiently exploring chemical spaces with high accuracy presents a significant challenge, particularly for the interface reaction system. This study introduces a workflow aimed at achieving this goal by incorporating the classical SOAP descriptor and practical PCA strategy to minimize redundancy and data requirements, while successfully capturing the features of complex potential energy surfaces. Specifically, the study focuses on interface combustion behaviors within promising alloy-based solid propellants. A neural network potential model tailored for modeling AlLi-AP interface reactions under varying conditions is constructed, showcasing excellent predictive capabilities in energy prediction, force estimation, and bond energies. A series of large-scale MD simulations reveal that Li doping significantly influences the initial combustion stage, enhancing reactivity and reducing thermal conductivity. Mass transfer analysis also highlights the considerably higher diffusion coefficient of Li compared to Al, with the former being three times greater. Consequently, the overall combustion process is accelerated by approximately 10%. These breakthroughs pave the way for virtual screening and the rational design of advanced propellant formulations and microstructures incorporating alloy-formula propellants.

Complex potential energy surfaces with projected CAP technique: Vibrational excitation of N2.

The projected complex absorbing potential (CAP) technique is one of the methods that allow one to extend the bound state methods for computing resonances' energies and widths. Here, we explore the accuracy of the potential energy curves generated with different electronic structure theory methods in combination with the projected CAP technique by considering resonant vibrational excitation (RVE) of N2 by electron impact as a model process. We report RVE cross sections computed using the boomerang model with potential energy curves obtained with CAP-based extended multistate complete active space perturbation theory (XMS-CASPT2) and equation of motion coupled-cluster method for electron attachment with single and double substitution (EOM-EA-CCSD) methods. We also compare potential energy curves computed with several electronic structure methods, including XMS-CASPT2, EOM-EA-CCSD, multireference configuration interaction with singles (MR-CIS) and singles and doubles (MR-CISD). A good agreement is observed between the experiment and simulated RVE cross sections obtained with the potential energy curves generated with XMS-CASPT2 and EOM-EA-CCSD methods, thus highlighting the potential of the projected CAP technique combined with accurate electronic structure methods for dynamical simulations of the processes that proceed through metastable electronic states.

Understanding the Vibrational Structure and Ultrafast Dynamics of the Metal Carbonyl Precatalyst [Mn(ppy)(CO)4

The solution phase structure, vibrational spectroscopy, and ultrafast relaxation dynamics of the precatalyst species [Mn(ppy)(CO)4] (1) in solution have been investigated using ultrafast two-dimensional infrared (2D-IR) spectroscopy. By comparing 2D-IR data with the results of anharmonic density functional theory (DFT) calculations, we establish an excellent agreement between measured and predicted inter-mode couplings of the carbonyl stretching vibrational modes of 1 that relates to the atomic displacements of axial and equatorial ligands in the modes and the nature of the molecular orbitals involved in M-CO bonding. Measurements of IR pump-probe spectra and 2D-IR spectra as a function of waiting time reveal the presence of ultrafast (few ps) intramolecular vibrational energy redistribution between carbonyl stretching modes prior to vibrational relaxation. The vibrational relaxation times of the CO-stretching modes of 1 are found to be relatively solvent-insensitive, suggestive of limited solvent-solute interactions in the ground electronic state. Overall, these data provide a detailed picture of the complex potential energy surface, bonding and vibrational dynamics of 1, establishing a fundamental basis for the next steps in understanding and modulating precatalyst behavior.

Theoretical kinetic investigation of the multichannel mechanism of O(3P) atmospheric oxidation reaction of but-3-enal

Several levels of theory such as Møller-Plesset MP2, G3, and CBS-QB3, have been used in order to investigate the complex and multichannel potential energy surface of the reaction of but-3-enal with the triplet oxygen atom. The results show that the O-addition channel is dominant. The different possible pathways of oxygen atom attack are thoroughly studied to better understand and explain the reaction mechanism. Regarding the oxidation of but-3-enal by triplet oxygen O(3P), it is shown that the major thermodynamic product is H3CC(O)CH2C(O)H (P3) being the most stable for the whole reaction. However, the most favored product kinetically is H2CC(OH)CH2C(O)H (P2). For the H-abstraction second possible pathway, the most favored product both kinetically and thermodynamically is found to be P8. The activation energy and calculated rate constants are consistent with the proposed addition mechanism.

Computational Investigations of the Detailed Mechanism of Reverse Intersystem Crossing in Inverted Singlet-Triplet Gap Molecules.

Inverted singlet-triplet gap (INVEST) materials have promising photophysical properties for optoelectronic applications due to an inversion of their lowest singlet (S1) and triplet (T1) excited states. This results in an exothermic reverse intersystem crossing (rISC) process that potentially enhances triplet harvesting, compared to thermally activated delayed fluorescence (TADF) emitters with endothermic rISCs. However, the processes and phenomena that facilitate conversion between excited states for INVEST materials are underexplored. We investigate the complex potential energy surfaces (PESs) of the excited states of three heavily studied azaphenalene INVEST compounds, namely, cyclazine, pentazine, and heptazine using two state-of-the-art computational methodologies, namely, RMS-CASPT2 and SCS-ADC(2) methods. Our findings suggest that ISC and rISC processes take place directly between the S1 and T1 electronic states in all three compounds through a minimum-energy crossing point (MECP) with an activation energy barrier between 0.11 to 0.58 eV above the S1 state for ISC and between 0.06 and 0.36 eV above the T1 state for rISC. We predict that higher-lying triplet states are not populated, since the crossing point structures to these states are not energetically accessible. Furthermore, the conical intersection (CI) between the ground and S1 states is high in energy for all compounds (between 0.4 to 2.0 eV) which makes nonradiative decay back to the ground state a relatively slow process. We demonstrate that the spin-orbit coupling (SOC) driving the S1-T1 conversion is enhanced by vibronic coupling with higher-lying singlet and triplet states possessing vibrational modes of proper symmetry. We also rationalize that the experimentally observed anti-Kasha emission of cyclazine is due to the energetically inaccessible CI between the bright S2 and the dark S1 states, hindering internal conversion. Finally, we show that SCS-ADC(2) is able to qualitatively reproduce excited state features, but consistently overpredict relative energies of excited state structural minima compared to RMS-CASPT2. The identification of these excited state features elaborates design rules for new INVEST emitters with improved emission quantum yields.

Identifying the charge states of carbon vacancies in 4H-SiC by ab initio metadynamics

4H Silicon carbide (4H-SiC) is widely recognized as a highly promising material for high-voltage and high-power electronic applications due to its exceptional properties. The performance of devices based on 4H-SiC is often weakened by the presence of carbon-related point defects, particularly carbon vacancies (VC). The defects of VC introduce deep-level traps (e.g., Z1/2 and EH6/7) that deteriorate device functionality. Experimental and theoretical studies on VC have led to some conflicting results about the charge states of VC, especially for the charge state ordering of EH6/7. We now employ ab initio metadynamics (META) to systematically investigate configuration space including the direction and magnitude of bond distortion and identify the most stable structures of VC. Eventually, the charge states of VC in 4H-SiC are identified. The Z1 (EH6) and Z2 (EH7) indicate transitions from acceptor (donor) levels of VC, located on the h and k sublattice sites, respectively. Z1 and Z2 demonstrate negative-U ordering, characterized by U values of −0.16 and −0.37 eV, respectively. Conversely, EH6 and EH7 display positive-U ordering, with U values of 0.16 and 0.08 eV, respectively. The current results provide insights into the properties of VC in 4H-SiC, highlighting the effectiveness of META in the exploration of complex potential energy surfaces associated with point defects in solids.

Mechanism of Asymmetric Homologation of Alkenylboronic Acids with CF3-Diazomethane via Borotropic Rearrangement.

Density functional theory calculations have been performed to investigate the mechanism for the BINOL-catalyzed asymmetric homologation of alkenylboronic acids with CF3-diazomethane. The reaction proceeds via a chiral BINOL ester of the alkenylboronic acid substrate. The calculations reveal a complex scenario for the formation of the chiral BINOL-alkenylboronate species, which is the key intermediate in the catalytic process. The aliphatic alcohol additive plays an important role in the reaction. This study provides a rationalization of the stereoinduction step of the reaction, and the enantioselectivity is mainly attributed to the steric repulsion between the CF3 group of the diazomethane reagent and the γ-substituent of the BINOL catalyst. The complex potential energy surface obtained by the calculations is analyzed by means of microkinetic simulations.

Benchmark ab initio characterization of the complex potential energy surfaces of the HOO- + CH3Y [Y = F, Cl, Br, I] reactions.

The α-effect is a well-known phenomenon in organic chemistry, and is related to the enhanced reactivity of nucleophiles involving one or more lone-pair electrons adjacent to the nucleophilic center. The gas-phase bimolecular nucleophilic substitution (SN2) reactions of α-nucleophile HOO- with methyl halides have been thoroughly investigated experimentally and theoretically; however, these investigations have mainly focused on identifying and characterizing the α-effect of HOO-. Here, we perform the first comprehensive high-level ab initio mapping for the HOO- + CH3Y [Y = F, Cl, Br and I] reactions utilizing the modern explicitly-correlated CCSD(T)-F12b method with the aug-cc-pVnZ [n = 2-4] basis sets. The present ab initio characterization considers five distinct product channels of SN2: (CH3OOH + Y-), proton abstraction (CH2Y- + H2O2), peroxide ion substitution (CH3OO- + HY), SN2-induced elimination (CH2O + HY + HO-) and SN2-induced rearrangement (CH2(OH)O- + HY). Moreover, besides the traditional back-side attack Walden inversion, the pathways of front-side attack, double inversion and halogen-bond complex formation have also been explored for SN2. With regard to the Walden inversion of HOO- + CH3Cl, the previously unaddressed discrepancies concerning the geometry of the corresponding transition state are clarified. For the HOO- + CH3F reaction, the recently identified SN2-induced elimination is found to be more exothermic than the SN2 channel, submerged by ∼36 kcal mol-1. The accuracy of our high-level ab initio calculations performed in the present study is validated by the fact that our new benchmark 0 K reaction enthalpies show excellent agreement with the experimental data in nearly all cases.

Global Neural Network Potential with Explicit Many-Body Functions for Improved Descriptions of Complex Potential Energy Surface.

The high dimensional machine learning potential (MLP) that has developed rapidly in the past decade represents a giant step forward in large-scale atomic simulation for complex systems. The long-range interaction and the poor description of chemical reactions are typical problems of high dimensional MLP, which are mainly caused by the poor structure discrimination of the atom-centered ML model. Herein, we propose a low-cost neural-network-based MLP architecture for fitting global potential energy surface data, namely, G-MBNN, that can offer improved energy and force resolution on a complex potential energy surface. In G-MBNN, a set of many-body energy terms based on the local atomic environment are explicitly included in computing the total energy─the total energy of the system is written as the sum of atomic energy and many-body energy contributions. These extra many-body energy terms are computationally low-cost and, importantly, can provide easy access to delicate energy terms in complex systems such as very short repulsion, long-range attractions, and sensitive angular-dependent covalent interactions. We implement G-MBNN in the LASP code and demonstrate the improved accuracy of the new framework in representative systems, including ternary-element energy materials LiCoOx, TiO2 with defects, and a series of organic reactions.

Fragment-based models for dissociation of strong acids in water: Electrostatic embedding minimizes the dependence on the fragmentation schemes.

Fragmentation methods such as MIM (Molecules-in-Molecules) provide a route to accurately model large systems and have been successful in predicting their structures, energies, and spectroscopic properties. However, their use is often limited to systems at equilibrium due to the inherent complications in the choice of fragments in systems away from equilibrium. Furthermore, the presence of charges resulting from any heterolytic bond breaking may increase the fragmentation error. We have previously suggested EE-MIM (Electrostatically Embedded Molecules-In-Molecules) as a method to mitigate the errors resulting from the missing long-range interactions in molecular clusters in equilibrium. Here, we show that the same method can be applied to improve the performance of MIM to solve the longstanding problem of dependency of the fragmentation energy error on the choice of the fragmentation scheme. We chose four widely used acid dissociation reactions (HCl, HClO4, HNO3, and H2SO4) as test cases due to their importance in chemical processes and complex reaction potential energy surfaces. Electrostatic embedding improves the performance at both one and two-layer MIM as shown by lower EE-MIM1 and EE-MIM2 errors. The EE-MIM errors are also demonstrated to be less dependent on the choice of the fragmentation scheme by analyzing the variation in fragmentation energy at the points with more than one possible fragmentation scheme (points where the fragmentation scheme changes). EE-MIM2 with M06-2X as the low-level resulted in a variation of less than 1 kcal/mol for all the cases and 1 kJ/mol for all but three cases, rendering our method fragmentation scheme-independent for acid dissociation processes.

Semiempirical Methods for Molecular Systems in Strong Magnetic Fields.

A general scheme is presented to extend semiempirical methods to include the effects of arbitrary strength magnetic fields, while maintaining computational efficiency. The approach utilizes three main modifications; a London atomic orbital (LAO) basis set is introduced, field-dependent kinetic energy corrections are added to the model Hamiltonian, and spin-Zeeman interaction energy terms are included. The approach is applied to the widely available density-functional tight-binding method GFN1-xTB. Considering the basis set requirements for the kinetic energy corrections in a magnetic field leads to two variants: a single-basis approach GFN1-xTB-M0 and a dual-basis approach GFN1-xTB-M1. The LAO basis in the latter includes the appropriate nodal structure for an accurate representation of the kinetic energy corrections. The variants are assessed by benchmarking magnetizabilities and nuclear magnetic resonance shielding constants calculated using weak magnetic fields. Remarkably, the GFN1-xTB-M1 approach also exhibits excellent performance for strong fields, || ≤ 0.2B0 (B0 = 2.3505 × 105 T), recovering exotic features such as the para- to dia-magnetic transition in the BH molecule and the preferred electronic configuration, molecular conformation, and orientation of benzene. At stronger field strengths, || > 0.2B0, a degradation in the quality of the results is observed. The utility of GFN1-xTB-M1 is demonstrated by performing conformer searches in a range of field strengths for the cyclooctatetraene molecule, with GFN1-xTB-M1 capturing the transition from tub to planar conformations at high field, consistent with much more computationally demanding current-density functional theory calculations. Magnetically induced currents are also shown to be well described for the benzene and infinitene molecules, the latter demonstrating the flexibility and computational efficiency of the approach. The GFN1-xTB-M1 approach is a useful tool for the study of structure, conformation, and dynamics of large systems in magnetic fields at the semiempirical level as well as for preoptimization of molecular structure in ab initio calculations, enabling more efficient exploration of complex potential energy surfaces and reactivity in the presence of external fields.

Machine learning meets quantum mechanics in catalysis

Over the past decade many researchers have applied machine learning algorithms with computational chemistry and materials science tools to explore properties of catalysts. There is a rapid increase in publications demonstrating the use of machine learning for rational catalyst design. In our perspective, targeted tools for rational catalyst design will continue to make significant contributions. However, the community should focus on developing high-throughput simulation tools that utilize molecular dynamics capabilities for thorough exploration of the complex potential energy surfaces that exist, particularly in heterogeneous catalysis. Catalyst-specific databases should be developed to contain enough data to represent the complex multi-dimensional space that defines structure-function relationships. Machine learning tools will continue to impact rational catalyst design; however, we believe that more sophisticated pattern recognition algorithms would yield better understanding of structure-function relationships for heterogeneous catalysis.

Machine-learning driven global optimization of surface adsorbate geometries

The adsorption energies of molecular adsorbates on catalyst surfaces are key descriptors in computational catalysis research. For the relatively large reaction intermediates frequently encountered, e.g., in syngas conversion, a multitude of possible binding motifs leads to complex potential energy surfaces (PES), however. This implies that finding the optimal structure is a difficult global optimization problem, which leads to significant uncertainty about the stability of many intermediates. To tackle this issue, we present a global optimization protocol for surface adsorbate geometries which trains a surrogate machine learning potential on-the-fly. The approach is applicable to arbitrary surface models and adsorbates and minimizes both human intervention and the number of required DFT calculations by iteratively updating the training set with configurations explored by the algorithm. We demonstrate the efficiency of this approach for a diverse set of adsorbates on the Rh(111) and (211) surfaces.

Robosample: A software using constrained molecular dynamics coupled with Gibbs sampling to investigate macromolecules.

An ongoing challenge for the simulation of macromolecules is their complex potential energy surface which is a result of high dimensionality and chemical complexity. The surface contains high energy barriers usually difficult to overcome by existing simulation methods. We present here Robosample, a software combining the advantages of robot mechanics for fast constrained simulations and of Gibbs sampling for achieving ergodicity. While energy calculation is the main computational bottleneck for macromolecule simulations, a big part of the work can be short-circuited by increasing the trajectory integration timestep. Usually this is achieved by modifying the system's mass matrix. The way this is implemented in Robosample is by using constraints in the framework of reduced coordinates. A good example is torsional dynamics, which transforms the system in Bond-Angle-Torsion coordinates and then ignores the motion of the bonds and the angles. Robosample sampling method is Hamiltonian Monte Carlo coupled with Gibbs sampling (GCHMC), with random variables divided among molecular degrees of freedom. With GCHMC, certain molecule segments are considered rigid and sub-sampled, while focusing on soft degrees of freedom as Gibbs blocks. Ergodicity is achieved by interleaving fully flexible Cartesian dynamics Monte Carlo proposals. The validity of the method was first proven on small systems correctly reproducing alanine dipeptide free energy landscape and showed to increase sampling efficiency of macrocycle containing small molecules which are notoriously difficult to simulate. The software was subsequently extended to include spherical and cylindrical robotic joints and proved to be efficient on medium sized biological molecules, such as glycans. Now we have reached a stage where we can successfully sample big macromolecules such as antibodies on graphical processing units.

PESPIP: Software to fit complex molecular and many-body potential energy surfaces with permutationally invariant polynomials.

We wish to describe a potential energy surface by using a basis of permutationally invariant polynomials whose coefficients will be determined by numerical regression so as to smoothly fit a dataset of electronic energies as well as, perhaps, gradients. The polynomials will be powers of transformed internuclear distances, usually either Morse variables, exp(-ri,j/λ), where λ is a constant range hyperparameter, or reciprocals of the distances, 1/ri,j. The question we address is how to create the most efficient basis, including (a) which polynomials to keep or discard, (b) how many polynomials will be needed, (c) how to make sure the polynomials correctly reproduce the zero interaction at a large distance, (d) how to ensure special symmetries, and (e) how to calculate gradients efficiently. This article discusses how these questions can be answered by using a set of programs to choose and manipulate the polynomials as well as to write efficient Fortran programs for the calculation of energies and gradients. A user-friendly interface for access to monomial symmetrization approach results is also described. The software for these programs is now publicly available.

Growth Pattern of Large Morse Clusters with Medium-Range Potentials.

Due to the extremely complex potential energy surfaces of large Morse clusters with medium-range potentials (i.e., ρ = 6 and 10), global optimization studies in the literature are limited to a cluster size (N) of ≤240. Starting from completely random structures, we successfully systematically studied Morse clusters with up to 700 atoms using our unbiased fuzzy global optimization (FGO) method. While all of the putative global minima reported previously have been efficiently obtained, new global minima with lower energies are identified for N values of 176, 258, 485, 561, 817, and 923 with ρ = 6 and for N values of 151, 202, 226, and 229 with ρ = 10. A detailed growth pattern and magic clusters are obtained. For the first time, we find that a central vacancy is present in Morse clusters containing 542, 543, 548, and 922 atoms with ρ = 6. FGO has achieved high performance in large clusters with different interatomic interaction ranges, thus showing great application potential in the global structure optimization of general clusters.

Automation of chemical kinetics: Status and challenges

Driven by synergic advancements in high performance computing and theory, the capability to estimate rate constants from first principles has evolved considerably recently. When this knowledge is coupled with a procedure to determine a list of all reactions relevant to describe the evolution of a reacting system, it becomes possible to envision a methodology to predict theoretically the reaction kinetics. However, if a thorough examination of all possible reaction channels is desired, the number of reactions for which a rate constant estimate is needed can become quite large. This determines the need for rate constant estimation automation. In the present work, the status of this rapidly evolving field is reviewed, with emphasis on recent advancements and present challenges. Thermochemistry is the field where automation is most advanced. Entropies, heat capacities, and enthalpies can be determined efficiently with accuracy comparable to experiments for most chemical species containing a limited number of atoms, while machine learning can be used to improve the computational predictions for large chemical species using reduced computational resources. Several approaches have been proposed to automatically investigate the reactivity over complex potential energy surfaces, while rate constants for elementary steps can be determined accurately for several reaction classes, such as abstraction, addition, beta-scission, and isomerization. Kinetic mechanisms can be automatically generated using methodologies that differ for level of complexity and required physical insight. Among the challenges that are still to be met are the estimation of rate constants for intrinsically multireference reaction classes, such as barrierless processes, the containment of the number of reactions to screen in mechanism development, and the integration of the existing automated software. It is suggested that the synergy between experiment and theory should evolve towards a stage where experiments are focused on the estimation of parameters where theoretical tools are least predictive, and vice versa.

Experimental study of the removal of excited state phosphorus atoms by H2O and H2: implications for the formation of PO in stellar winds

ABSTRACTThe reactions of the low-lying metastable states of atomic phosphorus, P(2D) and P(2P), with H2O and H2 were studied by the pulsed laser photolysis at 248 nm of PCl3, combined with laser-induced fluorescence detection of P(2D), P(2P), and PO. Rate coefficients between 291 and 740 K were measured, along with a yield for the production of PO from P(2D or 2P) + H2O of (35 ± 15) %. H2 reacts with both excited P states relatively efficiently; physical (i.e. collisional) quenching, rather than chemical reaction to produced PH + H, is shown to be the more likely pathway. A comprehensive phosphorus chemistry network is then developed using a combination of electronic structure theory calculations and a Master Equation treatment of reactions taking place over complex potential energy surfaces. The resulting model shows that at the high temperatures within two stellar radii of a MIRA variable AGB star in oxygen-rich conditions, collisional excitation of ground-state P(4S) to P(2D), followed by reaction with H2O, is a significant pathway for producing PO (in addition to the reaction between P(4S) and OH). The model also demonstrates that the PN fractional abundance in a steady (non-pulsating) outflow is underpredicted by about 2 orders of magnitude. However, under shocked conditions where sufficient thermal dissociation of N2 occurs at temperatures above 4000 K, the resulting N atoms convert a substantial fraction of PO into PN.