Skip to content

Commit 3e35fe7

Browse files
hczhaihczhai
committed
docs: spin-projected dmrg
1 parent a9cc0eb commit 3e35fe7

File tree

3 files changed

+264
-4
lines changed

3 files changed

+264
-4
lines changed

‎docs/source/index.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -152,6 +152,7 @@ Features
152152
tutorial/vibrational-hamiltonians
153153
tutorial/mpo-mps-quimb
154154
tutorial/mps-import-renormalizer
155+
tutorial/spin-projected-dmrg
155156

156157
.. raw:: latex
157158

‎docs/source/tutorial/spin-projected-dmrg.rst‎

Lines changed: 254 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,254 @@
1+
2+
.. _tutorial_spin_projected_dmrg:
3+
4+
Spin-Projected DMRG
5+
===================
6+
7+
Spin-projected DMRG (SP-DMRG) is a powerful technique for generating reliable initial guess
8+
Matrix Product States (MPS) for spin-adapted DMRG, particularly in systems with numerous
9+
competing broken-symmetry states. Due to its high computational cost, SP-DMRG is typically
10+
performed only at small bond dimensions. The resulting optimized MPS can then serve as a
11+
qualitatively reliable initial guess for subsequent, larger-scale optimization using
12+
spin-adapted DMRG (under SU2 symmetry).
13+
14+
Reference for the spin-projected DMRG algorithm:
15+
16+
* Li, Z., Chan, G. K.-L. Spin-Projected Matrix Product States: Versatile Tool for Strongly Correlated Systems. *Journal of Chemical Theory and Computation* 2017, **13**, 2681-2695. doi: `10.1021/acs.jctc.7b00270 <https://doi.org/10.1021/acs.jctc.7b00270>`_
17+
18+
The following example shows how to use spin-projected DMRG to generate the initial guess MPS.
19+
We study the three broken-symmetry states of the Fe4S4 active space model. The integral file
20+
can be found using
21+
22+
.. code-block:: bash
23+
wget -O Fe4S4.FCIDUMP https://raw.githubusercontent.com/zhendongli2008/Active-space-model-for-Iron-Sulfur-Clusters/main/Fe2S2_and_Fe4S4/Fe4S4/fe4s4
24+
25+
Exact MPO
26+
---------
27+
28+
In the first example, we use an exact MPO for the Hamiltonian, this can be done directly
29+
in the particle-number U1 symmetry mode. MPS can be initialized using a broken-symmetry
30+
determinant in the particle-number and projected spin symmetry mode, and then transformed
31+
to the particle-number U1 symmetry mode.
32+
33+
SP-DMRG is performed with particle-number U1 symmetry only, and the final MPS is transformed
34+
to the SU2 symmetry mode (``ket2, tag='KETX-0'``) which can be later loaded in the SU2
35+
symmetry mode to do spin-adapted DMRG with larger bond dimensions (not performed here).
36+
37+
.. code-block:: python3
38+
import numpy as np, sys
39+
import itertools
40+
from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes
41+
from pyblock2.algebra.io import MPSTools
42+
43+
istate = int(sys.argv[1])
44+
45+
driver = DMRGDriver(scratch="/tmp", symm_type=SymmetryTypes.SAnySZ, stack_mem=120 << 30, fp_codec_cutoff=0.0, n_threads=64)
46+
47+
bond_dims = [50] * 8 + [100] * 8
48+
noises = [1E-5] * (len(bond_dims) - 4) + [0] * 4
49+
thrds = [1E-7] * len(bond_dims)
50+
n_sweeps = len(bond_dims)
51+
52+
driver.read_fcidump(filename='Fe4S4.FCIDUMP', pg='d2h')
53+
driver.spin = 0
54+
twos = 0
55+
56+
npts = driver.get_spin_projection_npts(n_sites=driver.n_sites, n_elec=driver.n_elec, twos=twos)
57+
print("NPTS = %d" % npts)
58+
59+
driver.set_symmetry_groups("U1Fermi", "AbelianPG")
60+
driver.initialize_system(n_sites=driver.n_sites, n_elec=driver.n_elec, spin=driver.spin, orb_sym=driver.orb_sym)
61+
62+
mpo = driver.get_qc_mpo(h1e=driver.h1e, g2e=driver.g2e, ecore=driver.ecore, iprint=2, simple_const=True, add_ident=False)
63+
print("MPO = ", mpo.get_bond_dims())
64+
pmpo = driver.get_spin_projection_mpo(twos=twos, twosz=driver.spin, npts=npts, use_sz_symm=False, cutoff=1E-12, add_ident=True, iprint=1)
65+
66+
target = driver.target
67+
68+
driver.set_symmetry_groups("U1Fermi", "U1", "AbelianPG")
69+
driver.initialize_system(n_sites=driver.n_sites, n_elec=driver.n_elec, spin=driver.spin, orb_sym=driver.orb_sym)
70+
71+
n_sites = driver.n_sites
72+
73+
xdstr = [
74+
'22aaaaa2aaaa222222222222b2bbbbbbbb22',
75+
'22aaaaabbb2b222222222222a2aaabbbbb22',
76+
'222aaaab2bbb222222222222bbbbbaaaaa22',
77+
][istate]
78+
79+
print(istate, xdstr)
80+
81+
ket = driver.get_mps_from_csf_coefficients([xdstr], dvals=[1.0], tag='KET', dot=1)
82+
driver.align_mps_center(ket, ref=0)
83+
ket = driver.adjust_mps(ket, dot=2)[0]
84+
pket = driver.mps_change_symm(ket, 'PKET-0', target)
85+
86+
energy = driver.dmrg(mpo, pket, stacked_mpo=pmpo, metric_mpo=pmpo, context_ket=ket, n_sweeps=n_sweeps, bond_dims=bond_dims,
87+
noises=noises, thrds=thrds, lowmem_noise=True, twosite_to_onesite=None, tol=1E-12, cutoff=1E-24, iprint=2,
88+
dav_max_iter=400, dav_def_max_size=20)
89+
print('DMRG energy = %20.15f' % energy)
90+
91+
pmpo, mpo = None, None
92+
93+
ket = driver.adjust_mps(ket, dot=1)[0]
94+
driver.align_mps_center(ket, ref=0)
95+
96+
pyket = MPSTools.from_block2(ket)
97+
pyuket = MPSTools.trans_sz_to_su2(pyket, driver.basis, ket.info.target, target_twos=0)
98+
99+
driver.symm_type = SymmetryTypes.SU2
100+
driver.initialize_system(n_sites=driver.n_sites, n_elec=driver.n_elec, spin=driver.spin, orb_sym=driver.orb_sym)
101+
102+
impo = driver.get_identity_mpo()
103+
hmpo = driver.get_qc_mpo(h1e=driver.h1e, g2e=driver.g2e, ecore=driver.ecore, iprint=2, simple_const=True, add_ident=True, fast_no_orb_dep_op=True)
104+
105+
ket2 = MPSTools.to_block2(pyuket, driver.basis, tag='KETX-0')
106+
ket2.info.save_data(driver.scratch + "/%s-mps_info.bin" % ket2.info.tag)
107+
ket2.load_tensor(ket2.center)
108+
ket2.tensors[ket2.center].normalize()
109+
ket2.save_tensor(ket2.center)
110+
ket2.unload_tensor(ket2.center)
111+
norm = driver.expectation(ket2, impo, ket2)
112+
print('Norm = ', norm)
113+
114+
ket2.info.load_mutable()
115+
print('UMPS MAX BOND = ', ket2.info.get_max_bond_dimension())
116+
117+
energy = driver.expectation(ket2, hmpo, ket2, iprint=2)
118+
print('STATE %d Expt energy = %20.15f' % (istate, energy))
119+
120+
fe_idxs = [[2, 3, 4, 5, 6], [7, 8, 9, 10, 11], [24, 25, 26, 27, 28], [29, 30, 31, 32, 33]]
121+
122+
dm = driver.get_npdm(ket2, pdm_type=2, npdm_expr='((C+D)2+(C+D)2)0', mask=(0, 0, 1, 1), iprint=2, max_bond_dim=3000)
123+
dm = dm * (0.5 * -np.sqrt(3) / 2)
124+
fe_idxs = np.array([x for xx in fe_idxs for x in xx], dtype=int)
125+
dm = np.einsum('ijkl->ik', dm[fe_idxs, :][:, fe_idxs].reshape((4, 5, 4, 5)))
126+
127+
import matplotlib.pyplot as plt
128+
plt.matshow(dm, cmap='ocean_r')
129+
plt.gcf().set_dpi(300)
130+
plt.savefig("%02d-bip-spin-corr.png" % istate, dpi=300)
131+
132+
133+
Compressed MPO
134+
--------------
135+
136+
In the second example, we use a compressed Hamiltonian MPO, which can potentially save
137+
some computational cost. Note that to ensure that the Hamiltonian exactly preserves the
138+
total spin symmetry, the SVD compression needs to be done in the SU2 symmetry mode.
139+
After compression, the Hamiltonian MPO is transformed to lower symmetries.
140+
141+
.. code-block:: python3
142+
import numpy as np, sys
143+
import itertools
144+
from pyblock2.driver.core import DMRGDriver, SymmetryTypes, MPOAlgorithmTypes
145+
from pyblock2.algebra.io import MPSTools
146+
147+
istate = int(sys.argv[1])
148+
149+
driver = DMRGDriver(scratch="/tmp", symm_type=SymmetryTypes.SAnySU2, stack_mem=120 << 30, fp_codec_cutoff=0.0, n_threads=64)
150+
151+
bond_dims = [50] * 8 + [100] * 8
152+
noises = [1E-5] * (len(bond_dims) - 4) + [0] * 4
153+
thrds = [1E-7] * len(bond_dims)
154+
n_sweeps = len(bond_dims)
155+
156+
driver.read_fcidump(filename='Fe4S4.FCIDUMP', pg='d2h')
157+
driver.spin = 0
158+
twos = 0
159+
160+
npts = driver.get_spin_projection_npts(n_sites=driver.n_sites, n_elec=driver.n_elec, twos=twos)
161+
print("NPTS = %d" % npts)
162+
163+
driver.set_symmetry_groups("U1Fermi", "SU2", "SU2", "AbelianPG")
164+
driver.initialize_system(n_sites=driver.n_sites, n_elec=driver.n_elec, spin=driver.spin, orb_sym=driver.orb_sym)
165+
166+
umpo = driver.get_qc_mpo(h1e=driver.h1e, g2e=driver.g2e, ecore=driver.ecore, iprint=2, simple_const=True, add_ident=False,
167+
algo_type=MPOAlgorithmTypes.FastBlockedSVD, cutoff=1E-7, integral_cutoff=1E-12, fast_no_orb_dep_op=True)
168+
print("UMPO = ", umpo.get_bond_dims())
169+
170+
driver.set_symmetry_groups("U1Fermi", "U1", "AbelianPG")
171+
driver.initialize_system(n_sites=driver.n_sites, n_elec=driver.n_elec, spin=driver.spin, orb_sym=driver.orb_sym)
172+
173+
zmpo = driver.mpo_change_symm(umpo, add_ident=False)
174+
umpo = None
175+
print("ZMPO = ", zmpo.get_bond_dims())
176+
177+
driver.set_symmetry_groups("U1Fermi", "AbelianPG")
178+
driver.initialize_system(n_sites=driver.n_sites, n_elec=driver.n_elec, spin=driver.spin, orb_sym=driver.orb_sym)
179+
180+
mpo = driver.mpo_change_symm(zmpo, add_ident=True)
181+
zmpo = None
182+
print("MPO = ", mpo.get_bond_dims())
183+
pmpo = driver.get_spin_projection_mpo(twos=twos, twosz=driver.spin, npts=npts, use_sz_symm=False, cutoff=1E-12, add_ident=True, iprint=1)
184+
185+
driver.set_symmetry_groups("U1Fermi", "AbelianPG")
186+
driver.initialize_system(n_sites=driver.n_sites, n_elec=driver.n_elec, spin=driver.spin, orb_sym=driver.orb_sym)
187+
188+
target = driver.target
189+
190+
driver.symm_type = SymmetryTypes.SAnySZ
191+
driver.set_symmetry_groups("U1Fermi", "U1", "AbelianPG")
192+
driver.initialize_system(n_sites=driver.n_sites, n_elec=driver.n_elec, spin=driver.spin, orb_sym=driver.orb_sym)
193+
194+
n_sites = driver.n_sites
195+
196+
xdstr = [
197+
'22aaaaa2aaaa222222222222b2bbbbbbbb22',
198+
'22aaaaabbb2b222222222222a2aaabbbbb22',
199+
'222aaaab2bbb222222222222bbbbbaaaaa22',
200+
][istate]
201+
202+
print(istate, xdstr)
203+
204+
ket = driver.get_mps_from_csf_coefficients([xdstr], dvals=[1.0], tag='KET', dot=1)
205+
driver.align_mps_center(ket, ref=0)
206+
ket = driver.adjust_mps(ket, dot=2)[0]
207+
pket = driver.mps_change_symm(ket, 'PKET-0', target)
208+
209+
energy = driver.dmrg(mpo, pket, stacked_mpo=pmpo, metric_mpo=pmpo, context_ket=ket, n_sweeps=n_sweeps, bond_dims=bond_dims,
210+
noises=noises, thrds=thrds, lowmem_noise=True, twosite_to_onesite=None, tol=1E-12, cutoff=1E-24, iprint=2,
211+
dav_max_iter=400, dav_def_max_size=20)
212+
print('DMRG energy = %20.15f' % energy)
213+
214+
pmpo, mpo = None, None
215+
216+
ket = driver.adjust_mps(ket, dot=1)[0]
217+
driver.align_mps_center(ket, ref=0)
218+
219+
pyket = MPSTools.from_block2(ket)
220+
pyuket = MPSTools.trans_sz_to_su2(pyket, driver.basis, ket.info.target, target_twos=0)
221+
222+
driver.symm_type = SymmetryTypes.SU2
223+
driver.initialize_system(n_sites=driver.n_sites, n_elec=driver.n_elec, spin=driver.spin, orb_sym=driver.orb_sym)
224+
225+
impo = driver.get_identity_mpo()
226+
hmpo = driver.get_qc_mpo(h1e=driver.h1e, g2e=driver.g2e, ecore=driver.ecore, iprint=2, simple_const=True, add_ident=True, fast_no_orb_dep_op=True)
227+
228+
ket2 = MPSTools.to_block2(pyuket, driver.basis, tag='KETX-0')
229+
ket2.info.save_data(driver.scratch + "/%s-mps_info.bin" % ket2.info.tag)
230+
ket2.load_tensor(ket2.center)
231+
ket2.tensors[ket2.center].normalize()
232+
ket2.save_tensor(ket2.center)
233+
ket2.unload_tensor(ket2.center)
234+
norm = driver.expectation(ket2, impo, ket2)
235+
print('Norm = ', norm)
236+
237+
ket2.info.load_mutable()
238+
print('UMPS MAX BOND = ', ket2.info.get_max_bond_dimension())
239+
240+
energy = driver.expectation(ket2, hmpo, ket2, iprint=2)
241+
print('STATE %d Expt energy = %20.15f' % (istate, energy))
242+
243+
fe_idxs = [[2, 3, 4, 5, 6], [7, 8, 9, 10, 11], [24, 25, 26, 27, 28], [29, 30, 31, 32, 33]]
244+
245+
dm = driver.get_npdm(ket2, pdm_type=2, npdm_expr='((C+D)2+(C+D)2)0', mask=(0, 0, 1, 1), iprint=2, max_bond_dim=3000)
246+
dm = dm * (0.5 * -np.sqrt(3) / 2)
247+
fe_idxs = np.array([x for xx in fe_idxs for x in xx], dtype=int)
248+
dm = np.einsum('ijkl->ik', dm[fe_idxs, :][:, fe_idxs].reshape((4, 5, 4, 5)))
249+
250+
import matplotlib.pyplot as plt
251+
plt.matshow(dm, cmap='ocean_r')
252+
plt.gcf().set_dpi(300)
253+
plt.savefig("%02d-svd-spin-corr.png" % istate, dpi=300)
254+

