EFFECTIVE INTERATOMIC POTENTIALS BASED ON THE FIRST-PRINCIPLES MATERIAL DATABASE

Effective interatomic potentials are frequently utilized for large-scale simulations of materials. In this work, we generate an effective interatomic potential, with Niobium as an example, using the force-matching method derived from a material database which is created by the first-principle molecular dynamics. It is found that the potentials constructed in the present work are more transferable than other existing potential models. We further discuss how the first-principles material database should be organized for generation of additional potential.


INTRODUCTION
Molecular dynamics (MD) simulation is a useful tool for researching material properties, and various materials from metals to insulators have been studied by such simulations.The framework of effective potentials for metals was constructed on the basis of an embedded atom method (EAM) (Daw, 1984).Analytical EAM potentials have been generated for FCC metals.On the other hand, specialized potential types for semiconductors have been developed in order to reproduce their covalent bonds (Stillinger, 1985;Tersoff, 1986).In wake of this, the modified version of the EAM (MEAM) was derived to describe the covalent bonding of silicon and germanium (Baskes, 1987(Baskes, , 1989)).Moreover, the MEAM framework was applied to the BCC and HCP metals to reproduce realistic material properties (Baskes, 1992).For BCC metals, the Finnis-Sinclair model based on the tight-binding theory was much better than the EAM model (Finnis, 1984).Although effective interatomic potentials have been frequently utilized in limited situations, they need to be enhanced to simulate material properties under a much wider variety of conditions.
First-principles simulation based on the density functional theory (DFT) has been successful for calculating chemical and physical properties of solid materials, molecular systems, and liquids.First-principles MD simulations have been applied to discover dynamic properties of materials without any empirical parameters.In order to estimate equilibrium properties at finite temperatures from a MD simulation, a long simulation time with many atoms is required.Because of the lack of computer power, the simulation time of the first-principles MD, however, is presently restricted to the order of at least a pico-second, and the number of atoms in the simulation cell is limited to several hundred.To overcome this limitation, Ercolessi and Adams proposed a method to create effective interatomic potential applying information from first-principles MD by a force matching method (Ercolessi, 1994).The force matching method has been applied to metallic (Ercolessi, 1994;Li, 2003;Grochola 2005), covalent (Lenosky, 2000), and ionic systems (Umeno, 2002;Aguado, 2003).
In this paper, we report MD simulation results obtained using interatomic potentials constructed by the force-Data Science Journal, Volume 8, 20 May 2009 matching method on the basis of the first-principles forces and physical properties.We also discuss a first-principles material database for potential generation.Calculations on Niobium are presented as an example system.

The force matching method
Newtonian motion of an atomic system evolves by forces acting on atoms.Because the forces directly affect the dynamics of the atomic system, accurate forces are needed to simulate dynamics in a material.The Hellmann-Feynman forces calculated by the DFT method yield better results than those forces calculated by empirical methods.Because DFT calculations have a higher computational cost than the empirical ones, effective interatomic potentials reproducing the DFT forces are desired, which can be obtained by the force matching method.
We briefly describe the force-matching method (Ercolessi, 1994) for potentials of the EAM type as a typical example whose total energy is given by where φ is the two-body interatomic potential, ij r is the interatomic distance between atoms i and j, U is the embedded function, and f is the electron density (Daw, 1984).The effective interatomic forces are given by partially differentiating the total energy with respect to atomic coordinates.A first-principles MD simulation gives the forces acting on atoms in various atomic configurations.The mean square error between the EAM and firstprinciples force is A .P is the number of the properties, and l W is the weight factor for the square error of the property l.
Moreover, we add the penalty function to reduce corrugations of the spline functions.Finally, the potential fitting is done by minimizing the total objective function

Computational details
pseudopotential technique (Payne, 1992) to obtain physical properties of materials.Niobium is taken as an example.
We employ the generalized gradient approximation (GGA) for the exchange-correlation functional (Perdew, 1996) and an ultra-soft pseudopotential for the Niobium ion (Vanderbilt, 1990).The ground state of Niobium is a BCC structure, and hypothetical FCC structures are taken into account.For BCC (FCC) bulk Niobium, we use an energy cutoff of 30 Ry and a 10x10x10 (8x8x8) Monkhorst-Pack grid for k-space integrations.For other structures, we use the same energy cutoff and k-point samplings, which give almost equivalent accuracies to the bulk calculation.A database of chemical and physical properties of Niobium is created by first-principles calculations, and a database of Hellmann-Feynman forces acting on atoms at various configurations of solid and liquid states is generated by firstprinciples MD simulations.In each MD simulation, the volume is kept constant, and the Nosé-Hoover thermostat (Nosé, 1984;Hoover, 1985) keeps the temperature constant.
In the force-matching procedure, we use the forces obtained by the first-principles MD and material properties obtained by experiments and calculations.We perform the MD simulations for several Niobium structures, solid and liquid as listed in Table 1.The calculated material properties are lattice constants, bulk moduli, elastic constants, vacancy formation energy, the FCC-BCC energy difference, and pressures.The experimental and theoretical values are listed in Table 2, and the theoretical pressures of BCC solid Niobium calculated at the scaled lattice constants are shown in Table 3.We use the experimental value for the lattice constant of BCC solid and the theoretical value of FCC solid that is scaled with the ratio of experimental and theoretical values of BCC structure.The atomic configurations of the MD trajectories are also scaled in length.We employ a symplectic integrator for the classical MD simulation with the Nosé-Poincaré thermostat and the Andersen barostat (Aoki, 2008).The cubic simulation cell has 250 Niobium atoms.The external pressure is set to zero bar.The temperatures are set to a series of values up to 3300 K.

