eonc::pathintegral::RingPolymer Class Reference

Ring-polymer NVT step. More...

#include <PathIntegral.h>

Classes
struct	ModeGle

Public Member Functions
	RingPolymer (long nAtoms, std::vector< double > masses, std::vector< int > atomicNumbers, std::vector< char > free, Options opt)
void	setAllBeads (const double *q)
void	setBeads (const std::vector< VectorXd > &beads)
	One position vector of length 3 * nAtoms per bead.
const std::vector< VectorXd > &	beads () const
const std::vector< VectorXd > &	momenta () const
void	setMomenta (const std::vector< VectorXd > &momenta)
	One momentum vector of length 3 * nAtoms per bead; fixed coordinates are zeroed.
void	thermalMomenta ()
void	setHyperplane (const VectorXd &normal, const VectorXd &origin)
	Hold n · (q_centroid - origin) = 0.
void	step (Potential &pot, const double *box, bool record)
void	nveStep (Potential &pot, const double *box)
	Thermostat-free RPMD step: velocity Verlet with the free ring propagated exactly in normal modes.
double	kineticCv () const
double	meanForce () const
long	batches () const
VectorXd	centroid () const
VectorXd	centroidVelocity () const
void	resetAverages ()
Sample	sample (Potential &pot, const double *box, long equilibration, long production)

Private Member Functions
void	forces (Potential &pot, const double *box)
void	thermostat (double h)
void	kick (double h, bool dropParallel)
void	propagate (double h)
void	projectPosition ()
void	projectMomentum ()
void	toNormal (const std::vector< VectorXd > &src, std::vector< VectorXd > &dst) const
void	fromNormal (const std::vector< VectorXd > &src, std::vector< VectorXd > &dst) const
double	gauss ()
void	initGle ()

Private Attributes
Options	opt_
long	nAtoms_ {0}
long	nDof_ {0}
long	nBeads_ {0}
long	nFree_ {0}
std::vector< double >	mass_
std::vector< int >	atomicNumbers_
std::vector< char >	free_
std::vector< long >	freeIndex_
MatrixXd	modes_
VectorXd	omegaK_
std::vector< VectorXd >	q_
std::vector< VectorXd >	p_
std::vector< VectorXd >	f_
std::vector< VectorXd >	qnm_
std::vector< VectorXd >	pnm_
VectorXd	planeNormal_
VectorXd	planeOrigin_
bool	constrain_ {false}
bool	haveForces_ {false}
std::vector< ModeGle >	gle_
std::uint64_t	rng_
long	batches_ {0}
long	recorded_ {0}
double	kineticSum_ {0.0}
double	forceSum_ {0.0}

Detailed Description

Ring-polymer NVT step.

Physical forces on every bead are one forceBatch call, so a calculator group carries the beads together. The classical velocity Verlet step is not used.

Definition at line 92 of file PathIntegral.h.

Constructor & Destructor Documentation

RingPolymer()

eonc::pathintegral::RingPolymer::RingPolymer	(	long	nAtoms,
		std::vector< double >	masses,
		std::vector< int >	atomicNumbers,
		std::vector< char >	free,
		Options	opt )

