Investigation of Melting Dynamics of Hafnium Clusters - American

Subscriber access provided by University of Colorado Boulder

Article

Investigation of melting dynamics of hafnium clusters Wei Chun Ng, Thong Leng Lim, and Tiem Leong Yoon J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.6b00553 • Publication Date (Web): 08 Feb 2017 Downloaded from http://pubs.acs.org on February 12, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Investigation of Melting Dynamics of Hafnium Clusters Wei Chun Ng,a*‡ Thong Leng Lim†‡ and Tiem Leong Yoona‡ a

School of Physics, Universiti Sains Malaysia, 11800 USM Pulau Pinang, Malaysia.

AUTHOR INFORMATION Corresponding Author: Wei Chun Ng *

Email: [email protected]

Present Addresses: †

Faculty of Engineering and Technology, Multimedia University, Jln. Ayer Keroh Lama, 75450

Melaka, Malaysia. Author Contributions: The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. ‡These authors contributed equally, Wei Chun Ng, Thong Leng Lim and Tiem Leong Yoon respectively.

Keywords: Cluster melting, hafnium clusters, ground state structures, COMB potential, molecular dynamics (MD) simulation, density functional theory (DFT), low lying structures (LLS), chemical similarity, similarity index.

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 49

Abstract: Melting dynamics of Hafnium cluster is investigated using a novel approach based on the idea of chemical similarity index. Ground state configurations of small Hafnium clusters are first derived using Basin-Hopping and Genetic Algorithm in parallel tempering mode, employing COMB potential in energy calculator. These assumed ground state structure are verified by using Low Lying Structures (LLS) method. Melting process is carried out either by using direct heating method or prolonged simulated annealing. Melting point is identified by caloric curve. However, it is found that the global similarity index is much more superior in locating pre-melting and total melting point of Hafnium clusters.

1. Introduction Recent progress in nanotechnology has caused a surge in the interest of searching for a new generation of nanomaterials with exotic or desirable functionalities. Some of these newly established materials are nanoclusters and nanoalloys. Nanoclusters find their applications in crystal growth, adsorbents, thin film, catalysis, magnets design, medical use1 and other areas. The catalytic effect of nanoclusters is strongly related to their geometry, such as the core-shell structures which are commonly found in bimetallic nanoalloy. Son et al. 2 demonstrated the functionality of Ni/Pd core-shell nanoparticles in catalyzing the Sonogashira coupling reactions in a more economical way. Certain non-magnetic material displays a useful nature of magnetic behavior when they are in the form of a cluster. For example, Park and Cheon3 managed to synthesize a magnetic nanoalloy of cobalt-platinum via experiments, which is believed to be used in magnetic device applications. The nanoclusters, particularly gold metal clusters are of interest in the medical field due to their enhanced optical properties and inert nature of chemical reactions.4 The prediction of ground state structure of a cluster is a non-trivial task. Currently there are some widely used global optimizer in locating the ground state structure of atomic cluster, such as


2

Page 3 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulated annealing, genetic algorithm, basin hopping, PSO and evolutionary algorithm, to name a few. By coupling a global optimizer and interatomic potential evaluator, the global minimum configuration or ground state structure can be located before other properties are analysed. Surface melting is of paramount importance in the investigation of clusters. For atomic cluster cladded with three-dimensional structure, the surface atoms have fewer nearest neighbors and they are thus weakly bounded, and less constrained in thermal motion.5,8 Surface energy increases as the atomic cluster size decreases, thus bringing down the melting point of the cluster. According to Hanszen, melting is the temperature of equilibrium between a thin liquid layer and the solid sphere core. Experimental studies confirm that surface atoms bring down the melting temperature. Aguado9 argues that in nanometer region, when the size of the atomic cluster gets smaller, the surface to volume ratio increases and surface energy starts to show its contribution. Cluster melting is often studied through simulation. Aguado argued that a reliable molecular dynamics simulation has to be at least several nanoseconds long9 such that quantum chemistry method is rarely employed or only applied to small cluster due to high computational demand. By assuming the cluster to be studied can be modelled by Lennard-Jones potential10, Na cluster is studied through molecular dynamics by computing the short-time-average value for each calculated temperature9. Single atomic species of alkali metal atomic cluster displaying core-shell structure have also been reported in literature36,37. Size dependent melting is observed and further supported by experiments6. Bottani16 used Car-Parrinelo ab initio molecular dynamics simulation to study Si, Ge and Sn clusters up to 13 atoms. Chuang17 on the other hand, in his ab initio Langevin molecular dynamics simulation, revealed that for Sn6, Sn7, Sn10 and Sn13, all the melting temperatures are well above the bulk melting value. Steenbergen et. al.18 carried out ab initio


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 49

molecular dynamics on gallium cluster, with parallel tempering method for sizes between 7-12 atoms. Larger cluster behaves like bulk with shaper melting points when it melts. The melting point of 2.5 nm cluster size of gold is 930K12 while the bulk melting value is 1336K. The average latent heat of melting shows a size-scaling behavior similar to that of the melting point9 and a significant size dependence of the interfacial tension13,14 also being detected in experiments. 3

When the cluster size is relatively small, melting point, 𝑇𝑚 increases with √𝑛, and could show oscillation. The peaks in the oscillation can be explained by the existence of magic number in finite sized clusters.6 Duan et al.11 in their molecular simulation of pure Fe cluster showed that for large clusters such as Fe300, the surface can exist as molten phase while the core maintains the solid phase at the same time. The temperature range of melting is, however, narrow and precise. This coexistence of different phases on the surface and in the core proves that the melting of different shells at a different temperature is probable for pure element clusters. Logically, one would expect that only a core-shell clusters of different elements to undergo this kind of melting pattern. Melting of clusters is greatly affected by the nature of the PES that causes a large behavioral change. Duan et al.11 also concluded that, for small clusters, the coexistence is over time; whereas, for big clusters, the coexistence is over space. This paper is an effort to investigate the thermal property of hafnium in the form of nanoclusters. Further miniaturization of electronic components poses various challenges to silicon as the basic material and necessitates the search of next generation material in replacement of silicon. Among others, hafnium, in the form of various allotropes, might be an alternative material to silicon in overcoming the die shrinkage limit of the silicon transistors following the tests by Schlom et al. 5. However, the investigation of the properties of hafnium is essential for application of this material


4

Page 5 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in the industries. In particular, the thermal behavior of nanoscale hafnium allotropes has not been well studied in the literature. The outline of this paper is as follows. In section 2, the method used to optimize the structures of hafnium cluster, with COMB potential19,20 as interaction potential is described. A novel method for quantifying the melting transition is proposed, dubbed as global similarity index, which is based on chemical similarity principle. Section 3 presented input structures generation, simulation and post-processing steps. The application of global similarity index is extended to gain insights of the melting mechanism in Section 4.

2.

Method

2.1. Structural optimization Basin-Hopping (BH) is an unbiased global optimization method introduced by Wales and Doye21 to locate the ground state or global minimum of an atomic cluster. The complexity of the Potential Energy Surface (PES) determines the complication in obtaining the global minimum structure, which corresponds to the structure with the lowest energy. The objective of BH is to scan and later transform the PES into a simplified staircase topology. Henceforth, it is very efficient and easy in locating the minima. The transformed energy 𝐸̃ (𝑋) is defined by 𝐸̃ (𝑋) = 𝑚𝑖𝑛{𝐸(𝑋)}

(1)

where 𝐸(𝑋) represents a certain point on the PES with 𝑋 being the 3𝑛-dimensional position coordinates of the atoms, {𝑟1 , 𝑟2 , … , 𝑟𝑛 }. {𝑟𝑖 } carries the set of Cartesian coordinate of atom 𝑖 in the form of {𝑥𝑖 , 𝑦𝑖 , 𝑧𝑖 }. BH first generates an atomic configuration, where all atoms involved are 1/3

confined in a sphere of radius 𝑅𝑑 = 𝑟0 [1 + (3𝑛/4𝜋√2)

] with ro being the nearest-neighbor


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 49

