Monte Carlo simulation of joint density of states in one-dimensional
Lebwohl-Lasher model using Wang-Landau algorithm
Monte Carlo simulation using the Wang-Landau algorithm has been performed in an one-dimensional Lebwohl-Lasher model. Both one-dimensional and two-dimensional random walks have been carried out. The results are compared with the exact results which are available for this model.
PACS: 61.30.-v, 64.70.Md
Keywords: Monte Carlo, joint density of states
The Wang-Landau (WL) approach [wl] in Monte Carlo (MC) simulation, introduced in 2001, has since been applied to different areas of statistical physics. While proposing the algorithm the authors have demonstrated the application of the algorithm to systems with discrete energy levels (like the Ising system, Potts model, spin-glass etc.). Over the past few years several authors [lj1, lj2, jaya] have reported the application of the WL method to systems with continuous energy spectrum. The common feature of these simulations is the expected discretization of the energy range which has been investigated i.e. the division of the energy range into a number of bins and the subsequent application of the ideas of Wang and Landau. It may however be noted that in majority of these simulations a random walk in energy space alone was conducted. Relatively few papers [2dwalk, cont] have so far appeared where a two dimensional random walk has been performed in a continuous system.
It has generally been noted that even for a modest system size, the simulation of joint density of states, where the variables for instance are energy and the order parameter, (both quantities being continuously variable), the computational time necessary for the implementation of the WL algorithm is very large. It may further be noted that till now no work has been published regarding the errors involved in the WL-simulation in such systems. Also, no simulation has so far been reported where comparison has been made with the exact results available (as is usually done for discrete systems with the two dimensional Ising model). The reason for this may of course be attributed to the non-availability of results of exact calculation for any non-trivial system having a continuous energy spectrum.
The present work has been intended to fill the gap in the availability of adequate information in this direction to some extent. The system we have chosen for this purpose is an one-dimensional array of three-dimensional spins (1, 3, where is the system dimensionality and is the spin dimensionality) interacting with nearest neighbours via a potential , where is the second Legendre polynomial and is the angle between the nearest neighbour spins and , (the coupling constant in the interaction has been set to unity). This model, known as the Lebwohl-Lasher (LL) model [ll], is the lattice version of the Maier-Saupe (MS) model [ms1, ms2, ms3] which describes a nematic liquid crystal in the mean field approximation.
The one-dimensional LL model (1, 3), which of course does not have any ordered state and consequently can not exhibit any finite temperature order-disorder transition, has been solved exactly by Vuillermot and Romerio [exact1, exact2] in 1973, using a group theoretical method. We decided to choose this simple model to apply and test the performance of the WL algorithm for simulation of joint density of states so that a comparison can be made with the exact results available.
We have performed WL simulation using one-dimensional random walk, where the visits are confined to the energy space alone, and also a two-dimensional random walk where the two dimensional space spanned is the energy-order parameter space as well as energy-correlation function space. The partition function can directly be computed as a function of temperature from a knowledge of the density of states, , using the relation
where (the Boltzmann constant has been set to unity). In the two dimensional walk one computes where is the order parameter of the system or a correlation function which are defined in a following section. The partition function can be computed from a knowledge of :
The ensemble average of any function of at a temperature is given by,
We have computed from both 1-d and 2-d random walks and have compared the results obtained with those available from the exact results of Ref. [exact1, exact2].
2 The one dimensional Lebwohl-Lasher model and the exact results
The Hamiltonian of the Lebwohl-Lasher model is given by
where is the second order Legendre polynomial and is the angle between the nearest neighbour spins and . The spins are three dimensional and headless, i.e. the system has the O(3) as well as the local symmetry characteristic of a nematic liquid crystal. A vector order parameter is inadequate for the system and a traceless second rank tensor Q, as defined below, is used to describe the orientational order of the system [lub]. One uses,
where is the -th component of the unit vector , which points along the spin at the site . is the number of particles in the system. In the ordered state is non-zero. In a coordinate system with the Z-axis points along the direction of molecular alignment (director) the matrix is diagonal and for a uniaxial system,
where is the angle between a spin and the director.
MC simulations demonstrate that a three dimensional Lebwohl-Lasher model (=3, =3) exhibits a weakly first order transition, characteristic of a nematic-isotropic transition which is available from the Maier-Saupe model of a nematic in the mean field approximation. On the other hand, for lattice dimensionality =2 and 1, no true long range order is expected since Mermin-Wagner theorem [mermin] predicts a fluctuation destruction of long range order. The =2 LL model has been investigated by a number of authors [kunz, mondal] and the system shows a behaviour qualitatively similar to the two dimensional XY model. A quasi-long range order has been observed in this system and this is believed to be related to the existence of topological defects in the system [dutta].
The one-dimensional Lebwohl-Lasher model has been simulated by [zan] and has also been solved exactly [exact1, exact2]. The system is known to be disordered at all finite temperatures and critical behaviour is expected only at =0, which resembles an one-dimensional Ising model or the one dimensional Heisenberg model. The second rank spin-spin correlation function is defined as
where is the angle between two spins, lattice spacings apart. In the thermodynamic limit one would expect both and to vanish whereas in finite systems because of finite size effects both quantities may appear to have small but finite values.
Vuillermot and Romerio [exact1, exact2] presented an exact solution of the planar (=2) and spatial (=3) versions of the Lebwohl-Lasher model in one dimension (=1) for a nematic liquid crystal, without periodic boundary conditions. They also calculated the two-molecule correlation functions and have shown that these models do not exhibit any finite temperature order-disorder phase transition.
The partition function for the -particle system is given by
where =3/2T. is Dawson function [daw].
The dimensionless internal energy , the entropy and the specific heat are given by
The correlation function is given by
3 Computational details
In the model we have investigated, spins can take up any orientation in the three dimensional space and the orientation of each spin is stored in terms of the direction cosines (, , ). The starting configuration has always been chosen as a random one and to generate a new microstate, a randomly selected spin is chosen and each direction cosine of it is updated as (for i=1,2,3) where is a parameter to be chosen according to some criterion and is a random number between -1 to +1. To preserve the unit magnitude of the spins, (, , ) is always normalized.
The energy of the system in the LL model is a continuous variable and in one dimension (1) it can have any value between to /2. To have a discretization scheme for the implementation of the WL algorithm and for an one dimensional random walk in the energy space, we have chosen an energy range from (- to 0) and divided this energy range into bins each having a width .
We use where, is the number of micro-states corresponding to the -th bin for which the mid-point has the value . Initially we set all (i=1,M) to zero and the logarithm of the modification factor is taken as 1. Whenever a new microstate is generated by rotating a spin, the new system-energy and hence, the macrostate j is determined. Whether the move is accepted or not is decided according to the WL prescription [wl] for the probability
If the state j is accepted, we make and +1, where is the histogram count. Otherwise we make and + 1. This procedure is repeated for 10 MC sweeps (where one MC sweeps consists of attempted moves) and the flatness of the histogram is checked and the cycle is repeated till 90 flatness in the histogram is reached. This completes one iteration, following which we reduce the logarithm of the modification factor , reset the histogram, and the whole procedure is repeated. For each lattice size we have continued with the iterations till gets reduced to 10.
We have also calculated the quantity (Eq. (7)) which gives us the magnitude of the order parameter obtained from the largest eigenvalue of the ordering matrix defined in Eq. (5) and a two-dimensional random walk was performed in the (-) space for this purpose. This is necessary if one intends to determine quantities other than those like free energy, entropy, specific heat etc. which are directly related to energy, and is particularly useful, if one needs to compute, for instance, the variation of the order parameter in presence of an external field. It is also possible to calculate the order parameter from an one-dimensional walk in energy space alone, as has been demonstrated in ref. [jaya] by Jayasri et. al. for a liquid crystalline system (Here one uses the procedure of the so called histogram ‘unweighting’ and ‘reweighting’). But for a more accurate calculation of the order parameter (or the correlation function), it is perhaps a good idea, to generate a two dimensional walk in the (-) space (or in the (-) space) and check its flatness. For a given value of the energy of the system, the order parameter, has a distribution over a certain range of values. The whole range of is 0 to 1 and in order to perform the two dimensional random walk in the energy-order parameter (-) space we divide the two dimensional space into bins. We represent by the bin-width of the bins involving the parameter other than energy in the two-dimensional walk. Each microstate will now correspond to a macrostate labelled by the indices and and the acceptance probability given by Eq. (14) is now modified to
along with an appropriate modification of the procedure described after Eq. (14) for the two-dimensional random walk. Here, for instance, is the density of states for the -th energy and -th order parameter bin.
A two-dimensional random walk in the (-) space is a lot more expensive in terms of the CPU time than an one-dimensional walk in the energy space alone. The problem is particularly severe in a system with continuous energy and becomes worse as the lattice size increases. However, this ensures a much more uniform sampling of the order parameter bins that correspond to a particular energy bin and this improves the overall statistics of the work. It may be pointed out that it is impossible to arrive at a flat histogram in the (-) space if one attempts to visit the entire energy and order parameter ranges accessible to the system. For an one-dimensional walk one normally faces a problem in that, it takes a relatively long time to visit the lowest energy levels and this increases with the increase in system size. For a two-dimensional walk the possibility of uniformly visiting the entire rectangular (-) space is unphysical and one must have a prior knowledge of the range of -bins which are likely to be visited while the system energy has a given value . Our method of simulating the two dimensional random walk has resemblance to the work of Troster and Dellago [dellago], who have applied the WL algorithm to evaluate multidimensional integrals of sharply peaked functions. Our modified approach is elaborated in the following paragraph.
We have first mapped the (-) space which costed us 35 10 sweeps (to be called the pre-production run). The idea is to determine the minimum () and maximum () values of the -bins which are visited while the system energy is for i=1,M. We observe that there are always some -bins within the range (,) for each , where either no sampling or very low sampling takes place during the pre-production run. We therefore checked the histograms of the ) bins in the mapped region of the two-dimensional space and those which attain a 90 flatness during the pre-production run are marked with ‘1’ while other bins are marked ‘0’. This may be clarified as follows. We calculate the average histogram value for those bins which have been visited at least once, thus discarding the bins which are not visited at all. The flatness test (which needs each of the visited bins to have a histogram count at least equal to 90 of the average histogram) is then applied only to those bins and these are labelled with ‘1’. In the ‘production run’ part of the rest of the simulation we check the flatness of only those bins which were marked ‘1’ ignoring what is happening to the others. There is however, always a possibility, since the ‘production run’ generates many more microstates than that in the ‘pre-production run’, that larger areas in the (-) space would get included in the initial ‘visit-map’ or those bins, once marked ‘0’, would subsequently qualify for the label ‘1’. But it is impossible to improve upon the accuracy of the work indefinitely and we decided to stick to the map we obtained during a reasonable amount of the ‘pre-production run’, ignoring what is happening to the discarded bins.
In addition to the two-dimensional random walk in the (-) space we have also performed a number of other two-dimensional random walks. These involve the (-) space where is the correlation function defined in Eq. (8). We have done these only for the =160 lattice, for ranging from 2 to 40 and the ensemble averages of were evaluated for different temperatures using Eq. (3).