Skip to content

Commit

Permalink
sany: spin-dependent sgf
Browse files Browse the repository at this point in the history
  • Loading branch information
@hczhai
hczhai committed Jun 23, 2024
1 parent f987afc commit 9056137
Show file tree
Hide file tree
Showing 3 changed files with 173 additions and 21 deletions.
175 changes: 156 additions & 19 deletions pyblock2/driver/core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -333,6 +333,7 @@ def __init__(self, symm_type=SymmetryTypes.SU2):
self.b = b
self.symm_type = symm_type
self.qargs = None
self.hint = None
has_cpx = hasattr(b, "cpx")
has_sp = hasattr(b, "sp")
has_spcpx = has_sp and hasattr(b.sp, "cpx")
Expand Down Expand Up @@ -503,7 +504,7 @@ def init_phsu2(n, twos, nh=None):
self.VectorSX = b.VectorSGB
self.VectorVectorSX = b.VectorVectorSGB

def set_symmetry_groups(self, *args):
def set_symmetry_groups(self, *args, hint=None):
"""
Set the combination of symmetry sub-groups for ``symm_type = SAny``.
Expand All @@ -512,6 +513,8 @@ def set_symmetry_groups(self, *args):
List of names of (Abelian) symmetry groups. ``0 <= len(args) <= 6`` is required.
Possible sub-group names are "U1", "Z1", "Z2", "Z3", ..., "Z2055",
"U1Fermi", "Z1Fermi", "Z2Fermi", "Z3Fermi", ..., "Z2055Fermi", "LZ", and "AbelianPG".
hint : str or None
Hint for symmetry interpretation. Default is None.
"""
assert self.SXT == self.b.SAny and len(args) <= 6

Expand All @@ -531,6 +534,7 @@ def init_sany(*qargs):
return q

self.qargs = args
self.hint = hint
self.SX = init_sany


Expand Down Expand Up @@ -734,7 +738,7 @@ def set_symm_type(self, symm_type, reset_frame=True):
self.mpi.barrier()
self.frame.restart_dir = self.restart_dir

def set_symmetry_groups(self, *args):
def set_symmetry_groups(self, *args, hint=None):
"""
Set the combination of symmetry sub-groups for ``symm_type = SAny``.
Expand All @@ -743,8 +747,10 @@ def set_symmetry_groups(self, *args):
List of names of (Abelian) symmetry groups. ``0 <= len(args) <= 6`` is required.
Possible sub-group names are "U1", "Z1", "Z2", "Z3", ..., "Z2055",
"U1Fermi", "Z1Fermi", "Z2Fermi", "Z3Fermi", ..., "Z2055Fermi", "LZ", and "AbelianPG".
hint : str or None
Hint for symmetry interpretation. Default is None.
"""
self.bw.set_symmetry_groups(*args)
self.bw.set_symmetry_groups(*args, hint=hint)