distance. Potential energy of the said cluster is then calculated. Canonical Monte-Carlo is applied to search for the above mentioned PES. The BH method was adopted as part of a global search algorithm by Hsu and Lai22, in their fulllength algorithm called the PTMBHGA. The tendency of trapping in saddle point is solved by coupling the multi-canonical BH (MBH) together with Genetic Algorithm (GA) and parallel tempering Monte-Carlo (PT) method, hence the name PTMBHGA. Parallel tempering technique with various weight factors expand the searching pathway in the PES. It is said to be equivalent to the opening a multiple searching pathway, by running several independent Monte Carlos and multiple Basin Hopping in different machines concurrently. The deciding factor in PTMBHGA is the replacement of Boltzmann weighting factor with a multi-canonical weighting factor for facilitating the jumping over large basin and thus avoid the searching to be halted by trapping. Many who study cluster structures have come to realize that the global minimum obtained using classical interatomic potential does not always correspond to the result obtained from DFT. Oganov and Valle23 proposed that there are possibilities that the low-lying structures (LLS) might include the global minimum. Sometimes, the quantum refined global minimum happens to coincide with one of the higher energy local minimum instead of the global minimum of the empirical potential. The method of BH is applicable here due to the purpose of searching the minimum is aimed at locating a set of minima along the PES. During the search for the global minimum structure of a cluster, the effect of electronic charge transfer is not taken in account if less computationally intensive potentials are employed. Hence, the global minimum obtained from empirical interatomic potentials may not always be the one correspond to the real case, i.e. structure calculated quantum mechanically and thus assumed global minimum.


6

Page 7 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. here.

The use of LLS was introduced by Hartke24,25 in an attempt to improve the global minimum search based on ab initio methods such as DFT. The usage of LLS in the global minimum search is illustrated in Figure 1. The final aim is to obtain the true global minimum of a given cluster size at DFT level. In practice, the true global minimum is formidably difficult to be identified directly with limited computational resource. Hence, the overall strategy is to arrive at this true global minimum via an indirect route. In stage 1 of Figure 1, BH method is used to transform the PES into a simplified staircase topology which contains a collection of LLS in different basins. These LLS are indicated as green triangles in Figure 1. The LLS obtained in this stage are not to be minimized into the local minimum of each basin as only a minimal default setting for energy calculation is used (where no powerful minimization algorithm is called) so that the PES can be scanned through at a higher efficiency (at a fixed computational resource constraint). Stage 1 is a slightly faster process as compared to stage 2. In stage 2, the LLS from stage 1 are geometrically re-optimized (opt; indicated as down arrows in Figure 1) into their corresponding local minima using an expensive computational setting. These geometrically re-optimized structures are labelled by red diamonds in Figure 1 which are seen sitting at the local minimum of each basin. These red diamond are originally the LLS at stage 2 level. Geometrical re-optimization process in this strategy provides a mechanism to transform the LLS in stage 1 into the local minimum in each basin. Each basin contains a cluster structure with unique geometry. The lowest energy structure among these LLS is taken as the true global minimum structure.


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 49

2.2. Characterizing melting point for nanocluster Many experimental studies have established the following facts regarding the melting of a free cluster in comparison to its bulk counterpart, namely7 1. The melting temperature, 𝑇𝑚 of a cluster is generally lower than the bulk value, 𝑇𝑚𝑏𝑢𝑙𝑘 . 2. The latent heat of transformation is smaller than the bulk. 3. The melting stage does not occur at a fixed temperature, but begin with pre-melting, 𝑇𝑝𝑟𝑒 over a finite range of temperature. 4. The heat capacity of the finite-sized system can sometimes take a negative value. Statement one and two are analogous, since the melting temperature, 𝑇𝑚 and the latent heat of a cluster show similar fluctuations, while pre-melting is widely observed though experiments. Qiao et al.38 predicted that nanoparticles show pre-melting occurs at some temperature 𝑇𝑝𝑟𝑒 , lower than the actual melting point 𝑇𝑚 , and it begin from surface atoms. A negative heat capacity implies the energy is absorbed with decreasing cluster temperature. Heat absorption during melting is commonly understood in terms of latent heat of transformation, where the mean kinetic energy tends to remain constant. Schmidt and Haberland7 explained that a finite sized system tends to convert some of the kinetic energy into potential energy in order to avoid partially molten states. The melting data obtained via simulation can be used to study the dynamics of phase change. Comparison between commonly used post-processing methods are discussed by Lu et al.26 The most commonly used post-processing method is caloric curves, whereby the binding energy (or binding energy per atom) is plotted against temperature. Another route of measuring melting is through constant temperature specific heat capacity 𝑐𝑣 as a function of temperature, which is actually the fluctuations of the caloric curve,


8

Page 9 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


〈𝐸𝑡2 〉 𝑇 − 〈𝐸𝑡 〉2𝑇 𝑐𝑣 = 2𝑛𝑘𝐵 𝑇 2

(2)

where 𝐸𝑡 is the total energy of cluster, 𝑘𝐵 is the Boltzmann constant, 𝑛 the total number of atoms in the cluster and 〈 〉 𝑇 represent the thermal average at temperature 𝑇. It will be referred to as the 𝑐𝑣 hereafter. A typical 𝑐𝑣 curve appears in the form of a sharp peak at the melting temperature, 𝑇𝑚 . Lu et al.26 showed that there is a mismatch in the melting temperature of Co13 and Co14 clusters, if Lindemann index, 𝛿, is used as a means to gauge the melting process, as compared to 𝑐𝑣 curve. Lindemann index is 𝛿=

1 ∑ 𝛿𝑖 𝑛

(3)

𝑖

1

(〈𝑟𝑖𝑗2 〉 𝑇 − 〈𝑟𝑖𝑗 〉2𝑇 )2 1 𝛿𝑖 = ∑ 〈𝑟𝑖𝑗 〉 𝑇 𝑛−1

(4)

𝑗≠𝑖

where 𝑟𝑖𝑗 is the distance between the 𝑖th and 𝑗th atoms. Through this Lindemann index, structure is seen to deviate from its average form. To a bulk material, melting occurs when the root-meansquare deviation of bond length of the atoms is more than 10% of the equilibrium bond length distance in the ground state structure of the crystal. The dynamics during the cluster melting could also be studied by short-time averaged distance 〈𝜎𝑖 (𝑡)〉𝑠𝑡𝑎 , 〈𝜎𝑖 (𝑡)〉𝑠𝑡𝑎 = ∑|𝑟𝑖 (𝑡) − 𝑟𝑗 (𝑡)|

(5)

𝑗

where 𝑟𝑖 (𝑡) represent the position of the 𝑖th atom at time 𝑡, while 〈 〉sta denotes the average taken for a short interval of time steps, and is plotted against the time. Similar method has been used earlier by Aguado et al.27 on Na cluster. Computing the time-average value is troublesome but


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 49

potentially a solution to the thermal instability in the simulated annealing procedure. The problem of thermal stability is a big obstacle in simulation of cluster melting.

2.3. Chemical similarity and shape recognition The method of molecular shape comparison is used to comparing the geometry or spatial configuration of two or more molecular structures.28 The qualitative analysis is based on the similar property principle of molecule, which states that similar compound have similar properties. The selected structures for shape identification are normally represented by a binary molecular fingerprint of descriptors, also known as structural keys. Some examples include the size of the molecule, the number of bonds or type of bonding, the active functional groups and the pattern of target structure or substructure. The descriptors are generally classified into three types based on the dimensionality. The two-dimensional descriptors such as MACCS, MDL keys, and Daylight are mostly patented under their own company signature. Three-dimensional descriptors were difficult to compute until Ballester and Richards29 proposed the ultrafast shape recognition (USR) method based on four sets of distance distribution of the atoms in a molecule defined at different points of reference. The geometries of the atomic configurations are described by a total of three statistical moments, namely the mean, variance and skewness. It turns out that these USR descriptors are orientation independent, as the screened molecules do not need to be aligned to each other in order to obtain the similarity score. USR is said to perform at least three order of magnitude faster than other descriptors. The similarity index is denoted as 𝑆𝑞𝑖 ∈ (0,1), with the value of 1 represents “totally identical” while 0 means “vastly differed”. 𝑆𝑞𝑖 is defined as 12

1 𝑆𝑞𝑖 = (1 + ∑|𝑀𝑙𝑞 − 𝑀𝑙𝑖 |) 12

−1

(6)

𝑙=1


10

Page 11 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where the moments of shape descriptors 𝑀𝑞 and 𝑀𝑖 represent the query and the 𝑖th molecule.

3.

Pre-simulation preparations

