Reduced-variance orientational distribution functions from torque sampling

We introduce a method to sample the orientational distribution function in computer simulations. The method is based on the exact torque balance equation for classical many-body systems of interacting anisotropic particles in equilibrium. Instead of the traditional counting of events, we reconstruct the orientational distribution function via an orientational integral of the torque acting on the particles. We test the torque sampling method in two- and three-dimensions, using both Langevin dynamics and overdamped Brownian dynamics, and with two interparticle interaction potentials. In all cases the torque sampling method produces profiles of the orientational distribution function with better accuracy than those obtained with the traditional counting method. The accuracy of the torque sampling method is independent of the bin size, and hence it is possible to resolve the orientational distribution function with arbitrarily small angular resolutions.


I. INTRODUCTION
The spatial and orientational order in classical equilibrium many-body systems is the result of a delicate balance between forces and torques of internal, entropic (diffusive), and external origin.One-body distribution functions, obtained as statistical averages resolved in either space, orientation or both of these, are essential for the description and understanding of the organization of many-body systems at the microscopic level.For example, the density profile, which is an average over a statistical ensemble of the number of particles at a given position, provides information about the spatial structure of the many-body system.Traditionally, the density profile in computer simulations has been obtained by discretizing the simulation box and counting the number of particles in each element of the grid.Since the structure of the many-body system is the result of a force balance, an alternative to counting events in order to obtain the density profile consists of reconstructing it from the spatially resolved force contributions [1,2].The density profiles obtained via force-sampling methods have a reduced variance as compared to those obtained via the traditional counting method.Moreover, the density at a given position is constructed with information from the whole system.As a result the error in the density profile does not depend on the size of the elements of the grid [1,2].The density profile can therefore be resolved with arbitrarily high spatial resolution without increasing the computational cost.This is particularly useful for sampling two- [2] and three-dimensional [3] density profiles.Force-based estimators can be also used to improve the sampling of the radial distribution function [1,4,5], and that of the correlation functions required in the Green-Kubo expressions relevant for mobility profiles [6].
It is interesting to note that force-sampling methods can be derived from the general and versatile mapped averaging framework [7][8][9][10][11], in which approximate theo- * delasheras.daniel@gmail.com;www.danieldelasheras.comretical results are used to reformulate an ensemble average with reduced variance.Reduced-variance estimators were first introduced in classical and quantum Monte Carlo simulations [12,13].An account of reducedvariance estimators that make use of force sampling methods is given in a recent review [14].
Beyond constructing statistical estimators with low variance in equilibrium systems, the internal force can be used to derive force-based density functional theories [15,16], and it plays a fundamental role in the construction of exact sum rules using the symmetries of the system [17][18][19].Moreover, the use of the thermodynamic force can also improve the accuracy of adaptive resolution schemes [20] in which the simulation box is split in regions that can be treated with different levels of resolution [21][22][23].Another potential application of force sampling methods is to improve the convergence of Kirkwood-Buff [24] integrals in molecular simulations [25].
Moreover, the knowledge of the internal force field is not only beneficial in equilibrium systems.The adiabatic approximation, which substitutes the non-equilibrium internal forces by those in an equilibrium system, is at the core of popular dynamical theories such as dynamic density functional theory (DDFT) [26][27][28][29][30]. Sampling the internal forces in many-body non-equilibrium simulations and comparing them to those in equilibrium systems is therefore crucial to develop and test the accuracy of dynamical theories that go beyond the adiabatic approximation such as superadiabatic-DDFT [31] and power functional [32][33][34] theories.Knowledge of the nonequilibrium internal forces facilitates also the construction of the external force field that generates a desired dynamical response via custom flow methods [35,36], and serves to gain insight into physical processes such as the occurrence of viscous forces generated by the acceleration field [37].
In systems with translational and rotational degrees of freedom, such as liquid crystals, it is not only the forces but also the torques that are crucial in the determination of the equilibrium and non-equilibrium properties of the many-body system.The force balance equation is complemented and coupled with a torque balance equation.
Together, the force and the torque balance equations determine in equilibrium the positional and the orientational order of the system.
Here, we demonstrate that torque sampling, i.e. the analogue to force sampling in systems with orientational degrees of freedom, significantly improves the sampling of the orientational distribution function in computer simulations as compared to traditional counting methods.As a proof of concept, we sample the torques using several differing types of dynamics (overdamped Brownian and Langevin dynamics), dimensionality (two-and threedimensional systems), interparticle interaction potential (rectangular and Gay-Berne particles), type of orientational order (uniaxial and tetratic), and overall density.In all cases, torque sampling outperforms the traditional counting method.

