Crystal structure optimisation using an auxiliary equation of state

The standard operating procedure for computational investigations in solid-state chemistry is to begin with a crystal structure — obtained either from diffraction studies or through chemical analogy — and to optimise the lattice shape, volume, and internal positions to minimise all forces. It is from this equilibrium crystal structure (athermal for the majority of electronic-structure approaches) that the full range of material response functions (e.g., elastic, dielectric, magnetic) is calculated.¹ 

The optimisation of a crystal structure may require hundreds of self-consistent field iterations across a series of ionic configurations.² The most robust approach to optimisation is the calculation of an equation of state (EoS) for the material, relating the unit cell dimensions to energy and pressure.³ This is based on a series of calculations at differing volumes, where ideally the shape and internal positions are optimised at each point. The simplest case is a cubic lattice with high internal symmetry, e.g., rocksalt, where the only degree of freedom is the volume and computing the EoS reduces to a series of single-point calculations. For a triclinic cell, the lengths, angles, and internal positions in principle all require optimisation. While it is sometimes possible to directly optimise the cell volume by simultaneously minimising the stress tensor of the unit cell, this approach can run into artifacts when using plane-wave basis sets (i.e., Pulay forces) unless an iterative procedure is employed.⁴ 

It has become commonplace to use an “equilibrium” crystal geometry, determined using one exchange-correlation functional within density functional theory, for a “single-shot” higher-level calculation performed to give a more accurate electronic structure. This methodology has been applied to the calculation of properties as diverse as workfunctions, electronic bandgaps, optical properties, and defect formation energies.^5–12 The implicit assumption is that the qualitative behaviour is insensitive to small differences in the local structure. The approximation will fail where the electronic structure (chemical bonding) of a system is poorly described at the initial level of theory, e.g., the treatment of Mott insulators such as NiO within the local-density approximation (LDA).¹³ 

In this contribution, we outline a simple procedure for the rapid volume optimisation (RVO) of crystal structures. It takes advantage of the similarity in the pressure-volume relationship for a given material between different theoretical approaches, here being exchange-correlation (XC) functionals within density functional theory. Where an EoS is known for one functional, the equilibrium volume for another functional can be predicted with reasonable accuracy using a single calculation and further refined following an iterative procedure. The approach has particular utility for studies assessing material properties using a range of electronic-structure methods, and for studies employing methods with high computational cost (e.g., hybrid, meta-hybrid, and double-hybrid treatments of electron exchange and correlation). We validate the approach for four Zn and Pb chalcogenides and demonstrate its utility in describing the electronic structure and magnetic structure of one complex semiconductor (Cu₂ZnSnS₄) and one metal-organic framework (HKUST-1), respectively.

