# Constraining Local Magnetic Moments

Performing self-consistent calculatons with local magnetic moments constrained to given values can be a very useful tool. Among other things it can yield information about the longitudinal magnetic susceptibility.

Questaal can do this both at the DFT level and at the QS*GW* level. Both are accomplished via the program *lmf*. There are two kinds of ways to constrain the magnetic moment: constrain the global moment, or constrain individual local moments. Both are done through applying effective magnetic fields. The global scheme, first implemented by Moruzzi et al, can be found in many DFT codes. A global Zeeman field is added, whose value is determined to satisfy a contraint of a fixed global magnetic moment.

The second is less common, but more flexible. A site-local Zeeman field is applied to satisfy one or more constraints. One natural outcome of this approach is that the spatially resolved static longitudinal susceptibility can be resolved.

In this tutorial we apply the latter scheme to FeSe. This material system is of high current interest in large part because of its interesting many body effects (superconducting properties, nematicity, quantum critical point). The strong interplay between magnetism and structure makes this an ideal system to uncover some key physics by constraining local magnetic moments.

### Table of Contents

- Table of Contents
- Preliminaries
- Command summary
- Theory
- Introduction
- Getting Started
- Nonmagnetic case
- Magnetic case
- Constraining the magnetic moment in antiferrmagnetic FeSe

### Preliminaries

This tutorial assumes you have a familiarity with how to run the *lmf* code at the DFT level, as explained in this tutorial.

It also assumes you have compiled the Questaal executables and that they are in your path.

Note that everywhere you see the *lmf* command, you can run the same instruction in MPI mode, .e.g. a standard MPI launch with 8 processors would look replace `lmf`

with `mpirun -n 8 lmf`

. This tutorial uses `$mn`

as an abstraction for the MPI launcher.

This tutorial will do the following:

- Carry out a self-consistent, nonmagnetic calculation of FeSe at the experimental structure to establish that the force on Se is large
- Do the same calculation for antiferromagnetic FeSe, showing that the force is much improved
- Do another calculation for antiferromagnetic FeSe this time contraining the local Fe moment to an average moment obtained from a QS
*GW*calculation of a paramagnetic structure. With this constraint, the forces become small.

### Command summary

Start a fresh directory with *init.fese* in it. *ctrl* file:

```
blm --pbe --brief --mag --upcase --molstat --rsham~el=-.4,-1.2~npr=-50 --autobas~ehmx=-.4~eloc=-3.5~loc=5 --mix~nit=20 --gw~q0ps~wemax=0~emax=3~gcutb=3.0~gcutx=2.4~nk=6 --nk=8 ctrl.fese
lmchk ctrl.fese
```

Free nonmagnetic atomic density and basis set parameters

```
echo '-vnsp=1 --v8 --tbeh -vrsham=2 -vbeta=.1 -vnit=30' > switches-for-lm
lmfa ctrl.fese `cat switches-for-lm`
cp basp0.fese basp.fese
```

Self-consistent nonmagnetic PBE density

```
$mn lmf ctrl.fese `cat switches-for-lm` > out.nm
```

Free magnetic atomic density and basis set parameters

```
echo '-vnsp=2 --v8 --tbeh -vrsham=2 -vbeta=.1 -vnit=30' > switches-for-lm
lmfa ctrl.fese `cat switches-for-lm`
rm -f mixm.fese rst.fese
$mn lmf ctrl.fese `cat switches-for-lm` > out.afm
```

Constrained magnetic density

```
sed -i.bak 's/HAM/HAM BFIELD=2/' ctrl.fese
mcx -1:2 -e4 i 0 0 0 | sed s/.000000//g > bfield.fese
echo '-vnsp=2 --v8 --tbeh -vrsham=2 -vbeta=.1 -vnit=30 --fixmm~fullrho~ib=1~eq=-2~mom=2.355' > switches-for-lm
rm -f mixm.fese rst.fese
cp rst.pbe.afm rst.fese
$mn lmf ctrl.fese `cat switches-for-lm` > out.cnstm.afm
```

### Theory

Here we want to obtain a self-consistent density within density-functional theory (DFT), or the Quasiparticle Self-Consistent *GW*, or QS*GW*, Approximation), subject to the constraint that the magnetic moment (integrated magnetization density) at one or more site is of a prescribed value.

