Input File: Advanced Usage

Orbital Rotation

In this calculation we illustrate how to compute the ground state MPS in the given set of orbitals, find the (new) DMRG natural orbitals, transform integrals to new orbitals, transform the ground state MPS to new orbitals, and finally evaluate the energy of the transformed MPS in the new orbitals to verify the quality of the transformed MPS.

First, we compute the energy and 1-particle density matrix for the ground state using the following input file:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 0
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

onepdm
irrep_reorder

Note that we use the keyword irrep_reorder to reorder the orbitals so that orbitals belonging to the same point group irrep are grouped together. This can make the orbital rotation more local.

The DMRG occupation number (in original ordering) will be printed at the end of the calculation:

$ grep OCC dmrg-1.out
DMRG OCC =   1.957 1.625 1.870 1.870 0.361 0.098 0.098 0.006 0.008 0.008 0.008 0.013 0.014 0.014 0.011 0.006 0.006 0.006 0.005 0.005 0.002 0.002 0.002 0.001 0.001 0.001
$ grep Energy dmrg-1.out
DMRG Energy =  -75.728467269121097

Second, we use the keyword nat_orbs to compute the natural orbitals. The value of the keyword nat_orbs specifies the filename for storing the rotated integrals (FCIDUMP). If no value is associated with the keyword nat_orbs, the rotated integrals will not be computed. The keyword nat_orbs can only be used together with restart_onepdm or onepdm, since natural orbitals are found by diagonalizing 1-particle density matrix.

The following input file is used for this step (it can also be combined with the previous calculation):

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 0
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

restart_onepdm
nat_orbs C2.NAT.FCIDUMP
nat_km_reorder
nat_positive_def
irrep_reorder

Where the optional keyword nat_km_reorder can be used to remove the artificial reordering in the natural orbitals using Kuhn-Munkres algorithm. The optional keyword nat_positive_def can be used to avoid artificial rotation in the logarithm of the rotation matrix, by make the rotation matrix quasi-positive-definite, with “quasi” in the sense that the rotation matrix is not Hermitian. The two options may be good for weakly correlated systems, but have limited effects for highly correlated systems (but for highly correlated systems it is also recommended to be used).

The occupation number in natural orbitals will be printed at the end of the calculation:

$ grep OCC dmrg-2.out
DMRG OCC =   1.957 1.625 1.870 1.870 0.361 0.098 0.098 0.006 0.008 0.008 0.008 0.013 0.014 0.014 0.011 0.006 0.006 0.006 0.005 0.005 0.002 0.002 0.002 0.001 0.001 0.001
REORDERED OCC =   1.957 0.002 0.361 0.006 0.013 0.008 0.002 0.006 0.011 0.001 0.006 1.625 0.008 1.870 0.005 0.098 0.001 0.014 0.005 1.870 0.008 0.001 0.014 0.098 0.006 0.002
NAT OCC =   0.000465 0.003017 0.006424 0.007848 0.360936 1.968407 0.000081 0.000916 0.001991 0.004082 0.015623 1.628182 0.003669 0.008706 1.870680 0.000424 0.002862 0.110463 0.003667 0.008705 1.870678 0.000424 0.002862 0.110480 0.006422 0.001989

With the optional keyword nat_km_reorder there will be an extra line:

REORDERED NAT OCC =   1.968407 0.000465 0.360936 0.006424 0.007848 0.003017 0.001991 0.000081 0.004082 0.000916 0.015623 1.628182 0.008706 1.870680 0.003669 0.110463 0.000424 0.002862 0.003667 1.870678 0.008705 0.000424 0.002862 0.110480 0.006422 0.001989

The rotation matrix for natural orbitals, the logarithm of the rotation matrix, and the occupation number in natural orbitals are stored as nat_rotation.npy, nat_kappa.npy, nat_occs.npy in scartch folder, respectively. In this example, the rotated integral is stored as C2.NAT.FCIDUMP in the working directory.

Third, we load the MPS in the old orbitals and transform it into the new orbitals. This is done using time evolution. The keyword delta_t is used to set a time step and indicate that this is a time evolution calculation. The keyword orbital_rotation is used to indicate that the operator (exponentiated) applied into the MPS should be the orbital rotation operator (constructed from nat_kappa.npy saved in the previous step).

Typically, a large bond dimension should be used depending how non-local the orbital rotation operator is. The target_t for orbital rotation is automatically set to 1.

The following input file is used for this step:

sym d2h

nelec 8
spin 0
irrep 1

schedule
    0 1000 0 0
end

mps_tags BRA
orbital_rotation
delta_t 0.05
outputlevel 1
noreorder

Note that noreorder must be used for orbital rotation. The orbital reordering in previous step has already been taken into account.

The keyword te_type can be used to set the time-evolution algorithm. The default is rk4, which is the original time-step-targeting (TST) method. Another possible choice is tdvp, which is the time dependent variational principle with the projector-splitting (TDVP-PS) algorithm.

The output looks like the following:

