Constrained Molecular Dynamics Research Articles

Metal–N4 model single‐atom catalyst with electroneutral quadri‐pyridine macrocyclic ligand for CO2 electroreduction

AbstractMetal–N–C single‐atom catalysts, mostly prepared from pyrolysis of metal‐organic precursors, are widely used in heterogeneous electrocatalysis. Since metal sites with diverse local structures coexist in this type of material and it is challenging to characterize the local structure, a reliable structure–property relationship is difficult to establish. Conjugated macrocyclic complexes adsorbed on carbon support are well‐defined models to mimic the single‐atom catalysts. Metal–N4 site with four electroneutral pyridine‐type ligands embedded in a graphene layer is the most commonly proposed structure of the active site of single‐atom catalysts, but its molecular counterpart has not been reported. In this work, we synthesized the conjugated macrocyclic complexes with a metal center (Co, Fe, or Ni) coordinated with four electroneutral pyridinic ligands as model catalysts for CO2 electroreduction. For comparison, the complexes with anionic quadri‐pyridine macrocyclic ligand were also prepared. The Co complex with the electroneutral ligand expressed a turnover frequency of CO formation more than an order of magnitude higher than that of the Co complex with the anionic ligand. Constrained ab initio molecular dynamics simulations based on the well‐defined structures of the model catalysts indicate that the Co complex with the electroneutral ligand possesses a stronger ability to mediate electron transfer from carbon to CO2.

C−C Coupling in CO2 Electroreduction on Single Cu‐Modified Covalent Triazine Frameworks: A Static and Dynamic Density Functional Theory Study

AbstractThe electrochemical reduction of CO2 into C2+ products represents a promising solution to completing the carbon cycle, thereby fostering a sustainable energy supply. Single‐atom electrocatalysts (SAECs) have garnered significant attention as efficient electrocatalysts for the CO2 reduction reaction. Herein, we carried out a first‐principles study on the mechanism of C−C bond formation on single‐Cu‐atom‐modified covalent triazine frameworks (Cu‐CTFs), which are a promising platform for SAECs. Static density functional calculations indicated that the dimerization of CO, which is the main mechanism for C−C bond formation on bulk Cu metals, was not favorable for Cu‐CTFs because of the lack of adjacent Cu sites for co‐adsorption of CO molecules. Rather than CO dimerization, the reaction between adsorbed *CHO and CO to produce *COCHO has a relatively low reaction energy barrier. Constrained ab initio molecular dynamics analyses revealed that the C−C bond‐forming reaction proceeds via insertion of CO at the *CHO intermediate, which has a modest activation energy of 0.09 eV. Specifically, when CO molecule is constrained to be brought close to *CHO, CO insertion between *CHO and Cu occurs at a C−C distance of 1.8 Å. This insertion reaction is the transition step for this C−C bond formation.

Clumped isotope reordering kinetics in strontianite and witherite: Experiments and first-principles simulations

The carbonate clumped isotope (∆47) thermometer is a powerful tool for temperature reconstruction based on the thermodynamical 13C—18O bonding preference. Accurate interpretation of temperatures requires preservation of original carbonate ∆47 signatures, which are susceptible to resetting at elevated temperatures during burial. Kinetic models of carbonate 13C—18O bond reordering derived from experimental studies constrain the evaluation of ∆47 preservation and the reconstruction of thermal histories and initial temperatures for reordered carbonates. However, application of these models to different carbonates is limited because kinetic parameters vary with materials and are believed to be affected by factors such as crystal structure, cation, trace element composition, crystal defects, and internal water whose roles in reordering are still not well-understood. Here we present reordering data for strontianite and witherite from heating experiments and evaluate the effects of crystal structure, cation chemistry, cation vacancy and substitution, and internal water on reordering using constrained ab-initio molecular dynamics (cAIMD) simulations. We show that strontianite and witherite reorder significantly faster than other carbonates. Comparing experimental studies yields the following reordering rate sequence: strontianite > witherite > aragonite > calcite > dolomite. The general trend in Arrhenius activation energy for reordering is strontianite < witherite < aragonite < dolomite < calcite (with calcite values based on the average of individual materials). Our cAIMD simulations reveal a similar trend in free energy of activation for oxygen exchange in these carbonates. We suggest that orthorhombic carbonates are more susceptible to reordering than trigonal carbonates because of the shorter distance between carbonate groups and because the weaker M—O interaction in orthorhombic carbonates reduces the structural distortion and the energy cost during transition state formation. For carbonates sharing the same structure, the weak M—O bonding can be compensated by a large unit cell volume, leading to a non-systematic relationship between the free energy of activation and the cation radius. In simulations with lattice-bound water, the free energy of activation is greatly reduced compared with cation vacancy and substitution. The influence of fluid inclusions on measured reordering rates depends on the fraction of carbonates that can exchange with water as suggested by strontianite δ¹⁸O data. Our results serve as a starting point for a more sophisticated but flexible kinetic model containing material-specific factors (e.g., substituted ions, structure water, etc.) for reordering in abiotic and biotic carbonates.