@property
def basis(self):
Expand Down Expand Up @@ -932,7 +938,20 @@ def initialize_system(
[[0, 0, 1, 0], [0, 0, 0, -1], [0, 0, 0, 0], [0, 0, 0, 0]]
), # beta
}
if bw.qargs == ("U1Fermi", "AbelianPG"):
std_ops_sgf = {
"": np.array([[1, 0], [0, 1]]), # identity
"C": np.array([[0, 0], [1, 0]]), # +
"D": np.array([[0, 1], [0, 0]]), # -
}
if bw.qargs == ("U1Fermi", "AbelianPG") and bw.hint == "SGF":
site_basis, site_ops = [], []
for k in range(self.n_sites):
ipg = self.orb_sym[k]
basis = [(self.bw.SX(0, 0), 1), (self.bw.SX(1, ipg), 1)] # [0a]
site_basis.append(basis)
site_ops.append(std_ops_sgf)
self.ghamil = self.get_custom_hamiltonian(site_basis, site_ops)
elif bw.qargs == ("U1Fermi", "AbelianPG"):
site_basis, site_ops = [], []
for k in range(self.n_sites):
ipg = self.orb_sym[k]
Expand All @@ -944,6 +963,21 @@ def initialize_system(
site_basis.append(basis)
site_ops.append(std_ops)
self.ghamil = self.get_custom_hamiltonian(site_basis, site_ops)
elif bw.qargs == ("U1Fermi", "U1", "AbelianPG") and bw.hint == "SGF":
site_basis, site_ops = [], []
spin_sym = []
for k in range(self.n_sites):
ipg = self.orb_sym[k]
basis = [
(self.bw.SX(0, 0, 0), 1),
(self.bw.SX(1, 1 - (k % 2) * 2, ipg), 1),
] # [0a]
site_basis.append(basis)
site_ops.append(std_ops_sgf)
spin_sym.append({"C": 1 - (k % 2) * 2, "D": 1 - (1 - k % 2) * 2})
self.ghamil = self.get_custom_hamiltonian(
site_basis, site_ops, spin_dependent_ops="CD", spin_sym=spin_sym
)
elif bw.qargs == ("U1Fermi", "U1", "AbelianPG"):
site_basis, site_ops = [], []
for k in range(self.n_sites):
Expand Down Expand Up @@ -2226,7 +2260,14 @@ def deallocate(self):

return LZHamiltonian(self.vacuum, self.n_sites, self.orb_sym)

def get_custom_hamiltonian(self, site_basis, site_ops, orb_dependent_ops="cdCD"):
def get_custom_hamiltonian(
self,
site_basis,
site_ops,
orb_dependent_ops="cdCD",
spin_dependent_ops="",
spin_sym=None,
):
"""
Setting Hamiltonian in the general symmetry mode. ``SZ`` or ``SAny`` symmetry mode is required.
Expand All @@ -2246,6 +2287,10 @@ def get_custom_hamiltonian(self, site_basis, site_ops, orb_dependent_ops="cdCD")
List of operator names that can have point group irrep.
If point group or ``orb_sym`` is not used, this can be empty.
Default is "cdCD".
spin_dependent_ops : str
List of operator names that can have different spin in different sites. Default is empty.
spin_sym : None or list[dict[str, int]]
List of spin symmetries for site operators. Default is None.
Returns:
ghamil : CustomHamiltonian
Expand All @@ -2262,14 +2307,15 @@ def get_custom_hamiltonian(self, site_basis, site_ops, orb_dependent_ops="cdCD")
):

class CustomHamiltonian(GH):
def __init__(self, vacuum, n_sites, orb_sym):
def __init__(self, vacuum, n_sites, orb_sym, spin_sym=None):
GH.__init__(self)
self.opf = super_self.bw.bs.OperatorFunctions(
super_self.bw.brs.CG()
)
self.vacuum = vacuum
self.n_sites = n_sites
self.orb_sym = orb_sym
self.spin_sym = spin_sym
self.basis = super_self.bw.brs.VectorStateInfo(
[self.get_site_basis(m) for m in range(self.n_sites)]
)
Expand Down Expand Up @@ -2393,7 +2439,9 @@ def get_site_string_ops(self, m, ops):
xp = self.site_norm_ops[m][p]
q = xx.info.delta_quantum + xp.info.delta_quantum
mat = super_self.bw.bs.SparseMatrix(d_alloc)
mat.allocate(self.find_site_op_info(m, q))
dq = self.find_site_op_info(m, q)
assert dq is not None
mat.allocate(dq)
self.opf.product(0, xx, xp, mat)
xx = mat
ops[k] = self.site_norm_ops[m][k] = xx
Expand All @@ -2405,8 +2453,12 @@ def init_string_quanta(self, exprs, term_l, left_vacuum):
for norm_ops in self.site_norm_ops:
for k, v in norm_ops.items():
if k not in qs:
qs[k] = v.info.delta_quantum
qs[k].pg = 0
# must copy to prevent changing v.info.dq
qs[k] = v.info.delta_quantum[0]
if k in orb_dependent_ops:
qs[k].pg = 0
if k in spin_dependent_ops:
qs[k].twos = 0
return super_self.bw.VectorVectorSX(
[
super_self.bw.VectorSX(
Expand All @@ -2425,13 +2477,19 @@ def get_string_quanta(self, ref, expr, idxs, k):
"""Quantum number for string operators (orbital dependent part)."""
l, r = ref[k], ref[-1] - ref[k]
for j, (ex, ix) in enumerate(zip(expr, idxs)):
ipg = self.orb_sym[ix]
if ex not in orb_dependent_ops:
pass
elif j < k:
l.pg = l.pg ^ ipg
else:
r.pg = r.pg ^ ipg
if ex in orb_dependent_ops:
ipg = self.orb_sym[ix]
if j < k:
l.pg = l.pg ^ ipg
else:
r.pg = r.pg ^ ipg
if ex in spin_dependent_ops:
assert self.spin_sym is not None
ispin = self.spin_sym[ix][ex]
if j < k:
l.twos = l.twos + ispin
else:
r.twos = r.twos + ispin
return l, r

def get_string_quantum(self, expr, idxs):
Expand All @@ -2453,7 +2511,9 @@ def deallocate(self):
for bz in self.basis:
bz.deallocate()

return CustomHamiltonian(self.vacuum, self.n_sites, self.orb_sym)
return CustomHamiltonian(
self.vacuum, self.n_sites, self.orb_sym, spin_sym=spin_sym
)
else:
return NotImplemented

Expand Down Expand Up @@ -6994,6 +7054,13 @@ def get_mps_from_csf_coefficients(
mps = self.adjust_mps(mps, dot=dot)[0]
return mps

def get_spin_projection_npts(self, n_sites, n_elec, twos):
import numpy as np
if n_elec <= n_sites:
return int(np.ceil((n_elec / 2 + twos / 2 + 1) / 2))
else:
return int(np.ceil(((2 * n_sites - n_elec) / 2 + twos / 2 + 1) / 2))

def get_spin_projection_mpo(
self,
twos,
Expand Down Expand Up @@ -7048,7 +7115,10 @@ def get_spin_projection_mpo(
it = np.ones((1, 1, 1, 1))
pympo = None
for ixw, (xt, wt) in enumerate(zip(xts, wts)):
if not mpi_split or (mpi_split and self.mpi.rank == min(ixw, len(wts) - 1 - ixw) % self.mpi.size):
if not mpi_split or (
mpi_split
and self.mpi.rank == min(ixw, len(wts) - 1 - ixw) % self.mpi.size
):
ct = np.cos(xt / 2) * it
st, mt = np.sin(xt / 2) * it, -np.sin(xt / 2) * it
if SymmetryTypes.SZ in self.symm_type and use_sz_symm:
Expand Down Expand Up @@ -7105,7 +7175,74 @@ def get_spin_projection_mpo(
]
tensors.append(Tensor(blocks=blocks))
lqs = rqs
elif SymmetryTypes.SAny in self.symm_type:
elif (
SymmetryTypes.SAny in self.symm_type
and bw.qargs == ("U1Fermi", "AbelianPG")
and bw.hint == "SGF"
):
q, tensors = bw.SX(), []
for k, bz in enumerate(self.basis):
blocks = []
xq, yq = sorted(bz, key=lambda xq: xq.n)
if k % 2 == 0:
blocks.append(SubTensor(
reduced=np.array([1, -1]).reshape((1, 1, 1, 2)),
q_labels=(q, xq, xq, q),
))
blocks.append(SubTensor(
reduced=np.array([1, 1]).reshape((1, 1, 1, 2)),
q_labels=(q, yq, yq, q),
))
blocks.append(SubTensor(
reduced=np.array([1]).reshape((1, 1, 1, 1)),
q_labels=(q, yq, xq, (q + yq - xq)[0]),
))
blocks.append(SubTensor(
reduced=np.array([1]).reshape((1, 1, 1, 1)),
q_labels=(q, xq, yq, (q + xq - yq)[0]),
))
else:
cp, cm, st = (
(1 + np.cos(xt / 2)) / 2,
(1 - np.cos(xt / 2)) / 2,
np.sin(xt / 2),
)
blocks.append(SubTensor(
reduced=np.array([cp, -cm]).reshape((2, 1, 1, 1)),
q_labels=(q, xq, xq, q),
))
blocks.append(SubTensor(
reduced=np.array([cp, cm]).reshape((2, 1, 1, 1)),
q_labels=(q, yq, yq, q),
))
blocks.append(SubTensor(
reduced=np.array([st]).reshape((1, 1, 1, 1)),
q_labels=((q + yq - xq)[0], xq, yq, q),
))
blocks.append(SubTensor(
reduced=np.array([st]).reshape((1, 1, 1, 1)),
q_labels=((q + xq - yq)[0], yq, xq, q),
))
if k == 0:
blocks = [
SubTensor(
reduced=wt * x.reduced[0, ...],
q_labels=x.q_labels[1:],
)
for x in blocks
]
elif k == n_sites - 1:
blocks = [
SubTensor(
reduced=x.reduced[..., 0], q_labels=x.q_labels[:-1]
)
for x in blocks
]
tensors.append(Tensor(blocks=blocks))
elif SymmetryTypes.SAny in self.symm_type and bw.qargs == (
"U1Fermi",
"AbelianPG",
):
rt = np.array(
[
[np.cos(xt / 2), np.sin(xt / 2)],
Expand Down
18 changes: 16 additions & 2 deletions src/core/symmetry.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -147,12 +147,18 @@ template <int8_t L = 6> struct SAnyT {
return (int)values[k];
return 0;
}
int u1() const {
for (int8_t k = 0; k < L; k++)
if (types[k] == SAnySymmTypes::U1)
return (int)values[k];
return 0;
}
int twos() const {
for (int8_t k = 0; k < L; k++)
if (types[k] == SAnySymmTypes::SU2Fermi ||
types[k] == SAnySymmTypes::SU2)
return (int)values[k + 1];
return 0;
return u1();
}
int twos_low() const {
for (int8_t k = 0; k < L; k++)
Expand All @@ -175,13 +181,21 @@ template <int8_t L = 6> struct SAnyT {
break;
}
}
void set_u1(int n) {
for (int8_t k = 0; k < L; k++)
if (types[k] == SAnySymmTypes::U1) {
values[k] = n;
break;
}
}
void set_twos(int twos) {
for (int8_t k = 0; k < L; k++)
if (types[k] == SAnySymmTypes::SU2Fermi ||
types[k] == SAnySymmTypes::SU2) {
values[++k] = twos;
break;
return;
}
set_u1(twos);
}
void set_twos_low(int twos) {
for (int8_t k = 0; k < L; k++)
Expand Down
1 change: 1 addition & 0 deletions src/pybind/pybind_core.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3934,6 +3934,7 @@ template <typename S = void> void bind_symmetry(py::module &m) {
.def_static("init_sgb", &SAny::init_sgb, py::arg("n") = 0,
py::arg("pg") = 0)
.def_property("n", &SAny::n, &SAny::set_n)
.def_property("u1", &SAny::u1, &SAny::set_u1)
.def_property("twos", &SAny::twos, &SAny::set_twos)
.def_property("twos_low", &SAny::twos_low, &SAny::set_twos_low)
.def_property("pg", &SAny::pg, &SAny::set_pg)
Expand Down

0 comments on commit 9056137

Please sign in to comment.