In DFT there is no knob to apply this constraint directly; however, the magnetization and magnetic fields are conjugate pairs entering into the total energy, analogous to the charge and electrostatic potential. We can apply an external magnetic field in such a way as to satisfy the constraints. Here we are concerned with only the total magnetic moment within an augmentation sphere (local moment), so **M**(**r**) and **ᗷ**(**r**) get discretized into points *i* at nucleii where constraints are imposed. Even though **M** and **ᗷ** are vectors, at present the band code *lmf* only has the ability to vary the *z* component of **M**. Thus this tutorial treats *ᗷ _{i}* as a scalar (Zeeman) field and

*M*the conjugate magnetic moment to

_{j}*ᗷ*.

_{i}(Note: the following was adapted from Liqin Ke’s PhD thesis).

The magnetic susceptibility is the system’s magnetic response to an external magnetic field:

where are Cartesian indices.

The perturbed spin density can be written either in terms of full spin susceptibility or bare (non-interacting) spin susceptibility .

Equating these two expressions we obtain

This can be written in the compact form as

or equivalently as

where we have defined the Stoner matrix to be

If we assume to a local function of **r** and independent of frequency (as it is within the local-density approximation), this expression can be written more fully as

In collinear spin structures, transverse and longitudinal spin are decoupled. In such a case, matching component by component, we obtain

where the first term is purely transverse and the second one is purely longitudinal with respect to the local magnetization vector. The transverse term relies on the following relation within the LSDA:

where is the exchange-correlation energy per particle.

Taking *z* to be the spin quantization axis we arrive at Stoner parameters for the transverse and longitudinal parts

Similar expressions have been derived for the transeverse susceptibility in a QS*GW* framework. Most works apply to the transverse susceptibility (rotations of ). Here we are concerned with the longitudinal susceptibility (variations in amplitude). Unlike the transverse susceptibility, the longitudinal spin susceptibility couples to the charge susceptibility. This effect is crucially important in FeSe, as this tutorial shows.

##### Rigid Spin approximation

For a practical scheme within *lmf*, the preceding expressions need to be discretized. To see how to accomplish this, represent the relevant quantities in a basis set

In rigid spin approximation, we approximate the full basis with a single function per magnetic site, namely a normalized function proportional to the magnetization . Thus can either rotate (transverse case) or change amplitude (longitudinal case), subject to the constraint that the shape remains fixed. This is how we can discretize the full description. Assuming is static and a local function of **r**, it is obtained from

For constraining the moments, we require the longitudinal susceptibility.

### Introduction

As mentioned earlier, FeSe is of high current interest because of its many kinds of novel many-body effects. Because FeSe is paramagnetic, it is is commonly treated in density-functional theory as being nonmagnetic. It has long been known, however, that a nonmagnetic description of FeSe yields a poor description of the FeSe bond length. Fe resides in the basal plane in a checkerboard pattern, and the Se sites above the plane between Fe sites. The height of the Se above the plane (which controls the Fe-Se bond length and Fe-Se-Fe bond angle) crucially controls most key many-body effects in FeSe, for example the instability to superconductivity. Perhaps the careful measurement of this height is by the Cava group, which puts it 1.463 Å. How the Se height affects superconductivity is explained in some detail in this npj paper. See also this paper.

Even since FeSe was first discovered to be an unconventional superconductor (meaning superconductivity is not mediated by phonons), it has been known that the Se height has been difficult to predict in DFT. More recent work has shows that an antiferromagnetic description does a better job of predicting structure. Thus, so far as structure is concerned, an antiferromagnetic solution is a better proxy for the paramagnetic state than a nonmagnetic solution is.

This tutorial explores these observations, and shows that the local Fe moment has a significant effect on the structure. Further, it is found that within QS*GW*, the local Fe moment is further enhanced. (That is not demonstrated in this tutorial; we will just take the result). If we apply the PBE functional with the moment constrained to the QS*GW* result, we find we obtain an excellent description of the Se height. Thus, the inadequacy of DFT to predict the structure can be traced to its inadequate description of magnetism.

