DMRGSCF

In this section we explain how to use block2 (and optionally StackBlock) and pyscf for DMRGSCF (CASSCF with DMRG as the active space solver).

Preparation

One should first install the pyscf extension called dmrgscf, which can be obtained from https://github.com/pyscf/dmrgscf. If it is installed using pip, one also needs to create a file named settings.py under the dmrgscf folder, as follows:

$ pip install git+https://github.com/pyscf/dmrgscf
$ PYSCFHOME=$(pip show pyscf-dmrgscf | grep 'Location' | tr ' ' '\n' | tail -n 1)
$ wget https://raw.githubusercontent.com/pyscf/dmrgscf/master/pyscf/dmrgscf/settings.py.example
$ mv settings.py.example ${PYSCFHOME}/pyscf/dmrgscf/settings.py

Here we also assume that you have installed block2 either using pip or manually.

DMRGSCF (serial)

The following is an example python script for DMRGSCF using block2 running in a single node without MPI parallelism:

from pyscf import gto, scf, lib, dmrgscf
import os

dmrgscf.settings.BLOCKEXE = os.popen("which block2main").read().strip()
dmrgscf.settings.MPIPREFIX = ''

mol = gto.M(atom='C 0 0 0; C 0 0 1.2425', basis='ccpvdz',
    symmetry='d2h', verbose=4, max_memory=10000) # mem in MB
mf = scf.RHF(mol)
mf.kernel()

from pyscf.mcscf import avas
nactorb, nactelec, coeff = avas.avas(mf, ["C 2p", "C 3p", "C 2s", "C 3s"])
print('CAS = ', nactorb, nactelec)

mc = dmrgscf.DMRGSCF(mf, nactorb, nactelec, maxM=1000, tol=1E-10)
mc.fcisolver.runtimeDir = lib.param.TMPDIR
mc.fcisolver.scratchDirectory = lib.param.TMPDIR
mc.fcisolver.threads = int(os.environ.get("OMP_NUM_THREADS", 4))
mc.fcisolver.memory = int(mol.max_memory / 1000) # mem in GB

mc.canonicalization = True
mc.natorb = True
mc.kernel(coeff)

Note

Alternatively, to use StackBlock instead of block2 as the DMRG solver, one can change the line involving dmrgscf.settings.BLOCKEXE to:

dmrgscf.settings.BLOCKEXE = os.popen("which block.spin_adapted").read().strip()

Please see MPS Import/Export for the instruction for the installation of StackBlock.

Note

It is important to set a suitable mc.fcisolver.threads if you have multiple CPU cores in the node, to get high efficiency.

This will generate the following output:

$ grep 'CASSCF energy' cas1.out
CASSCF energy = -75.6231442712648

DMRGSCF (distributed parallel)

The following example is DMRGSCF in hybrid MPI (distributed) and openMP (shared memory) parallelism. For example, we can use 7 MPI processors and each processor uses 4 threads (so in total the calculation will be done with 28 CPU cores):

from pyscf import gto, scf, lib, dmrgscf
import os

dmrgscf.settings.BLOCKEXE = os.popen("which block2main").read().strip()
dmrgscf.settings.MPIPREFIX = 'mpirun -n 7 --bind-to none'

mol = gto.M(atom='C 0 0 0; C 0 0 1.2425', basis='ccpvdz',
    symmetry='d2h', verbose=4, max_memory=10000) # mem in MB
mf = scf.RHF(mol)
mf.kernel()

from pyscf.mcscf import avas
nactorb, nactelec, coeff = avas.avas(mf, ["C 2p", "C 3p", "C 2s", "C 3s"])
print('CAS = ', nactorb, nactelec)

mc = dmrgscf.DMRGSCF(mf, nactorb, nactelec, maxM=1000, tol=1E-10)
mc.fcisolver.runtimeDir = lib.param.TMPDIR
mc.fcisolver.scratchDirectory = lib.param.TMPDIR
mc.fcisolver.threads = 4
mc.fcisolver.memory = int(mol.max_memory / 1000) # mem in GB

