# Deriving ab initio model Hamiltonians for molecular crystals

## Abstract

Developing realistic and precise models of the electronic properties of organic molecular crystals is crucial for understanding the full range of strongly correlated phases that they exhibit. By using ab initio model construction methods, one can obtain unbiased non-interacting models of such systems from density functional theory, upon which one can base further (many-body) models. We will discuss the utility and advantages of ab initio model construction using Wannier orbitals. We will briefly review the approach, and then explain why it is so well suited to molecular crystals in particular. We discuss the ab initio construction of both non-interacting and interacting Hamitonians, and highlight recent examples where such first principles models lead to importantly different results than fitted models.

###### Contents:

## I Background

### i.1 Molecular Crystals

Molecular crystals, ordered periodic arrays of molecules, are known to exhibit a wide range of quantum mechanical phenomena, including unconventional superconductivity, quantum criticality, frustrated anti-ferromagnetism, and quantum spin liquid behaviour Chaikin et al. (1998); Powell and McKenzie (2011); Jacko et al. (2013a); Dressel (2007); Seo et al. (2004); Brown (2015); Jacko et al. (2013b); Seo et al. (2015). In some cases many of these phases can be found in a single material by tuning an external parameter, such as pressure or magnetic field. Often, one can control which phase is expressed by making subtle physical or chemical changes to the molecules Seo et al. (2004); Dressel (2007); Jacko et al. (2013b); Seo et al. (2015); Brown (2015). Along with the flexibility of the interactions within individual molecules, molecular crystals also have a range of intermolecular interactions. It is the subtle competition between the many intra- and inter-molecular interaction energies that brings the wide variety of phases seen in experiments so close together. These crystals tend to have a low effective dimension, and this likely contributes to the close competition between the phases Seo et al. (2004); Dressel (2007); Brown (2015); Coldea et al. (2010).

Molecular crystals are an exciting testing ground for finding and understanding new emergent states of matter. They often display competition between multiple emergent strongly correlated ground states Dressel (2007); Powell and McKenzie (2011). This, and the flexibility of organic chemistry, means that the emergent physics is often tuneable by subtle chemical and physical modifications (making slight variations around a core motif) Seo et al. (2004); Dressel (2007); Jacko et al. (2013b); Seo et al. (2015); Brown (2015). In one notable case, a superconducting state in an organic molecular crystal is destroyed by substituting some hydrogen atoms for deuterium, its heavy isotope Taniguchi et al. (1999). That this extremely subtle change can have such profound consequences is both exciting and intimidating. On one hand, it presents the inviting prospect of creating strongly correlated materials with technologically desireable properties; on the other, the level of detail required to correctly predict the phase of a material can be substatial.

Fig. 1 shows the phase diagram for a family of organic crystals called Fabre salts [salts of TMTTF (shown in Fig. 2) and an anion], which can be tuned through many different phases by applying physical pressure, or by slightly changing their anions, often thought of as a âchemical pressureâ. As the size of the anion decreases, the TMTTF molecules pack closer together, as they would under the application of physical pressure. This range of accessible phases indicates that the many competing interactions are very close in energy. The many competing energy scales in molecular crystals gives us access to phases with various interesting physical properties.

### i.2 Predictive versus postdictive modeling

There is an important philosophical point to be made here. The goal of science should be to make predictions about the nature of nature, and then to test those predictions. In the case of building models of materials, what is typically applied is a postdictive approach; one ‘knows’ that to have a model with the observed behaviour, it should be this lattice with that type of interaction (e.g. Heisenberg model on a triangular lattice, extended Hubbard model on a square lattice, and so on). Thus, there is limited information gained about the system, whether one finds the behaviour one was searching for in such a model or not. One one hand, one chose the model phenomenologically, so finding the correct phenomenology is not profound. On the other, not finding the expected behaviour could be due to any number of reasons from the profound to the trite.

Postdictive approaches can be useful, but one must go beyond them to gain a deeper understanding of the important commonalities and differences within a class of materials. For example, such an approach does not show much promise for describing all of the multitude of phases seen in the Fabre salts. What one ideally would like is a systematic way of constructing an effective many-body Hamiltonian from first principles. Constructing the non-interacting part from first principles is currently emminantly possible: By producing localised ‘Wannier’ orbitals (discussed in more detail later), one can use the results of density functional theory (DFT) to construct a tight-binding lattice without first assuming its form. Here we give a brief overview of the theoretical development of the concept and practical details of using Wannier orbitals. For a much more complete and mathematical account, turn to the excellent review of Marzari et al. Marzari et al. (2012).

On a related note, it is worth commenting on the distinction between ’ab initio’ and ‘first principles’; ab initio implies no empirical input, just calculations on the grounds of the many-electron Schrödinger equation using the fundamental constants of nature such as Planck’s constant, the charge of the electron, etc. On the other hand, first principles allows for empirical parameters. Both density functional theory and Wannier orbital construction are in principle ab initio, however particular implementations tend to include empirical parameterisations (specific density functionals, for example) that are properly considered first principles rather than ab initio.

### i.3 Development of Wannier Orbitals

In 1937 Gregory Wannier introduced the idea of constructing localised sets of wavefunctions by fourier transforming Bloch states Wannier (1937). For a Bloch wavefunction for band , , the corresponding Wannier function for band is

(1) |