Nuclear Equation of State and many-body phase-space correlations in the Constrained Molecular Dynamics

Many-body correlations characterizing the Constrained Molecular Dynamics (CoMD) are analyzed in the case of finite and zero-range effective microscopic interactions. The study begins by analyzing the case of infinite nuclear matter at zero temperature. A comparison with the predictions in the mean-field (MF) limit corresponding to different effective masses, highlights non-negligible differences regarding the produced Equation of State (EoS). A procedure is illustrated to determine the necessary corrections of the effective interaction parameters in the CoMD model so to reproduce the chosen EoS. The specific model calculations, the general feature of the discussed correlations gives a wider meaning to the resulting differences, which are in fact strongly related both to the Pauli principle constraint and to the localization effects related to wave-packets dynamics. Moving on to finite systems, preliminary results are shown in relation to the reaction mechanisms in the N64i+48Ca system described by the CoMD model. Therefore, the topic covered illustrates the effects produced by the correlations of the CoMD dynamics on the EoS and on some observables commonly studied in heavy ion collisions.

Many-body expansion for light nuclear systems

We extend the many-body expansion (MBE), previously applied to hydrogen bonded and molecular systems, to the light nuclear systems $^{3}\mathrm{H}$ and $^{3}\mathrm{He}$ by considering the nucleonic degrees of freedom as fundamental in the expansion. The analysis is based on the Pauli nucleonic dynamics (PND) model, a simple antisymmetrized dynamical code, inspired by the sophisticated constrained molecular dynamics (CoMD) model. The total energy of the $^{2}\mathrm{H}$ nucleus is calculated with this model at \ensuremath{-}2.312 MeV, which is within 4% of the experimental value of \ensuremath{-}2.225 MeV. The application of the MBE yields results for the three-body term in the $^{3}\mathrm{H}$ nucleus that is comparable with previous estimates, while it is reported for the first time for the nucleus of $^{3}\mathrm{He}$. The energies of $^{3}\mathrm{H}$ and $^{3}\mathrm{He}$ that include the sum of the one- and two-body terms with the model are $\ensuremath{-}6.97\ifmmode\pm\else\textpm\fi{}0.21$ MeV and $\ensuremath{-}6.19\ifmmode\pm\else\textpm\fi{}0.21$ MeV and the three-body terms, estimated from the difference of the sum of the one- and two-body terms from the experimentally measured energies, are $\ensuremath{-}1.51\ifmmode\pm\else\textpm\fi{}0.21$ MeV and $\ensuremath{-}1.53\ifmmode\pm\else\textpm\fi{}0.21$ MeV, respectively. The MBE for the $^{3}\mathrm{H}$ and $^{3}\mathrm{He}$ nuclei is qualitatively similar to the one previously reported for the water trimer: the two- and three-body terms are negative with the latter amounting to $\ensuremath{\sim}3%$ of the former. Additionally, the three-body terms correspond to about 16--23% of the total energy of the bound systems, a percentage that is also comparable to the one in the water trimer (17%). In this manner, the MBE analysis can be applied to light nuclear systems following the same protocol as the one that has been previously extensively used for hydrogen bonded molecular systems.

