# Basic Usage

In this documentation, we explain how to use `block2`

as an “executable”.
The input parameters are provided in a formatted input file.
The input file format used in `block2`

is highly compatible to the `StackBlock`

format.
For the most cases, the `StackBlock`

input/configuration file (“dmrg.conf”) can be directly understood by `block2`

,
but `block2`

also has some important extension for the keywords.

The information provided below is analogous to the corresponding `StackBlock`

documentation,
since the same input file format is used. However, the output format of `block2`

can be very different from that of `StackBlock`

.

## Preparation

If `block2`

is installed using `pip install block2`

, one can run a DMRG calculation using the following command:

```
block2main dmrg.conf > dmrg.out
```

Otherwise, for manual installation, please first compile the code according to
Installation with `cmake`

option `-DBUILD_LIB=ON`

(and other necessary options).
The following python script is used as the “block2 executable”:

```
${BLOCK2HOME}/pyblock2/driver/block2main
```

where `${BLOCK2HOME}`

is the `block2`

root directory. The `build`

directory under `block2`

root directory should be in `PYTHONPATH`

. You can add the following line in your environemnt
(such as `~/.bashrc`

) or submission script:

```
export PYTHONPATH=${BLOCK2HOME}/build:${PYTHONPATH}
```

Then you can run a DMRG calculation using the following command:

```
${BLOCK2HOME}/pyblock2/driver/block2main dmrg.conf > dmrg.out
```

where `dmrg.conf`

is the input file and `dmrg.out`

is the output file.

To run a DMRG calculation with MPI parallelization, please use the following command:

```
mpirun --bind-to core --map-by ppr:${SLURM_TASKS_PER_NODE}:node:pe=${OMP_NUM_THREADS} \
python -u ${BLOCK2HOME}/pyblock2/driver/block2main dmrg.conf > dmrg.out
```

where `${SLURM_TASKS_PER_NODE}`

is the number of mpi processes in each node.
`${OMP_NUM_THREADS}`

is the number of threads (CPU cores) used by each mpi process.
When executed in multiple nodes, a global scratch space (network file system) is required.

## Integral Generation

In the following we will use the C_{2} molecule to demonstrate the `block2`

features.
Integrals and orbitals should be supplied externally in the Molpro’s FCIDUMP format.
The integral file for C_{2} can be found in `${BLOCK2HOME}/data/C2.CAS.PVDZ.FCIDUMP.ORIG`

or
generated using the following script (only the RHF case is required):

```
from pyscf import gto, scf, mcscf
from pyblock2._pyscf.ao2mo import integrals as itg
from pyblock2.driver.core import DMRGDriver, SymmetryTypes
mol = gto.M(atom='C 0 0 0; C 0 0 1.2425', basis='ccpvdz', symmetry='d2h')
# RHF case (for spin-adapted / non-spin-adapted DMRG)
mf = scf.RHF(mol).run()
mc = mcscf.CASCI(mf, 26, 8)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf, mc.ncore, mc.ncas, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SU2)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
driver.write_fcidump(h1e, g2e, ecore=ecore, filename='./FCIDUMP', pg="d2h", h1e_symm=True)
# UHF case (for non-spin-adapted DMRG only)
mf = scf.UHF(mol).run()
mc = mcscf.UCASCI(mf, 26, 8)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_uhf_integrals(mf, mc.ncore[0], mc.ncas, g2e_symm=8)
driver = DMRGDriver(scratch="./tmp", symm_type=SymmetryTypes.SZ)
driver.initialize_system(n_sites=ncas, n_elec=n_elec, spin=spin, orb_sym=orb_sym)
driver.write_fcidump(h1e, g2e, ecore=ecore, filename='./FCIDUMP.UHF', pg="d2h", h1e_symm=True)
```

Alternatively, the integral file can be generated using the `pyscf/dmrgscf`

interface:

```
from pyscf import gto, scf, mcscf, dmrgscf
import os
dmrgscf.settings.BLOCKEXE = os.popen("which block2main").read().strip()
dmrgscf.settings.MPIPREFIX = ''
mol = gto.M(atom='C 0 0 0; C 0 0 1.2425', basis='ccpvdz', symmetry=1)
mf = scf.RHF(mol).run()
mc = mcscf.CASCI(mf, 26, 8)
mc.fcisolver = dmrgscf.DMRGCI(mol)
mc.canonicalization = False
dmrgscf.dryrun(mc)
```

Note

Please see DMRGSCF (pyscf) for the instruction for the installation of `pyscf/dmrgscf`

.

## Ground State Energy

