Matter.h

Go to the documentation of this file.

/*
** This file is part of eOn.
**
** SPDX-License-Identifier: BSD-3-Clause
**
** Copyright (c) 2010--present, eOn Development Team
** All rights reserved.
**
** Repo:
** https://github.com/TheochemUI/eOn
*/
#pragma once
#include "ConFileIO.h"
#include "Eigen.h"
#include "EonLogger.h"
#include "PbcSimd.h"
#include "Potential.h"
#include "StructureComparisonOptions.h"
#include "SurrogatePotential.h"
#include <array>
#include <cstdint>
#include <memory>
#include <span>
#include <string>
#include <vector>
 
// This is a forward declaration of BondBoost to avoid a circular dependency.
namespace eonc {
class BondBoost;
class Parameters;
 
// Position / difference MIC conventions (issue #176). Legacy matches historical
// applyPeriodicBoundary (fmod to [0,1) fractional). MinimumImage matches the
// existing eonc::pbc::apply path used for interatomic distances (frac in
// [-0.5,0.5)).
enum class PbcConvention {
  Legacy = 0,
  MinimumImage = 1,
};
 
namespace pbc {
 
// Minimum-image on a difference vector (fractional in [-0.5, 0.5)).
inline AtomMatrix apply(const AtomMatrix &diff, const Matrix3d &cell,
                        const Matrix3d &cellInverse) {
  // Fractional coordinates, minimum-image wrap to [-0.5, 0.5), back to
  // Cartesian. x - floor(x + 0.5), Highway when the build has it.
  AtomMatrix frac = diff * cellInverse;
  wrapMinimumImage(frac.data(), static_cast<size_t>(frac.size()));
  return frac * cell;
}
 
// Legacy position wrap: fractional coords into [0, 1) via fmod (historical
// Matter).
inline AtomMatrix applyLegacy(const AtomMatrix &coords, const Matrix3d &cell,
                              const Matrix3d &cellInverse) {
  AtomMatrix frac = coords * cellInverse;
  // Same values as fmod(x + 1, 1): (x + 1) - trunc(x + 1).
  wrapLegacyUnit(frac.data(), static_cast<size_t>(frac.size()));
  return frac * cell;
}
 
// Position wrap selecting convention (Legacy vs MinimumImage centering).
inline AtomMatrix applyPositions(const AtomMatrix &coords, const Matrix3d &cell,
                                 const Matrix3d &cellInverse,
                                 PbcConvention convention) {
  if (convention == PbcConvention::MinimumImage) {
    return apply(coords, cell, cellInverse);
  }
  return applyLegacy(coords, cell, cellInverse);
}
 
inline VectorXd applyV(const VectorXd &diffVector, const Matrix3d &cell,
                       const Matrix3d &cellInverse) {
  AtomMatrix pbcMatrix =
      apply(AtomMatrix::Map(diffVector.data(), diffVector.size() / 3, 3), cell,
            cellInverse);
  return VectorXd(VectorXd::Map(pbcMatrix.data(), diffVector.size()));
}
 
} // namespace pbc
 
/* Data describing an atomic structure. This class has been devised to handle
 * information about an atomic structure such as positions, velocities, masses,
 * etc. It also allow to associate a forcefield for the structure through a
 * pointer to function (potential()). The class can read and save data to a .con
 * file (atom2con() and con2atom()). It can also save to a .xyz file
 * (atom2xyz()).*/
 
class Matter {
public:
  ~Matter();
  Matter(std::shared_ptr<Potential> pot, const Parameters &params);
  Matter(const Matter &matter);                  // create a copy of matter
  const Matter &operator=(const Matter &matter); // copy the matter object
  Matter(Matter &&other) noexcept;
  Matter &operator=(Matter &&other) noexcept;
  bool compare(const Matter &matter, bool indistinguishable = false);
 
  double
  distanceTo(const Matter &matter); // the distance to the given matter object
  double perAtomNorm(const Matter &matter); // the maximum distance between two
                                            // atoms in the Matter objects
  void
  setPotential(std::shared_ptr<Potential> pot); // set potential function to use
  std::shared_ptr<Potential> getPotential();    // get potential function to use
  void resize(long int nAtoms);   // set or reset the number of atoms
  long int numberOfAtoms() const; // return the number of atoms
  Matrix3d getCell() const;
  void setCell(const Matrix3d &newCell);
  double getPosition(long int atom, int axis)
      const; // return the position of an atom along one of the axis
  void setPosition(
      long int atom, int axis,
      double position); // set the position of atom along axis to position
  void setVelocity(
      long int atom, int axis,
      double velocity); // set the velocity of atom along axis to velocity
  bool relax(bool quiet = false, bool writeMovie = false,
             bool checkpoint = false, std::string prefixMovie = std::string(),
             std::string prefixCheckpoint = std::string(),
             bool retainMovieFrames = false);