mc.canonicalization = True
mc.natorb = True
mc.kernel(coeff)

Note

To use MPI with block2, the block2 must be either (a) installed using pip install block2-mpi or (b) manually built with -DMPI=ON. Note that the block2 installed using pip install block2 cannot be used together with mpirun if there are more than one processors (if this happens, it will generate wrong results and undefined behavior).

If you have already pip install block2, you must first pip uninstall block2 then pip install block2-mpi.

Note

If you do not have the --bind-to option in the mpirun command, sometimes every processor will only be able to use one thread (even if you set a larger number in the script), which will decrease the CPU usage and efficiency.

This will generate the following output:

$ grep 'CASSCF energy' cas2.out
CASSCF energy = -75.6231442712753

CASSCF Reference

For this small (8, 8) active space, we can also compare the above DMRG results with the CASSCF result:

from pyscf import gto, scf, lib, mcscf
import os

mol = gto.M(atom='C 0 0 0; C 0 0 1.2425', basis='ccpvdz',
    symmetry='d2h', verbose=4, max_memory=10000) # mem in MB
mf = scf.RHF(mol)
mf.kernel()

from pyscf.mcscf import avas
nactorb, nactelec, coeff = avas.avas(mf, ["C 2p", "C 3p", "C 2s", "C 3s"])
print('CAS = ', nactorb, nactelec)

mc = mcscf.CASSCF(mf, nactorb, nactelec)
mc.fcisolver.conv_tol = 1E-10
mc.canonicalization = True
mc.natorb = True
mc.kernel(coeff)

This will generate the following output:

$ grep 'CASSCF energy' cas3.out
CASSCF energy = -75.6231442712446

State-Average with Different Spins

The following is an example python script for state-averaged DMRGSCF with singlet and triplet:

from pyscf import gto, scf, lib, dmrgscf, mcscf
import os

dmrgscf.settings.BLOCKEXE = os.popen("which block2main").read().strip()
dmrgscf.settings.MPIPREFIX = ''

mol = gto.M(atom='C 0 0 0; C 0 0 1.2425', basis='ccpvdz',
    symmetry='d2h', verbose=4, max_memory=10000) # mem in MB
mf = scf.RHF(mol)
mf.kernel()

from pyscf.mcscf import avas
nactorb, nactelec, coeff = avas.avas(mf, ["C 2p", "C 3p", "C 2s", "C 3s"])
print('CAS = ', nactorb, nactelec)

lib.param.TMPDIR = os.path.abspath(lib.param.TMPDIR)

solvers = [dmrgscf.DMRGCI(mol, maxM=1000, tol=1E-10) for _ in range(2)]
weights = [1.0 / len(solvers)] * len(solvers)

solvers[0].spin = 0
solvers[1].spin = 2

for i, mcf in enumerate(solvers):
    mcf.runtimeDir = lib.param.TMPDIR + "/%d" % i
    mcf.scratchDirectory = lib.param.TMPDIR + "/%d" % i
    mcf.threads = 8
    mcf.memory = int(mol.max_memory / 1000) # mem in GB

mc = mcscf.CASSCF(mf, nactorb, nactelec)
mcscf.state_average_mix_(mc, solvers, weights)

mc.canonicalization = True
mc.natorb = True
mc.kernel(coeff)

Note

The mc parameter in the function state_average_mix_ must be a CASSCF object. It cannot be a DMRGSCF object (will produce a runtime error).

This will generate the following output:

$ grep 'State ' cas4.out
State 0 weight 0.5  E = -75.6175232350073 S^2 = 0.0000000
State 1 weight 0.5  E = -75.298522666384  S^2 = 2.0000000

Unrestricted DMRGSCF