II. THEORY
We consider here classical systems of N identical interacting particles governed by either Langevin or overdamped Brownian dynamics.Exact one-body force and torque balance equations hold in equilibrium, and can be used to calculate one-body distribution functions from the forces and torques acting in the system.We start by revisiting the force balance equation in a many-body system with only translational degrees of freedom.

A. Force balance equation for isotropic particles
In many-body systems with only translational degrees of freedom, such as a system of isotropic particles (e.g. a fluid of Lennard-Jones particles), the exact one-body force density balance equation in equilibrium reads [38] The first term on the right hand side of Eq. ( 1) stems from the (ideal gas) diffusion, with k B being the Boltzmann constant, T is absolute temperature, ∇ is the derivative with respect to the spatial coordinate r, and ρ(r) is the one-body density distribution which is given by where the angles denote a statistical average over an equilibrium ensemble, δ(•) is the Dirac distribution, r i is the position of particle i, and the sum runs over all the particles in the system.The second term on the right hand side of Eq. ( 1) is the force density profile, given by where f i is the sum of the internal and the external forces acting on particle i in microstate r N = r 1 . . .r N with N particles.That is where ∇ i is the derivative with respect to r i , and u(r N ) is the total interparticle interaction potential.In equilibrium, the (imposed) external force f ext (r) must be conservative and hence with V ext (r) an imposed external potential.The force profile follows directly from the force density profile via normalization with the density profile, i.e. f (r) = F(r)/ρ(r).
The sum of the ideal, internal, and external force densities vanishes everywhere in space since the system is in equilibrium.Otherwise there would be a net flow of particles.In equilibrium, the exact force density balance equation, Eq. ( 1), holds in systems following either Newtonian dynamics, Langevin dynamics, or overdamped Brownian dynamics.
Both, the density profile ρ(r) and the force density profile F(r), can be easily sampled in computer simulations via Eqs.( 2) and (3), respectively.Sampling ρ(r) via Eq.( 2) is the traditional method of counting of events of particle occurrences at space points.
The exact force density balance equation, Eq. ( 1), can also be used to calculate the density profile ρ(r) via the forces instead of the direct traditional counting method.Inverting Eq. (1) results in with ρ 0 a constant and ∇ −1 the inverse ∇ operator.In effectively one-dimensional systems (e.g.planar geometry), the profiles depend only on one space coordinate and hence the ∇ −1 operator reduces to a simple spatial integral.Different approaches can be used to solve Eq. ( 6) in more general geometries [2].The unknown integration constant ρ 0 in Eq. ( 6) can be determined via normalization of the density where the integral is over the whole system volume.Results for the density profile calculated via force sampling, Eq. ( 6), carry a statistical uncertainty smaller than that of the standard counting method [2] since (i) force sampling avoids the inherent ideal gas fluctuations, and (ii) uses non-local information, the forces in the whole system, to determine the density profile at each space point.

