.. highlight:: bash .. _user_advanced: 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 hf_occ integral 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``. .. highlight:: python3 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 :math:`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') .. highlight:: text 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 :math:`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 :math:`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 :math:`|MPS_{out}\rangle = \exp (-t H) |MPS_{in}\rangle` (with multiple steps). When :math:`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 :math:`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 :math:`\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 :math:`t`. The "starting" time is considered as zero. Therefore, the number of time steps is computed as :math:`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 '' dmrg-3.out T = RE 0.00000 + IM 0.05000 = -75.726309692728165 = 0.999999608946318 T = RE 0.00000 + IM 0.10000 = -75.726336818185246 = 0.999994467614067 T = RE 0.00000 + IM 0.15000 = -75.726364807114123 = 0.999990200387707 T = RE 0.00000 + IM 0.20000 = -75.726389514836484 = 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 '' dmrg-4.out T = RE 0.00000 + IM 0.05000 = -75.637185795841717 = 0.999999661398567 T = RE 0.00000 + IM 0.10000 = -75.637212093724074 = 0.999995415040728 T = RE 0.00000 + IM 0.15000 = -75.637238086798163 = 0.999991630799571 T = RE 0.00000 + IM 0.20000 = -75.637260508028248 = 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 ````, but the energy is ``/``. 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 ```` printed by the time evolution (-75.726389514836484).