Competing collinear magnetic structures in superconducting FeSe by first principles quantum Monte Carlo calculations
Abstract
Resolving the interplay between magnetic interactions and structural properties in strongly correlated materials through a quantitatively accurate approach has been a major challenge in condensed matter physics. Here we apply highly accurate first principles quantum Monte Carlo (QMC) techniques to obtain structural and magnetic properties of the iron selenide (FeSe) superconductor under pressure. Where comparable, the computed properties are very close to the experimental values. Of potential ordered magnetic configurations, collinear spin configurations are the most energetically favorable over the explored pressure range. They become nearly degenerate in energy with bicollinear spin orderings at around 7 GPa, when the experimental critical temperature is the highest. On the other hand, ferromagnetic, checkerboard, and staggered dimer configurations become relatively higher in energy as the pressure increases. The behavior under pressure is explained by an accurate analysis of the charge compressibility and the orbital occupation as described by the QMC manybody wave function, which reveals how spin, charge and orbital degrees of freedom are strongly coupled in this compound. This remarkable pressure evolution suggests that stripelike magnetic fluctuations may be responsible for the enhanced in FeSe and that higher T is associated with nearness to a crossover between collinear and bicollinear ordering.
I Introduction
The quest for a microscopic theory of unconventional or hightemperature superconductivity is a major challenge in condensed matter physics. The discovery of ironbased superconductors in 2006 Kamihara et al. (2006) was an important contribution to the field since it added a second class of hightemperature unconventional superconductors to the experimental roster, along with the cuprate superconductors. Despite their different electronic structure, their phase diagrams have striking similaritiesNorman (2008); de’ Medici et al. (2014), particularly the proximity of the superconducting phase with an antiferromagnetic state. This behavior, along with other considerationsScalapino (2012); Wen et al. (2011); Dai (2015); Dai et al. (2012); Imai et al. (2009), makes it likely that spins and magnetism are important in determining the superconducting state.
FeSe is a particularly interesting example of the ironbased superconductors for several reasons. Its critical temperature is strongly dependent on pressureMizuguchi et al. (2008); Imai et al. (2009); Margadonna et al. (2009), reaching 37 K at 7 GPa. At ambient conditions, FeSe has a simple P4/nmm crystal structure with two inequivalent Fe and Se positions per unit cell, and it undergoes a distortion from tetragonal to orthorhombic symmetry by cooling it down below 91 K, while it becomes superconductors below 8 K at ambient pressure. An intriguing peculiarity of FeSe is that, at variance with most of the ironbased superconductors, it does not show any longrange magnetic order over the whole phase diagramMedvedev et al. (2009). In spite of this, very strong antiferromagnetic spin fluctuations have been revealed by neutron scattering experiments (see for example Ref. Rahn et al., 2015a) in the proximity of the superconducting phase. Their role in driving the nematic transition and their connection to superconductivity have been the subject of intense debate. All these aspects make it attractive for computational techniques to correlate microscopic electronic structure with the superconductivity and it is therefore one of the most studied ironbased superconductors. However, the precise calculation of the properties of this material remains challenging from first principles methods such as density functional theory (DFT) due to strong electron correlation.
For example, the PBE band structure is in poor agreement with experiments which report a considerably narrower bandwidthMaletz et al. (2014a); Audouard et al. (2015). Furthermore the FeSe lattice constants display an average error of 0.1 Å independently from the exchange correlation functional employed (see for instance Ref Caglieris et al. (2012) and Tab 1). Despite useful work using dynamical mean field theoryLiebsch and Ishida (2010); Yin, Z. P. et al. (2011); Aichhorn et al. (2010); Craco and Leoni (2010); Craco et al. (2014); Leonov et al. (2015) and GWTomczak et al. (2012); Hirayama et al. (2013); H. et al. (2015) methods, there is a strong need for high quality calculations that can better describe the electronic and crystal structure of these materials.
In this article, we describe the results of first principles quantum Monte Carlo simulations of the magnetic behavior of FeSe under pressure. The main method used in this article, fixed node diffusion Monte Carlo (FNDMC), has been shown recently to offer very accurate results on a number of challenging materials, including VOZheng and Wagner (2015), cuprates Wagner (2015); Foyevtsova et al. (2014) and other transition metal oxides as well as rare earths as cerium Devaux et al. (2015). Furthermore, a recent work Casula and Sorella (2013), based on quantum Monte Carlo techniques, successfully tackled the problem of pairing symmetry in FeSe itself.
We find that, compared to commonly used density functional theory calculations, the FNDMC calculations obtain more accurate lattice constants, bulk moduli, and band dispersion. By increasing the pressure, the difference in energy of ordered magnetic states with stripelike order goes to zero with pressure, while checkerboardlike magnetic states increase in energy. The convergence of the stripelike magnetic states is correlated with the increase in T in this material under pressure, which offers a tantalizing connection to spin fluctuations as a potential origin. Such behavior may be a calculable design principle for new unconventional superconducting materials.
Ii Methods
ii.1 Fixednode diffusion Monte Carlo
(a)  (b)  (c)  (d) 
Fixednode diffusion Monte Carlo Foulkes et al. (2001) is a stochastic approach to solving the exact firstprinciples imaginarytime Schrödinger equation. FNDMC is a variational (upper bound to the ground state) method with no adjustable parameters. FNDMC is a projectorbased approach, which projects out the groundstate wave function via repeated application of the operator on some trial wavefunction :
With no further approximation, this method is exact, however, it suffers from the wellknown fermion sign problem. This is circumvented in FNDMC with the fixednode approximation Anderson (1976); Ceperley and Alder (1984), which constrains the projected wave function to have the same nodes as the starting trial wave function. If the nodal structure of the trial wave function coincides with the groundstate one, the method remains exact. In practice, for accurate trial wavefunctions, this approximation introduces small errors that we will estimate.
We used a single SlaterJastrow trial wave function ansatz as follows:
(1) 
where , refer to electrons, refers to nuclei, and indicate spin up/down respectively. Since the FNDMC result is determined by the nodes of the determinants in Eq. 1, the orbitals determine the degree of the fixed node approximation. To test the effect of these orbitals, we use two approaches. The first optimizes a parameter in density functional theory to generate the orbitals, and is the less computationally demanding of the two. The parameter to be optimized in the Slater determinant is the amount of exact exchange , for which we find the optimal value near . This corresponds to the PBE0 functional, and this optimum is often the case for similar calculations on transition metal systems KolorenÄ et al. (2010). For consistency, we used for all our calculations, and we denote this method by DMC(PBE0). The details of the optimization of functional and of the parameter, , is reported in the Supplemental Material. The second QMC approach used in this work performs a complete variational optimization of the determinant orbitals within a relatively small basis set, and is more computationally demanding but in principle more accurate. These calculations are denoted throughout the paper with the label QMC(opt), where QMC will refer to VMC or DMC depending on the calculation. Further information on the orbital optimization procedure is provided in the Supplementary Material. In our tests for FeSe, while the DMC(opt) technique did obtain lower energies as expected, the energy differences were consistent between the two techniques, so most data is obtained from DMC(PBE0).
All DMC(PBE0) calculations were done within the opensource package QWalk Wagner et al. (2009), with orbitals generated by DFT calculations performed with the DFT code CRYSTAL R. Dovesi (2014). For the DMC(opt) method, we used the package TurboRVB Sorella ().
Our only approximation to the Hamiltonian is a HartreeFock pseudopotential designed specifically for quantum Monte Carlo calculations Burkatzki et al. (2007, 2008). The energy difference between the collinear and checkerboard magnetic state is consistent between an allelectron and pseudopotential PBE0 calculation within 0.01 eV. The convergence of the most important parameters in both our QMC methods is shown in Sec.A.1 of the Supplementary Material. For the FeSe crystal structure, the anion height above the iron planes is the only internal parameter of the compound in the tetragonal P4/nmm phase. This parameter represents a crucial ingredient to determine the magnetic behavior of FeSe, but its evaluation by firstprinciples methods is a challenging task, as detailed in Sec.III.2. The optimization of the Se height is carried out with two different procedures. For the DMC(PBE0) method, the relaxed value is obtained by fitting a total energy curve with a cubic function. In the VMC(opt), the optimized Se height value is obtained by a direct minimization of the ionic forces within the variational Monte Carlo framework Sorella and Capriotti (2010). By including the cell parameters in the minimization procedure, we are able to fully relax the crystal structure of FeSe at different magnetic orderings. We consider a minimization converged when both the forces and their error bars are lower than 10 Ha/a.u. per atom.
We found that effects due to the finite size of the simulation cell, or finite size errors (FSEs), constitute the major source of systematic error for both DMC(PBE0) and DMC(opt). We apply several techniques in order to reduce FSEs at best. All DMC(PBE0) calculations are twist averaged Lin et al. (2001) over 8 twist conditions on the 8 f.u. FeSe supercell; DMC(opt) results are instead obtained with a larger 16 f.u. supercell averaging over periodic and antiperiodic boundary conditions; further corrections are then applied to cure onebody and twobody FSEs. In both cases, we managed to reduce the impact of FSEs below the desired accuracy on energy calculations. A detailed explanation of the procedures used to control FSEs is given in the Supplemental Material.
Iii Results
iii.1 Trial wave functions and ground state
For the wave function in Eq. 1, there are many local minima both in preparing the KohnSham orbitals using density functional theory and in optimizing the orbitals directly. These minima correspond to different magnetic orderings of the Fe spins. The most relevant ones are presented in Fig. 1. We also included a type of paramagnetic state in which the up and down spin orbitals are constrained to be equal, but we found that state to be more than 0.5 eV/Fe higher in energy than any magnetically ordered wavefunction. The ground state thus seems to require large local moments on the Fe atoms.
While it is known experimentally that FeSe does not have longrange ordering Medvedev et al. (2009), the calculations here enforce periodic boundary conditions on a relatively small cell and thus cannot describe longrange fluctuations of the magnetic order that might be the cause of loss of longrange order. For the experimental crystal structure, the collinear magnetic ordering is the lowest in energy in our calculations and is observed to be the dominant shortrange order experimentally Rahn et al. (2015b). The energetic cost of introducing a “defect” into the magnetic order is quite small; we will discuss that aspect later. Both the DMC(opt) and DMC(PBE0) approaches result in a rather large magnetic moment on the Fe atom. For the collinear magnetic ordering we obtain a value of 3.4 for DMC(PBE0), and a slightly lower 3.1 for the fully optimized calculations. In both cases the magnetic moment is close to the atomic limit.
Between the two DMC approaches, the energy difference between different magnetic orderings is in agreement within stochastic errors, so there is good reason to believe that the cheaper DMC(PBE0) technique is accurate. In comparison to PBE calculations, which are the most common in the literature, the relative energies according DMC are quite different, including the lowest energy magnetic phase, which is the “staggered dimer” configuration in DFT Cao et al. (2015); Glasbrenner et al. (2015); Tresca et al. (2015), but turns out to be the collinear configuration in DMC. It appears that hybrid DFT calculations in the PBE0 approximation obtain reasonably good magnetic energy differences in comparison to DMC; since this functional also produced the orbitals that gave the lowest FNDMC energy, it may be capturing some of the correct physics for the magnetic properties of this material. However, the PBE0 functional predicts an insulating gap in agreement with previous workWu (2013) for FeSe for all magnetic orderings, in contrast to DMC and experiment.
iii.2 Crystal structure
Source  Magnetic Ord.  a  c  

