A New Many-Body Expansion Scheme for Atomic Clusters: Application to Nitrogen Clusters†

a. State Key Laboratory of Magnetic Resonance and Atomic and Molecular Physics, Wuhan Institute of Physics and Mathematics, Innovation Academy for Precision Measurement Science and Technology, Chinese Academy of Sciences, Wuhan 430071, China b. University of Chinese Academy of Sciences, Beijing 100049, China c. Institute of Chemical Materials, China Academy of Engineering Physics, Mianyang 621900, China d. Wuhan National Laboratory for Optoelectronics, Huazhong University of Science and Technology, Wuhan 430071, China

fragments, the energy of the system can be expressed in terms of the many-body expansion (MBE) method [8][9][10]. In practice, MBE is always truncated at a specific n-body level so that the energy of a system is approximated by numerous calculations of small fragments. The reliability of the MBE depends on whether high-order energy contribution can be neglected. The n-body interaction in the MBE method could be obtained from a pre-constructed potential energy surface or from quantum chemistry calculations.
MBE has been successfully applied to represent the interaction energy of a variety of systems, with which the physical and chemical properties of these systems were well investigated. Mayhall and Raghavachari proposed a many-overlapping-body (MOB) expansion and demonstrated that the MOB methodology significantly improves the two-body energies for four dendritic isomers of C 29 H 60 and 10 conformational isomers of a polypeptide molecule [11]. Paesani and his coworkers developed a "First Principles" water potential with flexible monomers by building upon the MBE of the interaction energy [12][13][14]. The calculated vibrationrotation tunneling spectrum, second virial coefficient, third virial coefficient, thermodynamic, and dynamical properties agree well with the available experimental data. Bowman and his coworkers reported the fulldimensional potential energy surfaces for water clusters and protonated water clusters, on which the predicted infrared spectra are in consistence with experimental measurements [15][16][17]. The cohesive energy of solid argon was computed by Schwerdtfeger et al. with an accuracy of about 5 J/mol by using wave-function-based methods within a many-body decomposition of the interaction energy [18,19]. The method of increments, a wavefunction-based ab initio correlation method, was proven to work well for solids [20]. Smits et al. studied the solid-to-liquid phase transition for krypton and xenon based on many-body relativistic coupled-cluster interaction potentials [21]. The calculated melting temperatures were in very good agreement with experimental values. Yao et al. combined the MBE with neural network, in which the n-body energies (up to threebody) were represented by trained neural networks [22]. Varandas and his coworkers constructed global molecular potentials for carbon clusters relying on the double many-body expansion approach with their PESs of ground-state singlet C 3 and triplet C 4 as input [23,24]. The convergence of the MBE of interaction potentials between atoms and molecules was analyzed by Schwerdtfeger and his coworkers in detail for different types of interactions [25]. They concluded that the MBE converges rapidly for long-range weak interactions, such as the interaction energy and the electron correlation energy of van der Waals-and hydrogen-bonded clusters, however, the n-body decomposition scheme converges very slowly in the case of metallic or covalent interactions. Therefore, the development of accurate interaction potentials for the metallic or covalent clusters is very challenging.
In recent years, machine learning has become a promising tool to construct interatomic potentials [26]. Impressively, Behler and Parrinello proposed a generalized neural-network method for constructing high-dimensional potential energy surfaces [27,28]. The potential energy E was given by the sum of local atomic energies E a of all atoms in the system. The atomic energy E i depends on the local chemical environment of atom a that is characterized by sets of atom-centered symmetry function descriptors. In this work, we propose a new MBE scheme to represent interatomic potentials. Like the high-dimensional neural network method, the potential energy is expressed as a sum of atomic energy. The atomic energy is, however, represented by low-order many-body expansion, instead of using individual neural network. The atom-based many-body expansion method is first applied to construct the interatomic potentials of nitrogen clusters.
Nitrogen is the most abundant element in atmosphere. Polynitrogen compound is one of the ideal candidates for high energy density materials (HEDM) [29][30][31][32][33][34][35][36]. Since a large amount of energy is released when polynitrogen compounds dissociate into N 2 molecules, they can be used as propellants and explosives. Some global PESs for N 3 and N 4 systems have been reported. Wang et al. proposed the first full-dimensional 4 A ′′ state PES for N 3 , in which the ab initio energies of sampled points were computed at the level of CCSD(T)/cc-pVQZ [37]. Galvão and Varandas used the double many-body expansion method to construct the ground state PES of N 3 , in which 1,592 sampled points are calculated by CCSD(T) and MRCI+Q extrapolated to complete basis set (CBS) limit [38]. Truhlar's group reported a series of global ground state PES for N 4 [39][40][41]. Impressively, Li et al. constructed the many-body component of the ground-state potential energy surface of N 4 by fitting a total of 21406 energy points by using neural networks with permutationally invariant polynomials and the mean unsigned deviation was reduced by a factor of three compared with least-squares methods [42]. Recently, the performance of the MBE approach for calculating potential energy of nitrogen clusters was assessed by our group [43]. We found that the MBE approach could give reasonable energy predictions for loose chain nitrogen clusters while it produced bad predictions for compact structures. The poor performance arose from the neglect of high-order (n≥5) n-body interactions that are significant in the bonding area. Considering the bad performance of traditional MBE method in calculating the energy of systems containing covalent bonds, we aim to develop a modified MBE method in this work to improve the prediction accuracy of energy.
When applying the atom-based many-body expansion method to calculate the potential energies of nitrogen clusters, the two-body, three-body and four-body interaction potentials are obtained by carrying out ab initio calculations. This paper is organized as follows. Section II briefly introduces the atom-based many-body expansion method, derived from the conventional MBE formalism. The results and discussions are presented in Section III. Finally, conclusions are given in Section IV.