where is any combination of the crystal lattice vectors with integer prefactors, is localised in the unit cell located at , and the integral runs over the first Brillouin zone (FBZ). These new wavefunctions have the advantages of atomic orbitals (such as locality) while enforcing orthogonality. This allows one to treat localised excitations of individual electrons in metallic materials on the same footing as the ‘bulk’ electrons (the delocalised Bloch states).

By the 1950’s these wavefunctions were widely known as Wannier functions, and of great use in understanding the physics of excitations in crystals. In 1953, George Koster introduced two new methods for defining Wannier functions without first having to solve the Schrödinger equation, and allowing one to use these orbitals to compute energy bands in crystals Koster (1953).

Walter Kohn put Wannier functions on a rigourous analytical grounding in 1959, showing that one can always find a unique, real, symmetry-preserving, and exponentially localised functions for a given single band Kohn (1959). In this work he showed (although not in so many words) that Wannier orbitals were the ideal basis for the recently-developed tight-binding method (closely related to the Hückel method used in chemistry) Hückel (1931); Slater and Koster (1954). (Kohn continued working on Wannier orbitals, and in the mid-90’s used the locality of Wannier functions, and the consequence that their interactions should decay exponentially, to propose a density functional theory method that scales linearly with the number of atoms Kohn (1993, 1995).) Jacques Des Cloizeaux further expanded the mathematical grounding of Wannier functions, and identified what would later become known as the disentangling problem: if bands overlap, it is difficult to construct Wannier functions for just one of those bands (requiring one to ‘disentangle’ the target band from the other bands it crosses) Cloizeaux (1964). This remains a general challenge in using Wannier functions to this day Marzari et al. (2012). Due to the practical difficulty of the disentangling problem, and the extra indeterminancy introduced in the disentangling proceedure, Wanner functions were not of signficant help in computational electronic structure theory until the 90’s, when approaches based on density functional theory were introduced Marzari and Vanderbilt (1997).

### i.4 Wannier orbitals in Density Functional Theory

The key breakthrough in the application of Wannier orbitals occured when Nicola Marzari and David Vanderbilt formulated a generalised approach for generating maximally localised Wannier functions for the case of multiple bands Marzari and Vanderbilt (1997). Not only that, they also described a numerical algorithm to produce such orbitals based on Bloch functions sampled on a mesh of points in -space, such as would be the output of a typical DFT code. This allowed for computations of Wannier orbitals in realistic, non-trivial cases. They also suggested the approach of using Wannier functions to construct effective model Hamiltonians for strongly correlated electron systems.

A few years later, Marzari and Vanderbilt, along with Ivo Souza, extended their original approach to allow for entangled bands. Together they introduced an efficient disentangling methodology Souza et al. (2001) that requires no additional information over a usual Wannier construction, just one additional assumption. That assumption is that the ‘character’ of the Wannier orbitals (the contributions from particular basis functions) should vary as smoothly and slowly as possible; this is enforced via minimising the change in character across the Brillouin zone.

These works laid the foundation for the wide-spread computation of Wannier orbitals in DFT codes. In 2008, the code wannier90 was released Mostofi et al. (2008). Developed by Arash Mostofi, Jonathan Yates, Young-Su Lee, along with Souza, Vanderbilt and Marzari, this code is now widely used, designed to interface with any DFT code to produce Wannier orbitals. It is now used in Wannier orbital construction in FPLO Koepernik and Eschrig (1999), WIEN2k Kuneš et al. (2010), Quantum ESPRESSO Giannozzi et al. (2009), ABINIT Gonze et al. (2009), and Fleur Freimuth et al. (2008), to list just a few of the more popular DFT codes for crystals.

## Ii Separation of energy scales in Molecular Crystals

Here we discuss the separation of energy scales that commonly occurs in molecular crystals and how this aids the construction of a minimal set of Wannier orbitals. Molecular crystals tend to have a separation of energy scales in their non-interacting states, while there is competition amongst many possible strongly correlated ground states in the full many-body treatment. Despite the advances made in disentangling procedures, it remains a highly challenging task to produce high quality, reliable Wannier orbitals from entangled bands. This is particularly important if one is concerned about capturing the fine features of the electronic structure, which can have significant effects on the many-body state, as I will discuss explicitly in the case of crystals based on Pd(dmit). Molecular crystals provide an exciting playground for applying Wannier orbital based techniques to their greatest potential. because one can bypass the difficulty and ambiguity of disentangling procedures, one can determine the significance of the fine features of the electronic structure in determining the rich phase diagrams of these materials.

The separation of energy scales in molecular crystals is straight-forward to understand: the molecules are held together by (strong) covalent bonds, while the crystal is held together by much weaker intermolecular forces; van der Waals, -stacking, and hydrogen bonding. The strong forces within a molecule produce a set of well spaced molecular orbitals (MOs), and these orbitals weakly couple between molecules, as illustrated in Fig. 2, producing bands that are narrow on the scale of the MO energy gaps.

We can understand this more concretely by considering a toy example: a 1D, two orbital tight-binding model. For simplicity we consider a chain of spinless fermions, with orbitals and on each site , with nearest neighbour interactions. The Hamiltonian is

(2) |