$ grep DW dmrg-3.out
Time elapsed =      2.263 | E =       0.0000000000 | Norm^2 =       0.9999999999 | DW = 1.76e-10
Time elapsed =      4.910 | E =      -0.0000000000 | Norm^2 =       0.9999999997 | DW = 1.43e-10
Time elapsed =      1.663 | E =      -0.0000000000 | Norm^2 =       0.9999999988 | DW = 4.46e-10
Time elapsed =      3.475 | E =       0.0000000000 | Norm^2 =       0.9999999983 | DW = 2.50e-10
... ...
Time elapsed =      3.011 | E =       0.0000000000 | Norm^2 =       0.9999999315 | DW = 1.04e-09
Time elapsed =      4.753 | E =       0.0000000000 | Norm^2 =       0.9999999284 | DW = 8.68e-10
Time elapsed =      1.786 | E =       0.0000000000 | Norm^2 =       0.9999999245 | DW = 1.07e-09
Time elapsed =      3.835 | E =       0.0000000000 | Norm^2 =       0.9999999213 | DW = 9.09e-10

Since in every time step an orthogonal transformation is applied on the MPS, the expectation value of the orthogonal transformation (printed as the energy expectation) calculated on the MPS should always be zero.

Note that largest discarded weight is 1.07e-09, and the norm of MPS is not far away from 1. So the transormation should be relatively accurate.

Finally, we calculate the energy expectation value using the transformed integral (C2.NAT.FCIDUMP) and the transformed MPS (stored in the scratch folder), using the following input file:

sym d2h
orbitals C2.NAT.FCIDUMP

nelec 8
spin 0
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

mps_tags BRA
restart_oh
restart_onepdm
noreorder

Note that noreorder must be used, since the MPS generated in the previous step is in unreordered natural orbitals. The keyword restart_oh will calculate the expectation value of the given Hamiltonian loaded from integrals on the MPS loaded from scartch folder.

We have the following output:

$ grep Energy dmrg-4.out
OH Energy =  -75.728457535820155

The difference compared to the energy generated in the first step DMRG Energy =  -75.728467269121097 is only 9.7E-6. One can increase the bond dimension in the evolution to make this closer to the value printed in the first step.

MPS Transform

The MPS can be copied and saved using another tag. For SU2 (spin-adapted) MPS, it can also be transformed to SZ (non-spin-adapted) MPS and saved using another tag.

Limitations:

• Total spin zero spin-adapted MPS can be transformed directly.

• For non-zero total spin, the spin-adapted MPS must be in singlet embedding format. See next section.

First, we compute the energy for the spin-adapted ground state using the following input file:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 0
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

irrep_reorder
mps_tags KET

The following script will read the spin-adapted MPS and tranform it to a non-spin-adapted MPS:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 0
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

irrep_reorder
mps_tags KET
restart_copy_mps ZKET
trans_mps_to_sz

Here the keyword restart_copy_mps indicates that the MPS will be copied, associated with a value indicating the new tag for saving the copied MPS. If the keyword trans_mps_to_sz is present, the MPS will be transformed to non-spin-adapted before being saved.

Finally, we calculate the energy expectation value using non-spin-adapted formalism and the transformed MPS (stored in the scratch folder), using the following input file:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 0
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

irrep_reorder
mps_tags ZKET
restart_oh
nonspinadapted

Some reference outputs for this example:

$ grep Energy dmrg-1.out
DMRG Energy =  -75.728467269121083
$ grep MPS dmrg-2.out
MPS =  KRRRRRRRRRRRRRRRRRRRRRRRRR 0 2
GS INIT MPS BOND DIMS =       1     3    10    35   120   263   326   500   500   500   500   500   500   500   500   500   500   500   498   500   407   219    94    32    10     3     1
$ grep 'MPS\|Energy' dmrg-3.out
MPS =  KRRRRRRRRRRRRRRRRRRRRRRRRR 0 2
GS INIT MPS BOND DIMS =       1     4    16    64   246   578   712  1114  1097  1102  1110  1121  1126  1130  1116  1111  1111  1107  1074  1103   895   444   186    59    16     4     1
OH Energy =  -75.728467269120898

We can see that the transformation from SU2 to SZ is nearly exact, and the required bond dimension for the SZ MPS is roughly two times of the SU2 bond dimension.

Singlet Embedding

For spin-adapted calculation with total spin not equal to zero, there can be some convergence problem even if in one-site algorithm. One way to solve this problem is to use singlet embedding. In StackBlock singlet embedding is used by default. In block2, by default singlet embedding is not used. If one adds the keyword singlet_embedding to the input file, the singlet embedding scheme will be used. For most total spin not equal to zero calculation, singlet embedding may be more stable. One cannot calculate transition density matrix between states with different total spins using singlet embedding. To do that one can translate the MPS between singlet embedding format and non-singlet-embedding format.

When total spin is equal to zero, the keyword singlet_embedding will not have any effect. If restarting a calculation, normally, the keyword singlet_embedding is not required since the format of the MPS can be automatically recognized.

For translating SU2 MPS to SZ MPS with total spin not equal to zero, the SU2 MPS must be in singlet embedding format.

First, we compute the energy for the spin-adapted with non-zero total spin using the following input file:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 2
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

