The effect of well anharmonicity

Here we follow the method of Haarhoff[164,166] to calculate an anharmonic correction term to the density of states. The energy levels for a Morse oscillator are given by

$$E=(n+\half)h\nu-(n+\half)^2{\left(h\nu\right)^2\over 4D_e}.$$ (3.17)

This quadratic can be solved for n. The root corresponding to a bound state is

$$n+\half={2D_e\over h\nu}\left(1-\sqrt{1-{E\over D_e}}\right).$$ (3.18)

Assuming n is continuous and differentiating with respect to E gives a classical density of states,

 $$\begin{eqnarray} \Omega(E)&=&{dn\over dE}={1\over h\nu\sqrt{1-\displaystyle{E\ov... ...E\over 2D_e} +{3\over 8}\left({E\over D_e}\right)^2+\cdots\right),\end{eqnarray}$$

(3.19)

where the square root has been expanded binomially. The expansion is valid if E/D_e is small. The first term in the series is the harmonic term. The full series will diverge as $E\rightarrow D_e$. This divergence is the correct behaviour for the Morse oscillator, but for the case of a cluster isomerization $\Omega(E)$ remains finite as $E\rightarrow \Delta$ from below, where $\Delta$ is the barrier height. As we are seeking a correction term for well anharmonicity, it seems reasonable to truncate equation 3.19 and examine the effect of the first term in this series. Laplace transformation then gives for the partition function

$$z={1\over h\nu}\left({1\over\beta}+{1\over 2\Delta\beta^2} \right).$$ (3.20)

The multi-dimensional partition function is then

$$Z=\prod_{j=1}^\kappa{1\over\beta h\nu_j}\left(1+{1\over 2\Delta_j\beta}\right).$$ (3.21)

Making the approximation that an average value of $1/\Delta$ can be used then gives

$$Z={1\over\prod_{j=1}^\kappa h\nu_j}\sum_{l=0}^\kappa {C^\kappa_l a^l\over\beta^{\kappa+l}},$$ (3.22)

where $a=\langle 1/2\Delta\rangle$ is an anharmonicity parameter and $C^\kappa_l$is a binomial coefficient. Inverse Laplace transforming gives

$$\Omega(E)={1\over\prod_{j=1}^\kappa h\nu_j}\sum_{l=0}^\kappa {C^\kappa_l a^l E^{\kappa+l-1}\over\Gamma(\kappa+l)}.$$ (3.23)

The total density of states is found by summing over all the minima, giving

$$\Omega(E)=\sum_{E^0_s<E}{n_s^*\over \prod^{\kappa}_{j=1} h\n... ...{C^\kappa_l a_s^l (E-E^0_s)^{\kappa+l-1}\over\Gamma(\kappa+l)}.$$ (3.24)

Converting the above equation to the $\gamma$ formulation using equation 3.7 gives

$$\Omega(E)\propto\sum_{E^0_s<E}\gamma(E')_s\left[\sum_{l=0}^\... ...a_l a_s^l (E'-E^0_s)^{\kappa+l-1}\over\Gamma(\kappa+l)}\right].$$ (3.25)

The temperature follows from equation 3.10  

$$T_\mu={\displaystyle{\sum_{E^0_s<E}\gamma(E')_s\left[\sum_{l... ..._l a_s^l(E'-E^0_s)^{\kappa+l-1}\over\Gamma(\kappa+l)}\right]}}.$$ (3.26)

This anharmonic correction has a similar form to that used by Chekmarev and Urmirzakov[165] but has been derived in a different way.

We consider two possible ways of calculating a_s from the potential energy surface of LJ₅₅. Firstly, we consider using the third derivatives of the potential. The dissociation energy of a one-dimensional Morse oscillator can be found from the second and third derivatives of the potential by

