eonc::Dynamics Class Reference

#include <Dynamics.h>

Public Member Functions
	Dynamics (Matter *matter, const DynamicsConfig &config)
	Dynamics (Matter *matter, const Parameters &parameters)
	~Dynamics ()
void	setTemperature (double temperature)
void	oneStep (int stepNumber=-1)
void	velocityVerlet ()
void	run ()
void	andersenCollision ()
void	setThermalVelocity ()
void	rescaleVelocity ()
void	noseHooverVerlet ()
	Nose-Hoover chain thermostat (Martyna-Klein-Tuckerman algorithm).
void	langevinVerlet ()
	Langevin dynamics (velocity-Verlet with friction and random forces).

Static Public Attributes
static const char	ANDERSEN [] = "andersen"
static const char	NOSE_HOOVER [] = "nose_hoover"
static const char	LANGEVIN [] = "langevin"
static const char	NONE [] = "none"

Private Attributes
long	nAtoms {0}
long	nFreeCoords {0}
Matter *	matter
DynamicsConfig	m_config
double	dt {0.0}
double	kB {0.0}
double	temperature {0.0}
double	vxi1 {0.0}
double	vxi2 {0.0}
double	xi1 {0.0}
double	xi2 {0.0}
eonc::log::Scoped	log

Detailed Description

Definition at line 54 of file Dynamics.h.

Constructor & Destructor Documentation

Dynamics() [1/2]

Dynamics::Dynamics	(	Matter *	matter,
		const DynamicsConfig &	config )

Definition at line 24 of file Dynamics.cpp.

    : matter{matter_in},
      m_config{config} {
  dt = m_config.time_step;
  nAtoms = matter->numberOfAtoms();
  nFreeCoords = matter->numberOfFreeAtoms() * 3;
  temperature = m_config.temperature;
  kB = m_config.kB;
  vxi1 = vxi2 = xi1 = xi2 = 0.0;
}

Dynamics() [2/2]

eonc::Dynamics::Dynamics ( Matter * matter, const Parameters & parameters )

inline

Definition at line 65 of file Dynamics.h.

      : Dynamics(matter, DynamicsConfig::fromParams(parameters)) {}

~Dynamics()

Dynamics::~Dynamics ( )

default

Member Function Documentation

andersenCollision()

void Dynamics::andersenCollision

(

)

Definition at line 139 of file Dynamics.cpp.

                                 {
  double alpha = m_config.andersen_alpha;
  double tCol = m_config.andersen_tcol;
  double pCol = 1.0 - std::exp(-m_config.time_step / tCol);
 
  AtomMatrix velocity = matter->getVelocities();
  auto mass = matter->getMasses();
 
  for (long i = 0; i < nAtoms; i++) {
    if (randomDouble() < pCol && !matter->getFixed(i)) {
      for (int j = 0; j < 3; j++) {
        double vOld = velocity(i, j);
        double vNew =
            std::sqrt(kB * temperature / mass[i]) * gaussRandom(0.0, 1.0);
        velocity(i, j) = std::sqrt(1.0 - alpha * alpha) * vOld + alpha * vNew;
      }
    }
  }
  matter->setVelocities(velocity);
}

langevinVerlet()

void Dynamics::langevinVerlet

(

)

Langevin dynamics (velocity-Verlet with friction and random forces).

