{'fraction_non_alphanumeric': 0.036337235040994495, 'fraction_numerical': 0.012962405482460023, 'mean_word_length': 4.509463414634146, 'pattern_counts': {'":': 1, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Importance Sampling in Rigid Body Diffusion Monte Carlo

We present an algorithm for rigid body diffusion Monte Carlo with importance sampling, which is based on a rigorous short-time expansion of the Green's function for rotational motion in three dimensions. We show that this short-time approximation provides correct sampling of the angular degrees of freedom, and provides a general way to incorporate importance sampling for all degrees of freedom. The full importance sampling algorithm significantly improves both calculational efficiency and accuracy of ground state properties, and allows rotational and bending excitations in molecular van der Waals clusters to be studied directly.


Introduction
Diffusion Monte Carlo (DMC) methods provide a powerful theoretical approach to the study of weakly bound clusters of atoms and molecules. The favorable polynomial scaling of DMC enables clusters with large numbers of degrees of freedom to be studied, provided that suitable schemes for sampling the multi-dimensional wavefunctions exist. In the case of van der Waals clusters containing one or more molecules, there often exist several different time scales, deriving on the one hand from the monomer internal degrees of freedom, and on the other hand from the 'external' monomer translations and rotations which are associated with the van der Waals modes between monomers. Typically, there is a large time scale separation between the high frequency internal vibrations of the monomer constituents, and lower frequency intermonomer van der Waals modes. This time scale separation is often invoked to motivate approximations in which one effectively freezes the high frequency degrees of freedom when studying the role of the low frequency modes. When applied to intramolecular degrees of freedom using Diffusion Monte Carlo, this approximation leads to an algorithm first proposed by Buch [1], commonly referred to as 'Rigid Body Diffusion Monte Carlo' (RBDMC), in which all monomers are treated as non-vibrating, rigid molecules which are free to rotate and translate in the presence of intermolecular interactions.
The first implementation of RBDMC for molecular clusters by Buch employed a short-time Green's function for the free rotation of each molecular monomer in its instantaneous principal axis frame [1]. A convenient Gaussian short-time approximation for such monomer rotational motion was assumed in that work. A number of applications of RBDMC have since been carried out, all of which use a combination of this angular sampling for monomer rotation in local monomer principal axis frames, together with the conventional sampling of translational motions of monomers in the laboratory frame [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. The rigid body rotations are carried out using either Euler angles, or quaternion algebra [16]. Considerable use has been made of this algorithm for water clusters, where the floppy intermolecular degrees of freedom are difficult to study with other methods for clusters larger than the water dimer [1,2,4,6,7,10,11,[14][15][16].
Nearly all of the RBDMC applications to date employ DMC in its simplest form, namely unbiased Monte Carlo sampling in which the ground state wavefunction is sampled directly without using any trial function. Exceptions to this include a calculation for HF-Ar n for which importance sampling of translational degrees of freedom was implemented [3], and a calculation for rotational excitations in clusters of SF 6 He n [17]. The latter calculation represents the first application of the complete importance sampling algorithm for rotations and translations which is described in detail in this work. Such biased sampling proves useful since unbiased sampling, while adequate for ground states of relatively strongly bound or fairly rigid clusters, rapidly becomes less efficient for weakly bound systems. Importance sampling also allows for easier computation of quantities other than the energy (e.g. structural quantities) that would require descendant weighting [18,19] in unbiased DMC. Furthermore, it is also more amenable to study of excited states than are unbiased sampling techniques.
Biased sampling in DMC is generally achieved with the importance sampling approach first proposed by Kalos and co-workers for translational motions [20]. The improved efficiency is dramatic for very weakly bound systems such as the clusters of helium and hydrogen, which are dominated by quantum effects.
Both pure helium clusters [21] and doped helium clusters [22] have been studied with this approach. Excited states can be obtained from fixed node calculations in which there is no explicit trial function but where nodal constraints are nevertheless imposed directly in an otherwise unbiased sampling procedure [23]. More efficient sampling of excited states within fixed node is often possible by using explicit trial functional forms. Exploration of nodal release, nodal optimization, spectral evolution and other approaches to go beyond the fixed node approximation is an active area of research [10,[24][25][26][27][28][29][30][31][32][33][34].
In this paper we first outline the standard rigid body diffusion Monte Carlo algorithm with unbiased sampling (Section 2). Our treatment includes a quantum mechanical derivation of the short-time Green's function for free rotation of a general asymmetric top (Section 2.1). For spherical tops, formal solutions for both the imaginary time, diffusive Green's function [35] and the real time Green's function [36] were established a long time ago. However for a general rigid body, i.e., an asymmetric top, formal solutions have been restricted to eigenfunction expansions [37,38], and no rigorous derivation of the expected Gaussian form for small angle diffusion has been presented. A formal solution for this is useful since it enables extension for rotation in the presence of external force fields (to which the quantum forces in biased diffusion Monte Carlo are formally equivalent) to be naturally and rigorously implemented. In Section 2.2 we show how the use of a rigid body diffusion Monte Carlo algorithm is motivated, with a detailed discussion of the factors influencing time step choice in a weakly bound system composed of both atoms and molecules. In Section 3 we then extend the rotational Green's function formalism to enable direct importance sampling of rotational degrees of freedom. Implementation of the full importance sampled RBDMC algorithm is described in Section 4. In Section 5 we provide some basic applications of this approach which illustrate its use. First, we demonstrate explicitly with detailed tests that such shorttime sampling of the angular degrees of freedom will cover the entire angular configuration space of a rigid body (Section 5.1). Then we show with specific examples taken from studies of two weakly bound molecule-helium systems, i.e., SF 6 4 He n and HCN 4 He n , that full importance sampled RBDMC considerably improves the efficiency of calculations for ground states (Section 5.2), as well as enabling calculation of Van der Waals excitations in such systems (Section 5.3).

Theory
The diffusion Monte Carlo method provides a stochastic approach to solution of the Schrödinger equation This is accomplished by first transforming the equation to imaginary time, and then solving the resulting relaxation equation with the use of a short time approximation to the Green's function forĤ. The formal solution to Eq. (2) is Expanding the initial wave function in a complete set of eigenfunctions of H,{|ψ n }: we obtain Choosing the energy scale such that E 0 = 0 will guarantee that the ground state solution is achieved asymptotically: The solution to Eq. (3) can be expressed in position representation by making use of the short-time Green's function: where Q refers to the set of coordinates specifying a single body, and Q to the entire set for the N-body system. In an all-atom calculation, these coordinates will consist of the 3N-dimensional cartesian coordinates of the N-body system, i.e., Q ≡ R. We shall use the notation R to refer to the three-dimensional (center of mass) cartesian coordinates of a single body, and R to the 3Ndimensional cartesian coordinates of the N-body system. If m of these bodies are molecules and have to be treated as rigid bodies, we have Q ≡ R, Ω, where Ω ≡ (ϕ, ϑ, χ) denotes the Euler angles specifying the orientation of the principal axis frame (PAF) of a monomer in the laboratory reference frame, Ω denotes the 3m-dimensional set of angular coordinates for m monomers, and R refers to the combined set of 3m cartesian coordinates of the monomer centers of mass, and the 3n atomic coordinates for the n = N − m atoms.
Eq. (8) can be propagated in imaginary time τ to find the asymptotic solution, Eq. (6), given an accurate short-time Green's function. The propagation is started from some initial guess for the wave function, Ψ( Q, τ = 0), which is represented as an ensemble of multi-dimensional configuration points or "walkers": Thus the algorithmical issue reduces to finding an accurate short time Green's function for a given HamiltonianĤ, together with an efficient sampling procedure for the multi-dimensional integral in Eq. (8).
The simplest diffusion Monte Carlo procedure samples the wave function Ψ( Q, τ ) directly. This unbiased, or non-importance sampled DMC has been extensively used in the past for applications to both electronic and nuclear structure and energetics [18,23,39]. We now summarize the conventional implementation of RBDMC using unbiased DMC, for Hamiltonians describing the nuclear motion of N-body atomic or molecular systems, witĥ E R is a reference energy whose significance will become apparent below. We assume for simplicity that the interaction potential is pair-wise additive, where r ij is the vector between the centers of mass R i and R j of monomers i and j, respectively.
Most conventional solutions implicitly utilize the second-order, split operator decomposition [40] of the Green's function: For Hamiltonians involving only translational kinetic energy, i.e., where D i =h 2 /2m i , the contribution e −T t ∆τ to the short time Green's function is given by the solution to the diffusion equation. In the position representation this is Note that this corresponds formally to rigid bodies moving with fixed relative orientations specified by the set of Euler angles {Ω i }. In general, this would not be a physically useful representation. However for a single molecule moving in an atomic cluster, e.g., in a rare gas cluster, this can sometimes provide a good representation, when the molecule is sufficiently heavy that the molecular rotational kinetic energy is small compared to the rare gas kinetic energy contributions [41]. The potential term in Eq. (12) is diagonal in the position representation: In the position representation, Eq. (12) then becomes a product of non-local single particle diffusive terms from the kinetic energy, and a local term from the potential energy. Depending on its sign, the latter can lead to exponential growth or decay of the amplitude Ψ( Q) in Eq. (8).
For Hamiltonians involving only rotational kinetic energy, we have for a general rigid body, i.e., an asymmetric top, whered pα = 1/2I pα . I pα is the principal axis moment of inertia for particle p relative to its principal axis α = x, y, z, andL pα is the corresponding component of the angular momentum operator in the pth PAF: Here φ ≡ (φ x , φ y , φ z ) specifies the angles of one-dimensional rotations around the molecular principal axes x, y, and z, respectively. The rotational Green's function for a general rigid body should formally be expressed as G r (Ω → Ω ′ , ∆τ ). However, applying the second-order short-time decomposition where leads immediately to the following short-time factorization into one-dimensional rotational Green's functions: Each of the one-dimensional rotational Green's functions can be evaluated by inserting a complete set of one-dimensional angular momentum eigenstates of the angular momentum along the corresponding molecular PAF. This is given by where with eigenvaluesL pα |L m pα = mh|L m pα , m = 0, ±1, ±2, ...
Then, e.g., for α = x, we readily obtain where we have used the scaled constant d pα =h 2d pα ≡h 2 /2I pα for particle p and principal axis α. We shall see below that in the short time limit, ∆τ → 0, this constant is identical to the rotational diffusion constant for this degree of freedom. Note that it is also equivalent to the molecular rotation constant B pα . As discussed in Section 2.2 below, the magnitude of this relative to the masses of the n = N − m atoms, of the m molecules, together with the strength of the interaction potentials, will determine the time step for a given imaginary time propagation.
It remains to be shown that Eq. (25) reduces to the expected Gaussian form in the small ∆τ limit. We do this by showing the reverse, i.e., by expanding the Gaussian form in eigenfunctions of the one-dimensional angular momentum operatorL = −ih∂/∂φ: Hence For ∆τ << 1, we can extend the limits of this integral to ±∞, allowing evaluation to yield Therefore we make the identification Replacing φ in Eq. (30) by (φ px − φ ′ px ) leads to the desired Gaussian form of the sum in Eq. (25), i.e., This yields the following short-time factorized limiting form for G r i : The constants d pα are now seen to act as one-dimensional diffusion constants for each of the rotational degrees of freedom in the p-th molecular PAF, as noted above.
We note that Eq. (32) is formally overcomplete in angle space, giving rise to a total integrated angular volume of (8π) 3 instead of the requisite 4π 2 . This is a consequence of the short-time factorization into three one-dimensional rotations, Eq. (19). For a spherical top this can be avoided when the Eulerian short-time limit described below is employed. A comparison of the two shorttime propagators for a spherical top is provided in Section 5, where we show that the formal violation of the angular configuration volume deriving from the product of one-dimensional rotations does not cause problems if the time step ∆τ is sufficiently small. The short-time factorization can be improved upon if necessary by making use of the split operator decomposition [40], but we have not found this necessary in any applications to date.
For the special case of a spherical top (d pα ≡ d p , α = x, y, z), the exact form of the Green' s function is well known [36], as is its short-time diffusive limit [35] G r p (Ω → Ω ′ , ∆τ ) = where ζ is the angle of rotation about an arbitrary axis. This short-time diffusive limit was already derived by Furry in 1957 [35]. It follows directly from taking the small ζ, small ∆τ limit to the imaginary time transcription of the infinitesimal propagator for a spherical top which was derived rigorously by Schulman (see Eq. (4.9) in Ref. [36], with Γ ≡ ζ, and d p ≡h 2 /2I). Eq. (33) corresponds to the Eulerian description of rotation of a rigid body as a rotation through a single angle ζ about an arbitrary axisn ζ [42]. It can be implemented in a Monte Carlo sampling procedure by first choosing a randomly oriented axisn ζ (uniform sampling of azimuthal angle φ on the interval [0, 2π) and of cos(θ) on the interval [−1, 0]), and then sampling the angle of rotation ζ about this axis from a Gaussian distribution of width 3 × 2d p ∆τ . The origin of the factor of 3 is described in Appendix A. Related Eulerian rotational sampling has been employed in classical Monte Carlo rigid body simulations, but with the Gaussian sampling of the rotation angle ζ replaced by uniform sampling inside a fixed interval [−∆ζ 0 /2, ∆ζ 0 /2] [43]. This satisfies the conditions of reversibility, detailed balance, and accessibility of all orientations, but does not correspond to the true Green's function, Eq. (33). It will therefore not necessarily yield the correct energies and expectation values in a quantum Monte Carlo calculation.
When the Hamiltonian involves both translational and rotational kinetic energy terms, the Green's functions G t and G r can be combined: A multi-dimensional configuration (commonly referred to as a 'walker'), now carries both translational and angular coordinates. Sampling of the rotational and translational degrees of freedom are carried out by independent diffusive moves in these coordinates, as described by the factorization in Eq. (34). However, there will be coupling of the translations with rotations via the potential energy term V ( R, Ω) in Eq. (12). Details of the stochastic evaluation of Eq. (15) is described in Section 4. Use of Eq. (34) is usually referred to as "rigid body diffusion Monte Carlo" (often abbreviated by its acronym RBDMC), because of the implicit representation of the full N-body kinetic energy by a sum of rigid body rotational contributions together with the sum of the center of mass translational terms.

Time step considerations
The time step parameter used during the Monte Carlo simulation is determined by two factors. First, as pointed out by Suhm and Watts [18], when all other factors are equal, the time step is controlled by the smallest mass. Indeed, for light particles, the diffusion coefficient is large and the time step has to be small enough so that the particle does not move too far in a single step. But this is not the only factor. The effect of the potential can also be important. For example, strongly bound particles must be sampled with a smaller displacement, and hence with a smaller time step, to avoid the majority of Monte Carlo configurational moves ending in destruction of the configuration.
In the general case of a mixture of light and heavy particles, the time step has therefore to be chosen not only according to the smallest mass value, but also taking the effect of the potential into account. The presence of heavy, strongly bound particles may require the use of a shorter time step. For example, in the case of the water dimer, Gregory and Clary [44] use time steps of 0.5 -5.0 au for the all-atom calculations. Those small values are needed because of the presence of high frequency internal modes of the monomers. When using RBDMC, they were able to increase the time step up to 30 -60 au, because the introduction of the rigid body approximation eliminated these strongly confined degrees of freedom. This decrease in computational time step is generally cited as the primary motivation for an RBDMC calculation.
The situation with regard to time steps is somewhat more complicated for helium clusters, where a wide range of time scales is operative. For pure helium clusters, the time steps used vary typically from 1000 to 10000 au [41,45]. Introduction of a non-rotating, rigid SF 6 impurity requires a time step reduction to 250 -500 au [41,46] even though the mass of SF 6 is much bigger than that of He. This reduction of time step is due to the fact that the helium-SF 6 interaction potential is considerably stronger than the He-He potential: the larger time step will lead to inefficient sampling of the He-SF 6 relative motion. When the molecular rotation is included, the situation becomes quite complicated as one now has the additional intrinsic time scale of the molecular rotation. Thus one needs to compare the corresponding diffusion coefficients d p and D p for each species p included in the simulation, and to assess both the relative roles of the masses/moments of inertia (measured via d and D), and the confining role of the interactions in all degrees of freedom. For the systems studied in this paper, namely SF 6 and HCN in clusters of 4 He, the defining parameters are given in Table 1. The time steps required to eliminate any time step bias for these two systems are of order 20 -40 au for SF 6 , 10 -30 au for HCN) in implementations where attempted moves involve all particles being moved together. Somewhat larger time steps can be used when attempted moves involve only single particles (e.g., ∼ 50 au for SF 6 in 4 He n ). The latter is preferred for all except very diffuse systems. For molecules interacting with helium, in general all-atom potentials are not available, and this is a second reason why RBDMC is essential here. Table 1 Diffusion coefficients and potential characteristics for SF 6 and HCN in 4 He clusters. For SF 6 , V sad corresponds to an adiabatic barrier for a helium atom to move between sites. It corresponds to the minimum of the potential along the C 2 symmetry axis. For HCN, V sad is defined as the maximum of min r V (r, θ) for θ varying from 0 to π. It corresponds to the minimum of the potential for θ = π.

Rigid Body Diffusion Monte Carlo with Importance sampling
We consider first the conventional form of importance sampling in cartesian coordinate systems. This starts with the introduction of a distribution function F ( R, τ ) which is biased by a trial function Ψ T ( R): Taking its imaginary time derivative and using Eqs. (2) and (10) withT ≡T t (Eq. (13)), leads to where E L ( R) is the "local energy", given by and F ( R) is referred to as the "quantum force", given For simplicity, we have set D p = D for all p = 1...N. The solution to Eq. (36) is then where the translational importance sampling Green's functionG t is formally defined byG The short-time solution forG t is well known [39]. It is usually derived within a Fokker-Planck formulation [47]. We derive here an alternative solution using operators which provides a introduction to the subsequent new derivation of the full translational and rotational importance sampled Green's function,G. G t is initially cast in the position representation, The local energy contributions can then be immediately evaluated: We now replace ∇ 2 by the operator −P 2 /h 2 (P = −ih ∇) and use the identity (see Appendix B for proof), to obtaiñ Inserting a complete set of momentum states {| p }, with p representing the eigenvalues of the momentum operatorP , allows the last factor in Eq. (44) to be rewritten as where we have used Substituting k = p/h, and evaluating Eq. (46) using the standard integral leads to the well known result, Now we expand the discussion to include the rotational kinetic energyT r , and to incorporate an orientation dependent interaction potential, V ( R, Ω). The biased distribution function becomes Identifying the quantum forces for one-dimensional angular motion around each of the principal axes, α = x, y, z, for each monomer p, and the angular momentum operatorsL p in the p-th PAF, Eq. (17), leads to the generalization of Eq. (35) for a general rigid body Hamiltonian Eq. (10), withT =T t +T r and interaction potentialV ( R, Ω): This can be simplified using Eq. (43) and its analog for angular motions (Appendix B), to arrive at The short-time Green's function forH is put into position representation using Eqs. (34) and (20), and then making use of Eqs. (12) and (15): Inserting a complete set of one-dimensional angular momentum eigenstates, Eq. (23), in each of the terms inside the product on the right hand side as before (Section 2) leads to, e.g., for α = x, with d pα =h 2d pα as before. This can be shown to reduce to a Gaussian form in the small ∆τ limit, via the same procedure as employed in Section 2. Thus, replacing φ in Eq. (30) now by the angular displacement including the contribution from the angular quantum force, e.g., for α = x, one has We thereby arrive at the full, short-time factorized, limiting form for the translational and rotational importance sampled Green's function for a set of interacting rigid bodies: This may be generalized to multiple values of D p , by replacing the first factor with 1 4π∆τ Details of the implementation of Eq. (58) are presented in the next Section.
We note that, for the specific case of diatomic molecules, Lewerenz has proposed another method to implement rigid rotor motion, namely the method of adiabatic constraints [48]. This method is attractively simple and easy to implement for ground state studies of linear molecules [48][49][50], although it becomes considerably more involved when applied to more general rigid bodies [18]. Importance sampling of rotational degrees of freedom would be problematic within such an approach however, since it is not clear how to combine this with the adiabatic constraints. In particular, it is not obvious at what stage the quantum forces should be introduced, i.e., before or after imposition of the adiabatic constraint. Recently, Sarsa et al have shown that classical constraint dynamics can be used to generally impose bond length and bond angle constraints, in both biased and unbiased DMC [34].

Mixed frame and fixed frame implementations
In order to demonstrate Monte Carlo solution of Eq. (8) with Eqs. (34) and (58), we consider a system composed of m rigid rotors and n atoms (n + m = N). The Monte Carlo propagation over a step ∆τ requires then 3(m + n) translational diffusion moves, 3m rotational diffusion moves, and the local exponential growth/decay operations. The translational Green's function, Eq. (15), consists of a product of a diffusive term deriving fromT t , the translational kinetic energy, and a branching term deriving from the potential energy, V . Two kinds of moves are thereby made on an initially established ensemble of multi-dimensional configuration points or "walkers", Eq. (9). The first kind of move is a diffusion from R j to new positions R ′ j , implemented by sampling the vector R j − R ′ j from a 3N-dimensional Gaussian distribution, of width 2D p ∆τ in each of the 3N cartesian coordinates. The second move is a "branching" which modifies the weight of each walker in the ensemble by the exponential factor e −∆τ [(V ( R ′ )+V ( R))/2−E R ] . This can be implemented by i) replicating/destroying configurations, or by ii) maintaining and updating a continuous weighting factor associated with each configuration, or by iii) a combination of weights and replication [47,49]. The role of the reference energy E R now becomes clear, i.e., it provides a means to reduce the fluctuations deriving from the branching term.
Propagation withG t , Eq.(48), is performed by a simple modification of the unbiased propagation described above, namely that the diffusive moves are modified to include the contribution from the quantum force F ( R). The coordinates of each particle are therefore updated according to, e.g., for R px , where ∆ px is a random number sampled from a Gaussian of width 2D p ∆τ . The associated modification of the branching term consists in using the local energy instead of the potential, as shown in the exponential factor e −∆τ [(E L ( R ′ )+E L ( R))/2−E R ] for the update of the weights and/or the replication of walkers.
The implementation for the rotational Green's function, unbiased G r and bi-asedG r , is similar. In the unbiased formalism, diffusive rotational moves consist in sampling the three rotation angles for each rotor from a 3m-dimensional Gaussian distribution of width 2d pα ∆τ . In the biased version, these rotation angles also include the quantum force : For the full Green's function, both translational and rotational diffusion moves are made. We now give a more detailed explanation of the various possible implementations of this in different frames of reference.

Mixed Frame
In the general case, i.e. for more than one rotor, it is necessary to use different frames for the rotation and the translation moves. The rotation must be performed in the principal axes frame, PAF of each rotor, the translations are done in the laboratory frame, SF. The evaluation of the potential and of the trial function have to be done in consistent frames. At each time step in the random walk, three rotations around the principal axes of the molecule are made. During the time step τ − ∆τ → τ , the PAF (PAF τ −∆τ ) is thus rotated into a new PAF: (PAF τ ). The rotation matrix describing the attempted move is :R whereR α (φ) represents a rotation of φ around the α-axis, and φ τ α denotes the angle moved around this axis in time τ . The corresponding relation between the coordinates in the new PAF and in the old can be derived using the following identity:R Here y ′ (z ′′ ) is the new orientation of the y (z) axis after the first (the first and second) rotation(s). This identity is easily verified by transforming each of the rotations appropriately [51], e.g.
The matrix representation of Eq.(63) is given by: If the initial orientation of the PAF is given by 3 Euler angles (ϕ 0 , ϑ 0 , χ 0 ) with ϕ 0 ∈ [0, π), ϑ 0 ∈ [0, π] and χ 0 ∈ [0, π), Eq. (66) can be used to compute the relation between the coordinates in the PAF and in the SF. Then In an unbiased "pure" DMC calculation, φ x , φ y and φ z are sampled from a Gaussian distribution which is centered at zero and has width √ 2d α ∆τ . In an importance sampled calculation, we have to add the quantum force, according to Eq. (61). Rotation around x modifies the orientation of the y principal axis of the rotor. The second rotation is thus made around the new orientation of this axis, followed by the third rotation around the subsequent new orientation of the z principal axis. Due to the small values of these angles, we neglect the non-commutation of these 3 rotations.
In the Mixed Frame formulation, the SF cartesian coordinates for the centers of mass of all particles are used at every time step in the scheme. The use of the P matrix allows us to compute the coordinates of the relative vectors between the atoms and the rotors, in the PAF of each rotor. These relative vectors are needed for evaluation of the potential term and of the trial wave function.
As noted above, the Mixed Frame formulation is needed whenever the system contains more than one rotor. It is also required whenever the Euler angles defining the orientation of the PAF in the SF are needed. These angles are defined by the P matrix. Knowledge of these angles is necessary for calculation of excited states within a fixed node approximation whenever the nodes are imposed in a SF frame, because the trial function is then an explicit function of these angles. In the particular case of a trial function that depends only on the second Euler angle , e.g. Ψ T = |1, 0, 0 ∝ cos ϑ = P 33 , computation and storage of only the last column of P is sufficient. Section 5.3 provides an illustrative example of calculation of internal bending excitations for He-HCN using this Mixed Frame formulation.

Fixed Molecular Frame
When the system under consideration contains only one rotor, the Monte Carlo propagation can be performed in the PAF of this rotor. Like in the previous section, at each time step in the random walk, three rotations around the principal axes of the molecule are made. The corresponding relations between the coordinates are given by equation (66).
In order to perform the translational moves, one continuously follows the molecular orientation. Then at each time step (τ ), the atoms and the center of mass of the rotor are defined by their coordinates in the molecular PAF τ . The matrix R (τ ) is used to compute the coordinates in the new PAF using their coordinates in the former one (PAF τ −∆τ ). The translational moves are then made according to Eq. (60), using the coordinates in the old PAF. Vectors describing the position of the atoms in the two PAF's are used to compute the atom-rotor interaction potential and the trial wave function, before and after the step.
Calculations for SF 6 -He system using the implementation described in this section leads to the same results as those obtained with the Fixed Laboratory Frame. One disadvantage of this Fixed Molecular Frame scheme is that it does not allow for easy computation of the Euler angles defining the orientation of the PAF with respect to the SF frame.

Fixed Laboratory Frame
In the case of a spherical top molecule, Eq. (17) is applicable in the space-fixed laboratory reference frame (SF) as well as in the PAF, because any choice of orthogonal frame diagonalizes the inertia tensor. The rotational moves can therefore be made around the axes of the SF frame (X, Y, Z). Although the rotation moves can be made without knowing the orientation of the rotor at a given time, it is necessary to keep track of the molecular orientation in the SF frame because computation of the interaction with the other rotors and atoms (Eq. (11)) requires knowledge of the orientation of the preassigned reference directions of the spherical top (e.g. the S-F bonds in the case of SF 6 He n clusters).
Rotational moves, performed according to Eq. (63), are combined with translational moves of the n atoms and of the m center of masses of the rotors, made in the cartesian coordinates of this Laboratory Frame (SF), within the conventional approach (Section 4 above). This fixed Laboratory Frame scheme has recently been implemented in a study of rotational excitations of SF 6 in helium clusters [17], and is also used in the SF 6 applications presented here.

Angular Sampling in Three Dimensions: SF 6 -He
The SF 6 He n system with n = 1 was used to test the rotational sampling algorithm. In particular, since for a spherical top we know the exact threedimensional rotational Green's function (Eq. (33) and Appendix A), we can check the accuracy of the short-time decomposition of this into three onedimensional rotational factors, Eq. (32). The calculations are carried out in the Fixed Laboratory Frame as described in section 4.3, and the SF 6 -He interaction is modeled using the anisotropic potential of Pack et al [52]. We Table 2 Unbiased DMC energies in K for the SF 6 -He system. The numbers in parentheses give one standard deviation statistical errors in the last digits. * The FBR/DVR value is converged, i.e., it is stable when the unsymmetrized basis set size is increased from 19890 to 33930.
FBR/DVR DMC DMC (no rotations) -37.14 * -37.4(2) -38.4 (3) implemented both sampling schemes for this system, i.e., one rotation about a randomly chosen axis, corresponding to Eq. (33), and three rotations about the space fixed X, Y , and Z axes, corresponding to Eq. (32). Figure 1 (dotted lines) shows the evolution of the distribution of Euler angles (orientation) obtained from the direct sampling of the three-dimensional rotational Green's function. It also shows (solid lines) the evolution of the Euler angle distributions obtained with the factorization into three independent rotations. Both of these methods yield similar distributions that exhibit the expected Gaussian form (see Appendix A). The computed ground state energies are summarized in Table 2. The value obtained with rotation incorporated is in good agreement with result obtained using basis set expansion methods. These reference calculations employed the Finite Basis Representation -Discrete Variable Representation (FBR-DVR) collocation scheme of Leforestier using a Wigner basis set [53]. We also show the energy obtained when the molecular rotation is omitted. This is ∼ 1 K lower, and is in good agreement with previous DMC studies made without rotation [41]. Since the energetic effect of rotation is very small for n = 1, i.e., the rotational energy increment in the SF 6 -He ground state is only ∼ 1 K, we prefer to investigate the accuracy of the two different short time Green's function approximations using an artificial SF 6 having rotational constant B 0 = 0.91 cm −1 , i.e., 10 times the true gas phase value. This increases the variation of the energy due to rotation and amplifies the effect of any approximations. Rotating by performing a single rotation about a randomly chosen axis yields a ground state energy (E 1rot = -33.9(6) K) that is consistent with both the value obtained by implementing the three rotations about fixed axes scheme (E 3rot = -33.8(2) K) and a FBR-DVR reference calculation (E F BR = -33.6 K). In contrast, incorrect sampling (omitting the a = 3 factor discussed in Appendix A) yields erroneous energies which are too low (E = -35.8 (7)). This omission is operationally equivalent to reducing the rotational constant by a factor of 3.

Importance sampling to overcome artificial dissociation in medium-sized weakly bound clusters
The necessity of importance sampling appears already in ground state studies of moderately sized, weakly bound clusters, e.g., SF 6 He n and HCN He n with n > 12. Both of these systems require importance sampling to avoid the unphysical dissociation of helium atoms which can result from the occasional diffusion of He out to large distances from the cluster. In unbiased DMC, once an atom diffuses very far from the weakly bound cluster it has only a very small probability to return diffusively to the region of stronger binding, while the energetic penalty for motion further away from the molecule becomes negligible. Moreover the configuration space for a motion toward the molecule is much smaller than the one for a motion away. This is therefore a metastable situation from which it is very difficult to return to configurations of strong binding. The result is that helium atoms become "lost" outside the boundaries of the simulation box, leading to false dissociation as well as to artificially high energies.
One of the attractive features of unbiased DMC is that it provides the fulldimensional wave function. This may be fit to obtain trial functions for importance sampling, as has been done by a number of authors [32,54,55]. For the molecule-helium clusters with n = 1, the helium wave function may be squared to obtain the density. However, for n > 1, three-dimensional representation of the n-helium wave function is not possible. In order to nevertheless provide a visual representation of the wave function, we shall extract an effective one-particle wave function by projection, binning independently the position of the n heliums in the molecular frame.
For the spherical top SF 6 molecule, these calculations are performed in the space-fixed, laboratory frame (SF), described in the previous section. For the linear HCN molecule, calculations are performed in the Mixed Frame scheme of Section 4.1. Rotation around the z-axis of the PAF is not allowed (B z = 0), and consequently, the rotation formalism presented there is reduced to two angles.
Simple trial wave functions can be used to avoid artificial dissociation. For one rigid molecule, a form of trial wave function possessing the correct permutation symmetry is given by the product of pair correlation terms [22]: In the examples presented here for N ≤ 20, the helium-helium correlation is relatively unimportant compared to the helium-molecule term, so for simplicity we have used Ψ He−He (R p , R q ) = 1. For SF 6 He 1 and SF 6 He 20 , the helium atoms lie within the first solvation shell. As long as we are interested primarily in the binding energies here, and not concerned with detailed structural information in the angular degrees of freedom, we can adequately describe the He-SF 6 correlations using simple radial wavefunctions, of the form where r pX = |R p − R X |, and X=SF 6 . For this demonstration with n = 20, we used a = 0.80 and c = 42000 in atomic units. These parameters were obtained from fitting results of unbiased simulations for similar n values. Care should be taken not to choose a trial wavefunction that is too restrictive (i.e. negligible amplitude in regions that may actually be important). Such choices of wavefunction can give rise to instabilities and to erroneous energies. If the radial trial function is too narrow, falling off too sharply in the repulsive regime close to the impurity, the convergence of the DMC scheme becomes more difficult and the energy can become stuck below the true value. We found this to be the case when a trial function fit to unbiased wave functions obtained for small n (n = 1 with a = 1.5, c = 6000 a.u.) was employed for importance sampling of a cluster with significantly larger n (n = 20). This emphasizes the importance of choosing a trial function which does not impose excessive restrictions on the sampling.
For HCN, we used the radial part obtained from the high quality anisotropic He-HCN trial function derived for use in the fixed node calculations of Section V.C. This is of the form where r pX = |R p − R X | with X=HCN, and parameters presented in Table 3. Table 3 Parameters (in au) used for the definition of the HCN-He trial wavefunction. The effective one-particle helium wave functions obtained by unbiased RB-DMC are shown in Figure 2a for SF 6 He 20 , and in Figure 3a for HCN He 20 . The "wings" at the outer edge of the binning domains are the signature of dissociating helium atoms. They appear here because we put all walkers having a position larger than the domain definition into the last bin at the edge of the domain. The average number of artificially dissociated atoms defined as those lying at a distance greater than 20 a.u. for both systems, is n diss ∼ 5 for HCN, and n diss ∼ 1 for SF 6 .  6 He 20 obtained using unbiased DMC and projecting the helium wave function into the molecular frame. This is arbitrarily normalized as if it were a density. (b) The mixed helium density < Ψ T |ρ|Ψ > for SF 6 He 20 obtained using a radial trial function (see text).
This dissociation phenomenon also has a pronounced effect on the energy. For clusters containing smaller numbers of helium atoms, unbiased RBDMC and importance sampled DMC lead to the same energies, as they should if both calculations are converged. For n = 20, convergence of the importance sampled DMC is straightforward, yielding a ground state energy of -237(4) K for HCN and -607 (11) K for SF 6 . However with unbiased DMC, in sharp contrast to this situation, full convergence of the energy is much more difficult -even impossible -for this number of helium atoms, because of the artificial dissociation of helium walkers, for both HCN and SF 6 . We can only obtain an upper estimate of the ground state energy as a result. For SF 6 this esti- Fig. 3. (a) One-particle effective ground state wave function for HCN 4 He 20 , obtained using unbiased DMC and projecting the helium wave function into the molecular frame. The HCN molecule lies on the z axis, with its center of mass at the origin and the H atom on the z > 0 side. The second axis (r) is the distance to that axis. The non-physical peaks at the edges of the box come from the binning procedure and show artificial dissociation of some helium atoms -see text. (b) Representation of the mixed density < Ψ T |ρ|Ψ > obtained using a radial trial function (see text).
mate is ∼ -564 K, and for HCN it is ∼ -196 K. In both cases the unbiased DMC value is considerably higher than the fully converged importance sampled value, reflecting the smaller number of truly bound helium atoms in the unbiased calculations. It is significant that this phenomenon is seen not only for the relatively weakly bound HCN molecule in 4 He n , but also for the much more strongly bound SF 6 molecule. Our analysis implies that this is a general phenomenon which will occur for any species beyond a minimum size which appears to be somewhat less than one solvation shell. Figures 2b and 3b show the mixed densities < Ψ T |ρ|Ψ > derived from the simple radial trial functions with importance sampled RBDMC, for SF 6 and for HCN, respectively. It is evident by comparison with Figures 2a and 3a that the use of a simple radial function which provides binding of the helium atoms to the embedded molecule very effectively removes the artificial dissociation problem, and allows us to fully converge a ground state energy calculation.
Although we do not evaluate structural quantities here other than this mixed density, we note that the use of importance sampling does allow a relatively straightforward access to computation of structural quantities using secondorder estimators [22]. Table 4 summarizes the energies and the number of artificially dissociated helium atoms for both clusters.

Excited States
Our third application addresses the computation of excited states for the He-HCN complex. Experiments in doped helium clusters show that the rotational spectrum of small molecules inside clusters of 4 He possesses the same symmetry as the corresponding spectra in the gas phase, but that the effective rotational constants are reduced [56,57]. For HCN this reduction is about 20%, [58,59] considerably less than the ∼ 80% reduction seen for the more strongly bound species such as SF 6 [60]. Table 1 shows that HCN is much more weakly bound to helium than SF 6 , but that in both cases the ground state energy of the n = 1 complex lies above the saddle point of the interaction potential. Analysis of such rotationally excited states within a fixed node approximation requires some estimate of the relevant nodal structure. Insight into this can be obtained from level assignments made according to simple models. The earlier study of SF 6 in 4 He N from our group calculated rotational levels corresponding to J ∼ j, where J is the total angular momentum, and j the molecular angular momentum [17]. For He-HCN, we evaluate here instead the excited state which has been assigned by Atkins and Hutson [61], and by Drucker et al [62], to the level |j = 1, ℓ = 1, J = 0 >. In this assignment scheme, J is the total rotational quantum number, j the quantum number of HCN, which is to a good approximation conserved in the weakly bound complex, and ℓ an "orbital" quantum number associated with the rotation of the helium around the HCN. For this excited state the fixed nodal approximation can be imposed in the molecular frame, as we describe below. Nevertheless, we employ here the Mixed Frame implementation of the rotational importance sampling algorithm, rather than the Fixed Molecular Frame version. For a linear molecule, while we must always perform the rotations in the molecular frame, nodal structures may be imposed in either the MF or the SF frame, leading to the natural choice of the Mixed Frame implementation as the most general in this case.
In order to obtain a trial wave function for the fixed node calculations, we first performed pure DMC runs which allow us to compute the ground state eigenfunction (see Figure 4). This was fit to an analytic form Ψ T (r, θ), where r and θ are the usual Jacobi coordinates of the He with respect to the center of mass and orientation of HCN, to obtain a high quality trial function. We use the following Legendre polynomial expansion: with {a i , c ki , k = 1, . . . , 5; i = 1, 2} the parameters of the fit presented in Table  3. Limiting the sums to ℓ ≤ 3 and i ≤ 2 does not noticeably affect the precision of the fit.
Such a high quality trial function is necessary for the fixed node calculations of the He-HCN complex described here. In order to construct a trial excited state function for a fixed node calculation, we combined this high quality ground state trial function with a simple function that imposes the nodal constraint. This is chosen here by comparison with the exact excited state function, which was obtained by repeating the collocation calculations of Aktins and Hutson [61] to compute and analyze the wavefunction for this level. The energy of this level is −8.5 K, and a contour plot of the corresponding wave function is shown on figure 5. The structure of this function implies a strong internal bending character to this excitation. The excited state function possesses a single nodal surface, which is seen to be approximately r-independent. This motivated us to employ a fixed nodal surface defined by cos(θ + χ), where χ is a parameter and θ is the internal Jacobi angle of the cluster i.e. the trial function used is thusΨ T (r, θ) = Ψ T (r, θ) cos(θ + χ). Since the calculation is carried out in the mixed frame formulation, this trial function is implicitly dependent upon the orientation of the rigid body in the space fixed frame, even though the nodal structure is imposed in the molecular frame. Comparison of the energy obtained for calculations restricted to each of the two sides of the nodal surface allows us to then obtain the optimal value of the parameter χ. This is arrived at when these two energies are equal. A similar procedure was used by McCoy and co-workers in an interactive sense to optimize nodal surfaces [32]. For both the ground state and excited state computation, we use an ensemble of approximately 5000-6000 walkers. We perform a block averaging in order to reduce the correlation between successive steps. After equilibration of the ensemble, we perform 1 run of 800 blocks, each of which consists of 150 time steps with ∆τ = 10 au. We reject all moves corresponding to a crossing of the nodal surface, keeping the corresponding walker at its former position [39].
The optimal value of χ was found to be χ = 15.65 degrees, which yielded a common energy of -8.4(2) K on both sides of the resulting nodal surface. This fixed node value is in excellent agreement with the value -8.48797 K (-5.89897 cm −1 ) obtained from the collocation method in Ref. [61], especially when the approximate nature of the nodal constraint is taken into account.
When studying clusters with larger numbers of helium atoms using unbiased RBDMC, we noticed that the θ dependence apparent in Figure 4 tends to smooth out [63]. This suggests that for larger clusters, it might be adequate for some purposes to reduce the trial function to the radial term only, e.g., for studying ground state energies (Section V.B).

Conclusions
We have derived an algorithm for Rigid Body Monte Carlo with importance sampling for all degrees of freedom, including both translation and rotation, and provided a rigorous derivation of the short time rotational Green's function and its associated quantum forces. Three possible implementations of the algorithm were presented, which differ according to the combination of frames used for rotational and translational moves during the propagation. The most general implementation scheme is the mixed frame implementation in which translational moves are made in the laboratory reference frame, while rotational moves of all rigid bodies are made in the principal axis frame of each monomer. This allows full importance sampling calculations of molecular clusters to be performed.
We have presented several applications of this importance sampled Rigid Body Diffusion Monte Carlo. First, we demonstrated the correctness of the shorttime factorization of the rotational Green's function by comparison of the resulting product of one-dimensional rotations with the well known threedimensional formulation for a spherical top. Second, we showed that impor-tance sampling of translational degrees of freedom is essential to avoid nonphysical dissociation of weakly bound helium atoms in doped quantum clusters of helium when more than a few helium atoms are present. Third, we showed how excited states can be easily accessed with this algorithm using the fixed node approximation and a trial wavefunction with an implicit dependence on the orientation of the rigid body. The same algorithm may now be applied to study excitations beyond the fixed node approximation, e.g., using the POITSE approach [31].

Acknowledgements
Financial and computational support from the National Science Foundation through NSF grants CHE-9318737, CHE-9616615, and the NSF NPACI program is gratefully acknowledged. We thank Kevin Schmidt for pointing out the need for the factor of three in the Gaussian sampling of angles for a spherical top (see Appendix A). We thank Claude Leforestier for permission to use his FBR-DVR Wigner Basis code [53].

Sampling of the Spherical Top Green's Function
In order to provide a sampling of the operator where d ≡ 1/2I, according to the configuration representation of Eq. (33), one must choose both an arbitrary axis of rotation, and an angle of rotation. In order to implement this, it is necessary to first specify the frame of reference of the operatorL, i.e., its axis of quantization for L z . This is evident on expansion of the operator in the angle space, whereΩ is a unit vector in an arbitrary direction, and a is a constant to be determined. Expanding both sides of Eq. (73) to order (∆τ ) 2 , we find 1 − dL 2 ∆τ = 1 − a∆τ 4π dΩ( L ·Ω) 2 .
The integral can be evaluated once the direction of quantization of L is established. Choosing this for simplicity asẑ, so that L ·Ω = L cos θ, and evaluating the integral over the spherical angular coordinates dΩ ≡ dφdθ, leads to the equality from which we conclude that a = 3. Sampling of Eq. (73) can then proceed by a) choosing a random vectorΩ, and then b) sampling the angle of rotation by choosing a value from a gaussian of width 3 × (2d∆τ ). This procedure provides access to the full angular space in Eq. (73). When we are dealing with small rotations, this approach gives rise to similar orientational distributions to the scheme in which 3 small rotations are performed about the 3 cartesian axes (see Figures 1). This similarity will continue to hold as long as the rotation ζ is small enough such that the single rotation about the arbitrary axisn ζ can be approximately factored into a product of three cartesian rotations, i.e., exp (iζĴ ·n ζ ) ≃ exp (iζJ x sin θ cos φ) exp (iζJ y sin θ sin φ) exp (iζJ z cos θ).(76) Then the corresponding one-dimensional rotation around, e.g., thex axis, is made with the small angle ζ sin θ cos φ.

