How to constrain magnetic moments in VASP

Some ab initio calculations necessitate constraining ion magnetic moments. For instance when determining the direction of anisotropy (in-plane, out-of-plane, etc.), or doing a 4-state mapping to calculate magnetic anisotropies (anisotropic J, DMI, single-ion anisotropy, etc.)¹I plan on writing a tutorial on how to use the 4-state mapping in the future.. In this tutorial, I will summarize what I learnt from the VASP documentation, and my personal experience, to hopefully give you a head start on your own calculations.

A naive starting point would be to use the MAGMOM tag to define the initial magnetic moments per ion. However, the final system’s ions may deviate from this initial configuration to achieve the lowest energy. To counter this, one could constrain the magnetic moment in a noncollinear calculation (see here for a collinear example; thanks to Liang in the comments below). Recall that a noncollinear calculation uses the vasp_ncl executable instead of vasp_std, and uses the tags LNONCOLLINEAR = .TRUE. and LSORBIT = .TRUE. (to turn on spin-orbit coupling), instead of ISPIN = 2 (collinear spin calculation). To the best of my knowledge, if you want to constrain moments in an otherwise collinear spin calculation (where you care only about spin up +1 or spin down -1, instead of a general vector), you will still need to use the noncollinear settings.

For noncollinear calculations, magnetic moments per ion should be specified as a triplet corresponding to a vector in Cartesian coordinates. For instance, if I have a 2-ion system and my first ion needs (1, 2, 0) , and the second needs (3, 1, 1) , the magnetic moment should be given as MAGMOM = 1 2 0 3 1 1. This is the same format required of the M_CONSTR tag as well (to be discussed later). If you are emulating a collinear spin calculation for a 3-atom system with spins up, down, down, you will use = 0 0 1 0 0 -1 0 0 -1 if you have out-of-plane anisotropy in a 2D material. If you want to quickly determine the direction of anisotropy (in-plane vs. out-of-plane, etc.), a simple test would be to check various configurations like = 1 0 0 1 0 0 1 0 0 (for x, in-plane), = 0 1 0 0 1 0 0 1 0 (for y, in-plane), = 0 0 1 0 0 1 0 0 1 (for z, out-of-plane), = 0.707107 0.707107 0 0.707107 0.707107 0 0.707107 0.707107 0 (for 45 degrees from x-axis in the in-plane xy-plane using vector \frac{1}{\sqrt{2}}(1, 1, 0) ), etc.

To actually constrain magnetic moments, you need to use the I_CONSTRAINED_M tag after setting the noncollinear calculation tags. I_CONSTRAINED_M = 1 constrains the direction only, = 2 constrains the size and direction, and = 4 constrains direction and sign. Note that you need at least VASP version 6.4.0 to use = 4. I could be wrong, but from my experience, for = 1, it seems like VASP does not differentiate between directions along the same parallel line. For example, between 0 0 1 and 0 0 -1. But = 4 differentiates between the two and will constrain exactly along the given direction including sign.

Following is the code I use in my noncollinear calculations. For exemplification, I use Nb_3 Cl_8 , where the Nb ions are taken as magnetic, while the Cl ions are not.

I_CONSTRAINED_M = 4
RWIGS = 1.270 1.111
LAMBDA = 9
MAGMOM = 1 0 0 1 0 0 1 0 0 24*0
M_CONSTR = 1 0 0 1 0 0 1 0 0 24*0

Above, note that I wrote 24*0 instead of writing 0 0 0 eight times for the Cl ions ( 8 \times 3 = 24 ). Also notice that I set MAGMOM and M_CONSTR to be the same. M_CONSTR specifies the desired moments, whereas MAGMOM gives the initial one. In my experience, it is important to set both tags to be the same to get the configuration you want.

Additionally, notice the tags RWIGS and LAMBDA. An annoying thing about these tags is that systematic convergence testing has to be done to determine their best values.

For RWIGS, the Wigner-Seitz radius (in Å), the VASP documentation gives a few suggestions for determining best values (it cannot be defined unambiguously for non-monoatomic systems). However, those suggestions can be time-consuming, and so what I sometimes do is use the RWIGS values given in the POTCAR file. Running grep -e RWIGS POTCAR in this case, I got:

   RWIGS  =    2.400; RWIGS  =    1.270    wigner-seitz radius (au A)
   RWIGS  =    2.810; RWIGS  =    1.487    wigner-seitz radius (au A)

Since RWIGS expects the values in Å, I chose 1.270 and 1.487 for Nb and Cl respectively. Another approach (listed in the documentation), is to use the atoms’ covalent radius as given in datasets such as this one. However, be warned that different sources may give different values for this. What ultimately matters is that the calculated moments make sense.