Definition at line 247 of file Dynamics.cpp.

                              {
  double gamma = m_config.langevin_friction;
  AtomMatrix pos = matter->getPositions();
  AtomMatrix vel = matter->getVelocities();
  AtomMatrix acc = matter->getAccelerations();
  AtomMatrix noise = acc; // same shape
  auto mass = matter->getMasses();
 
  // Generate friction + stochastic forces
  AtomMatrix friction = -gamma * vel;
  for (long i = 0; i < nAtoms; i++) {
    if (!matter->getFixed(i)) {
      for (int j = 0; j < 3; j++) {
        noise(i, j) = std::sqrt(4.0 * gamma * kB * temperature / dt / mass[i]) *
                      gaussRandom(0.0, 1.0);
      }
    }
  }
  acc += friction + noise;
 
  vel += acc * 0.5 * dt;
  pos += vel * dt;
  matter->setPositions(pos);
 
  // Second half-step
  acc = matter->getAccelerations();
  friction = -gamma * vel;
  for (long i = 0; i < nAtoms; i++) {
    if (!matter->getFixed(i)) {
      for (int j = 0; j < 3; j++) {
        noise(i, j) = std::sqrt(4.0 * gamma * kB * temperature / dt / mass[i]) *
                      gaussRandom(0.0, 1.0);
      }
    }
  }
  acc += friction + noise;
  vel += 0.5 * dt * acc;
  matter->setVelocities(vel);
}

noseHooverVerlet()

void Dynamics::noseHooverVerlet

(

)

Nose-Hoover chain thermostat (Martyna-Klein-Tuckerman algorithm).

Two chain variables (xi1, xi2) with velocities (vxi1, vxi2).

Definition at line 184 of file Dynamics.cpp.

                                {
  double dt2 = 0.5 * dt;
  double dt4 = 0.25 * dt;
  double dt8 = 0.125 * dt;
  double q1 = m_config.nose_mass;
  double q2 = q1;
  double Temp = kB * temperature;
 
  AtomMatrix vel = matter->getVelocities();
  AtomMatrix pos = matter->getPositions();
  double kinE = matter->getKineticEnergy();
 
  // Forward half-step for chain variables
  double g2 = (q1 * vxi1 * vxi1 - Temp);
  vxi2 += g2 * dt4;
  vxi1 *= std::exp(-vxi2 * dt8);
  double g1 = (2.0 * kinE - nFreeCoords * Temp) / q1;
  vxi1 += g1 * dt4;
  vxi1 *= std::exp(-vxi2 * dt8);
  xi1 += vxi1 * dt2;
  xi2 += vxi2 * dt2;
  double s = std::exp(-vxi1 * dt2);
  vel *= s;
  kinE *= s * s;
  vxi1 *= std::exp(-vxi2 * dt8);
  g1 = (2.0 * kinE - nFreeCoords * Temp) / q1;
  vxi1 += g1 * dt4;
  vxi1 *= std::exp(-vxi2 * dt8);
  g2 = (q1 * vxi1 * vxi1 - Temp) / q2;
  vxi2 += g2 * dt4;
 
  // Position + velocity Verlet step
  pos += vel * dt2;
  matter->setPositions(pos);
  AtomMatrix acc = matter->getAccelerations();
  vel += acc * dt;
  pos += vel * dt2;
  kinE = matter->getKineticEnergy();
 
  // Backward half-step for chain variables
  g2 = (q1 * vxi1 * vxi1 - Temp);
  vxi2 += g2 * dt4;
  vxi1 *= std::exp(-vxi2 * dt8);
  g1 = (2.0 * kinE - nFreeCoords * Temp) / q1;
  vxi1 += g1 * dt4;
  vxi1 *= std::exp(-vxi2 * dt8);
  xi1 += vxi1 * dt2;
  xi2 += vxi2 * dt2;
  s = std::exp(-vxi1 * dt2);
  vel *= s;
  kinE *= s * s;
  vxi1 *= std::exp(-vxi2 * dt8);
  g1 = (2.0 * kinE - nFreeCoords * Temp) / q1;
  vxi1 += g1 * dt4;
  vxi1 *= std::exp(-vxi2 * dt8);
  g2 = (q1 * vxi1 * vxi1 - Temp) / q2;
  vxi2 += g2 * dt4;
 
  matter->setPositions(pos);
  matter->setVelocities(vel);
}

oneStep()

void Dynamics::oneStep

(

int

stepNumber = -1

)