Commutation Relations
Expanding the standard commutator we arrive at Generalizing the one-dimensional function f (x) to the three-dimensional F ( R) = f x i + f y j + f z k, and evaluating Eq. (79) for each degree of freedom, leads then to the desired relation

Computational Scheme
Start with an initial ensemble of walkers (e.g. distributed according to the trial wavefunction). Compute the local energy (also used as the starting reference energy) for the initial ensemble and then propagate the ensemble in imaginary time by repeating the following steps for each walker (configuration).
(1) Move each host and impurity atom according to equation (60). This move involves a Gaussian random number and a quantum force. (2) Each move is then accepted or rejected with acceptance probability given by min(1, W ( Q ′ , Q)) where This involves having to recompute the quantum force after the particle has been moved.
(3) Rotate the impurity atom according to equation (61) (4) Again, this move is accepted or rejected using the above criterion. (5) Determine effective time step, ∆τ ef f . See Ref. [39]. One has to average the effective time step over the various kinds of moves since in general, each kind of move has a different diffusion constant. where ∆ is a random number that is uniformly distributed over (0,1). If M = 0, this walker is destroyed, otherwise M − 1 duplicates of the configuration are added to the ensemble. (8) To maintain a stable population size (N(τ )), update the reference energy according to where the population control parameter, α is chosen to be as small as possible (to avoid biasing the results) but large enough that the population size stays within acceptable limits. N * is the desired population size and is often chosen to be the starting population N(0) or N(τ ).
Of course, many variants of the above procedure will also work. In the HCN studies, for example, all of the translation moves and rotations are performed simultaneously and then accepted or rejected as a single move. When studying excited states using the fixed-node approximation, one also adds the additional constraint that moves crossing the nodal surface are automatically rejected.
The statistical error is estimated using block averaging, where the propagation is split into N b blocks of N s steps where N s is longer than the correlation length [47]. In the literature, various methods are used to compute averages and errors from this blocked data: (1) Take 1 data point per block, and compute the average and standard error from these N b data points (2) Same as above but use the block average values for the N b data points (3) Use all the data points in computing the mean. The standard error for each block, σ b , is computed from the N s values within each block. These are averaged to determine the reported error: In our examples, we report the largest of the above estimates (usually close in value).