The following input file can be used to compute the ground state energy:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
hf_occ integral
schedule default
maxM 500
maxiter 30
num_thrds 16
```

Note

Note that the integral file `C2.CAS.PVDZ.FCIDUMP.ORIG`

should be in the working direcotry.
By default, the orbitals will be reordered using the `fiedler`

method. One can optionally
add the keyword `noreorder`

to avoid orbital reordering.

`num_thrds`

indicates the number of OpenMP threads (shared-memory parallelism) to use.

`hf_occ integral`

has no effects in `block2`

, but it is required in `StackBlock`

.
If this line appear, `block2main`

will try to write some output files in a stackblock-compatible format.

By default, the calculation will be done in the spin-adapted mode, which is the most efficient.
One can optionally add the keyword `nonspinadapted`

to use the non-spin-adapted mode.

The keyword `prefix <scratch dir>`

can be used to set a folder for storing scratch files.
If running in a HPC supercomputer, it is highly recommended to use the high IO speed scratch space
(instead of the “home” storage) to achieve high performance.

Lines start with `!`

in the input file will be ignored. [1]

D_{2h} point group is enabled by `sym d2h`

.
The keywords `schedule default`

and `maxM`

sets the default sweep schedule and
the maximum number of renormalized states kept during the sweep, respectively.
`block2`

will then automatically set a sweep schedule as well as the defaults for various convergence thresholds.

The mps bond dimensions, sweep energies and the associated maximum discarded weights can be extracted by grepping the output `dmrg.out`

.

```
$ grep Bond dmrg.out
Sweep = 0 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-03 | Dav threshold = 1.00e-04
Sweep = 1 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-03 | Dav threshold = 1.00e-04
Sweep = 2 | Direction = forward | Bond dimension = 250 | Noise = 1.00e-03 | Dav threshold = 1.00e-04
Sweep = 3 | Direction = backward | Bond dimension = 250 | Noise = 1.00e-03 | Dav threshold = 1.00e-04
... ...
Sweep = 16 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-06
Sweep = 17 | Direction = backward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-06
Sweep = 0 | Direction = forward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-06
Sweep = 1 | Direction = backward | Bond dimension = 500 | Noise = 0.00e+00 | Dav threshold = 1.00e-06
$ grep DW dmrg.out
Time elapsed = 1.678 | E = -75.4879935448 | DW = 1.39e-05
Time elapsed = 2.936 | E = -75.6007921322 | DE = -1.13e-01 | DW = 9.88e-06
Time elapsed = 4.203 | E = -75.6367659659 | DE = -3.60e-02 | DW = 9.25e-05
Time elapsed = 5.750 | E = -75.6373954252 | DE = -6.29e-04 | DW = 3.91e-05
... ...
Time elapsed = 38.782 | E = -75.7283521752 | DE = -3.48e-05 | DW = 5.24e-06
Time elapsed = 41.169 | E = -75.7283676788 | DE = -1.55e-05 | DW = 5.28e-06
Time elapsed = 2.009 | E = -75.7283421257 | DW = 4.18e-17
Time elapsed = 4.158 | E = -75.7283421257 | DE = -2.84e-14 | DW = 2.47e-16
```

Note that in the last two sweeps (in default schedule) the 1-site algorithm is used. As a result, the discarded weights are nearly zero.

If you set `outputlevel 1`

in the input file, only essential information will be
printed and the `grep`

step can be skipped.

## Targeting States

You can target the states distinguished by the number of electrons `nelec`

,
the total spin `spin`

and the point-group symmetry of the state `irrep`

.

The following input file computes the energy for a single B_{1g} state in D_{2h} point group:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 4
hf_occ integral
schedule default
maxM 500
maxiter 30
```

Note

In D_{2h} point group, `irrep`

can be A_{1g} (1), B_{3u} (2),
B_{2u} (3), B_{1g} (4), B_{1u} (5), B_{2g} (6), B_{3g} (7), A_{1u} (8).

This will generate the following output:

```
$ grep DW dmrg.out
Time elapsed = 1.983 | E = -75.5422510106 | DW = 1.08e-05
Time elapsed = 3.580 | E = -75.6245880097 | DE = -8.23e-02 | DW = 9.97e-06
Time elapsed = 5.376 | E = -75.6366528654 | DE = -1.21e-02 | DW = 9.13e-05
Time elapsed = 7.172 | E = -75.6374064699 | DE = -7.54e-04 | DW = 4.03e-05
... ...
Time elapsed = 38.611 | E = -75.6389586629 | DE = -2.48e-05 | DW = 2.01e-06
Time elapsed = 40.981 | E = -75.6389699555 | DE = -1.13e-05 | DW = 2.05e-06
Time elapsed = 2.029 | E = -75.6389630224 | DW = 5.58e-15
Time elapsed = 4.106 | E = -75.6389632670 | DE = -2.45e-07 | DW = 2.40e-16
```

## State-Averaged Calculation

In the state-averaged DMRG algorithm, more than one state can be targeted in one calculation.
The states being calculated can have the same or different `nelec`

, `spin`

or `irrep`

.
Multiple values can be given for the above keywords. [1]
The number of states (roots) and the weight of each state can be specified using keywords
`nroots`

and `weights`

, respectively.
`block2`

will then try to find the low energy states within the space of targets formed
by all combintaions of the given values of `nelec`

, `spin`

and `irrep`

