Crystal structure optimisation using an auxiliary equation of state

Standard procedures for local crystal-structure optimisation involve numerous energy and force calculations. It is common to calculate an energy-volume curve, fitting an equation of state around the equilibrium cell volume. This is a computationally intensive process, in particular for low-symmetry crystal structures where each isochoric optimisation involves energy minimisation over many degrees of freedom. Such procedures can be prohibitive for non-local exchange-correlation functionals or other 'beyond' density functional theory electronic structure techniques, particularly where analytical gradients are not available. We present a simple approach for efficient optimisation of crystal structures based on a known equation of state. The equilibrium volume can be predicted from one single-point calculation, and refined with successive calculations if required. The approach is validated for PbS, PbTe, ZnS and ZnTe using nine density functionals, and applied to the quaternary semiconductor Cu$_{2}$ZnSnS$_{4}$ and the magnetic metal-organic framework HKUST-1.


I. INTRODUCTION
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) are calculated. 1 The optimisation of a crystal structure may require hundreds of self-consistent field iterations across a series of ionic configurations. 2 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. 3 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. 4 It has become commonplace to use an 'equilibrium' crystal geometry, determined using one exchangecorrelation 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 a) a.walsh@bath.ac.uk properties and defect formation energies. [5][6][7][8][9][10][11][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). 13 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 exchangecorrelation (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 electronicstructure 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 and magnetic structure of one complex semiconductor (Cu 2 ZnSnS 4 ) and one metal-organic framework (HKUST-1), respectively.

II. OUTLINE OF PROCEDURE
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) 14 is one of the most widely arXiv:1507.08815v2 [cond-mat.mtrl-sci] 12 Oct 2015 used electronic structure techniques for modelling solidstate materials. Most DFT codes provide optimisation algorithms for this purpose, e.g. the ISIF=4 setting in the Vienna Ab initio Simulations Package (VASP) 15 , the cell dofree='volume' setting in Quantum-Espresso 16 or the CVOLOPT setting in CRYSTAL. 17 A superficial resemblance is clear between E-V curves obtained with different exchange-correlation functionals, with similar shapes but different minima (Figure 1a). The extent of the similarity becomes apparent when using pressure-volume (P-V) curves (Figure 1b), where "pressure" refers to the scalar hydrostatic pressure on the periodic system. As this pressure P = − dE dV , the optimal geometries dE dV = 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: 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, 18 given their low cost and the availability of analytical gradients for the rapid calculation of forces. Comparative studies suggest that PBEsol 19 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.
3. Perform a calculation using a second approach (method B), e.g. hybrid DFT with the screened HSE06 functional 22 , using an estimated initial volume; this may be the equilibrium volume (V 0 ) for method A. Convert the resulting stress tensor to a hydrostatic pressure P 0 . (If no analytical stress tensor is available, use a finite difference as in Eqn. 1.) 5. Generate a unit cell with volume V 1 (e.g. by interpolating between the previous calculations with method A), and recalculate the desired properties with method B.

A. Accuracy of linear approximation
In this approach, a linear fit is used for the pressurevolume 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 0 when evaluated at the equilibrium volume V 0 .
As we have assumed this derivative to be constant, we combine with Eqn. (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 0 . 24 Taking its derivative form (Eqn. 18), we improve our error estimate: We can relate B 0 to the Murnaghan parameters k 0 , k 0 by finding the slope at V 0 :  We can relate B 0 to the Murnaghan parameters k 0 , k 0 0 by finding the slope at V 0 : 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 0 = 5. For smaller values of k 0 0 (i.e. closer to the fixed-bulk-modulus model), the errors are greatly reduced. In any case, the linear approximation appears to be su ciently accurate for a stable optimisation process.
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.

B. Dependence on accuracy of EoS
Returning to the simplistic EoS of Eqn. 10 we examine the residual pressure P 1 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 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.

A. Electronic Structure Calculations
Studies have been carried out on the binary chalcogenides PbS, PbTe, ZnS, and ZnTe as well as the quaternary semiconductor Cu 2 ZnSnS 4 and an organic-inorganic hybrid material HKUST-1.
All DFT calculations on the binary chalcogenides were carried out with VASP 25 using the two-atom primitive face-centred cubic unit cells. We employed projectoraugmented 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. 28 During electronic minimisation, the wavefunctions were optimised to an energy tolerance of 10 −8 eV. These parameters were found to be sufficient to converge the absolute total energies to within 1 meV atom −1 , 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. 29 As a baseline, we took the local-density approximation (LDA) with the Perdew-Zunger parameterisation of the correlation energy. 30 Calculations using the generalisedgradient (GGA) approximation were performed with the Perdew-Wang 91 (PW91) 31-33 and Perdew-Burke-Enzerhof (PBE) 34 functionals, plus the variant of PBE revised for solids (PBEsol). 19 To complement this set of functionals, we also tested the Grimme D2 dispersion correction to PBE. 35 "Meta-GGA" calculations were carried out using the Tao-Perdew-Staroverov-Scuseria (TPSS) functional 36 and the subsequent revision of Perdew et al. (revTPSS). 37 Finally, we tested two hybrids, viz. the popular HSE06 22 and B3LYP 38 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. 39 and 40 covering a range of approx. ±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 2 ZnSnS 4 (Section V B), energy-volume curves were formed from all-electron DFT calculations using the FHI-aims code. 41,42 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 2 ZnSnS 4 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 9 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. 43 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 Murnaghan fitting procedure.)