3.1. Generating structures In order to accurately locate the global minimum of a cluster, parallel tempering multi-canonical basin hopping plus genetic algorithm (PTMBHGA), developed by Lai et al.30, is adopted. The searching algorithm contains two global optimizers that complement each other. The basin hoping (BH) scans over the PES and generates 20 parent candidates around the minima of the PES. Local minimization is carried out before Genetic Algorithm (GA) discards the top five clusters with the highest total energies and generates new candidates to form the next-generation candidates. The process is then repeated. GA serves to complement BH by introducing new offspring located in different basins to avoid the trappings in saddle point. This BH-GA process will cease once the number of similar structure, all having same local minimum energy value, achieved a predefined value. Overall, the fitness value 𝑓𝑖 of the candidates is evaluated during all stages as a normalized control based on the equation 𝐹𝑖 20 ∑𝑗=1 𝐹𝑗

(7)

𝑉𝑚𝑎𝑥 − 𝑉𝑖 𝑉𝑚𝑎𝑥 − 𝑉𝑚𝑖𝑛

(8)

𝑓𝑖 =

𝐹𝑖 =

where 𝑉𝑖 is the potential energy of the cluster 𝑖 calculated using COMB potential. The PTMBHGA code comprised of two independent parts, the global-search algorithm part (comprised of GA and BH search algorithm) and an ‘energy calculator’ part. The COMB potential with charge optimization from the energy calculator LAMMPS31 is used in evaluating the energies


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 49

of the cluster generated by the GA-BH algorithm The PTMBHGA + LAMMPS hybrid package is dubbed PTMBHGA-MD hereafter.

3.2. The PTMBHGA-MD process line Candidate trial structures are generated by the PTMBHGA algorithm, followed by the calculations of their corresponding total energy employing the COMB potential implemented within the LAMMPS package. All the trial structures undergo relaxation using LAMMPS package at 𝑇 = 1 K,whereby each hafnium atomic cluster is allowed to equilibrate with minimal vibration. The resultant structure corresponds to the global minimum of the COMB PES. In the second stage, the clusters are locally and geometrically re-optimized with DFT package, Gaussian 0332 (dubbed G03), through implementation of two-electron integral computation for the purpose of high accuracy calculation. The quadratically convergent self-consistent field (SCF) procedure33 is activated in this stage whereas the basis set B3LYP/LanL2DZ is used as suggested by Sun et al.34. Kuntová et al.15 asserted that the similarity in geometry is the key to relate whether the structures of different PES lie in the same basin. Comparing the chemical similarities of the structures lying in empirical potential PES and that in DFT PES is a convenient way to computationally access the correspondence between two PES, specifically, whether both PES share a common basin. It is hypothesized that the structures obtained through PTMBHGA-MD lies within the basin containing the global minima of the DFT PES. Additional tests are carried out to justify the hypothesis by systematically comparing the structures obtained from the quantum refinement of PTMBHGAMD structures and the global minimum structures obtained from DFT approach.


12

Page 13 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1. here. Table 1 layout the comparison of structures obtained through PTMBHGA-MD and those underwent DFT geometrical re-optimization. Judging from the point groups, structures from both procedures are mostly the same or bear high similarity, thus justifying the use of the PTMBHGAMD method in generating structures of global minimum in nature for further melting purpose. For reference, as the index value closer to 1, i.e. 𝜉𝑖 → 1, the more similar both structures they are.

3.3. Melting Simulated annealing process refers to heating and quenching procedures in MD. In this work, the candidate structures are heated from 1 K to the desired target temperature and equilibrated at the same temperature for an extended period of time, during which the simulated results are sampled for statistical analysis. Simulated annealing process is performed using LAMMPS package, with the COMB potential being used to describe the interactions among hafnium atoms. Each hafnium cluster is the ground state structure obtained from the previous PTMBHGA-MD method. The free standing hafnium atomic cluster is positioned at the center of a simulation box, which is fixed with reflective walls at 100Å in each direction from the origin. The simulation timestep Δt is taken to be 0.5 fs. The temperature control is carried out using canonical (NVT) Nose-Hoover thermostat.


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 49

Relaxation at 0~1 K is carried out to give the hafnium cluster some good initial configuration of microstates. Next, the heating rate is selected to ensure the ergodicity of the system during melting, in which the size of the cluster is taken into. The choice of heating rate follows the criteria: 1. Optimum heating rate should give a 𝑐𝑣 curve that clearly shows a consistent range of melting point for clusters of all sizes (applied to both 𝑇𝑚 and 𝑇𝑝𝑟𝑒 ). 2. The computational time in heating the system to the desired target temperature is not over. 3. The pre-melting behavior can be captured. 4. The range of melting point is almost the same for the cluster of neighboring sizes.

Figure 2. here.

Figure 2 shows the resultant plot of melting and pre-melting points for Hf50 against the heating rate. Apparently, the first three statement converge once the heating rate come close to 7 × 1012 Ks −1 (the melting points become constant value, as shown in the marked area of Figure 2). It is found that, for smaller cluster size 𝑛, the difference in melting point tends to be larger compared to the neighboring size. But the forth statement is met even for higher heating rate, up to > 10 × 1012 Ks −1 . The heating rate is set at slightly lower than the converge point as indicates by the red arrow in Figure 2, which is 5 × 1012 Ks −1 (≡ 400 𝛥𝑡K −1). As a rule of thumb, the melting temperature, 𝑇𝑚 , is expected to increase with the number of atoms in the clusters, 𝑛. Due to the surface-to-body-atom ratio, 𝑇𝑚 will eventually saturate at certain size, whereby thermal behaviors starts to resemble those of bulk value. It is beyond the scope of this paper to go as big as 104 in size, but the size effect is being studied up to 𝑛 = 50. An additional cluster of size 𝑛 = 99 is designed to compare the results.


14

Page 15 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.4. Post-processing The widely accepted method to study the thermal properties of any heating process is by plotting the caloric curve and the 𝑐𝑣 vs. temperature curve. In this paper, both curves are plotted to obtain the melting temperature of hafnium clusters. The melting temperature is then compared to that obtained via the global similarity index, which will be introduced in the next subsection. There are two ways to obtain the caloric profile, which is through ‘direct heating’ and ‘prolonged annealing’. In the direct heating procedure, the cluster is heated directly to a high temperature with a slow heating rate. To ensure the melting dynamics is well represented by caloric profile, the heating rate used is 5 × 1012 Ks−1 . A slow hearing rate can minimize the error accumulated due to rapid velocity rescaling at every time step. In the second approach, each point on the caloric curve is taken from a single annealing process, which means separate runs for each interim temperature point on the curve. In this procedure, the candidate structure is heated to the desired interim temperature points and then allowed to equilibrate for a very long time to ensure the ergodicity of the microstate, specifically the temperature and the potential energy value. Both methods described above are used in this paper to see whether both approaches agree well with each other. The dynamical information of the simulation is traced by comparing the similarity index of each frame along the heating process. The comparisons are made with respect to the initial configuration at 1 K (after relaxation) prior to the heating process. The trajectories are recorded at every 300𝛥𝑡, which corresponds to 0.15 ps of simulation time. With the heating rate of 5 × 1012 Ks −1 , it is expected to have at least one record of the simulated system per unit temperature.


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 49

3.5. Global similarity index The global similarity index, 𝜉𝑖 proposed in this paper is a novel approach and derived based on the generic chemical similarity idea in order to detect the changes along the trajectory of the heating process. It is used as a measure to detect changes in the structures during the heating procedure. Compared to all the existing methods, global similarity index is capable of predicting detailed melting mechanism of the system. The functional form of 𝜉𝑖 is proposed to take the form of 𝑛

1 −1 𝜉𝑖 = ∑(𝑘𝑠,𝑖 + 1) 𝑛

(9)

𝑘𝑠,𝑖 = |√𝑑𝑠,𝑖 − √𝑑𝑠,0 |

(10)

𝑠=1

