Efficient Sampling of Thermal Averages of Interacting Quantum Particle Systems with Random Batches

An efficient sampling method, the pmmLang+RBM, is proposed to compute the quantum thermal average in the interacting quantum particle system. Benefiting from the random batch method (RBM), the pmmLang+RBM reduces the complexity due to the interaction forces per timestep from $O(NP^2)$ to $O(NP)$, where $N$ is the number of beads and $P$ is the number of particles. Although the RBM introduces a random perturbation of the interaction forces at each timestep, the long time effects of the random perturbations along the sampling process only result in a small bias in the empirical measure of the pmmLang+RBM from the target distribution, which also implies a small error in the thermal average calculation. We numerically study the convergence of the pmmLang+RBM, and quantitatively investigate the dependence of the error in computing the thermal average on the parameters including the batch size, the timestep, etc. We also propose an extension of the pmmLang+RBM, which is based on the splitting Monte Carlo method and is applicable when the interacting potential contains a singular part.


I. INTRODUCTION
Simulating complex chemical systems with quantum effects has been an appealing subject in computational physics and chemistry. In quantum systems, thermal properties are fully described by the canonical ensemble, and a considerable number of methods of calculating thermal averages are based on the path integral representation [1][2][3][4] , which reformulates the quantum system as a classical ring polymer system. In the past decades, the path integral Monte Carlo (PIMC) [5][6][7] and the path integral molecular dynamics (PIMD) [8][9][10] techniques have been developed and successfully applied to the calculation of quantum properties including reaction rates [11][12][13] , correlation functions 14,15 and quantum tunneling 16,17 .
In this paper we focus on long-range interacting quantum particle systems. The complexity to calculate the interaction forces of a P-particle system is O(P 2 ), and efficient computational methods are thus needed to simulate such big systems. Compared to short-range interactions which can be easily treated with cutoff 18 or data structures like the cell list 19,20 , the long-range nature of the interaction potential makes it more difficult to reduce the complexity due to the interaction forces.
So far a large variety of methods have been proposed to calculate the interaction forces for a particle system with electrostatic interactions. The most representative methods in this class are Ewald summation 21,22 for problems with periodic boundaries and fast multipole method (FMM) 23 for open systems, and the complexity of the interaction forces can be reduced to O(P log P) or even O(P). However, an obvious drawback of these methods is that they rely on the specific expression of the interaction potential, i.e., the Coulomb potential, and can hardly be applied to more general interacting systems.
Besides the methods which aim to calculate the interaction a) Electronic mail: abneryepku@pku.edu.cn b) Electronic mail: zhennan@bicmr.pku.edu.cn forces directly, there are also methods focusing on modifying the dynamics while maintaining the physical properties of interest. For example, the generalized Verlet algorithm 24 reduces the number of times of calculating interaction forces by introducing distance classes, and the random batch method (RBM) 25 simplifies the dynamics by random sampling at each iteration.
In this paper, the RBM is employed to reduce the complexity due to the interaction forces. The RBM is an efficient sampling method recently introduced for interacting particle systems, and has been applied to reduce the computational cost of certain examples in real time dynamics and ensemble average calculation 26,27 . In each timestep, the RBM randomly groups all particles into small batches, and the evolution of the whole system is replaced by the evolution in those batches respectively. Due to the random nature of the grouping in the RBM, the average force felt by each particle in the RBM is statistically the same as the force in the original dynamics. Since interaction forces are only calculated in the small batches, the complexity in a timestep is reduced from O(P 2 ) to O(P). Compared to other efficient sampling methods like the Ewald summation and the FMM, the RBM is easier to implement and more flexible to apply in complex interacting particle systems.
After including quantum effects, simulation of interacting particle systems meets additional challenges. In the context of the PIMD, suppose there are N beads in the reformulated ring polymer system. The complexity due to interaction forces in this reformulated system per timestep is O(NP 2 ), which is N times of the complexity in the classical case. In addition, the ring polymer system suffers from the stiffness of the ring polymer potential when the number of beads N is large. To resolve the stiffness, one can either precondition the dynamics by modifying the ring polymer mass matrix 28,29 or using the staging coordinates 30,31 , or employ other non-preconditioned methods 15,[32][33][34] .
In large interacting quantum systems, the compatibility between the RBM and the technique for resolving the stiffness must be taken into account. We choose to employ the pre-arXiv:2102.04688v4 [quant-ph] 4 Apr 2021 conditioned mass-modified Langevin dynamics (pmmLang) 29 rather than the staging coordinates to precondition the dynamics, since the pmmLang sticks to the use of the physical coordinates.
In this paper, we will show that the combination of the pmmLang and the RBM (denoted by the pmmLang+RBM) is a suitable sampling method in the large interacting quantum particle systems. The pmmLang+RBM costs only O(NP) complexity in a timestep to compute the interaction forces, and is able to obtain the accurate thermal average up to a small bias. Also, the pmmLang+RBM can be naturally combined with the splitting Monte Carlo method 26 to apply to systems with singular interacting potentials. Now we present our basic setup of this paper. In the Pparticle quantum system in R 3P with long-range interactions, we aim to calculate the thermal average of some observable operatorÂ. Suppose the potential of the system isV = V (q), with V (q) given by where q = (q 1 , · · · , q P ) ∈ R 3P is the position and momentum of the P particles, V (o) (q i ) is the external potential of q i in R 3 , and V (c) (q i − q j ) is the interacting potential between q i , q j . Assume each particle has mass m, then the Hamiltonian operator of the system iŝ H =p 2 2m +V (q), q, p ∈ R 3P (2) and the thermal average ofÂ is given by where β > 0 is the inverse temperature.
In the path integral representation, the quantum system (1) is reformulated as a classical ring polymer system with the potential function where β N = β /N and q = (q 1 , · · · , q N ) T ∈ R N×3P is the position of the ring polymer. In this way, the thermal average Â can be approximated by the ensemble average in the Boltzmann distribution π(q) ∝ e −β N U N (q) . The pmmLang+RBM we propose in this paper, which aims to sample the distribution π(q), can be implemented in the following steps: 1. Derive the pmmLang for the ring polymer system U N (q) to sample the Boltzmann distribution π(q).
2. Numerically integrate the sampling path of the pmm-Lang, where the interaction forces are efficiently computed by the RBM. Time averages of the weight functions are used to approximate the thermal average Â .
3. (optional) If the interacting potential V (c) (q) contains a singular part, combine the pmmLang+RBM with the splitting Monte Carlo method.
The pmmLang+RBM is simple and efficient to calculate the thermal average of the interacting quantum particle systems, and benefits from both the RBM and the pmmLang. Due to the RBM, the complexity of interaction forces in a timestep is reduced to O(NP), and the total complexity in a timestep is O(N log NP), which has satisfactory scaling properties for both N and P. When the interaction potential contains a singular part (e.g., the Lennard-Jones and the Morse potential), the pmmLang+RBM method can be extended by the use of the splitting Monte Carlo method, which lifts the constraint of using extremely small time steps. The paper is organized as follows. In Section II, we introduce the the PIMD for interacting particle systems and the difficulties in numerical simulation, then derive the pmm-Lang. In Section III, we introduce the random batch method and discuss of the error analysis of the pmmLang+RBM. In Section IV, we introduce the splitting Monte Carlo method and its combination with the pmmLang+RBM, the pmm-Lang+RBM+split. In Section V, we present the numerical results of the pmmLang+RBM and the pmmLang+RBM+split and report the error in the the calculation of thermal averages.