irrep_reorder
mps_tags KET

The above input file indicates that singlet embedding is not used. The output is:

$ grep 'MPS = ' dmrg-1.out
MPS =  CCRRRRRRRRRRRRRRRRRRRRRRRR 0 2 < N=8 S=1 PG=0 >
$ grep Energy dmrg-1.out
DMRG Energy =  -75.423916647509742

Here the printed target quantum number of the MPS indicates that it is a triplet.

We can add the keyword singlet_embedding to do a singlet embedding calculation:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 2
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

irrep_reorder
mps_tags SEKET
singlet_embedding

When singlet embedding is used, the output is:

$ grep 'MPS = ' dmrg-2.out
MPS =  CCRRRRRRRRRRRRRRRRRRRRRRRR 0 2 < N=10 S=0 PG=0 >
$ grep Energy dmrg-2.out
DMRG Energy =  -75.423879916245895

Here the printed target quantum number of the MPS indicates that it is a singlet (including some ghost particles).

One can use the keywords trans_mps_to_singlet_embedding and trans_mps_from_singlet_embedding combined with restart_copy_mps or copy_mps to translate between singlet embedding and normal formats.

The following script transforms the MPS from singlet embedding to normal format:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 2
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

irrep_reorder
mps_tags SEKET
restart_copy_mps TKET
trans_mps_from_singlet_embedding

We can verify that the transformed non-singlet-embedding MPS has the same energy as the singlet embedding MPS:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 2
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

irrep_reorder
mps_tags TKET
restart_oh

With the outputs:

$ grep 'MPS = ' dmrg-4.out
MPS =  KRRRRRRRRRRRRRRRRRRRRRRRRR 0 2 < N=8 S=1 PG=0 >
$ grep Energy dmrg-4.out
OH Energy =  -75.423879916245824

The following script will read the spin-adapted singlet embedding MPS and tranform it to a non-spin-adapted MPS:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 2
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

irrep_reorder
mps_tags SEKET
restart_copy_mps ZKETM2
trans_mps_to_sz
resolve_twosz -2
normalize_mps

Here the keyword resolve_twosz indicates that the transformed SZ MPS will have projected spin 2 * SZ = -2. For this case since 2 * S = 2, the possible values for resolve_twosz are -2, 0, 2. If the keyword resolve_twosz is not given, an MPS with ensemble of all possible projected spins will be produced (which is often not very useful). Getting one component of the SU2 MPS means that the SZ MPS will not have the same norm as the SU2 MPS. If the keyword normalize_mps is added, the transformed SZ MPS will be normalized. The keyword normalize_mps can only have effect when trans_mps_to_sz is present.

Finally, we calculate the energy expectation value using non-spin-adapted formalism and the transformed MPS (stored in the scratch folder), using the following input file:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin -2
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

irrep_reorder
mps_tags ZKETM2
restart_oh
nonspinadapted

Some reference outputs for this example:

$ grep MPS dmrg-6.out
MPS =  KRRRRRRRRRRRRRRRRRRRRRRRRR 0 2 < N=8 SZ=-1 PG=0 >
GS INIT MPS BOND DIMS =       1    12    48   192   601  1145  1398  1474  1476  1468  1466  1441  1356  1316  1255  1240  1217  1206  1198  1176   904   422   183    59    16     4     1
$ grep Energy dmrg-6.out
OH Energy =  -75.423879916245909

We can see that the transformation from SU2 to SZ is nearly exact. The other two components of the SU2 MPS will also have the same energy as this one.

CSF or Determinant Sampling

The overlap between the spin-adapted MPS and Configuration State Functions (CSFs), or between the non-spin-adapted MPS and determinants can be calculated. Since there are exponentially many CSFs or determinants (when the number of electrons is close to the number of orbitals), normally it only makes sense to sample CSFs or determinants with (absolute value of) the overlap larger than a threshold. The sampling is deterministic, meaning that all overlap above the given threshold will be printed.

The keyword sample or restart_sample can be used to sample CSFs or determinants after DMRG or from an MPS loaded from disk. The value associated with the keyword sample or restart_sample is the threshold for sampling.

Setting the threshold to zero is allowed, but this may only be useful for some very small systems.

Limitations: For non-zero total spin CSF sampling, the spin-adapted MPS must be in singlet embedding format. See the previous section.

The following is an example of the input file:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 0
irrep 1

hf_occ integral
schedule default
maxM 500
maxiter 30

irrep_reorder
mps_tags KET
sample 0.05

Some reference outputs for this example:

$ grep CSF dmrg-1.out
Number of CSF =         17 (cutoff =      0.05)
Sum of weights of sampled CSF =    0.909360149891891
CSF          0 20000000000202000002000000  =    0.828657540546610
CSF          1 20200000000002000002000000  =   -0.330323898091116
CSF          2 20+00000000+0200000-000-00  =   -0.140063445607095
CSF          3 20+00000000+0-0-0002000000  =   -0.140041987646036
... ...
CSF         16 200000000002000+0-02000000  =    0.050020205617060