The goal of local crystal-structure optimisation is to minimise all degrees of freedom (cell size, shape, and positions) with respect to the total energy of the system. It is convenient to employ an EoS based on an energy-volume (E-V) curve, where the remaining degrees of freedom (i.e., shape and positions) are minimised for each volume using standard numerical minimisation approaches (e.g., the conjugate-gradient method). Kohn-Sham density functional theory (DFT)¹⁴ is one of the most widely used electronic structure techniques for modelling solid-state materials. Most DFT codes provide optimisation algorithms for this purpose, e.g., the ISIF=4 setting in the Vienna Ab initio Simulations Package (VASP),¹⁵ the cell_dofree=`shape' setting in Quantum-Espresso¹⁶ or the CVOLOPT setting in CRYSTAL.¹⁷ 

A superficial resemblance is clear between E-V curves obtained with different exchange-correlation functionals, with similar shapes but different minima (Figure 1(a)). The extent of the similarity becomes apparent when using pressure-volume (P-V) curves (Figure 1(b)), where “pressure” refers to the scalar hydrostatic pressure on the periodic system. As this pressure $P=− d E d V$⁠, the optimal geometries $d E d V =0$ are now those intersecting the P = 0 line. We note that while these still differ depending on the chosen XC functional, the P-V curves have similar curvature, with the same approximate slopes about their zero-crossing points. From these we make our key assumption: the slope of one P-V curve may be used to estimate the crossing point of another.

For certain beyond-DFT calculation methods, the stress tensor is not computed directly. However, where the energy is available the hydrostatic pressure may always be estimated with a finite difference,

The procedure, outlined in Figure 2, is as follows.

1. Form a P-V curve using one description of the interatomic interactions (method A). This can be achieved by fitting an EoS to an energy-volume curve. If a system-specific set of classical potentials is available, this would be expected to perform very well as they are often fitted to the experimental lattice parameters and elastic properties. Within DFT, descriptions of electron exchange and correlation within the generalised-gradient approximation (GGA) are suitable,¹⁸ given their low cost and the availability of analytical gradients for the rapid calculation of forces. Comparative studies suggest that PBEsol¹⁹ gives especially good estimates for the lattice parameters and elastic properties of crystals.^20,21

2. Calculate the slope about P = 0 for method A. This will form our linear approximation,

   $m = d P A d V | P A = 0 .$

   (2)

3. Perform a calculation using a second approach (method B), e.g., hybrid DFT with the screened HSE06 functional,²² using an estimated initial volume; this may be the equilibrium volume (V₀) for method A. Convert the resulting stress tensor to a hydrostatic pressure P₀. (If no analytical stress tensor is available, use a finite difference as in Eq. (1).)

4. Estimate the corrected volume for method B,

   $V 1 = V 0 + P 0 m .$

   (3)

5. Generate a unit cell with volume V₁ (e.g., by interpolating between the previous calculations with method A) and recalculate the desired properties with method B.

6. Iterate steps 4 and 5 until P is acceptably low,

   $V n + 1 = V n + P n m , P n = f ( V n ) .$

   (4)

In this approach, a linear fit is used for the pressure-volume relationship,

This is by no means a conventional equation of state but may provide a useful approximation when close to the minimum volume. The standard definition of the bulk modulus,

yields the static bulk modulus B₀ when evaluated at the equilibrium volume V₀,

As we have assumed this derivative to be constant, we combine with Eq. (6) to give a physically meaningful expression of our assumption,

It is now straightforward to compare this approximate EoS with a more conventional form for solid materials, estimating an associated error (ϵ). The simplest case is a system with constant bulk modulus,

At a typical volume deviation of 5% (Table I),

where the pressure

and hence the fractional error $ϵ P$ is −2.5%.

Moving to an improved, while still relatively simple, EoS, the Murnaghan EoS adds a parameter, effectively giving a linear volume dependence to B₀.²⁴ Taking its derivative form (Eq. (18)), we improve our error estimate,

We can relate B₀ to the Murnaghan parameters $k 0 , k 0 ′$ by finding the slope at V₀,

Plotting these error estimates in Figure 3, we find that the error in pressure is less than 10% of the static bulk modulus for volumes 10% from the optimum, given a material that obeys the Murnaghan EoS with a typical value $k 0 ′ =5$⁠. For smaller values of $k 0 ′$ (i.e., closer to the fixed-bulk-modulus model), the errors are greatly reduced. In any case, the linear approximation appears to be sufficiently accurate for a stable optimisation process.

Returning to the simplistic EoS of Eq. (10),

we examine the residual pressure P₁ following a single step of RVO from an initial volume V_i,

We note that ln(x) ≈ x − 1 for x close to 1, and hence the residual pressure is approximately linear with respect to the error in initial volume estimate. The term $B 0 , B V 0 , A B 0 , A V 0 , B$ indicates a smaller linear dependence on the similarity of the EoS for method A and method B. Moving to the Murnaghan EoS,

Again, the pressure is dominated by the initial position, with a smaller contribution from the difference between EoS “stiffness” and volume minima. The non-linearity of this relationship follows the non-linearity of the true EoS through the exponent $k 0 , B ′$⁠. We conclude therefore that the performance of the first step is equally sensitive to percentage differences in equilibrium volume and bulk modulus between method A and method B. Convergence is impacted by the non-linearity of the EoS but not by how accurately this non-linearity is reproduced by method A.

Studies have been carried out on the binary chalcogenides PbS, PbTe, ZnS, and ZnTe as well as the quaternary semiconductor Cu₂ZnSnS₄ and an organic-inorganic hybrid material HKUST-1.

All DFT calculations on the binary chalcogenides were carried out with VASP²⁵ using the two-atom primitive face-centred cubic unit cells. We employed projector-augmented wave (PAW) frozen-core potentials^26,27 treating the outermost s and p electrons of S, Te, and Pb and the outermost s, p, and d electrons of Zn explicitly as valence. For consistency, we used the LDA PAW potential set. The PAW potentials are highly transferable and tests showed a very weak dependence of the resulting optimised electronic and crystal structures. A plane-wave kinetic-energy cutoff of 550 eV was employed in all these calculations, and the Brillouin zone was sampled with an 8 × 8 × 8 Γ-centered Monkhorst-Pack mesh.²⁸ During electronic minimisation, the wavefunctions were optimised to an energy tolerance of 10⁻⁸ eV. These parameters were found to be sufficient to converge the absolute total energies to within 1 meV atom⁻¹, and the stress tensors to well within 1 kbar (0.1 GPa).

The simplicity of the binary systems allowed us to test a wide range of functionals, spanning different “levels” of approximations to the exchange-correlation potential.²⁹ As a baseline, we took the LDA with the Perdew-Zunger parameterisation of the correlation energy.³⁰ Calculations using the GGA were performed with the Perdew-Wang 91 (PW91)^18,31,32 and Perdew-Burke-Enzerhof (PBE)³³ functionals, plus the variant of PBE revised for solids (PBEsol).¹⁹ To complement this set of functionals, we also tested the Grimme D2 dispersion correction to PBE.³⁴ Meta-GGA calculations were carried out using the Tao-Perdew-Staroverov-Scuseria (TPSS) functional³⁵ and the subsequent revision of Perdew et al. (revTPSS).³⁶ Finally, we tested two hybrids, viz., the popular HSE06²² and B3LYP³⁷ functionals. For each material and functional, we calculated an E-V curve by adjusting the lattice parameter to yield 21 volumes about the experimental 300 K lattice parameters reported in Refs. 38 and 39 covering a range of approximately. ±5%. We note that, as a result of the high symmetry of these systems, the lattice parameter is the only degree of freedom, and thus relaxation of the cell shape and internal positions was not required.

For Cu₂ZnSnS₄ (Section V B), energy-volume curves were formed from all-electron DFT calculations using the FHI-aims code.^40,41 These calculations employed numerically tabulated atom-centered basis functions (the “tight” defaults were used, which correspond to expected convergence of <10 meV per atom) and evenly spaced k-point grids. Additional hybrid (HSE06) DFT calculations and primitive-cell optimisations used VASP with the PAW-PBE potential set and a 500 eV cutoff for the plane wave basis set. All calculations on Cu₂ZnSnS₄ sampled the Brillouin zone with 6 × 6 × 6 Γ-centered k-point grids.

For the Cu-based metal-organic framework HKUST-1, calculations were again performed with the VASP code, considering only the point Γ in reciprocal space due to the large size of unit cell. Owing to the presence of the open-shell Cu(II) ion (d⁹ configuration) all calculations were spin-polarised, and a range of magnetic structures were tested as discussed in Section V C. The PBEsol and HSEsol functionals were used along with the PAW-PBE potential set. Here “HSEsol” refers to a modification of HSE06, with PBEsol replacing PBE as the local exchange-correlation functional.⁴² Due to the complexity of the crystal structure, only three energy-volume points were included in the EoS and a single iteration of RVO was performed to recover the ground-state HSEsol structure. A slightly different procedure was followed in this case: a quadratic E-V curve was fitted to the three PBEsol points. The initial HSEsol calculation was carried out at the fully optimised PBEsol point, and the E-V curve was followed assuming a constant pressure offset to estimate the equilibrium volume for HSEsol. (This application was the first chronologically, and led to the development of the fitting procedure based on the Murnaghan EoS.) Calculation data and E-V curves are made available in the supplementary material.⁶⁸ 

The RVO method was implemented and tested with code written in Python 2.7.3, using the standard library and Numpy/Scipy/Matplotlib.^43–45 (The code is freely available; details in the supplementary material.⁶⁸) Non-linear fitting to the Murnaghan EoS uses the curve_fit routine in the Scipy Python library, which accesses Minpack, an open-source Fortran library.⁴⁴ This implements least-squares fitting with the Levenberg-Marquardt algorithm.⁴⁶ Initial guesses of 50.0 eV Å⁻³ and 5.0 were used for the k′ and $k 0 ′$ parameters, respectively.

For PbS there is a significant range in equilibrium lattice parameters for different exchange-correlation functionals, corresponding to a maximum volume difference of over 10%, between the LDA and B3LYP calculations (Figure 1). Values are tabulated in Table I and compared to a recent low-temperature study by Knight.²³ The computed values are similar, but slightly different, to the computational results reported by Hummer et al.⁴⁷ Direct bandgaps were also calculated at each volume expansion, at the special k-point X for the lead compounds and at Γ for the zinc compounds. (Note that PbS and PbTe have another, smaller direct bandgap at the L point.) It can be seen in Figure 4 that over the lattice-parameter expansion and contraction of up to 5%, the bandgaps vary by around 1 eV, with the direction of change depending on the structure type and chemistry. In this case, using a LDA-predicted geometry for a “single-shot” B3LYP calculation would lead to a difference in bandgap of ∼0.4 eV compared to that at the equilibrium geometry for B3LYP.

An iterative application of the RVO procedure was then simulated from the data. The Murnaghan EoS (Eq. (36)) in its integrated form (Eq. (37)) was fitted to each E-V curve from DFT calculations. This allowed energy and pressure to be calculated for arbitrary volumes without carrying out additional DFT calculations,

The quality of these fits was sufficient for this exercise, with RMS fitting errors of <1 meV. Fitting parameters and full data are included in ESI. For each “test functional”-“reference functional” pair, the minimum volume (corresponding to the fitting parameter V₀) of the reference functional was taken as the initial volume guess, and an external pressure calculation modelled by evaluating the pressure at this volume using the EoS for the test functional. This refined pressure and volume was then used as the basis for further iterations. The external pressure over successive iterations is shown for PbS in Figure 5 for each combination of functionals; convergence is rapid with the residual pressure dropping almost logarithmically with subsequent steps, typically by a factor of ∼10³ in three steps.

In general terms, a direct optimisation with method B will take N_opt,B steps, each requiring an average computing time t_B, to converge to the equilibrium volume. Constructing an EoS for RVO using method A requires N_EoS,A optimisations, which, as for method B, take N_opt,A steps of time t_A. We note that in most cases N_opt,B will be larger than N_opt,A, since the direct optimisation with method B must adjust the internal coordinates, cell shape, and volume, while the optimisation with method A needs only to optimise the internal coordinates and the shape. Subsequent application of N_{RV O} iterations of the algorithm then requires 1 + N_{RV O} single-point calculations using method B, each again requiring t_B time. RVO is expected to be more efficient than a direct optimisation with method B if the following inequality holds:

The cubic systems considered in this section, for which the cell volume is the only degree of freedom, represent a special case where N_opt,A = 1. We assume that energy gradients are available with method B and that the optimisation algorithm would converge in three steps, i.e., N_opt,B = 3. This is reasonable if a good estimate of the starting volume is available, such as a room-temperature lattice constant. The inequality simplifies to

it can be seen that RVO will outperform a direct optimisation if a suitable pressure is obtained after one iteration while $t A < t B N EoS , A$⁠.

As a concrete example, we compared a direct optimisation of PbS with HSE06 to an optimisation with RVO using PBEsol and HSE06 as method A and method B, respectively. The initial structure for both optimisations was the experimentally measured room temperature volume, and an 11-point EoS for RVO was computed about this value using PBEsol. The direct optimisation used a quasi-Newtonian algorithm as implemented by VASP (with the input tag "IBRION=1"). Both sets of calculations were performed on a dual-CPU Intel Xeon workstation with 12 physical cores and 64 Gb RAM, allowing the timings to be compared directly. The comparison is given in Table II.

In this test, a single-point calculation with HSE06 was on average 150 times more expensive than a calculation with PBEsol; this is both due to the higher complexity of non-local hybrid functionals compared to semi-local GGA methods, and to the different scaling properties with the number of k-points used to sample the Brillouin zone. With the force convergence criterion of 10⁻² eV Å⁻¹ for direct optimisation, the pressure was reduced to −0.1 kbar from an initial pressure of 3.42 kbar in three steps, taking 4.75 h on our test hardware. A single iteration of RVO yielded a pressure of 0.15 kbar in 4.18 h, while a second iteration yielded 0.03 kbar in a total time of 6 h.

It can be seen from the data in Table II that the direct optimisation takes on average less time per force calculation than RVO; the procedure implemented in VASP re-uses calculated wavefunctions to speed up the convergence of the second and third steps. In this case, we found that one of the conjugate-gradient electronic-minimisation cycles during the first single-point calculation for the RVO algorithm took some 1500 s longer than both the other steps in this series and the first step of the quasi-Newtonian volume optimisation, despite the latter being notionally an identical calculation. This artefact contributes significantly to the cost of the single-iteration RVO calculation.

Nonetheless, even for this relatively simple test case, useful savings in computing time could potentially be obtained in practice with RVO. Given the poor scaling of computational cost with system size when using advanced electronic-structure methods, we would expect more substantial savings for more complex unit cells. This would also be true in systems where direct optimisation requires the minimisation of a larger number of degrees of freedom, leading a larger number of steps with method B, which is the subject of the following case studies.

Cu₂ZnSnS₄ (CZTS) is an attractive light-absorbing material for thin-film photovoltaics, with a direct bandgap and consisting of earth-abundant components, which has attracted significant experimental^49–53 and computational^54–58 research effort. In the search for new materials for solar energy conversion, the prediction of accurate bandgaps from first-principles is a serious challenge and CZTS represents a suitable case for probing the effect of crystal structure.

An initial structure for CZTS in the kesterite phase, optimised with PBEsol, was obtained during previous work.⁵⁹ This was reduced from a conventional 16-atom body-centered-tetragonal cell with $I ̄ 4$ symmetry to the corresponding 8-atom primitive cell using Spglib.⁶⁰ A set of seven structures was obtained for both cells by multiplying each lattice vector by a scale factor from 0.97 to 1.03 and performing a local optimisation of the atomic positions, this forming the “method A” energy-volume curves. In addition to this isotropic scaling, an additional set of structures were calculated including optimisation of the cell shape (i.e., the tetragonal c/a ratio) for each volume point. The iterative RVO procedure was followed in order to minimise the pressure and obtain a more accurate electronic structure using the HSE06 functional. The results are given in Table III; pressure minimisation was rapid in all cases, decreasing logarithmically with each step.

We note that the resulting lattice parameters from these calculations, especially those using the primitive cell, are very close to both the experimental lattice parameter a = 5.427 Å and theoretical lattice parameter a = 5.448 Å reported by Paier et al. following a conventional optimisation procedure with a variant of the HSE functional.⁶¹ We also note a bandgap shift of almost 0.1 eV when the E-V curve was provided by a non-isotropic set of primitive lattices.

This case also highlights the importance of internal structure optimisation. After two steps of optimisation using the E-V curve further calculations were carried out, employing the HSE06 exchange-correlation functional, where the internal atomic positions were relaxed while fixing the unit cell. As shown in the table, these lead to an increase in the absolute pressure, but also a considerable improvement in the bandgap estimation compared to experimental measurements. Previous electronic structure studies have shown that the bandgap of CZTS is highly sensitive to the S positions, which correspond to deviations away from the ideal tetrahedral coordination environment.⁵⁵ In the ideal kesterite crystal structure, the metal nuclei all occupy high symmetry Wyckoff positions (2a and 2c by Cu, 2d by Zn, and 2b by Sn). However, the sulfur anions occupy the lower symmetry 8g positions, with x, y, and z displacement parameters. The change of ∼0.3 eV in the bandgap following further optimisation (Table III) emphasises the importance of internal relaxations for quantitative studies of electronic structure.

In 1999, Williams and co-workers isolated Cu₃(btc)₂ (HKUST-1).⁶² Since then this material has been widely studied in the field of metal-organic frameworks (MOFs), with possible applications in catalysis, ionic, and electrical conductivity, photovoltaics, batteries, and gas capture.⁶³ First-principles calculations of MOF properties have traditionally posed challenges for computational chemists because they combine large unit cells with complex organic and inorganic components.

HKUST-1 features an additional layer of complexity: it is composed of Cu-Cu “paddlewheel” inorganic regions, where each Cu(II) atom is associated with an unpaired electron and each paddlewheel is antiferromagnetically (AFM) coupled in the ground state configuration.^64,65 The magnetic interactions are highly sensitive to the Cu-Cu separation. Moreover, previous studies^66,67 have shown that the ionisation potentials and bandgaps of porous materials are sensitive to cell pressure and volume, similar to some of their inorganic counterparts. HKUST-1 represents an extreme case, where deviations from the equilibrium Cu-Cu distance result in spin flipping and formation of a ferromagnetic (FM) state, which impacts the electronic structure.

Typically, PBEsol-optimised structures agree with low-temperature experimental measurements of MOFs to within 1%. This is the case here, and PBEsol also reproduces the correct AFM state. However, a single-point HSE06 calculation on the PBEsol structure yields an incorrect FM ground-state, as shown in Fig. 6, with an associated HSE06 cell pressure of −13.98 kbar (Table IV). In order to recover an accurate HSE06 crystal structure (and the associated correct magnetic structure), a single iteration of RVO was required. Notably, there is not only a magnetic difference but also a significant difference in predicted electronic bandgap and workfunction.

The RVO approach presented here uses information from an inexpensive energy-volume curve to obtain a useful estimate of the optimal unit cell volume for a different level of theory. Our focus was in bridging between different exchange-correlation functionals within density functional theory, but a measured bulk modulus or classical interatomic potential could also be used to construct the reference energy-volume data. In sensitive systems the volume change can lead to qualitative differences in the electronic and magnetic properties. The results depend on the initial volume estimate and are relatively insensitive to the accuracy of the E-V curve. The RVO method is expected to be competitive with conventional optimisation approaches for simple symmetric unit cells, as demonstrated for rocksalt structured PbS. For materials such as Cu₂ZnSnS₄ that are sensitive to the atomic positions within the unit cell, RVO may form part of the optimisation approach but direct optimisation of internal positions is still needed. More significantly, it allows for the improved estimation of properties for large unit cells as demonstrated for HKUST-1, where conventional optimisation methods may be infeasible. As the only property used is the hydrostatic pressure, it is possible to employ calculation methods which return a total energy without analytical gradients by evaluating the energy of a single finite difference. In this case, an improvement over the “single-shot” may be obtained with two additional high-level calculations and an inexpensive E-V curve; a fourth high-level calculation would give an estimate of the convergence. We expect that in many cases this will prove an economical approach for the application of state-of-the-art electronic structure calculations in the solid state.

The authors thank J. M. Frost for useful discussions. We acknowledge our use of the ARCHER computing facility through our membership of the UK’s HPC Materials Chemistry Consortium, which is funded by EPSRC Grant No. EP/L000202. J.M.S. is funded by an EPSRC Programme Grant (No. EP/K004956/1). A.J.J. is funded by the EPSRC Doctoral Training Centre in Sustainable Chemical Technologies (Grant Nos. EP/G03768X/1 and EP/L016354/1). A.W. acknowledges support from the Royal Society and the ERC (Grant No. 277757).