.

Note

In `StackBlock`

, state-averaged calculation can only be done for states with the same
`nelec`

, `spin`

and `irrep`

. In `block2`

, targetting multiple `nelec`

, `spin`

or `irrep`

may cause the calculation hard to converge to the lowest energy states. Typically,
one needs larger `nroots`

than the number of states actually needed, to make sure that
the low energy states are converged.

For normal non-state-averaged calculation, namely, when `nroots`

is 1, you can also target
multiple `nelec`

, `spin`

or `irrep`

.

The following input file performs state-averged DMRG for two A_{1g} states in D_{2h} point group:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
nroots 2
weights 0.5 0.5
hf_occ integral
schedule default
maxM 500
maxiter 30
```

This will generate the following output:

```
$ grep DW dmrg.out
Time elapsed = 3.257 | E[ 2] = -75.5019604920 -75.4800275143 | DW = 1.54e-05
Time elapsed = 5.109 | E[ 2] = -75.5980474127 -75.5776457885 | DE = -9.76e-02 | DW = 1.98e-05
Time elapsed = 6.854 | E[ 2] = -75.6711500018 -75.6363593637 | DE = -5.87e-02 | DW = 1.86e-04
Time elapsed = 8.635 | E[ 2] = -75.6717525884 -75.6368970346 | DE = -5.38e-04 | DW = 1.35e-04
Time elapsed = 45.946 | E[ 2] = -75.7279558636 -75.6386525742 | DE = -3.41e-05 | DW = 2.49e-05
Time elapsed = 48.491 | E[ 2] = -75.7279954715 -75.6386699048 | DE = -1.73e-05 | DW = 1.67e-05
Time elapsed = 2.215 | E[ 2] = -75.7279403993 -75.6386251036 | DW = 1.77e-05
Time elapsed = 4.338 | E[ 2] = -75.7279224367 -75.6386152528 | DE = 9.85e-06 | DW = 8.35e-06
```

## State-Specific Calculation

### Orthogonalization Approach

The state-specific calculation can be done as a restart calculation which assumes that a previous
state-averaged DMRG calculation has been converged. The state-specific DMRG calculation then reads the MPS
from scratch folder and refines them for each root separately.
The state-specific DMRG calculation can be done with any of `onedot`

, `twodot`

or `twodot_to_onedot`

(default)
keywords. [1]

Note

In `StackBlock`

, state-specific calculation can only be done with `onedot`

.

A state-specific DMRG calculation for two A_{1g} states in D_{2h} point group consists of two steps.

First, using the input file given in the previous section to obtain the state-averaged MPSs (in the scratch folder).

Second, the state-specific DMRG calculation can be performed by setting the keyword

`statespecific`

. The MPSs from the previous DMRG calculation will be read from the scratch folder. The following input file can be used for this step:sym d2h orbitals C2.CAS.PVDZ.FCIDUMP.ORIG nelec 8 spin 0 irrep 1 nroots 2 weights 0.5 0.5 statespecific hf_occ integral schedule default maxM 500 maxiter 30

This will generate the following output:

```
$ grep Energy dmrg.out
DMRG Energy for root 0 = -75.728342642601376
DMRG Energy for root 1 = -75.638959372610813
```

Sometimes, the orthogonalization approach can be unstable and when computing the exciated state
it may fall back to the ground state. Adding the keyword `onedot`

for the second step can alleviate this problem.

### Level Shift Approach

The second step of the above can also be done with the level shift approach, by changing Hamiltonian from \(\hat{H}\) to \(\hat{H} + \sum_i w_i |\phi_i\rangle \langle \phi_i|\). Normally, the weights \(w_i\) are positive and they should be larger than the energy gap.

The following input file can be used for the second step:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
nroots 2
weights 0.5 0.5
statespecific
proj_weights 5 5
hf_occ integral
schedule default
maxM 500
maxiter 30
```

This will generate the following output:

```
$ grep Energy dmrg.out
DMRG Energy for root 0 = -75.728341047222145
DMRG Energy for root 1 = -75.638958637510370
```

### Without State-Average

The excited MPS and energies can also be obtained without performing a state-averaged calculation as the first step. Instead, we can do several DMRG, and each time projecting out MPSs from all previous DMRG.

Note

It is recommended to use `noreorder`

or fixed manual orbital reordering for this approach.
Otherwise, one should carefully check that the orbital reordering in all DMRG calculations are the same.

We first get the ground state using the following input file `dmrg-1.conf`

:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
schedule default
maxM 500
maxiter 30
mps_tags KET1
```

After this is finished, we compute the first excited state using the following input file `dmrg-2.conf`

:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
schedule default
maxM 500
maxiter 30
mps_tags KET2
proj_mps_tags KET1
proj_weights 5
```

Then we compute the second excited state using the following input file `dmrg-3.conf`

:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
schedule default
maxM 500
maxiter 30
mps_tags KET3
proj_mps_tags KET1 KET2
proj_weights 5 5
```

And so on.

This will generate the following output:

```
$ grep Energy dmrg-*.out
dmrg-1.out:DMRG Energy = -75.728342508616663
dmrg-2.out:DMRG Energy = -75.638961566176221
dmrg-3.out:DMRG Energy = -75.629597871820607
dmrg-4.out:DMRG Energy = -75.467766576734363
dmrg-5.out:DMRG Energy = -75.350470798772307
dmrg-6.out:DMRG Energy = -75.312672909521751
```

### Mixed with State-Average

The above approach can also be used together with the state-average approach. Namely, we can first compute the two lowest
states, then we compute the next three lowest states, by projecting out the two lowest states.
The MPS to be projected must not be in state-averaged format, so we need to use the `split_states`

keyword to break
state-averaged MPS into individual MPSs, so that they can be used for projection in the subsequent calculations.

Currently, this type of state-average calulcation cannot be used together with multiple targets.

We first get the two lowest states using the following input file `dmrg-1.conf`

:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
nroots 2
weights 0.5 0.5
schedule default
maxM 500
maxiter 30
mps_tags KET
copy_mps
split_states
```

After this is finished, we compute the next three states using the following input file `dmrg-2.conf`

:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
nroots 3
weights 0.5 0.5 0.5
schedule default
maxM 500
maxiter 30
mps_tags EXKET
proj_mps_tags KET-0 KET-1
proj_weights 5 5
copy_mps
split_states
```

After this is finished, we compute the next one state using the following input file `dmrg-3.conf`

:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
schedule default
maxM 500
maxiter 30
mps_tags EXXKET
proj_mps_tags KET-0 KET-1 EXKET-0 EXKET-1 EXKET-2
proj_weights 5 5 5 5 5
```

This will generate the following output:

```
$ grep DW dmrg-1.out | tail -1
Time elapsed = 5.461 | E[ 2] = -75.7279224622 -75.6386156808 | DE = 9.32e-06 | DW = 8.33e-06
$ grep DW dmrg-2.out | tail -1
Time elapsed = 13.165 | E[ 3] = -75.6290377907 -75.4669665917 -75.3494878435 | DE = 8.63e-07 | DW = 8.45e-05
$ grep DW dmrg-3.out | tail -1
Time elapsed = 8.651 | E = -75.3126745298 | DE = -6.24e-07 | DW = 3.79e-15
```

## n-Particle Reduced Density Matrix

The 1-, 2-, 3-, and 4-particle DMRG reduced density matrix for a particular state can be calculated using
the keywords `onepdm`

, `twopdm`

, `threepdm`

and `fourpdm`

.
The reduced density matrix calculation can be done with either `onedot`

or `twodot`

keywords. [1]

Note

Most of the time, only `onedot`

density matrix calculation makes sense, since the MPS should not change
during the sweep.

Density matrices of the \(n\)-th state are calculated and stored in a `numpy`

binary file
named `1pdm-n-n.npy`

, `2pdm-n-n.npy`

, `3pdm-n-n.npy`

, etc. (in the scratch folder), respectively,
starting with `n = 0`

.
If there is only one root, the files are named `1pdm.npy`

, `2pdm.npy`

, `3pdm.npy`

, etc. respectively.

The following input file computes the energy and 2-particle density matrix for the ground state:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
schedule default
maxM 500
maxiter 30
twopdm
num_thrds 16
```

The 2-particle density matrix file can be loaded using the following python script:

```
>>> import numpy as np
>>> _2pdm = np.load('./nodex/2pdm.npy')
>>> print(_2pdm.shape)
(3, 26, 26, 26, 26)
```

The following input file computes the energy and 2-particle density matrix for two
state-averaged A_{1g} states:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
nroots 2
weights 0.5 0.5
schedule default
maxM 500
maxiter 30
twopdm
num_thrds 16
```

The 2-particle density matrix file for the first state can be loaded using the following python script:

```
>>> import numpy as np
>>> n = 0
>>> _2pdm = np.load('./nodex/2pdm-%d-%d.npy' % (n, n))
>>> print(_2pdm.shape)
(3, 26, 26, 26, 26)
```

The 1-particle density matrix (in both the non-spin-adapted and spin-adapted mode) is stored as an array with the shape \([2,n,n]\),
where `n`

is the number of spatial orbitals, and the two components with indicies \([:,a,b]\) are for
\(\langle a^\dagger_{a\alpha} a_{b\alpha} \rangle\), and \(\langle a^\dagger_{a\beta} a_{b\beta} \rangle\),
respectively.

The 2-particle density matrix (in both the non-spin-adapted and spin-adapted mode) is stored as an array with the shape \([3,n,n,n,n]\), where the three components with indicies \([:,a,b,c,d]\) are for \(\langle a^\dagger_{a\alpha} a^\dagger_{b\alpha} a_{c\alpha} a_{d\alpha} \rangle\), \(\langle a^\dagger_{a\alpha} a^\dagger_{b\beta} a_{c\beta} a_{d\alpha} \rangle\), and \(\langle a^\dagger_{a\beta} a^\dagger_{b\beta} a_{c\beta} a_{d\beta} \rangle\), respectively.