A. Ring polymer representation and Langevin sampling
In the path integral representation, the thermal average (3) is approximated as where β N = β /N, and is the position of the ring polymer. In the coordinate notation q i k , the subscript k indicates the bead index in the ring polymer representation, and the superscript i indicates the particle index in the physical system. Besides, U N (q) is the total potential of the ring polymer system defined in (4), is the partition function, and W N (q) is the weight corresponding to the observable operatorÂ. Standard theory of the path integral 5 shows that the approximation (5) is exact as the number of beads N → ∞. Define the Boltzmann distribution as then (5) can be equivalently written as i.e., Â is approximated as the ensemble average of the weight function W N (q) in the distribution π(q).
In this paper, the observable operators of interest are the position-dependent operatorÂ = A(q) and the kinetic energy operatorÂ =p 2 /(2m), where A(q) is an analytic function in R 3P . For the position-dependent operator, the weight W N (q) is simply For the kinetic energy operator, the weight W N (q) is chosen as the virial estimator 35 is the center of the ring polymer.
To compute the ensemble average W N (q) π , we employ the path integral molecular dynamics (PIMD) technique to sample the distribution π(q). By introducing an auxiliary momentum variable p ∈ R N×3P , the PIMD couples the the Hamiltonian dynamics and a thermostat scheme to preserve the invariant distribution. In practice, the widely used thermostats include the Andersen thermostat 36 , the Nosé-Hoover thermostat 37 , the Langevin thermostat, etc. In this paper we focus on the last one to develop an efficient sampling algorithm with random batches, and our underlying dynamics is thus the second-order Langevin dynamics.
For heuristic purposes, we present in the following the PIMD with the Langevin thermostat. Introduce the positive definite mass matrix M ∈ R N×N (maybe different from the physical mass m) and extend U N (q) to the Hamiltonian where ·, · F is the Frobenius inner product in R N×3P . By adding damping and diffusion terms in the Hamiltonian dynamics of (12), we obtain the second-order Langevin dynamics in the PIMD, where γ > 0 is the friction constant and B is standard Brownian motion in R N×3P .
The invariant distribution of (13) is whose marginal distribution in q is exactly π(q) as in (7). From the ergodicity of the Langevin dynamics, W N (q) π can be computed by the infinite time average of W N (q(t)) with the classical trajectory (q(t), p(t)) propogated by (13), i.e., In practice, one employs the discrete time trajectory q( j∆t) numerically solved by (13) to approximate the time integral in (15), i.e., where ∆t is the timestep. In the large interacting particle system, the most costly part of the Langevin dynamics (13) is computing the gradient ∇U N (q), which involves heavy calculation of the interaction forces. As we shall elaborate in the next part, the complexity due to interaction forces is O(NP 2 ).

B. Simulation bottleneck in Langevin sampling
In this paper we are mainly concerned with the efficiency of calculating the thermal average (3) for a large quantum interacting particle system, where the numerical challenges mainly originate from two aspects of reasons. On the one hand, given a quantum interacting particle system, the number of beads N in the path integral representation needs to be sufficiently large to approximate the thermal average Â . On the other hand, in order to accurately calculate a wide class of physical quantities, it is often desirable to simulate in large size systems, for example, the liquid water, since the number of particles P should be large enough to incorporate the correct scientific phenomenon.
When N and P are large, one of the major difficulties in the numerical simulation of (13) is that the computational cost per timestep is extremely heavy. To perform a complexity analysis, we write the total ∇U N (q) as where L ∈ R N×N is the second order difference matrix The gradient ∇U N (q) is thus given by where each ∇V (q k ) involves the full interaction forces in the system (1). In fact, and the interaction force felt by the i-th particle q i k is Therefore, the complexity of calculating each ∇V (q k ) is O(P 2 ), and the complexity of calculating the full gradient It is worth noting that the O(P 2 ) computational cost per timestep is the common difficulty in the simulation of interacting particle systems, regardless of the thermostat scheme. As long as one uses molecular dynamics approaches, heavy calculation of the interaction forces is inevitable. In the PIMD where an N-bead ring polymer is used, there are N duplicates of the original interacting system, and the computational cost becomes O(NP 2 ).
Except for the interaction forces ∇V (c) (q i k −q j k ), calculation of the weight W N (q) can also have O(NP 2 ) time complexity in a single timestep. For the position-dependent operatorÂ = A(q), if the analytic function A(q) is given by then the complexity to compute the weight function (9) is O(NP 2 ). Also, the virial estimator (10) for the kinetic energy involves the full gradient {∇V (q k )} N k=1 , thus the complexity is O(NP 2 ).
The O(NP 2 ) complexity in the calculation of the interaction forces and the weight is a huge impediment on efficient sampling, especially for large P. Hence we seek effective means to accelerate the simulation of the PIMD Langevin dynamics (13) by reducing the computational cost per timestep to O(NP). In the classical case, the Ewald summation 21 and FMM 23 have been successfully applied to reduce the complexity of interaction forces, but they are complicated to implement, cumbersome in the ring polymer representation, and only efficacious for specific interacting potentials, for example, the Coulomb potential. Therefore, we aim to develop a method to reduce the O(NP 2 ) complexity, which is easy-touse, consistent with the PIMD framework, and applicable for any interaction potential V (c) (q). For this reason, we propose a sampling algorithm motivated by the recently proposed random batch method (RBM) 25 .
In the RBM, the group of P particles is randomly divided into small batches, where the division is chosen independently in different timesteps. During one timestep, the P particles are restricted to interact within their own batches, and in this way the total complexity of interaction forces is reduced from O(NP 2 ) to O(NP). A more detailed description of incorporating the RBM in Langevin sampling is presented in Section III.