A. Interaction many-body expansion (IMBE)
The many-body expansion of the energy of a Natomic system is written as, where E (n) denotes individual n-body contribution. The explicit forms of E (n) are [44], . . .
The one-body energy term E (1) is the sum of the energies of all separated atoms. The two-body energy term E (2) is the sum of two-body interaction energy V 2 (i, j) that is defined as the difference between the two-body energy E 2 (i, j) and the one-body energy of the atom i and j. Similarly, the three-body interaction V 3 (i, j, k) is defined as the difference between the three-body energy E 3 (i, j, k) and the sum of onebody energy Similarly, one can obtain n-body energy term with n up to N .
As has been pointed by Hermann et al. [25], the convergence of MBE is very slow for metallic or covalent clusters and the many-body potential is difficult to be constructed. Here, we reformulate the MBE and test the convergence of the new method. The method is based on the many-body expansion of the interaction of monomers and is named as interaction many-body expansion (IMBE). For the system with two atoms, the MBE energy is, Because the interaction V 2 is contributed from two atoms, it can be divided into two parts. Combined with the one-body energy, the energy of the system is expressed as, For a system including three atoms, the energy is given by For a system containing four atoms, the energy is expressed as From the definition of V k (1, 2, · · · , k) (k≥2), the value of V k is close to zero if the distance between any two atoms is large enough. Thus, the number of terms in summation of interaction could be significantly reduced if neglecting these vanishing terms. For a system including N atoms, the general energy expression is Such a reformulation of the MBE is similar to the calculation of atomic charge in the Mullikan population analysis. The energy of a system is the sum of individual atomic energy E(k). E(k) includes both the energy of noninteraction atom and the interaction of the atom with its surroundings. Similar to the MBE, the interaction V k is obtained from ab initio calculations.

B. Non-overlapping and overlapping schemes
In practical application of the MBE method, the expansion is always truncated at a certain order. The same principle is also suitable for the IMBE approach. Here a tetra-atomic cluster is used to demonstrate the truncation scheme in the IMBE approach. For example, the energy E(1) is exactly expressed as, Supposing that the expansion is truncated up to threebody interaction, two fragmentation schemes could be applied.