Definition at line 320 of file PathIntegral.cpp.

    : opt_(std::move(opt)),
      nAtoms_(nAtoms),
      nDof_(3 * nAtoms),
      nBeads_(opt_.beads),
      atomicNumbers_(std::move(atomicNumbers)),
      free_(std::move(free)),
      rng_(opt_.seed == 0 ? 1 : opt_.seed) {
  if (nAtoms_ < 1 || nBeads_ < 1) {
    throw std::invalid_argument("path integral needs atoms and beads");
  }
  if (static_cast<long>(masses.size()) != nAtoms_ ||
      static_cast<long>(atomicNumbers_.size()) != nAtoms_ ||
      static_cast<long>(free_.size()) != nDof_) {
    throw std::invalid_argument("path integral mass, number or mask size");
  }
  if (!(opt_.temperature > 0.0) || !(opt_.kB > 0.0) || !(opt_.hbar > 0.0) ||
      !(opt_.dt > 0.0) || !(opt_.pileTau > 0.0) || !(opt_.pileScale > 0.0)) {
    throw std::invalid_argument(
        "path integral temperature, timestep and damping must be positive");
  }
  if (opt_.springs == Springs::Eco && opt_.thermostat == Thermostat::Piglet) {
    throw std::invalid_argument(
        "economised springs cannot be combined with a normal-mode GLE");
  }
  mass_.assign(static_cast<size_t>(nDof_), 0.0);
  nFree_ = 0;
  for (long i = 0; i < nAtoms_; ++i) {
    if (!(masses[static_cast<size_t>(i)] > 0.0)) {
      throw std::invalid_argument("path integral atom has no mass");
    }
    for (int axis = 0; axis < 3; ++axis) {
      const long a = 3 * i + axis;
      mass_[static_cast<size_t>(a)] = masses[static_cast<size_t>(i)];
      if (free_[static_cast<size_t>(a)]) {
        freeIndex_.push_back(a);
        ++nFree_;
      }
    }
  }
  if (nFree_ < 1) {
    throw std::invalid_argument("path integral has no free coordinate");
  }
  modes_ = normalModeMatrix(nBeads_);
  const double omegan =
      static_cast<double>(nBeads_) * opt_.kB * opt_.temperature / opt_.hbar;
  if (opt_.springs == Springs::Eco) {
    const double xmax =
        opt_.ecoOmegaMax * opt_.hbar / (opt_.kB * opt_.temperature);
    omegaK_ = omegan * ecoEigenvalues(nBeads_, xmax);
  } else {
    omegaK_ = omegan * trotterEigenvalues(nBeads_);
  }
  q_.assign(static_cast<size_t>(nBeads_), VectorXd::Zero(nDof_));
  p_.assign(static_cast<size_t>(nBeads_), VectorXd::Zero(nDof_));
  f_.assign(static_cast<size_t>(nBeads_), VectorXd::Zero(nDof_));
  qnm_.assign(static_cast<size_t>(nBeads_), VectorXd::Zero(nDof_));
  pnm_.assign(static_cast<size_t>(nBeads_), VectorXd::Zero(nDof_));
  if (opt_.thermostat == Thermostat::Piglet) {
    initGle();
  }
  thermalMomenta();
}

Member Function Documentation

batches()

long eonc::pathintegral::RingPolymer::batches ( ) const

inlinenodiscard

Definition at line 120 of file PathIntegral.h.

{ return batches_; }

beads()

const std::vector< VectorXd > & eonc::pathintegral::RingPolymer::beads ( ) const

inlinenodiscard

Definition at line 102 of file PathIntegral.h.

{ return q_; }

centroid()

VectorXd eonc::pathintegral::RingPolymer::centroid ( ) const

nodiscard

Definition at line 585 of file PathIntegral.cpp.

                                     {
  VectorXd c = VectorXd::Zero(nDof_);
  for (long bead = 0; bead < nBeads_; ++bead) {
    c += q_[static_cast<size_t>(bead)];
  }
  c /= static_cast<double>(nBeads_);
  return c;
}

centroidVelocity()

VectorXd eonc::pathintegral::RingPolymer::centroidVelocity ( ) const

nodiscard

Definition at line 594 of file PathIntegral.cpp.

                                             {
  std::vector<VectorXd> pnm(static_cast<size_t>(nBeads_),
                            VectorXd::Zero(nDof_));
  toNormal(p_, pnm);
  VectorXd v = VectorXd::Zero(nDof_);
  const double scale = std::sqrt(static_cast<double>(nBeads_));
  for (long a = 0; a < nDof_; ++a) {
    const double m = mass_[static_cast<size_t>(a)];
    if (m > 0.0) {
      v[a] = pnm[0][a] / (scale * m);
    }
  }
  return v;
}