C. Choice of the preconditioning method
Before continuing the discussion on the RBM, we take a short detour to address another ubiquitous numerical issue in simulating the PIMD when N is large. The ring polymer potential U N (q) leads to a stiff term in the Langevin dynamics (13), which prevents the use of large time steps in integrating the sampling path. To get a glimpse of the stiffness, note that the condition number of the matrix L ∈ R N×N defined in (18) is hence if the mass matrix is simply chosen as the timestep ∆t needs to be small as O(1/N) to integrate the highest oscillation mode of (13), which means the Langevin dynamics (13) shows stiffness when N is large. The numerical stiffness can be resolved by certain preconditioning methods, which is often realized by introducing a proper change of coordinates or choosing a proper mass matrix. For example, the staging coordinates transformation 30 given bỹ has shown to be a powerful technique in relaxing the use of time steps. However, the staging coordinates is not compatible with the RBM unless one repeatedly uses the staging coordinates transformation and its inverse (see appendix A for a detailed discussion), which are inconvenient for large systems. In fact, a preconditioning method which can be directly implemented in the physical coordinates is preferable to be combined with the RBM. For this reason, we adopt the preconditioned mass-modified Langevin dynamics (pmmLang) 29 as the preconditioning method, which we present in the following. In the Langevin dynamics (13), we choose the mass matrix in such a specific form where the regularization parameter α > 0 is introduced to make sure L α is positive definite. Let be the velocity, and define the modified potential then (13) becomes The obtained Langevin dynamics (29) is exactly the pmm-Lang we aim to derive. Compared to the staging coordinates, the pmmLang can be conveniently implemented in the physical coordinates and is naturally compatible with the RBM. A detailed derivation of the pmmLang and the discussion on its implementation is given in Appendix A. Given the regularity parameter α > 0, the invariant distribution of the pmmLang (29) is (30) Since the right hand sides of (14) and (30) are exactly the same except for the velocity transformation v = (L α ) −1 p, the pmmLang (29) produces the correct invariant distribution (L α ) −1 π(q) in the marginal of q. If we could neglect the ∇U α (q) term in (29), the Hamiltonian part of the pmmLang is simply thus there is no stiffness in (29). With further analysis, the pmmLang (29) can be shown to be a successful preconditioning method 29 in the PIMD.
It is worth pointing out that the pmmLang has also been studied in some recent work ? ? , where the authors obtained a quantitative analysis of the convergence rate with its dependence on the regularization constant α. A qualitative corollary of their result is that α being close to zero or too large both slowdowns the simulation efficiency of the pmmLang. To our point of view, if α is close to zero, the entries of the matrix (L α ) −1 in (29) will be large as O(1/α), hence the magnitude of (L α ) −1 ∇U N (q) restricts the timestep for numerical stability. For mathematical well-posedness, α should be chosen to ensure that the modified potential U α (q) is confined, i.e., lim q→∞ U α (q) = +∞.
While it is not a trivial task to choose the optimal parameter α in the pmmLang, there are other non-preconditioned numerical integrators for the PIMD Langevin dynamics which are parameter free 15,34 . Nevertheless, how to choose α is not the focus of this work, and we employ the pmmLang as a preconditioning method in the PIMD under certain simplified conditions given below. We assume the external potential V (o) (q) in the quantum system (1) is harmonic, i.e., bluefor some constant α 0 > 0 we have We choose regularization parameter α > 0 to coincide with the constant α 0 in (33), thus the modified potential U α (q) is simply the sum of all interacting potentials, and the pmmLang (29) can be equivalently written as For simplicity, the incorporation of the RBM will be discussed only for (35), under the assumption of (33). Clearly, the RBM can be applied without (33), and there will be additional terms in U α (q) except the interacting potentials in (34). These additional terms contributes only O(NP) complexity in a single timestep and does not change the major difficulties, since the O(NP 2 ) complexity is only due to the interaction forces ∇V (c) (q i k − q j k ).

III. RANDOM BATCH METHOD IN LANGEVIN SAMPLING
A. RBM for the pmmLang In the simulation of large interacting particle systems, high computational cost per timestep has always been an impediment to efficient sampling. With the setup in the previous sections, to integrate the sampling trajectory for a single timestep, it requires a cost of O(NP 2 ) complexity to compute all the interaction forces V (c) (q i − q j ). Since the calculation of the interaction forces is rather expensive when P is large, reducing the complexity due to interaction forces is crucial to attain high efficiency in the PIMD simulation of interacting particle systems.
The recently proposed random batch method (RBM) provides a simple approach to fulfill such a reduction of cost in force evaluation. The RBM avoids calculation of the full interaction by randomly dividing the group of P particles into small-size batches, and only allowing interactions within each small batch with adjusted interaction strength. Since for each timestep, a new random division is conducted, after sufficiently many times, each particle will not only has an equal probability to interact with any other particles in the group, but also will have plenty of random interactions with the rest of the particles. This is why we expect the RBM to produce the correct statistical properties of the interacting system in the long time simulation. One does not necessarily need to trace the interactions between a typical particle with every other particle to have an accurate prediction of the statistical properties, but rather, the interactions with other particles may vastly cancel each other and result in an average field to the typical particle. For heuristic purpose, we present the RBM for the pmmLang (35) in the following.
During each timestep, we first randomly divide the group of P particles into n batches C 1 , · · · , C q , where each batch is of size p and n = P/p. The batch size p should be far less than P to avoid massive calculation of interaction forces, and be at least two to capture the pairwise interactions ∇V (c) (q i − q j ). Then we approximate the full interaction of the P-particle sys-tem with the interaction forces within each small batch C l , and adjust the interaction strength to ensure statistical consistency. To be specific, the interaction force felt by the i-th particle With the use of RBM, each q i is assigned to a batch C l for some l = 1, · · · , q, and hence the approximation of the interaction force (36) within C l is Recall that the division C 1 , · · · , C q is randomly generated for each timestep, thus for fixed q i , the batch C l which contains q i is also random. Therefore, (37) can be interpreted as a (p−1)term sample approximation of a (P−1)-term summation (36). Also, the approximation (37) is unbiased, i.e., (38) In this way, the random-batch approximated pmmLang within the batch C l is given by and the RBM for the pmmLang (denoted by the pmm-Lang+RBM) in a single timestep ∆t just evolves (39) for all batches C l , l = 1, · · · , n, as shown in Algorithm 1. In consecutive timesteps, the previous divisions are discarded and new sets of divisions are randomly chosen. Hence a fixed particle q i ∈ R N×3 can be assigned to different batches C l in different timesteps, and only interacts with the particles in the current batch.
In Algorithm 1, to generate a random division C 1 , · · · , C q , one may use the Fisher-Yates shuffle 38 with O(P) complexity to obtain a random permutation of 1, · · · , P and then divide it into n batches in order. Note that there are p(p−1) 2 pairs of interactions within each batch C l , the complexity of Algorithm 35 due to interaction forces is hence the batch size p should be small to attain high efficiency in the pmmLang+RBM. In particular, if p is chosen as small integers such as 2 and 4, the complexity due to interaction forces is only O(NP). In this way, we can employ the pmm-Lang+RBM (Algorithm 1), which is much more efficient than the original pmmLang (35), to compute the ensemble average W N (q) π .
Finally, we point out that the RBM has an analogy to the stochastic gradient descent (SGD) [39][40][41][42] , the stochastic optimization method widely used in machine learning. In machine learning problems, the size of the data set can be extremely large, and the SGD approximates the full gradient with the average calculated from samples in a randomly chosen small batch. Although the RBM and the SGD shares the idea of full gradient approximation, the high dimensional natures of their underlying problems are different. In an optimization problem where the SGD applies, the dimension of the parameter space does not increases with respect to the size of the data set. However, in an interacting particle system where the RBM applies, the dimension of the coordinate space growth linearly with the number of particles P.
Algorithm 1: RBM for the pmmLang (35) in a timestep ∆t Randomly divide the P particles into n batches C 1 , · · · , C n of size p, where n = P/p. for l = 1, · · · , n do Evolve the pmmLang within the batch C l in a timestep ∆t:

B. Error analysis in PIMD with the RBM
In this section, we aim to investigate the error in calculating the ensemble average W N (q) π due to the use of the RBM. By the ergodicity of the pmmLang (35), W N (q) π can be accurately computed from the time average as the number of samples J → ∞, where {q(t)} t 0 is the pmm-Lang trajectory evolved by (35) and ∆t is the timestep specified in Algorithm 1. Correspondingly, the pmmLang+RBM uses the following time average to estimate W N (q) π , where {q(t)} t 0 is the pmm-Lang+RBM trajectory evolved by Algorithm 1. Therefore, the error analysis in the RBM is actually asking whether the time average (42) is a good estimator of W N (q) π for sufficiently large J.
In general, exploring the effects of the consecutive uses of random-batch approximations (37) along the sampling process is difficult. To our knowledge, the strong error in a contracting dynamics 25 and the weak error in a short period 43 have been rigorously justified, whereas none of those theoretical results apply to the RBM in the long time sampling process.
If the adoption of the RBM results in a small strong error in the long time simulation, i.e., the pmmLang+RBM trajectoryq(t) almost coincides with the pmmLang trajectory q(t) all the time, we can certainly deduce that the RBM also leads to a small error in the ensemble average calculation. However, this is not true. As we shall show in the following, it is likely that the pmmLang+RBM trajectoryq(t) drifts significantly apart from q(t) in a short time, but stillq(t) provides a faithful and accurate approximation of W N (q) π in the long time simulation. We shall further rationalize the error induced by the use of the RBM in the PIMD in the following.

Strong error analysis
The strong error characterizes the deviation of the numerically computed trajectory from the exact one in form of the mean squared error 44 . The strong error of the pmm-Lang+RBM is then defined as where we assume the pmmLang+RBM trajectoryq(t) and the pmmLang trajectory q(t) are driven by the same Brownian motion B(t) and the same initial state. In addition, the expectation is taken over all possible Brownian motions and choices of random batches in the time interval [0,t]. When (35)(39) are integrated exactly, the deviation ofq(t) from q(t) is solely due to the use of the random batches. According to Algorithm 1, the pmmLang+RBM trajectorỹ q(t) coincides with q(t) when the batch size p = P or in the limit ∆t → 0. In fact, the batch size p = P implies the interaction forces are calculated accurately, while in the limit ∆t → 0 each particle is supposed to be driven by the average effects. Nevertheless, such parameters are impractical to apply in the pmmLang+RBM, because choosing either p = P or ∆t → 0 greatly increases the total computation cost and contradicts with our original intention to use the RBM. To be specific, the cost due to interaction forces of the pmmLang+RBM in the which is extremely large when p = P or ∆t → 0. When the parameters p, ∆t are chosen such that the pmm-Lang+RBM becomes an efficient sampling method, i.e., the batch size p is small and the timestep ∆t is relatively large, we cannot expect the pmmLang+RBM to produce the accurate pmmLang trajectory q(t). In contrast, the strong error e(t) shortly grows large in a short time period , as we show in the following example.
In the Coulomb interacting system (the potential function is defined in (68) in Section V), we show in Figure 1 how the pmmLang+RBM trajectoryq(t) deviates from q(t) and the strong error e(t) grows large along the sampling process. It can be seen from Figure 1 that the pmmLang+RBM trajectorỹ q(t) drifts apart from q(t) at about t ≈ 5, and the strong error e(t) grows large at t ≈ 1, no matter how small the timestep ∆t is. Therefore, there will always be large deviation in the pmmLang+RBM trajectoryq(t) from the pmmLang trajectory q(t), unless one chooses p = P or ∆t → 0. As we shall show in the next, requiring a small path-wise error is totally unnecessary in thermal average calculation.