As a workflow strategy, I would start by using the RWIGS values from POTCAR, and use those to do LAMBDA testing. Then I would use the optimal value of LAMBDA to check a few other values of RWIGS, and keep fine-tuning from there, also making sure the integrated magnetic moment is consistent with expectations (see the paragraph after the next one).

As for LAMBDA, this parameter corresponds to the energy penalty added to the system. To constrain magnetic moments along a specific direction, energy has to be added to force the moments to constrain. The manual for I_CONSTRAINED_M states that the total DFT energy should be arrived at by gradually increasing LAMBDA until the energy penalty contribution becomes very small. You can find the value of this contribution (which I gather is given in units of eV), in the last block at the end of the OSZICAR file, where they give: E_p = …. The actual DFT total energy is the energy in OUTCAR (e.g., energy(sigma->0)) minus this E_p. From experience, this is because the energy in OUTCAR already includes the penalty (I could not find the manual stating this, but some simple tests of varying LAMBDA were able to verify this claim). I would read that manual page very thoroughly to be aware of the many subtleties.

For LAMBDA convergence, what I usually do is test a range of values like [0, 1, 10, 25, 50, 75, 100, 125, 150] and keep funneling into a working range based on values that don’t give me an error. In my experience, values that are too high give an error. Keep decreasing the range, e.g., to [1, 10, 25], and then [1, 3, 5, 7, 10, 13, 15, 17, 20], to [4, 5, 6, 7, 8, 9, 10, 11, 12], etc., until you find an optimal value. The goal is to find LAMBDA that is large enough that E_p is very small, while giving the final magnetic moments in the desired direction. You can verify the latter by inspecting the final occurrences of magnetization (x), magnetization (y) and magnetization (z) in OUTCAR. You need to set LORBIT = 11 for the code to display this information in OUTCAR. These magnetizations should be consistent with the moments defined by the MAGMOM and M_CONSTR tags in INCAR. Note that if you are not constraining the size/magnitude, your magnitudes can differ. For instance, for = 0 0 1 0 0 1 0 0 1 24*0, I got:

 magnetization (z)

# of ion       s       p       d       tot
------------------------------------------
    1        0.002   0.002   0.212   0.217
    2        0.002   0.002   0.213   0.217
    3        0.002   0.002   0.212   0.216
    4       -0.000  -0.001   0.001   0.001
    5       -0.000  -0.001   0.001   0.001
    .....

and ~0 for all ions in magnetization (x) and magnetization (y). The ~0.21 magnitude is reasonable in this case, because we expect a single spin-half to be shared by each Nb atom in this case (for this specific material system). This means each Nb atom should roughly be \frac{1}{2}/3 = 1/6 \approx 0.17 in magnitude. One could make this ~0.21 come very close to 0.17 by changing RWIGS. In this case, RWIGS would have to be made smaller so that the total area of integration becomes smaller and doesn’t capture overlaps. But please be careful following this reasoning, as there are several other subtleties involving M_int and MW_int in the I_CONSTRAINED_M page. However, if you used I_CONSTRAINED_M = 2 (to constrain size too), you should expect the final moment size to be close to the ones you defined in MAGMOM and M_CONSTR.

Finally, you now have all the ingredients you need to do a constrained magnetic moment calculation. Please note that there might be a lot of issues with convergence, and the internet has many potential solutions. It can be a tricky problem, as exemplified by one page on the VASP manual stating “praying” as a potential workaround. One strategy I found quite helpful when testing out different moment directions systematically, was to first do a spin-unpolarized calculation ISPIN = 1, and then use its CHGCAR along with ICHARG = 1 to start from that charge density. It seems to help with convergence, speed, and accuracy. But, before using this CHGCAR, please make sure your unpolarized calculation detects the expected symmetry using the OUTCAR file (see this tutorial I wrote here a while ago).

Finally, and in a sudden-jump-fashion, I wanted to share a simple Python script I use to help me write out the MAGMOM and M_CONSTR tags for my INCAR files. For large systems, they can be cumbersome (and dangerous!) to write out manually. So, hopefully the following will be of use!
I exemplify using a 4 \times 4 supercell for Nb_3 Cl_8 , which has 48 Nb ions and 128 Cl. Below, I give the first 3 Nb ions 1 0 0, the next 3 Nb ions -1 0 0, the remaining 42 Nb ions 0 1 0, and all the Cl ions 0 0 0:

def generate_string(a, b, c, n1, n2, n3):
    result = []
    result.extend(a * n1)
    result.extend(b * n2)
    result.extend(c * n3)
    return "M_CONSTR = " + " ".join(map(str, result)) + " 384*0" # hardcoding 128 * 3 = 384 for Cl

n1 = 3; n2 = 3; n3 = 42 # number of ions

a = [1, 0, 0]; b = [-1, 0, 0]; c = [0, 1, 0] # desired magnetic moments for Nb
print(generate_string(a, b, c, n1, n2, n3))