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 Pythonbased 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 DMRGCASCI and DMRGCASSCF, MPSPT can be called through BlockPySCF interface. Using Block with PySCF, systems around 50activeorbital DMRGCASSCF or 30activeorbital 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 selectingthenreordering 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 meanfield 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 meanfield calculation on the molecule. The
meanfield 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 = 1e6
; 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 DMRGCASCI and
DMRGCASSCF methods. But the rest of the input code should be the same to the
regular CASCI/CASSCF calculation. You need the molecule and the meanfield
objects to create the DMRGCASCI/DMRGCASSCF object mc
. You can adjust the
parameters in mc
object to control the DMRGCASCI/DMRGCASSCF 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_adapted1.5.3.gz
(compiled with Boost1.55, OpenMPI1.10.3, MKL11) and the MPIdisabled version
block.spin_adapted1.5.3serial.gz
(compiled with Boost1.55, MKL11).
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"
BLOCKEXE_COMPRESS_NEVPT = "/path/to/serially/compiled/Block/block.spin_adapted"
BLOCKSCRATCHDIR = "/path/to/scratch"
MPIPREFIX = "mpirun" # or srun for SLURM system
You need at least set BLOCKEXE
for DMRGCASCI and DMRGCASSCF 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 DMRGCASCI and DMRGCASSCF 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
wavefunction 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 DMRGCASSCF method
whenever is possible.
Control Block program through PySCF wrapper¶
Parallelism¶
Block1.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 Block1.5, Block code introduces
threading level parallelism, more specifically, the OpenMP threading. To enable
the multithreading 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 multithreading with the multiprocessing model (mpirun np) potentially offers higher performance and better scaling for DMRG parallelism. It is recommended to enable the multithreading 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
#SBATCH ntaskspernode=4
#SBATCH cpuspertask=4
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 multithreading
automatically. You still need explicitly set mc.fcisolver.num_thrds = 4
in
the PySCF input script to turn on the multithreading 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, 1e5, 1e7, 1e7, 1e7]
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 1e5 0.0001
20 2000 1e7 0.0001
24 2000 1e7 0.0
30 2000 1e7 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 twodot algorithm
can effectively help DMRG solver moving out of the local minimum. Attribute
twodot_to_onedot
indicates when to switch to the onedot 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 1step DMRGCASSCF 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 includewarmup local_2sitenonspinadaptedfiedlerSee Keywords for the details of Block code keywords

configFile
¶ By default, keywords are written to file dmrg.conf

dmrg_switch_tol
¶ In 1step DMRGCASSCF, 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 1e3.

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 Block1.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 sameBLOCKSCRATCHDIR
settings, it’s easy to cause name conflicts on the scratch directory, especially when two DMRGCASSCF calculations are run on the same node.

spin
¶ 2S (= nelec_alpha  nelec_beta). If the argument
nelec
ofDMRGCI.kernel()
function is a twoitem 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 onedot 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 1particle density matrix from the directory indicated by
scratchDirectory
.

make_rdm12
(state, norb, nelec)¶ Given state ID, read its 1particle and 2particle density matrices from the directory indicated by
scratchDirectory
. Note the 2particle 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 3particle density matrices from the directory indicated by
scratchDirectory
. Note the 2particle 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 3particle 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 1particle density matrix from the directory indicated by
scratchDirectory
.