Definition at line 41 of file Dynamics.cpp.

                                     {
  if (m_config.thermostat_kind == ANDERSEN) {
    andersenCollision();
    velocityVerlet();
  } else if (m_config.thermostat_kind == NOSE_HOOVER) {
    noseHooverVerlet();
  } else if (m_config.thermostat_kind == LANGEVIN) {
    langevinVerlet();
  } else if (m_config.thermostat_kind == NONE) {
    velocityVerlet();
  }
 
  if (stepNumber != -1) {
    if (stepNumber == 1) {
      QUILL_LOG_DEBUG(log, "{} {:8s} {:10s} {:12s} {:12s} {:10s}\n",
                      "[Dynamics]", "Step", "KE", "PE", "TE", "KinT");
    }
    double kinE = matter->getKineticEnergy();
    double potE = matter->getPotentialEnergy();
    double kinT = 2.0 * kinE / nFreeCoords / kB;
 
    if (stepNumber % m_config.write_movies_interval == 0) {
      QUILL_LOG_DEBUG(log, "{} {:8} {:10.4} {:12.4} {:12.4} {:10.2}\n",
                      "[Dynamics]", stepNumber, kinE, potE, kinE + potE, kinT);
    }
  }
}

rescaleVelocity()

void Dynamics::rescaleVelocity

(

)

Definition at line 175 of file Dynamics.cpp.

                               {
  AtomMatrix velocity = matter->getVelocities();
  double kinE = matter->getKineticEnergy();
  double kinT = 2.0 * kinE / nFreeCoords / kB;
  matter->setVelocities(velocity * std::sqrt(temperature / kinT));
}

run()

void Dynamics::run

(

)

Definition at line 82 of file Dynamics.cpp.

                   {
  double sumT = 0.0, sumT2 = 0.0;
 
  setThermalVelocity();
 
  if (m_config.thermostat_kind != NONE) {
    QUILL_LOG_DEBUG(log,
                    "{} Running NVT molecular dynamics: {:8.2f} K for {} "
                    "steps ({:.4e} s)\n",
                    "[Dynamics]", temperature, m_config.steps,
                    1e-15 * m_config.time_step * m_config.timeUnit *
                        m_config.steps);
  } else {
    QUILL_LOG_DEBUG(log, "{} Running NVE molecular dynamics: {} steps\n",
                    "[Dynamics]", m_config.steps);
  }
 
  if (m_config.write_movies) {
    if (!eonc::io::io_ok(matter->matter2con("dynamics", false))) {
      QUILL_LOG_WARNING(log, "Failed to write dynamics movie header frame");
    }
  }
 
  QUILL_LOG_DEBUG(log, "{} {:8} {:10} {:12} {:12} {:10}\n", "[Dynamics]",
                  "step", "KE", "PE", "TE", "kinT");
 
  for (long step = 0; step <= m_config.steps; step++) {
    oneStep();
 
    double kinE = matter->getKineticEnergy();
    double potE = matter->getPotentialEnergy();
    double kinT = 2.0 * kinE / nFreeCoords / kB;
    sumT += kinT;
    sumT2 += kinT * kinT;
 
    if (step % m_config.write_movies_interval == 0) {
      QUILL_LOG_DEBUG(log, "{} {} {} {} {} {}\n", "[Dynamics]", step, kinE,
                      potE, kinE + potE, kinT);
    }
 
    if (m_config.write_movies && (step % m_config.write_movies_interval == 0)) {
      if (!eonc::io::io_ok(matter->matter2con("dynamics", true))) {
        QUILL_LOG_WARNING(log, "Failed to append dynamics movie frame");
      }
    }
  }
 
  double avgT = sumT / static_cast<double>(m_config.steps);
  double varT = sumT2 / static_cast<double>(m_config.steps) - avgT * avgT;
  double stdT = std::sqrt(varT);
  QUILL_LOG_DEBUG(log,
                  "{} Temperature : Average = {:.2f} ; StdDev = {:.2f} ; "
                  "Factor = {:.2f}\n",
                  "[Dynamics]", avgT, stdT,
                  varT / avgT / avgT * nFreeCoords / 2.0);
}