### Getting Started

Copy the contents of the box below to file *init.fese*

```
HEADER FeSe, from Nitsche et al, Inorg. Chem., 2012, 51 (13), p 7370
LATTICE
% const a=3.779
ALAT={a} UNITS=A
PLAT= 1.0000000 0.0000000 0.0000000
0.0000000 1.0000000 0.0000000
0.0000000 0.0000000 1.4583488
SPEC ATOM=Feup Z=26 MMOM=0,0,2
ATOM=Fedn Z=26 MMOM=0,0,-2
SITE
ATOM=Feup X= 0.0000000 0.0000000 0.0000000
ATOM=Fedn X= 0.5000000 0.5000000 0.0000000
ATOM=Se X= 0.0000000 0.5000000 0.2655000
ATOM=Se X= 0.5000000 0.0000000 0.7345000
```

As noted above, many-body effects depend crucially on the height of the Se, and also in a non-negligible way on the *c*/*a* ratio. The *init* file above takes the experimentally determined height from Cava’s group. (**X** specifies the positions in units of the reciprocal lattice vectors, as explained on this page). The tutorial below uses the PBE functional to demonstrate the strong interplay between magnetism and this height.

Make the *ctrl* file (input file) with the following instruction

```
blm --pbe --brief --mag --upcase --molstat --rsham~el=-.4,-1.2~npr=-50 --autobas~ehmx=-.4~eloc=-3.5~loc=5 --mix~nit=20 --gw~q0ps~wemax=0~emax=3~gcutb=3.0~gcutx=2.4~nk=6 --nk=8 ctrl.fese
```

The switches to *blm* are explained in the command line documentation. Note in particular, `--pbe`

sets tags to use the PBE functional, `--molstat`

tags for molecular relaxation, `--rsham`

to use the real-space, tight-binding form of the basis set, `--gw`

sets some tags should this file be used for subsequent QS*GW* calculations (not done in this tutorial), `--nk=8`

sets the number of *k* divisions. A finer mesh would change results slightly.

First, check the structure, and confirm that the Se height is very close to Cava value.

```
lmchk ctrl.fese
```

Look at locations of the Fe and Se nucleii: the Fe lie in the basal plane *z*=0, while the Se lie above, at *z*=0.3872 in units of the lattice constant. To convert to Å, note the lattice constant is `alat`=7.141 a.u., or 3.779 Å. Thus the Se height is 0.3872×7.141×.5292 = 1.463 Å.

### Nonmagnetic case

The PBE functional yields a larger and better magnetic moment than straight LDA theory; that is important here.

Make the free atomic density and basis set parameters

```
echo '-vnsp=1 --v8 --tbeh -vrsham=2 -vbeta=.1 -vnit=30' > switches-for-lm
lmfa ctrl.fese `cat switches-for-lm`
cp basp0.fese basp.fese
```

The first line sets switches, mostly for the band code *lmf*. `-vnsp=1` tells it to perform a nonmagnetic calculation; `–v8 –tbeh -vrsham=2` turn on the tight-binding mode; `-vbeta=.1` sets a relatively low mixing of the output density. (The density has both fast and slow degrees of freedom, especially when magnetism is involved, and it is safer to mix output densities into the trial density rather slowly.)

Make the density self-consistent

```
rm -f mixm.fese rst.fese
$mn lmf ctrl.fese `cat switches-for-lm` > out.nm
```

The first instruction is not necessary unless you are rerunning a prior calculation.

Here are some useful things to check.

That the Harris-Foulkes and Hohenberg-Kohn-Sham energies coincide. You can see them from the last lines of the output file

*out.nm*, and can confirm tye are nearly the same.That the

*k*convergence is adequate. The simplest is to check the size of Bloechl correction to the tetrahedron integration of the band sum. Find this with`grep 'Bloechl' out.nm`

. You should see that the correction is about 1.5 mRy, which is quite good enough for this purpose. A more careful check would be to rerun the calculation with a finer*k*-mesh.Inspect the forces at self-consistency:

`grep 'Maximum Harris force' out.nm`

This tells you the largest force on the system, as it cycles to self-consistency. The table should look like this:

Forces, with eigenvalue correction ib estatic eigval total 1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3 0.00 0.00 654.63 0.00 0.00 -718.36 0.00 0.00 -63.73 4 0.00 0.00 -654.63 0.00 0.00 718.36 0.00 0.00 63.73 Maximum Harris force = 63.7 mRy/au (site 3) Max eval correction = 2.3

63.7 mRy/au is a very large force, which means that the the minimum total energy for this functional corresponds to a different Se height. Note that the force is negative for atom 3, which means it wants to relax towards the Fe in the basal plane.

*lmf* will relax the structure, and it is interesting to see what the equilibrium height is. It is not essential to this tutorial, so you can omit this step.

The *ctrl* file has these lines

% ifdef minx DYN MSTAT[MODE=6 HESS=t XTOL=.001 GTOL=0 STEP=.010 NKILL=0] NIT={minx} % endif

which tell *lmf* to relax the structure according to the forces it generates. `MODE=6` tells *lmf* to use a Broyden scheme; the value of `minx` supplies an upper bound to the number of relaxation steps it will take before giving up.

Do the following

```
rm -f mixm.fese rst.fese hssn.fe
$mn lmf ctrl.fese `cat switches-for-lm` --rs=3 -vminx=10 --wpos=pos > out.nm.relax &
```

You can see how the relaxation proceeded with

```
grep RELAX out.nm.relax
```

You should see

RELAX: completed Broyden iter 1; max shift=-0.0100 |g|=0.0901 RELAX: completed Broyden iter 2; max shift=-0.0100 |g|=0.0659 RELAX: completed Broyden iter 3; max shift=-0.0100 |g|=0.0398 RELAX: completed Broyden iter 4; max shift=-0.0039 |g|=0.0111 RELAX: converged to tolerance; max shift=9.26e-5 |g|=2.72e-4

Also check the maximum force with `grep RELAX out.nm.relax`

. You should see that it is very small at the last step. Relaxed positions are written to file *pos.fese*. To see the height, do

```
lmchk ctrl.fese --rpos=pos
```

You should find that the Se height is 1.336 Å. The many-body properties of this system are profoundly different when the Se height is 1.34 Å as opposed to 1.46 Å.

### Magnetic case

Let’s repeat the calculation, this time allowing (anti)ferromagntrism to develop. Note the *init* file assumes a magnetic moment of ±2 *μ*_{B} on the Fe *d* orbitals. This moment will be included in the generation of the free-atomic density (*lmfa*), and it gives a reasonable starting trial moment for *lmf*. It will determined self-consistently.

Redo the nonmagnetic calculation but this time for a spin polarized case

```
echo '-vnsp=2 --v8 --tbeh -vrsham=2 -vbeta=.1 -vnit=30' > switches-for-lm
lmfa ctrl.fese `cat switches-for-lm`
rm -f mixm.fese rst.fese
$mn lmf ctrl.fese `cat switches-for-lm` > out.afm
cp rst.fese rst.pbe.afm
```

The last step preserves the restart file for later steps in the tutorial.

You can see how the moment develops as it iterates to self-consistency. Do:

```
grep -A 1 'true mm' out.afm
```

At final iteration you should see something similar to

mkrout: Qtrue sm,loc local true mm smooth mm local mm ovlp 1 6.649772 2.393983 4.255789 2.025467 0.479954 1.545512 1.957406 contr. to mm extrapolated for r>rmt: 0.094986 est. true mm = 2.120452 2 6.649773 2.393983 4.255789 -2.025467 -0.479954 -1.545512 -1.957406 contr. to mm extrapolated for r>rmt: -0.094986 est. true mm =-2.120452

The PBE local moment turns out to be very close to ±2 *μ*_{B} for the experimental structure.

Comparing the total energies of the AFM and NM cases, you can see that the AFM case is lower in energy by about 25 mRy. (The total energy can be found at the end of the output file.)

Inspect the forces at self-consistency: `grep 'Maximum Harris force' out.afm`

. You should see the maximum force is about 18 mRy/au — much better than the nonmagnetic structure. Curiously though, this AFM structure has a very small gap (`grep gap out.afm`

), contrary to experiment. Thus the AFM state describes the total energy better, but the energy bands are poorly described.

If you relax this structure, you can observe that the equilibrium Se height is 1.417 Å. This is much closer to the experimental 1.46 Å, but it continues to be too small.

### Constraining the magnetic moment in antiferrmagnetic FeSe

A QS*GW* calculation of the paramagnetic phase predicts the average absolute value of the local moment to be 2.355 *μ*_{B}. (Such is a tour de force calculation, and here we only use the result. An enhancement of the moment in the change AFM→paramagnetic is *not* found with the PBE functional; and in this crucial respect many-body perturbation theory and density-functional theory are quite different.) 2.355 *μ*_{B} is substantially larger than the 2.03 *μ*_{B} PBE predicts for the AFM phase. If we allow 2.355 *μ*_{B} to be more representative of the PM moment, we can perform a PBE calculation subject to the constraint that the local moment matches the mean QS*GW* paramagnetic moment.

Switch `--fixmm~options`

turns on Questaal’s constrained local moments facility; its options are explained in the command line arguments documentation.

In preparation for using this switch, you must do three things.

- Tag
**BFIELD=2**must be added to the**HAM**category. Add this tag with your text editor, or just do`sed -i.bak 's/HAM/HAM BFIELD=2/' ctrl.fese`

Confirm that the

**HAM**category in original*ctrl.fese*now reads**HAM BFIELD=2**. - You must create file
*bfield.fese*, specifying the site-dependent external field. Even though the field will be determined by the constraint in this context,*lmf*requires this file be present whenever**BFIELD**is turned on.

Create the file with e.g.`mcx -1:2 -e4 i 0 0 0 | sed s/.000000//g > bfield.fese`

Inspect the file. It uses Questaal’s standard format for ASCII data, with the caveat that

*lmf*reads this array in sparse format. The first column specifies the site, and columns 2-4 specify the*x*-,*y*-, and*z*- components of the field (all zero in this case). At present*lmf*allows only collinear moments along the*z*axis, and the*x*-,*y*- components are just placeholders. The*z*component will be determined in the course of the calculation. - You must specify which sites are to be constrained, and what the magnetic moment is for each site. The most general way is to supply a file
*mmom.ext*for each moment you want to constrain.*mmom.fese*would look like following:

% rows 2 cols 4 1 0 0 2.355 2 0 0 -2.355

*mmom.fese* and *bfield.fese* have the same structure, with columns 1-4 specifying the site, and the *x*-, *y*-, and *z*- components of the moment or field.

Here we have only two moments, and moreover the the up- and down- local moments should be constrained to be equal and opposite. There only one independent degree of freedom. Instead of specifying constraints through file *mmom.fese* as described above, you can alternatively impose constraints by modifiers to `--fixmm`

. Moreover, you can specify that moments on certain sites are linked to other sites. The modifier `~ib=1~eq=-2~mom=2.355`

, accomplishes what we want: it says that site 1 has moment 2.355 and site 2 has the same magnitude but opposite sign — the minus sign flags that the sign should be opposite. (In practice the code does not impose a separate constraint, but assigns the *ᗷ* field of site 2 to be the same magnitude as site 1, but with opposite sign. Since the constraints are specified on the command-line, *mmom.fese* is not needed here.

```
echo '-vnsp=2 --v8 --tbeh -vrsham=2 -vbeta=.1 -vnit=30 --fixmm~fullrho~ib=1~eq=-2~mom=2.355' > switches-for-lm
rm -f mixm.fese rst.fese
cp rst.pbe.afm rst.fese
$mn lmf ctrl.fese `cat switches-for-lm` > out.cnstm.afm
```

The `~fullrho`

modifier tells *lmf* to use the full magnetic moment inside the augmentation sphere. This corresponds to the column ‘true mm’ in the printout of magnetic moments (2.025) in the table shown above for the unconstrained moment. Without this modifier, *lmf* would choose the moment from the ‘ovlp’ column in this table.

Inspect this table with the constraint applied. The last occurence in *out.cnstm.afm* should read something similar to

mkrout: Qtrue sm,loc local true mm smooth mm local mm ovlp 1 6.637305 2.386850 4.250455 2.355607 0.569491 1.786116 2.276871 contr. to mm extrapolated for r>rmt: 0.113828 est. true mm = 2.469435 2 6.637305 2.386850 4.250455 -2.355607 -0.569491 -1.786116 -2.276871 contr. to mm extrapolated for r>rmt: -0.113828 est. true mm =-2.469435

At self-consistency, ‘true mm’ is ±2.3556 *μ*_{B}, equal to the constraint (2.355) imposed within the default tolerance.

Next look at how *lmf* imposed the constraint. Do: `grep cnstm: out.cnstm.afm`

and you should see a long sequence of lines.

The first lines print information about the setup conditions

sucnstm: 0 constraints imported from moms file sucnstm: 1 constraint on moments (m from full density) sucnstm: 1 dependent sites: 2 sucnstm: maximum number of cycles = 3 (grad + 2 lin-resp) mtol=1e-3 btol=1e-4

They tell us that one constraint is imposed, and there is one other magnetic site dependent on it. The last line says that the constraint algorithm will make three passes, or cycles. The first cycle needs to build the inverse susceptibility , and for a general case it is time consuming as 2*N* band passes are required (*N* being the number of independent constraints). It does this by taking derivatives *dM _{i}* /

*dB*via finite difference, by varying

_{j}*B*numerically at each constraining site.

Next appear the lines

cnstm: constrain 1 loc mmom: maximum 3 cycles (grad + 2 lin-resp) Mtol=1e-3 dB=5e-3 dBmx=0.03 cnstm: begin cycle 1 (max 3), grad mode, 1 constrained moment diffm = 0.3300 cnstm: begin iterations for constraint 1 of 1 (channel 1). cnstm: iter 1 cnst 1 cycle 1 m-mtarg -0.32997 new B 0.002500 dB 0.002500 ir -1 (1st point) cnstm: continue search (iter 1, glob 1) for constraint 1 channel 1 cnstm: iter 2 cnst 1 cycle 1 m-mtarg -0.37242 new B -0.002500 dB -0.005000 ir -2 (2nd point) cnstm: continue search (iter 2, glob 2) for constraint 1 channel 1 cnstm: iter 3 cnst 1 cycle 1 m-mtarg -0.28848 new B 0.000000 dB 0.002500 ir 0 (completed) cnstm: end cycle 1 (max 3), grad mode diffm = 0.2885

The first cycle generates ; that is what is revealed in the lines above.

Once in hand, it can be assumed to be held fixed and several additional cycles can be carried out to solve the nonlinear coupled equations (only one equation here). In the present setup the total number of cycles is capped at 3, so it carries two additional cycles in this configuration.

cnstm: band pass for linear response cycle 1 : make M from linear response estimate of B cnstm: begin cycle 3 (max 3), linear response mode, 1 constrained moment diffm = 0.02791 cnstm: band pass for linear response cycle 2 : make M from linear response estimate of B cnstm: begin cycle 4 (max 3), grad mode, 1 constrained moment diffm = 0.002747 cnstm: cycle count reached limit (3) ... ending search cnstm: terminate cycles, deviation of M from target 0.002747 after 3 band passes

You can see the deviation of *M* from the target value drops by one order of magnitude in each cycle.

At this point the algorithm has reached the maximum number of cycles. It gives up searching for *M* and returns to the outer loop that iterates to charge self-consistency. The charge and potential are updated. You should find this line in *out.cnstm.afm* after the lines shown above <pre> mixrho: sought 2 iter, 12586 elts from file mixm; read 0. RMS DQ=8.78e-3 </pre> The change in the density is significant: it indicates that there is a strong cross-coupling between charge and longitudinal spin degrees of freedom.

*lmf* starts a new pass, running through three moment-searching cycles again. Now the starting is much closer so has a much easier time iterating to the target *M*. In the second loop you can see the following in the search for *M*.

... cnstm: end cycle 1 (max 3), grad mode diffm = 0.04042 cnstm: moments converged in 2 cycles, 3 calls: deviation 1.638e-5 from target within tol (1e-3) cnstm: terminate cycles, deviation of M from target 1.638e-5 after 3 band passes

*lmf* exits the inner loop and the density is again updated. The inner cycle is repeated for each outer one. After 15 outer cycles the density is self-consistent, including an appropriate *B* to satisfy the constraint that *M* alignes with the target value

mixrho: sought 2 iter, 12586 elts from file mixm; read 0. RMS DQ=8.78e-3 mixrho: sought 2 iter, 12586 elts from file mixm; read 1. RMS DQ=8.01e-3 last it=8.78e-3 mixrho: sought 2 iter, 12586 elts from file mixm; read 2. RMS DQ=7.24e-3 last it=8.01e-3 ... mixrho: sought 2 iter, 12586 elts from file mixm; read 0. RMS DQ=1.11e-5 last it=1.24e-5

The last printout from the inner cycle reads

cnstm: terminate cycles, deviation of M from target 6.067e-4 after 0 band passes

##### Compare interacting and bare susceptibility

We can infer both the bare and interacting susceptibility for this system. Note that *lmf* generates a file *chi0.fese.h5*. To see what is in the file, do `h5ls chi0.fese.h5`

.

In particular the inverse of the bare susceptibility is stored. In general it is a matrix, but here it is only a 1×1 matrix because there is only one degree of freedom. Read with

```
h5dump -d /chi0inv chi0.fese.h5
```

It should be about −0.074. You can also estimate this number from the ratio [Δ*M*/Δ*B*]^{−1} as determined from the last iteration of the first large cycle. Look for last occurence of `bnow` in *out.cnstm.afm* before the close of the first charge loop. You should find this