Accurate and efficient constrained molecular dynamics of polymers using Newton's method and special purpose code

In molecular dynamics simulations we can often increase the time step by imposing constraints on bond lengths and bond angles. This allows us to extend the length of the time interval and therefore the range of physical phenomena that we can afford to simulate. We examine the existing algorithms and software for solving nonlinear constraint equations in parallel and we explain why it is necessary to advance the state-of-the-art. We present ILVES-PC, a new algorithm for imposing bond constraints on proteins accurately and efficiently. It solves the same system of differential algebraic equations as the celebrated SHAKE algorithm, but ILVES-PC solves the nonlinear constraint equations using Newton's method rather than the nonlinear Gauss-Seidel method. Moreover, ILVES-PC solves the necessary linear systems using a specialized linear solver that exploits the structure of the protein. ILVES-PC can rapidly solve constraint equations as accurately as the hardware will allow. The run-time of ILVES-PC is proportional to the number of constraints. We have integrated ILVES-PC into GROMACS and simulated proteins of different sizes. Compared with SHAKE, we have achieved speedups of up to 4.9× in single-threaded executions and up to 76× in shared-memory multi-threaded executions. Moreover, ILVES-PC is more accurate than P-LINCS algorithm. Our work is a proof-of-concept of the utility of software designed specifically for the simulation of polymers.

The Giant Monopole and Dipole Resonances within the Constrained Molecular Dynamics (CoMD) Framework

The Constrained Molecular Dynamics (CoMD) model is used to describe the collective motion of various nuclear systems. A CoMD-inspired phenomenology for the GDR is developed. In addition, the dependence of the GDR upon the effective interaction parameters is studied. Furthermore, both the monopole and dipole main and soft modes of 68Ni are reliably reproduced. We conclude that a hard EoS with K=308 MeV increases the GDR energy, without altering the GMR energy. Thus, this EoS gives rather consistent results in both the monopole and dipole giant resonances.

Investigating the possible existence of hyper-heavy nuclei in a neutron-star environment

The synthesis of hyper-heavy elements is investigated under conditions simulating neutron star environment. The Constrained Molecular Dynamics (CoMD) approach is used to simulate low energy collisions of extremely n-rich nuclei. A new type of the fusion barrier due to a "neutron wind" is observed when the effect of neutron star environment (screening of Coulomb interaction) is introduced implicitly. When introducing also a background of surrounding nuclei, the nuclear fusion becomes possible down to temperatures of 10$^{8}$ K and synthesis of extremely heavy and n-rich nuclei appears feasible. A possible existence of hyper-heavy nuclei in a neutron star environment could provide a mechanism of extra coherent neutrino scattering or an additional mechanism, resulting in X-ray burst or a gravitational wave signal and, thus, becoming another crucial process adding new information to the suggested models on neutron star evolution.

Thermalization of nuclear matter in heavy-ion collisions at Fermi energies

We analyze the time evolution of the kinetic properties of nuclear matter produced in heavy-ion collisions at Fermi energies. The collision system is simulated using Constrained Molecular Dynamics (CoMD) transport calculations whose output is the isospin, position, and momentum of the nucleons. Focusing on central 35A⋅MeV 40Ca+40Ca collisions we utilize this information to extract localized momentum distributions in volume elements of 8fm3 and time steps of 5fm/c. We then parameterize the single-particle momentum distributions with thermally motivated fit functions in the local rest frame of each cell. While the transverse-momentum distributions are well reproduced by thermal ones, the longitudinal ones carry a marked imprint of the initial nuclear motion which we capture by introducing a centroid motion into our fit functions. In particular, we find that Fermi distributions yield significantly better fits than Boltzmann ones, a consequence of the Pauli blocking implemented in CoMD. From the fits we extract the time dependence of the thermodynamic and collective properties of the excited nuclear medium. We find that the transverse temperature gradually rises to about 6MeV, which is accompanied by a dissipation of the initial centroid motion of the incoming nuclei which vanishes at about 100fm/c after initial impact. We are therefore able to track the transition of beam energy into random kinetic energy for nucleons, suggesting a three-dimensional equilibration of energy in the late stages of the collision.