setTemperature()

void Dynamics::setTemperature

(

double

temperature

)

Definition at line 37 of file Dynamics.cpp.

                                                   {
  temperature = temperature_in;
}

setThermalVelocity()

void Dynamics::setThermalVelocity

(

)

Definition at line 160 of file Dynamics.cpp.

                                  {
  AtomMatrix velocity = matter->getVelocities();
  auto mass = matter->getMasses();
 
  for (long i = 0; i < nAtoms; i++) {
    if (!matter->getFixed(i)) {
      for (int j = 0; j < 3; j++) {
        velocity(i, j) =
            std::sqrt(kB * temperature / mass[i]) * gaussRandom(0.0, 1.0);
      }
    }
  }
  matter->setVelocities(velocity);
}

velocityVerlet()

void Dynamics::velocityVerlet

(

)

Definition at line 69 of file Dynamics.cpp.

                              {
  AtomMatrix positions = matter->getPositions();
  AtomMatrix velocities = matter->getVelocities();
  AtomMatrix accInit = matter->getAccelerations();
 
  positions += dt * velocities + 0.5 * dt * dt * accInit;
  matter->setPositions(positions);
 
  AtomMatrix accFinal = matter->getAccelerations();
  velocities += 0.5 * dt * (accInit + accFinal);
  matter->setVelocities(velocities);
}

Member Data Documentation

ANDERSEN

const char Dynamics::ANDERSEN = "andersen"

static

Definition at line 57 of file Dynamics.h.

dt

double eonc::Dynamics::dt {0.0}

private

Definition at line 86 of file Dynamics.h.

{0.0};

kB

double eonc::Dynamics::kB {0.0}

private

Definition at line 87 of file Dynamics.h.

{0.0};

LANGEVIN

const char Dynamics::LANGEVIN = "langevin"

static

Definition at line 59 of file Dynamics.h.

log

eonc::log::Scoped eonc::Dynamics::log

private

Definition at line 90 of file Dynamics.h.

m_config

DynamicsConfig eonc::Dynamics::m_config

private

Definition at line 84 of file Dynamics.h.

matter

Matter* eonc::Dynamics::matter

private

Definition at line 83 of file Dynamics.h.

nAtoms

long eonc::Dynamics::nAtoms {0}

private

Definition at line 81 of file Dynamics.h.

{0}, nFreeCoords{0};

nFreeCoords

long eonc::Dynamics::nFreeCoords {0}

private

Definition at line 81 of file Dynamics.h.

{0}, nFreeCoords{0};

NONE

const char Dynamics::NONE = "none"

static

Definition at line 60 of file Dynamics.h.

NOSE_HOOVER

const char Dynamics::NOSE_HOOVER = "nose_hoover"

static

Definition at line 58 of file Dynamics.h.

temperature

double eonc::Dynamics::temperature {0.0}

private

Definition at line 88 of file Dynamics.h.

{0.0};

vxi1

double eonc::Dynamics::vxi1 {0.0}

private

Definition at line 89 of file Dynamics.h.

{0.0}, vxi2{0.0}, xi1{0.0}, xi2{0.0};

vxi2

double eonc::Dynamics::vxi2 {0.0}

private

Definition at line 89 of file Dynamics.h.

{0.0}, vxi2{0.0}, xi1{0.0}, xi2{0.0};

xi1

double eonc::Dynamics::xi1 {0.0}

private

Definition at line 89 of file Dynamics.h.

{0.0}, vxi2{0.0}, xi1{0.0}, xi2{0.0};

xi2

double eonc::Dynamics::xi2 {0.0}

private

Definition at line 89 of file Dynamics.h.

{0.0}, vxi2{0.0}, xi1{0.0}, xi2{0.0};

---

The documentation for this class was generated from the following files:

• /home/runner/work/eOn/eOn/include/eon/Dynamics.h
• /home/runner/work/eOn/eOn/client/Dynamics.cpp