From f43487d0c76b9c41780015bea450c58074b8ea0f Mon Sep 17 00:00:00 2001 From: "tobias.hainer" Date: Fri, 25 Sep 2026 10:19:35 +0200 Subject: [PATCH 1/3] Add thermodynamic integration ensemble for NEP --- src/integrate/ensemble.cuh | 3 +- src/integrate/ensemble_ti_nep.cu | 334 ++++++++++++++++++++++++++++++ src/integrate/ensemble_ti_nep.cuh | 87 ++++++++ src/integrate/integrate.cu | 4 + 4 files changed, 427 insertions(+), 1 deletion(-) create mode 100644 src/integrate/ensemble_ti_nep.cu create mode 100644 src/integrate/ensemble_ti_nep.cuh diff --git a/src/integrate/ensemble.cuh b/src/integrate/ensemble.cuh index c5027ada9..e2de203bd 100644 --- a/src/integrate/ensemble.cuh +++ b/src/integrate/ensemble.cuh @@ -57,7 +57,8 @@ enum class EnsembleType TI_RS, TI_AS, WALL_HARMONIC, - TI_LIQUID + TI_LIQUID, + TI_NEP }; inline bool is_standard_nvt(const EnsembleType type) diff --git a/src/integrate/ensemble_ti_nep.cu b/src/integrate/ensemble_ti_nep.cu new file mode 100644 index 000000000..b389407d4 --- /dev/null +++ b/src/integrate/ensemble_ti_nep.cu @@ -0,0 +1,334 @@ +/* + Copyright 2017 Zheyong Fan and GPUMD development team + This file is part of GPUMD. + GPUMD is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + GPUMD is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with GPUMD. If not, see . +*/ + +/*----------------------------------------------------------------------------80 +Nonequilibrium thermodynamic integration between two NEP potentials. +The first potential in run.in (lambda = 0) drives the regular force +calculation and the second potential (lambda = 1) is evaluated here. +------------------------------------------------------------------------------*/ + +#include "ensemble_ti_nep.cuh" +#include "utilities/gpu_macro.cuh" + +namespace +{ +static __global__ void gpu_add_nep2_force( + int number_of_atoms, + double lambda, + double* fx, + double* fy, + double* fz, + const double* fx_nep2, + const double* fy_nep2, + const double* fz_nep2) +{ + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < number_of_atoms) { + fx[i] = (1 - lambda) * fx[i] + lambda * fx_nep2[i]; + fy[i] = (1 - lambda) * fy[i] + lambda * fy_nep2[i]; + fz[i] = (1 - lambda) * fz[i] + lambda * fz_nep2[i]; + } +} + +static __global__ void initialize_properties( + int N, double* g_fx, double* g_fy, double* g_fz, double* g_pe, double* g_virial) +{ + int n1 = blockIdx.x * blockDim.x + threadIdx.x; + if (n1 < N) { + g_fx[n1] = 0.0; + g_fy[n1] = 0.0; + g_fz[n1] = 0.0; + g_pe[n1] = 0.0; + for (int d = 0; d < 9; ++d) { + g_virial[n1 + d * N] = 0.0; + } + } +} + +static __global__ void gpu_get_pe_sum(const int N, double* pe) +{ + //<<<1, 1024>>> + int tid = threadIdx.x; + int batch, n; + int number_of_batches = (N - 1) / 1024 + 1; + __shared__ double s_data[1024]; + s_data[tid] = 0.0; + + for (batch = 0; batch < number_of_batches; batch++) { + n = tid + batch * 1024; + if (n < N) + s_data[tid] += pe[n]; + } + __syncthreads(); + for (int offset = blockDim.x >> 1; offset > 0; offset >>= 1) { + if (tid < offset) + s_data[tid] += s_data[tid + offset]; + __syncthreads(); + } + if (tid == 0) + pe[0] = s_data[0]; +} + +} // namespace + +Ensemble_TI_Nep::Ensemble_TI_Nep(const std::vector& tokens) +{ + const int num_params = tokens.size(); + temperature_coupling = 100; + int i = 2; + while (i < num_params) { + if (tokens[i] == "tswitch") { + if (i + 1 >= num_params) + PRINT_INPUT_ERROR("Missing value for tswitch keyword."); + auto_switch = false; + if (!is_valid_int(tokens[i + 1], &t_switch)) + PRINT_INPUT_ERROR("Wrong inputs for t_switch keyword."); + i += 2; + } else if (tokens[i] == "tequil") { + if (i + 1 >= num_params) + PRINT_INPUT_ERROR("Missing value for tequil keyword."); + auto_switch = false; + if (!is_valid_int(tokens[i + 1], &t_equil)) + PRINT_INPUT_ERROR("Wrong inputs for t_equil keyword."); + i += 2; + } else if (tokens[i] == "temp") { + if (i + 1 >= num_params) + PRINT_INPUT_ERROR("Missing value for temp keyword."); + if (!is_valid_real(tokens[i + 1], &temperature)) + PRINT_INPUT_ERROR("Wrong inputs for temp keyword."); + i += 2; + } else if (tokens[i] == "tperiod") { + if (i + 1 >= num_params) + PRINT_INPUT_ERROR("Missing value for tperiod keyword."); + if (!is_valid_real(tokens[i + 1], &temperature_coupling)) + PRINT_INPUT_ERROR("Wrong inputs for t_period keyword."); + i += 2; + } else { + PRINT_INPUT_ERROR("Unknown keyword."); + } + } + if ((t_switch < 0 && t_equil >= 0) || (t_switch >= 0 && t_equil < 0)) { + PRINT_INPUT_ERROR( + "Error: Please specify either both t_switch and t_equil, or neither (to let the program " + "auto-determine)"); + } + printf( + "Thermostat: target temperature is %f k, t_period is %f timesteps.\n", + temperature, + temperature_coupling); + type = EnsembleType::NVT_LAN; + c1 = exp(-0.5 / temperature_coupling); + c2 = sqrt((1 - c1 * c1) * K_B * temperature); +} + +void Ensemble_TI_Nep::init( + const int number_of_steps, const Atom& atom, const GPU_Vector& thermo) +{ + if (auto_switch) { + t_switch = (int)(number_of_steps * 0.4); + t_equil = (int)(number_of_steps * 0.1); + } else + printf("The number of steps should be set to %d!\n", 2 * (t_equil + t_switch)); + printf( + "Nonequilibrium thermodynamic integration: t_switch is %d timestep, t_equil is %d timesteps.\n", + t_switch, + t_equil); + output_file = my_fopen("ti_nep.csv", "w"); + fprintf(output_file, "lambda,dlambda,pe,pe_nep2\n"); + int N = atom.number_of_atoms; + + curand_states.resize(N); + int grid_size = (N - 1) / 128 + 1; + initialize_curand_states<<>>(curand_states.data(), N, rand()); + GPU_CHECK_KERNEL + + thermo_cpu.resize(thermo.size()); + potential_nep2.resize(N); + force_nep2.resize(3 * N); + virial_nep2.resize(9 * N); +} + +void Ensemble_TI_Nep::find_thermo( + const Box& box, + const std::vector& group, + const Atom& atom, + GPU_Vector& thermo) +{ + Ensemble::find_thermo( + box.get_volume(), + group, + atom.mass, + atom.potential_per_atom, + atom.velocity_per_atom, + atom.virial_per_atom, + thermo); + thermo.copy_to_host(thermo_cpu.data()); + pe = thermo_cpu[1]; +} + +double Ensemble_TI_Nep::get_pe_nep2(const int number_of_atoms) +{ + // reduces in place; potential_nep2 is re-initialized every step + double temp; + gpu_get_pe_sum<<<1, 1024>>>(number_of_atoms, potential_nep2.data()); + GPU_CHECK_KERNEL + potential_nep2.copy_to_host(&temp, 1); + return temp; +} + +Ensemble_TI_Nep::~Ensemble_TI_Nep(void) { close_output_file(false); } + +void Ensemble_TI_Nep::finalize_run(const Atom& atom, const Box& box) +{ + FILE* yaml_file = my_fopen("ti_nep.yaml", "w"); + fprintf(yaml_file, "F_diff: %f\n", F_diff); + fprintf(yaml_file, "T: %f\n", temperature); + + close_output_file(true); + fclose(yaml_file); + + printf("\n"); + printf("-----------------------------------------------------------------------\n"); + printf("Helmholtz free energy difference (NEP1 - NEP2): %f eV/atom.\n", F_diff); + printf("This value is stored in ti_nep.yaml.\n"); + printf("-----------------------------------------------------------------------\n"); +} + +void Ensemble_TI_Nep::close_output_file(const bool print_message) +{ + if (output_file != nullptr) { + if (print_message) { + printf("Closing ti_nep output file...\n"); + } + fclose(output_file); + output_file = nullptr; + } +} + +void Ensemble_TI_Nep::add_nep2_force(Atom& atom) +{ + int N = atom.number_of_atoms; + gpu_add_nep2_force<<<(N - 1) / 128 + 1, 128>>>( + N, + lambda, + atom.force_per_atom.data(), + atom.force_per_atom.data() + N, + atom.force_per_atom.data() + 2 * N, + force_nep2.data(), + force_nep2.data() + N, + force_nep2.data() + 2 * N); + GPU_CHECK_KERNEL +} + +void Ensemble_TI_Nep::initialize_before_first_step( + const double, + const int number_of_steps, + const std::vector&, + Box&, + Atom& atom, + GPU_Vector& thermo) +{ + init(number_of_steps, atom, thermo); +} + +void Ensemble_TI_Nep::find_lambda( + const int step, + const Box& box, + const std::vector& group, + Atom& atom, + GPU_Vector& thermo) +{ + find_thermo(box, group, atom, thermo); + bool need_output = false; + + const int t = step - t_equil; + pe_nep2 = get_pe_nep2(atom.number_of_atoms); + + const double r_switch = 1.0 / t_switch; + + if ((t >= 0) && (t <= t_switch)) { + lambda = switch_func(t * r_switch); + dlambda = dswitch_func(t * r_switch); + need_output = true; + } else if ((t >= t_equil + t_switch) && (t <= (t_equil + 2 * t_switch))) { + lambda = switch_func(1.0 - (t - t_switch - t_equil) * r_switch); + dlambda = -dswitch_func(1.0 - (t - t_switch - t_equil) * r_switch); + need_output = true; + } + + if (need_output) { + fprintf( + output_file, + "%e,%e,%e,%e\n", + lambda, + dlambda, + pe / atom.number_of_atoms, + pe_nep2 / atom.number_of_atoms); + F_diff += 0.5 * (pe - pe_nep2) * fabs(dlambda) / atom.number_of_atoms; + } +} + +void Ensemble_TI_Nep::compute2( + const double time_step, + const int step, + const int number_of_steps, + const std::vector& group, + Box& box, + Atom& atoms, + GPU_Vector& thermo, + Force& force) +{ + if (!potentials_checked) { + if (force.get_number_of_potentials() < 2) { + PRINT_INPUT_ERROR("ensemble ti_nep requires two potentials (NEP1 and NEP2) in run.in."); + } + if (force.get_potential(1).nep_model_type == 3) { + PRINT_INPUT_ERROR("ensemble ti_nep does not support a temperature-dependent second NEP."); + } + potentials_checked = true; + } + + const int N = atoms.number_of_atoms; + initialize_properties<<<(N - 1) / 128 + 1, 128>>>( + N, + force_nep2.data(), + force_nep2.data() + N, + force_nep2.data() + N * 2, + potential_nep2.data(), + virial_nep2.data()); + GPU_CHECK_KERNEL + + force.get_potential(1).compute( + box, atoms.type, atoms.position_per_atom, potential_nep2, force_nep2, virial_nep2); + + find_lambda(step, box, group, atoms, thermo); + add_nep2_force(atoms); + + Ensemble_LAN::compute2(time_step, step, number_of_steps, group, box, atoms, thermo, force); +} + +double Ensemble_TI_Nep::switch_func(double t) +{ + double t2 = t * t; + double t5 = t2 * t2 * t; + return ((70.0 * t2 * t2 - 315.0 * t2 * t + 540.0 * t2 - 420.0 * t + 126.0) * t5); +} + +double Ensemble_TI_Nep::dswitch_func(double t) +{ + double t2 = t * t; + double t4 = t2 * t2; + return ((630 * t2 * t2 - 2520 * t2 * t + 3780 * t2 - 2520 * t + 630) * t4) / t_switch; +} diff --git a/src/integrate/ensemble_ti_nep.cuh b/src/integrate/ensemble_ti_nep.cuh new file mode 100644 index 000000000..7450192a3 --- /dev/null +++ b/src/integrate/ensemble_ti_nep.cuh @@ -0,0 +1,87 @@ +/* + Copyright 2017 Zheyong Fan and GPUMD development team + This file is part of GPUMD. + GPUMD is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + GPUMD is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with GPUMD. If not, see . +*/ + +#pragma once +#include "ensemble_lan.cuh" +#include "force/force.cuh" +#include "langevin_utilities.cuh" +#include "model/box.cuh" +#include "utilities/common.cuh" +#include "utilities/error.cuh" +#include "utilities/read_file.cuh" +#include +#include +#include + +class Ensemble_TI_Nep : public Ensemble_LAN +{ +public: + Ensemble_TI_Nep(const std::vector& tokens); + ~Ensemble_TI_Nep(void) override; + + void finalize_run(const Atom& atom, const Box& box) override; + + void initialize_before_first_step( + const double time_step, + const int number_of_steps, + const std::vector& group, + Box& box, + Atom& atom, + GPU_Vector& thermo) override; + + void compute2( + const double time_step, + const int step, + const int number_of_steps, + const std::vector& group, + Box& box, + Atom& atoms, + GPU_Vector& thermo, + Force& force) override; + + void find_thermo( + const Box& box, + const std::vector& group, + const Atom& atom, + GPU_Vector& thermo); + double get_pe_nep2(const int number_of_atoms); + void add_nep2_force(Atom& atom); + void init(const int number_of_steps, const Atom& atom, const GPU_Vector& thermo); + void find_lambda( + const int step, + const Box& box, + const std::vector& group, + Atom& atom, + GPU_Vector& thermo); + double switch_func(double t); + double dswitch_func(double t); + +protected: + FILE* output_file = nullptr; + double lambda = 0, dlambda = 0; + int t_equil = -1, t_switch = -1; + double pe, pe_nep2; + double F_diff = 0; // Helmholtz free energy difference F(NEP1) - F(NEP2) (eV/atom) + bool auto_switch = true; + bool potentials_checked = false; + // per-atom properties computed with the second potential + GPU_Vector potential_nep2; + GPU_Vector force_nep2; + GPU_Vector virial_nep2; + std::vector thermo_cpu; + +private: + void close_output_file(const bool print_message); +}; diff --git a/src/integrate/integrate.cu b/src/integrate/integrate.cu index 633b0c067..6fdbf5f4b 100644 --- a/src/integrate/integrate.cu +++ b/src/integrate/integrate.cu @@ -34,6 +34,7 @@ The driver class for the various integrators. #include "ensemble_ti.cuh" #include "ensemble_ti_as.cuh" #include "ensemble_ti_liquid.cuh" +#include "ensemble_ti_nep.cuh" #include "ensemble_ti_rs.cuh" #include "ensemble_ti_spring.cuh" #include "ensemble_ttm.cuh" @@ -419,6 +420,9 @@ void Integrate::parse_ensemble( } else if (tokens[1] == "ti_liquid") { type = EnsembleType::TI_LIQUID; ensemble_ = std::make_unique(tokens); + } else if (tokens[1] == "ti_nep") { + type = EnsembleType::TI_NEP; + ensemble_ = std::make_unique(tokens); } else { PRINT_INPUT_ERROR("Invalid ensemble type."); } From 3b008b0ccde976fc92fd7c3595c3efa6a4679487 Mon Sep 17 00:00:00 2001 From: "tobias.hainer" Date: Tue, 29 Sep 2026 16:57:26 +0200 Subject: [PATCH 2/3] Changed the precision of reported free energy --- src/integrate/ensemble_ti_nep.cu | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/integrate/ensemble_ti_nep.cu b/src/integrate/ensemble_ti_nep.cu index b389407d4..8a1f5c232 100644 --- a/src/integrate/ensemble_ti_nep.cu +++ b/src/integrate/ensemble_ti_nep.cu @@ -193,7 +193,7 @@ Ensemble_TI_Nep::~Ensemble_TI_Nep(void) { close_output_file(false); } void Ensemble_TI_Nep::finalize_run(const Atom& atom, const Box& box) { FILE* yaml_file = my_fopen("ti_nep.yaml", "w"); - fprintf(yaml_file, "F_diff: %f\n", F_diff); + fprintf(yaml_file, "F_diff: %.16e\n", F_diff); fprintf(yaml_file, "T: %f\n", temperature); close_output_file(true); @@ -201,7 +201,7 @@ void Ensemble_TI_Nep::finalize_run(const Atom& atom, const Box& box) printf("\n"); printf("-----------------------------------------------------------------------\n"); - printf("Helmholtz free energy difference (NEP1 - NEP2): %f eV/atom.\n", F_diff); + printf("Helmholtz free energy difference (NEP1 - NEP2): %.16e eV/atom.\n", F_diff); printf("This value is stored in ti_nep.yaml.\n"); printf("-----------------------------------------------------------------------\n"); } @@ -271,7 +271,7 @@ void Ensemble_TI_Nep::find_lambda( if (need_output) { fprintf( output_file, - "%e,%e,%e,%e\n", + "%.16e,%.16e,%.16e,%.16e\n", lambda, dlambda, pe / atom.number_of_atoms, From 28ff74677c01c2a72a25b74ac5de60d6a1830fe2 Mon Sep 17 00:00:00 2001 From: "tobias.hainer" Date: Mon, 5 Oct 2026 09:32:33 +0200 Subject: [PATCH 3/3] Add documentation for ensemble ti_nep Co-Authored-By: Claude Opus 5.5 --- doc/bibliography.rst | 12 ++++ .../input_parameters/ensemble_ti_nep.rst | 65 +++++++++++++++++++ doc/gpumd/input_parameters/index.rst | 1 + 3 files changed, 78 insertions(+) create mode 100644 doc/gpumd/input_parameters/ensemble_ti_nep.rst diff --git a/doc/bibliography.rst b/doc/bibliography.rst index 74132432a..c4b5af2dc 100644 --- a/doc/bibliography.rst +++ b/doc/bibliography.rst @@ -10,6 +10,12 @@ Bibliography | Applied Mathematics Research eXpress **2014**, 332–352 (2014) | DOI: `10.1093/amrx/abu005 `_ +.. [Berger2026] + | Ethan Berger, Tobias Hainer, Tobias Möslinger, Paul Erhart, and Julia Wiktor + | *Stability of Polarons and Oxygen Vacancies in BiVO₄ Revealed by Machine-Learning-Driven Simulations* + | PRX Energy **5**, 033004 (2026) + | DOI: `10.1103/xwxv-6s55 `_ + .. [Bernetti2020] | Mattia Bernetti and Giovanni Bussi | *Pressure control using stochastic cell rescaling* @@ -160,6 +166,12 @@ Bibliography | Computational Materials Science **175**, 109584 (2020) | DOI: `10.1016/j.commatsci.2020.109584 `_ +.. [Hainer2025] + | Tobias Hainer, Ethan Berger, Esmée Berger, Olof Hildeberg, Paul Erhart, and Julia Wiktor + | *Thermal Stabilization of Defect Charge States and Finite-Temperature Charge Transition Levels* + | arXiv:2512.15463 (2025) + | DOI: `10.48550/arXiv.2512.15463 `_ + .. [Hoover1996] | William G. Hoover and Brad Lee Holian | *Kinetic moments method for the canonical ensemble distribution* diff --git a/doc/gpumd/input_parameters/ensemble_ti_nep.rst b/doc/gpumd/input_parameters/ensemble_ti_nep.rst new file mode 100644 index 000000000..be770cc03 --- /dev/null +++ b/doc/gpumd/input_parameters/ensemble_ti_nep.rst @@ -0,0 +1,65 @@ +.. _kw_ensemble_ti_nep: + +:attr:`ensemble` (TI_NEP) +========================= + +This keyword is used to set up a nonequilibrium thermodynamic integration integrator between two :term:`NEP` models. Please check [Freitas2016]_ for more details on the nonequilibrium switching scheme. + +The system is driven by the Hamiltonian + +.. math:: + + H(\lambda) = (1-\lambda) U_1 + \lambda U_2, + +where :math:`U_1` and :math:`U_2` are the potential energies from the first (NEP1) and the second (NEP2) :ref:`potential ` specified in the :attr:`run.in` file. +Lambda is switched from 0 to 1 (forward) and back from 1 to 0 (backward), and the Helmholtz free energy difference :math:`F_{\rm NEP1} - F_{\rm NEP2}` is obtained from the average of the forward and backward work. +This is an NVT ensemble with a Langevin thermostat, so the free energy difference is evaluated at the volume of the simulation box. + +Syntax +------ + +The parameters can be specified as follows:: + + ensemble ti_nep temp tperiod tequil tswitch + +- :attr:``: The temperature of the simulation. +- :attr:``: This parameter is optional, and defaults to ``100``. It determines the period of the thermostat in units of the timestep. It determines how strongly the system is coupled to the thermostat. +- :attr:``: The number timesteps to equilibrate the system at :math:`\lambda=0` and at :math:`\lambda=1`. +- :attr:``: The number timesteps to vary lambda from 0 to 1 (and from 1 to 0). + +If you do not specify :attr:`` and :attr:``, they will be automatically set in a 1:4 ratio. +If you specify them, you must specify both, and the number of steps in the :ref:`run keyword ` should be 2 × (:attr:`` + :attr:``). + +Exactly two NEP potentials are needed in the :attr:`run.in` file. +The first one is NEP1 (:math:`\lambda=0`) and the second one is NEP2 (:math:`\lambda=1`). +Only the forces are mixed during the switching; the potential energies of NEP1 and NEP2 are kept separately, because both are needed to compute the free energy difference. +The potential energy and pressure written by :ref:`dump_thermo ` are therefore those of NEP1 during the whole run. + +Example +------- + +.. code-block:: rst + + potential nep1.txt + potential nep2.txt + velocity 300 + ensemble ti_nep temp 300 tequil 10000 tswitch 40000 + run 100000 + +This command equilibrates for 10000 timesteps with NEP1, switches lambda from 0 to 1 in 40000 timesteps, equilibrates for 10000 timesteps with NEP2, and switches lambda back from 1 to 0 in 40000 timesteps. + +Example: switching between defect states +---------------------------------------- + +This ensemble has been used to compute the free energy difference between two charge states of a point defect [Hainer2025]_ and between polaronic states [Berger2026]_. +One NEP model is trained for one charge state and then relabeled for the other charge state; the original and the relabeled model are given as NEP1 and NEP2. +The resulting :attr:`F_diff` is the raw Helmholtz free energy difference between the two charge states, which goes into the expression for the formation free energy of the defect. + +Output file +----------- + +This command will produce a csv file :attr:`ti_nep.csv`. The columns are lambda, dlambda, potential energy of NEP1 and potential energy of NEP2 (eV/atom). The file is only written during the switching. + +This command will also produce a yaml file :attr:`ti_nep.yaml`, which contains the Helmholtz free energy difference :attr:`F_diff` (:math:`F_{\rm NEP1} - F_{\rm NEP2}`, eV/atom) and the temperature. + +Both files are overwritten by a subsequent :attr:`ti_nep` run in the same directory, so use separate directories for separate calculations. diff --git a/doc/gpumd/input_parameters/index.rst b/doc/gpumd/input_parameters/index.rst index 64a395267..f9b68e302 100644 --- a/doc/gpumd/input_parameters/index.rst +++ b/doc/gpumd/input_parameters/index.rst @@ -32,6 +32,7 @@ Below you can find a listing of keywords for the ``run.in`` input file. ensemble_ti_as ensemble_ti_rs ensemble_ti + ensemble_ti_nep ensemble_shock_nemd ensemble_msst ensemble_nphug