Diffusion quantum Monte Carlo study of martensitic phase transition energetics: The case of phosphorene

Phosphorene is the two-dimensional counterpart of layered black phosphorus, and a detailed investigation of this particular material is of great interest for various technological applications because the material is predicted to exhibit unique electronic and optical properties.^1,2 The most commonly discussed phase is black phosphorene, and the material has already been synthesized and has shown good performance in a few applications for solar energy conversion.^3,4 In additional to the black phosphorene phase, more than ten different phases of phosphorene have been predicted to exist in literature.^1,2,5 Among these different phases of phosphorene, the so-called blue phosphorene phase is predicted by density functional theory (DFT) to have a similar honeycomb structure and ground state energy as the black phosphorene phase yet possessing distinctively different electronic properties.² Despite these theoretical predictions, however, it has yet to be experimentally synthesized. In this work, we investigate the possibility of martensitic phase transition from the common black phosphorene phase to blue phosphorene phase by employing diffusion Monte Carlo (DMC) calculation⁶ to obtain accurate reaction energetics along the transition pathway.

Several works have explored the possibility of diffusionless phase transitions between different phases because such phase transition can be facilitated through an application of external strain on the material in experiments.² Determination of the relevant reaction coordinate and obtaining accurate reaction energetics along the reaction coordinate are two important requirements in theoretical investigations of such transitions. Generally, minimum energy path (MEP) defines the reaction coordinate and transition states of infrequent transition events such as chemical reactions when the entropic contribution is not dominant, and various numerical methods have been proposed over the years for locating the MEP.⁷ These include chain-of-states methods such as the nudged elastic band (NEB) method⁸ and the string method⁹ as well as other types of numerical approaches, which are based on local curvatures of the potential energy surface such as the dimer method¹⁰ and the activation-relaxation technique¹¹ for locating transition states. Unlike for chemical reactions, however, phase transitions involve not only changes of atom coordinates but also changes in the lattice vectors. Thus, the lattice degrees of freedom need to be included in the search space for locating the MEP. For the particular case of the phosphorene phase transition in this work, we are concerned with a martensitic (diffusionless) phase transformation in which the phase changes involve cooperative/concerted movement of all atoms without involving long-range diffusion of atoms. In this work we employ the variable cell NEB method¹² to locate the MEP in the extended configuration space, which includes lattice coordinates so that the reaction coordinate can be defined also for the phase transition.

Another important aspect of the present study is the electronic structure calculation that is used to obtain the energetics of the transformation for large extended systems. For theoretical investigations of extended systems, DFT calculations have become quite common because of their balanced accuracy and computational affordability. However, while geometries appear to be well described, reaction energetics is often poorly described by DFT calculations when employed with most common exchange-correlation (XC) approximations.¹³ In this regard, a higher-level of first-principles approaches is desired for accurately characterizing reaction energetics. Among several first-principles many-body approaches, diffusion quantum Monte Carlo (DMC) method is promising for large systems because of its favorable scaling of the computational cost with the system size. When systems do not show strong correlation behaviors (e.g., no chemical bond formation/dissociation), DMC’s accuracy often rivals that of CCSD(T) method with a large basis set even with the fixed node approximation using only one Slater determinant.^14–17 At the same time, application of explicitly many-body methods such as DMC calculation to extended systems like phosphorene is not straightforward because of the many-body finite-size error that is associated with using the periodic boundary conditions to describe infinitely large systems.¹⁸ There are already extensive literatures on how the many-body finite-site error can be corrected in the context of DMC,^19–23 and the topic remains an important area of further pursuit. The finite-size error in DMC is very different from the finite-size error in DFT, which is entirely due to a finite sampling of Brillouin zone (BZ). Detailed discussion of the origin of the finite-size error in DMC can be found, for example, in Ref. 23. In our present study on the martensitic phase transition of the extended phosphorene, we examined both the extrapolation approach²⁴ as well as the DFT-based correction scheme by Kwee et al.,²⁵ along with employing the twist-averaging approach for Brillouin zone integration.²⁶ 

The variable cell nudged elastic band (VCNEB) method¹² is used to locate the minimum energy path (MEP) for the phase transition between blue and black phosphorenes (Figure 1). VCNEB is a nudged elastic band (NEB) method that allows unit cell vectors to be varied in addition to the atomic positions.¹² This builds on the NEB framework by including unit cell vectors in the configuration space during the search for the MEP. The configuration space can be expressed as

