Challenges
Datasets
Workspaces
Papers
Discussions
Leaderboard
Log inSign up
Challenges
Datasets
Workspaces
Papers
Discussions
Leaderboard
Blog
Job Board
Q3AS

© 2026 Aqora Quantum S.A.S.

TermsPrivacyLegal Notice

The Journal of Chemical Physics

A staging Monte Carlo algorithm for sampling off-diagonal density matrix elements via open-chain path integrals

Alan Robledo, Mark E. Tuckerman

Published May 15, 2026• 0 citations
Publisher

Abstract

We introduce a computationally simple algorithm for sampling open-chain distributions within the framework of imaginary-time Feynman path integration. The present method is based on the staging algorithm introduced by Pollock and Ceperley [Phys. Rev. B 30, 2555 (1984)] originally developed for computing position-dependent observables. Here, we sample off-diagonal elements of the density matrix, formulated as a distribution describing a linear polymer-like chain of beads, each connected via nearest-neighbor springs to calculate momentum-dependent quantities. This is achieved using a Monte Carlo scheme that ensures efficient and unbiased sampling of all beads along the chain from the free-particle distribution via a staging transformation; we refer to this approach as staging open path integral Monte Carlo (OPIMC). The proposed algorithm is straightforward to implement, as it only involves sampling Gaussian distributions through a transformation defined by a set of recursion relations, followed by a standard Metropolis acceptance/rejection step. The staging OPIMC method accurately reproduces end-to-end and momentum distributions for quantum systems ranging from coupled harmonic oscillators to liquid water.

DOI
10.1063/5.0331949

Related papers

Exact ground state Monte Carlo method for Bosons without importance sampling

M. Rossi, Marco Nava, L. Reatto +1 more

Generally "exact" quantum Monte Carlo computations for the ground state of many bosons make use of importance sampling. The importance sampling is based either on a guiding function or on an initial variational wave function. Here we investigate the need of importance sampling in the case of path integral ground state (PIGS) Monte Carlo. PIGS is based on a discrete imaginary time evolution of an initial wave function with a nonzero overlap with the ground state, which gives rise to a discrete path which is sampled via a Metropolis-like algorithm. In principle the exact ground state is reached in the limit of an infinite imaginary time evolution, but actual computations are based on finite time evolutions and the question is whether such computations give unbiased exact results. We have studied bulk liquid and solid (4)He with PIGS by considering as initial wave function a constant, i.e., the ground state of an ideal Bose gas. This implies that the evolution toward the ground state is driven only by the imaginary time propagator, i.e., there is no importance sampling. For both phases we obtain results converging to those obtained by considering the best available variational wave function (the shadow wave function) as initial wave function. Moreover we obtain the same results even by considering wave functions with the wrong correlations, for instance, a wave function of a strongly localized Einstein crystal for the liquid phase. This convergence is true not only for diagonal properties such as the energy, the radial distribution function, and the static structure factor, but also for off-diagonal ones, such as the one-body density matrix. This robustness of PIGS can be traced back to the fact that the chosen initial wave function acts only at the beginning of the path without affecting the imaginary time propagator. From this analysis we conclude that zero temperature PIGS calculations can be as unbiased as those of finite temperature path integral Monte Carlo. On the other hand, a judicious choice of the initial wave function greatly improves the rate of convergence to the exact results.

arXiv:0907.4430•The Journal of Chemical Physics•Oct 19, 2009•62 citations

Algorithms for quantum molecular dynamics: from matrix product states to path integrals

Dmitri Iouchtchenko

This thesis describes several novel approaches in quantum molecular dynamics for obtaining properties of molecular systems in different regimes. We investigate ground state properties of chains of linear rotors with dipole–dipole interactions via the density matrix renormalization group (DMRG), by deriving the appropriate form of the interaction operator and implementing it in ITensor. This provides us with further evidence of a quantum phase transition in this system. We also improve the sampling of Gaussian mixture distributions for finite temperature path integral Monte Carlo (PIMC) of vibronic Hamiltonians. To do this, we replace random sampling by quasi-random sampling, and improve sampling distributions by optimizing their parameters. Finally, we introduce estimators and integrators for constrained free energy simulations in path integral molecular dynamics (PIMD). This method is applied to the study of a water dimer, for which we obtain a quantum potential of mean force.

