hkl_symmetry.hh

// \file hkl_symmetry.hh

#ifndef HKL_SYMMETRY_HEADER
#define HKL_SYMMETRY_HEADER

#include <vector>
// Clipper
#include <clipper/clipper.h>
using clipper::Message;
using clipper::Message_fatal;
using clipper::Message_info;

#include "hkl_datatypes.hh"
#include "lattice.hh"

namespace scala {
  //--------------------------------------------------------------
  std::string SGnameHtoR(const std::string& sgname, const char& HorR = ' ');
  //--------------------------------------------------------------

  class SpaceGroup : public clipper::Spacegroup
                     // Clipper Spacegroup with symops in defined order
                     // ie same as in constructor
  {
  public:
    // Constructors
    SpaceGroup() : Nsymp(0) {}
    SpaceGroup(std::vector<clipper::Symop>& symops); 
    SpaceGroup(const std::string& spgname);          
    SpaceGroup(const int& SpgNumber);                
    SpaceGroup(const clipper::Spacegroup& ClpSG);    

    void init(std::vector<clipper::Symop>& symops);
    void init(const std::string& spgname);         
    void init(const int& SpgNumber);               
    void init(const clipper::Spacegroup& ClpSG);   

    bool Null() const {return (Nsymp == 0);}

    clipper::Symop Symop(const int& symN) const
    {return csymops.at(symN);}  
    clipper::Symop RotSymop(const int& symN) const
    {return rotsymops.at(symN);} 
    clipper::Symop InvRotSymop(const int& symN) const
    {return invrotsymops.at(symN);}  

    clipper::Symop_code RotSymopCode(const int& symN) const
    {return clipper::Symop_code(rotsymops.at(symN));}
    clipper::Symop_code InvRotSymopCode(const int& symN) const
    {return clipper::Symop_code(invrotsymops.at(symN));}

    SpaceGroup PattersonGroup() const;
    SpaceGroup PointGroup() const;

    char LatType() const {return lattype;}
    SpaceGroup NewLatticePointGroup(const char& Lattype) const;

    std::vector<clipper::Symop_code> SymopCodes() const;

    bool IsSymopIdentity(const int& symN) const; 

    std::string formatAllSymops_as_xyz() const; 
    std::string formatAllSymops_as_hkl() const; 

    std::string Symbol_hm() const {return spacegroupname;}
    std::string Symbol_hall() const;
    int Spacegroup_number() const {return spacegroupnumber;}
    int CCP4_Spacegroup_number() const {return CCP4spacegroupnumber;}

    clipper::HKL put_in_asu(const clipper::HKL& hkl, int& isym) const;
    clipper::HKL get_from_asu(const clipper::HKL& hkl, const int& isym) const;

    friend bool operator == (const SpaceGroup& a,const SpaceGroup& b);
    friend bool operator != (const SpaceGroup& a,const SpaceGroup& b);

    void ChangeBasis(const scala::ReindexOp& reindex);
    
  private:
    int Nsymp;
    char lattype;  // 'H' for hexagonal-setting rhombohedral
    std::string spacegroupname;
    int spacegroupnumber;
    int CCP4spacegroupnumber;

    // The following operators are in constructor-order
    // (unlike clipper::Spacegroup which may reorder them
    // Note that if [Ri] is the rotation part of the i'th symmetry operator,
    // then in reciprocal space we have indices actual h(i) and reduced h(red)
    //   h(red) = h(i)T   [Ri]
    //   h(i)   = h(red)T [Ri]^-1 
    std::vector<clipper::Symop> csymops;      // full space group operators
    std::vector<clipper::Symop> rotsymops;    // rotation operators [Ri]
    std::vector<clipper::Symop> invrotsymops; // inverse rotation operators [Ri]^-1

    // Private functions
    void init();
    // Return lattice type extracted from space group name & symops
    void SetLatType();
    // CCP4 symop from clipper symop
    CSym::ccp4_symop MakeCCP4symop(const clipper::Symop& S) const;
    // Use CCP4 library routine to get name from operators
    std::string CCP4spaceGroupName() const;
    // Use CCP4 library routine to get numbers from operators
    void CCP4spaceGroupNumber();
    // Clipper symop from CCP4 symop
    clipper::Symop MakeClippersymop(const CSym::ccp4_symop& S) const;
    // Make symop string from CCP4 space group
    std::string AllSymopsfromCCP4(const CSym::CCP4SPG* ccp4sg) const;
  };  // end class SpaceGroup
  //--------------------------------------------------------------
  class SymElement {
  public:
    int Nfold;  // order of rotation
    clipper::Vec3<int> iaxis;  // rotation axis direction vector
    std::vector<int> symops;  // list of symop numbers representing this element
    SymElement() : Nfold(0) {};  // null constructor
    void print() const;
  };
  //--------------------------------------------------------------
  class hkl_symmetry
  //
  //
  //  Constructor
  //    hkl_symmetry(const scala::SpaceGroup& SG)
  //    hkl_symmetry(const std::string SpgName);
  //    hkl_symmetry(const int& SpgNumber);
  //    hkl_symmetry(const clipper::Spacegroup& ClpSG);
  //
  //  Reset symmetry
  //    void set_symmetry()
  //
  //  Constructor & set_symmetry setup tables to identify symmetry
  //  elements (N-fold rotations) and list corresponding symmetry
  //  operators
  //  Also sets up the operation corresponding to pairs of ISYM codes
  //  ( ISYM code = SymopNumber*2 + 1(hkl) or +2(-h-k-l))
  //
  //  An n-fold rotation axis corresponds to n operators including I
  //  I is a 1-fold rotation
  //
  //  int get_symelmt(isym1, isym2) returns symmetry element number
  //                               (from 0) of primitive operator
  //                                relating a pair ofobservations
  //                                with isym1 & isym2
  //
  {
  public:
    // Constructors
    hkl_symmetry();
    hkl_symmetry(const clipper::Spacegroup& cspgp); 
    hkl_symmetry(const scala::SpaceGroup& cspgp); 
    hkl_symmetry(const std::string SpgName); 
    hkl_symmetry(const int& SpgNumber); 

    bool IsNull() const {return (Nsymp == 0);}
    void set_reindex(const ReindexOp& op) {ReindexToZeroFrame = op;}

    SpaceGroup GetSpaceGroup() const {return spaceGroup;}
    clipper::Spacegroup ClipperGroup() const {return spaceGroup;}

    Hkl put_in_asu(const Hkl& hkl, int& isym) const; 
    clipper::HKL put_in_asu(const clipper::HKL& hkl, int& isym) const;
    Hkl get_from_asu(const Hkl& hkl, const int& isym) const;
    bool is_centric(const Hkl& hkl) const; 
    float epsilon(const Hkl& hkl) const;  

    bool IsCentro() const {return centro;}  
    bool IsChiral() const {return chiral;}  

    bool LatticePresent(const Hkl& hkl) const; 
    void ChangeBasis(const scala::ReindexOp& reindex);

    // *****
    // ***  Symmetry element stuff
    int Nelement() const {return Nelement_;}
    bool IsElementIdent(const int& kelement) const;
    clipper::Vec3<int>  AxisDirection(const SymElement& elmt) const;
    std::string format_element(const int& kelement) const;
    std::string XML_element(const int& kelement) const;
    void print_element(const int& idx) const; 
    void print_elements() const; 


    int NopInElement(const int& kelement) const;
    std::vector<double> SymopInElement(const int& j,
                                       const int& kelement) const;

    int get_symelmt(const int& isym1, const int& isym2)  const;


    std::vector<int> get_symelmt(const Hkl& h1, const Hkl& h2) const;

    // Return spacegroup names etc
    std::string symbol_Hall() const;        
    std::string symbol_xHM(const char& HorR=' ') const;  
    std::string symbol() const {return spgname;} 

    char lattice_type() const {return LatType;}
    int NsymP() const {return Nsymp;}
    int Nsym() const {return spaceGroup.num_symops();}

    std::vector<clipper::Symop> RotSymopsInElement(const int& kelement) const;
    std::vector<clipper::Symop> PrimRotSymops() const;


    CrystalSystem CrysSys() const {return cryssys;}

    std::vector<int> CellConstraint() const;


    bool equals_rt(const hkl_symmetry& other) const; 
    bool equals_r(const hkl_symmetry& other) const; 

    bool operator==(const hkl_symmetry& other) const {return equals_r(other);}
    bool operator!=(const hkl_symmetry& other) const {return !equals_r(other);}

    bool equalElement(const int& jel, const hkl_symmetry& other, const int& jel2) const;

  private:
    int Nsymp;
    SpaceGroup spaceGroup; //  the actual space group
    std::string spgname;  // name from constructor or xHM name
    char LatType;
    int Nelement_;
    CrystalSystem cryssys;
    bool centro;
    bool chiral;

    std::vector<clipper::Symop_code> inv_rot_symcodes;
    std::vector<clipper::Symop_code> rot_symcodes;
    std::vector<SymElement> elements; // symops for each sub-set
    std::vector<int> element_index;  // element number for this symop (Nsymp)
    std::vector<int> symelmt;   // a 2D array Nelement x Nelement

    // This is the reciprocal space operator to convert back to
    // initial frame
    // Reindex operator H 
    // Applies to index h such that h'T = hT H
    ReindexOp ReindexToZeroFrame; // initially identity

    // private functions
    bool ElementEqual(const SymElement& a, const SymElement& b) const;
    void set_symmetry(); // set up symmetry stuff on initialisation
    // Get crystal system
    CrystalSystem CrysSysfromSGnumber(const int& SpaceGroupNumber);
    // return true if space group is chiral
    bool ChiralTest();
  }; 
  //--------------------------------------------------------------
}
#endif

---