where ε_ij is the elements in the finite strain tensor $ε¯¯$ which is defined through $h=1+ε¯¯h0$⁠. h is the matrix of lattice vectors, and h₀ is a reference matrix of lattice vectors, which was taken to be that of the blue phosphorene. In the VCNEB calculation, the configuration space is represented by a total of 9 + 3N_atom variables. The convergence threshold for the VCNEB calculation was set to 0.005 eV/Å, and we use a total of 47 images to represent the path between the initial and final states for the phase transition. We performed the VCNEB calculation by interfacing the USPEX code²⁷ with Quantum Espresso code²⁸ such that the underlying electronic structure is derived from DFT calculation. The Perdew-Burke-Ernzerhof (PBE) exchange-correlation (XC) functional²⁹ was employed, and Vanderbilt ultrasoft pseudopotentials³⁰ were used for describing core-valence electron interactions. The planewave cutoff for Kohn-Sham wavefunctions and the smooth charge density cutoff were 20 Ry and 200 Ry, respectively. The Monkhorst-Pack k-point grid of 8 × 8 × 1 was used for sampling the Brillouin zone using the unit cell with four atoms. For stationary states along the located MEP, further DFT calculations were performed with PBE0 hybrid XC functional for additional comparison. In this case, Troullier-Martin norm-conserving pseudopontials³¹ were used with the planewave cutoff for Kohn-Sham wavefunctions of 80 Ryd.

Fixed-node diffusional Monte Carlo (DMC) calculations were performed using the QWALK code³² to obtain accurate energetics at the stationary points along the MEP located. Fermion nodes here were given by Slater determinant of Kohn-Sham orbitals, which were obtained with PBE functional using the CRYSTAL code.³³ 3s² 3p³ electrons of phosphorous atoms are treated explicitly as valence electrons using the BFD pseudopotential,³⁴ where the Ne core is treated within the pseudopotential. For the basis sets, we use a modified Burkatzki-Filippi-Dolg (BFD) basis, and its mathematical form is provided in supplementary material. Locality approximation is used for dealing with the non-local part of the pseudopotential.³⁴ Importance sampling was introduced with the trial wavefunction, which was given by the Slater determinant multiplied by one-body and two-body Jastrow correlation factor. Variational quantum Monte Carlo (QMC) calculation was first used with the variance minimization³⁵ to obtain the trial wavefunctions. As implemented in the QWALK code, the DMC algorithm by Umrigar, Nightingale, and Runge³⁶ was used with a modification to the branching part for keeping the number of walkers constant.³² We used the imaginary time step of 0.01 a.u. to perform the DMC calculations, and the convergence test is presented in supplementary material. In the recent 2016 work, Zen et al. have shown that DMC energies could suffer the size consistency problem when the time step is not sufficiently small.³⁷ The authors proposed an alternative DMC algorithm that significantly reduces the size-consistent problem, allowing the use of larger time steps. We, however, did not implement this new algorithm in this work but instead used the rather small time step (0.01 a.u.) with the above-mentioned modified Umrigar-Nightingale-Runge (UNR) algorithm.

Independent-particle methods such as DFT can employ the Bloch theorem to eliminate the finite-size error by obtaining the macroscopic limit via Brillouin zone (BZ) integration. Finite-size error is one of the major sources of error in quantum Monte Carlo calculation of extended systems.⁶ There are two sources contributing to the finite-size error: one-body and two-body errors. The one-body contribution to the error is related to BZ sampling and can be largely removed using twist averaging, an approach discussed in Ref. 26. We used 4 k-points (2 × 2 × 1) with a 64-atom simulation cell for the “twist.” We show the convergence test in Figure 2 for the blue phase and also for an intermediate transition state structure along the MEP because the latter was found to be metallic. The test shows that the 4 k-points are sufficient for the BZ sampling even for the smaller 16-atom simulation cell.

The two-body error is more problematic, and it originates from the artificial periodicity in the long-range electron-electron Coulomb interaction when extended systems are simulated using periodic boundary conditions. The most straight-forward approach (and most rigorous, in principle if there are no size-consistency errors associated with the DMC algorithm as discussed in Ref. 37 and above) is an extrapolation using increasingly larger simulation cells, as was utilized in a recent work for a two-dimensional material by Wagner and co-workers.²⁴ Assuming that phosphorenes can be considered as a two-dimensional system, one expects the total energy to scale as N^−5/4, where N is the system size according to the work by Drummond et al.²³ Fitting DMC energies as a linear function of N^−5/4, where N is the number of unit cells in the supercell, we can extrapolate the linear fit to the limit of N^−5/4 = 0 (i.e., an infinitely large simulation cell).

Alternatively, Kwee et al. proposed a scheme that is based on DFT calculations using a modified finite-size local density approximation (LDA) to XC functional for estimating the two-body error,²⁵ and it has thus far has been applied successfully to bulk systems.^38,39 A computational advantage of this approach is that finite-size correction can be approximated using DFT calculations rather than performing a series of more expensive QMC calculations for the extrapolation. Whereas the LDA XC functional was originally parameterized for an infinite jellium using QMC results,⁴⁰ this finite-size LDA XC functional is instead parameterized to fit finite size QMC calculations of jellium, thus the finite size error is intentionally retained in the modified XC functional itself. The Kohn-Sham single-particle Hamiltonian can be expressed as