where 𝑑𝑠,𝑖 and 𝑑𝑠,0 represent the sorted distance of atoms relative to the average positions (center of mass) in the cluster for the 𝑖th and 0th frame (the input structure for the simulated annealing process), while 𝑛 corresponds to the number of atoms, which is equals to the number of pairs of 𝑘𝑠,𝑖 . The value of 𝜉𝑖 = 1 corresponds to totally identicalness and 𝜉𝑖 → 0 for vast difference. Instead of just a single, conventionally defined center of mass, this work used three different generalized mean positions to act as the center of mass of a cluster, namely the arithmetic mean, the harmonic mean, and the quadratic mean (a.k.a. root mean square). Each of these means is viewed as a center of reference (COR) for the immediate state of clusters. The idea of COR works like the center of mass, where the whole cluster is represented by a single point mass. With this generalization, there are now three COR for each frame of cluster. The generalized mean, also known as the power mean, can be expressed in the following form 𝑛

1 𝑀𝑝 (𝑥1 , … , 𝑥𝑛 ) = ( ∑ 𝑥𝑖𝑝 ) 𝑛

1 𝑝

(11)

𝑖=1


16

Page 17 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where 𝑝 must be a non-zero real number, and 𝑥1 , … , 𝑥𝑛 are variables. When 𝑝 takes on the value 𝑝 = 1, −1, 2, we have arithmetic mean, harmonic mean, and quadratic mean respectively. The cluster will be quantified by three different sets of 𝑘𝑠,𝑖 , each based on a different COR. The definitions of Eq.(9) and (10) are hence generalized to include an COR index to indicate which COR is being referred to when evaluating 𝜉𝑖 , 𝑛

1 −1 𝐶𝑂𝑅 = ∑(𝑘𝑠,𝑖 + 1) 𝑛

(12)

𝐶𝑂𝑅 𝐶𝑂𝑅 𝐶𝑂𝑅 𝑘𝑠,𝑖 = |√𝑑𝑠,𝑖 − √𝑑𝑠,0 |

(13)

𝜉𝑖𝐶𝑂𝑅

𝑠=1

In each COR, there are 𝑛 sorted atomic distances, giving a total of 3𝑛 pairs of 𝑘𝑠,𝑖 . The distances of atoms from the mean position are sorted in ascending order from the shortest to the longest. The difference in Eq. (12) and (13) are measured with respect to the distance of atoms in the 1 K frame of the same sorting sequence. After the value of 𝜉𝑖 is calculated for each COR, averaging of 𝜉𝑖 over these COR is then obtained, via 𝜉𝑖̅ =

1 𝑁𝐶𝑂𝑅

𝑛

∑ 𝐶𝑂𝑅

𝜉𝑖𝐶𝑂𝑅

1 1 −1 𝐶𝑂𝑅 = ∑ ∑(𝑘𝑠,𝑖 + 1) 3 𝑛 𝐶𝑂𝑅

𝑠=1

𝑛

1 −1 𝐶𝑂𝑅 = ∑ ∑(𝑘𝑠,𝑖 + 1) 3𝑛

(14)

𝐶𝑂𝑅 𝑠=1

where 𝑁𝐶𝑂𝑅 = 3, the number of COR used in this paper. It is found that taking the square of Eq. (14) improves the sensitivity of 𝜉𝑖 in detecting similarity in the clusters. Since a 𝜉𝑖 ranges from 0 → 1, squaring Eq. (14) does not change the normalization property of 𝜉𝑖 . In fact, by squaring Eq. (14), the variation in 𝜉𝑖 2 becomes numerically more pronounced and sensitive due to the parabolic behavior of squaring a number within the range of


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 49

(0,1). In practice, the value of 𝜉𝑖 tends to vary only very slightly. Hence, squaring the 𝜉𝑖 is a convenient way to blow up the visibility of the variation effect. It is instructive to compare Eq.(14) to the USR defined by Ballester and Richards29 12

1 𝑆𝑞𝑖 = (1 + ∑|𝑀𝑙𝑞 − 𝑀𝑙𝑖 |) 12

−1

𝑙=1

𝑆𝑞𝑖 as defined in the definition of USR has a very different working principle from that of the global similarity index. 𝑆𝑞𝑖 has the general form of 𝑆𝑞𝑖 ~

1 1 1 + 12 ∑12 𝑙=1(⋯ )𝑙

, (⋯ )𝑙 ≡ |𝑀𝑙𝑞 − 𝑀𝑙𝑖 |

(15)

whereas 𝜉𝑖 has the general form 𝑛

1 1 𝐶𝑂𝑅 𝐶𝑂𝑅 𝜉𝑖 ~ ∑∑ , (⋯ )𝑠 ≡ |√𝑑𝑠,𝑖 − √𝑑𝑠,0 | (⋯ ) 3𝑛 1+ 𝑠

(16)

𝐶𝑂𝑅 𝑠=1

The fraction of

1 12

in 𝑆𝑞𝑖 is sitting inside the reciprocal, whereas the fraction

1 3n

in 𝜉𝑖 is placed

𝐶𝑂𝑅 outside the reciprocal. This will ensure that 𝑘𝑠,𝑖 is not averaged before the index 𝑠 is summed

over. As it was originally proposed by Ballester and Richards29, the 𝑀𝑙 terms as appear in 𝑆𝑞𝑖 have to undergo a prescribed normalization procedure so that they give values with the same order of magnitude for all 𝑙 (𝑀𝑙 tends to give values with different order of magnitude for different 𝑙 if not 𝐶𝑂𝑅 normalized). On the other hand, normalization procedure is not required for 𝜉𝑖 as each 𝑘𝑠,𝑖 , being

the differences in the distances of atoms from the COR, has the same order of magnitude. The main idea of the global similarity index is to compare the configuration of one image of the molecule (cluster) to another. Imagine two images of identical cluster with one of them is obtained by performing a linear transformation (including translation, rotation, reflection and rescaling) on the other. The global similarity index must be able to discern the fact that these two images, despite


18

Page 19 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


having different coordinates, are in fact the same cluster, and return a unique value of 𝜉𝑖 = 1. In general, when two images are being compared, the transformational relation between them have to be factored in while evaluating the global similarity index. The measurement protocol of the similarity index 𝜉𝑖 can be further studied by computing the fluctuation of 𝜉𝑖 over the course of simulation. In general, 𝜉𝑖 , which is related to the microscopic property, is not directly measurable via experiment. However, its fluctuation could be related to a macroscopic observable. As in the case of 𝑐𝑣 , it is a measure of fluctuation in the binding energy. The fluctuation of 𝜉𝑖 remains a prospect for future study. The fluctuation equation is considered to be proportional to the variance of 𝜉𝑖 as 𝑆𝜉𝑖 ∝ 𝜎 2 (𝜉𝑖 )

(17)

The 𝑆𝜉𝑖 curve shows a sudden change in the similarity index 𝜉𝑖 in the form of sharp peaks at the particular moment when major transition occurs in the configuration.

4.

Application and discussion

The global similarity index 𝜉𝑖 is used to study the dynamics of hafnium clusters during the simulated annealing processes. The value of 𝜉𝑖 in each instance 𝑖 is a measure of similarity of the cluster structure in that instance as compared to the ground state structure (the structure during the initial 1 K equilibration) in the 0th frame. During the early part of the heating-up process, only little variation in 𝜉𝑖 is detected, indicating that only minor change to the geometry happens in the cluster. Up to this point, 𝜉𝑖 tends to measure an almost constant value close to 1. The abrupt change to the global similarity index will only happen during (or close to) the melting temperature. Visually, the melting transition happens with a total distortion in the structure spanning the whole cluster. Following the sequential evolution of 𝜉𝑖 allows one to locate the sudden change in the


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 49

structure of the clusters during a heating process. After the melting point, the whole cluster is expected to be fluid-like and each atom tends to roam across the simulation box randomly. Thus, the global similarity index is expected to be constantly close to 0. Note that, 𝜉𝑖 stays almost constant during most of the simulated annealing process, except during the range of temperature in which the melting process occurs. Thus, the fluctuation of similarity index 𝑆𝜉𝑖 becomes useful to locate the temperature where these transition happened. In addition to the melting transition, patterns of distortion to the clusters during the simulated annealing process are also detectable by the global similarity index based on the fluctuation curve 𝑆𝜉𝑖 . To begin, we discuss the MD heating simulation of Hf20 as an illustrative example. The caloric and 𝑐𝑣 curves of Hf20 generated by prolonged annealing process are presented in Figure 3. Each point on the caloric curve has a corresponding counterpart on the 𝑐𝑣 curve at the same temperature. It is easy to understand the connection between the two graphs based on the thermodynamic interpretation for the constant volume specific heat capacity, cv=(∂E/∂T)v, where the subscript v denotes the derivatives taken under constant volume thermosetting. However, here we use the fluctuations in the caloric curve to obtain 𝑐𝑣 curve.