DFTPBE  paramagnetic  3.6802  6.1663  2.6023  1.3862 
DFTPBE  collinear  3.8007  6.2363  2.6966  1.4568 
VMC  paramagnetic  3.71(1)  5.49(1)  2.62(1)  1.437(5) 
VMC  collinear  3.72(1)  5.68(1)  2.63(1)  1.56(1) 
experimental Louca et al. (2010)  T 7 K  3.7646(1)  5.47920(9)  1.4622  
experimental Kumar et al. (2010)  T 8 K  3.7685(1)  5.5194(9)  2.6647(3)  1.5879  
experimental Margadonna et al. (2009)  T 300 K  3.7724(1)  5.5217(1)  1.4759 
Obtaining the correct crystal structure for FeSe is a major challenge, since the layers interact through nonbonded interactions. The lattice parameter in particular is affected by Van der Waals interactions and electron correlation plays a key role in determining the inplane physics. The exotic behavior of FeSe superconducting properties under pressure gives another clue on the importance of structural variations in its description. A firstprinciples prediction of the lattice parameters is thus an important test of the description of this physics. Since the DMC calculations are computationally costly, we limited our study to the tetragonal phase of FeSe. Because the low temperature orthorhombic distortion is smallMargadonna et al. (2009), one might expect that its effect on the overall electronic structure is also small. We leave such considerations to another paper.
The equilibrium lattice parameters of FeSe are presented in Table 1. As mentioned in the previous sections, these results are obtained with a direct optimization of FeSe cell parameters with the VMC(opt) method. The inplane FeSe properties should be well captured by QMC since the lattice parameter is in close agreement with experimental results (within 4 ) independently of the chosen magnetic configuration. Both collinear and paramagnetic wavefunctions shows also a general improvement with respect to DFT concerning the lattice parameter. This provides evidence of the accuracy in treating Van der Waals interactions with the QMC wavefunction, mainly achieved with the Jastrow factorBenali et al. (2014); Sorella et al. (2007). The evaluation of the interplane distance might be affected by the dispersion along the zaxis, which we did not take into account in our supercell which always contains only one Fe plane. We check this dependency by performing a test structural relaxation on a FeSe supercell with 16 Fe atoms in two planes and 8 Fe atoms with only one plane considered in the supercell. We find that the difference between the parameter obtained in the two configurations is negligible.
The final internal parameter represents the height of the selenium anion above the plane, and it has been experimentally demonstrated Okabe et al. (2010) to be of key importance in determining superconducting properties of ironbased superconductors in general. We collect all our calculations of , as well as some experimental results, in Fig. 2. We find that both the magnetic state and the accuracy of the calculation have an important effect on the prediction of this parameter. At approximately the same level of finite size error, our two DMC calculations agree very closely, determining that fixed node and basis set error is likely to be unimportant. However, we found that is surprisingly sensitive to finite size effects, both in the inplane and out of plane directions. Given the supercells that we studied, we found a variation in of approximately 0.05 Å, depending on the twisted boundary conditions and supercell. With experimental lattice parameters, our best estimate for is thus 1.54(5)Å, which is quite close to the experimental range. As we shall see later, the properties of FeSe depend sensitively on , so to account for this uncertainty, we will consider properties as a function of selenium height as well as pressure.
By fitting an equation of state previously used by Anton et al. Anton and Schmidt (1997) to our DMC(PBE0) energies as a function of volume, we extract the bulk modulus and the pressure dependence on volume , shown in Fig 3. The collection of ambientpressure bulkmoduli results is reported in the inset of Fig. 3, in units of GPa. For all these calculations, experimental lattice constants Kumar et al. (2010) have been used. and the bulk modulus show a strong dependence on the magnetic order.
While has scatter between experiments, they are more consistent in the bulk modulus, so we base our comparisons of the theoretical calculations on the latter quantity. The DMC(PBE0) calculation demonstrates excellent agreement with all three experiments if the collinear magnetic ordering is imposed, but it is less close to experiment for the other magnetic orderings. Our PBE0 calculations are also in somewhat good agreement with DMC(PBE0), except a notable disagreement for the ferromagnetic ordering. On the other hand, PBE bulk moduli are significantly lower than both experiment and the other calculations, generally predicting bulk moduli between 7 and 10 GPa, depending only slightly on the magnetic ordering. Since the collinear ordering is also the lowest energy for DMC(PBE0), for the remainder of this article, we use the collinear equation of state to estimate the pressures that correspond to the volumes used in the calculations.
iii.3 Interaction of structure and magnetism
Fig 4 shows the interaction between pressure, magnetic ordering, and selenium height. As has been found before Lischner et al. (2015), the magnetic energies depend strongly on the selenium height, and this dependence changes with pressure. For a given pressure, there are two points on the selenium height curves that are of substantial interest. The first is the minimum energy (solid vertical line), which as can be seen in Fig 2, does not change very much with pressure(or volume) in our calculations. The second is the crossing point (dashed vertical line) between the collinear and bicollinear magnetic orderings, which are competing lowenergy states. This crossing point depends on the pressure, and approaches the minimum energy point at higher pressures(lower volumes), as shown on the rightmost plot in the figure.
Another interesting feature visible in Fig 4 is that the checkerboard magnetic ordering intersects the bicollinear and collinear magnetic orders at zero pressure (78.4 Å) and large , but there is a shift of the checkerboard curve to higher energies once pressure is applied. The underlying physics of this effect will be discussed in Section III.5.