forces()

void eonc::pathintegral::RingPolymer::forces ( Potential & pot, const double * box )

private

Definition at line 558 of file PathIntegral.cpp.

                                                          {
  double zeroBox[9] = {};
  const double *cell = box != nullptr ? box : zeroBox;
  std::vector<const double *> pos(static_cast<size_t>(nBeads_));
  std::vector<const int *> nrs(static_cast<size_t>(nBeads_));
  std::vector<double *> frc(static_cast<size_t>(nBeads_));
  std::vector<const double *> boxes(static_cast<size_t>(nBeads_), cell);
  for (long bead = 0; bead < nBeads_; ++bead) {
    pos[static_cast<size_t>(bead)] = q_[static_cast<size_t>(bead)].data();
    nrs[static_cast<size_t>(bead)] = atomicNumbers_.data();
    frc[static_cast<size_t>(bead)] = f_[static_cast<size_t>(bead)].data();
  }
  std::vector<double> energies(static_cast<size_t>(nBeads_), 0.0);
  std::vector<double> variances(static_cast<size_t>(nBeads_), 0.0);
  pot.forceBatch(nBeads_, nAtoms_, pos.data(), nrs.data(), frc.data(),
                 energies.data(), variances.data(), boxes.data());
  ++batches_;
  for (long bead = 0; bead < nBeads_; ++bead) {
    for (long a = 0; a < nDof_; ++a) {
      if (!free_[static_cast<size_t>(a)]) {
        f_[static_cast<size_t>(bead)][a] = 0.0;
      }
    }
  }
  haveForces_ = true;
}

fromNormal()

void eonc::pathintegral::RingPolymer::fromNormal ( const std::vector< VectorXd > & src, std::vector< VectorXd > & dst ) const

private

Definition at line 546 of file PathIntegral.cpp.

                                                             {
  MatrixXd packed(nDof_, nBeads_);
  for (long k = 0; k < nBeads_; ++k) {
    packed.col(k) = src[static_cast<size_t>(k)];
  }
  const MatrixXd out = packed * modes_;
  for (long j = 0; j < nBeads_; ++j) {
    dst[static_cast<size_t>(j)] = out.col(j);
  }
}

gauss()

double eonc::pathintegral::RingPolymer::gauss ( )

private

Definition at line 487 of file PathIntegral.cpp.

                          {
  rng_ = rng_ * 6364136223846793005ULL + 1ULL;
  const double u1 =
      std::max((rng_ >> 11) * (1.0 / 9007199254740992.0), 1.0e-16);
  rng_ = rng_ * 6364136223846793005ULL + 1ULL;
  const double u2 = (rng_ >> 11) * (1.0 / 9007199254740992.0);
  return std::sqrt(-2.0 * std::log(u1)) * std::cos(2.0 * kPi * u2);
}

initGle()

void eonc::pathintegral::RingPolymer::initGle ( )

private

Definition at line 386 of file PathIntegral.cpp.

                          {
  if (nBeads_ == 1) {
    return;
  }
  if (opt_.gleFile.empty()) {
    throw std::invalid_argument("normal-mode GLE needs a matrix file");
  }
  const GleFile file = readGle(opt_.gleFile);
  long first = 0;
  if (file.nModes == nBeads_) {
    first = 1;
  } else if (file.nModes != nBeads_ - 1) {
    throw std::invalid_argument(
        "GLE matrix count must be the bead count or one less");
  }
  const double h = 0.5 * opt_.dt;
  gle_.resize(static_cast<size_t>(nBeads_ - 1));
  for (long k = 1; k < nBeads_; ++k) {
    const long src = first + (k - 1);
    ModeGle mode;
    mode.drift = file.drift[static_cast<size_t>(src)];
    mode.covariance = file.covariance[static_cast<size_t>(src)];
    mode.propagate = (-mode.drift * h).exp();
    const MatrixXd cov =
        opt_.kB * (mode.covariance - mode.propagate * mode.covariance *
                                         mode.propagate.transpose());
    mode.noise = factorCovariance(cov);
    mode.extended = MatrixXd::Zero(file.dim, nFree_);
    gle_[static_cast<size_t>(k - 1)] = std::move(mode);
  }
}