Effect of Charged Surfaces, High Concentrated and Localized High Concentrated Electrolytes on Lithium Ion Solvation Complex Evolution Near the Electrode Surface

The lack of understanding of the Li-ion deposition mechanism remains due to its complexity. However, the migration of the lithium ion from the electrolyte towards the Li metal deposition is heavily influenced by the evolution of its solvation complex along this pathway from bulk to surface. In this work, the desolvation and deposition dynamics near charged surfaces are studied to understand the electron transfer to the solvation complex that leads to the ion deposition We focus on the structural changes of the primary solvation shell as well as the energy barriers that the cation encounters through the diffusion pathway that is defined by the electrolyte and all its components including anions, additives, and solvents. Constrained molecular dynamics and cluster calculations with Density Functional Theory (DFT) are used to study the solvation complex changes and the thermodynamics behind those changes. The relevance of this work lies in identifying the role of each of the electrolyte component to facilitate or impede the ion mobility and the structural changes of the ions solvation complex towards a better understanding of ionic diffusion and deposition.

Amphiphile conformation impacts aggregate morphology and solution structure across multiple lengthscales

Although the self-assembly of amphiphiles is well-studied in aqueous solutions, much less is understood about the fundamental driving forces and structure property relationships in non-polar media. In recent work [Journal of Physical Chemistry B, 2020, 124, 10822.] the authors have studied a series of malonomide-based amphiphiles that are relevant to liquid-liquid extraction. That work demonstrated that aggregation is largely driven by local dipole-dipole interactions between molecules. Here, we build upon this observation to develop a more detailed understanding of how the balance of dipole-dipole interactions (controlled by conformation) and molecular architecture influences the morphology of the aggregates across lengthscales. Using constrained molecular dynamics about key degrees of freedom, we demonstrate that the conformation of N,N-dimethyl,N,N-dioctylhexylethoxy malonamide (DMDOHEMA) and N,N-dimethyl,N,N-dibutyltetradecyl malonamide (DMDBTDMA) has a significant impact upon self-association - where appropriate conformational sampling is essential. To quantify the aggregate morphology, several graph theoretic and persistent homology based properties are determined. The former examines the patterns of intermolecular interactions within clusters, while the latter examines the 3-dimensional spatial distribution. Based upon these analyses, we find that the morphology of aggregates, particularly at higher malonamide concentration, depends on a balance of dipole alignment and alkyl tail sterics. Dipole alignment encourages linear patterns of the intermolecular interactions within aggregates, while the alkyl tail steric interactions between the malonamide result in noticeably less linear aggregates for DMDOHEMA than DMDBTDMA. This is reflected in the spatial distribution, where more holes or voids exist between extractants within the DMDOHEMA that distribute within the solution in more of a “swiss cheese” arrangement as opposed to the more filamentous distribution of DMBDTDMA. This study links conformation and molecular structure to the morphology of amphiphile assemblies, and serves as a basis for ongoing study of multicomponent amphiphile solutions with polar and other solutes, and how these impact aggregation phenomena.

Constrained Molecular Dynamic Simulation of the Potential Mean Force of Lithium Bromide Ion Pairs in Acetonitrile

Molecular dynamic simulations of Li+, and Br− ions in acetonitrile were carried out. The simulated structural properties were compared to experimental data. The solvent potentials of Li+-Br−, Li+-Li+, and Br−-Br− were evaluated using constrained molecular dynamics (CMD) simulations, to determine the solvent contribution to the total force acting on the solute and estimate the liquid arrangements according to the potential of mean force (PMF) values. The PMF of friction kernels and the passage across the Li+-Br− barrier was studied using the Grote–Hynes theory. The union-separation development happens in the polarization confining system.

Coulomb field correction due to virtual e+e− production in heavy ion collisions

The correction to the Coulomb energy due to virtual production of e+e− pairs, which is on the order of one percent of the Coulomb energy at nuclear scales is discussed. The effects of including a pair-production term in the semi-empirical mass formula and the correction to the Coulomb barrier for a handful of nuclear collisions using the Bass and Coulomb potentials are studied. With an eye toward future work using Constrained Molecular Dynamics (CoMD) model, we also calculate the correction to the Coulomb energy and force between protons after folding with a Gaussian spatial distribution.

