# DMRG for electronic structure calculations¶

Block program supports two executing modes: running standalone through command line or as a plugin of other quantum chemistry package. The Python-based quantum chemistry program package PySCF provides a simple solution to run Block program. It is the recommended way to use Block program in most scenario.

In PySCF, DMRG program is mainly used as a replacement of Full CI solver for large active space CASCI or CASSCF problem. On top of DMRG-CASCI and DMRG-CASSCF, MPS-PT can be called through Block-PySCF interface. Using Block with PySCF, systems around 50-active-orbital DMRG-CASSCF or 30-active-orbital MPSPT can be studied in a regular basis.

## CASCI/CASSCF in PySCF¶

PySCF is a collection Python modules for electronic structure simulation and theory developing. In this section, we briefly review the usage of PySCF. More usage details of PySCF package can be found in PySCF online documents http://www.pyscf.org. If you have PySCF installed and setup correctly (see http://www.pyscf.org/install.html), you can create a Python script for CASCI and CASSCF calculation:

$cat c2_cas.py from pyscf import gto, scf, mcscf mol = gto.M(atom="N 0 0 0; N 0 0 1", basis="ccpvdz") mf = scf.RHF(mol).run() ncas = 6 nelec_cas = 6 mc = mcscf.CASCI(mf, ncas, nelec_cas) mc.kernel() mc = mcscf.CASSCF(mf, ncas, nelec_cas) mc.kernel()  The Python script c2_cas.py can be executed by Python interpreter in command line: $ python c2_cas.py
converged SCF energy = -108.929838385609
CASCI E = -108.980200822148  E(CI) = -11.9360559617961  S^2 = 0.0000000
CASSCF energy = -109.044401900068


Based on the given active space size and number of correlated electrons, the CASCI/CASSCF solver by default takes the highest occupied and lowest unoccupied orbitals to form the active space. To change the active space, you need prepare a set of orbitals and reorder the orbitals to place the required active orbitals in the HOMO and LUMO space. You can feed the reordered orbitals to function mc.kernel(orbs) as the initial guess. CASCI/CASSCF solver will take the “HOMO/LUMO” orbitals from orbs as the active space. It is inconvenient to prepare the active space through this selecting-then-reordering procedure. To simplify this procedure, PySCF package provides some helper functions, such as sort_mo(), sort_mo_by_irrep(), dmet_cas.guess_cas(), atomic_valence() [1]. In the following example, we selected 4 $$\pi$$ orbitals and 1 $$\sigma$$ orbital and 1 $$\sigma^*$$ orbital from mean-field molecular orbitals to form the active space using the helper function sort_mo_by_irrep().

 1 2 3 4 5 6 7 8  from pyscf import gto, scf, mcscf mol = gto.M(atom="N 0 0 0; N 0 0 1", basis="ccpvdz", symmetry=True) mf = scf.RHF(mol).run() mc = mcscf.CASSCF(mf, 6, 6) ncas_by_irreps = {'E1ux':2, 'E1uy':2, 'A1g':1, 'A1u':1} orbs = mc.sort_mo_by_irrep(mf.mo_coeff, ncas_by_irreps) mc.kernel(orbs) 

In the above example, you should read the input as a Python program. Line 2 creates a molecule and applied mean-field calculation on the molecule. The mean-field results are saved in mf object so that you can access them later. For example, in line 5, HF MOs mf.mo_coeff are passed to function mc.sort_mo_by_irrep(). mc.sort_mo_by_irrep() read the configs from ncas_by_irreps and return the reordered orbitals orbs which is then fed to mc.kernel() function as the initial guess. mc is the CASSCF object created by mcscf.CASSCF function. More options can be specified for mc object to control the calculation. For example, you can set the convergence tolerance mc.conv_tol = 1e-6; require more computation details to be printed in the output with mc.verbose=5; call mc.analyze() to print out the population analysis of the CASCI/CASSCF results. The FCI solver of CASCI/CASSCF object is handled by the attribute mc.fcisolver. You can control the number of roots to compute by setting mc.fcisolver.nroots = 3, or change the symmetry of the correlated wave function with mc.fcisolver.wfnsym = 'A1u':

from pyscf import gto, scf, mcscf
mol = gto.M(atom="N 0 0 0; N 0 0 1", basis="ccpvdz", symmetry=True)
mf = scf.RHF(mol).run()

mc = mcscf.CASCI(mf, 6, 6)
ncas_by_irreps = {'E1ux':2, 'E1uy':2, 'A1g':1, 'A1u':1}
orbs = mcscf.sort_mo_by_irrep(mc, mf.mo_coeff, ncas_by_irreps)
mc.fcisolver.nroots = 3
mc.fcisolver.wfnsym = 'a1u'
mc.kernel(orbs)
mc.analyze()


Replacing mc.fcisolver with DMRG solver leads to the DMRG-CASCI and DMRG-CASSCF methods. But the rest of the input code should be the same to the regular CASCI/CASSCF calculation. You need the molecule and the mean-field objects to create the DMRG-CASCI/DMRG-CASSCF object mc. You can adjust the parameters in mc object to control the DMRG-CASCI/DMRG-CASSCF calculation and adjust DMRG configs through the mc.fcisolver object. More CASCI/CASSCF parameters are documented in http://www.pyscf.org/mcscf.html

## Setup Block in PySCF package¶

First you need prepare the Block executable binary. You can either compile it from source code [2] or download the precompiled binary block.spin_adapted-1.5.3.gz (compiled with Boost-1.55, OpenMPI-1.10.3, MKL-11) and the MPI-disabled version block.spin_adapted-1.5.3-serial.gz (compiled with Boost-1.55, MKL-11). Next, you need setup the Block runtime environment in PySCF. In the config file /path/to/pyscf/future/dmrgscf/settings.py (see also the template /path/to/pyscf/future/dmrgscf/settings.py.template), you need specify:

BLOCKEXE = "/path/to/Block/block.spin_adapted"
BLOCKSCRATCHDIR = "/path/to/scratch"
MPIPREFIX = "mpirun"  # or srun for SLURM system


You need at least set BLOCKEXE for DMRG-CASCI and DMRG-CASSCF methods. BLOCKSCRATCHDIR is the directory where to store temporary data and the DMRG wave function.

Note

Usually, the size of DMRG wave function is very large. Be sure that the disk which BLOCKSCRATCHDIR pointed to has enough space.

In the input script, you can replace the mc.fcisolver by DMRGCI object to call Block program in CASCI/CASSCF calculation:

from pyscf import gto, scf, mcscf
from pyscf import dmrgscf
mol = gto.M(atom="N 0 0 0; N 0 0 1", basis="ccpvdz", symmetry=True)
mf = scf.RHF(mol).run()

mc = mcscf.CASCI(mf, 6, 6)
ncas_by_irreps = {'E1ux':2, 'E1uy':2, 'A1g':1, 'A1u':1}
orbs = mcscf.sort_mo_by_irrep(mc, mf.mo_coeff, ncas_by_irreps)
mc.fcisolver = dmrgscf.DMRGCI(mol)
mc.fcisolver.nroots = 3
mc.fcisolver.wfnsym = 'a1u'
mc.kernel(orbs)
mc.analyze()


Generally speaking, this simple replacement of mc.fcisolver is enough to call the DMRG-CASCI and DMRG-CASSCF methods in your calculation. The rest settings of the mc object are all the same to the regular CASCI/CASSCF. When mc.kernel() is finished, the CASCI/CASSCF results such as orbital coefficients, natural occupancy etc. are held in mc object. But the DMRG wave-function is not. It is stored in the directory specified by the attribute DMRGCI.scratchDirectory or BLOCKSCRATCHDIR (the default value) in the config pyscf/future/dmrgscf/settings.py.

To make the embedded DMRG solver work more efficiently in CASSCF optimization, one needs carefully tune the DMRG parameters and dynamically update the parameters during the CASSCF optimization. It requires more codes in the interface to let CASSCF and DMRG talk to each other. We provided a shortcut function DMRGSCF() in the dmrgscf module to handle this functionality:

from pyscf import gto, scf, mcscf
from pyscf import dmrgscf
mol = gto.M(atom="N 0 0 0; N 0 0 1", basis="ccpvdz", symmetry=True)
mf = scf.RHF(mol).run()

mc = dmrgscf.DMRGSCF(mf, 6, 6)
ncas_by_irreps = {'E1ux':2, 'E1uy':2, 'A1g':1, 'A1u':1}
orbs = mcscf.sort_mo_by_irrep(mc, mf.mo_coeff, ncas_by_irreps)
mc.fcisolver.wfnsym = 'a1u'
mc.kernel(orbs)


We recommend to use dmrgscf.DMRGSCF() as the entry of DMRG-CASSCF method whenever is possible.

## Control Block program through PySCF wrapper¶

### Parallelism¶

Block-1.1.1 or older version support MPI level parallelization. The MPI parallelization parameters are controlled by the variable MPIPREFIX in pyscf/future/dmrgscf/settings.py or the attribute mpiprefix of DMRGCI object. For example, if you want to run Block using 4 processors on 2 nodes with Infiniband as the communication layer, you can specify in the input script:

mc = dmrgscf.DMRGSCF(mf, 6, 6)
mc.fcisolver.mpiprefix = 'mpirun -np 4 -npernode --mca btl self,openib'
mc.kernel()


If you are using SLURM system for job manager, you can put MPIPREFIX = 'srun' in the settings.py

To efficiently use memory, starting from Block-1.5, Block code introduces threading level parallelism, more specifically, the OpenMP threading. To enable the multi-threading feature in Block, you need specify the attribute num_thrds in DMRGCI object to indicate the maximum number of threads to be used by each MPI process:

mc = dmrgscf.DMRGSCF(mf, 6, 6)
mc.fcisolver.num_thrds = 4
mc.kernel()


By default, Block code uses 1 thread in each process. Using the multi-threading with the multi-processing model (mpirun -np) potentially offers higher performance and better scaling for DMRG parallelism. It is recommended to enable the multi-threading feature if your block program is newer than version 1.5.

On SLURM job system, the hybrid parallelism settings are controlled by SLURM runtime environment variables. You can control the parallel model by either configuring the resources through the #SBATCH flags or setting the \$SLURM_XXX variables in the SLURM script. For example, the following slurm script allocated in total 32 CPUs which are distributed in 8 processes on 2 nodes:

#SBATCH --nodes=2

python c2_cas.py


Specifying mc.fcisolver.mpiprefix = 'srun' will use SLURM to lanuch the Block program which will be executed on 2 nodes with 4 processes on each node. Note Block program does not detect the environment and setup the multi-threading automatically. You still need explicitly set mc.fcisolver.num_thrds = 4 in the PySCF input script to turn on the multi-threading for Block program.

### Bond dimension and sweep scheduler¶

Depending on the system, you may need change the DMRG bond dimension to improve the accuracy or balance the accuracy and efficiency. The default bond dimension is 1000. You can change the bond dimension by setting fcisolver.maxM:

from pyscf import gto, scf, mcscf
from pyscf import dmrgscf
mol = gto.M(atom="N 0 0 0; N 0 0 1", basis="ccpvdz", symmetry=True)
mf = scf.RHF(mol).run()

mc = dmrgscf.DMRGSCF(mf, 6, 6)
mc.fcisolver.maxM = 50
mc.kernel()


Generally, other default scheduler implemented in the PySCF wrapper should work fine in most systems. You can adjust the sweep schedule through the DMRGCI object:

dmrgsolver.scheduleSweeps = [0, 4, 8, 12, 16, 20, 24, 30]
dmrgsolver.scheduleMaxMs  = [200, 400, 800, 1200, 2000, 2000, 2000, 2000]
dmrgsolver.scheduleTols   = [0.0001, 0.0001, 0.0001, 0.0001, 1e-5, 1e-7, 1e-7, 1e-7]
dmrgsolver.scheduleNoises = [0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0, 0.0]
dmrgsolver.twodot_to_onedot = 34
dmrgsolver.maxIter = 50


The first four attributes which prefixed with schedule will be converted to the schedule section in the Block config file:

schedule
0  200   0.0001  0.0001
4  400   0.0001  0.0001
8  800   0.0001  0.0001
12 1200  0.0001  0.0001
16 2000  1e-5    0.0001
20 2000  1e-7    0.0001
24 2000  1e-7    0.0
30 2000  1e-7    0.0
end


In the early stage of Block sweep, the wave function is easy to stuck at local minimum. Although less efficient and accurate, applying the two-dot algorithm can effectively help DMRG solver moving out of the local minimum. Attribute twodot_to_onedot indicates when to switch to the one-dot algorithm which is efficient and stable to converge.

### DMRGCI functions and attributes¶

class DMRGCI

The interface of Block and PySCF. The class exposes the Block keywords to PySCF so that the Block code can be run and controlled in Python script.

approx_maxIter

In 1-step DMRG-CASSCF algorithm, the number of sweeps during the approximate FCI/DMRG updating. Default is 4

block_extra_keyword

It allows you to input Block keywords which were not exposed in DMRGCI class. Some commonly used keywords include

warmup local_2site
fiedler

See Keywords for the details of Block code keywords

configFile

By default, keywords are written to file dmrg.conf

dmrg_switch_tol

In 1-step DMRG-CASSCF, when the orbital gradients is smaller than this value, the DMRG calculation starts to read the solution from previous step as the initial guess (to reduce the computational cost). Default is 1e-3.

executable

Default is settings.BLOCKEXE

integralFile

The file to store FCIDUMP. Default is FCIDUMP

maxIter

Max number of sweeps

memory

The maximum memory (in GB) to use. Default is 2 GB. When you enabled multi threading and had large bond dimensions maxM, you might need more memory to hold the intermediates. Generally, large memory is helpful to improve efficiency.

New in version Block-1.5: (stackblock)

mpiprefix

Default is settings.MPIPREFIX

nroots

Number of states to solver simultaneously.

num_thrds

Number of OpenMP threads to be used in each MPI process. Default is 1.

outputFile

Block output. Default is dmrg.out

outputlevel

0 (less output) to 3 (very noise). Default is 2.

restart

Whether to read the wave function from the temporary directory (specified by scratchDirectory) as the initial guess.

Note

Block code does not check whether the system of the existed wave funciton matches the one in study. A mismatched DMRG wave function (from wrong DMRGCI.scratchDirectory) may lead to wrong solution or cause DMRG program crash.

runtimeDir

Where to put files dmrg.conf, dmrg.out etc temporarily. Default is current directory (where you execute python).

scratchDirectory

The directory where to store the intermediates and wave functions. Default is settings.BLOCKSCRATCHDIR.

Note

Be sure mc.fcisolver.scratchDirectory is properly assigned. Since all DMRGCI object by default uses the same BLOCKSCRATCHDIR settings, it’s easy to cause name conflicts on the scratch directory, especially when two DMRG-CASSCF calculations are run on the same node.

spin

2S (= nelec_alpha - nelec_beta). If the argument nelec of DMRGCI.kernel() function is a two-item list to represent the number of alpha and beta electrons, the Block program will use the given alpha and beta electron numbers to determine the spin. Otherwise, Block program takes this value as the spin of the system.

twodot_to_onedot

When to switch to one-dot algorithm.

weights

In state average calculation, the weight assocated to each state.

wfnsym

