Density-Matrix Quantum Monte Carlo Methods

Introduction

The Density-Matrix Quantum Monte Carlo (DMQMC) [1] methods stochastically sample the N-body thermal density matrix \(\hat{\rho} = e^{-\beta \hat{H}}\). This sampling allows exact thermodynamic properties to be computed for many-body systems at finite temperature.

The density matrix \(\hat{\rho}\) is represented in a discrete basis of Slater determinants. The matrix elements are stochastically sampled using a population of signed “psi-particles” or “psips”. The initial population of psips is initialized at \(\beta = 0\). As \(\beta\) is increased, the psips are propagated throughout the density matrix following a process of spawning, death, cloning, and annihilation events. Each of these actions has a probability determined by the Hamiltonian and current state of the density matrix as sampled by the psips.

The relationship \(\partial\hat{\rho}/\partial\beta\) between the density matrix and the inverse temperature can be constrained by either the symmetric or asymmetric Bloch equations. Both approaches are defined by their own classes in pydmqmc: SymmetricBlochDMQMC and AsymmetricBlochDMQMC. The way to use these classes is the same and will be covered below.

The pydmqmc library also offers Interaction Picture DMQMC (IP-DMQMC) [2]. As the name suggests, in this variation the density matrix \(\hat{\rho}\) is sampled in the interaction picture. This improves sampling—and statistical accuracy–for weakly correlated systems. The InteractionPictureDMQMC class enables use of IP-DMQMC.

Summary of Available Classes

Click on any entry for more details about the API or follow the guidelines below for setting up and running a DMQMC simulation.

• SymmetricBlochDMQMC - Performs DMQMC using the symmetric Bloch equation.

• AsymmetricBlochDMQMC - Performs DMQMC using the asymmetric Bloch equation.

• InteractionPictureDMQMC - Performs DMQMC in the interaction picture.

Setting up a DMQMC Simulation

Before starting with a DMQMC system, you must instantiate your desired System object. Once your system has been defined, you can create an instance of your desired DMQMC method. For this example, we’ll use SymmetricBlochDMQMC but this procedure is the same for all DMQMC classes listed above:

import pydmqmc.systems
from pydmqmc.methods import SymmetricBlochDMQMC

# load your system file with the appropriate class
sys = pydmqmc.systems...(...)

# use your system to instantiate your desired DMQMC method
mtd = SymmetricBlochDMQMC(sys)

Once your desired DMQMC method has been instantiated, you need to run the setup() method. This method is used to define the final inverse temperature \(\beta\) for the simulation, the framework for defining the initial density matrix, and the list of quantities to report during the iteration. You may need to specify additional parameters based on the initialization method chosen. An example of a simulation setup is below:

mtd.setup(
    final_beta = 0.75,
    initialization = "deterministic",
    report_quants = ["trace", "energy-expectation"]
)

For more information on reporting quantities as the simulation runs, see Iteration Progress Reports.

The setup() method can be called any number of times so long as the run() method has not been invoked. This allows you to change the final beta(), initial density matrix, or list of reported quantities before running the simulation. Neither attribute can be changed once the simulation has been run to prevent loss of data.

Density Matrix Initialization

The SymmetricBlochDMQMC and AsymmetricBlochDMQMC methods both offer the same set of frameworks for the initial density matrix: the identity matrix ("deterministic"), randomly placing a number of psips along the diagonal following a uniform distribution ("random-uniform"), and a user-supplied initial diagonal ("fixed").

The InteractionPictureDMQMC method also offers "deterministic" and "random-uniform" options but uses the Hartree-Fock thermal weights instead of integer psips. It also offers A "random-thermal" framework where the diagonal is sampled with probabilities proportional to the thermal weights and the "random-grand-canoncial" framework that uses the grand canonical density matrix to inform sampling.

Use the links above to see the appropriate API documentation for more details on these options and the additional parameters required.

Once the density matrix has been set up, you can view it with the density_matrix attribute:

print(mtd.density_matrix)

You can also double check the final beta() that was specified with the final_beta attribute:

print(mtd.final_beta)

Running a DMQMC Simulation

