DMRGSCF (OpenMOLCAS)
In this section we explain how to use block2 and OpenMOLCAS for CASSCF and CASPT2 with DMRG as the active space solver.
Preparation
First, make sure block2 is installed correctly (either compiled manually or installed using pip,
and for pip the version of block2 should be >=0.5.1rc17),
so that the command which block2main can print a valid file path to block2main.
For example, the required block2 can be installed using pip as:
pip install block2>=0.5.1rc17 --extra-index-url=https://block-hczhai.github.io/block2-preview/pypi/
Then we need to compile an OpenMOLCAS with the block2 interface.
The source code of the required OpenMOLCAS code can be found in https://github.com/hczhai/OpenMolcas, which is
a slightly modified version of https://github.com/quanp/OpenMolcas (the OpenMOLCAS interface for the block 1.5 and StackBlock).
To activate the blcok2 interface, run cmake for this OpenMOLCAS with the option -DBLOCK2=ON. The detailed procedure is as follows:
git clone https://github.com/hczhai/OpenMolcas
export MOLCASHOME=$PWD/OpenMolcas
cd OpenMolcas
mkdir build
cd build
CC=gcc CXX=g++ FC=gfortran MKLROOT=/usr/local cmake .. -DCMAKE_INSTALL_PREFIX=../install -DLINALG=MKL -DOPENMP=ON -DBLOCK2=ON
make -j 10
make install
Remember to change the MKLROOT variable in the above example for your case.
Then one can run OpenMolcas using the following command:
MOLCAS=$MOLCASHOME/install MOLCAS_WORKDIR=/content/tmp pymolcas test.in
Where test.in is an OpenMolcas input file.
Sometimes you may need to add the --not-here option to pymolcas if it cannot find the molcas executable.
DMRGSCF
The following is an example input file for DMRGSCF for a O2 triplet state (see DMRGSCF (pyscf) for the similar calculation using pyscf):
&GATEWAY
Title
O2 Molecule
Coord
2
O 0 0 -0.6035
O 0 0 0.6035
Basis set
CC-PVDZ
&SEWARD
&SCF
Spin = 1
&RASSCF
Spin
3
Symmetry
4
nActEl
2 0 0
Inactive
3 1 1 0 2 0 0 0
Ras2
0 0 0 0 0 1 1 0
CIROOT = 1 1 ; 1
&RASSCF
Spin
3
Symmetry
4
nActEl
8 0 0
Inactive
2 0 0 0 2 0 0 0
Ras2
1 1 1 0 1 1 1 0
CIROOT = 1 1 ; 1
CISOlver = BLOCK
DMRG = 1000
Note that the first RASSCF is actually a ROHF mean-field calculation.
The same calculation in pyscf is:
from pyscf import gto, scf, mcscf, mrpt, dmrgscf, lib, symm
from pyblock2._pyscf.ao2mo import integrals as itg
import os
mol = gto.M(atom='O 0 0 0; O 0 0 1.207', basis='cc-pvdz', spin=2, symmetry='d2h', cart=False, verbose=4)
mf = scf.RHF(mol).run(conv_tol=1E-20)
ncas, n_elec, spin, ecore, h1e, g2e, orb_sym = itg.get_rhf_integrals(mf, g2e_symm=8)
print(orb_sym)
print(mf.mo_occ)
orb_sym_name = [symm.irrep_id2name(mol.groupname, ir) for ir in orb_sym]
print(orb_sym_name)
mc = mcscf.CASSCF(mf, 6, 8)
mc.fcisolver.conv_tol = 1e-14
mc.canonicalization = True
mc.natorb = True
mc.run()
From the pyscf output we can see the occupation number and orbtial irreps are :
[0, 5, 0, 5, 0, 6, 7, 2, 3, 5, 5, 6, 7, 0, 2, 3, 0, 5, 6, 7, 0, 1, 4, 5, 0, 2, 3, 5] # XOR irreps
[2. 2. 2. 2. 2. 2. 2. 1. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] # occ
['Ag', 'B1u', 'Ag', 'B1u', 'Ag', 'B2u', 'B3u', 'B2g', 'B3g', 'B1u', 'B1u', 'B2u', 'B3u', 'Ag', 'B2g', 'B3g', 'Ag', 'B1u', 'B2u', 'B3u', 'Ag', 'B1g', 'Au', 'B1u', 'Ag', 'B2g', 'B3g', 'B1u']
The MOLCAS ordering of irreps of D2h is:
ag b3u b2u b1g b1u b2g b3g au
This information can help us setting the Inactive and Ras2 in the MOLCAS inputfile.
From the pyscf output we have:
$ grep 'converged SCF energy' pyscf.out
converged SCF energy = -149.608181589162
$ grep 'CASSCF energy' pyscf.out
CASSCF energy = -149.708657770064
From the openMOLCAS output we have:
$ grep ':: RASSCF' o2.out
:: RASSCF root number 1 Total energy: -149.60818159
:: RASSCF root number 1 Total energy: -149.70865773
Note that in the openMOLCAS output, the first line is actually the SCF (ROHF) energy,
and the second line is the CASSCF energy. So they are consistent.
DMRG-cu-CASPT2
The following is an example input file for CASPT2 calculation after DMRGSCF for a O2 triplet state. In this example, the cumulant approximation of 4PDM is used for CASPT2. Note that the IPEA shift = 0.25 is used by default.
&GATEWAY
Title
O2 Molecule
Coord
2
O 0 0 -0.6035
O 0 0 0.6035
Basis set
CC-PVDZ
&SEWARD
&SCF
Spin = 1
&RASSCF
Spin
3
Symmetry
4
nActEl
2 0 0
Inactive
3 1 1 0 2 0 0 0
Ras2
0 0 0 0 0 1 1 0
CIROOT = 1 1 ; 1
&RASSCF
Spin
3
Symmetry
4
nActEl
8 0 0
Inactive
2 0 0 0 2 0 0 0
Ras2
1 1 1 0 1 1 1 0
CIROOT = 1 1 ; 1
CISOlver = BLOCK
DMRG = 1000
3RDM
NO4R
&CASPT2
MULT = 1 1
CUMU
The keyword NO4R is required in the RASSCF section to avoid spending time on computing 4pdms.
This will generate the following output:
$ grep ':: CASPT2' o2.out
:: CASPT2 Root 1 Total energy: -149.97055932
DMRG-CASPT2
The following is an example input file for CASPT2 calculation after DMRGSCF for a O2 triplet state. In this example, the exact 4PDM is computed and used.
&GATEWAY
Title
O2 Molecule
Coord
2
O 0 0 -0.6035
O 0 0 0.6035
Basis set
CC-PVDZ
&SEWARD
&SCF
Spin = 1
&RASSCF
Spin
3
Symmetry
4
nActEl
2 0 0
Inactive
3 1 1 0 2 0 0 0
Ras2
0 0 0 0 0 1 1 0
CIROOT = 1 1 ; 1
&RASSCF
Spin
3
Symmetry
4
nActEl
8 0 0
Inactive
2 0 0 0 2 0 0 0
Ras2
1 1 1 0 1 1 1 0
CIROOT = 1 1 ; 1
CISOlver = BLOCK
DMRG = 1000
3RDM
&CASPT2
BLOCK
MULT = 1 1
In the above example, we use the keyword BLOCK to replace the old keyword CUMU
so that the cumulant approximation is not used.
Note
By default there will be frozen orbitals in the CASPT2 treatment. One can add
FROZEN
0 0 0 0 0 0 0 0
in the CASPT2 section in the above example to avoid frozen orbitals.
This will generate the following output:
$ grep ':: CASPT2' o2.out
:: CASPT2 Root 1 Total energy: -149.96959847
State-Average
The following is an example input file for state-averaged DMRGSCF for three states, and then the CASPT2 treatment of each of the three states. In this example, the exact 4PDM is computed and used.
&GATEWAY
Title
O2 Molecule
Coord
2
O 0 0 -0.6035
O 0 0 0.6035
Basis set
CC-PVDZ
&SEWARD
&SCF
Spin = 1
&RASSCF
Spin
3
Symmetry
4
nActEl
2 0 0
Inactive
3 1 1 0 2 0 0 0
Ras2
0 0 0 0 0 1 1 0
CIROOT = 1 1 ; 1
&RASSCF
Spin
3
Symmetry
4
nActEl
8 0 0
Inactive
2 0 0 0 2 0 0 0
Ras2
1 1 1 0 1 1 1 0
CIROOT = 3 3 1
CISOlver = BLOCK
DMRG = 1000
3RDM
&CASPT2
BLOCK
MULT = 1 1
&CASPT2
BLOCK
MULT = 1 2
&CASPT2
BLOCK
MULT = 1 3
From the output we have:
$ grep ':: RASSCF' o2.out
:: RASSCF root number 1 Total energy: -149.60818159
:: RASSCF root number 1 Total energy: -149.69063345
:: RASSCF root number 2 Total energy: -149.09370540
:: RASSCF root number 3 Total energy: -148.86158577
$ grep ':: CASPT2' o2.out
:: CASPT2 Root 1 Total energy: -149.96175902
:: CASPT2 Root 1 Total energy: -149.39685470
:: CASPT2 Root 1 Total energy: -149.13012648