bnow = -0.021460850551 mnow = 2.352252638950 mtrg = 2.355000000000

The magnetic moment without constraint is 2.025, as noted previously. Thus [Δ*M*/Δ*B*]^{−1} = [(2.3522-2.025)/-0.021]^{−1} = −0.064. That Δ*M*/Δ*B* and δ*M*/δ*B* differ by 15% is a measure of nonlinearity: is not independent of *M*.

Now dind the last occurence `bnow` after self-consistency has been reached. You should find something like

bnow = -0.011353139979 mnow = 2.355606682064 mtrg = 2.355000000000

so that χ^{−1} = [(2.3522-2.025)/-0.011]^{−1} = -0.034.

If we define an effective Stoner parameter through , we see that . This is not the way Stoner parameters are usually defined. In this context, χ^{−1} should be a large matrix capturing many neighbors, and should be large enough that χ^{−1} is small enough for the largest connecting vectors. To a good approximation is a site-diagonal matrix (it is exactly so in density-functional theory). If a sizeable supercell were built and inferred approximately in the manner shown here, it should be larger. Typically is of order 1 eV in 3*d* transition metals.

##### Effect of Constraint on Total Energy and Forces

You can read the total energy of unconstrained self-consistent from the end of *out.afm*. You should see something similar to

it 27 of 30 ehf= -14811.5605405 ehk= -14811.5605305

