Density Matrix Quantum Monte Carlo¶
dmqmc {
sys = system,
qmc = { ... },
dmqmc = { ... },
ipdmqmc = { ... },
operators = { ... },
rdm = { ... },
restart = { ... },
reference = { ... },
qmc_state = qmc_state,
}
- Returns:
- a qmc_state object.
a lua table containing the sampling probabilities found if
find_weightsis set totrue. This can be passed directly to theweightsoption of a subsequent DMQMC calculation and/or manipulated inside the lua script. Iffind_weightsisfalse, only the qmc_state object is returned.
dmqmc performs a density matrix quantum Monte Carlo (DMQMC) calculation on a system.
Unlike Coupled Cluster Monte Carlo and Full Configuration Interaction Quantum Monte Carlo, where quantities are averaged inside each report loop, any quantities in DMQMC are evaluated at the first iteration of the report loop only. This is because different iterations represent different temperatures in DMQMC, and so averaging over a report loop would average over different temperatures, which is not the desired behaviour.
Note
Density Matrix Quantum Monte Carlo is currently rather experimental. In particular, it is not implemented for all systems yet and some options are only implemented for specific systems. In particular, DMQMC is only implemented for the Heisenberg model, the UEG, the real and momentum-space Hubbard model, and for molecular systems. The evaluation of operators other than the total energy, such as correlation functions and entanglement measures, is currently only possible for the Heisenberg model. The calculation of the reduced density matrices from DMQMC is also only supported for the Heisenberg model (for both temperature-dependent and ground state RDM calculations).
Options¶
systype: system object.
Required.
The system on which to perform the calculation. Must be created via a system function.
qmctype: lua table.
Required.
Further options that are common to all implemented QMC algorithms. See qmc options.
dmqmctype: lua table.
Optional.
Further options to control the DMQMC algorithm. See dmqmc options.
ipdmqmctype: lua table.
Optional.
If set, even to an empty table, then interaction picture DMQMC [Malone15] is performed. The table can contain further options to control the IP-DMQMC algorithm. See ipdmqmc options.
operatorstype: lua table.
Optional.
Further options to select the operators for which expectation values are evaluated. See operators options.
rdmtype: lua table.
Optional.
Further options to select which (if any) reduced density matrices and corresponding operators are to be evaluated. See rdm options.
restarttype: lua table.
Optional.
Further options to control restarting the calculation from a previous calculation. See restart options.
referencetype: lua table.
Optional.
Further options to select the reference state used. See reference options.
qmc_statetype: qmc_state object.
Optional.
Output of a previous calculation to resume.
Warning
The qmc_state object must have been returned by a previous DMQMC calculation. The validity of this is not checked. The system must also be unchanged.
Warning
This destroys the qmc_state object and so it cannot be re-used in subsequent QMC calculations.
dmqmc options¶
symmetric_blochtype: boolean.
Optional. Default: true.
Use the symmetrized form of the Bloch equation, \(\frac{d\hat{\rho}}{d\beta}=-\frac{1}{2}\{\hat{H},\hat{\rho}\}\), to propagate the density matrix when true. Otherwise the non-symmetrized form, \(\frac{d\hat{\rho}}{d\beta}=-\hat{\rho} \hat{H}\) of the Bloch equation is used. The symmetrize option only works with the symmetric version of the Bloch equation.
Note
The use of symmetric or asymmetric propagation in the DMQMC methods can impact the behavior of the sign problem as well the convergence with respect to beta loops. For more information see [Petras21].
replica_trickstype: boolean.
Optional. Default: false.
Perform replica simulations (i.e. evolve two independent DMQMC simulations concurrently) if true. This allows calculation of unbiased estimators that are quadratic in the density matrix.
fermi_temperaturetype: boolean.
Optional. Default: false.
Rescale tau so that the simulation runs in timesteps of \(\Delta\tau / T_F\) where \(T_F\) is the Fermi temperature. This is so results are at dimensionless inverse temperatures of \(\Theta^{-1} =T_F/T\). This option is only valid for systems with a well defined Fermi energy.
all_sym_sectorstype: boolean.
Optional. Default: false.
Sample states with all symmetries of the system instead of just those which conserve the symmetry of the reference state.
all_spin_sectorstype: boolean.
Optional. Default: false.
Sample states with all spin polarisations of the system instead of just those which conserve the spin polarisation of the reference state.
beta_loopstype: integer.
Optional. Default: 100.
The number of loops over the desired temperature range (each starting from \(T=\infty\) and performing the desired number of iterations) to perform. Each beta loop samples the initial conditions independently.
Note
Estimators must be averaged at each temperature from different beta loops. As each beta loop is independent, this can be done in separate calculations in an embararassingly parallel fashion.
final_betatype: float.
Optional. Default: 0.0.
Sets the final inverse temperature the density matrix is propagated to. When not provided in the input, the number of reports and Monte Carlo cycles controls the final temperature instead. If specified while using the interaction picture, the interaction picture and Bloch equation are used in a piecewise fashion to sample a range of temperatures from target_beta to final_beta. [VanBenschoten22]
sampling_weightstype: vector of floats.
Optional. Default: none.
Specify factors used to alter the spawning probabilities in the DMQMC importance sampling procedure. See PRB, 89, 245124 (2014) for an explanation, in particular section IV and appendix B.
The length of the vector of floats should be equal to the maximum number of excitations from any determinant in the space. For a chemical system with \(N\) electrons and more than \(2N\) spin orbitals, this would be equal to \(N\). For a Heisenberg model with \(N\) spins in the \(M_s=0\) spin sector, this should be equal to \(N/2\) (each pair of opposite spins flipped is one excitation).
vary_weightstype: integer.
Optional. Default: 0
The number of iterations over which to introduce the weights in the importance sampling scheme (see PRB, 89, 245124 (2014)). If not set then the full weights will be used from the first iteration. Otherwise, the weights will be increased by a factor of \((W_{\gamma})^{\beta/\beta_{target}}\) each iteration, where \(W_{\gamma}\) is the final weight of excitation level \(\gamma\) and \(\beta_{target}\) is the beta value to vary the weights until (equal to the value specified by this option, multiplied by the time step size).
find_weightstype: boolean.
Optional. Default: false.
Run a simulation to attempt to find appropriate weights for use in the DMQMC importance sampling procedure. This algorithm will attempt to find weights such that the population of psips is evenly distributed among the various excitation levels when the ground state is reached (at large beta values). The algorithm should be run for several beta loops until the weights settle down to a roughly constant value.
The weights are output at the end of each beta loop.
This option should be used together with the find_weights_start option, which is used to specify at which iteration the ground state is reached and therefore when averaging of the excitation distribution begins.
This option cannot be used together with the excit_dist option. The find_weights option averages the excitation distribution in the ground state, whereas the excit_dist option accumulates and prints out the excitation distribution at every report loop.
Warning
This feature is found to be unsuccessful for some larger lattices (for example, 6x6x6, for the Heisenberg model). The weights output should be checked. Increasing the number of psips used may improve the weights calculated.
find_weights_starttype: integer.
Optional. Default: 0.
The iteration number at which averaging of the excitation distribution begins, when using the find_weights option.
symmetrizetype: boolean.
Optional. Default: false.
Explicitly symmetrize the density matrix, thus only sampling one triangle of the matrix. This can yield significant improvements in stochastic error in some cases.
initiator_leveltype: integer.
Optional. Default: -1.
Set all density matrix elements at excitation level initiator_level and below to be initiator determinants. An initiator_level of -1 indicates that no preferential treatment is given to density matrix elements and the usual initiator approximation is imposed, 0 indicates that the diagonal elements are initiators, etc.
This is experimental and the user should identity when convergence has been reached.
piecewise_shifttype: float.
Optional. Default: 0.
Sets the value of the simulation shift when the propagator change occurs at the target_beta when running the piecewise interaction picture method.
walker_scale_factortype: integer.
Optional. Default: 1.
Scales the walker population on the initial trial density matrix by a constant factor. The simulations target population is scaled as well.
Warning
This feature is experimental, and results should be tested for accuracy.
ipdmqmc options¶
target_betatype: float.
Optional. Default: 1.0.
The inverse temperature to propagate the density matrix to. If fermi_temperature is set to True then target_beta is interpreted as the inverse reduced temperature \(\tilde{\beta} = 1/\Theta = T_F/T\), where \(T_F\) is the Fermi temperature. Otherwise target_beta is taken to be in atomic units.
Note
If final_beta is set to a value greater than target_beta, the interaction picture will be used until the target_beta has been reached. Thereafter, the Bloch equation will be used to sample continously until the final_beta has been reached.
initial_matrixtype: string.
Optional. Default: ‘hartree_fock’.
Possible values: ‘free_electron’, ‘hartree_fock’.
Initialisation of the density matrix at \(\tau=0\). ‘free_electron’ samples the free electron density matrix, i.e. \(\hat{\rho} = \sum_i e^{-\beta \sum_j \varepsilon_j \hat{n}_j} |D_i\rangle\langle D_i|\), where \(\varepsilon_j\) is the single-particle eigenvalue and \(\hat{n}_j\) the corresponding number operator. ‘hartree_fock’ samples a ‘Hartree–Fock’ density matrix defined by \(\hat{\rho} = \sum e^{-\beta H_{ii}} |D_i\rangle\langle D_i|\), where \(H_{ii} = \langle D_i|\hat{H}|D_i\rangle\).
It is normally best to use the hartree-fock option as this removes cloning/death on the diagonal if the shift is fixed at zero. This requires slightly more work when also using the grand_canonical_initialisation, but this is negligeable.
grand_canonical_initialisationtype: boolean.
Optional. Default: false.
Use the grand canonical partition function to initialise the psip distribution. The default behaviour will randomly distribute particles among the determinants requiring a non-zero value of metropolis_attempts to be set for the correct distribution to be reached.
skip_gci_reference_checktype: boolean.
Optional. Default: false.
When performing grand_canonical_initialisation, we check that \(H_{ii}\) is not lower in energy than \(H_{00}\). If a lower energy \(H_{ii}\) is found this can cause many spawns to occur with a weight lower than 1.0 which is undesirable, and so the simulation exits with information to update the reference. Setting this flag to true will ignore the lower energy \(H_{ii}\).
Warning
It is recommended that the orbital single particle eigenvalues in the FCIDUMP are recalculated with the new reference.
metropolis_attemptstype: integer.
Optional. Default: 0.
Number of Metropolis moves to perform (per particle) on the initial distribution. It is up to the user to determine if the desired distribution has been reached, i.e. by checking if results are independent of metropolis_attempts.
symmetric_interaction_picturetype: boolean.
Optional. Default: true.
Use symmetric version of ip-dmqmc where now \(\hat{f}(\tau) = e^{-\frac{1}{2}(\beta-\tau)\hat{H}^0}e^{-\tau\hat{H}}e^{-\frac{1}{2}(\beta-\tau)\hat{H}^0}\).
Warning
This feature is experimental and only tested for the 3D uniform electron gas.
count_diagonal_occupationstype: boolean.
Optional. Default: false.
When performing grand_canonical_initialisation, instead of accumulating the number of walkers being added to the trace count the number of diagonal elements that are occupied. The original grand_canonical_initialisation would count the number of successful occupations which could lead to substantially more particles being added then the provided initial population. Generally only applicable when initial_matrix is set to ‘hartree_fock’.
operators options¶
renyi2type: boolean.
Optional. Default: false.
Calculate the Renyi-2 entropy of the entire system. Requires
replica_tricksto be enabled.energytype: boolean.
Optional. Default: false.
Calculate the thermal expectation value of the Hamiltonian operator.
energy2type: boolean.
Optional. Default: false.
Calculate the thermal expectation value of the Hamiltonian operator squared. Only available for the Heisenberg model.
staggered_magnetisationtype: boolean.
Optional. Default: false.
Calculate the thermal expectation value of the staggered magnetisation operator. Only available for the Heisenberg model and with bipartite lattices.
excit_disttype: boolean.
Optional. Default: false.
Calculate the fraction of psips at each excitation level, where the excitation level is the number of excitations separating the two states labelling a given density matrix element. This fraction is then output to the data table at each report loop, and so the temperature-dependent excitation distribution is printed out.
This option should not be used with the find_weights option, which averages the excitation distribution within the ground state.
correlationtype: 2D vector of integers.
Optional. Default: false.
Calculate the spin-spin correlation function between the two specified lattice sites, \(i\) and \(j\), which is defined as the thermal expectation value of:
\[\hat{C}_{ij} = \hat{S}_{xi}\hat{S}_{xj} + \hat{S}_{yi}\hat{S}_{yj} + \hat{S}_{zi}\hat{S}_{zj}.\]Only available for the Heisenberg model.
potential_energytype: boolean
Optional. Default: false
Evaluate the bare Coulomb energy. Only available for the UEG.
kinetic_energytype: boolean
Optional. Default: false
Evaluate the kinetic energy. Only available for the UEG.
H0_energytype: boolean
Optional. Default: false
Evaluate the thermal expectation value of the zeroth order Hamiltonian where \(\hat{H} = \hat{H}^0 + \hat{V}\). See initial_matrix option. Only available when using the ip-dmqmc algorithm.
HI_energyEvaluate the expectation value of the interaction picture Hamiltonian where
\[\hat{H}_I\left(\frac{1}{2}(\beta-\tau)\right) = e^{\frac{1}{2}(\beta-\tau)\hat{H}^0}\hat{H}e^{-\frac{1}{2}(\beta-\tau)\hat{H}^0}.\]mom_disttype: float
Optional. Default: 0.0
Evaluate the (spin averaged) momentum distribution in kspace, i.e., \(\langle \hat{n}_{\mathbf{k}} \rangle\), up to a maximum wavevector defined by kmax which is a multiple of the Fermi wavevector. The momentum distribution will be printed out at unique kpoints which have the same kinetic energy. Results can be extracted from the analysed (i.e. by using the finite_temp_analysis script in the tools/dmqmc (see tutorial for more information)) dmqmc output using the extract_momentum_correlation.py script in the tools/dmqmc directory.
Only currently implemented for the UEG.
structure_factortype: float
Optional. Default: 0.0
Evaluate the static structure factor:
\[S_{\sigma\sigma'}(q) = \frac{N_{\sigma}\delta_{\sigma\sigma'}}{N} + \frac{1}{N} \sum_{kp} \left\langle c^{\dagger}_{k+q\sigma}c^{\dagger}_{p-q\sigma'} c_{p\sigma'}c_{k\sigma}\right\rangle\]up to a maximum wavevector defined by qmax which is a multiple of the Fermi wavevector. The static structure factor will be printed out at unique kpoints which have the same kinetic energy. Note that in the output file we actually print out \(S(q)-1\), \(S_{\uparrow\uparrow}(q)+S_{\downarrow\downarrow}(q)-1\) and \(S_{\uparrow\downarrow}(q)+S_{\downarrow\uparrow}(q)\), where \(S(q) = \sum_{\sigma\sigma'}S_{\sigma\sigma'}\). Results can be extracted from the analysed (i.e. by using the finite_temp_analysis script in the tools/dmqmc (see tutorial for more information)) dmqmc output using the extract_momentum_correlation.py script in the tools/dmqmc directory. The extraction script takes care of the factors of 1.
Currently only implemented for the UEG.
rdm options¶
Note that the use of RDMs is currently only available with the Heisenberg model.
rdmstype: table of 1D vectors.
Required.
Each vector corresponds to the subsystem of a reduced density matrix as a list of the basis function indices in the subsystem. For example:
rdms = { { 1, 2 } }
specifies one RDM containing basis functions with indices 1 and 2, and
rdms = { { 1, 2 }, { 3, 4} }
specifies two RDMs, with the first containing basis functions with indices 1 and 2, and the second basis functions 3 and 4.
Either
instantaneousorground_statemust be enabled to set the desired mode of evaluating the RDM (but both options cannot be used together).instantaneoustype: boolean.
Optional. Default: false.
Calculate the RDMs at each temperature based upon the instantaneous psip distribution.
Cannot be used with the ground_state option (either ground_state or instantaneous RDMs can be calculated, but not both concurrently).
ground_statetype: boolean.
Optional. Default: false.
Accumulate the RDM once the ground state (as specified by
ground_state_start) is reached. This has two limitations: only one RDM can be accumulated in a calculation and the subsystem should be at most half the size of the system (which is always sufficient for ground-state calculations).Cannot be used with the instantaneous option (either ground_state or instantaneous RDMs can be calculated, but not both concurrently).
spawned_state_sizetype: integer.
Required if
instantaneousis true. Ignored otherwise.Maximum number of states (i.e. reduced density matrix elements) to store in the “spawned” list, which limits the number of unique RDM elements that each processor can set. Should be a sizeable fraction of
state_size(see qmc options) and depends on the size of the subsystem compared to the full space.Note
This is a per processor quantity. It is usually safe to assume that each processor has approximately the same number of states.
ground_state_starttype: integer.
Optional. Default: 0.
Monte Carlo cycle from which the RDM is to be accumulated in each beta loop. Relevant only if
ground_stateis set to true and, as such, should be set to an iteration (which is a measure of temperature) such that the system has reached the ground state.concurrencetype: boolean.
Optional. Default: false.
Calculate the unnormalised concurrence and the trace of the reduced density matrix at the end of each beta loop. The normalised concurrence can be calculated from this using the
average_entropy.pyscript.Valid for
ground_stateonly; temperature-dependent concurrence is not currently implemented.renyi2type: boolean.
Optional. Default: false.
Calculate the Renyi-2 entropy of each subsystem. More accurately, the quantity output to the data table is \(S^n_2 = \sum_{ij} (\rho^n_{ij})^2\), (which differs from the Renyi-2 entropy by a minus sign and a logarithm) where \(\rho^n\) is the reduced density matrix of the \(n\)-th subsystem. The temperature-dependent estimate of the Renyi-2 entropy can then be obtained using the
finite_temp_analysis.pyscript.Valid for
instantaneousonly; ground-state Renyi-2 averaged over a single beta loop is not currently implemented. Requiresreplica_tricksto be enabled in order to obtained unbiased estimates.von_neumanntype: boolean.
Optional. Default: false.
Calculate the unnormalised von Neumann entropy and the trace of the reduced density matrix at the end of each beta loop. The normalised von Neumann entropy can be calculated from this using the
average_entropy.pyscript.Valid for
ground_stateonly; temperature-dependent von Neumann entropy is not currently implemented.writetype: boolean.
Optional. Default: false.
Print out the ground-state RDM to a file at the end of each beta loop. The file contains the trace of the RDM in the first line followed by elements of the upper triangle of the RDM labelled by their index.
Valid for
ground_stateonly.ref_projected_energytype: boolean.
Optional. Default: false.
Calculate the numerator and denominator for the projected energy as well as the total walker population for the reference row (or column) of the density matrix. Currently only available for read in systems.