$$D_e={(V'')^3\over (V''')^2}.$$ (3.27)

If off-diagonal elements in the multidimensional case are ignored, the above equation can be used to calculate the barrier heights associated with each normal mode. Using analytical third derivatives of the LJ potential, the Cartesian third derivatives for the LJ₅₅ icosahedron were calculated and transformed to obtain the diagonal elements in the Hessian eigenvector basis. This scheme underestimates the anharmonicity and gives values for a_s which when substituted into equation 3.26 have an insignificant effect on $\Omega(E)$, because the off-diagonal elements, which far outnumber the diagonal elements, should not be neglected. However, there is no obvious way to calculate the effect of the off-diagonal elements and for a system such as LJ₅₅ the transformation of the third derivatives into the Hessian eigenvector basis is too computationally expensive for the off-diagonal elements.

The second method considered was to use the barrier heights of our transition state sample. An average barrier height was calculated for each minimum, but for some minima it led to a gross overestimation of the anharmonicity and unphysical caloric curves. The barrier heights are not a good quantitative measure of the anharmonicity of an individual minimum. This method was therefore not considered further.

Therefore, instead of trying to obtain a_s from the PES, we considered how it could be obtained by comparison with the simulation results. If a_s is assumed to be the same for all minima, a value of a can be found which leads to the reproduction of the Van der Waals loop in the caloric curve at the observed temperature and energy, but the temperature difference between the two turning points is still too small. However, one would expect the anharmonicity to depend upon the energy of the minimum: the higher energy minima are likely to have a greater anharmonicity than the solid-like minima.

  

Table 3.5: Values of the anharmonicity parameter, a_i, for LJ₅₅ in the six energy ranges given in Table 3.1. Values for I, II, III and VI were found by fitting the temperatures for these regions to the simulation results.

In §3.3, we saw an example of how the properties of different regions of configuration space could be obtained using suitable order parameters. Indeed, for LJ₅₅, the caloric curves associated with minima in different energy ranges were found. Using the superposition method, it is easy to calculate the properties of a particular region of configuration space simply by restricting the summation to those minima with suitable values of an order parameter, i.e. 

$$\Omega_i=\sum_{E^0_s<E, s\in i} \Omega_s,$$ (3.28)

where $\Omega_i$ is the densities of states associated with region i. From the definition, $1/kT_i=(\partial \ln \Omega_i/\partial E)_{N,V}$, where T_i is the temperature of region i, it follows that  

$$T_i={\displaystyle{\sum_{E^0_s<E, s\in i}\gamma(E')_s\left[\... ..._l a_s^l(E'-E^0_s)^{\kappa+l-1}\over\Gamma(\kappa+l)}\right]}}.$$ (3.29)

We assigned values of the anharmonicity parameter, a_i, to minima in the six different energy ranges (Table 3.1). Values were chosen for regions I, II, III and VI which reproduced the simulation results of Figure 3.4. The caloric curves for each region accurately fitted the simulation curves showing that our anharmonic correction term has an appropriate form. Values of a_i for regions IV and V were chosen that were intermediate between the values for regions III and VI, but as these regions only make a small contribution to $\Omega(E)$ the exact value chosen will only have a very small effect on the overall caloric curve. The values of a_i are given in Table 3.5. As would be expected the anharmonicity increases with the potential energy of the minima.

  

Figure 3.7: Microcanonical caloric curves for LJ₅₅ using the $\gamma$ formula including anharmonic terms with minima samples A (solid line) and B (dashed line). The calculated caloric curves are compared to simulation (dashed line with diamonds).

The microcanonical caloric curve was calculated using these values of a_i. Figure 3.7 shows that for sample A the calculated curve is in remarkable agreement with the simulation results. The calculated Van der Waals loop now has the correct depth, because the effect of the greater anharmonicity of the liquid-like state is to increase the difference in temperature between the two branches of the caloric curve. The success of this method can be understood from the equation  

$${1\over T}=\sum_i{p_i\over T_i},$$ (3.30)

where p_i is the probability that the cluster is in region i. The probabilities at E', p_i(E') are fixed by the quench frequencies,  

$$p_i(E')={\sum_{s\in i}\gamma_s(E')\over\sum_s\gamma_s(E')}$$ (3.31)

Furthermore, the p_i calculated from sample A are in good agreement with simulation results over a wide range of energy (see §3.6). The values of a_i have been chosen to reproduce the simulation temperatures, T_i, and so it is unsurprising that the overall temperature, T, is very accurate. Sample B produces worse results because there are fewer quenches to regions 1-3, and so the quench frequencies are subject to larger statistical errors than for sample A.

For LJ₁₃, the STA temperature distribution is bimodal. The high temperature peak corresponds to solid-like clusters associated with the icosahedral global minimum, and the low temperature peak to liquid-like clusters. The residence times in the solid-like and liquid-like states are much shorter than for LJ₅₅ and so the information provided by short-time averaging for LJ₁₃ is less well-resolved. We partitioned the minima distribution into these two regions, and assigned values of a_i to each. The values of a_i were higher than for LJ₅₅ and the curvature of the T_i curves differed from the simulation results. The apparently greater anharmonicity can be explained by the larger number of surface atoms for LJ₁₃. The incorrect curvature suggests that the energy dependence of the density of states is inappropriate and so the 2nd-order Morse correction term to the density of states was included. The resulting T_i curves fitted the simulation results much more accurately. The partition function for a single minimum including this term is

$$\begin{eqnarray} Z&=&\prod_{j=1}^\kappa{1\over\beta h\nu_j}\left(1+{a\over\beta}... ...{\kappa-l} D^\kappa_{l, m} {3^m a^{l+2m}\over\beta^{\kappa+l+2m}},\end{eqnarray}$$

(3.32)


(3.33)

where

$$D^\kappa_{l, m}={\kappa!\over l!m!(\kappa-l-m)!}.$$ (3.34)

Inverse Laplace transforming and summing over all minima gives for the total density of states,

$$\Omega(E)=\sum_{E^0_s<E}{n_s^*\over \prod^{\kappa}_{j=1} h\n... ...a_s^{l+2m} (E-E^0_s)^{\kappa+l+2m-1} \over\Gamma(\kappa+l+2m)}.$$ (3.35)

In an MD simulation one calculates the kinetic temperature, T_K. For LJ₅₅ the difference between the thermodynamic and kinetic temperatures is negligible. For LJ₁₃ it is still small but not insignificant, and so T_K rather than $T_\mu$ is fitted to the simulation results. Bixon and Jortner found both temperatures produced very similar results in their model calculations of the microcanonical caloric curve[152]. T_K has been obtained for the harmonic superposition method by Franke using the equipartition theorem[154], however this approach cannot be extended to include anharmonicity. We obtain T_K by a different method. Firstly we note that

 $$\begin{eqnarray} \langle E_K\rangle &=&\int_0^E p(E_K) E_K dE_K \\ &=& {1 \over \Omega(E)} \int_0^E E_K \Omega_c(E-E_K) \Omega _K(E_K) dE_K,\end{eqnarray}$$

(3.36)


(3.37)

where $\Omega_K(E_K)$ is the kinetic density of states, $\Omega_c(E_c)$ is the configurational density of states, and the potential energy E_c is given by E_c=E-E_K. Deconvoluting $\Omega_K(E_K)$ from $\Omega(E)$ gives

$$\Omega_c(E_c)={1\over (2\pi m)^{\kappa/2} \prod^\kappa_{j=1}... ...+2m} (E_c-E^0_s)^{\kappa/2+l+2m-1} \over\Gamma(\kappa/2+l+2m)}.$$ (3.38)

Integrating equation 3.37, summing over all minima and substituting into equation 3.3 gives

$$T_K={\displaystyle{\sum_{E^0_s<E} {n_s^*\over\prod^{\kappa}_... ...s^{l+2m} (E-E^0_s)^{\kappa+l+2m-1} \over\Gamma(\kappa+l+2m)}}}.$$ (3.39)

 

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

C 13.7768 95

D 18.3268 295

 

$$a_i/\epsilon^{-1}$$ Region

I 0.470

II 0.625

Sample C (Table 3.6) was used to calculate the caloric curve of LJ₁₃ for the $\gamma$ formulation because it was produced from an MD run in the coexistence region and so should have more accurate quench statistics for the low energy minima, and sample D was used for the n^* formulation because it is from the liquid-like region and so has a larger sample of minima. The values of a_i assigned to solid-like and liquid-like clusters are given in Table 3.7. From Figure 3.8 it can be seen that the caloric curve given by the n^* formulation agrees very well with simulation, because it is possible to obtain a near-complete set of minima for LJ₁₃. (Tsai and Jordan found 1328 minima in a comprehensive search of the LJ₁₃ PES[162]. Many of these, though, are too high in energy to be thermodynamically significant.) The $\gamma$ formulation, however, has too high a transition temperature. This error arises because the probability that the cluster is in the solid-like state derived from quenching is higher than the corresponding probability derived from the STA temperature distributions. In the $\gamma$ formulation the probabilities, p_i, at E' are fixed by the quench statistics (equation 3.31), and so this constraint leads to the higher transition temperature. The difference between the probabilities derived from quenching and the STA temperature distributions may be because there are regions of the PES which are in the basin of attraction of the icosahedron but which have thermodynamic properties similar to a liquid-like well, or because the short-time averaging does not distinguish between the solid-like and liquid-like states with the same resolution as for LJ₅₅. This difference does not stem from the quench method since we obtained similar quench statistics with steepest-descent[170,171], conjugate gradient[87] and eigenvector-following[88] methods.

  

Figure 3.8: Microcanonical caloric curves for LJ₁₃ using the n^* formula (solid line), the $\gamma$ formula (dashed line) and from simulation (dashed line with diamonds). Both the n^* and $\gamma$ formula curves include anharmonic terms.