An effective Hamiltonian for the simulation of open quantum molecular systems

We discuss the derivation of an effective Hamiltonian for open quantum many-particle systems. The aim is to define an operator that can be used for (molecular) simulations where, through the exchange of energy and matter with the surrounding environment (reservoir), the number of particles, n, becomes a variable of the problem. The Hamiltonian is formally derived from the Von Neumann equation; specifically, we derive an n-hierarchy of equations for the density matrix, ρ^n , for near equilibrium situations. Such a hierarchy, in case of stationary equilibrium, delivers the standard grand canonical density matrix as it would be expected. We report that a similar Hamiltonian was conjectured, from empirical considerations, in the field of superconductivity. Thus our result also provide a formal basis for this long-standing hypothesis. Finally, an application is discussed for Path Integral simulations of molecular systems.


Introduction
The simulation of open quantum many-particle systems that exchange matter with the environment is becoming a pressing request by modern research and associated technology (see e.g.[1] and references therein).In current research a systematic description of an open quantum system derived from the Von Neumann equation.i.e. from first principles, is achieved by the Lindblad equation and related approaches [2][3][4].However, the equation is used mostly for * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.problems of transport (see e.g.[2][3][4][5][6][7] and references therein) and, at the best of our knowledge, does not describe the situation of open systems in equilibrium with a reservoir where the number of particles naturally fluctuates in a grand canonical fashion [2,3,8].On the other hand, an effective Hamiltonian for open quantum systems in equilibrium with a particle's reservoir has been conjectured long time ago in the field of superconductivity; such a conjecture is based on physical intuition, but so far has not been explicitly derived from the Von Neumann equation (see e.g.[9,10]).Systems in equilibrium are of high interest in current research thus, the key question to address is whether it is possible to explicitly derive an n-dependent operator, with n being the variable number of particles, physically consistent and mathematically well defined, so that one can treat the evolution of a system that exchanges energy and matter with a bath at equilibrium or near equilibrium and as a consequence is characterized by particle number fluctuations that are not induced by any external source.For near equilibrium it is meant that the deviation from equilibrium can be eventually added to the system as a (small) perturbation As an example of system of interest, a relevant class of problems that falls in the category is that of numerical algorithms for ab initio molecular dynamics for open molecular systems [1,11,12].In such systems, within the Born-Oppenheimer approximation [13], that is under the hypothesis of nuclear and electronic scale separation, nuclei are considered classical objects while electrons are treated at quantum mechanical level; these latter require a proper quantum treatment as open system and a corresponding numerical prescription for calculations.This field is, in our view, still in its infancy, since most of the current open-system ab initio algorithms are built on empirical basis and thus are often numerically not reliable [14] due to conceptual inconsistencies that lead to incorrect statistical sampling [1,15].This means that the construction of a solid theoretical framework would assure physical consistency and computational robustness to the underlying numerical methods, as it occurred for molecular dynamics of classical and semi-classical open systems with a system-reservoir exchange of particles [16][17][18][19][20].
Using basic statistical mechanics arguments, for the ab initio treatment of open molecular systems the following condition has been proposed: the electronic chemical potential of the open system must have the same value of the electronic chemical potential of reference (i.e. of a large system) [12]; the results of this work provides an explicit formal justification to the idea described above.Furthermore, the current model sets the formal basis for moving beyond equilibrium by adding, on the equilibrium set up, possible physical perturbations to the molecular system.
Actually, in molecular simulations of open systems where the quantum effects of spatial delocalization of atoms in space are included, the numerical algorithm makes use of the effective Hamiltonian derived in this work.This specific example of practical application will be explicitly discussed in the paper.From the conceptual point of view, the key aspect at the basis of this work, is that for situations of equilibrium or near equilibrium, the fundamental principles of quantum statistical mechanics state that the density matrix of a quantum system takes the form of a grand canonical operator (see e.g.[21]).In this context, the necessary request to assure physical consistency for an equation with system-reservoir exchange of particles is that it should deliver the grand canonical density matrix in the limiting case of full statistical equilibrium.As anticipated above, interestingly, from empirical considerations within the framework of statistical mechanics, N.N.Bogoljubov in the treatment of superconductivity introduced a scalar parameter (identified a posteriori as the chemical potential of the system) associated to the particle number operator with the task of fixing the particle's number variation and conjectured an effective Hamiltonian for open systems with similar form as the one that we will derive here [9,10].However, in our work, differently from the empirical derivation of [9,10], we propose a formal derivation of the effective Hamiltonian starting from the Von Neumann equation for the density matrix of the entire system.Next, having in mind systems typically treated in molecular simulations, we integrate out the degrees of freedom of the bath region and derive an n-hierarchy of equations for the density matrix of the subsystem of interest, ρn .This procedure is conceptually equivalent to the procedure followed for the Liouville equation in the classical case [16,17] developed within the framework of molecular dynamics [22].The novelty of our proposal lies in the explicit (step-by-step) formal derivation of the equation from first principles, having in mind a molecular system for numerical simulations.The relevant implication is that the procedure employed allows us to explicitly define the physical conditions of validity of the model; this aspect may be rather relevant in a computational set up.It must be noticed that, in our knowledge, the principle that the grand canonical density matrix should be automatically and explicitly obtained from the Von Neumann equation in the limiting case of equilibrium has not been applied to any of the currently used quantum master equations where the number of particles are a variable of the problem.