where is the energy difference between the orbitals, and is the inter-site hopping for orbital . In the case of molecular crystals, the orbitals are molecular orbitals of single molecules. The energy difference between the molecular orbitals () comes from the inter-atomic hopping within a molecule, typically a -type overlap. A typical energy scale for this difference between molecular orbitals in an organic molecule is a few eV (see for example Hinze and Beveridge (1971)). The inter-site hopping comes from the overlap of molecular orbitals on different molecules, and is exponentially suppresed by distance. These energy scales are on the order of 10 - 100 meV for nearest neighbour overlaps (see for example Jacko et al. (2013b)). Thus it is often the case in such systems that ; the bands resulting from each molecular orbital are narrow enough and well-separated enough that they do not overlap in energy. Thus, depending on filling, one can consider just one orbital or the other as the foundation for an effective model Hamiltonian.

In applying the Wannier construction procedure outlined above, we have glossed over the details of limiting the Fourier transform window to some small energy range. In the simplest case of a single band system, it is clear that this originally-infinite window can be truncated to be exactly the bandwidth of the single band without any loss of generality. However, in multi-band systems, extracting a subset of bands can become difficult. To understand why, let us consider again the 1D, two orbital model.

For , there are gapless excitations possible between these bands; the two bands have weight in an overlapping energy range. Thus, even in this simplest case without band crossings or interactions, the presence of a second band makes picking out the states of the first band non-trivial. When these bands cross or hybridise, this further complicates the procedure. This problem is very difficult to solve in general and is known as the disentangling problem.

Now, the separation of energy scales in molecular crystals comes in to play. As discussed above, because of the often quite different energy scales of inter- and intra-molecular interactions, it is typical to find well isolated sets of bands in the band structure of a molecular crystal, as illustrated in Fig. 3. This property means that one can bypass the difficulty and ambiguity of projective disentangling procedures. Thus minimal input is needed into the WO construction procedure in molecular crystals; one just inputs how many orbitals you would like, spanning what energy window.

## Iii ab initio Model Construction

It is important when modeling these systems that the models we use be as accurate and unbiased as possible; starting with preconceptions of how the model ’shouldâ look can limit what one finds.

### iii.1 Tight-Binding Models

DFT can give us useful information about the non-interacting electronic properties of a system. To utilise this information, we will construct a tight-binding model from the DFT using a rigorous Wannier orbital construction technique. This procedure creates localised orbitals that accurately represent the electronic properties found by DFT for the frontier electrons, those that determine the low-temperature physics. As discussed above, the separation of energy scales makes Wannier orbital construction straightforward in molecular crystals. Once one has local Wannier orbitals one can construct a first principles tight-binding model.

#### First principles versus Fitting

One might ask why is all of this effort justified? Why not simply write down a perfectly good tight binding model and fit it to a first principles band structure? (As is often done, see Kandpal et al. (2009); Scriven and Powell (2012); Jeschke et al. (2012); Seo et al. (2015) for just a sample.) For very simple systems, where the relevant tight-binding model is clear, and only has a few parameters, a first principles method is probably not justified as the band structure is well reproduced by fitting methods Kandpal et al. (2009); Scriven and Powell (2012); Jeschke et al. (2012); Seo et al. (2015). However, in the much more common situation of a somewhat ambiguous tight-binding model with an unknown number of relevant parameters, a first principles approach is ideal, as I will discuss further below. In fact, in some systems where a quite simple model seemed obvious, it has been shown that a somewhat more nuanced model does a better job of capturing the important many-body physics of the material (as discussed in the context of Pd(dmit), below).

John von Neumann is famously claimed to have said “With four parameters I can fit an elephant, and with five I can make him wiggle his trunk.” Dyson (2004) (this is more-or-less true, as shown in Fig. 4, after Ref. Mayer et al. (2010)); this captures the essence of the problem of fitting energy bands to a tight-binding model. With enough parameters in the fit, it is difficult to have a *bad* fit of the band structure; a model with many parameters might reproduce the dispersion without having any connection to a realistic microscopic description of the system. As such, a good fit provides almost no information, especially not about the quality of the microscopic model to which you are fitting. Worse, when there are many parameters, very different values can produce similarly good fits by whatever optimisation metric you are using. Often, these different sets of values have importantly different physical consequences (for example changing the electronic dimerisation of a chain, or the localisation of charge). These differences can lead to importantly different many-body ground states, as will be discussed further below.

To demonstrate this, I will discuss a particular case of producing a tight-binding model for an organic molecular crystal by fitting and via Wannier orbitals. TMTTFAsF is one of the Fabre salts, a family of organic charge-transfer salts with a rich phase diagram, given in Fig. 1. These flat organic molecules -stack into one dimensional chains along the crystallographic direction (shown in Fig. 5), with some inter-chain coupling in the plane, and minimal coupling in the direction (where the anionic AsF layers introduce a large spacing between TMTTF layers). Thus this system is largely 1D electronically, with some 2D nature introduced by next-nearest neighbour hopping and further terms Nogami et al. (2005); Jacko et al. (2013b); Brown (2015).

In this family of materials, once one decides to include any 2D hopping terms, one encounters the problem that there are many terms of similar magnitude; to avoid neglecting terms of the order one keeps, one needs to add many parameters to the tight-binding model. These many degrees of freedom cause problems for fitting procedures; one finds many local minima with similar optimiziation functions values, but very different parameters, with important physical consequences. This is the heart of the problem; that such fits are examples of sloppy modelsTranstrum et al. (2015): changes in one parameter can be almost completely masked by compensatory changes in other parameters.