Figure 3. here, a) and b).

Melting transition is expected to happen when the potential energy of the cluster changed drastically during the melting point, 𝑇𝑚 . This temperature corresponds to the latent heat of fusion. Based on the definition of cv=(∂E/∂T)v, drastic change of 𝐸 value is expected to return a relatively larger value of 𝑐𝑣 in the vicinity of 𝑇𝑚 . Thus, melting transition corresponds to a large peak in 𝑐𝑣 curve. Apparently, from Figure 3, the caloric curve is less effective in giving the exact temperature


20

Page 21 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for both the melting point and pre-melting point as compared to the 𝑐𝑣 curve. There are many random peaks in the 𝑐𝑣 curve (Figure 3b), the first apparent peak is the pre-melting, while the highest peak is the melting point. The melting point corresponds to the most obvious transition in the caloric curve (Figure 3a). The green arrow indicates the pre-melting temperature 𝑇𝑝𝑟𝑒 = 1400 K while the red arrow indicates the melting temperature 𝑇𝑚 = 1850 K in both graphs. In addition to the prolonged annealing method, an independent method to obtain melting behaviour of the clusters is also attempted. The purposes for conducting this procedure on top of the prolonged method are two-fold: (1) it is used as an independent check to the results obtained via the prolonged method, and (2) to generate the evolutional history of global similarity index of the clusters, which carries the dynamical information throughout the evolutionary history of a MD simulation. The trajectory file from the LAMMPS simulation contain the frames of the atoms that reflect the heating process as recorded in the corresponding log file. However, in the prolonged annealing simulation, in which the target temperatures are generally set to a very high value, e.g. 𝑇 = 3000 K, a large fluctuation in the energy value is usually present throughout the equilibration steps. This large fluctuation is mainly due to the finite size effects of the cluster system35. The system is too small in the vicinity of the simulation box and the fact that the cluster in study is not periodic has caused thermal fluctuations in the system studied. Due to the large fluctuation in the temperature, the trajectory of the atoms also fluctuates rigorously. Over an extended equilibration time step, the trajectory of each individual hafnium atom in the cluster will move all over the simulation box and tends to be ergodic. The ergodicity is not desired as the average position of each atom is seen to overlap in the confined space of simulation box. Therefore, global similarity index cannot be used in the prolonged heating procedure to quantify the heating process. In


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 49

contrast, the direct heating method, due to the imposition of the 300∆𝑡 time averaging protocol, the fluctuations of all recorded data are confined within the 300∆𝑡 time window.


On the other hand, the direct heating method is also able to predict the melting temperature 𝑇𝑚 and pre-melting temperature 𝑇𝑝𝑟𝑒 by plotting the caloric curve and the 𝑐𝑣 curve. In this respect, the caloric and 𝑐𝑣 curve are known as ‘indicators’, meaning molecular dynamics quantities from which thermodynamically interesting transitions can be monitored. Figure 4 shows the caloric curve and 𝑐𝑣 curve of Hf10 via direct heating method. The 𝑐𝑣 curve shows a lot of peaks (Figure 4b). The first apparent peak is the pre-melting and the tallest peak is the melting point. The melting point corresponds to the largest transition in caloric curve (Figure 4a). By using the same color code as for Figure 3, green arrow indicates the pre-melting temperature 𝑇𝑝𝑟𝑒 = 1350 K, and the red arrow indicates the melting point 𝑇𝑚 = 2200 K. Both curves provide complementary information to help locating the transition temperatures. Only one unique transition temperature is reported based on these two indicators, after comparing and contrasting both curves to estimate the most precise value for the transition temperatures. Due to its higher sensitivity, the transition temperatures reported in this study are in all cases read off from the 𝑐𝑣 curve. As it turns out, indication of transition in the caloric curve is relatively less distinctive (hence higher in uncertainty), hence it is used only as a supplementary confirmation to the transition value reported from the 𝑐𝑣 curve. As a comparison to the cluster system, a bulk hafnium system is presented under the same temperature range and heating rate. Figure 5 shows the caloric curve and 𝑐𝑣 curve of hafnium crystal via direct heating method. The melting point of the hafnium crystal is in the range of


22

Page 23 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


𝑇𝑚𝑏𝑢𝑙𝑘 (𝐻𝑓)~2600 𝐾 to 2650 K. As compared to the cluster samples shown in Figure 3 and Figure 4, the simulation with hafnium crystal clearly shows that fluctuations in energy value are greatly reduced.

Figure 5. here, a) and b)

The identification of transition temperatures can also be made based on another independent indicator, which is the global similarity index. The trajectory file used to obtain the similarity index profile is based on the direct heating method. The dynamical evolution of global similarity index as a function of temperature is shown in Figure 6 for Hf13 as a sample candidate. Figure 6a shows the similarity index 𝜉𝑖 profile for each frames during the heating process as compared to the ground state structure (the 0th frame), derived using Eq. 12. Figure 6b is the fluctuation (variance) of similarity index 𝑆𝜉𝑖 corresponds to the same trajectory calculated by using Eq. 17.


The features of thermodynamically interest in the example of Figure 6 is now interpreted. Melting and pre-melting features can be identified from both the 𝑆𝜉𝑖 vs 𝑇 and 𝜉𝑖 vs 𝑇 curves, albeit with different sensitivity. The pre-melting transition is not obvious in the 𝑆𝜉𝑖 curve due to the suppression of the overall scale in 𝑆𝜉𝑖 by the high peak at 𝑇 ≈ 2050 K. The pre-melting transition can merely become visible by zooming in to the 𝑆𝜉𝑖 curve at around 𝑇 ≈ 1050 K, as displayed by the inset of Figure 6b. As is observed throughout the MD calculations, it is found that usually the signal of pre-melting or melting in the 𝑆𝜉𝑖 curve is only moderate unless the transitions in the


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 49

clusters occurs in a really abrupt manner. In practice, both the 𝑆𝜉𝑖 and 𝜉𝑖 curves have to be simultaneously deployed so that the complementary features in both curves are able to pinpoint a transition with much a reliable precision. As in the case of Figure 6, the pre-melting point correspond to the first obvious drop in the similarity index (indicated by the green arrow) at 𝑇𝑝𝑟𝑒 = 1600 K. On the other hand, the melting point is visibly apparent in both graphs (indicated by the red arrow) at 𝑇𝑚 = 2050 K. Three different indicators are used in this study to identify melting transitions in two independent heating procedures. In the prolonged heating method, only two indicators are used to identify melting transitions, namely the caloric and 𝑐𝑣 curves. In the direct heating method, the caloric, the 𝑐𝑣 curves and global similarity index are used as indicators. In this section, the transition temperatures as a function of cluster size 𝑛 are reported.

Table 2. here.

Table 2 shows the results of melting and pre-melting temperatures in both heating procedures. The first two columns are the melting temperatures obtained by prolonged annealing methods. The third and fourth columns are directly extracted from a direct heating simulation elevated from 1 K to 3000 K at the optimum rate of 5 × 1012 Ks −1 . These transition temperatures are determined from the 𝑐𝑣 curve but not the caloric curve. Meanwhile the last two are transition temperatures deduced from the graph of fluctuations in the global similarity index, 𝑆𝜉𝑖 . Recall that global similarity index is only applicable in the direct heating process. In fact, the prolonged simulated annealing method involves compilation of many heating processes equilibrated at different target temperatures. Every single data point in the caloric curve is the result of statistically sampling a


24

Page 25 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cluster’s energy over a lengthy period of equilibration at the fixed temperature. The pre-melting temperature, 𝑇𝑝𝑟𝑒 and the exact melting temperature, 𝑇𝑚 , are plotted in the same graph in Figure 7 for each method. For the sake of easy comparison, the melting temperature 𝑇𝑚 vs 𝑛 obtained from the three different approaches listed in Table 2 are compiled into a single graph as shown in Figure 8. Results for the clusters Hf60 through Hf99 are used to extrapolate the graph for bigger cluster sizes.

Figure 7. here, a), b) and c).

