diff --git a/docs/source/user/interfaces.rst b/docs/source/user/interfaces.rst index e03c4b1e..dd4c8d19 100644 --- a/docs/source/user/interfaces.rst +++ b/docs/source/user/interfaces.rst @@ -7,7 +7,7 @@ Interfaces Input File ---------- -Like many quantum chemistry packages, ``block2`` can be used by reading parameters and instrutions from a formatted input file. +Like many quantum chemistry packages, ``block2`` can be used by reading parameters and instructions from a formatted input file. This interface is ``StackBlock`` compatible. See :ref:`user_basic`, :ref:`user_advanced`, and :ref:`user_keywords`. diff --git a/pyblock2/driver/core.py b/pyblock2/driver/core.py index 164d62d8..d733c24d 100644 --- a/pyblock2/driver/core.py +++ b/pyblock2/driver/core.py @@ -2541,7 +2541,8 @@ def init_site_ops(self): ] blocks.append((q_map[i][1], q_map[j][1], mat)) - assert dqs is not None + if dqs is None: + dqs = [self.vacuum] assert len(dqs) >= 1 dq = dqs[0] @@ -5418,6 +5419,7 @@ def get_npdm( cutoff=1e-24, iprint=0, max_bond_dim=None, + fermionic_ops=None, ): """ Compute the N-Particle Density Matrix (NPDM) for the given MPS. @@ -5473,6 +5475,11 @@ def get_npdm( max_bond_dim : None or int If not None, the MPS bond dimension will be restricted during the sweeps. Default is None. + fermionic_ops : None or str + If not None, the given set of operator names will be treated as Fermion + operators (for computing signs for swapping operators). + Default is None, and operators like "cdCD" will be treated as Fermion + operators. Returns: dms : np.ndarray[flat|complex] or list[np.ndarray[flat|complex]] @@ -5529,6 +5536,7 @@ def get_npdm( self.mpi.barrier() if SymmetryTypes.SU2 in bw.symm_type: + assert fermionic_ops is None if npdm_expr is not None and "%s" not in npdm_expr: op_str = npdm_expr else: @@ -5556,6 +5564,8 @@ def get_npdm( perms = bw.b.VectorSpinPermScheme( [ bw.b.SpinPermScheme.initialize_sz(pdm_type * 2, cd, True) + if fermionic_ops is None else + bw.b.SpinPermScheme.initialize_sany(pdm_type * 2, cd, fermionic_ops) for cd in op_str ] ) @@ -5568,6 +5578,9 @@ def get_npdm( [ bw.b.SpinPermScheme.initialize_sz( pt * 2, cd, True, mask=bw.b.VectorUInt16(xm) + ) if fermionic_ops is None else + bw.b.SpinPermScheme.initialize_sany( + pt * 2, cd, fermionic_ops, mask=bw.b.VectorUInt16(xm) ) for cd, xm, pt in zip(op_str, mask, pts) ] @@ -5580,6 +5593,9 @@ def get_npdm( [ bw.b.SpinPermScheme.initialize_sz( pt * 2, cd, True, mask=bw.b.VectorUInt16(mask) + ) if fermionic_ops is None else + bw.b.SpinPermScheme.initialize_sany( + pt * 2, cd, fermionic_ops, mask=bw.b.VectorUInt16(mask) ) for cd, pt in zip(op_str, pts) ] @@ -5595,6 +5611,8 @@ def get_npdm( perms = bw.b.VectorSpinPermScheme( [ bw.b.SpinPermScheme.initialize_sz(pdm_type * 2, cd, True) + if fermionic_ops is None else + bw.b.SpinPermScheme.initialize_sany(pdm_type * 2, cd, fermionic_ops) for cd in op_str ] ) @@ -5607,6 +5625,9 @@ def get_npdm( [ bw.b.SpinPermScheme.initialize_sz( pt * 2, cd, True, mask=bw.b.VectorUInt16(xm) + ) if fermionic_ops is None else + bw.b.SpinPermScheme.initialize_sany( + pt * 2, cd, fermionic_ops, mask=bw.b.VectorUInt16(xm) ) for cd, xm, pt in zip(op_str, mask, pts) ] @@ -5619,10 +5640,43 @@ def get_npdm( [ bw.b.SpinPermScheme.initialize_sz( pt * 2, cd, True, mask=bw.b.VectorUInt16(mask) + ) if fermionic_ops is None else + bw.b.SpinPermScheme.initialize_sany( + pt * 2, cd, fermionic_ops, mask=bw.b.VectorUInt16(mask) ) for cd, pt in zip(op_str, pts) ] ) + elif SymmetryTypes.SAny in bw.symm_type: + assert npdm_expr is not None + op_str = [npdm_expr] if isinstance(npdm_expr, str) else npdm_expr + fermionic_ops = "cdCD" if fermionic_ops is None else fermionic_ops + if mask is None: + perms = bw.b.VectorSpinPermScheme( + [ + bw.b.SpinPermScheme.initialize_sany(len(cd), cd, fermionic_ops) + for cd in op_str + ] + ) + elif len(mask) != 0 and not isinstance(mask[0], int): + assert len(mask) == len(op_str) + perms = bw.b.VectorSpinPermScheme( + [ + bw.b.SpinPermScheme.initialize_sany( + len(cd), cd, fermionic_ops, mask=bw.b.VectorUInt16(xm) + ) + for cd, xm in zip(op_str, mask) + ] + ) + else: + perms = bw.b.VectorSpinPermScheme( + [ + bw.b.SpinPermScheme.initialize_sany( + len(cd), cd, fermionic_ops, mask=bw.b.VectorUInt16(mask) + ) + for cd in op_str + ] + ) if iprint >= 1: print("npdm string =", op_str) diff --git a/src/core/spin_permutation.hpp b/src/core/spin_permutation.hpp index ddb6ac25..16531f7b 100644 --- a/src/core/spin_permutation.hpp +++ b/src/core/spin_permutation.hpp @@ -221,6 +221,22 @@ struct SpinPermTensor { z[pcnt[x[i]]] = xops[i], perm[i] = pcnt[x[i]]++; return make_pair(z, permutation_parity(perm) ? -1 : 1); } + static pair auto_sort_string(const vector &x, + const string &xops, + const string &fermionic_ops) { + vector px = x; + string z = xops; + int xsign = 1; + for (int l = (int)px.size() - 1; l > 0; l--) + for (int p = 0; p < l; p++) + if (px[p] > px[p + 1]) { + swap(px[p], px[p + 1]), swap(z[p], z[p + 1]); + if (fermionic_ops.find(z[p]) != string::npos && + fermionic_ops.find(z[p + 1]) != string::npos) + xsign *= -1; + } + return make_pair(z, xsign); + } static vector dot_product(const SpinPermTensor &a, const SpinPermTensor &b) { assert(a.data.size() == b.data.size()); @@ -929,6 +945,7 @@ struct SpinPermPattern { const vector &mask = vector()) { if (x.size() == 0) return x; + bool mask_sorted = is_sorted(mask.begin(), mask.end()); vector> pp; for (auto &ix : x) if (pp.size() == 0 || ix != pp.back().first) @@ -971,9 +988,16 @@ struct SpinPermPattern { r[ir + k] = g[j + kk++]; if (mask.size() != 0) { bool skip = false; - for (int k = 1; k < (int)x.size(); k++) - skip = skip || (mask[k] == mask[k - 1] && - r[ir + k] != r[ir + k - 1]); + if (mask_sorted) { + for (int k = 1; k < (int)x.size(); k++) + skip = skip || (mask[k] == mask[k - 1] && + r[ir + k] != r[ir + k - 1]); + } else { + for (int k = 0; k < (int)x.size(); k++) + for (int l = k + 1; l < (int)x.size(); l++) + skip = skip || (mask[k] == mask[l] && + r[ir + k] != r[ir + l]); + } if (skip) continue; } @@ -1179,7 +1203,7 @@ struct SpinPermScheme { left_vacuum = r.left_vacuum; } static SpinPermScheme - initialize_sz(int nn, string spin_str, bool is_fermion = true, + initialize_sz(int nn, const string &spin_str, bool is_fermion = true, const vector &mask = vector()) { using T = SpinPermTensor; using R = SpinPermRecoupling; @@ -1214,6 +1238,43 @@ struct SpinPermScheme { r.mask = mask; return r; } + static SpinPermScheme + initialize_sany(int nn, const string &spin_str, + const string &fermionic_ops = "cdCD", + const vector &mask = vector()) { + using T = SpinPermTensor; + using R = SpinPermRecoupling; + SpinPermPattern spat(nn, mask); + vector mptr; + SpinPermScheme r; + r.index_patterns.resize(spat.count()); + r.data.resize(spat.count()); + size_t k = 0; + for (size_t i = 0; i < spat.count(); i++) { + vector irr = spat[i]; + r.index_patterns[k] = irr; + vector rr = SpinPermPattern::all_reordering(irr, mask); + int nj = irr.size() == 0 ? 1 : (int)(rr.size() / irr.size()); + for (int jj = 0; jj < nj; jj++) { + vector indices(rr.begin() + jj * irr.size(), + rr.begin() + (jj + 1) * irr.size()); + vector perm = + SpinPermTensor::find_pattern_perm(indices); + r.data[k][perm] = vector>(); + vector> &rec_formula = r.data[k].at(perm); + auto pis = SpinPermTensor::auto_sort_string(indices, spin_str, + fermionic_ops); + rec_formula.push_back(make_pair((double)pis.second, pis.first)); + } + k += nj != 0; + } + r.index_patterns.resize(k); + r.data.resize(k); + r.is_su2 = false; + r.left_vacuum = 0; + r.mask = mask; + return r; + } static SpinPermScheme initialize_su2_old(int nn, string spin_str, bool is_npdm = false) { using T = SpinPermTensor; diff --git a/src/core/tensor_functions.hpp b/src/core/tensor_functions.hpp index b652bc8d..ac8d79eb 100644 --- a/src/core/tensor_functions.hpp +++ b/src/core/tensor_functions.hpp @@ -1020,9 +1020,13 @@ template struct TensorFunctions { vector rmx(rpat.size(), 0); uint64_t mxx = r_step; for (int k = (int)perm.size() - 1; k >= 0; k--) { - if (mask.size() != 0 && k != 0 && - mask[k] == mask[k - 1]) - continue; + if (mask.size() != 0) { + bool ok = true; + for (int kk = 0; kk < k; kk++) + ok = ok && mask[k] != mask[kk]; + if (!ok) + continue; + } if (perm[k] < lmx.size()) lmx[perm[k]] = mxx; else diff --git a/src/dmrg/sweep_algorithm.hpp b/src/dmrg/sweep_algorithm.hpp index 82ea0d04..8ac7a23d 100644 --- a/src/dmrg/sweep_algorithm.hpp +++ b/src/dmrg/sweep_algorithm.hpp @@ -6541,10 +6541,14 @@ struct Expect { for (int i = 0; i < (int)scheme->perms.size(); i++) { int n_op = (int)scheme->perms[i]->index_patterns[0].size(); if (scheme->perms[i]->mask.size() != 0) { - n_op = 1; - for (int k = 1; k < scheme->perms[i]->mask.size(); k++) - n_op += scheme->perms[i]->mask[k] != - scheme->perms[i]->mask[k - 1]; + n_op = 0; + for (int k = 0; k < (int)scheme->perms[i]->mask.size(); k++) { + bool ok = true; + for (int j = 0; j < k; j++) + ok = ok && scheme->perms[i]->mask[k] != + scheme->perms[i]->mask[j]; + n_op += (int)ok; + } } vector shape(n_op, n_physical_sites); r[i] = make_shared>(shape); diff --git a/src/pybind/pybind_core.hpp b/src/pybind/pybind_core.hpp index 5f39912d..6d61da3e 100644 --- a/src/pybind/pybind_core.hpp +++ b/src/pybind/pybind_core.hpp @@ -2246,7 +2246,14 @@ template void bind_io(py::module &m) { .def_static("T", &SpinPermTensor::T) .def_static("permutation_parity", &SpinPermTensor::permutation_parity) .def_static("find_pattern_perm", &SpinPermTensor::find_pattern_perm) - .def_static("auto_sort_string", &SpinPermTensor::auto_sort_string) + .def_static("auto_sort_string", + (pair(*)(const vector &, + const string &, const string &)) & + SpinPermTensor::auto_sort_string) + .def_static( + "auto_sort_string", + (pair(*)(const vector &, const string &)) & + SpinPermTensor::auto_sort_string) .def_static("mul", (SpinPermTensor(*)(const SpinPermTensor &, const SpinPermTensor &, int16_t, const SU2CG &)) & @@ -2360,6 +2367,10 @@ template void bind_io(py::module &m) { py::arg("nn"), py::arg("spin_str"), py::arg("is_fermion") = true, py::arg("mask") = vector()) + .def_static("initialize_sany", &SpinPermScheme::initialize_sany, + py::arg("nn"), py::arg("spin_str"), + py::arg("fermionic_ops") = string("cdCD"), + py::arg("mask") = vector()) .def_static("initialize_su2_old", &SpinPermScheme::initialize_su2_old, py::arg("nn"), py::arg("spin_str"), py::arg("is_npdm") = false)