.. _dmqmc: Density Matrix Quantum Monte Carlo ================================== .. code-block:: lua 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_weights`` is set to ``true``. This can be passed directly to the ``weights`` option of a subsequent DMQMC calculation and/or manipulated inside the lua script. If ``find_weights`` is ``false``, only the qmc_state object is returned. ``dmqmc`` performs a density matrix quantum Monte Carlo (DMQMC) calculation on a system. Unlike :ref:`ccmc` and :ref:`fciqmc`, 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 ------- ``sys`` type: system object. Required. The system on which to perform the calculation. Must be created via a system function. ``qmc`` type: lua table. Required. Further options that are common to all implemented QMC algorithms. See :ref:`qmc_table`. ``dmqmc`` type: lua table. Optional. Further options to control the DMQMC algorithm. See :ref:`dmqmc_table`. ``ipdmqmc`` type: 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 :ref:`ipdmqmc_table`. ``operators`` type: lua table. Optional. Further options to select the operators for which expectation values are evaluated. See :ref:`operators_table`. ``rdm`` type: lua table. Optional. Further options to select which (if any) reduced density matrices and corresponding operators are to be evaluated. See :ref:`rdm_table`. ``restart`` type: lua table. Optional. Further options to control restarting the calculation from a previous calculation. See :ref:`restart_table`. ``reference`` type: lua table. Optional. Further options to select the reference state used. See :ref:`reference_table`. ``qmc_state`` type: 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_table: dmqmc options ------------- ``symmetric_bloch`` type: boolean. Optional. Default: true. Use the symmetrized form of the Bloch equation, :math:`\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, :math:`\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_tricks`` type: 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_temperature`` type: boolean. Optional. Default: false. Rescale tau so that the simulation runs in timesteps of :math:`\Delta\tau / T_F` where :math:`T_F` is the Fermi temperature. This is so results are at dimensionless inverse temperatures of :math:`\Theta^{-1} =T_F/T`. This option is only valid for systems with a well defined Fermi energy. ``all_sym_sectors`` type: 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_sectors`` type: 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_loops`` type: integer. Optional. Default: 100. The number of loops over the desired temperature range (each starting from :math:`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_beta`` type: 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_weights`` type: 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 :math:`N` electrons and more than :math:`2N` spin orbitals, this would be equal to :math:`N`. For a Heisenberg model with :math:`N` spins in the :math:`M_s=0` spin sector, this should be equal to :math:`N/2` (each pair of opposite spins flipped is one excitation). ``vary_weights`` type: 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 :math:`(W_{\gamma})^{\beta/\beta_{target}}` each iteration, where :math:`W_{\gamma}` is the final weight of excitation level :math:`\gamma` and :math:`\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_weights`` type: 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_start`` type: integer. Optional. Default: 0. The iteration number at which averaging of the excitation distribution begins, when using the **find_weights** option. ``symmetrize`` type: 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_level`` type: 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_shift`` type: 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_factor`` type: 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_table: ipdmqmc options --------------- ``target_beta`` type: 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 :math:`\tilde{\beta} = 1/\Theta = T_F/T`, where :math:`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_matrix`` type: string. Optional. Default: 'hartree_fock'. Possible values: 'free_electron', 'hartree_fock'. Initialisation of the density matrix at :math:`\tau=0`. 'free_electron' samples the free electron density matrix, i.e. :math:`\hat{\rho} = \sum_i e^{-\beta \sum_j \varepsilon_j \hat{n}_j} |D_i\rangle\langle D_i|`, where :math:`\varepsilon_j` is the single-particle eigenvalue and :math:`\hat{n}_j` the corresponding number operator. 'hartree_fock' samples a 'Hartree--Fock' density matrix defined by :math:`\hat{\rho} = \sum e^{-\beta H_{ii}} |D_i\rangle\langle D_i|`, where :math:`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_initialisation`` type: 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_check`` type: boolean. Optional. Default: false. When performing **grand_canonical_initialisation**, we check that :math:`H_{ii}` is not lower in energy than :math:`H_{00}`. If a lower energy :math:`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 :math:`H_{ii}`. .. warning:: It is recommended that the orbital single particle eigenvalues in the FCIDUMP are recalculated with the new reference. ``metropolis_attempts`` type: 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_picture`` type: boolean. Optional. Default: true. Use symmetric version of ip-dmqmc where now :math:`\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_occupations`` type: 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_table: operators options ----------------- ``renyi2`` type: boolean. Optional. Default: false. Calculate the Renyi-2 entropy of the entire system. Requires ``replica_tricks`` to be enabled. ``energy`` type: boolean. Optional. Default: false. Calculate the thermal expectation value of the Hamiltonian operator. ``energy2`` type: boolean. Optional. Default: false. Calculate the thermal expectation value of the Hamiltonian operator squared. Only available for the Heisenberg model. ``staggered_magnetisation`` type: 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_dist`` type: 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. ``correlation`` type: 2D vector of integers. Optional. Default: false. Calculate the spin-spin correlation function between the two specified lattice sites, :math:`i` and :math:`j`, which is defined as the thermal expectation value of: .. math:: \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_energy`` type: boolean Optional. Default: false Evaluate the bare Coulomb energy. Only available for the UEG. ``kinetic_energy`` type: boolean Optional. Default: false Evaluate the kinetic energy. Only available for the UEG. ``H0_energy`` type: boolean Optional. Default: false Evaluate the thermal expectation value of the zeroth order Hamiltonian where :math:`\hat{H} = \hat{H}^0 + \hat{V}`. See **initial_matrix** option. Only available when using the ip-dmqmc algorithm. ``HI_energy`` Evaluate the expectation value of the interaction picture Hamiltonian where .. math:: \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_dist`` type: float Optional. Default: 0.0 Evaluate the (spin averaged) momentum distribution in kspace, i.e., :math:`\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_factor`` type: float Optional. Default: 0.0 Evaluate the static structure factor: .. math:: 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 :math:`S(q)-1`, :math:`S_{\uparrow\uparrow}(q)+S_{\downarrow\downarrow}(q)-1` and :math:`S_{\uparrow\downarrow}(q)+S_{\downarrow\uparrow}(q)`, where :math:`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_table: rdm options ----------- Note that the use of RDMs is currently only available with the Heisenberg model. ``rdms`` type: 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: .. code-block:: lua rdms = { { 1, 2 } } specifies one RDM containing basis functions with indices 1 and 2, and .. code-block:: lua 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 ``instantaneous`` or ``ground_state`` must be enabled to set the desired mode of evaluating the RDM (but both options cannot be used together). ``instantaneous`` type: 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_state`` type: 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_size`` type: integer. Required if ``instantaneous`` is 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 :ref:`qmc_table`) and depends on the size of the subsystem compared to the full space. .. todo - should allow -ve numbers to specify the MB usage instead (as in the main spawned_state_size?) .. 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_start`` type: integer. Optional. Default: 0. Monte Carlo cycle from which the RDM is to be accumulated in each beta loop. Relevant only if ``ground_state`` is 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. ``concurrence`` type: 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.py`` script. Valid for ``ground_state`` only; temperature-dependent concurrence is not currently implemented. ``renyi2`` type: boolean. Optional. Default: false. Calculate the Renyi-2 entropy of each subsystem. More accurately, the quantity output to the data table is :math:`S^n_2 = \sum_{ij} (\rho^n_{ij})^2`, (which differs from the Renyi-2 entropy by a minus sign and a logarithm) where :math:`\rho^n` is the reduced density matrix of the :math:`n`-th subsystem. The temperature-dependent estimate of the Renyi-2 entropy can then be obtained using the ``finite_temp_analysis.py`` script. Valid for ``instantaneous`` only; ground-state Renyi-2 averaged over a single beta loop is not currently implemented. Requires ``replica_tricks`` to be enabled in order to obtained unbiased estimates. ``von_neumann`` type: 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.py`` script. Valid for ``ground_state`` only; temperature-dependent von Neumann entropy is not currently implemented. ``write`` type: 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_state`` only. ``ref_projected_energy`` type: 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.