The 3-particle density matrix in the spin-adapted mode is stored as the spin-traced format with the shape \([n,n,n,n,n,n]\), defined as

The 3-particle density matrix in the non-spin-adapted mode is stored as an array with the shape \([4,n,n,n,n,n,n]\), where the four components with indicies \([:,a,b,c,d,e,f]\) are for \(\langle a^\dagger_{a\alpha} a^\dagger_{b\alpha} a^\dagger_{c\alpha} a_{d\alpha}a_{e\alpha} a_{f\alpha} \rangle\), \(\langle a^\dagger_{a\alpha} a^\dagger_{b\alpha} a^\dagger_{c\beta} a_{d\beta}a_{e\alpha} a_{f\alpha} \rangle\), \(\langle a^\dagger_{a\alpha} a^\dagger_{b\beta} a^\dagger_{c\beta} a_{d\beta}a_{e\beta} a_{f\alpha} \rangle\), and \(\langle a^\dagger_{a\beta} a^\dagger_{b\beta} a^\dagger_{c\beta} a_{d\beta}a_{e\beta} a_{f\beta} \rangle\), respectively.

The 4-particle density matrix in the spin-adapted mode is stored as the spin-traced format with the shape \([n,n,n,n,n,n,n,n]\), defined as

The 4-particle density matrix in the non-spin-adapted mode is stored as an array with the shape \([5,n,n,n,n,n,n,n,n]\), where the five components with indicies \([:,a,b,c,d,e,f,g,h]\) are for \(\langle a^\dagger_{a\alpha} a^\dagger_{b\alpha} a^\dagger_{c\alpha} a^\dagger_{d\alpha} a_{e\alpha} a_{f\alpha} a_{g\alpha} a_{h\alpha} \rangle\), \(\langle a^\dagger_{a\alpha} a^\dagger_{b\alpha} a^\dagger_{c\alpha} a^\dagger_{d\beta} a_{e\beta} a_{f\alpha} a_{g\alpha} a_{h\alpha} \rangle\), \(\langle a^\dagger_{a\alpha} a^\dagger_{b\alpha} a^\dagger_{c\beta} a^\dagger_{d\beta} a_{e\beta} a_{f\beta} a_{g\alpha} a_{h\alpha} \rangle\), \(\langle a^\dagger_{a\alpha} a^\dagger_{b\beta} a^\dagger_{c\beta} a^\dagger_{d\beta} a_{e\beta} a_{f\beta} a_{g\beta} a_{h\alpha} \rangle\), and \(\langle a^\dagger_{a\beta} a^\dagger_{b\beta} a^\dagger_{c\beta} a^\dagger_{d\beta} a_{e\beta} a_{f\beta} a_{g\beta} a_{h\beta} \rangle\), respectively.

In the general spin orbital mode (with the keyword `use_general_spin`

), the 1-, 2-, 3-, and 4-particle density matrices are stored
with the shape \([1,n,n]\), \([1,n,n,n,n]\), \([1,n,n,n,n,n,n]\), and \([1,n,n,n,n,n,n,n,n]\) respectively,
where `n`

is the number of spin orbitals. The content is the expectation value for \(\langle a^\dagger_{a} a_{b} \rangle\),
\(\langle a^\dagger_{a} a^\dagger_{b} a_{c} a_{d} \rangle\), \(\langle a^\dagger_{a} a^\dagger_{b} a^\dagger_{c} a_{d} a_{e} a_{f} \rangle\),
and \(\langle a^\dagger_{a} a^\dagger_{b} a^\dagger_{c} a^\dagger_{d} a_{e} a_{f} a_{g} a_{h} \rangle\), respectively.

## n-Particle Transition Reduced Density Matrix

The 1-, 2-, 3- and 4-particle DMRG transition density matrix can be calculated using
the keywords `tran_onepdm`

, `tran_twopdm`

, `tran_threepdm`

and `tran_fourpdm`

.

Transition density matrices between the \(m\)-th (bra) and \(n\)-th (ket) states are calculated and stored in a `numpy`

binary file
named `1pdm-m-n.npy`

, `2pdm-m-n.npy`

, etc. (in the scratch folder), respectively, starting with `m = n = 0`

.

The following input file computes the 2-particle transition density matrix for two
state-averaged A_{1g} states:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
nroots 2
weights 0.5 0.5
schedule default
maxM 500
maxiter 30
tran_twopdm
num_thrds 16
```

Note

There can be a overall undetermined +1/-1 factor in Transition density matrices due to the relative phase in two MPSs.

The following input file computes the state-specific 2-particle transition density matrix for two
refined A_{1g} states:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
nroots 2
weights 0.5 0.5
statespecific
schedule default
maxM 500
maxiter 30
tran_twopdm
num_thrds 16
```

