fts/perturbation/EinsteinCrystalPerturbation.h

#ifndef RP_EINSTEIN_CRYSTAL_PERTURBATION_H
#define RP_EINSTEIN_CRYSTAL_PERTURBATION_H
 
#include <util/containers/DArray.h>  // member
#include <iostream>
 
namespace Pscf {
namespace Rp {
 
   using namespace Util;

   template <int D, class T>
   class EinsteinCrystalPerturbation : public T::Perturbation
   {
 
   public:
 
      // Protected constructor and destructor (see below).

      virtual void readParameters(std::istream& in);

      virtual void setup();

      virtual double hamiltonian(double unperturbedHamiltonian);

      virtual void incrementDc(DArray<typename T::RField>& dc);

      virtual double df();

      virtual void saveState();

      virtual void restoreState();
 
   protected:

      EinsteinCrystalPerturbation(typename T::Simulator& simulator);

      ~EinsteinCrystalPerturbation() = default;
 
      // Alias for base class.
      using PerturbationT = typename T::Perturbation;
 
      // Inherited protected members (selected).
      using PerturbationT::lambda_;
      using PerturbationT::simulator;
      using PerturbationT::system;
 
   private:
 
      using RFieldT = typename T::RField;
 
      // Parameters used in Einstein crystal integration
      DArray<double> epsilon_;
 
      // Reference w field
      DArray< RFieldT > w0_;
 
      // Eigenvector components of the reference w fields
      DArray< RFieldT > wc0_;
 
      // Work space
      RFieldT dw_;
 
      // Current Einstein crystal Hamiltonian
      double ecHamiltonian_;
 
      // Current unperturbed Hamiltonian
      double unperturbedHamiltonian_;
 
      // Saved Einstein crystal Hamiltonian
      double stateEcHamiltonian_;
 
      // Saved unperturbed Hamiltonian
      double stateUnperturbedHamiltonian_;
 
      // Reference field file name
      std::string fieldFileName_;
 
      // Have epsilon values been set?
      bool hasEpsilon_;
 
      // Compute eigenvector components of the reference field.
      void computeWcReference();
 
   };
 
}
}
#endif