B. Implementation
The RVO method was implemented and tested with code written in Python 2.7.3, using the standard library and Numpy/Scipy/Matplotlib. 44

Simulated application across a range of methods
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 K. S. Knight. 23 The computed values are similar, but slightly different, to the computational results reported by Hummer et al. 48 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 an 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 (Eqn. 37) in its integrated form (Eqn. 38) 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 0 ) 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 3 in three steps.

Comparison with a standard optimisation procedure
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 RVO iterations of the algorithm then requires 1+N RVO singlepoint 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 FIG. 5: Residual pressures following six iterations of volume optimisation of PbS using different "reference"/"test" combinations of exchange-correlation functionals. The "reference" refers to the functional used for the EoS ("method A"), while the groups of bars correspond to the functionals used in simulated single-point calculations ("method B"). Note that the pressure is presented on a logarithmic scale.
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 eleven-point EoS for RVO was computed about this value using PBEsol. The direct optimisation used a quasi-Newton 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 −2 eV A −1 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 hrs on our test hardware. A single iteration of RVO yielded a pressure of 0.15 kbar in 4.18 hrs, while a second iteration yielded 0.03 kbar in a total time of 6 hrs.
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 electronicminimisation cycles during the first single-point calculation for the RVO algorithm took some 1500s longer than both the other steps in this series and the first step of the quasi-Newton volume optimisation, despite the latter being notionally an identical calculation. This artefact contributes significantly to the cost of the single-iteration RVO calculation. II: Comparison of a direct HSE06 volume optimisation and one and two iterations of the RVO algorithm in determining the equilibrium volume of PbS. For each step, the total time for each step is recorded alongside the cell volume and pressure after the cycle where appropriate. For the direct optimisation, the timings of the three steps are printed alongside the total for the complete calculation, so the latter includes additional operations such as setup time and is slightly longer than the sum of the three electronic minimisations.

Algorithm
Step 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.

B. Quaternary Sulfide Cu 2 ZnSnS 4
Cu 2 ZnSnS 4 (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 50-54 and computational 55-59 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. 60 This was reduced from a conventional 16-atom body-centered-tetragonal cell withĪ4 symmetry to the corresponding 8-atom primitive cell using Spglib. 61 A set of seven structures was obtained for both cells by multiplying each lattice vector by a scale factor from 0.97-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 a=5.448Å reported by Paier et al. following a conventional optimisation procedure with a variant of the HSE functional. 62 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 show that the bandgap of CZTS is highly sensitive to the S positions, which correspond to deviations away from the ideal tetrahedral coordination environment. 56 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)  In 1999, Williams and co-workers isolated Cu 3 (btc) 2 (HKUST-1). 63 Since then this material has been widely studied in the field metal-organic frameworks (MOFs), with possible applications in catalysis, ionic and electrical conductivity, photovoltaics, batteries and gas capture. 64 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. 65,66 The magnetic interactions are highly sensitive to the Cu-Cu separation. Moreover, previous studies 67, 68 have III: Results from application of RVO to Cu 2 ZnSnS 4 , using the HSE06 functional and a PBEsol-derived E-V curve. The unit of pressure P is kbar (10 8 Pa) and volumes V are given inÅ 3 . a is the lattice parameter inÅ; these are calculated as a mean over the a and b vectors (crystal symmetry is not enforced in FHI-AIMS). E g is the electronic bandgap in eV taken from the Kohn-Sham eigenvalues at the Γ-point. The methods in parentheses refer to the process by which the E-V curve was generated; isotropic expansion energies with FHI-aims and volume-conserving relaxations with VASP. Iteration "2*" is the data from a final set of calculations, where the internal atomic positions are relaxed while maintaining the unit cell from iteration 2.  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 lowtemperature 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.

VI. CONCLUSIONS
The rapid volume optimisation (RVO) approach presented here uses information from an inexpensive energyvolume 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 energyvolume 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 2 ZnSnS 4 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.

VII. SUPPLEMENTARY DATA
Python code providing a reference implementation of this method is available at http://github.com/ wmd-bath/rvo; a snapshot as of this publication is available at DOI:10.5281/zenodo.31940. This includes a program to generate the plots in this publication from the binary chalcogenide data. Calculation data has been deposited online with Figshare at the DOI: 10.6084/m9.figshare.1468388. Raw output files from the hybrid DFT calculations are made available for CZTS and HKUST-1, and energy-volume curves are available for all the systems. The full set of graphs and fitting parameters for PbS, PbTe, ZnS, and ZnTe are also included as supplementary data with this paper.