[Not-So-Simple-DMRG]: A High-Performance Massively Parallel DMRG library to solve interacting spin-S (Heisenberg) and fermions (Hubbard) models in 1D. The software exploits consevation of Abelian Quantum Numbers (e.g.
The structure of this code is largely inspired by the excellent simple-DMRG project: GitHub and Zenodo, but somehow diverged in its very nature, thus the name.
-
- Milestone 1 Basic development and testing.
- Milestone 2 Quantum Number conservation.
- Milestone 3 Code reshuffle, SuperBlock modules.
- Milestone 4 Finite DMRG and Measure.
-
Milestone 5 Direct
$H_{SB}\vec{v}$ product. Serial. - Milestone 6 Further development
- Milestone 7 Distributed parallel setup. MPI-DMRG.
-
Milestone 8 Profiling and
BlockI/O restart. - Milestone 9 Implement Periodic Boundary Conditions.
- Milestone 10 Fix Symmetry Fragmentation using MPI sub-communicators.
- Milestone 11 Lazy direct MVP operator filtering.
- Milestone 12 Memory management and further optimizations. Fix restart and checkpointing.
- Milestone 13 Observables and multi-point static correlation functions. Structure factors and topological marker.
- Milestone 14 Change in the QN target definition.
- Milestone 15 Wave-function prediction and PRIMME eigensolver interface.
The code is based on:
-
MPI
This DMRG library is available in the form of a static Fortran library (libdmrg.a) and the related Fortran module DMRG.
A standard installation from source is available through CMake, via the conventional out-of-source method.
Clone the repo
git clone https://github.com/aamaricci/Lattice_DMRG
Optionally export the correct Fortran compiler:
export FC=mpif90/gfortran/ifort
From the repository directory (cd Lattice_DMRG) configure and compile the library:
mkdir build
cd build
cmake ..
make -j
The building step can be further configured passing suitable variables to CMake:
-DBUILD_TYPE=RELEASE [/DEBUG/TESTING/AGGRESSIVE]. Set build mode flags.-DUSE_MPI= yes(True)/no(False). Enable MPI support at compile time.-DVERBOSE=yes(True)/no(False). Verbose CMake output. Superseded bymake VERBOSE=yes/no.
WARNING BUILD_TYPE=AGGRESSIVE includes many deep level debug options which might not compile on some systems or breakdown compilation at linking step.
System-wide installation is completed after the build step using
make install
Please follow the instructions on the screen to complete installation on your environment. The library can be loaded using one of the following, automatically generated, files :
- A generated environment module , installed to
$HOME/.modules.d/dmrg/${PREFIX}and available after module environment setup (instruction on the screen) - A generated bash script installed to
${PREFIX}/bin/dmrg_ 8000 config_user.sh, to be sourced for static loading. - A generated pkg-config fileinstalled to
~/.pkg-config.d/dmrg.pc
For ease of use a specific and automatically generated recap message is printed after installation.
- Extend Spin iDMRG code to Electronic case from the Spin problem.
- Solve the fermion sign problem.
- Develop test code with recursive block construction and digonalization.
- Develop code without Quantum Number (QN) conservation.
- Test code against ED solution and exact behavior for the non-interacting regime.
- Develop the infrastructure to deal with QN conservation.
- Develop the construction of the states of a given QN sector.
- Code recursive blocking withing a specific QN sector. Compare with ED.
- Develop full iDMRG algorithm.
- Test code against iDMRG without QN and ED.
- Wrap the SuperBlock construction into a dedicated module.
- Encapsulate and generalize the code structure.
- Merge with Spin DMRG code.
- Add layer to save rotation/truncation matrices
- Develop finite DMRG algorithm
- Measure local observables
- Measure entanglement
- Measure nearest neighbor correlations
- Implement a better strategy for the SuperBlock matrix-vector product
$H_{sb}|\psi\rangle$ , exploiting the tensor product structure of$H_{sb}= H_L\otimes 1_R + 1_L\otimes H_R + H_{LR}$ .
With this update we introduce two distinct method to evaluate and diagonalize the SuperBlock (SB) Hamiltonian$H_{sb}$ , according to the value of the parametersparse_H=T,F. IfTthe matrix$H_{sb}$ is evaluated as a sparse matrix performing the kroenecker products above. IfFeach non-trivial part of the the matrix$H_{sb}$ , i.e.$H_L$ ,$H_R$ and$H_{LR}$ are directly applied to the SB vector$|\psi\rangle$ .To improve execution each term is decomposed into its blocks structure corresponding to conserved quantum numbers$\vec{Q}=[Q_1,\dots,Q_N]$ . In addition, exploiting the fact that each SB vector can be decomposed into a matrix form of Left and Right states,$|\psi\rangle = \sum_{lr} c_{lr} |l\rangle|r\rangle$ , we can operate directly as follows:
$$H_L\otimes 1_R \circ |\psi\rangle \rightarrow H_L \circ \sum_l c_{lr}|l\rangle$$
$$1_L\otimes H_R \circ |\psi\rangle \rightarrow H_R \circ \sum_r c_{lr}|r\rangle$$
$$H_{LR} \circ |\psi\rangle \rightarrow \sum_p \sum_{Q_i} A_p(Q_i)\otimes B_p(Q_i) |\psi\rangle$$
where the matrices$A_p$ ,$B_p$ are the connecting operators (i.e.$c_{i\sigma}$ or$S_a(i)$ ) restricted between two blocks with quantum numbers$Q_i$ and$Q_j=Q_i-1_{\sigma/a}$ .
- Further development of finite DMRG algorithm.
- Massive distributed MPI parallelization of the DMRG algorithm.
So far different ideas have been be explored in the literature. Some focused on the parallelization of the DMRG algorithm itself, by distributing workload from different sites/MPS across nodes. Other have focused on accelarating the parallel execution of the tensor product for$H_{sb}$ forsparse_H=T, see for instance Ref.[1]. Here we take a slightly different approach. A simple profiling shows that more than 80% of the code executin time is spent in the solution of the eigenvalue problem for the superblock$H_{sb}|v\rangle = E|v\rangle$ . As we have shown in milestone 5, the super-block Hamiltonian$H_{sb}$ can be decomposed into different parts, each having the form of a suitable tensor product. This form is analog to the one assumed by Hubbard-like models with spin-conservation for which a massively parallel Exact Diagonalization algorithm is available, see Ref.[2]-[3]. Using this approach we constructed a simple yet efficient distributed MPI DMRG scheme as follows. We represent the generic superblock vector$|v\rangle \in {\cal H}_{sb}$ as a block matrix$V=diag[V_q]$ with$q=1,\dots,N_q$ the quantum number index and$N_q$ the total number of quantum number sectors contributing to the current superblock. Each block$V_q$ corresponds to a sub-matrix with$D_R(q)$ rows and$D_L(q)$ columns.
In the current MPI impementation of the DMRG we distribute across the$M$ nodes a share of$Q=D_L(q)/M$ columns. As we will explain in the following, this enables to perform application of each term in$H_{sb}$ locally in the nodes$p$ memory with only minimal MPI communication overhead. For, we consider the functionvector_transpose_MPIinDMRG_GLOBAL. This function exploits the MPI functionMPI_All2AllVto perform the parallel transposition of a given matrix$A$ which is distributed column-wise. In addition we note that each term appearing in the expression for$H_{sb}$ corresponds to a given block with quantum number$q$ , as by dictated by the symmetries of the problem:$H_{sb}=diag[H_{sb}(q)]$ . Thus we can analyze the matrix-vector product (MVP) of each term with a given quantum number one-by-one:
-
$1_L\otimes H_R$ : when applied to the block$V_q$ this term involves multiplication of the (small) matrix$H_R$ rows with columns of$V_q$ with indices$j_q(p)$ , with$p$ the node index. This multiplication is local in the memory of the node and does not required communication. -
$H_L\otimes 1_R$ : when applied to the block$V_q$ this term involves multiplication of the (small) matrix$H_L$ rows with the rows of$V_q$ with indices$i_q(p)$ . As the rows are distributed among the nodes a global communication is required. To optimize this step we operate as follows:- MPI-transpose the block
$V_q\rightarrow V^T_q$ , - Apply the matrix
$H_L$ which now operates locally in the nodes memory. - MPI-transpose back the result
$(H_LV^T_q)^T\rightarrow H_LV$ .
- MPI-transpose the block
-
$H_{LR} = \sum_a \sum_{k} A_a(k)\otimes B_a(k)$ : This term requires a more involved treatment. First we observe that each term in the double sum applies as:$$A_a(k)\otimes B_a(k) |v_q\rangle = A_a(k)\otimes B_a(k) \cdot vec(V_q)$$ where the operator$vec({\cdot})$ takes the vector on input and transforms it into a matrix columns-wise. By a well-known properties of the tensor products the last term is equal to:$$vec{(B_a(q,k)\cdot V_k \cdot A_a(k,q)^T)}$$ which is an apparently involved double matrix-matrix product (MMP). However, as indicated in the last expression, because of the block structure of the operators$A_a$ and$B_a$ (indicated by the index$k$ ) some restrictions applies to this product: only the off-diagonal block components of$A_a$ and$B_a$ which ensures the final result of the MMP contributes to the specific$q$ quantum number are possible. Recalling that each$V_k$ is distributed column-wise we perform the first product in parallel as for$1_L\otimes H_R$ :$C_a(q,k)= B_a(q,k)\cdot V_k$ . The resulting dense matrix$C_a(q,k)$ is distributed column-wise by construction. We are then left with a final MMP:
$$vec{(C_a(q,k)\cdot A_a(k,q)^T )} = vec{( [A_a(q,k)\cdot C^T_a(k,q)]^T )}$$
which, as for the term$H_L\otimes 1_R$ requires using MPI-transpose two times: i) to get$C^T_a(k,q)$ and ii) to transpose the final result to be accumulated in the outgoing vector. The overall cost of communication is minimized relying exclusively on theMPI_All2All-type of communication, thus unlocking massively parallel scaling of the MVP at the heart of the super-block diagonalization.
- Write/Read
Blockobjects to restart DMRG. - Implement accurate timing of the different DMRG steps beyond the actual form.
- Development of PBC strategies using double step 2-site increment.
- Adapt the measurement strategy to the PBC
See graphics upon running to get a glimpse of how PBC are implemented. Schematically the difference is as follows:
With OBC we grow L.block to the right L=[o-o...o-o]-o and R.block to the left R=o-[o-o...o-o], then we form the superblock by connecting L.block with R.block and solve:
OBC: [o-o..L..o-o]-o-o-[o-o..R..o-o]
In PBC we need to close the chain links. To avoid having links between largely renormalized operators sitting at one end of the chain, we now grow each block in two directions (n, p).
So L.block grows at first step as [o-o...o-o]-o and the next step as o-[o-o...o-o-o], the R.block does the same in inverted order. The SuperBlock is formed by contacting the two blocks at each end at every step:
PBC 1:
.[o-o..L..o-o]-o-o-[o-o..R..o-o].
^-------------------------------^
PBC 2:
.o-[o-o..L..o-o-o]-[o-o-o..R..o-o]-o.
^-----------------------------------^
In this milestone we fixed the issue of symmetry fragmentation in the parallel implementation of the DMRG. The problem is that, as we distribute the superblock vector blocks
To solve this issue we implemented a strategy based on MPI sub-communicators, which allows to split the global communicator into smaller communicators, one per symmetry sectors, restricted to the maximum required number of cores fixed by the conditions:
$$
N_{active}(q) = \max( \min(D_L(q),N_{cores}), \min(D_R(q),N_{cores}) )
$$
This way we can ensure that each node has a more balanced workload and we can optimize the communication patterns accordingly. The implementation of this strategy required just few modifications to the existing code, essentially concentrated in the MPI matrix-vector product and vector_transpose_MPI function, where the communication is performed.
- Implement lazy direct MVP operator filtering to reduce the memory footprint of large SuperBlock calculations.
The direct matrix-vector product introduced in milestone 5 and parallelized in milestone 7 avoids constructing the full sparse SuperBlock Hamiltonian. However, in its cached form it still precomputes and stores, for every SuperBlock quantum-number sector
To address this regime we introduced an optional lazy direct MVP mode, controlled by the input variable DIRECT_H_LAZY. The relevant choices are now:
-
SPARSE_H = T: build and apply the sparse SuperBlock Hamiltonian$H_{sb}$ . -
SPARSE_H = F,DIRECT_H_LAZY = F: use the direct MVP with cached filtered operator blocks. This is the fastest direct mode, but stores all$H_L(q)$ ,$H_R(q)$ ,$A_a(q,k)$ and$B_a(q,k)$ blocks. -
SPARSE_H = F,DIRECT_H_LAZY = T: use the direct MVP with lazy sector filtering. Only the original enlarged block operators and light sector maps are kept persistently; the filtered operator blocks required by a given MVP term are constructed on the fly, applied, and immediately released.
In the lazy workflow the SuperBlock setup stores the sector state lists and inverse maps associated with the left and right block bases. The MVP then obtains the required restricted operators through local filtering operations such as: $$ O(q,k) = O_{\mathrm{full}}\big|_{{\cal H}(q)\rightarrow{\cal H}(k)} $$ only at the moment when that contribution is applied to the vector. Thus the memory scaling is governed by the original block operators plus the largest temporary filtered blocks needed for a single term, rather than by the sum over all terms and all SuperBlock sectors. The price is additional CPU time due to repeated filtering during the iterative eigensolver.
For fermionic models the left boundary hopping operators are prepared directly in the contracted form
The practical goal of DIRECT_H_LAZY = T is therefore to trade repeated sparse filtering work for a substantially smaller persistent memory footprint. This is useful when the calculation is memory-bound rather than CPU-bound, and especially when additional MPI ranks or threads can compensate for the slower MVP while enabling larger retained-state spaces.
For MPI builds compiled with profiling, HXV_DETAIL_PROFILE = T adds time_hxv_filter, time_hxv_local, and time_hxv_mpi to full_profile_iDMRG.dmrg when the lazy direct MVP is used. The last timer covers the entire vector_transpose_MPI call, including its setup, MPI collectives, local transpose, and synchronization. Each reported timer is the maximum over ranks; their values need not add up exactly to time_hxv_direct. The option defaults to F and does not change the Hamiltonian application.
LANC_TOLERANCE defaults to 1d-10. It is an ARPACK relative Ritz
convergence tolerance; GS_DEG_THRESHOLD remains a separate absolute energy
cutoff for selecting the states included in the ground-state density matrix.
For a nearly degenerate manifold, compare the low-energy subspace and the
subsequent DMRG trajectory before relaxing the solver tolerance further.
LANC_CHECK_RESIDUAL=T optionally measures the absolute norm
||H*psi_i-E_i*psi_i|| after each solve, at the cost of one extra H*v per
returned state.
Wave-function prediction and the PRIMME interface are summarized in Milestone 15.
- Memory management and further optimizations.
- Fix restart and checkpointing.
This version of the code (5.0.0) includes a complete revision of the memory footprint. Using several input variables we can now control how the results of the calculations are stored in the memory and/or the disk. The relevant variables are:
-
SAVE_BLOCK=T/F: Logical flag to trigger block storage. DEBUG enforces T. Default: T. This is used to save the block objects to file at each step of the DMRG, which can be used to restart the calculation from that given step. The block files are stored with the name prefixBLOCK_FILEplus a suffix. -
SAVE_ALL_BLOCKS=T/F: Logical flag to save all blocks (T) or just the last (F). DEBUG enforces T. Default: F
This is activated to save the blocks at any step to a different file. Mostly for debugging purposes, as it can lead to a large number of files and a large memory footprint. -
BLOCK_UMAT_CACHE=T/F: Logical flag to keep all rotation matrices in the block memory. Default=F. This is used to keep (T) or not (F) the rotation matrices in the block memory. While keeping the rotation matrices in memory can speed up the finite DMRG sweeps, it can lead to memory overflow in large calculations. -
SAVE_UMAT=T/F:Logical flag to save rotation matrices to file. Default=T. Rotation matrices are used to measure in DMRG, thus they are stored in the memory or to file according to the value ofBLOCK_UMAT_CACHE. IfTthe rotation matrices are stored in file with the name prefixUMAT_FILEplus a suffix. -
SAVE_MEASURE_STATE=T/F: Logical flag to save SuperBlock state for post-processing measures. Default=T The superblock calculations results are stored to file with the name prefixMEASURE_FILEplus a suffix. This is used to save the SuperBlock state at each step, which can be used to perform measurements in post-processing without the need to run the DMRG again. This is particularly useful for large calculations, where the DMRG can be time consuming, while the measurements can be performed in a much faster post-processing step. -
BLOCK_FILE=block: Name prefix of the stored block file, used to restart DMRG. -
UMAT_FILE=umat:Name suffix of the stored rotation matrices, used to measure in DMRG. -
MEASURE_FILE=measure:Name prefix of the stored SuperBlock measurement state.
On top of that this version includes a number of optimizations and bug fixes which removed linear growing of memory footprint with the number of DMRG steps. Larger calculations can now be performed with a reasonable memory footprint, while the restart and checkpointing functionality is now fully working.
- Multi-point static correlation functions, Structure factors and string observables.
- Restartable post-processing with OBC/PBC-specific structure-factor paths.
The following table lists every public interface intended for external
drivers. In the signatures below, brackets denote optional arguments and
-> denotes a returned value or an intent(out) argument.
Notation:
K = operator key,
O = sparse_matrix,
R = real(8),
Z = complex(8),
I = integer,
L = logical,
V = real/complex according to the build.
The measurement session is initialized with Init_Measure_DMRG and released with End_Measure_DMRG. If a restart structure is present the session uses the stored rotations, sector maps and ground-state data.
Main routines:
| Interface | Inputs | Output |
|---|---|---|
Init_Measure_DMRG |
[msg] |
Opens/rebuilds the measurement state |
End_Measure_DMRG |
— | Releases the measurement state |
Measure_DMRG |
Op or Ops(:), [pos(:),file] |
[avOp(:)] or [avOp(:,:)] |
Measure_Op_DMRG |
Op,pos |
R = <Op(pos)> |
A sequence of lower-level calls can be enclosed in Init_Measure_DMRG / End_Measure_DMRG to reuse the loaded rotations, sector
maps and ground-state data. High-level interfaces initialize lazily when needed.
Measure_DMRG measures one or several local opeators (as sparse matrices) over pos. Without pos it uses the complete chain. The scalar-operator overload returns avOp(Npos), while the operator-array overload returns avOp(Noperators,Npos).
If file is present, the same values are appended to
the suffixed output file.
A number of generic correlation and ordered-product measurements are available. User can pass the key of the operator, or the operator itself as a sparse matrix.
The key overloads are preferred: operator matrix, QN shift and fermionic
grading are read consistently from LIST_OPERATORS. In the matrix overload,
dqA and dqB follow O |q> -> |q+dq(O)>.
type='fermionic' marks an odd operator and activates the required
Jordan--Wigner strings; parity-even operators may use bosonic, psign, or
the corresponding stored type. connected=.true. subtracts
<A(posA)><B(posB)>. Products preserve the supplied factor order, including
repeated positions.
Products preserve the supplied factor order, including repeated positions.
The following table lists the public interfaces and their overloads:
| Interface | Overload | Signature |
|---|---|---|
Measure_Corr_DMRG |
keys | (keyA,keyB,posA,posB,[connected]) -> V |
| matrices | (OpA,dqA,OpB,dqB,posA,posB,[typeA,typeB,connected]) -> V |
|
Measure_Product_DMRG |
keys | (keys(:),positions(:)) -> V |
| matrices | (Ops(:),dqs(Qdim,M),types(:),positions(:)) -> V |
|
Measure_SpinSpin_DMRG |
— | (posA,posB) -> <S(posA).S(posB)> |
Measure_DensityDensity_DMRG |
— | (posA,posB,[connected]) -> V(Nso,Nso) |
Structure factors and nonlocal observables are measured by specific routines, in real or reciprocal space. The following table lists the public interfaces and their overloads:
| Interface | Overload | Required inputs and output |
|---|---|---|
Measure_Structure_DMRG |
keys | (keyA,keyB,q(:),Fq(:)->, options) |
| matrices | (OpA,dqA,OpB,dqB,q(:),Fq(:)->, options) |
|
Measure_Structure_Aq_DMRG |
key | (key,q(:),Fq(:)->,[positions,connected]) |
| matrix | (Op,dq,q(:),Fq(:)->,[positions,connected,type]) |
|
Measure_Polarization_DMRG |
key/matrix | (key or Op,[positions,x,length]) -> Z |
Measure_String_DMRG |
— | (Op,i,j,[theta]) -> Z |
Measure_Parity_DMRG |
— | (Op,i,j,[reference,theta]) -> Z |
Common optional arguments of Measure_Structure_DMRG are:
| Option | Type | Meaning |
|---|---|---|
positions |
I(:) |
Selected sites; default: complete chain |
connected |
L |
Subtract one-point contribution |
Cij |
Z(:,:) output |
Complete matrix; exact mode only |
mode |
character | exact (default), auto, bulk, periodic/pbc |
same_operator |
L |
Caller assertion A=B |
hermitian |
L |
Caller assertion A=A^dagger |
commuting |
L |
Caller assertion [A_i,A_j]=0, i/=j |
rmax |
I |
Maximum bulk distance |
norigins |
I |
Origins sampled in bulk/periodic mode |
window_weights |
L |
Apply `1- |
Cr |
Z(:) output |
Distance-resolved correlations |
The matrix overload of Measure_Structure_DMRG also accepts
[typeA:character,typeB:character], with the same fermionic-grading convention
as Measure_Corr_DMRG.
mode='auto' selects periodic only when PBCdmrg is true, positions is
the complete ordered ring, every momentum is q=2*pi*n/M, and options that
require the exact pair matrix are absent. It otherwise falls back to exact
and writes the selected route to LOGfile. This permits one common OBC/PBC
post-processing call without silently applying wrapped-distance assumptions to
an OBC window or to an incomplete PBC selection.
The Measure_Structure_Aq_DMRG key overload additionally accepts
[positions,connected]; its matrix overload accepts
[positions,connected,type]. It requires a Hermitian, parity-even,
QN-conserving local operator and computes <Aq^dagger Aq>/M without assuming
translation invariance.
Measure_Polarization_DMRG accepts [positions:I(:),x:R(:),length:R].
positions selects the local factors, x supplies their coordinates, and
length is the period used in the phase. Measure_String_DMRG defaults to
theta=pi; Measure_Parity_DMRG defaults to reference=2 and theta=pi.
| Interface | Signature |
|---|---|
Measure_FermionBond_DMRG |
(posA,posB,Hij(:,:)) -> Eij |
Measure_SpinBond_DMRG |
(posA,posB,Hij(:,:)) -> Eij |
Measure_SpinExchangeEnergy_DMRG |
(Hij(:,:),[Jij->]) -> Espin |
Measure_LocalEnergy_DMRG |
([Hi->,H0loc,Hint,Hshift,E0loc->,Eint->,Eshift->]) -> Eloc |
Measure_Energy_DMRG |
(Hij,Ebond->,Eloc->,Etotal->,[E0loc->,Eint->,Eshift->,H0loc,Hint,Hshift]) |
For Measure_LocalEnergy_DMRG and Measure_Energy_DMRG, each optional scalar
component must be supplied together with its matching operator: (E0loc,H0loc),
(Eint,Hint), or (Eshift,Hshift). Hi, when requested, is a sparse diagonal
container of site-resolved local energies. Jij stores each measured spin-bond
energy once in its upper triangle.
| Interface | Signature |
|---|---|
Structure_Factor_DMRG |
(C(:,:),q(:),[x(:)]) -> Fq(:) |
Correlation_Distance_DMRG |
(C,positions,r->,Cr->,counts->) |
Build_Op_DMRG |
(Op,pos,[set_basis]) -> Oi |
Advance_Op_DMRG |
(Op,pos,[nstep]) -> Oi |
Advance_Corr_DMRG |
(Op,pos,[nstep]) -> Oi |
Average_Op_DMRG |
(Oi,pos) -> <Oi> |
Write_DMRG |
(file,val or vals(:)/vals(:,:),[x]) |
Structure_Factor_DMRG Fourier-transforms an already available correlation
matrix and performs no DMRG contractions. Correlation_Distance_DMRG groups
the entries of C by absolute physical separation and returns distances,
averages and pair counts. The Build/Advance/Average routines expose the
internal operator-propagation stages for advanced code; ordinary analysis
drivers should normally use Measure_Op_DMRG, Measure_Corr_DMRG, or the
structure-factor interfaces. Write_DMRG appends the current DMRG suffix to
file and writes scalar, vector or column-oriented matrix data.
Measure_Structure_DMRG evaluates
| Path | DMRG contractions | Assumption/use |
|---|---|---|
mode='exact' |
|
Exact pair sum; optionally returns Cij
|
mode='bulk' |
approximately norigins*rmax
|
OBC bulk sampling; tests distance homogeneity |
mode='periodic' |
norigins*M |
Full PBC ring, |
mode='auto' |
dispatches above | Uses periodic only when its PBC conditions are satisfied; otherwise exact
|
Measure_Structure_Aq_DMRG |
approximately Nq*M block propagations |
Exact collective-operator route, no translation assumption |
Existing calls to Measure_Structure_DMRG keep the default mode='exact' which analyzes all pairs i and j.
For identical Hermitian operators same_operator=.true., hermitian=.true., commuting=.true. to measure only
one triangle, including the diagonal. These flags are caller assertions.
The code checks operator identity but does not infer the remaining properties.
The exact result is no
F58B
rmalized by the number of selected sites. C_AB(i,j) is
allocated only if requested (optional input). Otherwise the Fourier sum is accumulated directly.
call Measure_Structure_DMRG(Sz, dq0, Sz, dq0, q, Sq, &
positions=bulk_sites, mode='bulk', rmax=30, norigins=2, &
same_operator=.true., hermitian=.true., commuting=.true., &
connected=.true., Cr=Cr)bulk_sites is the list of sites from which pairs i and j are obtained.
It must be consecutive and ascending. Bulk mode measures central
pairs for each distance, including zero; it assumes correlations in the
selection depend mainly on distance. This is an approximation in OBC and
must be checked against different origins, windows and distances. Merely
setting symmetry flags does not assert translation invariance.
rmax defaults to half the selected window (at most M-1); norigins defaults
to one. Origins are adjacent and centered for each distance; near the maximum
distance their number is capped at the available pairs. Both parameters refer
to the selected window. Without all three symmetry flags, positive and negative
distances are measured separately. Cr has bounds -rmax:rmax; bulk mode
rejects Cij, because no complete correlation matrix has been measured.
By default mode=bulk returns the unweighted truncated sum
F(q)=sum_r exp(i*q*r)*Cr(r), with no extra division by window size.
window_weights=.true. inserts 1-abs(r)/M to estimate the finite-window
pair sum with the same normalization as exact mode. Sampling every available
origin (norigins>=M) and setting rmax=M-1 with these weights reproduces
the exact sum, but restores quadratic measurement cost.
rmax, norigins, window_weights, and Cr are bulk-only options.
A truncated or sampled bulk transform need not be nonnegative, even when the
exact same-operator structure factor is. PBC Hamiltonians still use physical
non-wrapping separations here: no periodic distance reduction is inferred.
For a translation-invariant PBC correlator, mode='periodic' instead measures
the wrapped row of the correlation matrix:
call Measure_Structure_DMRG(keyA,keyB,q,Sq,mode='periodic',&
same_operator=.true.,hermitian=.true.,commuting=.true.)It evaluates F(q)=sum(r=0:M-1) exp(i*q*r)*C(r). The default
norigins=1 needs only M DMRG contractions, independently of the number of
requested momenta. Increasing norigins averages uniformly distributed rows
around the ring; norigins=M is exactly equivalent to the complete pairwise
sum divided by M and does not assume translation invariance. This periodic
Fourier convention requires the ring momenta q=2*pi*n/M; other momenta are
rejected. Periodic mode also requires PBCdmrg and the complete ordered ring
in positions (or no positions argument). It rejects Cij, rmax, and
window_weights; optional Cr is returned with bounds 0:M-1.
Measure_Structure_Aq_DMRG computes the same finite-selection structure factor
as Measure_Structure_DMRG(...,mode='exact'), without measuring individual
pairs or assuming translation symmetry:
call Measure_Structure_Aq_DMRG(Sz, dq0, q, Sq, connected=.true.)
call Measure_Structure_Aq_DMRG(density, dq0, q, Nq, positions=selected_sites)
! Or use a registered local operator key:
call Measure_Structure_Aq_DMRG('observable_key', q, Fq, connected=.true.)The direct route folds momenta related by q -> q + 2*pi*n or q -> -q
and performs the DMRG contractions only for the inequivalent values in
[0,pi]. The returned array still follows the order and length of the input
q. This is exact for the Hermitian, parity-even local operators accepted by
this interface and nearly halves the cost of a symmetric [-pi,pi] grid.
For M selected sites, connected=.true. subtracts
Coordinates are physical site numbers, even for irregular selections or PBC. There are no distance cutoffs or window weights: all selected pairs are included through the collective operator. The result is a complex array with zero imaginary part for the supported Hermitian observables.
This version supports parity-even Hermitian local operators that conserve all quantum numbers, including total local Sz and density. The matrix API checks Hermiticity and the supplied zero QN shift; the caller must give correct quantum-number/parity metadata. The key API checks stored metadata on every selected site and also permits site-dependent local matrices. QN-changing and parity-odd operators are not supported by this interface.
The implementation accumulates
- Change in the QN target definition.
- Input variables
QN_DENSITYandQN_OFFSETreplace the oldDMRG_QNandQNTYPErecords.
The new input variable QN_DENSITY sets the background charge per site and QN_OFFSET adds a fixed
total change:
Q = int(L * QN_DENSITY) + QN_OFFSET.
Both vectors have QNDIM entries and default to zero. For fermions density includes all orbitals, while offsets are not multiplied by NORB or by the local spin.
Example: for one-orbital fermions, half filling plus one up electron is:
QN_DENSITY = 0.5, 0.5
QN_OFFSET = 1.0, 0.0
Quarter filling uses density 0.25,0.25 and zero offsets. While for half-filled two-orbital Hubbard or BHZ models:
QN_DENSITY = 1.0, 1.0
QN_OFFSET = 0.0, 0.0
Non-commensurate lengths are accepted, e.g. rho=1/3 at L=8 gives two particles. Truncation is unchanged, including numerical effects: 0.33333333 * 6 gives one.
For both spin-1/2 and spin-1 chains, density zero and offset one select total Sz=1 as reversing a spin-1/2 changes Sz by one. An offset of 0.5 is impossible on four spin-1/2 sites, whose total Sz is integer. Changing the target selects another sector, not a specific excited wavefunction.
Fermion normal mode uses (N_up,N_down), nonsu2 their sum, and superc
their difference, equal to 2*Sz.
There a number of checks on the target. The requested total QN is compared to the local limits, which are set by the number of sites and orbitals. With equal local QN limits, the final target is checked before growth or restart loading. Five up electrons on four one-orbital sites, or Sz=0.5 on four spin-1/2 sites, fail immediately. During growth only, targets exceeding those limits are temporarily reduced and logged: density 0.5 plus offset 3 uses
four up electrons at L=4, then the requested six at L=6. The final target is never corrected. No fixed cutoff on the offset is imposed.
Same site family does not mean same QN limits. Different fields are fine; spin-1/2 and spin-1 both use type SPIN and physical Sz, but their limits differ. For mixed limits, early final checks and temporary corrections are disabled; superblock construction checks the actual states. Three spin-1/2 sites plus one spin-1 site allow Sz=0.5 but not Sz=0. Small mixed-spin examples agree with exact diagonalization; arbitrary ordering, measurements and mixed-chain restarts remain unvalidated. Spin/fermion mixtures are rejected.
Custom bases with missing QN combinations also require the superblock check. Truncation can cause a missing sector too: blocks retaining only Sz=0,1 cannot form total Sz=-1; the code reports this rather than changing the target.
The old QNTYPE and DMRG_QN records are removed. Migrate local inputs to density old_DMRG_QN * NORB with zero offset, and global inputs to zero density
with the previous effective target as offset. Restart formats are unchanged;
continuation requires matching the old targets. mainSep25 keeps the old code.
- Correlated wave-function prediction for OBC/PBC iDMRG.
- Serial and distributed PRIMME interfaces with block initial guesses.
- Optional Ritz-residual validation and eigensolver diagnostics.
LANC_METHOD selects the iterative solver used when the superblock is above
LANC_DIM_THRESHOLD. Small sectors still use dense LAPACK automatically.
| Capability | LANC_METHOD='arpack' |
LANC_METHOD='primme' |
|---|---|---|
| Default | Yes | No |
| Serial callback | scalar H*v |
block H*V adapter |
| MPI callback | distributed PARPACK-style vector | distributed PRIMME MatMat |
| Predicted initial state | one vector | orthonormal block |
Prediction with LANC_NEIGEN>1 |
disabled by input validation | supported |
| Returned diagnostics | NumOp |
NumOp, iterations, converged states, exit status |
| Failure handling | solver error | nonzero PRIMME status is reduced across ranks and is fatal |
The PRIMME path uses PRIMME_DEFAULT_MIN_MATVECS through SciFortran and
supports real/complex, serial/MPI builds. The existing superblock H*v
kernels remain the source of the matrix action; the PRIMME adapters only expose
the required block callback and do not assemble a new Hamiltonian. PRIMME
returns iteration/convergence/status counters to sb_diag; a nonzero exit
status is logged together with the number of converged states.
| Input | Default | Meaning |
|---|---|---|
LANC_METHOD |
arpack |
arpack or primme |
LANC_NEIGEN |
2 |
Requested eigenstates per superblock sector |
LANC_TOLERANCE |
1d-10 |
Relative iterative-solver tolerance |
LANC_NITER |
512 |
Maximum iterative-solver steps |
LANC_NCV_FACTOR |
10 |
ARPACK/solver search-space factor |
LANC_NCV_ADD |
0 |
Additive search-space size |
LANC_CHECK_RESIDUAL |
F |
Evaluate ` |
LANC_TARGET_PREDICTION |
F |
Enable correlated iDMRG initial guesses |
LANC_CHECK_RESIDUAL=T adds one extra Hamiltonian application per returned
state and reports the absolute Ritz residual. It is intended for validation;
GS_DEG_THRESHOLD remains an independent absolute energy threshold used to
select the ground-state manifold for the density matrix.
LANC_TARGET_PREDICTION=T is available for iDMRG only, with either OBC or
PBC. The first growth steps use the solver's ordinary initialization until the
previous center wavefunctions and left/right truncation isometries are
available. Later steps:
- project the previous target states through the new left/right isometries;
- join the two half states through an SVD pseudoinverse of the previous center;
- stream one local-site pair at a time into the current superblock ordering;
- orthonormalize the valid target columns before passing them to PRIMME.
PBC prediction explicitly follows the alternating n/p link ordering,
which inserts physical sites at opposite block edges. The streamed contraction
avoids a dense four-index wavefunction tensor. Rank zero constructs the trial
block; MPI runs scatter each valid column to the normal distributed
superblock-vector layout. If an MPI sector requires the serial empty-row
fallback, PRIMME can still receive the complete root-owned prediction block;
ARPACK retains its unseeded fallback there.
ARPACK accepts a predicted state only for LANC_NEIGEN=1; input parsing logs
the incompatibility and disables prediction otherwise. PRIMME accepts up to
LANC_NEIGEN independent predicted columns. Prediction is not currently
enabled for finite-DMRG sweeps.
LANC_METHOD=primme
LANC_NEIGEN=2
LANC_TARGET_PREDICTION=T
LANC_TOLERANCE=1d-10
LANC_CHECK_RESIDUAL=T
Reference BHZ tests reported the following solver-level reductions while preserving the converged energy within the quoted numerical precision:
| Geometry | Setup | H*v without/with prediction |
Energy agreement |
|---|---|---|---|
| OBC | one target, LANC_TOLERANCE=1d-10 |
150 -> 60 on eligible steps |
2e-14 per site |
| PBC | LDMRG=10, MDMRG=20, 4 MPI ranks |
50-55 -> 30-35 |
3.2e-14 per site |
These are small reference cases, not general performance guarantees. Peak
memory and wall-time scaling for large retained-state spaces, especially
around MDMRG=500, still require production benchmarking.
Here are some results for the Heisenberg model:
In the top panels we show the groundstate energy
Finally, in the bottom-right panel we report the spatial distribution of the local magnetization for a Spin 1 chain with open boundary conditions, showing Spin 1/2 edge modes.
Results for the Hubbard model are available here:
$$
H = -t\sum \sum_{i\sigma} c^\dagger_{i\sigma} c_{i+1\sigma} +
U \sum_{i} n_{i\uparrow} n_{i\downarrow}
$$
In the top-left panel we compare the energy per site
The regression tests perform complete DMRG runs from scratch and compare the resulting energies, entanglement entropies, local observables, nearest-neighbor correlations, and long-range correlations against trusted reference data. They currently cover the spin-$1/2$ and one-band Hubbard drivers.
After configuring and building the project, run the complete test suite from the build directory with:
cmake ..
cmake --build .
ctest --output-on-failureWhen using a Makefile generator, the last two commands can equivalently be written as:
cmake ..
make
make testIndividual DMRG regression tests can be selected with:
ctest -R dmrg_spin_regression --output-on-failure
ctest -R dmrg_hubbard_regression --output-on-failureUse ctest -N to list the registered tests and ctest -V to display their complete output. Test runs are performed in isolated working directories under:
build/test/work/Spin1d_SU2
build/test/work/Hubbard1d
The directories under test/reference contain the manual drivers, input files, and trusted numerical tables used to generate the reference data. These drivers are neither compiled nor executed by CMake. Once validated, the model input file and the resulting *.check tables must be copied manually to the corresponding directory under test/src. CMake stages the approved files from test/src into the isolated test working directories. After updating them, rerun cmake .. from the build directory before executing the tests again.
The code contains a detailed profiling of timing and parallel data communication for the crucial parts of the algorithm. The profiling output is regulated with the CMake configuration variable PROFILE.
The output is written to the file named full_profile_, with a suffix determined by the DMRGtype input variable. The profile results are appended to the file step-by-step separated by the pattern "# STEP: <LEFT.length>".
This one-line bash script can be used to split the profile output into several file, named profile_step_<counter>:
awk '/# STEP:/{close(f); f="profile_step_" ++c;next} {print>f;}' full_profile_iDMRG.out
There are a number of known issues with this code which, mostly for time reasons, we did not solve completely. Please report to the author or open an issue in GitHub.
[1] Large-Scale Implementation of the Density Matrix Renormalization Group Algorithm, J.Vance, MHPC Thesis (2017).
[2] EDIpack: A parallel exact diagonalization package for quantum impurity problems, A.Amaricci et al., Computer Physics Communications, 273, 108261 (April 2022).
[3] Massively parallel exact diagonalization of strongly correlated systems, A.Dolfen, Master Thesis, IFF Julich (2006).
[4] Density-matrix algorithms for quantum renormalization groups, S.White, Phys. Rev. B 48, 10345 (1993).
[5] Entanglement-informed distributed wavefunction approach to scalable quantum many-body systems, A.Amaricci, arXiv:2605.07621 (2026).
--
This program is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation, either version 3 of the License, or any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU LGPL for more details.
The software is provided as it is and can be read and copied, in agreement with the Terms of Service of GITHUB.
You should have received a copy of the GNU LGPL along with this program. If not, see http://www.gnu.org/licenses/.