•Jan 13, 2021•0 citations

Path-integral diffusion Monte Carlo: Calculation of observables of many-body systems in the ground state

Balázs Hetényi, Eran Rabani, B. J. Berne

We propose a new method to calculate ground state position dependent observables in quantum many-body systems. The method, which we call the path-integral diffusion Monte Carlo (PI-DMC) method, is essentially a combination of path-integral Monte Carlo (PIMC) and diffusion Monte Carlo (DMC) methods. The distribution resulting from a DMC simulation is further propagated in imaginary time by PIMC sampling. Tests of the new method for simple cases such as the harmonic oscillator, a double well, and a set of ten coupled harmonic oscillators show that the results for several observables converge rapidly to the exact result.

The Journal of Chemical Physics•Apr 1, 1999•21 citations

Iterative Monte Carlo formulation of real-time correlation functions

Cristian O. Baltaretu, Nancy Makri

We present an iterative Monte Carlo path integral methodology for evaluating thermally averaged real-time correlation functions. Standard path integral Monte Carlo methods are used to sample paths along the imaginary time contour. Propagation of the density matrix is performed iteratively on a grid composed of the end points of the sampled paths. Minimally oscillatory propagators are constructed using energy filtering techniques. A single propagation yields the values of the correlation function at all intermediate time points. Model calculations suggest that the method yields accurate results over several oscillation periods and the statistical error grows slowly with increasing propagation time.

The Journal of Chemical Physics•Oct 25, 2010•7 citations

Path-integral formulation of bosonic Markovian open quantum dynamics with Monte Carlo stochastic trajectories using the Glauber–Sudarshan P, Wigner, and Husimi Q functions and hybrids

Toma Yoneya, Kazuya Fujimoto, Yuki Kawaguchi

The Monte Carlo trajectory sampling of stochastic differential equations based on the quasiprobability distribution functions, such as the Glauber–Sudarshan P, Wigner, and Husimi Q functions, enables us to investigate bosonic open quantum many-body dynamics described by the Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) equation. In this method, the Monte Carlo samplings for the initial distribution and stochastic noises incorporate quantum fluctuations, and thus, we can go beyond the mean-field approximation. However, description using stochastic differential equations is possible only when the corresponding Fokker–Planck equation has a positive-semidefinite diffusion matrix. In this work, we analytically derive the stochastic differential equations for arbitrary Hamiltonian and jump operators based on the path-integral formula, independently of the derivation of the Fokker–Planck equation. In the course of the derivation, we formulate the path-integral representation of the GKSL equation by using the -ordered quasiprobability distribution function, which systematically describes the aforementioned quasiprobability distribution functions by changing the real parameter . The essential point of this derivation is that we employ the Hubbard–Stratonovich transformation in the path integral, and its application is not always feasible. We find that the feasible condition of the Hubbard–Stratonovich transformation is identical to the positive-semidefiniteness condition of the diffusion matrix in the Fokker–Planck equation. In the benchmark calculations, we confirm that the Monte Carlo simulations of the obtained stochastic differential equations well reproduce the exact dynamics of physical quantities and non-equal time correlation functions of numerically solvable models, including the Bose–Hubbard model. This work clarifies the applicability of the approximation and gives systematic and simplified procedures to obtain the stochastic differential equations to be numerically solved.

arXiv:2508.01991•Annals of Physics•Dec 16, 2025•0 citations

Path integral Monte Carlo in the angular momentum basis for a chain of planar rotors

Estêvão V. B. de Oliveira, Muhammad Shaeer Moeed, Pierre–Nicholas Roy