kick()

void eonc::pathintegral::RingPolymer::kick ( double h, bool dropParallel )

private

Definition at line 681 of file PathIntegral.cpp.

                                                  {
  VectorXd removal = VectorXd::Zero(nDof_);
  if (dropParallel && constrain_) {
    VectorXd fc = VectorXd::Zero(nDof_);
    for (long bead = 0; bead < nBeads_; ++bead) {
      fc += f_[static_cast<size_t>(bead)];
    }
    fc /= static_cast<double>(nBeads_);
    double along = 0.0;
    for (long a : freeIndex_) {
      along += planeNormal_[a] * fc[a];
    }
    removal = along * planeNormal_;
  }
  for (long bead = 0; bead < nBeads_; ++bead) {
    for (long a : freeIndex_) {
      p_[static_cast<size_t>(bead)][a] +=
          (f_[static_cast<size_t>(bead)][a] - removal[a]) * h;
    }
  }
}

kineticCv()

double eonc::pathintegral::RingPolymer::kineticCv ( ) const

nodiscard

Definition at line 609 of file PathIntegral.cpp.

                                    {
  double k = 0.5 * static_cast<double>(nFree_) * opt_.kB * opt_.temperature;
  if (!haveForces_) {
    return k;
  }
  const VectorXd c = centroid();
  double virial = 0.0;
  for (long bead = 0; bead < nBeads_; ++bead) {
    for (long a : freeIndex_) {
      virial += (q_[static_cast<size_t>(bead)][a] - c[a]) *
                f_[static_cast<size_t>(bead)][a];
    }
  }
  k += -0.5 / static_cast<double>(nBeads_) * virial;
  return k;
}

meanForce()

double eonc::pathintegral::RingPolymer::meanForce ( ) const

nodiscard

Definition at line 626 of file PathIntegral.cpp.

                                    {
  if (!constrain_ || recorded_ < 1) {
    return 0.0;
  }
  return forceSum_ / static_cast<double>(recorded_);
}

momenta()

const std::vector< VectorXd > & eonc::pathintegral::RingPolymer::momenta ( ) const

inlinenodiscard

Definition at line 103 of file PathIntegral.h.

{ return p_; }

nveStep()

void eonc::pathintegral::RingPolymer::nveStep	(	Potential &	pot,
		const double *	box )

Thermostat-free RPMD step: velocity Verlet with the free ring propagated exactly in normal modes.

Refused while a hyperplane is set.

Definition at line 801 of file PathIntegral.cpp.

                                                           {
  if (constrain_) {
    throw std::logic_error("path integral NVE step with a hyperplane set");
  }
  const double half = 0.5 * opt_.dt;
  if (!haveForces_) {
    forces(pot, box);
  }
  kick(half, false);
  propagate(opt_.dt);
  forces(pot, box);
  kick(half, false);
}

projectMomentum()

void eonc::pathintegral::RingPolymer::projectMomentum ( )

private

Definition at line 741 of file PathIntegral.cpp.

                                  {
  if (!constrain_) {
    return;
  }
  VectorXd sum = VectorXd::Zero(nDof_);
  for (long bead = 0; bead < nBeads_; ++bead) {
    sum += p_[static_cast<size_t>(bead)];
  }
  // p_0 = sum / sqrt(P); remove lambda n from p_0.
  const double scale = std::sqrt(static_cast<double>(nBeads_));
  double num = 0.0;
  double den = 0.0;
  for (long a : freeIndex_) {
    const double m = mass_[static_cast<size_t>(a)];
    num += planeNormal_[a] * (sum[a] / scale) / m;
    den += planeNormal_[a] * planeNormal_[a] / m;
  }
  if (den > 0.0) {
    const double shift = num / den / scale;
    for (long bead = 0; bead < nBeads_; ++bead) {
      for (long a : freeIndex_) {
        p_[static_cast<size_t>(bead)][a] -= shift * planeNormal_[a];
      }
    }
  }
}

