The Harmonic Superposition Method

 The total energy density of states associated with a single minimum on the PES is in the harmonic approximation[161]

$$\Omega(E)={(E-E^0)^{\kappa-1}\over \Gamma(\kappa)\prod_{j=1}^\kappa h\nu_j}\theta(E-E^0),$$ (3.4)

where E⁰ is the potential energy of the minimum, $\theta$ is the Heaviside step function, and $\kappa$,the number of vibrational degrees of freedom, is 3N-6. (As we will be comparing our results to those for stationary, non-rotating clusters, we do not include translational or rotational contributions to the density of states.) To calculate the density of states for the whole system, all the minima on the PES need to be considered. We make a superposition approximation and sum the density of states over all the minima low enough in energy to contribute. This approximation is equivalent to assuming that the phase space hyperellipsoids associated with each minimum do not overlap. This gives  

$$\Omega(E)=\sum_{E^0_s<E}{n_s^*(E-E^0_s)^{\kappa-1} \over \Gamma(\kappa)\prod^{\kappa}_{j=1} h\nu^s_j},$$ (3.5)

where the sum is over all the geometrically distinct minima on the surface; n^*_s, the number of permutational isomers of minimum s, is given by n_s^*=2N!/h_s where h_s is the order of the point group of s.

The difficulty with equation 3.5 is that for all but the very smallest clusters this sum involves an impractically large number of minima. Hoare and McInnes[64], and more recently Tsai and Jordan[162] have enumerated lower bounds to the number of geometric isomers for LJ clusters from 6 to 13 atoms. This number rises exponentially with N. Extrapolating this trend gives for LJ₅₅ an estimate of 10²¹ geometric isomers. In such a case, as it is not possible to obtain a complete set of minima, a representative sample is needed. A large set of minima can be obtained by systematic quenching from a high energy MD trajectory. However, this gives a greater proportion of the low energy minima than of the high energy minima. Consequently, if the sample is used in equation 3.5 it is likely to underestimate the density of states due to the high energy minima, and so be inaccurate at high energies.

A method is needed which corrects for the incomplete nature of the sample of minima. This correction can be achieved by weighting the density of states for each known minimum by g_s, the number of minima of energy E⁰_s for which the minimum s is representative. Hence,  

$$\Omega(E)=\sum_{E^0_s<E}{g_sn_s^*(E-E^0_s)^{\kappa-1} \over \Gamma(\kappa)\prod^{\kappa}_{j=1} h\nu^s_j},$$ (3.6)

where the sum is now over a representative sample of minima. The effect of g_s can be incorporated by using the quench statistics. If the system is ergodic and the MD run is performed at constant energy, the number of quenches to a minimum, $\gamma$, is assumed to be proportional to the density of states of the set of g_s minima, i.e. $\gamma(E')_s\propto g_s \Omega(E')_s$. Hence,

  $$\begin{eqnarray} \Omega(E)&\propto &\sum_{E^0_s<E} \gamma \left(E' \right)_s {\O... ...eft(E' \right)_s \left({E-E^0_s \over E'-E^0_s}\right)^{\kappa-1},\end{eqnarray}$$

(3.7)


(3.8)

where E' is the energy of the MD run. This reweighting technique is analogous to histogram Monte Carlo, but instead of determining the configurational density of states from the canonical potential energy distribution, g, effectively a density of minima, is found from the microcanonical probability distribution for quenching to a minimum.

If all the low energy minima are known, the n^* formula (equation 3.5) will be accurate at low energies. Therefore, the proportionality constant in the above equation can be found by matching it to the low energy form of the n^* formula. For LJ₁₃ and LJ₅₅, the term due to the icosahedron is dominant at low energies, and other terms in the sum can be neglected. Comparing the first terms of equations 3.5 and 3.8 gives for the proportionality constant, c  

$$c={n_0^*(E'-E^0_0)^{\kappa-1} \over \gamma \left(E' \right)_0 \Gamma(\kappa)\prod^{\kappa}_{j=1} h\nu^0_j},$$ (3.9)

where E⁰₀ is the energy of the global minimum.

  

Figure 3.5: Microcanonical caloric curves for LJ₅₅. The solid lines were calculated from the $\gamma$ formula, the dashed lines were calculated from the n^* formula, and the dashed line with diamonds shows the simulation results. The sample of minima used in the calculations is marked on the graph.

A critical test of these formulae for $\Omega(E)$ is the predicted microcanonical caloric curve, which for LJ₅₅ has a Van der Waals loop[135]. Using the thermodynamic definition of the microcanonical temperature, $T_\mu$,  

$${1\over kT_\mu}=\left({\partial\ln\Omega\over\partial E}\rig... ...\over\Omega}\left({\partial\Omega\over\partial E}\right)_{N,V},$$ (3.10)

an expression for $T_\mu$ can be derived[153]. For LJ₅₅ we have two samples of minima produced by systematic quenching[153], details of which are given in Table 3.4. Sample A is from a MD run at an energy in the upper end of the coexistence region, and B at an energy just into the liquid-like region. The results for samples A and B using the n^* and $\gamma$ formulations are given in Figure 3.5. It can be seen the n^* formula fails badly, predicting only an inflection in the caloric curve which is too high in energy. This failure is because the contribution to $\Omega(E)$ from the higher energy minima is underestimated. The $\gamma$ formula is much more successful, reproducing the Van der Waals loop at the observed energy. That the harmonic superposition method produces a caloric curve with the correct features shows, as Bixon and Jortner suggested[152], that the distribution of minima is crucial in determining the basic form of the caloric curve. However, the Van der Waals loop is too shallow and lies at too high a temperature. The discrepancy is due to the harmonic approximation. The temperature rises linearly with energy for a single harmonic well. The actual wells of the cluster, however, are flatter than assumed in the harmonic approximation, especially around the transition state regions. Consequently, the cluster spends more time in these high potential energy, low temperature regions of the PES, and so the true temperature is lower than that given by the harmonic approximation.

 

$$E'/\epsilon$$ Number of minima

A 64.7485 989

B 70.2485 1153

For smaller clusters it was found that the performance of the n^* formula improved[155]. This improvement occurs because there are fewer minima on the PES of a smaller cluster and so the set of minima obtained from quenching is more complete.

The harmonic superposition method has three main possible sources of error. The first is associated with the statistical accuracy of the quench frequencies. These errors can be mostly eliminated by having a long enough quench run to ensure ergodicity, and by choosing an appropriate energy for the MD run so that the relevant regions of phase space are all significantly sampled. When studying the thermodynamics of melting it is most appropriate to choose E' to lie in the coexistence region, as in the case of sample A, so that quenches to solid-like and liquid-like states are frequent. The second possible source of error is the assumption that the phase space volumes for each minima can be summed independently, i.e. the hyperellipsoids in phase space do not overlap. If overlap occurred $\Omega(E)$ would be overestimated. Of course, above an energy threshold the true phase space volumes of each minimum are interconnected, but this interconnection is normally due to the anharmonic extension of the phase volumes to form necks along the transition state valleys. The third possible source of error is the harmonic approximation. Near the bottom of the well this is a reasonable assumption, but as the energy is increased some parts of the well become increasingly flat. Consequently the harmonic approximation causes $\Omega(E)$ to be underestimated. Comparison of the caloric curves from simulation and the harmonic superposition method (Figure 3.5) shows that the harmonic superposition method does indeed underestimate $\Omega(E)$ and so the harmonic approximation is likely to be the main source of error.