Computational Study of the Structure and Thermodynamic Properties of Ammonium Chloride Clusters Using a Parallel J-Walking Approach

The thermodynamic and structural properties of (NH$_4$Cl)$_n$ clusters, n=3-10 are studied. Using the method of simulated annealing, the geometries of several isomers for each cluster size are examined. Jump-walking Monte Carlo simulations are then used to compute the constant-volume heat capacity for each cluster size over a wide temperature range. To carry out these simulations a new parallel algorithm is developed using the Parallel Virtual Machine (PVM) software package. Features of the cluster potential energy surfaces, such as energy differences among isomers and rotational barriers of the ammonium ions, are found to play important roles in determining the shape of the heat capacity curves.

thermodynamic properties of the clusters such as C V , the constant-volume heat capacity.
Heat capacities for Lennard-Jones clusters [1,3] have been calculated in this manner, and anomalies in C V such as peaks and shoulders have been found to correspond to the onset of isomerization. These isomerization transitions have been interpreted as signatures of cluster analogs of phase changes for some cluster sizes. [3] In addition, Lopez and Freeman [13] have calculated heat capacities for Ni-Pd clusters, observing a "melting" transition and a low-temperature anomaly corresponding to an order-disorder transition analogous to those seen in bulk bimetallic alloys. [14] Metropolis Monte Carlo simulations [15] in the temperature region corresponding to the onset of isomerization suffer from the inability to sample all available configuration space ergodically. The deficiency of Metropolis sampling, which cannot be effectively overcome by increasing the step size, is its inability to move among the statistically important regions of configuration space that are separated by significant transition state barriers. One way to overcome the problems of Metropolis sampling to use the the jump-walking (or J-walking) move strategy. [1,2] J-walking combines the small step size of Metropolis Monte Carlo with occasional jumps to configurations belonging to a higher temperature ergodic distribution.
The higher temperature distribution contains information about the potential energy surface of the system, and the jumps to configurations belonging to separated but statistically important regions of configuration space are thus attempted with sufficient frequency.
Additional complexities in cluster behavior can occur in molecular aggregates. For example, in the case of ammonium chloride studied here, the effects from the rotations of the ammonium ions significantly contribute to the thermodynamic properties. The rotational potential energy barriers are, in many instances, smaller than barriers between distinct minima, causing the rotational effects to be reflected in the heat capacity at lower temperatures than the isomerization effects.
In the current work we report two principal outcomes. First, we develop a parallel Jwalking algorithm that enables us to carry out a Monte Carlo simulation efficiently in a multi-processor computing environment. Second, the parallel J-walking algorithm is applied to the study of the thermodynamic properties of (NH 4 Cl) n , n = 3-10 clusters. The contents of the remainder of this paper are as follows. In Section II we describe the computational methods, including the model potential and a discussion of the simulated annealing methods used to locate the ammonium chloride cluster isomers. This is followed by a discussion of the J-walking Monte Carlo simulation technique used in computing the constant-volume heat capacity, and the implementation of Parallel Virtual Machine (PVM) software [16] into our Monte Carlo algorithm. Section II concludes with a discussion of the methods used in locating transition states. In Section III we present our results. Isomers of the ammonium chloride clusters found by simulated annealing are shown, along with a few structures of transition states. Then, a constant-volume heat capacity curve for each cluster size is presented, along with a discussion of the features of each curve. Finally, Section IV contains our concluding remarks and directions for future work.

