.. _program_listing_file_src_lattice.cpp:

Program Listing for File lattice.cpp
====================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_lattice.cpp>` (``src/lattice.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   /* This file is part of brille.
   
   Copyright © 2019,2020 Greg Tucker <greg.tucker@stfc.ac.uk>
   
   brille is free software: you can redistribute it and/or modify it under the
   terms of the GNU Affero General Public License as published by the Free
   Software Foundation, either version 3 of the License, or (at your option)
   any later version.
   
   brille is distributed in the hope that it will be useful, but
   WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
   or FITNESS FOR A PARTICULAR PURPOSE.
   See the GNU Affero General Public License for more details.
   
   You should have received a copy of the GNU Affero General Public License
   along with brille. If not, see <https://www.gnu.org/licenses/>.            */
   #include "lattice.hpp"
   #include "hall_symbol.hpp"
   #include "approx.hpp"
   
   using namespace brille;
   
   Lattice::Lattice(const double* latmat, const int h){
     double l[3]={0,0,0}, a[3]={0,0,0};
     latmat_to_lenang(latmat,3,1,l,a); // (3,1) -> assume row-ordered matrix
     this->set_len_pointer(l,1);
     this->set_ang_pointer(a,1, AngleUnit::radian);
     this->volume=this->calculatevolume();
     this->check_hall_number(h);
   }
   Lattice::Lattice(const double* lengths, const double* angles, const int h, const AngleUnit au){
     this->set_len_pointer(lengths,1);
     this->set_ang_pointer(angles,1, au);
     this->volume=this->calculatevolume();
     this->check_hall_number(h);
   }
   Lattice::Lattice(const double la, const double lb, const double lc, const double al, const double bl, const double cl, const int h):
     len{{la,lb,lc}}, ang{{al,bl,cl}} {
     // this->set_len_scalars(la,lb,lc);
     // this->set_ang_scalars(al,bl,cl, AngleUnit::not_provided);
     this->check_ang(AngleUnit::not_provided);
     this->volume = this->calculatevolume();
     this->check_hall_number(h);
   }
   Lattice::Lattice(const double* latmat, const std::string& itname, const std::string& choice){
     double l[3]={0,0,0}, a[3]={0,0,0};
     latmat_to_lenang(latmat,3,1,l,a);
     this->set_len_pointer(l,1);
     this->set_ang_pointer(a,1, AngleUnit::radian);
     this->volume=this->calculatevolume();
     this->check_IT_name(itname, choice);
   }
   Lattice::Lattice(const double *lengths, const double *angles, const std::string& itname, const std::string& choice, const AngleUnit au){
     this->set_len_pointer(lengths,1);
     this->set_ang_pointer(angles,1, au);
     this->volume=this->calculatevolume();
     this->check_IT_name(itname, choice);
   }
   Lattice::Lattice(const double la, const double lb, const double lc, const double al, const double bl, const double cl, const std::string& itname, const std::string& choice):
   len({{la,lb,lc}}), ang({{al,bl,cl}}) {
     // this->set_len_scalars(la,lb,lc);
     // this->set_ang_scalars(al,bl,cl, AngleUnit::not_provided);
     this->check_ang(AngleUnit::not_provided);
     this->volume = this->calculatevolume();
     this->check_IT_name(itname, choice);
   }
   // void Lattice::set_len_scalars(const double a, const double b, const double c){
   //   this->len[0] = a;
   //   this->len[1] = b;
   //   this->len[2] = c;
   // }
   // void Lattice::set_ang_scalars(const double a, const double b, const double g, const AngleUnit angle_unit){
   //   this->ang[0] = a;
   //   this->ang[1] = b;
   //   this->ang[2] = g;
   //   this->check_ang(angle_unit);
   // }
   void Lattice::check_ang(const AngleUnit angle_unit){
     AngleUnit au = angle_unit;
     if (au == AngleUnit::not_provided){
       double minang = (std::numeric_limits<double>::max)();
       double maxang = (std::numeric_limits<double>::lowest)();
       for (int i=0; i<3; ++i){
         if (this->ang[i] < minang) minang = this->ang[i];
         if (this->ang[i] > maxang) maxang = this->ang[i];
       }
       if (minang < 0.) throw std::runtime_error("Unexpected negative inter-facial cell angle");
       au = (maxang < brille::pi) ? AngleUnit::radian : AngleUnit::degree;
     }
     if (au != AngleUnit::radian){
       double conversion = brille::pi/((AngleUnit::degree == au) ? 180.0 : 1.0);
       for (int i=0;i<3;i++) this->ang[i] *= conversion;
     }
   }
   void Lattice::check_hall_number(const int h){
     Spacegroup spg(h); // if h is invalid the next three lines might fail
     this->bravais = spg.get_bravais_type();
     this->spgsym = spg.get_spacegroup_symmetry();
     this->ptgsym = spg.get_pointgroup_symmetry();
   }
   void Lattice::check_IT_name(const std::string& itname, const std::string& choice){
     int hall_number = brille::string_to_hall_number(itname, choice);
     if (hall_number > 0 && hall_number < 531){
       this->check_hall_number(hall_number);
     } else {
       // maybe itname is actually a non-standard Hall symbol?
       HallSymbol hs(itname);
       Symmetry gens;
       if (hs.validate()){
         gens = hs.get_generators();
         bravais = hs.getl();
       } else {
         // last-ditch effort: maybe x,y,z notation Seitz matrices were passed?
         gens.from_ascii(itname);
         bravais = Bravais::P; // without any extra information this must be primitive
       }
       this->set_spacegroup_symmetry(gens);
     }
   }
   double Lattice::unitvolume() const{
     // The volume of a parallelpiped with unit length sides and our body angles
     double c[3];
     for (int i=0;i<3;++i) c[i]=std::cos(this->ang[i]);
     double sos=0, prd=2;
     for (int i=0;i<3;++i){
       sos+=c[i]*c[i];
       prd*=c[i];
     }
     return std::sqrt( 1 - sos + prd );
   }
   double Lattice::calculatevolume(){
     // we could replace this by l[0]*l[1]*l[2]*this->unitvolume()
     double tmp=1;
     double *a = this->ang.data();
     double *l = this->len.data();
     tmp *= sin(( a[0] +a[1] +a[2])/2.0);
     tmp *= sin((-a[0] +a[1] +a[2])/2.0);
     tmp *= sin(( a[0] -a[1] +a[2])/2.0);
     tmp *= sin(( a[0] +a[1] -a[2])/2.0);
     tmp = sqrt(tmp);
     tmp *= 2*l[0]*l[1]*l[2];
     this->volume = tmp;
     if (std::isnan(tmp)){
       std::string msg = "Invalid lattice unit cell ";
       if (std::isnan(this->unitvolume())){
         msg += "angles [";
         for (int i=0; i<3; ++i) msg += " " + std::to_string(a[0]/brille::pi*180.);
         msg += " ]";
       } else {
         msg += "lengths [";
         for (int i=0; i<3; ++i) msg += " " + std::to_string(l[i]);
         msg += " ]";
       }
       throw std::domain_error(msg);
     }
     return tmp;
   }
   Lattice Lattice::inner_star() const {
     std::array<double,3> lstar, astar, cosang, sinang;
     for (size_t i=0; i<3u; ++i){
       cosang[i] = std::cos(this->ang[i]);
       sinang[i] = std::sin(this->ang[i]);
     }
     for (size_t i=0; i<3u; ++i){
       size_t j = (i+1)%3;
       size_t k = (i+2)%3;
       lstar[i] = 2*brille::pi*sinang[i]*this->len[j]*this->len[k]/this->volume;
       astar[i] = std::acos((cosang[j]*cosang[k]-cosang[i])/(sinang[j]*sinang[k]));
     }
     return Lattice(lstar, astar, this->bravais, this->spgsym, this->ptgsym, this->basis);
   }
   void Lattice::get_metric_tensor(double * mt) const {
     const double *a = this->ang.data();
     const double *l = this->len.data();
   
     double cosa, cosb, cosc;
     cosa = cos(a[0]);
     cosb = cos(a[1]);
     cosc = cos(a[2]);
   
     mt[0] = l[0]*l[0];
     mt[1] = l[1]*l[0]*cosc;
     mt[2] = l[2]*l[0]*cosb;
   
     mt[3] = l[0]*l[1]*cosc;
     mt[4] = l[1]*l[1];
     mt[5] = l[2]*l[1]*cosa;
   
     mt[6] = l[0]*l[2]*cosb;
     mt[7] = l[1]*l[2]*cosa;
     mt[8] = l[2]*l[2];
   }
   void Lattice::get_covariant_metric_tensor(double *mt) const {
     this->get_metric_tensor(mt);
   }
   void Lattice::get_contravariant_metric_tensor(double *mt) const {
     double *tmp = nullptr;
     tmp = new double[9]();
     this->get_metric_tensor(tmp);
     brille::utils::matrix_inverse(mt,tmp);
     delete[] tmp;
   }
   
   std::vector<double> Lattice::get_metric_tensor() const {
     std::vector<double> m(9);
     this->get_metric_tensor(m.data());
     return m;
   }
   std::vector<double> Lattice::get_covariant_metric_tensor() const {
     std::vector<double> m(9);
     this->get_covariant_metric_tensor(m.data());
     return m;
   }
   std::vector<double> Lattice::get_contravariant_metric_tensor() const {
     std::vector<double> m(9);
     this->get_contravariant_metric_tensor(m.data());
     return m;
   }
   
   bool Lattice::issame(const Lattice& lat) const{
     return brille::approx::vector(3, this->ang.data(), lat.ang.data()) && brille::approx::vector(3, this->len.data(), lat.len.data());
   }
   
   bool Lattice::isapprox(const Lattice& lat) const {
     return this->ispermutation(lat)==0 ? false : true;
   }
   int Lattice::ispermutation(const Lattice& lat) const {
     double a[3], l[3];
     int i, j, ap[3]{0,2,1}; // for anti-permutations
     for (j=0; j<3; ++j){
       for (i=0; i<3; ++i){
         a[i] = this->ang[(i+j)%3];
         l[i] = this->len[(i+j)%3];
       }
       if (brille::approx::vector(3, a, lat.ang.data()) && brille::approx::vector(3, l, lat.len.data()))
         return j+1;
     }
     for (j=0; j<3; ++j){
       for (i=0; i<3; ++i){
         a[i] = this->ang[ap[(i+j)%3]];
         l[i] = this->len[ap[(i+j)%3]];
       }
       if (brille::approx::vector(3, a, lat.ang.data()) && brille::approx::vector(3, l, lat.len.data()))
         return -j-1;
     }
     return 0;
   }
   
   
   void Lattice::print(){
     printf("(%g %g %g) " ,this->len[0], this->len[1], this->len[2]);
     printf("(%g %g %g)\n",this->ang[0]/brille::pi*180, this->ang[1]/brille::pi*180, this->ang[2]/brille::pi*180);
   }
   std::string lattice2string(const Lattice& l, const std::string& lenunit="", const std::string& angunit="°"){
     std::string repr;
     repr = "(" + std::to_string(l.get_a()) + " "
                + std::to_string(l.get_b()) + " "
                + std::to_string(l.get_c()) + ")" + lenunit
          +" (" + std::to_string(l.get_alpha()/brille::pi*180) + " "
                + std::to_string(l.get_beta()/brille::pi*180)  + " "
                + std::to_string(l.get_gamma()/brille::pi*180) + ")" + angunit;
     return repr;
   }
   std::string Lattice::string_repr(){ return lattice2string(*this); }
   std::string Direct::string_repr(){return lattice2string(*this,"Å");}
   std::string Reciprocal::string_repr(){return lattice2string(*this,"Å⁻¹");}
   
   Reciprocal Direct::star() const {
     return Reciprocal(this->inner_star());
   }
   Direct Reciprocal::star() const {
     return Direct(this->inner_star());
   }
   
   void Direct::get_lattice_matrix(double *latmat) const{
     this->get_lattice_matrix(latmat, 3u, 1u);
   }
   template<class I> void Direct::get_lattice_matrix(double *latmat, std::vector<I>& strides) const {
     this->get_lattice_matrix(latmat, static_cast<size_t>(strides[0]/sizeof(double)), static_cast<size_t>(strides[1]/sizeof(double)));
   }
   void Direct::get_lattice_matrix(double *latmat, const size_t c, const size_t r) const{
     // define the lattice basis vectors using the same convention as spglib
     double c0=std::cos(this->ang[0]), c1=std::cos(this->ang[1]), c2=std::cos(this->ang[2]);
     double s2=std::sin(this->ang[2]);
   
     double xhat[3]={1,0,0}, yhat[3]={c2,s2,0};
     double zhat[3]={c1, (c0-c2*c1)/s2, this->unitvolume()/s2};
   
     for (int i=0;i<3;++i){
       // latmat[i*c+0*r] = this->len[0]*xhat[i];
       // latmat[i*c+1*r] = this->len[1]*yhat[i];
       // latmat[i*c+2*r] = this->len[2]*zhat[i];
       latmat[0*c+i*r] = this->len[0]*xhat[i]; // ̂x direction, first row
       latmat[1*c+i*r] = this->len[1]*yhat[i]; // ̂y direction, second row
       latmat[2*c+i*r] = this->len[2]*zhat[i]; // ̂z direction, third row
     }
   }
   void Direct::get_xyz_transform(double* toxyz) const{
     this->get_xyz_transform(toxyz, 3u, 1u);
   }
   template<class I> void Direct::get_xyz_transform(double* x, std::vector<I>& s) const {
     this->get_xyz_transform(x, static_cast<size_t>(s[0]/sizeof(double)), static_cast<size_t>(s[1]/sizeof(double)));
   }
   void Direct::get_xyz_transform(double *toxyz, const size_t c, const size_t r) const {
     // there are infinite possibilities for your choice of axes.
     // the original spglib used x along a and y in the (a,b) plane
     // here we're going with x along astar and y in the (astar, bstar) plane -- this is the "B-matrix"
     double B[9], t[9];
     this->star().get_B_matrix(B, 3u, 1u);
     // we want toxyz to be 2*pi*transpose(inverse(B));
     // use t as a buffer
     brille::utils::matrix_inverse(t,B);
     brille::utils::matrix_transpose(t); // in-place transpose using swap
     for (int i=0; i<3; ++i) for (int j=0; j<3; ++j) toxyz[i*c+j*r] = 2*brille::pi*t[3*i+j];
   }
   std::vector<double> Direct::get_xyz_transform() const {
     std::vector<double> mat(9);
     this->get_xyz_transform(mat.data());
     return mat;
   }
   void Direct::get_inverse_xyz_transform(double* fromxyz) const{
     this->get_inverse_xyz_transform(fromxyz, 3u, 1u);
   }
   template<class I> void Direct::get_inverse_xyz_transform(double* x, std::vector<I>& s) const {
     this->get_inverse_xyz_transform(x, static_cast<size_t>(s[0]/sizeof(double)), static_cast<size_t>(s[1]/sizeof(double)));
   }
   void Direct::get_inverse_xyz_transform(double *fromxyz, const size_t c, const size_t r) const{
     double toxyz[9], tmp[9];
     this->get_xyz_transform(toxyz,3u,1u);
     if(!brille::utils::matrix_inverse(tmp, toxyz))
       throw std::runtime_error("xyz_transform matrix has zero determinant");
     for (size_t i=0; i<3u; ++i) for (size_t j=0; j<3u; ++j) fromxyz[i*c+j*r] = tmp[i*3+j];
   }
   std::vector<double> Direct::get_inverse_xyz_transform() const {
     std::vector<double> mat(9);
     this->get_inverse_xyz_transform(mat.data());
     return mat;
   }
   void Reciprocal::get_lattice_matrix(double *latmat) const{
     this->get_lattice_matrix(latmat, 3u, 1u);
   }
   template<class I> void Reciprocal::get_lattice_matrix(double* x, std::vector<I>& strides) const {
     this->get_lattice_matrix(x, static_cast<size_t>(strides[0]/sizeof(double)), static_cast<size_t>(strides[1]/sizeof(double)));
   }
   void Reciprocal::get_lattice_matrix(double *latmat, const size_t c, const size_t r) const{
     Direct d = this->star();
     double m[9], tmp[9];
     d.get_lattice_matrix(m, 3u, 1u);
     brille::utils::matrix_inverse(tmp,m);
     brille::utils::matrix_transpose(tmp);
     for (int i=0; i<9; ++i) tmp[i]*=2*brille::pi;
     for (size_t i=0; i<3u; ++i) for (size_t j=0; j<3u; ++j) latmat[c*i+r*j] = tmp[3*i+j];
   }
   void Reciprocal::get_B_matrix(double* B) const {
     this->get_B_matrix(B, 3u, 1u);
   }
   template<class I> void Reciprocal::get_B_matrix(double* B, std::vector<I>& strides) const {
     this->get_B_matrix(B, static_cast<size_t>(strides[0]/sizeof(double)), static_cast<size_t>(strides[1]/sizeof(double)));
   }
   void Reciprocal::get_B_matrix(double *B, const size_t c, const size_t r) const {
     //Calculate the B-matrix as in Acta Cryst. (1967). 22, 457
     // http://dx.doi.org/10.1107/S0365110X67000970
     Direct d = this->star();
   
     double asb, asg;
     asb = sin(this->get_beta());
     if (asb<0) asb*=-1.0;
     asg = sin(this->get_gamma());
     if (asg<0) asg*=-1.0;
     // Be careful about indexing B. Make sure each vector goes in as a column, not a row!
     // if you mix this up, you'll only notice for non-orthogonal space groups
   
     // a-star along x -- first column (constant row index = 0)
     B[0*c+0*r] = this->get_a();
     B[1*c+0*r] = 0.0;
     B[2*c+0*r] = 0.0;
     // b-star in the x-y plane -- second column
     B[0*c+1*r] = this->get_b()*cos(this->get_gamma());
     B[1*c+1*r] = this->get_b()*asg;
     B[2*c+1*r] = 0.0;
     // and c-star -- third column
     B[0*c+2*r] = this->get_c()*cos(this->get_beta());
     B[1*c+2*r] = -1.0*(this->get_c())*asb*cos(d.get_alpha());
     B[2*c+2*r] = 2*brille::pi/d.get_c();
   }
   void Reciprocal::get_xyz_transform(double *toxyz) const { this->get_xyz_transform(toxyz, 3u, 1u); }
   template<class I> void Reciprocal::get_xyz_transform(double* toxyz, std::vector<I>& strides) const {
     this->get_xyz_transform(toxyz, static_cast<size_t>(strides[0]/sizeof(double)), static_cast<size_t>(strides[1]/sizeof(double)));
   }
   void Reciprocal::get_xyz_transform(double *toxyz, const size_t c, const size_t r) const {
     // there are infinite possibilities for your choice of axes.
     // the original spglib used x along a and y in the (a,b) plane
     // here we're going with x along astar and y in the (astar, bstar) plane -- this is the "B-matrix"
     this->get_B_matrix(toxyz, c, r);
   }
   std::vector<double> Reciprocal::get_xyz_transform() const {
     std::vector<double> mat(9);
     this->get_xyz_transform(mat.data());
     return mat;
   }
   void Reciprocal::get_inverse_xyz_transform(double *fromxyz) const{
     this->get_inverse_xyz_transform(fromxyz, 3u, 1u);
   }
   template<class I> void Reciprocal::get_inverse_xyz_transform(double* x, std::vector<I>& strides) const {
     this->get_inverse_xyz_transform(x, static_cast<size_t>(strides[0]/sizeof(double)), static_cast<size_t>(strides[1]/sizeof(double)));
   }
   void Reciprocal::get_inverse_xyz_transform(double *fromxyz, const size_t c, const size_t r) const {
     double toxyz[9], tmp[9];
     this->get_xyz_transform(toxyz, 3u, 1u);
     if(!brille::utils::matrix_inverse(tmp, toxyz))
       throw std::runtime_error("xyz_transform matrix has zero determinant");
     for (size_t i=0; i<3u; ++i) for(size_t j=0; j<3u; ++j) fromxyz[i*c+j*r] = tmp[i*3+j];
   }
   std::vector<double> Reciprocal::get_inverse_xyz_transform() const {
     std::vector<double> mat(9);
     this->get_inverse_xyz_transform(mat.data());
     return mat;
   }
   
   // We have to define these separately from Lattice since Lattice.star() doesn't exist.
   bool Direct::isstar(const Reciprocal& latt) const{
     // two lattices are the star of each other if the star of one is the same as the other
     // we need to check both ways in case rounding errors in the inversion cause a problem
     return ( this->issame( latt.star() ) || latt.issame( this->star() ) );
   }
   bool Reciprocal::isstar(const Direct& latt) const{
     // two lattices are the star of each other if the star of one is the same as the other
     // we need to check both ways in case rounding errors in the inversion cause a problem
     return ( this->issame( latt.star() ) || latt.issame( this->star() ) );
   }
   
   //bool Direct::issame(const Reciprocal) const {printf("call to Direct::issame(Reciprocal)\n"); return false;}
   bool Direct::isstar(const Direct&) const {return false;}
   //bool Reciprocal::issame(const Direct) const {printf("call to Reciprocal::issame(Direct)\n"); return false;}
   bool Reciprocal::isstar(const Reciprocal&) const {return false;}
   
   void Direct::print(){
     printf("(%g %g %g)A " ,this->len[0], this->len[1], this->len[2]);
     printf("(%g %g %g)\n",this->ang[0]/brille::pi*180, this->ang[1]/brille::pi*180, this->ang[2]/brille::pi*180);
   }
   void Reciprocal::print(){
     printf("(%g %g %g)/A " ,this->len[0], this->len[1], this->len[2]);
     printf("(%g %g %g)\n",this->ang[0]/brille::pi*180, this->ang[1]/brille::pi*180, this->ang[2]/brille::pi*180);
   }
   
   
   Direct Direct::primitive(void) const{
     PrimitiveTransform P(this->bravais);
     if (P.does_anything()){
       double plm[9], lm[9];
       this->get_lattice_matrix(lm); // now returns *row* vectors!
       // The transformation matrix P gives us the primitive basis column-vector
       // matrix Aₚ from the standard basis column-vector matrix Aₛ by
       // Aₚ = Aₛ P. But here we have Aₛᵀ and want Aₚᵀ:
       // Aₚᵀ = (Aₛ P)ᵀ = Pᵀ Aₛᵀ
       // So we need the transpose of P from the PrimitiveTransform object:
       brille::utils::multiply_matrix_matrix<double,double,double,3>(plm,P.get_Pt().data(),lm);
       return Direct(plm);
     }
     return *this;
   }
   Reciprocal Reciprocal::primitive(void) const{
     // We could try to use the fact that for column vectors Bₚ = Bₛ (P⁻¹)ᵀ
     // but it's probably safer to continue using the reciprocal of the primitive
     // of the reciprocal of this lattice.
     return this->star().primitive().star();
   }