GbShifts.cpp

Go to the documentation of this file.

//
// Created by Nikhil Chandra Admal on 2/4/24.
//
#include "../../include/Lattices/GbShifts.h"
#include "../../include/Utilities/randomInteger.h"
 
namespace oILAB {
template <int dim>
GbShifts<dim>::GbShifts(const Gb<dim> &gb,
                        const ReciprocalLatticeVector<dim> &axis,
                        const std::vector<LatticeVector<dim>> &gbCslVectors,
                        const double &bhalfMax)
    : gb(gb), axis(axis), gbCslVectors(gbCslVectors),
      bShiftPairs(getbShiftPairs(gb, gbCslVectors, bhalfMax)) {
  std::cout << "--------------------GBShifts class construction "
               "---------------------------"
            << std::endl;
  std::cout << "GB CSL vectors = " << std::endl;
  for (const auto &elem : gbCslVectors)
    std::cout << elem.cartesian().transpose() << std::endl;
  std::cout << std::endl;
 
  std::cout << "GB reciprocal CSL vectors = " << std::endl;
  Eigen::Matrix<double, dim, dim - 1> gbCslBasis;
  for (int i = 0; i < dim - 1; ++i)
    gbCslBasis.col(i) = gbCslVectors[i].cartesian();
  Eigen::Matrix<double, dim, dim - 1> gbCslReciprocalBasis =
      gbCslBasis.completeOrthogonalDecomposition().pseudoInverse().transpose();
  std::cout << gbCslReciprocalBasis.transpose() << std::endl;
  std::cout << std::endl;
 
  std::cout << "Maximum b < "
            << 2 * bhalfMax * gb.bc.A.latticeBasis.col(0).norm() << std::endl;
  std::cout << std::endl;
 
  VectorDimD normal;
  if (dim == 3)
    normal = gbCslVectors[0].cross(gbCslVectors[1]).cartesian().normalized();
  else
    normal = gbCslVectors[0].cross().cartesian().normalized();
 
  std::cout << "Exploring the following translation-shift pairs:" << std::endl;
  for (const auto &[b, s] : bShiftPairs) {
    std::cout << "b = " << b.cartesian().transpose();
    std::cout << "; s = " << s.transpose() << std::endl;
    // if (abs(s.dot(normal)) > FLT_EPSILON)
    if (abs(s.dot(normal)) > 1e-6)
      throw std::runtime_error(
          "GBShifts construction failed - shifts are not parallel to the GB.");
 
    Eigen::MatrixXd shiftCoordinates = gbCslReciprocalBasis.transpose() * s;
    if ((shiftCoordinates.array() < -FLT_EPSILON).any() ||
        (shiftCoordinates.array() > 1 + FLT_EPSILON).any()) {
      std::cout << "Shift coordinates = " << shiftCoordinates.transpose()
                << std::endl;
      throw std::runtime_error(
          "GB shifts are not in the area spanned by the GB CSL vectors.");
    }
  }
  std::cout << "----------------------------" << std::endl;
  std::cout << std::endl;
 
    }
 
    template<int dim>
    std::vector<std::pair<LatticeVector<dim>, typename GbShifts<dim>::VectorDimD>> GbShifts<dim>::getbShiftPairs(const Gb<dim>& gb,
                                                                                                                 const std::vector<LatticeVector<dim>>& gbCslVectors,
                                                                                                                 const double& bhalfMax)
    {
        std::vector<std::pair<LatticeVector<dim>, VectorDimD>> output;
 
        assert(gbCslVectors.size()==dim-1);
        auto nC= gb.bc.getReciprocalLatticeDirectionInC(gb.nB.reciprocalLatticeVector());
 
        // form a CSL cell for modulo operations
        auto gbPlaneParallelCslBasis= gb.bc.csl.planeParallelLatticeBasis(nC,true);
        std::vector<LatticeVector<dim>> cslSubLatticeVectors;
        cslSubLatticeVectors.push_back(gbPlaneParallelCslBasis[0].latticeVector());
        cslSubLatticeVectors.push_back(gbCslVectors[0]);
        cslSubLatticeVectors.push_back(gbCslVectors[1]);
        auto cslPoints= gb.bc.csl.box(cslSubLatticeVectors);
 
        // form a T lattice cell for exploring translations
        auto nT= gb.getReciprocalLatticeDirectionInT(gb.bc.getReciprocalLatticeDirectionInD(gb.nA.reciprocalLatticeVector()).reciprocalLatticeVector());
        auto planeParallelBasisT= gb.T.planeParallelLatticeBasis(nT,true);
        std::vector<LatticeVector<dim>> latticeVectorsT;
 
        double latticeConstant= gb.bc.A.latticeBasis.col(0).norm();
        for(int i=0; i<dim; ++i)
        {
            // scale the lattice vectors of T based on the bhalfMax parameter
            //int factor= floor(bhalfMax*latticeConstant/planeParallelBasisT[i].latticeVector().cartesian().norm()+FLT_EPSILON);
            int factor= floor(35*bhalfMax*latticeConstant/planeParallelBasisT[i].latticeVector().cartesian().norm()+FLT_EPSILON);
            factor= (factor>0 ? factor : 1);
            latticeVectorsT.push_back(factor * planeParallelBasisT[i].latticeVector());
        }
        auto points= gb.T.box(latticeVectorsT,"T.txt");
 
 
        VectorDimD shiftT, shiftC;
        shiftT << -0.5, -0.5, -0.5;
        shiftC << -0.5, -FLT_EPSILON, -FLT_EPSILON;
        for(auto& point : points) {
            //if (point.cartesian().norm() > bhalfMax*gb.bc.A.latticeBasis.col(0).norm())
            //    continue;
            LatticeVector<dim>::modulo(point, latticeVectorsT, shiftT);
            if (point.cartesian().norm() > bhalfMax*gb.bc.A.latticeBasis.col(0).norm())
                continue;
            auto cslShift = LatticeVector<dim>((gb.bc.LambdaA * gb.basisT * point).eval(), gb.bc.dscl);
            for(const auto& cslPoint : cslPoints) {
                // only for those CSL points on the boundary
                if (gb.bc.getLatticeVectorInA(cslPoint).dot(gb.nA)!=0) continue;
                VectorDimD cslShiftCentered = cslShift.cartesian() + cslPoint.cartesian() - point.cartesian() / 2;
                LatticeVector<dim>::modulo(cslShiftCentered, cslSubLatticeVectors, shiftC);
                output.push_back(std::make_pair(point, cslShiftCentered));
            }
        }
        return output;
    }
    /*
    template<int dim>
    std::vector<LatticeVector<dim>> GbShifts<dim>::getGbCslVectors(const Gb<dim>& gb, const ReciprocalLatticeVector<dim>& axis)
    {
        std::vector<LatticeVector<dim>> output;
        LatticeVector<dim> axisA(gb.bc.A.latticeDirection(axis.cartesian()).latticeVector());
        LatticeVector<dim> axisC(gb.bc.getLatticeDirectionInC(axisA).latticeVector());
        for(int i=0; i<dim-1; ++i) {
            if (i==0) output.push_back(gb.getPeriodVector(axis));
            if (i==1) output.push_back(axisC);
        }
        return output;
    }
     */
 
    //template class GbShifts<2>;
    template class GbShifts<3>;
 
    } // namespace oILAB