Figure 6 shows the results of 5000 runs of a fitting procedure applied to the band structure of TMTTFAsF jes (); Jacko et al. (2013b); a pre-defined tight binding lattice is input, and the values of the 8 different hopping integrals (c.f. Fig. 8) are optimised to provide the best fit to the band structure (quantified by the least-squares error over the set of points the bands are sampled on). Each run starts with a (different) randomised set of fitting parameters and iterated on. Due to the nature of this method, each of the 5000 fits is basically unique. Rounding each to the nearest 0.1 meV, there are 26 unique fits; of the 5000 results, 74 have the minimum objective function. A given run has a 1.5% chance (1/67) of finding this ‘best’ solution! There are many more fits with a just slightly higher value of the objective function, and these fits contain contradictory physical information. For example, these different fits make different predictions about the electronic dimerisation of the system; whether the electronic dimers are on the structural dimers or not. This is an important difference, and can change as a function of pressure for example Jacko et al. (2013b). While most of the parameters have positive and negative equivalents, a few parameters are very precisely determined, and with a fixed phase. This means that the relative phases of the hopping integrals cannot be absorbed by a gauge transformation; these different sets of parameters have different physical meanings. The number of parameters used in the fit also effects the outcome. Removing some parameters from the inputed model will cause the other values to change. Damningly, in this system the optimal tight-binding fit is quite different to the set of hopping integrals found from Wannier orbitals, as shown in Fig. 6. The Wannier and fitted parameters make contraditary predictions about the electronic dimerisation of the system (as seen in the magnitude of the first two ’s). In the limit of including all Wannier overlaps out to infinite distance, the set of Wannier tight-binding parameters will exactly reproduce the band structure when the bands are not entagled.

Overall, it is hard to take fitting procedures too seriously with more than a couple of parameters in the fit. With small numbers of parameters, fitting becomes more stable, and in certain systems a few parameters is enough to acurately describe the band structure Kandpal et al. (2009); Jeschke et al. (2012); Seo et al. (2015). Even then, a fitted model should be considered an effective model that has potentially lost important detail in ‘integrating out’ the full set of parameters. In a sense these are ‘variational Hamiltonians’; they are optimised by some metric, but there is no assurance that they represent the underlying microscopic physics.

By producing Wannier orbitals for molecular crystals, and computing a set of ’s from those, one finds a single set of parameters that is reliable and robust; one can believe them just as much as one believes the other results of the DFT computation. In the general case of Wannier orbital construction, there is ambiguity involved in disentangling bands to produce the desired set of Wanniers Marzari et al. (2012). When the bands one cares about are well-separated from the bulk, this ambiguity is gone. Not only that, but one can gain knowledge by looking at the Wannier orbitals themselves. In the case of TMTTFAsF, the Wannier orbitals are localised to single TMTTF molecules, as shown in Fig. 7. Here, the Wannier orbital is very much like the HOMO (highest occupied molecular orbital) of a single TMTTF molecule in vacuum. Having this real-space orbital allows us to do many further computations based on the DFT, as well as producing a single robust parameter set for a tight-binding model (Fig. 8). This Wannier based model construction technique is being used more and more in molecular systems Nakamura et al. (2009, 2012); Jacko et al. (2013b, a); Altmeyer et al. (2015).

### iii.2 Including Interactions

Correctly parameterising many-body effects is an important and challenging task; it is the competition between energy scales that leads to the interesting physics in many systems, so small relative changes in large parameters can have large effects Jacko et al. (2010). Often, these parameters are estimated without careful consideration of the assumptions involved. For instance, if one considers a Hubbard model on a dimer with an inter-monomer hopping , on-monomer Hubbard repulsion , and inter-monomer Hubbard repulsion ; in the limit , , the effective Hubbard repulsion in the dimer orbitals is Komatsu et al. (1996); McKenzie (1998). This assumption is often used, since it allows one to estimate many-body parameters from straightforward band structure calculations, or molecular Hückel calculations. However, it is not well justified. It has since been shown that in the more realistic case of , that Scriven and Powell (2009). Thus, one still needs to be able to correctly compute many-body parameters to estimate the dimer parameters (even before considering screening). None-the-less, this approximation continues to be used (for example Jeschke et al. (2012); Yoshimi et al. (2012); Tsumuraya et al. (2013)), often without stating the strong underlying assumptions.

One might think that, given these real-space Wannier orbitals for a particular system, it must be straightforward to calculate the many-body Coulomb integrals and parameterise a Hubbard model. However, computing these terms by simply evaluating the Coloumb energy for each orbital neglects screening (equivalently, relaxation of the bulk states). Screening can easily suppress the Hubbard by an order of magnitude Nakamura et al. (2006).

In classical electromagnetics, there are many techniques for determining the response of a bulk to a perturbing field (analogous to the case here of computing the screening/relaxation of a doubly occupied orbital). The discrete dipole approximations (also called the coupled dipole approximation) is one such technique. In this approximation, one discretises the bulk as a set of polarisable dipoles, and self-consistently solves their response to the perturbing field and to each other DeVoe (1964). This method is very much like a technique applied to molecular crystals to compute screened Coloumb parameters. By representing each molecule by a set of polarisable (classical) dipoles, and placing a perturbing charge on one lattice site, one can compute the correction to the Hubbard repulsion due to the polarisation of the rest of the molecules in the crystal Cano-Cortés et al. (2010). This technique, though promising, has only been applied to a single molecular crystal (TTF-TCNQ), with no new applications apparent in the 5 years since the original publication.