where the first three terms are kinetic energy, ionic potential energy, and Hartree term, respectively, and the last term is the XC potential. This term is calculated as $Vxc∞r=δ(nrϵxc∞n)/δn(r)$ where $nr$ is the charge density, and $ϵxc∞n$ is obtained from the QMC calculation of the infinite size jellium. In the finite-size error LDA functional by Kwee et al., the XC energy has the form of $ϵxcFSn≡ϵxrs,L+εc(rs,L)$⁠, where r_s specifies the density of the jellium via $4πrs33≡1n$ and L denotes the linear size of the supercell (for example, in cubic cell, L = V^1/3). The finite size error then can be estimated using DFT calculations as

where the $E∞$ is the DFT total energy (which is converged with respect to k-point integration) using the original LDA XC functional, and E^FS(L) is the total energy using the modified finite-size LDA XC functional with the same k points as in the QMC calculation. This correction can then be used to correct the QMC calculated value to estimate the QMC energy in the limit of an infinite cell as

where $EQMCFS$ is the finite-size QMC result, and $EQMC∞$ is QMC result in the limit of infinite size. In our calculations, we treat phosphorenes as strictly a two-dimensional system, and we take the linear size L to be the square root of the area of the plane of the simulation cell that is parallel to the phosphorene (xy plane). We checked the dependence of $EQMC∞$ on the size of the vacuum region along the z direction and found it to be negligibly small (less than 0.0005 eV/atom).

We compared the above-mentioned extrapolation approach and the correction scheme using the modified finite-size XC functional within DFT calculations (Kwee, Zhang, Krakauer (KZK) correction). Supercells with 2 × 2, 3 × 3, and 4 × 4 cell sizes in blue phosphorene phase were considered here, and they correspond to N = 4, 9, and 16 of the N^−5/4 extrapolation. Figure 3 shows QMC total energies with and without the KZK correction as a function of N^−5/4. The extrapolation with the uncorrected total energies using the least-square fitting gives −178.9190 eV/atom in the limit of N^−5/4 = 0. The KZK-corrected total energy gives a similar value −178.9119 eV/atom when the largest simulation cell size with N = 4 × 4 = 16 is used, resulting in the difference of only 0.0071 eV/atom. Given the very small difference, which is certainly below the necessary accuracy for our work, we chose the computationally less demanding correction scheme using the KZK functional for calculating $EQMC∞$ using the simulation cell of N = 4 × 4 = 16 for the rest of the paper.

Previously, Zhu et al. estimated an upper bound of the transition energy barrier to be approximately 0.47 eV/atom using DFT calculation with PBE XC functional. The calculation assumed the out-of-plane displacement of one atom to be the reaction coordinate for the phase transition.² Here, we employed the variable cell nudged elastic band (VCNEB) method to locate the minimum energy path (MEP) between black and blue phosphorene phases as shown in Figure 4. The transition involves significant changes to atom positions as well as to lattice vectors as shown in Figure 5. The transition from black phosphorene to the first transition state (TS1) has the energy change of 0.308 eV/atom, according to DFT-PBE calculation. The cell parameter a changes from 3.30 Å to 3.48 Å, while the cell parameter b increases from 4.58 Å to 5.15 Å. The local minimum (LM) is only 0.024 eV/atom lower in energy compared to TS1 and 0.052 eV/atom lower in energy compared to the second transition state (TS2). These values indicate that the local minimum is highly unstable even at room temperature. The cell parameters in the region of the LM vary significantly while the total energy remains essential unchanged suggesting a structure with a small Young’s modulus for the LM structure. From TS1 to TS2, the cell parameter a decreases slightly from 3.48 Å to 3.27 Å, and the cell parameter b increases significantly from 5.15 Å to 5.90 Å. The energy decreases by 0.336 eV/atom from TS2 to the blue phosphorene phase. In this part of the transition the lattice parameters remain almost the same changing by only 5% from the cell parameters of TS2 while significant atomic displacements accompany the transition. DFT-PBE calculation yields the overall energy barrier of 0.336 eV/atom for the phase transition.