B. Torque balance equation for anisotropic particles
For anisotropic particles, the one-body density distribution depends not only on the space coordinates r but also on the orientation, which is denoted here by the unit vector û: In addition to the exact equilibrium one-body force density balance equation, there exists an exact one-body torque density balance equation: Here, F(r, û) and T(r, û) are the force density and the torque density, respectively.Both, F(r, û) and T(r, û), contain external and internal (inter-particle) contributions and they depend in general on position and orientation.As before, ∇ is the gradient operator acting on the position, and R is the orientational counterpart operator acting on the orientation û, i.e.
with ∇ û the derivative with respect to the Cartesian coordinates of û.
The one-body torque density is accessible in computer simulations via (12) with ûN = û1 . . .ûN and ûi = (sin θ i cos ϕ i , sin θ i sin ϕ i , cos θ i ) being the orientation of particle i.Here, θ i and ϕ i are the polar and azimuthal angles of particle i, respectively.The torque on particle i is t i , given by (13) with Ri = ûi × ∇ ûi .Note that both the total interparticle potential u(r N , ûN ) and the external potential V ext (r, û) are allowed to carry a dependence on the particle orientation.The one-body torque density is therefore the sum of internal and external contributions with Using eq. ( 8) the external contribution is simply Further details regarding the derivation of the onebody torque density balance in equilibrium are given in Appendix A.
In general, the force and the torque density balance equations are linked via the one-body density distribution.Here, we focus only on the role of the torque balance equation.For this we consider in what follows systems that are homogeneous in space and therefore cases in which the force balance equation does not play any role.In such systems ρ(r, û) = ρ b f (û), with ρ b being the bulk density, and f (û) being the orientational distribution function.That is, f (û)dû is the probability of finding a particle with orientation û within a solid angle dû.The orientational distribution function is therefore normalized such that Using the traditional sampling method, the orientational distribution function can be sampled in computer simulations as The prefactor 1/N ensures the proper normalization of the orientational distribution function.
For spatially homogeneous systems, the one-body torque density balance equation (10) Isolating the ideal gas term of Eq. ( 21) and integrating appropriately we obtain an expression for the orientational distribution function Here, R−1 is formally the inverse operator of R and f 0 is an integration constant that ensures the proper normalization of the orientational distribution function.For two-dimensional systems and uniaxial three-dimensional systems, the inverse operator R−1 reduces to a simple angular integral, see Appendix B.
Obtaining the orientational distribution function via the one-body torque density balance has the advantage of treating the ideal gas part explicitly and hence, it avoids the corresponding fluctuations present in the counting method.The only source of statistical inaccuracies is in the sampled one-body torque density which is integrated over in order to obtain the orientational distribution function.As it turns out, this process reduces the statistical noise significantly.Orientational distribution function sampled with the counting method (orange) and the torque sampling method (blue) for three different sampling times: (a) 10τ , (b) 10 3 τ , and (c) 10 5 τ .The black dashed line is the "true" equilibrium profile obtained by sampling over 10 7 τ and taking the arithmetic mean of the counting and the torque sampling methods.The bin size is 10 −4 π.The inset in (c) is a close view of the encircled region.(d) A characteristic snapshot of the system: N = 64 particles with rectangular shape subject to a weak external potential that orients the particles along the x−axis.A log − log plot of the error ∆ as a function of the sampling time (e) and the error ∆ as a function of the bin size (f) using the counting (orange squares) and the torque sampling (blue circles) method.In panel (e) the bin size is fixed to 10 −4 π and the dashed lines are linear fits.In panel (f) the sampling time is fixed to 10 2 τ and the solid lines are guides for the eye.

III. RESULTS
As a proof of concept, we test the validity of the torque sampling method with two different systems: (i) twodimensional rectangular particles following Langevin dynamics and (ii) three-dimensional Gay-Berne particles following overdamped Brownian dynamics.

