Program Listing for File symmetry_group.hpp

↰ Return to documentation for file (src_ed/symmetry_group.hpp)

#ifndef symmetry_group_h
#define symmetry_group_h

#include <iostream>
#include <vector>
#include <typeinfo>
#include <cstdint>

#include "qcm_ED.hpp"
#include "matrix.hpp"
#include "binary_state.hpp"
#include "sector.hpp"
#include "parser.hpp"
#include "index.hpp"
#include "symmetric_orbital.hpp"
#include "block_matrix.hpp"
#include "index_pair.hpp"

using namespace std;

#define MAX_GENERATORS 4

struct rep_map
{
  binary_state b;
  int phase;
  int length;
};


struct symmetry_group
{
  bool has_complex_irrep;
  bool bath_irrep;
  int N;
    int n_sites;
    int g;
  int L;
    vector<vector<int>> generator;
    vector<vector<int>> e;
    matrix<int> phi;
    matrix<Complex> chi;
    matrix<Complex> S;
  matrix<double> S_real;
  vector<bool> complex_irrep;
  vector<int> site_irrep_dim;
    vector<int> site_irrep_dim_cumul;
    matrix<int> tensor;
    vector<int> conjugate;
  vector<int> phaseR;
  vector<Complex> phaseC;


  symmetry_group(int _N, int _n_sites, const vector<vector<int>> &gen, bool bath_irrep);
  bool is_valid_element(const vector<int>& v, bool full=false);
  bool is_identity(const vector<int>& v, bool full=false);
  vector<int> product(const vector<int>& x, const vector<int>& y);
  vector<int> inverse(const vector<int>& x);
  vector<int> identity();
  int sign(const vector<int>& x);
  uint32_t apply_single(const vector<int>& x, const uint32_t b);
  pair<uint32_t, int> apply(const vector<int>& x, const uint32_t b);
  template<typename T>
  vector<T> apply_to_vector(const vector<int>& v, const vector<T>& x);
  pair<int, int> map_orbital(const vector<int>& v, int a);
  bool diff(const vector<int>& x, const vector<int>& y, bool full=false);
  void build();
  pair<binary_state, int> apply(const int &elem, const binary_state &b);
  pair<binary_state, int> apply(const int &elem, const int &irrep, const binary_state &b);
  rep_map Representative(const binary_state &b, int irrep);
  void to_site_basis(int r, vector<Complex> &x, vector<Complex> &y, int n_mixed);
  void to_site_basis(int r, vector<double> &x, vector<double> &y, int n_mixed);
  void to_site_basis(block_matrix<Complex> &B, matrix<Complex> &G, int n_mixed);
  bool sector_is_valid(const sector &sec);
  sector shift_sector(const sector &sec, int pm, int spin, int _irrep);
  bool operator!=(const symmetry_group &y);

  template<typename HilbertField>
  inline HilbertField phaseX(int i);

  template<typename HilbertField>
  void check_invariance(const map<index_pair,HilbertField> &elements, const string& name, bool interaction);
  vector<vector<symmetric_orbital>> build_symmetric_orbitals(int mixing);
};


std::ostream & operator<<(std::ostream &s, const symmetry_group &x);


// implementation


template<typename HilbertField>
void symmetry_group::check_invariance(const map<index_pair,HilbertField> &elements, const string& name, bool interaction)
{
  if(g==1) return;

  for(int i=0; i<generator.size(); i++){
    for(auto &x : elements){
      auto A = map_orbital(generator[i], x.first.r);
      A.second *= -1; // sign because c_a^\dg is a creation operator
      auto B = map_orbital(generator[i], x.first.c);
      int f = (A.second + B.second + 4*g)%(2*g); // positive modulo
      index_pair xp(A.first, B.first);
      auto y = elements.find(xp);
      if(y == elements.end() and interaction) y = elements.find(index_pair(B.first, A.first));
      if(y == elements.end()){
        ostringstream sout;
        sout << "Term " << x.first << " of parameter " << name << " is incompatible with symmetry generator " << i+1 << ", which maps it into " << xp;
        qcm_ED_throw(sout.str());
      }
      Complex ratio = x.second*phaseC[f]/y->second;
      if(abs(ratio-1.0) > 1e-8){
        ostringstream sout;
        sout << "Term " << x.first << " of parameter " << name << " is incompatible with symmetry generator " << i+1 << ", which maps it into " << xp << " but with value " << phaseC[f]*y->second << " instead of " << x.second;
        qcm_ED_throw(sout.str());
      }
    }
  }
}

template<>
inline double symmetry_group::phaseX<double>(int i){return phaseR[i];}

template<>
inline Complex symmetry_group::phaseX<Complex>(int i){return phaseC[i];}

#endif