Weak error analysis
The weak error measures the difference of the approximate and the exact values of the ensemble average 44 . In the PIMD, the weight function of interest is W N (q), which takes different forms for different observables. To be specific, the weak error of the pmmLang+RBM is defined as where the expectation is taken over all possible Brownian motions and choices of random batches, where the observable operator is the kinetic energy, and the corresponding weight function of interest is the virial estimator (10). To give a first impression of the weak error, we plot in Figure 2 the expectations EW N (q(t)) and EW N (q(t)) with various timesteps in the Coulomb interacting system. Unlike the strong error (43), the curve of EW N (q(t)) associated with the pmmLang+RBM is very close to the curve of EW N (q(t)) associated with the pmmLang, and the weak error (44) remains small all along the sampling process. It can be seen from Figure 2 that the ring polymer system quickly goes into equilibrium about t ≈ 2.5, for both the pmmLang and the pmmLang+RBM dynamics, hence the introduction of the batch force approximation (37) does not influence the convergence mechanism of the original pmmLang (35). Also, the stochastic error in computing W N (q(t)) due to the use of the force approximation (37) cancels over all possible choices of random batches, which finally results in a small bias on the average EW N (q(t)), even if the batch size p is small and the timestep ∆t is relatively large. However, in the long time simulation it's impossible to use EW N (q(t)) or EW N (q(t)) to estimate W N (q) π , since the calculation of the expectations in (44) requires a large number of sampling trajectories. In practice, it is more feasible to use the time averages (41)(42) to compute W N (q) π , which are ob-tained by numerically integrating one sampling path, respectively.
If the numerical error is neglected, we know that due to the ergodicity of the Langevin dynamics, the time average (41) generated by the pmmLang converges to the correct limit as the number of samples J → ∞, i.e.
We can actually view such a convergence from the perspective of the relative entropy. Recall that we aim to sample the Boltzmann distribution π(q) as in (7) and for any distribution f (q), we can use the relative entropy 45,46 to measure how much the distribution f deviates from the target distribution π. The relative entropy is always nonnegative, and it is zero only when the two distributions are identical. Note that we can define a family of the empirical distributions from the samples generated by the pmmLang, which are given by Thus, the convergence (45) can be interpreted as The diminishing of the relative entropy manifests the time irreversibility and the ergodicity of the sampling trajectory. That is, the thermostat effect of the Langevin dynamics brings in dissipation to the relative entropy with respect to the invariant measure, such that the empirical measure µ J (q) converges to the target Boltzmann distribution π(q) in the weak sense. However, when the random batches are used along a sampling path, the empirical distributions are repeated perturbed, mostly likely away from the Boltzmann distribution. In fact, when the random batches are used for each time, although the force can be viewed as an unbiased approximation from the total interaction force, the sampling path deviates from the Langevin dynamics without the use of random batches with a deterministic bias for each specific choice of the random divisions. Let us denote the empirical measure of the pmm-Lang+RBM byμ J (q), which is defined bỹ Because the random divisions are chosen independently, we cannot expectμ J (q) to converge as J → ∞. However, as we shall demonstrate with the following numerical example, the accumulation of the bias for each use of the random batches does not add up to an unbounded error. On the contrary, the perturbations of the interaction forces along the sampling path only results in a small bias inμ J (q) from the target distribution π(q).
We present a numerical experiment to observe the difference betweenμ J (q) and π(q) for sufficiently large J. Note that π(q) is defined over the high-dimensional coordinate space R N×3P , it's intractable to directly compute D(μ J ||π), hence we instead numerically simulate the distributions of certain observables associated withμ J (q) and π(q). In this example, we choose the observable operator to be the position q 1 1 ∈ R 3 and the pairwise distance |q 1 1 −q 2 1 | ∈ R. Recall that the superscript denotes the index of particles, while the subscript denotes the index of beads. We plot the relative entropy betweenμ J (q) and π(q) in the observables along the sampling process in Figure 3. As a comparison, the relative entropy for µ J (q), the empirical measure of the pmmLang, is also plotted.
We observe from Figure 3 that the empirical measure of the pmmLang+RBM converges to π(q) as fast as the pmmLang when the sampling time T < 10 4 . However, as T > 10 4 , there is a significant slowdown in the convergence of the pmm-Lang+RBM, especially in the bottom figure where the observable is the pairwise distance |q 1 1 − q 2 1 | ∈ R. Finally, there is a certain bias of the empirical measureμ J (q) from the target distribution π(q). The slowndown in the bottom panel of Figure 3 is more noticeable than the top one, because the RBM modifies the interaction forces and directly impacts the calculation of the pairwise distance. Meanwhile, the RBM has relatively less influence on the marginal distribution of a single particle.
The convergence ofμ J (q) to π(q) also implies the convergence of the time average (42) to W N (q) π , which can be verified in the following example. We plot the time averages (41)(42) of the pmmLang and the pmmLang+RBM in Figure 4, where the observable operator is the kinetic energy. It can be seen from Figure 4 that although (42) deviates from (41) at the first phase of the sampling process, (42) gradually converges to a fixed limit as time evolves. Furthermore, the bias of (41) from Â is small and diminishes if we shrink the timestep.
To understand the the approximation property of the pmm-Lang+RBM, we present a heuristic explanation on the decay of the relative entropy D(μ J ||π). Our key observation lies in the fact that the use of the random batches does not affect the thermostat part of the Langevin dynamics. And thus, the pmmLang+RBM can be viewed as an integrated process, where the thermostat part is constantly driving the empirical distribution approaching the target measure and the use of random batches gives rise to a sequence of random perturbations which presents the convergences of the empirical distribution. The numerical results strong suggests that the thermostat effect is dominating the random perturbations, and as a consequence, although the empirical distributionμ J (q) does not converge, it stays in a relatively small vicinity of the target distribution π(q) for J sufficiently large. To further analyze the dominance of the thermostat in such an integrated process, we continue the discussion in two distinct phases of sampling.
In the first phase of sampling, the samples generated by the Langevin dynamics are statistically dependent on the choice of the initial state, and within a fairly short time, the sampling path becomes uncorrelated with the initial state. This period roughly corresponds to the upper plot in Figure 3. During this phase of sampling, we observe that the time average of pmm-Lang+RBM is actually not very close to that of pmmLang, but the trends of both running averages are similar: after a short time, both of them start to fluctuate within a same small neighborhood. We remark that in practice, the samples from the first phase are "burned-in" anyway due to their low qualify, but it is crucial that the consecutive random perturbations are suppressed by the thermostatting mechanism of the Langevin dynamics such that the running average of the pmmLang+RBM is able to approach a small vicinity of the true value by the end of the first phase. In probability language, it means the difference between the probability distribution function of the stochastic process pmmLang+RBM and the target distribution significantly reduces within a short time period, after which the pmmLang+RBM is able to produce the correct samples from the target distribution.
In the second phase of sampling, the Langevin dynamics effectively produces a vast amount of samples of the target distribution such that the empirical measure is converging to the target Boltzmann distribution π(q), which is manifested by the decaying of the relative entropy. When the random batches are used, the generated position samplesq( j∆t) from pmm-Lang+RBM are not, unfortunately, unbiased samples from π(q). Whereas, the thermostatting effect that dissipates the relative entropy with respect to π(q) still cause that the difference between the empirical measureμ J (q) and the Boltzmann distribution π(q) diminishes in time. This argument is further confirmed by the numerical tests as shown in Figure 3, where we observe that with relative entropy of D(μ J ||π) from pmm-Lang+RBM decays exponentially in time, although the con-vergence speed may be slightly reduced comparing with the decaying of D(µ J ||π). This implies the large number of small biases induced by the use of the RBM in long time sampling only add to a small bias in the empirical measure. Furthermore, we conclude that the use of random batches only lead to small random perturbations to the sample values of the observables, and thus by the law of large numbers, we expect a small bias in the thermal average calculations.
We summarize the weak error analysis for pmmLang+RBM in the following mathematical statement: there exists a small constant C p,∆t depending only on the batch size p, and the timestep ∆t such that for J 1, In other words, the fluctuation-dissipation relationship for the pmmLang+RBM can be established in an approximate sense, when the RBM mechanism is incorporated. The rigorous justification of this statement is not yet complete and beyond the scope of this paper.

C. Efficient calculation of the weight function
In Section II, we have mentioned that the calculation of the weight W N (q) also has O(NP 2 ) complexity, if the observable operator is the kinetic energy or in the form of a summation over pairwise contributions. Using again the idea of random batches, we can make an unbiased approximation of W N (q) and reduce the complexity of calculating W N (q) from O(NP 2 ) to O(NP).
For the kinetic energyÂ =p 2 /(2m), rewrite the virial estimator (10) as where the modified potential U α (q) defined in (34) is the sum of all interacting potentials V (c) (q i k − q j k ). We can directly obtain an unbiased approximation of ∇U α (q) (and thus W N (q)) from the random-batch approximation (37). For the positiondependent operatorÂ = A(q) with A(q) given in the form we randomly pick a batch C of size p from the group of P particles (a full division is not required now) and approximate A(q) as which is unbiased. Hence the weight W N (q) can be efficiently computed by Finally, we emphasize that such approximation of the weight W N (q) does not change the original dynamics.

