Langevin thermostat

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"
@@ -68,6 +69,7 @@
 #include "ensemble/CavityEnsemble.h"
 
 #include "thermostats/VelocityScalingThermostat.h"
+#include "thermostats/TemperatureObserver.h"
 #include "thermostats/TemperatureControl.h"
 
 #include "utils/FileUtils.h"
@@ -115,6 +117,7 @@ Simulation::Simulation()
 	_rand(8624),
 	_longRangeCorrection(nullptr),
 	_temperatureControl(nullptr),
+	_temperatureObserver(nullptr),
 	_FMM(nullptr),
 	_timerProfiler(),
 #ifdef TASKTIMINGPROFILE
@@ -150,6 +153,8 @@ Simulation::~Simulation() {
 	_longRangeCorrection = nullptr;
 	delete _temperatureControl;
 	_temperatureControl = nullptr;
+	delete _temperatureObserver;
+	_temperatureObserver = nullptr;
 	delete _FMM;
 	_FMM = nullptr;
 
@@ -182,6 +187,9 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) {
 			_integrator = new Leapfrog();
 		} else if (integratorType == "LeapfrogRMM") {
 			_integrator = new LeapfrogRMM();
+		} else if (integratorType == "Langevin") {
+			_integrator = new Langevin();
+			_thermostatType = LANGEVIN_THERMOSTAT;
 		} else {
 			Log::global_log-> error() << "Unknown integrator " << integratorType << std::endl;
 			Simulation::exit(1);
@@ -521,17 +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);
-                    }
-                }
-				else
-				{
+					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 {
 					Log::global_log->warning() << "Unknown thermostat " << thermostattype << std::endl;
 					continue;
 				}
@@ -810,6 +826,10 @@ void Simulation::prepare_start() {
 		_cellProcessor = new LegacyCellProcessor( _cutoffRadius, _LJCutoffRadius, _particlePairsHandler);
 	}
 
+	if(dynamic_cast<Langevin*>(_integrator) != nullptr) {
+		_integrator->init();
+	}
+
 	if (_FMM != nullptr) {
 
 		double globalLength[3];
@@ -898,6 +918,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;
@@ -1201,13 +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);
 	}
-	// <-- TEMPERATURE_CONTROL
-
-
 
 	advanceSimulationTime(_integrator->getTimestepLength());
 

src/Simulation.h

@@ -54,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;
@@ -226,6 +228,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 +396,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;
 

src/integrators/Langevin.cpp

@@ -0,0 +1,190 @@
+//
+// 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 <random>
+
+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);
+
+	_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() {
+	_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) {
+	addLangevinContribution(particleContainer);
+
+	std::map<int, unsigned long> N;
+	std::map<int, unsigned long> rotDOF;
+	std::map<int, double> summv2;
+	std::map<int, double> sumIw2;
+
+	if (domain->severalThermostats()) {
+		#if defined(_OPENMP)
+		#pragma omp parallel
+		#endif
+		{
+			std::map<int, unsigned long> N_l;
+			std::map<int, unsigned long> rotDOF_l;
+			std::map<int, double> summv2_l;
+			std::map<int, double> 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) {
+	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 * _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 &region = _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]);
+		}
+	}
+}