In the DMRGCI interface, the wave function symmetry ID follows the PySCF convention (see http://www.pyscf.org/symm.html). But Block code follows Molpro convention. A mapping between two symmetry ID is invoked in the DMRGCI initialization function. It is recommended to put the label of wave function (such as ‘A1g’, ‘B2u’) here to avoid the ambiguity.

make_rdm1(state, norb, nelec)

Given state ID, read its 1-particle density matrix from the directory indicated by scratchDirectory.

make_rdm12(state, norb, nelec)

Given state ID, read its 1-particle and 2-particle density matrices from the directory indicated by scratchDirectory. Note the 2-particle density matrix is reordered to match the 2e integrals of chemists’ notation, dm2[p,q,r,s] $$= \langle p^\dagger r^\dagger s q\rangle$$.

make_rdm123(state, norb, nelec)

Given state ID, read 1, 2 and 3-particle density matrices from the directory indicated by scratchDirectory. Note the 2-particle density matrix is reordered to match the 2e integrals of chemists’ notation. dm2[p,q,r,s] = $$\langle p^\dagger r^\dagger s q\rangle$$; The 3-particle density matrix takes the similar convention, dm3[p,q,r,s,t,u] $$= \langle p^\dagger r^\dagger t^\dagger u s q\rangle$$.

trans_rdm1(statebra, stateket, norb, nelec)

Given the state ID of bra and ket, read the 1-particle density matrix from the directory indicated by scratchDirectory.

trans_rdm12(statebra, stateket, norb, nelec)

Given the state ID of bra and ket, read the 1-particle and 2-particle density matrices from the directory indicated by scratchDirectory. Note the 2-particle density matrix is reordered to match the 2e integrals of chemists’ notation, dm2[p,q,r,s] $$= \langle p^\dagger r^\dagger s q\rangle$$.

kernel(h1e, eri, norb, nelec, fciRestart=None, ecore=0)

The kernel function to call Block program. “eri” is the array of 2-electron integrals (ij|kl). 8-fold permutation symmetry is required. The function returns the total energy and the state ID which is corresponding to the wave-function files in scratchDirectory. If multiple roots were required, the function returns two lists. The first list is the energy of each state. The second is a list of state ID.

DMRGSCF(mf, norb, nelec)

Shortcut function to setup CASSCF with the DMRG solver. The DMRG solver is properly initialized in this function so that the 1-step algorithm can applied efficiently in DMRG-CASSCF method.

Examples:

>>> mol = gto.M(atom='N 0 0 0; N 0 0 1')
>>> mf = scf.RHF(mol).run()
>>> mc = DMRGSCF(mf, 4, 4)
>>> mc.kernel()
-74.414908818611522


## State-average and state-specific DMRG-CASCI/DMRG-CASSCF¶

State-average and state-specific calculations were also supported in the DMRG-CASCI/DMRG-CASSCF through the Block-PySCF interface. The usage is the same to that in regular CASCI/CASSCF calculation. mc.state-average_() function provides the average over the multiple solutions over a single fcisolver:

mc = dmrgscf.DMRGSCF(mf, 6, 6)
# half-half average over ground state and first excited state
mc.state_average_([0.5, 0.5])
mc.kernel()


In this example, DMRGSCF() replaced the fcisolver with the DMRGCI object. Two DMRG states with the same spin and spatial (point group) symmetry are computed and half-half averaged. The two states are saved on the disk indicated by mc.fcisolver.scratchDirectory. In many calculations, one would require the state-average for states with different spin or spatial symmetry. Multiple FCI/DMRG solvers need to be created and each solver should handle one particular symmetry. Function mcscf.state_average_mix_() offers this functionality to mix different solvers in a single fcisolver object:

from pyscf import gto, scf, mcscf, dmrgscf
mol = gto.M(atom="N 0 0 0; N 0 0 1", basis="ccpvdz", symmetry=1, verbose=4)
mf = scf.RHF(mol).run()

mc = dmrgscf.DMRGSCF(mf, 6, 6)
weights = [.5, .25, .25]
solver1 = dmrgscf.DMRGCI(mol)
solver1.scratchDirectory = '/scratch/solver1'
solver1.nroots = 1
solver1.wfnsym = 'a1g'
solver1.spin = 2  # nelec_alpha - nelec_beta
solver2 = dmrgscf.DMRGCI(mol)
solver2.scratchDirectory = '/scratch/solver2'
solver2.nroots = 2
solver2.wfnsym = 'a1u'
mcscf.state_average_mix_(mc, [solver1, solver2], weights)
mc.kernel()


In this example, one solver for a triplet state of A1g symmetry and another solver for two singlet states of A1u symmetry are combined into one faked solver and assigned to fcisolver by state_average_mix_(). If the fake solver needs to handle solvers of different spin symmetry, you need explicitly assign the spin attribute to the solver. For first solver solver1, solver1.spin = 2 indciates that the number of alpha electrons is 2 more than the number of beta electrons. The kernel() function of fake solver mc.fcisolver will return 3 states in a list [0, 0, 1]. The number in the list represents the state ID in each solver. The first state (the first 0 in the list) is obtained from solver1. Its wave-function and density matrices can be found in /scratch/solver1. The second and third elements of [0, 0, 1] are the states obtained from solver2. The relevant wave functions and density matrices are all stored in /scratch/solver2.

Note

Block program stores the wave function in scratchDirectory. You must assign different scratchDirectory for different DMRG solvers. If two Block wave function are put in the same scratchDirectory, the solver may crash or produce wrong solution.

State-specific DMRG-CASSCF is the other common calculation one would take. Setting up state-specific DMRG-CASSCF object is the same to the regular CASSCF code. By calling mc.state_specific_() function with state ID: 0 for ground state, 1 for first excited state ..., you can optimize the target state with DMRG-CASSCF:

# Optimize the first excited state
mc = dmrgscf.DMRGSCF(mf, 6, 6)
mc.state_specific_(state=1)
mc.kernel()


The mc.state_specific_() function can be applied with DMRG-CASCI object as well. However, a straightforward solution for DMRG-CASCI is to compute multiple states simultaneously with attribute nroots:

mc = mcscf.CASCI(mf, 6, 6)
mc.fcisolver = dmrgscf.DMRGCI(mol)
mc.fcisolver.nroots = 5
mc.kernel()


In PySCF source code, you can find more examples of state-average and state-specific calculations.

## DMRG-NEVPT2¶

For Block 1.1.1 version or older, the standard DMRG-NEVPT2 calculation can be carried out on top of the DMRG-CASCI or DMRG-CASSCF calculation:

from pyscf import gto, scf, dmrgscf, mrpt
mol = gto.M(atom="N 0 0 0; N 0 0 1", basis="ccpvdz")
mf = scf.RHF(mol).run()

mc = dmrgscf.DMRGSCF(mf, 6, 6).run()
mrpt.NEVPT(mc).run()

mc = mcscf.CASCI(mf, 6, 6)
mc.fcisolver = dmrgscf.DMRGCI(mol)
mc.run()
mrpt.NEVPT(mc).run()


The standard DMRG-NEVPT2 method requires the 4-particle density matrix. Computing and storing the 4-particle density matrix is extremely demanding. It limits the system size to at most 26 orbitals. Starting from Block 1.1 version, we implemented an effective approximation based on compressed MPS-perturber technique which can significantly reduce the computation cost. The MPS-perturber NEVPT2 implementation requires the MPI4Py library and the serial version of Block program. You need set in the config file /path/to/pyscf/future/dmrgscf/settings.py the variable BLOCKEXE_COMPRESS_NEVPT:

BLOCKEXE_COMPRESS_NEVPT = "/path/to/serially/compiled/Block/block.spin_adapted-serial"


Note

The wavefunction structure from different Block versions are incompatible. If BLOCKEXE for zeroth order wavefunction is set to Block-1.1, BLOCKEXE_COMPRESS_NEVPT for PT should also be Block-1.1. Similarly, Block-1.5 (stackblock) PT code only compatible with the zeroth order wavefunction of Block-1.5 (stackblock).

Now you can use compress_approx() function to initialize a compressed pertuber NEVPT2 method. In the compress_approx() function, we precomputed the most demanding intermediates and stored them on disk:

from pyscf import gto, scf, dmrgscf, mrpt
mol = gto.M(atom="N 0 0 0; N 0 0 1", basis="ccpvdz")
mf = scf.RHF(mol).run()
mc = dmrgscf.dmrgci.DMRGSCF(mf, 6, 6).run()

mrpt.NEVPT(mc).compress_approx().run()


Note

The compressed NEVPT2 algorithm is also very demanding, especially on the memory usage. It can support up to about 35 orbitals in Block-1.5. Please refer to the Benchmark for approximate costs.

If the excitation energy is of interest, we can use DMRG-NEVPT2 to compute the energy of excited state based on the multiple-root CASCI calculations:

mc = mcscf.CASCI(mf, 6, 6)
mc.fcisolver = dmrgscf.DMRGCI(mol)
mc.fcisolver.nroots = 2
mc.kernel()
mrpt.NEVPT(mc, root=0).compress_approx(maxM=100).run()
mrpt.NEVPT(mc, root=1).compress_approx(maxM=100).run()


In the above example, two NEVPT2 calculations are called separately for two states which are indicated by the argument root=*. If the DMRG-NEVPT2 calculations are called based on the state-average DMRG-CASSCF calculation, you should be very careful with scratchDirectory for the DMRG wave function that NEVPT2 perturbation is applied on. In the multiple-solver state-average DMRG-CASSCF calculation, you need assign the right fcisolver and state ID to the mc object before passing it to mrpt.NEVPT() method.:

mc = dmrgscf.DMRGSCF(mf, 6, 6)
weights = [.5, .25, .25]
solver1 = dmrgscf.DMRGCI(mol)
solver1.scratchDirectory = '/scratch/solver1'
solver1.nroots = 1
solver1.wfnsym = 'a1g'
solver1.spin = 2  # nelec_alpha - nelec_beta
solver2 = dmrgscf.DMRGCI(mol)
solver2.scratchDirectory = '/scratch/solver2'
solver2.nroots = 2
solver2.wfnsym = 'a1u'
mcscf.state_average_mix_(mc, [solver1, solver2], weights)
mc.kernel()

mc.fcisolver = solver1
mrpt.NEVPT(mc, root=0).compress_approx(maxM=100).run()

mc.fcisolver = solver2
mrpt.NEVPT(mc, root=1).compress_approx(maxM=100).run()


## Case study¶

#!/usr/bin/env python
#
# Contributors:
#       Zhendong Li <zhendongli2008@gmail.com>
#       Qiming Sun <osirpt.sun@gmail.com>
#

from functools import reduce
import numpy
import scipy.linalg
from pyscf import gto, scf, mcscf
from pyscf import tools

mol = gto.Mole()
mol.verbose = 4
mol.output = 'hs.out'
mol.atom = '''
Mo  7.34411020207581      1.17495097908005      6.72284263905920
Fe  7.84036274632279      2.40948832380662      3.90857987198295
S   8.11413397508734      3.34683967317511      5.92473122721237
Cl  9.42237962692288      2.83901882053830      2.40523971787167
S   7.63129189448338      0.24683725427736      4.48256715460659
Cl  5.78420653505383      3.15381896731458      3.13969003482939
N   7.05276738605521      2.42445066370842      8.66076404425459
N   6.64167403727998     -1.00707407196440      7.11600799320614
N   5.24002742536262      1.70306900993116      5.97156233521481
H   6.27407522563538      2.37009344884271      9.32021452836747
H   7.93656914286549      2.22405698642280      9.14675757456406
H   7.12313637861828      3.37423478174186      8.26848891472229
H   4.53157107313027      2.02429015953190      6.63557341863725
H   5.36325579034589      2.43796505637839      5.24458486826946
H   4.86812692298093      0.90155604764634      5.45231738540969
H   5.82209673287966     -1.30027608533573      7.65147593202357
H   6.56861368828978     -1.33871396574670      6.14286445056596
H   7.48831993436433     -1.42577562013418      7.51985443522153
S   9.40594188780477      0.42545761808747      7.87277304102829
H   8.82966334944139     -0.10099345030206      8.99111747895267
'''
mol.basis = 'tzp-dkh'
mol.charge = -1
mol.spin = 8
mol.build()

#
# X2C correction for relativistic effects
#
# first pass, to generate initial guess
#
mf = scf.sfx2c(scf.UHF(mol))
mf.chkfile = 'hs.chk'
mf.level_shift = 0.1
mf.conv_tol = 1e-2
mf.kernel()
#
# second pass to converge SCF calculation
#
mf = scf.newton(mf)
mf.conv_tol = 1e-12
mf.kernel()

##################################################
#
# Analyze SCF results and make MCSCF initial guess
#
##################################################

# This parameter to control the call to matplotlib
ifplot = False

def sqrtm(s):
e, v = numpy.linalg.eigh(s)
return numpy.dot(v*numpy.sqrt(e), v.T.conj())

def lowdin(s):
e, v = numpy.linalg.eigh(s)
return numpy.dot(v/numpy.sqrt(e), v.T.conj())

##################################################
#
# 1. Read UHF-alpha/beta orbitals from chkfile
#
##################################################
fname = 'hs.chk'
chkfile = fname
mo_coeff = mf["mo_coeff"]
ma = mo_coeff[0]
mb = mo_coeff[1]
nb = ma.shape[1]
nalpha = (mol.nelectron+mol.spin)/2
nbeta  = (mol.nelectron-mol.spin)/2
print('Nalpha = %d, Nbeta %d, Sz = %d, Norb = %d' % (nalpha, nbeta, mol.spin, nb))

#=============================
# DUMP from chkfile to molden
#=============================
#
# One can view the orbitals in many visualization tool like Jmol, IBOviewer
#
moldenfile = fname+'0.molden'
tools.molden.from_chkfile(moldenfile, chkfile)
if 0:
# Jmol script to generate orbital images.  Run this jmol script in command line
#       jmol hs.spt
# It writes out images for 10 HOMO and 10 LUMO orbitals.
tools.molden.from_mo(mol, fname+'_alpha.molden', ma)
tools.molden.from_mo(mol, fname+'_beta.molden', mb)
jmol_script = 'hs.spt'
fspt = open(jmol_script,'w')
fspt.write('''
initialize;
set background [xffffff];
set frank off
set autoBond true;
set bondTolerance 0.5;
set forceAutoBond false;
#cd /home/abc/pyscf/examples/dmrg
''')
fspt.write('rotate -30 y;\n')
fspt.write('rotate 20 x;\n')
for i in range(nalpha-10,nalpha+10):
fspt.write('isoSurface MO %d fill noMesh noDots;\n' % (i+1))
fspt.write('#color isoSurface translucent 0.6 [x00ff00];\n')
fspt.write('write JPG 90 "%s-alpha-%d.jpg";\n' % (jmol_script, (i+1)))

fspt.write('rotate -30 y;\n')
fspt.write('rotate 20 x;\n')
for i in range(nbeta-10,nbeta+10):
fspt.write('isoSurface MO %d fill noMesh noDots;\n' % (i+1))
fspt.write('#color isoSurface translucent 0.6 [x0000ff];\n')
fspt.write('write JPG 90 "%s-beta-%d.jpg";\n' % (jmol_script, (i+1)))
fspt.close()

##################################################
#
# 2. Sanity check, using eg orthogonality
#
##################################################
ova = mol.intor_symmetric("cint1e_ovlp_sph")
diff = reduce(numpy.dot,(mo_coeff[0].T,ova,mo_coeff[0])) - numpy.identity(nb)
print(numpy.linalg.norm(diff))
diff = reduce(numpy.dot,(mo_coeff[1].T,ova,mo_coeff[1])) - numpy.identity(nb)
print(numpy.linalg.norm(diff))

#=============================
# Natural orbitals
# Lowdin basis X=S{-1/2}
# psi = chi * C
#     = chi' * C'
#     = chi*X*(X{-1}C')
#=============================
pTa = numpy.dot(ma[:,:nalpha],ma[:,:nalpha].T)
pTb = numpy.dot(mb[:,:nbeta],mb[:,:nbeta].T)
pT = pTa+pTb
pT = 0.5*pT
# Lowdin basis
s12 = sqrtm(ova)
s12inv = lowdin(ova)
pT = reduce(numpy.dot,(s12,pT,s12))
print('idemponency of DM: %s' % numpy.linalg.norm(pT.dot(pT)-pT))
enorb = mf["mo_energy"]
print('\nCMO_enorb:')
print(enorb)
if ifplot:
import matplotlib.pyplot as plt
plt.plot(range(nb),enorb[0],'ro')
plt.plot(range(nb),enorb[1],'bo')
plt.show()

#
# Non-orthogonal cases: FC=SCE
# Fao = SC*e*C{-1} = S*C*e*Ct*S
# OAO basis:
# F = Xt*Fao*X = S1/2*C*e*Ct*S1/2
#
fa = reduce(numpy.dot,(ma,numpy.diag(enorb[0]),ma.T))
fb = reduce(numpy.dot,(mb,numpy.diag(enorb[1]),mb.T))
fav = (fa+fb)/2
fock_sf = fOAO = reduce(numpy.dot,(s12,fav,s12))
#
# Small level shift on density matrix to break occupation degeneracy in natual
# orbitals
#
shift = 1e-7
pTshift = pT + shift*fOAO
#
# 'natural' occupations and orbitals
#
eig,coeff = scipy.linalg.eigh(pTshift)
eig = 2*eig
print('Natual occupancy %s ' % eig)
eig[abs(eig)<1.e-14]=0.0
if ifplot:
import matplotlib.pyplot as plt
plt.plot(range(nb),eig,'ro')
plt.show()
#
# Rotate back to AO representation and check orthogonality
#
coeff = numpy.dot(s12inv,coeff)
ova = mol.intor_symmetric("cint1e_ovlp_sph")
diff = reduce(numpy.dot,(coeff.T,ova,coeff)) - numpy.identity(nb)
print('CtSC-I',numpy.linalg.norm(diff))

##################################################
#
# 3. Search for active space
#
##################################################

#
# 3.1 Transform the entire MO space into core, active, and external space
# based on natural occupancy
#
# Expectation value of natural orbitals <i|F|i>
fexpt = reduce(numpy.dot,(coeff.T,ova,fav,ova,coeff))
enorb = numpy.diag(fexpt)
# Sort by occupancy
index = numpy.argsort(-eig)
enorb = enorb[index]
nocc  = eig[index]
coeff = coeff[:,index]
#
# Reordering and define active space according to thresh
#
thresh = 0.01
active = (thresh <= nocc) & (nocc <= 2-thresh)
print('\nNatural orbitals:')
print('Offdiag(F) = %s' % numpy.linalg.norm(fexpt - numpy.diag(enorb)))
for i in range(nb):
print('orb:',i,active[i],nocc[i],enorb[i])
actIndices = numpy.where(active)[0]
print('active orbital indices %s' % actIndices)
print('Num active orbitals %d' % len(actIndices))
cOrbs = coeff[:,:actIndices[0]]
aOrbs = coeff[:,actIndices]
vOrbs = coeff[:,actIndices[-1]+1:]
nb = cOrbs.shape[0]
nc = cOrbs.shape[1]
na = aOrbs.shape[1]
nv = vOrbs.shape[1]
print('core orbs:',cOrbs.shape)
print('act  orbs:',aOrbs.shape)
print('vir  orbs:',vOrbs.shape)
assert nc+na+nv == nb

#
# 3.2 Localizing core, active, external space separately, based on certain
# local orbitals.
#
# We now dump out UHF natual orbitals and localized orbitals to help identify
# active space.
#
# dump UHF natrual orbital
#
tools.molden.from_mo(mol, fname+'_uno.molden', coeff)

#=============================
# localized orbitals
#=============================
iflocal  = False
if iflocal:
# We implemented different localization later
from pyscf.tools import localizer
loc = localizer.localizer(mol,ma[:,:mol.nelectron/2],'boys')
loc.verbose = 10
new_coeff = loc.optimize()
loc = localizer.localizer(mol,ma[:,mol.nelectron/2:],'boys')
new_coeff2 = loc.optimize()
lmo = numpy.hstack([new_coeff,new_coeff2])
tools.molden.from_mo(mol, fname+'lmo.molden', lmo)

#
# Test orthogonality because occasionally localization procedure may break the
# orbital orthogonality (when AO functions are close to linear dependent).
#
cOrbsOAO = numpy.dot(s12,cOrbs)
aOrbsOAO = numpy.dot(s12,aOrbs)
vOrbsOAO = numpy.dot(s12,vOrbs)
print('Ortho-cOAO',numpy.linalg.norm(numpy.dot(cOrbsOAO.T,cOrbsOAO)-numpy.identity(nc)))
print('Ortho-aOAO',numpy.linalg.norm(numpy.dot(aOrbsOAO.T,aOrbsOAO)-numpy.identity(na)))
print('Ortho-vOAO',numpy.linalg.norm(numpy.dot(vOrbsOAO.T,vOrbsOAO)-numpy.identity(nv)))

#==========================================
# Now try to get localized molecular orbitals (SCDM)
#==========================================
def scdm(coeff, overlap, aux):
#
# Argument coeff is a set of orthogonal orbitals |O> (eg occupied HF
# orbitals); aux is a set of localized orbitals.  One can define a subset |B>
# of aux, which has the closest overlap to the coeff space.
# The (orthogonalized) resultant local orbitals |B> can be considered as the
# localized coeff |O>
#
#       |B> = |O><O|aux>, in which det(<O|aux>) is maximized;
#       return lowdin(|B>)
#
no = coeff.shape[1]
ova = reduce(numpy.dot,(coeff.T, overlap, aux))
# ova = no*nb
q,r,piv = scipy.linalg.qr(ova, pivoting=True)
# piv[:no] defines the subset of aux which has the largest overlap to coeff space
bc = ova[:,piv[:no]]

ova = numpy.dot(bc.T,bc)
s12inv = lowdin(ova)
cnew = reduce(numpy.dot,(coeff,bc,s12inv))
return cnew
#
# Various choices for the localized orbitals
# * the non-orthogonal AOs
#       aux=numpy.identity(nb)
# * Lowdin orthogonalized AOs
aux = s12inv
# * Meta-lowdin orthogonalized AOs
#       from pyscf import lo
#       aux = lo.orth.orth_ao(mol,method='meta_lowdin',pre_orth_ao=lo.orth.pre_orth_ao(mol))
# * ...
#
ova = mol.intor_symmetric("cint1e_ovlp_sph")
clmo = scdm(cOrbs, ova, aux)  # local "AOs" in core space
almo = scdm(aOrbs, ova, aux)  # local "AOs" in active space
vlmo = scdm(vOrbs, ova, aux)  # local "AOs" in external space

#
# 3.3 Sorting each space (core, active, external) based on "orbital energy" to
# prevent high-lying orbitals standing in valence space.
#
# Get <i|F|i>
def psort(ova, fav, coeff):
# pT is density matrix, fav is Fock matrix
# OCC-SORT
pTnew = 2.0*reduce(numpy.dot,(coeff.T,s12,pT,s12,coeff))
nocc  = numpy.diag(pTnew)
index = numpy.argsort(-nocc)
ncoeff = coeff[:,index]
nocc   = nocc[index]
enorb = numpy.diag(reduce(numpy.dot,(ncoeff.T,ova,fav,ova,ncoeff)))
return ncoeff, nocc, enorb

# E-SORT
mo_c, n_c, e_c = psort(ova, fav, clmo)
mo_o, n_o, e_o = psort(ova, fav, almo)
mo_v, n_v, e_v = psort(ova, fav, vlmo)
#
# coeff is the local molecular orbitals
#
coeff = numpy.hstack((mo_c, mo_o, mo_v))

#
# Test orthogonality for the localize MOs as before
#
diff = reduce(numpy.dot,(coeff.T,ova,coeff)) - numpy.identity(nb)
print('diff=',numpy.linalg.norm(diff))
tools.molden.from_mo(mol, fname+'_scdm.molden', coeff)

#
# Population analysis to confirm that our LMO (coeff) make sense
#
#==========================================
# lowdin-pop of the obtained LMOs in OAOs
#==========================================
lcoeff = s12.dot(coeff)
# Orthogonality test
diff = reduce(numpy.dot,(lcoeff.T,lcoeff)) - numpy.identity(nb)
print('diff=',numpy.linalg.norm(diff))

print('\nLowdin population for LMOs:')

pthresh = 0.02
labels = mol.ao_labels(None)
ifACTONLY = False #True
nelec = 0.0
nact = 0.0
for iorb in range(nb):
vec = lcoeff[:,iorb]**2
idx = list(numpy.argwhere(vec>pthresh))
if ifACTONLY == False:
if iorb < nc:
print(' iorb_C=',iorb,' occ=',n_c[iorb],' fii=',e_c[iorb])
nelec += n_c[iorb]
elif iorb >= nc and iorb < nc+na:
print(' iorb_A=',iorb,' occ=',n_o[iorb-nc],' faa=',e_o[iorb-nc])
nelec += n_o[iorb-nc]
else:
print(' iorb_V=',iorb,' occ=',n_v[iorb-nc-na],' fvv=',e_v[iorb-nc-na])
nelec += n_v[iorb-nc-na]
for iao in idx:
print('    iao=',labels[iao],' pop=',vec[iao])
else:
if iorb >= nc and iorb < nc+na:
print(' iorb_A=',iorb,' faa=',e_o[iorb-nc])
for iao in idx:
print('    iao=',labels[iao],' pop=',vec[iao])
print('nelec=',nelec)

#
# 3.4 select 'active' orbitals
#
# By reading the orbital images with Jmol, we characterized some of the local
# orbitals
#
a1 = [80,82,83,84,85,86] # S-3p = 7
o1 = [ 2]*6  # approximate occupancies, to help obtain the electrons in active space
a2 = [87,88,89,90,91,92,93,94,95,96] # Fe-Mo (3d,4d) = 10
o2 = [ 1]*8+[0]*2  # approximate occupancies
a3 = [97,98,99,101,103,105] # Mo-s + Fe4d = 6
o3 = [0]*6  # approximate occupancies

#
# There are many different choices for active space, here we just demonstrate
# one which is consists of Fe 3d, Mo 4d and S 3p orbitals
#
#==========================
# select 'active' orbitals
#==========================
caslst = a1+a2
norb = len(caslst)
ne_act = sum(o1) + sum(o2)
s = 1 # 0,1,2,3,4, High-spin case ms = s
ne_alpha = ne_act/2 + s
ne_beta  = ne_act/2 - s
nalpha = ne_alpha
nbeta = ne_beta
norb = len(caslst)
print('norb/nacte=',norb,[nalpha,nbeta])

##################################################
#
# 4. DMRG-CASSCF and DMRG-NEVPT2
#
##################################################
#
# Adjust the MPI schedular and scratch directory if needed.
# NOTE the DMRG-NEVPT2 is expensive, it requires about 8 GB memory per processor
#
#from pyscf.dmrgscf import settings
#settings.MPIPREFIX = 'srun'
#settings.BLOCKSCRATCHDIR = '/scratch'

from pyscf.dmrgscf import DMRGCI, DMRGSCF
from pyscf import mrpt

#
# Redirect output to another file
#
mol.build(verbose=7, output = 'hs_dmrg.out')

mf = scf.sfx2c1e(scf.RHF(mol))
mc = DMRGSCF(mf, norb, [nalpha,nbeta])
mc.chkfile = 'hs_mc.chk'
mc.max_memory = 30000
mc.fcisolver.maxM = 1000
mc.fcisolver.tol = 1e-6
orbs = mc.sort_mo(caslst, coeff, base=0)
mc.mc1step(orbs)

#
# CASCI-NEVPT2
#
# If DMRG-CASSCF was finished without any problems (eg convergence, wall time
# limits on cluster etc),  one can simply continue with DMRG-NEVPT2
#       mrpt.NEVPT(mc).kernel()
#
# But it's highly possible that one needs to restore the calculation from
# previous work.  The following is an example to restore the calculation.
# Assuming DMRG-CASSCF has converged and the DMRG temporary files were
# deleted, we just need the DMRG-CASCI calculation with the converged MCSCF
# orbitals to get the DMRG wavefunction.
#
mc = mcscf.CASCI(mf, norb, [nalpha,nbeta])
mc.chkfile = 'hs_mc.chk'
mc.fcisolver = DMRGCI(mol,maxM=500, tol =1e-8)
#
# Tune DMRG parameters.  It's not necessary in most scenario.
#
#mc.fcisolver.outputlevel = 3
#mc.fcisolver.scheduleSweeps = [0, 4, 8, 12, 16, 20, 24, 28, 30, 34]
#mc.fcisolver.scheduleMaxMs  = [200, 400, 800, 1200, 2000, 4000, 3000, 2000, 1000, 500]
#mc.fcisolver.scheduleTols   = [0.0001, 0.0001, 0.0001, 0.0001, 1e-5, 1e-6, 1e-7, 1e-7, 1e-7, 1e-7 ]
#mc.fcisolver.scheduleNoises = [0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0001, 0.0, 0.0, 0.0, 0.0]
#mc.fcisolver.twodot_to_onedot = 38
#mc.fcisolver.maxIter = 50
mc.casci(mo)

#
# DMRG-NEVPT2
#
mrpt.NEVPT(mc).kernel()

#
# There is also a fast DMRG-NEVPT2 implementation.  See also the example
# pyscf/examples/dmrg/02-dmrg_nevpt2.py
#
mrpt.NEVPT(mc).compress_approx().kernel()

##################################################
#
# Don't forget to clean up the scratch.  DMRG calculation can produce large
# amount of temporary files.
#
##################################################


## Run Block standalone¶

Block program can be run standalone without the PySCF environments. In PySCF-1.3, the DMRG interface provides dry run mode to generate the Block input config dmrg.conf and the integral file FCIDUMP.:

from pyscf import gto, scf, dmrgscf
mf = gto.M(atom="N 0 0 0; N 0 0 1", basis="ccpvdz").apply(scf.RHF).run()
mc = dmrgscf.DMRGSCF(mf, 6, 6)
dmrgscf.dryrun(mc)


You can execute Block program in command line:

mpirun -n 2 block.spin_adapted dmrg.conf > dmrg.out


See more examples in Chapter Block program as a standalone solver.

Footnotes

 [1] cite paper.