When there are more than 50 determinants, only the first 50 with largest weights will be printed. The complete list of determinants and coefficients are stored in sample-dets.npy and sample-vals.npy in the scratch folder, respectively.

So the restricted Hartree-Fock determinant/CSF has a very large coefficient (0.83).

To verify this, we can also directly compress the ground-state MPS to bond dimension 1, to get the CSF with the largest coefficient. Note that the compression method may converge to some other CSFs if there are many determinants with similar coefficients.

MPS Compression

MPS compression can be used to compress or fit a given MPS to a different (larger or smaller) bond dimension.

The following is an example of the input file for the compression (which will load the MPS obtailed from the previous ground-state DMRG):

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

nelec 8
spin 0
irrep 1

hf_occ integral
schedule
0  250  0 0
2  125  0 0
4   62  0 0
6   31  0 0
8   15  0 0
10   7  0 0
12   3  0 0
14   1  0 0
end
maxiter 16

compression
overlap
read_mps_tags KET
mps_tags BRA

irrep_reorder

Here the keyword compression indicates that this is a compression calculation. When the keyword overlap is given, the loaded MPS will be compressed, otherwise, the result of H|MPS> will be compressed. The tag of the input MPS is given by read_mps_tags, and the tag of the output MPS is given by mps_tags.

Some reference outputs for this example:

$ grep 'Compression overlap' dmrg-2.out
Compression overlap =    0.828657540546619

We can see that the value obtained from compression is very close to the sampled value. But when a lower bound of the overlap is known, the sampling method should be more reliable and efficient for obtaining the CSF with the largest weight.

If the CSF or determinat pattern is required, one can do a quick sampling on the compressed MPS using the keyword restart_sample 0.

If the given MPS has a very small bond dimension, or the target (output) MPS has a very large bond dimension (namely, “decompression”), one should use the keyword random_mps_init to allow a better random initial guess for the target MPS. Otherwise, the generated output MPS may be inaccurate.

LZ Symmetry

For diatomic molecules or model Hamiltonian with translational symmetry (such as 1D Hubbard model in momentum space), it is possible to utilize additional K space symmetry. To support the K space symmetry, the code must be compiled with the option -DUSE_KSYMM=ON (default).

One can add the keyword k_symmetry in the input file to use this additional symmetry. Point group symmetry can be used together with k symmetry. Therefore, even for system without K space symmetry, the calculation can still run as normal when the keyword k_symmetry is added. Note, however, the MPS or MPO generated from an input file with/without the keyword k_symmetry, cannot be reloaded with an input file without/with the keyword k_symmetry.

For molecules, the integral file (FCIDUMP file) must be generated in a special way so that the K/LZ symmetry can be used. the following python script can be used to generate the integral with \(C_2 \otimes L_z\) symmetry:

import numpy as np
from functools import reduce
from pyscf import gto, scf, ao2mo, symm, tools, lib
from block2 import FCIDUMP, VectorUInt8, VectorInt

# adapted from https://github.com/hczhai/pyscf/blob/1.6/examples/symm/33-lz_adaption.py
# with the sign of lz
def lz_symm_adaptation(mol):
    z_irrep_map = {} # map from dooh to lz
    g_irrep_map = {} # map from dooh to c2
    symm_orb_map = {} # orbital rotation
    for ix in mol.irrep_id:
        rx, qx = ix % 10, ix // 10
        g_irrep_map[ix] = rx & 4
        z_irrep_map[ix] = (-1) ** ((rx & 1) == ((rx & 4) >> 2)) * ((qx << 1) + ((rx & 2) >> 1))
        if z_irrep_map[ix] == 0:
            symm_orb_map[(ix, ix)] = 1
        else:
            if (rx & 1) == ((rx & 4) >> 2):
                symm_orb_map[(ix, ix)] = -np.sqrt(0.5) * ((rx & 2) - 1)
            else:
                symm_orb_map[(ix, ix)] = -np.sqrt(0.5) * 1j
            symm_orb_map[(ix, ix ^ 1)] = symm_orb_map[(ix, ix)] * 1j

    z_irrep_map = [z_irrep_map[ix] for ix in mol.irrep_id]
    g_irrep_map = [g_irrep_map[ix] for ix in mol.irrep_id]
    rev_symm_orb = [np.zeros_like(x) for x in mol.symm_orb]
    for iix, ix in enumerate(mol.irrep_id):
        for iiy, iy in enumerate(mol.irrep_id):
            if (ix, iy) in symm_orb_map:
                rev_symm_orb[iix] = rev_symm_orb[iix] + symm_orb_map[(ix, iy)] * mol.symm_orb[iiy]
    return rev_symm_orb, z_irrep_map, g_irrep_map

