Engineering Papers⌕ Search

SEARCH · Engineering Papers

Results for “equilibrium modeling”

Search indexed NASA NTRS and DOE OSTI research on propulsion, heat transfer, battery materials and energy systems. Follow report and document links to the original sources.

Quote a phrase for an exact phrase match. Source license links do not imply unrestricted reuse.

At least 541 records · Page 30

Phase-field modeling of radiation-induced segregation for multicomponent alloys

Structural alloys under irradiation are known to undergo radiation-induced solute redistribution (RIS) at grain boundaries, leading to detrimental effects such as intergranular corrosion and stress-assisted cracking. To better understand the phenomenon, improved models of RIS applicable to concentrated, multicomponent alloys, and mesoscale microstructures are needed. In this talk, we present a novel grand-potential-based phase-field model to account for the complete set of multicomponent kinetic and thermodynamic couplings between atoms and point defects in the Onsager transport equations. We demonstrate multiscale modeling capability by deriving the Onsager coefficient matrix from atomistic-based Kinetic Monte Carlo simulations. Model predictions and validations of RIS and the effect of defect production, grain boundary sink strength and density are demonstrated for a model FCC FeCrNi system. We also demonstrate the novel capability to describe RIS in the presence of equilibrium segregation described using a density-based CALPHAD thermodynamic framework. This modeling approach overcomes certain limitations in conventional RIS models and provides a step closer towards multiscale modeling.

73 NUCLEAR PHYSICS AND RADIATION PHYSICS↗

Sea-surface temperature pattern effects have slowed global warming and biased warming-based constraints on climate sensitivity

The observed rate of global warming since the 1970s has been proposed as a strong constraint on equilibrium climate sensitivity (ECS) and transient climate response (TCR)—key metrics of the global climate response to greenhouse-gas forcing. Using CMIP5/6 models, we show that the inter-model relationship between warming and these climate sensitivity metrics (the basis for the constraint) arises from a similarity in transient and equilibrium warming patterns within the models, producing an effective climate sensitivity (EffCS) governing recent warming that is comparable to the value of ECS governing long-term warming under CO 2 forcing. However, CMIP5/6 historical simulations do not reproduce observed warming patterns. When driven by observed patterns, even high ECS models produce low EffCS values consistent with the observed global warming rate. The inability of CMIP5/6 models to reproduce observed warming patterns thus results in a bias in the modeled relationship between recent global warming and climate sensitivity. Correcting for this bias means that observed warming is consistent with wide ranges of ECS and TCR extending to higher values than previously recognized. These findings are corroborated by energy balance model simulations and coupled model (CESM1-CAM5) simulations that better replicate observed patterns via tropospheric wind nudging or Antarctic meltwater fluxes. Because CMIP5/6 models fail to simulate observed warming patterns, proposed warming-based constraints on ECS, TCR, and projected global warming are biased low. The results reinforce recent findings that the unique pattern of observed warming has slowed global-mean warming over recent decades and that how the pattern will evolve in the future represents a major source of uncertainty in climate projections.

58 GEOSCIENCES↗

Analytical population dynamics underlying harmful algal blooms triggered by prey avoidance

A modified version of the NPZD ecosystem model is used to analytically examine the effects of predation avoidance, a possible mechanism for triggering harmful algal blooms (HAB). To resolve HAB development caused by predation avoidance, an additional phytoplankton functional group is considered, one that has slower nutrient uptake and better predation avoidance characteristics than the non-harmful phytoplankton group used in traditional NPZD models. Because the two phytoplankton groups (one non-harmful and one HAB) compete for only one resource within the same system, steady state (equilibrium) conditions cannot occur without the presence of zooplankton; only the non-harmful phytoplankton group, which defeats the HAB group in the resource competition, can survive in the equilibrium. The presence of sufficient zooplankton effectively acts to replenish the nutrient pool by consuming the non-harmful phytoplankton. When this occurs, two equilibrium states are found: one with both phytoplankton groups coexisting, and one that only includes the HAB group. Further, the condition required for equilibrium is that the total nitrogen within the system should be larger than a threshold determined by model coefficients. The threshold and feasibility of the equilibrium are sensitive to the relative HAB predation avoidance coefficient. If the coefficient is larger than the ratio of net growth rates between the HAB and non-harmful phytoplankton group, the threshold becomes infinite, and an equilibrium is not feasible. The time scale for the system to reach an equilibrium state that includes a HAB group is determined asymptotically. The dependence of a threshold condition as a controlling factor may explain the regime shift of dominant species causing HABs. The ecosystem model is fully implemented into the Regional Ocean Modeling System and applied to an idealized coastal embayment (with depths and geometry taken from San Francisco Bay) to show numerically the dominance of prey avoidance dynamics in a natural shallow water environment that includes advection and diffusion. The analytical results improve strategies for HAB modeling and provide guidance for setting model coefficients necessary to resolve a HAB event.