An alternative approach to computing screened Coulomb parameters from first principles has gained prominence recently, the constrained random phase approximation (cRPA) Aryasetiawan et al. (2004). This technique is also based around computing the polarisation of the system, but in this case, the quantum mechanical polarisation function in the random phase approximation. Here we discuss RPA, cRPA, and its application to molecular crystals in more detail. Practically, it also relies on having Wannier orbitals for a few relevant bands, and so like the tight-binding model construction it is particularly suitable to molecular crystals.

#### Constrained Random Phase Approximation

The random phase approximation (RPA) was introduced by David Bohm and David Pines in the 1950’s to include the effects of screening into models of electron gases Bohm and Pines (1951); Pines and Bohm (1952); Bohm and Pines (1953); Pines (1953). Murray Gell-Mann and Keith Brueckner placed this approximation on a firmer footing, showing that the RPA can be derived from a self-consistent series of leading order Feynman diagrams Gell-Mann and Brueckner (1957), an example of which is illustrated in Fig. 9.

Given a basis of occupied and unoccupied states, one can compute an RPA polarisation function

(3) |

where () runs over the occupied (unoccupied) single particle states. With this polarisation function one can compute the screening effects of this system. This is exactly what one wants if computing the effects of an impurity in such a system, for example. However, if one wants to include many-body effects in a lattice model, then this constitutes an overcounting of the effect.

In 2004, Aryasetiawan et al. introduced a new, precise method for constructing sets of screened effective model parameters for strongly correlated lattice models Aryasetiawan et al. (2004). The constrained random phase approximation (constrained RPA or cRPA) is a systematic way of accounting for screening in the many-body parameters computed for some basis orbitals. The system is divided into two fragments; the active subspace (often labelled ‘d’), the space spanned by the orbitals of interest; and the rest of the bands (labelled ‘r’). On a conceptual level, this procedure computes the effects of transitions involving the ‘r’ subspace with RPA (by computing the polarisation function due to these transitions), while leaving the transitions within the ‘d’ subspace to be dealt with in the many-body model that results Aryasetiawan et al. (2004). This proceedure allows one to generate all the terms resulting from the Coulomb interaction, on- and off-site repulsive and magnetic interactions. Practically, one constrains the sums in Eq. 3 to exclude transitions within the active subspace.

The partitioning idea at the core of cRPA works best in the same situation that the Wannier orbital proceedure itself works best: a set of relevant bands well separated from the bulk. In the situation of entangled bands; where the natural basis one would like to use mixes with bands due to other states; one can apply the Wannier disentangling proceedure to construct a disentangled basis Miyake et al. (2009). In inorganic system, there are difficulties in disentangling the target bands from the bulk. Nonetheless this approach was quickly applied to transition metal systems and simple transition metal oxides Nakamura et al. (2006); Miyake and Aryasetiawan (2008); Nakamura et al. (2008); Miyake et al. (2009).

This approach has been applied to only a few organic crystals Nakamura et al. (2009, 2012). In those cases where it has, it finds sometimes importantly different parameter values. In the ET charge transfer salt -(ET) Cu (CN), cRPA predicts a value of for the dimer about twice as large as was estimated from a Hückel analysis of a dimer (using an optial conductivity estimate of the monomer value, ), vs Komatsu et al. (1996); Nakamura et al. (2009). In a simple Hubbard model, this would place this material well into the insulating phase, contrary to the observed metallic behaviour. The cRPA analysis also showed that off-site terms are significant, , meaning that to properly understand the system one must consider an extended Hubbard model Nakamura et al. (2009). While the optical conductivity estimate for the monomer is quite reliable, the assumptions in using this to estimate are not. This Wannier-based approach gives us a reliable first principles estimate of all the Hamiltonian parameters on the same footing.

## Iv First priciples approach finds important differences

Here we discuss particular examples to show that using a first principles approach can give importantly different results and insights than a standard fitting approach; be it caused by subtleties of parameter variations or qualitatively different lattices.

### iv.1 EtMeSb[Pd(dmit)]: Fine details matter

To demonstrate the importance of finding a robust set of model parameters we will turn to the example of EtMeSb[Pd(dmit)], a spin-liquid candidate material and part of a family of organic molecular crystals with a rich phase diagram; as well as the spin-liquid phase, these materials have Mott insulating, superconducting, spin density wave and valence bond solid phases Powell and McKenzie (2011); Kobayashi et al. (1991); Seya et al. (1995); Tamura and Kato (2002); Itou et al. (2008). Constructing a coherent picture of this family of materials and their many phases is highly challenging. This effort has been hindered by the fact that, in the usual development of microscopic models, many approximations are made without fully understanding their consequences Jacko et al. (2013a).