projectPosition()

void eonc::pathintegral::RingPolymer::projectPosition ( )

private

Definition at line 725 of file PathIntegral.cpp.

                                  {
  if (!constrain_) {
    return;
  }
  const VectorXd c = centroid();
  double sigma = 0.0;
  for (long a : freeIndex_) {
    sigma += planeNormal_[a] * (c[a] - planeOrigin_[a]);
  }
  for (long bead = 0; bead < nBeads_; ++bead) {
    for (long a : freeIndex_) {
      q_[static_cast<size_t>(bead)][a] -= sigma * planeNormal_[a];
    }
  }
}

propagate()

void eonc::pathintegral::RingPolymer::propagate ( double h )

private

Definition at line 703 of file PathIntegral.cpp.

                                    {
  toNormal(q_, qnm_);
  toNormal(p_, pnm_);
  for (long a : freeIndex_) {
    const double m = mass_[static_cast<size_t>(a)];
    qnm_[0][a] += pnm_[0][a] / m * h;
    for (long k = 1; k < nBeads_; ++k) {
      const double omega = omegaK_[k];
      const double c = std::cos(omega * h);
      const double s = std::sin(omega * h);
      const double pk = pnm_[static_cast<size_t>(k)][a];
      const double qk = qnm_[static_cast<size_t>(k)][a];
      pnm_[static_cast<size_t>(k)][a] = c * pk - m * omega * s * qk;
      qnm_[static_cast<size_t>(k)][a] = c * qk + s * pk / (omega * m);
    }
  }
  fromNormal(qnm_, q_);
  fromNormal(pnm_, p_);
}

resetAverages()

void eonc::pathintegral::RingPolymer::resetAverages

(

)

Definition at line 633 of file PathIntegral.cpp.

                                {
  recorded_ = 0;
  kineticSum_ = 0.0;
  forceSum_ = 0.0;
}

sample()

Sample eonc::pathintegral::RingPolymer::sample	(	Potential &	pot,
		const double *	box,
		long	equilibration,
		long	production )

Definition at line 815 of file PathIntegral.cpp.

                                                                {
  if (equilibration < 0 || production < 1) {
    throw std::invalid_argument("path integral sample length");
  }
  for (long step = 0; step < equilibration; ++step) {
    this->step(pot, box, false);
  }
  resetAverages();
  const long batchesBefore = batches_;
  for (long step = 0; step < production; ++step) {
    this->step(pot, box, true);
  }
  Sample out;
  out.kineticCv = kineticSum_ / static_cast<double>(recorded_);
  out.meanForce = meanForce();
  out.batches = batches_ - batchesBefore;
  return out;
}

setAllBeads()

void eonc::pathintegral::RingPolymer::setAllBeads

(

const double *

q

)

Definition at line 418 of file PathIntegral.cpp.

                                             {
  if (q == nullptr) {
    throw std::invalid_argument("path integral positions are missing");
  }
  for (long bead = 0; bead < nBeads_; ++bead) {
    for (long a = 0; a < nDof_; ++a) {
      q_[static_cast<size_t>(bead)][a] = q[a];
    }
  }
  if (constrain_) {
    projectPosition();
  }
  haveForces_ = false;
}

setBeads()

void eonc::pathintegral::RingPolymer::setBeads

(

const std::vector< VectorXd > &

beads

)

One position vector of length 3 * nAtoms per bead.

The centroid is projected onto the hyperplane when one is set.