54 ENVIRONMENTAL SCIENCES↗

Fourier-like Thermal Relaxation of Nanoscale Explosive Hot Spots

Hot spots are local regions of high temperature that are widely considered to govern explosive initiation. Hot spot dynamics rests on a delicate balance between heat generation due to chemical reactions and heat loss through thermal conduction, making accurate determinations of the conductivity under extreme conditions a key component of predictive explosive models. We develop here an approach to directly determine the thermal transport properties of explosive hot spots with realistic initial structures through a combination of molecular dynamics (MD) and diffusive heat equation (HEq) modeling. Effective thermal conductivity values are determined by fitting HEq models to MD predictions of long timescale hot spot relaxation. The approach is applied to model hot spots in the molecular crystalline explosive 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) for a range of shock strengths and two limiting cases for impact orientation. Isotropic and anisotropic HEq models yield similar results, despite TATB exhibiting some of the largest and most anisotropic thermal conductivity values for explosive near normal conditions. The conductivity is found to be a strong function of density, which parametrically captures dependence on temperature, pressure, and material state. The associated root-mean-square errors of the fitted HEq models are approximately 5% of MD predicted final equilibrium temperatures. Here, the conductivity values determined here for TATB hot spots are considerably larger than those used in a prior hot spot criticality study, which may significantly impact predictions for critical hot spot sizes. The approach provides a convenient foundation for determining the effective thermal conductivity for hot spot problems in other explosives and directly yields information on reasonable approximations that might be taken in higher-level models for those materials.

36 MATERIALS SCIENCE↗

Characteristics of locational uncertainty marginal price for correlated uncertainties of variable renewable generation and demands

With the rapid increase of variable renewable energy sources in power systems, how to manage and price the uncertainty of renewable resources’ power outputs is becoming an urgent issue. Current market designs considering the uncertainties are mainly based on the probabilistic scenario set of demand and renewable energy resources power outputs. This consideration makes market designs vulnerable to three significant challenges when put into practice. First, the accurate probability distribution of renewable generation is hard to obtain in real-time. Second, it is challenging to clear the market timely with many scenarios to guarantee accuracy. Third, generation cost recovery cannot be guaranteed for some scenarios. To overcome these challenges, this paper proposes a locational uncertainty marginal price model to price the uncertainty explicitly based on a scenario-free stochastic market-clearing model. Instead of using the probabilistic scenario set, the uncertainty of renewable energy sources and loads is modeled with distributionally-robust chance constraints. The correlation of uncertainties can be endogenously modeled in both the market-clearing and the locational uncertainty marginal price formation. Furthermore, this paper proves that generation cost recovery, revenue adequacy, and partial market equilibrium can be achieved using the locational uncertainty marginal price model. Numerical results from both the small and large systems simulations validate that the generation cost recovery is maintained no matter the generation participates in uncertainty mitigation or not. The transmission congestion surplus is also allocated appropriately among loads, renewable energy sources, and financial transmission right owners.

29 ENERGY PLANNING, POLICY, AND ECONOMY↗

Phase Diagram and Structure Map of Binary Nanoparticle Superlattices from a Lennard-Jones Model

A first-principles prediction of the binary nanoparticle phase diagram assembled by solvent evaporation has eluded theoretical approaches. In this paper, we show that a binary system interacting through the Lennard-Jones (LJ) potential contains all experimental phases in which nanoparticles are effectively described as quasi hard spheres. We report a phase diagram consisting of 53 equilibrium phases, whose stability is quite insensitive to the microscopic details of the potentials, thus giving rise to some type of universality. Furthermore, we show that binary lattices may be understood as consisting of certain particle clusters, i.e., motifs, that provide a generalization of the four conventional Frank–Kasper polyhedral units. Our results show that metastable phases share the very same motifs as equilibrium phases. Finally, we discuss the connection with packing models, phase diagrams with repulsive potentials, and the prediction of likely experimental superlattices.

37 INORGANIC, ORGANIC, PHYSICAL, AND ANALYTICAL CH↗

