OrthorhombicBoundary.h

#ifndef SIMP_ORTHORHOMBIC_BOUNDARY_H
#define SIMP_ORTHORHOMBIC_BOUNDARY_H

/*
* Util Package - C++ Utilities for Scientific Computation
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/

#include "OrthoRegion.h"                // base class
#include <util/crystal/LatticeSystem.h> // member
#include <util/containers/FSArray.h>    // member template
#include <util/space/Vector.h>          // member template argument
#include <util/space/IntVector.h>       // inline functions
#include <util/space/Dimension.h>       // member template argument
#include <util/global.h>                // asserts in inline functions

#include <cmath>
#include <iostream>

class OrthorhombicBoundaryTest;

namespace Util {
   class Random;
}

namespace Simp
{

   using namespace Util;

   class OrthorhombicBoundary : private OrthoRegion
   {

   public:

      OrthorhombicBoundary();

      // Use default copy constructor.

      void setOrthorhombic(const Vector &lengths);

      void setTetragonal(double a, double bc);

      void setCubic(double a);

      template <class Archive>
      void serialize(Archive& ar, const unsigned int version);



      void shift(Vector &r) const;

      void shift(Vector &r, IntVector& shift) const;

      void applyShift(Vector &r, int i, int t) const;

      void shiftGen(Vector &r) const;

      void shiftGen(Vector &r, IntVector& shift) const;



      double distanceSq(const Vector &r1, const Vector &r2) const;

      double distanceSq(const Vector &r1, const Vector &r2, IntVector& shift) 
      const;

      double distanceSq(const Vector &r1, const Vector &r2, Vector &dr) const;


     
      bool isMinImageGen(const Vector &dr);

      bool isMinImageCart(const Vector &dr);



      void transformCartToGen(const Vector& Rc, Vector& Rg) const;

      void transformGenToCart(const Vector& Rg, Vector& Rc) const;



      LatticeSystem latticeSystem();

      const Vector& lengths() const;

      double length(int i) const;

      double minLength() const;

      double volume() const;

      const Vector& bravaisBasisVector(int i) const;

      const Vector& reciprocalBasisVector(int i) const;

      void randomPosition(Random &random, Vector &r) const;

      bool isValid();


   private:

      FSArray<Vector, Dimension>  bravaisBasisVectors_;

      FSArray<Vector, Dimension>  reciprocalBasisVectors_;

      Vector invLengths_;

      double minLength_;

      LatticeSystem lattice_;

   // friends:

      friend class ::OrthorhombicBoundaryTest;

      friend std::istream& operator >> (std::istream& in, 
                                        Simp::OrthorhombicBoundary& boundary);

      friend std::ostream& operator << (std::ostream& out, 
                                        const Simp::OrthorhombicBoundary& boundary);

      void reset();

   };

   // Friend function declarations

   std::istream& operator >> (std::istream& in, OrthorhombicBoundary& boundary);

   std::ostream& operator << (std::ostream& out, 
                              const OrthorhombicBoundary& boundary);

   // Inline member function definitions

   /* 
   * Return Vector of lengths by const reference.
   */
   inline const Vector& OrthorhombicBoundary::lengths() const 
   {  return lengths_; }

   /* 
   * Get length = maximum - minimum in direction i.
   */
   inline double OrthorhombicBoundary::length(int i) const 
   {  return lengths_[i]; }

   /* 
   * Return the maximum validity range of the distances.
   */
   inline double OrthorhombicBoundary::minLength() const
   {  return minLength_; }

   /* 
   * Return region volume.
   */
   inline double OrthorhombicBoundary::volume() const
   {  return volume_; }

   /* 
   * Return Bravais lattice basis vector number i.
   */
   inline const Vector& OrthorhombicBoundary::bravaisBasisVector(int i) const
   {  return bravaisBasisVectors_[i]; }

   /* 
   * Return reciprocal lattice basis vector number i.
   */
   inline const Vector& OrthorhombicBoundary::reciprocalBasisVector(int i) const
   {  return reciprocalBasisVectors_[i]; }

   /* 
   * Shift Cartesian Vector r to primitive unit cell.
   */
   inline void OrthorhombicBoundary::shift(Vector& r) const
   {
      for (int i = 0; i < Dimension; ++i) {
         if( r[i] >= maxima_[i] ) {
            r[i] = r[i] - lengths_[i];
            assert(r[i] < maxima_[i]);
         } else
         if ( r[i] <  minima_[i] ) {
            r[i] = r[i] + lengths_[i];
            assert(r[i] >= minima_[i]);
         }
      }
   }

   /* 
   * Shift Vector r to periodic cell, R[axis] < r[axis] < maxima_[axis].
   */
   inline void OrthorhombicBoundary::shift(Vector& r, IntVector& shift) const
   {
      for (int i = 0; i < Dimension; ++i) {
         if (r[i] >= maxima_[i]) {
            r[i] = r[i] - lengths_[i];
            ++(shift[i]);
            assert(r[i] < maxima_[i]);
         } else
         if (r[i] <  minima_[i]) {
            r[i] = r[i] + lengths_[i];
            --(shift[i]);
            assert(r[i] >= minima_[i]);
         }
      }
   }

   /*
   * Shift Cartesian Vector r by multiple t of a Bravais lattice vector.
   */
   inline void OrthorhombicBoundary::applyShift(Vector &r, int i, int t) const
   {  r[i] += t*lengths_[i]; }

   /* 
   * Shift generalized Vector r to primitive unit cell.
   */
   inline void OrthorhombicBoundary::shiftGen(Vector& r) const
   {
      for (int i = 0; i < Dimension; ++i) {
         if( r[i] >= 1.0 ) {
            r[i] = r[i] - 1.0;
            assert(r[i] < 1.0);
         } else
         if ( r[i] <  0.0 ) {
            r[i] = r[i] + 1.0;
            assert(r[i] >= 0.0);
         }
      }
   }

   /* 
   * Shift generalized Vector r to primitive cell, 0 < r[axis] < 1.0.
   */
   inline 
   void OrthorhombicBoundary::shiftGen(Vector& r, IntVector& shift) const
   {
      for (int i = 0; i < Dimension; ++i) {
         if (r[i] >= 1.0) {
            r[i] = r[i] - 1.0;
            ++(shift[i]);
            assert(r[i] < 1.0);
         } else
         if (r[i] <  0.0) {
            r[i] = r[i] + 1.0;
            --(shift[i]);
            assert(r[i] >= 0.0);
         }
      }
   }

   /* 
   * Calculate squared distance by minimum image convention.
   */
   inline 
   double OrthorhombicBoundary::distanceSq(const Vector &r1, const Vector &r2, 
                               IntVector& translate) const
   {
      double dr;
      double norm = 0.0;
      for (int i = 0; i < Dimension; ++i) {
         dr = r1[i] - r2[i];
         if ( fabs(dr) > halfLengths_[i] ) {
            if (dr > 0.0) {
               dr -= lengths_[i];
               translate[i] = -1;
            } else {
               dr += lengths_[i];
               translate[i] = +1;
            }
            assert(fabs(dr) <= halfLengths_[i]);
         } else {
            translate[i] = 0;
         }
         norm += dr*dr;
      }
      return norm;
   }

   /* 
   * Return squared distance and separation with minimum image convention.
   */
   inline 
   double OrthorhombicBoundary::distanceSq(const Vector &r1, const Vector &r2) const
   {
      double dr;
      double norm = 0.0;
      for (int i = 0; i < Dimension; ++i) {
         dr = r1[i] - r2[i];
         if (fabs(dr) > halfLengths_[i]) {
            if (dr > 0.0) {
               dr -= lengths_[i];
            } else {
               dr += lengths_[i];
            }
            assert(fabs(dr) <= halfLengths_[i]);
         }
         norm += dr*dr;
      }
      return norm;
   }

   /* 
   * Calculate squared distance between two positions, and separation vector,
   * using the minimum image convention.
   */
   inline 
   double OrthorhombicBoundary::distanceSq(const Vector &r1, const Vector &r2, 
                                           Vector &dr) const
   {
      for (int i = 0; i < Dimension; ++i) {
         dr[i] = r1[i] - r2[i];
         if (fabs(dr[i]) > halfLengths_[i]) {
            if (dr[i] > 0.0) {
               dr[i] -= lengths_[i];
            } else {
               dr[i] += lengths_[i];
            }
            assert(fabs(dr[i]) <= halfLengths_[i]);
         }
      }
      return ( dr[0]*dr[0] + dr[1]*dr[1] + dr[2]*dr[2] );
   }

   /*
   * Return actual lattice system.
   */
   inline LatticeSystem OrthorhombicBoundary::latticeSystem()
   { return lattice_; }

   /*
   * Transform Cartesian Vector Rc to generalized Vector Rg.
   */
   inline void 
   OrthorhombicBoundary::transformCartToGen(const Vector& Rc, Vector& Rg) 
   const
   {
      for (int i = 0; i < Dimension; ++i) {
         Rg[i] = Rc[i] * invLengths_[i];
      }
   }
      
   /*
   * Transform Cartesian Vector Rc to generalized Vector Rg.
   */
   inline void 
   OrthorhombicBoundary::transformGenToCart(const Vector& Rg, Vector& Rc) 
   const
   {
      for (int i = 0; i < Dimension; ++i) {
         Rc[i] = Rg[i] * lengths_[i];
      }
   }

   /*
   * Is a generalized separation vector a minimimum image of itself?
   */
   inline
   bool OrthorhombicBoundary::isMinImageGen(const Vector &dr)
   {
      for (int i = 0; i < Dimension; ++i) {
         if (fabs(dr[i]) > 0.5) {
            return false;
         }
      }
      return true;
   }

   /*
   * Is Cartesian separation vector dr a minimimum image of itself?
   */
   inline
   bool OrthorhombicBoundary::isMinImageCart(const Vector &dr)
   {
      for (int i = 0; i < Dimension; ++i) {
         if (fabs(dr[i]) > halfLengths_[i]) {
            return false;
         }
      }
      return true;
   }

   // Member function template

   /*
   * Serialize an OrthorhombicBoundary to/from an archive.
   */
   template <class Archive> void 
   OrthorhombicBoundary::serialize(Archive& ar, const unsigned int version)
   {
      OrthoRegion::serialize(ar, version);
      ar & bravaisBasisVectors_;
      ar & reciprocalBasisVectors_;
      ar & invLengths_;
      ar & minLength_;
      serializeEnum(ar, lattice_, version);
      if (Archive::is_loading()) {
         isValid();
      }
   }

}
 
#ifdef UTIL_MPI
#include <util/mpi/MpiSendRecv.h>
#include <util/mpi/MpiTraits.h>

namespace Util
{

   using namespace Simp; 

   template <>
   void send<Simp::OrthorhombicBoundary>(MPI::Comm& comm, 
             Simp::OrthorhombicBoundary& data, int dest, int tag);

   template <>
   void recv<Simp::OrthorhombicBoundary>(MPI::Comm& comm, 
             Simp::OrthorhombicBoundary& data, int source, int tag);

   template <>
   void bcast<Simp::OrthorhombicBoundary>(MPI::Intracomm& comm, 
              Simp::OrthorhombicBoundary& data, int root);

   template <>
   class MpiTraits<Simp::OrthorhombicBoundary>
   {
   public:
      static MPI::Datatype type;         
      static bool hasType;               
   };

}
#endif // ifdef  UTIL_MPI

#endif // ifndef UTIL_ORTHORHOMBIC_BOUNDARY_H