Skip to content

Commit

Permalink
get_npdm flexible mask
Browse files Browse the repository at this point in the history
  • Loading branch information
@hczhai
hczhai committed Sep 2, 2024
1 parent cdd0f99 commit e696b80
Show file tree
Hide file tree
Showing 6 changed files with 148 additions and 14 deletions.
2 changes: 1 addition & 1 deletion docs/source/user/interfaces.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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`.

Expand Down
56 changes: 55 additions & 1 deletion pyblock2/driver/core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]

Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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]]
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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
]
)
Expand All @@ -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)
]
Expand All @@ -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)
]
Expand All @@ -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
]
)
Expand All @@ -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)
]
Expand All @@ -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)
Expand Down
69 changes: 65 additions & 4 deletions src/core/spin_permutation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<string, int> auto_sort_string(const vector<uint16_t> &x,
const string &xops,
const string &fermionic_ops) {
vector<uint16_t> 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<double> dot_product(const SpinPermTensor &a,
const SpinPermTensor &b) {
assert(a.data.size() == b.data.size());
Expand Down Expand Up @@ -929,6 +945,7 @@ struct SpinPermPattern {
const vector<uint16_t> &mask = vector<uint16_t>()) {
if (x.size() == 0)
return x;
bool mask_sorted = is_sorted(mask.begin(), mask.end());
vector<pair<uint16_t, uint16_t>> pp;
for (auto &ix : x)
if (pp.size() == 0 || ix != pp.back().first)
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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<uint16_t> &mask = vector<uint16_t>()) {
using T = SpinPermTensor;
using R = SpinPermRecoupling;
Expand Down Expand Up @@ -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<uint16_t> &mask = vector<uint16_t>()) {
using T = SpinPermTensor;
using R = SpinPermRecoupling;
SpinPermPattern spat(nn, mask);
vector<double> 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<uint16_t> irr = spat[i];
r.index_patterns[k] = irr;
vector<uint16_t> 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<uint16_t> indices(rr.begin() + jj * irr.size(),
rr.begin() + (jj + 1) * irr.size());
vector<uint16_t> perm =
SpinPermTensor::find_pattern_perm(indices);
r.data[k][perm] = vector<pair<double, string>>();
vector<pair<double, string>> &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;
Expand Down
10 changes: 7 additions & 3 deletions src/core/tensor_functions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1020,9 +1020,13 @@ template <typename S, typename FL> struct TensorFunctions {
vector<uint64_t> 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
Expand Down
12 changes: 8 additions & 4 deletions src/dmrg/sweep_algorithm.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<MKL_INT> shape(n_op, n_physical_sites);
r[i] = make_shared<GTensor<FLX>>(shape);
Expand Down
13 changes: 12 additions & 1 deletion src/pybind/pybind_core.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2246,7 +2246,14 @@ template <typename S = void> 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<string, int>(*)(const vector<uint16_t> &,
const string &, const string &)) &
SpinPermTensor::auto_sort_string)
.def_static(
"auto_sort_string",
(pair<string, int>(*)(const vector<uint16_t> &, const string &)) &
SpinPermTensor::auto_sort_string)
.def_static("mul", (SpinPermTensor(*)(const SpinPermTensor &,
const SpinPermTensor &, int16_t,
const SU2CG &)) &
Expand Down Expand Up @@ -2360,6 +2367,10 @@ template <typename S = void> void bind_io(py::module &m) {
py::arg("nn"), py::arg("spin_str"),
py::arg("is_fermion") = true,
py::arg("mask") = vector<uint16_t>())
.def_static("initialize_sany", &SpinPermScheme::initialize_sany,
py::arg("nn"), py::arg("spin_str"),
py::arg("fermionic_ops") = string("cdCD"),
py::arg("mask") = vector<uint16_t>())
.def_static("initialize_su2_old", &SpinPermScheme::initialize_su2_old,
py::arg("nn"), py::arg("spin_str"),
py::arg("is_npdm") = false)
Expand Down

0 comments on commit e696b80

Please sign in to comment.