Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Langevin thermostat #296

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Langevin thermostat #296

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
10e6c5d
Langevin Thermo/Integrator implemented
Snapex2409 Mar 6, 2024
c5e0029
Target Temp is now reached
Snapex2409 Mar 8, 2024
0e2a9d6
Fixed review items
Snapex2409 Mar 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 30 additions & 8 deletions src/Simulation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -115,6 +116,7 @@ Simulation::Simulation()
_rand(8624),
_longRangeCorrection(nullptr),
_temperatureControl(nullptr),
_temperatureObserver(nullptr),
Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved
_FMM(nullptr),
_timerProfiler(),
#ifdef TASKTIMINGPROFILE
Expand Down Expand Up @@ -150,6 +152,8 @@ Simulation::~Simulation() {
_longRangeCorrection = nullptr;
delete _temperatureControl;
_temperatureControl = nullptr;
delete _temperatureObserver;
_temperatureObserver = nullptr;
Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved
delete _FMM;
_FMM = nullptr;

Expand Down Expand Up @@ -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 {
Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved
Log::global_log-> error() << "Unknown integrator " << integratorType << std::endl;
Simulation::exit(1);
}
Expand Down Expand Up @@ -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 {
Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved
Log::global_log->warning() << "Unknown thermostat " << thermostattype << std::endl;
continue;
}
Expand Down Expand Up @@ -810,6 +824,10 @@ void Simulation::prepare_start() {
_cellProcessor = new LegacyCellProcessor( _cutoffRadius, _LJCutoffRadius, _particlePairsHandler);
}

if(dynamic_cast<Langevin*>(_integrator) != nullptr) {
_integrator->init();
}

Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved
if (_FMM != nullptr) {

double globalLength[3];
Expand Down Expand Up @@ -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);
}

Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved
// initializing plugins and starting plugin timers
for (auto& plugin : _plugins) {
Log::global_log->info() << "Initializing plugin " << plugin->getPluginName() << std::endl;
Expand Down Expand Up @@ -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) {
Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved
if (_thermostatType ==VELSCALE_THERMOSTAT) {
Log::global_log->debug() << "Velocity scaling" << std::endl;
if (_domain->severalThermostats()) {
Expand Down Expand Up @@ -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);
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not 100% sure why this temperature Observer is exactly needed for the moment...

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  • It was initially created due to my misunderstanding of the formula. I thought we needed the current temperature in a specific region. However, that is not the case.
  • But since we still require a method of disabling all other thermostat functionality, it seemed viable to keep the Temperature Observer as a "thermostat" that does nothing. The other option would be using an existing thermostat and setting certain values to zero, which seemed to me to be prone to misconfiguration.


advanceSimulationTime(_integrator->getTimestepLength());

Expand Down
7 changes: 7 additions & 0 deletions src/Simulation.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@

// plugins
#include "plugins/PluginFactory.h"
#include "thermostats/TemperatureObserver.h"

#if !defined (SIMULATION_SRC) or defined (IN_IDE_PARSER)
class Simulation;
Expand Down Expand Up @@ -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; }

Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved
/** 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 */
Expand Down Expand Up @@ -391,6 +395,9 @@ class Simulation {
/** Temperature Control (Slab Thermostat) */
TemperatureControl* _temperatureControl;

/** No Thermostat - only measures temp in selected regions */
TemperatureObserver* _temperatureObserver;

Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved
/** The Fast Multipole Method object */
bhfmm::FastMultipoleMethod* _FMM;

Expand Down
190 changes: 190 additions & 0 deletions src/integrators/Langevin.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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);

_xi = 0;
Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved
xmlconfig.getNodeValue("friction", _xi);
mardyn_assert(_xi > 0);
}

void Langevin::init() {
std::vector<std::pair<std::array<double, 3>, std::array<double, 3>>> 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);
}
}

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<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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the following shall be OpenMP-parallelized I miss some worksharing directives, e.g., #omp for ...

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I understand the regionIterator correctly, then no further OMP directives are needed, as everything is handled by the iterator.

{
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) {
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 {};

// 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;
}
61 changes: 61 additions & 0 deletions src/integrators/Langevin.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
//
// Created by Alex Hocks on 05.03.24.
//

#ifndef MARDYN_LANGEVIN_H
#define MARDYN_LANGEVIN_H


#include "Integrator.h"

#include <array>
#include <vector>

/**
* 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.
Comment on lines +15 to +20
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would be good to have a reference of a paper/some material to make this easier to follow.

* */
class Langevin : public Integrator {
public:
void readXML(XMLfileUnits &xmlconfig) override;
Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved

void init() override;

void eventNewTimestep(ParticleContainer *moleculeContainer, Domain *domain) override;

void eventForcesCalculated(ParticleContainer *moleculeContainer, Domain *domain) override;
private:
using d3 = std::array<double, 3>;
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;
Snapex2409 marked this conversation as resolved.
Show resolved Hide resolved

//! @brief time step length halved
double _dt_half;

//! @brief regions in which friction is not set to 0
std::vector<box_t> _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);
};


#endif //MARDYN_LANGEVIN_H
Loading
Loading