The determination of quartic force fields for use in vibrational second-order perturbation (VPT2) calculations, currently available in numerous electronic structure packages, becomes very expensive as the size of the molecule increases, especially if high-level coupled-cluster theory is used. Machine-learned potentials (MLPs) for large molecules and clusters offer a viable alternative to obtaining the quartic force field (QFF). Here, we report Fortran and Python software to determine the QFF and perform VPT2 calculations of energies from the MLPs. We describe this software and then apply it to H2O and protonated oxalate as the test cases. The Fortran software is applied to 21-atom aspirin using a fast MLP reported by us. Despite the fact that there are 32,509 unique cubic force constants for aspirin, the computer time to calculate them using this MLP is trivial, i.e., around 1 min. The new software provides an efficient way to calculate quantum anharmonic energies, using the established VPT2 methodology, for machine learned potentials of large molecules.
Recent research on photovoltaic/thermal (PV/T) collectors has focused on two key strategies to enhance performance: geometric modifications of the thermal flow channel (such as fins, baffles, and ribbed structures) and the integration of advanced materials like phase change materials (PCMs) and porous media to improve heat transfer and overall efficiency. In this direction, this numerical study investigates the performance enhancement of a photovoltaic/thermal (PV/T) solar collector through the integration of a solid layer along the lower wall of the airflow channel, coupled with a porous medium. The solid layer is introduced to accelerate the airflow and intensify convective heat transfer, thereby improving the thermal management of the photovoltaic cells. To solve the governing transport equations, an in-house computational code was developed in the Fortran programming language based on the finite volume method, coupled with the SIMPLER algorithm. The effects of solid-layer thickness and length, porous-layer thickness, and Darcy number are systematically investigated under a constant Reynolds number (Re = 500) and a uniform heat flux of 1000 W/m². The obtained results show that increasing the solid-layer thickness significantly enhances airflow acceleration and leads to a pronounced reduction in PV cell temperature of up to 33 °C. Extending the solid-layer length further improves the convective cooling process and increases both electrical and thermal efficiencies. When combined with a sufficiently permeable porous layer, additional performance gains are achieved, particularly at high Darcy numbers (Da = 10- 1). Compared to a conventional PV/T collector, the optimized configuration demonstrates enhancements of up to 60% in thermal efficiency and 28% in electrical efficiency.
Solid rocket propellants based on hydroxyl-terminated polybutadiene (HTPB), in which HTPB acts as the polymeric binder and fuel matrix, are widely used in aerospace propulsion. During storage, transport, and service, these composite energetic materials are exposed to sustained mechanical loads as well as environmental variations, which may induce time-dependent inelastic deformation. Such creep deformation can alter the grain geometry, affect combustion stability, and reduce the structural reliability of rocket motors. In this work, room-temperature tensile creep tests were conducted on an HTPB-based solid propellant under different stress levels. Several viscoelastic and power-law constitutive models were compared, and a composite time-hardening creep model was established to describe the experimental strain-time response. The model was further implemented in Abaqus through a Fortran user subroutine for finite element simulation. The results provide a useful basis for creep deformation assessment, formulation optimization, and structural reliability analysis of HTPB-based propellants.
Computational vibrational spectroscopy beyond the harmonic approximation relies on the molecular potential and ideally dipole and possibly higher moments of charge distributions. In the past decade, there has been a paradigm shift in generating highly accurate Machine-Learned potentials (MLPs). These are precise fits to thousands of electronic energies, using modern methods of regression. With such MLPs, it is possible to combine these with a variety of post-harmonic quantum methods ranging from perturbation theory to full variational calculations. After a short review of these methods, we focus on vibrational self-consistent field and configuration interaction (VSCF + VCI) calculations, as implemented in the code MULTIMODE. Two applications of this software to complex parts of the infrared spectra of formic acid dimer and the protonated oxalate anion are presented. Two new interfaces to MULTIMODE are then given. One is a Python-based GUI to enable user-friendly input to MULTIMODE. The second interface, PyFort, which is written in Fortran, uses MLPs written in Python in MULTIMODE via a C wrapper. Demonstrations of this are given for a PhysNet potential of Meuwly and co-workers for protonated oxalate anion (C2O4H-) and for the "universal" force field MACE-OFF of Csányi and co-workers. MULTIMODE VSCF + VCI vibrational energies of C2O4H- using the PhysNet MLP agree well with those using a permutationally invariant potential, trained on the datasets used to train the PhysNet MLP. A test of the MACE-OFF interface is done for H2CO. The PyFort software for both these examples is provided in the supplementary material.
In this paper, we present the recent advances in the computation of the Dirac-Kohn-Sham (DKS) method of the BERTHA code. We show here that the simple underlined structure of the FORTRAN code also favors efficient porting of the code to GPUs, leading to a particularly efficient hybrid CPU/GPU implementation (OpenMP/OpenACC), where the most computationally intensive part for DKS matrix evaluation (three-center two-electron integrals evaluated via the McMurchie-Davidson scheme) is efficiently offloaded to the GPU via compiler directives based on the OpenACC programming model. This scheme in combination with the use of a linear algebra library optimized for GPUs (cuBLAS, cuSOLVER) significantly accelerates the DKS calculations. In addition, the low-level integral kernel developed here at FORTRAN level was used to port our real-time DKS (RT-TDDKS) implementation based on Python (PyBERTHART) for the utilization of the GPU. The results obtained on the new Tier-0 EuroHPC supercomputer (LEONARDO) of the CINECA Supercomputing Centre with a single NVIDIA A100 card are very satisfactory. We achieve a speedup up to 30 for Au16 in a single-point DKS energy calculation and up to 10 for the Au8 systems in an RT-TDDKS calculation, compared to our OpenMP (i.e., CPU only) parallel implementation (with 32 cores). The approach presented here is very general and, to our knowledge, represents the first port of a Python API to GPUs based on a FORTRAN kernel for the evaluation of two-electron integrals. The implementation is currently limited to the use of a single GPU accelerator, but future paths to an actual exascale implementation are discussed.
A novel fourth-order, L-stable, exponential time differencing Runge-Kutta type scheme is developed to solve non-linear systems of reaction-diffusion equations, particularly with non-smooth data. The new scheme, ETDRK4RDP, is constructed by approximating the matrix exponentials in the ETDRK4 scheme with a fourth-order, L-acceptable, non-Padé rational function having real and distinct poles. The L-acceptable rational approximation ensures efficient damping of spurious oscillations arising from non-smooth initial and/or boundary conditions. The real and distinct poles of the rational function eliminate the presence of complex arithmetic in the solution of linear systems and thus make the scheme more attractive for parallelization in low-level languages like Fortran and C++. We verify empirically that the new ETDRK4RDP scheme is fourth-order accurate for several reaction diffusion systems with Dirichlet and Neumann boundary conditions and show it to be more efficient than competing exponential time differencing schemes when implemented in parallel, with up to five times speed-up in CPU time.
As a novel type of high-energy-density, environmentally friendly, and low-sensitivity energetic materials (EMs), cyclo-pentazolate salts are being extensively studied. However, their detonation mechanism remains unclear. This study developed a neural network potential (NNP) to simulate the shock-induced detonation process of NH3OH+N5-, a representative salt of the pentazolate anion (N5-). The well-trained NNP exhibits high precision comparable to DFT, as well as high efficiency. The NNP-based large-scale molecular dynamics (MD) simulations for NH3OH+N5- produced an ideal C-J detonation velocity of 9.4 km/s, which is in agreement with the value estimated by the Cheetah 7.0 program (9.93 km/s). The simulation demonstrates that the proton transfer from NH3OH+ to N5- is the initial reaction, while the primary decomposition pathway of N5- is a ring-opening reaction, or the bimolecular reactions with its initial decomposition intermediate azide anion N3- resulting in the formation of N8 ring. Quantum chemical calculations show that these pathways possess low activation barriers. The influence of nuclear quantum effects on shock-induced chemical reactions was also studied, which shows that nuclear quantum corrections not only improve the accuracy of predicted ideal detonation velocity but also improve temperature in simulations, which results in the different reaction mechanism of shock-induced detonation reaction of NH3OH+N5-, facilitating the ring-opening reaction of N5- ring and preventing its reaction with N3. This study enhances the understanding of the detonation mechanism of cyclo-pentazolate salts. In this work, NNP potential was trained by the DeePMD-kit package and homemade FORTRAN code. The density functional theory (DFT) calculation of structural energies and atomic forces, as well as ab initio molecular dynamics (AIMD), was conducted using the Vienna Ab initio Simulation Package (VASP) software. The PAW method and the GGA-PBE functional were adopted. The shock wave response MD simulations were conducted by LAMMPS with the multiscale shock technique (MSST) and QB-MSST methods. Quantum chemical calculations were carried out at the M06-2X/TZVP level using the Gaussian 09 program.
The oscillatory motion of graphene oxide-Go nanoparticles in water lubricating Maxwell nanofluid flow for heat and mass transfer enhancement in steady and fluctuating regime is important significance of this study. The aim of this work indicates the temperature distribution, nanoparticle concentration rate and velocity field around stretching radiating-cylinder in drilling systems. The nonlinear radiating energy, entropy generation, mixed convection, buoyancy ratio and oscillatory effects are assumed for heat and mass performance. The partial differential based mathematical expressions are developed to estimate the values of current analysis. The oscillatory stokes conditions, complex variables, and primitive transformations are applied to develop steady and fluctuating results. The computational outputs are secured using very efficient methods like finite difference and Gaussian elimination. The graphical outputs are displayed through TecPlot-360 and FORTRAN. The fluid velocity field, energy, and concentration outputs are executed with the help of various physical factors. The steady friction and steady thermal rate are depicted and are used in transient algorithm to depict the fluctuating friction rate and oscillatory thermal-mass transport. The magnitude of fluid velocity, fluid temperature and fluid concentration enhances as Maxwell parameter is enhanced. The variation in fluid velocity amplitude and fluid temperature increases as radiating energy is enhanced. For high parametric range of Maxwell number, the steady skin friction and mass transfer increases but heat transfer decreases. At smaller Eckert number, the oscillations in transient skin friction, transient heat and mass transfer are enhanced. At higher Maxwell parameter, buoyancy and Schmidt number, the large amplitude in heat and mass oscillations is observed.
We present a software to calculate phase-resolved resonant vibrational sum-frequency generation (vSFG) susceptibility χ(2)(ω) of water and hydroxyls at planar interfaces, e.g., air/water or solid/liquid or (bio)membrane/liquid interfaces of aqueous solutions. The released code (i) reads several formats of molecular trajectories, both from ab initio (AIMD) and classical MD (CMD), (ii) calculates instantaneous surfaces to allow flexible interfaces, (iii) is written in FORTRAN, parallelized by OpenMP and optimized for memory usage, (iv) allows processing of systems up of tens of thousand atoms and for unlimited simulation time, and (v) includes many tunable processing parameters. The code and its documentation are available via GitHub. Flexible models of water and surface hydroxyl (if evaluated) (CMD or AIMD) must be used. The derivatives of the polarizability tensors and dipole moments with the change of O-H distance must be calculated externally by ab initio methods and provided as input data. We present the impact of various parameters of the MD simulations (simulation length, nonbonded interaction cutoff, size of the system, and thermostat relaxation time) as well as of the processing code (filter relaxation, cutoff of cross-terms) and provide representative results for air/water, charged quartz (101)/aqueous solution, and neutral α-alumina (0001)/aqueous solution interfaces. Further extensions are planned to distinguish signals from specific O-H or C-H bonds of interfacial molecules.
Incorporation of machine learning (ML) techniques into atomic-scale modeling has proven to be an extremely effective strategy to improve the accuracy and reduce the computational cost of simulations. It also entails conceptual and practical challenges, as it involves combining very different mathematical foundations as well as software ecosystems that are very well developed in their own right but do not share many commonalities. To address these issues and facilitate the adoption of ML in atomistic simulations, we introduce two dedicated software libraries. The first one, metatensor, provides multi-platform and multi-language storage and manipulation of arrays with many potentially sparse indices, designed from the ground up for atomistic ML applications. By combining the actual values with metadata that describes their nature and that facilitates the handling of geometric information and gradients with respect to the atomic positions, metatensor provides a common framework to enable data sharing between ML software-typically written in Python-and established atomistic modeling tools-typically written in Fortran, C, or C++. The second library, metatomic, provides an interface to store an atomistic ML model and metadata about this model in a portable way, facilitating the implementation, training, and distribution of models, and their use across different simulation packages. We showcase a growing ecosystem of tools, including low-level libraries, training utilities, and interfaces with existing software packages, that demonstrate the effectiveness of metatensor and metatomic in bridging the gap between traditional simulation software and modern ML frameworks.
The Community Multiscale Air Quality (CMAQ) model simulates atmospheric phenomena, including advection, diffusion, gas-phase chemistry, aerosol physics and chemistry, and cloud processes. Gas-phase chemistry is often a major computational bottleneck due to its representation as large systems of coupled nonlinear stiff differential equations. We leverage the parallel computational performance of graphics processing unit (GPU) hardware to accelerate the numerical integration of these systems in CMAQ's CHEM module. Our implementation, dubbed CMAQ-CUDA, in reference to its use in the Compute Unified Device Architecture (CUDA) general purpose GPU (GPGPU) computing solution, migrates CMAQ's Rosenbrock solver from Fortran to CUDA Fortran. CMAQ-CUDA accelerates the Rosenbrock solver such that simulations using the chemical mechanisms RACM2, CB6R5, and SAPRC07 require only 51%, 50%, or 35% as much time, respectively, as CMAQv5.4 to complete a chemistry time step. Our results demonstrate that CMAQ is amenable to GPU acceleration and highlight a novel Rosenbrock solver implementation for reducing the computational burden imposed by the CHEM module.
Pipelines in soft soil are prone to deformation and failure under traffic loads. Therefore, it is highly important to accurately characterize the dynamic response of pipelines under traffic loads and reasonably evaluate their operational status. First, for rigid pipelines, a Dload subroutine is written in the FORTRAN language to accurately characterize traffic loads, and a 3D numerical analysis model of the rigid pipe‒soil system is established using ABAQUS software to simulate the dynamic response of rigid pipelines in soft soil under traffic loads. The simulation results are validated against data from field tests. Second, an improved particle swarm optimization (PSO) algorithm is introduced to optimize the weights and thresholds of a backpropagation neural network (BPNN). An improved PSO-BPNN method for predicting the dynamic response of pipelines is proposed, and the accuracy and applicability of the method are verified. Finally, using the prediction model, a reliability analysis is conducted on the dynamic response of rigid pipelines in soft soil under traffic loads. The results show that compared with smaller-diameter pipelines, larger-diameter pipelines exhibit lower axial stress and vertical displacement, with more concentrated distributions. During pipeline construction, larger-diameter pipelines should be chosen whenever possible to reduce the adverse impact of factors such as traffic loads on the dynamic response of pipelines. These research results provide a new theoretical basis and technical support for enhancing the reliability of rigid pipelines in soft soil and conducting in-depth safety assessments of pipelines under traffic loads.
To evaluate the effect of suction direction on zonular fibers' mechanical behavior using finite element modeling. Laboratory study. 3D finite element model experimental study. A 3D finite element model was developed, including nucleus covered by cortex and capsule. Similar to cataract surgery, a circular rapture was considered at the top of the capsule to apply suction pressure. Finally, zonular fibers were modeled as a continuum body using 3D solid elements (C3D8R). A custom FORTRAN subroutine was implemented to enforce tension-only behavior, mimicking the physiological characteristics of the zonular fibers, which are resistant to tensile but not compressive loads. This method allows for a realistic simulation of zonular mechanics during cortical aspiration. The suction processes in tangential and normal angles were simulated. The resulting relative displacements between the cortex and capsule, as a criterion of cortex separation, vs the resulting maximum zonule displacements, were recorded in each model. The cortex-capsule relative displacement vs maximum zonule displacement indicated a diagram slope of 0.09 for tangential applying pressure and 0.02 in the case of applying normal pressure. The results illustrated that zonular fibers were less tensile in a specific magnitude of cortex-capsule relative displacement under tangential applying pressure, which means that eye zonular fibers have a lower risk of failure until the separation of the lens cortex.
PyPWDFT is a Python software designed for performing plane-wave density functional theory (DFT) calculations. It can perform large-scale DFT calculations using only a single process on a single node, including local density functional for 10,000 atoms and nonlocal hybrid functional for 4096 atoms. Our benchmark test results demonstrate that PyPWDFT achieves performance comparable to that of Fortran/C++ codes, despite being developed in a native Python environment. In addition, it requires only NumPy, SciPy, and CuPy, enabling CPU-GPU heterogeneous computing, achieving a two-order-of-magnitude speedup compared to single-threaded CPU execution. Due to its excellent cross-platform compatibility, medium-scale DFT calculations can be performed through a graphical user interface on personal computers and Windows systems using consumer-grade GPUs, such as the NVIDIA GeForce RTX 4090. The computational efficiency is comparable to that of professional-grade GPUs such as the NVIDIA V100. The efficient performance, scalability to handle large-scale systems, high numerical accuracy, and different interfaces for molecular dynamics collectively underscore the considerable potential of PyPWDFT to develop into versatile DFT software.
Understanding the interactions of phenylboronic acid with its surrounding water molecules is essential for several applications in solvated systems. In the present work, we investigated the microhydration of the phenylboronic acid (PBA) and calculated its hydration free energy using the cluster continuum solvation model. Microhydration of PBA has not been investigated previously in the literature. It requires the structures of PBA to be surrounded by n explicit water molecules (PBA-W n ). The results show that the B(OH) 2 unit of phenylboronic acid forms clusters with water molecules that are similar to those of neutral water clusters. The QTAIM analysis shows that the structures of phenylboronic acid-water clusters are stabilized by strong OH ⋯ O and weak CH ⋯ O hydrogen bonds. In addition to QTAIM analysis, NBO analysis was also performed on the most stable configurations to better understand the delocalization of electron density from donor to proper acceptor within the compound. In addition, we found that the most stable structures dominate the population of the clusters for temperatures from 20 to 400 K. Finally, using the structures of the microhydrated phenylboronic acid, we estimated the free energy of hydration and the enthalpy of hydration of PBA. At room temperature, the phenylboronic acid's free energy and enthalpy of hydration are respectively evaluated to be - 72.1 kcal/mol and - 85.5k cal/mol. Assessment of temperature effects on the free energy and the enthalpy of hydration shows that the enthalpy is temperature-independent, while the free energy increases linearly with temperature. Initial configurations of PBA-W n have been generated using classical molecular dynamics and subsequently optimized using the level of theory, ω B97X-D/def2-TZVP. Optimizations, frequency calculations, and NBO analysis are performed using the Gaussian 16 suite of programs. On the most stable configurations, we have performed the quantum theory of atoms in molecules (QTAIM) analysis to get insights into the hydrogen bond network of PBA-W n . QTAIM is performed using AIMAll software. Thermodynamic properties as a function of temperature are evaluated using a homemade FORTRAN code-named TEMPO.
Nonadiabatic molecular dynamics (NAMD) simulations are crucial for revealing the underlying mechanisms of photochemical and photophysical processes. Typical NAMD simulation software packages rely on on-the-fly ab initio electronic structure and nonadiabatic coupling calculations, and thus become challenging when dealing with large complex systems. We here introduce a new Simulation Package for non-Adiabatic Dynamics in Extended systems (SPADE), which is designed to address the limitations of traditional surface hopping methods in dealing with these problems. By design, SPADE enables the users to define arbitrary quasi-diabatic Hamiltonians through parametrized functions and incorporates a variety of algorithms (e.g., global flux hopping probabilities, complex crossing and decoherence corrections), which can realize efficient and reliable NAMD simulations without using nonadiabatic couplings at all. All the employed methods and expressions for diabatic Hamiltonian matrix elements can be flexibly set through the input files. SPADE is mainly written in Fortran based on a modular design and has a great capacity for further implementation of new methods. SPADE can be used to simulate both model and atomistic systems as long as proper Hamiltonians are provided. As demonstrations, a series of representative models are studied to show the main features and capabilities.
Controlling water hammer pressure is essential, necessitates a transient surge analysis to identify critical pressure points along a pipeline system. A pressurized air vessel is a pressure control device used to control both positive and negative pressure fluctuations. This study investigates three key parameters that affect the sizing of the pressurized air vessel: orifice diameter (the throttling aperture), the vessel diameter, and water volume fraction ratio. A mathematical model, developed using the FORTRAN programming language and based on the unsteady one-dimensional momentum and continuity equations, determines the optimal sizing of these parameters. These equations are solved using the method of characteristics, and the pressurized air vessel is mathematically modelled as a quasi-one-dimensional flow system. An experimental test rig, equipped with a rapid closing solenoid valve and pressure sensors, is used to validate the mathematical model results. Both the experimental and numerical results demonstrate the effectiveness of the pressurized air vessel to dampen water hammer pressure. The findings indicate that the throttling action has a significant effect on the required size of the pressurized air vessel. This study presents a novel approach that provides quantitative insights into key parameters that affect the performance of the pressurized air vessel by using the combined modelling and experimental validation. The orifice diameter is the most influential parameter on the water hammer head, vessel air head, and water level inside the vessel.
Architectural forms in practical engineering projects exhibit significant diversity and complexity, often exceeding the representational capacity of traditional curves or even B-spline curves. This limitation creates substantial challenges in establishing accurate initial models.Given a set of key points on a B-spline and the introduction of a reference plane, the internal and external boundaries were determined using a B-spline curved surface formed by the key points and a parametric B-spline curve on the reference plane, where the reference plane was divided using Delaunay triangulation. Then, an initial structure with a complicated boundary was obtained. To address the challenges of computational efficiency and mesh distortion, the strain energy sensitivity to the key points was derived using the relationship between the strain energy and the key points. A new method was established to minimize the strain energy while improving the computation efficiency. The optimization model, sensitivity analysis, optimum algorithm, and mechanics analysis program were all implemented in FORTRAN language, with two free-form continuous shell structure surfaces to demonstrate the correctness and effectiveness of the method. The failure pattern, displacement contours and load-displacement curves were studied, for cases with optimized typical positive and negative Gaussian curvature. Finally, the ultimate load before and after optimization, as well as the effect of concrete strength, shell thickness and reinforcement location on the ultimate load, were studied using the finite element analysis software ABAQUS.
The development of ordered structures upon immediately quenching from 473 K to different crystallization temperatures (Tc = 300 K, 320 K, 340 K, and 360 K) of three n-alkanes (n-eicosane (C20H42), n-tetracontane (C40H82) and n-octacontane (C80H162) to accommodate the formation of bilayer, monolayer and the folded-chain configurations, respectively), was quantitatively analyzed by evaluating chain/bond order parameters, interaction energies, conformational statistics, local dynamics, structural pair correlation function, and X-ray scattering profiles. The formation of ordered structures seems to be best at Tc = 340 K. Generally, systems that tend to form the monolayer and folded-chain configurations exhibit the most and the least ordered structures, respectively. Chains in monolayer structures tend to have a larger fraction of trans state, higher anisotropy of bond orientation, slower monomer dynamics, and more densely packed structures than chains in bilayer structures, while longer alkanes with chain-folded structures exhibit the least characteristics of these properties. The nucleation temperature (Tn) can be affected by their structural configuration, with the highest Tn for the monolayer structure. Nevertheless, the melting temperature (Tm) tends to depend solely on the molecular weights of these alkanes, not on their structural configuration. Monte Carlo (MC) simulation of the coarse-grained (CG) models of three n-alkanes (one CG bead equivalent to an ethylene unit) on the second nearest neighbor diamond (2nnd) lattice with the same periodic box dimension of 5 nm in each size. The energetics of CG models were composed of the Rotational Isomeric State (RIS) model and the Lennard-Jones (LJ) potential energies to represent their intra- and intermolecular interactions, respectively. Structure developments within 100 million Monte Carlo steps (MCS) trajectories were monitored, and data analysis was based on snapshots collected at intervals of 10,000 MCS. All simulations and data analysis were performed using in-house FORTRAN codes with the g77 compiler.
The field of chemical reaction dynamics has evolved considerably since its inception, driven by advances in computational power and theoretical methodologies. While ab initio molecular dynamics (AIMD) simulations offer high accuracy by computing forces directly from electronic structure theory, its high computational cost limits its applicability to small systems and short time scales. Machine learning (ML) potentials and exascale computing offer new ways of developing potential energy surfaces (PESs) for chemical reaction dynamics simulations of longer times scales and larger systems. In this work, we introduce VENUSpy, a Python-based reimplementation and extension of the classical VENUS code, designed to interface with existing ML-based potentials and the exascale-ready quantum chemistry package NWChemEx. While NWChemEx is in development, VENUSpy demonstrates interfacing with its top-level classes and objects for use in reaction dynamics simulations. VENUSpy preserves the large library of initial sampling and classical trajectory propagation of the original, Fortran-based VENUS, while expanding its versatility through integration with modern Python tools. This modular framework enables ML/ab initio hybrid dynamics simulations, which offers distinct advantages to AIMD simulations and MLMD simulations, as well as fosters rapid development of new methodologies for studying complex reactive systems.