IV. SPLITTING MONTE CARLO METHOD
In the previous section, we have shown that the pmm-Lang+RBM is an efficient approximate integrator of the pmm-Lang (35). However, when the interaction potential V (c) (q) is singular (e.g., Lennard-Jones or Morse potential), extremely small timesteps ∆t are needed to integrate the dynamics near the singular kernel, no matter if the RBM is used.
In the simulation of the pmmLang, the singular kernel of V (c) (q) limits the timestep when the collision occurs, i.e., there are two particles q i , q j close enough to each other. Although the collision is a rare event, the requirement for small timesteps largely increases the computational cost. To overcome the singularity of the interaction potential, the pmm-Lang (35) should be modified to avoid making use of gradients of singular potentials, be compatible with the random batch method, and do not change the invariant distribution (30) of the pmmLang (35).
In this section we introduce a modified version of the pmm-Lang, which is based on the splitting Monte Carlo method 26 . Briefly speaking, we split the potential U α (q) into the smooth part U 1 (q) and the singular part U 2 (q). At each timestep, a proposal point q * is generated by evolving the pmmLang driven by U 1 (q), and is accepted or rejected according to the potential difference U 2 (q * ) −U 2 (q). Therefore, the splitting Monte Carlo method can be seen as a variant of the Metropolis-adjusted Langevin algorithm (MALA) 47 .
To illustrate the splitting Monte Carlo method, we split the interaction potential V (c) (q) into where V (c) 2 (q) is short-ranged and captures the singular part of V (c) (q), and thus V (c) then U α (q) is split into the sum of the smooth potential U 1 (q) and the singular potential U 2 (q). On the one hand, the pmm-Lang driven by the smooth U 1 (q) is whose invariant distribution is

Algorithm 3: RBM with splitting Monte Carlo for the pmmLang (35) in a timestep ∆t
Randomly divide the P particles into n batches C 1 , · · · , C n of size p, where n = P/p. for l = 1, · · · , n do Evolve the pmmLang within the batch C l in a timestep ∆t: On the other hand, when the proposal q * is accepted with probability a(q, q * ) = min{1, e −β N (U 2 (q * )−U 2 (q)) }, the Metropolis algorithm for the singular potential U 2 (q) samples its distribution Therefore, the key to assure that the coupled dynamics of the pmmLang (58) driven by U 1 (q) and the Metropolis algorithm for U 2 (q) is the obvious relation Since the pmmLang is the second-order Langevin dynamics, special care should be taken to preserve the detailed balance. According to 48 , the detailed balance of the second-order Langevin dynamics is given by where T ∆t ((q, v), (q , v )) is the transition probability density from (q, v) to (q , v ) in a timestep ∆t. The velocities v, v appearing in the right hand side of (63) are flipped as −v, −v , thus in the Metropolis algorithm for U 2 (q), the velocity v should also be flipped if the proposal q * is rejected. To sum up, the splitting Monte Carlo method for the pmmLang (denoted by the pmmLang+split) is given in Algorithm 2. The correctness of Algorithm 2 is guaranteed by the following theorem.
Theorem 1 Let T ∆t ((q, v), (q , v )) be the transition probability density of the splitting Monte Carlo method (Algorithm 2) in a timestep ∆t, which satisfies then the detailed balance holds The detailed balance (65) implies π(q, v) is the invariant distribution of the splitting Monte Carlo method (Algorithm 2), In the large interacting particle system, the computational cost of Algorithm 2 within a single timestep originates from   1 (q i k − q j k ) and the potential difference U 2 (q * ) −U 2 (q). Note that U 2 (q) defined in (56)(57) is the sum of all short-ranged interacting potentials V (c) 2 (q i k − q j k ), thus can be efficiently calculated by cutoff or using data structures such as the cell list 19 . In this paper we employ the cutoff method to compute U 2 (q), i.e., the particle pair (q i k , q j k ) is counted in the summation (56)(57) only when their distance |q i k − q j k | is less than the cutoff distance. Therefore, for large interacting particle systems, the majority of the computational cost is still the calculation of the interaction forces ∇V . Now we employ the RBM to reduce the complexity due to the interaction forces ∇V . At each iteration, we randomly divide the P particles to n = P/p batches, where each batch C is of size p. The pmmLang driven by the smooth potential U 1 (q) within the batch C is then By coupling the random-batch approximated pmmLang (67) driven by U 1 (q) and the Metropolis-Hastings algorithm for U 2 (q), we obtain the RBM with splitting Monte Carlo for the pmmLang (denoted by the pmmLang+RBM+split, presented in Algorithm 3).
To show how the splitting Monte Carlo method accelerates the sampling efficiency, we compare the performance of the pmmLang and the pmmLang+split in the numerical example below. In the mixed Coulomb-Lennard-Jones system (the potential function and the splitting scheme are defined in (69)(70)(71) in Section V), we plot in Figure 5 the time averages computed by the pmmLang and the pmmLang+split with various timesteps. Different inverse temperatures β = 1, 4 are tested respectively.
We observe from Figure 5 that the time average of the pmm-Lang is very sensitive to the timestep ∆t, while in the pmm-Lang+split we can adopt relatively large timesteps to obtain the correct thermal average, which greatly improves the efficiency of the simulation.
Despite of the satisfactory performance of the pmm-Lang+split in this example, the numerical efficiency of this Monte Carlo-type method can be influenced by the rejection rate along the sampling process. When the number of beads N or the number of particles P is large, it is likely that the highdimensional nature of the dynamics will slow down the simu-lation since the the rejection rate is high. A detailed numerical investigation of the rejection rate with different parameters is presented in Section V.

V. NUMERICAL TESTS
A. Examples of the interacting particle system In Section III we have proposed and analyzed the efficient sampling method for the quantum interacting particle systems, the pmmLang+RBM. To further explore this method, we present in this section more numerical results, where different parameters and interaction potentials are tested. The efficiency and the error in the computing thermal average Â are primarily used to quantify the numerical performance of the pmmLang+RBM. In the quantum system (1), we choose the interacting potential V (c) (q) as either the Coulomb potential or the mixed Coulomb-Lennard-Jones potential where κ = 1, σ = 0.3 and r = |q|. The kernel of (69) is provided by the Lennard-Jones potential, and is much more singular than the Coulomb potential. Hence we employ the splitting Monte Carlo method introduced in Section IV to and split the mixed potential (69) into the sum of and V (c) 2 is short-ranged. The graphs of the Coulomb potential (68) and the mixed potential (69) with its splitting scheme (70)(71) are shown in Figure 6. In the quantum system (1), the external potential is chosen to be harmonic (33), where the parameter α is chosen as where P is the number of particles. Here the choice (72) is to assure that all the particles are evenly distributed in R 3 . While placing P particles in the harmonic potential (33), the size of the potential well should be O(P − 1 3 ) to ensure the particles have O(1) pairwise distance, thus α in (33)   The position-dependent observableÂ = A(q) is chosen in the form of (53). For the Coulomb system (68), choose a(q) = V (c) (q), then is the average interacting potential of the system. For the mixed system (69), choose a(q) = e −θ |q| 2 with θ = 0.1, then B. Outline of the numerical simulation method We outline the procedure to compute the thermal average Â defined in (42) in this paper. By choosing sufficiently large N, Â is approximated by the ensemble aver-age W N (q) π . Then the pmmLang and the pmmLang+RBM estimates W N (q) π by the time averages (41) (42).
We stress that the calculation of the time averages only requires a single sampling path, thus (41)(42) are both random variables. That is to say, the values of (41)(42) depend on the Browian motion and the choices of random divisions along the samlping process. In the numerical experiments below, only one typical sampling path is used to compute the time avrages (41) (42).
We use the BAOAB scheme 49,50 to integrate the pmm-Lang dynamics numerically. By combining the velocity Verlet method 50 with the Langevin thermostat, we obtain the BAOAB scheme for the pmmLang (29) in a timestep ∆t, where η ∈ R N×3P is a random variable and each column of η obeys the Gaussian distribution N(0, (L α ) −1 ). Similar BAOAB schemes can be derived for the pmmLang within the batch (39) or the pmmLang in the splitting Monte Carlo method (58)(67).
Since the spectrum of L α is known, we cam sample from the Gaussian distribution N(0, (L α ) −1 ) using the fast Fourier transform (FFT), whose complexity is only O(N log N) (see Appendix A for a detailed mathematical formulation). Therefore, the computational cost due to algebraic operations is O(N log NP). In the interacting particle system with N not relatively large, the complexity from the interaction forces is still the dominant difficulty in the numerical simulation.
Finally, we point out that the performance of the BAOAB scheme for the pmmLang may depend on N, although the continuous-time pmmLang dynamics has been shown to have a dimension-independent convergence rate 29 . Designing a dimension-independent integrator for the pmmLang is beyond the scope of this article.