Heavy ion induced charge exchange reactions: A preliminary study within the Constrained Molecular Dynamics model

We present a preliminary study of charge exchange heavy ion induced reactions based on the constrained molecular dynamics (CoMD) model. The purpose is to test the capability of the model in predicting the occurrence of single charge exchange (SCE) and double charge exchange (DCE) exit channels for three different entrance channels at the same laboratory incident energy. The nuclear reaction dynamics and nuclear interaction within the CoMD approach are the only ingredients that have given, at this stage, promising results for SCE and DCE cross section calculations. The obtained results suggest an upgrade and possible future employment of the model for studies relating to the production of exotic nuclei through charge exchange reactions or DCE reactions and their connection with neutrinoless double beta decay.

Two-Channel Model for Electron Transfer in a Dye-Catalyst-Dye Supramolecular Complex for Photocatalytic Water Splitting.

To improve the performance of dye‐sensitized photoelectrochemical cell (DS‐PEC) devices for splitting water, the tailoring of the photocatalytic four‐photon water oxidation half‐reaction represents a principle challenge of fundamental significance. In this study, a Ru‐based water oxidation catalyst (WOC) covalently bound to two 2,6‐diethoxy‐1,4,5,8‐diimide‐naphthalene (NDI) dye functionalities provides comparable driving forces and channels for electron transfer. Constrained ab initio molecular dynamics simulations are performed to investigate the photocatalytic cycle of this two‐channel model for photocatalytic water splitting. The introduction of a second light‐harvesting dye in the Ru‐based dye‐WOC‐dye supramolecular complex enables two separate parallel electron‐transfer channels, leading to a five‐step catalytic cycle with three intermediates and two doubly oxidized states. The total spin S=1 is conserved during the catalytic process and the system with opposite spin on the oxidized NDI proceeds from the Ru=O intermediate to the final Ru‐O2 intermediate with a triplet molecular 3O2 ligand that is eventually released into the environment. The in‐depth insight into the proposed photocatalytic cycle of the two‐channel model provides a strategy for the development of novel high‐efficiency supramolecular complexes for DS‐PEC devices with buildup and conservation of spin multiplicity along the reaction coordinate as a design principle.

A Constrained Molecular Dynamics (CoMD) study of nuclear near-ground-state properties

The Constrained Molecular Dynamics (CoMD) model is used to describe the properties of nuclear systems near the ground state. A procedure for global optimization of the initial configurations of the nuclei is developed. In addition, the neutron skins of various nuclear systems are calculated. Finally, the GDR and GMR spectra of medium-mass nuclear systems are studied. The effect of the model parameters to the spectra is explored. We conclude that an increased compressibility of K = 308 MeV results in increased GDR energy and decreased skin, while the total energy and the GMR energy remain almost unaltered.

Two-stage description of charge distributions in 56Fe+p spallation reactions at 0.3-1.5 GeV/A

Spallation reactions are studied in the framework of a two-stage model. The intra-nuclear cascade (INC) stage of the reaction is simulated with the code ISABEL and a Constrained Molecular Dynamics (CoMD) code. The de-excitation of the highly excited pre-fragments is described with the multi-sequential binary decay code MECO, based on a generalized Weisskopf-Ewing formalism. Emission of nucleons, gamma rays and IMFs in their ground, excited bound and unbound states is considered. Calculated cross sections are compared with the experimental charge distributions of 56Fe + p spallation reaction products studied at GSI with the fragment separator FRS in the energy-range 0.3 − 1.5 GeV/A. At all bombarding energies, a good description of the experimental data is obtained with ISABEL-MECO using a global set of parameters. A preliminary calculation with CoMD-MECO at 1GeV/A is discussed.

Correlation between time and angular alignment in molecular dynamics simulations of heavy ion collisions