  [[nodiscard]] const std::vector<readcon::ConFrame> &movieFrames() const {
    return movie_frames_;
  }
  void clearMovieFrames() { movie_frames_.clear(); }
  [[nodiscard]] std::vector<readcon::ConFrame> takeMovieFrames() {
    return std::move(movie_frames_);
  }
 
  AtomMatrix pbc(const AtomMatrix &diff) const;
  VectorXd pbcV(const VectorXd &diff) const;
 
  size_t getPotentialCalls() const;
  const AtomMatrix &getPositions() const; // return coordinates of atoms
  AtomMatrix getPositionsCopy() const;    // return a modifiable copy
  VectorXd getPositionsV() const;
  AtomMatrix
  getPositionsFree() const; // return coordinates of free atoms in array pos
  VectorXd getPositionsFreeV() const;
  void
  setPositions(const AtomMatrix &pos); // update Matter with the new positions
                                       // of the free atoms given in array pos
  void setPositionsV(const VectorXd &pos);
  void setPositionsFree(
      const AtomMatrix &pos); // update Matter with the new positions of the
                              // free atoms given in array pos
  void setPositionsFreeV(const VectorXd &pos);
 
  AtomMatrix getVelocities() const;
  void setVelocities(const AtomMatrix &v);
  void setBiasForces(const AtomMatrix &bf);
  void setBiasPotential(BondBoost *bondBoost);
  BondBoost *getBiasPotential() const;
  void assignKeepingBias(const Matter &other);
  void setForces(const AtomMatrix &f);
  AtomMatrix getAccelerations();
 
  const AtomMatrix &getForces() const;
  const AtomMatrix &getForcesRaw() const;
  AtomMatrix getBiasForces();
  VectorXd getForcesV() const;
  // Const so callers can read free-atom forces without mutating Matter (#171).
  AtomMatrix getForcesFree() const;
  VectorXd getForcesFreeV() const;
 
  PbcConvention getPbcConvention() const { return pbcConvention; }
  void setPbcConvention(PbcConvention convention) {
    pbcConvention = convention;
  }
 
  double getMass(long int atom) const; // return the mass of the atom specified
  void setMass(long int atom, double mass); // set the mass of an atom
  void setMasses(const VectorXd &massesIn); // set the mass of an atom
  long getAtomicNr(
      long int atom) const; // return the atomic number of the atom specified
  void setAtomicNr(long int atom,
                   long atomicNr);           // set the atomic number of an atom
  VectorXi getAtomicNrs() const;             // Get the vector of atomic numbers
  VectorXi getAtomicNrsFree() const;         // Get the vector of atomic numbers
  void setAtomicNrs(const VectorXi &atmnrs); // set the vector of atomic numbers

  int getFixed(long int atom) const;
  int getFixed(long int atom, int axis) const;
  [[nodiscard]] std::array<bool, 3> getFixedMask(long int atom) const;
  void setFixed(long int atom, int isFixed);
  void setFixed(long int atom, int axis, int isFixed);
  void setFixedMask(long int atom, std::array<bool, 3> mask);
  double getEnergyVariance() const;
  double getPotentialEnergy() const;
  Matrix3d cauchyStress();

  [[nodiscard]] bool needsForceUpdate() const { return recomputePotential; }

  double *forcesData();
  void setEvaluation(const AtomMatrix &forces, double energy) {
    restoreFileForces(forces, true, energy);
  }

  void setComputedPotential(double energy, double variance);
  double getKineticEnergy() const;
  double getMechanicalEnergy() const;
 
  double distance(long index1, long index2)
      const; // return the distance between two atoms in same configuration
  double pdistance(long index1, long index2, int axis) const;
  double distance(const Matter &matter, long index)
      const; // the distance between the same atom in two cofigurations
 
  long int
  numberOfFreeAtoms() const; // return the number of free (or movable) atoms
  long int numberOfFixedAtoms() const; // return the number of fixed atoms
 
  long
  getForceCalls() const; // return how many force calls that have been performed
  void resetForceCalls(); // zeroing the value of force calls
 
  double maxForce(void) const;

  [[nodiscard]] double maxFreeAtomForce(const AtomMatrix &rows) const;

  [[nodiscard]] bool getWriteConForces() const noexcept;
 