A. Two-dimensional system of rectangular particles
We consider a two-dimensional system of particles with rectangular shape undergoing Langevin dynamics (implemented according to Ref. [39]).The interaction between two particles is modeled via a purely repulsive potential φ(r) = σ r 12 .Here, r is the minimum distance between the two particles, σ is our length scale, and is our energy scale.The potential acts only between the two closest points (one on each particle) located on the particles' perimeter.The interparticle potential generates both an internal force and an internal torque.Details about the calculation of the forces and the torques, as well as about the integration of the equations of motion are given in the Appendix C We study a system of N = 64 rectangular particles with length L/σ = 10 and width D/σ = 2 in a square box of length 100σ and periodic boundary conditions.We set the temperature to k B T / = 1 and the integration time step to ∆t/τ = 10 −3 with τ = σ m/ and m the mass of one particle.We sample every 10 ∆t.Since the system is very diluted, the equilibrium bulk state is isotropic.We induce orientational order via the external potential V ext (ϕ)/ = −0.5 cos 2 ϕ, with the angle ϕ measured anticlockwise with respect to the x−axis.A characteristic snapshot of the system is shown in Fig. 1.
We initialize the particles randomly and equilibrate the system with a simulation lasting 10 3 τ .After equilibration we sample the orientational distribution function via the counting and the torque sampling methods.The results are shown in Fig. 1 for three different sampling times: 10τ panel (a), 10 3 τ panel (b), and 10 5 τ panel (c).Due to the head-tail symmetry of the particles we represent the orientational distribution function in the interval ϕ ∈ [0, π] only.Torque sampling provides at each time a profile which is closer to the "true" equilibrium profile than the one provided by the counting method.The "true" equilibrium profile f eq (ϕ) is defined here as the arithmetic mean of the profiles obtained with the counting and the torque sampling methods in a long simulation (total simulation time 10 7 τ ).For all sampling times the statistical noise in the profiles using the counting method is significantly larger than that using the torque sampling method.To quantify the accuracy of each method, we define an error parameter as the integrated square difference between the "true" equilibrium profile and the sampled profile Here, f s is the profile sampled using the counting or the torque sampling methods.As can be seen in Fig. 1(e) the error of the torque sampling method is for all sampling times below the error of the counting method.For this particular bin size (10 −4 π) one has to sample about 10 times longer using the counting method than using the torque sampling method to reach the same accuracy.In Fig. 1(f) we investigate the effect of varying the bin size at a fixed sampling time (10 2 τ ).By decreasing the bin size we increase the level of detail with which we resolve the orientational distribution function.However, decreasing the bin size obviously increases the number of bins and, as a direct consequence, the error in the traditional counting method also increases.Note that in the counting method the number of events that contribute to each bin is proportional to the bin size.On the other hand, the error in the torque sampling method is essentially independent of the bin size.The error does not increase by decreasing the bin size because the orientational distribution function is not determined by the local number of events.Instead, at each orientation the orientational distribution function is obtained via an orientational integral over the torque density.Analogue behaviour occurs also when sampling the density profile using the force sampling method in systems with only translational degrees of freedom [1,2,14].
Tetratic order.Instead of sampling the complete, angle-resolved, orientational distribution function, it is common to sample only a reduced set of orientational order parameters (moments of the distribution).However, having access to the complete orientational distribution function can help to fully understand the type of order in the system.To illustrate this, we investigate a densely packed system of N = 290 particles with length L/σ = 4 and width D/σ = 2 in a square box of length 75σ.The equilibration time was 10 4 τ .Due to their small length-towidth aspect ratio, the particles form in bulk at moderate densities a tetratic phase [40][41][42].In the tetratic phase the particles are equally likely oriented along two directions perpendicular to each other.We add an external potential of the form V ext (ϕ)/ = −0.5 sin 2 (ϕ − ϕ 0 ) with ϕ 0 /π = 1/4 and set the temperature to k B T / = 1.The external potential breaks the symmetry of the tetratic phase by favoring the orientation along the bottom-right to top-left diagonal of the square simulation box.
A snapshot of the system is shown in Fig. 2(a).The particles are colored according to their orientation.The resulting orientational distribution function is shown in Fig. 2(b) for a short sampling time of 1τ and in Fig. 2(c) for a sampling time of 10 3 τ .Clearly more particles are aligned along the bottom-right to top-left diagonal (ϕ/π = 0.75) than along the other diagonal (ϕ/π = 0.25) due to the external potential.In this example, the uniaxial order parameter or even the combination of both the uniaxial and the tetratic order parameters would not give enough information about the orientational order in the system.
The distributions sampled with torque sampling are always smoother than those sampled with the counting method.However, torque sampling sometimes produces artifacts for very short sampling times (of the order of 1τ ), like the negative values around ϕ/π = 0 shown in Fig. 2(b).It might be possible to eliminate these artifacts by either using a combination of linear estimators [43] or the mapped averaging framework [10].The artifacts are at least partially due to local angular currents orginated by fluctuations that do not vanish (on average) due to the short sampling times.The occurrence of these angular currents is apparent when comparing the orientational distribution functions sampled at short, Fig. 2(b), and long, Fig. 2(c), sampling times (Cf. the evolution of the value of the orientational distribution functions at the peaks).For longer sampling times, Fig. 2(c), the angular current averages to zero for all orientations, and the distribution function calculated with torque sampling is free of artifacts.The profile obtained with torque sampling is more precise than that obtained via counting.Even at very long sampling times, e.g. 10 5 τ in 2(c), torque sampling outperforms counting.This is particularly clear when looking at the numerical angular derivative of the distribution function, see inset of Fig. 2(c).
Sampling the torques is not only useful to improve the sampling of the orientational distribution function but it also helps to understand the underlying physics.As an illustration, we show in Fig. 3 the components of the torque balance equation in the system with tetratic ordering and an external potential.The torques point along the z−direction.That is, positive (negative) torques tend to rotate the particles anticlockwise (clockwise), increasing (decreasing) therefore the value of ϕ.The diffusive torque (blue) always favors an isotropic state by trying to remove the inhomogeneities in the orientational distribution function.In the current configuration, the diffusive torque tries to orient the particles away from the diagonals.The behaviour of the internal torque depends on several factors such as the interparticle potential, the temperature, and the density.In the current example, the internal torque (red) favors tetratic ordering by trying to align the particles along the diagonals.The imposed external torque (yellow) tries to orient the particles along the bottom-right to top-left diagonal (ϕ/π = 0.75) and it also tries to orient the particles away from the other diagonal at ϕ/π = 0.25.As dictated by the torque balance equation, the sum of all three components (diffusive, internal, and external) vanishes since the system is in equilibrium, see Fig. 3 (black line).