Prediction of alkaline earth metal ion adsorption on goethite for various background electrolytes with the CD-MUSIC model

As water scarcity drives the use of more saline water sources, contaminant fate and transport models must capture the impact of high concentrations of alkaline earth metal ions (AEMs) and background electrolytes in these more complex waters. By utilizing macroscopic adsorption data from various electrolyte systems, a Charge Distribution – Multisite Complexation (CD-MUSIC) model, capable of incorporating electrolyte adsorption, was able to accurately simulate the adsorption behavior of alkaline earth metal ions onto goethite. The modeling effort was guided by previous spectroscopic and surface complexation modeling of alkaline earth metal adsorption and built on previous CD-MUSIC modeling that accounted for changes in crystal face contributions to the surface site density as a function of specific surface area. Further, the model was constrained to consider only two dominant surface complex species for each metal ion adsorption reaction. These two species were selected from 44 possible species through objective curve fitting of single-solute macroscopic adsorption data. While most of the alkaline earth metal surface complexes formed outer-sphere complexes at the goethite surface, an inner-sphere species was utilized for Mg 2+ . With the surface complex species and equilibrium constants obtained from this study, the calibrated model successfully predicted alkaline earth metal ion adsorption over a wide range of solution and surface conditions; the model predictions encompassed a wide range of pH (5–11), solute/solid ratio (1.37 × 10 -5 – 8.33 × 10 -4 mol -solute /g -solid ), ionic strengths (0.01 M – 0.7 M), and background electrolytes (Na + , Cs + , Rb + , Cl - , and NO 3 - ) using the same crystal face contribution methodology for site density, capacitance values, and surface acidity constants adopted for proton and cadmium adsorption in previous work (Han and Katz, 2019). Model simulations for a range of background water chemistries demonstrated the potential for Mg 2+ to reduce Cd 2+ adsorption to goethite in model seawater and oil- and gas-produced waters.

42 ENGINEERING↗

Design of novel refractory equiatomic multi-principal elemental alloys based on Mo-Nb-Ti system for Gen IV reactor applications

Excellent irradiation damage resistance demonstrated by multi-principal elemental alloys (MPEAs) has sparked significant interest among researchers, prompting exploration into their vast compositional space, to validate their suitability for nuclear applications. A combined approach of thermodynamic and empirical parameters calculations alongside CALPHAD (CALculation of PHAse Diagrams) for phase formation predictions enable high-throughput material selection for sophisticated applications like nuclear, overcoming laborious and time-consuming experiments. Key thermodynamic and empirical parameters for eight novel equiatomic MPEAs, based on seven low thermal neutron cross section refractory elements, for predicting phase formation were calculated, and equilibrium and non-equilibrium simulations in CALPHAD were employed to comprehensively model the systems. Pseudo binary phase diagram simulations showed that Zr, V or equiatomic CrV additions to the base MoNbTi alloy (MoNbTi-Zr, MoNbTi-V and MoNbTi-CrV alloys) favor the formation of isomorphous body-centered cubic (BCC) phase at high temperatures, while Cr, Al, equiatomic ZrV, or equiatomic CrAl additions (MoNbTi-Cr, MoNbTi-Al, MoNbTi-ZrV or MoNbTi-CrAl alloys) limit the solubility of them. Equilibrium CALPHAD simulations at 750 oC were consistent with XRD results on MoNbTi, MoNbTiZr and MoNbTiCr alloys, and partially for others. Notably, elemental segregation observed in the backscattered electron (BSE) scanning electron microscopy (SEM) images of the alloys was accurately simulated through non-equilibrium Scheil solidification calculations in CALPHAD, further verified by experiments. The precipitation of TiCr2 Laves phase in Cr containing MoNbTiCr and MoNbTiCrAl was accurately predicted while discrepancies were noted in MoNbTiCrV. The equilibrium simulations also provided insights into phase compositions at specific temperatures offering a pathway for tailoring the desired microstructure and properties of these systems. Empirical parameters calculations successfully predicted random solid solution in the base MoNbTi alloy, and with an exception in MoNbTiV and MoNbTiAl, predicted intermetallic precipitation in the rest, especially, Laves phase precipitation in Cr containing alloys.

36 - MATERIALS SCIENCE↗

On the relationship between the multi-region relaxed variational principle and resistive inner-layer theory