As mentioned in Section 3.5, direct heating of the structure to a high temperature allows one to yield detailed dynamics of the clusters. It can be seen from Figure 8 that this approach overall predicts a higher melting point than that from the prolonged annealing method, but the prediction still agrees well with expectation as compared to the bulk melting of hafnium, which is 𝑇𝑚𝑏𝑢𝑙𝑘 (𝐻𝑓)~2504 𝐾 according to the periodic table, and 𝑇𝑚𝑏𝑢𝑙𝑘 (𝐻𝑓)~2600 𝐾 according to current work as shown in Figure 5. Nevertheless, both simulation approaches agreed that the melting temperature of hafnium cluster increases with the size 𝑛.

Figure 8. here.

One of the most important results obtained from this study is that global similarity index, 𝜉𝑖 , has the capacity to predict the dynamics of the heating process from the MD trajectory of clusters. In general, the pre-melting effect take place at the first significant peak of 𝑆𝜉𝑖 , whereby the exterior of the cluster is usually heavily distorted, but the constituent atoms are still bounded relatively


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 49

close to each other. As the heating continues, the cluster tends to tear apart into fragmented state and start drifting in random direction, entering the fluidic phase. In order to portray the distance between atoms in the clusters, the dynamic bond length is set to be around 3.127Å. The following subsections focus only to a few specific sizes of hafnium clusters.

4.1. Cluster Hf30 From Figure 9, the first significant change of Hf30 cluster in the heating process occurs around 𝑇 = 1000 K to 1300 K, indicated by the sudden drop of 𝜉𝑖 . The outermost layer atoms become further apart from each other but remain intact with the ‘core’ atoms. The pre-melting takes place at the first peak of the 𝑆𝜉𝑖 , at around 𝑇𝑝𝑟𝑒 = 1750 K with the ‘shell’ atoms begin to drift apart from the cluster. The highest peak of 𝑆𝜉𝑖 , corresponds to a temperature around 𝑇𝑚 = 1850 K, showing a total breakdown of the cluster. A final screenshot at 𝑇 = 2150 K shows the remnant cluster ‘core’ after many surface atoms have drifted away.


4.2. Cluster Hf50 For the case of Hf50 as shown in Figure 10, the 𝜉𝑖 curve displays some unique features, which are the sudden drop at around 𝑇 = 1850 K and a sudden rise at around 𝑇 = 1900 K. This behavior is observed to be the effect of a single hafnium atom leaving the cluster drifting to the wall of the simulation box and bouncing back to the cluster. We could see this behavioral change as a result of kinetic energy accumulated within a single atom that causes it to escape from the cluster. In


26

Page 27 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


order to ensure that this artifact is not some kind of random error in the simulation, the overall simulated annealing procedure is re-run with slightly different initial condition (different random initial velocity). The same feature still shows up when different initial condition is imposed. The sharp change is also detected by the fluctuation plot against temperature, showing two sharp peaks at about the same temperature. The left one corresponded to a sudden drop in 𝜉𝑖 while the right one corresponds to a sudden rise of 𝜉𝑖 value. The true melting is at around 𝑇𝑚 = 2150 K when the cluster breaks down completely and drifts apart.


4.3. Cluster Hf99 Cluster size Hf99 is used to describe the complete trace of dynamics during the single annealing process from 𝑇 = 0 K to 3000 K, as shown in Figure 11. The atomic vibration slowly increases in magnitude upon velocity rescaling. As observed from the 𝜉𝑖 curve, every significant drop of the 𝜉𝑖 value represents significant changes to the cluster geometry, until 𝑇𝑝𝑟𝑒 = 1900 K where small groups of dimers and trimers begin to tear apart from the core structure, marking the surface premelting effect. At 𝑇 = 2000 K, the cluster geometry is completely distorted and becomes nonspherical. The distortion of the whole cluster continues till 𝑇𝑚 = 2300 K before more atoms being drifting away into the vacuum (empty spaces) of the simulation box.



27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

5.

Page 28 of 49

Conclusions

Melting transition of hafnium clusters for the size of clusters up to 𝑛 = 99 is reported. The melting point of hafnium cluster shows gradual increment with the increase in the number of atoms, 𝑛, but remains lower than the bulk value from periodic table, which is 𝑇𝑚𝑏𝑢𝑙𝑘 (𝐻𝑓)~2504 𝐾. While simulation predicted the bulk melting point at some higher value of 𝑇𝑚𝑏𝑢𝑙𝑘 (𝐻𝑓)~2600 𝐾. In addition to that, hafnium clusters experience surface atom pre-melting effect at some temperature 𝑇𝑝𝑟𝑒 lower than that of their melting point 𝑇𝑚 . Besides the conventional methods of caloric curve and 𝑐𝑣 curve, the global similarity index 𝜉𝑖 successfully predicted the temperature for melting and pre-melting transition in hafnium clusters. The global similarity index also enabled the tracking of unique structural changes happened during the heating process in the direct heating procedure.

ACKNOWLEDGMENT This work was supported financially by FRGS Grant (FASA2/2013) from the Ministry of Higher Education, Malaysia. S.K. Lai and Peter Yen are gratefully acknowledged for their original PTMBHGA algorithm.

REFERENCES 1. Ferrando, R.; Jellinek, J.; Johnston, R. L., Nanoalloys: from theory to applications of alloy clusters and nanoparticles. Chemical reviews 2008, 108, 845-910.


28

Page 29 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2. Son, S. U.; Jang, Y.; Park, J.; Na, H. B.; Park, H. M.; Yun, H. J.; Lee, J.; Hyeon, T., Designed synthesis of atom-economical Pd/Ni bimetallic nanoparticle-based catalysts for sonogashira coupling reactions. Journal of the American Chemical Society 2004, 126, 5026-5027. 3. Park, J.-I.; Cheon, J., Synthesis of “solid solution” and “core-shell” type cobalt-platinum magnetic nanoparticles via transmetalation reactions. Journal of the American Chemical Society 2001, 123, 5743-5746. 4. Giasuddin, A.; Jhuma, K.; Haq, A. M., Use of gold nanoparticles in diagnostics, surgery and medicine: a review. Bangladesh Journal of Medical Biochemistry 2013, 5, 56-60. 5. Berry, R. S., When the Melting and Freezing Points Are Not the Same. Scientific American 1990. 6. Schmidt, M.; Kusche, R.; von Issendorff, B.; Haberland, H., Irregular variations in the melting point of size-selected atomic clusters. Nature 1998, 393, 238-240. 7. Schmidt, M.; Haberland, H., Phase transitions in clusters. Comptes Rendus Physique 2002, 3, 327-340. 8. Couchman, P., The Lindemann hypothesis and the size dependence of melting temperatures. II. Philosophical Magazine A 1979, 40, 637-643. 9. Aguado, A.; Jarrold, M. F., Melting and freezing of metal clusters. Annual review of physical chemistry 2011, 62, 151-172. 10. Wales, D. J.; Berry, R. S., Melting and freezing of small argon clusters. The Journal of chemical physics 1990, 92, 4283-4295.


29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 49

11. Duan, H.; Ding, F.; Rosén, A.; Harutyunyan, A. R.; Curtarolo, S.; Bolton, K., Size dependent melting mechanisms of iron nanoclusters. Chemical physics 2007, 333, 57-62. 12. Koga, K.; Ikeshoji, T.; Sugawara, K.-i., Size-and temperature-dependent structural transitions in gold nanoparticles. Physical review letters 2004, 92, 115507. 13. Solliard, C.; Flueli, M., Surface stress and size effect on the lattice parameter in small particles of gold and platinum. Surface Science 1985, 156, 487-494. 14. Weissmüller, J., Comment on “Lattice contraction and surface stress of fcc nanocrystals”. The Journal of Physical Chemistry B 2002, 106, 889-890. 15. Kuntová, Z.; Rossi, G.; Ferrando, R., Melting of core-shell Ag-Ni and Ag-Co nanoclusters studied via molecular dynamics simulations. Physical Review B 2008, 77, 205431. 16. Bottani, C.; Bassi, A. L.; Tanner, B.; Stella, A.; Tognini, P.; Cheyssac, P.; Kofman, R., Melting in metallic Sn nanoparticles studied by surface Brillouin scattering and synchrotron-x-ray diffraction. Physical Review B 1999, 59, R15601. 17. Chuang, F.-c.; Wang, C.; Öğüt, S.; Chelikowsky, J. R.; Ho, K., Melting of small Sn clusters by ab initio molecular dynamics simulations. Physical Review B 2004, 69, 165408. 18. Steenbergen, K. G.; Gaston, N., First-principles melting of gallium clusters down to nine atoms: structural and electronic contributions to melting. Physical Chemistry Chemical Physics 2013, 15, 15325-15332. 19. Shan, T.-R.; Devine, B. D.; Kemper, T. W.; Sinnott, S. B.; Phillpot, S. R., Chargeoptimized many-body potential for the hafnium/hafnium oxide system. Physical Review B 2010, 81, 125328.