B. Three-dimensional Gay-Berne fluid
We further test the method in a three-dimensional system of N = 500 Gay-Berne particles [44] confined in a box of size lengths L x /σ 0 = 10, L y /σ 0 = 10, and L z /σ 0 = 25 with periodic boundary conditions.We use the parameters σ 0 and 0 of the Gay-Berne potential as our length and energy scales, respectively.All details about the interparticle potential are presented in Appendix D. We set the length-to-width ratio of the particles to three.The particles follow overdamped Brownian dynamics.Time is measured in units of τ 0 = γσ 2 0 / 0 , with γ the translational friction coefficient against the implicit solvent.The particles are subject to an external potential V ext (θ)/ 0 = −0.5 cos 2 (θ), with θ the polar angle.Hence, the external potential favors uniaxial alignment of the particles along the z-axis.The temperature is set to k B T / 0 = 0.5.For details regarding the implementation of the overdamped Brownian dynamics see Appendix E The orientational distribution functions obtained via torque sampling and counting are shown for different sampling times in Fig. 4. Again, torque sampling provides profiles with better accuracy than counting.The differences between both methods are more acute for small values of the polar angle.The area of the bins on the unit sphere decreases close to the poles.Therefore, less events contribute to each bin, which produces large fluctuations of the profile obtained with the counting method.However, the profile obtained with torque sampling is unaffected by this problem since the error is independent of the bin size.

