{'fraction_non_alphanumeric': 0.029332289015077343, 'fraction_numerical': 0.012845114548658703, 'mean_word_length': 4.946902654867257, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Development of a Transferable Reactive Force Field of P/H Systems: Application to the Chemical and Mechanical Properties of Phosphorene

ReaxFF provides a method to model reactive chemical systems in large-scale molecular dynamics simulations. Here, we developed ReaxFF parameters for phosphorus and hydrogen to give a good description of the chemical and mechanical properties of pristine and defected black phosphorene. ReaxFF for P/H is transferable to a wide range of phosphorus and hydrogen containing systems including bulk black phosphorus, blue phosphorene, edge-hydrogenated phosphorene, phosphorus clusters and phosphorus hydride molecules. The potential parameters were obtained by conducting unbiased global optimization with respect to a set of reference data generated by extensive ab initio calculations. We extend ReaxFF by adding a 60{\deg} correction term which significantly improves the description of phosphorus clusters. Emphasis has been put on obtaining a good description of mechanical response of black phosphorene with different types of defects. Compared to nonreactive SW potential [1], ReaxFF for P/H systems provides a huge improvement in describing the mechanical properties the pristine and defected black phosphorene and the thermal stability of phosphorene nanotubes. A counterintuitive phenomenon is observed that single vacancies weaken the black phosphorene more than double vacancies with higher formation energy. Our results also show that mechanical response of black phosphorene is more sensitive to defects for the zigzag direction than for the armchair direction. Since ReaxFF allows straightforward extensions to the heterogeneous systems, such as oxides, nitrides, ReaxFF parameters for P/H systems build a solid foundation for the reactive force field description of heterogeneous P systems, including P-containing 2D van der Waals heterostructures, oxides, etc.


Introduction
In recent years, two-dimensional (2D) materials have attracted much interest because of their fascinating electronic, 2,3 mechanical, 4 optoelectronic, 5,6 and chemical 7,8 properties. The epic discovery of graphene opened up the possibility of isolating and studying the intriguing properties of a whole family of 2D materials including the 2D insulator boron nitride (BN), 9-11 2D semiconductor molybdenum disulfide 9,12,13 and recently, 2D phosphorus, i.e. phosphorene. 14,15 Single layer black phosphorus, i.e. phosphorene, was obtained in experiments in 2014. 15 Because of its tunable band gap and high carrier mobility, phosphorene holds great potential in electronic and optoelectronic applications.
Over the past decade, tremendous success has been achieved in the synthesis of 2D materials.
However, the cycles of synthesis, characterization and test for 2D materials are slow and costly, which inspired the development of computational tools to design new 2D materials [16][17][18] and to provide guidance for the fabrication of 2D devices. [19][20][21][22] Although ab initio methods (such as density functional theory, DFT) provide accurate description of the electronic structure of 2D crystals, they are limited to small systems (several hundreds of atoms) with short time scales (picoseconds). To the contrary, molecular dynamics simulations powered by force fields are able to reach much larger scale with much longer time. To date, several force fields have been developed for black phosphorus.
A valence force field (VFF) for black phosphorus was first proposed in 1982 23 and used to study the elastic properties in black phosphorus. More recently, Jiang et al. 1

developed a
Stillinger-Weber (SW) potential for phosphorene based on the VFF model by fitting parameters to experimental phonon spectrum. In the SW potential, the energy parameters were taken from the VFF model, and geometrical parameters were derived analytically from the equilibrium state of individual potential terms. While both VFF model and SW potential have been used to describe phonons and elastic deformations, they are not suitable to describe states far from equilibrium. 24 Moreover, the SW potential strongly underestimated the Young's modulus of black phosphorene in the zigzag direction. 25 Due to its nonreactive nature, SW potential also had difficulty describing phosphorene with defects. An improved force field which balances accuracy and computational efficiency is therefore highly desirable. In 2001, van Duin et al. developed a reactive force field (ReaxFF) for hydrocarbons. 26 ReaxFF is a bond order interaction model, capable of handling bond breaking and forming with associated changes in atomic hybridization.
In this work, we develop a ReaxFF parameter set for P and H to describe the chemical and mechanical properties of pristine and defected black phosphorene. ReaxFF for P/H is transferable to a wide range of phosphorus and hydrogen containing systems including bulk black phosphorus, blue phosphorene, hydrogenated phosphorene, phosphorus clusters and phosphorus hydride molecules. The ReaxFF parameters for P/H were fitted to a set of reference data generated by extensive ab initio calculations. The proposed ReaxFF for P/H provides a distinctive improvement in describing the thermomechanical properties the pristine and defected black phosphorene, as well as that of the phosphorene nanotubes (PNTs) over the SW potential.
The ReaxFF parameters for P/H presented here provide a first step in the development of a reactive force field description for the heterogeneous P systems.

DFT Calculations
The fitting data used for P/H systems was obtained from DFT calculations performed with the Cambridge series of total-energy package (CASTEP). 34,35 For these calculations, ultrasoft pseudopotentials were used to describe the core electrons and the electron exchange-correlation effects was described by the Perdew-Burke-Ernzerhof (PBE) 36 generalized gradient approximation. In this work, the empirical dispersion correction schemes proposed by Grimme (D2) 37 was used in combination with PBE functional. In computing the energies of phosphorus clusters, phosphorus hydride molecules and phosphorene with defects and adatoms, spin polarization was used to account for the energy contributions from magnetization. Periodic boundary conditions were used for all the calculations, with monolayer structures represented by a periodic array of slabs separated by a 15 Å thick vacuum region. A large 5 × 7 supercell of black phosphorene was adopted to study the effect of defects and adatoms. A plane wave cutoff of 520 eV was used to determine the self-consistent charge density. For condensed phases, Brillouin zone integrations were performed with Monkhorst-Pack 38 mesh with 0.02 Å -1 k-point spacing. For cluster calculations, a cubic supercell of 20 Å (to ensure the interactions between clusters in adjacent cells is negligible) was used with the clusters or molecules placed at the center of the cell with Brillouin zone sampled at the Γ point. All geometries were optimized by CASTEP using the conjugate gradient method (CG) with convergence tolerances of a total energy within 5.0 × 10 −6 eV atom −1 , maximum Hellmann-Feynman force within 0.01 eV Å −1 , maximum ionic displacement within 5.0 × 10 −5 Å, and maximum stress within 0.01 GPa. For black phosphorene, the stress-strain responses in the armchair and zigzag directions were calculated using the method described in the references 39

ReaxFF
The ReaxFF model [26][27][28] is a bond-order interaction model. For ReaxFF, the interatomic potential describes chemical reactions through a bond-order framework, in which the bond order is directly calculated from interatomic distances. Within the bond order framework, the electronic interactions (i.e. the driving force of the chemical bonding) are treated implicitly, allowing the method to simulate chemical reactions without expensive quantum chemical calculations. Typical empirical force field (EFF) potentials adopt empirical equations to describe the bond stretching, bond bending, and bond torsion events, with additional expressions to handle the van der Waals (vdW) and Coulomb interactions. These EFF potentials require a userspecified connectivity table, while ReaxFF is able to calculate the atom connectivity on the fly, which distinguishes ReaxFF from conventional EFF potentials since the breaking and forming of bonds can be captured during MD simulations.
For a ReaxFF description of P/H systems, the bond energies ( ) are corrected with overcoordination penalty energies ( ). Energy contributions from valence angle ( ) and torsion angle ( ) are included. Dispersion interactions are described by the combination of the original van der Waals term ( ) and low-gradient vdW correction term ( ) 31 . Energy contribution from Coulomb interactions ( ) are taken into account between all atom pairs, where the atomic charges are calculated based on connectivity and geometry using the Electron Equilibration Method (EEM). 45 All energy terms except the last three are bond-order dependent and a detailed description of them (except 60 ) can be found in Refs. 26,27,31 . The total energy is the summation of these energy pieces, shown by The stability of P4 cluster and the instability of larger phosphorus clusters has been an ongoing puzzle for several decades. 46 The phosphorus is often expected to favor valence angles near 101°. 47 If this is true, the strain energy of bonds in P4 cluster (with 60° valence angles) should make it unstable. In QM calculations, this problem was resolved by including the effect of d-orbitals. 48 In order to address the stability of the P4 cluster (and other phosphorus clusters with valence angles near 60°), we added a 60° angle correction term to Eq. 1 * .
* Note that, for the simulation of P/H systems with the 60° angle correction, one needs to use the force field file with 60° angle correction and recompile the LAMMPS package with our modified source file, reaxc_valence_angles.cpp. We verified that the 60° angle correction term would only affect the properties of P/H systems with valence angles near 60°. Therefore, for the simulation of condensed phases (either pristine or defected) and phosphorus hydride molecules, the original software of LAMMPS package can be used with the force field file without 60° angle correction. Because the lgvdw term was included in the ReaxFF, the "pair_style" command in the input file of LAMMPS should be: pair_style reax/c NULL lgvdw yes.
In section 3.2.2, it is demonstrated that the accuracy of cluster formation energies was significantly improved by the addition of the 60° angle correction term.

DFT Training of Force Field
The ReaxFF parameters for P/H systems were optimized using a modified version of the evolutionary algorithms (EA) software suite OGOLEM, 52,53 which is able to globally optimize ReaxFF parameter sets with high parallel efficiency. Based on DFT calculations for bulk black phosphorus, pristine and defected black phosphorene, blue phosphorene, phosphorus hydride molecules and phosphorus clusters, ReaxFF parameters were generated for P-P and P-H bond energies, P-P-P, H-P-P and H-P-H valence angle energies and for H-P-P-P and H-P-P-H torsion energies.
The parametrization of ReaxFF for P/H systems consisted of following steps:

Parameterization and Validation of ReaxFF
Our final fitted, global optimized ReaxFF for P/H systems is given in Tables 1-7. The potential form is given in Eq. 1 (a detailed description of all terms can be found in Refs. 26,27,31 ).
Unless otherwise stated, all ReaxFF results in the following discussion refer to our global optimized ReaxFF parameter set.

Relative stabilities for bulk black phosphorus, black and blue phosphorene
For ReaxFF to accurately describe phosphorus in the condensed phase, descriptions for different crystalline phases should be included in the DFT training set. Relative stabilities of the black phosphorus crystal as a function of unit cell volumes and relative stabilities of the both black and blue phosphorene as a function of unit cell in-plane areas were calculated. In general, the ReaxFF model gives a good description of lattice parameters of all three crystal phases (see Table 8) and shows a good consistency of the crystal structures of these crystal phases (see Fig.   1).
In particular, the DFT and ReaxFF results of cohesive energies are compared to SW results and experimental data in Table 9. The equilibrium cohesive energy of bulk phosphorus used in the fitting procedure was the experimental data 54 of -3.26 eV rather than the value computed from DFT (-3.43 eV). ReaxFF predicts a black phosphorus cohesive energy of -2.91 eV. The cohesive energy of black phosphorene calculated by ReaxFF is -2.84 eV, which slightly underestimates the DFT result of -3.35 eV. Still, ReaxFF provides a much better prediction of cohesive energy of phosphorene than that of SW potential, which underestimates the cohesive energy of phosphorene by an order of magnitude. ReaxFF are able to correctly reproduce the relative order of stability of three crystal phases (shown in Table 9). In Fig. 2(a) and Fig. 2 the results from ReaxFF correctly describe the relative stabilities of bulk black phosphorus for a broad range of cell volume, as well as that of black phosphorene for a broad range of cell area. In the training set, not all the data can be fitted equally well. For blue phosphorene (Fig. 2(c)), ReaxFF slightly overestimates the in-plane area of the unit cell, leading to a small offset of the energy profile of the relative stability. Given that no existing force field can describe the properties of blue phosphorene, the present ReaxFF may represent a major step forward.

Relative stabilities of phosphorus clusters
For ReaxFF to provide accurate description of phosphorus in clusters, the geometries and formation energies of P clusters of sizes 4,5,6 and 8 atoms are included in the training set. The formation energies per atom of clusters, , defined by where is the energy of the relaxed phosphorus cluster with n atoms, ℎ (bulk) is the cohesive energy of the bulk black phosphorus. As can be seen from Fig. 3, ReaxFF is capable of providing a very good description of the geometries of P clusters. Table 10 shows that the cluster formation energies per atom calculated by ReaxFF with 60° correction agree well with the DFT results. It is intriguing that a simple 60° angle correction term is able to provide such a notable improvement in terms of cluster formation energies.

Potential energy curves for phosphorus hydride molecules
Data for selected phosphorus hydride molecules was also included in the training set to train the P-H interactions and to enhance the transferability of the ReaxFF for P/H systems. To include DFT data for P-H, P-H bonds, dissociation profiles were determined from DFT calculations for phosphine, P2H2 and P2H4 molecules. The bond dissociation profiles were generated from the equilibrium geometries of these molecules by changing the bond length from the equilibrium value while allowing other structural parameters to relax, which are shown in Fig. 4 (a-c). Only the lowest-energy states (singlet, triplet or quintet depending on geometry) were included in bond dissociation profiles. The DFT and ReaxFF curves are shown in Fig. 4 (a-c).
To include DFT data for P-P-P, H-P-P and H-P-H valence angles, P3H5, P2H2 and phosphine molecules were used, respectively. Following the same procedure of constructing the bond dissociation profiles, P3H5, P2H2 and phosphine molecules were geometry optimized to create reference states. Afterwards the valence angles were modified while other structural parameters were optimized. The resulting angle distortion curves are shown in Fig. 4 (d-f).
Energy profiles for torsion distortion of H-P-P-H in the P2H4 molecule and of H-P-P-P in the P4H2 molecule were also included in the training set. The torsion distortion curves were generated from the equilibrium geometries of these molecules by changing the relevant torsion angle from the equilibrium value while allowing other structural parameters to relax, which are shown in Fig. 4 (g-h).
In Fig. 4 (a, d, f, g, h), it is visible that the interactions between phosphorus and hydrogen atoms in phosphorus hydride molecules are well reproduced with ReaxFF. For the interactions between phosphorus atoms in phosphorus hydride molecules (see Fig. 4 (b, c, e)), agreement between the ReaxFF and DFT results are not perfect, because the crystal phases of phosphorus were prioritized over the phosphorus hydride molecules in ReaxFF. The depth of the ReaxFF potential well in Fig. 4(b) is shallow, in order to offset the errors in cohesive energy for bulk black phosphorus (cf. Table 9) and the ultimate strength of black phosphorene in zigzag direction (cf. Fig. 7).

Defects for black phosphorene
Properties and applications of 2D materials are strongly affected by defects, 55 where is the energy of the defected phosphorene (geometry optimized) with n phosphorus atoms, ℎ (black) is the energy per atom of the black phosphorene. Fig. 5 shows that ReaxFF performs very well in predicting defects geometries of all three types with respect to DFT calculations, whereas the SW potential fails to predict the structure of all three type of defects. From Table 11, ReaxFF provides a good description of the defect formation energy of single vacancy and double vacancy in phosphorene, as well as the relative stability between single vacancy and double vacancy. The formation energy of Stone-Wales defect is overestimated by 36% by ReaxFF. By comparison, for SW potential, the formation energies of single and double vacancy are seriously underestimated (see Table 11) and the Stone-Wales defect is unstable (see Fig. 5), leading to an erroneous 0 eV formation energy. Compared to SW potential, ReaxFF provides a significant improvement in describing different types of defects in phosphorene.

Adatoms for black phosphorene: a transferability test
Due to its 2D nature, the large surface area to volume ratio of a black phosphorene nanosheet leads to a high chemical activity to foreign atoms. Thus the accurate description of where / ℎ is the total energy of phosphorene with/without adatoms and is the energy of the isolated atom. Fig. 6 shows that ReaxFF agrees very well with DFT calculations for predicting the adsorption structures of P and H adatoms. By contrast, the SW potential overestimates the bond length between P adatom and upper P atoms in black phosphorene. Without P-H interactions, SW potential is not capable to describe the H adatoms for black phosphorene. In Table 12, it can be seen that ReaxFF provides a good description of the binding energy of P adatom and slightly overestimates the binding energy of H adatom.
However, the SW potential underestimates the binding energy of P adatom by an order of magnitude. Overall, ReaxFF provides a good description of P and H adatoms on black phosphorene. Since the structures and formation energies of P and H adatoms for black phosphorene were not included in the training set, these results indicate a good transferability of the ReaxFF for P/H systems.

ReaxFF
In Table 13, the Young's modulus of black phosphorene in armchair and zigzag directions calculated by ReaxFF and SW potential are compared to DFT results. ReaxFF performs fairly well in reproducing the Young's modulus of black phosphorene in both directions, while SW potential underestimates the Young's modulus of black phosphorene in both directions. Fig. 7 (a) shows the stress-strain curves of black phosphorene in zigzag and armchair directions calculated with DFT, ReaxFF and SW potential. For zigzag direction, ReaxFF is able to capture the modulus change as the strain increases, providing a reasonable agreement in ultimate strength and failure strain. However, the SW potential severely underestimates the ultimate strength and failure strain in the zigzag direction. For armchair direction, ReaxFF overpredicts the failure strain while SW potential underpredicts it. The ultimate strength of black phosphorene in the armchair direction is slightly overestimated by ReaxFF, while it is severely underestimated by SW potential. ReaxFF yields a smaller failure strain at 300 K than 1.0 K for both the zigzag and armchair directions (see Fig. 7 (b)). Generally, ReaxFF gives a much better representation of the mechanical response of pristine black phosphorene over the SW potential.

Effect of Defects on the Mechanical Response of Black Phosphorene
Stress-strain curves of defected black phosphorene in the armchair and zigzag directions calculated with ReaxFF at 1 K are shown in Fig. 8 (a) and Fig. 8 (b), respectively. For armchair direction, black phosphorene with single vacancies shows a larger reduction in the failure strain than black phosphorene with double vacancies (keeping defect density the same), even though the double vacancy has a higher formation energy than the single vacancy. The reduction in the failure strain induced by Stone-Wales defect is in between that of single and double vacancy.
The Young's modulus in the armchair direction is more or less unaffected by all three types of defects. For zigzag direction, all three types of defects reduce the failure strain by about 50%.
Only minor reduction in the Young's modulus in the zigzag direction is induced by all three types of defects. Thus, the mechanical response of black phosphorene in the zigzag direction is more sensitive to defects than that for the armchair direction.
To understand these phenomenon, the structural deformation and stress distribution of defected black phosphorene under tension (ε=0.13) in the armchair ( Fig. 9 (a-c)) and zigzag ( Fig.   9 (d-f)) directions were analyzed. For armchair direction ( Fig. 9 (a-c)), stress at the single vacancy is more concentrated than that of double vacancy and Stone-Wales defect, due to the ReaxFF agrees fairly well with DFT results. 57 This clearly shows that ReaxFF for P/H systems provides a robust tool to study the effect of defects on the mechanical response of black phosphorene on a much larger space and time scale compared to DFT.

Thermal stability of phosphorene nanotubes
Similar to carbon nanotubes, the electrical and optical properties of the one-dimensional phosphorus nanotube (PNT) are chirality dependent and can be tuned by strain and size, [58][59][60][61][62][63] shedding light on its potential applications in transistors, strain sensors and photodetectors. Thus the accurate description of the properties of PNTs is important for ReaxFF. Two types of PNTs were designed by wrapping up a phosphorene sheet along the zigzag and armchair directions, i.e.
(m, 0) zigzag PNTs and (0, n) armchair PNTs. 61,64 Fig. 10 shows that compared to SW potential, ReaxFF provides a more accurate description of the cohesive energies change of the zigzag PNTs and armchair PNTs with respect to their sizes. SW potential underpredicts the cohesive energies of PNTs by an order of magnitude, indicating that SW potential could seriously underestimate the thermal stability of PNTs. 65 The phase diagrams for thermal stability of the zigzag PNTs and the armchair PNTs with varying temperatures and wrapping vectors of the nanotube are shown in Fig. 10 (c) and Fig. 10 (d), respectively. It is seen that SW potential strongly underpredicts the thermal stability of PNTs, compared to ReaxFF. Guan et al. 60

Concluding Remarks
We present a reactive force field (ReaxFF) for phosphorus and hydrogen, which gives an accurate description of the chemical and mechanical properties of pristine and defected black phosphorene. A 60° correction term is added which significantly improves the description of phosphorus clusters. ReaxFF for P/H is transferable to a wide range of P/H systems including bulk black phosphorus, blue phosphorene, phosphorus clusters and phosphorus hydride molecules. Emphasis has been placed on obtaining a good description of mechanical response of black phosphorene with different types of defects. Compared to SW potential, ReaxFF for P/H systems provides a notable improvement in describing the cohesive energy, mechanical response of pristine and defected black phosphorene and the thermal stability of phosphorene nanotubes.
We observe a counterintuitive phenomenon that single vacancies weaken the black phosphorene more than relatively more unstable double vacancies. It is shown that the mechanical response of black phosphorene is more sensitive to defects in the zigzag direction than the armchair direction. Straightforward extensions to the heterogeneous systems, including oxides, nitrides, etc., enable the ReaxFF parameters for P/H systems to build a solid foundation for the simulation of a wide range of P-containing materials.

Supporting Information
The ReaxFF file for P/H systems with or without 60° angle correction and the modified source file (reaxc_valence_angles.cpp) are available. Fig. 1. Crystal structures of bulk black phosphorus, black phosphorene and blue phosphorene calculated by DFT and ReaxFF.         this table and Tables 2-6 can be found in Refs. 26,27,31 . For H, parameters from Ref. 66 were used.  Table 4 Bond energy and bond-order parameters for the P-P, For the H-H bond, parameters from Ref. 66 were used.  The parameters of 60° angle correction for P-P-P are designed to improve the description of phosphorus clusters with ReaxFF, explained in section 2.2. Table 7 Torsion angle parameters