30

Page 31 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


20. Shan, T.-R.; Devine, B. D.; Hawkins, J. M.; Asthagiri, A.; Phillpot, S. R.; Sinnott, S. B., Second-generation charge-optimized many-body potential for Si/SiO 2 and amorphous silica. Physical Review B 2010, 82, 235302. 21. Wales, D. J.; Doye, J. P., Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms. The Journal of Physical Chemistry A 1997, 101, 5111-5116. 22. Hsu, P.; Lai, S., Structures of bimetallic clusters. The Journal of chemical physics 2006, 124, 044711. 23. Oganov, A. R.; Valle, M., How to quantify energy landscapes of solids. The Journal of chemical physics 2009, 130, 104504. 24. Hartke, B., Global geometry optimization of clusters guided by N-dependent model potentials. Chemical physics letters 1996, 258, 144-148. 25. Hartke, B., Global geometry optimization of small silicon clusters at the level of density functional theory. Theoretical Chemistry Accounts 1998, 99, 241-247. 26. Lu, S.; Zhang, J.; Duan, H., Melting behaviors of CoN (N= 13, 14, 38, 55, 56) clusters. Chemical Physics 2009, 363, 7-12. 27. Aguado, A.; López, J. M.; Alonso, J. A.; Stott, M., Melting in large sodium clusters: An orbital-free molecular dynamics study. The Journal of Physical Chemistry B 2001, 105, 2386-2392.


31


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 49

28. Grant, J. A.; Gallardo, M.; Pickup, B. T., A fast method of molecular shape comparison: A simple application of a Gaussian description of molecular shape. Journal of computational chemistry 1996, 17, 1653-1666. 29. Ballester, P. J.; Richards, W. G. Ultrafast shape recognition for similarity search in molecular databases. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2007; The Royal Society: 2007; Vol. 463; pp 13071321. 30. Lai, S.; Hsu, P.; Wu, K.; Liu, W.; Iwamatsu, M., Structures of metallic clusters: Mono-and polyvalent metals. The Journal of chemical physics 2002, 117, 10715-10725. 31. Plimpton, S., Fast parallel algorithms for short-range molecular dynamics. Journal of computational physics 1995, 117, 1-19. 32. Frisch, A., Gaussian 03. Gaussian: 2004. 33. Bacskay, G. B., A quadratically convergent Hartree—Fock (QC-SCF) method. Application to closed shell systems. Chemical Physics 1981, 61, 385-404. 34. Sun, X.; Du, J.; Zhang, P.; Jiang, G., A Systemic DFT Study on Several 5d-Electron Element Dimers: Hf2, Ta2, Re2, W2, and Hg2. Journal of Cluster Science 2010, 21, 619636. 35. Robbins, M. O.; Grest. G. S.; Kremer, K., Effect of finite system size on thermal fluctuations: Implications for melting. Physical Review B 1990, 42, 5579-5585.


32

Page 33 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


36. Dugourd, Ph.; Rayane, D.; Labastie, P.; Vezin, B.; Chevaleyre, J.; Broyer, M., Measurements of lithium cluster ionization potentials. Chemical Physics Letters 1992, 197, 433-437. 37. Knight, W. D.; Walt A. de Heer; Clemenger, K.; Saunders, W. A., Electronic shell structure in potassium clusters. Solid State Communications 1985, 53, 445-446. 38. Qiao, Z.; Feng, H.; Zhou, J., Molecular dynamics simulations on the melting of gold nanoparticles. Phase Transitions 2014, 87, 59-70.


33


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 49

List of figures:

Figure 1. A schematic sketch indicating the strategy to obtain true global minimum by the way of sampling LLS at a coarse level search with BH method without structural optimization (stage 1), and subsequently undergo a refined geometrical re-optimization of these LLS using DFT method (stage 2). The doted profile in stage 1 represent the simplified staircase topology of the PES.

Figure 2. The melting point and pre-melting point of Hf50 with various heating rate. The circle region marks the convergence of melting point lower than certain heating rate (denoted as H.R. in the horizontal axis).


34

Page 35 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a) a)

b) b)

Figure 3. a) The caloric curve and b) 𝑐𝑣 curve of Hf20 obtained via prolonged annealing process. The green arrow indicates the preFigure 4. a) The caloric curve and b) 𝑐𝑣 curve of Hf10 obtained via melting temperature at 𝑇𝑝𝑟𝑒 = 1400 K, and the red arrow direct heating process. The green arrow indicates the preindicates the melting point at 𝑇𝑚 = 1850 K.

melting temperature at 𝑇𝑝𝑟𝑒 = 1350 𝐾, and the red arrow indicates the melting point at 𝑇𝑚 = 2200 𝐾.

a)

b)

Figure 5. a) The caloric curve and b) 𝑐𝑣 curve for the heating of bulk hafnium crystal via direct heating process. The predicted melting point is around 𝑇𝑚𝑏𝑢𝑙𝑘 (𝐻𝑓)~2600 𝐾 to 2650 K.


35


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 49

a)

b)

Figure 6. a) The similarity index 𝜉𝑖 and b) fluctuation of similarity index 𝑆𝜉𝑖 of Hf13 obtain via a direct heating process. The green arrow indicates the pre-melting temperature at 𝑇𝑝𝑟𝑒 = 1600 𝐾, and the red arrow indicates the melting point at 𝑇𝑚 = 2050 𝐾.


36

Page 37 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a)

b)

c)

Figure 7. Plotting together the estimated pre-melting temperature, 𝑇𝑝𝑟𝑒 and the exact melting point, 𝑇𝑚 of Hf clusters of various size 𝑛 for a) prolonged simulated annealing, b) direct heating process, and c) the global similarity index.


37


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 49

Figure 8. The estimated melting point of the hafnium cluster against the cluster size 𝑛, based on three different approaches.

a)

b)

Figure 9. a) Similarity index 𝜉𝑖 curve and b) fluctuation of the similarity index 𝑆𝜉𝑖 of Hf30. The screenshots show the configuration of the cluster Hf30 during that particular temperature. Hafnium atoms in opposite colour tone represent those which are broke free from the cluster core.


38

Page 39 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a)

b)

Figure 10. a) Similarity index 𝜉𝑖 curve and b) fluctuation of the similarity index 𝑆𝜉𝑖 of Hf50. The screenshots shows the configurations of the cluster Hf50 during that particular temperature. The hafnium atoms which are drifted and breaking free from the cluster core are highlighted in inverted color tone.


39


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 49

a)

b)

Figure 11. a) Similarity index 𝜉𝑖 curve and b) fluctuation of the similarity index 𝑆𝜉𝑖 of Hf99. The screenshots shows the configurations of the cluster Hf99 during that particular temperature. The contrast colour of certain hafnium atom showing they are breaking free from the cluster core.


40

Page 41 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


List of tables:

Table 1: Comparing the clusters obtained via COMB potential after 1 K relaxation and those via DFT geometrical re-optimization with B3LYP basis set. Average bond length Point group Similarity 𝑛 index, 𝜉𝑖 MD DFT MD DFT 2 2.97095 2.572 C∞ C∞ 1.0 3 3.00117 2.68649 D3h D3h 0.999195 4 3.01735 2.78485 Td Td 0.999208 5 3.02978 2.82413 D3h D3h 0.977741 6 2.95319 2.87518 Oh Oh 0.974136 7 2.89142 2.88128 C2v D5h 0.923361 8 2.92052 2.89056 C2v C2v 0.960306 9 2.95442 2.91814 D3h D3h 0.952129 10 2.93217 2.90191 C2 D2 0.903499 11 2.95215 2.90484 D4h C2v 0.929598 12 2.90624 2.81625 D5d D5h 0.97928 13 3.00836 2.95388 Ih Ih 0.973953 14 2.94974 2.94608 C2v C2v 0.959454 15 2.96477 2.91166 Cs D6d 0.928096 16 2.92923 2.95073 Cs C2v 0.926347 17 2.93245 2.93065 C4v Td 0.935931 18 2.97061 2.93476 D2 D2 0.955958 19 2.99142 2.97752 D5h D5h 0.969007 20 2.94963 2.94803 C2 D2h 0.961776