We introduce a path integral Monte Carlo approach that uses the angular momentum representation for the description of interacting rotor systems. Such a choice of representation allows the calculation of momentum properties without having to break the paths. The discrete nature of the momentum basis also allows the use of rejection-free Gibbs sampling techniques. To illustrate the method, we study the collective behavior of N confined planar rotors with dipole-dipole interactions, a system known to exhibit a quantum phase transition from a disordered to an ordered state at zero temperature. Ground state properties are obtained using the path integral ground state method. We propose a Bond-Hamiltonian decomposition for the high temperature density matrix factorization of the imaginary time propagator. We show that cluster-loop type moves are necessary to overcome ergodicity issues and to achieve efficient Markov chain updates. Ground state energies and angular momentum properties are computed and compared with density matrix renormalization group benchmark results. In particular, the derivative of the kinetic energy with respect to the interaction strength estimator is presented as a successful order parameter for the detection of the quantum phase transition.

arXiv:2506.14977•The Journal of Chemical Physics•Oct 9, 2025•2 citations
We introduce a computationally simple algorithm for sampling open-chain distributions within the framework of imaginary-time Feynman path integration. The present method is based on the staging algorithm introduced by Pollock and Ceperley [Phys. Rev. B 30, 2555 (1984)] originally developed for computing position-dependent observables. Here, we sample off-diagonal elements of the density matrix, formulated as a distribution describing a linear polymer-like chain of beads, each connected via nearest-neighbor springs to calculate momentum-dependent quantities. This is achieved using a Monte Carlo scheme that ensures efficient and unbiased sampling of all beads along the chain from the free-particle distribution via a staging transformation; we refer to this approach as staging open path integral Monte Carlo (OPIMC). The proposed algorithm is straightforward to implement, as it only involves sampling Gaussian distributions through a transformation defined by a set of recursion relations, followed by a standard Metropolis acceptance/rejection step. The staging OPIMC method accurately reproduces end-to-end and momentum distributions for quantum systems ranging from coupled harmonic oscillators to liquid water.
Generally "exact" quantum Monte Carlo computations for the ground state of many bosons make use of importance sampling. The importance sampling is based either on a guiding function or on an initial variational wave function. Here we investigate the need of importance sampling in the case of path integral ground state (PIGS) Monte Carlo. PIGS is based on a discrete imaginary time evolution of an initial wave function with a nonzero overlap with the ground state, which gives rise to a discrete path which is sampled via a Metropolis-like algorithm. In principle the exact ground state is reached in the limit of an infinite imaginary time evolution, but actual computations are based on finite time evolutions and the question is whether such computations give unbiased exact results. We have studied bulk liquid and solid (4)He with PIGS by considering as initial wave function a constant, i.e., the ground state of an ideal Bose gas. This implies that the evolution toward the ground state is driven only by the imaginary time propagator, i.e., there is no importance sampling. For both phases we obtain results converging to those obtained by considering the best available variational wave function (the shadow wave function) as initial wave function. Moreover we obtain the same results even by considering wave functions with the wrong correlations, for instance, a wave function of a strongly localized Einstein crystal for the liquid phase. This convergence is true not only for diagonal properties such as the energy, the radial distribution function, and the static structure factor, but also for off-diagonal ones, such as the one-body density matrix. This robustness of PIGS can be traced back to the fact that the chosen initial wave function acts only at the beginning of the path without affecting the imaginary time propagator. From this analysis we conclude that zero temperature PIGS calculations can be as unbiased as those of finite temperature path integral Monte Carlo. On the other hand, a judicious choice of the initial wave function greatly improves the rate of convergence to the exact results.
This thesis describes several novel approaches in quantum molecular dynamics for obtaining properties of molecular systems in different regimes. We investigate ground state properties of chains of linear rotors with dipole–dipole interactions via the density matrix renormalization group (DMRG), by deriving the appropriate form of the interaction operator and implementing it in ITensor. This provides us with further evidence of a quantum phase transition in this system. We also improve the sampling of Gaussian mixture distributions for finite temperature path integral Monte Carlo (PIMC) of vibronic Hamiltonians. To do this, we replace random sampling by quasi-random sampling, and improve sampling distributions by optimizing their parameters. Finally, we introduce estimators and integrators for constrained free energy simulations in path integral molecular dynamics (PIMD). This method is applied to the study of a water dimer, for which we obtain a quantum potential of mean force.
We propose a new method to calculate ground state position dependent observables in quantum many-body systems. The method, which we call the path-integral diffusion Monte Carlo (PI-DMC) method, is essentially a combination of path-integral Monte Carlo (PIMC) and diffusion Monte Carlo (DMC) methods. The distribution resulting from a DMC simulation is further propagated in imaginary time by PIMC sampling. Tests of the new method for simple cases such as the harmonic oscillator, a double well, and a set of ten coupled harmonic oscillators show that the results for several observables converge rapidly to the exact result.
We present an iterative Monte Carlo path integral methodology for evaluating thermally averaged real-time correlation functions. Standard path integral Monte Carlo methods are used to sample paths along the imaginary time contour. Propagation of the density matrix is performed iteratively on a grid composed of the end points of the sampled paths. Minimally oscillatory propagators are constructed using energy filtering techniques. A single propagation yields the values of the correlation function at all intermediate time points. Model calculations suggest that the method yields accurate results over several oscillation periods and the statistical error grows slowly with increasing propagation time.
The Monte Carlo trajectory sampling of stochastic differential equations based on the quasiprobability distribution functions, such as the Glauber–Sudarshan P, Wigner, and Husimi Q functions, enables us to investigate bosonic open quantum many-body dynamics described by the Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) equation. In this method, the Monte Carlo samplings for the initial distribution and stochastic noises incorporate quantum fluctuations, and thus, we can go beyond the mean-field approximation. However, description using stochastic differential equations is possible only when the corresponding Fokker–Planck equation has a positive-semidefinite diffusion matrix. In this work, we analytically derive the stochastic differential equations for arbitrary Hamiltonian and jump operators based on the path-integral formula, independently of the derivation of the Fokker–Planck equation. In the course of the derivation, we formulate the path-integral representation of the GKSL equation by using the -ordered quasiprobability distribution function, which systematically describes the aforementioned quasiprobability distribution functions by changing the real parameter . The essential point of this derivation is that we employ the Hubbard–Stratonovich transformation in the path integral, and its application is not always feasible. We find that the feasible condition of the Hubbard–Stratonovich transformation is identical to the positive-semidefiniteness condition of the diffusion matrix in the Fokker–Planck equation. In the benchmark calculations, we confirm that the Monte Carlo simulations of the obtained stochastic differential equations well reproduce the exact dynamics of physical quantities and non-equal time correlation functions of numerically solvable models, including the Bose–Hubbard model. This work clarifies the applicability of the approximation and gives systematic and simplified procedures to obtain the stochastic differential equations to be numerically solved.
We introduce a path integral Monte Carlo approach that uses the angular momentum representation for the description of interacting rotor systems. Such a choice of representation allows the calculation of momentum properties without having to break the paths. The discrete nature of the momentum basis also allows the use of rejection-free Gibbs sampling techniques. To illustrate the method, we study the collective behavior of N confined planar rotors with dipole-dipole interactions, a system known to exhibit a quantum phase transition from a disordered to an ordered state at zero temperature. Ground state properties are obtained using the path integral ground state method. We propose a Bond-Hamiltonian decomposition for the high temperature density matrix factorization of the imaginary time propagator. We show that cluster-loop type moves are necessary to overcome ergodicity issues and to achieve efficient Markov chain updates. Ground state energies and angular momentum properties are computed and compared with density matrix renormalization group benchmark results. In particular, the derivative of the kinetic energy with respect to the interaction strength estimator is presented as a successful order parameter for the detection of the quantum phase transition.