The transition density matrices between states with different point group irreducible representations are also available by simply
adding the keyword `tran_twopdm`

after the corresponding multi-target state-averaged calculation. [1]

## Restart DMRG Energy Calculation

DMRG energy calculations can be restarted, using the MPS (stored in scratch folder) generated in the previous calculation,
by specifying the keyword `fullrestart`

.
If the previous calulcation stopped during the middle of a sweep, it will be restarted from the middle of a sweep.

Alternatively, the user can also set a directory for storing MPS after each sweep using the keyword `restart_dir`

. [1]
When restarting, the MPS data and `mps_info.bin`

in the scratch folder should be copied from the `restart_dir`

to the
scartch folder of the restarting calculation.

The keyword `restart_dir_per_sweep`

can be used to save a copy of MPS for each sweep. The MPS from different sweeps will
be put into different folders (by adding suffix to the given direcotry).

You may need to change the (custom) scheudle in the input file so that the sweeps (with smaller bond dimension) finished in previous calculations will not be repeated, when you are restarting an interrupted calculation.

The following input file restarts an interrupted calculation:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
hf_occ integral
schedule default
maxM 500
maxiter 30
fullrestart
```

## Load MPS for Density Matrix Calculation

The density matrix and transition density matrix calculation can be carried out separately, by restarting from an existing MPS, state-averged MPSs or state-specific MPSs (stored in scartch folder from a previous DMRG energy calculation).

Assuming a previous ground-state energy calculation has been finished, the following input file computes the 2-particle density matrix for the ground-state (loaded from scratch folder):

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
hf_occ integral
schedule default
maxM 500
maxiter 30
restart_twopdm
```

Assuming a previous state-averaged energy calculation has been finished,
the following input file computes the 2-particle transition density matrix for two
state-averaged A_{1g} states (loaded from scratch folder):

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
nroots 2
weights 0.5 0.5
hf_occ integral
schedule default
maxM 500
maxiter 30
restart_tran_twopdm
```

Now we explain how to compute 2-particle transition density matrix for bra and ket states
belonging to different point group irreducible representations.
We consider the A_{1g} (bra) and B_{3u} (ket) states.

The following input file computes the energy for a single B_{3u} state in D_{2h} point group.
The keyword `mps_tags`

can be used to assign a tag to the mps for later reference: [1]

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 2
hf_occ integral
schedule default
maxM 500
maxiter 30
mps_tags KET
```

The following input file computes the energy for a single A_{1g} state in D_{2h} point group:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
hf_occ integral
schedule default
maxM 500
maxiter 30
mps_tags BRA
```

The output looks like the following:

```
$ grep Energy dmrg-1.out
DMRG Energy = -75.675393353797631
$ grep Energy dmrg-2.out
DMRG Energy = -75.728342388135175
```

The following input file computes the 2-particle transition density matrix for the two states:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
mps_tags BRA KET
hf_occ integral
schedule default
maxM 500
maxiter 30
restart_tran_twopdm
```

Note that in the above input file, keywords such as `nelec`

, `spin`

, `irrep`

, and `nroots`

will be unimportant.
The keyword `mps_tags`

lists the tags for the MPSs that should be loaded. [1]

## Diagonal 2-Particle Density Matrix

Since the full two-particle density matrix calculation can be expensive for some systems,
it is possible to calculate only the diagonal parts, which is much cheaper, using the keywords
`restart_diag_twopdm`

or `diag_twopdm`

. [1]
The time cost for diagonal 2pdm is roughly 2 times of the cost of 1pdm.

Note that `diag_twopdm`

implies `onepdm`

and `correlation`

. The diagonal 2pdm is defined as:

The computed diagonal 2pdm will be stored as `e_pqqp.npy`

and `e_pqpq.npy`

in scratch folder.

If one also computed the full 2pdm using the keyword `twopdm`

or `restart_twopdm`

,
we can verify that its diagonal part matches the `e_pqqp.npy`

and `e_pqpq.npy`

obtained here:

```
>>> import numpy as np
>>> _2pdm = np.load('./nodex/2pdm.npy')
>>> print(_2pdm.shape)
(3, 26, 26, 26, 26)
>>> _e_pqqp = np.load('./nodex/e_pqqp.npy')
>>> _e_pqpq = np.load('./nodex/e_pqpq.npy')
>>> _2pdm_spat = _2pdm[0] + 2 * _2pdm[1] + _2pdm[2]
>>> _2pdm_spat_pqqp = np.einsum('pqqp->pq', _2pdm_spat)
>>> _2pdm_spat_pqpq = np.einsum('pqpq->pq', _2pdm_spat)
>>> print(np.linalg.norm(_e_pqqp - _2pdm_spat_pqqp))
3.28666776770176e-14
>>> print(np.linalg.norm(_e_pqpq - _2pdm_spat_pqpq))
1.6947732597975102e-14
```

## Custom Sweep Schedule