We show that the variational energy principle of the multi-region relaxed magnetohydrodynamic (MRxMHD) model can be used to predict finite-pressure linear tearing instabilities. In this model, the plasma volume is sliced into sub-volumes separated by 'ideal interfaces', and in each volume the magnetic field relaxes to a Taylor state, where the pressure gradient $\nabla p = 0$. The MRxMHD model is implemented in the Stepped-Pressure Equilibrium Code (SPEC) so that the equilibrium solution in each region is computed while preserving the force balance across the interfaces. As SPEC computes the Hessian matrix (a discretized stability matrix), the stability of an MRxMHD equilibrium can also be computed with SPEC. In this article, using SPEC, we investigate the effect of local pressure gradients and the $\nabla p = 0$ in the vicinity of the resonant surface of a tearing mode. For low-beta plasma, we have been able to illustrate a relationship between the resistive singular-layer theory (Coppi et al 1966 Nucl. Fusion6 101; Glasser et al 1975 Phys. Fluids18 875–88) and the MRxMHD model. Within the singular layer, the volume-averaged magnetic helicity and the flux-averaged toroidal flux are shown to be the invariants for the linear tearing modes in SPEC simulations. Our technique to compute MRxMHD stability is first tested numerically in a cylindrical tokamak and its application in toroidal geometry is demonstrated. We demonstrate an agreement between the stability boundary obtained with SPEC simulation and the resistive inner-layer theories.

70 PLASMA PHYSICS AND FUSION TECHNOLOGY↗

Toroidal modeling of energetic passing particle drift kinetic effects on tearing mode stability

Abstract Drift kinetic effects of the neutral beam injection induced passing energetic particles (EPs) on the linear stability of the n = 1 tearing mode (TM) (with the dominant poloidal harmonic of m = 2) are numerically investigated utilizing the MARS-K code (Liu et al 2008 Phys. Plasmas 15 112503), in a tokamak plasma with finite equilibrium pressure and anisotropic thermal transport. In the low plasma pressure regime, it is found that co- (counter-) passing EPs stabilize (destabilize) the TM, agreeing with previous studies. However, as the plasma pressure increases beyond a critical value, it is found that co-passing EPs also destabilize the mode. An in-depth analysis reveals that the net effect of co-passing EPs is a result of competition between the stabilizing contribution from the non-adiabatic drift kinetic terms and the destabilizing contribution associated with adiabatic terms, with the latter becoming more dominant at higher equilibrium pressure. Non-perturbative magnetohydrodynamic-kinetic hybrid modeling also finds that co- and counter-passing EPs modify the TM eigenfunction differently, with the counter-passing EPs enhancing the sideband harmonics. Furthermore, effects of the plasma resistivity and toroidal rotation, as well as that of the equilibrium distribution of EPs in the particle pitch angle space, are also investigated, showing asymmetric results on the TM stability between the co- and counter-passing EPs. The first order finite orbit width correction is found to be stabilizing with co-passing EPs and destabilizing with counter-passing particles. Finally, drift resonances between passing EPs and the TM induce finite frequency to the mode and generate finite net torques inside the plasma, due to the neoclassical toroidal viscosity and the Reynolds stress associated with 3D perturbations.

Physics↗

Heavy flavor quenching and flow: the roles of initial condition, pre-equilibrium evolution, and in-medium interaction

Within an advanced Langevin-hydrodynamics framework coupled to a hybrid fragmentation-coalescence hadronization model, we study heavy flavor quenching and flow in relativistic heavy-ion collisions. We investigate how the initial heavy quark spectrum, the in-medium energy loss and hadronization mechanisms of heavy quarks, the evolution profile of the pre-equilibrium stage, the medium flow, and the temperature dependence of heavy quark diffusion coefficients influence the suppression and elliptic flow of heavy mesons at the RHIC and the LHC. Our results show that the different modeling of initial conditions, pre-equilibrium evolution, and in-medium interactions can individually yield uncertainties of approximately 10-40% in D meson suppression and flow at a low transverse momentum. Furthermore, we also find that proper combinations of collisional versus radiative energy loss, coalescence versus fragmentation in hadronization, and the inclusion of medium flow are the most important factors for describing the suppression and elliptic flow of heavy mesons.

72 PHYSICS OF ELEMENTARY PARTICLES AND FIELDS↗

Generalizable machine learning potentials for quantum-accurate predictions of non-equilibrium behavior in 2D materials

