MonoclinicBoundary.h

#ifndef SIMP_MONOCLINIC_BOUNDARY_H
#define SIMP_MONOCLINIC_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 <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 methods
#include <util/space/Dimension.h>       // member template argument
#include <util/global.h>                // asserts in inline methods

#include <cmath>
#include <iostream>

class MonoclinicBoundaryTest;

namespace Util { 
   class Random; 
   #ifdef UTIL_MPI
   template <class T> void send(MPI::Comm& comm, T& data, int dest, int tag);
   template <class T> void recv(MPI::Comm& comm, T& data, int source, int tag);
   template <class T> void bcast(MPI::Intracomm& comm, T& data, int root);
   #endif
}

namespace Simp
{

   using namespace Util;

   class MonoclinicBoundary 
   {

   public:

      MonoclinicBoundary();

      void setMonoclinic(const Vector &lengths, const double d);

      void setOrthorhombic(const Vector &lengths);

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



      void shift(Vector &r) const;

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

      void shiftGen(Vector &r) const;

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

      void applyShift(Vector &r, int i, int t) 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;



      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:

      Vector minima_;

      Vector maxima_;

      Vector l_;

      double d_;

      Vector halfL_;

      Vector invL_;

      Vector lengths_;

      double volume_;

      double c3_;

      double e_;

      double minLength_;

      FSArray<Vector, Dimension>  bravaisBasisVectors_;

      FSArray<Vector, Dimension>  reciprocalBasisVectors_;

      LatticeSystem lattice_;

   // friends:

      friend class ::MonoclinicBoundaryTest;

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

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

      #ifdef UTIL_MPI
      friend void Util::send<>(MPI::Comm& comm, Simp::MonoclinicBoundary& data, 
                               int dest, int tag);

      friend void Util::recv<>(MPI::Comm& comm, Simp::MonoclinicBoundary& data, 
                               int source, int tag);

      friend void Util::bcast<>(MPI::Intracomm& comm, Simp::MonoclinicBoundary& data, 
                                int root);
      #endif

      void reset();

   };

   // Inline method definitions

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

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

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

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

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

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

   /* 
   * Shift Vector r to periodic cell, R[axis] < r[axis] < maxima_[axis].
   */
   inline void MonoclinicBoundary::shift(Vector& r) const
   {
       if(r[0] >= maxima_[0]) {
          r[0] -= l_[0];
          assert(r[0] < maxima_[0]);
       } else
       if (r[0] <  minima_[0]) {
          r[0] += l_[0];
          assert(r[0] >= minima_[0]);
       }
       if(r[1] >= maxima_[1]) {
          r[1] -= l_[1];
          r[2] -= d_;
          assert(r[1] < l_[1]);
       } else
       if (r[1] <  minima_[1]) {
          r[1] += l_[1];
          r[2] += d_;
          assert(r[1] >= minima_[1]);
       }
       double u2 = r[2] + c3_*r[1];
       if(u2 >= maxima_[2]) {
          r[2] -= l_[2];
          assert(r[2] + c3_*r[1] >= minima_[2]);
       } else
       if (u2 <  minima_[2]) {
          r[2] += l_[2];
          assert(r[2] + c3_*r[1] >= minima_[2]);
       }
   }

   /* 
   * Shift Vector r to periodic cell, R[axis] < r[axis] < maxima_[axis].
   */
   inline void MonoclinicBoundary::shift(Vector& r, IntVector& shift) const
   {
       if (r[0] >= maxima_[0]) {
          r[0] -= l_[0];
          ++(shift[0]);
          assert(r[0] < maxima_[0]);
       } else
       if (r[0] <  minima_[0]) {
          r[0] += l_[0];
          --(shift[0]);
          assert(r[0] >= minima_[0]);
       }
       if (r[1] >= maxima_[1]) {
          r[1] -= l_[1];
          r[2] -= d_;
          ++(shift[1]);
          assert(r[1] < l_[1]);
       } else
       if (r[1] <  minima_[1]) {
          r[1] += l_[1];
          r[2] += d_;
          --(shift[1]);
          assert(r[1] >= minima_[1]);
       }
       double u2 = r[2] + c3_*r[1];
       if (u2 >= maxima_[2]) {
          r[2] -= l_[2];
          ++(shift[2]);
          assert(r[2] + c3_*r[1] < maxima_[2]);
       } else
       if (u2 <  minima_[2]) {
          r[2] += l_[2];
          --(shift[2]);
          assert(r[2] + c3_*r[1] >= minima_[2]);
       }
   }

   /*
   * Shift Cartesian Vector r by multiple t of a Bravais lattice vector.
   */
   inline void MonoclinicBoundary::applyShift(Vector &r, int i, int t) const
   { 
      switch (i) {
      case 0: 
         r[0] += t*l_[0];
         break;
      case 1: 
         r[1] += t*l_[1];
         r[2] += t*d_;
         break;
      case 2: 
         r[2] += t*l_[2];
         break;
      }
   }
   