Fig 5 shows a cut through the data in Fig 4 along the minimum energy (subfigure a), and the experimentally determined (subfigure b). Along this cut we evaluated many magnetic orderings to establish a set of trends, and checked finite size errors by considering an 8 f.u. cell and 16 f.u. cell with twist averaging. Further information on finite size corrections are available in the Supplemental Material. Under pressure, the checkerboard, ferromagnetic, and staggered dimer magnetic orderings rise in energy compared to the lowest energy collinear ordering. On the other hand, the stripelike orderings, including the bicollinear and collinear orderings with defects converge with applied pressure.
From Fig 5 (bottom panels) is apparent the failure of PBE in capture this trend in FeSe energetics under pressure. Even with lattice constants fixed to experimental ones, the PBE energies of magnetically ordered states are quite different from the FNDMC energies. In agreement with recent work, PBE does predict the staggered dimer as ground state. Despite the failure of PBE0 in describing the conducting behavior of FeSe, the magnetic energies are reasonably close to the DMC results.
Given the data available to us, we can determine some properties that are robust to the finite size errors and uncertainty in in our calculations. The first is that the relative energetics of magnetic orders changes strongly as a function and pressure. In FNDMC and PBE0, which would a priori be expected to be more accurate, the collinear and bicollinear orders become degenerate as a function of pressure for reasonable values of . According to FNDMC, this effect is robust against variations, depending mainly on the change in the relative magnetic energies as a function of pressure.
The energetic cost of reversing a single spin in the collinear ordered state, labeled “collinear, flip 1” in Fig 5a, follows the bicollinear energy quite closely. Because this cost decreases with pressure, we can surmise that magnetic fluctuations become more energetically available as pressure is increased.
iii.4 Optical excitations and magnetism
The direct optical gap was calculated by promoting the highest energy orbital in the Slater determinant part of the trial wavefunction to the next excited state orbital. This constructs a wave function ansatz for an electronhole excitation. The results are shown in Fig 6. The resulting DMC(PBE0) energy relative to the DMC(PBE0) ground state is our estimation of the gap. Interestingly, the DMC(PBE0) gap is within statistical uncertainties of 0 despite the fact PBE0 estimates a rather large gap, regardless of magnetic ordering. Experimentally Shimojima et al. (), the gap is no more than 80 meV at any point, which is consistent with our results for the bicollinear and collinear magnetic ordering. Only the checkerboard state is gapped according to DMC(PBE0).
The charge degrees of freedom are therefore coupled to the spin degrees of freedom. According to these calculations, in FeSe there is a coupling between the mobility of charge and the spin ordering. In the remainder of the paper, we will correlate these properties with those of the ground state for different spin orderings.
iii.5 Interaction of charge and orbitals with magnetism