‎docs/source/user/references.rst‎

Lines changed: 9 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -112,18 +112,23 @@ Multi-Reference Correlation Theories
112112
* Sharma, S., Alavi, A. Multireference Linearized Coupled Cluster Theory for Strongly Correlated Systems Using Matrix Product States. *The Journal of Chemical Physics* 2015, **143**, 102815. doi: `10.1063/1.4928643 <https://doi.org/10.1063/1.4928643>`_
113113
* Sharma, S., Jeanmairet, G., Alavi, A. Quasi-Degenerate Perturbation Theory Using Matrix Product States. *The Journal of Chemical Physics* 2016, **144**, 034103. doi: `10.1063/1.4939752 <https://doi.org/10.1063/1.4939752>`_
114114

115-
* Larsson, H. R., Zhai, H., Gunst, K., Chan, G. K. L. Matrix product states with large sites. *Journal of Chemical Theory and Computation* 2022, **18**, 749-762. doi: `10.1021/acs.jctc.1c00957 <https://doi.org/10.1021/acs.jctc.1c00957>`_
115+
* Larsson, H. R., Zhai, H., Gunst, K., Chan, G. K.-L. Matrix product states with large sites. *Journal of Chemical Theory and Computation* 2022, **18**, 749-762. doi: `10.1021/acs.jctc.1c00957 <https://doi.org/10.1021/acs.jctc.1c00957>`_
116116