Machine learning interatomic potentials (ML-IAPs) are emerging as transformative tools in materials modeling, promising quantum-level accuracy at a fraction of the computational cost. However, their ability to generalize beyond equilibrium configurations and to reliably capture defect- and temperature-driven behavior remains underexplored. Here, we develop and benchmark two state-of-the-art ML-IAPs, Spectral Neighbor Analysis Potential (SNAP) and Allegro, on a comprehensive dataset for monolayer MoSe₂. Using density functional theory (DFT) as the reference, we evaluate their performance in capturing stress–strain behavior, phase transition energetics, defect evolution, edge stability, and fracture toughness. Allegro, a deep equivariant neural network potential, surpasses both SNAP and the classical Tersoff potential in accuracy, efficiency, and transferability. Importantly, both ML potentials accurately reproduce experimental fracture measurements and ab initio predictions of inversion domain formation—phenomena well beyond their training sets. Our findings establish ML-IAPs as viable replacements for traditional force fields in the study of non-equilibrium mechanical phenomena, enabling large-scale, high-fidelity simulations in 2D materials and beyond. In conclusion, this work provides a broadly applicable framework for data-driven modeling of structural and functional transformations under extreme conditions.

2D materials↗

An Improved Model For Determining Salinity Recharge Time for Deep Borehole Disposal - 20405

In a previous publication (K. P. Travis, D. Burley and F. G. F. Gibb, WM2017 Conference, Phoenix Arizona paper no. 17480) we introduced a numerical model to estimate the time it would take for a pressure perturbation to subside following its creation from the construction of a deep geological borehole in water-saturated rock. The model was based on the notion of a point source of momentum - a solution of the time-dependent pressure-diffusion equation. Using the model, we estimated that it would take on the order of 10 k years for physical equilibrium to become re-established following the sinking of a 5 km borehole in granite. The significance of such a model is that it places an upper bound on the lifetime required of a borehole sealing system - a necessary input to a borehole post-closure safety case assessment. No engineered sealing system has ever been devised which is capable of retaining its sealing properties for 10 half-lives of long-lived radioisotopes in spent fuel or high-level waste. One of the key advantages of disposing of nuclear waste via Deep Borehole Disposal (DBD) is the natural sealing provided by density stratified groundwater. Upon sinking of a borehole, and subsequent filling with fresh water, brine or drilling mud, this natural barrier may be temporarily damaged. Over time, fresh brine from the far-field will flow towards or away from the hole (driven by a pressure gradient) and will eventually re-establish the original salinity gradient. It follows that the engineered seals need last only as long as the time required for this salinity gradient to reset itself. One of the limitations of our previous model was the use of a static boundary condition on the borehole wall. The model used a boundary pressure which varied quadratically with depth (arising from differences between the pressure of a column of fresh water in the borehole and that of a column of brine in the host rock). However as brine replaces fresh water in the borehole (driven by a pressure gradient), the boundary function must change with time. We now introduce an improved model which takes this time dependent boundary condition into account. The model also takes into account the time taken for a mixture of brine and fresh water to re-establish chemical equilibrium through the process of diffusion. The paper contains the mathematical details of our new iterative model as well as results showing the time taken to reach steady state, and concentration profiles for the components of the brine-filled borehole as a function of time together with a discussion on the implications of the results for developing a post-closure safety assessment for DBD. (authors)

12 MANAGEMENT OF RADIOACTIVE AND NON-RADIOACTIVE W↗

Cosmological consequences of first-order general-relativistic viscous fluid dynamics

In this work, we investigate the out-of-equilibrium dynamics of viscous fluids in a spatially flat Friedmann-LemaîtreRobertson-Walker cosmology using the most general causal and stable viscous energy-momentum tensor defined at first order in spacetime derivatives. In this new framework a pressureless viscous fluid having equilibrium energy density ρ can evolve to an asymptotic future solution in which the Hubble parameter approaches a constant while ρ → 0, even in the absence of a cosmological constant (i.e., Λ = 0). Thus, while viscous effects in this model drive an accelerated expansion of the universe, the equilibrium energy density itself vanishes, leaving behind only the acceleration. This behavior emerges as a consequence of causality in first-order theories of relativistic fluid dynamics and it is fully consistent with Einstein’s equations.

79 ASTRONOMY AND ASTROPHYSICS↗

Narrow operator models of stellarator equilibria in Fourier Zernike basis