Molecular dynamics simulations
We performed classical MD simulations with a symplictic integrator (Aoki, 2008) for Niobium using the forcematching EAM potentials.Figure 2 shows that the calculated thermal lattice expansion from the EAM FM-B model agrees very well with experimental results over the whole range of temperatures.The melting point is about 3000 K, which is in rough agreement with the experimental value of 2750 K.It should be noted that the EAM FM-A model cannot reproduce the experimental result because of the lack of the information of the solid state.The lattice constants of both models extracted at 0 K agree well with the experimental value of 0.3295 nm.The lattice constant of the EAM FM-A model has a minimum at 200 K.It is considered that the low electron density at short interatomic distances results in a low many-body repulsive force.Because of the lower electron density as compared to the FM-B model as discussed in the previous section, the lattice contraction below 200 K for the FM-A model would be caused by the lack of repulsive force.Notably, the thermal expansion coefficients above 500 K by the EAM FM-A and FM-B models are better than by the extended Finnis-Sinclair potential (Dai, 2006).

The first-principles material database
In the force-matching procedure, highly reliable materials databases containing structural properties data and force data are very important.In our calculation, we used the experimental lattice constant of the BCC ground state structure and the scaled theoretical lattice constant of the hypothesis FCC structure, without using either the theoretical or experimental cohesive energies.Because normally, calculated lattice constants and cohesive energies are not exactly equivalent to the experimental ones, theoretical values from different methods, using values from more advanced theoretical scheme is recommended.Our database can be further improved with theoretical lattice constants and cohesive energies that are calculated by first-principles methods that treat electron correlations more accurately, such as quantum Monte Carlo methods.
The force database was constructed from the first-principles MD simulations.However, the DFT forces within the GGA approximation do not exactly reproduce the experimental dynamical properties.Therefore, we scaled the atomic configurations into lengths that fit the experimental lattice constant.Scaling factors in length and energy can be determined by a trial and error manner that fits the radial distribution function using a trial potential based on the scaled-force database to the experimental value (Grochola, 2005).It would be preferable to have the scaling factors be determined by beyond-DFT methods because it is not easy to obtain good experimental data of pure materials in various cases.
In metals with localized electrons such as Cerium, we also need force and property data that take into account the valence fluctuation effects.These effects can be treated by both extended DFT methods based on the hybrid exchange-correlation functional combined with the density functional and an orbital-dependent functional (Hay,  2006; Da Silva; 2007) and with the effective coulomb interaction correction, namely the DFT+U method (Fabris,  2005; Loschen 2008).
In order to apply the force-matching method to various materials, the force database calculated by the first-principles method in high accuracy, by both DFT and beyond-DFT, is very important because it is almost impossible to obtain experimental force data.A general concept in our present scheme is shown in Figure 3.

CONCLUSION
We constructed an effective interatomic potential for Niobium by the force-matching method based on a firstprinciples material database.We then simulated the thermal expansion of Niobium by classical MD using the forcematching EAM potential.It was shown that our force-matching potential produced more reasonable simulation results than what is obtained by the extended Finnis-Sinclair potential.We further discussed the preferable construction of the database.The first principles force database in the future will become one of the most important fundamental data sources to improve calculation of effective interatomic potentials by the force-matching method.

FF
is the first-principles force of the atom i at the configuration k, and ( ) α ik is the EAM force at the atomic position in the first-principles MD trajectory.M is the number of configurations, and N k is the number of atoms in the configuration k. α is the set of parameters that define the EAM-potential functionsφ , U , and f .We do not use the analytic functions frequently used in empirical fittings but alternatively use spline functions.The flexibility of spline functions makes it feasible to realize transferability of the potential.The force fitting is done by minimizing the objective function chemical and physical properties given by the EAM potential Results for two different EAM potentials of Niobium obtained by the force-matching method are shown in Figures1 (a), (b), and (c).One, denoted as FM-A, was generated based on the MD trajectories of liquid only, while another one, denoted as FM-B, was created by the full trajectories of our database.The minimum of the pair-wise potential is located at 0.33 nm.The embedded function in Figure1(c) has a minimum at unity of the electron density.The electron density functions, shown in Figure1(b), deviate at small interatomic distance.The FM-B potential results in a larger electron density than the FM-A potential at interatomic distances lower than 0.28 nm because the firstprinciples MD simulation of a solid state generates force data at a smaller interatomic distance than that of liquid state.The difference comes from whether the potential used the force data of the solid states in addition to the force data of the liquid state.

Figure 1 .
Figure 1.(a)The pair-wise potential functions, (b) the electron density functions, and (c) the embedded functions of the force-matching EAM potentials (EAM FM-A and -B) for Niobium

Figure 2 .
Figure 2. Thermal expansions of Niobium obtained by the EAM FM-A and -B potentials (long-and short-dashed lines, respectively) and the extended Finnis-Sinclair potential (dotted line).The experimental curve is also indicated in the solid line.

Figure 3 .
Figure 3. Building up the first-principles (FP) database for materials design

Table 1 .
Niobium structures for the first-principles MD simulation, and the number of configurations used in the force-matching procedure for each structure

Table 3 .
Theoretical pressures of BCC solid Niobium.Lattice constants are given in nm, and pressures are in GPa.