The sweep schedule defines number of the renormalized states \(M\) kept ,
the convergence threshold for Davidson algorithm (in the unit of norm^{2}), and the noise
(in the unit of norm^{2}) in successive DMRG sweeps.
For finer control over the sweeps, customized sweep schedule should be used.

The following input file computes the ground state energy using a custom sweep schedule:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
hf_occ integral
schedule
0 100 1E-4 1E-3
4 250 1E-4 1E-3
8 400 1E-5 1E-4
10 600 1E-6 1E-5
12 800 1E-7 1E-6
14 1000 1E-8 1E-7
16 1000 1E-8 0E+0
end
twodot_to_onedot 18
maxiter 100
sweep_tol 1E-9
```

In the above input file, `twodot_to_onedot`

specifies the sweep at which the switch is made from
a 2-site to a 1-site DMRG algorithm (counting from 0). `maxiter`

gives the maximum number of sweep
iterations to be performed. `sweep_tol`

gives the final tolerance on the DMRG energy, and is analogous
to an energy convergence threshold in other quantum chemistry methods.

In the above input file, between `schedule`

and `end`

each line has four values. They are corresponding
to starting sweep iteration (counting from zero), MPS bond dimension, tolerance for the Davidson iteration,
and noise, respectively. Starting sweep iteration is the sweep iteration in which the given parameters
in the line should take effect.

This will generate the following output:

```
$ grep DW dmrg.out
Time elapsed = 1.686 | E = -74.1599100997 | DW = 4.86e-05
Time elapsed = 3.332 | E = -74.6555553068 | DE = -4.96e-01 | DW = 7.28e-05
Time elapsed = 4.461 | E = -75.6224601188 | DE = -9.67e-01 | DW = 1.55e-04
Time elapsed = 5.648 | E = -75.6302268887 | DE = -7.77e-03 | DW = 1.26e-04
Time elapsed = 7.491 | E = -75.6347292246 | DE = -4.50e-03 | DW = 6.46e-05
Time elapsed = 10.732 | E = -75.6367873793 | DE = -2.06e-03 | DW = 2.96e-05
Time elapsed = 13.383 | E = -75.6372588510 | DE = -4.71e-04 | DW = 1.01e-04
Time elapsed = 16.138 | E = -75.6375874124 | DE = -3.29e-04 | DW = 3.83e-05
Time elapsed = 20.541 | E = -75.6687725683 | DE = -3.12e-02 | DW = 8.76e-06
Time elapsed = 26.404 | E = -75.7265879915 | DE = -5.78e-02 | DW = 9.21e-06
Time elapsed = 36.001 | E = -75.7282887562 | DE = -1.70e-03 | DW = 3.43e-06
Time elapsed = 47.351 | E = -75.7283943399 | DE = -1.06e-04 | DW = 3.04e-06
Time elapsed = 64.673 | E = -75.7284858001 | DE = -9.15e-05 | DW = 1.24e-06
Time elapsed = 86.412 | E = -75.7285031554 | DE = -1.74e-05 | DW = 1.21e-06
Time elapsed = 118.443 | E = -75.7285302492 | DE = -2.71e-05 | DW = 4.82e-07
Time elapsed = 158.894 | E = -75.7285335786 | DE = -3.33e-06 | DW = 5.44e-07
Time elapsed = 176.071 | E = -75.7285376489 | DE = -4.07e-06 | DW = 5.73e-07
Time elapsed = 191.672 | E = -75.7285377336 | DE = -8.46e-08 | DW = 5.76e-07
Time elapsed = 10.790 | E = -75.7285342605 | DW = 1.47e-16
Time elapsed = 21.186 | E = -75.7285342992 | DE = -3.87e-08 | DW = 3.21e-14
Time elapsed = 31.924 | E = -75.7285343224 | DE = -2.32e-08 | DW = 3.07e-17
Time elapsed = 42.348 | E = -75.7285343375 | DE = -1.51e-08 | DW = 8.17e-15
Time elapsed = 53.073 | E = -75.7285343475 | DE = -9.98e-09 | DW = 4.35e-17
Time elapsed = 63.362 | E = -75.7285343571 | DE = -9.58e-09 | DW = 6.64e-16
Time elapsed = 73.965 | E = -75.7285343630 | DE = -5.87e-09 | DW = 3.96e-17
Time elapsed = 84.094 | E = -75.7285343661 | DE = -3.17e-09 | DW = 1.14e-16
Time elapsed = 94.525 | E = -75.7285343678 | DE = -1.71e-09 | DW = 1.34e-16
Time elapsed = 104.658 | E = -75.7285343721 | DE = -4.29e-09 | DW = 2.45e-16
Time elapsed = 114.925 | E = -75.7285343746 | DE = -2.44e-09 | DW = 1.38e-16
Time elapsed = 124.710 | E = -75.7285343763 | DE = -1.76e-09 | DW = 3.03e-16
Time elapsed = 135.115 | E = -75.7285343763 | DE = 5.68e-14 | DW = 2.24e-17
```

## Sweep Energy Extrapolation

In practice the sweep energy converges almost linearly as a function of the “maximum discarded weight”. Therefore, it is convenient to use the “maximum discarded weight” quantity as an estimate of the error of the DMRG calculation. It is recommended to use the 2-site algorithm for energy extrapolation since the 2-site DMRG wavefunction provides additional variational freedom over the 1-site DMRG wavefunction. A strong deviation from a linear function (e.g. a plateau behavior followed by a sudden drop of the energy as a function of discarded weight) indicates that the DMRG was stuck in a local minimum.

One can use restart a converged DMRG calculation with a “reverse schedule” to generate data for energy extrapolation. This can guarentee that the energy for each different MPS bond dimension is fully converged and not representing any local minima.

The following input file restarts the previous calculation using a custom reverse sweep schedule:

```
sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nelec 8
spin 0
irrep 1
hf_occ integral
twodot
schedule
0 800 1E-8 0E+0
4 600 1E-8 0E+0
8 400 1E-8 0E+0
12 200 1E-8 0E+0
end
maxiter 16
sweep_tol 0.0
fullrestart
```

This will generate the following output (`dmrg-2.out`

):

```
$ grep DW dmrg-2.out
Time elapsed = 12.597 | E = -75.7285358881 | DW = 1.75e-06
Time elapsed = 23.720 | E = -75.7285188420 | DE = 1.70e-05 | DW = 1.42e-06
Time elapsed = 33.955 | E = -75.7285186195 | DE = 2.23e-07 | DW = 1.35e-06
Time elapsed = 44.842 | E = -75.7285186529 | DE = -3.34e-08 | DW = 1.34e-06
Time elapsed = 52.432 | E = -75.7285113908 | DE = 7.26e-06 | DW = 4.98e-06
Time elapsed = 59.530 | E = -75.7284626837 | DE = 4.87e-05 | DW = 3.66e-06
Time elapsed = 66.036 | E = -75.7284622858 | DE = 3.98e-07 | DW = 3.49e-06
Time elapsed = 73.045 | E = -75.7284623697 | DE = -8.39e-08 | DW = 3.47e-06
Time elapsed = 77.523 | E = -75.7284421278 | DE = 2.02e-05 | DW = 1.71e-05
Time elapsed = 81.396 | E = -75.7282631341 | DE = 1.79e-04 | DW = 1.11e-05
Time elapsed = 85.001 | E = -75.7282618298 | DE = 1.30e-06 | DW = 1.02e-05
Time elapsed = 88.824 | E = -75.7282620286 | DE = -1.99e-07 | DW = 1.02e-05
Time elapsed = 91.267 | E = -75.7282077342 | DE = 5.43e-05 | DW = 1.04e-04
Time elapsed = 93.148 | E = -75.7270840401 | DE = 1.12e-03 | DW = 5.65e-05
Time elapsed = 95.144 | E = -75.7270844505 | DE = -4.10e-07 | DW = 5.01e-05
Time elapsed = 96.921 | E = -75.7270854757 | DE = -1.03e-06 | DW = 4.85e-05
```

Sweep energy extrapolation can be plotted using the following python script:

```
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats
fname = 'dmrg-2.out'
out = open(fname, 'r').readlines()
eners, dws = [], []
for l in out:
if "DW" in l:
eners.append(float(l.split()[7]))
dws.append(float(l.split()[-1]))
eners, dws = eners[3::4], dws[3::4]
reg = scipy.stats.linregress(dws, eners)
x_reg = np.array([0, 1E-4])
emin, emax = min(eners), max(eners)
de = emax - emin
plt.plot(x_reg, reg.intercept + reg.slope * x_reg, '--', linewidth=1, color='#5FA8AB')
plt.plot(dws, eners, 'o', color='#38686A', markerfacecolor='white', markersize=5)
plt.text(2E-6, emax, "$E(M=\\infty) = %.6f \pm %.6f \\mathrm{\\ Hartree}$" %
(reg.intercept, abs(reg.intercept - emin) / 5), color='#38686A', fontsize=12)
plt.text(2E-6, emax - de * 0.1, "$R^2 = %.6f$" % (reg.rvalue ** 2),
color='#38686A', fontsize=12)
plt.xlim((0, 5E-5))
plt.ylim((emin - de * 0.1, emax + de * 0.1))
plt.xlabel("Largest Discarded Weight")
plt.ylabel("Sweep Energy (Hartree)")
plt.subplots_adjust(left=0.16, bottom=0.1, right=0.95, top=0.95)
plt.savefig("extra.png", dpi=600)
```

Alternatively, the keyword `extrapolation`

can be added to the previous script,
so that the extrapolation energy will be printed and the figure named `extrapolation.png`

will be saved in the scartch folder.

The script will generate the following figure:

In the above script, we have used the largest discarded weights and associated sweep energies in the last sweep iteration of each bond dimension (\(M = 800, 600, 400, 200\)) to make linear regression. The extrapolated DMRG sweep energy is -75.728567 Hartree.