One can also perform Unrestricted CASSCF (UCASSCF) with block2 using a UHF reference. Currently this is not directly supported by the pyscf/dmrgscf package, but here we can add some small modifications. The following is an example:

from pyscf import gto, scf, lib, dmrgscf, mcscf, fci
import os

dmrgscf.settings.BLOCKEXE = os.popen("which block2main").read().strip()
dmrgscf.settings.MPIPREFIX = ''

mol = gto.M(atom='C 0 0 0; C 0 0 1.2425', basis='ccpvdz',
    symmetry=False, verbose=4, max_memory=10000) # mem in MB
mf = scf.UHF(mol)
mf.kernel()

def write_uhf_fcidump(DMRGCI, h1e, g2e, n_sites, nelec, ecore=0, tol=1E-15):

    import numpy as np
    from pyscf import ao2mo
    from subprocess import check_call
    from block2 import FCIDUMP, VectorUInt8

    if isinstance(nelec, (int, np.integer)):
        na = nelec // 2 + nelec % 2
        nb = nelec - na
    else:
        na, nb = nelec

    assert isinstance(h1e, tuple) and len(h1e) == 2
    assert isinstance(g2e, tuple) and len(g2e) == 3

    mh1e_a = h1e[0][np.tril_indices(n_sites)]
    mh1e_b = h1e[1][np.tril_indices(n_sites)]
    mh1e_a[np.abs(mh1e_a) < tol] = 0.0
    mh1e_b[np.abs(mh1e_b) < tol] = 0.0

    g2e_aa = ao2mo.restore(8, g2e[0], n_sites)
    g2e_bb = ao2mo.restore(8, g2e[2], n_sites)
    g2e_ab = ao2mo.restore(4, g2e[1], n_sites)
    g2e_aa[np.abs(g2e_aa) < tol] = 0.0
    g2e_bb[np.abs(g2e_bb) < tol] = 0.0
    g2e_ab[np.abs(g2e_ab) < tol] = 0.0

    mh1e = (mh1e_a, mh1e_b)
    mg2e = (g2e_aa, g2e_bb, g2e_ab)

    cmd = ' '.join((DMRGCI.mpiprefix, "mkdir -p", DMRGCI.scratchDirectory))
    check_call(cmd, shell=True)
    if not os.path.exists(DMRGCI.runtimeDir):
        os.makedirs(DMRGCI.runtimeDir)

    fd = FCIDUMP()
    fd.initialize_sz(n_sites, na + nb, na - nb, 1, ecore, mh1e, mg2e)
    fd.orb_sym = VectorUInt8([1] * n_sites)
    integral_file = os.path.join(DMRGCI.runtimeDir, DMRGCI.integralFile)
    fd.write(integral_file)
    DMRGCI.groupname = None
    DMRGCI.nonspinAdapted = True
    return integral_file

def make_rdm12s(DMRGCI, state, norb, nelec, **kwargs):

    import numpy as np

    if isinstance(nelec, (int, np.integer)):
        na = nelec // 2 + nelec % 2
        nb = nelec - na
    else:
        na, nb = nelec

    file2pdm = "2pdm-%d-%d.npy" % (state, state) if DMRGCI.nroots > 1 else "2pdm.npy"
    dm2 = np.load(os.path.join(DMRGCI.scratchDirectory, "node0", file2pdm))
    dm2 = dm2.transpose(0, 1, 4, 2, 3)
    dm1a = np.einsum('ikjj->ki', dm2[0]) / (na - 1)
    dm1b = np.einsum('ikjj->ki', dm2[2]) / (nb - 1)

    return (dm1a, dm1b), dm2

dmrgscf.dmrgci.writeIntegralFile = write_uhf_fcidump
dmrgscf.DMRGCI.make_rdm12s = make_rdm12s

