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"
static const char	PILE [] = "pile"
static const char	PIGLET [] = "piglet"

Private Member Functions
void	nhcChainHalfStep (AtomMatrix &vel, double &kinE)
	One Martyna-Klein-Tuckerman chain half-step.
void	runPathIntegral ()
	Ring-polymer trajectory. Beads do not enter velocityVerlet.

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 70 of file Dynamics.h.

Constructor & Destructor Documentation

Dynamics() [1/2]

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

Definition at line 29 of file Dynamics.cpp.

    : matter{matter_in},
      m_config{config} {
  if (!matter) {
    throw std::invalid_argument("Dynamics: null Matter");
  }
  dt = m_config.time_step;
  nAtoms = matter->numberOfAtoms();
  // Unfixed axes only. A partly fixed atom is not three degrees of freedom.
  nFreeCoords = 0;
  for (long i = 0; i < nAtoms; ++i) {
    for (int axis = 0; axis < 3; ++axis) {
      if (!matter->getFixed(i, axis)) {
        ++nFreeCoords;
      }
    }
  }
  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 83 of file Dynamics.h.

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

~Dynamics()

eonc::Dynamics::~Dynamics ( )

default

Member Function Documentation

andersenCollision()

void eonc::Dynamics::andersenCollision

(

)

Definition at line 224 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 (eonc::rng::randomDouble() < pCol && !matter->getFixed(i)) {
      for (int j = 0; j < 3; j++) {
        double vOld = velocity(i, j);
        const double vNew = (mass[i] > 0.0 && kB > 0.0 && temperature > 0.0)
                                ? std::sqrt(kB * temperature / mass[i]) *
                                      eonc::rng::gaussRandom(0.0, 1.0)
                                : 0.0;
        velocity(i, j) = std::sqrt(1.0 - alpha * alpha) * vOld + alpha * vNew;
      }
    }
  }
  matter->setVelocities(velocity);
}

langevinVerlet()

void eonc::Dynamics::langevinVerlet

(

)

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

Matter::getFixed(atom) is true only when every axis is fixed, so noise and the position kick are applied per free axis.

Definition at line 332 of file Dynamics.cpp.

                              {
  const double gamma = m_config.langevin_friction;
  AtomMatrix pos = matter->getPositions();
  AtomMatrix vel = matter->getVelocities();
  AtomMatrix acc = matter->getAccelerations();
  AtomMatrix noise = AtomMatrix::Zero(nAtoms, 3);
  const auto mass = matter->getMasses();
 
  // Zero a frozen axis by assignment. Multiplying by the free mask leaves
  // a NaN in place, and setPositions stores every component.
  auto holdFixed = [&](AtomMatrix &m) {
    for (long i = 0; i < nAtoms; i++) {
      for (int j = 0; j < 3; j++) {
        if (matter->getFixed(i, j)) {
          m(i, j) = 0.0;
        }
      }
    }
  };
  auto addNoise = [&]() {
    noise.setZero();
    for (long i = 0; i < nAtoms; i++) {
      for (int j = 0; j < 3; j++) {
        if (matter->getFixed(i, j)) {
          continue;
        }
        noise(i, j) = std::sqrt(4.0 * gamma * kB * temperature / dt / mass[i]) *
                      eonc::rng::gaussRandom(0.0, 1.0);
      }
    }
  };
 
  addNoise();
  AtomMatrix friction = -gamma * vel;
  holdFixed(friction);
  holdFixed(acc);
  acc += friction + noise;
 
  vel += acc * 0.5 * dt;
  holdFixed(vel);
  pos += vel * dt;
  matter->setPositions(pos);
 
  acc = matter->getAccelerations();
  addNoise();
  friction = -gamma * vel;
  holdFixed(friction);
  holdFixed(acc);
  acc += friction + noise;
  vel += 0.5 * dt * acc;
  holdFixed(vel);
  matter->setVelocities(vel);
}

nhcChainHalfStep()

void eonc::Dynamics::nhcChainHalfStep ( AtomMatrix & vel, double & kinE )

private

One Martyna-Klein-Tuckerman chain half-step.

G2 is always (Q1 vxi1^2 - kT) / Q2.