# copied from https://github.com/hczhai/pyscf/blob/1.6/pyscf/symm/addons.py#L29
# with the support for complex orbitals
def label_orb_symm(mol, irrep_name, symm_orb, mo, s=None, check=True, tol=1e-9):
    nmo = mo.shape[1]
    if s is None:
        s = mol.intor_symmetric('int1e_ovlp')
    s_mo = np.dot(s, mo)
    norm = np.zeros((len(irrep_name), nmo))
    for i, csym in enumerate(symm_orb):
        moso = np.dot(csym.conj().T, s_mo)
        ovlpso = reduce(np.dot, (csym.conj().T, s, csym))
        try:
            s_moso = lib.cho_solve(ovlpso, moso)
        except:
            ovlpso[np.diag_indices(csym.shape[1])] += 1e-12
            s_moso = lib.cho_solve(ovlpso, moso)
        norm[i] = np.einsum('ki,ki->i', moso.conj(), s_moso).real
    norm /= np.sum(norm, axis=0)  # for orbitals which are not normalized
    iridx = np.argmax(norm, axis=0)
    orbsym = np.asarray([irrep_name[i] for i in iridx])

    if check:
        largest_norm = norm[iridx,np.arange(nmo)]
        orbidx = np.where(largest_norm < 1-tol)[0]
        if orbidx.size > 0:
            idx = np.where(largest_norm < 1-tol*1e2)[0]
            if idx.size > 0:
                raise ValueError('orbitals %s not symmetrized, norm = %s' %
                                (idx, largest_norm[idx]))
            else:
                raise ValueError('orbitals %s not strictly symmetrized.',
                            np.unique(orbidx))
    return orbsym

mol = gto.M(
    atom=[["C", (0, 0, 0)],
          ["C", (0, 0, 1.2425)]],
    basis='ccpvdz',
    symmetry='dooh')

mol.symm_orb, z_irrep, g_irrep = lz_symm_adaptation(mol)
mf = scf.RHF(mol)
mf.run()

h1e = mf.mo_coeff.conj().T @ mf.get_hcore() @ mf.mo_coeff
print('h1e imag = ', np.linalg.norm(h1e.imag))
assert np.linalg.norm(h1e.imag) < 1E-14
e_core = mol.energy_nuc()
h1e = h1e.real.ravel()
_eri = ao2mo.restore(1, mf._eri, mol.nao)
g2e = np.einsum('pqrs,pi,qj,rk,sl->ijkl', _eri,
    mf.mo_coeff.conj(), mf.mo_coeff, mf.mo_coeff.conj(), mf.mo_coeff, optimize=True)
print('g2e imag = ', np.linalg.norm(g2e.imag))
assert np.linalg.norm(g2e.imag) < 1E-14
print('g2e symm = ', np.linalg.norm(g2e - g2e.transpose((1, 0, 3, 2))))
print('g2e symm = ', np.linalg.norm(g2e - g2e.transpose((2, 3, 0, 1))))
print('g2e symm = ', np.linalg.norm(g2e - g2e.transpose((3, 2, 1, 0))))
g2e = g2e.real.ravel()

fcidump_tol = 1E-13
na = nb = mol.nelectron // 2
n_mo = mol.nao
h1e[np.abs(h1e) < fcidump_tol] = 0
g2e[np.abs(g2e) < fcidump_tol] = 0

orb_sym_z = label_orb_symm(mol, z_irrep, mol.symm_orb, mf.mo_coeff, check=True)
orb_sym_g = label_orb_symm(mol, g_irrep, mol.symm_orb, mf.mo_coeff, check=True)
print(orb_sym_z)

fcidump = FCIDUMP()
fcidump.initialize_su2(n_mo, na + nb, na - nb, 1, e_core, h1e, g2e)

orb_sym_mp = VectorUInt8([tools.fcidump.ORBSYM_MAP['D2h'][i] for i in orb_sym_g])
fcidump.orb_sym = VectorUInt8(orb_sym_mp)
print('g symm error = ', fcidump.symmetrize(VectorUInt8(orb_sym_g)))

fcidump.k_sym = VectorInt(orb_sym_z)
fcidump.k_mod = 0
print('z symm error = ', fcidump.symmetrize(fcidump.k_sym, fcidump.k_mod))

fcidump.write('FCIDUMP')

Note that, if only the LZ symmetry is required, one can simply set orb_sym_g[:] = 0.

The following input file can be used to perform the calculation with \(C_2 \otimes L_z\) symmetry:

sym d2h
orbitals FCIDUMP
k_symmetry
k_irrep 0

nelec 12
spin 0
irrep 1

hf_occ integral
schedule
0  500 1E-8 1E-3
4  500 1E-8 1E-4
8  500 1E-9 1E-5
12 500 1E-9 0
end
maxiter 30

Where the k_irrep can be used to set the eigenvalue of LZ in the target state. Note that it can be easier for the Davidson procedure to get stuck in local minima with high symmetry. It is therefore recommended to use a custom schedule with larger noise and smaller Davidson threshold.

Some reference outputs for this input file:

$ grep 'Time elapsed' dmrg-1.out | tail -1
Time elapsed =     73.529 | E =     -75.7291544157 | DE = -6.31e-07 | DW = 1.28e-05
$ grep 'DMRG Energy' dmrg-1.out
DMRG Energy =  -75.729154415733063

When there are too many orbitals, and the default warmup fci initial guess is used, the initial MPS can have very large bond dimension (especially when the LZ symmetry is used, since LZ is not a finite group) and the first sweep will take very long time.