41


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 49

Table 2. The melting and pre-melting temperatures obtained from three different approaches. The first four columns are obtained from caloric curves and 𝑐𝑣 curves. All the values of temperatures are in the unit of Kelvin (K).

𝑛 6 10 14 17 18 20 22 24 26 30 34 38 39 40 42 46 50 60 70 80 90 99

Prolonged annealing 𝑇𝑝𝑟𝑒 1250 1500 1550 1200 1450 1400 1600 1350 1500 1650 1500 1550 1600 1350 1650 1600 1450 1500 1650 1800 1750 1650

simulated Direct heating process 𝑇𝑚 1750 1750 1800 1900 1700 1850 1850 1800 1850 1750 1800 1850 1900 1900 1950 1850 1900 1900 2200 2200 2200 2050

𝑇𝑝𝑟𝑒 1250 1350 1200 1250 1450 1500 1500 1700 1950 1650 1950 1650 1700 1400 1800 1800 1600 1550 1650 1800 1900 2150

𝑇𝑚 2000 2200 1900 2250 2200 2150 2150 2050 2300 2250 2300 2500 2350 2100 2100 2250 2250 1950 2300 2500 2350 2500

Fluctuation of the similarity index, 𝑆𝜉𝑖 curves 𝑇𝑝𝑟𝑒 𝑇𝑚 1250 1950 1850 2250 1700 2200 1700 2100 1850 2100 1750 2150 1750 2150 1800 2100 2050 2250 1750 1850 1850 2100 1850 2000 1600 2000 1800 2150 1800 1950 1800 2200 1800 2150 1850 2000 1950 2200 1950 2250 1850 2100 1900 2300


42

Page 43 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Novelty of the work: A method to quantify the chemical similarity of two structures and applied further to locate the melting point of clusters.

List of figures:

Figure 1. A schematic sketch indicating the strategy to obtain true global minimum by the way of sampling LLS at a coarse level search with BH method without structural optimization (stage 1), and subsequently undergo a refined geometrical re-optimization of these LLS using DFT method (stage 2). The doted profile in stage 1 represent the simplified staircase topology of the PES.

Figure 2. The melting point and pre-melting point of Hf50 with various heating rate. The circle region marks the convergence of melting point lower than certain heating rate (denoted as H.R. in the horizontal axis).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 49

a) a)

b) b)

Figure 3. a) The caloric curve and b) 𝑐𝑣 curve of Hf20 obtained via prolonged annealing process. The green arrow indicates the premelting temperature at 𝑇𝑝𝑟𝑒 = 1400 K , and the red arrow indicates the melting point at 𝑇𝑚 = 1850 K.

a)

Figure 4. a) The caloric curve and b) 𝑐𝑣 curve of Hf10 obtained via direct heating process. The green arrow indicates the premelting temperature at 𝑇𝑝𝑟𝑒 = 1350 𝐾 , and the red arrow indicates the melting point at 𝑇𝑚 = 2200 𝐾.

b)

Figure 5. a) The caloric curve and b) 𝑐𝑣 curve for the heating of bulk hafnium crystal via direct heating process. The predicted melting point is around 𝑇𝑚𝑏𝑢𝑙𝑘 (𝐻𝑓)~2600 𝐾 to 2650 K.


Page 45 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a)

a)

b) b)

c) Figure 6. a) The similarity index 𝜉𝑖 and b) fluctuation of similarity index 𝑆𝜉𝑖 of Hf13 obtain via a direct heating process. The green arrow indicates the pre-melting temperature at 𝑇𝑝𝑟𝑒 = 1600 𝐾, and the red arrow indicates the melting point at 𝑇𝑚 = 2050 𝐾.

Figure 7. Plotting together the estimated pre-melting temperature, 𝑇𝑝𝑟𝑒 and the exact melting point, 𝑇𝑚 of Hf clusters of various size 𝑛 for a) prolonged simulated annealing, b) direct heating process, and c) the global similarity index.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 49

a)

Figure 8. The estimated melting point of the hafnium cluster against the cluster size 𝑛, based on three different approaches.

b)

Figure 9. a) Similarity index 𝜉𝑖 curve and b) fluctuation of the similarity index 𝑆𝜉𝑖 of Hf30. The screenshots show the configuration of the cluster Hf30 during that particular temperature. Hafnium atoms in opposite colour tone represent those which are broke free from the cluster core.


Page 47 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a)

a)

b)

b)

Figure 10. a) Similarity index 𝜉𝑖 curve and b) fluctuation of the similarity index 𝑆𝜉𝑖 of Hf50. The screenshots shows the configurations of the cluster Hf50 during that particular temperature. The hafnium atoms which are drifted and breaking free from the cluster core are highlighted in inverted color tone.

Figure 11. a) Similarity index 𝜉𝑖 curve and b) fluctuation of the similarity index 𝑆𝜉𝑖 of Hf99. The screenshots shows the configurations of the cluster Hf99 during that particular temperature. The contrast colour of certain hafnium atom showing they are breaking free from the cluster core.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 48 of 49

List of tables: Table 2: The melting and pre-melting temperatures obtained from three different approaches. The first four columns are obtained from caloric curves and 𝑐𝑣 curves. All the values of temperatures are in the unit of Kelvin (K). 𝑛

Prolonged simulated annealing 𝑇𝑝𝑟𝑒

6 10 14 17 18 20 22 24 26 30 34 38 39 40 42 46 50 99

1250 1500 1550 1200 1450 1400 1600 1350 1500 1650 1500 1550 1600 1350 1650 1600 1450 1650

𝑇𝑚 1750 1750 1800 1900 1700 1850 1850 1800 1850 1750 1800 1850 1900 1900 1950 1850 1900 2050

Direct heating process 𝑇𝑝𝑟𝑒 1250 1350 1200 1250 1450 1500 1500 1700 1950 1650 1950 1650 1700 1400 1800 1800 1600 2150

Table 1: Comparing the clusters obtained via COMB potential after 1 K relaxation and those via DFT geometrical re-optimization with B3LYP basis set. Average bond length Point group Similarity 𝑛 index, 𝜉𝑖 MD DFT MD DFT 2 2.97095 2.572 C∞ C∞ 1.0 3 3.00117 2.68649 D3h D3h 0.999195 4 3.01735 2.78485 Td Td 0.999208 5 3.02978 2.82413 D3h D3h 0.977741 6 2.95319 2.87518 Oh Oh 0.974136 7 2.89142 2.88128 C2v D5h 0.923361 8 2.92052 2.89056 C2v C2v 0.960306 9 2.95442 2.91814 D3h D3h 0.952129 10 2.93217 2.90191 C2 D2 0.903499 11 2.95215 2.90484 D4h C2v 0.929598 12 2.90624 2.81625 D5d D5h 0.97928 13 3.00836 2.95388 Ih Ih 0.973953 14 2.94974 2.94608 C2v C2v 0.959454 15 2.96477 2.91166 Cs D6d 0.928096 16 2.92923 2.95073 Cs C2v 0.926347 17 2.93245 2.93065 C4v Td 0.935931 18 2.97061 2.93476 D2 D2 0.955958 19 2.99142 2.97752 D5h D5h 0.969007 20 2.94963 2.94803 C2 D2h 0.961776

𝑇𝑚 2000 2200 1900 2250 2200 2150 2150 2050 2300 2250 2300 2500 2350 2100 2100 2250 2250 2500


Fluctuation of the similarity index, 𝑆𝜉𝑖 curves 𝑇𝑝𝑟𝑒 1250 1850 1700 1700 1850 1750 1750 1800 2050 1750 1850 1850 1600 1800 1800 1800 1800 1900

𝑇𝑚 1950 2250 2200 2100 2100 2150 2150 2100 2250 1850 2100 2000 2000 2150 1950 2200 2150 2300

Page 49 of 49


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Investigation of Melting Dynamics of Hafnium Clusters - American

Recommend Documents