We performed DMC calculations at the five stationary structures along the MEP (see Figure 4). In particular, the DFT predictions of the near degeneracy between the black and blue phases as well as the existence of the metastable structure need to be verified. In DMC calculations of periodic systems, a large error may come from the finite size effects. We utilized two approaches to calculate and to reduce the finite size error as discussed in Theoretical methods section. Consistent with the DFT result, DMC calculation confirms that black and blue phosphorene phases have almost the same total energies where black and blue phosphorene structures yield −179.0611 ± 0.0016 eV/atom and −179.0647 ± 0.0015 eV/atom, respectively. Interestingly, the metastable LM structure is located energetically lower than the two TS structures accordingly to DMC calculation as well, supporting the DFT prediction. However, DMC predicts a much greater energy barrier than DFT for the overall transition. The DMC calculation yields an energy barrier of 0.4237 ± 0.0027 eV/atom for the transition from blue phosphorene to TS1, which is about 1.4 times the value predicted by DFT (0.308 eV/atom). The energy difference between TS1 and LM is 0.0283 ± 0.0024 eV/atom compared to 0.024 eV/atom in DFT. The energy difference between TS2 and LM is 0.0793 ± 0.0029 eV/atom compared to the DFT result of 0.052 eV/atom. The energy change from TS2 to the black phosphorene phase (0.4783 ± 0.0030 eV/atom) is again about 1.4 times the DFT value (0.336 eV/atom).

In recent years, PBE0 XC approximation⁴¹ within DFT calculations has become quite popular for investigating extended systems due to its superior performance compared to standard GGA approximations such as PBE. In terms of computational cost, PBE0 is approximately one-to-two orders of magnitude more expensive than PBE using planewaves as basis set but still much more affordable than performing DMC calculations. A number of studies on molecular systems have shown that the GGA functionals with some amount of exact exchange describe the reaction barriers significantly better as compared to PBE.^42,43 In the present case of phosphorene martensitic phase transition, there is also a sizable improvement with PBE0 in the direction toward DMC calculation result as shown in Table I. For both TS1 and TS2 energy barriers, PBE0 calculation predicts the values that are in-between DFT-PBE and DMC calculations. We also observe the relative energy of the local minimum is improved while the overall energy difference between the black and blue phases is essentially unchanged.

Most experimental studies on phosphorenes are focused on the black phosphorene phase^3,4 while we confirmed here the earlier DFT prediction² that the blue phosphorene phase is equally stable energetically by employing accurate quantum Monte Carlo calculations. Having shown that DFT calculations can be trusted for the qualitative description of the martensitic phase transition, we study the extent to which an application of external strain can change the transition energetics. Since the lattice vectors of these two phases differ significantly along the b-vector direction (see Figure 5), we examined to what extent an application of external strain might induce a more favorable transition energetics for the blue phosphorene. In their DFT work, Hu and Dong have previously studied the effect of applied strain and showed that the blue phosphorene phase can be made energetically more stable than the black phase.⁴⁴ Such an application of strain is experimentally conceivable as demonstrated for graphenes and nanowires as they are suspended over a trench.⁴⁵ We mimic such an experimental setting, by constraining the b-vector to be 100%, 110%, and 120% of the blue phosphorene (of all 47 images in the VCNEB calculation) so that we could examine how the MEP as well as the reaction energetics change. The structures along all the MEPs do not change substantially when compared to the case with no strains, except that a new local minimum appears between the black phosphorene and TS1 structures for the 120% case as shown in Figure 6. In all cases under the strain, the blue phosphorene phase is now approximately 0.1 eV/atom lower in energy than the black phosphorene phase. At the same time, the energy barrier (in DFT-PBE calculations) changes substantially from 0.34 eV/atom to 0.25 eV/atom to 0.12 eV/atom to 0.05 eV/atom as the b-vector is constrained to larger values.

In this work, we presented a theoretical study of martensitic phase transition of a two-dimensional phosphorene between its so-called black and blue phases. The minimum energy path was first determined using the variable-cell nudged elastic band method based on the atomic forces and stress tensors from DFT calculation so that the reaction coordinate can be defined without assuming a particular transition mechanism. The phase transition was found to consist of two distinct stages separated such that the initial atomic displacements are accompanied by a lattice deformation, while the reorganization occurring in the second stage is nearly exclusively due to the atomic displacements alone. DMC method was then used to examine the energetics of the phase transition, which was found at DFT level of theory to have a rather diffusive barrier⁴⁶ characterized by two transitions states. The DMC calculations support this interesting finding in addition to confirming that the two phases are essentially degenerate. Having demonstrated that DFT calculations predict the qualitatively correct behavior of the martensitic phase transition energetics, we also showed that the application of mechanical strain considerably alters the stability of the two phases of phosphorene, and the blue phase can be made much more stable than the black phase, with accompanied lowering of the energy barrier for the phase transition.

See supplementary material for the convergence test of DMC energy with respect to the imaginary time step and the basis set for phosphorous atoms.

We thank Lucas K. Wagner (UIUC) for helpful discussions related to quantum Monte Carlo calculations. This material is based upon work supported by the National Science Foundation under Grant No. DGE-1144081. We thank Lawrence Livermore National Laboratory for computational resources.