   /* 
   * Shift generalized Vector r to primitive unit cell.
   */
   inline void MonoclinicBoundary::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 MonoclinicBoundary::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 minimum-image squared distance.
   */
   inline 
   double MonoclinicBoundary::distanceSq(const Vector &r1, const Vector &r2, 
                                         IntVector& translate) const
   {

      double dx = r1[0] - r2[0];
      if ( fabs(dx) > halfL_[0] ) {
         if (dx > 0.0) {
            dx -= l_[0];
            translate[0] = -1;
         } else {
            dx += l_[0];
            translate[0] = +1;
         }
         assert(fabs(dx) <= halfL_[0]);
      } else {
         translate[0] = 0;
      }

      double dy = r1[1] - r2[1];
      double dz = r1[2] - r2[2];
      double u2 = dz + c3_*dy;
      if ( fabs(dy) > halfL_[1] ) {
         if (dy > 0.0) {
            dy -= l_[1];
            dz -= d_;
            translate[1] = -1;
         } else {
            dy += l_[1];
            dz += d_;
            translate[1] = +1;
         }
         assert(fabs(dy) <= halfL_[1]);
      } else {
         translate[1] = 0;
      }

      if ( fabs(u2) >  halfL_[2]) {
         if (u2 > 0.0) {
            dz -= l_[2];
            translate[2] = -1;
         } else {
            dz += l_[2];
            translate[2] = +1;
         }
         assert(fabs(dz + c3_*dy) <= halfL_[2]);
      } else {
         translate[2] = 0;
      }
      return dx*dx+dy*dy+dz*dz;
   }

   /* 
   * Compute minimum image squared distance and separation.
   */
   inline 
   double MonoclinicBoundary::distanceSq(const Vector &r1, 
                                         const Vector &r2) const
   {

      double dx = r1[0] - r2[0];
      if ( fabs(dx) > halfL_[0] ) {
         if (dx > 0.0) {
            dx -= l_[0];
         } else {
            dx += l_[0];
         }
         assert(fabs(dx) <= halfL_[0]);
      }

      double dy = r1[1] - r2[1];
      double dz = r1[2] - r2[2];
      double u2 = dz + c3_*dy;
      if ( fabs(dy) > halfL_[1] ) {
         if (dy > 0.0) {
            dy -= l_[1];
            dz -= d_;
         } else {
            dy += l_[1];
            dz += d_;
         }
         assert(fabs(dy) <= halfL_[1]);
      }

      if ( fabs(u2) >  halfL_[2]) {
         if (u2 > 0.0) {
            dz -= l_[2];
         } else {
            dz += l_[2];
         }
         assert(fabs(dz+c3_*dy) <= halfL_[2]);
      }

      return dx*dx+dy*dy+dz*dz;
   }

   /* 
   * Calculate squared distance and separation vector by mininum convention.
   */
   inline 
   double MonoclinicBoundary::distanceSq(const Vector &r1, const Vector &r2, 
                                         Vector &dr) const
   {

      double dx = r1[0] - r2[0];
      if (fabs(dx) > halfL_[0]) {
         if (dx > 0.0) {
            dx -= l_[0];
         } else {
            dx += l_[0];
         }
         assert(fabs(dx) <= halfL_[0]);
      }

      double dy = r1[1] - r2[1];
      double dz = r1[2] - r2[2];
      double u2 = dz + c3_*dy;

      if (fabs(dy) > halfL_[1]) {
         if (dy > 0.0) {
            dy -= l_[1];
            dz -= d_;
         } else {
            dy += l_[1];
            dz += d_;
         }
         assert(fabs(dy) <= halfL_[1]);
      }

      if (fabs(u2) >  halfL_[2]) {
         if (u2 > 0.0) {
            dz -= l_[2];
         } else {
       dz += l_[2];
         }
         assert(fabs(dz+c3_*dy) <= halfL_[2]);
      }

      dr[0] = dx;
      dr[1] = dy;
      dr[2] = dz;

      return (dx*dx + dy*dy + dz*dz);
   }

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

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

   inline LatticeSystem MonoclinicBoundary::latticeSystem()
   { return lattice_; }

   inline 
   void MonoclinicBoundary::transformCartToGen(const Vector& Rc, Vector& Rg) 
   const
   {
      Rg[0] = Rc[0] * invL_[0];
      Rg[1] = Rc[1] * invL_[1];
      Rg[2] = (Rc[2] + c3_ *Rc[1]) * invL_[2];
   }
      
   inline 
   void MonoclinicBoundary::transformGenToCart(const Vector& Rg, Vector& Rc) 
   const
   {
      Rc[0] = Rg[0]*l_[0];
      Rc[1] = Rg[1]*l_[1];
      Rc[2] = Rg[2]*l_[2] + Rg[1]*d_;
   }

   /*
   * Serialize an OrthorhombicBoundary to/from an archive.
   */
   template <class Archive> void 
   MonoclinicBoundary::serialize(Archive& ar, const unsigned int version)
   {
      ar & l_;
      ar & d_;
      if (Archive::is_loading()) {
         reset();
         isValid();
      }
   }

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

namespace Util
{

   using namespace Simp;

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

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

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

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

}
#endif // ifdef  UTIL_MPI
#endif // ifndef UTIL_MONOCLINIC_BOUNDARY_H