Non-overlapping scheme
The scheme divides the surround atoms into several non-overlapping fragments. As shown in the left panel of FIG. 1, the interaction of the 1st atom with its surrounding atoms is approximated as the sum of interactions between the 1st atom with the fragment 23 and the 4th atom. This fragmentation scheme is denoted as 1(23/4) and the corresponding energy is expressed as The omitted terms are, As aforementioned, these interactions will be negligible if the distance R 24 or R 34 is large enough.

Overlapping scheme
When the interaction among the atoms 1, 3, and 4 cannot be neglected, the overlapping scheme is applied (the right panel of FIG. 1). This scheme is denoted as 1(23/34). It is noted that the atom 3 is included in both 123 and 134 fragments and the interaction between the atom 1 and 3 is counted two times. Thus, one 1−3 interaction should be removed from the energy expression, The omitted terms are Compared with the traditional MBE method, further truncations are implemented in the IMBE method, which depends on what fragmentation scheme are employed. If the MBE and IMBE are truncated at a given order and all terms up to the order are included in the IMEB method, both give the same result.

III. RESULTS AND DISCUSSION
The IMBE approach is first employed to evaluate the energy of nitrogen clusters. In this work, 19 nitrogen atomics clusters are studied by both the MBE and IMBE to assess the performance of the IMBE approach. These clusters are obtained from Ref. [45] and their structure are presented in FIG. 2. The energies of these clusters are calculated at the level of B3LYP/aug-cc-pVTZ. All ab initio calculations are performed by Gaussian 09 software [46]. In order to analyze the convergence of n-body expansion, the energies of fragments for each cluster are calculated at the same level as well. Because the energy of cluster or fragment is related to its spin multiplicity, the calculations with different spin multiplicities are carried out and the lowest energy is used in the following MBE/IMBE analysis. In the strict sense, the spin of each fragment should follow the Wigner-Witmer spin-spatial correlation rules [47,48], which make the Fragment-based approach give the correct dissociation products. This is an important prerequisite for accurately simulating molecular dynamics. However, as to the IMBE method, the central atom is always shared by different fragments in both non-overlapping and overlapping schemes. Therefore, it is infeasible to correlate the spin of each fragment with the term manifold of the constituent cluster.

A. Fragmentation rules
As introduced in Section II, the IMBE method decomposes the energy of cluster as the sum of the atomic energy E(k). The atomic energy E(k) includes the energy of isolated atom and the interaction energy of the atom with its surrounding atoms. In this work, the many-body expansion is truncated to the fourth order. The surrounding atoms are divided into a number of fragments so that the number of atoms in each fragment is not larger than four. Obviously, there exist many ways to divide the surrounding atoms into fragments and their effects on the evaluation of the cluster energy are unclear. This subsection will explore the basic fragmentation rules within the IMBE approach by inspecting four nitrogen clusters. The four nitrogen clusters are N 4 (T d 1 A 1 ), N 5 (D 2d 2 B 1 ), N 6 (C 2h 1 A g ), and N 6 (D 3h 1 A ′ 1 ). The energies from different fragmentation schemes are calculated for each cluster, for which the schemes with smaller errors are listed in Table I.  Table I, the fragmentation schemes with small errors have some common characters: (i) The distances between each twoatom in each fragment should be as small as possible. Small interatomic distance means large interaction, which should be taken into account in the fragmentation scheme; (ii) If some surrounding atoms are close enough and the number of these atoms excesses three, the overlapping fragmentation scheme would give better results.