Autocorrelation
First, we aim to investigate the impacts of the RBM on the variance of the estimator, by computing the autocorrelation (abbreviated by AC) of the time averages (41)(42) numerically.
Recall that in Section III, we propose the RBM approach to reduce the complexity due to the weight function W N (q), if the observable operator is in specific forms. Thus in practice, there are two ways to implement the pmmLang+RBM: one is to use the RBM in the dynamics (Algorithm 35) but compute W N (q) as it is, and another is to use the RBM in both the dynamics and calculation of the weight function. To analyze the effects of the RBM, we will test both versions of the pmmLang+RBM. In the Coulomb system, we compute the pmmLang and the pmmLang+RBM time averages (41)(42), where the observable operator of interest is the kinetic energy and both versions of the pmmLang+RBM are considered. The autocorrelation of the time averages is shown in Figure 7. In Figure 7 and Table I, the autocorrelation of the pmmLang+RBM decays faster than the pmmLang. Since the random divisions at different timesteps are independent, it is likely that the weak correlation of the interaction forces in different timesteps leads to a small autocorrelation in the pmmLang+RBM. Still, the randomness in the batch force approximation (37) enlarged the effective variance of the estimator. Compared to the pmmLang of the version 1, the version 2 computes the weight function W N (q) approximately, thus has additional variance in the estimator.
Since the random approximation to the weight function W N (q) is necessary to reduce the computational cost per  timestep, in the following tests we always employ the pmm-Lang+RBM of the version 2. That is to say, we use the RBM in both the dynamics and the calculation of the weight function.

Convergence with the number of beads
We study the convergence of pmmLang+RBM with the number of beads N. In the PIMD representation, the ensemble average W N (q) π converges to the exact thermal average Â as the number of beads N → ∞, hence it's necessary to check if the pmmLang+RBM time averages (42) also possess this convergence property.
In the Coulomb interacting system, we plot the in Figure  8 the time averages (41)(42) of the pmmLang and the pmm-Lang+RBM, where the observable operator of interest is the kinetic energy. Different numbers of beads are tested to study the convergence of the pmmLang and the pmmLang+RBM. When the number of particles P = 16, we record the effective variance of (41) (42) in Table II. In Figure 8, we observe that pmmLang+RBM time averages (42) converge as the number of beads N enlarges. Besides, Table II shows that the variance of the estimator is not sensitive to the number of beads N, thus it's safe increase N in the pmmLang+RBM to obtain a more accurate approximation of the thermal average.
We also note that, even if N is sufficiently large, there is a bias of the pmmLang+RBM time average (42) from Â . Moreover, the bias becomes significant when the number of particles P grows large or the batch size p is small as 2. This bias though can be reduced with fixed batch size p by decreas-  ing the timestep ∆t, which are to be confirmed in the tests of the next subsection. Overall, the pmmLang+RBM time average (42) is still an accurate approximation of Â . When N = 128, the relative error of (42) from Â is no more than 2.5%. This bias though can be reduced with fixed batch size p by decreasing the timestep ∆t, which are to be confirmed in the follow-up tests. For simplicity, in the following we will take the number of beads N = 16 in both the pmmLang and the pmmLang+RBM.

Error in the calculation of the thermal average
In this part we test the error of time average (42) of the pmmLang+RBM with respect to varying timesteps ∆t. In the Coulomb interacting system, we employ the pmmLang and the pmmLang+RBM to compute the time averages (41)(42), where the observable of interest is the position-dependent one with A(q) given in (73). In Figure 9, we plot the time averages (41)(42) with different timesteps. Different numbers of particles are used to test the sampling methods. The relative error of the time averages for β = 4 is shown in Table III. In Figure 9, we observe that the bias of the pmmLang+RBM time average (42) from Â diminishes as the timestep ∆t approaches 0. Also, the bias associated with the batch size p = 4 is much less than p = 2. The numerical results above confirm our arguments on the constant C p,∆t in the weak error analysis (50).     Finally, we report the time compexity of the pmmLang and the pmmLang+RBM in Table IV. When the number of particles P is large, the pmmLang+RBM is much more efficient than the pmmLang, while the relative error always still keeps small.

Rejection rates
We focus on the mixed Coulomb-Lennard-Jones system, where the pairwise interacting potential V (c) (q) is very singular. As we have demonstrated in Section IV, we shall use the splitting Monte Carlo method to avoid using too small timesteps. In this case, the rejection rate is an important index to assess the numerical efficiency of the method. Below we compute the rejection rates of the pmmLang+split and the pmmLang+RBM+split with different parameters including N, P, p and ∆t.
In Table V and VI, we fix the number of beads N = 16 and the number of particles P = 16 respectively and compute the rejection rate of the pmmLang+split and the pmm-Lang+RBM+split. The inverse temperature β = 4, the mass m = 1, and different batch sizes p and timesteps ∆t are used.
Under the proper scaling α = P − 2 3 specified in (72), Table  V shows that the rejection rate gently grows with the number of particles P, and can be significantly reduced by shrinking the timestep. For example, when N = 16 and P = 32, the timestep ∆t = 1/16 makes sure the rejection rate not larger than 25%. Additionally, the application of the RBM does not change the rejection rate too much.
However, when we fix the number of particles P, Table VI shows that the rejection rate rapidly deteriorates with the number of beads N, even in the pmmLang+split without the use of random batches. This is due to the growing N makes a particle q i ∈ R N×3 more easily to collide with other particles. How to properly design the splitting Monte Carlo method to overcome the effects of the singular potential within the PIMD framework remains a topic to study in the future.