Definition at line 275 of file Dynamics.cpp.

                                                             {
  const double dt2 = 0.5 * dt;
  const double dt4 = 0.25 * dt;
  const double dt8 = 0.125 * dt;
  const double q1 = m_config.nose_mass;
  const double q2 = q1;
  const double Temp = kB * temperature;
  if (!(q1 > 0.0)) {
    throw std::invalid_argument("thermostat.nose_mass must be positive");
  }
 
  // Martyna, Klein, Tuckerman JCP 97, 2635 (1992): G2 = (Q1 v_ξ1² − kT) / Q2.
  auto g2 = [&]() { return (q1 * vxi1 * vxi1 - Temp) / q2; };
  auto g1 = [&]() { return (2.0 * kinE - nFreeCoords * Temp) / q1; };
 
  vxi2 += g2() * dt4;
  vxi1 *= std::exp(-vxi2 * dt8);
  vxi1 += g1() * dt4;
  vxi1 *= std::exp(-vxi2 * dt8);
  xi1 += vxi1 * dt2;
  xi2 += vxi2 * dt2;
  const double s = std::exp(-vxi1 * dt2);
  vel *= s;
  kinE *= s * s;
  vxi1 *= std::exp(-vxi2 * dt8);
  vxi1 += g1() * dt4;
  vxi1 *= std::exp(-vxi2 * dt8);
  vxi2 += g2() * dt4;
}

noseHooverVerlet()

void eonc::Dynamics::noseHooverVerlet

(

)

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

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

Definition at line 307 of file Dynamics.cpp.

                                {
  const double dt2 = 0.5 * dt;
  AtomMatrix vel = matter->getVelocities();
  AtomMatrix pos = matter->getPositions();
  double kinE = matter->getKineticEnergy();
 
  nhcChainHalfStep(vel, kinE);
 
  pos += vel * dt2;
  matter->setPositions(pos);
  AtomMatrix acc = matter->getAccelerations();
  vel += acc * dt;
  pos += vel * dt2;
  matter->setPositions(pos);
  matter->setVelocities(vel);
  kinE = matter->getKineticEnergy();
 
  nhcChainHalfStep(vel, kinE);
 
  matter->setVelocities(vel);
}

oneStep()

void eonc::Dynamics::oneStep

(

int

stepNumber = -1

)

Definition at line 57 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();
  } else if (m_config.thermostat_kind == PILE ||
             m_config.thermostat_kind == PIGLET) {
    throw std::invalid_argument(
        "path-integral dynamics is a ring-polymer trajectory, not a "
        "velocity Verlet step");
  }
 
  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();
    const double kinT =
        (nFreeCoords > 0 && kB > 0.0) ? 2.0 * kinE / nFreeCoords / kB : 0.0;
 
    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 eonc::Dynamics::rescaleVelocity

(

)

Definition at line 264 of file Dynamics.cpp.

                               {
  AtomMatrix velocity = matter->getVelocities();
  double kinE = matter->getKineticEnergy();
  const double kinT =
      (nFreeCoords > 0 && kB > 0.0) ? 2.0 * kinE / nFreeCoords / kB : 0.0;
  if (!(kinT > 0.0) || !(temperature > 0.0)) {
    return;
  }
  matter->setVelocities(velocity * std::sqrt(temperature / kinT));
}

run()

void eonc::Dynamics::run

(

)

Definition at line 161 of file Dynamics.cpp.

                   {
  if (m_config.thermostat_kind == PILE || m_config.thermostat_kind == PIGLET) {
    runPathIntegral();
    return;
  }
  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();
    const double kinT =
        (nFreeCoords > 0 && kB > 0.0) ? 2.0 * kinE / nFreeCoords / kB : 0.0;
    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");
      }
    }
  }
 
  const double nstat = static_cast<double>(std::max(m_config.steps, 1L));
  double avgT = sumT / nstat;
  double varT = sumT2 / nstat - 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);
}

runPathIntegral()

void eonc::Dynamics::runPathIntegral ( )

private

Ring-polymer trajectory. Beads do not enter velocityVerlet.

