From 10e6c5d264369af71e1f960d43c959ed40df9198 Mon Sep 17 00:00:00 2001 From: AlexHocks Date: Wed, 6 Mar 2024 16:13:02 +0100 Subject: [PATCH 1/3] Langevin Thermo/Integrator implemented --- src/Simulation.cpp | 38 +++-- src/Simulation.h | 7 + src/integrators/Langevin.cpp | 177 ++++++++++++++++++++++++ src/integrators/Langevin.h | 55 ++++++++ src/thermostats/TemperatureObserver.cpp | 150 ++++++++++++++++++++ src/thermostats/TemperatureObserver.h | 142 +++++++++++++++++++ 6 files changed, 561 insertions(+), 8 deletions(-) create mode 100644 src/integrators/Langevin.cpp create mode 100644 src/integrators/Langevin.h create mode 100644 src/thermostats/TemperatureObserver.cpp create mode 100644 src/thermostats/TemperatureObserver.h diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 53be31d563..8bda1146f4 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -41,6 +41,7 @@ #include "integrators/Integrator.h" #include "integrators/Leapfrog.h" #include "integrators/LeapfrogRMM.h" +#include "integrators/Langevin.h" #include "plugins/PluginBase.h" #include "plugins/PluginFactory.h" @@ -115,6 +116,7 @@ Simulation::Simulation() _rand(8624), _longRangeCorrection(nullptr), _temperatureControl(nullptr), + _temperatureObserver(nullptr), _FMM(nullptr), _timerProfiler(), #ifdef TASKTIMINGPROFILE @@ -150,6 +152,8 @@ Simulation::~Simulation() { _longRangeCorrection = nullptr; delete _temperatureControl; _temperatureControl = nullptr; + delete _temperatureObserver; + _temperatureObserver = nullptr; delete _FMM; _FMM = nullptr; @@ -182,7 +186,9 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) { _integrator = new Leapfrog(); } else if (integratorType == "LeapfrogRMM") { _integrator = new LeapfrogRMM(); - } else { + } else if (integratorType == "Langevin") { + _integrator = new Langevin(); + } else { Log::global_log-> error() << "Unknown integrator " << integratorType << std::endl; Simulation::exit(1); } @@ -530,8 +536,16 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) { Simulation::exit(-1); } } - else - { + else if (thermostattype == "TemperatureObserver") { + if (_temperatureObserver == nullptr) { + _temperatureObserver = new TemperatureObserver(); + _temperatureObserver->readXML(xmlconfig); + } else { + Log::global_log->error() << "Instance of TemperatureObserver already exists!" << std::endl; + Simulation::exit(-1); + } + } + else { Log::global_log->warning() << "Unknown thermostat " << thermostattype << std::endl; continue; } @@ -810,6 +824,10 @@ void Simulation::prepare_start() { _cellProcessor = new LegacyCellProcessor( _cutoffRadius, _LJCutoffRadius, _particlePairsHandler); } + if(dynamic_cast(_integrator) != nullptr) { + _integrator->init(); + } + if (_FMM != nullptr) { double globalLength[3]; @@ -898,6 +916,11 @@ void Simulation::prepare_start() { if(nullptr != _temperatureControl) _temperatureControl->prepare_start(); // Has to be called before plugin initialization (see below): plugin->init(...) + if(_temperatureObserver != nullptr) { + _temperatureObserver->init(); + _temperatureObserver->step(_moleculeContainer); + } + // initializing plugins and starting plugin timers for (auto& plugin : _plugins) { Log::global_log->info() << "Initializing plugin " << plugin->getPluginName() << std::endl; @@ -1172,7 +1195,7 @@ void Simulation::simulateOneTimestep() // scale velocity and angular momentum // TODO: integrate into Temperature Control - if ( !_domain->NVE() && _temperatureControl == nullptr) { + if ( !_domain->NVE() && _temperatureControl == nullptr && _temperatureObserver == nullptr) { if (_thermostatType ==VELSCALE_THERMOSTAT) { Log::global_log->debug() << "Velocity scaling" << std::endl; if (_domain->severalThermostats()) { @@ -1204,10 +1227,9 @@ void Simulation::simulateOneTimestep() } else if ( _temperatureControl != nullptr) { // mheinen 2015-07-27 --> TEMPERATURE_CONTROL _temperatureControl->DoLoopsOverMolecules(_domainDecomposition, _moleculeContainer, _simstep); - } - // <-- TEMPERATURE_CONTROL - - + } else if (_temperatureObserver != nullptr) { + _temperatureObserver->step(_moleculeContainer); + } advanceSimulationTime(_integrator->getTimestepLength()); diff --git a/src/Simulation.h b/src/Simulation.h index 93055edfd2..f97de44a3b 100644 --- a/src/Simulation.h +++ b/src/Simulation.h @@ -15,6 +15,7 @@ // plugins #include "plugins/PluginFactory.h" +#include "thermostats/TemperatureObserver.h" #if !defined (SIMULATION_SRC) or defined (IN_IDE_PARSER) class Simulation; @@ -226,6 +227,9 @@ class Simulation { /** Get pointer to the molecule container */ ParticleContainer* getMoleculeContainer() { return _moleculeContainer; } + /** Get pointer to the temperature observer */ + TemperatureObserver* getTemperatureObserver() { return _temperatureObserver; } + /** Set the number of time steps to be performed in the simulation */ void setNumTimesteps( unsigned long steps ) { _numberOfTimesteps = steps; } /** Get the number of time steps to be performed in the simulation */ @@ -391,6 +395,9 @@ class Simulation { /** Temperature Control (Slab Thermostat) */ TemperatureControl* _temperatureControl; + /** No Thermostat - only measures temp in selected regions */ + TemperatureObserver* _temperatureObserver; + /** The Fast Multipole Method object */ bhfmm::FastMultipoleMethod* _FMM; diff --git a/src/integrators/Langevin.cpp b/src/integrators/Langevin.cpp new file mode 100644 index 0000000000..526b56d05b --- /dev/null +++ b/src/integrators/Langevin.cpp @@ -0,0 +1,177 @@ +// +// Created by alex on 05.03.24. +// + +#include "Langevin.h" +#include "utils/Logger.h" +#include "utils/mardyn_assert.h" +#include "utils/xmlfileUnits.h" +#include "particleContainer/ParticleContainer.h" +#include "Simulation.h" +#include "Domain.h" + +#include + +void Langevin::readXML(XMLfileUnits &xmlconfig) { + _timestepLength = 0; + xmlconfig.getNodeValueReduced("timestep", _timestepLength); + Log::global_log->info() << "Timestep: " << _timestepLength << std::endl; + mardyn_assert(_timestepLength > 0); + + _xi = 0; + xmlconfig.getNodeValue("friction", _xi); + mardyn_assert(_xi > 0); +} + +void Langevin::init() { + std::vector, std::array>> regions; + if(_simulation.getTemperatureObserver() == nullptr) return; // we will check again later + + _simulation.getTemperatureObserver()->getRegions(regions); + for(auto& [low, high] : regions) { + _stochastic_regions.emplace_back(low, high); + } + + _dt_half = _timestepLength / 2; +} + +void Langevin::eventForcesCalculated(ParticleContainer *particleContainer, Domain *domain) { + for(std::size_t index = 0; index < _stochastic_regions.size(); index++) { + auto& region = _stochastic_regions[index]; + double T = _simulation.getTemperatureObserver()->getTemperature(index); + T = _simulation.getEnsemble()->T(); + #if defined(_OPENMP) + #pragma omp parallel + #endif + for(auto it = particleContainer->regionIterator(std::data(region.low), + std::data(region.high), + ParticleIterator::ONLY_INNER_AND_BOUNDARY); it.isValid(); ++it) { + d3 v_old = it->v_arr(); + d3 v_rand = sampleRandomForce(it->mass(), T); + d3 v_delta {}; + + for(int d = 0; d < 3; d++) { + v_delta[d] = - _dt_half * _xi * v_old[d] + v_rand[d]; + } + + it->vadd(v_delta[0], v_delta[1], v_delta[2]); + } + } + + + std::map N; + std::map rotDOF; + std::map summv2; + std::map sumIw2; + + if (domain->severalThermostats()) { + #if defined(_OPENMP) + #pragma omp parallel + #endif + { + std::map N_l; + std::map rotDOF_l; + std::map summv2_l; + std::map sumIw2_l; + + for (auto tM = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); tM.isValid(); ++tM) { + int cid = tM->componentid(); + int thermostat = domain->getThermostat(cid); + tM->upd_postF(_dt_half, summv2_l[thermostat], sumIw2_l[thermostat]); + N_l[thermostat]++; + rotDOF_l[thermostat] += tM->component()->getRotationalDegreesOfFreedom(); + } + + #if defined(_OPENMP) + #pragma omp critical (thermostat) + #endif + { + for (auto & it : N_l) N[it.first] += it.second; + for (auto & it : rotDOF_l) rotDOF[it.first] += it.second; + for (auto & it : summv2_l) summv2[it.first] += it.second; + for (auto & it : sumIw2_l) sumIw2[it.first] += it.second; + } + } + } + else { + #if defined(_OPENMP) + #pragma omp parallel + #endif + { + unsigned long Ngt_l = 0; + unsigned long rotDOFgt_l = 0; + double summv2gt_l = 0.0; + double sumIw2gt_l = 0.0; + + for (auto i = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); i.isValid(); ++i) { + i->upd_postF(_dt_half, summv2gt_l, sumIw2gt_l); + mardyn_assert(summv2gt_l >= 0.0); + Ngt_l++; + rotDOFgt_l += i->component()->getRotationalDegreesOfFreedom(); + } + + #if defined(_OPENMP) + #pragma omp critical (thermostat) + #endif + { + N[0] += Ngt_l; + rotDOF[0] += rotDOFgt_l; + summv2[0] += summv2gt_l; + sumIw2[0] += sumIw2gt_l; + } + } // end pragma omp parallel + } + for (auto & thermit : summv2) { + domain->setLocalSummv2(thermit.second, thermit.first); + domain->setLocalSumIw2(sumIw2[thermit.first], thermit.first); + domain->setLocalNrotDOF(thermit.first, N[thermit.first], rotDOF[thermit.first]); + } +} + +void Langevin::eventNewTimestep(ParticleContainer *particleContainer, Domain *domain) { + for(std::size_t index = 0; index < _stochastic_regions.size(); index++) { + auto& region = _stochastic_regions[index]; + double T = _simulation.getTemperatureObserver()->getTemperature(index); + T = _simulation.getEnsemble()->T(); + #if defined(_OPENMP) + #pragma omp parallel + #endif + for(auto it = particleContainer->regionIterator(std::data(region.low), + std::data(region.high), + ParticleIterator::ONLY_INNER_AND_BOUNDARY); it.isValid(); ++it) { + d3 v_old = it->v_arr(); + d3 v_rand = sampleRandomForce(it->mass(), T); + d3 v_delta {}; + + for(int d = 0; d < 3; d++) { + v_delta[d] = - _dt_half * _xi * v_old[d] + v_rand[d]; + } + + it->vadd(v_delta[0], v_delta[1], v_delta[2]); + } + } + + #if defined(_OPENMP) + #pragma omp parallel + #endif + { + for (auto i = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); i.isValid(); ++i) { + i->upd_preF(_timestepLength); + } + } + + +} + +Langevin::d3 Langevin::sampleRandomForce(double m, double T) { + static thread_local std::mt19937 generator; // keeping default seed for reproducibility + std::normal_distribution normal{0.0, 1.0}; + d3 r_vec {}; + + double scale = std::sqrt(_timestepLength * T * _xi / m); + for(int d = 0; d < 3; d++) { + r_vec[d] = scale * normal(generator); + } + + return r_vec; +} diff --git a/src/integrators/Langevin.h b/src/integrators/Langevin.h new file mode 100644 index 0000000000..6f69fc351a --- /dev/null +++ b/src/integrators/Langevin.h @@ -0,0 +1,55 @@ +// +// Created by Alex Hocks on 05.03.24. +// + +#ifndef MARDYN_LANGEVIN_H +#define MARDYN_LANGEVIN_H + + +#include "Integrator.h" + +#include +#include + +/** + * New time integrator that deviates from the standard Verlet integration scheme. + * Equations of motion are changed to follow Langevin equations. + * This is a stochastic integration scheme, which is incorporated as random "kicks" from the connected heat bath. + * */ +class Langevin : public Integrator { +public: + void readXML(XMLfileUnits &xmlconfig) override; + + void init() override; + + void eventNewTimestep(ParticleContainer *moleculeContainer, Domain *domain) override; + + void eventForcesCalculated(ParticleContainer *moleculeContainer, Domain *domain) override; +private: + using d3 = std::array; + struct box_t { + box_t() = default; + box_t(const d3& l, const d3& h) : low(l), high(h) {} + d3 low, high; + }; + + //! @brief friction strength + double _xi; + + //! @brief time step length halved + double _dt_half; + + //! @brief regions in which friction is not set to 0 + std::vector _stochastic_regions; + + /** + * Sample from Gaussian with technically 0 mean and sigma**2 = 2 * m * friction * k_b * T_target / delta_t. + * Is already adapted to match the integration scheme. + * @param m mass + * @param T temp target + * */ + d3 sampleRandomForce(double m, double T); +}; + + +#endif //MARDYN_LANGEVIN_H diff --git a/src/thermostats/TemperatureObserver.cpp b/src/thermostats/TemperatureObserver.cpp new file mode 100644 index 0000000000..b3477bc104 --- /dev/null +++ b/src/thermostats/TemperatureObserver.cpp @@ -0,0 +1,150 @@ +// +// Created by alex on 06.03.24. +// + +#include "TemperatureObserver.h" +#include "Simulation.h" +#include "ensemble/EnsembleBase.h" + +void TemperatureObserver::readXML(XMLfileUnits &xmlconfig) { + std::size_t num_regions = 0; + XMLfile::Query query = xmlconfig.query("ActiveRegions/region"); + num_regions = query.card(); + _region_data.resize(num_regions); + + std::string oldpath = xmlconfig.getcurrentnodepath(); + XMLfile::Query::const_iterator rIt; + std::size_t index = 0; + for(rIt = query.begin(); rIt; rIt++) { + xmlconfig.changecurrentnode(rIt); + d3 low {}, high{}; + xmlconfig.getNodeValue("lowX", low[0]); + xmlconfig.getNodeValue("lowY", low[1]); + xmlconfig.getNodeValue("lowZ", low[2]); + + xmlconfig.getNodeValue("highX", high[0]); + xmlconfig.getNodeValue("highY", high[1]); + xmlconfig.getNodeValue("highZ", high[2]); + + _region_data[index].first.low = low; + _region_data[index].first.high = high; + index++; + } + xmlconfig.changecurrentnode(oldpath); + + _max_samples = 100; + xmlconfig.getNodeValue("samples", _max_samples); +} + +void TemperatureObserver::init() { + for(auto& [region, region_cache] : _region_data) { + region_cache.history.init(_max_samples); + region_cache.current_temp = 0.0; + } +} + +void TemperatureObserver::step(ParticleContainer* particleContainer) { + for(auto& [region, region_cache] : _region_data) { + region_cache.current_temp = measureTemp(region.low, region.high, particleContainer, region_cache.history); + } +} + +void TemperatureObserver::getRegions(std::vector, std::array>>& regions) { + for(auto& [region, region_cache] : _region_data) { + regions.emplace_back(region.low, region.high); + } +} + +double TemperatureObserver::getTemperature(std::size_t index) { + mardyn_assert(index < _region_data.size() && index >= 0); + return _region_data[index].second.current_temp; +} + +double TemperatureObserver::measureTemp(const std::array &low, const std::array &high, + ParticleContainer *particleContainer, BinData& binData) { + std::array v_mean = computeMeanVelocityStep(low, high, particleContainer, binData); + double m = _simulation.getEnsemble()->getComponent(0)->m(); + double v = 0; + std::size_t n = 0; +#if defined(_OPENMP) +#pragma omp parallel reduction(+:v, n) +#endif + for (auto it = particleContainer->regionIterator(std::data(low), std::data(high), + ParticleIterator::ALL_CELLS); it.isValid(); ++it) { + for(int d = 0; d < 3; d++) { + v += std::pow(it->v(d) - v_mean[d], 2); + } + n += 1; + } + + return v * m / (3 * static_cast(n)); +} + +std::array TemperatureObserver::computeMeanVelocityStep(const std::array &low, const std::array &high, + ParticleContainer *particleContainer, BinData& binData) { + //compute mean velocity per dim using SMC method, see: Karimian et al. 2011 + std::array v_avg{0}; + double* v_raw = std::data(v_avg); + std::size_t count = 0; +#if defined(_OPENMP) +#pragma omp parallel reduction(+:v_raw[:3], count) +#endif + for (auto it = particleContainer->regionIterator(std::data(low), std::data(high), + ParticleIterator::ALL_CELLS); it.isValid(); ++it) { + for(int d = 0; d < 3; d++) { + v_raw[d] += it->v(d); + } + count += 1; + } + //TODO check this + if(count == 0) count = 1; + binData.mol_counts.insert(count); + binData.velocities.insert(v_avg); + + //compute CAM average first + { + std::size_t steps = binData.mol_counts.size(); + auto d = static_cast(steps<=1?1:steps-1); + + std::array v_cam{ 0 }; + for(std::size_t i = 0; i < steps; i++) { + const auto& vec = binData.velocities.get(i); + for(int dim = 0; dim < 3; dim++) { + v_cam[dim] += vec[dim]; + } + } + for(int dim = 0; dim < 3; dim++) { + v_cam[dim] /= d; + } + + std::size_t nks = 0; + for(std::size_t i = 0; i < steps; i++) { + nks += binData.mol_counts.get(i); + } + double nks_d = static_cast(nks) / d; + + for(int dim = 0; dim < 3; dim++) { + v_cam[dim] /= nks_d; + } + binData.v_cam.insert(v_cam); + } + + //compute SAM over CAM + { + std::size_t steps = binData.v_cam.size(); + auto d = static_cast(steps<=1?1:steps-1); + + std::array v_scm{ 0 }; + for(std::size_t i = 0; i < steps; i++) { + const auto& vec = binData.v_cam.get(i); + for(int dim = 0; dim < 3; dim++) { + v_scm[dim] += vec[dim]; + } + } + for(int dim = 0; dim < 3; dim++) { + v_scm[dim] /= d; + } + + return v_scm; + } +} diff --git a/src/thermostats/TemperatureObserver.h b/src/thermostats/TemperatureObserver.h new file mode 100644 index 0000000000..269aa5648c --- /dev/null +++ b/src/thermostats/TemperatureObserver.h @@ -0,0 +1,142 @@ +// +// Created by alex on 06.03.24. +// + +#ifndef MARDYN_TEMPERATUREOBSERVER_H +#define MARDYN_TEMPERATUREOBSERVER_H + + +#include "utils/xmlfileUnits.h" +#include "utils/mardyn_assert.h" + +#include +#include + + +class Domain; +class ParticleContainer; +class Ensemble; + +/** + * Ring buffer that stores up to n entries. Will override oldest entry once n+1 elements are inserted. + * */ +template +struct NBuffer { + NBuffer() = default; + [[maybe_unused]] explicit NBuffer(std::size_t n) : _n(n), _begin(-1), _end(0) { + _data.resize(n); + } + void init(std::size_t n) { + _n = n; + _begin = -1; + _end = 0; + _data.resize(n); + } + void insert(const T& t) { + if(_begin == -1) { + _data[0] = t; + _begin = 0; + _end = 1; + } + else if (_begin == _end) { + _data[_end] = t; + _begin = (_begin+1) % _n; + _end = _begin; + } + else { + _data[_end] = t; + _end = (_end+1) % _n; + } + } + const T& get(std::size_t index) { + mardyn_assert(index < _n); + return _data[(_begin + index) % _n]; + } + std::size_t size() { + if(_begin == -1) { + return 0; + } + else if (_begin == _end) { + return _n; + } + else { + return _end - _begin; + } + } + std::size_t _n{}; + long _begin{}; + long _end{}; + std::vector _data; +}; + +/** + * Data struct for computing SAM and CAM averages + * */ +struct BinData { + NBuffer mol_counts; + NBuffer> velocities; + NBuffer> v_cam; + + void init(std::size_t n) { + mol_counts.init(n); + velocities.init(n); + v_cam.init(n); + } +}; + +/** + * Measures the temperature in selected regions. See Karimian et al. 2011 and 2014 + * */ +class TemperatureObserver { +public: + /** + * Init all buffers + * */ + void init(); + void readXML(XMLfileUnits &xmlconfig); + /** + * Measure temp in all regions, updates caches and stores current temp + * */ + void step(ParticleContainer* particleContainer); + + /** + * @returns all regions, for which the temp is being measured. + * */ + void getRegions(std::vector, std::array>>& regions); + + /** + * Gets the temperature for the selected region. + * Regions can be queried by accessing getRegions(). + * @returns T * k_b + * */ + double getTemperature(std::size_t index); +private: + std::size_t _max_samples; + + using d3 = std::array; + struct box_t { + box_t() = default; + box_t(const d3& l, const d3& h) : low(l), high(h) {} + d3 low, high; + }; + struct box_data_t { + BinData history{ }; + double current_temp{ }; + }; + + //! @brief each region is defined by a box in 3d space (box_t) and a data cache (box_data_t) + std::vector> _region_data; + + /** + * Based on: Details about pressure calculation in molecular dynamic analysis by Karimian et al. 2014 + * returns T * k_b + * */ + double measureTemp(const std::array& low, const std::array& high, ParticleContainer *particleContainer, BinData& binData); + + /** + * Computes the SAM over the CAM average of the mean velocity for the specified region. (SCM method) + * */ + std::array computeMeanVelocityStep(const std::array& low, const std::array& high, ParticleContainer *particleContainer, BinData& binData); +}; + +#endif //MARDYN_TEMPERATUREOBSERVER_H From c5e0029c868a98e3dbcc744268a09e29eca23b55 Mon Sep 17 00:00:00 2001 From: AlexHocks Date: Fri, 8 Mar 2024 14:37:50 +0100 Subject: [PATCH 2/3] Target Temp is now reached --- src/integrators/Langevin.cpp | 21 +++++++++++++++++---- src/integrators/Langevin.h | 6 ++++++ 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/src/integrators/Langevin.cpp b/src/integrators/Langevin.cpp index 526b56d05b..559a333fb8 100644 --- a/src/integrators/Langevin.cpp +++ b/src/integrators/Langevin.cpp @@ -13,6 +13,7 @@ #include void Langevin::readXML(XMLfileUnits &xmlconfig) { + _checkFailed = false; _timestepLength = 0; xmlconfig.getNodeValueReduced("timestep", _timestepLength); Log::global_log->info() << "Timestep: " << _timestepLength << std::endl; @@ -25,14 +26,23 @@ void Langevin::readXML(XMLfileUnits &xmlconfig) { void Langevin::init() { std::vector, std::array>> regions; - if(_simulation.getTemperatureObserver() == nullptr) return; // we will check again later + _dt_half = _timestepLength / 2; + + if(_simulation.getTemperatureObserver() == nullptr && _checkFailed) { + Log::global_log->warning() << "Langevin Integrator used, but no Temperature Observer defined as thermostat." + << " Will behave identical to leapfrog integrator." << std::endl; + } + if(_simulation.getTemperatureObserver() == nullptr && !_checkFailed) { _checkFailed = true; return; } // we will check again later _simulation.getTemperatureObserver()->getRegions(regions); + if(regions.empty()) { + Log::global_log->warning() << "Langevin Integrator used, but no regions defined in Temperature Observer." + << " Will behave identical to leapfrog integrator." << std::endl; + return; + } for(auto& [low, high] : regions) { _stochastic_regions.emplace_back(low, high); } - - _dt_half = _timestepLength / 2; } void Langevin::eventForcesCalculated(ParticleContainer *particleContainer, Domain *domain) { @@ -168,7 +178,10 @@ Langevin::d3 Langevin::sampleRandomForce(double m, double T) { std::normal_distribution normal{0.0, 1.0}; d3 r_vec {}; - double scale = std::sqrt(_timestepLength * T * _xi / m); + // additional factor 2 here + // according to book about Langevin integration not needed, but according to Langevin's equations of motion it does exist + // by using it we actually reach the target temperature + double scale = std::sqrt(2 * _timestepLength * T * _xi / m); for(int d = 0; d < 3; d++) { r_vec[d] = scale * normal(generator); } diff --git a/src/integrators/Langevin.h b/src/integrators/Langevin.h index 6f69fc351a..cef1b07efe 100644 --- a/src/integrators/Langevin.h +++ b/src/integrators/Langevin.h @@ -15,6 +15,9 @@ * New time integrator that deviates from the standard Verlet integration scheme. * Equations of motion are changed to follow Langevin equations. * This is a stochastic integration scheme, which is incorporated as random "kicks" from the connected heat bath. + * + * This is to be used with the TemperatureObserver, which disables the normal velocity scaling and defines + * the regions in which the Langevin Thermostat should be active. * */ class Langevin : public Integrator { public: @@ -42,6 +45,9 @@ class Langevin : public Integrator { //! @brief regions in which friction is not set to 0 std::vector _stochastic_regions; + //! @brief We need to check at least twice if a TemperatureObserver exists + bool _checkFailed; + /** * Sample from Gaussian with technically 0 mean and sigma**2 = 2 * m * friction * k_b * T_target / delta_t. * Is already adapted to match the integration scheme. From 0e2a9d686b5a73cedb280e7f0e5d449193ecfc69 Mon Sep 17 00:00:00 2001 From: AlexHocks Date: Mon, 11 Mar 2024 12:15:12 +0100 Subject: [PATCH 3/3] Fixed review items --- src/Simulation.cpp | 71 +++--- src/Simulation.h | 11 +- src/integrators/Langevin.cpp | 320 ++++++++++++------------ src/integrators/Langevin.h | 82 +++--- src/thermostats/TemperatureObserver.cpp | 251 ++++++++++--------- src/thermostats/TemperatureObserver.h | 227 +++++++++-------- 6 files changed, 507 insertions(+), 455 deletions(-) diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 8bda1146f4..3afdf81418 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -69,6 +69,7 @@ #include "ensemble/CavityEnsemble.h" #include "thermostats/VelocityScalingThermostat.h" +#include "thermostats/TemperatureObserver.h" #include "thermostats/TemperatureControl.h" #include "utils/FileUtils.h" @@ -116,7 +117,7 @@ Simulation::Simulation() _rand(8624), _longRangeCorrection(nullptr), _temperatureControl(nullptr), - _temperatureObserver(nullptr), + _temperatureObserver(nullptr), _FMM(nullptr), _timerProfiler(), #ifdef TASKTIMINGPROFILE @@ -152,8 +153,8 @@ Simulation::~Simulation() { _longRangeCorrection = nullptr; delete _temperatureControl; _temperatureControl = nullptr; - delete _temperatureObserver; - _temperatureObserver = nullptr; + delete _temperatureObserver; + _temperatureObserver = nullptr; delete _FMM; _FMM = nullptr; @@ -187,8 +188,9 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) { } else if (integratorType == "LeapfrogRMM") { _integrator = new LeapfrogRMM(); } else if (integratorType == "Langevin") { - _integrator = new Langevin(); - } else { + _integrator = new Langevin(); + _thermostatType = LANGEVIN_THERMOSTAT; + } else { Log::global_log-> error() << "Unknown integrator " << integratorType << std::endl; Simulation::exit(1); } @@ -527,25 +529,25 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) { } } else if(thermostattype == "TemperatureControl") { - if (nullptr == _temperatureControl) { - _temperatureControl = new TemperatureControl(); - _temperatureControl->readXML(xmlconfig); - } else { - Log::global_log->error() << "Instance of TemperatureControl allready exist! Programm exit ..." - << std::endl; - Simulation::exit(-1); - } - } + if (nullptr == _temperatureControl) { + _temperatureControl = new TemperatureControl(); + _temperatureControl->readXML(xmlconfig); + } else { + Log::global_log->error() << "Instance of TemperatureControl allready exist! Programm exit ..." + << std::endl; + Simulation::exit(-1); + } + } else if (thermostattype == "TemperatureObserver") { - if (_temperatureObserver == nullptr) { - _temperatureObserver = new TemperatureObserver(); - _temperatureObserver->readXML(xmlconfig); - } else { - Log::global_log->error() << "Instance of TemperatureObserver already exists!" << std::endl; - Simulation::exit(-1); - } - } - else { + if (_temperatureObserver == nullptr) { + _temperatureObserver = new TemperatureObserver(); + _temperatureObserver->readXML(xmlconfig); + } else { + Log::global_log->error() << "Instance of TemperatureObserver already exists!" << std::endl; + Simulation::exit(-1); + } + } + else { Log::global_log->warning() << "Unknown thermostat " << thermostattype << std::endl; continue; } @@ -824,9 +826,9 @@ void Simulation::prepare_start() { _cellProcessor = new LegacyCellProcessor( _cutoffRadius, _LJCutoffRadius, _particlePairsHandler); } - if(dynamic_cast(_integrator) != nullptr) { - _integrator->init(); - } + if(dynamic_cast(_integrator) != nullptr) { + _integrator->init(); + } if (_FMM != nullptr) { @@ -916,10 +918,10 @@ void Simulation::prepare_start() { if(nullptr != _temperatureControl) _temperatureControl->prepare_start(); // Has to be called before plugin initialization (see below): plugin->init(...) - if(_temperatureObserver != nullptr) { - _temperatureObserver->init(); - _temperatureObserver->step(_moleculeContainer); - } + if(_temperatureObserver != nullptr) { + _temperatureObserver->init(); + _temperatureObserver->step(_moleculeContainer); + } // initializing plugins and starting plugin timers for (auto& plugin : _plugins) { @@ -1195,7 +1197,7 @@ void Simulation::simulateOneTimestep() // scale velocity and angular momentum // TODO: integrate into Temperature Control - if ( !_domain->NVE() && _temperatureControl == nullptr && _temperatureObserver == nullptr) { + if ( !_domain->NVE() && _temperatureControl == nullptr) { if (_thermostatType ==VELSCALE_THERMOSTAT) { Log::global_log->debug() << "Velocity scaling" << std::endl; if (_domain->severalThermostats()) { @@ -1224,12 +1226,13 @@ void Simulation::simulateOneTimestep() } + else if (_thermostatType == LANGEVIN_THERMOSTAT) { + _temperatureObserver->step(_moleculeContainer); + } } else if ( _temperatureControl != nullptr) { // mheinen 2015-07-27 --> TEMPERATURE_CONTROL _temperatureControl->DoLoopsOverMolecules(_domainDecomposition, _moleculeContainer, _simstep); - } else if (_temperatureObserver != nullptr) { - _temperatureObserver->step(_moleculeContainer); - } + } advanceSimulationTime(_integrator->getTimestepLength()); diff --git a/src/Simulation.h b/src/Simulation.h index f97de44a3b..0ea9d1990d 100644 --- a/src/Simulation.h +++ b/src/Simulation.h @@ -15,7 +15,6 @@ // plugins #include "plugins/PluginFactory.h" -#include "thermostats/TemperatureObserver.h" #if !defined (SIMULATION_SRC) or defined (IN_IDE_PARSER) class Simulation; @@ -55,10 +54,12 @@ class LongRangeCorrection; class Homogeneous; class Planar; class TemperatureControl; +class TemperatureObserver; class MemoryProfiler; // by Stefan Becker const int VELSCALE_THERMOSTAT = 1; +const int LANGEVIN_THERMOSTAT = 2; namespace bhfmm { class FastMultipoleMethod; @@ -227,8 +228,8 @@ class Simulation { /** Get pointer to the molecule container */ ParticleContainer* getMoleculeContainer() { return _moleculeContainer; } - /** Get pointer to the temperature observer */ - TemperatureObserver* getTemperatureObserver() { return _temperatureObserver; } + /** Get pointer to the temperature observer */ + TemperatureObserver* getTemperatureObserver() { return _temperatureObserver; } /** Set the number of time steps to be performed in the simulation */ void setNumTimesteps( unsigned long steps ) { _numberOfTimesteps = steps; } @@ -395,8 +396,8 @@ class Simulation { /** Temperature Control (Slab Thermostat) */ TemperatureControl* _temperatureControl; - /** No Thermostat - only measures temp in selected regions */ - TemperatureObserver* _temperatureObserver; + /** No Thermostat - only measures temp in selected regions */ + TemperatureObserver* _temperatureObserver; /** The Fast Multipole Method object */ bhfmm::FastMultipoleMethod* _FMM; diff --git a/src/integrators/Langevin.cpp b/src/integrators/Langevin.cpp index 559a333fb8..9e43acf345 100644 --- a/src/integrators/Langevin.cpp +++ b/src/integrators/Langevin.cpp @@ -13,178 +13,178 @@ #include void Langevin::readXML(XMLfileUnits &xmlconfig) { - _checkFailed = false; - _timestepLength = 0; - xmlconfig.getNodeValueReduced("timestep", _timestepLength); - Log::global_log->info() << "Timestep: " << _timestepLength << std::endl; - mardyn_assert(_timestepLength > 0); - - _xi = 0; - xmlconfig.getNodeValue("friction", _xi); - mardyn_assert(_xi > 0); + _checkFailed = false; + _timestepLength = 0; + xmlconfig.getNodeValueReduced("timestep", _timestepLength); + Log::global_log->info() << "Timestep: " << _timestepLength << std::endl; + mardyn_assert(_timestepLength > 0); + + _gamma = 0; + xmlconfig.getNodeValue("friction", _gamma); + mardyn_assert(_gamma > 0); + + std::size_t num_regions = 0; + XMLfile::Query query = xmlconfig.query("ActiveRegions/region"); + num_regions = query.card(); + _stochastic_regions.resize(num_regions); + + std::string oldpath = xmlconfig.getcurrentnodepath(); + XMLfile::Query::const_iterator rIt; + std::size_t index = 0; + for(rIt = query.begin(); rIt; rIt++) { + xmlconfig.changecurrentnode(rIt); + d3 low {}, high{}; + xmlconfig.getNodeValue("lower/x", low[0]); + xmlconfig.getNodeValue("lower/y", low[1]); + xmlconfig.getNodeValue("lower/z", low[2]); + + xmlconfig.getNodeValue("upper/x", high[0]); + xmlconfig.getNodeValue("upper/y", high[1]); + xmlconfig.getNodeValue("upper/z", high[2]); + + _stochastic_regions[index].low = low; + _stochastic_regions[index].high = high; + index++; + } + xmlconfig.changecurrentnode(oldpath); } void Langevin::init() { - std::vector, std::array>> regions; - _dt_half = _timestepLength / 2; - - if(_simulation.getTemperatureObserver() == nullptr && _checkFailed) { - Log::global_log->warning() << "Langevin Integrator used, but no Temperature Observer defined as thermostat." - << " Will behave identical to leapfrog integrator." << std::endl; - } - if(_simulation.getTemperatureObserver() == nullptr && !_checkFailed) { _checkFailed = true; return; } // we will check again later - - _simulation.getTemperatureObserver()->getRegions(regions); - if(regions.empty()) { - Log::global_log->warning() << "Langevin Integrator used, but no regions defined in Temperature Observer." - << " Will behave identical to leapfrog integrator." << std::endl; - return; - } - for(auto& [low, high] : regions) { - _stochastic_regions.emplace_back(low, high); - } + _dt_half = _timestepLength / 2; + + if (_simulation.getTemperatureObserver() == nullptr && !_checkFailed) { + _checkFailed = true; + return; + } // we will check again later + + if (_simulation.getTemperatureObserver() == nullptr && _checkFailed) { + Log::global_log->warning() << "Langevin Integrator used, but no Temperature Observer defined as thermostat." + << " Behaviour is undefined." + << " Use Temperature Observer to disable other thermostat functionality!" << std::endl; + } } void Langevin::eventForcesCalculated(ParticleContainer *particleContainer, Domain *domain) { - for(std::size_t index = 0; index < _stochastic_regions.size(); index++) { - auto& region = _stochastic_regions[index]; - double T = _simulation.getTemperatureObserver()->getTemperature(index); - T = _simulation.getEnsemble()->T(); - #if defined(_OPENMP) - #pragma omp parallel - #endif - for(auto it = particleContainer->regionIterator(std::data(region.low), - std::data(region.high), - ParticleIterator::ONLY_INNER_AND_BOUNDARY); it.isValid(); ++it) { - d3 v_old = it->v_arr(); - d3 v_rand = sampleRandomForce(it->mass(), T); - d3 v_delta {}; - - for(int d = 0; d < 3; d++) { - v_delta[d] = - _dt_half * _xi * v_old[d] + v_rand[d]; - } - - it->vadd(v_delta[0], v_delta[1], v_delta[2]); - } - } - - - std::map N; - std::map rotDOF; - std::map summv2; - std::map sumIw2; - - if (domain->severalThermostats()) { - #if defined(_OPENMP) - #pragma omp parallel - #endif - { - std::map N_l; - std::map rotDOF_l; - std::map summv2_l; - std::map sumIw2_l; - - for (auto tM = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); tM.isValid(); ++tM) { - int cid = tM->componentid(); - int thermostat = domain->getThermostat(cid); - tM->upd_postF(_dt_half, summv2_l[thermostat], sumIw2_l[thermostat]); - N_l[thermostat]++; - rotDOF_l[thermostat] += tM->component()->getRotationalDegreesOfFreedom(); - } - - #if defined(_OPENMP) - #pragma omp critical (thermostat) - #endif - { - for (auto & it : N_l) N[it.first] += it.second; - for (auto & it : rotDOF_l) rotDOF[it.first] += it.second; - for (auto & it : summv2_l) summv2[it.first] += it.second; - for (auto & it : sumIw2_l) sumIw2[it.first] += it.second; - } - } - } - else { - #if defined(_OPENMP) - #pragma omp parallel - #endif - { - unsigned long Ngt_l = 0; - unsigned long rotDOFgt_l = 0; - double summv2gt_l = 0.0; - double sumIw2gt_l = 0.0; - - for (auto i = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); i.isValid(); ++i) { - i->upd_postF(_dt_half, summv2gt_l, sumIw2gt_l); - mardyn_assert(summv2gt_l >= 0.0); - Ngt_l++; - rotDOFgt_l += i->component()->getRotationalDegreesOfFreedom(); - } - - #if defined(_OPENMP) - #pragma omp critical (thermostat) - #endif - { - N[0] += Ngt_l; - rotDOF[0] += rotDOFgt_l; - summv2[0] += summv2gt_l; - sumIw2[0] += sumIw2gt_l; - } - } // end pragma omp parallel - } - for (auto & thermit : summv2) { - domain->setLocalSummv2(thermit.second, thermit.first); - domain->setLocalSumIw2(sumIw2[thermit.first], thermit.first); - domain->setLocalNrotDOF(thermit.first, N[thermit.first], rotDOF[thermit.first]); - } + addLangevinContribution(particleContainer); + + std::map N; + std::map rotDOF; + std::map summv2; + std::map sumIw2; + + if (domain->severalThermostats()) { + #if defined(_OPENMP) + #pragma omp parallel + #endif + { + std::map N_l; + std::map rotDOF_l; + std::map summv2_l; + std::map sumIw2_l; + + for (auto tM = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); tM.isValid(); ++tM) { + int cid = tM->componentid(); + int thermostat = domain->getThermostat(cid); + tM->upd_postF(_dt_half, summv2_l[thermostat], sumIw2_l[thermostat]); + N_l[thermostat]++; + rotDOF_l[thermostat] += tM->component()->getRotationalDegreesOfFreedom(); + } + + #if defined(_OPENMP) + #pragma omp critical (thermostat) + #endif + { + for (auto &it: N_l) N[it.first] += it.second; + for (auto &it: rotDOF_l) rotDOF[it.first] += it.second; + for (auto &it: summv2_l) summv2[it.first] += it.second; + for (auto &it: sumIw2_l) sumIw2[it.first] += it.second; + } + } + } else { + #if defined(_OPENMP) + #pragma omp parallel + #endif + { + unsigned long Ngt_l = 0; + unsigned long rotDOFgt_l = 0; + double summv2gt_l = 0.0; + double sumIw2gt_l = 0.0; + + for (auto i = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); i.isValid(); ++i) { + i->upd_postF(_dt_half, summv2gt_l, sumIw2gt_l); + mardyn_assert(summv2gt_l >= 0.0); + Ngt_l++; + rotDOFgt_l += i->component()->getRotationalDegreesOfFreedom(); + } + + #if defined(_OPENMP) + #pragma omp critical (thermostat) + #endif + { + N[0] += Ngt_l; + rotDOF[0] += rotDOFgt_l; + summv2[0] += summv2gt_l; + sumIw2[0] += sumIw2gt_l; + } + } // end pragma omp parallel + } + for (auto &thermit: summv2) { + domain->setLocalSummv2(thermit.second, thermit.first); + domain->setLocalSumIw2(sumIw2[thermit.first], thermit.first); + domain->setLocalNrotDOF(thermit.first, N[thermit.first], rotDOF[thermit.first]); + } } void Langevin::eventNewTimestep(ParticleContainer *particleContainer, Domain *domain) { - for(std::size_t index = 0; index < _stochastic_regions.size(); index++) { - auto& region = _stochastic_regions[index]; - double T = _simulation.getTemperatureObserver()->getTemperature(index); - T = _simulation.getEnsemble()->T(); - #if defined(_OPENMP) - #pragma omp parallel - #endif - for(auto it = particleContainer->regionIterator(std::data(region.low), - std::data(region.high), - ParticleIterator::ONLY_INNER_AND_BOUNDARY); it.isValid(); ++it) { - d3 v_old = it->v_arr(); - d3 v_rand = sampleRandomForce(it->mass(), T); - d3 v_delta {}; - - for(int d = 0; d < 3; d++) { - v_delta[d] = - _dt_half * _xi * v_old[d] + v_rand[d]; - } - - it->vadd(v_delta[0], v_delta[1], v_delta[2]); - } - } - - #if defined(_OPENMP) - #pragma omp parallel - #endif - { - for (auto i = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); i.isValid(); ++i) { - i->upd_preF(_timestepLength); - } - } + addLangevinContribution(particleContainer); + + #if defined(_OPENMP) + #pragma omp parallel + #endif + { + for (auto i = particleContainer->iterator(ParticleIterator::ONLY_INNER_AND_BOUNDARY); i.isValid(); ++i) { + i->upd_preF(_timestepLength); + } + } } Langevin::d3 Langevin::sampleRandomForce(double m, double T) { - static thread_local std::mt19937 generator; // keeping default seed for reproducibility - std::normal_distribution normal{0.0, 1.0}; - d3 r_vec {}; - - // additional factor 2 here - // according to book about Langevin integration not needed, but according to Langevin's equations of motion it does exist - // by using it we actually reach the target temperature - double scale = std::sqrt(2 * _timestepLength * T * _xi / m); - for(int d = 0; d < 3; d++) { - r_vec[d] = scale * normal(generator); - } - - return r_vec; + static thread_local std::mt19937 generator; // keeping default seed for reproducibility + std::normal_distribution normal{0.0, 1.0}; + d3 r_vec{}; + + // additional factor 2 here + // according to book about Langevin integration not needed, but according to Langevin's equations of motion it does exist + // by using it we actually reach the target temperature + double scale = std::sqrt(2 * _timestepLength * T * _gamma / m); + for (int d = 0; d < 3; d++) { + r_vec[d] = scale * normal(generator); + } + + return r_vec; +} + +void Langevin::addLangevinContribution(ParticleContainer *particleContainer) { + for (std::size_t index = 0; index < _stochastic_regions.size(); index++) { + auto ®ion = _stochastic_regions[index]; + double T = _simulation.getEnsemble()->T(); + #if defined(_OPENMP) + #pragma omp parallel + #endif + for (auto it = particleContainer->regionIterator(std::data(region.low), + std::data(region.high), + ParticleIterator::ONLY_INNER_AND_BOUNDARY); it.isValid(); ++it) { + d3 v_old = it->v_arr(); + d3 v_rand = sampleRandomForce(it->mass(), T); + d3 v_delta{}; + + for (int d = 0; d < 3; d++) { + v_delta[d] = -_dt_half * _gamma * v_old[d] + v_rand[d]; + } + + it->vadd(v_delta[0], v_delta[1], v_delta[2]); + } + } } diff --git a/src/integrators/Langevin.h b/src/integrators/Langevin.h index cef1b07efe..5061786893 100644 --- a/src/integrators/Langevin.h +++ b/src/integrators/Langevin.h @@ -21,40 +21,62 @@ * */ class Langevin : public Integrator { public: - void readXML(XMLfileUnits &xmlconfig) override; + /** + * XML Format: + * + * DOUBLE + * DOUBLE + * + * + * DOUBLE DOUBLE DOUBLE + * DOUBLE DOUBLE DOUBLE + * + * + * + * */ + void readXML(XMLfileUnits &xmlconfig) override; - void init() override; + void init() override; - void eventNewTimestep(ParticleContainer *moleculeContainer, Domain *domain) override; + void eventNewTimestep(ParticleContainer *moleculeContainer, Domain *domain) override; + + void eventForcesCalculated(ParticleContainer *moleculeContainer, Domain *domain) override; - void eventForcesCalculated(ParticleContainer *moleculeContainer, Domain *domain) override; private: - using d3 = std::array; - struct box_t { - box_t() = default; - box_t(const d3& l, const d3& h) : low(l), high(h) {} - d3 low, high; - }; - - //! @brief friction strength - double _xi; - - //! @brief time step length halved - double _dt_half; - - //! @brief regions in which friction is not set to 0 - std::vector _stochastic_regions; - - //! @brief We need to check at least twice if a TemperatureObserver exists - bool _checkFailed; - - /** - * Sample from Gaussian with technically 0 mean and sigma**2 = 2 * m * friction * k_b * T_target / delta_t. - * Is already adapted to match the integration scheme. - * @param m mass - * @param T temp target - * */ - d3 sampleRandomForce(double m, double T); + using d3 = std::array; + + struct box_t { + box_t() = default; + + box_t(const d3 &l, const d3 &h) : low(l), high(h) {} + + d3 low, high; + }; + + //! @brief friction strength + double _gamma; + + //! @brief time step length halved + double _dt_half; + + //! @brief regions in which friction is not set to 0 + std::vector _stochastic_regions; + + //! @brief We need to check at least twice if a TemperatureObserver exists + bool _checkFailed; + + /** + * Sample from Gaussian with technically 0 mean and sigma**2 = 2 * m * friction * k_b * T_target / delta_t. + * Is already adapted to match the integration scheme. + * @param m mass + * @param T temp target + * */ + d3 sampleRandomForce(double m, double T); + + /** + * Adds the parts of the equation, that only exist in the Langevin Equation's of motion + * */ + void addLangevinContribution(ParticleContainer* particleContainer); }; diff --git a/src/thermostats/TemperatureObserver.cpp b/src/thermostats/TemperatureObserver.cpp index b3477bc104..34f8143e92 100644 --- a/src/thermostats/TemperatureObserver.cpp +++ b/src/thermostats/TemperatureObserver.cpp @@ -7,144 +7,145 @@ #include "ensemble/EnsembleBase.h" void TemperatureObserver::readXML(XMLfileUnits &xmlconfig) { - std::size_t num_regions = 0; - XMLfile::Query query = xmlconfig.query("ActiveRegions/region"); - num_regions = query.card(); - _region_data.resize(num_regions); - - std::string oldpath = xmlconfig.getcurrentnodepath(); - XMLfile::Query::const_iterator rIt; - std::size_t index = 0; - for(rIt = query.begin(); rIt; rIt++) { - xmlconfig.changecurrentnode(rIt); - d3 low {}, high{}; - xmlconfig.getNodeValue("lowX", low[0]); - xmlconfig.getNodeValue("lowY", low[1]); - xmlconfig.getNodeValue("lowZ", low[2]); - - xmlconfig.getNodeValue("highX", high[0]); - xmlconfig.getNodeValue("highY", high[1]); - xmlconfig.getNodeValue("highZ", high[2]); - - _region_data[index].first.low = low; - _region_data[index].first.high = high; - index++; - } - xmlconfig.changecurrentnode(oldpath); - - _max_samples = 100; - xmlconfig.getNodeValue("samples", _max_samples); + std::size_t num_regions = 0; + XMLfile::Query query = xmlconfig.query("ActiveRegions/region"); + num_regions = query.card(); + _region_data.resize(num_regions); + + std::string oldpath = xmlconfig.getcurrentnodepath(); + XMLfile::Query::const_iterator rIt; + std::size_t index = 0; + for (rIt = query.begin(); rIt; rIt++) { + xmlconfig.changecurrentnode(rIt); + d3 low{}, high{}; + xmlconfig.getNodeValue("lower/x", low[0]); + xmlconfig.getNodeValue("lower/y", low[1]); + xmlconfig.getNodeValue("lower/z", low[2]); + + xmlconfig.getNodeValue("upper/x", high[0]); + xmlconfig.getNodeValue("upper/y", high[1]); + xmlconfig.getNodeValue("upper/z", high[2]); + + _region_data[index].first.low = low; + _region_data[index].first.high = high; + index++; + } + xmlconfig.changecurrentnode(oldpath); + + _max_samples = 100; + xmlconfig.getNodeValue("samples", _max_samples); } void TemperatureObserver::init() { - for(auto& [region, region_cache] : _region_data) { - region_cache.history.init(_max_samples); - region_cache.current_temp = 0.0; - } + for (auto &[region, region_cache]: _region_data) { + region_cache.history.init(_max_samples); + region_cache.current_temp = 0.0; + } } -void TemperatureObserver::step(ParticleContainer* particleContainer) { - for(auto& [region, region_cache] : _region_data) { - region_cache.current_temp = measureTemp(region.low, region.high, particleContainer, region_cache.history); - } +void TemperatureObserver::step(ParticleContainer *particleContainer) { + for (auto &[region, region_cache]: _region_data) { + region_cache.current_temp = measureTemp(region.low, region.high, particleContainer, region_cache.history); + } } -void TemperatureObserver::getRegions(std::vector, std::array>>& regions) { - for(auto& [region, region_cache] : _region_data) { - regions.emplace_back(region.low, region.high); - } +void TemperatureObserver::getRegions(std::vector, std::array>> ®ions) { + for (auto &[region, region_cache]: _region_data) { + regions.emplace_back(region.low, region.high); + } } double TemperatureObserver::getTemperature(std::size_t index) { - mardyn_assert(index < _region_data.size() && index >= 0); - return _region_data[index].second.current_temp; + mardyn_assert(index < _region_data.size() && index >= 0); + return _region_data[index].second.current_temp; } double TemperatureObserver::measureTemp(const std::array &low, const std::array &high, - ParticleContainer *particleContainer, BinData& binData) { - std::array v_mean = computeMeanVelocityStep(low, high, particleContainer, binData); - double m = _simulation.getEnsemble()->getComponent(0)->m(); - double v = 0; - std::size_t n = 0; -#if defined(_OPENMP) -#pragma omp parallel reduction(+:v, n) -#endif - for (auto it = particleContainer->regionIterator(std::data(low), std::data(high), - ParticleIterator::ALL_CELLS); it.isValid(); ++it) { - for(int d = 0; d < 3; d++) { - v += std::pow(it->v(d) - v_mean[d], 2); - } - n += 1; - } - - return v * m / (3 * static_cast(n)); + ParticleContainer *particleContainer, BinData &binData) { + std::array v_mean = computeMeanVelocityStep(low, high, particleContainer, binData); + double m = _simulation.getEnsemble()->getComponent(0)->m(); + double v = 0; + std::size_t n = 0; + #if defined(_OPENMP) + #pragma omp parallel reduction(+:v, n) + #endif + for (auto it = particleContainer->regionIterator(std::data(low), std::data(high), + ParticleIterator::ALL_CELLS); it.isValid(); ++it) { + for (int d = 0; d < 3; d++) { + v += std::pow(it->v(d) - v_mean[d], 2); + } + n += 1; + } + + return v * m / (3 * static_cast(n)); } -std::array TemperatureObserver::computeMeanVelocityStep(const std::array &low, const std::array &high, - ParticleContainer *particleContainer, BinData& binData) { - //compute mean velocity per dim using SMC method, see: Karimian et al. 2011 - std::array v_avg{0}; - double* v_raw = std::data(v_avg); - std::size_t count = 0; -#if defined(_OPENMP) -#pragma omp parallel reduction(+:v_raw[:3], count) -#endif - for (auto it = particleContainer->regionIterator(std::data(low), std::data(high), - ParticleIterator::ALL_CELLS); it.isValid(); ++it) { - for(int d = 0; d < 3; d++) { - v_raw[d] += it->v(d); - } - count += 1; - } - //TODO check this - if(count == 0) count = 1; - binData.mol_counts.insert(count); - binData.velocities.insert(v_avg); - - //compute CAM average first - { - std::size_t steps = binData.mol_counts.size(); - auto d = static_cast(steps<=1?1:steps-1); - - std::array v_cam{ 0 }; - for(std::size_t i = 0; i < steps; i++) { - const auto& vec = binData.velocities.get(i); - for(int dim = 0; dim < 3; dim++) { - v_cam[dim] += vec[dim]; - } - } - for(int dim = 0; dim < 3; dim++) { - v_cam[dim] /= d; - } - - std::size_t nks = 0; - for(std::size_t i = 0; i < steps; i++) { - nks += binData.mol_counts.get(i); - } - double nks_d = static_cast(nks) / d; - - for(int dim = 0; dim < 3; dim++) { - v_cam[dim] /= nks_d; - } - binData.v_cam.insert(v_cam); - } - - //compute SAM over CAM - { - std::size_t steps = binData.v_cam.size(); - auto d = static_cast(steps<=1?1:steps-1); - - std::array v_scm{ 0 }; - for(std::size_t i = 0; i < steps; i++) { - const auto& vec = binData.v_cam.get(i); - for(int dim = 0; dim < 3; dim++) { - v_scm[dim] += vec[dim]; - } - } - for(int dim = 0; dim < 3; dim++) { - v_scm[dim] /= d; - } - - return v_scm; - } +std::array +TemperatureObserver::computeMeanVelocityStep(const std::array &low, const std::array &high, + ParticleContainer *particleContainer, BinData &binData) { + //compute mean velocity per dim using SMC method, see: Karimian et al. 2011 + std::array v_avg{0}; + double *v_raw = std::data(v_avg); + std::size_t count = 0; + #if defined(_OPENMP) + #pragma omp parallel reduction(+:v_raw[:3], count) + #endif + for (auto it = particleContainer->regionIterator(std::data(low), std::data(high), + ParticleIterator::ALL_CELLS); it.isValid(); ++it) { + for (int d = 0; d < 3; d++) { + v_raw[d] += it->v(d); + } + count += 1; + } + // setting to 1 in case there were no molecules to not divide by zero; will not change result as v is still zero + if (count == 0) count = 1; + binData.mol_counts.insert(count); + binData.velocities.insert(v_avg); + + //compute CAM average first + { + std::size_t steps = binData.mol_counts.size(); + auto d = static_cast(steps <= 1 ? 1 : steps - 1); + + std::array v_cam{0}; + for (std::size_t i = 0; i < steps; i++) { + const auto &vec = binData.velocities.get(i); + for (int dim = 0; dim < 3; dim++) { + v_cam[dim] += vec[dim]; + } + } + for (int dim = 0; dim < 3; dim++) { + v_cam[dim] /= d; + } + + std::size_t nks = 0; + for (std::size_t i = 0; i < steps; i++) { + nks += binData.mol_counts.get(i); + } + double nks_d = static_cast(nks) / d; + + for (int dim = 0; dim < 3; dim++) { + v_cam[dim] /= nks_d; + } + binData.v_cam.insert(v_cam); + } + + //compute SAM over CAM + { + std::size_t steps = binData.v_cam.size(); + auto d = static_cast(steps <= 1 ? 1 : steps - 1); + + std::array v_scm{0}; + for (std::size_t i = 0; i < steps; i++) { + const auto &vec = binData.v_cam.get(i); + for (int dim = 0; dim < 3; dim++) { + v_scm[dim] += vec[dim]; + } + } + for (int dim = 0; dim < 3; dim++) { + v_scm[dim] /= d; + } + + return v_scm; + } } diff --git a/src/thermostats/TemperatureObserver.h b/src/thermostats/TemperatureObserver.h index 269aa5648c..5fb5465182 100644 --- a/src/thermostats/TemperatureObserver.h +++ b/src/thermostats/TemperatureObserver.h @@ -14,7 +14,9 @@ class Domain; + class ParticleContainer; + class Ensemble; /** @@ -22,66 +24,68 @@ class Ensemble; * */ template struct NBuffer { - NBuffer() = default; - [[maybe_unused]] explicit NBuffer(std::size_t n) : _n(n), _begin(-1), _end(0) { - _data.resize(n); - } - void init(std::size_t n) { - _n = n; - _begin = -1; - _end = 0; - _data.resize(n); - } - void insert(const T& t) { - if(_begin == -1) { - _data[0] = t; - _begin = 0; - _end = 1; - } - else if (_begin == _end) { - _data[_end] = t; - _begin = (_begin+1) % _n; - _end = _begin; - } - else { - _data[_end] = t; - _end = (_end+1) % _n; - } - } - const T& get(std::size_t index) { - mardyn_assert(index < _n); - return _data[(_begin + index) % _n]; - } - std::size_t size() { - if(_begin == -1) { - return 0; - } - else if (_begin == _end) { - return _n; - } - else { - return _end - _begin; - } - } - std::size_t _n{}; - long _begin{}; - long _end{}; - std::vector _data; + NBuffer() = default; + + [[maybe_unused]] explicit NBuffer(std::size_t n) : _n(n), _begin(-1), _end(0) { + _data.resize(n); + } + + void init(std::size_t n) { + _n = n; + _begin = -1; + _end = 0; + _data.resize(n); + } + + void insert(const T &t) { + if (_begin == -1) { + _data[0] = t; + _begin = 0; + _end = 1; + } else if (_begin == _end) { + _data[_end] = t; + _begin = (_begin + 1) % _n; + _end = _begin; + } else { + _data[_end] = t; + _end = (_end + 1) % _n; + } + } + + const T &get(std::size_t index) { + mardyn_assert(index < _n); + return _data[(_begin + index) % _n]; + } + + std::size_t size() { + if (_begin == -1) { + return 0; + } else if (_begin == _end) { + return _n; + } else { + return _end - _begin; + } + } + + std::size_t _n{}; + long _begin{}; + long _end{}; + std::vector _data; }; /** * Data struct for computing SAM and CAM averages * */ struct BinData { - NBuffer mol_counts; - NBuffer> velocities; - NBuffer> v_cam; - - void init(std::size_t n) { - mol_counts.init(n); - velocities.init(n); - v_cam.init(n); - } + NBuffer mol_counts; + NBuffer> velocities; + NBuffer> v_cam; + + void init(std::size_t n) { + mol_counts.init(n); + velocities.init(n); + v_cam.init(n); + } }; /** @@ -89,54 +93,75 @@ struct BinData { * */ class TemperatureObserver { public: - /** - * Init all buffers - * */ - void init(); - void readXML(XMLfileUnits &xmlconfig); - /** - * Measure temp in all regions, updates caches and stores current temp - * */ - void step(ParticleContainer* particleContainer); - - /** - * @returns all regions, for which the temp is being measured. - * */ - void getRegions(std::vector, std::array>>& regions); - - /** - * Gets the temperature for the selected region. - * Regions can be queried by accessing getRegions(). - * @returns T * k_b - * */ - double getTemperature(std::size_t index); + /** + * Init all buffers + * */ + void init(); + + /** + * XML Format: + * + * DOUBLE + * + * + * DOUBLE DOUBLE DOUBLE + * DOUBLE DOUBLE DOUBLE + * + * + * + * */ + void readXML(XMLfileUnits &xmlconfig); + + /** + * Measure temp in all regions, updates caches and stores current temp + * */ + void step(ParticleContainer *particleContainer); + + /** + * @returns all regions, for which the temp is being measured. + * */ + void getRegions(std::vector, std::array>> ®ions); + + /** + * Gets the temperature for the selected region. + * Regions can be queried by accessing getRegions(). + * @returns T * k_b + * */ + double getTemperature(std::size_t index); + private: - std::size_t _max_samples; - - using d3 = std::array; - struct box_t { - box_t() = default; - box_t(const d3& l, const d3& h) : low(l), high(h) {} - d3 low, high; - }; - struct box_data_t { - BinData history{ }; - double current_temp{ }; - }; - - //! @brief each region is defined by a box in 3d space (box_t) and a data cache (box_data_t) - std::vector> _region_data; - - /** - * Based on: Details about pressure calculation in molecular dynamic analysis by Karimian et al. 2014 - * returns T * k_b - * */ - double measureTemp(const std::array& low, const std::array& high, ParticleContainer *particleContainer, BinData& binData); - - /** - * Computes the SAM over the CAM average of the mean velocity for the specified region. (SCM method) - * */ - std::array computeMeanVelocityStep(const std::array& low, const std::array& high, ParticleContainer *particleContainer, BinData& binData); + std::size_t _max_samples; + + using d3 = std::array; + + struct box_t { + box_t() = default; + + box_t(const d3 &l, const d3 &h) : low(l), high(h) {} + + d3 low, high; + }; + + struct box_data_t { + BinData history{}; + double current_temp{}; + }; + + //! @brief each region is defined by a box in 3d space (box_t) and a data cache (box_data_t) + std::vector> _region_data; + + /** + * Based on: Details about pressure calculation in molecular dynamic analysis by Karimian et al. 2014 + * returns T * k_b + * */ + double measureTemp(const std::array &low, const std::array &high, + ParticleContainer *particleContainer, BinData &binData); + + /** + * Computes the SAM over the CAM average of the mean velocity for the specified region. (SCM method) + * */ + std::array computeMeanVelocityStep(const std::array &low, const std::array &high, + ParticleContainer *particleContainer, BinData &binData); }; #endif //MARDYN_TEMPERATUREOBSERVER_H