trans_rdm12
(statebra, stateket, norb, nelec)¶ Given the state ID of bra and ket, read the 1particle and 2particle density matrices from the directory indicated by
scratchDirectory
. Note the 2particle 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 2electron integrals (ijkl). 8fold permutation symmetry is required. The function returns the total energy and the state ID which is corresponding to the wavefunction 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 1step algorithm can applied efficiently in DMRGCASSCF 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
Stateaverage and statespecific DMRGCASCI/DMRGCASSCF¶
Stateaverage and statespecific calculations were also supported in the
DMRGCASCI/DMRGCASSCF through the BlockPySCF interface. The usage is
the same to that in regular CASCI/CASSCF calculation.
mc.stateaverage_()
function provides the average over the
multiple solutions over a single fcisolver
:
mc = dmrgscf.DMRGSCF(mf, 6, 6)
# halfhalf 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 halfhalf averaged. The two states are
saved on the disk indicated by mc.fcisolver.scratchDirectory
.
In many calculations, one would require the stateaverage 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 wavefunction 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.
Statespecific DMRGCASSCF is the other common calculation one would take.
Setting up statespecific DMRGCASSCF 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 DMRGCASSCF:
# 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 DMRGCASCI object as
well. However, a straightforward solution for DMRGCASCI 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 stateaverage and statespecific calculations.
DMRGNEVPT2¶
For Block 1.1.1 version or older, the standard DMRGNEVPT2 calculation can be carried out on top of the DMRGCASCI or DMRGCASSCF 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 DMRGNEVPT2 method requires the 4particle density matrix.
Computing and storing the 4particle 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 MPSperturber
technique which can significantly reduce the computation cost. The MPSperturber
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_adaptedserial"
Note
The wavefunction structure from different Block versions are incompatible. If BLOCKEXE for zeroth order wavefunction is set to Block1.1, BLOCKEXE_COMPRESS_NEVPT for PT should also be Block1.1. Similarly, Block1.5 (stackblock) PT code only compatible with the zeroth order wavefunction of Block1.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 Block1.5. Please refer to the Benchmark for approximate costs.
If the excitation energy is of interest, we can use DMRGNEVPT2 to compute the energy of excited state based on the multipleroot 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 DMRGNEVPT2
calculations are called based on the stateaverage DMRGCASSCF calculation, you
should be very careful with scratchDirectory
for the DMRG wave function
that NEVPT2 perturbation is applied on. In the multiplesolver stateaverage
DMRGCASSCF 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 = 'tzpdkh'
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 = 1e2
mf.kernel()
#
# second pass to converge SCF calculation
#
mf = scf.newton(mf)
mf.conv_tol = 1e12
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 UHFalpha/beta orbitals from chkfile
#
##################################################
fname = 'hs.chk'
chkfile = fname
mol, mf = scf.chkfile.load_scf(chkfile)
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.nelectronmol.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 bondRadiusMilliAngstroms 66;
set bondTolerance 0.5;
set forceAutoBond false;
#cd /home/abc/pyscf/examples/dmrg
''')
fspt.write('load %s_beta.molden;\n' % fname)
fspt.write('rotate 30 y;\n')
fspt.write('rotate 20 x;\n')
for i in range(nalpha10,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 "%salpha%d.jpg";\n' % (jmol_script, (i+1)))
fspt.write('load %s_alpha.molden;\n' % fname)
fspt.write('rotate 30 y;\n')
fspt.write('rotate 20 x;\n')
for i in range(nbeta10,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 "%sbeta%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()
#
# Nonorthogonal 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 = 1e7
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.e14]=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('CtSCI',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 <iFi>
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 <= 2thresh)
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('OrthocOAO',numpy.linalg.norm(numpy.dot(cOrbsOAO.T,cOrbsOAO)numpy.identity(nc)))
print('OrthoaOAO',numpy.linalg.norm(numpy.dot(aOrbsOAO.T,aOrbsOAO)numpy.identity(na)))
print('OrthovOAO',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><Oaux>, in which det(<Oaux>) 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 nonorthogonal AOs
# aux=numpy.identity(nb)
# * Lowdin orthogonalized AOs
aux = s12inv
# * Metalowdin 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 highlying orbitals standing in valence space.
#
# Get <iFi>
def psort(ova, fav, coeff):
# pT is density matrix, fav is Fock matrix
# OCCSORT
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
# ESORT
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
#
#==========================================
# lowdinpop 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[iorbnc],' faa=',e_o[iorbnc])
nelec += n_o[iorbnc]
else:
print(' iorb_V=',iorb,' occ=',n_v[iorbncna],' fvv=',e_v[iorbncna])
nelec += n_v[iorbncna]
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[iorbnc])
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] # S3p = 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] # FeMo (3d,4d) = 10
o2 = [ 1]*8+[0]*2 # approximate occupancies
a3 = [97,98,99,101,103,105] # Mos + 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, Highspin 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. DMRGCASSCF and DMRGNEVPT2
#
##################################################
#
# Adjust the MPI schedular and scratch directory if needed.
# NOTE the DMRGNEVPT2 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 = 1e6
orbs = mc.sort_mo(caslst, coeff, base=0)
mc.mc1step(orbs)
#
# CASCINEVPT2
#
# If DMRGCASSCF was finished without any problems (eg convergence, wall time
# limits on cluster etc), one can simply continue with DMRGNEVPT2
# 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 DMRGCASSCF has converged and the DMRG temporary files were
# deleted, we just need the DMRGCASCI calculation with the converged MCSCF
# orbitals to get the DMRG wavefunction.
#
mc = mcscf.CASCI(mf, norb, [nalpha,nbeta])
mc.chkfile = 'hs_mc.chk'
mo = scf.chkfile.load('hs_mc.chk', "mcscf/mo_coeff")
mc.fcisolver = DMRGCI(mol,maxM=500, tol =1e8)
#
# 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, 1e5, 1e6, 1e7, 1e7, 1e7, 1e7 ]
#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)
#
# DMRGNEVPT2
#
mrpt.NEVPT(mc).kernel()
#
# There is also a fast DMRGNEVPT2 implementation. See also the example
# pyscf/examples/dmrg/02dmrg_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 PySCF1.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. 
[2]  Please contact “Sandeep Sharma” <sanshar@gmail.com> or “Garnet Chan” <gkc1000@gmail.com> for Block source code. 