QMeCha: quantum Monte Carlo package for fermions in embedding environments

TL;DR

QMeCha introduces an open-access quantum Monte Carlo package designed to treat fermions embedded in semi-quantum environments comprising classical charges and quantum Drude oscillators, enabling explicit modeling of dispersion, polarization, and electrostatics in large, mixed systems. The framework combines flexible fermionic wavefunctions (Pfaffian, AGP, Slater) with comprehensive Jastrow factors and dedicated drudon wavefunctions, optimized via Variational Monte Carlo and Diffusion Monte Carlo with correlated sampling and size-consistent improvements. The paper details the Hamiltonian formalism, embedding strategies (El-QDO), and a wide range of basis sets and pseudopotentials, together with modular, scalable code structure implemented in Fortran 2008 and demonstrated on sizable benchmarks. Three primary applications—vdW-rich macromolecules, electron-positron systems, and quantum embedding—show QMeCha’s potential as a reference tool for benchmark studies and for developing advanced multi-scale methods, with future prospects for backflow and force evaluations to broaden its utility.

Abstract

We present the first open access version of the QMeCha (Quantum MeCha) code, a quantum Monte Carlo (QMC) package developed to study many-body interactions between different types of quantum particles, with a modular and easy-to-expand structure. The present code has been built to solve the Hamiltonian of a system that can include nuclei and fermions of different mass and charge, e.g. electrons and positrons, embedded in an environment of classical charges and quantum Drude oscillators. To approximate the ground state of this many-particle operator, the code features different wavefunctions. For the fermionic particles, beyond the traditional Slater determinant, QMeCha also includes Geminal functions such as the Pfaffian, and presents different types of explicit correlation terms in the Jastrow factors. The classical point charges and quantum Drude oscillators, described through different variational ansätze, are used to model a molecular environment capable of explicitly describing dispersion, polarization, and electrostatic effects experienced by the nuclear and fermionic subsystem. To integrate these wavefunctions, efficient variational Monte Carlo and diffusion Monte Carlo protocols have been developed, together with a robust wavefunction optimization procedure that features correlated sampling. In conclusion, QMeCha is a massively parallel code introduced here to explore quantum correlation effects in mixed systems with thousands of fermions and bosonic particles, beyond what was previously accessible to other reference methods.

Figures (9)

• Figure 1: a. Schematic representation of a cluster of 30 water moleculesrak+19jcp from https://sites.uw.edu/wdbase/database-of-water-clusters/ in which the central water dimer is represented through the fermionic Hamiltonian (b) while the rest of the molecules are substituted by a model obtained from point charges ($Q_{H}$ and $Q_M$) and quantum Drude Oscillators (QDOs) (c). In b also a positron is displayed to indicate that the code can handle also different types of fermions. In c the quadratic potential between the QDO center and the drudon is represented as a spring. All the other interactions are described by Coulomb potentials.
• Figure 2: Schematic representation of the total wavefunction presented in the code (eq. \ref{['eq:tot_wf']}). The Pfaffian includes both Singlet (SG) and Triplet (TG) geminal correlations. Both SG and TG reduce to the Slater determinant (SD) with a particular set of constraints (Section \ref{['ssec:eleorposwf']}). These wavefunction can be used to describe the fermionic behavior of the separate populations of electrons and positrons (eq. \ref{['eq:elecpos_wvf2']}). The cusp functions in the Jastrow factor are used for all sets of particles in the Hamiltonian, nuclei, fermions, QDO centers and drudons.
• Figure 3: Schematic representation of the main structure of the QMeCha code.
• Figure 4: Weak scaling tests for the system of 30 water molecules all described as electronic systems. The original geometry is taken from ref. rak+19jcp (https://sites.uw.edu/wdbase/database-of-water-clusters/) and is shown in Fig. \ref{['fig:fig1']}a.
• Figure 5: a. Difference between the binding energies (BE) of the dimers computed with FN-DMC and DFT (PBE0+MBD) and the LNO-CCSD(T) calculations. For the LNO-CCSD(T) calculations we also report the estimated uncertainty that includes the basis set extrapolation and the error due to the orbital localization (see ref. Puleva2025 for further details), indicated with the different $\sigma,2\sigma,3\sigma$ curves. These are only used as a guide for the eye to underline the dimers for which there is the largest discrepancy between FN-DMC and LNO-CCSD(T). b. Relative error (in percentage) of the FN-DMC and DFT (PBE0+MBD) calculations defined with respect to the LNO-CCSD(T) ones as $| (\text{BE}-\text{BE}_{\text{LNO-CCSD(T)}} )/ \text{BE}_{\text{LNO-CCSD(T)}}|$. In this case the errors on the FN-DMC and PBE0+MBD calculations are obtained through the propagation of error that also includes the uncertainty in the LNO-CCSD(T) results. On the left we report the average relative absolute errors for FN-DMC and PBE0+MBD respectively in blue and red. In panels c, d, e, f we display the structures of the conformers F1I2, L2I2, SF1I2 and SF3I3 respectively, that show the larges discrepancies between FN-DMC and LNO-CCSD(T) results, as shown in panel a. Data is taken from refs. Puleva2025QUIDrep.