From the energetic properties, we note two classes of magnetic order in FeSe: ones which are stripelike, and ones which are not stripelike. The stripelike orderings converge in energy with pressure, while the checkerboard and staggered dimer pattern increases in energy relative to those orderings. Similarly, the gap calculated in DMC(PBE0) distinguishes between different orderings, with metallic character in the stripelike ordering. In this section, we will make the following observations:

Delocalized minority spin electrons are most affected by the spin ordering.

The oneparticle orbitals are occupied differently depending on the spin ordering.

The above two effects change the effect of correlation in FeSe.
These effects combine to give a cartoon picture of the physics that explains the difference in pressure behavior between the magnetic orders.
Delocalization.
To address delocalization, we evaluate the charge compressibility of the Fe sites: , where is the number of electrons within a Voronoi polehedron around the th Fe site of spin . Larger values of the compressibility indicate more delocalized electrons. For a Fe atom with net spin, the electrons are labeled majority electrons and the minority, and vice versa for Fe atoms with net spin. In Fig 7a, these results are presented. For all magnetic orders, the majority spin is very similar and shows a low charge compressibility, while the minority spin is different between different magnetic orders, and its charge compressibility is larger. From this, the minority electrons are more delocalized and their localization is affected by the magnetic order. For the stripelike orders, the minority electrons are the most delocalized.
One particle orbitals.
In Fig 7d, we present the orbital occupation of the orbitals in different spin orderings. For stripelike orderings, the , , and orbitals are occupied, in agreement with ARPES results Maletz et al. (2014b). On the other hand, the orbital is occupied for the checkerboard ordering. This gives a simple explanation for the delocalization patterns: The checkerboard pattern causes the minority spin to occupy the out of plane orbital, which would rise to an insulating state if it were the ground state. This idea can be confirmed by checking the offdiagonal oneparticle density matrix elements between Fe atoms with parallel and antiparallel net spins, in Fig 7e. The atomic orbitals are more hybridized between parallel spin Fe atoms for the stripelike orders. The charge degrees of freedom, which are mainly the minority spins from the Fe, interact strongly with the magnetic ordering. This effect also interacts with the net magnetic moment and onsite correlations (Fig 7b and Fig 7c).
iii.6 A cartoon picture of FeSe
A simple picture based on Hund’s coupling can explain the energetics and other properties presented in the results section. Hund’s rules dictate that for an atom with a partially filled shell, we expect the electrons to have total spin that maximizes the multiplicity . This is consistent with our computed magnetic moment, which find that the majority channel is mostly filled, bringing the moment to around 3.13.4 . The spin occupation of the states in a reference iron is diagrammatically shown in the top row of Fig. 8. Also due to Hund’s coupling, the electron that is most likely to hop to nearby iron atoms would be the electron in the minority channel, to keep a large . As illustrated in Fig. 8, this minority channel is already filled for neighboring irons that are antiparallel, so only majority spin electrons can hop to those atoms. Conversely, minority electrons can hop to neighboring parallel irons, since that spin channel is not filled. Thus, irons with parallel spins allow the minority electrons to delocalize, therefore decrease the kinetic energy. As seen in Fig. 7d the magnetic ordering affects the occupation of the states, hence affects the labeling of the states in Fig. 8, but the basic idea is unchanged.
While the minority spins require at least some parallel iron magnetic moments, the large localized magnetic moments also interact antiferromagnetically, leading to a competition between these two mechanisms. As a compromise, antiferromagnetic configurations with ferromagnetic chains emerge as the lowest energy configurations.
This picture unifies many of the observations from our calculations. The checkerboard state is distinguished from the other states by its lack of parallel nearest neighbors, similar to how the ferromagnetic state is distinguished by its lack of antiparallel neighbors. These two extremes are higher in energy, and are disfavored as pressure increases the importance of FeFe interactions. Because the checkerboard has no parallel nearest neighbor, its iron states are localized on the irons, leading to a low charge variance, and states that primarily occupy the orbital. All stripelike states have a combination of antiparallel and parallel nearest neighbors, and allow the electrons to delocalize along the irons chains, leading to higher correlations, higher variance, and more FeFe hybridization. Although the staggered dimer ordering is energetically competitive at low pressures, its energy, charge variance, and magnetic moment are similar to checkerboard, and at high pressures, becomes energetically unfavorable just as the checkerboard ordering does. Although the staggered dimer does allow some delocalization between the dimered parallel spins, the itinerant spins are still trapped on the dimers, and therefore this state’s energetics follow the checkerboard behavior at higher pressures.
This competition of interactions sets up a fine balance between many qualitatively different magnetic configurations. Parameters in the structure can tilt this balance one way or another, leading to a strong magnetostructural coupling. This is evident both from the strong magnetic dependence of the bulk modulus in Fig. 3, as well as in Fig. 4, where can exchange the ground state configuration between at least two magnetic configurations. This logic can be straightforwardly applied to iron telluride (FeTe), the non superconducting parent compound of FeSe. This material has the ground state magnetic ordering at a around Å, which implies that the should be decreased to force a crossover. Therefore FeTe would superconduct if it were put in tensile stress, as it has been observed Han et al. (2010).
Iv Conclusion
In summary, we have shown that QMC calculations can obtain an accurate description of the electronic structure of FeSe. The lattice constants, bulk modulus, and bandwidth are all very close to the experimental values and significantly improve over DFT calculations. Our results are substantiated by the agreement between different and complementary QMC techniques employed. We showed that they yield sufficiently small statistical and systematic errors to study the relative energetics of different magnetic orders, which behave differently from those predicted by DFT. The largest error in the calculations appear to be due to finite size supercells, which we checked to be small enough that the trends presented here are preserved.
As an outcome of the highaccuracy calculations, we have determined that collinear and bicollinear motifs become close in energy as pressure increases, while the checkerboard motif increases in energy with pressure. This behavior is correlated with delocalization of the minority electrons on the highspin Fe atom. Collinear and bicollinear motifs allow for more delocalization, which increases in importance as the pressure is increased. This delocalization effect is strong enough to change the occupation of atomic orbitals in FeSe depending on the magnetic ordering, so it is larger than the crystal field splitting of the orbitals. The spontaneous breaking of C symmetry (or more properly S symmetry)Casula and Sorella (2013); Wang et al. (2015), is a result of this physics. Magnetic configurations which contain spin chains are thus favored over the whole considered range of pressure.
From the above results, we can see that the magnetic degrees of freedom are strongly coupled with the charge and orbital degrees of freedom. In a similar way, since the relative magnetic energies are dependent on the hopping of minority electrons from site to site, they are also strongly dependent on the structure. There is thus both spincharge and spinstructural coupling in this system. As one of us showed recentlyWagner and Abbamonte (2014); Wagner (2015), the cuprates also show strong magnetostructural and magnetocharge coupling. One might speculate that both of these effects are necessary for high T and it may prove fruitful to look for similar effects in proposed new superconductors.
Acknowledgements.
We would like to thank the many people who have given useful suggestions and comments throughout this work. They include B.B. would like to thank the NSF Graduate Research Fellowship for funding. L.K.W. was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program under Award Number FG0212ER46875. Computational resources were provided through the INCITE PhotoSuper and SuperMatSim programs. S.S, M.C. and M.D. acknowledge computational resources provided through the HPCI System Research Project (Nos. hp120174 and hp140092) on K computer at RIKEN Advanced Institute for Computational Science, and on the HOKUSAI GreatWave computer under the project G15034. M. C. thanks the GENCI french program to provide additional computer time through the Grant No. 2014096493.Appendix A Supplemental Material
a.1 Convergence and validation