B. Performance
This subsection compares the MBE and IMBE energy for twelve nitrogen clusters. Because the high-order ex-  Note: the number before the first slash denotes the central atom and the number(s) in between the slash or after the slash represent the atoms in a fragment.
pansions in both the MBE and IMBE approaches are expensive, the expansion is limited up to the fourth order. The ab initio calculations are carried out to obtain the energies of the twelve clusters and the two-, three-, four-body fragment energies that are used to evaluate the MBE and IMBE energies. For the IMBE energy, the fragmentation schemes are determined by utilizing the rules obtained in the above subsection, which are listed in Table II.
The MBE and IMBE errors with respect to the DFT results are given in Table III. The binding energy for the chosen nitrogen clusters ranges from −779.1 kcal/mol to −369.8 kcal/mol. As expected, the performance of the MBE method is very poor. The MBE energy deviates seriously from the DFT energy, with the largest difference being 7740.8 kcal/mol. In contrast, the IMBE calculation produces satisfactory results. The difference is mostly less than 10 kcal/mol. The largest deviation occurs for the cluster N 10 (D 5h 1 A ′ 1 ), being 23.2 kcal/mol. Compared with the other clusters, the cluster N 10 (D 5h 1 A ′ 1 ) has a very compact structure, indicating that the higher-order interaction is possibly nonnegligible. Overall, the IMBE method provides a noticeable improvement over the traditional MBE method. Compared with the MBE method, the threebody and (or) four-body terms are further truncated in the IMBE method. These truncated terms appear to be able to compensate the errors introduced by the neglection of higher-order (n>4) interactions.
For the clusters involving covalent bonds, the MBE error has been shown to depend strongly on the size of clusters. Thus, it is infeasible to apply the traditional MBE method to large clusters with strong interactions. FIG. 3 shows the size and structure dependence of the MBE and IMBE errors with respect to the DFT calculations. For clarity, the error in this figure is defined as the ratio of the absolute difference between the IMBE (or MBE) energy and the DFT energy to the DFT energy. The errors of isomers for a specified number of nitrogen atoms are averaged in FIG. 3(a). Clearly, the     8 23.2 a The binding energy is the difference between the DFT energy and the sum of the energy of non-interaction nitrogen atom. b The MBE error is the difference between the DFT energy and the MBE energy. c The IMBE error is the difference between the DFT energy and the IMBE energy.
MBE error is positively correlated with the size of the cluster. As the number of nitrogen atoms increases, the error rises from 24.7% to 948.4%. In contrast, the IMBE error does not have a clear size dependence, which oscillates around 1%. This indicates the IMBE method is suitable to be applied to larger systems.
To explore the structure dependence of error, the structure of each cluster is described by the ratio of the maximum rotational inertia to the minimum one in FIG. 3(b)  error appears to decline with the increase of I max /I min , in sharp contrast to the MBE error.
As aforementioned, the lowest energy of a cluster or a fragment with different spin multiplicities is used in the IMBE analysis. This treatment, however, cannot guarantee the IMBE approach to give correct dissociation products. To test the practical performance of IMBE approach, FIG. 4 compares the IMBE energy with the ab initio energy for the N 6 (C 2h 1 A g ), N 7 (C 2v 2 B 1 a) and N 8 (C 2h 1 A g ) clusters along possible dissociating paths.
The same fragmentation rules are applied in the IMBE calculations. Clearly, the IMBE curve varies smoothly from the strong interaction region to the dissociation region. The change of the IMBE energy agrees reason-ably well with ab initio calculations. Although there exist some visible discrepancies between the IMBE and ab initio energies in the bonding region, the treatment could give correct energies for dissociating fragments.

IV. CONCLUSION
The traditional many-body expansion (MBE) approach is infeasible to be applied to atomic clusters including strong interactions. In this work, we develop a new approach, namely the interaction many-body expansion (IMBE) approach, to calculate the binding energy of covalent or metal atomic clusters. With the same spirit as the high-dimensional neural network method, the potential energy is defined as the sum of atomic energies. The atomic energy is represented by low-order many-body expansion. To assess the performance of the IMBE approach, it is utilized to calculate the energies of nitrogen clusters. The IMBE method has made particularly noticeable improvements over the MBE method. The energy differences between DFT and IMBE calculations for the nitrogen clusters are around 10 kcal/mol, much smaller than those from MBE calculations. In addition, the IMBE error does not have a clear size and structure dependence, implying that the IMBE method can be applied to evaluate larger clusters. Although the performance of the IMBE method in the bond-breaking and bond-forming area has not been fully assessed, the conclusion is thought to be reasonable at least for stationary structures. We plan to apply the IMBE method to metal clusters in the near future.