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 QSGW 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
  • Rigid Spin approximation
• Introduction
• Getting Started
• Nonmagnetic case
• Magnetic case
• Constraining the magnetic moment in antiferrmagnetic FeSe
  • Compare interacting and bare susceptibility
  • Effect of Constraint on Total Energy and Forces

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 QSGW 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 QSGW, 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_j the conjugate magnetic moment to ᗷ_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: $$\delta B_{\mathrm{tot}}^{\alpha}(\mathbf{r}') = \delta B_{\mathrm{ext}}^{\alpha}(\mathbf{r}') +\sum_{\beta} \frac{\delta B_{\mathrm{xc}}^{\alpha}(\mathbf{r}')}{\delta M^\beta(\mathbf{r})} \delta M^\beta(\mathbf{r})$$

where $\alpha,\beta$ are Cartesian indices.

The perturbed spin density $\delta M^{\alpha}$ can be written either in terms of full spin susceptibility $\chi$ or bare (non-interacting) spin susceptibility $\chi^{0}$. $$\delta M^{\alpha} = \sum_{\beta} \chi^{\alpha \beta} \delta B^{\beta,\mathrm{ext}}$$ $$\delta M^{\alpha} = \sum_{\beta} {(\chi^0)}^{\alpha \beta} \delta B^{\beta,\mathrm{tot}}$$

Equating these two expressions we obtain $$\sum_{\beta} \chi^{\alpha \beta} \delta B^{\beta,\mathrm{ext}} = \sum_{\beta} {(\chi^0)}^{\alpha \beta} \large[ \delta B^{\beta,{\mathrm{ext}}} +\sum_{\gamma} \frac{\delta B^{\beta,{\mathrm{xc}}}}{\delta M^\gamma} \sum_{\delta} \chi^{\gamma \delta} \delta B^{\delta,\mathrm{ext}} \large]$$

This can be written in the compact form as $$\chi=\chi^{0}+\chi^{0} \, I^{\mathrm{xc}}$$

or equivalently as $$\chi^{-1} = \left( \chi^0 \right)^{-1} - {I^{\mathrm{xc}}}$$

where we have defined the Stoner matrix to be $$I^{\beta \gamma}= \frac{\delta B^{\beta,\mathrm{xc}}}{\delta M^\gamma} .$$

If we assume $I$ 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 $$\chi(\mathbf{r},\mathbf{r}',\omega)= \chi^0(\mathbf{r},\mathbf{r}',\omega) + \int d\mathbf{r}'' \chi^0(\mathbf{r},\mathbf{r}'',\omega) I^{\mathrm{xc}}(\mathbf{r}'') \chi(\mathbf{r}'',\mathbf{r}',\omega)$$

In collinear spin structures, transverse and longitudinal spin are decoupled. In such a case, matching component by component, we obtain $$I^{\beta \gamma}= \frac{\delta B^{\beta,{\mathrm{xc}}}} {\delta M^\gamma}= \frac{B^{\mathrm{xc}}}{M} \left(\delta_{\beta \gamma}-\frac{M^\beta M^\gamma}{M^2} \right)+\frac{M^\beta M^\gamma}{M^2} \frac{\partial B^{\mathrm{xc}}}{\partial M}$$

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: $$\mathbf{B}^{\mathrm{xc}}=-2 \frac{\mathbf{M}}{M} \frac{\partial}{\partial M}[n \epsilon^{\mathrm{xc}}]= \frac{\mathbf{M}}{M} \, B^{\mathrm{xc}}$$

where $\epsilon^{\mathrm{xc}}$ 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 $$I^{z,{\mathrm{xc}}}=\frac{\partial B^{\mathrm{xc}}}{\partial M} \qquad \qquad I^{+-,{\mathrm{xc}}}=\frac{B^\mathrm{xc}}{M}$$

Similar expressions have been derived for the transeverse susceptibility in a QSGW framework. Most works apply to the transverse susceptibility (rotations of $\mathbf{M}$). 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 $\mathcal{B}$ $${M}(\mathbf{r}) = \sum_{i} {m}_i \mathcal{B}_{i}(\mathbf{r})$$ $$\chi(\mathbf{r},\mathbf{r}',\omega) = \sum_{ij} \chi_{ij} (\omega) \mathcal{B}_{i}(\mathbf{r}) \mathcal{B}_{j}(\mathbf{r}')$$ $$I(\mathbf{r},\mathbf{r}',\omega) = \sum_{ij} I_{ij} (\omega) \mathcal{B}_{i}(\mathbf{r}) \mathcal{B}_{j}(\mathbf{r}')$$

In rigid spin approximation, we approximate the full basis $\mathcal{B}_{i}(\mathbf{r})$ with a single function per magnetic site, namely a normalized function proportional to the magnetization $M_{i}(\mathbf{r})$. Thus $M_{i}(\mathbf{r})$ 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 $I$ is static and a local function of r, it is obtained from $$I_{ii}^{z,{\mathrm{xc}}} = \left< \mathcal{B}_i \right| \frac{\partial B_i^{\mathrm{xc}}}{\partial M_i} \left| \mathcal{B}_i \right> \qquad \qquad I_{ii}^{+-,{\mathrm{xc}}} = \left< \mathcal{B}_i \right| \frac{B_i^{\mathrm{xc}}}{M_i} \left| \mathcal{B}_i \right>, \qquad \qquad \mathcal{B_i}(\mathbf{r}) = M_{i}(\mathbf{r})/ \smallint_R d\mathbf{r} M_{i}(\mathbf{r})$$

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 QSGW, 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 QSGW 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 QSGW 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.

1. 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.

2. 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.

3. 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.

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

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 &

grep RELAX out.nm.relax

 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

lmchk ctrl.fese --rpos=pos

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 QSGW 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 QSGW 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 $(\chi^{0})^{-1}$, and for a general case it is time consuming as 2N band passes are required (N being the number of independent constraints). It does this by taking derivatives dM_i / dB_j via finite difference, by varying 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 $(\chi^{0})^{-1}$; 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 $M$ 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 $\chi^0$ 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 $\chi_0^{-1}$ with

h5dump -d  /chi0inv chi0.fese.h5

It should be about −0.074. You can also estimate this number from the ratio [ΔM/ΔB]⁻¹ 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]⁻¹ = [(2.3522-2.025)/-0.021]⁻¹ = −0.064. That ΔM/ΔB and δM/δB differ by 15% is a measure of nonlinearity: $\chi_0$ 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 χ⁻¹ = [(2.3522-2.025)/-0.011]⁻¹ = -0.034.

If we define an effective Stoner parameter through $\chi^{-1} = \chi_{0}^{-1} + I$, we see that $I\,{\sim}\,0.4\,\mathrm{Ry}\,{\sim}\,0.5\,\mathrm{eV}$. This is not the way Stoner parameters are usually defined. In this context, χ⁻¹ should be a large matrix capturing many neighbors, and should be large enough that χ⁻¹ is small enough for the largest connecting vectors. To a good approximation $I$ is a site-diagonal matrix (it is exactly so in density-functional theory). If a sizeable supercell were built and $I$ inferred approximately in the manner shown here, it should be larger. Typically $I$ is of order 1  eV in 3d 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 QSGW), 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 QSGW is larger than PBE).