Quantum grand canonical ensemble
In statistical mechanics, an open system of particles, S, that exchanges energy and matter in equilibrium with a reservoir B, is characterized by a grand canonical density matrix: ρn = 1 QGC exp(− 1 kT ( Ĥn − µn)), where Ĥn is the Hamiltonian operator at n particles in the system S, where n is a variable quantity and n the number operator that counts the number of particles in a given quantum state.Finally, µ is the chemical potential, that is a quantity that regulates the exchange of particles with the bath and Q GC a normalization factor.In essence, µ can be seen as a rate of change of the free energy of the system w.r.t. the change of number of particles, and n, the particle number operator [21].This latter, in a Fock space formalism, is defined as [23]

Quantum equation in the limiting case of equilibrium: grand canonical operator
The approach employed in this work to derive a quantum equation for the density matrix of a (sub-)system embedded in a large reservoir of energy and particles starts from the Von Neumann equation [24] for the density matrix ρ of a (large) total system with a given Hamiltonian Ĥ: In the following, instead of a pure dynamical view, as done for example in the derivation of the Lindblad equation, we consider a statistical view, that is consider a partitioning of the total system with N particles in a subsystem S with n particles and a reservoir B with N − n particles, as done for the Liouville equation in the case of classical systems [16,17].Next we trace out the degrees of freedom of B in the Von Neumann equation, considering all the partitioning of N particles in n and N − n sets.The result is a hierarchy of equations for ρn (t) = N! (N−n)!n! ρS (t), which in the case of stationary state of equilibrium must deliver the grand canonical density matrix: ρn = 1 QGC exp(− 1 kT ( Ĥn − µn)) as a solution.For hierarchy here is meant that one obtains a set of equations, one for each value of n, so that the n-dimensional operator ρn represents the n-fluctuating ensemble for the open system.

Partitioning the total system into a system of interest and reservoir
We consider the following Hamiltonian for the total system: Ĥ = ĤS ⊗ 1 B + 1 S ⊗ ĤB + Ĥint ; where the notation 1 S and 1 B represent the identity operator of the Hilbert space associated to the domain S and B respectively.In essence, we consider the total Hamiltonian partitioned in the Hamiltonian of particles in S, the Hamiltonian of the particles in B and the Hamiltonian of interaction between the particles of S and B.
We assume || Ĥint ||≪|| ĤS || always, this is the so-called surface-to-volume ratio approximation used in the statistical derivation of grand canonical ensemble that can be found in the most popular textbooks of statistical mechanics (see e.g.[21]).In essence, one neglects the direct interaction between S and B and considers the total Hamiltonian partitioned between the part of S and the part of B: n , thus the particle exchange between S and B happens according to the balance of free energy between S and B. The free energy of S and the free energy of B are characterized by the particle-particle interactions between particles only in S and particles only in B respectively, and by a particle-particle interaction between particles in S and particles in B (direct interactions).Under the hypothesis that S is large enough, the direct particle-particle interaction between S and B is negligible compared to the effect of all the other interactions, e.g.surface-to-volume ratio hypothesis (see e.g.[21]).The density matrix of the whole system needs to be defined according to the partitioning in a system S with n (instantaneous) particles and a bath B with N − n (instantaneous) particles, this means that, for indistinguishable particles, one must consider all possible partitioning (realizations) of N in n particles for S and N − n particles for B, i.e.

