Solid state calculations

In this tutorial we will run periodic boundary conditions CCMC calculation on a diamond crystal in STO-3G basis with 2x1x1 sampling of the Brillouin zone. Familiarity with CCMC and FCIQMC tutorials is assumed. The input and output files can be found under the documentation/manual/tutorials/calcs/ccmc_solids subdirectory of the source distribution.

First of all, we need the one- and two- electron integrals from an external source. We will use PySCF software package [Sun18] to perform preliminary Hartree-Fock calculation and generate the integrals. PySCFDump script can be used to save the integrals in the FCIDUMP format readable by HANDE [1] .

To correctly address exchange divergence, additional exchange integrals are needed. Those are written in a FCIDUMP_X file by PySCFDump. Theoretical details of this procedure will be elaborated on as part of an upcoming paper.

Cell object can be conveniently prepared using build_cell function of the pyscfdump.helpers module and ASE library. Set of basis functions, kinetic energy cutoff and pseudopotential used are specified.

The run_khf function of the pyscfdump.scf module is used to run the HF calculation. Number of k-points in each dimension is specified. By default, the Monkhorst-Pack grid is used. If gamma is set to true, the grid will be shifted to include the Gamma point. Exchange divergence treatment scheme should also be chosen. The function returns the converged HF calculation object and a list of scaled k-points (i.e. (0.5,0.5,0.5) is at the very corner of the Brillouin zone).

Finally, fcidump function of the pyscfdump.pbcfcidump module is used to dump the integrals to a file. Name of the resulting file, SCF calculation object, number of k points, list of scaled k points and boolean variable MP (true if the grid is Monkhorst-Pack - i.e. not shifted to include the gamma point for even grids) must be provided.

from ase.lattice import bulk
import pyscfdump.scf as scf
from pyscfdump.helpers import build_cell
from pyscfdump.pbcfcidump import fcidump

A2B = 1.889725989	#angstrom to bohr conversion

a=3.567		#lattice parameter
ke=1000		#kinetic energy cutoff in units of rydberg
basis='sto-3g'	#basis set choice
nmp = [2,1,1]	#k point sampling

#prepare the cell object
ase_atom = bulk('C', 'diamond', a=a*A2B)
cell = build_cell(ase_atom, ke=ke, basis=basis, pseudo='gth-pade')

#run the HF calculation
kmf,scaled_kpts = scf.run_khf(cell, nmp=nmp, exxdiv='ewald', gamma=True)

#dump the integrals
fcidump('fcidumpfile',kmf,nmp,scaled_kpts,False)

Now we are ready to run the CCMC calculation. System definition is read in from the FCIDUMP files. There are two points to notice

  • path to files with additional exchange integrals is specified as ex_int_file
  • to properly exploit translational symmetry in the crystal lattice, the orbitals and hence the integrals must be complex. The complex mode of HANDE is enabled by setting complex = true. Consequently, numbers of particles on excitors also have both real and imaginary part, as well as projected energy.

In a non-initiator calculation, we try setting the time step tau as big as possible before too many blooms happen. A shoulder plot should be used to determine target_population. The heat_bath excitation generator [Holmes16] (adapted to HANDE as described in [Neufeld19]) is often a good choice in small systems. The input script used in this tutorial is:

sys = read_in {
    int_file = "fcidumpfile",
    complex = true,
    ex_int_file = "fcidumpfile_X",
}
ex_l=2

ccmc {
    sys = sys,
    qmc = {
        tau = 0.02,
        rng_seed = 13086,
        mc_cycles = 10,
        init_pop = 200,
        nreports = 50000,
        target_population = 1e4,
        state_size = -800,
        spawned_state_size = -500,
        excit_gen = "heat_bath",
        real_amplitudes = true,
    },
    ccmc = {
	even_selection = true,
        full_non_composite=true,
    },
    reference = {
        ex_level = ex_l,
    },
    blocking = {
        blocking_on_the_fly = true,
        auto_shift_damping = true,
    },
    restart = {
        write = true,
    },
}

The Hilbert space of the system we are dealing with is quite small and so is the target population. The general rule is that in order to use MPI parallelism, each process should contain at least \(10^{5}\) occupied excitors [Spencer18]. Having fewer excitors on each process is both inefficient and in extreme cases can lead to biased results. This is why we will use only one process here. However, use of openMP shared memory threads is recommended in order to make full use of the available resources.

$ hande.x ccmc.lua > diamond_ccmc.out

We can now plot the population

(Source code, png, hires.png, pdf)

../../_images/solids-1.png

and correlation energy:

(Source code, png, hires.png, pdf)

../../_images/solids-2.png

It is worth noting that the projected energy is in fact a complex quantity, whose imaginary part evaluates to zero in a non-trivial way [Booth13]. For the plot we calculate the result by only using real parts of both \(\sum_j H_{0j} N_j(\tau)\) and \(N_0(\tau)\), using the fact that

\(E(\tau) = \Re\left(\frac{\sum_j H_{0j} N_j(\tau)}{N_0(\tau)}\right) = \frac{\Re\left(\sum_j H_{0j} N_j(\tau)\right)}{\Re\left(N_0(\tau)\right)}\)

where the second equality holds provided that imaginary part of \(E(\tau)\) is zero.

The reblocking analysis uses magnitudes rather than real parts as this prevents problems with potential changes of phase during the calculation. Neverthless, for a well behaved calculation such as the one presented here, it is found that reblocking using real parts would give identical results.

In any case, the shift remains a strictly real measure of the correlation energy [2].

To analyse the calculation we can use reblock_hande.py script:

$ reblock_hande.py --quiet diamond_ccmc.out

which results in:

Block from: 20790
Recommended statistics from optimal block size:

                       Block from  # H psips \sum H_0j N_j      N_0        Shift Proj. Energy
diamond_ccmc_2.out 2.07900000e+04  33940(10)      -1015(2)  3328(6)  -0.30485(8)  -0.30495(5)

Footnotes

[1]The pyscfdump module code will be released shortly. For preliminary access contact Alex Thom at ajwt3@cam.ac.uk .
[2]As discussed in [Booth13] for FCIQMC - the CCMC case is exactly analogous.