crystal.h

#ifndef CRYSTAL_H
#define CRYSTAL_H
 
#include "eigen.h"
#include <string>
#include <vector>
#include "context.h"
 
struct SymmetryOperation {
  Eigen::Matrix3d rotation;
  Eigen::Vector3d translation;
};
 
class Crystal {
protected:
  static Eigen::Matrix3d calcReciprocalCell(const Eigen::Matrix3d &directUnitCell);
 
  void setDirectUnitCell(const Eigen::Matrix3d &directUnitCell_);
 
  Eigen::Matrix3d directUnitCell;
  Eigen::Matrix3d reciprocalUnitCell;
  double volumeUnitCell;
  int numAtoms;
  int numSpecies;
  int dimensionality;
 
  // these vectors/matrices  running over the number of atoms
  Eigen::MatrixXd atomicPositions;         // size (numAtoms,3)
  Eigen::VectorXi atomicSpecies;           // size (numAtoms)
  std::vector<std::string> atomicNames;    // size (numAtoms)
  Eigen::VectorXd atomicMasses;            // size (numAtoms)
  Eigen::VectorXd atomicIsotopeCouplings;  // size (numAtoms)
  Eigen::Tensor<double, 3> bornCharges;    // size (numAtoms,3,3)
  Eigen::Matrix3d dielectricMatrix;        // size (3,3)
 
  // vectors running over the number of species
  std::vector<std::string> speciesNames;   // size (numSpecies)
  Eigen::VectorXd speciesMasses;           // size (numSpecies)
  Eigen::VectorXd speciesIsotopeCouplings; // size (numSpecies)
 
  // Untested for now
  std::vector<SymmetryOperation> symmetryOperations;
  int numSymmetries;
 
public:
  Crystal(Context &context, Eigen::Matrix3d &directUnitCell_, Eigen::MatrixXd &atomicPositions_,
          Eigen::VectorXi &atomicSpecies_,
          std::vector<std::string> &speciesNames_,
          Eigen::VectorXd &speciesMasses_, Eigen::Tensor<double, 3>& bornCharges_,
          Eigen::Matrix3d &dielectricMatrix);
 
  Crystal();
 
  Crystal(const Crystal &obj);
 
  Crystal &operator=(const Crystal &obj);
 
  void print() const;
 
  //  Setter and getter for all the variables above
 
  const Eigen::Matrix3d &getDirectUnitCell();
 
  const Eigen::Matrix3d &getReciprocalUnitCell();
 
  const int &getNumAtoms();
 
  const Eigen::Tensor<double, 3> getBornEffectiveCharges();
 
  const inline Eigen::Matrix3d getDielectricMatrix() { return dielectricMatrix; }
 
  double getVolumeUnitCell(int dimensionality_ = 3) const;
 
  const std::vector<SymmetryOperation> &getSymmetryOperations();
 
  int getNumSymmetries() const;
 
  const Eigen::MatrixXd &getAtomicPositions();
 
  const Eigen::VectorXi &getAtomicSpecies();
 
  const std::vector<std::string> &getAtomicNames();
 
  const Eigen::VectorXd &getAtomicMasses();
 
  const std::vector<std::string> &getSpeciesNames();
 
  const Eigen::VectorXd &getSpeciesMasses();
 
  int getDimensionality() const;
 
  int getNumSpecies() const;
 
  const Eigen::VectorXd &getAtomicIsotopeCouplings();
 
  Eigen::Vector3d crystalToCartesian(const Eigen::Vector3d &point);
 
  Eigen::Vector3d cartesianToCrystal(const Eigen::Vector3d &point);
 
  std::tuple<Eigen::MatrixXd, Eigen::VectorXd>
  buildWignerSeitzVectors(const Eigen::Vector3i &grid,
                          const int &superCellFactor = 2);
 
  std::tuple<Eigen::MatrixXd, Eigen::Tensor<double, 3>>
  buildWignerSeitzVectorsWithShift(const Eigen::Vector3i &grid,
                                   const Eigen::MatrixXd &shift,
                                   const int &superCellFactor = 2);
 
  void generateSymmetryInformation(Context &context);
 
 
};
 
#endif