One way to solve this is to limit the LZ to a finite group, using modular arithmetic. We can limit LZ to Z4 or Z2. The efficiency gain will be smaller, but the convergence may be more stable. The keyword k_mod can be used to set the modulus. When k_mod = 0, it is the original infinite LZ group.

The following input file can be used to perform the calculation with \(C_2 \otimes Z_4\) symmetry:

sym d2h
orbitals FCIDUMP
k_symmetry
k_irrep 0
k_mod 4

nelec 12
spin 0
irrep 1

hf_occ integral
schedule
0  500 1E-8 1E-3
4  500 1E-8 1E-4
8  500 1E-9 1E-5
12 500 1E-9 0
end
maxiter 30

Some reference outputs for this input file:

$ grep 'Time elapsed' dmrg-2.out | tail -1
Time elapsed =    111.491 | E =     -75.7292222457 | DE = -8.17e-08 | DW = 1.28e-05
$ grep 'DMRG Energy' dmrg-2.out
DMRG Energy =  -75.729222245693876

Similarly, setting k_mod 2 gives the following output:

$ grep 'Time elapsed' dmrg-3.out | tail -1
Time elapsed =    135.394 | E =     -75.7314583188 | DE = -3.97e-07 | DW = 1.49e-05
$ grep 'DMRG Energy' dmrg-3.out
DMRG Energy =  -75.731458318751280

Initial Guess with Occupation Numbers

Once can use warmup occ initial guess to solve the initial guess problem, where another keywrod occ should be used, followed by a list of (fractional) occupation numbers separated by the space character, to set the occupation numbers. The occupation numbers can be obtained from a DMRG calculation using the same integral with/without K symmetry (or some other methods like CCSD and MP2). If onepdm is in the input file, the occupation numbers will be printed at the end of the output.

The following input file will perform the DMRG calculation using the same integral without the K symmetry (but with C2 symmetry):

sym d2h
orbitals FCIDUMP

nelec 12
spin 0
irrep 1

hf_occ integral
schedule
0  500 1E-8 1E-3
4  500 1E-8 1E-4
8  500 1E-9 1E-5
12 500 1E-9 0
end
maxiter 30
onepdm

Some reference outputs for this input file:

$ grep 'Time elapsed' dmrg-1.out | tail -2 | head -1
Time elapsed =    190.549 | E =     -75.7314655815 | DE = -1.88e-07 | DW = 1.53e-05
$ grep 'DMRG Energy' dmrg-1.out
DMRG Energy =  -75.731465581478815
$ grep 'DMRG OCC' dmrg-1.out
DMRG OCC =   2.000 2.000 1.957 1.626 1.870 1.870 0.360 0.098 0.098 0.006 0.008 0.008 0.008 0.013 0.014 0.014 0.011 0.006 0.006 0.006 0.005 0.005 0.002 0.002 0.002 0.001 0.001 0.001

The following input file will perform the DMRG calculation using the K symmetry, but with initial guess generated from occupation numbers:

sym d2h
orbitals FCIDUMP
k_symmetry
k_irrep 0
warmup occ
occ 2.000 2.000 1.957 1.626 1.870 1.870 0.360 0.098 0.098 0.006 0.008 0.008 0.008 0.013 0.014 0.014 0.011 0.006 0.006 0.006 0.005 0.005 0.002 0.002 0.002 0.001 0.001 0.001
cbias 0.2

nelec 12
spin 0
irrep 1

hf_occ integral
schedule
0  500 1E-8 1E-3
4  500 1E-8 1E-4
8  500 1E-9 1E-5
12 500 1E-9 0
end
maxiter 30

Here cbias is the keyword to add a constant bias to the occ, so that 2.0 becomes 2.0 - cbias, and 0.098 becomes 0.098 + cbias. Without the bias it is also easy to converge to a local minima.

Some reference outputs for this input file:

$ grep 'Time elapsed' dmrg-3.out | tail -1
Time elapsed =     55.938 | E =     -75.7244716369 | DE = -5.25e-07 | DW = 7.45e-06
$ grep 'DMRG Energy' dmrg-3.out
DMRG Energy =  -75.724471636942383

Here the calculation runs faster because the better initial guess, but the energy becomes worse.

Time Evolution

Now we give an example on how to do time evolution. The computation will apply \(|MPS_{out}\rangle = \exp (-t H) |MPS_{in}\rangle\) (with multiple steps). When \(t\) is a real floating point value, we will do imaginary time evolution of the MPS (namely, optimizing to ground state or finite-temperature state). When \(t\) is a pure imaginary value, we will do real time evolution of the MPS (namely, solving the time dependent Schrodinger equation).

To get accurate results, the time step has to be sufficiently small. The keyword delta_t is used to set a time step \(\Delta t\) and indicate that this is a time evolution calculation. The keyword target_t is used to set a target “stopping” time, namely, the \(t\). The “starting” time is considered as zero. Therefore, the number of time steps is computed as \(nsteps = t / \Delta t\) and printed.