Neutron-proton equilibration is a process which has been used to study the density dependence of the symmetry energy term in the nuclear equation-of-state. This study utilizes constrained molecular dynamics (CoMD) simulations of $^{70}\mathrm{Zn}\phantom{\rule{4pt}{0ex}}+\phantom{\rule{4pt}{0ex}}^{70}\mathrm{Zn}$ with collision energies of 35 and 45 MeV/nucleon. An algorithm is used which searches through CoMD events and identifies the PLF* after it separates from the target and determines its lifetime, $\mathrm{\ensuremath{\Delta}}t$. It also determines the fragments that the PLF* breaks apart into and determines their angular alignment. This technique gives an opportunity to explore how the average alignment of dynamically produced fragments, ${\ensuremath{\langle}\ensuremath{\alpha}\ensuremath{\rangle}}_{\mathrm{dyn}}$, evolves with PLF* lifetime. An approximately linear relationship was determined with $d{\ensuremath{\langle}\ensuremath{\alpha}\ensuremath{\rangle}}_{\mathrm{dyn}}/d\mathrm{\ensuremath{\Delta}}t=0.98\ifmmode\pm\else\textpm\fi{}0.08\phantom{\rule{0.16em}{0ex}}\mathrm{rad}/\mathrm{zs}$ and $1.06\ifmmode\pm\else\textpm\fi{}0.09\phantom{\rule{0.16em}{0ex}}\mathrm{rad}/\mathrm{zs}$ for the 35 and $45\phantom{\rule{0.16em}{0ex}}\mathrm{MeV}/\mathrm{nucleon}$, respectively, indicating a correlation with magnitude consistent with classically determined values which were used for prior experimental studies.

Computational method for highly constrained molecular dynamics of rigid bodies: Coarse-grained simulation of auxetic two-dimensional protein crystals.

The increasing number of protein-based metamaterials demands reliable and efficient theoretical and computational methods to study the physicochemical properties they may display. In this regard, we develop a simulation strategy based on Molecular Dynamics (MD) that addresses the geometric degrees of freedom of an auxetic two-dimensional protein crystal. This model consists of a network of impenetrable rigid squares linked through massless rigid rods. Our MD methodology extends the well-known protocols SHAKE and RATTLE to include highly non-linear holonomic and non-holonomic constraints, with an emphasis on collision detection and response between anisotropic rigid bodies. The presented method enables the simulation of long-time dynamics with reasonably large time steps. The data extracted from the simulations allow the characterization of the dynamical correlations featured by the protein subunits, which show a persistent motional interdependence across the array. On the other hand, non-holonomic constraints (collisions between subunits) increase the number of inhomogeneous deformations of the network, thus driving it away from an isotropic response. Our work provides the first long-timescale simulation of the dynamics of protein crystals and offers insights into promising mechanical properties afforded by these materials.

Robosample: A rigid-body molecular simulation program based on robot mechanics

BackgroundCompared with all-atom molecular dynamics (MD), constrained MD methods allow for larger time steps, potentially reducing computational cost. For this reason, there has been continued interest in improving constrained MD algorithms to increase configuration space sampling in molecular simulations. MethodsHere, we introduce Robosample, a software package that implements high-performance constrained dynamics algorithms, originally developed for robotics, and applies them to simulations of biomolecular systems. As in the gMolmodel package developed by Spiridon and Minh in 2017, Robosample uses Constrained Dynamics Hamiltonian Monte Carlo (CDHMC) as a Gibbs sampling move - a type of Monte Carlo move where a subset of coordinates is allowed to change. In addition to the previously described Cartesian and torsional dynamics moves, Robosample implements spherical and cylindrical joints that can be distributed along the molecule by the user. ResultsIn alanine dipeptide simulations, the free energy surface is recovered by mixing fully flexible with torsional, cylindrical, or spherical dynamics moves. Ramachandran dynamics, where only the two key torsions are mobile, accelerate the slowest transition by an order of magnitude. We also show that simulations of a complex glycan cover significantly larger regions of the configuration space when mixed with constrained dynamics. Major conclusionsRobosample is a tool of choice for efficient conformational sampling of large biomolecules. General significanceRobosample is intended as a reliable and user-friendly simulation package for fast biomolecular sampling that does not require extensive expertise in mechanical engineering or in the statistical mechanics of reduced coordinates.