GleThermostat Class Reference

Colored-noise (generalized Langevin) thermostat. More...

#include <GleThermostat.h>

Public Member Functions
	GleThermostat (const MatrixXd &a_drift, double kbt, double dt_half, long n_dof)
	Build from the drift matrix, target kB*T (eV), the half-step dt (internal units), and the number of free degrees of freedom.
void	apply (VectorXd &vel, const VectorXd &masses3N, const std::function< double()> &gauss)
	One OU half-step on the mass-scaled velocities; gauss supplies standard-normal draws from the job's deterministic stream.
bool	valid () const
	True when construction produced a usable propagator.

Static Public Member Functions
static MatrixXd	loadDriftMatrix (const std::string &path)
	Load a drift matrix from a gle4md-layout text file ('#' starts a comment).

Private Attributes
bool	m_valid {false}
MatrixXd	m_T
	exp(-A dt/2)
MatrixXd	m_S
	noise factor, S S^T = kbt (I - T T^T)
MatrixXd	m_Z
	extended state, (1 + n_aux) x n_dof

Detailed Description

Colored-noise (generalized Langevin) thermostat.

Ceriotti, Bussi & Parrinello, "Colored-noise thermostats a la carte", J. Chem. Theory Comput. 6, 1170 (2010): each Cartesian degree of freedom carries n_aux auxiliary momenta, and the extended momentum vector evolves under the exact Ornstein-Uhlenbeck map

(p, s) <- T (p, s) + S xi, T = exp(-A dt/2), S S^T = C - T C T^T, C = kB T * I (canonical sampling),

applied as half-steps around a velocity Verlet integrator. The frequency-dependent friction shaped by A thermalises every mode of a stiff spectrum at comparable efficiency, where white noise critically damps only one band.

The drift matrix A is an operator input ((n_aux+1) x (n_aux+1) plain-text rows, inverse internal time units, the layout produced by the gle4md generator). A 1x1 matrix [gamma] degenerates to standard white-noise Langevin, so the thermostat carries no hidden per-system tuning of its own.

All degrees of freedom propagate together: with the extended state held as a (1 + n_aux) x n_dof matrix Z (row 0 the mass-scaled velocities, rows 1.. the auxiliaries), one half-step is Z <- T Z + S Xi with Xi standard normal.

Definition at line 48 of file GleThermostat.h.

Constructor & Destructor Documentation

GleThermostat()

eonc::GleThermostat::GleThermostat	(	const MatrixXd &	a_drift,
		double	kbt,
		double	dt_half,
		long	n_dof )

Build from the drift matrix, target kB*T (eV), the half-step dt (internal units), and the number of free degrees of freedom.

Definition at line 60 of file GleThermostat.cpp.

                                                         {
  const long dim = a_drift.rows();
  if (dim == 0 || a_drift.cols() != dim || !(kbt > 0.0) || !(dt_half > 0.0) ||
      n_dof <= 0) {
    return;
  }
  m_T = (-a_drift * dt_half).exp();
  // Fluctuation-dissipation: S S^T = kbt (I - T T^T). The right-hand
  // side is symmetric positive semi-definite; factor through its
  // eigendecomposition with negative rounding clamped so a
  // deterministic (dt -> 0) or overdamped block cannot poison the
  // Cholesky.
  const MatrixXd cov =
      kbt * (MatrixXd::Identity(dim, dim) - m_T * m_T.transpose());
  Eigen::SelfAdjointEigenSolver<MatrixXd> es(cov);
  if (es.info() != Eigen::Success) {
    return;
  }
  m_S = es.eigenvectors() *
        es.eigenvalues().cwiseMax(0.0).cwiseSqrt().asDiagonal();
  m_Z = MatrixXd::Zero(dim, n_dof);
  m_valid = true;
}

Member Function Documentation

apply()

void eonc::GleThermostat::apply	(	VectorXd &	vel,
		const VectorXd &	masses3N,
		const std::function< double()> &	gauss )

One OU half-step on the mass-scaled velocities; gauss supplies standard-normal draws from the job's deterministic stream.

Definition at line 85 of file GleThermostat.cpp.

                                                              {
  if (!m_valid || vel.size() != m_Z.cols()) {
    return;
  }
  const VectorXd sqrtM = masses3N.cwiseSqrt();
  m_Z.row(0) = (vel.cwiseProduct(sqrtM)).transpose();
  MatrixXd xi(m_Z.rows(), m_Z.cols());
  for (long i = 0; i < xi.rows(); ++i) {
    for (long j = 0; j < xi.cols(); ++j) {
      xi(i, j) = gauss();
    }
  }
  m_Z = m_T * m_Z + m_S * xi;
  vel = m_Z.row(0).transpose().cwiseQuotient(sqrtM);
}

loadDriftMatrix()

MatrixXd eonc::GleThermostat::loadDriftMatrix ( const std::string & path )

static

Load a drift matrix from a gle4md-layout text file ('#' starts a comment).

Returns an empty (0 x 0) matrix on parse failure.

Definition at line 22 of file GleThermostat.cpp.

                                                             {
  std::ifstream in(path);
  if (!in) {
    return MatrixXd();
  }
  std::vector<std::vector<double>> rows;
  std::string line;
  while (std::getline(in, line)) {
    const auto hash = line.find('#');
    if (hash != std::string::npos) {
      line = line.substr(0, hash);
    }
    std::istringstream ss(line);
    std::vector<double> row;
    double v;
    while (ss >> v) {
      row.push_back(v);
    }
    if (!row.empty()) {
      rows.push_back(std::move(row));
    }
  }
  const auto n = static_cast<long>(rows.size());
  if (n == 0) {
    return MatrixXd();
  }
  MatrixXd a(n, n);
  for (long i = 0; i < n; ++i) {
    if (static_cast<long>(rows[i].size()) != n) {
      return MatrixXd();
    }
    for (long j = 0; j < n; ++j) {
      a(i, j) = rows[i][j];
    }
  }
  return a;
}

valid()

bool eonc::GleThermostat::valid ( ) const

inline

True when construction produced a usable propagator.

Definition at line 65 of file GleThermostat.h.

{ return m_valid; }

Member Data Documentation

m_S

MatrixXd eonc::GleThermostat::m_S

private

noise factor, S S^T = kbt (I - T T^T)

Definition at line 70 of file GleThermostat.h.

m_T

MatrixXd eonc::GleThermostat::m_T

private

exp(-A dt/2)

Definition at line 69 of file GleThermostat.h.

m_valid

bool eonc::GleThermostat::m_valid {false}

private

Definition at line 68 of file GleThermostat.h.

{false};

m_Z

MatrixXd eonc::GleThermostat::m_Z

private

extended state, (1 + n_aux) x n_dof

Definition at line 71 of file GleThermostat.h.

---

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

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