N!
(N−n)!n! .It follows that a consistent definition of the total density matrix is (see also [16]): with ρSn being the n-realization of the density matrix that represents S and ρBN−n the corresponding N − n-realization of the density matrix that represents B.
We then define: At this point we make a first step in the modeling, that is we employ the Born-Markov approximation [25]: In essence we assume that the bath is large enough and thermodynamically in a stationary state where the exchange of particle with S can be ignored (for the properties of B).One can formally see it also in the following manner: one can write ρB = 1 QB γB , where Q B is the partition function of B, that is: Q B = Tr B γB and γB is the density matrix operator of B not normalized, for example in the canonical form: γB = e −βHB .It is known from statistical mechanics that −kT ln Q B = F B which is the free energy of B. Since we assumed that B is large enough and in equilibrium so that its properties do not significantly change with time, then the exchange of particles with S can be ignored (w.r.t. the properties of B), thus F B can be considered constant and thus independent of N − n, it follows: Tr B ρB (0) = 1.The followed approach corresponds to the standard procedure of statistical mechanics when the grand canonical partition function is derived from statistical considerations (see e.g.[21]).Proceeding further with the model, since this is a statistical approach where in S the particle number n can take values from zero to N, one has to count all such possible realizations in the ensemble and thus: n=0,N Trρ n = 1, which implies that Trρ = n=0,N Tr S Tr B N! (N−n)!n! ρSn ρBN−n = 1, as it occurs in a grand-canonical ensemble [16,21].
In the next to avoid a heavy formalism we will indicate S for the system and B for the bath without carrying further the particle number label n and N − n.
At this point, in order to reduce equation ( 1) only to the degrees of freedom of S, we integrate both sides of the equation w.r.t. the degrees of freedom of B, that is we calculate the trace over B: Due to additivity of the Lie bracket, the r.h.s of equation ( 4) becomes Tr B [ Ĥn We have: For the next term we have: It can be notices that Tr B ( ĤN−n ρB ) and Tr B (ρ B ĤN−n ) are nothing else than the ensemble averaged energy of B, which we denote by ⟨E B ⟩.This quantity depends on the number of particles of S and B (i.e.n and N − n).Here comes the crucial point of our modeling proposal, as explained below.⟨E B (N − n)⟩ is an unknown function of the variable n, if we consider the situation of equilibrium between S and B, then we can reasonably expand the energy of B in Taylor series at the first order: The term ⟨E B (N − n)⟩ is not known explicitly, the term ⟨E B (N)⟩ is not crucial since this is a constant that can be considered as the zero of the energy scale, instead − ∂⟨EB(N−n)⟩ ∂n | n≪N is actually an explicitly known quantity that characterizes the thermodynamics of the system.It corresponds, by definition, to the chemical potential µ, at constant entropy and volume for B (since we consider B in stationary equilibrium, see also [26]).It must be noticed, as expressed in the formalism above, that the expansion in n must not be seen as n → 0, but as an expansion in the limit N≫n, so that effectively n is negligible compared to N. This hypothesis implies that the effects of a variation of n on the physics of B is negligible, that is the macroscopic properties of B are determined as if B was a canonical ensemble with, effectively, a fixed number of particles N (see e.g.[21]).In practical terms, the validity of equation ( 5) strongly relies on considering a large B for which the variation of its macroscopic physical quantities, due to a variation of n, is negligible w.r.t. the numerical accuracy chosen for the calculations of physical quantities in the open system.
Moreover, the term − ∂⟨EB(N−n)⟩ ∂n n = µn must be interpreted within the framework of quantum statistical mechanics, that is, n in equation ( 5) must be considered as an operator, n, which, at given thermodynamic conditions consistent with µ, counts the number of particles of the corresponding quantum state of S. Next, we notice that while B is considered large enough so that the change of number of particles, from the physical point of view does not affect its properties (as it was at fixed number of particles), the same cannot be said for S. In fact, a change of number of particles carries sizable fluctuations in the statistical properties of S since ⟨E B (N − n)⟩ enters into equation (4) 'acting on' ρS .Thus for situations close to equilibrium equation ( 4) is reduced to: and for the exact equilibrium with d ρn dt = 0, one would recover the grand canonical density matrix for S: ρn = 1 . Note that equation ( 6) leads to a solution for ρn that preserves the trace as: The correspondence between classical and quantum statistical mechanics used above to expand ⟨E B (N − n)⟩ and write it in operator form is the main proposal of our model; it is inspired by the procedure usually used in literature to define statistical and thermodynamic quantities in quantum mechanics starting from the known classical case [21].Once the model above is accepted, near equilibrium one may add any source of energy of (small) perturbation (e.g.particle current, heat, external fields) on the r.h.s. of equation (6).The meaning of 'near equilibrium' varies from system to system and problem to problem, a general criterion would be to consider perturbations that lead to a linear response with an input/perturbation of energy, with respect to the energy of equilibrium of the system, within a certain threshold, so that the physical hypothesis of 'small perturbation' holds (e.g.10% of energy input/output, see e.g.[27]).The extended equation would then allow for a larger set of possibilities compared to the physical situations treatable by other approaches, e.g. the standard Lindblad equation.

Empirical derivation of Ĥ − µn v.s. our first principles derivation
As anticipated above, an effective Hamiltonian similar to our result, ĤGC = Ĥ − µn was conjectured by Bogoljubov in the framework of superconductivity for open quantum systems.The argument used for advancing such a proposal is based on the physical intuition that the result obtained from the canonical ensemble can be carried over to the grand canonical ensemble by extending the canonical Hamiltonian with a term that regulates the variation of number of particles at the rate expressed by the chemical potential, µ, as known from thermodynamics.That is, it seems reasonable to define a grand canonical effective Hamiltonian: ĤGC = Ĥ − µn so that for similarity to the canonical density matrix, ρC = e −β Ĥ, one would obtain: ρGC = e −β ĤGC = e −β( Ĥ−µn) , that is the well known result of statistical mechanics regarding the grand canonical density matrix.It must be added that in literature the results of quantum statistical mechanics are derived through a formal similarity to the classical results, that is, the canonical or grand canonical density matrix are a formal extensions of the classical results to the Hilbert space (see e.g.[28] and references therein).In this sense, the Hamiltonian ĤGC = Ĥ − µn, in [9,10], is based on an empirical formal extension of a classical result, it is not derived from a rigorous first principles equation of quantum mechanics (Von Neumann equation).
The approach shown here instead is not empirical in such sense, in fact we start from the Von Neumann equation, that is an equation that expresses the very nature of a quantum system, divide the system in an open system and a large reservoir, make physical hypothesis about the reservoir and integrate its degrees of freedom in the equation.The relevant consequence is that the chemical potential is not empirically 'imposed' a priori, but is formally obtained as a result of the integration of the reservoir (as the constant of the first order expansion in the number of particles for the system-reservoir coupling).This means that we do not impose the grand canonical density matrix as a constraints of our model, instead we verify that our Hamiltonian automatically leads to the grand canonical density matrix as a natural solution of the Von Neumann equation for an open system in a stationary state.
In essence, our procedure embeds a formal justification of the empirical model of [9,10].Furthermore it clarifies that Ĥ − µn is a first order approximation within the system-reservoir coupling and the error is in the truncation terms of the Taylor series w.r.t. the number of particles of the open system.
It is important to notice that while, for technical reasons, ĤGC = Ĥ − µn, may be popular in fields such as theory of superconductivity, for quantum molecular simulations this Hamiltonian has not been considered yet, except in this work (see the next section).In fact in simulations, the variation of number of particles is implemented through a stochastic process within a Monte Carlo approach (see e.g.[29]), where particles are inserted or deleted according to the chemical potential.Instead, ĤGC = Ĥ − µn, as derived by us, allows for a fully dynamic picture, without the need of separating the purely dynamic part, generated by Ĥ, with a stochastic part regulated by µ.As a practical conclusion, one of the major advantage of the extended Hamiltonian is that it allows for the possibility of designing very efficient Molecular Dynamics algorithms for open systems [1,12,22].

Example of a practical application: path integral molecular dynamics of open systems
The Path Integral formalism of Feynman can be used in molecular dynamics (MD) to account for the quantum delocalization in space of atoms/atomic nuclei (quantum nuclear effects).The essence of the idea (for the specific details see the appendix or appendix C and D of [11]) is that the quantum Hamiltonian of an atomic or molecular system can be mapped onto an effective classical Hamiltonian where the atoms in MD, or the nuclei in ab initio electronic structure calculations, are represented as ring polymers with beads connected via springs.The three dimensional fluctuation in space of such polymer rings describes the quantum atomic or nuclear delocalization in space.Moreover, the sampling of configurations in space, via classical trajectories produced by the MD, results in the statistical sampling of the corresponding atomic/nuclear quantum density matrix.As a consequence any physical quantity averaged over the MD trajectory corresponds to the quantum average of the quantity in question.Of specific interest for us is the Hamiltonian used for doing such simulation in the case of open systems that exchange energy and matter with a reservoir.In fact, in such simulations [30] the effective classical Hamiltonian contains an additional one-particle potential term that acts (mostly) at the border between the molecular system and the reservoir and regulates the system-reservoir exchange of particles.The system-reservoir exchange occurs in a way that the resulting statistics reproduces the Grand Canonical distribution at the (correct) chemical potential of the system, this latter corresponds approximately to the statistical average of the one-particle potential.In formal terms, the MD of the polymer rings (i.e. the effective classical simulation on which the quantum Hamiltonian has been mapped) contains the standard path integral Hamiltonian plus an additional, one-particle potential term: , with ⟨ϕ i (x)⟩ ≈ µ and n the instantaneous number of particles in the system [31][32][33][34].The corresponding MD simulation, as explained above, samples effectively the quantum, grand canonical, density matrix of the molecular system.
As a matter of fact the classical Hamiltonian H PI − n i =1 ϕ i (x), since it samples the grand canonical distribution, consistently maps on an effective classical problem the Hamiltonian of an open quantum system in equilibrium, that is, it effectively maps the Hamiltonian of equation ( 6) Ĥn − µn in the stationary case of equilibrium.
In the path integral MD the term, n i =1 ϕ i (x), as underlined above, regulates the net number of particles that according to the statistical rate (µ, the chemical potential) is allowed to enter into/exit from the system so that the statistics satisfies the grand canonical distribution of equilibrium, which, effectively, is the quantum grand canonical density matrix: In a quantum view of the molecular system, the equivalent of the action of n i =1 ϕ i is the action of the operator µn, i.e. counting the instantaneous number of particles of the system according to the rate µ so that it satisfies the grand canonical density matrix at equilibrium (i.e.[ Ĥn − µn, ρn ] = 0).As a consequence, the effective classical Hamiltonian in the path integral MD, H PI − n i =1 ϕ i (x), when mapped back to the quantum Hamiltonian is physically consistent with the expression: Ĥ − µn of equation ( 6).Furthermore, since ⟨ϕ i (x)⟩ ≈ µ, one has that ⟨ n i =1 ϕ i (x)⟩ ≈ µ⟨n⟩ and accordingly in quantum terms one would have that ⟨ψ|µn|ψ⟩ = µ⟨n⟩, which enforces the statistical correspondence between the classical operator used in path integral simulations of open systems and the quantum operator derived in this work.In essence, the result of this paper gives the needed conceptual legitimization to an otherwise purely empirical numerical approach.Until now, due to the lack of conceptual justification for the effective Hamiltonian used, the application of the open system approach to Path Integral Molecular Dynamics has often required a case-by-case numerical validation.The numerical validation was done running additional and rather expensive full quantum simulations of large systems [19,20,30,35,36].The conceptual justification of the effective Hamiltonian shown in this paper implies that a specific case-by-case numerical validation is not required.

Conclusions
We have addressed the problem of system-reservoir exchange of energy and particles for quantum systems with the aim of deriving an effective Hamiltonian for molecular simulations.We start from the Von Neumann equation for the density matrix of the whole (large) system, and model it according to the partitioning of the total system in a system of n variable particles and a large stationary reservoir; furthermore we focus on the situation of equilibrium or near equilibrium.Next, by integrating out the degrees of freedom of the environment we obtain an equation for the density matrix of the subsystem of interest carrying terms of interaction with the environment/reservoir.The result is a n-hierarchy of equations for the density matrix of the system ρn , where n is the number of particles governed by an effective Hamiltonian, Ĥn − µn.
Since n is a variable quantity it follows that each ρn represents one of the possible statistical realizations; in the limit of stationary equilibrium we recover a known solution, that is the grand canonical density matrix as it would be expected for naturally fluctuating open system in equilibrium.The procedure employed for the derivation of the hierarchy of equations is similar to the procedure employed for a system of classical particles [16,17]; such a procedure was used as a theoretical platform to design algorithms of (classical) molecular simulation [22].We have pointed out that our derivation can be also considered as the explicit formal derivation of a result that was empirically conjectured within the field of superconductivity and that, in addition, now can be extended to the larger field of quantum molecular simulations.At numerical level, the theoretical model proposed here already gives sufficiently solid ground to computational methods for studying open molecular systems.In order to support such a claim we discussed how the model proposed here forms the theoretical basis of algorithms that are successfully used in Path Integral simulations of open molecular systems.In perspective, the proposed model can be used as a theoretical basis for ab initio simulations of open systems where molecules explicitly carry electrons.Existing algorithms for such systems are based on empirical considerations only [1,12].This field, despite its potential relevance, is, in our view, still in its infancy and the search for a solid formal ground with the associated conceptual development of the subject is highly desired.

Appendix. Path integral formalism and its application in molecular dynamics
In this section, we report the basic information about the the path integral formalism within the framework of quantum statistical mechanics and its implementation in Molecular Dynamics to simulate quantum spatial delocalization in space of atoms of nuclei.A more complete treatment can be found in [37].
The starting point is the Hamiltonian of a single particle subject to a potential U: with K and Û the kinetic and potential energy operators, respectively, fulfilling the commutation relation: [ K, Û] ̸ = 0.The element in the space representation of the density matrix are written as: K and Û do not commute, thus to the aim of feasible calculation the Trotter theorem is used.For two non-commuting operators, A and B, the Trotter theorem reads: where P is known as the 'Trotter number'.The application of the theorem to the density matrix element gives: P − 1 insertions of the identity operator, Î = ´dx|x⟩⟨x|, between each Ω gives: As a consequence, the calculation of the matrix elements is simplified: Concerning the kinetic operator, one introduces the identity operator in momentum space: whose integration delivers the following expression: and, by substituting in the expression of the density matrix, one obtains: The partition function of the system is calculated as the trace of the density matrix within a certain interval [0, L] (the linear size of the box containing the particle): in explicit terms: it can be concluded that the 'effective' Hamiltonian corresponds to: with ω P = √ P βℏ and x P+1 = x 1 .The result above corresponds (and it is known in literature as) the discretized path integral (quantum) representation of the partition function of a single particle.As a matter of fact, it can be seen as the partition function of a polymer ring with P beads harmonically linked to each other via nearest neighbor connections, with ω P being the coupling strength.The extension N interacting particles is, in principle, not trivial since the statistics/symmetry (fermionic or bosonic) must be included.However, in first approximation, independently on whether the particles are fermions or bosons, the basic quantum effects of particle spatial delocalization do not need the introducing of the particle symmetry.
As a consequence, for a Hamiltonian of N particles in d-dimensions: Ĥ = A natural consequence of the formalism is that for two particles (atoms) in their bead representation, the atom-atom interaction implies that only the beads with the same index k do interact with each other.The density quantum density matrix (or the partition function) obtained above can be sampled using a an effective fictitious dynamics as that produced by classical Molecular Dynamics, as reported in the subsection below.

A.1. Path integral formalism in molecular dynamics
The effective fictitious dynamics for MD simulations requires to have momenta associated to each particle (atom), which are not present in the effective ring polymer Hamiltonian.This problem can be solved by writing: ( mP 2π βℏ 2 ) = ´dp 1 . ..dpP (−β ) that is by artificially adding P-Gaussian integrals corresponding to fictitious momentum variables p 1 , . .., p P where m ′ = mP (2π ℏ) 2 , is an arbitrary mass parameter.Here we want only to show the general idea so we can restrict our attention to a system of equivalent atoms, thus m j = m a for the j-atom and m ′ j = m ′ for the fictitious mass.The extension to multi-atoms components is indeed straightforward.
The partition function is written as: to be used for the MD simulation.
In this way we have mapped the quantum problem on an effective sampling in space of interacting polymer rings via Molecular Dynamics, that is, while the dynamics is fictitious, the space sampling of the conformations is statistically equivalent to the quantum statistics of interacting atoms which are not rigid spheres but space-delocalized entities (the fluctuating shape of each ring).