UnitCellBase.h

#ifndef PSCF_UNIT_CELL_BASE_H
#define PSCF_UNIT_CELL_BASE_H
 
/*
* PSCF - Polymer Self-Consistent Field Theory
*
* Copyright 2016 - 2022, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
#include <util/containers/FArray.h>
#include <util/containers/FSArray.h>
#include <util/containers/FMatrix.h>
#include <pscf/math/IntVec.h>
#include <pscf/math/RealVec.h>
#include <util/math/Constants.h>
 
namespace Pscf
{
 
   using namespace Util;
 
   template <int D>
   class UnitCellBase
   {
   public:
 
      UnitCellBase();
 
      ~UnitCellBase();
 
      void setParameters(FSArray<double, 6> const & parameters);
 
      virtual double ksq(IntVec<D> const & k) const;
 
      virtual double dksq(IntVec<D> const & vec, int n) const;
 
      int nParameter() const;
 
      FSArray<double, 6> parameters() const;
 
      double parameter(int i) const;
 
      const RealVec<D>& rBasis(int i) const;
 
      const RealVec<D>& kBasis(int i) const;
 
      double drBasis(int k, int i, int j) const;
 
      double dkBasis(int k, int i, int j) const;
 
      double drrBasis(int k, int i, int j) const;
 
      double dkkBasis(int k, int i, int j) const;
 
      bool isInitialized() const
      {  return isInitialized_; }
 
   protected:
 
      FArray<RealVec<D>, D> rBasis_;
 
      FArray<RealVec<D>, D> kBasis_;
 
      FArray<FMatrix<double, D, D>, 6> drBasis_;
 
      FArray<FMatrix<double, D, D>, 6> dkBasis_;
 
      FArray<FMatrix<double, D, D>, 6> drrBasis_;
 
      FArray<FMatrix<double, D, D>, 6> dkkBasis_;
 
      FArray<double, 6> parameters_;
 
      int nParameter_;
 
      bool isInitialized_;
 
      void setLattice();
 
   private:
 
      void initializeToZero();
 
      virtual void setBasis() = 0;
 
      void computeDerivatives();
 
   };
 
   /*
   * Get the number of Parameters in the unit cell.
   */
   template <int D>
   inline
   int UnitCellBase<D>::nParameter() const
   {  return nParameter_;  }
 
   /*
   * Get the unit cell parameters.
   */
   template <int D>
   inline
   FSArray<double, 6> UnitCellBase<D>::parameters() const
   {
      FSArray<double, 6> parameters;
      for (int i = 0; i < nParameter_; ++i) {
         parameters.append(parameters_[i]);
      }
      return parameters;
   }
 
   /*
   * Get a single parameter of the unit cell.
   */
   template <int D>
   inline
   double UnitCellBase<D>::parameter(int i) const
   {  return parameters_[i]; }
 
   /*
   * Get Bravais basis vector i.
   */
   template <int D>
   inline
   const RealVec<D>& UnitCellBase<D>::rBasis(int i) const
   {  return rBasis_[i];  }
 
   /*
   * Get reciprocal basis vector i.
   */
   template <int D>
   inline
   const RealVec<D>& UnitCellBase<D>::kBasis(int i) const
   {  return kBasis_[i];  }
 
   /*
   * Get component j of derivative of r basis vector i w/respect to param k.
   */
   template <int D>
   inline
   double UnitCellBase<D>::drBasis(int k, int i, int j) const
   {  return drBasis_[k](i,j);  }
 
   /*
   * Get component j of derivative of r basis vector i w/respect to param k.
   */
   template <int D>
   inline
   double UnitCellBase<D>::dkBasis(int k, int i, int j) const
   {  return dkBasis_[k](i, j);  }
 
   /*
   * Get the derivative of ki*kj with respect to parameter k.
   */
   template <int D>
   inline
   double UnitCellBase<D>::dkkBasis(int k, int i, int j) const
   {  return dkkBasis_[k](i, j);  }
 
   /*
   * Get the derivative of ri*rj with respect to parameter k.
   */
   template <int D>
   inline
   double UnitCellBase<D>::drrBasis(int k, int i, int j) const
   {  return drrBasis_[k](i, j);  }
 
   // Non-inline member functions
 
   /*
   * Constructor.
   */
   template <int D>
   UnitCellBase<D>::UnitCellBase()
    : nParameter_(0),
      isInitialized_(false)
   {}
 
   /*
   * Destructor.
   */
   template <int D>
   UnitCellBase<D>::~UnitCellBase()
   {}
 
   /*
   * Set all the parameters in the unit cell.
   */
   template <int D>
   void UnitCellBase<D>::setParameters(FSArray<double, 6> const& parameters)
   {
      UTIL_CHECK(parameters.size() == nParameter_);
      isInitialized_ = false;
      for (int i = 0; i < nParameter_; ++i) {
         parameters_[i] = parameters[i];
      }
      setLattice();
   }
   
   /*
   * Get square magnitude of reciprocal basis vector.
   */
   template <int D>
   double UnitCellBase<D>::ksq(IntVec<D> const & k) const
   {
      RealVec<D> g(0.0);
      RealVec<D> p;
      for (int i = 0; i < D; ++i) {
         p.multiply(kBasis_[i], k[i]);
         g += p;
      }
      double value = 0.0;
      for (int i = 0; i < D; ++i) {
         value += g[i]*g[i];
      }
      return value;
   }
 
   /*
   * Get magnitude of derivative of square of reciprocal basis vector.
   */
   template <int D>
   double UnitCellBase<D>::dksq(IntVec<D> const & vec, int n) const
   {
      double element = 0.0;
      double value = 0.0;
 
      for (int p = 0; p < D; ++p){
         for (int q = 0; q < D; ++q){
            element = dkkBasis(n, p, q);
            value += vec[p]*vec[q]*element;
         }
      }
 
      return value;
   }
 
   // Protected member functions
 
   /*
   * Initialize internal arrays to zero.
   */
   template <int D>
   void UnitCellBase<D>::initializeToZero()
   {
      // Initialize all elements to zero
      int i, j, k;
      for (i = 0; i < D; ++i) {
         for (j = 0; j < D; ++j) {
            rBasis_[i][j] = 0.0;
            kBasis_[i][j] = 0.0;
         }
      }
      for (k = 0; k < 6; ++k){
         for (i = 0; i < D; ++i) {
            for (j = 0; j < D; ++j) {
               drBasis_[k](i,j) = 0.0;
               dkBasis_[k](i,j) = 0.0;
               drrBasis_[k](i,j) = 0.0;
               dkkBasis_[k](i,j) = 0.0;
            }
         }
      }
   }
 
   /*
   * Compute quantities involving derivatives.
   */
   template <int D>
   void UnitCellBase<D>::computeDerivatives()
   {
      // Compute dkBasis
      int p, q, r, s, t;
      for (p = 0; p < nParameter_; ++p) {
         for (q = 0; q < D; ++q) {
            for (r = 0; r < D; ++r) {
 
               // Loop over free indices s, t
               for (s = 0; s < D; ++s) {
                  for (t = 0; t < D; ++t) {
                     dkBasis_[p](q,r)
                       -= kBasis_[q][s]*drBasis_[p](t,s)*kBasis_[t][r];
                  }
               }
               dkBasis_[p](q,r) /= 2.0*Constants::Pi;
 
            }
         }
      }
 
      // Compute drrBasis and dkkBasis
      for (p = 0; p < nParameter_; ++p) {
         for (q = 0; q < D; ++q) {
            for (r = 0; r < D; ++r) {
               for (s = 0; s < D; ++s) {
                  drrBasis_[p](q,r) += rBasis_[q][s]*drBasis_[p](r,s);
                  drrBasis_[p](q,r) += rBasis_[r][s]*drBasis_[p](q,s);
                  dkkBasis_[p](q,r) += kBasis_[q][s]*dkBasis_[p](r,s);
                  dkkBasis_[p](q,r) += kBasis_[r][s]*dkBasis_[p](q,s);
 
               }
            }
         }
      }
 
   }
 
   /*
   * Set all lattice parameters.
   */
   template <int D>
   void UnitCellBase<D>::setLattice()
   {
      initializeToZero();
      setBasis();
      computeDerivatives();
      isInitialized_ = true;
   }
 
}
#endif