Once the density matrix and target inverse temperature have been specified with the setup() method, the run() method can be used to execute the actual DMQMC simulation. The run() method follows the same basic structure for all of the DMQMC classes; see their APIs for more detail. At minimum, all that’s needed to run a simulation is to specify three parameters that control the stability of the iterations as detailed in the next section:

mtd.run(
    dbeta=0.001,
    cycles_per_shift=10,
    shift_dampening=0.05
)

Integration Control

DMQMC simulations start from an inverse temperature of \(\beta = 0\) and increase it gradually to the final \(\beta\) specified during setup. Currently this is accomplished using a fixed step size \(\Delta\beta\) which is set by the required parameter dbeta.

In DMQMC, an energy shift \(S\) is applied to the Hamiltonian to keep the normalization approximately constant. This shift is updated every \(n\) cycles where \(n\) is set by the required cycles_per_shift parameter. How much the shift is allowed to vary with each update is controlled by the required shift_dampening parameter, which maps to the variable \(\zeta\) in Eqn 16 of the DMQMC paper [1]. The shift can either be specified for the Hamiltonian as a whole or independently for each row of Hamiltonian if shift_by_rows is True.

Since DMQMC is integrating \(\partial\hat{\rho}/\partial\beta\) to obtain the final density matrix \(\hat{\rho}\), you can specify the integration method with the update_method parameter. By default this parameter is set to “euler” but any of the integration methods recognized by pydmqmc.methods.Iterative.parse_method() are supported.

Note

Claire hasn’t actually tested that RK4 with DMQMC, though she’s tested both parts independently.

Psip Spawning Control

Note

This text is Claire’s best effort to understand these approximation methods based on conversations with Will. Corrections are welcome.

Several of the optional parameters for run() control how psips are spawned. The basic rules for spawning are as follows.

For a given site \(\rho_{ij}\) attempting to spawn a psip at \(\rho_{ik}\):

1. If \(\rho_{ik} \neq 0\), always let \(\rho_{ij}\) spawn psips.

2. If \(\rho_{ik} = 0\), let \(\rho_{ij}\) spawn if \(|\rho_{ij}| \ge \mathtt{n\_add}\).

3. If neither of the above are satisfied but another \(\rho_{pq}\), where \(p \neq i\) and \(q \neq j\), spawns psips at \(\rho_{ik}\) with the same sign, then \(\rho_{ij}\) may spawn as well.

4. If none of the above are satisfied, nullify the spawn from \(\rho_{ij}\) to \(\rho_{ik}\).

5. If the resulting \(\rho_{ik} \le \mathtt{spawn\_cutoff}\), nullify the spawn.

Warning

Claire isn’t sure rule 3 is currently implemented in the code.

The initiator approximation is represented by step 2 and may be disabled by setting n_add = None (\(n_\mathrm{add} = 0\)). This allows all spawns from \(\rho_{ik}\).

Similarly, the lower limit on psip spawning represented by step 5 can be disabled by setting spawn_cutoff = 0.

The ilevel parameter adds the initiator level approximation, which modifies the above rules to include another rule before 4: if the number of excitations between the determinant labels for \(i\) and \(j\) is less than or equal to the supplied level, allow the spawn at \(\rho_{ik}\). Currently, only initiator level zero is supported.

Saving Simulation Results

Both the final density matrix and the iteration report (the set of specified quantities that are calculated periodically as the simulation runs) can be saved to disk using one simple command:

mtd.save_data("my_sim")

This will create two files: my_sim_density_matrix.csv that contains the final density matrix and my_sim_report.csv that contains the iteration report. See the save_data() API for more information on the available options.

References

[1] (1,2)

N. S. Blunt, T. W. Rogers, J. S. Spencer, and W. M. C. Foulkes. Density-matrix quantum monte carlo method. Phys. Rev. B, 89:245124, Jun 2014. URL: https://link.aps.org/doi/10.1103/PhysRevB.89.245124, doi:10.1103/PhysRevB.89.245124.

[2]

Fionn D. Malone, N. S. Blunt, James J. Shepherd, D. K. K. Lee, J. S. Spencer, and W. M. C. Foulkes. Interaction picture density matrix quantum monte carlo. The Journal of Chemical Physics, 143(4):044116, 07 2015. URL: https://doi.org/10.1063/1.4927434, arXiv:https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4927434/15498760/044116_1_online.pdf, doi:10.1063/1.4927434.