The typical approach in EtMeSb[Pd(dmit)] and the related family of materials is to focus on a dimer model, where the dimers of Pd(dmit) sit on a triangular lattice. Parameters are either fit or mapped to this non-interacting model before many-body effects are considered (see for example Refs. Itou et al. (2008); Scriven and Powell (2012)). It has since been shown that a fully-anisotropic triangular lattice (FATL; ) better represents the electronic structure Tsumuraya et al. (2013); Jacko et al. (2013a); Manna et al. (2014). Further, it was shown that a FATL allows one to reproduce the observed many-body properties, predicting a spin-liquid ground state for reasonable parameter values in EtMeSb[Pd(dmit)], while the model does not Jacko et al. (2013a). Fig. 10 shows the phase diagram for the Hubbard model (as a function of ) on the isotropic triangular lattice, triangular lattice, and fully anisotropic triangular lattice (FATL), each with tight-binding parameters consistent with EtMeSb[Pd(dmit)] (computed with variational quantum Monte Carlo). First principles estimates predict Nakamura et al. (2012). The FATL enters the spin-liquid phase at this point, while the and isotropic lattices would predict an insulating phase, with this value of very far from the critical value. Generally, the extra anisotropy seems to destablise the insulating phase relative to the metallic and spin-liquid phases. It is worth noting that these variational quantum Monte Carlo results are not definitive; however, if nothing else, they are indicitive of the important consiquences that even slight parameter changes can have.

Such highly anisotropic models have since become increasingly used in investigations of this family of materials Seo et al. (2015). Having reliable and believable predictions of the degrees of anisotropy in these materials (were the fine variation in parameter values can be attributed to physics and not a quirk of the particular fit one is applying) will be vital for building an understanding of the whole class of materials.

### iv.2 -(Bedt-Ttf) salts: Long range terms

The BEDT-TTF (bis(ethylenedithio)-tetrathiafulvalene, or ET) organic charge transfer salts are a family of quasi-2D crystals that exhibit a wide range of strongly correlated phases (such as non-BCS -wave superconductivity) Mayaffre et al. (1995); Kino and Fukuyama (1996); Vorontsov and Vekhter (2007); Powell and McKenzie (2011); Milbradt et al. (2013); Brown (2015). Understanding this wide range of phases requires a good effective model and good model parameters. It was in these materials that the shortcomings of the approximation (discussed above) were made clear, showing that it leads to a systematic underestimate of Scriven and Powell (2009). Once a realistically large value of is used (computed with cRPA and found to be a 50% - 100% increase over previous estimates), a straightforward Hubbard model of the dimer lattice does not capture anything but the Mott insulating phase Nakamura et al. (2009). These cRPA parameter estimates also showed that the nearest neighbour inter-site Coulomb interactions are significant (), and that they decay slowly with distance Nakamura et al. (2009).

By applying first principles model building techniques, it was found that describing the phases of the ET salts requires models like the extended Hubbard model with significant and long-ranged inter-site interactions. This kind of model, although it has more parameters, it has no more free parameters. Additionally, the inclusion of long-range Coulomb terms has important implications for the energetics of ordered phases Nakamura et al. (2009).

### iv.3 MoS(dmit): An unexpected lattice

In the previous examples, we showed how details in the model parameters are found to have significant consequences on the predicted ground state. We now turn to a system where, by using a first principles method, one finds a totally different lattice model than any previously considered for this system.

MoS(dmit) is a single component molecular crystal that was designed to be metallic. However, it was found to be an activated insulator with an activation energy of 34 meV Llusar et al. (2004). Further, it was found to have no sign of any magnetic order down to very low temperatures ( Llusar et al. (2004)); a possible sign of a spin-liquid state. Based on the apparent 1D physicical properties, its crystal structure, and initial bandstructure calculations, MoS(dmit) was modeled with a one-dimensional lattice Llusar et al. (2004); Janani et al. (2014a, b). This is the 1D ‘triangular necklace’ lattice, illustrated in Fig. 11. The ground state of the Hubbard model on this lattice at filling (as appropriate for MoS(dmit)) is found to be in the Haldane phase, consistant with the experimental evidence Janani et al. (2014a, b).

However, Wannier orbital tight-binding model construction based on density functional theory for MoS(dmit) predicts that at the single electron level, this system is actually 2D with coupling between the 2D layers. The lattice of the 2D layers is an unusual decorated honeycomb lattice, the ‘kagomene’ lattice (illustrated in Fig. 11); interpolating between the graphene (honeycomb) and kagomé lattices Jacko et al. (2015). These lattices have quite different properties, and provide quite different pictures of the physics of MoS(dmit). The kagomene lattice has been studied theoretically before Yao and Kivelson (2007); Wen et al. (2010); Rüegg et al. (2010); Tikhonov and Feigel’man (2010); Yao et al. (2007), but never seen in a real system. This first priciples approach found a layered kagomene lattice in MoS(dmit) quite unexpectedly, demonstrating the novel insights this appoach can yield. The microscopic picture produced is quite different from the phenomenological model.

The one dimensional behaviour can be understood on the grounds of the kagomene lattice: just like kagomé, this lattice has exactly localised states Bergman et al. (2008), illustrated in Fig. 12. Once the 2D kagomene lattice is extended into 3D, these localised states become 1D bands. These emergent 1D states are topological; their degeneracy depends on the boundary conditions of the lattice. Thus, dispite the hopping integrals having similar magnitudes in every direction (in fact, slightly smaller in the stacking direction), one recovers the 1D behaviour and gains some important insights about the potential topological properties of this system.