IV. CONCLUSION
Reduced-variance estimators can be constructed using force sampling methods [14] to measure e.g. the density profile and the radial distribution function in computer simulations with better accuracy than the traditional counting method.We have shown here that in equilibrium systems of interacting anisotropic particles, reduced-variance estimators can be also constructed via torque sampling.By sampling the torques and using the exact torque balance equation of equilibrium manybody systems, we have developed a method to accurately reconstruct the orientational distribution function.Although the cases that we have studied here are arguably toy models, they do cover a wide range of situations, including two-vs three-dimensional systems, dilute vs dense systems, uniaxial vs tetratic orientational order, and Langevin vs overdamped Brownian dynamics.In all cases, torque sampling has outperformed counting.
Force sampling works equally well in Brownian dynamics, molecular dynamics, and Monte Carlo simulations [2].Hence, although we have used here Brownian and Langevin dynamics, the torque sampling method is expected to also outperform the counting method in molecular dynamics and Monte Carlo simulations.
For small bin sizes, the statistical error for the counting method diverges, while the error for the torque sampling method is independent of the size of the bin.Hence, torque sampling can be particularly useful in cases where a large number of bins might be required such as for example when investigating biaxial nematics in threedimensional systems [45][46][47].
We have restricted our study to cases without positional order such that force and torque balance equations are decoupled.There exist several fully inhomogeneous standard situations accessible in computer simulations [48][49][50][51] in which both the density profile and the orientational distribution profile depend on the position coordinate, i.e. ρ(r, û) = ρ(r)f (r, û).These include, among others, the formation of stable bulk phases with both positional and orientational order [52][53][54], confinement [55][56][57], sedimentation [58][59][60], formation of topological defects [61][62][63][64], and nucleation [65] in liquid crystals.The force balance equation and the torque balance equation are then coupled and jointly determine the spatial and the orientational order of the system.The combination of force and torque sampling should be in such cases substantially better than counting which requires filling a multidimensional histogram in both positions and orientations.
The formulation of the torque sampling method presented here cannot be directly applied to hard particle models [66], in which forces arise only due to particle collisions.However, the mapped averaging framework is applicable in hard particle models [67].Using the torque balance equation as input for the mapped aver-aging framework might result in reduced-variance estimators for the orientational distribution function.Exact sum rules involving the torques follow from the symmetries of the system [17] and might be also useful in the derivation of reduced-variance estimators in computer simulations of anisotropic particles.
The forces between individual colloidal particles are also accessible experimentally [68].It might therefore be possible to use force and torque sampling methods for the determination of distribution functions in experimental systems.

