#ifndef qbasis_h #define qbasis_h #define MKL_Complex16 std::complex #ifndef lapack_int #define lapack_int MKL_INT #endif #ifndef lapack_complex_double #define lapack_complex_double MKL_Complex16 #endif #include #include #include #include #include #include #include #include #include #include #include "mkl.h" namespace qbasis { enum class which_sym: int { /// full basis (no translational symmetry) full, /// representative basis of translational symmetry repr, /// variational basis vrnl }; } namespace qbasis { // -------------part 0: global variables, forward declarations ---------------- // ---------------------------------------------------------------------------- /** @file qbasis.h * \var const double machine_prec * \brief machine precision */ extern const double machine_prec; extern const double opr_precision; extern const double sparse_precision; extern const double lanczos_precision; /** @file qbasis.h * \var bool enable_ckpt * \brief use checkpoint and restart in Lanczos, for long runs */ extern bool enable_ckpt; /** @file qbasis.h * \fn void initialize(const bool &enable_ckpt_) * \brief initialize global variables & print out info */ void initialize(const bool &enable_ckpt_=false); template class multi_array; template void swap(multi_array&, multi_array&); /** \brief Multi-dimensional array */ template class multi_array { friend void swap <> (multi_array&, multi_array&); public: multi_array(): dim_(0), size_(0) {} multi_array(const std::vector &linear_size_input); multi_array(const std::vector &linear_size_input, const T &element); multi_array(const multi_array &old): // copy constructor dim_(old.dim_), size_(old.size_), linear_size_(old.linear_size_), data(old.data) {} multi_array(multi_array &&old) noexcept : // move constructor dim_(old.dim_), size_(old.size_), linear_size_(std::move(old.linear_size_)), data(std::move(old.data)) {} multi_array& operator=(multi_array old) { swap(*this, old); return *this; } multi_array& operator=(const T &element); ~multi_array() {} uint32_t dim() const { return dim_; } uint64_t size() const { return size_; } std::vector linear_size() const { return linear_size_; } T& index(const std::vector &pos); const T& index(const std::vector &pos) const; private: uint32_t dim_; uint64_t size_; std::vector linear_size_; std::vector data; }; typedef multi_array MltArray_uint32; typedef multi_array MltArray_double; typedef multi_array,std::vector>> MltArray_PairVec; class basis_prop; class mbasis_elem; template class wavefunction; template class opr; template class opr_prod; template class mopr; class lattice; template class csr_mat; template class model; // Operations on basis void basis_props_split(const std::vector &parent, std::vector &sub1, std::vector &sub2); bool q_bosonic(const std::vector &props); void swap(mbasis_elem&, mbasis_elem&); bool operator<(const mbasis_elem&, const mbasis_elem&); bool operator==(const mbasis_elem&, const mbasis_elem&); bool operator!=(const mbasis_elem&, const mbasis_elem&); /** \brief DEPRECATED! */ bool trans_equiv(const mbasis_elem&, const mbasis_elem&, const std::vector &props, const lattice&); // computational heavy, use with caution void zipper_basis(const std::vector &props_parent, // defined in Weisse's PRE 87, 043305 (2013) const std::vector &props_sub_a, const std::vector &props_sub_b, const mbasis_elem &sub_a, const mbasis_elem &sub_b, mbasis_elem &parent); void unzipper_basis(const std::vector &props_parent, const std::vector &props_sub_a, const std::vector &props_sub_b, const mbasis_elem &parent, mbasis_elem &sub_a, mbasis_elem &sub_b); // generate states compatible with given symmetry template void enumerate_basis(const std::vector &props, std::vector &basis, std::vector> conserve_lst = {}, std::vector val_lst = {}); /** @file qbasis.h * \fn void sort_basis_normal_order(std::vector &basis) * \brief sort basis with a simple "<" comparison */ void sort_basis_normal_order(std::vector &basis); /** @file qbasis.h * \fn void sort_basis_Lin_order(const std::vector &props, std::vector &basis) * \brief sort basis according to Lin Table convention (Ib, then Ia) */ void sort_basis_Lin_order(const std::vector &props, std::vector &basis); // generate Lin Tables for a given basis void fill_Lin_table(const std::vector &props, const std::vector &basis, std::vector &Lin_Ja, std::vector &Lin_Jb); /** \brief (sublattice) for a given list of full basis, find the reps according to translational symmetry. * Note: any state = T(disp2rep) * |rep> */ void classify_trans_full2rep(const std::vector &props, const std::vector &basis_all, const lattice &latt, const std::vector &trans_sym, std::vector &reps, std::vector &belong2rep, std::vector> &dist2rep); // (sublattice) for a given list of reps, find the corresponding translation group void classify_trans_rep2group(const std::vector &props, const std::vector &reps, const lattice &latt, const std::vector &trans_sym, const std::vector>,uint32_t>> &groups, std::vector &omega_g, std::vector &belong2group); // tabulate the maps for (ga,gb,ja,jb) -> (i,j), and (ga,gb,j) -> w void classify_Weisse_tables(const std::vector &props_parent, const std::vector &props_sub, const std::vector &basis_sub_repr, const lattice &latt_parent, const std::vector &trans_sym, const std::vector &belong2rep, const std::vector> &dist2rep, const std::vector &belong2group, const std::vector>,uint32_t>> &groups_parent, const std::vector>,uint32_t>> &groups_sub, MltArray_PairVec &Weisse_e_lt, MltArray_PairVec &Weisse_e_eq, MltArray_PairVec &Weisse_e_gt, MltArray_uint32 &Weisse_w_lt, MltArray_uint32 &Weisse_w_eq, MltArray_uint32 &Weisse_w_gt); // ^{-1} double norm_trans_repr(const std::vector &props, const mbasis_elem &repr, const lattice &latt_parent, const std::pair>,uint32_t> &group_parent, const std::vector &momentum); // Operations on wavefunctions template void swap(wavefunction&, wavefunction&); template wavefunction operator+(const wavefunction&, const wavefunction&); template wavefunction operator*(const wavefunction&, const T&); template wavefunction operator*(const T&, const wavefunction&); template wavefunction operator*(const mbasis_elem&, const T&); template wavefunction operator*(const T&, const mbasis_elem&); template T inner_product(const wavefunction&, const wavefunction&); // -------- operators ----------- template void swap(opr&, opr&); template bool operator==(const opr&, const opr&); template bool operator!=(const opr&, const opr&); template bool operator<(const opr&, const opr&); // only compares site, orbital and fermion, for sorting purpose template opr operator*(const T&, const opr&); template opr operator*(const opr&, const T&); template opr normalize(const opr&, T&); // sum_{i,j} mat[i,j]^2 == dim; the 1st nonzero element (in memory) be real positive template void swap(opr_prod&, opr_prod&); template bool operator==(const opr_prod&, const opr_prod&); template bool operator!=(const opr_prod&, const opr_prod&); template bool operator<(const opr_prod&, const opr_prod&); // compare only length, and if each lhs.mat_prod < rhs.mat_prod template opr_prod operator*(const opr_prod&, const opr_prod&); template opr_prod operator*(const opr_prod&, const opr&); template opr_prod operator*(const opr&, const opr_prod&); template opr_prod operator*(const opr_prod&, const T&); template opr_prod operator*(const T&, const opr_prod&); template opr_prod operator*(const opr&, const opr&); // cast up template void swap(mopr&, mopr&); template bool operator==(const mopr&, const mopr&); template bool operator!=(const mopr&, const mopr&); template mopr operator+(const mopr&, const mopr&); template mopr operator+(const mopr&, const opr_prod&); template mopr operator+(const opr_prod&, const mopr&); template mopr operator+(const mopr&, const opr&); template mopr operator+(const opr&, const mopr&); template mopr operator+(const opr_prod&, const opr_prod&); // cast up template mopr operator+(const opr_prod&, const opr&); // cast up template mopr operator+(const opr&, const opr_prod&); // cast up template mopr operator+(const opr&, const opr&); // cast up template mopr operator-(const mopr&, const mopr&); template mopr operator-(const mopr&, const opr_prod&); template mopr operator-(const opr_prod&, const mopr&); template mopr operator-(const mopr&, const opr&); template mopr operator-(const opr&, const mopr&); template mopr operator-(const opr_prod&, const opr_prod&); // cast up template mopr operator-(const opr_prod&, const opr&); // cast up template mopr operator-(const opr&, const opr_prod&); // cast up template mopr operator-(const opr&, const opr&); // cast up template mopr operator*(const mopr&, const mopr&); template mopr operator*(const mopr&, const opr_prod&); template mopr operator*(const opr_prod&, const mopr&); template mopr operator*(const mopr&, const opr&); template mopr operator*(const opr&, const mopr&); template mopr operator*(const mopr&, const T&); template mopr operator*(const T&, const mopr&); // opr * | orb0, orb1, ..., ORB, ... > = | orb0, orb1, ..., opr*ORB, ... >, fermionic sign has to be computed when traversing orbitals // wavefunction = opr * wavefunction (overwritten) template void oprXphi(const opr&, const std::vector&, wavefunction&, const bool &append = false); // wavefunction = opr * mbasis_elem template void oprXphi(const opr&, const std::vector&, wavefunction&, mbasis_elem, const bool &append = false); // wavefunction = opr * wavefunction0 template void oprXphi(const opr&, const std::vector&, wavefunction&, wavefunction phi0, const bool &append = false); // wavefunction = opr_prod * wavefunction (overwritten) template void oprXphi(const opr_prod&, const std::vector&, wavefunction&); // wavefunction = opr_prod * mbasis_elem template void oprXphi(const opr_prod&, const std::vector&, wavefunction&, mbasis_elem, const bool &append = false); // wavefunction = opr_prod * wavefunction0 template void oprXphi(const opr_prod&, const std::vector&, wavefunction&, wavefunction phi0, const bool &append = false); // wavefunction = mopr * mbasis_elem template void oprXphi(const mopr&, const std::vector&, wavefunction&, mbasis_elem, const bool &append = false); // wavefunction = mopr * wavefunction0 template void oprXphi(const mopr&, const std::vector&, wavefunction&, wavefunction phi0, const bool &append = false); // mopr * {a list of mbasis} -->> {a new list of mbasis} template void gen_mbasis_by_mopr(const mopr&, std::list&, const std::vector&, std::vector> conserve_lst = {}, std::vector val_lst = {}); void rm_mbasis_dulp_trans(const lattice&, std::list &, const std::vector&); // csr matrix template void swap(csr_mat&, csr_mat&); // divide into two identical sublattices, if Nsites even. To be used in the divide and conquer method lattice divide_lattice(const lattice &parent); int basis_disk_read(const std::string &filename, std::vector &basis); int basis_disk_write(const std::string &filename, const std::vector &basis); // -------------------- part 1: basis of the wave functions ------------------ // ---------------------------------------------------------------------------- // ------------ basic info of a particular basis ----------------- struct extra_info { uint8_t Nmax; // maximum number of particles // more items can be filled here... }; /** \brief The common properties of a given basis */ class basis_prop { public: basis_prop() = default; basis_prop(const uint32_t &n_sites, const uint8_t &dim_local_, const std::vector &Nf_map = std::vector(), const bool &dilute_ = false); /** \brief constructor of basis_prop, by proving number of sites and name of basis * * Implemented choices of basis name (s): * - spin-1/2 * - spin-1 * - dimer * - electron * - tJ * - spinless-fermion */ basis_prop(const uint32_t &n_sites, const std::string &s, const extra_info &ex = extra_info{0}); /** \brief question if corresponding to fermionic basis */ bool q_fermion() const { return (! Nfermion_map.empty()); } /** \brief split (site by site) into basis_prop for two basis */ void split(basis_prop &sub1, basis_prop &sub2) const; uint8_t dim_local; ///< local (single-site, single-orbital) dimension < 256 uint8_t bits_per_site; ///< <= 8 uint8_t bits_ignore; ///< for each orbital (with many sites), there are a few bits ignored uint16_t num_bytes; ///< for multi-orbital system, sum of num_bytes < 65536 uint32_t num_sites; std::vector Nfermion_map; ///< Nfermion_map[i] corresponds to the number of fermions of state i std::string name; ///< store the name of the basis bool dilute; ///< if dilute, bit-rep is not a good representation }; /** \brief Class for representing a quantum basis using bits * * Fundamental class for basis elements. * For given number of sites, and several orbitals, store the vectors of bits */ class mbasis_elem { friend void swap(mbasis_elem&, mbasis_elem&); friend bool operator<(const mbasis_elem&, const mbasis_elem&); friend bool operator==(const mbasis_elem&, const mbasis_elem&); friend bool trans_equiv(const mbasis_elem&, const mbasis_elem&, const std::vector &props, const lattice&); template friend class wavefunction; template friend void oprXphi(const opr&, const std::vector&, wavefunction&, const bool&); template friend void oprXphi(const opr&, const std::vector&, wavefunction&, mbasis_elem, const bool&); template friend void oprXphi(const opr_prod&, const std::vector&, wavefunction&, mbasis_elem, const bool&); template friend T inner_product(const wavefunction &, const wavefunction &); friend int basis_disk_read(const std::string&, std::vector&); friend int basis_disk_write(const std::string&, const std::vector&); public: /** \brief default constructor */ mbasis_elem() : mbits(nullptr) {} /** \brief constructor with its properties */ explicit mbasis_elem(const std::vector &props); /** \brief copy constructor */ mbasis_elem(const mbasis_elem& old); /** \brief move constructor */ mbasis_elem(mbasis_elem &&old) noexcept; /** \brief copy/move assignment constructor */ mbasis_elem& operator=(mbasis_elem old) { swap(*this, old); return *this; } /** \brief destructor */ ~mbasis_elem(); // ---------- status changes ----------- /** \brief read out a state from a given site and given orbital */ uint8_t siteRead(const std::vector &props, const uint32_t &site, const uint32_t &orbital) const; /** \brief write to a given site and given orbital */ mbasis_elem& siteWrite(const std::vector &props, const uint32_t &site, const uint32_t &orbital, const uint8_t &val); /** \brief reset all bits to 0 for a particular orbital */ mbasis_elem& reset(const std::vector &props, const uint32_t &orbital); /** \brief reset all bits to 0 in all orbitals */ mbasis_elem& reset(); /** \brief change mbasis_elem to the next available state, for a particular orbital */ mbasis_elem& increment(const std::vector &props, const uint32_t &orbital); /** \brief change mbasis_elem to the next available state */ mbasis_elem& increment(const std::vector &props); /** \brief print the basis in the bit representation */ void prt_bits(const std::vector &props) const; /** \brief print the basis in numbers, each corresponding to a state on one site (vacuum not printed) */ void prt_states(const std::vector &props) const; /** \brief generate a file to plot the basis on a lattice */ void plot_states(const std::vector &props, const lattice &latt, const std::string &filename) const; // ----------- basic inquiries ---------- /** \brief question if the vacuum state on a particular site */ [[nodiscard]] bool q_zero_site(const std::vector &props, const uint32_t &site) const; /** \brief question if the vacuum state on a particular orbital */ [[nodiscard]] bool q_zero_orbital(const std::vector &props, const uint32_t &orbital) const; /** \brief question if the vacuum state */ [[nodiscard]] bool q_zero() const; /** \brief question if every site occupied by the highest state, for a given orbital */ [[nodiscard]] bool q_maximized(const std::vector &props, const uint32_t &orbital) const; /** \brief question if every site occupied by the highest state */ [[nodiscard]] bool q_maximized(const std::vector &props) const; /** \brief question if every site occupied by the same state, for a given orbital */ [[nodiscard]] bool q_same_state_all_site(const std::vector &props, const uint32_t &orbital) const; /** \brief question if every site occupied by the same state */ [[nodiscard]] bool q_same_state_all_site(const std::vector &props) const; // get a label // preferred size of work: 2*num_sites uint64_t label(const std::vector &props, const uint32_t &orbital, std::vector &work) const; // preferred size of work1: 2*num_sites, work2: 2*num_orbs uint64_t label(const std::vector &props, std::vector &work1, std::vector &work2) const; // preferred size of work: 2*num_sites void label_sub(const std::vector &props, const uint32_t &orbital, uint64_t &label1, uint64_t &label2, std::vector &work) const; // preferred size of work1: 2*num_sites, work2: 4*num_orbs void label_sub(const std::vector &props, uint64_t &label1, uint64_t &label2, std::vector &work1, std::vector &work2) const; // return a vector of length dim_local (for orbital), reporting # of each state [[nodiscard]] std::vector statistics(const std::vector &props, const uint32_t &orbital) const; // a direct product of the statistics of all orbitals, size: dim_orb1 * dim_orb2 * ... [[nodiscard]] std::vector statistics(const std::vector &props) const; /** \brief output the center of the non-vacuum states, in units of the lattice basis \f$ \vec{a}_i \f$ */ [[nodiscard]] std::vector center_pos(const std::vector &props, const lattice &latt) const; /** \brief transform the basis according to the given plan (single given orbital). * sgn = 0 or 1 denoting if extra minus sign generated by transforming fermions */ mbasis_elem& transform(const std::vector &props, const std::vector &plan, int &sgn, const uint32_t &orbital); /** \brief transform the basis according to the given plan (all orbitals transform in same way): site1 -> site2 */ mbasis_elem& transform(const std::vector &props, const std::vector &plan, int &sgn); /** \brief transform the basis according to the given plan (different orbs transform in different ways): (site1, orb1) -> (site2, orb2) * \param plan: outer vector (orbital), middle vector (site), inner pair(first=site, second=orbital) */ mbasis_elem& transform(const std::vector &props, const std::vector>> &plan, int &sgn); /* mbasis_elem& translate(const std::vector &props, const lattice &latt, const std::vector &disp, int &sgn, const uint32_t &orbital); mbasis_elem& translate(const std::vector &props, const lattice &latt, const std::vector &disp, int &sgn); */ /** \brief (OBC assumed; state 0 assumed to be vacuum) old state = T(-disp_vec) * centralized state */ mbasis_elem& translate2center_OBC(const std::vector &props, const lattice &latt, std::vector &disp_vec); // ------------ measurements ------------ [[nodiscard]] double diagonal_operator(const std::vector &props, const opr &lhs) const; [[nodiscard]] std::complex diagonal_operator(const std::vector &props, const opr> &lhs) const; [[nodiscard]] double diagonal_operator(const std::vector &props, const opr_prod &lhs) const; [[nodiscard]] std::complex diagonal_operator(const std::vector &props, const opr_prod> &lhs) const; [[nodiscard]] double diagonal_operator(const std::vector &props, const mopr &lhs) const; [[nodiscard]] std::complex diagonal_operator(const std::vector &props, const mopr> &lhs) const; private: /** store an array of basis elements, for multi-orbital site (or unit cell) * the first 2 bytes are used to store the total number of bytes used by this array * * in terms of bits when perfoming "<" comparison (e.g. bits_per_site = 2): * e.g. 0 1 0 1 0 0 0 0 1 1 0 1 0 0 0 0 1 1 1 0 * ^ ^ * / \ * bits[1] bits[0] * \ / * site[0] * Note1: for arrangement of orbitals, they are similar: * ... orb[2] orb[1] orb[0] * Note2: The wavefunction is always defined as: * |alpha_0, beta_1, gamma_2, ... > = alpha_0^\dagger beta_1^\dagger gamma_2^\dagger ... |GS> * (where alpha_i^\dagger is creation operator of state alpha on site i) */ uint8_t* mbits; }; // -------------- class for wave functions --------------- // Use with caution, may hurt speed when not used properly template class wavefunction { friend class model; friend void swap <> (wavefunction &, wavefunction &); friend wavefunction operator+ <> (const wavefunction&, const wavefunction&); friend void oprXphi <> (const opr&, const std::vector&, wavefunction&, const bool&); friend void oprXphi <> (const opr&, const std::vector&, wavefunction&, mbasis_elem, const bool&); friend void oprXphi <> (const opr&, const std::vector&, wavefunction&, wavefunction, const bool&); friend void oprXphi <> (const opr_prod&, const std::vector&, wavefunction&, mbasis_elem, const bool&); friend void oprXphi <> (const opr_prod&, const std::vector&, wavefunction&, wavefunction, const bool&); friend void oprXphi <> (const mopr&, const std::vector&, wavefunction&, mbasis_elem, const bool&); friend void oprXphi <> (const mopr&, const std::vector&, wavefunction&, wavefunction, const bool&); friend T inner_product <> (const wavefunction &, const wavefunction &); friend void gen_mbasis_by_mopr <> (const mopr &, std::list&, const std::vector&, std::vector> conserve_lst, std::vector val_lst); public: // constructor from props explicit wavefunction(const std::vector &props, const MKL_INT &capacity = 64); // constructor from an element explicit wavefunction(const mbasis_elem &old, const MKL_INT &capacity = 64); explicit wavefunction(mbasis_elem &&old, const MKL_INT &capacity = 64); // copy constructor wavefunction(const wavefunction &old) : bgn(old.bgn), end(old.end), total_bytes(old.total_bytes), ele(old.ele) {} // move constructor wavefunction(wavefunction &&old) noexcept : bgn(old.bgn), end(old.end), total_bytes(old.total_bytes), ele(std::move(old.ele)) {} // copy assignment constructor and move assignment constructor, using "swap and copy" wavefunction& operator=(wavefunction old) { swap(*this, old); return *this; } // destructor ~wavefunction() = default; // ----------- basic inquiries ---------- MKL_INT size() const { MKL_INT capacity = static_cast(ele.size()); return ( (end + capacity - bgn) % capacity ); } // check if zero [[nodiscard]] bool q_empty() const { return (size() == 0); } // check if full [[nodiscard]] bool q_full() const { return (size() + 1 == static_cast(ele.size())); } // check if sorted [[nodiscard]] bool q_sorted() const; // check if sorted and there are no dulplicated terms [[nodiscard]] bool q_sorted_fully() const; // for \sum_i \alpha_i * element[i], return \sum_i |\alpha_i|^2 [[nodiscard]] double amplitude() const; // ---------------- print --------------- void prt_bits(const std::vector &props) const; void prt_states(const std::vector &props) const; // ----------- element access ----------- std::pair& operator[](MKL_INT n); const std::pair& operator[](MKL_INT n) const; // ------------ arithmetics ------------- wavefunction& reserve(const MKL_INT& capacity_new); // enlarge capacity wavefunction& clear(); // clear but memory not released wavefunction& operator+=(mbasis_elem ele_new); wavefunction& operator+=(std::pair ele_new); wavefunction& operator+=(wavefunction rhs); wavefunction& operator*=(const T &rhs); // simplify wavefunction& simplify(); private: // store an array of basis elements, and their corresponding coefficients // note: there should not be any dulplicated elements // using cycled queue data structure (wasting 1 element to avoid end==bgn confusion: empty or full) MKL_INT bgn; MKL_INT end; int total_bytes; std::vector> ele; // add one basis element (promised not to overlap) wavefunction& add(const mbasis_elem &rhs); // copy from one basis element (promised not to overlap) wavefunction& copy(const mbasis_elem &rhs); // add one basis element (promised not to overlap) wavefunction& add(const std::pair &rhs); // copy from one basis element (promised not to overlap) wavefunction& copy(const std::pair &rhs); // add another wavefunction (promised not to overlap) wavefunction& add(const wavefunction &rhs); // copy from another wavefuction (promised not to overlap) wavefunction& copy(const wavefunction &rhs); }; // ------------------------------ part 2: operators -------------------------- // ---------------------------------------------------------------------------- /** \brief An operator on a given site and orbital * (fundamental class for operators) */ template class opr { friend void swap <> (opr&, opr&); friend bool operator== <> (const opr&, const opr&); friend bool operator!= <> (const opr&, const opr&); friend bool operator< <> (const opr&, const opr&); friend opr operator* <> (const T&, const opr&); friend opr operator* <> (const opr&, const T&); friend opr normalize <> (const opr&, T&); friend class opr_prod; friend class mopr; friend class mbasis_elem; friend void oprXphi <> (const opr&, const std::vector&, wavefunction&, const bool&); friend void oprXphi <> (const opr&, const std::vector&, wavefunction&, mbasis_elem, const bool&); friend void oprXphi <> (const opr&, const std::vector&, wavefunction&, wavefunction, const bool&); public: /** \brief default constructor */ opr() : mat(nullptr) {} /** \brief constructor from diagonal elements */ opr(const uint32_t &site_, const uint32_t &orbital_, const bool &fermion_, const std::vector &mat_); /** \brief constructor from a matrix */ opr(const uint32_t &site_, const uint32_t &orbital_, const bool &fermion_, const std::vector> &mat_); /** \brief copy constructor */ opr(const opr &old); /** \brief move constructor */ opr(opr &&old) noexcept; /** \brief copy/move assignment constructor */ opr& operator=(opr old) { swap(*this, old); return *this; } /** \brief destructor */ ~opr(); /** \brief print details of opr */ void prt() const; // ----------- basic inquiries ---------- /** \brief question if it is zero operator */ [[nodiscard]] bool q_zero() const; /** \brief question if it is identity operator */ [[nodiscard]] bool q_diagonal() const { return diagonal; } /** \brief question if it is identity operator */ [[nodiscard]] bool q_identity() const; /** \brief return site index */ [[nodiscard]] uint32_t pos_site() const { return site; } /** \brief return orbital index */ [[nodiscard]] uint32_t pos_orb() const { return orbital; } // ------------ arithmetics ------------- /** * \brief \f$ \sqrt { \sum_{i,j} |M_{i,j}|^2 } \f$, * where \f$ M \f$ is the matrix representation of the operator */ [[nodiscard]] double norm() const; /** \brief simplify the structure if possible */ opr& simplify(); /** \brief invert the sign */ opr& negative(); /** \brief take Hermitian conjugate */ opr& dagger(); /** \brief change site index */ opr& change_site(const uint32_t &site_); /** \brief change the site index of operator * * Say plan[i]=j, it means site i becomes site j after the transformation. * If opr has site index i, then after transformation, the site index becomes j. */ opr& transform(const std::vector &plan); /** \brief opr = opr + rhs */ opr& operator+=(const opr &rhs); /** \brief opr = opr - rhs */ opr& operator-=(const opr &rhs); /** \brief opr = opr * rhs */ opr& operator*=(const opr &rhs); /** \brief opr = opr * rhs */ opr& operator*=(const T &rhs); private: uint32_t site; // site No. uint32_t orbital; // orbital No. uint32_t dim; // number of rows(columns) of the matrix bool fermion; // fermion or not bool diagonal; // diagonal in matrix form T *mat; // matrix form, or diagonal elements if diagonal }; // -------------- class for operator products ---------------- // note: when mat_prod is empty, it represents identity operator, with coefficient coeff // all matrices in this class should have the same type (may decrease effiency, think later how we can improve) template class opr_prod { friend void swap <> (opr_prod&, opr_prod&); friend bool operator== <> (const opr_prod&, const opr_prod&); friend bool operator!= <> (const opr_prod&, const opr_prod&); friend bool operator< <> (const opr_prod&, const opr_prod&); friend opr_prod operator* <> (const opr_prod&, const opr_prod&); friend opr_prod operator* <> (const opr_prod&, const opr&); friend opr_prod operator* <> (const opr&, const opr_prod&); friend opr_prod operator* <> (const opr_prod&, const T&); friend opr_prod operator* <> (const T&, const opr_prod&); friend opr_prod operator* <> (const opr&, const opr&); friend class mopr; friend class mbasis_elem; friend void oprXphi <> (const opr_prod&, const std::vector&, wavefunction&); friend void oprXphi <> (const opr_prod&, const std::vector&, wavefunction&, mbasis_elem, const bool&); friend void oprXphi <> (const opr_prod&, const std::vector&, wavefunction&, wavefunction, const bool&); public: // default constructor opr_prod() = default; // constructor from one fundamental operator explicit opr_prod(const opr &ele); // copy constructor opr_prod(const opr_prod &old): coeff(old.coeff), mat_prod(old.mat_prod) {} // move constructor opr_prod(opr_prod &&old) noexcept : coeff(std::move(old.coeff)), mat_prod(std::move(old.mat_prod)) {} // copy assignment constructor and move assignment constructor, using "swap and copy" opr_prod& operator=(opr_prod old) { swap(*this, old); return *this; } // destructor ~opr_prod() = default; void prt() const; // ----------- basic inquiries ---------- // question if it is zero operator [[nodiscard]] bool q_zero() const; // question if each opr is diagonal [[nodiscard]] bool q_diagonal() const; // question if it is proportional to identity operator [[nodiscard]] bool q_prop_identity() const; [[nodiscard]] uint32_t len() const; opr& operator[](uint32_t n); const opr& operator[](uint32_t n) const; // ------------ arithmetics ------------- // invert the sign opr_prod& negative(); // Hermitian conjugate (not sorted after operation) opr_prod& dagger(); // sort according to {site, orbital}. If exchanging fermions with odd times, invert sign. // replace in future, if used heavily opr_prod& bubble_sort(); opr_prod& transform(const std::vector &plan); opr_prod& operator*=(opr rhs); opr_prod& operator*=(opr_prod rhs); // in this form to avoid self-assignment opr_prod& operator*=(const T &rhs); private: T coeff; std::list> mat_prod; // each opr should be normalized }; // -------------- class for a combination of operators ---------------- // a linear combination of products of operators template class mopr { friend class model; friend void swap <> (mopr&, mopr&); friend bool operator== <> (const mopr&, const mopr&); friend bool operator!= <> (const mopr&, const mopr&); friend mopr operator+ <> (const mopr&, const mopr&); friend mopr operator+ <> (const mopr&, const opr_prod&); friend mopr operator+ <> (const opr_prod&, const mopr&); friend mopr operator+ <> (const mopr&, const opr&); friend mopr operator+ <> (const opr&, const mopr&); friend mopr operator+ <> (const opr_prod&, const opr_prod&); friend mopr operator+ <> (const opr_prod&, const opr&); friend mopr operator+ <> (const opr&, const opr_prod&); friend mopr operator+ <> (const opr&, const opr&); friend mopr operator- <> (const mopr&, const mopr&); friend mopr operator- <> (const mopr&, const opr_prod&); friend mopr operator- <> (const opr_prod&, const mopr&); friend mopr operator- <> (const mopr&, const opr&); friend mopr operator- <> (const opr&, const mopr&); friend mopr operator- <> (const opr_prod&, const opr_prod&); friend mopr operator- <> (const opr_prod&, const opr&); friend mopr operator- <> (const opr&, const opr_prod&); friend mopr operator- <> (const opr&, const opr&); friend mopr operator* <> (const mopr&, const mopr&); friend mopr operator* <> (const mopr&, const opr_prod&); friend mopr operator* <> (const opr_prod&, const mopr&); friend mopr operator* <> (const mopr&, const opr&); friend mopr operator* <> (const opr&, const mopr&); friend mopr operator* <> (const mopr&, const T&); friend mopr operator* <> (const T&, const mopr&); friend void oprXphi <> (const mopr&, const std::vector&, wavefunction&, mbasis_elem, const bool&); friend void oprXphi <> (const mopr&, const std::vector&, wavefunction&, wavefunction, const bool&); friend void gen_mbasis_by_mopr <> (const mopr &, std::list&, const std::vector&, std::vector> conserve_lst, std::vector val_lst); public: // default constructor mopr() = default; // constructor from one fundamental operator explicit mopr(const opr &ele); // constructor from operator products explicit mopr(const opr_prod &ele); // copy constructor mopr(const mopr &old): mats(old.mats) {} // move constructor mopr(mopr &&old) noexcept : mats(std::move(old.mats)) {} // copy assignment constructor and move assignment constructor, using "swap and copy" mopr& operator=(mopr old) { swap(*this, old); return *this; } // destructor ~mopr() = default; void prt() const; // ----------- basic inquiries ---------- [[nodiscard]] bool q_zero() const { return mats.empty(); } // question if each opr_prod is diagonal [[nodiscard]] bool q_diagonal() const; [[nodiscard]] uint32_t size() const { return static_cast(mats.size()); } opr_prod& operator[](uint32_t n); const opr_prod& operator[](uint32_t n) const; // ------------ arithmetics ------------- mopr& simplify(); // invert the sign mopr& negative(); // Hermitian conjugate mopr& dagger(); mopr& transform(const std::vector &plan); mopr& operator+=(opr rhs); mopr& operator+=(opr_prod rhs); mopr& operator+=(mopr rhs); mopr& operator-=(opr rhs); mopr& operator-=(opr_prod rhs); mopr& operator-=(mopr rhs); mopr& operator*=(opr rhs); mopr& operator*=(opr_prod rhs); mopr& operator*=(mopr rhs); mopr& operator*=(const T &rhs); private: // the outer list represents the sum of operators, inner data structure taken care by operator products std::list> mats; }; // -------------------------- part 3: sparse matrices ------------------------ // ---------------------------------------------------------------------------- // Note: sparse matrices in this code are using zero-based convention // By default, all diagonal elements are stored, even if they are zero (to be compatible with pardiso) template struct lil_mat_elem { T val; MKL_INT col; }; template class lil_mat { friend class csr_mat; public: // default constructor lil_mat() = default; // constructor with the Hilbert space dimension explicit lil_mat(const MKL_INT &n, bool sym_ = false); // add one element void add(const MKL_INT &row, const MKL_INT &col, const T &val); // explicitly destroy, free space void destroy(const MKL_INT &row); void destroy(); // destructor ~lil_mat() = default;; MKL_INT dimension() const { return dim; } MKL_INT num_nonzero() const { return nnz; } bool q_sym() { return sym; } // print void prt() const; void use_full_matrix() { sym = false; } private: MKL_INT dim; // dimension of the matrix MKL_INT nnz; // number of non-zero entries bool sym; // if storing only upper triangle std::vector>> mat; }; // 3-array form of csr sparse matrix format, zero based template class csr_mat { friend void swap <> (csr_mat&, csr_mat&); public: MKL_INT dim; MKL_INT nnz; // number of non-zero entries bool sym; // if storing only upper triangle T *val; MKL_INT *ja; MKL_INT *ia; sparse_matrix_t handle; // default constructor csr_mat() : dim(0), nnz(0), sym(false), val(nullptr), ja(nullptr), ia(nullptr), handle(nullptr) {} // copy constructor csr_mat(const csr_mat &old); // move constructor csr_mat(csr_mat &&old) noexcept; // copy assignment constructor and move assignment constructor, using "swap and copy" csr_mat& operator=(csr_mat old) { swap(*this, old); return *this; } // explicitly destroy, free space void destroy(); // destructor ~csr_mat(); // print void prt() const; MKL_INT dimension() const { return dim; } // construcotr from lil_mat, and if sym_ == true, use only the upper triangle // then destroy the lil_mat explicit csr_mat(lil_mat &old); // matrix vector product // y = H * x + y void MultMv2(const T *x, T *y) const; // y = H * x void MultMv(const T *x, T *y) const; std::vector to_dense() const; }; // ------------------------------part 4: Lanczos ------------------------------ // ---------------------------------------------------------------------------- // m step of Lanczos (for iram: m = k + np; for simple Lanczos, m depends on the convergence speed) // on entry, assuming k steps of Lanczos already performed // v of length m+1, hessenberg matrix of size m*m (m-step Lanczos) // after decomposition, mat * v[0:m-1] = v[0:m-1] * hessenberg + b[m] * v[m] * e_m^T, // where e_m has only one nonzero element: e[0:m-2] == 0, e[m-1] = 1 // // maxit: maximum allowed Lanczos steps, m < maxit // maxit also serves as the leading dimension of hessenberg matrix // a[j] = hessenberg[j+maxit], diagonal of hessenberg matrix // b[j] = hessenberg[j] // a[0] b[1] -> note: b[0] not used // b[1] a[1] b[2] // b[2] a[2] b[3] // b[3] a[3] b[4] // .. .. .. b[k-1] // b[k-1] a[k-1] // on entry: // a[0], ..., a[k-1] in hessenberg matrix // b[0], ..., b[k-1], b[k] in hessenberg matrix // on exit: // a[0], ..., a[m-1] in hessenberg matrix // b[0], ..., b[m-1], b[m] in hessenberg matrix // // if purpose == "iram": // on entry: // v[0], ..., v[k] stored in v // on exit: // v[0], ..., v[m] stored in v // // if purpose == "sr_val0" (smallest eigenvalue), or purpose == "dnmcs": // on entry: // v[k-1], v[k] stored in v. If k%2==0, stored as {v[k],v[k-1]}; else stored as {v[k-1],v[k]} // on exit: // v[m-1], v[m] stored in v. If m%2==0, stored as {v[m],v[m-1]}; else stored as {v[m-1],v[m]} // // if purpose == "sr_vec0" (smallest eigenvector): // same as "sr_val0", but with one extra column storing the eigenvector template void lanczos(MKL_INT k, MKL_INT np, const MKL_INT &maxit, MKL_INT &m, const MKL_INT &dim, const MAT &mat, T v[], double hessenberg[], const std::string &purpose); // Iterative sparse solver using conjugate gradient method // given well converged eigenvalue E0, and a good initital guess v[0], solves (H - E0) * v[j] = 0 // on entry: assuming m-step finished, v contains the initital guess v[m] // r[0] = -(H-E0)*v[0], p[0] = r[0] // on exit: v rewritten by the converged solution template void eigenvec_CG(const MKL_INT &dim, const MKL_INT &maxit, MKL_INT &m, const MAT &mat, const T &E0, double &accu, T v[], T r[], T p[], T pp[]); // compute eigenvalues and eigenvectors of hessenberg matrix // on entry, hessenberg should have leading dimension maxit // on exit, ritz of size m, s of size m*m // order = "sm", "lm", "sr", "lr", where 's': small, 'l': large, 'm': magnitude, 'r': real part void hess_eigen(const double hessenberg[], const MKL_INT &maxit, const MKL_INT &m, const std::string &order, std::vector &ritz, std::vector &s); // implicitly restarted Arnoldi method // nev: number of eigenvalues needed // ncv: length of each individual lanczos process // 2 < nev + 2 <= ncv template void iram(const MKL_INT &dim, MAT &mat, T v0[], const MKL_INT &nev, const MKL_INT &ncv, const MKL_INT &maxit, const std::string &order, MKL_INT &nconv, double eigenvals[], T eigenvecs[]); // feast method to get eigenvalues inside [E_low, E_high] // nev: number of eigenvalues needed // nconv: number of eigenvalues converged template void call_feast(csr_mat &mat, MKL_INT nev, const double &E_low, const double &E_high, MKL_INT max_loops, MKL_INT &nconv, double eigenvals[], T eigenvecs[]); // ----------------------------- part 5: Lattices ---------------------------- // ---------------------------------------------------------------------------- /** In general, the position of any site can be represented as: * \f[ * \vec{r} = \sum_{i=1}^{D} m_i^0 \vec{a}_i + \sum_{i=1}^D M_i \vec{A}_i + \vec{d}_s * = \sum_{i=1}^{D} m_i \vec{a}_i + \vec{d}_s, * \f] * where \f$ m_i^0 \f$, \f$ m_i \f$, and \f$ M_i \f$ are integers, * and \f$ \vec{d}_s \f$ are the position shift of sublattices. * Note: \f$ m_i^0 \f$ is located in the 1st supercell defined by * \f[ * \sum_{i=1}^D m_i^0 \vec{a}_i = \sum_{i=1}^D x_i \vec{A}_i, * \f] * where \f$ 0 \le x_i < 1 \f$. * In other words, any site has two coordinate labels, one within the supercell, * another as the supercell coordinate. * Most of time, the supercell coordinate is set at origin; but they can be effective in some other cases. */ class lattice { friend lattice divide_lattice(const lattice &parent); public: lattice() = default; /** \brief constructor from particular requirements. e.g. square, triangular... */ lattice(const std::string &name, const std::vector &L_, const std::vector &bc_, bool auto_dim_spec = true); /** \brief constructor from manually designed lattices, specified from input files */ explicit lattice(const std::string &filename, bool auto_dim_spec = true); /** \brief if all directions untilted, return true */ bool q_tilted() const; bool q_dividable() const; /** \brief for given \f$ m_i \f$, return its corresponding \f$ m_i^0 \f$ in the 1st supercell, and the supercell index \f$ M_i \f$ */ void coor2supercell0(const int *coor, int *coor0, int *M) const; /** \brief with given \f$ m_i \f$, output the site index */ void coor2site(const std::vector &coor, const int &sub, uint32_t &site, std::vector &work) const; // coordinates <-> site indices // first find a direction (dim_spec) which has even size, if not successful, use the following: // 1D: site = i * num_sub + sub // 2D: site = (i + j * L[0]) * num_sub + sub // 3D: site = (i + j * L[0] + k * L[0] * L[1]) * num_sub + sub // otherwise, the dim_spec should be counted first void coor2site_old(const std::vector &coor, const int &sub, uint32_t &site) const; /** \brief for any site, output the coordinate \f$ m_i \f$ (NOT \f$ m_i^0 \f$). */ void site2coor(std::vector &coor, int &sub, const uint32_t &site) const; /** \brief with given \f$ m_i \f$, output cartisian coordinates */ void coor2cart(const std::vector &coor, std::vector &cart) const; /** \brief with given \f$ m_i \f$ and sublattice index, output cartisian coordinates */ void coor2cart(const std::vector &coor, std::vector &cart, const int &sub) const; /** with given \f$ \vec{k}=\sum_i k_i \vec{b}_i \f$, * return \f$ \vec{K} = \sum_i K_i \vec{B}_i \f$ and \f$ \tilde{k} = \sum_i \tilde{k}_i \vec{B}_i \f$, * which satisfies \f$ \vec{k} = \vec{K} + \tilde{k} \f$. * Here \f$ K_i \f$ are integers, and \f$ \tilde{k}_i \in [0,1) \f$. */ void k2superBZ(const double *k, int *K, double *ktilde) const; /** \brief on return, the vector "plan" contains the positions of each site after translation T(disp) */ void translation_plan(std::vector &plan, const std::vector &disp, std::vector &scratch_coor, std::vector &scratch_work) const; /** \brief return a vector containing the positions of each site after rotation (2D-only for now): * x -> x', where (x' - x0) = R (x - x0) */ std::vector rotation_plan(const uint32_t &origin, const double &angle) const; /** \brief return a vector containing the positions of each site after reflection of a given plane (2D-only for now): * x -> x', where the plane is given by a*x + b*y + c == 0 */ std::vector reflection_plan(const std::vector &abc) const; /** \brief combine two plans \f$P_\text{total} = P_1 P_2\f$. Length of P = number of orbitals, * and length of P[0] = number of sites. Each pair.first = new site, Each pair.second = new orbital */ std::vector>> plan_product(const std::vector>> &P1, const std::vector>> &P2) const; // inverse of a transformation std::vector>> plan_inverse(const std::vector>> &old) const; uint32_t dimension_spec() const { return dim_spec; } /** \brief return the center of the lattice */ std::vector center_pos() const { return center; } uint32_t Lx() const { return L[0]; } uint32_t Ly() const { if (L.size() > 1) return L[1]; else return 1; } uint32_t Lz() const { if (L.size() > 2) return L[2]; else return 1; } // obtain all possible divisors of a lattice, for the divide and conquer method // the returned value is a list of lists: {{divisors for Lx}, {divisors for Ly}, ...} std::vector> divisor_v1(const std::vector &trans_sym) const; // {{1,1,1}, {1,1,2}, {1,1,5}, {1,2,1}, {1,2,2}, {1,2,5},...}, i.e., combine results from v1 to a single list //std::vector> divisor_v2(const std::vector &trans_sym) const; // enumerate all possible commensurate magnetic bravis basis, ordered by the unit cell size // note: equally, it gives all the possible translational subgroups (without dulplicates) // in the following format: // 1D: { pair( {{1}}, size=1 ), // pair( {{2}}, size=2 ), // ... // } // 2D: { pair( {{1,0},{0,1}}, size=1 ), // pair( {{1,0},{0,1}}, size=1 ), // ..., // } // 3D: { pair( {{1,0,0},{0,1,0},{0,0,1}}, size=1 ), // pair( {{2,0,0},{0,1,0},{0,0,1}}, size=2 ), // ..., // } std::vector>,uint32_t>> trans_subgroups(const std::vector &trans_sym) const; std::vector> a; //!< real space basis {\f$ \vec{a}_0,\, \vec{a}_1,\,\ldots \f$} std::vector> b; //!< reciprocal space basis \f$ \vec{b}_i \f$ std::vector> A; //!< superlattice basis, in units of \f$ \vec{a}_i \f$ std::vector Amat; //!< superlattice basis in matrix format (column j -> A_j) std::vector Bmat; //!< reciprocal superlattice basis in matrix format (column i -> B_i), units: \f$ \vec{b}_i \f$ std::vector> pos_sub; //!< position shift \f$ \vec{d}_s \f$ of each sublattice, in units of \f$ \vec{a}_i \f$ std::vector bc; //!< boundary condition. Only explicitly used when all in 1st supercell, otherwise too complicated uint32_t dim; //!< dimension of lattice uint32_t num_sub; //!< number of sublattices uint32_t Nsites; //!< total number of sites in the 1st supercell std::vector L; //!< linear size in each dimension private: // coordinates <-> site indices // first find a direction (dim_spec) which has even size, if not successful, use the following: // 1D: site = i * num_sub + sub // 2D: site = (i + j * L[0]) * num_sub + sub // 3D: site = (i + j * L[0] + k * L[0] * L[1]) * num_sub + sub // otherwise, the dim_spec should be counted first void site2coor_old(std::vector &coor, int &sub, const uint32_t &site) const; std::vector tilted; //!< \brief true if {A1,A2,...} unparallel to {a1,a2,...} uint32_t dim_spec; //!< \brief the code starts labeling sites from a dimension which has even # of sites (if applicable) std::vector,int>> site2coor_map; //!< store the coordinate \f$ m_i \f$ (NOT \f$ m_i^0 \f$), and sublattice index for each site std::vector> site2super_map; //!< store the supercell index \f$ M_i \f$ for each site (empty vector if all in 1st supercell) std::vector,uint32_t>> coor2site_map; //!< for given \f$ m_i^0 \f$ (in the 1st supercell) and sublattice index, return the site label std::vector center; //!< center of mass of the lattice, in units of \f$ \vec{a}_i \f$ // one more variable here, denoting the divide and conquer partition // if empty(false), then force to store the matrix when working with translational symmetry }; // ---------------part 6: Routines to construct Hamiltonian ------------------- // ---------------------------------------------------------------------------- template class model { public: bool matrix_free{}; ///< if generating matrix on the fly std::vector props, props_sub_a, props_sub_b; mopr Ham_diag; ///< diagonal part of H mopr Ham_off_diag; ///< offdiagonal part of H mopr Ham_vrnl; ///< used for generating Trugman's basis MKL_INT nconv{}; // controls which sector of basis to be active // by default sec_full = 0 (e.g. Sz=0 ground state sector of Heisenberg model); // when needed, sec_full will be switched to 1 to activate another sector (e.g. Sz=1 sector), // such setting can avoid messing up the code when calculating correlation functions which_sym sec_sym{}; ///< choose from {full, repr, vrnl} uint32_t sec_mat{}; ///< which sector the matrix is relevant. std::vector dim_full; std::vector dim_repr; std::vector dim_vrnl; std::vector> momenta; std::vector> momenta_vrnl; /** \brief full basis without translation sym */ std::vector> basis_full; /** \brief basis with translation sym */ std::vector> basis_repr; /** \brief variational basis for Trugman's method */ std::vector> basis_vrnl; /** \brief ground state representative for Trugman's method */ mbasis_elem gs_vrnl; /** \brief basis for half lattice, used for building Weisse Table */ std::vector basis_sub_full; /** \brief reps for half lattice */ std::vector basis_sub_repr; // Lin tables, for both full basis and translation basis std::vector> Lin_Ja_full; std::vector> Lin_Jb_full; std::vector> Lin_Ja_repr; std::vector> Lin_Jb_repr; /** \brief 1 / */ std::vector> norm_repr; /** \brief 1 / = omega_g, for the variational vacuum state */ std::vector gs_norm_vrnl; /** \brief omega_g, orbital size of the variational vacuum state */ uint32_t gs_omegaG_vrnl{}; /** \brief ground state energy of the variational vacuum state */ double gs_E0_vrnl{}; /** \brief ground state momentum */ std::vector gs_momentum_vrnl; std::vector> HamMat_csr_full; std::vector> HamMat_csr_repr; std::vector> HamMat_csr_vrnl; std::vector eigenvals_full; std::vector eigenvecs_full; std::vector eigenvals_repr; std::vector> eigenvecs_repr; std::vector eigenvals_vrnl; std::vector> eigenvecs_vrnl; // ---------------- deprecated -------------------- std::vector> basis_belong_deprec; // size: dim_target_full, store the position of its repr std::vector>> basis_coeff_deprec; // size: dim_target_full, store the coeff std::vector> basis_repr_deprec; // ---------------- deprecated -------------------- model() = default; explicit model(lattice latt, const uint32_t &num_secs = 5, const double &fake_pos_ = 100.0); ~model() = default; void prt_Ham_diag() const { Ham_diag.prt(); } void prt_Ham_offdiag() const { Ham_off_diag.prt(); } void add_orbital(const uint32_t &n_sites, const uint8_t &dim_local_, const std::vector &Nf_map = std::vector(), const bool &dilute_ = false) { props.emplace_back(n_sites,dim_local_,Nf_map,dilute_); basis_props_split(props, props_sub_a, props_sub_b); } void add_orbital(const uint32_t &n_sites, const std::string &s, const extra_info &ex = extra_info{0}) { props.emplace_back(n_sites, s, ex); basis_props_split(props, props_sub_a, props_sub_b); } uint32_t local_dimension() const; void add_Ham(const opr &rhs); void add_Ham(const opr_prod &rhs); void add_Ham(const mopr &rhs); void add_Ham_vrnl(const opr &rhs); void add_Ham_vrnl(const opr_prod &rhs); void add_Ham_vrnl(const mopr &rhs); void switch_sec_mat(const uint32_t &sec_mat_); void fill_Weisse_table(); // check if translational symmetry satisfied void check_translation(); // naive way of enumerating all possible basis state void enumerate_basis_full(std::vector> conserve_lst = {}, std::vector val_lst = {}, const uint32_t &sec_full = 0); // Need to build Weiss Tables before enumerating representatives void enumerate_basis_repr(const std::vector &momentum, std::vector> conserve_lst = {}, std::vector val_lst = {}, const uint32_t &sec_repr = 0); // build the variational basis to run Trugman's method void build_basis_vrnl(const std::list &initial_list, const mbasis_elem &gs, const std::vector &momentum_gs, const std::vector &momentum, const uint32_t &iteration_depth, std::vector> conserve_lst = {}, std::vector val_lst = {}, const uint32_t &sec_vrnl = 0); // momentum has to be in format {m,n,...} corresponding to (m/L1) b_1 + (n/L2) b_2 + ... void basis_init_repr_deprecated(const std::vector &momentum, const uint32_t &sec_full = 0, const uint32_t &sec_repr = 0); // generate the Hamiltonian using basis_full void generate_Ham_sparse_full(const uint32_t &sec_full = 0, const bool &upper_triangle = true); // generate the Hamiltonian using basis_repr // a few artificial diagonal elements above 100, corresponding to zero norm states void generate_Ham_sparse_repr(const uint32_t &sec_repr = 0, const bool &upper_triangle = true); // generate the Hamiltonian using basis_vrnl void generate_Ham_sparse_vrnl(const uint32_t &sec_vrnl = 0, const bool &upper_triangle = true); // a few artificial diagonal elements above 100, corresponding to zero norm states void generate_Ham_sparse_repr_deprecated(const uint32_t &sec_full = 0, const uint32_t &sec_repr = 0, const bool &upper_triangle = true); // generate the Hamiltonian using basis_repr // generate a dense matrix of the Hamiltonian std::vector> to_dense(const uint32_t &sec_mat_ = 0); /** \brief y = H * x + y (matrix generated on the fly) */ void MultMv2(const T *x, T *y) const; /** \brief y = H * x (matrix generated on the fly) */ void MultMv(const T *x, T *y) const; // Note: in this function, (nev, ncv, maxit) have different meanings comparing to IRAM! // 1 <= nev <= 2, nev-1 <= ncv <= nev // nev = 1, calcualte up to ground state energy // nev = 2, calculate up to 1st excited state energy // ncv = 1, calculate up to ground state eigenvector // ncv = 2, calculate up to 1st excited state eigenvector void locate_E0_lanczos(const which_sym &sec_sym_, const MKL_INT &nev = 1, const MKL_INT &ncv = 1, MKL_INT maxit = 1000); /** \brief calculate the lowest eigenstates using IRAM * nev, ncv, maxit following ARPACK definition */ void locate_E0_iram(const which_sym &sec_sym_, const MKL_INT &nev = 2, const MKL_INT &ncv = 6, MKL_INT maxit = 0); /** \brief calculate the highest eigenstates using IRAM * nev, ncv, maxit following ARPACK definition. * for repr and vrnl, there may be a few artificial eigenvalues above fake_pos (default to 100), corresponding to zero norm states. */ void locate_Emax_iram(const which_sym &sec_sym_, const MKL_INT &nev = 2, const MKL_INT &ncv = 6, MKL_INT maxit = 0); /** \brief calculate eigenstates in [E_low, E_high]. Note: for repr and vrnl, there may be a few artificial eigenvalues above fake_pos (default to 100), corresponding to zero norm states. * \param nev number of estimated number of eigenstates in the requested energy range */ void locate_Es_feast(const which_sym &sec_sym_, const MKL_INT &nev, const double &E_low, const double &E_high, MKL_INT max_loop = 20); /** \brief return dim_full */ std::vector dimension_full() const { return dim_full; } /** \brief return dim_repr */ std::vector dimension_repr() const { return dim_repr; } /** \brief return dim_vrnl */ std::vector dimension_vrnl() const { return dim_vrnl; } /** \brief return E0 */ double energy_min() const { return E0; } /** \brief return Emax */ double energy_max() const { return Emax; } /** \brief return gap */ double energy_gap() const { return gap; } /** \brief lhs | phi >, where | phi > is an input state */ void moprXvec_full(const mopr &lhs, const uint32_t &sec_old, const uint32_t &sec_new, const T* vec_old, T* vec_new) const; /** \brief lhs | phi >, where | phi > is an eigenstate */ void moprXvec_full(const mopr &lhs, const uint32_t &sec_old, const uint32_t &sec_new, const MKL_INT &which_col, T* vec_new) const; /** transform a vector in the full space according to the plan */ void transform_vec_full(const std::vector &plan, const uint32_t &sec_full, const T* vec_old, T* vec_new) const; /** transform an eigenvector in the full space according to the plan */ void transform_vec_full(const std::vector &plan, const uint32_t &sec_full, const MKL_INT &which_col, T* vec_new) const; /** \brief project a state into a given momentum sector */ void projectQ_full(const std::vector &momentum, const uint32_t &sec_full, const T* vec_old, T* vec_new) const; /** \brief project an eigenstate into a given momentum sector */ void projectQ_full(const std::vector &momentum, const uint32_t &sec_full, const MKL_INT &which_col, T* vec_new) const; /** \brief < phi | lhs | phi > */ T measure_full_static(const mopr &lhs, const uint32_t &sec_full, const MKL_INT &which_col) const; /** \brief < phi | ... * lhs2 * lhs1 * lhs0 | phi >, where sec_old has to be given for each lhs_i */ T measure_full_static(const std::vector> &lhs, const std::vector &sec_old_list, const MKL_INT &which_col) const; /** \brief calculate dynamical structure factors * * \f[ * G_A(z) = \langle \phi | A_q^\dagger (z-H)^{-1} A_q | \phi \rangle, * \f] * where \f$ z = \omega + E_0 + i \eta \f$, and * norm = \f$ \sqrt{\langle \phi | A_q^\dagger A_q | \phi \rangle} \f$. * \f[ * G_A(z) = \frac{\langle \phi | A_q^\dagger A_q | \phi \rangle}{z-a_0 - * \frac{b_1^2}{z-a_1 - * \frac{b_2^2}{z-a_2 - \cdots}}} * \f] * * on exit: \f$ a_i \f$, \f$ b_i \f$ and norm are given. */ void measure_full_dynamic(const mopr &Aq, const uint32_t &sec_old, const uint32_t &sec_new, const MKL_INT &maxit, MKL_INT &m, double &norm, double hessenberg[]) const; /** \brief \f$ A_q | \phi \rangle \f$ * * Requirements: after operation of Aq, the new state is still an eigenstate of translation operator. * e.g. \f$ A_q = \frac{1}{\sqrt{N}} \sum_r e^{i q \cdot r} A_r \f$, * the old state has momentum k, the new state has momentum k+q */ void moprXvec_repr(const mopr &Aq, const uint32_t &sec_old, const uint32_t &sec_new, const T* vec_old, T* vec_new) const; /** \brief \f$ A_q | \phi \rangle \f$, where \f$ | \phi \rangle \f$ is an eigenstate */ void moprXvec_repr(const mopr &lhs, const uint32_t &sec_old, const uint32_t &sec_new, const MKL_INT &which_col, T* vec_new) const; void transform_vec_repr(const std::vector &plan, const uint32_t &sec_full, const MKL_INT &which_col, std::vector &vec_new) const; /** \brief < phi | lhs | phi > */ T measure_repr_static(const mopr &lhs, const uint32_t &sec_repr, const MKL_INT &which_col) const; /** \brief Calculate dynamical structure factors (in translational symmetric basis). * For details, see measure_full_dynamic. */ void measure_repr_dynamic(const mopr &Aq, const uint32_t &sec_old, const uint32_t &sec_new, const MKL_INT &maxit, MKL_INT &m, double &norm, double hessenberg[]) const; /** \f[ * A_q | G(Q_0) \rangle = \sum_i \frac{p_i}{N} | \varphi_i (Q_0 + q) \rangle, * \f] * where \f$ | G(Q_0) \rangle \f$ is the variational ground state. * * Note: in this input, will use * \f[ * B_q \equiv \sqrt{N} A_q = \sum_i A_{r_i} e^{i q \cdot r_i}, * \f] * since we cannot directly work with \f$ N \rightarrow \infty \f$. * (summation only over the box) * * Also, we output \f$ p_i \f$ instead of \f$ p_i/N \f$ for the same reason. */ void moprXgs_vrnl(const mopr &Bq, const uint32_t &sec_vrnl, T* vec_new) const; /** \f[ * A_q | vec \rangle = \sum_i \frac{p_i}{\sqrt{N}} | \varphi_i (k+q) \rangle + p_G | G(k+q) \rangle, * \f] * where \f$ | \varphi_k \rangle \f$ is a state with momentum k * * Note: in this input, will use * \f[ * B_q \equiv \sqrt{N} A_q = \sum_i A_{r_i} e^{i q \cdot r_i}, * \f] * since we cannot directly work with \f$ N \rightarrow \infty \f$. * (summation only over the box) * * Also, we output \f$ p_i \f$ instead of \f$ p_i/\sqrt{N} \f$ for the same reason. */ void moprXvec_vrnl(const mopr &Bq, const uint32_t &sec_old, const uint32_t &sec_new, const T* vec_old, T* vec_new, T &pG) const; void moprXvec_vrnl(const mopr &Bq, const uint32_t &sec_old, const uint32_t &sec_new, const MKL_INT &which_col, T* vec_new, T &pG) const; /** \brief < phi | lhs | phi >, where lhs is translational invariant */ T measure_vrnl_static_trans_invariant(const mopr &lhs, const uint32_t &sec_vrnl, const MKL_INT &which_col) const; /** Calculate dynamical structure factors (in Trugman basis). * For details, see measure_full_dynamic. * Note: the input \f$ B_q = \sqrt{N} A_q \f$ (summation only over the box) */ void measure_vrnl_dynamic(const mopr &Bq, const uint32_t &sec_vrnl, const MKL_INT &maxit, MKL_INT &m, double &norm, double hessenberg[]) const; /** Build the matrix for calculating the observables in the Wannier state * * \f[ * \mu_{k_1 k_2} = \langle \phi(k_1) | B_{k_1 - k_2} | \phi(k_2) \rangle , * \f] * where * \f[ * B_q = \sqrt{N} A_q = \sum_i e^{i q \cdot r_i} A_{r_i}. * \f] * Note: \f$ \phi(k_1) \f$ and \f$ \phi(k_2) \f$ should be in the same band (same quantum numbers and well-defined dispersion). * Thus, this function ONLY works if you really have a well defined band. * * Each pair denotes the position \f$ r_i \f$ (in cartesin coordinates) and the operator \f$ A_{r_i} \f$. */ void WannierMat_vrnl(const std::vector,mopr>> &Ar_list, const uint32_t &sec_vrnl, const std::vector> &momenta_list, std::vector> &matrix_mu, const std::function&, const uint32_t&)> &locate_state); // later add conserved quantum operators and corresponding quantum numbers? private: double Emax{}; double E0{}; double E1{}; double gap{}; double fake_pos{}; lattice latt_parent; lattice latt_sub; std::vector trans_sym; // if translation allowed in each dimension bool dim_spec_involved{}; // if translation allowed in partitioned direction // Weisse Tables for translation symmetry // Note: different from Weisse's paper, here we use Weisse_w to store the (parent) group label, instead of omega_g std::vector belong2rep_sub; std::vector> dist2rep_sub; std::vector>,uint32_t>> groups_parent; std::vector>,uint32_t>> groups_sub; std::vector omega_g_sub; std::vector belong2group_sub; MltArray_PairVec Weisse_e_lt, Weisse_e_eq, Weisse_e_gt; MltArray_uint32 Weisse_w_lt, Weisse_w_eq, Weisse_w_gt; void ckpt_lczsE0_init(bool &E0_done, bool &V0_done, bool &E1_done, bool &V1_done, std::vector &v); void ckpt_lczsE0_updt(const bool &E0_done, const bool &V0_done, const bool &E1_done, const bool &V1_done); }; // --------------------------- Kernel polynomial ---------------------------- // ---------------------------------------------------------------------------- /** \brief use Lanczos to determine the upper and lower bound of the eigenvalues of the matrix. */ template void energy_scale(const MKL_INT &dim, const MAT &mat, T v[], double &lo, double &hi, const double &extend = 0.1, const MKL_INT &iters = 128); // --------------------------- Miscellaneous stuff ---------------------------- // ---------------------------------------------------------------------------- std::string date_and_time(); inline int round2int(const double &x) { return static_cast(x < 0.0 ? x - 0.5 : x + 0.5); } inline double conjugate(const double &rhs) { return rhs; } inline std::complex conjugate(const std::complex &rhs) { return std::conj(rhs); } // calculate base^index, in the case both are integers template T2 int_pow(const T1 &base, const T1 &index); // given two arrays: num & base, get the result of: // num[0] + num[1] * base[0] + num[2] * base[0] * base[1] + num[3] * base[0] * base[1] * base[2] + ... // e.g. (base = {2,2,2,2,2}) // 0 0 1 0 1 // ^ ^ // | | // num[1] num[0] // 1 + 0 * 2 + 1 * 2^2 + 0 * 2^3 + 0 * 2^4 template T2 dynamic_base(const std::vector &nums, const std::vector &base); // deprecated template void dynamic_base_vec2num(const MKL_INT &len, const T1* base, const T1* vec, T2 &num); // the other way around template std::vector dynamic_base(const T2 &total, const std::vector &base); // nums + 1 template std::vector dynamic_base_plus1(const std::vector &nums, const std::vector &base); // check if maximized template bool dynamic_base_maximized(const std::vector &nums, const std::vector &base); // check overflow template bool dynamic_base_overflow(const std::vector &nums, const std::vector &base); template bool is_sorted_norepeat(const std::vector &array); // note: end means the position which has already passed the last element template T2 binary_search(const std::vector &array, const T1 &val, const T2 &bgn, const T2 &end); // return the number of exchanges happened during the bubble sort template int bubble_sort(std::vector &array, const int &bgn, const int &end); // b1 // a0 + --------------- // b2 // a1 + ---------- // a2 + ... template T continued_fraction(T a[], T b[], const MKL_INT &len); // b0 not used template void vec_swap(const MKL_INT &n, T *x, T *y); /** \beief fill x with zeros */ template void vec_zeros(const MKL_INT &n, T *x); /** \brief fill x with Lehmer 16807 random numbers, * where seed == 0 reserved for filling 1/n to each element */ template void vec_randomize(const MKL_INT &n, T *x, const uint32_t &seed = 1); template int vec_disk_read(const std::string &filename, MKL_INT n, T *x); template int vec_disk_write(const std::string &filename, MKL_INT n, T *x); int basis_disk_read(const std::string &filename, std::vector &basis); int basis_disk_write(const std::string &filename, const std::vector &basis); } #endif /* qbasis_h */