Definition at line 104 of file Dynamics.cpp.

                               {
  auto pot = matter->getPotential();
  if (!pot) {
    throw std::invalid_argument("path-integral dynamics has no potential");
  }
  const long atoms = matter->numberOfAtoms();
  std::vector<double> masses(static_cast<size_t>(atoms));
  std::vector<int> numbers(static_cast<size_t>(atoms));
  std::vector<char> free(static_cast<size_t>(3 * atoms), 1);
  const VectorXi z = matter->getAtomicNrs();
  for (long i = 0; i < atoms; ++i) {
    masses[static_cast<size_t>(i)] = matter->getMass(i);
    numbers[static_cast<size_t>(i)] = z[i];
    for (int axis = 0; axis < 3; ++axis) {
      if (matter->getFixed(i, axis)) {
        free[static_cast<size_t>(3 * i + axis)] = 0;
      }
    }
  }
  pathintegral::Options opt;
  opt.beads = m_config.path_beads;
  opt.temperature = temperature;
  opt.kB = kB;
  opt.hbar = tunneling::kHbar;
  opt.dt = dt;
  opt.springs = m_config.path_springs == "eco" ? pathintegral::Springs::Eco
                                               : pathintegral::Springs::Trotter;
  opt.thermostat = m_config.thermostat_kind == PIGLET
                       ? pathintegral::Thermostat::Piglet
                       : pathintegral::Thermostat::Pile;
  opt.pileTau = m_config.path_pile_tau;
  opt.pileScale = m_config.path_pile_scale;
  opt.ecoOmegaMax = m_config.path_eco_omega_max;
  opt.gleFile = m_config.path_gle_file;
  opt.seed = m_config.path_seed;
  pathintegral::RingPolymer ring(atoms, masses, numbers, free, opt);
  const AtomMatrix pos = matter->getPositions();
  ring.setAllBeads(pos.data());
  const Matrix3d box =
      matter->getPeriodic() ? matter->getCell() : Matrix3d::Zero().eval();
  for (long step = 0; step < m_config.steps; ++step) {
    ring.step(*pot, box.data(), false);
  }
  const VectorXd centroid = ring.centroid();
  const VectorXd velocity = ring.centroidVelocity();
  AtomMatrix out = pos;
  AtomMatrix vel = matter->getVelocities();
  for (long i = 0; i < atoms; ++i) {
    for (int axis = 0; axis < 3; ++axis) {
      out(i, axis) = centroid[3 * i + axis];
      vel(i, axis) = velocity[3 * i + axis];
    }
  }
  matter->setPositions(out);
  matter->setVelocities(vel);
}

setTemperature()

void eonc::Dynamics::setTemperature

(

double

temperature

)

Definition at line 53 of file Dynamics.cpp.

                                                   {
  temperature = temperature_in;
}

setThermalVelocity()

void eonc::Dynamics::setThermalVelocity

(

)

Definition at line 247 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) = (mass[i] > 0.0 && kB > 0.0 && temperature > 0.0)
                             ? std::sqrt(kB * temperature / mass[i]) *
                                   eonc::rng::gaussRandom(0.0, 1.0)
                             : 0.0;
      }
    }
  }
  matter->setVelocities(velocity);
}

velocityVerlet()

void eonc::Dynamics::velocityVerlet

(

)

Definition at line 91 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 eonc::Dynamics::ANDERSEN = "andersen"

static

Definition at line 73 of file Dynamics.h.

dt

double eonc::Dynamics::dt {0.0}

private

Definition at line 110 of file Dynamics.h.

{0.0};

kB

double eonc::Dynamics::kB {0.0}

private

Definition at line 111 of file Dynamics.h.

{0.0};

LANGEVIN

const char eonc::Dynamics::LANGEVIN = "langevin"

static

Definition at line 75 of file Dynamics.h.

log

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

private

Definition at line 114 of file Dynamics.h.

m_config

DynamicsConfig eonc::Dynamics::m_config

private

Definition at line 108 of file Dynamics.h.

matter

Matter* eonc::Dynamics::matter

private

Definition at line 107 of file Dynamics.h.

nAtoms

long eonc::Dynamics::nAtoms {0}

private

Definition at line 105 of file Dynamics.h.

{0}, nFreeCoords{0};

nFreeCoords

long eonc::Dynamics::nFreeCoords {0}

private

Definition at line 105 of file Dynamics.h.

{0}, nFreeCoords{0};

NONE

const char eonc::Dynamics::NONE = "none"

static

Definition at line 76 of file Dynamics.h.

NOSE_HOOVER

const char eonc::Dynamics::NOSE_HOOVER = "nose_hoover"

static

Definition at line 74 of file Dynamics.h.

PIGLET

const char eonc::Dynamics::PIGLET = "piglet"

static

Definition at line 78 of file Dynamics.h.

PILE

const char eonc::Dynamics::PILE = "pile"

static

Definition at line 77 of file Dynamics.h.

temperature

double eonc::Dynamics::temperature {0.0}

private

Definition at line 112 of file Dynamics.h.

{0.0};

vxi1

double eonc::Dynamics::vxi1 {0.0}

private

Definition at line 113 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 113 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 113 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 113 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