Source code for ensembler.conditions.thermostats

import warnings

import numpy as np
from scipy import constants as const
from scipy.stats import maxwell

from ensembler.conditions._basicCondition import _conditionCls
from ensembler.util.ensemblerTypes import systemCls as systemType


[docs]class thermostat(_conditionCls): """ Thermostat This is the parent class of Thermostats. The apply function of this class is not implemented and needs to be overwritten by each subclass. """ _currentTemperature: float def __init__(self, system: systemType, tau: int, verbose:bool=False): super().__init__(system=system, tau=tau, verbose=verbose)
[docs]class andersonThermostat(thermostat): """ -UnderConstuction- andersenThermostat This thermostat was developed by anderson 1980 and is also called stochastic collisions method. This method is optimal for single particle simulations as it tries to simulate a particle that collides with virtual particles from a temperature bath. At each collision, the velocity and the interval of the next collision is randomly reassigned and reselected from a Maxwell-Boltzmann distribution. This applying the thermostat leads to an at least microcanonical Simulation. Simulated energies should be gaussian distributed. [Molecular Modelling Principles and Applications, A. R. Leach, second edition] FIX MAXWWELL POSITIONING! """ # collision parameters kb: float = const.k * const.Avogadro a: float = 1 # dimensionless constant k: float = 1 # thermal konductivity N_dens: float = 0.1 # particle density N: float = 10 # particle # temperaturebath: temperature: float # desired temperature area temperature_noise_range: float _new_temperature: float _lambda: float = 1 # scaling factor for velocities """Under COnstructurion"""
[docs] def __init__(self, temperature: float = 298, temperature_noise_range: float = 25, MConstraintsDims: int = 1, system: systemType = None, tau:int=1, kb=const.k * const.Avogadro, a: float = 1, k: float = 1, N_dens: float = 0.005, N: float = 1, verbose: bool = False): """ __This Thermostat is underconstruction! __ Parameters ---------- temperature : float, optional desired temperature temperature_noise_range : float, optional noise range MConstraintsDims : system : system, optional a system, that shall be thermostated tau : int, optional every n step apply thermostat kb : float, optional boltzman constant a : float, optional collision parameter k : float, optional collision parameter N_dens : float, optional virtual particle density N : int, optional virtual particle number verbose : bool, optional I can be loud and noisy! """ warnings.warn("__Under construction___!") super().__init__(system=system, tau=tau, verbose=verbose) # Collision parameters self.kb = kb self.a = a # dimensionless constant self.k = k # thermal konductivity self.N_dens = N_dens # particle density self.N = N # particle # Temperature Scalingparameters self.temperature_noise_range = temperature_noise_range self.M = MConstraintsDims if (system != None): self.temperature = self.system.temperature else: self.temperature = temperature
def _collision(self) -> bool: p_collision = (2 * self.a * self.k) / (3 * self.kb * self.N_dens ** (1 / 3) * self.N ** (2 / 3)) options = [True, False] collision = np.random.choice(a=options, p=(p_collision, 1 - p_collision)) return collision def _rescale_velocities(self): orig_vels = self.system._currentVelocities new_vels = self._lambda * orig_vels if (self._lambda * orig_vels) else 0.0001 # do not allow 0 vel self.system._currentVelocities = new_vels def _calculate_scaling_factor(self): # pick new temperature randomly self._new_temperature = maxwell.rvs(loc=self.temperature - self.temperature_noise_range, scale=self.temperature_noise_range) self._currentTemperature = ( self.system._currentTemperature if (self.system._currentTemperature != 0) else 0.000001) self.system._currentTemperature = self._new_temperature # scaling factor self._lambda = self._new_temperature / self._currentTemperature # (1+(self.dt/self.tau)*((self.system.temperature/T_t)-1))**0.5
[docs] def apply_coupled(self): if (self._collision()): self._calculate_scaling_factor() self._rescale_velocities() if (self.verbose): print("THERMOSTAT: get to temp: ", self.system._currentTemperature, "\n" 'THERMOSTAT: tot_kin: ', self.system.calculate_total_kinetic_energy(), "\n" "THERMOSTAT: lambda: ", self._lambda, "\n" "THERMOSTAT: current_Velocity: ", self.system._currentVelocities, "\n" "\n") else: pass
[docs]class berendsenThermostate(thermostat):
[docs] def __init__(self, tau: float, dt: float, MConstraintsDims: int = 1, system: systemType = None, verbose:bool=False): """ __under contsruction! __ This thermostat is not tested, no guarantee for correctness! reference: Molecular dynamics with coupling to an external bath; H.J.C. Berendsen Parameters ---------- tau : int, optional apply every n steps. dt : float, optional time step MConstraintsDims : int, optional number of dimensions system : system, optional a target system to be thermostated verbose : bool, optional More output? You want more output? """ super().__init__(system=system, tau=tau, verbose=verbose) self._lambda: float = 1 # scaling factor of velocities self._current_temperatur = 1 self.dt = dt self.M = MConstraintsDims
[docs] def apply_coupled(self): self._calculate_current_temperature() self._calculate_scaling_factor() self._rescale_velocities() if (self.verbose): print("THERMOSTAT: get to temp: ", self.system.temperature, "\n" 'THERMOSTAT: tot_kin: ', self.system.calculate_total_kinetic_energy(), "\n" "THERMOSTAT: curr temp: ", self._current_temperatur, "\n" "THERMOSTAT: lambda: ", self._lambda, "\n" "THERMOSTAT: current_Velocity: ", self.system._currentVelocities, "\n" "\n")
def _rescale_velocities(self): orig_vels = self.system._currentVelocities new_vels = abs(self._lambda) * orig_vels if (self._lambda * orig_vels) else 0.0001 # do not allow 0 vel self.system._currentVelocities = new_vels
[docs] def _calculate_current_temperature(self): """ autofunction: _calculate current Temperature (eq. 32,33) :return: """ # M = constraints N = self.system.nparticles * const.Avogadro self._current_temperatur = (2 / (3 * N - self.M - 3)) * self.system.calculate_total_kinetic_energy() * N self.system._currentTemperature = self._current_temperatur
[docs] def _calculate_scaling_factor(self): """ autofunction: _calculate_scaling_factor (eq.34) :return: """ T_t = (self._current_temperatur if (self._current_temperatur != 0) else 0.000001) self._lambda = (1 + (self.dt / self.tau) * ((self.system.temperature / T_t) - 1)) ** 0.5
""" class nhIntegrator(samplers): Nosehover-leapfrog def scaleVel(self, sys): freetemp = 2.0 / const.gas_constant / 1000.0 * sys.mu * sys.new_velocity ** 2 # t+0.5Dt self.oldxi = self.xi # t-0.5Dt self.xi += self.dt / (self.tau * self.tau) * (freetemp / sys.temp - 1.0) # t+0.5t scale = 1.0 - self.xi * self.dt return scale def step(self, sys): sys.pos = sys.newpos # t sys.force = -sys.pot.dhdpos(sys.lam, sys.pos) # t sys.oldvel = sys.new_velocity # t - 0.5 Dt sys.new_velocity += sys.force / sys.mu * self.dt # t+0.5t sys.new_velocity *= self.scaleVel(sys) sys.vel = 0.5 * (sys.oldvel + sys.new_velocity) sys.newpos += self.dt * sys.new_velocity # t+Dt sys.veltemp = sys.mu / const.gas_constant/1000.0 * sys.vel ** 2 # t sys.updateEne() return [sys.pos, sys.vel, sys.veltemp, sys.totkin, sys.totpot, sys.totene, sys.lam, sys.dhdlam] def __init__(self, xi=0.0, tau=0.1, dt=1e-1): self.xi = xi self.oldxi = xi self.tau = tau self.dt = dt raise NotImplementedError("This "+__class__+" class is not implemented") class hmcIntegrator(samplers): def step(self, sys): accept = 0 oldene = sys.totene oldpos = sys.pos # t oldvel = sys.vel sys.initVel() # t-0.5Dt for i in range(self.steps): sys.pos += self.dt * sys.new_velocity # t force = -sys.pot.dhdpos(sys.lam, sys.pos) # t sys.vel = (oldvel + sys.new_velocity) / 2.0 sys.veltemp = sys.mu / const.gas_constant / 1000.0 * sys.vel ** 2 # t sys.updateEne() oldvel = sys.new_velocity sys.new_velocity += force / sys.mu * self.dt # t+0.5t accept = 0 if sys.totene < oldene: accept = 1 else: if np.random.rand() <= np.exp(-1 / (const.gas_constant / 1000.0* sys.temp) * (sys.totene - oldene)): accept = 1 if accept == 0: sys.pos = oldpos sys.vel = oldvel sys.veltemp = sys.mu / const.gas_constant * sys.vel ** 2 # t sys.updateEne() return [sys.pos, sys.vel, sys.veltemp, sys.totkin, sys.totpot, sys.totene, sys.lam, sys.dhdlam] def __init__(self, steps=5, dt=1e-1): self.steps = steps self.dt = dt raise NotImplementedError("This "+__class__+" class is not implemented") """