From 1e6049294a25eac95a8a67618c234c88e0ba52e5 Mon Sep 17 00:00:00 2001 From: monst Date: Tue, 5 May 2026 15:25:39 +0900 Subject: [PATCH 1/9] feat: add HF cryogenic plasma etching model for SiO2 Implements HFCryoEtching, a physics-based process model for pure HF gas cryogenic plasma etching of SiO2. The model captures temperature-dependent surface chemistry via Frenkel-Arrhenius desorption and Arrhenius reaction kinetics, replacing chemical passivation layers with thermal suppression of lateral etching. Key additions: - psHFCryoParameters.hpp: parameter struct with Desorption (nu0, E_des), Reaction (A_r, E_a), SiO2, and Ion sub-structs; k_des(T) and k_r(T) helper methods compute Arrhenius rates at runtime. - psHFCryoEtching.hpp: two-particle model (HFCryoIon, HFCryoEtchant) with a surface model that solves steady-state HF coverage and computes etch rate from chemical, ion-enhanced, and sputter contributions. - viennaps.hpp: include psHFCryoEtching.hpp alongside other models. Co-Authored-By: Claude Sonnet 4.6 --- include/viennaps/models/psHFCryoEtching.hpp | 373 ++++++++++++++++++ .../viennaps/models/psHFCryoParameters.hpp | 65 +++ include/viennaps/viennaps.hpp | 1 + 3 files changed, 439 insertions(+) create mode 100644 include/viennaps/models/psHFCryoEtching.hpp create mode 100644 include/viennaps/models/psHFCryoParameters.hpp diff --git a/include/viennaps/models/psHFCryoEtching.hpp b/include/viennaps/models/psHFCryoEtching.hpp new file mode 100644 index 00000000..1f2ea7a2 --- /dev/null +++ b/include/viennaps/models/psHFCryoEtching.hpp @@ -0,0 +1,373 @@ +#pragma once + +#include +#include + +#include "../process/psProcessModel.hpp" +#include "../psUnits.hpp" + +#include "psHFCryoParameters.hpp" +#include "psIonModelUtil.hpp" + +namespace viennaps { + +using namespace viennacore; + +// HF cryogenic plasma etching of SiO2. +// +// Physics: +// SiO2 + 4HF -> SiF4(g) + 2H2O +// +// Coverage (steady-state): +// phi_HF = Gamma_HF / (Gamma_HF + k_des(T) + k_r(T) + Y_ie * Gamma_ion) +// +// k_des(T) = nu0 * exp(-E_des / kB*T) [Frenkel-Arrhenius desorption] +// k_r(T) = A_r * exp(-E_a / kB*T) [Arrhenius reaction rate] +// +// Etch rate: +// v = -(1/rho) * (k_r(T)*phi_HF + Y_sp*Gamma_ion + Y_ie*Gamma_ion*phi_HF) +// +// No passivation layer is needed: anisotropy is provided by ion bombardment +// and temperature-suppressed chemical etching on sidewalls. + +namespace impl { + +// ── Ion ────────────────────────────────────────────────────────────────────── +template +class HFCryoIon final + : public viennaray::Particle, NumericType> { +public: + explicit HFCryoIon(const HFCryoParameters &pParams) + : params(pParams), + A_energy(1. / (1. + params.Ions.n_l * + (M_PI_2 / params.Ions.inflectAngle - 1.))), + sqrt_E_th_ie(std::sqrt(params.SiO2.Eth_ie)) {} + + void surfaceCollision(NumericType rayWeight, const Vec3D &rayDir, + const Vec3D &geomNormal, + const unsigned int primID, const int /*materialId*/, + viennaray::TracingData &localData, + const viennaray::TracingData *, + RNG &) override { + const NumericType cosTheta = + -DotProduct(rayDir, geomNormal); + const NumericType angle = + std::acos(std::max(std::min(cosTheta, NumericType(1)), NumericType(0))); + + // Angular factor: ion-enhanced etching drops off at glancing angles + NumericType f_ie = 1.; + if (cosTheta < 0.5) + f_ie = std::max(NumericType(3) - NumericType(6) * angle / NumericType(M_PI), + NumericType(0)); + + const NumericType sqrtE = std::sqrt(E); + + // Physical sputtering yield + const NumericType Y_sp = + params.SiO2.A_sp * + std::max(sqrtE - std::sqrt(params.SiO2.Eth_sp), NumericType(0)); + + // Ion-enhanced chemical etching yield + const NumericType Y_ie = + params.SiO2.A_ie * + std::max(sqrtE - sqrt_E_th_ie, NumericType(0)) * f_ie; + + // ionSputterFlux + localData.getVectorData(0)[primID] += Y_sp * rayWeight; + // ionEnhancedFlux + localData.getVectorData(1)[primID] += Y_ie * rayWeight; + } + + std::pair> + surfaceReflection(NumericType rayWeight, const Vec3D &rayDir, + const Vec3D &geomNormal, + const unsigned int /*primId*/, const int /*materialId*/, + const viennaray::TracingData *, + RNG &rng) override { + const NumericType cosTheta = + -DotProduct(rayDir, geomNormal); + const NumericType incAngle = + std::acos(std::max(std::min(cosTheta, NumericType(1)), NumericType(0))); + + const NumericType newEnergy = updateEnergy( + rng, E, incAngle, A_energy, params.Ions.inflectAngle, params.Ions.n_l); + + if (newEnergy > params.SiO2.Eth_ie) { + E = newEnergy; + auto dir = viennaray::ReflectionConedCosine( + rayDir, geomNormal, rng, + M_PI_2 - std::min(incAngle, params.Ions.minAngle)); + return {NumericType(0), dir}; + } + return VIENNARAY_PARTICLE_STOP; + } + + void initNew(RNG &rng) override { + E = initNormalDistEnergy(rng, params.Ions.meanEnergy, + params.Ions.sigmaEnergy); + } + + NumericType getSourceDistributionPower() const override { + return params.Ions.exponent; + } + + [[nodiscard]] std::vector getLocalDataLabels() const override { + return {"ionSputterFlux", "ionEnhancedFlux"}; + } + +private: + const HFCryoParameters ¶ms; + const NumericType A_energy; + const NumericType sqrt_E_th_ie; + NumericType E = 0.; +}; + +// ── HF etchant ─────────────────────────────────────────────────────────────── +template +class HFCryoEtchant final + : public viennaray::Particle, NumericType> { +public: + explicit HFCryoEtchant(const HFCryoParameters &pParams) + : params(pParams) {} + + void surfaceCollision(NumericType rayWeight, const Vec3D &, + const Vec3D &, const unsigned int primID, + const int /*materialId*/, + viennaray::TracingData &localData, + const viennaray::TracingData *globalData, + RNG &) override { + // Effective sticking: HF only adsorbs on uncovered surface + const NumericType phi_HF = globalData->getVectorData(0)[primID]; + const NumericType S_eff = + params.gamma_HF * std::max(NumericType(1) - phi_HF, NumericType(0)); + + localData.getVectorData(0)[primID] += rayWeight * S_eff; + } + + std::pair> + surfaceReflection(NumericType /*rayWeight*/, const Vec3D &, + const Vec3D &geomNormal, + const unsigned int primID, const int /*materialId*/, + const viennaray::TracingData *globalData, + RNG &rng) override { + const NumericType phi_HF = globalData->getVectorData(0)[primID]; + const NumericType S_eff = + params.gamma_HF * std::max(NumericType(1) - phi_HF, NumericType(0)); + + auto dir = viennaray::ReflectionDiffuse(geomNormal, rng); + return {S_eff, dir}; + } + + NumericType getSourceDistributionPower() const override { return 1.; } + + [[nodiscard]] std::vector getLocalDataLabels() const override { + return {"etchantFlux"}; + } + +private: + const HFCryoParameters ¶ms; +}; + +// ── Surface model ───────────────────────────────────────────────────────────── +template +class HFCryoSurfaceModel : public SurfaceModel { +public: + using SurfaceModel::coverages; + using SurfaceModel::surfaceData; + + explicit HFCryoSurfaceModel(const HFCryoParameters &pParams) + : params(pParams) {} + + void initializeCoverages(unsigned numGeometryPoints) override { + if (coverages == nullptr) + coverages = viennals::PointData::New(); + else + coverages->clear(); + + std::vector zeros(numGeometryPoints, 0.); + coverages->insertNextScalarData(zeros, "HFCoverage"); + } + + void initializeSurfaceData(unsigned numGeometryPoints) override { + if (!Logger::hasIntermediate()) + return; + + if (surfaceData == nullptr) + surfaceData = viennals::PointData::New(); + else + surfaceData->clear(); + + std::vector zeros(numGeometryPoints, 0.); + surfaceData->insertNextScalarData(zeros, "chemicalRate"); + surfaceData->insertNextScalarData(zeros, "ionEnhancedRate"); + surfaceData->insertNextScalarData(zeros, "sputterRate"); + } + + SmartPointer> + calculateVelocities(SmartPointer> rates, + const std::vector> &coordinates, + const std::vector &materialIds) override { + updateCoverages(rates, materialIds); + + const auto numPoints = rates->getScalarData(0)->size(); + std::vector etchRate(numPoints, 0.); + + const auto ionSputterFlux = rates->getScalarData("ionSputterFlux"); + const auto ionEnhancedFlux = rates->getScalarData("ionEnhancedFlux"); + const auto phi_HF = coverages->getScalarData("HFCoverage"); + + std::vector *chRate = nullptr, *ieRate = nullptr, + *spRate = nullptr; + if (Logger::hasIntermediate()) { + chRate = surfaceData->getScalarData("chemicalRate"); + ieRate = surfaceData->getScalarData("ionEnhancedRate"); + spRate = surfaceData->getScalarData("sputterRate"); + chRate->resize(numPoints); + ieRate->resize(numPoints); + spRate->resize(numPoints); + } + + const NumericType kr = params.k_r(); + const NumericType unitConversion = + units::Time::convertSecond() / units::Length::convertNanometer(); + + bool stop = false; + +#pragma omp parallel for reduction(|| : stop) + for (size_t i = 0; i < numPoints; ++i) { + if (coordinates[i][D - 1] < params.etchStopDepth || stop) { + stop = true; + continue; + } + + if (!MaterialMap::isMaterial(materialIds[i], Material::SiO2)) + continue; + + const NumericType chemical = + kr * phi_HF->at(i); + const NumericType ionEnhanced = + phi_HF->at(i) * ionEnhancedFlux->at(i) * params.ionFlux; + const NumericType sputter = + ionSputterFlux->at(i) * params.ionFlux; + + etchRate[i] = -(1. / params.SiO2.rho) * + (chemical + ionEnhanced + sputter) * unitConversion; + + if (Logger::hasIntermediate()) { + chRate->at(i) = chemical; + ieRate->at(i) = ionEnhanced; + spRate->at(i) = sputter; + } + } + + if (stop) { + std::fill(etchRate.begin(), etchRate.end(), NumericType(0)); + VIENNACORE_LOG_INFO("Etch stop depth reached."); + } + + return SmartPointer>::New(std::move(etchRate)); + } + + void updateCoverages(SmartPointer> rates, + const std::vector & /*materialIds*/) override { + const auto numPoints = rates->getScalarData(0)->size(); + const auto etchantFlux = rates->getScalarData("etchantFlux"); + const auto ionEnhancedFlux = rates->getScalarData("ionEnhancedFlux"); + + auto phi_HF = coverages->getScalarData("HFCoverage"); + phi_HF->resize(numPoints); + + const NumericType k_des = params.k_des(); + const NumericType k_r = params.k_r(); + +#pragma omp parallel for + for (size_t i = 0; i < numPoints; ++i) { + // Adsorption flux (scaled) + const NumericType Gamma_HF = etchantFlux->at(i) * params.etchantFlux; + // Ion-enhanced removal of adsorbed HF + const NumericType Y_ie_ion = + ionEnhancedFlux->at(i) * params.ionFlux; + + // Steady-state coverage: + // phi = Gamma_HF / (Gamma_HF + k_des + k_r + Y_ie_ion) + const NumericType denom = Gamma_HF + k_des + k_r + Y_ie_ion; + phi_HF->at(i) = (denom > NumericType(1e-30)) + ? std::min(Gamma_HF / denom, NumericType(1)) + : NumericType(0); + } + } + +private: + const HFCryoParameters ¶ms; +}; + +} // namespace impl + +// ── Public process model ────────────────────────────────────────────────────── +template +class HFCryoEtching final : public ProcessModelCPU { +public: + HFCryoEtching() { initializeModel(); } + + // Convenience constructor + HFCryoEtching(NumericType ionFlux, NumericType etchantFlux, + NumericType temperature, // K + NumericType meanEnergy, // eV + NumericType sigmaEnergy, // eV + NumericType E_des = 0.25, // eV (Frenkel desorption) + NumericType E_a = 0.10, // eV (reaction activation) + NumericType etchStopDepth = + std::numeric_limits::lowest()) { + params.ionFlux = ionFlux; + params.etchantFlux = etchantFlux; + params.temperature = temperature; + params.Ions.meanEnergy = meanEnergy; + params.Ions.sigmaEnergy = sigmaEnergy; + params.Desorption.E_des = E_des; + params.Reaction.E_a = E_a; + params.etchStopDepth = etchStopDepth; + initializeModel(); + } + + explicit HFCryoEtching(const HFCryoParameters &pParams) + : params(pParams) { + initializeModel(); + } + + void setParameters(const HFCryoParameters &pParams) { + params = pParams; + initializeModel(); + } + + HFCryoParameters &getParameters() { return params; } + +private: + void initializeModel() { + if (units::Length::getUnit() == units::Length::UNDEFINED || + units::Time::getUnit() == units::Time::UNDEFINED) { + VIENNACORE_LOG_ERROR("Units have not been set."); + } + + auto ion = std::make_unique>(params); + auto etchant = + std::make_unique>(params); + + auto surfModel = + SmartPointer>::New(params); + auto velField = + SmartPointer>::New(); + + this->setSurfaceModel(surfModel); + this->setVelocityField(velField); + this->setProcessName("HFCryoEtching"); + this->particles.clear(); + this->insertNextParticleType(ion); + this->insertNextParticleType(etchant); + } + + HFCryoParameters params; +}; + +PS_PRECOMPILE_PRECISION_DIMENSION(HFCryoEtching) + +} // namespace viennaps diff --git a/include/viennaps/models/psHFCryoParameters.hpp b/include/viennaps/models/psHFCryoParameters.hpp new file mode 100644 index 00000000..08c7d179 --- /dev/null +++ b/include/viennaps/models/psHFCryoParameters.hpp @@ -0,0 +1,65 @@ +#pragma once + +#include "../psConstants.hpp" +#include + +namespace viennaps { + +template struct HFCryoParameters { + // Fluxes in units of (1e15 / cm² / s) + NumericType ionFlux = 1.; + NumericType etchantFlux = 1.0e3; + + // Substrate temperature (K) + NumericType temperature = 200.; // cryogenic, e.g. 173-223 K + + // HF sticking probability on bare SiO2 + NumericType gamma_HF = 0.9; + + NumericType etchStopDepth = std::numeric_limits::lowest(); + + // Frenkel-Arrhenius desorption: k_des = nu0 * exp(-E_des / kB*T) + struct DesorptionType { + NumericType nu0 = 1.0e13; // attempt frequency (1/s) + NumericType E_des = 0.25; // desorption activation energy (eV) + } Desorption; + + // Arrhenius reaction rate: k_r = A_r * exp(-E_a / kB*T) + struct ReactionType { + NumericType A_r = 3.0e2; // pre-exponential factor (1e15 cm⁻²s⁻¹) + NumericType E_a = 0.10; // reaction activation energy (eV) + } Reaction; + + // SiO2 material properties + struct SiO2Type { + NumericType rho = constants::SiO2::rho; // 1e22 atoms/cm³ + NumericType Eth_sp = 18.; // sputtering threshold (eV) + NumericType A_sp = constants::SiO2::A_sp; + NumericType Eth_ie = 12.; // ion-enhanced etching threshold (eV) + NumericType A_ie = 2.0; + } SiO2; + + // Ion properties + struct IonType { + NumericType meanEnergy = 100.; // eV + NumericType sigmaEnergy = 10.; // eV + NumericType exponent = 300.; // angular distribution exponent + NumericType inflectAngle = constants::Ion::inflectAngle; + NumericType n_l = constants::Ion::n_l; + NumericType minAngle = constants::Ion::minAngle; + } Ions; + + // Compute k_des(T) in units consistent with fluxes (1/s) + NumericType k_des() const { + return Desorption.nu0 * + std::exp(-Desorption.E_des / (constants::kB * temperature)); + } + + // Compute k_r(T) in units of (1e15 cm⁻²s⁻¹) + NumericType k_r() const { + return Reaction.A_r * + std::exp(-Reaction.E_a / (constants::kB * temperature)); + } +}; + +} // namespace viennaps diff --git a/include/viennaps/viennaps.hpp b/include/viennaps/viennaps.hpp index 8e6b34f7..be5623aa 100644 --- a/include/viennaps/viennaps.hpp +++ b/include/viennaps/viennaps.hpp @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include From 5a70b52824bb7331da0ebd92b88514d2aa2cd647 Mon Sep 17 00:00:00 2001 From: monst Date: Tue, 5 May 2026 15:58:05 +0900 Subject: [PATCH 2/9] feat: add HFCryoEtching example Adds a minimal 2D trench example for the HFCryoEtching model demonstrating: - SiO2 trench geometry setup - HFCryoParameters configuration (temperature, Arrhenius E_des / E_a) - Process execution with ray tracing and Runge-Kutta advection - Temperature sweep output showing k_r(T) and k_des(T) at 150-300 K Co-Authored-By: Claude Sonnet 4.6 --- examples/HFCryoEtching/CMakeLists.txt | 3 + examples/HFCryoEtching/HFCryoEtching.cpp | 102 +++++++++++++++++++++++ 2 files changed, 105 insertions(+) create mode 100644 examples/HFCryoEtching/CMakeLists.txt create mode 100644 examples/HFCryoEtching/HFCryoEtching.cpp diff --git a/examples/HFCryoEtching/CMakeLists.txt b/examples/HFCryoEtching/CMakeLists.txt new file mode 100644 index 00000000..2248dc39 --- /dev/null +++ b/examples/HFCryoEtching/CMakeLists.txt @@ -0,0 +1,3 @@ +project(HFCryoEtching LANGUAGES CXX) + +viennaps_add_example(${PROJECT_NAME} "${PROJECT_NAME}.cpp") diff --git a/examples/HFCryoEtching/HFCryoEtching.cpp b/examples/HFCryoEtching/HFCryoEtching.cpp new file mode 100644 index 00000000..399616a8 --- /dev/null +++ b/examples/HFCryoEtching/HFCryoEtching.cpp @@ -0,0 +1,102 @@ +#include +#include +#include +#include + +using namespace viennaps; + +int main() { + using NumericType = double; + constexpr int D = 2; + + Logger::setLogLevel(LogLevel::INFO); + omp_set_num_threads(4); + + // Units + units::Length::setUnit(units::Length::Nanometer); + units::Time::setUnit(units::Time::Second); + + // Geometry: SiO2 trench with mask + NumericType gridDelta = 1.0; // nm + NumericType xExtent = 60.0; // nm + NumericType yExtent = 60.0; // nm + NumericType trenchWidth = 20.0; // nm + NumericType maskHeight = 20.0; // nm + + auto domain = SmartPointer>::New( + gridDelta, xExtent, yExtent); + + // Fill substrate with SiO2, mask with Mask material + MakeTrench(domain, + trenchWidth, + 0.0, // trench depth (0 = flat start) + 0.0, // taper angle + maskHeight, + 0.0, // mask taper angle + false) // periodic boundary + .apply(); + + // Replace default Si substrate material with SiO2 + domain->setMaterial(Material::SiO2); + + // HF cryo etching parameters + HFCryoParameters params; + params.ionFlux = 1.0; // 1e15 /cm²/s + params.etchantFlux = 1.0e3; // 1e15 /cm²/s + params.temperature = 200.; // K (cryogenic ~-73°C) + params.gamma_HF = 0.9; // HF sticking probability on bare SiO2 + + // Frenkel-Arrhenius desorption + params.Desorption.nu0 = 1.0e13; // 1/s + params.Desorption.E_des = 0.25; // eV + + // Arrhenius reaction rate + params.Reaction.A_r = 3.0e2; // 1e15 cm⁻²s⁻¹ + params.Reaction.E_a = 0.10; // eV + + // Ion properties + params.Ions.meanEnergy = 100.; // eV + params.Ions.sigmaEnergy = 10.; // eV + params.Ions.exponent = 300.; // angular distribution + + auto model = SmartPointer>::New(params); + + // Print computed Arrhenius rates at the set temperature + std::cout << "Temperature : " << params.temperature << " K\n"; + std::cout << "k_des(T) : " << params.k_des() << " /s\n"; + std::cout << "k_r(T) : " << params.k_r() << " (1e15 cm⁻²s⁻¹)\n"; + + // Process + NumericType processTime = 30.; // seconds + + Process process(domain, model, processTime); + + RayTracingParameters rayParams; + rayParams.raysPerPoint = 3000; + process.setParameters(rayParams); + + AdvectionParameters advParams; + advParams.spatialScheme = SpatialScheme::WENO_3RD_ORDER; + advParams.temporalScheme = TemporalScheme::RUNGE_KUTTA_2ND_ORDER; + process.setParameters(advParams); + + // Save initial surface + domain->saveSurfaceMesh("HFCryo_initial.vtp"); + + std::cout << "Running HF cryo etching simulation...\n"; + process.apply(); + + // Save final surface + domain->saveSurfaceMesh("HFCryo_final.vtp"); + std::cout << "Done. Output: HFCryo_final.vtp\n"; + + // --- Temperature sweep: compare 150 K vs 200 K vs 250 K --- + std::cout << "\n--- Temperature effect on k_r ---\n"; + for (NumericType T : {150., 200., 250., 300.}) { + params.temperature = T; + std::cout << "T=" << T << " K k_r=" << params.k_r() + << " k_des=" << params.k_des() << "\n"; + } + + return 0; +} From 527ddeed98e9ac02b8f289763a1e49e165e6d88d Mon Sep 17 00:00:00 2001 From: monst Date: Wed, 6 May 2026 17:52:17 +0900 Subject: [PATCH 3/9] feat: extend HFCryoEtching with physisorption, surface diffusion, and 4-case example Physical model additions to psHFCryoEtching / psHFCryoParameters: - Two-state adsorption: theta_phys (physisorbed) + theta_chem (chemisorbed) with steady-state denominators including all loss pathways - Pure chemical etching term k_r_direct * theta_phys (thermal, no ions) - Steady-state surface diffusion PDE via Thomas algorithm on 1D surface chain - ModelConfig flags: useTemperatureDependence, usePhysisorption, useSurfaceDiffusion - DirectReactionType parameters and effective_k_* helper methods - Fix surfaceReflection: return {1-S_eff, dir} for correct Knudsen transport New example HFCryo4Cases: - 4 cases x 4 temperatures (150/200/250/300 K) = 16 simulations - VTP output with Temperature_K and CaseID scalars for ParaView Co-Authored-By: Claude Sonnet 4.6 --- examples/HFCryo4Cases/CMakeLists.txt | 3 + examples/HFCryo4Cases/HFCryo4Cases.cpp | 167 +++++++ examples/HFCryoEtching/HFCryoEtching.cpp | 177 ++++--- include/viennaps/models/psHFCryoEtching.hpp | 433 +++++++++++++----- .../viennaps/models/psHFCryoParameters.hpp | 63 ++- 5 files changed, 648 insertions(+), 195 deletions(-) create mode 100644 examples/HFCryo4Cases/CMakeLists.txt create mode 100644 examples/HFCryo4Cases/HFCryo4Cases.cpp diff --git a/examples/HFCryo4Cases/CMakeLists.txt b/examples/HFCryo4Cases/CMakeLists.txt new file mode 100644 index 00000000..806a3511 --- /dev/null +++ b/examples/HFCryo4Cases/CMakeLists.txt @@ -0,0 +1,3 @@ +project(HFCryo4Cases LANGUAGES CXX) + +viennaps_add_example(${PROJECT_NAME} "${PROJECT_NAME}.cpp") diff --git a/examples/HFCryo4Cases/HFCryo4Cases.cpp b/examples/HFCryo4Cases/HFCryo4Cases.cpp new file mode 100644 index 00000000..4de95500 --- /dev/null +++ b/examples/HFCryo4Cases/HFCryo4Cases.cpp @@ -0,0 +1,167 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +using namespace viennaps; + +static constexpr double gridDelta = 1.0; +static constexpr double xExtent = 80.0; +static constexpr double yExtent = 2000.0; // 700 → 2000 (깊은 트렌치 수용) +static constexpr double trenchWidth = 10.0; +static constexpr double maskHeight = 80.0; // 50 → 80 (두꺼운 마스크) +static constexpr double processTime = 400.0; // 90 → 400 (긴 공정 시간) + +// Save mesh with Temperature_K and CaseID as point scalars. +// In ParaView: Color by "Temperature_K" → temperature gradient (150-300K) +// Color by "CaseID" → 4 distinct case colors (0-3) +static void saveColoredMesh(SmartPointer> domain, + const std::string &filename, + double T, int caseId) { + auto mesh = viennals::Mesh::New(); + viennals::ToMultiSurfaceMesh(domain->getLevelSets(), mesh).apply(); + + const size_t n = mesh->nodes.size(); + mesh->getPointData().insertNextScalarData( + std::vector(n, T), "Temperature_K"); + mesh->getPointData().insertNextScalarData( + std::vector(n, double(caseId)), "CaseID"); + + viennals::VTKWriter(mesh, filename).apply(); +} + +static void runCase(const std::string &filename, + HFCryoParameters params, + int caseId) { + auto domain = SmartPointer>::New(gridDelta, xExtent, yExtent); + MakeTrench(domain, trenchWidth, 0.0, 0.0, + maskHeight, 0.0, false, Material::SiO2).apply(); + + auto model = SmartPointer>::New(params); + Process process(domain, model, processTime); + + RayTracingParameters rayParams; + rayParams.raysPerPoint = 5000; // 3000 → 5000 (깊은 트렌치 정확도↑) + process.setParameters(rayParams); + + AdvectionParameters advParams; + advParams.spatialScheme = SpatialScheme::WENO_3RD_ORDER; + advParams.temporalScheme = TemporalScheme::RUNGE_KUTTA_2ND_ORDER; + process.setParameters(advParams); + + process.apply(); + + saveColoredMesh(domain, filename + ".vtp", params.temperature, caseId); + std::cout << "done -> " << filename << ".vtp\n"; +} + +int main() { + Logger::setLogLevel(LogLevel::WARNING); + omp_set_num_threads(4); + + units::Length::setUnit(units::Length::NANOMETER); + units::Time::setUnit(units::Time::SECOND); + + // ── Base parameters ───────────────────────────────────────────────────────── + HFCryoParameters base; + base.ionFlux = 3.0; // 1.0 → 3.0 (이온 강화 → 수직 식각↑) + base.etchantFlux = 5.0e3; // 1e3 → 5e3 (HF 공급↑) + base.gamma_HF = 0.9; + base.Desorption.nu0 = 1.0e8; + base.Desorption.E_des = 0.20; + base.Reaction.A_r = 3.0e2; + base.Reaction.E_a = 0.10; + base.DirectReaction.A_r = 1.0e4; + base.DirectReaction.E_a = 0.25; + base.IonActivation.A_act = 1.0; + base.Diffusion.D0 = 0.0; + base.Diffusion.omega = 0.25; + base.Ions.meanEnergy = 150.0; // 100 → 150 eV (더 높은 이온 에너지) + base.Ions.sigmaEnergy = 10.0; + base.Ions.exponent = 600.0; // 300 → 600 (이온빔 더 수직, 이방성↑) + base.Config.T_ref = 300.0; + + // ── Case definitions ───────────────────────────────────────────────────────── + struct CaseInfo { + std::string tag; + std::string desc; + bool tempDep; + bool phys; + bool diff; + double D0; + }; + + const CaseInfo cases[4] = { + {"C0_none", "nothing (constant k_des, single theta, no diffusion)", + false, false, false, 0.0}, + {"C1_temp", "+Arrhenius k_des (single theta, no diffusion)", + true, false, false, 0.0}, + {"C2_phys", "+physisorption (theta_phys+theta_chem, no diffusion)", + true, true, false, 0.0}, + {"C3_diff", "+surface diffusion (D0=1e3)", + true, true, true, 1.0e3}, + }; + + const double Tlist[4] = {150., 200., 250., 300.}; + + // ── Print rate table ───────────────────────────────────────────────────────── + std::cout << "=== HF Cryo 4-Case x 4-Temperature Comparison ===\n"; + std::cout << "Geometry: " << trenchWidth << "nm trench, " + << processTime << "s process\n\n"; + + std::cout << "Effective k_des per case/temperature:\n"; + std::cout << std::setw(20) << "" << std::setw(10) << "150K" + << std::setw(10) << "200K" << std::setw(10) << "250K" + << std::setw(10) << "300K\n"; + std::cout << std::string(60, '-') << "\n"; + for (const auto &ci : cases) { + std::cout << std::setw(20) << ci.tag; + for (double T : Tlist) { + auto p = base; + p.temperature = T; + p.Config.useTemperatureDependence = ci.tempDep; + std::cout << std::setw(10) << std::setprecision(1) << std::fixed + << p.effective_k_des(); + } + std::cout << "\n"; + } + std::cout << "\n"; + + // ── Run all 16 simulations ─────────────────────────────────────────────────── + for (int c = 0; c < 4; ++c) { + const auto &ci = cases[c]; + std::cout << "=== Case " << c << ": " << ci.desc << " ===\n"; + + for (double T : Tlist) { + std::cout << " T=" << std::setw(3) << (int)T << "K ... "; + std::cout.flush(); + + auto p = base; + p.temperature = T; + p.Config.useTemperatureDependence = ci.tempDep; + p.Config.usePhysisorption = ci.phys; + p.Config.useSurfaceDiffusion = ci.diff; + p.Diffusion.D0 = ci.D0; + + std::string fname = ci.tag + "_T" + std::to_string((int)T) + "K"; + runCase(fname, p, c); + } + std::cout << "\n"; + } + + std::cout << "=== Done: 16 VTP files generated ===\n\n"; + std::cout << "ParaView: open all *.vtp files\n"; + std::cout << " Color by 'Temperature_K' -> blue(150K) to red(300K)\n"; + std::cout << " Color by 'CaseID' -> 4 discrete case colors\n"; + + return 0; +} diff --git a/examples/HFCryoEtching/HFCryoEtching.cpp b/examples/HFCryoEtching/HFCryoEtching.cpp index 399616a8..31f2508d 100644 --- a/examples/HFCryoEtching/HFCryoEtching.cpp +++ b/examples/HFCryoEtching/HFCryoEtching.cpp @@ -1,3 +1,7 @@ +#include +#include +#include + #include #include #include @@ -5,71 +9,20 @@ using namespace viennaps; -int main() { - using NumericType = double; - constexpr int D = 2; +static void runCase(const std::string &label, + double T, double gridDelta, + double xExtent, double yExtent, + double trenchWidth, double maskHeight, + double processTime, + HFCryoParameters params) { + params.temperature = T; - Logger::setLogLevel(LogLevel::INFO); - omp_set_num_threads(4); + auto domain = SmartPointer>::New(gridDelta, xExtent, yExtent); + MakeTrench(domain, trenchWidth, 0.0, 0.0, + maskHeight, 0.0, false, Material::SiO2).apply(); - // Units - units::Length::setUnit(units::Length::Nanometer); - units::Time::setUnit(units::Time::Second); - - // Geometry: SiO2 trench with mask - NumericType gridDelta = 1.0; // nm - NumericType xExtent = 60.0; // nm - NumericType yExtent = 60.0; // nm - NumericType trenchWidth = 20.0; // nm - NumericType maskHeight = 20.0; // nm - - auto domain = SmartPointer>::New( - gridDelta, xExtent, yExtent); - - // Fill substrate with SiO2, mask with Mask material - MakeTrench(domain, - trenchWidth, - 0.0, // trench depth (0 = flat start) - 0.0, // taper angle - maskHeight, - 0.0, // mask taper angle - false) // periodic boundary - .apply(); - - // Replace default Si substrate material with SiO2 - domain->setMaterial(Material::SiO2); - - // HF cryo etching parameters - HFCryoParameters params; - params.ionFlux = 1.0; // 1e15 /cm²/s - params.etchantFlux = 1.0e3; // 1e15 /cm²/s - params.temperature = 200.; // K (cryogenic ~-73°C) - params.gamma_HF = 0.9; // HF sticking probability on bare SiO2 - - // Frenkel-Arrhenius desorption - params.Desorption.nu0 = 1.0e13; // 1/s - params.Desorption.E_des = 0.25; // eV - - // Arrhenius reaction rate - params.Reaction.A_r = 3.0e2; // 1e15 cm⁻²s⁻¹ - params.Reaction.E_a = 0.10; // eV - - // Ion properties - params.Ions.meanEnergy = 100.; // eV - params.Ions.sigmaEnergy = 10.; // eV - params.Ions.exponent = 300.; // angular distribution - - auto model = SmartPointer>::New(params); - - // Print computed Arrhenius rates at the set temperature - std::cout << "Temperature : " << params.temperature << " K\n"; - std::cout << "k_des(T) : " << params.k_des() << " /s\n"; - std::cout << "k_r(T) : " << params.k_r() << " (1e15 cm⁻²s⁻¹)\n"; - - // Process - NumericType processTime = 30.; // seconds - - Process process(domain, model, processTime); + auto model = SmartPointer>::New(params); + Process process(domain, model, processTime); RayTracingParameters rayParams; rayParams.raysPerPoint = 3000; @@ -80,23 +33,95 @@ int main() { advParams.temporalScheme = TemporalScheme::RUNGE_KUTTA_2ND_ORDER; process.setParameters(advParams); - // Save initial surface - domain->saveSurfaceMesh("HFCryo_initial.vtp"); - - std::cout << "Running HF cryo etching simulation...\n"; process.apply(); - // Save final surface - domain->saveSurfaceMesh("HFCryo_final.vtp"); - std::cout << "Done. Output: HFCryo_final.vtp\n"; + std::string tag = label + "_T" + std::to_string((int)T) + "K"; + domain->saveSurfaceMesh(tag + ".vtp", true); + domain->saveHullMesh(tag + "_hull.vtp"); + std::cout << "done -> " << tag << ".vtp\n"; +} + +int main() { + Logger::setLogLevel(LogLevel::WARNING); + omp_set_num_threads(4); + + units::Length::setUnit(units::Length::NANOMETER); + units::Time::setUnit(units::Time::SECOND); + + const double gridDelta = 1.0; + const double xExtent = 80.0; + const double yExtent = 700.0; + const double trenchWidth = 10.0; + const double maskHeight = 50.0; + const double processTime = 90.0; + + // Base parameters shared by both cases + HFCryoParameters params; + params.ionFlux = 1.0; + params.etchantFlux = 1.0e3; + params.gamma_HF = 0.9; + params.Desorption.nu0 = 1.0e8; + params.Desorption.E_des = 0.20; + params.Reaction.A_r = 3.0e2; + params.Reaction.E_a = 0.10; + params.IonActivation.A_act = 1.0; + params.Ions.meanEnergy = 100.0; + params.Ions.sigmaEnergy = 10.0; + params.Ions.exponent = 300.0; + + // Print D_s per temperature + std::cout << "=== HF Cryo Etching | Surface Diffusion Comparison ===\n\n"; + std::cout << "T(K) k_des k_r k_r_direct D_s(D0=1e3)\n"; + std::cout << "------------------------------------------------------\n"; + double Tlist[4] = {150.0, 200.0, 250.0, 300.0}; + for (int i = 0; i < 4; ++i) { + params.temperature = Tlist[i]; + params.Diffusion.D0 = 1.0e3; + params.Diffusion.omega = 0.25; + std::cout << Tlist[i] << "K " + << params.k_des() << " " + << params.k_r() << " " + << params.k_r_direct() << " " + << params.D_s() << "\n"; + } + std::cout << "\n"; + + // Save initial geometry + { + auto domain = SmartPointer>::New(gridDelta, xExtent, yExtent); + MakeTrench(domain, trenchWidth, 0.0, 0.0, + maskHeight, 0.0, false, Material::SiO2).apply(); + domain->saveSurfaceMesh("HFCryo_initial.vtp", true); + std::cout << "Initial geometry saved.\n\n"; + } - // --- Temperature sweep: compare 150 K vs 200 K vs 250 K --- - std::cout << "\n--- Temperature effect on k_r ---\n"; - for (NumericType T : {150., 200., 250., 300.}) { - params.temperature = T; - std::cout << "T=" << T << " K k_r=" << params.k_r() - << " k_des=" << params.k_des() << "\n"; + // ── Case A: No surface diffusion (D0=0) ────────────────────────────────── + std::cout << "=== Case A: No Surface Diffusion (D0=0) ===\n"; + params.Diffusion.D0 = 0.0; + params.Diffusion.omega = 0.25; + for (int i = 0; i < 4; ++i) { + std::cout << " T=" << Tlist[i] << "K ... "; + std::cout.flush(); + runCase("NoDiff", Tlist[i], gridDelta, xExtent, yExtent, + trenchWidth, maskHeight, processTime, params); } + // ── Case B: With surface diffusion (D0=1e3) ─────────────────────────────── + std::cout << "\n=== Case B: With Surface Diffusion (D0=1e3) ===\n"; + params.Diffusion.D0 = 1.0e3; + params.Diffusion.omega = 0.25; + for (int i = 0; i < 4; ++i) { + std::cout << " T=" << Tlist[i] << "K ... "; + std::cout.flush(); + runCase("Diff", Tlist[i], gridDelta, xExtent, yExtent, + trenchWidth, maskHeight, processTime, params); + } + + std::cout << "\n=== Done ===\n"; + std::cout << "Compare in ParaView:\n"; + std::cout << " NoDiff_T150K.vtp vs Diff_T150K.vtp\n"; + std::cout << " NoDiff_T200K.vtp vs Diff_T200K.vtp\n"; + std::cout << "Color by MaterialIds to see Mask / SiO2\n"; + return 0; } diff --git a/include/viennaps/models/psHFCryoEtching.hpp b/include/viennaps/models/psHFCryoEtching.hpp index 1f2ea7a2..5332a6ca 100644 --- a/include/viennaps/models/psHFCryoEtching.hpp +++ b/include/viennaps/models/psHFCryoEtching.hpp @@ -1,5 +1,9 @@ #pragma once +#include +#include +#include + #include #include @@ -15,20 +19,21 @@ using namespace viennacore; // HF cryogenic plasma etching of SiO2. // -// Physics: -// SiO2 + 4HF -> SiF4(g) + 2H2O -// -// Coverage (steady-state): -// phi_HF = Gamma_HF / (Gamma_HF + k_des(T) + k_r(T) + Y_ie * Gamma_ion) +// Two-state adsorption + surface diffusion model: // -// k_des(T) = nu0 * exp(-E_des / kB*T) [Frenkel-Arrhenius desorption] -// k_r(T) = A_r * exp(-E_a / kB*T) [Arrhenius reaction rate] +// Physisorbed HF (theta_phys) — solved via 1D diffusion PDE along surface: +// D_s * d²theta/ds² + Gamma_HF*(1-theta) - k_des*theta +// - k_act*Gamma_ion*theta = 0 +// BC: Dirichlet at chain ends (open surface = local Langmuir equilibrium) // -// Etch rate: -// v = -(1/rho) * (k_r(T)*phi_HF + Y_sp*Gamma_ion + Y_ie*Gamma_ion*phi_HF) +// Chemisorbed HF (theta_chem) — ion-activated, zero on sidewalls (Lill 2023): +// theta_chem = k_act * Gamma_ion * theta_phys / k_r (bottom/non-sidewall) +// theta_chem = 0 (sidewall, per paper) // -// No passivation layer is needed: anisotropy is provided by ion bombardment -// and temperature-suppressed chemical etching on sidewalls. +// Etch rate: +// v = -(Y_ie * Gamma_ion * theta_phys +// + k_r * theta_chem +// + Y_sp * Gamma_ion) / rho namespace impl { @@ -45,16 +50,14 @@ class HFCryoIon final void surfaceCollision(NumericType rayWeight, const Vec3D &rayDir, const Vec3D &geomNormal, - const unsigned int primID, const int /*materialId*/, + const unsigned int primID, const int, viennaray::TracingData &localData, const viennaray::TracingData *, RNG &) override { - const NumericType cosTheta = - -DotProduct(rayDir, geomNormal); + const NumericType cosTheta = -DotProduct(rayDir, geomNormal); const NumericType angle = std::acos(std::max(std::min(cosTheta, NumericType(1)), NumericType(0))); - // Angular factor: ion-enhanced etching drops off at glancing angles NumericType f_ie = 1.; if (cosTheta < 0.5) f_ie = std::max(NumericType(3) - NumericType(6) * angle / NumericType(M_PI), @@ -62,30 +65,25 @@ class HFCryoIon final const NumericType sqrtE = std::sqrt(E); - // Physical sputtering yield const NumericType Y_sp = params.SiO2.A_sp * std::max(sqrtE - std::sqrt(params.SiO2.Eth_sp), NumericType(0)); - // Ion-enhanced chemical etching yield const NumericType Y_ie = params.SiO2.A_ie * std::max(sqrtE - sqrt_E_th_ie, NumericType(0)) * f_ie; - // ionSputterFlux localData.getVectorData(0)[primID] += Y_sp * rayWeight; - // ionEnhancedFlux localData.getVectorData(1)[primID] += Y_ie * rayWeight; } std::pair> - surfaceReflection(NumericType rayWeight, const Vec3D &rayDir, + surfaceReflection(NumericType, const Vec3D &rayDir, const Vec3D &geomNormal, - const unsigned int /*primId*/, const int /*materialId*/, + const unsigned int, const int, const viennaray::TracingData *, RNG &rng) override { - const NumericType cosTheta = - -DotProduct(rayDir, geomNormal); + const NumericType cosTheta = -DotProduct(rayDir, geomNormal); const NumericType incAngle = std::acos(std::max(std::min(cosTheta, NumericType(1)), NumericType(0))); @@ -103,8 +101,7 @@ class HFCryoIon final } void initNew(RNG &rng) override { - E = initNormalDistEnergy(rng, params.Ions.meanEnergy, - params.Ions.sigmaEnergy); + E = initNormalDistEnergy(rng, params.Ions.meanEnergy, params.Ions.sigmaEnergy); } NumericType getSourceDistributionPower() const override { @@ -112,7 +109,7 @@ class HFCryoIon final } [[nodiscard]] std::vector getLocalDataLabels() const override { - return {"ionSputterFlux", "ionEnhancedFlux"}; + return {"ionSputterFlux", "ionActivationFlux"}; } private: @@ -132,30 +129,33 @@ class HFCryoEtchant final void surfaceCollision(NumericType rayWeight, const Vec3D &, const Vec3D &, const unsigned int primID, - const int /*materialId*/, + const int, viennaray::TracingData &localData, const viennaray::TracingData *globalData, RNG &) override { - // Effective sticking: HF only adsorbs on uncovered surface - const NumericType phi_HF = globalData->getVectorData(0)[primID]; + const NumericType theta_phys = globalData->getVectorData(0)[primID]; + const NumericType theta_chem = globalData->getVectorData(1)[primID]; + const NumericType occupied = + std::min(theta_phys + theta_chem, NumericType(1)); const NumericType S_eff = - params.gamma_HF * std::max(NumericType(1) - phi_HF, NumericType(0)); - + params.gamma_HF * std::max(NumericType(1) - occupied, NumericType(0)); localData.getVectorData(0)[primID] += rayWeight * S_eff; } std::pair> - surfaceReflection(NumericType /*rayWeight*/, const Vec3D &, + surfaceReflection(NumericType, const Vec3D &, const Vec3D &geomNormal, - const unsigned int primID, const int /*materialId*/, + const unsigned int primID, const int, const viennaray::TracingData *globalData, RNG &rng) override { - const NumericType phi_HF = globalData->getVectorData(0)[primID]; + const NumericType theta_phys = globalData->getVectorData(0)[primID]; + const NumericType theta_chem = globalData->getVectorData(1)[primID]; + const NumericType occupied = + std::min(theta_phys + theta_chem, NumericType(1)); const NumericType S_eff = - params.gamma_HF * std::max(NumericType(1) - phi_HF, NumericType(0)); - + params.gamma_HF * std::max(NumericType(1) - occupied, NumericType(0)); auto dir = viennaray::ReflectionDiffuse(geomNormal, rng); - return {S_eff, dir}; + return {NumericType(1) - S_eff, dir}; // reflected weight = 1 - sticking prob } NumericType getSourceDistributionPower() const override { return 1.; } @@ -185,13 +185,13 @@ class HFCryoSurfaceModel : public SurfaceModel { coverages->clear(); std::vector zeros(numGeometryPoints, 0.); - coverages->insertNextScalarData(zeros, "HFCoverage"); + coverages->insertNextScalarData(zeros, "physCoverage"); // index 0 + coverages->insertNextScalarData(zeros, "chemCoverage"); // index 1 } void initializeSurfaceData(unsigned numGeometryPoints) override { if (!Logger::hasIntermediate()) return; - if (surfaceData == nullptr) surfaceData = viennals::PointData::New(); else @@ -199,63 +199,150 @@ class HFCryoSurfaceModel : public SurfaceModel { std::vector zeros(numGeometryPoints, 0.); surfaceData->insertNextScalarData(zeros, "chemicalRate"); - surfaceData->insertNextScalarData(zeros, "ionEnhancedRate"); surfaceData->insertNextScalarData(zeros, "sputterRate"); } + // Local (point-wise) steady-state — no diffusion. + // Called first; result is used as Dirichlet BC at chain endpoints. + void updateCoverages(SmartPointer> rates, + const std::vector &) override { + const auto numPoints = rates->getScalarData(0)->size(); + const auto etchantFlux = rates->getScalarData("etchantFlux"); + const auto ionActFlux = rates->getScalarData("ionActivationFlux"); + + auto theta_phys = coverages->getScalarData("physCoverage"); + auto theta_chem = coverages->getScalarData("chemCoverage"); + theta_phys->resize(numPoints); + theta_chem->resize(numPoints); + + const NumericType k_des = params.effective_k_des(); + const NumericType k_r = params.effective_k_r(); + const NumericType k_r_direct = params.effective_k_r_direct(); + const NumericType A_act = params.IonActivation.A_act; + +#pragma omp parallel for + for (size_t i = 0; i < numPoints; ++i) { + const NumericType Gamma_HF = etchantFlux->at(i) * params.etchantFlux; + + if (!params.Config.usePhysisorption) { + // Single coverage: no chemisorption distinction + const NumericType denom = Gamma_HF + k_des; + theta_phys->at(i) = (denom > NumericType(1e-30)) + ? std::min(Gamma_HF / denom, NumericType(1)) + : NumericType(0); + theta_chem->at(i) = NumericType(0); + } else { + // Two-state: physisorbed + ion-activated chemisorbed + const NumericType k_act_ion = ionActFlux->at(i) * params.ionFlux * A_act; + const NumericType denom = Gamma_HF + k_des + k_act_ion + k_r_direct; + theta_phys->at(i) = (denom > NumericType(1e-30)) + ? std::min(Gamma_HF / denom, NumericType(1)) + : NumericType(0); + const NumericType num_chem = k_act_ion * theta_phys->at(i); + theta_chem->at(i) = (k_r > NumericType(1e-30)) + ? std::min(num_chem / k_r, + NumericType(1) - theta_phys->at(i)) + : NumericType(0); + } + } + } + SmartPointer> calculateVelocities(SmartPointer> rates, const std::vector> &coordinates, const std::vector &materialIds) override { + // Step 1: local coverage (no diffusion) — also serves as Dirichlet BC values updateCoverages(rates, materialIds); const auto numPoints = rates->getScalarData(0)->size(); - std::vector etchRate(numPoints, 0.); - const auto ionSputterFlux = rates->getScalarData("ionSputterFlux"); - const auto ionEnhancedFlux = rates->getScalarData("ionEnhancedFlux"); - const auto phi_HF = coverages->getScalarData("HFCoverage"); + const auto etchantFlux = rates->getScalarData("etchantFlux"); + const auto ionActFlux = rates->getScalarData("ionActivationFlux"); - std::vector *chRate = nullptr, *ieRate = nullptr, - *spRate = nullptr; + const NumericType k_des = params.effective_k_des(); + const NumericType k_r = params.effective_k_r(); + const NumericType k_r_direct = params.effective_k_r_direct(); + const NumericType A_act = params.IonActivation.A_act; + const NumericType Ds = params.D_s(); + + std::vector Gamma_HF_vec(numPoints); + std::vector k_act_ion_vec(numPoints); + for (size_t i = 0; i < numPoints; ++i) { + Gamma_HF_vec[i] = etchantFlux->at(i) * params.etchantFlux; + k_act_ion_vec[i] = ionActFlux->at(i) * params.ionFlux * A_act; + } + + auto theta_phys = coverages->getScalarData("physCoverage"); + auto theta_chem = coverages->getScalarData("chemCoverage"); + + // Step 2: build surface chain + sidewall flags + const auto chain = buildSurfaceChain(coordinates); + const auto sidewall = buildSidewallFlags(chain, coordinates); + + // Step 3: surface diffusion (only when enabled and D0 > 0) + if (params.Config.useSurfaceDiffusion && Ds > NumericType(0)) + applySurfaceDiffusion(*theta_phys, Gamma_HF_vec, k_act_ion_vec, + chain, coordinates, Ds, k_des, k_r_direct); + + // Step 4: recompute theta_chem + if (params.Config.usePhysisorption) { + // Two-state: zero on sidewalls (Lill 2023), steady-state elsewhere + for (size_t i = 0; i < numPoints; ++i) { + if (sidewall[i]) { + theta_chem->at(i) = NumericType(0); + continue; + } + const NumericType num_chem = k_act_ion_vec[i] * theta_phys->at(i); + theta_chem->at(i) = (k_r > NumericType(1e-30)) + ? std::min(num_chem / k_r, + NumericType(1) - theta_phys->at(i)) + : NumericType(0); + } + } else { + // Single coverage: no chemisorption + std::fill(theta_chem->begin(), theta_chem->end(), NumericType(0)); + } + + // Step 5: etch rate + const auto ionSputterFlux = rates->getScalarData("ionSputterFlux"); + const auto ionActivationFlux = rates->getScalarData("ionActivationFlux"); + const NumericType unitConversion = + units::Time::convertSecond() / units::Length::convertNanometer(); + + std::vector etchRate(numPoints, 0.); + + std::vector *chRate = nullptr, *spRate = nullptr; if (Logger::hasIntermediate()) { chRate = surfaceData->getScalarData("chemicalRate"); - ieRate = surfaceData->getScalarData("ionEnhancedRate"); spRate = surfaceData->getScalarData("sputterRate"); chRate->resize(numPoints); - ieRate->resize(numPoints); spRate->resize(numPoints); } - const NumericType kr = params.k_r(); - const NumericType unitConversion = - units::Time::convertSecond() / units::Length::convertNanometer(); - bool stop = false; - #pragma omp parallel for reduction(|| : stop) for (size_t i = 0; i < numPoints; ++i) { if (coordinates[i][D - 1] < params.etchStopDepth || stop) { stop = true; continue; } - if (!MaterialMap::isMaterial(materialIds[i], Material::SiO2)) continue; - const NumericType chemical = - kr * phi_HF->at(i); - const NumericType ionEnhanced = - phi_HF->at(i) * ionEnhancedFlux->at(i) * params.ionFlux; - const NumericType sputter = - ionSputterFlux->at(i) * params.ionFlux; + // Ion-enhanced etching: Y_ie * Gamma_ion * theta_phys (or single theta) + // Pure chemical: k_r_direct * theta_phys (only when usePhysisorption, + // otherwise the concept of "physisorbed" HF is not defined) + const NumericType ionEnhanced = + ionActivationFlux->at(i) * params.ionFlux * theta_phys->at(i); + const NumericType pureChemical = k_r_direct * theta_phys->at(i); + const NumericType sputter = ionSputterFlux->at(i) * params.ionFlux; etchRate[i] = -(1. / params.SiO2.rho) * - (chemical + ionEnhanced + sputter) * unitConversion; + (ionEnhanced + pureChemical + sputter) * + unitConversion; if (Logger::hasIntermediate()) { - chRate->at(i) = chemical; - ieRate->at(i) = ionEnhanced; + chRate->at(i) = ionEnhanced + pureChemical; spRate->at(i) = sputter; } } @@ -268,37 +355,181 @@ class HFCryoSurfaceModel : public SurfaceModel { return SmartPointer>::New(std::move(etchRate)); } - void updateCoverages(SmartPointer> rates, - const std::vector & /*materialIds*/) override { - const auto numPoints = rates->getScalarData(0)->size(); - const auto etchantFlux = rates->getScalarData("etchantFlux"); - const auto ionEnhancedFlux = rates->getScalarData("ionEnhancedFlux"); +private: + const HFCryoParameters ¶ms; - auto phi_HF = coverages->getScalarData("HFCoverage"); - phi_HF->resize(numPoints); + // Nearest-neighbor chain ordering surface points as a 1D curve. + // Starts from the topmost point (open to plasma). + std::vector + buildSurfaceChain(const std::vector> &points) const { + const size_t N = points.size(); + if (N == 0) return {}; + + std::vector visited(N, false); + std::vector chain; + chain.reserve(N); + + size_t start = 0; + for (size_t i = 1; i < N; ++i) + if (points[i][D - 1] > points[start][D - 1]) + start = i; + + chain.push_back(start); + visited[start] = true; + + NumericType thresholdSq = NumericType(9); + + for (size_t step = 1; step < N; ++step) { + const size_t curr = chain.back(); + NumericType bestDistSq = std::numeric_limits::max(); + size_t best = N; + + for (size_t j = 0; j < N; ++j) { + if (visited[j]) continue; + NumericType distSq = NumericType(0); + for (int d = 0; d < D; ++d) { + NumericType dd = points[j][d] - points[curr][d]; + distSq += dd * dd; + } + if (distSq < bestDistSq) { + bestDistSq = distSq; + best = j; + } + } - const NumericType k_des = params.k_des(); - const NumericType k_r = params.k_r(); + if (best == N || bestDistSq > thresholdSq) + break; -#pragma omp parallel for - for (size_t i = 0; i < numPoints; ++i) { - // Adsorption flux (scaled) - const NumericType Gamma_HF = etchantFlux->at(i) * params.etchantFlux; - // Ion-enhanced removal of adsorbed HF - const NumericType Y_ie_ion = - ionEnhancedFlux->at(i) * params.ionFlux; - - // Steady-state coverage: - // phi = Gamma_HF / (Gamma_HF + k_des + k_r + Y_ie_ion) - const NumericType denom = Gamma_HF + k_des + k_r + Y_ie_ion; - phi_HF->at(i) = (denom > NumericType(1e-30)) - ? std::min(Gamma_HF / denom, NumericType(1)) - : NumericType(0); + if (step == 1) + thresholdSq = bestDistSq * NumericType(4); + + chain.push_back(best); + visited[best] = true; } + + return chain; } -private: - const HFCryoParameters ¶ms; + // Classify each point in chain as sidewall (true) or not (false). + // Sidewall: local chain tangent is predominantly in height direction (D-1). + std::vector + buildSidewallFlags(const std::vector &chain, + const std::vector> &points) const { + const size_t numPoints = points.size(); + const size_t M = chain.size(); + std::vector sidewall(numPoints, false); + + for (size_t j = 0; j < M; ++j) { + // Centered tangent (or one-sided at endpoints) + Vec3D tang{0, 0, 0}; + if (j == 0) + for (int d = 0; d < D; ++d) + tang[d] = points[chain[1]][d] - points[chain[0]][d]; + else if (j == M - 1) + for (int d = 0; d < D; ++d) + tang[d] = points[chain[M-1]][d] - points[chain[M-2]][d]; + else + for (int d = 0; d < D; ++d) + tang[d] = points[chain[j+1]][d] - points[chain[j-1]][d]; + + NumericType len = NumericType(0); + for (int d = 0; d < D; ++d) len += tang[d] * tang[d]; + len = std::sqrt(len); + if (len < NumericType(1e-10)) continue; + + // Sidewall: tangent mostly vertical (height component > 0.7) + const NumericType vertFrac = std::abs(tang[D - 1]) / len; + sidewall[chain[j]] = (vertFrac > NumericType(0.7)); + } + + return sidewall; + } + + // Solve steady-state surface diffusion PDE on 1D chain via Thomas algorithm. + // + // Boundary conditions (Dirichlet): + // theta[0] = local Langmuir value at open-surface end (chain start) + // theta[M-1] = local Langmuir value at open-surface end (chain end) + // + // Interior (j = 1..M-2): + // -alpha*theta[j-1] + (2*alpha + k_des + Gamma[j] + k_act[j])*theta[j] + // - alpha*theta[j+1] = Gamma[j] + void applySurfaceDiffusion( + std::vector &theta_phys, + const std::vector &Gamma_HF_vec, + const std::vector &k_act_ion_vec, + const std::vector &chain, + const std::vector> &points, + NumericType Ds, NumericType k_des, NumericType k_r_direct) const { + const size_t M = chain.size(); + if (M < 3) return; + + // Average arc-length spacing + NumericType avgDs = NumericType(0); + for (size_t j = 0; j + 1 < M; ++j) { + NumericType distSq = NumericType(0); + for (int d = 0; d < D; ++d) { + NumericType dd = points[chain[j + 1]][d] - points[chain[j]][d]; + distSq += dd * dd; + } + avgDs += std::sqrt(distSq); + } + avgDs /= NumericType(M - 1); + if (avgDs < NumericType(1e-10)) return; + + const NumericType alpha = Ds / (avgDs * avgDs); + + // Save Dirichlet values at chain endpoints (local Langmuir equilibrium) + const NumericType bc_left = theta_phys[chain[0]]; + const NumericType bc_right = theta_phys[chain[M - 1]]; + + // Build tridiagonal system + std::vector a(M, -alpha); + std::vector b(M); + std::vector c(M, -alpha); + std::vector rhs(M); + + // Dirichlet BC at endpoints + a[0] = NumericType(0); b[0] = NumericType(1); + c[0] = NumericType(0); rhs[0] = bc_left; + a[M - 1] = NumericType(0); b[M - 1] = NumericType(1); + c[M - 1] = NumericType(0); rhs[M - 1] = bc_right; + + // Interior nodes: all theta_phys loss terms in diagonal + for (size_t j = 1; j + 1 < M; ++j) { + const size_t idx = chain[j]; + const NumericType G = Gamma_HF_vec[idx]; + const NumericType K = k_act_ion_vec[idx]; + b[j] = NumericType(2) * alpha + k_des + G + K + k_r_direct; + rhs[j] = G; + } + + // Thomas algorithm — forward sweep + std::vector c_prime(M), d_prime(M); + c_prime[0] = c[0] / b[0]; + d_prime[0] = rhs[0] / b[0]; + + for (size_t j = 1; j < M; ++j) { + const NumericType denom = b[j] - a[j] * c_prime[j - 1]; + const NumericType inv = (std::abs(denom) > NumericType(1e-30)) + ? NumericType(1) / denom + : NumericType(0); + c_prime[j] = c[j] * inv; + d_prime[j] = (rhs[j] - a[j] * d_prime[j - 1]) * inv; + } + + // Back substitution + theta_phys[chain[M - 1]] = + std::max(NumericType(0), std::min(d_prime[M - 1], NumericType(1))); + + for (size_t j = M - 2; ; --j) { + theta_phys[chain[j]] = + std::max(NumericType(0), + std::min(d_prime[j] - c_prime[j] * theta_phys[chain[j + 1]], + NumericType(1))); + if (j == 0) break; + } + } }; } // namespace impl @@ -309,26 +540,6 @@ class HFCryoEtching final : public ProcessModelCPU { public: HFCryoEtching() { initializeModel(); } - // Convenience constructor - HFCryoEtching(NumericType ionFlux, NumericType etchantFlux, - NumericType temperature, // K - NumericType meanEnergy, // eV - NumericType sigmaEnergy, // eV - NumericType E_des = 0.25, // eV (Frenkel desorption) - NumericType E_a = 0.10, // eV (reaction activation) - NumericType etchStopDepth = - std::numeric_limits::lowest()) { - params.ionFlux = ionFlux; - params.etchantFlux = etchantFlux; - params.temperature = temperature; - params.Ions.meanEnergy = meanEnergy; - params.Ions.sigmaEnergy = sigmaEnergy; - params.Desorption.E_des = E_des; - params.Reaction.E_a = E_a; - params.etchStopDepth = etchStopDepth; - initializeModel(); - } - explicit HFCryoEtching(const HFCryoParameters &pParams) : params(pParams) { initializeModel(); @@ -348,10 +559,8 @@ class HFCryoEtching final : public ProcessModelCPU { VIENNACORE_LOG_ERROR("Units have not been set."); } - auto ion = std::make_unique>(params); - auto etchant = - std::make_unique>(params); - + auto ion = std::make_unique>(params); + auto etchant = std::make_unique>(params); auto surfModel = SmartPointer>::New(params); auto velField = diff --git a/include/viennaps/models/psHFCryoParameters.hpp b/include/viennaps/models/psHFCryoParameters.hpp index 08c7d179..0df96c45 100644 --- a/include/viennaps/models/psHFCryoParameters.hpp +++ b/include/viennaps/models/psHFCryoParameters.hpp @@ -25,11 +25,43 @@ template struct HFCryoParameters { } Desorption; // Arrhenius reaction rate: k_r = A_r * exp(-E_a / kB*T) + // Chemisorbed HF reacting with SiO2 (ion-assisted pathway) struct ReactionType { NumericType A_r = 3.0e2; // pre-exponential factor (1e15 cm⁻²s⁻¹) NumericType E_a = 0.10; // reaction activation energy (eV) } Reaction; + // Direct (pure) chemical etching: physisorbed HF reacts with SiO2 thermally, + // no ion activation needed. Higher barrier than chemisorption pathway. + // Negligible at cryo (150K~200K), significant at room temperature (300K). + struct DirectReactionType { + NumericType A_r = 1.0e4; // pre-exponential factor + NumericType E_a = 0.25; // activation energy (eV) > chemisorption E_a + } DirectReaction; + + // Ion activation: ions convert physisorbed HF -> chemisorbed HF + // k_act_ion = A_act * ionEnhancedFlux * ionFlux [same units as Gamma_HF] + struct IonActivationType { + NumericType A_act = 1.0; // dimensionless scaling factor + } IonActivation; + + // Model configuration: select which physical features are active. + // Enables progressive comparison: none → +temp → +phys → +diffusion + struct ModelConfigType { + bool useTemperatureDependence = true; // Frenkel-Arrhenius/Arrhenius for rates + bool usePhysisorption = true; // two-state model (theta_phys + theta_chem) + bool useSurfaceDiffusion = true; // 1D diffusion PDE along surface chain + NumericType T_ref = 300.; // fixed reference T when !useTemperatureDependence + } Config; + + // Surface diffusion of physisorbed HF along the surface + // D_s(T) = D0 * exp(-omega * E_des / kB*T) + // E_diff = omega * E_des (corrugation ratio from Lill 2023, Table IV: omega~0.2-0.3) + struct DiffusionType { + NumericType D0 = 1.0e3; // nm²/s, pre-exponential + NumericType omega = 0.25; // corrugation ratio + } Diffusion; + // SiO2 material properties struct SiO2Type { NumericType rho = constants::SiO2::rho; // 1e22 atoms/cm³ @@ -49,16 +81,33 @@ template struct HFCryoParameters { NumericType minAngle = constants::Ion::minAngle; } Ions; - // Compute k_des(T) in units consistent with fluxes (1/s) + // Raw Arrhenius/Frenkel-Arrhenius at actual temperature NumericType k_des() const { - return Desorption.nu0 * - std::exp(-Desorption.E_des / (constants::kB * temperature)); + return Desorption.nu0 * std::exp(-Desorption.E_des / (constants::kB * temperature)); } - - // Compute k_r(T) in units of (1e15 cm⁻²s⁻¹) NumericType k_r() const { - return Reaction.A_r * - std::exp(-Reaction.E_a / (constants::kB * temperature)); + return Reaction.A_r * std::exp(-Reaction.E_a / (constants::kB * temperature)); + } + NumericType k_r_direct() const { + return DirectReaction.A_r * std::exp(-DirectReaction.E_a / (constants::kB * temperature)); + } + NumericType D_s() const { + const NumericType E_diff = Diffusion.omega * Desorption.E_des; + return Diffusion.D0 * std::exp(-E_diff / (constants::kB * temperature)); + } + + // Config-aware rates: use T_ref when !useTemperatureDependence + NumericType effective_k_des() const { + const NumericType T = Config.useTemperatureDependence ? temperature : Config.T_ref; + return Desorption.nu0 * std::exp(-Desorption.E_des / (constants::kB * T)); + } + NumericType effective_k_r() const { + const NumericType T = Config.useTemperatureDependence ? temperature : Config.T_ref; + return Reaction.A_r * std::exp(-Reaction.E_a / (constants::kB * T)); + } + NumericType effective_k_r_direct() const { + const NumericType T = Config.useTemperatureDependence ? temperature : Config.T_ref; + return DirectReaction.A_r * std::exp(-DirectReaction.E_a / (constants::kB * T)); } }; From c7f9af268315af8e81fbb1034a7450c8b7609f73 Mon Sep 17 00:00:00 2001 From: monst Date: Tue, 19 May 2026 11:52:03 +0900 Subject: [PATCH 4/9] refactor: simplify HFCryo4Cases to single full-model run at 4 temperatures Replace 4-case comparison example with a cleaner single simulation that enables all three physical features (Arrhenius temperature dependence, physisorption, surface diffusion) and sweeps 150/200/250/300 K. Outputs HFCryo_T{T}K.vtp files colored by Temperature_K for ParaView. Co-Authored-By: Claude Sonnet 4.5 --- examples/HFCryo4Cases/HFCryo4Cases.cpp | 155 ++++++++----------------- 1 file changed, 47 insertions(+), 108 deletions(-) diff --git a/examples/HFCryo4Cases/HFCryo4Cases.cpp b/examples/HFCryo4Cases/HFCryo4Cases.cpp index 4de95500..3c21e559 100644 --- a/examples/HFCryo4Cases/HFCryo4Cases.cpp +++ b/examples/HFCryo4Cases/HFCryo4Cases.cpp @@ -16,41 +16,59 @@ using namespace viennaps; static constexpr double gridDelta = 1.0; static constexpr double xExtent = 80.0; -static constexpr double yExtent = 2000.0; // 700 → 2000 (깊은 트렌치 수용) +static constexpr double yExtent = 2000.0; static constexpr double trenchWidth = 10.0; -static constexpr double maskHeight = 80.0; // 50 → 80 (두꺼운 마스크) -static constexpr double processTime = 400.0; // 90 → 400 (긴 공정 시간) - -// Save mesh with Temperature_K and CaseID as point scalars. -// In ParaView: Color by "Temperature_K" → temperature gradient (150-300K) -// Color by "CaseID" → 4 distinct case colors (0-3) -static void saveColoredMesh(SmartPointer> domain, - const std::string &filename, - double T, int caseId) { +static constexpr double maskHeight = 80.0; +static constexpr double processTime = 400.0; + +static void saveMesh(SmartPointer> domain, + const std::string &filename, double T) { auto mesh = viennals::Mesh::New(); viennals::ToMultiSurfaceMesh(domain->getLevelSets(), mesh).apply(); const size_t n = mesh->nodes.size(); mesh->getPointData().insertNextScalarData( - std::vector(n, T), "Temperature_K"); - mesh->getPointData().insertNextScalarData( - std::vector(n, double(caseId)), "CaseID"); + std::vector(n, T), "Temperature_K"); viennals::VTKWriter(mesh, filename).apply(); } -static void runCase(const std::string &filename, - HFCryoParameters params, - int caseId) { +static void runCase(double T) { + std::cout << " T=" << std::setw(3) << (int)T << "K ... "; + std::cout.flush(); + auto domain = SmartPointer>::New(gridDelta, xExtent, yExtent); MakeTrench(domain, trenchWidth, 0.0, 0.0, - maskHeight, 0.0, false, Material::SiO2).apply(); + maskHeight, 0.0, false, Material::SiO2).apply(); + + // ── Parameters: all 3 features ON ─────────────────────────────────────────── + HFCryoParameters params; + params.temperature = T; + params.ionFlux = 3.0; + params.etchantFlux = 5.0e3; + params.gamma_HF = 0.9; + params.Desorption.nu0 = 1.0e8; + params.Desorption.E_des = 0.20; + params.Reaction.A_r = 3.0e2; + params.Reaction.E_a = 0.10; + params.DirectReaction.A_r = 1.0e4; + params.DirectReaction.E_a = 0.25; + params.IonActivation.A_act = 1.0; + params.Diffusion.D0 = 1.0e3; // surface diffusion ON + params.Diffusion.omega = 0.25; + params.Ions.meanEnergy = 150.0; + params.Ions.sigmaEnergy = 10.0; + params.Ions.exponent = 600.0; + params.Config.useTemperatureDependence = true; // feature 1: Arrhenius + params.Config.usePhysisorption = true; // feature 2: physisorption + params.Config.useSurfaceDiffusion = true; // feature 3: surface diffusion + params.Config.T_ref = 300.0; auto model = SmartPointer>::New(params); Process process(domain, model, processTime); RayTracingParameters rayParams; - rayParams.raysPerPoint = 5000; // 3000 → 5000 (깊은 트렌치 정확도↑) + rayParams.raysPerPoint = 5000; process.setParameters(rayParams); AdvectionParameters advParams; @@ -60,8 +78,9 @@ static void runCase(const std::string &filename, process.apply(); - saveColoredMesh(domain, filename + ".vtp", params.temperature, caseId); - std::cout << "done -> " << filename << ".vtp\n"; + const std::string fname = "HFCryo_T" + std::to_string((int)T) + "K.vtp"; + saveMesh(domain, fname, T); + std::cout << "done -> " << fname << "\n"; } int main() { @@ -71,97 +90,17 @@ int main() { units::Length::setUnit(units::Length::NANOMETER); units::Time::setUnit(units::Time::SECOND); - // ── Base parameters ───────────────────────────────────────────────────────── - HFCryoParameters base; - base.ionFlux = 3.0; // 1.0 → 3.0 (이온 강화 → 수직 식각↑) - base.etchantFlux = 5.0e3; // 1e3 → 5e3 (HF 공급↑) - base.gamma_HF = 0.9; - base.Desorption.nu0 = 1.0e8; - base.Desorption.E_des = 0.20; - base.Reaction.A_r = 3.0e2; - base.Reaction.E_a = 0.10; - base.DirectReaction.A_r = 1.0e4; - base.DirectReaction.E_a = 0.25; - base.IonActivation.A_act = 1.0; - base.Diffusion.D0 = 0.0; - base.Diffusion.omega = 0.25; - base.Ions.meanEnergy = 150.0; // 100 → 150 eV (더 높은 이온 에너지) - base.Ions.sigmaEnergy = 10.0; - base.Ions.exponent = 600.0; // 300 → 600 (이온빔 더 수직, 이방성↑) - base.Config.T_ref = 300.0; - - // ── Case definitions ───────────────────────────────────────────────────────── - struct CaseInfo { - std::string tag; - std::string desc; - bool tempDep; - bool phys; - bool diff; - double D0; - }; - - const CaseInfo cases[4] = { - {"C0_none", "nothing (constant k_des, single theta, no diffusion)", - false, false, false, 0.0}, - {"C1_temp", "+Arrhenius k_des (single theta, no diffusion)", - true, false, false, 0.0}, - {"C2_phys", "+physisorption (theta_phys+theta_chem, no diffusion)", - true, true, false, 0.0}, - {"C3_diff", "+surface diffusion (D0=1e3)", - true, true, true, 1.0e3}, - }; + std::cout << "=== HF Cryo Etching: full model (Arrhenius + physisorption + surface diffusion) ===\n"; + std::cout << "Geometry: " << trenchWidth << "nm trench, " + << maskHeight << "nm mask, " << processTime << "s process\n\n"; const double Tlist[4] = {150., 200., 250., 300.}; - // ── Print rate table ───────────────────────────────────────────────────────── - std::cout << "=== HF Cryo 4-Case x 4-Temperature Comparison ===\n"; - std::cout << "Geometry: " << trenchWidth << "nm trench, " - << processTime << "s process\n\n"; - - std::cout << "Effective k_des per case/temperature:\n"; - std::cout << std::setw(20) << "" << std::setw(10) << "150K" - << std::setw(10) << "200K" << std::setw(10) << "250K" - << std::setw(10) << "300K\n"; - std::cout << std::string(60, '-') << "\n"; - for (const auto &ci : cases) { - std::cout << std::setw(20) << ci.tag; - for (double T : Tlist) { - auto p = base; - p.temperature = T; - p.Config.useTemperatureDependence = ci.tempDep; - std::cout << std::setw(10) << std::setprecision(1) << std::fixed - << p.effective_k_des(); - } - std::cout << "\n"; - } - std::cout << "\n"; - - // ── Run all 16 simulations ─────────────────────────────────────────────────── - for (int c = 0; c < 4; ++c) { - const auto &ci = cases[c]; - std::cout << "=== Case " << c << ": " << ci.desc << " ===\n"; - - for (double T : Tlist) { - std::cout << " T=" << std::setw(3) << (int)T << "K ... "; - std::cout.flush(); - - auto p = base; - p.temperature = T; - p.Config.useTemperatureDependence = ci.tempDep; - p.Config.usePhysisorption = ci.phys; - p.Config.useSurfaceDiffusion = ci.diff; - p.Diffusion.D0 = ci.D0; - - std::string fname = ci.tag + "_T" + std::to_string((int)T) + "K"; - runCase(fname, p, c); - } - std::cout << "\n"; - } - - std::cout << "=== Done: 16 VTP files generated ===\n\n"; - std::cout << "ParaView: open all *.vtp files\n"; - std::cout << " Color by 'Temperature_K' -> blue(150K) to red(300K)\n"; - std::cout << " Color by 'CaseID' -> 4 discrete case colors\n"; + for (double T : Tlist) + runCase(T); + + std::cout << "\n=== Done: 4 VTP files generated ===\n"; + std::cout << "ParaView: open HFCryo_T*.vtp, Color by 'Temperature_K'\n"; return 0; } From 64c16829c4558a9f24ead0cbc68747b14b0d81ac Mon Sep 17 00:00:00 2001 From: monst Date: Tue, 19 May 2026 16:45:05 +0900 Subject: [PATCH 5/9] feat: add surface diffusion OFF/ON comparison to HFCryo4Cases Run all four temperatures (150-300 K) with surface diffusion both disabled and enabled, and report the resulting trench etch depth for each case. VTP outputs are tagged HFCryo_diffOFF/diffON_T*K.vtp. Co-Authored-By: Claude Opus 4.7 (1M context) --- examples/HFCryo4Cases/HFCryo4Cases.cpp | 80 ++++++++++++++++++-------- 1 file changed, 55 insertions(+), 25 deletions(-) diff --git a/examples/HFCryo4Cases/HFCryo4Cases.cpp b/examples/HFCryo4Cases/HFCryo4Cases.cpp index 3c21e559..f5a83670 100644 --- a/examples/HFCryo4Cases/HFCryo4Cases.cpp +++ b/examples/HFCryo4Cases/HFCryo4Cases.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -21,27 +22,29 @@ static constexpr double trenchWidth = 10.0; static constexpr double maskHeight = 80.0; static constexpr double processTime = 400.0; -static void saveMesh(SmartPointer> domain, - const std::string &filename, double T) { +// Writes the surface mesh and returns the etch depth: how far the deepest +// point of the etch front sits below the original substrate surface (y = 0). +static double saveMeshAndDepth(SmartPointer> domain, + const std::string &filename, double T) { auto mesh = viennals::Mesh::New(); viennals::ToMultiSurfaceMesh(domain->getLevelSets(), mesh).apply(); const size_t n = mesh->nodes.size(); mesh->getPointData().insertNextScalarData( std::vector(n, T), "Temperature_K"); - viennals::VTKWriter(mesh, filename).apply(); -} -static void runCase(double T) { - std::cout << " T=" << std::setw(3) << (int)T << "K ... "; - std::cout.flush(); + double minY = 0.0; + for (size_t i = 0; i < n; ++i) + minY = std::min(minY, mesh->nodes[i][1]); + return -minY; +} +static double runCase(double T, bool useDiffusion) { auto domain = SmartPointer>::New(gridDelta, xExtent, yExtent); MakeTrench(domain, trenchWidth, 0.0, 0.0, maskHeight, 0.0, false, Material::SiO2).apply(); - // ── Parameters: all 3 features ON ─────────────────────────────────────────── HFCryoParameters params; params.temperature = T; params.ionFlux = 3.0; @@ -54,14 +57,14 @@ static void runCase(double T) { params.DirectReaction.A_r = 1.0e4; params.DirectReaction.E_a = 0.25; params.IonActivation.A_act = 1.0; - params.Diffusion.D0 = 1.0e3; // surface diffusion ON + params.Diffusion.D0 = 1.0e3; params.Diffusion.omega = 0.25; params.Ions.meanEnergy = 150.0; params.Ions.sigmaEnergy = 10.0; params.Ions.exponent = 600.0; - params.Config.useTemperatureDependence = true; // feature 1: Arrhenius - params.Config.usePhysisorption = true; // feature 2: physisorption - params.Config.useSurfaceDiffusion = true; // feature 3: surface diffusion + params.Config.useTemperatureDependence = true; + params.Config.usePhysisorption = true; + params.Config.useSurfaceDiffusion = useDiffusion; // toggled per run params.Config.T_ref = 300.0; auto model = SmartPointer>::New(params); @@ -78,9 +81,10 @@ static void runCase(double T) { process.apply(); - const std::string fname = "HFCryo_T" + std::to_string((int)T) + "K.vtp"; - saveMesh(domain, fname, T); - std::cout << "done -> " << fname << "\n"; + const std::string tag = useDiffusion ? "diffON" : "diffOFF"; + const std::string fname = + "HFCryo_" + tag + "_T" + std::to_string((int)T) + "K.vtp"; + return saveMeshAndDepth(domain, fname, T); } int main() { @@ -90,17 +94,43 @@ int main() { units::Length::setUnit(units::Length::NANOMETER); units::Time::setUnit(units::Time::SECOND); - std::cout << "=== HF Cryo Etching: full model (Arrhenius + physisorption + surface diffusion) ===\n"; - std::cout << "Geometry: " << trenchWidth << "nm trench, " - << maskHeight << "nm mask, " << processTime << "s process\n\n"; - const double Tlist[4] = {150., 200., 250., 300.}; - for (double T : Tlist) - runCase(T); - - std::cout << "\n=== Done: 4 VTP files generated ===\n"; - std::cout << "ParaView: open HFCryo_T*.vtp, Color by 'Temperature_K'\n"; - + std::cout << "=== HF Cryo Etching: surface diffusion OFF vs ON ===\n"; + std::cout << "Geometry: " << trenchWidth << "nm trench, " << maskHeight + << "nm mask, " << processTime << "s process\n\n"; + + double depthOff[4], depthOn[4]; + + std::cout << "--- Surface diffusion OFF ---\n"; + for (int i = 0; i < 4; ++i) { + std::cout << " T=" << std::setw(3) << (int)Tlist[i] << "K ... " + << std::flush; + depthOff[i] = runCase(Tlist[i], false); + std::cout << "depth = " << std::fixed << std::setprecision(1) + << depthOff[i] << " nm\n"; + } + + std::cout << "\n--- Surface diffusion ON ---\n"; + for (int i = 0; i < 4; ++i) { + std::cout << " T=" << std::setw(3) << (int)Tlist[i] << "K ... " + << std::flush; + depthOn[i] = runCase(Tlist[i], true); + std::cout << "depth = " << std::fixed << std::setprecision(1) + << depthOn[i] << " nm\n"; + } + + std::cout << "\n=== Summary: trench etch depth (nm) ===\n"; + std::cout << " T[K] diffOFF diffON gain\n"; + for (int i = 0; i < 4; ++i) { + const double gain = + depthOff[i] > 1e-9 ? (depthOn[i] / depthOff[i] - 1.0) * 100.0 : 0.0; + std::cout << " " << std::setw(4) << (int)Tlist[i] << " " << std::setw(8) + << std::setprecision(1) << depthOff[i] << " " << std::setw(8) + << depthOn[i] << " " << std::setw(7) << std::setprecision(1) + << gain << "%\n"; + } + + std::cout << "\n=== Done: 8 VTP files (HFCryo_diffOFF/ON_T*.vtp) ===\n"; return 0; } From bf3f42dabc22e4800d6dc9328cf5b0c18c95e84f Mon Sep 17 00:00:00 2001 From: monst Date: Tue, 19 May 2026 22:27:58 +0900 Subject: [PATCH 6/9] fix: surface chain in HFCryoEtching now traverses the full trench contour buildSurfaceChain started the nearest-neighbor walk from the topmost point (interior of the flat mask top) with a tight gap threshold, so the walk covered only one mask-top strip and never descended into the trench. Surface diffusion was therefore applied to the wrong part of the surface and had no effect on the etch front. Start from a true endpoint of the open surface (leftmost point) and widen the gap threshold so the chain follows mask top -> sidewall -> bottom -> sidewall -> mask top. Co-Authored-By: Claude Opus 4.7 (1M context) --- include/viennaps/models/psHFCryoEtching.hpp | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/include/viennaps/models/psHFCryoEtching.hpp b/include/viennaps/models/psHFCryoEtching.hpp index 5332a6ca..cf2ac840 100644 --- a/include/viennaps/models/psHFCryoEtching.hpp +++ b/include/viennaps/models/psHFCryoEtching.hpp @@ -359,7 +359,11 @@ class HFCryoSurfaceModel : public SurfaceModel { const HFCryoParameters ¶ms; // Nearest-neighbor chain ordering surface points as a 1D curve. - // Starts from the topmost point (open to plasma). + // Starts from a true endpoint of the open surface (leftmost point) so the + // greedy walk traverses the whole contour: mask top -> trench sidewall -> + // bottom -> opposite sidewall -> mask top. Starting in the interior (e.g. + // the topmost point) makes the walk cover only one flat strip and miss the + // trench entirely. std::vector buildSurfaceChain(const std::vector> &points) const { const size_t N = points.size(); @@ -371,13 +375,15 @@ class HFCryoSurfaceModel : public SurfaceModel { size_t start = 0; for (size_t i = 1; i < N; ++i) - if (points[i][D - 1] > points[start][D - 1]) + if (points[i][0] < points[start][0] || + (points[i][0] == points[start][0] && + points[i][D - 1] > points[start][D - 1])) start = i; chain.push_back(start); visited[start] = true; - NumericType thresholdSq = NumericType(9); + NumericType thresholdSq = std::numeric_limits::max(); for (size_t step = 1; step < N; ++step) { const size_t curr = chain.back(); @@ -400,8 +406,11 @@ class HFCryoSurfaceModel : public SurfaceModel { if (best == N || bestDistSq > thresholdSq) break; + // Calibrate gap threshold from the first (typical) spacing. Generous + // factor so corners and into-trench turns survive, but a jump across + // the open trench mouth still breaks the chain. if (step == 1) - thresholdSq = bestDistSq * NumericType(4); + thresholdSq = bestDistSq * NumericType(36); chain.push_back(best); visited[best] = true; From e88403158bc5a7a55536097eabc4e0442752fa8c Mon Sep 17 00:00:00 2001 From: monst Date: Tue, 19 May 2026 23:35:04 +0900 Subject: [PATCH 7/9] fix: replace HFCryoEtching surface diffusion with a stable mesh-free solver The surface-diffusion step solved a steady-state 1D PDE on a nearest-neighbor surface chain with Dirichlet boundaries at the chain endpoints. Chain-ordering defects turned the discrete Laplacian into a non-smoothing operator that injected geometry-reconstruction noise; at large diffusion coefficients this destabilized the etch front and a D0 scan produced chaotic, non-monotonic etch depths (-53% to +39% vs diffusion-off). Replace it with a mesh-free reaction-diffusion solve on the surface point cloud, iterated with Jacobi. Each update is a positive, diagonally dominant contraction, so it is unconditionally stable and strictly smoothing, with no chain reconstruction and no imposed boundary values; open mask-top points act as natural diffusion sources. The D0 scan is now monotonic and saturating, and surface diffusion produces a stable ARDE-mitigation effect at physical D0. Co-Authored-By: Claude Opus 4.7 (1M context) --- include/viennaps/models/psHFCryoEtching.hpp | 149 +++++++++++--------- 1 file changed, 79 insertions(+), 70 deletions(-) diff --git a/include/viennaps/models/psHFCryoEtching.hpp b/include/viennaps/models/psHFCryoEtching.hpp index cf2ac840..e811edd2 100644 --- a/include/viennaps/models/psHFCryoEtching.hpp +++ b/include/viennaps/models/psHFCryoEtching.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include #include @@ -282,7 +283,7 @@ class HFCryoSurfaceModel : public SurfaceModel { // Step 3: surface diffusion (only when enabled and D0 > 0) if (params.Config.useSurfaceDiffusion && Ds > NumericType(0)) applySurfaceDiffusion(*theta_phys, Gamma_HF_vec, k_act_ion_vec, - chain, coordinates, Ds, k_des, k_r_direct); + coordinates, Ds, k_des, k_r_direct); // Step 4: recompute theta_chem if (params.Config.usePhysisorption) { @@ -454,89 +455,97 @@ class HFCryoSurfaceModel : public SurfaceModel { return sidewall; } - // Solve steady-state surface diffusion PDE on 1D chain via Thomas algorithm. + // Surface diffusion of physisorbed HF, solved as a steady-state reaction- + // diffusion balance directly on the surface point cloud (mesh-free): // - // Boundary conditions (Dirichlet): - // theta[0] = local Langmuir value at open-surface end (chain start) - // theta[M-1] = local Langmuir value at open-surface end (chain end) + // Ds * Lap(theta) + Gamma*(1 - theta) + // - (k_des + k_act_ion + k_r_direct) * theta = 0 // - // Interior (j = 1..M-2): - // -alpha*theta[j-1] + (2*alpha + k_des + Gamma[j] + k_act[j])*theta[j] - // - alpha*theta[j+1] = Gamma[j] + // The Laplacian Lap(theta_i) = sum_j w_ij (theta_j - theta_i) is taken over + // neighbouring surface points (w_ij = 1/d_ij^2). The balance is solved by + // Jacobi iteration; each update is a weighted average of the local Langmuir + // value and the neighbour average, so it is a positive, diagonally dominant + // contraction — unconditionally stable and strictly smoothing. Open mask-top + // points keep a high local coverage and act as natural diffusion sources, so + // no reconstructed 1D chain and no imposed boundary values are needed. The + // earlier chain + Dirichlet formulation injected geometry-reconstruction + // noise that destabilised the etch front at large Ds. void applySurfaceDiffusion( std::vector &theta_phys, const std::vector &Gamma_HF_vec, const std::vector &k_act_ion_vec, - const std::vector &chain, const std::vector> &points, NumericType Ds, NumericType k_des, NumericType k_r_direct) const { - const size_t M = chain.size(); - if (M < 3) return; + const size_t N = points.size(); + if (N < 3) + return; - // Average arc-length spacing - NumericType avgDs = NumericType(0); - for (size_t j = 0; j + 1 < M; ++j) { - NumericType distSq = NumericType(0); + auto distSq = [&](size_t i, size_t j) { + NumericType s = NumericType(0); for (int d = 0; d < D; ++d) { - NumericType dd = points[chain[j + 1]][d] - points[chain[j]][d]; - distSq += dd * dd; + const NumericType dd = points[i][d] - points[j][d]; + s += dd * dd; } - avgDs += std::sqrt(distSq); - } - avgDs /= NumericType(M - 1); - if (avgDs < NumericType(1e-10)) return; - - const NumericType alpha = Ds / (avgDs * avgDs); - - // Save Dirichlet values at chain endpoints (local Langmuir equilibrium) - const NumericType bc_left = theta_phys[chain[0]]; - const NumericType bc_right = theta_phys[chain[M - 1]]; - - // Build tridiagonal system - std::vector a(M, -alpha); - std::vector b(M); - std::vector c(M, -alpha); - std::vector rhs(M); - - // Dirichlet BC at endpoints - a[0] = NumericType(0); b[0] = NumericType(1); - c[0] = NumericType(0); rhs[0] = bc_left; - a[M - 1] = NumericType(0); b[M - 1] = NumericType(1); - c[M - 1] = NumericType(0); rhs[M - 1] = bc_right; - - // Interior nodes: all theta_phys loss terms in diagonal - for (size_t j = 1; j + 1 < M; ++j) { - const size_t idx = chain[j]; - const NumericType G = Gamma_HF_vec[idx]; - const NumericType K = k_act_ion_vec[idx]; - b[j] = NumericType(2) * alpha + k_des + G + K + k_r_direct; - rhs[j] = G; - } + return s; + }; - // Thomas algorithm — forward sweep - std::vector c_prime(M), d_prime(M); - c_prime[0] = c[0] / b[0]; - d_prime[0] = rhs[0] / b[0]; - - for (size_t j = 1; j < M; ++j) { - const NumericType denom = b[j] - a[j] * c_prime[j - 1]; - const NumericType inv = (std::abs(denom) > NumericType(1e-30)) - ? NumericType(1) / denom - : NumericType(0); - c_prime[j] = c[j] * inv; - d_prime[j] = (rhs[j] - a[j] * d_prime[j - 1]) * inv; - } + // Neighbour radius from a robust estimate of the point spacing. + std::vector nnSq(N, std::numeric_limits::max()); +#pragma omp parallel for + for (size_t i = 0; i < N; ++i) + for (size_t j = 0; j < N; ++j) + if (j != i) + nnSq[i] = std::min(nnSq[i], distSq(i, j)); + std::vector sortedNN(nnSq); + std::sort(sortedNN.begin(), sortedNN.end()); + const NumericType medianNNSq = sortedNN[N / 2]; + if (medianNNSq < NumericType(1e-20)) + return; + const NumericType radiusSq = medianNNSq * NumericType(9); // (3x spacing)^2 - // Back substitution - theta_phys[chain[M - 1]] = - std::max(NumericType(0), std::min(d_prime[M - 1], NumericType(1))); + // Neighbour lists with inverse-square-distance weights. + std::vector> nbrIdx(N); + std::vector> nbrW(N); + std::vector wSum(N, NumericType(0)); +#pragma omp parallel for + for (size_t i = 0; i < N; ++i) { + for (size_t j = 0; j < N; ++j) { + if (j == i) + continue; + const NumericType d2 = distSq(i, j); + if (d2 > radiusSq || d2 < NumericType(1e-20)) + continue; + const NumericType w = NumericType(1) / d2; + nbrIdx[i].push_back(j); + nbrW[i].push_back(w); + wSum[i] += w; + } + } - for (size_t j = M - 2; ; --j) { - theta_phys[chain[j]] = - std::max(NumericType(0), - std::min(d_prime[j] - c_prime[j] * theta_phys[chain[j + 1]], - NumericType(1))); - if (j == 0) break; + // Jacobi iteration of the steady-state balance. The update is a contraction + // (all coefficients positive, diagonally dominant) so it always converges. + std::vector next(theta_phys); + constexpr int maxIter = 1000; + for (int it = 0; it < maxIter; ++it) { +#pragma omp parallel for + for (size_t i = 0; i < N; ++i) { + const NumericType G = Gamma_HF_vec[i]; + const NumericType losses = k_des + k_act_ion_vec[i] + k_r_direct; + NumericType nbrSum = NumericType(0); + for (size_t n = 0; n < nbrIdx[i].size(); ++n) + nbrSum += nbrW[i][n] * theta_phys[nbrIdx[i][n]]; + const NumericType denom = G + losses + Ds * wSum[i]; + const NumericType val = (denom > NumericType(1e-30)) + ? (G + Ds * nbrSum) / denom + : theta_phys[i]; + next[i] = std::max(NumericType(0), std::min(val, NumericType(1))); + } + NumericType maxChange = NumericType(0); + for (size_t i = 0; i < N; ++i) + maxChange = std::max(maxChange, std::abs(next[i] - theta_phys[i])); + theta_phys.swap(next); + if (maxChange < NumericType(1e-5)) + break; } } }; From 2f4736e5c151eccaf0ffff41e5ad48a2eb366384 Mon Sep 17 00:00:00 2001 From: monst Date: Wed, 20 May 2026 01:43:46 +0900 Subject: [PATCH 8/9] docs: correct HFCryoEtching model header to match the implementation The header comment described an etch rate of k_r*theta_chem and a 1D-chain diffusion PDE, neither of which the code implements. Update it to the actual model: mesh-free surface diffusion, the clamped chemisorbed coverage, theta_chem's role as site-blocking in the etchant sticking (not a separate etch term), and the implemented etch rate (Y_ie*Gamma_ion*theta_phys + k_r_direct*theta_phys + Y_sp*Gamma_ion). Co-Authored-By: Claude Opus 4.7 (1M context) --- include/viennaps/models/psHFCryoEtching.hpp | 25 ++++++++++++--------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/include/viennaps/models/psHFCryoEtching.hpp b/include/viennaps/models/psHFCryoEtching.hpp index e811edd2..f1072f3f 100644 --- a/include/viennaps/models/psHFCryoEtching.hpp +++ b/include/viennaps/models/psHFCryoEtching.hpp @@ -22,19 +22,22 @@ using namespace viennacore; // // Two-state adsorption + surface diffusion model: // -// Physisorbed HF (theta_phys) — solved via 1D diffusion PDE along surface: -// D_s * d²theta/ds² + Gamma_HF*(1-theta) - k_des*theta -// - k_act*Gamma_ion*theta = 0 -// BC: Dirichlet at chain ends (open surface = local Langmuir equilibrium) +// Physisorbed HF (theta_phys): local Langmuir balance of adsorption, +// desorption and reaction, optionally smoothed by surface diffusion. The +// diffusion is a steady-state reaction-diffusion balance solved mesh-free +// on the surface point cloud (see applySurfaceDiffusion). // -// Chemisorbed HF (theta_chem) — ion-activated, zero on sidewalls (Lill 2023): -// theta_chem = k_act * Gamma_ion * theta_phys / k_r (bottom/non-sidewall) -// theta_chem = 0 (sidewall, per paper) +// Chemisorbed HF (theta_chem): ion-activated, zero on sidewalls (Lill 2023): +// theta_chem = min(k_act * Gamma_ion * theta_phys / k_r, 1 - theta_phys) +// theta_chem = 0 on sidewalls (no ion bombardment, per paper) +// theta_chem is NOT a separate etch term. It occupies surface sites and so +// lowers the HF sticking probability in the etchant ray tracing: +// S_eff = gamma_HF * (1 - theta_phys - theta_chem). // -// Etch rate: -// v = -(Y_ie * Gamma_ion * theta_phys -// + k_r * theta_chem -// + Y_sp * Gamma_ion) / rho +// Etch rate (as implemented in calculateVelocities): +// v = -(Y_ie * Gamma_ion * theta_phys // ion-enhanced etching +// + k_r_direct * theta_phys // direct thermal reaction +// + Y_sp * Gamma_ion) / rho // physical sputtering namespace impl { From c0caee0e2546abe3e1e91e4ced7e8f9d1326652b Mon Sep 17 00:00:00 2001 From: filipovic Date: Fri, 22 May 2026 08:43:52 +0200 Subject: [PATCH 9/9] Format HF cryo etching model --- examples/HFCryo4Cases/HFCryo4Cases.cpp | 70 ++++++----- examples/HFCryoEtching/HFCryoEtching.cpp | 83 ++++++------- include/viennaps/models/psHFCryoEtching.hpp | 117 ++++++++++-------- .../viennaps/models/psHFCryoParameters.hpp | 60 +++++---- 4 files changed, 175 insertions(+), 155 deletions(-) diff --git a/examples/HFCryo4Cases/HFCryo4Cases.cpp b/examples/HFCryo4Cases/HFCryo4Cases.cpp index f5a83670..d23d71f9 100644 --- a/examples/HFCryo4Cases/HFCryo4Cases.cpp +++ b/examples/HFCryo4Cases/HFCryo4Cases.cpp @@ -15,11 +15,11 @@ using namespace viennaps; -static constexpr double gridDelta = 1.0; -static constexpr double xExtent = 80.0; -static constexpr double yExtent = 2000.0; +static constexpr double gridDelta = 1.0; +static constexpr double xExtent = 80.0; +static constexpr double yExtent = 2000.0; static constexpr double trenchWidth = 10.0; -static constexpr double maskHeight = 80.0; +static constexpr double maskHeight = 80.0; static constexpr double processTime = 400.0; // Writes the surface mesh and returns the etch depth: how far the deepest @@ -30,8 +30,8 @@ static double saveMeshAndDepth(SmartPointer> domain, viennals::ToMultiSurfaceMesh(domain->getLevelSets(), mesh).apply(); const size_t n = mesh->nodes.size(); - mesh->getPointData().insertNextScalarData( - std::vector(n, T), "Temperature_K"); + mesh->getPointData().insertNextScalarData(std::vector(n, T), + "Temperature_K"); viennals::VTKWriter(mesh, filename).apply(); double minY = 0.0; @@ -41,31 +41,33 @@ static double saveMeshAndDepth(SmartPointer> domain, } static double runCase(double T, bool useDiffusion) { - auto domain = SmartPointer>::New(gridDelta, xExtent, yExtent); - MakeTrench(domain, trenchWidth, 0.0, 0.0, - maskHeight, 0.0, false, Material::SiO2).apply(); + auto domain = + SmartPointer>::New(gridDelta, xExtent, yExtent); + MakeTrench(domain, trenchWidth, 0.0, 0.0, maskHeight, 0.0, false, + Material::SiO2) + .apply(); HFCryoParameters params; - params.temperature = T; - params.ionFlux = 3.0; - params.etchantFlux = 5.0e3; - params.gamma_HF = 0.9; - params.Desorption.nu0 = 1.0e8; - params.Desorption.E_des = 0.20; - params.Reaction.A_r = 3.0e2; - params.Reaction.E_a = 0.10; - params.DirectReaction.A_r = 1.0e4; - params.DirectReaction.E_a = 0.25; - params.IonActivation.A_act = 1.0; - params.Diffusion.D0 = 1.0e3; - params.Diffusion.omega = 0.25; - params.Ions.meanEnergy = 150.0; - params.Ions.sigmaEnergy = 10.0; - params.Ions.exponent = 600.0; + params.temperature = T; + params.ionFlux = 3.0; + params.etchantFlux = 5.0e3; + params.gamma_HF = 0.9; + params.Desorption.nu0 = 1.0e8; + params.Desorption.E_des = 0.20; + params.Reaction.A_r = 3.0e2; + params.Reaction.E_a = 0.10; + params.DirectReaction.A_r = 1.0e4; + params.DirectReaction.E_a = 0.25; + params.IonActivation.A_act = 1.0; + params.Diffusion.D0 = 1.0e3; + params.Diffusion.omega = 0.25; + params.Ions.meanEnergy = 150.0; + params.Ions.sigmaEnergy = 10.0; + params.Ions.exponent = 600.0; params.Config.useTemperatureDependence = true; - params.Config.usePhysisorption = true; - params.Config.useSurfaceDiffusion = useDiffusion; // toggled per run - params.Config.T_ref = 300.0; + params.Config.usePhysisorption = true; + params.Config.useSurfaceDiffusion = useDiffusion; // toggled per run + params.Config.T_ref = 300.0; auto model = SmartPointer>::New(params); Process process(domain, model, processTime); @@ -75,13 +77,13 @@ static double runCase(double T, bool useDiffusion) { process.setParameters(rayParams); AdvectionParameters advParams; - advParams.spatialScheme = SpatialScheme::WENO_3RD_ORDER; + advParams.spatialScheme = SpatialScheme::WENO_3RD_ORDER; advParams.temporalScheme = TemporalScheme::RUNGE_KUTTA_2ND_ORDER; process.setParameters(advParams); process.apply(); - const std::string tag = useDiffusion ? "diffON" : "diffOFF"; + const std::string tag = useDiffusion ? "diffON" : "diffOFF"; const std::string fname = "HFCryo_" + tag + "_T" + std::to_string((int)T) + "K.vtp"; return saveMeshAndDepth(domain, fname, T); @@ -107,8 +109,8 @@ int main() { std::cout << " T=" << std::setw(3) << (int)Tlist[i] << "K ... " << std::flush; depthOff[i] = runCase(Tlist[i], false); - std::cout << "depth = " << std::fixed << std::setprecision(1) - << depthOff[i] << " nm\n"; + std::cout << "depth = " << std::fixed << std::setprecision(1) << depthOff[i] + << " nm\n"; } std::cout << "\n--- Surface diffusion ON ---\n"; @@ -116,8 +118,8 @@ int main() { std::cout << " T=" << std::setw(3) << (int)Tlist[i] << "K ... " << std::flush; depthOn[i] = runCase(Tlist[i], true); - std::cout << "depth = " << std::fixed << std::setprecision(1) - << depthOn[i] << " nm\n"; + std::cout << "depth = " << std::fixed << std::setprecision(1) << depthOn[i] + << " nm\n"; } std::cout << "\n=== Summary: trench etch depth (nm) ===\n"; diff --git a/examples/HFCryoEtching/HFCryoEtching.cpp b/examples/HFCryoEtching/HFCryoEtching.cpp index 31f2508d..f591344f 100644 --- a/examples/HFCryoEtching/HFCryoEtching.cpp +++ b/examples/HFCryoEtching/HFCryoEtching.cpp @@ -9,17 +9,17 @@ using namespace viennaps; -static void runCase(const std::string &label, - double T, double gridDelta, - double xExtent, double yExtent, - double trenchWidth, double maskHeight, - double processTime, +static void runCase(const std::string &label, double T, double gridDelta, + double xExtent, double yExtent, double trenchWidth, + double maskHeight, double processTime, HFCryoParameters params) { params.temperature = T; - auto domain = SmartPointer>::New(gridDelta, xExtent, yExtent); - MakeTrench(domain, trenchWidth, 0.0, 0.0, - maskHeight, 0.0, false, Material::SiO2).apply(); + auto domain = + SmartPointer>::New(gridDelta, xExtent, yExtent); + MakeTrench(domain, trenchWidth, 0.0, 0.0, maskHeight, 0.0, false, + Material::SiO2) + .apply(); auto model = SmartPointer>::New(params); Process process(domain, model, processTime); @@ -29,7 +29,7 @@ static void runCase(const std::string &label, process.setParameters(rayParams); AdvectionParameters advParams; - advParams.spatialScheme = SpatialScheme::WENO_3RD_ORDER; + advParams.spatialScheme = SpatialScheme::WENO_3RD_ORDER; advParams.temporalScheme = TemporalScheme::RUNGE_KUTTA_2ND_ORDER; process.setParameters(advParams); @@ -48,26 +48,26 @@ int main() { units::Length::setUnit(units::Length::NANOMETER); units::Time::setUnit(units::Time::SECOND); - const double gridDelta = 1.0; - const double xExtent = 80.0; - const double yExtent = 700.0; - const double trenchWidth = 10.0; - const double maskHeight = 50.0; - const double processTime = 90.0; + const double gridDelta = 1.0; + const double xExtent = 80.0; + const double yExtent = 700.0; + const double trenchWidth = 10.0; + const double maskHeight = 50.0; + const double processTime = 90.0; // Base parameters shared by both cases HFCryoParameters params; - params.ionFlux = 1.0; - params.etchantFlux = 1.0e3; - params.gamma_HF = 0.9; - params.Desorption.nu0 = 1.0e8; - params.Desorption.E_des = 0.20; - params.Reaction.A_r = 3.0e2; - params.Reaction.E_a = 0.10; - params.IonActivation.A_act = 1.0; - params.Ions.meanEnergy = 100.0; - params.Ions.sigmaEnergy = 10.0; - params.Ions.exponent = 300.0; + params.ionFlux = 1.0; + params.etchantFlux = 1.0e3; + params.gamma_HF = 0.9; + params.Desorption.nu0 = 1.0e8; + params.Desorption.E_des = 0.20; + params.Reaction.A_r = 3.0e2; + params.Reaction.E_a = 0.10; + params.IonActivation.A_act = 1.0; + params.Ions.meanEnergy = 100.0; + params.Ions.sigmaEnergy = 10.0; + params.Ions.exponent = 300.0; // Print D_s per temperature std::cout << "=== HF Cryo Etching | Surface Diffusion Comparison ===\n\n"; @@ -75,46 +75,45 @@ int main() { std::cout << "------------------------------------------------------\n"; double Tlist[4] = {150.0, 200.0, 250.0, 300.0}; for (int i = 0; i < 4; ++i) { - params.temperature = Tlist[i]; - params.Diffusion.D0 = 1.0e3; + params.temperature = Tlist[i]; + params.Diffusion.D0 = 1.0e3; params.Diffusion.omega = 0.25; - std::cout << Tlist[i] << "K " - << params.k_des() << " " - << params.k_r() << " " - << params.k_r_direct() << " " - << params.D_s() << "\n"; + std::cout << Tlist[i] << "K " << params.k_des() << " " << params.k_r() + << " " << params.k_r_direct() << " " << params.D_s() << "\n"; } std::cout << "\n"; // Save initial geometry { - auto domain = SmartPointer>::New(gridDelta, xExtent, yExtent); - MakeTrench(domain, trenchWidth, 0.0, 0.0, - maskHeight, 0.0, false, Material::SiO2).apply(); + auto domain = + SmartPointer>::New(gridDelta, xExtent, yExtent); + MakeTrench(domain, trenchWidth, 0.0, 0.0, maskHeight, 0.0, false, + Material::SiO2) + .apply(); domain->saveSurfaceMesh("HFCryo_initial.vtp", true); std::cout << "Initial geometry saved.\n\n"; } // ── Case A: No surface diffusion (D0=0) ────────────────────────────────── std::cout << "=== Case A: No Surface Diffusion (D0=0) ===\n"; - params.Diffusion.D0 = 0.0; + params.Diffusion.D0 = 0.0; params.Diffusion.omega = 0.25; for (int i = 0; i < 4; ++i) { std::cout << " T=" << Tlist[i] << "K ... "; std::cout.flush(); - runCase("NoDiff", Tlist[i], gridDelta, xExtent, yExtent, - trenchWidth, maskHeight, processTime, params); + runCase("NoDiff", Tlist[i], gridDelta, xExtent, yExtent, trenchWidth, + maskHeight, processTime, params); } // ── Case B: With surface diffusion (D0=1e3) ─────────────────────────────── std::cout << "\n=== Case B: With Surface Diffusion (D0=1e3) ===\n"; - params.Diffusion.D0 = 1.0e3; + params.Diffusion.D0 = 1.0e3; params.Diffusion.omega = 0.25; for (int i = 0; i < 4; ++i) { std::cout << " T=" << Tlist[i] << "K ... "; std::cout.flush(); - runCase("Diff", Tlist[i], gridDelta, xExtent, yExtent, - trenchWidth, maskHeight, processTime, params); + runCase("Diff", Tlist[i], gridDelta, xExtent, yExtent, trenchWidth, + maskHeight, processTime, params); } std::cout << "\n=== Done ===\n"; diff --git a/include/viennaps/models/psHFCryoEtching.hpp b/include/viennaps/models/psHFCryoEtching.hpp index f1072f3f..b7e66903 100644 --- a/include/viennaps/models/psHFCryoEtching.hpp +++ b/include/viennaps/models/psHFCryoEtching.hpp @@ -64,8 +64,9 @@ class HFCryoIon final NumericType f_ie = 1.; if (cosTheta < 0.5) - f_ie = std::max(NumericType(3) - NumericType(6) * angle / NumericType(M_PI), - NumericType(0)); + f_ie = + std::max(NumericType(3) - NumericType(6) * angle / NumericType(M_PI), + NumericType(0)); const NumericType sqrtE = std::sqrt(E); @@ -73,9 +74,9 @@ class HFCryoIon final params.SiO2.A_sp * std::max(sqrtE - std::sqrt(params.SiO2.Eth_sp), NumericType(0)); - const NumericType Y_ie = - params.SiO2.A_ie * - std::max(sqrtE - sqrt_E_th_ie, NumericType(0)) * f_ie; + const NumericType Y_ie = params.SiO2.A_ie * + std::max(sqrtE - sqrt_E_th_ie, NumericType(0)) * + f_ie; localData.getVectorData(0)[primID] += Y_sp * rayWeight; localData.getVectorData(1)[primID] += Y_ie * rayWeight; @@ -83,9 +84,8 @@ class HFCryoIon final std::pair> surfaceReflection(NumericType, const Vec3D &rayDir, - const Vec3D &geomNormal, - const unsigned int, const int, - const viennaray::TracingData *, + const Vec3D &geomNormal, const unsigned int, + const int, const viennaray::TracingData *, RNG &rng) override { const NumericType cosTheta = -DotProduct(rayDir, geomNormal); const NumericType incAngle = @@ -105,7 +105,8 @@ class HFCryoIon final } void initNew(RNG &rng) override { - E = initNormalDistEnergy(rng, params.Ions.meanEnergy, params.Ions.sigmaEnergy); + E = initNormalDistEnergy(rng, params.Ions.meanEnergy, + params.Ions.sigmaEnergy); } NumericType getSourceDistributionPower() const override { @@ -159,7 +160,8 @@ class HFCryoEtchant final const NumericType S_eff = params.gamma_HF * std::max(NumericType(1) - occupied, NumericType(0)); auto dir = viennaray::ReflectionDiffuse(geomNormal, rng); - return {NumericType(1) - S_eff, dir}; // reflected weight = 1 - sticking prob + return {NumericType(1) - S_eff, + dir}; // reflected weight = 1 - sticking prob } NumericType getSourceDistributionPower() const override { return 1.; } @@ -172,7 +174,8 @@ class HFCryoEtchant final const HFCryoParameters ¶ms; }; -// ── Surface model ───────────────────────────────────────────────────────────── +// ── Surface model +// ───────────────────────────────────────────────────────────── template class HFCryoSurfaceModel : public SurfaceModel { public: @@ -210,19 +213,19 @@ class HFCryoSurfaceModel : public SurfaceModel { // Called first; result is used as Dirichlet BC at chain endpoints. void updateCoverages(SmartPointer> rates, const std::vector &) override { - const auto numPoints = rates->getScalarData(0)->size(); + const auto numPoints = rates->getScalarData(0)->size(); const auto etchantFlux = rates->getScalarData("etchantFlux"); - const auto ionActFlux = rates->getScalarData("ionActivationFlux"); + const auto ionActFlux = rates->getScalarData("ionActivationFlux"); auto theta_phys = coverages->getScalarData("physCoverage"); auto theta_chem = coverages->getScalarData("chemCoverage"); theta_phys->resize(numPoints); theta_chem->resize(numPoints); - const NumericType k_des = params.effective_k_des(); - const NumericType k_r = params.effective_k_r(); + const NumericType k_des = params.effective_k_des(); + const NumericType k_r = params.effective_k_r(); const NumericType k_r_direct = params.effective_k_r_direct(); - const NumericType A_act = params.IonActivation.A_act; + const NumericType A_act = params.IonActivation.A_act; #pragma omp parallel for for (size_t i = 0; i < numPoints; ++i) { @@ -237,16 +240,17 @@ class HFCryoSurfaceModel : public SurfaceModel { theta_chem->at(i) = NumericType(0); } else { // Two-state: physisorbed + ion-activated chemisorbed - const NumericType k_act_ion = ionActFlux->at(i) * params.ionFlux * A_act; + const NumericType k_act_ion = + ionActFlux->at(i) * params.ionFlux * A_act; const NumericType denom = Gamma_HF + k_des + k_act_ion + k_r_direct; theta_phys->at(i) = (denom > NumericType(1e-30)) ? std::min(Gamma_HF / denom, NumericType(1)) : NumericType(0); const NumericType num_chem = k_act_ion * theta_phys->at(i); - theta_chem->at(i) = (k_r > NumericType(1e-30)) - ? std::min(num_chem / k_r, - NumericType(1) - theta_phys->at(i)) - : NumericType(0); + theta_chem->at(i) = + (k_r > NumericType(1e-30)) + ? std::min(num_chem / k_r, NumericType(1) - theta_phys->at(i)) + : NumericType(0); } } } @@ -255,24 +259,25 @@ class HFCryoSurfaceModel : public SurfaceModel { calculateVelocities(SmartPointer> rates, const std::vector> &coordinates, const std::vector &materialIds) override { - // Step 1: local coverage (no diffusion) — also serves as Dirichlet BC values + // Step 1: local coverage (no diffusion) — also serves as Dirichlet BC + // values updateCoverages(rates, materialIds); const auto numPoints = rates->getScalarData(0)->size(); const auto etchantFlux = rates->getScalarData("etchantFlux"); - const auto ionActFlux = rates->getScalarData("ionActivationFlux"); + const auto ionActFlux = rates->getScalarData("ionActivationFlux"); - const NumericType k_des = params.effective_k_des(); - const NumericType k_r = params.effective_k_r(); + const NumericType k_des = params.effective_k_des(); + const NumericType k_r = params.effective_k_r(); const NumericType k_r_direct = params.effective_k_r_direct(); const NumericType A_act = params.IonActivation.A_act; - const NumericType Ds = params.D_s(); + const NumericType Ds = params.D_s(); std::vector Gamma_HF_vec(numPoints); std::vector k_act_ion_vec(numPoints); for (size_t i = 0; i < numPoints; ++i) { - Gamma_HF_vec[i] = etchantFlux->at(i) * params.etchantFlux; + Gamma_HF_vec[i] = etchantFlux->at(i) * params.etchantFlux; k_act_ion_vec[i] = ionActFlux->at(i) * params.ionFlux * A_act; } @@ -280,7 +285,7 @@ class HFCryoSurfaceModel : public SurfaceModel { auto theta_chem = coverages->getScalarData("chemCoverage"); // Step 2: build surface chain + sidewall flags - const auto chain = buildSurfaceChain(coordinates); + const auto chain = buildSurfaceChain(coordinates); const auto sidewall = buildSidewallFlags(chain, coordinates); // Step 3: surface diffusion (only when enabled and D0 > 0) @@ -297,10 +302,10 @@ class HFCryoSurfaceModel : public SurfaceModel { continue; } const NumericType num_chem = k_act_ion_vec[i] * theta_phys->at(i); - theta_chem->at(i) = (k_r > NumericType(1e-30)) - ? std::min(num_chem / k_r, - NumericType(1) - theta_phys->at(i)) - : NumericType(0); + theta_chem->at(i) = + (k_r > NumericType(1e-30)) + ? std::min(num_chem / k_r, NumericType(1) - theta_phys->at(i)) + : NumericType(0); } } else { // Single coverage: no chemisorption @@ -308,7 +313,7 @@ class HFCryoSurfaceModel : public SurfaceModel { } // Step 5: etch rate - const auto ionSputterFlux = rates->getScalarData("ionSputterFlux"); + const auto ionSputterFlux = rates->getScalarData("ionSputterFlux"); const auto ionActivationFlux = rates->getScalarData("ionActivationFlux"); const NumericType unitConversion = units::Time::convertSecond() / units::Length::convertNanometer(); @@ -336,14 +341,13 @@ class HFCryoSurfaceModel : public SurfaceModel { // Ion-enhanced etching: Y_ie * Gamma_ion * theta_phys (or single theta) // Pure chemical: k_r_direct * theta_phys (only when usePhysisorption, // otherwise the concept of "physisorbed" HF is not defined) - const NumericType ionEnhanced = + const NumericType ionEnhanced = ionActivationFlux->at(i) * params.ionFlux * theta_phys->at(i); const NumericType pureChemical = k_r_direct * theta_phys->at(i); - const NumericType sputter = ionSputterFlux->at(i) * params.ionFlux; + const NumericType sputter = ionSputterFlux->at(i) * params.ionFlux; etchRate[i] = -(1. / params.SiO2.rho) * - (ionEnhanced + pureChemical + sputter) * - unitConversion; + (ionEnhanced + pureChemical + sputter) * unitConversion; if (Logger::hasIntermediate()) { chRate->at(i) = ionEnhanced + pureChemical; @@ -371,7 +375,8 @@ class HFCryoSurfaceModel : public SurfaceModel { std::vector buildSurfaceChain(const std::vector> &points) const { const size_t N = points.size(); - if (N == 0) return {}; + if (N == 0) + return {}; std::vector visited(N, false); std::vector chain; @@ -395,7 +400,8 @@ class HFCryoSurfaceModel : public SurfaceModel { size_t best = N; for (size_t j = 0; j < N; ++j) { - if (visited[j]) continue; + if (visited[j]) + continue; NumericType distSq = NumericType(0); for (int d = 0; d < D; ++d) { NumericType dd = points[j][d] - points[curr][d]; @@ -440,15 +446,17 @@ class HFCryoSurfaceModel : public SurfaceModel { tang[d] = points[chain[1]][d] - points[chain[0]][d]; else if (j == M - 1) for (int d = 0; d < D; ++d) - tang[d] = points[chain[M-1]][d] - points[chain[M-2]][d]; + tang[d] = points[chain[M - 1]][d] - points[chain[M - 2]][d]; else for (int d = 0; d < D; ++d) - tang[d] = points[chain[j+1]][d] - points[chain[j-1]][d]; + tang[d] = points[chain[j + 1]][d] - points[chain[j - 1]][d]; NumericType len = NumericType(0); - for (int d = 0; d < D; ++d) len += tang[d] * tang[d]; + for (int d = 0; d < D; ++d) + len += tang[d] * tang[d]; len = std::sqrt(len); - if (len < NumericType(1e-10)) continue; + if (len < NumericType(1e-10)) + continue; // Sidewall: tangent mostly vertical (height component > 0.7) const NumericType vertFrac = std::abs(tang[D - 1]) / len; @@ -473,12 +481,12 @@ class HFCryoSurfaceModel : public SurfaceModel { // no reconstructed 1D chain and no imposed boundary values are needed. The // earlier chain + Dirichlet formulation injected geometry-reconstruction // noise that destabilised the etch front at large Ds. - void applySurfaceDiffusion( - std::vector &theta_phys, - const std::vector &Gamma_HF_vec, - const std::vector &k_act_ion_vec, - const std::vector> &points, - NumericType Ds, NumericType k_des, NumericType k_r_direct) const { + void applySurfaceDiffusion(std::vector &theta_phys, + const std::vector &Gamma_HF_vec, + const std::vector &k_act_ion_vec, + const std::vector> &points, + NumericType Ds, NumericType k_des, + NumericType k_r_direct) const { const size_t N = points.size(); if (N < 3) return; @@ -555,7 +563,8 @@ class HFCryoSurfaceModel : public SurfaceModel { } // namespace impl -// ── Public process model ────────────────────────────────────────────────────── +// ── Public process model +// ────────────────────────────────────────────────────── template class HFCryoEtching final : public ProcessModelCPU { public: @@ -580,12 +589,12 @@ class HFCryoEtching final : public ProcessModelCPU { VIENNACORE_LOG_ERROR("Units have not been set."); } - auto ion = std::make_unique>(params); - auto etchant = std::make_unique>(params); + auto ion = std::make_unique>(params); + auto etchant = + std::make_unique>(params); auto surfModel = SmartPointer>::New(params); - auto velField = - SmartPointer>::New(); + auto velField = SmartPointer>::New(); this->setSurfaceModel(surfModel); this->setVelocityField(velField); diff --git a/include/viennaps/models/psHFCryoParameters.hpp b/include/viennaps/models/psHFCryoParameters.hpp index 0df96c45..51bee24c 100644 --- a/include/viennaps/models/psHFCryoParameters.hpp +++ b/include/viennaps/models/psHFCryoParameters.hpp @@ -20,46 +20,49 @@ template struct HFCryoParameters { // Frenkel-Arrhenius desorption: k_des = nu0 * exp(-E_des / kB*T) struct DesorptionType { - NumericType nu0 = 1.0e13; // attempt frequency (1/s) - NumericType E_des = 0.25; // desorption activation energy (eV) + NumericType nu0 = 1.0e13; // attempt frequency (1/s) + NumericType E_des = 0.25; // desorption activation energy (eV) } Desorption; // Arrhenius reaction rate: k_r = A_r * exp(-E_a / kB*T) // Chemisorbed HF reacting with SiO2 (ion-assisted pathway) struct ReactionType { - NumericType A_r = 3.0e2; // pre-exponential factor (1e15 cm⁻²s⁻¹) - NumericType E_a = 0.10; // reaction activation energy (eV) + NumericType A_r = 3.0e2; // pre-exponential factor (1e15 cm⁻²s⁻¹) + NumericType E_a = 0.10; // reaction activation energy (eV) } Reaction; // Direct (pure) chemical etching: physisorbed HF reacts with SiO2 thermally, // no ion activation needed. Higher barrier than chemisorption pathway. // Negligible at cryo (150K~200K), significant at room temperature (300K). struct DirectReactionType { - NumericType A_r = 1.0e4; // pre-exponential factor - NumericType E_a = 0.25; // activation energy (eV) > chemisorption E_a + NumericType A_r = 1.0e4; // pre-exponential factor + NumericType E_a = 0.25; // activation energy (eV) > chemisorption E_a } DirectReaction; // Ion activation: ions convert physisorbed HF -> chemisorbed HF // k_act_ion = A_act * ionEnhancedFlux * ionFlux [same units as Gamma_HF] struct IonActivationType { - NumericType A_act = 1.0; // dimensionless scaling factor + NumericType A_act = 1.0; // dimensionless scaling factor } IonActivation; // Model configuration: select which physical features are active. // Enables progressive comparison: none → +temp → +phys → +diffusion struct ModelConfigType { - bool useTemperatureDependence = true; // Frenkel-Arrhenius/Arrhenius for rates - bool usePhysisorption = true; // two-state model (theta_phys + theta_chem) - bool useSurfaceDiffusion = true; // 1D diffusion PDE along surface chain - NumericType T_ref = 300.; // fixed reference T when !useTemperatureDependence + bool useTemperatureDependence = + true; // Frenkel-Arrhenius/Arrhenius for rates + bool usePhysisorption = true; // two-state model (theta_phys + theta_chem) + bool useSurfaceDiffusion = true; // 1D diffusion PDE along surface chain + NumericType T_ref = + 300.; // fixed reference T when !useTemperatureDependence } Config; // Surface diffusion of physisorbed HF along the surface // D_s(T) = D0 * exp(-omega * E_des / kB*T) - // E_diff = omega * E_des (corrugation ratio from Lill 2023, Table IV: omega~0.2-0.3) + // E_diff = omega * E_des (corrugation ratio from Lill 2023, Table IV: + // omega~0.2-0.3) struct DiffusionType { - NumericType D0 = 1.0e3; // nm²/s, pre-exponential - NumericType omega = 0.25; // corrugation ratio + NumericType D0 = 1.0e3; // nm²/s, pre-exponential + NumericType omega = 0.25; // corrugation ratio } Diffusion; // SiO2 material properties @@ -67,15 +70,15 @@ template struct HFCryoParameters { NumericType rho = constants::SiO2::rho; // 1e22 atoms/cm³ NumericType Eth_sp = 18.; // sputtering threshold (eV) NumericType A_sp = constants::SiO2::A_sp; - NumericType Eth_ie = 12.; // ion-enhanced etching threshold (eV) + NumericType Eth_ie = 12.; // ion-enhanced etching threshold (eV) NumericType A_ie = 2.0; } SiO2; // Ion properties struct IonType { - NumericType meanEnergy = 100.; // eV - NumericType sigmaEnergy = 10.; // eV - NumericType exponent = 300.; // angular distribution exponent + NumericType meanEnergy = 100.; // eV + NumericType sigmaEnergy = 10.; // eV + NumericType exponent = 300.; // angular distribution exponent NumericType inflectAngle = constants::Ion::inflectAngle; NumericType n_l = constants::Ion::n_l; NumericType minAngle = constants::Ion::minAngle; @@ -83,13 +86,16 @@ template struct HFCryoParameters { // Raw Arrhenius/Frenkel-Arrhenius at actual temperature NumericType k_des() const { - return Desorption.nu0 * std::exp(-Desorption.E_des / (constants::kB * temperature)); + return Desorption.nu0 * + std::exp(-Desorption.E_des / (constants::kB * temperature)); } NumericType k_r() const { - return Reaction.A_r * std::exp(-Reaction.E_a / (constants::kB * temperature)); + return Reaction.A_r * + std::exp(-Reaction.E_a / (constants::kB * temperature)); } NumericType k_r_direct() const { - return DirectReaction.A_r * std::exp(-DirectReaction.E_a / (constants::kB * temperature)); + return DirectReaction.A_r * + std::exp(-DirectReaction.E_a / (constants::kB * temperature)); } NumericType D_s() const { const NumericType E_diff = Diffusion.omega * Desorption.E_des; @@ -98,16 +104,20 @@ template struct HFCryoParameters { // Config-aware rates: use T_ref when !useTemperatureDependence NumericType effective_k_des() const { - const NumericType T = Config.useTemperatureDependence ? temperature : Config.T_ref; + const NumericType T = + Config.useTemperatureDependence ? temperature : Config.T_ref; return Desorption.nu0 * std::exp(-Desorption.E_des / (constants::kB * T)); } NumericType effective_k_r() const { - const NumericType T = Config.useTemperatureDependence ? temperature : Config.T_ref; + const NumericType T = + Config.useTemperatureDependence ? temperature : Config.T_ref; return Reaction.A_r * std::exp(-Reaction.E_a / (constants::kB * T)); } NumericType effective_k_r_direct() const { - const NumericType T = Config.useTemperatureDependence ? temperature : Config.T_ref; - return DirectReaction.A_r * std::exp(-DirectReaction.E_a / (constants::kB * T)); + const NumericType T = + Config.useTemperatureDependence ? temperature : Config.T_ref; + return DirectReaction.A_r * + std::exp(-DirectReaction.E_a / (constants::kB * T)); } };