For the constrained case (*out.cnst.afm*) you should find

it 15 of 30 ehf= -14811.5571212 ehk= -14811.5571277

As it should, imposing the constraint weakens the binding energy.

However, the PBE functional itself has limited accuracy. In particular we can see how the forces at the experimental geometry compare in the unconstrained and constrained cases.

Reading from the output of the unconstrained case you should see a table like this one

Forces, with eigenvalue correction ib estatic eigval total 1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3 0.00 0.00 670.51 0.00 0.00 -688.48 0.00 0.00 -17.97 4 0.00 0.00 -670.51 0.00 0.00 688.48 0.00 0.00 17.97 Maximum Harris force = 18 mRy/au (site 3) Max eval correction = 0.2

while in the constrained case it should be similar to

Forces, with eigenvalue correction ib estatic eigval total 1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3 0.00 0.00 675.18 0.00 0.00 -678.70 0.00 0.00 -3.51 4 0.00 0.00 -675.18 0.00 0.00 678.70 0.00 0.00 3.51 Maximum Harris force = 3.51 mRy/au (site 4) Max eval correction = 0.9

The forces in the latter case are quite small. This shows that if average magnetic moment matched the paramagnetic case (determined from QS*GW*), the structure is well described. This indicates nonlocal spin fluctuations are essential in determining the structure. (It is also true that the average moment from QS*GW* is larger than PBE).