As a phenomenological model, the necklace lattice does a good job of reproducing the observed magnetic properties of MoS(dmit) Janani et al. (2014b). On the other hand, the kagomene model provides a natural explanation for the quasi-1D behaviour, and highlights the interesting topological flat bands analogous to those seen in the kagomé lattice Jacko et al. (2015). In addition, one can find a lattice closely related to the necklace model as a limiting behaviour of an interacting model in the kagomene lattice, and the many-body behaviour of this model is very similar to the necklace model comm (). One can naturally find new terms to extend the necklace model in a consistent way by introducing higher-order terms in these limits, for example including the chiral next-nearest neighbour terms Jacko et al. (2015).

## V Summary

Over the last decade Wannier orbitals have become an important tool for predictive physics. By constructing Wannier orbitals for frontier bands, we can derive effective models that avoid our biases. These models are robust and reliable, and allow us to make detailed comparisons between materials, and start to extract some general behaviours. They can also find models that we might never have expected to see, leading us to new insights. By moving to this kind of assumption-free methodology for model construction, we can move to a truly predictive approach. We can avoid the dangers of relying on variational Hamiltonians, and allow ourselves to find truely unexpected things.

## Vi Acknowledgments

I would like to thank Ben Powell, Ross McKenzie, Roser Valentí, Harald Jeschke, and Klaus Koepernik for many interesting discussions on these topics. I thank Ben Powell and Ross McKenzie for their helpful comments on this manuscript. I was supported by the Australian Research Council (ARC) through Grant No. DP130100757.

### References