mc = mcscf.UCASSCF(mf, 8, 8)
mc.fcisolver = dmrgscf.DMRGCI(mol, maxM=1000, tol=1E-7)
mc.fcisolver.runtimeDir = lib.param.TMPDIR
mc.fcisolver.scratchDirectory = lib.param.TMPDIR
mc.fcisolver.threads = int(os.environ["OMP_NUM_THREADS"])
mc.fcisolver.memory = int(mol.max_memory / 1000) # mem in GB

mc.canonicalization = True
mc.natorb = True
mc.kernel()

Note

In the above example, mf is the UHF object and mc is the UCASSCF object. It is important to ensure that both of them are with unrestricted orbitals. Otherwise the calculation may be done with only restricted orbitals. DMRGSCF wrapper cannot be used for this example.

Note

Due to limitations in pyscf/UCASCI, currently the point group symmetry is not supported in UCASSCF/UCASCI with DMRG solver. pyscf/avas does not support creating active space with unrestricted orbtials so here we did not use avas. The above example will not work with StackBlock (the compatibility with StackBlock will be considered in future).

This will generate the following output:

$ grep 'UCASSCF energy' cas5.out
UCASSCF energy = -75.6231442541606

UCASSCF Reference

We compare the above DMRG results with the UCASSCF result using the FCI solver:

mc = mcscf.UCASSCF(mf, 8, 8)
mc.fcisolver.conv_tol = 1E-10
mc.canonicalization = True
mc.natorb = True
mc.kernel(coeff)

This will generate the following output:

$ grep 'UCASSCF energy' cas6.out
UCASSCF energy = -75.6231442706386

DMRGSCF Nuclear Gradients and Geometry Optimization

The following is an example python script for computing DMRGSCF nuclear gradients and geometry optimization using block2:

from pyscf import gto, scf, lib, dmrgscf
import os

dmrgscf.settings.BLOCKEXE = os.popen("which block2main").read().strip()
dmrgscf.settings.MPIPREFIX = ''

mol = gto.M(atom='C 0 0 0; C 0 0 1.2425', basis='ccpvdz',
    symmetry='d2h', verbose=4, max_memory=10000) # mem in MB
mf = scf.RHF(mol)
mf.kernel()

from pyscf.mcscf import avas
nactorb, nactelec, coeff = avas.avas(mf, ["C 2p", "C 3p", "C 2s", "C 3s"])
print('CAS = ', nactorb, nactelec)

mc = mcscf.CASSCF(mf, nactorb, nactelec)
mc.fcisolver = dmrgscf.DMRGCI(mol, maxM=1000, tol=1E-10)
mc.fcisolver.runtimeDir = lib.param.TMPDIR
mc.fcisolver.scratchDirectory = lib.param.TMPDIR
mc.fcisolver.threads = int(os.environ.get("OMP_NUM_THREADS", 4))
mc.fcisolver.memory = int(mol.max_memory / 1000) # mem in GB

mc.canonicalization = True
mc.natorb = True
mc.kernel(coeff)

grad = mc.nuc_grad_method().kernel()

mol_eq = mc.nuc_grad_method().optimizer(solver='geomeTRIC').kernel()
print(mol_eq.atom_coords())

This will generate the following output (the nuclear gradient at the initial geometry and the optimized geometry):

$ grep -A 4 'SymAdaptedCASSCF gradients' cas7.out
--------------- SymAdaptedCASSCF gradients ---------------
        x                y                z
0 C     0.0000000000     0.0000000000     0.0388202961
1 C     0.0000000000     0.0000000000    -0.0388202961
----------------------------------------------
$ tail -n 3 cas7.out
cycle 3: E = -75.6240204052  dE = -5.51573e-07  norm(grad) = 9.37108e-05
[[ 0.          0.         -1.19709701]
[ 0.          0.          1.19709701]]

Note

Currently, gradients for UCASSCF is not supported in pyscf. The geometry optimization part requires an additional module called geomeTRIC, which can be installed via pip install geometric.