Definition at line 433 of file PathIntegral.cpp.

                                                           {
  if (static_cast<long>(beads.size()) != nBeads_) {
    throw std::invalid_argument("path integral bead count mismatch");
  }
  for (long bead = 0; bead < nBeads_; ++bead) {
    if (beads[static_cast<size_t>(bead)].size() != nDof_) {
      throw std::invalid_argument("path integral bead has the wrong length");
    }
    q_[static_cast<size_t>(bead)] = beads[static_cast<size_t>(bead)];
  }
  if (constrain_) {
    projectPosition();
  }
  haveForces_ = false;
}

setHyperplane()

void eonc::pathintegral::RingPolymer::setHyperplane	(	const VectorXd &	normal,
		const VectorXd &	origin )

Hold n · (q_centroid - origin) = 0.

n is normalised on the free coordinates. origin has length 3 * nAtoms.

Definition at line 466 of file PathIntegral.cpp.

                                                        {
  if (normal.size() != nDof_ || origin.size() != nDof_) {
    throw std::invalid_argument("hyperplane vectors have the wrong length");
  }
  planeNormal_ = VectorXd::Zero(nDof_);
  planeOrigin_ = origin;
  double norm = 0.0;
  for (long a : freeIndex_) {
    planeNormal_[a] = normal[a];
    norm += normal[a] * normal[a];
  }
  if (!(norm > 0.0)) {
    throw std::invalid_argument("hyperplane normal has no free component");
  }
  planeNormal_ /= std::sqrt(norm);
  constrain_ = true;
  projectPosition();
  projectMomentum();
}

setMomenta()

void eonc::pathintegral::RingPolymer::setMomenta

(

const std::vector< VectorXd > &

momenta

)

One momentum vector of length 3 * nAtoms per bead; fixed coordinates are zeroed.

Definition at line 449 of file PathIntegral.cpp.

                                                               {
  if (static_cast<long>(momenta.size()) != nBeads_) {
    throw std::invalid_argument("path integral bead count mismatch");
  }
  for (long bead = 0; bead < nBeads_; ++bead) {
    if (momenta[static_cast<size_t>(bead)].size() != nDof_) {
      throw std::invalid_argument(
          "path integral momentum has the wrong length");
    }
    for (long a = 0; a < nDof_; ++a) {
      p_[static_cast<size_t>(bead)][a] =
          free_[static_cast<size_t>(a)] ? momenta[static_cast<size_t>(bead)][a]
                                        : 0.0;
    }
  }
}

step()

void eonc::pathintegral::RingPolymer::step	(	Potential &	pot,
		const double *	box,
		bool	record )

Definition at line 768 of file PathIntegral.cpp.

                                                                     {
  const double half = 0.5 * opt_.dt;
  thermostat(half);
  projectMomentum();
  forces(pot, box);
  kick(half, true);
  projectMomentum();
  propagate(half);
  propagate(half);
  projectPosition();
  forces(pot, box);
  if (record) {
    kineticSum_ += kineticCv();
    if (constrain_) {
      VectorXd fc = VectorXd::Zero(nDof_);
      for (long bead = 0; bead < nBeads_; ++bead) {
        fc += f_[static_cast<size_t>(bead)];
      }
      fc /= static_cast<double>(nBeads_);
      double along = 0.0;
      for (long a : freeIndex_) {
        along += planeNormal_[a] * fc[a];
      }
      forceSum_ += along;
    }
    ++recorded_;
  }
  kick(half, true);
  projectMomentum();
  thermostat(half);
  projectMomentum();
}

thermalMomenta()

void eonc::pathintegral::RingPolymer::thermalMomenta

(

)