- P. M. Chaikin, E. I. Chashechkina, I. J. Lee, and M. J. Naughton, Journal of Physics: Condensed Matter 10, 11301 (1998).
- B. J. Powell and R. H. McKenzie, Rep. Prog. Phys. 74, 056501 (2011).
- A. C. Jacko, L. F. Tocchio, R. Valentí, and H. O. Jeschke, Phys. Rev. B 88, 155139 (2013a).
- M. Dressel, Naturwissenschaften 94, 527 (2007).
- H. Seo, C. Hotta, and H. Fukuyama, Chemical Reviews 104, 5005 (2004).
- S. Brown, Physica C: Superconductivity and its Applications 514, 279 (2015).
- A. C. Jacko, H. Feldner, E. Rose, F. Lissner, M. Dressel, R. Valentí, and H. O. Jeschke, Phys. Rev. B 87, 155139 (2013b).
- H. Seo, T. Tsumuraya, M. Tsuchiizu, T. Miyazaki, and R. Kato, J. Phys. Soc. Jap. 84, 044716 (2015).
- R. Coldea, D. A. Tennant, E. M. Wheeler, E. Wawrzynska, D. Prabhakaran, M. Telling, K. Habicht, P. Smeibidl, and K. Kiefer, Science 327, 177 (2010).
- H. Taniguchi, A. Kawamoto, and K. Kanoda, Phys. Rev. B 59, 8424 (1999).
- D. Jérome, Science 252, 1509 (1991).
- N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
- G. H. Wannier, Phys. Rev. 52, 191 (1937).
- G. F. Koster, Phys. Rev. 89, 67 (1953).
- W. Kohn, Phys. Rev. 115, 809 (1959).
- E. Hückel, Zeit. für Physik 70, 204 (1931).
- J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).
- W. Kohn, Chem. Phys. Lett. 208, 167 (1993).
- W. Kohn, Int. J. Quant. Chem. 56, 229 (1995).
- J. D. Cloizeaux, Phys. Rev. 135, A685 (1964).
- N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
- I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001).
- A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comp. Phys. Comm. 178, 685 (2008).
- K. Koepernik and H. Eschrig, Phys. Rev. B 59, 1743 (1999).
- J. Kuneš, R. Arita, P. Wissgott, A. Toschi, H. Ikeda, and K. Held, Comp. Phys. Comm. 181, 1888 (2010).
- P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al., J. Phys.: Cond. Matt. 21, 395502 (2009).
- X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken, F. Bottin, P. Boulanger, F. Bruneval, D. Caliste, R. Caracas, M. CÃ´tÃ©, et al., Comp. Phys. Comm. 180, 2582 (2009).
- F. Freimuth, Y. Mokrousov, D. Wortmann, S. Heinze, and S. Blügel, Phys. Rev. B 78, 035120 (2008).
- J. Hinze and D. L. Beveridge, Journal of the American Chemical Society 93, 3107 (1971).
- H. C. Kandpal, I. Opahle, Y.-Z. Zhang, H. O. Jeschke, and R. Valentí, Phys. Rev. Lett. 103, 067004 (2009).
- E. P. Scriven and B. J. Powell, Phys. Rev. Lett. 109, 097206 (2012).
- H. O. Jeschke, M. de Souza, R. Valentí, R. S. Manna, M. Lang, and J. A. Schlueter, Phys. Rev. B 85, 035125 (2012).
- F. Dyson, Nature (London) 427, 297 (2004).
- J. Mayer, K. Khairy, and J. Howard, American Journal of Physics 78, 648 (2010).
- Y. Nogami, T. Ito, K. Yamamoto, N. Irie, S. Horita, T. Kambe, N. Nagao, K. Oshima, N. Ikeda, and T. Nakamura, J. Phys. IV France 131, 39 (2005).
- M. K. Transtrum, B. B. Machta, K. S. Brown, B. C. Daniels, C. R. Myers, and J. P. Sethna, The Journal of Chemical Physics 143, 010901 (2015).
- Tight binding fitting code written by H. O. Jeschke, used with permission.
- K. Nakamura, Y. Yoshimoto, T. Kosugi, R. Arita, and M. Imada, Journal of the Physical Society of Japan 78, 083710 (2009).
- K. Nakamura, Y. Yoshimoto, and M. Imada, Phys. Rev. B 86, 205117 (2012).
- M. Altmeyer, R. Valentí, and H. O. Jeschke, Phys. Rev. B 91, 245137 (2015).
- A. C. Jacko, B. J. Powell, and R. H. McKenzie, J. Chem. Phys. 133, 124314 (2010).
- T. Komatsu, N. Matsukawa, T. Inoue, and G. Saito, Journal of the Physical Society of Japan 65, 1340 (1996).
- R. H. McKenzie, Comments Cond. Mat. Phys. 18, 309 (1998).
- E. Scriven and B. J. Powell, Phys. Rev. B 80, 205107 (2009).
- K. Yoshimi, H. Seo, S. Ishibashi, and S. E. Brown, Phys. Rev. Lett. 108, 096402 (2012).
- T. Tsumuraya, H. Seo, M. Tsuchiizu, R. Kato, and T. Miyazaki, Journal of the Physical Society of Japan 82, 033709 (2013).
- K. Nakamura, R. Arita, Y. Yoshimoto, and S. Tsuneyuki, Phys. Rev. B 74, 235113 (2006).
- H. DeVoe, J. Chem. Phys. 41, 393 (1964).
- L. Cano-Cortés, A. Dolfen, J. Merino, and E. Koch, Physica B 405, S185 (2010).
- F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann, and A. I. Lichtenstein, Phys. Rev. B 70, 195104 (2004).
- D. Bohm and D. Pines, Phys. Rev. 82, 625 (1951).
- D. Pines and D. Bohm, Phys. Rev. 85, 338 (1952).
- D. Bohm and D. Pines, Phys. Rev. 92, 609 (1953).
- D. Pines, Phys. Rev. 92, 626 (1953).
- M. Gell-Mann and K. A. Brueckner, Phys. Rev. 106, 364 (1957).
- RPA bubble diagram produced by wikimedia user ‘Condmatstrel’, and released under the Creative Commons Attribution 3.0 license.
- T. Miyake, F. Aryasetiawan, and M. Imada, Phys. Rev. B 80, 155134 (2009).
- T. Miyake and F. Aryasetiawan, Phys. Rev. B 77, 085122 (2008).
- K. Nakamura, R. Arita, and M. Imada, Journal of the Physical Society of Japan 77, 093711 (2008).
- A. Kobayashi, H. Kobayashi, A. Miyamoto, R. Kato, R. A. Clark, and A. E. Underhill, Chem. Lett. 20, 2063 (1991).
- K. Seya, Y. Kobayashi, T. Nakamura, T. Takahashi, Y. Osako, H. Kobayashi, R. Kato, A. Kobayashi, and H. Iguchi, Synth. Met. 70, 1043 (1995).
- M. Tamura and R. Kato, J. Phys.: Condens. Matt. 14, L729 (2002).
- T. Itou, A. Oyamada, S. Maegawa, M. Tamura, and R. Kato, Phys. Rev. B 77, 104413 (2008).
- R. S. Manna, M. de Souza, R. Kato, and M. Lang, Phys. Rev. B 89, 045113 (2014).
- H. Mayaffre, P. Wzietek, D. Jérome, C. Lenoir, and P. Batail, Phys. Rev. Lett. 75, 4122 (1995).
- H. Kino and H. Fukuyama, Journal of the Physical Society of Japan 65, 2158 (1996).
- A. B. Vorontsov and I. Vekhter, Phys. Rev. B 75, 224502 (2007).
- S. Milbradt, A. A. Bardin, C. J. S. Truncik, W. A. Huttema, A. C. Jacko, P. L. Burn, S.-C. Lo, B. J. Powell, and D. M. Broun, Phys. Rev. B 88, 064501 (2013).
- R. Llusar, S. Uriel, C. Vicent, J. M. Clemente-Juan, E. Coronado, C. J. Gómez-GarcÄ±Â´a, B. Braïda, and E. Canadell, J. Am. Chem. Soc. 126, 12076 (2004).
- C. Janani, J. Merino, I. P. McCulloch, and B. J. Powell, Phys. Rev. Lett. 113, 267204 (2014a).
- C. Janani, J. Merino, I. P. McCulloch, and B. J. Powell, Phys. Rev. B 90, 035120 (2014b).
- A. C. Jacko, C. Janani, K. Koepernik, and B. J. Powell, Phys. Rev. B 91, 125140 (2015).
- H. Yao and S. A. Kivelson, Phys. Rev. Lett. 99, 247203 (2007).
- J. Wen, A. Rüegg, C.-C. J. Wang, and G. A. Fiete, Phys. Rev. B 82, 075125 (2010).
- A. Rüegg, J. Wen, and G. A. Fiete, Phys. Rev. B 81, 205115 (2010).
- K. S. Tikhonov and M. V. Feigel’man, Phys. Rev. Lett. 105, 067207 (2010).
- N. Yao, C. Laumann, A. Gorshkov, H. Weimer, L. Jiang, J. Cirac, P. Zoller, and M. Lukin, Nat. Comms. 4, 1585 (2007).
- D. L. Bergman, C. Wu, and L. Balents, Phys. Rev. B 78, 125104 (2008).
- H. Nourse, in preparation.