DATA AVAILABILITY
All the data supporting the findings are available from the corresponding author upon reasonable request.statistical averages • within the Fokker-Planck formalism can be computed as where the integrals cover the complete space of configurations.The velocity v i and the angular velocity ω i of particle i are given by Here, γ and γ r are the translational and rotational friction constants against the implicit solvent, respectively, and τ ext is an external torque field.The current J and the angular current J ω are Multiplying equations (A2) and by δ(r−r i )δ(û−û i ), summing over all particles i, and applying the average in Eq. (A1) yields directly γJ(r, û, t) = −k B T ∇ρ(r, û, t) + F(r, û, t), (A6) The above force (A6) and torque (A7) density balance equations hold in full non-equilibrium situations.In equilibrium, the equations simplify further since: (i) the time dependence drops from the density profile, the force density F, and the torque density T, (ii) the external force and the external torque are conservative, and (iii) both J and J ω vanish.Hence, in equilibrium Eqs.(A6) and (A7) reduce to Eqs. ( 9) and (10), respectively.The equilibrium force and torque balance equations do not change if the particles obey Langevin or molecular dynamics instead of Brownian dynamics but the derivation is slightly different.To derive the force and the torque balance equations in Langevin dynamics or molecular dynamics, one needs to time derive the current and the angular current, both of which also vanish in equilibrium: Here, the average • is again performed over the complete configuration space which in molecular dynamics includes integrals over the linear and angular momenta in addition to those over the positions and the orientations of the particles.Incorporating the time derivative inside the averages in Eqs.(A8) and (A9) results in the force and torque balance equations.In equilibrium, the integrals over the linear and the angular momenta can be carried out explicitly.
Appendix B: Torque sampling for single angular dependencies We derive here the expressions for the orientational distribution function as an angular integral over the torques in the The rotational operator can be written as with e ϕ and e θ being the unit vectors on the unit sphere in the azimuthal and in the polar directions, respectively.Two-dimensional system.In the two-dimensional system of rectangular particles, the orientational distribution function depends only on the azimuthal angle f = f (ϕ).Hence, using θ = π/2, the rotational operator, eq.(B1), simplifies to R = e z ∂/∂ϕ, and the torque density balance equation ( 21) is then with T also directed along the e z direction.The orientational distribution function can be then reconstructed with the sampled torques via with f 0 a normalization constant such that dϕf (ϕ) = 1.Using Eq. ( 17) it follows that which inserted in Eq. (B2) can be used to first solve the homogeneous equation analytically and then treat the internal torque density as an inhomogeneity.We find that both approaches give similar results.
Three-dimensional uniaxial system.We consider now the three dimensional system of Gay-Berne particles with a uniaxial distribution function, i.e. f = f (θ) and an external potential V ext (θ) that depends also only on the polar angle.Hence the rotational operator, Eq. (B1), is R = e ϕ ∂/∂θ, which inserted in the torque density balance equation yields where T = (T int (θ) + T ext (θ))e ϕ , from which we obtain The normalization constant f 0 is here such that dûf (θ) = 1.Again, it is possible to first analytically solve the homogeneous differential equation of (B5) by writing the external torques explicitly and treat the internal part as the inhomogeneous part.
Appendix C: Interparticle interaction between two rectangles and Langevin dynamics Two rectangles interact via a purely repulsive pairpotential of the form with r being the minimum distance between the two rectangles.Depending on the positions and the orientations of the particles, there are two possible scenarios, see Fig. 5: (i) the minimum distance is between a corner of one particle and a point located on an edge of the other particle, or (ii) the minimum distance is between two corners.We introduce a cut-off distance of r c = 2L + 3σ between the centers of mass of the particles.The potential generates a contact force between the two particles.The effect of the contact force between the two closest points is equivalent to apply both a force and a torque on the center of masses of the particles.The force acting on the center of mass of particle i due to particle j is given by and the torque acting on particle i due to particle j is given by Here, r i and ûi are the position of the center of mass and the orientation of particle i, respectively.The vector r c i joins the center of mass of particle i with the point of application of the force.
To calculate the minimum distance between two particles we calculate all possible corner-corner and corneredge distances and select the minimum of all of them.
Verlet-type integration algorithm for Langevin dynamics.We calculate the trajectories following the integration scheme for Langevin dynamics presented in Figure 5.The minimum distance between two rectangles r is either between a corner and an edge (a) or between two corners (b).The effect of the contact force acting on the points of minimum distance (top panels) is equivalent to the effect of the same force and a torque acting on the center of mass (bottom panels).
Ref. [39] for isotropic particles.The translational equations of motion for particle i are given by ṙi =v i , (C4) Here r i , v i , and f i are the position, the velocity, and the total force (internal plus external) of particle i (the overdot indicates time derivative), m is the mass of one particle, γ is the translational friction coefficient, and f rand i is a delta-correlated Gaussian random force (described in detail in Appendix E).These equations are integrated with the following Verlet-type scheme [39] r The equations of motion for the angular degrees of freedom in our two-dimensional system of rectangular particles following Langevin dynamics are given by φi = ω i , (C8) Here ϕ i , ω i , and t i are the azimuthal angle, the angular velocity, and the torque (internal and external) of particle i, γ r is the rotational friction coefficient, t rand i is a random torque (see details in Appendix E), and I = m(L 2 + D 2 )/12 is the moment of inertia around the axis normal to the particle that passes through the center of mass (the particles are assumed to have a homogeneous mass distribution).The torques and the angular velocity point along the z−direction (normal to the particles).Equations (C8) and (C9) have the same mathematical structure as Eqs.(C4) and (C5).We therefore apply the same integration scheme as for the positional degrees of freedom, replacing the mass m by the moment of inertia I and the translational friction γ by the rotational friction γ r .For simplicity we set γ r = γσ 2 .The value of the friction constants does not affect the equilibrium properties.
Appendix D: Gay-Berne potential We use the same implementation of the Gay-Berne potential as that in Ref. [44].The interaction potential between two particles is Here, ξ takes the values χ or χ and the parameters 0 and σ 0 set the energy and the length scales.We select a length-to-width ratio l = 3, and set the energy ratio between the side-by-side and the tip-to-tip configurations to d = 5.
Appendix E: Overdamped Brownian dynamics with orientational degree of freedom The equations of motion of particle i are Here, the angles denote an average of the noise, 1l is the identity matrix, ûû indicates the dyadic product, and δ ik is the Kronecker delta.The angular velocity is Using the vector triple product and the fact that ûi • ui = 0 due to ûi • ûi = 1, it follows directly that ui = ω i × ûi .(E6) Hence, the equations of motion can be integrated in time using the standard Euler algorithm via r i (t + ∆t) = r i (t) + ∆t γ −∇ i u(r N , ûN ) − ∇ i V ext (r i , ûi ) + η i (t) (E7) ûi (t + ∆t) = ûi (t) + ∆t γ r − Ri u r N , ûN − Ri V ext (r i , ûi ) × ûi (t) + Γ i (t).(E8) Here, η i is a Gaussian random displacement with zero mean and standard deviation 2k B T ∆t/γ, and is a random rotation.Here, U 1 i and U 2 i are Gaussian random numbers with zero mean and unit width, and ŵ1 i = e x × ûi /|e x × ûi |, ŵ2 i = ŵ1 i × ûi , and u i are orthonormal to each other.We renormalize ûi after each time step such that it remains a unit vector.
We arbitrarily relate the rotational friction coefficient to the translational friction coefficient via γ r = γσ 2 .Also, we use a single translational friction coefficient γ for displacements parallel and perpendicular to the main axis of the particle.The values of the friction coefficients do not play any role in the equilibrium distribution functions.The Euler integration time step is ∆t = 10 −4 τ and we sample every 10 −2 τ .Although we have used here a simple Euler scheme to integrate the equations of motion, it would be useful to extend the recently developed adaptive Brownian dynamics [70] algorithm to systems with orientational degrees of freedom to speed up the simulations.

