Case study 2: Optimizations and benchmarking - steventopper/TransRot GitHub Wiki
We are grateful to William Wang (The Lawrenceville School, NJ) for completing and sharing this case study for the TransRot community.
Optimizations of small noble gas clusters
In Part 2, we use the OPLS-AA Lennard-Jones parameters established for noble gas atoms to carry out optimizations of small clusters. Simulations using the Lennard-Jones parameters given in the previous wiki post for the noble gases were conducted in TransRot to determine minimum energy structures of pure noble gas clusters having 2 to 6 atoms. This was done in order to compare the cluster structures and their energies to predictions by Hoare and Pal, which are tabulated in the Cambridge Energy Landscape Database. The energies and point groups of these clusters are given in Table 1.
Table 1. Literature Values for Energies of Minimum-Energy Cluster Structures
| Cluster size | Energy ($\epsilon$) | Point Group |
|---|---|---|
| 2 atoms | -1.000000 | $D_{\infty h}$ |
| 3 atoms | -3.000000 | $D_{3h}$ |
| 4 atoms | -6.000000 | $T_{d}$ |
| 5 atoms | -9.103852 | $D_{5h}$ |
| 6 atoms | -12.712062 | $O_{h}$ |
Runtime configurations and the dbase.txt file are provided below. Four simulations were conducted for each single-atom cluster on a personal computer. The initial annealing temperature was generally set to 150K, based on the melting point of xenon. Only helium clusters had their initial annealing temperature set to 25K, which was necessary due to helium's far weaker intermolecular forces. The initial box size parameter was set to $10\text{Å}$ for all systems; the code automatically increased the box size to $15\text{Å}$ for the larger 6-atom clusters during the course of the run. For comparison, TransRot outputs were converted into units of $\epsilon$ by dividing each final cluster energy by the original $\epsilon$ value in kcal/mole.
Figure 1: Annealed $\mathrm{Xe}_6$ cluster structure visualized in Avogadro
The results of the study are shown below. Each annealing calculation took about 3 seconds or less on a single core. We see that the correct lowest-energy basin was located for all of the systems considered here, with percent deviations of 0.045% or less in all cases. Note that since TransRot does not currently implement a final gradient-based optimization, we do not expect to resolve the final structure with 0.0% deviation.
Table 2. Comparison of TransRot Optimized Energies to Literature Values
| Element | Cluster size | Energy (kcal/mol) | Energy ($\epsilon$) | % Deviation |
|---|---|---|---|---|
| He | 2 | -0.020000 | -1.00000 | 0.0000 |
| 3 | -0.060000 | -2.99996 | 0.0001 | |
| 4 | -0.119991 | -5.999562 | 0.0073 | |
| 5 | -0.182056 | -9.102807 | 0.0115 | |
| 6 | -0.254128 | -12.706402 | 0.0445 | |
| Ne | 2 | -0.069000 | -1.000000 | 0.0000 |
| 3 | -0.207000 | -2.999996 | 0.0001 | |
| 4 | -0.413976 | -5.999650 | 0.0058 | |
| 5 | -0.628053 | -9.102218 | 0.0180 | |
| 6 | -0.876846 | -12.707920 | 0.0326 | |
| Ar | 2 | -0.233900 | -1.000000 | 0.0000 |
| 3 | -0.701700 | -2.999998 | 0.0001 | |
| 4 | -1.403319 | -5.999655 | 0.0057 | |
| 5 | -2.129259 | -9.103289 | 0.0062 | |
| 6 | -2.972871 | -12.710008 | 0.0162 | |
| Kr | 2 | -0.317000 | -1.000000 | 0.0000 |
| 3 | -0.951000 | -2.999999 | 0.0001 | |
| 4 | -1.901978 | -5.999932 | 0.0011 | |
| 5 | -2.885681 | -9.103094 | 0.0083 | |
| 6 | -4.028944 | -12.709601 | 0.0194 | |
| Xe | 2 | -0.433000 | -1.000000 | 0.0000 |
| 3 | -1.299000 | -2.999999 | 0.0000 | |
| 4 | -2.597944 | -5.999870 | 0.0022 | |
| 5 | -3.941667 | -9.103158 | 0.0076 | |
| 6 | -5.503664 | -12.710541 | 0.0120 |
Figure 2. Configuration information used in this case study (config.txt)
/Sawtooth Annealing Configuration Options
//Use 2 or more spaces between label and value
Max Temperature (K): //150.0K (all) 25.0K (He)
Moves per Point: 20000
Points per Tooth: 10
Points Added per Tooth: 0
Number of Teeth: 4
Temperature Decrease per Tooth: 0.8
Max Translation Distance (Angstroms): 0.25
Magwalk Translation Multiplication Factor: 10.00
Magwalk Translation Probability: 0.1
Max Rotation (Radians): 0.523599
Magwalk Rotation Probability: 0.1
Length of Cubic Space: 10.00
Max Failures During Propagation: 20
Use Input.xyz (true/false): false
0K Finale (true/false): false
Static Temperature (true/false): false
Write Energies When Static Temperature (true/false): true
Override Interaction Parameters (true/false): false
Write Acceptance Ratios (true/false): true
//Below are specifications for desired cluster size and noble gas i.e., He6
Possible discovery of local minima during annealing
In certain trials, energy values at the end of the last tooth occasionally oscillated between the minimum energy value and values clearly above the minimum energy values. For example, in one simulation of a 6-atom neon cluster, the following behavior was observed:
...
Initialization and propagation done in 0.027s
Cluster movement constrained within cube with side length 15.0
Starting Energy: -0.012837376980522992
Energy at end of tooth 1: -0.8451735560158807
Energy at end of tooth 2: -0.7956055590276574
Energy at end of tooth 3: -0.848612785295177
Energy at end of tooth 4: -0.876717108162562
Annealing done in 2.798s.
While the energies at the end of tooth 1, tooth 3, and tooth 4 are relatively close to the energy of the global minimum-energy configuration, the energy at the end of tooth 2 is not. The simulation quenches to an anomalously high-energy structure at the end of a tooth when a temperature drop occurs too quickly, leaving the cluster trapped on an energy shoulder as opposed to falling into an energy well.
One might guess that the energy at the end of tooth 2 is a local minimum. However, the 6-atom Lennard-Jones cluster is known (see Daniel B. Faken's doctoral thesis) to have exactly two higher-energy isomers, both with energies only about 0.409 $\epsilon$, or 0.0282 kcal/mol, above the global minimum. The energy of the structure at the end of tooth 2 is about 0.081 kcal/mol above the final value from tooth 4. This structure might represent an unresolved local minimum, but without a gradient-based optimization of the structure it is difficult to tell for certain.
This behavior is actually desirable if the goal is to locate not only the global minimum structure, but also local minima / isomers / conformers (structures which are higher in energy than the global minimum-energy structure). As long as the system is not too large, TransRot can be used for this purpose as well by carrying out multiple quenches ("teeth") with slightly faster annealing rates, with magwalking turned off and the uppermost temperature scaling factor ($\alpha$, or "Temperature decrease per tooth" in the config.txt file) set equal to 1.
End of Part 2
-WW and RQT