Definition at line 496 of file PathIntegral.cpp.

                                 {
  for (long bead = 0; bead < nBeads_; ++bead) {
    p_[static_cast<size_t>(bead)].setZero();
  }
  toNormal(p_, pnm_);
  const double tSim = static_cast<double>(nBeads_) * opt_.kB * opt_.temperature;
  for (long k = 0; k < nBeads_; ++k) {
    const bool gleMode = opt_.thermostat == Thermostat::Piglet && k > 0;
    if (gleMode) {
      continue;
    }
    for (long a : freeIndex_) {
      pnm_[static_cast<size_t>(k)][a] =
          gauss() * std::sqrt(mass_[static_cast<size_t>(a)] * tSim);
    }
  }
  if (opt_.thermostat == Thermostat::Piglet) {
    for (long k = 1; k < nBeads_; ++k) {
      ModeGle &mode = gle_[static_cast<size_t>(k - 1)];
      const MatrixXd factor = factorCovariance(opt_.kB * mode.covariance);
      MatrixXd noise(mode.extended.rows(), mode.extended.cols());
      for (long r = 0; r < noise.rows(); ++r) {
        for (long c = 0; c < noise.cols(); ++c) {
          noise(r, c) = gauss();
        }
      }
      mode.extended = factor * noise;
      for (long col = 0; col < nFree_; ++col) {
        const long a = freeIndex_[static_cast<size_t>(col)];
        pnm_[static_cast<size_t>(k)][a] =
            mode.extended(0, col) * std::sqrt(mass_[static_cast<size_t>(a)]);
      }
    }
  }
  fromNormal(pnm_, p_);
  haveForces_ = false;
}

thermostat()

void eonc::pathintegral::RingPolymer::thermostat ( double h )

private

Definition at line 639 of file PathIntegral.cpp.

                                     {
  toNormal(p_, pnm_);
  const double tSim = static_cast<double>(nBeads_) * opt_.kB * opt_.temperature;
  auto langevin = [&](long k, double tau) {
    const double damp = std::exp(-h / tau);
    const double noise = std::sqrt(tSim * (1.0 - damp * damp));
    for (long a : freeIndex_) {
      const double sm = std::sqrt(mass_[static_cast<size_t>(a)]);
      double pms = pnm_[static_cast<size_t>(k)][a] / sm;
      pms = damp * pms + noise * gauss();
      pnm_[static_cast<size_t>(k)][a] = pms * sm;
    }
  };
  langevin(0, opt_.pileTau);
  for (long k = 1; k < nBeads_; ++k) {
    if (opt_.thermostat == Thermostat::Piglet) {
      ModeGle &mode = gle_[static_cast<size_t>(k - 1)];
      for (long col = 0; col < nFree_; ++col) {
        const long a = freeIndex_[static_cast<size_t>(col)];
        mode.extended(0, col) = pnm_[static_cast<size_t>(k)][a] /
                                std::sqrt(mass_[static_cast<size_t>(a)]);
      }
      MatrixXd noise(mode.extended.rows(), mode.extended.cols());
      for (long r = 0; r < noise.rows(); ++r) {
        for (long c = 0; c < noise.cols(); ++c) {
          noise(r, c) = gauss();
        }
      }
      mode.extended = mode.propagate * mode.extended + mode.noise * noise;
      for (long col = 0; col < nFree_; ++col) {
        const long a = freeIndex_[static_cast<size_t>(col)];
        pnm_[static_cast<size_t>(k)][a] =
            mode.extended(0, col) * std::sqrt(mass_[static_cast<size_t>(a)]);
      }
    } else {
      const double tau = 1.0 / (2.0 * opt_.pileScale * omegaK_[k]);
      langevin(k, tau);
    }
  }
  fromNormal(pnm_, p_);
}

toNormal()

void eonc::pathintegral::RingPolymer::toNormal ( const std::vector< VectorXd > & src, std::vector< VectorXd > & dst ) const

private

Definition at line 534 of file PathIntegral.cpp.

                                                           {
  MatrixXd packed(nDof_, nBeads_);
  for (long j = 0; j < nBeads_; ++j) {
    packed.col(j) = src[static_cast<size_t>(j)];
  }
  const MatrixXd out = packed * modes_.transpose();
  for (long k = 0; k < nBeads_; ++k) {
    dst[static_cast<size_t>(k)] = out.col(k);
  }
}