Figure 1 .
Figure1.Orientational distribution function sampled with the counting method (orange) and the torque sampling method (blue) for three different sampling times: (a) 10τ , (b) 10 3 τ , and (c) 10 5 τ .The black dashed line is the "true" equilibrium profile obtained by sampling over 10 7 τ and taking the arithmetic mean of the counting and the torque sampling methods.The bin size is 10 −4 π.The inset in (c) is a close view of the encircled region.(d) A characteristic snapshot of the system: N = 64 particles with rectangular shape subject to a weak external potential that orients the particles along the x−axis.A log − log plot of the error ∆ as a function of the sampling time (e) and the error ∆ as a function of the bin size (f) using the counting (orange squares) and the torque sampling (blue circles) method.In panel (e) the bin size is fixed to 10 −4 π and the dashed lines are linear fits.In panel (f) the sampling time is fixed to 10 2 τ and the solid lines are guides for the eye.

Figure 2 .
Figure 2. (a) Characteristic snapshot of a Langevin dynamics simulation of N = 290 rectangular particles in a square box of side length L box /σ = 75 subject to an external potential that favors particle orientations along the bottom-right to top-left diagonal of the box.The particles are colored according to their orientation ϕ, measured with respect to the x−direction, see colorbar.Orientational distribution f (ϕ) obtained via counting (orange lines) and torque sampling (blue lines) using a bin size of 10 −3 π and for two sampling times: 1τ (b) and 10 5 τ (c).The inset in (c) shows the numerical angular derivative of the orientational distribution function f (ϕ) = df /dϕ using the central difference.

Figure 3 .
Figure 3. Components of the torque balance equation (normalized with the bulk density ρ b ) as a function of the angle in the tetratic configuration with an external field shown in Fig. 2. The torques point in the z-direction.Positive (negative) torques try to rotate the particles anticlockwise (clockwise), as indicated at selected angles by the color arrows.The external potential favors particle alignments along the bottom right to top left diagonal (ϕ/π = 0.75) of the simulation box.The bottom left to top right diagonal is located at ϕ/π = 0.25.Shown are the internal torque density (red), the diffusive torque density (blue), the external torque density (yellow), and the total torque density (black).

Figure 4 .
Figure 4. Overdamped Brownian dynamics simulation of N = 500 Gay-Berne particles in a three-dimensional box with periodic boundary conditions.Orientational distribution function f as a function of the polar angle θ obtained via counting (orange line) and torque sampling (blue line) using a bin size of ∆θ = 10 −3 π/2 for different sampling times: 10 3 τ0 (a), 10 4 τ0 (b), and 10 5 τ0 (c).The inset in panel (c) is a close view of the region of small polar angles as indicated.