Numerical computation of the ideal magnetohydrodynamic (MHD) equilibrium magnetic field is at the base of stellarator optimisation and provides the starting point for solving more sophisticated partial differential equations like transport or turbulence models. Conventional approaches solve for a single stationary point of the ideal MHD equations, which is fully defined by three invariants and the numerical scheme employed by the solver. We present the first numerical approach that can solve for a continuous distribution of equilibria with fixed boundary and rotational transform, varying only the pressure invariant. This approach minimises the force residual by optimising parameters of multilayer perceptrons that map from a scalar pressure multiplier to the Fourier Zernike basis as implemented in the modern stellarator equilibrium solver DESC.

fusion plasma↗

Noise reduction in X-ray photon correlation spectroscopy with convolutional neural networks encoder–decoder models

Abstract Like other experimental techniques, X-ray photon correlation spectroscopy is subject to various kinds of noise. Random and correlated fluctuations and heterogeneities can be present in a two-time correlation function and obscure the information about the intrinsic dynamics of a sample. Simultaneously addressing the disparate origins of noise in the experimental data is challenging. We propose a computational approach for improving the signal-to-noise ratio in two-time correlation functions that is based on convolutional neural network encoder–decoder (CNN-ED) models. Such models extract features from an image via convolutional layers, project them to a low dimensional space and then reconstruct a clean image from this reduced representation via transposed convolutional layers. Not only are ED models a general tool for random noise removal, but their application to low signal-to-noise data can enhance the data’s quantitative usage since they are able to learn the functional form of the signal. We demonstrate that the CNN-ED models trained on real-world experimental data help to effectively extract equilibrium dynamics’ parameters from two-time correlation functions, containing statistical noise and dynamic heterogeneities. Strategies for optimizing the models’ performance and their applicability limits are discussed.

36 MATERIALS SCIENCE↗

Modeling of a 43 GHz microplasma within a photonic crystal

The electron density (n e ) and argon metastable density (1 s 5 ) of a 43GHz microplasma are obtained using the zero-dimensional plasma kinetics solver (ZDPlasKin®) for an incident power of 160 and 1000mW and argon pressures of 10–600Torr (1.3 ×10 3 to 8.0 × 10 4 Pa). These simulations are compared with previously published experimental data. To determine the self-consistent electric field in the modeled plasma, the three-dimensional millimeter wave fields are computed as a function of electron density using ANSYS EM19.2, HFSS®. This electromagnetic field model is coupled to ZDPlasKin such that any increase in the simulated plasma density correctly attenuates the simulated electric field within the microplasma. The electron density is found to be sensitive to argon gas temperature, so a two-zone temperature model was needed to obtain agreement with experimental measurements. The temperature in the central core of the microplasma was determined by a previous experimental study. That temperature was used as an input to the model for the simulation of volume recombination losses. The outer regions of the microplasma are assumed to be in equilibrium with the walls (300K). This second temperature was used in the model to determine diffusion losses. The modeled electron and metastable densities are of the order 10 20 and 10 18 m –3 , respectively. This is in good agreement with those measured experimentally as long as the two-zone temperature model is used. In the absence of a hot gas core, the modeled three-body recombination rates are excessive and the simulation severely under-predicts the electron density and over-estimates the metastable density. We conclude that the millimeter wave microplasma has a hot core (2500K at 600 Torr) that rarifies the argon gas and effectively reduces three-body recombination. Furthermore, this allows one to achieve high electron density on the order of 10 20 m –3 with only 100mW of wave power.

70 PLASMA PHYSICS AND FUSION TECHNOLOGY↗

Phase diagram of the TIP4P/Ice water model by enhanced sampling simulations

In this work, we studied the phase diagram for the TIP4P/Ice water model using enhanced sampling molecular dynamics simulations. Our approach is based on the calculation of ice–liquid free energy differences from biased coexistence simulations that reversibly sample the melting and growth of layers of ice. We computed a total of 19 melting points for five different ice polymorphs, which are in excellent agreement with the melting lines obtained from the integration of the Clausius–Clapeyron equation. For proton-ordered and fully proton-disordered ice phases, the results are in very good agreement with previous calculations based on thermodynamic integration. For the partially proton-disordered ice III, we find a large increase in stability that is in line with previous observations using direct coexistence simulations for the TIP4P/2005 model. This issue highlights the robustness of the approach employed here for ice polymorphs with diverse degrees of proton disorder. Our approach is general and can be applied to the calculation of other complex phase diagrams.

71 CLASSICAL AND QUANTUM MECHANICS, GENERAL PHYSIC↗