Member Data Documentation

atomicNumbers_

std::vector<int> eonc::pathintegral::RingPolymer::atomicNumbers_

private

Definition at line 149 of file PathIntegral.h.

batches_

long eonc::pathintegral::RingPolymer::batches_ {0}

private

Definition at line 174 of file PathIntegral.h.

{0};

constrain_

bool eonc::pathintegral::RingPolymer::constrain_ {false}

private

Definition at line 161 of file PathIntegral.h.

{false};

f_

std::vector<VectorXd> eonc::pathintegral::RingPolymer::f_

private

Definition at line 156 of file PathIntegral.h.

forceSum_

double eonc::pathintegral::RingPolymer::forceSum_ {0.0}

private

Definition at line 177 of file PathIntegral.h.

{0.0};

free_

std::vector<char> eonc::pathintegral::RingPolymer::free_

private

Definition at line 150 of file PathIntegral.h.

freeIndex_

std::vector<long> eonc::pathintegral::RingPolymer::freeIndex_

private

Definition at line 151 of file PathIntegral.h.

gle_

std::vector<ModeGle> eonc::pathintegral::RingPolymer::gle_

private

Definition at line 171 of file PathIntegral.h.

haveForces_

bool eonc::pathintegral::RingPolymer::haveForces_ {false}

private

Definition at line 162 of file PathIntegral.h.

{false};

kineticSum_

double eonc::pathintegral::RingPolymer::kineticSum_ {0.0}

private

Definition at line 176 of file PathIntegral.h.

{0.0};

mass_

std::vector<double> eonc::pathintegral::RingPolymer::mass_

private

Definition at line 148 of file PathIntegral.h.

modes_

MatrixXd eonc::pathintegral::RingPolymer::modes_

private

Definition at line 152 of file PathIntegral.h.

nAtoms_

long eonc::pathintegral::RingPolymer::nAtoms_ {0}

private

Definition at line 144 of file PathIntegral.h.

{0};

nBeads_

long eonc::pathintegral::RingPolymer::nBeads_ {0}

private

Definition at line 146 of file PathIntegral.h.

{0};

nDof_

long eonc::pathintegral::RingPolymer::nDof_ {0}

private

Definition at line 145 of file PathIntegral.h.

{0};

nFree_

long eonc::pathintegral::RingPolymer::nFree_ {0}

private

Definition at line 147 of file PathIntegral.h.

{0};

omegaK_

VectorXd eonc::pathintegral::RingPolymer::omegaK_

private

Definition at line 153 of file PathIntegral.h.

opt_

Options eonc::pathintegral::RingPolymer::opt_

private

Definition at line 143 of file PathIntegral.h.

p_

std::vector<VectorXd> eonc::pathintegral::RingPolymer::p_

private

Definition at line 155 of file PathIntegral.h.

planeNormal_

VectorXd eonc::pathintegral::RingPolymer::planeNormal_

private

Definition at line 159 of file PathIntegral.h.

planeOrigin_

VectorXd eonc::pathintegral::RingPolymer::planeOrigin_

private

Definition at line 160 of file PathIntegral.h.

pnm_

std::vector<VectorXd> eonc::pathintegral::RingPolymer::pnm_

private

Definition at line 158 of file PathIntegral.h.

q_

std::vector<VectorXd> eonc::pathintegral::RingPolymer::q_

private

Definition at line 154 of file PathIntegral.h.

qnm_

std::vector<VectorXd> eonc::pathintegral::RingPolymer::qnm_

private

Definition at line 157 of file PathIntegral.h.

recorded_

long eonc::pathintegral::RingPolymer::recorded_ {0}

private

Definition at line 175 of file PathIntegral.h.

{0};

rng_

std::uint64_t eonc::pathintegral::RingPolymer::rng_

private

Definition at line 173 of file PathIntegral.h.

---

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

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