In this section we present an extensive study of the convergence of the main parameters involved in our QMC calculations. Within DMC(PBE0), the only parameter to be optimized in the Slater determinant is the amount of exact exchange . The optimization of FNDMC energy as a function of is presented in Fig 10. As already mentioned in the main manuscript, we notice that the best FNDMC is generally obtained with of exact exchange for all magnetic configurations; this corresponds to the PBE0 density functional. The convergence of the FNDMC energy with other parameters of QMC methods is reported in Fig 9. For the DMC(PBE0) method, in Fig. 9b, we present energy convergence as a function of Fe and Se basis set exponents by showing the behavior of the exponent of our most diffuse gaussian basis exponent as it gets more diffuse. Benchmark calculations of DMC(PBE0) against the time step used in the projection are presented in Fig 9d. We employed a time step of 0.01 for all DMC(PBE0) calculations. On the other hand, lattice regularized FNDMC algorithm employed for DMC(opt) method (see next section) suffers from the lattice step error in Laplacian discretization. Convergence with lattice step is shown for the collinear configuration in Fig 9f. We used a lattice step of a.u. for all DMC(opt) calculations.
We turn now our attention to finite size errors (FSE) which represent the main source of error in our QMC calculations. We performed several DFT test calculations to determine the impact of onebody FSE. For both 8 f.u. and 16 f.u. supercells (used for DMC(PBE0) and DMC(opt) methods respectively) a 4x4x4 kpoint grid is enough to obtain results converged within 1 meV independently from the density functional used, as shown in 9a. The same kpoint grid in the PBE0 wavefunction is also sufficient to converge FNDMC energies within the same threshold. For the DMC(PBE0) calculations, we twist average over a set of 8 twist conditions Lin et al. (2001), for unit cells ranging from 4 to 16 f.u., expanding the supercell in the  plane, and direction (adding an additional layer). The resulting finitesize extrapolations are depicted in Fig 9c. The finite size extrapolation in the direction is the checkerboard line above the other checkerboard line. The true infinite size limit will likely lie in between these two extrapolations. Although the finite size effects are relevant, they do not alter any conclusions of the main text. Going from low to high pressure, the finite size errors amplify the change in energy differences between checkerboard and collinear. Extrapolating with more than two points is prohibitively expensive for the bicollinear state, and so how finite size effects affects the pressure dependence of the energy differences between collinear and bicollinear is not clear, as the extrapolations are within error bars. However, it is certainly clear that bicollinear both remains very close in energy to the collinear state—even in the extrapolation—whereas checkerboard’s energy certainly rises well above the other states.
For DMC(opt), all calculations have been done with a 16 f.u. supercell. Structural optimization was performed with periodic boundary conditions only. All other calculations were instead averaged between periodic (PBC) and fully antiperiodic (APBC) boundary conditions. Further finite size corrections to DMC energies are obtained by adding onebody corrections estimated from fully converged DFTLDA calculations and twobody corrections evaluated within the KZK approach Kwee et al. (2008); Sorella et al. (2011). In Figure 9e is reported the convergence of DMC(opt) energy in the paramagnetic phase as a function of the system size after applying the corrections.
a.2 QMC full orbital optimization
This section reports a detailed explanation of the QMC(opt) method. The wavefunction optimization procedure is carried out entirely with the variational Monte Carlo scheme. The initial Slater determinant is taken from a DFT calculation in the local density approximation (LDA) obtained with a builtin DFT package. The Slater determinant is expanded over a singleparticle, localized gaussian basis set where is the basis set index and is the nuclear index. The number of gaussian functions employed for each angular momentum channel is chosen to attain the best compromise between basis set size and DFT(LDA) total energy. The basis set we use is for Fe and for Se. Following the notation, the round brackets contain the primitive basis. The set of variational parameters involved in the orbital optimization can be highly reduced by building atomic basis contractions, i.e. linear combinations of the primitive orbitals. The contracted basis employed for Fe and Se is contained in the square brackets. The procedure for finding the optimal number of contracted orbitals will be explained in detail in a manuscript is going to be published Sorella et al. (2015). The Jastrow factor is multiplied to the DFT(LDA) determinant to build the final QMC ansatz. For the orbital optimization scheme we choose a complex parameterization of this term in order to include both charge and spin dependent terms. The argument of the exponential function in Eq. 1 reads:
(2) 
where () is the number (position) of electrons and () is the number (position) of atoms. Both terms on righthand side are expanded over the same localized gaussian basis set .
The first term on the righthand side of Eq. 2 is the electronnucleus term:
(3) 
The term explicitly depends on two electron positions and it accounts for both charge and spin fluctuations. It has the following form:
(4) 
where the homogeneous term fulfills the electronelectron cusp conditions for unlikespin particles; the second term in the righthand side of Eq. 4 describes chargecharge interactions only, while the last term explicitly includes spin correlations and it is therefore suitable to describe spinspin fluctuations. We verified that the inclusion of spinspin term yields a much ( eV/Fe) lower energy than considering the chargecharge term only. This Jastrow has therefore been used for all QMC(opt) calculations presented in the main manuscript.
The matrix elements , , are taken as variational parameters and optimized during the energy minimization procedure.
The full wave function optimization requires several steps. At first the variational parameters of the Jastrow factor are optimized while the Slater determinant is kept fixed at DFT(LDA) level. Afterwards the Slater determinant orbitals parameters (both the contracted orbitals coefficients and the exponents of the gaussian basis functions) are included in the minimization procedure. All the variational parameters of the JastrowSlater wavefunction are optimized by means of the stochastic reconfiguration method with Hessian accelerator, also known as “linear method” Sorella (2001, 2005); Umrigar et al. (2007).
Thanks to the adjoint algorithmic differentiation technique introduced in Ref. Sorella and Capriotti, 2010, energy derivatives with respect to ionic positions and lattice cell parameters can be performed in an efficient way (computational scaling as for total energy evaluation), by treating those quantities on the same footing as the other wavefunction variational parameters. This opens the way to obtain relaxed structures of periodic systems within a fully ab initio variational Monte Carlo approach, as presented in the main manuscript.
The energy gradients with respect to both wavefunction parameters and nuclear positions have always finite variance during the minimization procedure thanks to an efficient reweighting scheme Attaccalite and Sorella (2008); Zen et al. (2013) which cures the divergence of variance in the proximity of the nodes.
In order to compute accurate groundstate energies, once the trial wavefunction is fully optimized, we use a particular implementation of the FNDMC, the socalled Lattice Regularized diffusion Monte Carlo (LRDMC), introduced in Ref Casula et al., 2005. This method provides a solution to the wellknown issue of standard DMC when nonlocal pseudopotential are present as for our FeSe calculations. The problem is solved by introducing an effective Hamiltonian whose Laplacian (kinetic energy) term is evaluated on a discrete lattice. Within LRDMC, the total energy is a variational upper bound to the true ground state energy also in presence of a nonlocal pseudopotential. The tradeoff is represented by the Laplacian discretization error which is assessed in the previous section.
a.3 Bulk moduli and magnetic moments
We report in this section the complete set of results concerning the bulk moduli (Table 2) and the Fe magnetic moments (Table 3) in FeSe. QMC calculations are compared with PBE and PBE0 outcome.
source  ordering  Bulk Mod. (GPa) 

