Diamond having unique physicochemical properties is considered as one of the most promising materials for next generation semiconductors.^[1] It is well-known that p-type diamond have been successfully obtained by boron (B) doping,^[2,3] however, the fabrication of n-type diamond with a shallow donor level is extremely difficult. In order to achieve the shallow doping for n-type diamond, kinds of donors have been attempted, such as lithium,^[4,5] nitrogen,^[6–9] phosphorus,^[10–13] oxygen,^[14] and sulfur (S),^[15–17] etc.

Doping with S was reported to be an effective method to obtaining the n-type diamond,^[15,18] then a number of computational approaches have been employed to estimate the donor model and donor level value.^[19–31] One possibility of the origin of the donor levels may be substitutional S (S_sub),^[1] and the most stable configuration of S_sub in diamond is generally considered to S atom and four carbon (C) neighbors with symmetry,^[26–29] and the configuration is considered as the less stable structure.^[27,28] Then, it has been achieved kinds of donor level values between −0.15 eV to −1.63 eV for S_sub in diamond with symmetry.^[19–29] However, nearly all of the supercells were relaxed at a fixed volume using structural optimizations to simulate the potential field of diamond, and most of supercells calculated only contained 64 atoms to minimize the computational cost.^[30] It is assumed that the donor center structure of S_sub and the donor level value are qualitatively affected by the supercell size (namely, the doping concentration) and optimization methods.

Since the molecular dynamic (MD) simulations can exhibit the changing of the structure parameter of S_sub in diamond at room temperature as the function of time, it has been successfully used to predict the new structure of ammonium azide under high pressure,^[32] consider the tanglesome rotation of ammonium cation in ammonium nitrate^[33] and understand the kinetic and thermodynamic properties of the closing and opening of a base pair.^[34] In this paper, both the ab initio MD simulations and structural optimizations are considered for the structural model of S_sub in diamond with the larger supercells. We find that the symmetry is the stable structure for S_sub in diamond, and the minimum supercell is 216-atom.

All calculations performed in this paper are based upon ab initio density functional theory simulations, as implemented in Vienna Ab-initio Simulation Package (VASP) code.^[35] The projector-augmented wave (PAW)^[36,37] method is adopted and the cutoff energy of 520 eV is set. The Perdew–Burke–Ernzerhof (PBE) function based on the generalized gradient approximation (GGA)^[37] is employed to describe the exchange and correlation effects of electrons. The MD simulations are implemented both in the isobaric isothermal (NPT) ensemble at 300 K, and atmospheric pressure and canonical (NVT) ensemble at 300 K. The optimized lattice parameters of pristine diamonds with 64-, 216-, and 512-atom supercells are 7.146 (3.573 × 2), 10.7189 (3.573 × 3), and 14.292 (3.573 × 4) Å, which are agree well with the experimental single cell lattice value of 3.567 Å,^[38] and then one C atom is replaced with a S atom for MD simulations. In MD simulations, the Brillouin zone is sampled using the Monkhorst–Pack scheme, with mesh of 3 × 3 × 3, 2 × 2 × 2, and 1 × 1 × 1 for 64-, 216-, and 512-atom supercells, respectively. The time steps of 1.0 fs and total simulation times of at least 10 ps are used, and all the MD simulations achieve thermalisation after 3 ps by examining the system pressure and temperature. After MD simulations, we optimize the averaged structures and calculate the total energies of the three supercells with k mesh of 7 × 7 × 7, 5 × 5 × 5, and 3 × 3 × 3, respectively.

First, the MD simulations both in the NVT and NPT ensembles at 300 K are performed to investigate the S_sub structure. The main difference between the NVT and NPT ensembles is that the volume, side length and side angle of supercell for the former are restricted constant, while those in the latter are allowed to optimize. Figure 1 show the bond lengths and bond angles between the S dopant and four adjacent C atoms with the 64-atom supercell. For NVT MD simulation steps shown in Figs. 1(a) and 1(c), after the simulation reaches equilibrium, one S–C bond length increases to 1.95 Å, and the other three bonds keep at ∼1.72 Å. The lattice parameters in the last 1.0 ps are averaged, and a symmetry donor center is obtained by further geometry relaxation. Additionally, the longer S–C bond length and lager C–S–C bond angle appear alternately, indicating that the system hops from one configuration to the other three equivalent structures. It is consistent with previous results, which considered that pattern was more stable.^[27,28] As a contrast, for NPT MD simulation in Figs. 1(b) and 1(d), the four S–C bond length values are almost equal as the time increases, and a symmetry donor center is obtained by averaging over the lattice parameters after 9.0 ps. Moreover, the hoping of S atoms in the donor center of supercells can also be observed. It is concluded that ( ) symmetry is more stable for NVT (NPT) ensemble for S_sub in diamond with the 64-atom supercell. Considering the opposite results of NPT and NPV MD simulations, the supercell is increased in the following.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. The distance between S and four adjacent C atoms (d_SCi, i = 1, 2, 3, 4) as a function of the MD simulation time in (a) NVT and (b) NPT ensembles, and the six bond angles in donor center ( CiSCj, i, j = 1, 2, 3, and 4) as a function of the MD simulation time in (c) NVT and (d) NPT ensembles. The calculations are performed with 64-atom supercell.