If delta_t is too big, the time step error will be large. If delta_t is small, for fixed target time we have to do more time steps, with MPS bond dimension truncation happening after each sweep. So if delta_t is too small, the accumulated bond dimension truncation error will be large. Some meaningful time steps may be 0.01 to 0.1.

Real Time Evolution

First, we do a state-averaged calculation for the lowest two states using the following input file:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nroots 2

hf_occ integral
schedule default
maxM 500
maxiter 30

noreorder

Note that the orbital reordering is disabled. The output:

$ grep elapsed dmrg-1.out | tail -1
Time elapsed =      5.762 | E[  2] =     -75.7268133875    -75.6376794953 | DE = -8.89e-08 | DW = 6.38e-05
$ grep Final dmrg-1.out
Final canonical form =  LLLLLLLLLLLLLLLLLLLLLLLLLJ 25

The energy of the MPS at the last site is actually -75.72629673 and -75.63717415, which are slightly different from the above values.

Second, we can use the following input file to load the state-averaged MPS and then split it into individual MPSs:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG
nroots 2

hf_occ integral
schedule default
maxM 500
maxiter 30

restart_copy_mps
split_states
trans_mps_to_complex
noreorder

Note that here nroots must be the same as the previous case (or smaller, but larger than one), otherwise the state-averaged MPS cannot be correctly loaded. The state-averaged MPS has the default tag KET. We use calculation type keyword restart_copy_mps to do this transformation. The new keyword split_states indicates that we want to split the MPS, this keyword should only be used together with restart_copy_mps. The extra keyword trans_mps_to_complex will further make the MPS a complex MPS. This is required for real time evolution, where delta_t can be imaginary.

For imaginary time evolution and real delta_t and real target_t, everything will be real during the time evolution, so normally we do not need this extra keyword trans_mps_to_complex (but if you add it it is also okay).

The output looks like :

$ tail -7 dmrg-2.out
----- root =   0 /   2 -----
    final tag = KET-CPX-0
    final canonical form = LLLLLLLLLLLLLLLLLLLLLLLLLT
----- root =   1 /   2 -----
    final tag = KET-CPX-1
    final canonical form = LLLLLLLLLLLLLLLLLLLLLLLLLT
MPI FINALIZE: rank 0 of 1

By default, the tranformed MPS will have tags KET-0, KET-1 etc, if it is real, or KET-CPX-0, KET-CPX-1 etc if it is complex. If you set a custom tag, for example, when the input is like restart_copy_mps SKET, the tranformed MPS will have tags SKET-0, SKET-1, etc, no matter it is real or complex.

Third, we use the following script to do real time evolution:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

hf_occ integral
schedule
0 500 0 0
end
maxiter 10

read_mps_tags KET-CPX-0
mps_tags BRA
delta_t 0.05i
target_t 0.20i
complex_mps
noreorder

Note that a custom sweep schedule has to be used, to set the bond dimension to 500 (for example). The keyword maxiter and noise in the sweep schedule are ignored.

For every time step, there can be multiple sweeps, called “sub sweeps”. The total number of sweeps is n_sweeps = nsteps * n_sub_sweeps. The keyword n_sub_sweeps can be used to set the number of sub sweeps. Default value is 2.

For real time evolution, delta_t and target_t should be pure imaginary values. But they can also be general complex values. When doing imaginary time evolution, delta_t and target_t should be all real.

The tag of the input MPS (old MPS) is given by read_mps_tags. The tag of the output MPS (new MPS) is given by mps_tags. The two tags cannot be the same. They should (better) not have common prefix. For example, KET and KET-1 may not be used together, as -1 may be used by the code internally which will lead to confusion.

For this example, target_t is four times delta_t, so we will have 4 steps. Each time step has 2 sweeps. In total there will be 8 sweeps. The output is the result of applying \exp(-0.2i H) to the input.

Whenever a complex MPS is used, the keyword complex_mps should be used, otherwise the code will load the MPS incorrectly.

The output :

$ grep 'final' dmrg-3.out
    mps final tag = BRA
    mps final canonical form = MRRRRRRRRRRRRRRRRRRRRRRRRR
$ grep '<E>' dmrg-3.out
T = RE    0.00000 + IM    0.05000 <E> =  -75.726309692728165 <Norm^2> =    0.999999608946318
T = RE    0.00000 + IM    0.10000 <E> =  -75.726336818185246 <Norm^2> =    0.999994467614067
T = RE    0.00000 + IM    0.15000 <E> =  -75.726364807114123 <Norm^2> =    0.999990200387707
T = RE    0.00000 + IM    0.20000 <E> =  -75.726389514836484 <Norm^2> =    0.999986418355937

Here we see that the expectation value is printed after each time step. The energy is roughly conserved (similar to the DMRG output -75.72629673), and the norm is roughly one. Decreasing the time step may give more accurate results.

We can do the same for the excited state:

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

hf_occ integral
schedule
0 500 0 0
end
maxiter 10

read_mps_tags KET-CPX-1
mps_tags BRAEX
delta_t 0.05i
target_t 0.20i
complex_mps
noreorder

The output :

