PathIntegral.h

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** i-PI Copyright (C) 2014-2015 i-PI developers
**
** Permission is hereby granted, free of charge, to any person obtaining
** a copy of this software and associated documentation files (the
** "Software"), to deal in the Software without restriction, including
** without limitation the rights to use, copy, modify, merge, publish,
** distribute, sublicense, and/or sell copies of the Software, and to
** permit persons to whom the Software is furnished to do so, subject to
** the following conditions:
**
** The above copyright notice and this permission notice shall be
** included in all copies or substantial portions of the Software.
**
** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
** NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
** BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
** ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
** CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
** SOFTWARE.
**
** SPDX-License-Identifier: MIT
*/
#pragma once
 
#include "eon/Eigen.h"
#include "eon/Potential.h"
 
#include <cstdint>
#include <string>
#include <vector>
 
namespace eonc::pathintegral {

enum class Springs { Trotter, Eco };

enum class Thermostat { Pile, Piglet };
 
struct Options {
  long beads{8};
  double temperature{1.0};
  double kB{1.0};
  double hbar{1.0};
  double dt{0.005};
  Springs springs{Springs::Trotter};
  Thermostat thermostat{Thermostat::Pile};
  double pileTau{0.2};
  double pileScale{1.0};
  double ecoOmegaMax{0.0};
  std::string gleFile;
  std::uint64_t seed{1};
};
 
struct Sample {
  double kineticCv{0.0};
  double meanForce{0.0};
  long batches{0};
};

void requireTrotterSprings(const std::string &springs, const char *use);

VectorXd trotterEigenvalues(long nBeads);
VectorXd ecoEigenvalues(long nBeads, double xmax);

MatrixXd normalModeMatrix(long nBeads);

class RingPolymer {
public:
  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);
  [[nodiscard]] const std::vector<VectorXd> &beads() const { return q_; }
  [[nodiscard]] const std::vector<VectorXd> &momenta() const { return p_; }
  void setMomenta(const std::vector<VectorXd> &momenta);
  void thermalMomenta();

  void setHyperplane(const VectorXd &normal, const VectorXd &origin);
 
  void step(Potential &pot, const double *box, bool record);
  void nveStep(Potential &pot, const double *box);
 
  [[nodiscard]] double kineticCv() const;
  [[nodiscard]] double meanForce() const;
  [[nodiscard]] long batches() const { return batches_; }
  [[nodiscard]] VectorXd centroid() const;
  [[nodiscard]] VectorXd centroidVelocity() const;
 
  void resetAverages();
 
  Sample sample(Potential &pot, const double *box, long equilibration,
                long production);
 
private:
  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();
 
  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};
 
  struct ModeGle {
    MatrixXd drift;
    MatrixXd covariance;
    MatrixXd propagate;
    MatrixXd noise;
    MatrixXd extended;
  };
  std::vector<ModeGle> gle_;
 
  std::uint64_t rng_;
  long batches_{0};
  long recorded_{0};
  double kineticSum_{0.0};
  double forceSum_{0.0};
};
 
} // namespace eonc::pathintegral