A. MODEL POTENTIAL
In the Monte Carlo simulations carried out in this work, a pair potential dominated by Coulombic interactions is used. The form of the potential is where r represents the entire set of coordinates for the system, r ij is the distance between particles i and j, q Cl = −1.0, q N = −0.4, q H = 0.35, and the remaining parameters are given in Table 1. The potential is derived from potentials of Klein et al. [17] and Pettitt and Rossky [18,19]. The parameters not available from other sources are obtained with the standard combination rules. Because the present simulations are classical and the internal vibrations of the ammonium ions are expected to be quantum mechanical with high frequencies, the ammonium ions in the potential of Eq. (1) are assumed to be rigid tetrahedra. With the potential in Eq. (1), in the CsCl phase at 0K the lattice constant for bulk ammonium chloride is found to be 3.79Å with a cohesive energy of -720 kJ mol −1 . These numbers can be compared with the experimental lattice constant (3.868Å) [20] and the experimental cohesive energy at 298K (-697 kJ mol −1 ). [21] Because clusters have finite vapor pressures in constant temperature simulations, an external constraining potential [22] has been included about the center of mass of each cluster. For this constraining potential we have chosen the same form used elsewhere [23] where r i is the coordinate of particle i, R cm is the coordinate of the center of mass of the cluster, κ has units of energy and R c is a parameter that defines the radius of the constraining potential. In the current calculations R c is taken to be 15 Bohr for the trimer, 20 Bohr for the tetramer, 25 Bohr for n = 5 − 7, and 30 Bohr for n = 8 − 10. In this work κ is taken to be unity.
In the simulated annealing calculations that are used in finding structures of the isomers, it is computationally more convenient to include all degrees of freedom of the system. The internal vibrations of the ammonium ions are described with a harmonic force-field potential of the form where atom 1 is the nitrogen, and atoms 2-5 are the hydrogens, θ ij is the angle formed by hydrogen atoms i and j and the nitrogen atom, θ 0 is the H−N−H angle in a perfect The search for the isomers of the NH 4 Cl clusters has been carried out using the method of simulated annealing. This method has been used previously to locate the isomers of other systems. The basic idea of simulated annealing is to take a system at high temperature and gradually cool the system until the global minimum or a local minimum on the potential energy surface is attained. If the cooling is performed adiabatically, the system will find the global minimum. Using a finite cooling rate allows the system to become trapped in local minima, yielding information about the various isomers.
In this work we use a Brownian dynamics approach. [13,26,27] Starting with chloride and ammonium ions randomly distributed in a fixed volume, the system is propagated in time according to the Langevin equation. As the temperature is decreased, the kinetic energy is drained from the system. When the temperature finally reaches 0 K, the friction term in the Langevin equation drains the remaining kinetic energy from the system, leaving it in some local minimum.
For the smaller clusters the number of isomers is manageable, and the lowest energy isomer for each cluster size can be identified from the Brownian dynamics simulations.
For the larger clusters, however, the number of isomers becomes overwhelming, making it increasingly difficult to identify the lowest energy isomer. We use the J-walking Monte Carlo simulations described below to confirm that the lowest energy isomer is correctly identified.
Because the J-walking Monte Carlo simulations are expected to be fully ergodic, the cluster will be in its lowest energy isomer as the simulation temperature approaches 0 K.
where V 2 and V 2 are the average of the square and the square of the average of the potential energy, respectively, T is the temperature, n is the number of particles, and k B is the Boltzmann constant.
The non-ergodicity of the random walk can be alleviated using a move strategy called J-walking. J-walking combines ordinary Monte Carlo moves with jump attempts to configurations in an ergodic distribution at a higher temperature. We discuss briefly J-walking here and refer the reader to original papers for detailed descriptions. [ In practice, J-walking can be implemented in two ways, each with its own set of drawbacks. The first approach uses tandem walkers, one at a high temperature where Metropolis sampling is ergodic, and one or multiple walkers at lower temperatures. The configurations generated by the high temperature walker are used by the lower temperature walkers for attempting J-walking moves. This tandem approach has been used infrequently because correlations inherent in the Metropolis walks can introduce systematic errors in the J-walking results.
The second approach writes the configurations from the simulation at the J-walking temperature to an external file and accesses configurations randomly from the external file while carrying out a simulation at the lower temperature. Correlation errors are avoided with external configuration files because of two features in this approach. First, as discussed in Reference [3], by writing configurations to an external file infrequently (about once every 40-100 moves), the correlations are partially broken. Second, because correlations still persist between configurations separated by 40-100 moves, it is also necessary to access the external files randomly. The difficulty with this approach is the storage requirements for the external distributions. The large storage requirements have limited the application of the method only to small systems.

D. PARALLEL STRATEGY
As we discussed above, J-walking can be implemented using either tandem walkers or previously generated external distributions. The best features of these two approaches can be combined into a single J-walking algorithm with the use of multiple processors and the Parallel Virtual Machine (PVM) software package. [16] PVM enables processes running on the same or different machines to exchange information while they are executing. We incorporate the PVM subroutines into the Monte Carlo computer code, where these subroutines are used to send and receive configuration geometries and potential energies of the clusters.
Using a multiple-processor computing environment, we can execute a separate Monte Carlo walk on each processor, with PVM enabling these Monte Carlo walks to communicate with each other. Instead of generating external distributions and storing them before the actual simulation, we generate the required distributions "on the fly" and pass them to the lower-temperature walkers. The walker accepts or rejects the configuration it receives based on Eq. (5) of Reference [1].
Our computational scheme consists of two kinds of Monte Carlo processes, the configuration generating processes and the computing processes. The generating processes are designed to be a source of ergodic distributions at particular temperatures. Each set of computing processes carries out Monte Carlo walks at temperatures below that of the generating process, periodically jumping to a configuration provided by the generating process to maintain ergodicity.
A diagram of a model PVM J-walking simulation is shown in Fig. 1. Each box represents a process in a parallel machine executing a Metropolis simulation at a particular temperature. The set of boxes on the left-hand side of Fig. 1 represents the generating processes, executing walks at temperatures T 1 , T 2 and T 3 . T 1 is assumed to be sufficiently high that Metropolis Monte Carlo is ergodic without modification. The set of boxes on the right-hand side of Fig. 1 represent the computing processes. The necessity of having the generating processes running at different generating temperatures is precisely analogous to generating external distributions at these temperatures in a J-walking simulation using the serial algorithm. [1] As in the serial J-walking method, a new distribution is needed when the jump acceptance rate falls below 10-15%. The generating processes at T 2 maintain their ergodicity just like the computing processes by making jump attempts to configurations contained in the T 1 processes. In turn, the T 2 processes serve as sources of ergodic distributions for computing processes executing at temperatures between T 2 and T 3 . Additionally the generating processes at T 3 are made ergodic by periodic jumps to the configurations contained in the T 2 processes. If still lower temperature generating processes are needed, this procedure can be continued. It is important to note that four processes are executed at each generating temperature. The four downward arrows from T 1 to T 2 and from T 2 to T 3 denote that the generating processes at each temperature are independent of each other.
The boxes on the right-hand side of Fig. 1 represent the computing processes. The thick arrows pointing from the generating processes to the computing processes indicates that each of the four generating processes is feeding configurations to each computing process.
The reason for having four generating processes at each temperature instead of just one is the following. The bottleneck of this algorithm is the frequency at which a generating process can pass a configuration to a computing process or to a lower-temperature generating process. This frequency must be kept sufficiently low to allow the Monte Carlo simulation at the generating temperature ample opportunity to alter the configuration of the system.
Previous J-walking studies of Lennard-Jones clusters have shown that a configuration can be saved once every 40-100 Monte Carlo passes, depending on the size of the cluster. [3] It has also been determined that a computing process should attempt jumps to the ergodic distribution as often as once every 10-15 passes. We can maintain the 10-15 pass criterion in the computing process by having each computing process access four independent generating processes. The result is that each generating process is accessed once every 60 moves, but since there are four, the computing process can attempt a jump to one of the four generating processes once every 15 moves. In tests performed on (NH 4 Cl) 9 , storing a configuration every 100 passes instead of every 60 did not affect the calculated heat capacity. The computing time is reduced by a factor of four in the case where the number of processors is not limited.
As in the serial algorithm there are correlations in the Metropolis walks of each generating process that must be broken. To break the correlations, each configuration-generating process keeps an array of configurations in memory. When a configuration is transmitted to a lower-temperature process, it is a configuration randomly chosen from this array, not the current configuration of the walker. The current configuration of the walker then replaces the configuration just passed from the array to another processor. The minimum size for this array has been found to be 2500 configurations from tests on the (NH 4 Cl) 9 cluster.
In all calculations presented below, a 5000-configuration array is used. Compared to array sizes of at least 50000 configurations used in the serial algorithm, the arrays in the parallel method are small and do not inhibit applications of the method to large systems.
An issue that must be addressed is the initiation of the 5000-configuration array at each generating temperature T 1 , T 2 , etc. We have examined two approaches. One approach is to create the arrays during the computation. The startup time required to populate the configuration arrays sequentially can be prohibitively long. For example, the computing processes between T 1 and T 2 and the generating processes at T 2 must wait until the 5000configuration arrays at T 1 are created. The other approach, which we have found to be more useful, is to create the small 5000-configuration distributions using serial J-walking coefficients are then determined, yielding a trajectory with a particular total energy. The details of finding the Fourier coefficients are given elsewhere. [28,29] The trajectory with the lowest total energy is usually found to pass through or near the transition state, identified by having exactly one negative eigenvalue in its Hessian matrix.
The other method used in locating the transition states is the eigenmode-following technique initially proposed by Hildenbrand [30] and developed further by other workers. [31] More recently, Tsai and Jordan [6] and Wales et al. [9] have employed this method to locate transition states of atomic and molecular clusters. Briefly, a search is conducted by starting at a local minimum and following one of the normal modes by maximizing the potential energy along that particular mode while minimizing the potential energy along the remaining normal modes, until a transition state is reached.
For the case of (NH 4 Cl) 4 and (NH 4 Cl) 9 we are concerned with rotational barriers of the ammonium ions. The double-ended trajectory method for finding transition states is ideally suited to these motions because of their simplicity. For the more complex transition states of the (NH 4 Cl) 3 cluster connecting distinct isomers we use the eigenmode-following method.
The double-ended trajectory method is unable to locate the transition states between the isomers of the trimer with a low number of Fourier expansion coefficients. To apply the double-ended trajectory approach to a case with a complex rearrangement, a large number of Fourier coefficients is required to describe the path of the system between the two minima, making the method computationally expensive. Rather than pursuing the transition state search with the double-ended method, we have opted to use the eigenmode-following method in the transition state searches in (NH 4 Cl) 3 .

A. STRUCTURES
The ammonium chloride monomer is known [32] to be a van der Waals complex rather than an ionic species and cannot be treated using the model potential of Section II.A. We do not investigate (NH 4 Cl) 2 because it has only a single isomer and is not expected to exhibit interesting thermodynamic behavior. We assume that all clusters in this work are ionic. The cluster size at which there is a transition from van der Waals to ionic bonding represents an interesting and complex electronic structure problem that is beyond the scope of this work.
We begin our discussion of the ammonium chloride cluster structures with (NH 4 Cl) 3 , for which we have found three isomers with our potential model. These three isomers, located with the simulated annealing procedure described in the previous Section, are shown in Figs.  The relative stability of the even-number clusters versus the odd-number clusters can be compared by examining the gain in binding energy of the lowest energy isomer, ∆V = V n−1 − V n as a function of the cluster size n, where V n is the potential energy of the lowest energy isomer of size, n. The trend observed in Fig. 5 clearly shows that the gain in potential energy is greater when going from an odd to an even cluster, except for (NH 4 Cl) 9 .
The deviation of (NH 4 Cl) 9 from the observed trend can be understood by examining its geometry, shown in Fig. 4(e). The (NH 4 Cl) 9 cluster forms a rock-salt structure that has more in common with the even-number clusters than with the other odd-number clusters.

B. HEAT CAPACITIES
The constant-volume heat capacity curve for (NH 4 Cl) 3 is shown in Fig. 6. The error bars in this and all subsequent heat capacity curves represent two standard deviations. Also, only the potential energy contribution to the heat capacity is shown in this and all subsequent heat capacity curves, i.e. we set There are no additional structural features from the constant kinetic energy contribution in Eq. (4). The curve starts at approximately 10.5, the equipartition value for the trimer, undergoes a slow increase in the 0-80 K range, then rises rapidly to a peak at 140 K, and finally the heat capacity declines gradually over the remaining temperature range. The presence of an early peak in the heat capacity is a direct consequence of the small differences in energy among the isomers.
We have found that during the simulation at 50 K all configurations are in the potential energy wells belonging to isomers (a) and (e) in in most of the cases studied in this paper. As we shall discuss in the context of (NH 4 Cl) 9 , the temperature region where C v is flat is the region where the NH + 4 ions are nearly freely rotating. Using the method of double-ended classical trajectories, the barrier for rotation of a NH + 4 ion about the N-H bond pointing outward from the cluster has been determined to be 4300 K.
The heat capacity curve for (NH 4 Cl) 5 is shown in Fig. 8(a). It has qualitative similarities to the (NH 4 Cl) 4 heat capacity shown in Fig. 7, with a melting peak in the 1000-1100 K region. The melting peak is, however, significantly smaller relative to the T = 0 K heat capacity than in (NH 4 Cl) 4 . This is a trend that we observe for all the odd-number clusters, with the exception of (NH 4 Cl) 9 ; namely the ratio of the maximum heat capacity relative to the T = 0 K value is significantly larger for the even-number clusters than for odd-number clusters.
Unlike the tetramer, (NH 4 Cl) 5 has several low-energy isomers. Isomerization in (NH 4 Cl) 5 is seen at significantly lower temperatures than in the tetramer. The isomerization at the lower temperatures in (NH 4 Cl) 5 occurs among isomers with compact structures, two of which are shown in Fig. 9 (a) and (b). Figures 9(c) and 9(d) show examples of higherenergy open structures of (NH 4 Cl) 5 . The open-type structures are at least 12000 K above the lowest-energy isomer and become accessible in the vicinity of the 1000-1100 K melting peak. In contrast, only the lowest energy isomer of (NH 4 Cl) 4 has a compact structure and the remaining isomers are all high energy open structures (see Fig. 3). Figure 8(b) shows the heat capacity curve for (NH 4 Cl) 6 , and the lowest energy isomer is shown in Fig. 3(b). The heat capacity curve is similar to that of (NH 4 Cl) 4 . The region between 500 K and 800 K is remarkably flat, and the melting peak at 1200 K is sharp compared with (NH 4 Cl) 5 . (NH 4 Cl) 6 has at least one other low-energy isomer (not shown).
These low-energy isomers of (NH 4 Cl) 6 are 2051 K apart and isomerization between them is seen at temperatures as low as 400 K. The remaining isomers found in this work are open-type structures and are in the range of 12000-14000 K above the lowest energy isomer, becoming accessible at temperatures near the melting peak.
The lowest energy isomer for the (NH 4 Cl) 7 cluster is shown in Fig. 3(c). The heat capacity curve, shown in Fig. 8(c), is the only curve in this study that clearly exhibits two distinct peaks. The first peak occurs around 600 K and is in the same region where other cluster sizes (most notably 4,6,8, and 9) show a flat region. In this case, however, the peak at 600 K is a result of isomerization between the global minimum and a local minimum approximately 2000 K above the global minimum. The second peak, centered at about 1200 K, is a melting peak indicating the onset of isomerization to open-type structures.
To investigate further the origin of the 600 K peak in the heat capacity of (NH 4 Cl) 7 , we have performed quench studies of the configurations obtained from the simulations at temperatures of 425, 600, and 800 K. At each temperature, configurations were taken each 750 Monte Carlo passes, and 400 configurations at each temperature were quenched to the nearest local minimum using Brownian dynamics. At 425K, well before the peak in the heat capacity, almost all 400 configurations from an ergodic distribution quench to the lowest energy isomer. The dominance of the lowest energy isomer is reflected in the large peak in the histogram in Fig. 10(a) at 0 K and the small peaks representing higher energy isomers (the energies of the isomers in Fig. 10 are given relative to the lowest energy isomer defined to be 0 K). At 600 K (Fig. 10(b)) slightly more than half of the configurations quench to the lowest energy isomer, with the majority of remaining configurations distributed between two isomers 830 K and 2194 K above the lowest energy isomer. The time-scale for isomerization at 600 K is still slow, corroborated by examining configurations of the (NH 4 Cl) 7 cluster saved during the simulation. The infrequent isomerization events result in large fluctuations in the potential energy, leading to a peak in the heat capacity at 600 K. Finally, at 800 K (see Fig. 10(c)) two isomers are seen to dominate the histogram with roughly equal populations.
Rapid isomerization near 800 K leads to the observed decrease in the heat capacity. The second peak at 1200 K is the melting peak.
To determine the features of (NH 4 Cl) 7 that are unique in this series, we compare the energetics of the clusters. (NH 4 Cl) 4 has only one compact structure, the global minimum (see Fig. 3), and the rest of the isomers are at least 13000 K higher in energy. (NH 4 Cl) 5 has a few low-energy isomers (see Fig. 9), but there is a large gap between the group of low-energy isomers and the group of open-type higher energy isomers. (NH 4 Cl) 6 has just one other low energy isomer in addition to the global minimum, and the energy difference between the two low energy isomers in (NH 4 Cl) 6 is roughly the same as the difference between the two isomers of (NH 4 Cl) 7 that dominate the the histogram in Fig. 10(c). However, in contrast to the case of (NH 4 Cl) 7 , (NH 4 Cl) 6 does not have any other low-energy isomers. The melting peak in (NH 4 Cl) 6 , as in the other even-number clusters, is significantly larger compared to the T = 0 K heat capacity than in the odd-number clusters. Any decline in the heat capacity accompanying the onset of rapid isomerization between the two low-lying isomers of (NH 4 Cl) 6 is masked in the 800 K region by the melting peak. Conversely, the melting peak in the heat capacity of (NH 4 Cl) 7 is small and does not mask the drop in the heat capacity.
The structure (Fig. 3(d)) and the heat capacity of (NH 4 Cl) 8 (Fig. 7(d)) are similar to those of (NH 4 Cl) 6 (see Figs. 3(b) and 7(b), respectively). The slope of the curve in the 500-800 K "shoulder" region is not as flat as in the case of (NH 4 Cl) 6 . Similar to (NH 4 Cl) 6 , there are a few low-energy isomers with compact structures accessible at lower temperatures, and the majority of isomers are higher-energy open-type structures that become accessible in the vicinity of the melting peak.
The (NH 4 Cl) 9 cluster, shown in Fig. 3(e), has more in common with the even-number clusters than it does with the odd-number clusters. The structure of (NH 4 Cl) 9 is slab-like, reminiscent of rock-salt structure of NaCl. The heat capacity curve, shown in Fig. 7 shows a flat region in the 400 -800 K region and a large melting peak at about 1100 K. An interesting feature of the (NH 4 Cl) 9 cluster is the low rotational barrier of the ammonium ion in the center of the cluster. The center ammonium ion has 12 identical potential energy minima as it rotates by 360 degrees, and the barrier to rotation between each pair of minima is only 48 K. Although we can see the center ammonium ion hopping between different minima at temperatures as low as 5 K, the contribution of this motion to the heat capacity is small and masked by the steady rise in the heat capacity from the anharmonic vibrational motions in the rest of the cluster.
We have claimed that the main contributing factor to the appearance of a temperature region where the slope of C v is small is the onset of free rotation of the ammonium ions. As a qualitative test of this claim we carry out two additional Metropolis Monte Carlo simulations in this temperature range for (NH 4 Cl) 9 . In the first simulation, only the rotational Monte Carlo moves are allowed, with the translational moves of the ammonium and chloride ions excluded from the simulation. In the second simulation, the translational moves of the ammonium and chloride ions are included, but the rotational moves of the ammonium ions are not allowed. It must be stressed that this calculation has been done to demonstrate qualitatively the effects of the onset of free rotation on the heat capacity. The two curves are shown in Fig. 11. The upper curve is the heat capacity resulting from the translational moves, and the lower curve includes only rotations. The upper (translational) curve increases through the entire temperature range in Fig. 11. The lower (rotational) curve rises to a maximum at about 600 K and then undergoes a decrease. The increase in the rotational contribution to the heat capacity coincides with the ammonium ions undertaking infrequent hops between the equivalent orientations. Above 600 K, however, the motions of the ammonium ions are more representative of free rotation rather than hops, which is consistent with the decreasing rotational contribution to the heat capacity. We have confirmed this by examining the configurations of the (NH 4 Cl) 9 cluster at several temperatures throughout the flat region of the heat capacity. The rotational barrier of a corner ammonium ion is 2908 K, and the rotational barrier of a side ammonium ion is 1948 K. Both of these values are consistent with the onset of free rotation at 600 K. The sum of the rotational and translational contributions is nearly flat as in the full C v curve shown in Fig. 8(e).
The heat capacity for the remaining cluster in the study, (NH 4 Cl) 10 , is shown in Fig. 8(f) and the lowest-energy isomer is shown in Fig. 3(f). The heat capacity rises steadily, reaches its maximum at 700 K and then slowly decreases with increasing temperature. There is a change in slope in the 1000 -1200 K region, where clusters of size 4-9 exhibit a well-defined peak.

CONCLUDING REMARKS
In this paper we have presented a study of structures and temperature-dependent thermodynamic properties of (NH 4 Cl) n clusters, n=3-10. The structures have been determined using the approach of simulated annealing, and the constant-volume heat capacities have been computed using a newly developed parallel PVM J-walking algorithm. We can compare our results to the structures of (NaCl) n clusters determined by Phillips et al. [7] For (NaCl) 3 , only one isomer is found, compared to three for (NH 4 Cl) 3 . The structure of (NaCl) 3 is a planar ring, similar to the highest energy isomer of (NH 4 Cl) 3 , shown in Fig. 2(e). A close examination of the (NH 4 Cl) 3 isomer shown in Fig. 2(c) together with the absence of such a structure in (NaCl) 3 suggests that the (NH 4 Cl) 3 isomer shown in Fig. 2(c) is stabilized by the attractive hydrogen-chloride interactions. The geometries of the lowest energy isomers for (NaCl) n , n = 4,7,8,10 are similar to the geometries of the lowest energy isomers that we find for the respective (NH 4 Cl) n clusters. For the remaining isomers, n = 5,6, and 9, the geometries of the the second lowest isomers of (NaCl) n are similar to those of the lowest energy isomers of (NH 4 Cl) n . In the work of Diefenbach and Martin [33], lowest energy isomers are identified for for several alkali halide clusters. They find that the lowest energy isomer geometries are dependent on both the constituent ions and the potential model used. It is noteworthy that unlike bulk ammonium chloride crystal, for which the cesium chloride structure is thermodynamically stable, we have found no cesium chloride type structures in the clusters. The cluster size at which the cesium chloride structure begins to appear represents an interesting problem.
A major thrust of this paper is the development of the parallel J-walking algorithm with the use of Parallel Virtual Machine (PVM) package. The ability to maintain the J-walking distributions dynamically at each required temperature makes it possible to perform an ergodic Monte Carlo simulation spanning the entire temperature range of interest. The PVM J-walking algorithm is designed to take full advantage of a multi-processor computing environment. Ergodic simulations spanning the entire temperature range of interest require the same amount of time as a single Metropolis Monte Carlo simulation, provided that a separate processor can be assigned to each process in a PVM algorithm. PVM is designed to be most efficient when the passing of information among the different processes is infrequent.
In that respect, PVM is ideally suited for J-walking simulations because the most frequent jumps are made only once every fifteen passes in the computing processes. The PVM J-walking algorithm does not require the storage of large external configuration files needed in the serial algorithm, and ergodic Monte Carlo simulations can be performed on significantly larger systems than possible previously.
The constant-volume heat capacities for each cluster size have been determined, revealing two important temperature regions. The first region, seen most prominently in the (NH 4 Cl) n , n=4, 6, 8, and 9 clusters, is the flat region in the 500 -800 K vicinity. This flat region in the heat capacity is apparently a consequence of competing contributions from rotational motion of the ammonium ions and the anharmonic vibrational motion in the clusters. The other feature of the heat capacity curves for the majority of the clusters is the prominent maximum in the 1000-1200 K range that we have called a melting peak. This peak coincides with isomerization to open-type configurations. Cluster geometries at temperatures below 1000 K are normally dominated by compact structures belonging to potential energy wells of low lying isomers.
Future work on ammonium chloride clusters will be directed toward a quantum mechanical treatment of this system. J-walking has already been extended to quantum systems by incorporating it into the Fourier path integral Monte Carlo method. [2] Realistic inclusion of all degrees of freedom in this system cannot be accomplished using classical mechanics because of the intrinsic quantum-mechanical nature of the internal vibrations of the ammonium ions. As an example of the potential importance of quantum contributions, it can be anticipated that the quantum-mechanical treatment will have a significant effect on the heat capacity curve of the (NH 4 Cl) 3 cluster because the differences in zero-point energies of the isomers can be expected to rearrange the ordering and change the energy spacings between the isomers. Furthermore, significant isotope effects are known in the bulk system, implying the importance of quantum effects in the clusters. [34] TABLES  [17] c Reference [18] d Reference [19] e Units of energy in Hartree and units of length in Bohr          Carlo walk whose nearest local minimum on the potential energy surface has potential energy V . Figure 11. Upper curve shows the constant-volume heat capacity for (NH 4 Cl) 9 in the absence of rotational moves, and the lower curve is the heat capacity with the translational motion of ammonium and chloride ions excluded. The sum of the two curves changes little explaining the plateau features in Fig. 8.