$ grep 'final' dmrg-4.out
    mps final tag = BRAEX
    mps final canonical form = MRRRRRRRRRRRRRRRRRRRRRRRRR
$ grep '<E>' dmrg-4.out
T = RE    0.00000 + IM    0.05000 <E> =  -75.637185795841717 <Norm^2> =    0.999999661398567
T = RE    0.00000 + IM    0.10000 <E> =  -75.637212093724074 <Norm^2> =    0.999995415040728
T = RE    0.00000 + IM    0.15000 <E> =  -75.637238086798163 <Norm^2> =    0.999991630799571
T = RE    0.00000 + IM    0.20000 <E> =  -75.637260508028248 <Norm^2> =    0.999988252849994

The energy is close to the DMRG value -75.63717415.

For imaginary time evolution, since the propagator is not unitary, the norm will increase exponentially. You may use the extra keyword normalize_mps to normalize MPS after each time step. The norm will still be computed and printed, but it will not be accumulated.

Finally, we can verify the energy at T = 0.0 and T = 0.2 and compute the overlap for these states. The overlap between the all four states can be computed using the following input :

sym d2h
orbitals C2.CAS.PVDZ.FCIDUMP.ORIG

hf_occ integral
schedule
0 500 0 0
end
maxiter 10

mps_tags KET-CPX-0 BRA KET-CPX-1 BRAEX
restart_tran_oh
complex_mps
overlap
noreorder

The output is:

$ grep 'OH' dmrg-5.out
OH Energy    0 -    0 = RE    1.000000000000002 + IM    0.000000000000000
OH Energy    1 -    0 = RE   -0.845792004408687 + IM   -0.533433527528264
OH Energy    1 -    1 = RE    0.999986418355938 + IM    0.000000000000000
OH Energy    2 -    0 = RE   -0.000000000000000 + IM    0.000000000000000
OH Energy    2 -    1 = RE   -0.000000827506956 + IM   -0.000000742303613
OH Energy    2 -    2 = RE    1.000000000000004 + IM    0.000000000000000
OH Energy    3 -    0 = RE    0.000001731091412 + IM   -0.000000316659748
OH Energy    3 -    1 = RE   -0.000001122421894 + IM    0.000002348984005
OH Energy    3 -    2 = RE   -0.836158473098047 + IM   -0.548435696470209
OH Energy    3 -    3 = RE    0.999988252849993 + IM    0.000000000000000

Here in the output each MPS gets a number, according to the order of tags in mps_tags. We have 0 (KET-CPX-0), 1 (BRA), 2 (KET-CPX-1) and 3 (BRAEX).

Note that state 1 (not normalized) is time evolved from state 0 (normalized). We see that the overlap <1|1> is exactly 1. To get the overlap between the normalized states, we have:

< normlized(0) | normlized(1) >
= <0|1> / sqrt(<0|0> * <1|1>)
= (-0.845792004408687 -0.533433527528264j) / sqrt( 0.999986418355938 * 1.000000000000002)
= -0.8457977480901698 -0.5334371500173138j

The absolute value and the angle of this complex overlap is :

   np.abs( -0.8457977480901698 -0.5334371500173138j ) =  0.9999645112167714
np.angle ( -0.8457977480901698 -0.5334371500173138j ) = -2.578911293480138

The absolute value is close to one. So the time evolution simply introduced a complex phase factor for the state, as expected. The complex phase factor can be computed as the remainder of E t divided by 2 pi:

-75.72638951483646 * 0.2 % (2 * np.pi) - 2 * np.pi = -2.5789072886081197

Which is close to the printed value.

Also note that the overlap between the ground state and the excited state <2|0> is exactly zero. The corresponding overlap between the time evolved states <3|1> is slightly different from zero, mainly due to the time step error and truncation error.

We can also get the energy expetation, by removing the keyword overlap:

$ grep 'OH' dmrg-6.out
OH Energy    0 -    0 = RE  -75.726296730204453 + IM    0.000000000000000
OH Energy    1 -    0 = RE   64.049088006450049 + IM   40.394772180607831
OH Energy    1 -    1 = RE  -75.725361025967970 + IM   -0.000000000000007
OH Energy    2 -    0 = RE    0.000000000000008 + IM    0.000000000000000
OH Energy    2 -    1 = RE    0.000061050951670 + IM    0.000056012958492
OH Energy    2 -    2 = RE  -75.637174152353893 + IM    0.000000000000000
OH Energy    3 -    0 = RE   -0.000132735557064 + IM    0.000024638559206
OH Energy    3 -    1 = RE    0.000086585167013 + IM   -0.000178008928209
OH Energy    3 -    2 = RE   63.244928578558032 + IM   41.482021915322555
OH Energy    3 -    3 = RE  -75.636371985782972 + IM    0.000000000000000

Note that here not all states are normalized, the printed value is not directly the energy. The printed value is <A|H|B>, but the energy is <A|H|B>/<A|B>. So the printed value should be divided by the square of the norm of the MPS (see previous output). For example, for state 1 we have :

-75.725361025967970 / 0.999986418355938 = -75.72638951483646

Which is the same as the number <E> printed by the time evolution (-75.726389514836484).