  // I/O delegates to eonc::io free functions (IoStatus for bindings).
  [[nodiscard]] io::IoStatus writeTibble(std::string filename) {
    return io::writeTibble(*this, filename);
  }
  [[nodiscard]] io::IoStatus con2matter(std::string filename) {
    return io::con2matter(*this, filename);
  }
  [[nodiscard]] io::IoStatus
  con2matter(const readcon::ConFrame &frame,
             io::ConFrameMetadata *out_metadata = nullptr) {
    return io::con2matter(*this, frame, out_metadata);
  }
  [[nodiscard]] io::IoStatus convel2matter(std::string filename) {
    return io::convel2matter(*this, filename);
  }
  [[nodiscard]] io::IoStatus
  matter2con(std::string filename, bool append = false,
             const io::ConFrameMetadata *metadata = nullptr) {
    return io::matter2con(*this, filename, append, metadata);
  }
  [[nodiscard]] io::IoStatus matter2convel(std::string filename) {
    return io::matter2convel(*this, filename);
  }
  [[nodiscard]] io::IoStatus matter2xyz(std::string filename,
                                        bool append = false) {
    return io::matter2xyz(*this, filename, append);
  }
 
  AtomMatrix getFree() const;
  VectorXd getFreeV() const;
  Eigen::Matrix<double, Eigen::Dynamic, 1> getMasses() const;

  [[nodiscard]] std::int64_t getAtomIndex(long int atom) const;
  void setAtomIndex(long int atom, std::int64_t index);

  [[nodiscard]] long mapFileRow(long file_row) const {
    if (file_row < 0 || file_row >= nAtoms ||
        fileToMatter.size() != static_cast<size_t>(nAtoms)) {
      return file_row;
    }
    return fileToMatter[static_cast<size_t>(file_row)];
  }
  void setFileToMatter(std::vector<long> map) { fileToMatter = std::move(map); }

  [[nodiscard]] const std::array<std::string, 5> &getHeaderCon() const {
    return headerCon;
  }
  void setHeaderCon(const std::array<std::string, 5> &headers) {
    headerCon = headers;
  }
  void setHeaderConLine(size_t i, std::string line) {
    headerCon.at(i) = std::move(line);
  }
 
  [[nodiscard]] bool getPeriodic() const noexcept {
    return usePeriodicBoundaries;
  }
  void setPeriodic(bool periodic) {
    usePeriodicBoundaries = periodic;
    recomputePotential = true;
    recomputeMaskedForces = true;
  }

  void applyPeriodicBoundaryIfEnabled() {
    if (usePeriodicBoundaries) {
      applyPeriodicBoundary();
    }
  }
 
private:
  // Friend declarations for eonc::io free functions that need private access
  // con2matter still needs private headerCon + recompute flags for faithful
  // force/energy restore without an extra potential evaluation.
  friend io::IoStatus io::con2matter(Matter &, const readcon::ConFrame &,
                                     io::ConFrameMetadata *);
 
  eonc::log::Scoped m_log;
  std::shared_ptr<Potential>
      potential; // pointer to function calculating the energy and forces
  bool usePeriodicBoundaries; // boolean telling periodic boundaries are used

  void assertIsolatedMoleculeLayoutSafe() const;
  PbcConvention pbcConvention;     // position-wrap MIC convention (issue #176)
  mutable bool recomputePotential; // boolean indicating if the potential energy
                                   // and forces need to be recalculated
  mutable long
      forceCalls; // keep track of how many force calls have been performed
 
  // CON file header lines (indices 0-4 map to old headerCon1,2,4,5,6)
  std::array<std::string, 5> headerCon;
 
  void computePotential() const;
  void applyPeriodicBoundary();
  void applyPeriodicBoundary(double &component, int axis);
  void applyPeriodicBoundary(AtomMatrix &diff);
 
  // Narrowed from Parameters: only fields Matter actually reads
  bool removeNetForce{true};
  StructureComparisonOptions structComp;
  // Full Parameters pointer retained solely for relax() delegation
  const Parameters *parameters;
  long nAtoms;
  struct Impl;
  std::unique_ptr<Impl> impl_;
  BondBoost *biasPotential{nullptr};
  std::vector<long> fileToMatter; // CON file row -> Matter row
  void restoreFileForces(const AtomMatrix &fileForces, bool trustEnergy,
                         double energy);
  mutable std::vector<int> freeIndices; // cached indices of free atoms
  mutable bool recomputeFreeMask{true};
  mutable bool recomputeMaskedForces{true};
  mutable double energyVariance;
  std::vector<readcon::ConFrame> movie_frames_;
  mutable double potentialEnergy;
};

void evaluateTogether(Potential &pot, std::span<Matter *const> systems);
 
} // namespace eonc