117117
Determinant Coefficients
118118
------------------------
119119

120-
* Lee, S., Zhai, H., Sharma, S., Umrigar, C. J., Chan, G. K. Externally Corrected CCSD with Renormalized Perturbative Triples (R-ecCCSD (T)) and the Density Matrix Renormalization Group and Selected Configuration Interaction External Sources. *Journal of Chemical Theory and Computation* 2021, **17**, 3414-3425. doi: `10.1021/acs.jctc.1c00205 <https://doi.org/10.1021/acs.jctc.1c00205>`_
120+
* Lee, S., Zhai, H., Sharma, S., Umrigar, C. J., Chan, G. K.-L. Externally Corrected CCSD with Renormalized Perturbative Triples (R-ecCCSD (T)) and the Density Matrix Renormalization Group and Selected Configuration Interaction External Sources. *Journal of Chemical Theory and Computation* 2021, **17**, 3414-3425. doi: `10.1021/acs.jctc.1c00205 <https://doi.org/10.1021/acs.jctc.1c00205>`_
121121

122122
Perturbative DMRG
123123
-----------------
124124

125-
* Guo, S., Li, Z., Chan, G. K. L. Communication: An efficient stochastic algorithm for the perturbative density matrix renormalization group in large active spaces. *The Journal of chemical physics* 2018, **148**, 221104. doi: `10.1063/1.5031140 <https://doi.org/10.1063/1.5031140>`_
126-
* Guo, S., Li, Z., Chan, G. K. L. A perturbative density matrix renormalization group algorithm for large active spaces. *Journal of chemical theory and computation* 2018, **14**, 4063-4071. doi: `10.1021/acs.jctc.8b00273 <https://doi.org/10.1021/acs.jctc.8b00273>`_
125+
* Guo, S., Li, Z., Chan, G. K.-L. Communication: An efficient stochastic algorithm for the perturbative density matrix renormalization group in large active spaces. *The Journal of chemical physics* 2018, **148**, 221104. doi: `10.1063/1.5031140 <https://doi.org/10.1063/1.5031140>`_
126+
* Guo, S., Li, Z., Chan, G. K.-L. A perturbative density matrix renormalization group algorithm for large active spaces. *Journal of chemical theory and computation* 2018, **14**, 4063-4071. doi: `10.1021/acs.jctc.8b00273 <https://doi.org/10.1021/acs.jctc.8b00273>`_
127+
128+
Spin-Projected DMRG
129+
------------------
130+
131+
* Li, Z., Chan, G. K.-L. Spin-Projected Matrix Product States: Versatile Tool for Strongly Correlated Systems. *Journal of Chemical Theory and Computation* 2017, **13**, 2681-2695. doi: `10.1021/acs.jctc.7b00270 <https://doi.org/10.1021/acs.jctc.7b00270>`_
127132

128133
Orbital Reordering
129134
------------------

0 commit comments

Comments
 (0)