Margadonna Margadonna et al. (2009)  –  30.7(1.1) 
Millican Millican et al. (2009)  –  31 
Kumar Kumar et al. (2010)  –  30.9(3) 
DMC(PBE0)  bicollinear  26.4(8) 
DMC(PBE0)  checkerboard  43.1(8) 
DMC(PBE0)  collinear  31.2(7) 
DMC(PBE0)  collinear, flip 1  27(1) 
DMC(PBE0)  ferromagnetic  43(1) 
PBE0  bicollinear  27.6(1) 
PBE0  checkerboard  40(1) 
PBE0  collinear  29.7(2) 
PBE0  collinear, flip 1  31.8(4) 
PBE0  ferromagnetic  28.3(4) 
PBE  bicollinear  7.03(3) 
PBE  checkerboard  8.9(2) 
PBE  collinear  7.3(1) 
PBE  staggered dimer  8.8(1) 
Source  Magnetic Ord.  Mag. Mom. () 

DMC(PBE0)  bicollinear  3.518 
DMC(PBE0)  checkerboard  3.500 
DMC(PBE0)  collinear  3.443 
DMC(opt)  checkerboard  3.134(6) 
DMC(opt)  collinear  3.014(8) 
PBE0  bicollinear  3.492 
PBE0  checkerboard  3.473 
PBE0  collinear  3.433 
PBE  bicollinear  2.653 
PBE  checkerboard  2.332 
PBE  collinear  2.561 
PBE  stag. dimer  2.543 
References
 Y. Kamihara, H. Hiramatsu, M. Hirano, R. Kawamura, H. Yanagi, T. Kamiya, and H. Hosono, Journal of the American Chemical Society 128, 10012 (2006), ISSN 00027863, 15205126, URL http://pubs.acs.org/doi/abs/10.1021/ja063355c.
 M. Norman, Physics 1 (2008), ISSN 19432879, URL http://link.aps.org/doi/10.1103/Physics.1.21.
 L. de’ Medici, G. Giovannetti, and M. Capone, Phys. Rev. Lett. 112, 177001 (2014), URL http://link.aps.org/doi/10.1103/PhysRevLett.112.177001.
 D. J. Scalapino, Reviews of Modern Physics 84, 1383 (2012), URL http://link.aps.org/doi/10.1103/RevModPhys.84.1383.
 J. Wen, G. Xu, G. Gu, J. M. Tranquada, and R. J. Birgeneau, Reports on Progress in Physics 74, 124503 (2011), ISSN 00344885, 13616633, URL http://stacks.iop.org/00344885/74/i=12/a=124503?key=crossref.7793474e4a4f48d2ddbe92eb3c906f37.
 P. Dai, Reviews of Modern Physics 87, 855 (2015), URL http://link.aps.org/doi/10.1103/RevModPhys.87.855.
 P. Dai, J. Hu, and E. Dagotto, Nature Physics 8, 709 (2012), ISSN 17452473, 17452481, URL http://www.nature.com/doifinder/10.1038/nphys2438.
 T. Imai, K. Ahilan, F. Ning, T. McQueen, and R. Cava, Physical Review Letters 102 (2009), ISSN 00319007, 10797114, URL http://link.aps.org/doi/10.1103/PhysRevLett.102.177005.
 Y. Mizuguchi, F. Tomioka, S. Tsuda, T. Yamaguchi, and Y. Takano, Applied Physics Letters 93, 152505 (2008), ISSN 00036951, URL http://link.aip.org/link/APPLAB/v93/i15/p152505/s1&Agg=doi.
 S. Margadonna, Y. Takabayashi, Y. Ohishi, Y. Mizuguchi, Y. Takano, T. Kagayama, T. Nakagawa, M. Takata, and K. Prassides, Physical Review B 80, 064506 (2009), URL http://link.aps.org/doi/10.1103/PhysRevB.80.064506.
 S. Medvedev, T. M. McQueen, I. A. Troyan, T. Palasyuk, M. I. Eremets, R. J. Cava, S. Naghavi, F. Casper, V. Ksenofontov, G. Wortmann, et al., Nature Materials 8, 630 (2009).
 M. C. Rahn, R. A. Ewings, S. J. Sedlmaier, S. J. Clarke, and A. T. Boothroyd, Phys. Rev. B 91, 180501 (2015a).
 J. Maletz, V. B. Zabolotnyy, D. V. Evtushinsky, S. Thirupathaiah, A. U. B. Wolter, L. Harnagea, A. N. Yaresko, A. N. Vasiliev, D. A. Chareev, A. E. Böhmer, et al., Phys. Rev. B 89, 220506 (2014a), URL http://link.aps.org/doi/10.1103/PhysRevB.89.220506.
 A. Audouard, F. Duc, L. Drigo, P. Toulemonde, S. Karlsson, P. Strobel, and A. Sulpice, EPL (Europhysics Letters) 109, 27003 (2015), URL http://stacks.iop.org/02955075/109/i=2/a=27003.
 F. Caglieris, F. Ricci, G. Lamura, A. Martinelli, A. Palenzona, I. Pallecchi, A. Sala, G. Profeta, and M. Putti, Sci. Tech. Adv. Mater. 13, 054402 (2012).
 A. Liebsch and H. Ishida, Phys. Rev. B 82, 155106 (2010), URL http://link.aps.org/doi/10.1103/PhysRevB.82.155106.
 Yin, Z. P., Haule, K., and Kotliar, G., Nature Materials 10, 932â935 (2011).
 M. Aichhorn, S. Biermann, T. Miyake, A. Georges, and M. Imada, Physical Review B 82, 064504 (2010), URL http://link.aps.org/doi/10.1103/PhysRevB.82.064504.
 L. Craco and S. Leoni, EPL (Europhysics Letters) 92, 67003 (2010), URL http://stacks.iop.org/02955075/92/i=6/a=67003.
 L. Craco, M. S. Laad, and S. Leoni, J. Phys.: Conf. Ser. 487, 012017 (2014), URL http://stacks.iop.org/17426596/487/i=1/a=012017.
 I. Leonov, S. L. Skornyakov, V. I. Anisimov, and D. Vollhardt, Phys. Rev. Lett. 115, 106402 (2015), URL http://link.aps.org/doi/10.1103/PhysRevLett.115.106402.
 J. M. Tomczak, M. van Schilfgaarde, and G. Kotliar, Phys. Rev. Lett. 109, 237010 (2012), URL http://link.aps.org/doi/10.1103/PhysRevLett.109.237010.
 M. Hirayama, T. Miyake, and M. Imada, Phys. Rev. B 87, 195144 (2013), URL http://link.aps.org/doi/10.1103/PhysRevB.87.195144.
 M. H., T. M., T. Miyake, and M. Imada, J. Phys. Soc. Jap. 84, 093703 (2015).
 H. Zheng and L. K. Wagner, Physical Review Letters 114, 176401 (2015), URL http://link.aps.org/doi/10.1103/PhysRevLett.114.176401.
 L. K. Wagner, Physical Review B 92, 161116 (2015), URL http://link.aps.org/doi/10.1103/PhysRevB.92.161116.
 K. Foyevtsova, J. T. Krogel, J. Kim, P. R. C. Kent, E. Dagotto, and F. A. Reboredo, Phys. Rev. X 4, 031003 (2014), URL http://link.aps.org/doi/10.1103/PhysRevX.4.031003.
 N. Devaux, M. Casula, F. Decremps, and S. Sorella, Phys. Rev. B 91, 081101 (2015).
 M. Casula and S. Sorella, Phys. Rev. B 88, 155125 (2013), URL http://link.aps.org/doi/10.1103/PhysRevB.88.155125.
 W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Reviews of Modern Physics 73, 33 (2001), URL http://link.aps.org/doi/10.1103/RevModPhys.73.33.
 J. B. Anderson, J. Chem. Phys. 65 (1976).
 D. M. Ceperley and B. J. Alder, J. Chem. Phys. 81 (1984).
 J. KolorenÄ, S. Hu, and L. Mitas, Physical Review B 82, 115108 (2010), URL http://link.aps.org/doi/10.1103/PhysRevB.82.115108.
 L. K. Wagner, M. Bajdich, and L. Mitas, Journal of Computational Physics 228, 3390 (2009), ISSN 00219991, URL http://www.sciencedirect.com/science/article/pii/S0021999109000424.
 A. E. C. M. Z.W. B. C. S. C. L. M. M. F. M. D. L. P. P. D. Y. N. M. C. M. R. B. K. R. Dovesi, R. Orlando, Int. J. Quantum Chem. 114, 1287 (2014).
 S. Sorella, Turborvb: Quantum monte carlo software for electronic structure calculations, http://people.sissa.it/~sorella/web/index.html.
 M. Burkatzki, C. Filippi, and M. Dolg, The Journal of Chemical Physics 126, 234105 (2007), ISSN 00219606, URL http://jcp.aip.org/resource/1/jcpsa6/v126/i23/p234105_s1.
 M. Burkatzki, C. Filippi, and M. Dolg, The Journal of Chemical Physics 129, 164115 (2008), ISSN 00219606, URL http://jcp.aip.org/resource/1/jcpsa6/v129/i16/p164115_s1?isAuthorized=no.
 S. Sorella and L. Capriotti, J. Chem. Phys. 133, 234111 (2010).
 C. Lin, F. H. Zong, and D. M. Ceperley, Phys. Rev. E 64, 016702 (2001), URL http://link.aps.org/doi/10.1103/PhysRevE.64.016702.
 R. S. Kumar, Y. Zhang, S. Sinogeikin, Y. Xiao, S. Kumar, P. Chow, A. L. Cornelius, and C. Chen, The Journal of Physical Chemistry B 114, 12597 (2010), eprint http://dx.doi.org/10.1021/jp1060446, URL http://dx.doi.org/10.1021/jp1060446.
 J. N. Millican, D. Phelan, E. L. Thomas, J. B. LeÃ£o, and E. Carpenter, Solid State Communications 149, 707 (2009), ISSN 00381098, URL http://www.sciencedirect.com/science/article/pii/S0038109809000854.
 M. C. Rahn, R. A. Ewings, S. J. Sedlmaier, S. J. Clarke, and A. T. Boothroyd, Phys. Rev. B 91, 180501 (2015b), URL http://link.aps.org/doi/10.1103/PhysRevB.91.180501.
 H.Y. Cao, S. Chen, H. Xiang, and X.G. Gong, Phys. Rev. B 91, 020504 (2015).
 J. K. Glasbrenner, I. I. Mazin, H. O. Jeschke, P. J. Hirschfeld, R. M. Fernandes, and R. Valentí, Nature Physics 1, 953 (2015).
 C. Tresca, F. Ricci, and G. Profeta, 2D Materials 2, 015001 (2015).
 W. Wu, Solid State Communications 161, 23 (2013), ISSN 00381098, URL http://www.sciencedirect.com/science/article/pii/S0038109813001087.
 P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al., Journal of Physics: Condensed Matter 21, 395502 (19pp) (2009), URL http://www.quantumespresso.org.
 D. Louca, K. Horigane, A. Llobet, R. Arita, S. Ji, N. Katayama, S. Konbu, K. Nakamura, T.Y. Koo, P. Tong, et al., Phys. Rev. B 81, 134524 (2010).
 A. Benali, L. Shulenburger, N. A. Romero, J. Kim, and O. A. von Lilienfeld, J. Chem. Theory Comp. 10, 3417 (2014).
 S. Sorella, M. Casula, and D. Rocca, J. Chem. Phys. 127, 014105 (2007).
 H. Okabe, N. Takeshita, K. Horigane, T. Muranaka, and J. Akimitsu, Phys. Rev. B (R) 81, 205119 (2010).
 H. Anton and P. C. Schmidt, Intermetallics 5, 449 (1997), ISSN 09669795, URL http://www.sciencedirect.com/science/article/pii/S0966979597000174.
 J. Lischner, T. Bazhirov, A. H. MacDonald, M. L. Cohen, and S. G. Louie, Phys. Rev. B 91, 020502 (2015).
 T. Shimojima, Y. Suzuki, T. Sonobe, A. Nakamura, M. Sakano, J. Omachi, K. Yoshioka, M. KuwataGonokami, K. Ono, H. Kumigashira, et al., 90, 121111 (????), URL http://link.aps.org/doi/10.1103/PhysRevB.90.121111.
 J. Maletz, V. B. Zabolotnyy, D. V. Evtushinsky, S. Thirupathaiah, A. U. B. Wolter, L. Harnagea, A. N. Yaresko, A. N. Vasiliev, D. A. Chareev, A. E. Böhmer, et al., Phys. Rev. B 89, 220506 (2014b), URL http://link.aps.org/doi/10.1103/PhysRevB.89.220506.
 Y. Han, W. Y. Li, L. X. Cao, X. Y. Wang, B. Xu, B. R. Zhao, Y. Q. Guo, and J. L. Yang, Physical Review Letters 104, 017003 (2010), URL http://link.aps.org/doi/10.1103/PhysRevLett.104.017003.
 Y. Wang, T. Berlijn, P. J. Hirschfeld, D. J. Scalapino, and T. A. Maier, Phys. Rev. Lett. 114, 107002 (2015), URL http://link.aps.org/doi/10.1103/PhysRevLett.114.107002.
 L. K. Wagner and P. Abbamonte, Physical Review B 90, 125129 (2014), ISSN 10980121, 1550235X, URL http://link.aps.org/doi/10.1103/PhysRevB.90.125129.
 H. Kwee, S. Zhang, and H. Krakauer, Phys. Rev. Lett. 100, 126404 (2008).
 S. Sorella, M. Casula, L. Spanu, and A. Dal Corso, Phys. Rev. B 83, 075119 (2011).
 S. Sorella, N. Devaux, M. Dagrada, G. Mazzola, and M. Casula, to be published in The Journal of Chemical Physics (2015).
 S. Sorella, Phys. Rev. B 64, 024512 (2001).
 S. Sorella, Phys. Rev. B 71, 241103 (2005).
 C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Hennig, Phys. Rev. Lett. 98, 110201 (2007).
 C. Attaccalite and S. Sorella, Phys. Rev. Lett. 100, 114501 (2008).
 A. Zen, Y. Luo, S. n, and L. Guidoni, J. Chem. Theory Comp. 9, 4332 (2013).
 M. Casula, C. Filippi, and S. Sorella, Phys. Rev. Lett. 95, 100201 (2005).