Figure 2(a)–2(d) exhibit the S–C bond length and C–S–C bond angle in 216-atom supercell as a function of MD simulation time. Both the NVT and NPT ensemble simulations suggest that the and donor center appears alternately. Furthermore, the S–C bond length and C–S–C bond angle as a function of MD simulation steps in the 512-atom supercells in Figs. 3(a)–3(d) show that at the most of NVT simulation time, the donor center performs as symmetry, and during almost the entire NPT simulation, the donor center is symmetry. Thus, for the NVT simulation, the symmetry donor center transforms into structure, when the supercell is larger than 216, namely, symmetry donor center for S_sub in diamond with a larger supercell is more possible. However, for NPT ensemble calculations, symmetry always appears with the 64- and 512-atom supercells, while the and structures hop alternatively with the 216-atom supercell.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. The distance between S and four adjacent C atoms (d_SCi, i = 1, 2, 3, 4) as a function of the MD simulation time in (a) NVT and (b) NPT ensembles, and the six bond angles in donor center ( CiSCj, i, j = 1, 2, 3, and 4) as a function of the MD simulation time in (c) NVT and (d) NPT ensembles. The calculations are performed with 216-atom supercell.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. The distance between S and four adjacent C atoms (d_SCi, i = 1, 2, 3, 4) as a function of the MD simulation time in (a) NVT and (b) NPT ensembles, and the six bond angles in donor center ( CiSCj, i, j = 1, 2, 3, and 4) as a function of the MD simulation time in (c) NVT and (d) NPT ensembles. The calculations are performed with 512-atom supercell.

To explain it, the three lattice angles of S_sub in diamond with 64-, 216-, and 512-atom supercells and pristine diamond with 512-atom supercell during NPT MD simulations are given in Fig. 4. For 64-atom supercell, the averaged lattice angle of 91.0656° or 88.9465° is far away from that of cubic structure. Especially, the averaged angle values of 216- and 512-atom supercell are nearly 90°, which are close to the cubic supercell. Thus, the result obtained from NPT MD calculations with 64-atom supercell is not reliable. Consequently, from the NPT MD results, we also obtain the same conclusion that the minimum supercell to obtain converged values of S_sub in diamond is 216-atom. The different symmetry of S_sub in diamond are obtained for NVT and NPT ensembles with 64-, 216-atom supercell, while the same symmetry is realized with the 512-atom supercell. This is because that during the MD simulation with limited number of atoms, there are still interaction between the boundary atoms and S dopant. Thus, the different restrictions for the atomic location in the boundary between the NVT and NPT ensembles result in the different geometry of impurity centers. When the simulated supercell is large enough, the influence of the impurity atom at the supercell boundary atoms is negligible, and thus the difference between the NVT and the NPT ensemble is also negligible.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. The lattice angle of Ssub in diamond with (a) 64-, (b) 216-, and (c) 512-atom supercells and (d) pristine diamond with 512-atom supercell as a function of NPT MD simulation time.

Then, the structures under a series of density functional theory (DFT) structural optimizations are obtained to further confirm the stable S_sub structure. First, the ionic positions, the lattice shape and parameters of and structures are fully relaxed, namely, the pressure is zero. The relative energies for fully relaxed structures show that the structure is more stable for 64-, 216-, and 512-atom supercells (Table 1). Moreover, the lattice parameters are closer to those of cubic, when the supercell is larger than 216-atom. Then, the relative atomic positions in and structures are kept invariant, and the lattice parameters are adjusted to be cubic, only the ionic positions are relaxed. Four representative lattice parameters of 7.206, 7.203, 7.146, and 7.134 Å are selected for the supercells, corresponding to the structure with zero pressure, structure with zero pressure, relaxed pristine diamond, and experimental value of pristine diamond.^[39] As shown in Table 2, in the 64-atom cubic structures, the energy of donor center is lower than that of , just as discussed in literature.^[27,28] However, the pressures of structures are higher than those of with the same lattice parameters, denoting that the conclusion need to be further versified. Further calculations with 216- and 512-atom supercell indicate that the symmetry is more stable. Moreover, the deviation of lattice angles has a limited impact on the energy difference between and supercell at 512-atom condition, but it significantly reduces the energy difference at 216-atom condition.

Table 1.

Table 1.

Table 1. The lattice parameters of fully relaxed supercells with and donor centers with 64-, 216-, and 512-atom supercells. . 64-atom Angle/(°) 91.0702 89.9997 89.9997 90.3020 90.3018 90.3017 Side length/Å 7.119 7.244 7.244 7.205 7.205 7.205 Relative energy/eV 0.2243 216-atom Angle/(°) 89.7028 90.0000 90.0000 89.9465 90.0000 90.0000 Side length/Å 10.707 10.761 10.761 10.745 10.761 10.761 Relative energy/eV 0.1803 512- atom Angle/(°) 90.1213 90.0003 90.0002 89.9800 90.0215 89.9787 Side length/Å 14.285 14.315 14.316 14.307 14.307 14.307 Relative energy/eV 0.1291