Error in the calculation of the thermal average
Next, we test the error of the time average computed by the pmmLang+RBM+split with respect to varying timesteps ∆t. In the mixed Coulomb-Lennard-Jones system, we employ the pmmLang+split and the pmmLang+RBM+split to compute the time averages, where the observable of interest is the position-dependent one with A(q) given in (74). In Figure  10, we plot the time averages computed by the pmLang+split and the pmmLang+RBM+split the with different timesteps. Different numbers of particles are used to test the sampling methods. The relative error of the time averages for β = 4 is in Table VII. It can be seen from Figure 10 and Table  VII that the error estimation (50) also holds for the pmm-Lang+RBM+split. That is to say, the bias of the time average computed by the pmmLang+RBM+split diminishes as the timestep ∆t → 0 or the batch size p increases. Finally, we report the time complexity of the pmm-Lang and the pmmLang+RBM in Table VIII. Comparing the results in Table IV and Table VIII, we observe that the pmmLang+RBM+split is slightly slower than the pmm-Lang+RBM, but is still much more efficient than the pmm-    Lang+split. Therefore, the additional calculation of the potential different U 2 (q * ) − U 2 (q) in the splitting Monte Carlo method does not increases the total computational cost too much.

VI. CONCLUSION
We have proposed the pmmLang+RBM, an efficient sampling method of the quantum interacting particle system in the PIMD framework. The pmmLang+RBM properly combines the preconditioned mass-modified Langevin dynamics (pmm-Lang) and the random batch method (RBM) to resolve the stiffness of the ring polymer and reduce the complexity due to interaction forces. In the pmmLang+RBM, the computational cost due to interaction forces in a timestep is reduced from O(NP 2 ) to O(NP), where N is the number of beads in the ring polymer and P is the number of particles.
In the extensive numerical tests, the pmmLang+RBM shows fine performance in the calculation of the thermal average. The pmmLang+RBM shares a similar convergence mechanism with the original pmmLang, and results in a small bias from the target distribution even for large N and P. Nevertheless, a rigorous error estimation of the pmmLang+RBM remains to be studied further.
Under the circumstances of singular interacting potentials (e.g., the Lennard-Jones potential), we introduce the pmm-Lang+split and the pmmLang+RBM+split to avoid using extremely small timesteps. When the pmmLang+split has a low rejection rate, the pmmLang+RBM+split greatly reduces the computational cost per timestep, keeps small error in computing the thermal average, and does not increase the rejection rate too much. It will be valuable to explore how to control the rejection rate when N or P is large.
Appendix A: Discussion on the preconditioning methods

Design of the preconditioning methods
In this section we briefly introduce the theory of the preconditioning methods in the PIMD, which aim to resolve the stiffness in the ring polymer. As we have shown in Section II, the stiffness originates from the stiffness matrix L ∈ R N×N , which is defined by Now let's focus on the Hamiltonian dynamics of the ring polymer system, dq = −M −1 pdt, dp = −Lqdt, where M ∈ R N×N is the positive definite mass matrix. The first step of preconditioning is to use a coordinate transformation q = Dq to decompose (A1) into different modes. To be specific, using the transformatioñ q = D −1 q,p = D T p, we can rewrite (A1) as dq = −(D T MD) −1p dt, dp = −(D T LD)qdt, where we require D T MD and D T LD to be both diagonal. In this, a preconditioning method in the PIMD is designed via the following steps: 1. Find a transformation matrix D ∈ R N×N such that 2. For some suitable constants µ 1 , · · · , µ N > 0, choose the mass matrix M = D −T diag{µ 1 , · · · , µ N }D −1 (A5) The transformed Hamiltonian dynamics (A3) is then decomposed into different modes, where the frequency of the k-th mode is and the stiffness is resolved if the frequencies ω k are uniform for k = 1, · · · , N. In the following we introduce two specific preconditioning methods, the staging coordinates and the preconditioned mass-modified Langevin dynamics (pmmLang).

Staging coordinates
The staging coordinates transformation is given bỹ q 1 = q 1 ,q k = q k − (k − 1)q k+1 + q 1 k , k = 2, · · · , N (A7) which directly yields the transformation matrix D. Using this transformation, the ring polymer potential becomes where m 1 = 0 and m k = k k−1 for k = 2, · · · , N. Therefore L ∈ R N×N is diagonalized as in (A4), where In the staging coordinates method, the constants µ 1 , · · · , µ N are chosen to be thus the mass matrix M ∈ R N×N is given by where E 11 ∈ R N×N is the matrix with only (1, 1) entry equal to 1. If one attempts to use the staging coordinates in the physical coordinates q, p, then it will be a tough task to deal with the algebra of M in (A11).

Preconditioned mass-modified Langevin dynamics
The pmmLang is based on the spectral decomposition of L, and the corresponding transformation matrix D is orthogonal. For simplicity, assume the number of beads N is even. Then L ∈ R N×N is orthogonally diagonalized as in (A4), where The constants µ 1 , · · · , µ N are simply chosen to be µ k = λ k + α, k = 1, · · · , N (A12) where the regularization parameter α > 0 is a fixed constant, and the corresponding mass matrix is M = L α := L + αI.
The advantage of the pmmLang is that it is simple can be directly applied in the spatial coordinates q. On the contrary, the mass matrix of the staging coordinates is complicated in the spatial coordinates q, and is inconvenient to combine with the RBM, hence we choose the pmmLang as the preconditioning method in the interacting particle system.
Finally, we discuss the complexity due to algebraic operations in the pmmLang. In the BAOAB scheme in a timestep, one needs to solve the linear system and compute the matrix-vector multiplication to obtain the Gaussian random variable η ∼ N(0, (L α ) −1 ). Since L α is tridiagonal, the complexity of the linear system is O(NP). With the use of the spectrum of L given above, the matrix-vector multiplication (L α ) − 1 2 ξ can be calculated by the fast Fourier transform, and the complexity is O(N log NP). In conclusion, the complexity due to algebraic operations in the pmmLang in a timestep is O(N log NP).

Appendix B: Splitting Monte Carlo method
In this section we establish the detailed balance for the splitting Monte Carlo method, and prove that the corresponding second-order Langevin dynamics preserves the desired Boltzmann distribution. To simplify our arguments, consider the target Boltzmann distribution where q, v ∈ R d is the position and velocity, M ∈ R d×d is the positive definite mass matrix and β > 0 is the inverse temperature. The Langevin dynamics which preserving the distribution π(q, v) is given by where B is the standard Brownian motion in R d , and γ > 0 is the friction constant. The detailed balance for (B2) reads 48 : Theorem 2 Let T ((q, v), (q , v )) be the transition probability density of the Langevin dynamics (B2) in time t, then the detailed balance holds: π(q, v)T ((q, v), (q , v )) = π(q , −v )T ((q , −v ), (q, −v)), (B3) which implies π(q, v) is the invariant distribution of (B2): R 2d π(q, v)T ((q, v), (q , v ))dqdv = π(q , v ). (B4) When the potential function U(q) is singular, we split it into where U 1 (q) is smooth and U 2 (q) is short-ranged and singular. Then we obtain the splitting Monte Carlo method method (Algorithm 4), which preserves the Boltzmann distribution π(q, v) and avoids the gradient of the singular part of U 1 (q).
The detailed balance of this method is given in the following theorem.