Table 1. The lattice parameters of fully relaxed supercells with and donor centers with 64-, 216-, and 512-atom supercells. .

Table 2.

Table 2.

Table 2. The parameters of relaxed cubic supercells with and donor centers with 64-, 216-, and 512-atom supercells. ( pressure ( ). . Side length/Å Relative energy /eV Pressure (kB, Pressure (kB, /kB 64- atom 7.206 0 –0.1514 –5.58 0.00 –5.58 7.203 0 –0.1498 0.00 5.59 –5.59 7.146 0 –0.1185 103.23 109.08 –5.85 7.134 0 –0.1119 126.18 132.06 –5.88 216-atom 10.745 0 0.0686 –2.00 0.00 –2.00 10.743 0 0.0693 0.00 1.99 –1.99 10.719 0 0.0844 29.94 31.96 –2.02 10.701 0 0.0923 52.24 54.29 –2.05 512-atom 14.306 0 0.1181 –0.94 0.00 –0.94 14.305 0 0.1184 0.00 0.88 –0.88 14.292 0 0.1230 12.24 13.15 –0.91 14.268 0 0.1309 34.36 35.21 –0.85

Table 2. The parameters of relaxed cubic supercells with and donor centers with 64-, 216-, and 512-atom supercells. ( pressure ( ). .

The band structures of 512-atom with symmetry in neutral and positive charge state (named as S_ and S_ are shown in Figs. 5(a) and 5(b), respectively. The former is spin-degenerated, while partial bands near the Fermi level for the latter are spin-splitting. For the case of S_ (Fig. 5(a)), the donor level is about 1.1 eV with respect to the conduction band minimum (CBM). For S_ (Fig. 5(b)), the occupied spin-up donor level is about 1.6 eV below the CBM, and the spin-down donor level is un-occupied. It denotes that the S_sub in diamond is a deep donor. The partial charge density (PCD) for donor level at Fermi surface in the inset denotes that the lone pair electrons and the electrons in the nearest S–C bond have a major contribution to the donor level, as the discussion of substitutional phosphorus in diamond.^[12] In order to have the further insight into the bonds near the S atom, the partial structures of electron localization function (ELF), which has widely been used to reveal the bonding, lone electron pairs, atomic shell structure,^[40] is shown in Fig. 6(a). At the same time, two-dimensional contour plots of ELF for the plane crossing the S, C₁, and C₂, and S, C₃, and C₄ atoms are shown in Figs. 6(b) and 6(c), respectively. One can see that the S atom is bonded with four neighbor C atoms, with the remaining two electrons forming a lone-pair orbital between the C₃ and C₄ atoms. In detail, the bonding orbital between the S–C₃ and S–C₄ are repulsed from the lone-pair electron. Furthermore, we present the charge density distribution and its two-dimensional contour plots of the plane crossing S, C₁, and C₂, and S, C₃, and C₄ in Figs. 6(d)–6(f). The charge localization between C₃ and S (or C₄ and S) is more evident than that between C₁ and S, indicating that the lone pair electrons have a major effect on the two nearest S–C₃ and S–C₄ bonds, and further impact the electron distribution between two C atoms those in the same plane with S, C₃ and C₄ (labeled in Fig. 6(e)), while in other direction, the impact of electron distribution is neglected (Fig. 6(f)).

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. The band structures of Ssub in diamond with 512-atom supercell with symmetry in (a) neutral and (b) positive charge states. The insert shows the isosurface plot for the donor level near Fermi energy. The brown balls are C atom, and red is S atom.

Fig. 6.

Figure Option ViewDownloadNew Window

Fig. 6. Partial structures of [(a), (b), (c)] ELF and [(d), (e), (f)] charge density of Ssub in diamond in 512-atom supercell with symmetry. 2D ELF contour plots of the plane crossing (b) S, C3, and C4 atoms, and (c) S, C1, and C2 atoms. 2D charge density contour plots of the plane crossing (e) S, C3, and C4 atoms, and (f) S, C1, and C2 atoms. The brown balls are C atom, and red is S atom. Blue to red indicates the gradually increased charge localization.

In conclusion, the stable donor center structure of S_sub in diamond is studied by a series of MD simulations and structural optimizations using the 64-, 216-, and 512-atom supercells. The center is the more stable configuration in the cubic supercell with 64 atoms, while the center is more energetically stable with a larger supercell. Our results indicate that 216-atom supercell is the minimum supercell to obtain the accurate structure of S_sub in diamond. The S_sub with donor center has a deep donor level, which mainly comes from the lone pair electrons located between two S–C bonds. It is expected that more studies on S-doped diamond and co-doped diamond with S and other donor impurities are needed.

The calculations were performed in the High Performance Computing Center (HPCC) of Jilin University.