Constrain the size and direction of the magnetic moments

[liang_zhang]

Dear all,

I am trying to constrain the size and direction of the magnetic moments of O2 molecule to its triplet state. But I found the penalty energy does not decrease (instead increase) when increasing LAMBDA in the range [1,10,25,50]. The results show that it's hard to constrain the size of the magnetic moments?

Here is my INCAR:
SYSTEM = O2
ISTART = 0
ICHARG = 1
GGA=RP
ISMEAR = 0
SIGMA = 0.01
ALGO = Fast
ENCUT = 400
EDIFF = 1E-5
LREAL = Auto
NELM = 200

LNONCOLLINEAR = .TRUE.
MAGMOM=1 0 0 1 0 0
LSORBIT = .TRUE.
LORBIT = 11
I_CONSTRAINED_M = 2
RWIGS = 0.820
LAMBDA = 50
M_CONSTR = 1 0 0 1 0 0
AMIX = 0.2
BMIX = 0.0001
AMIX_MAG = 0.8
BMIX_MAG = 0.0001
LCHARG = .TRUE.

Here is my POSCAR:
O2 molecule in a box
1.00000000000000
7.5000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 8.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 8.9000000000000004
O
2
Selective dynamics
Direct
0.0000000000000000 0.0000000002719517 -0.0016986474597765 T T T
0.0000000000000000 -0.0000000002719517 0.1373615688080896 T T T

The K-point is 1x1x1.

Here is the output of OSZICAR:
E_p = 0.20616E+01 lambda = 0.500E+02
<lVp>= -0.24593E+02
DBL = 0.26655E+02
ion MW_int M_int
1 0.856 -0.000 0.000 1.391 -0.000 0.000
2 0.856 0.000 -0.000 1.391 0.000 -0.000
RMM: 24 -0.164229872423E+01 -0.28394E-05 -0.18771E-06 74 0.841E-03
1 F= -.16422987E+01 E0= -.16422987E+01 d E =-.101177-114 mag= 1.9997 0.0000 -0.0000

E_p = 0.20616E+01 lambda = 0.500E+02
ion lambda*MW_perp
1 -0.71788E+01 -0.28455E-19 0.00000E+00
2 -0.71788E+01 -0.41028E-19 0.13553E-17

Please correct me if there are any errors. Thanks in advance!

Best wishes,
Liang

[fabien_tran1]

Hi,

In the OSZICAR there is " mag= 1.9997 0.0000 -0.0000", i.e., a triplet state is obtained. The discussion at https://www.vasp.at/forum/viewtopic.php?t=19187 about the penalty energy may be relevant.

[liang_zhang]

Hi,

I also tried only constrain the direction of the magnetic moments using I_CONSTRAINED_M = 1 (not 2). The result seems to be reasonable, the E_p is really small:

ion MW_int M_int
1 0.439 -0.000 0.000 0.828 -0.000 0.000
2 0.439 -0.000 0.000 0.828 0.000 0.000
RMM: 13 -0.957636086353E+01 0.10245E-05 -0.19022E-05 73 0.352E-02
1 F= -.95763609E+01 E0= -.95763609E+01 d E =-.136553E-21 mag= 1.9999 -0.0000 0.0000

E_p = 0.45920E-40 lambda = 0.100E+01
ion lambda*MW_perp
1 0.00000E+00 -0.92644E-22 0.00000E+00
2 0.00000E+00 -0.52940E-22 -0.67763E-20

So I am not sure what causes E_p to be large in the calculations of I_CONSTRAINED_M =2. From Wiki https://www.vasp.at/wiki/index.php/I_CONSTRAINED_M, E_p should be made vanishingly small, is that right?

I further calculate the frequency (please see below) of the triplet O2 using I_CONSTRAINED_M =2. But the results show two unreasonably large frequencies (870 and 862 cm-1) in addition to the O-O stretch vibration frequency of 1771 cm-1. This means there must be something wrong with the calculation of I_CONSTRAINED_M =2, am I right?

1 f = 53.101034 THz 333.643639 2PiTHz 1771.259848 cm-1 219.608133 meV
2 f = 26.110033 THz 164.054175 2PiTHz 870.936949 cm-1 107.982371 meV
3 f = 25.843543 THz 162.379770 2PiTHz 862.047806 cm-1 106.880258 meV
4 f/i= 0.007422 THz 0.046637 2PiTHz 0.247588 cm-1 0.030697 meV
5 f/i= 0.022141 THz 0.139115 2PiTHz 0.738538 cm-1 0.091567 meV
6 f/i= 0.023756 THz 0.149264 2PiTHz 0.792416 cm-1 0.098247 meV

Thanks for your discusssion.

[fabien_tran1]

I have done calculations on O2 using your INCAR and POSCAR files. I could see that even without using the constrained local moments method, the triplet state is obtained (as indicated in OUTCAR, the total magnetic moment in the unit cell is 2, i.e., one unpaired electron at each O atom). Thus, no need to use the constrained local moments method. This should also explain why the penalty energy increases when LAMBDA increases: the difference between the target and calculated moments (M_I-M_I^0)^2 stays the same (actually very small) when LAMBDA is increased, making E_p=LAMBDA*(M_I-M_I^0)^2 larger.

[liang_zhang]

Thanks for your clarification. I got it now.

Sorry, I have one more question.
While the total magnetic moment is indeed ~2 (1.999), however, the sum of the local ion magnetization is only 1.631 (the magnetization of each ion is 0.816, not 1). This mismatch was found by both spin-polarised and constrained local moments methods. For constrained local moments calculations, I also tried to increase RWIGS, but that seemed to have no effect on the value of the local magnetization.

Can you please explain to me what's the reason of the mismatch and why RWIGS has no effect on the local magnetization value? Thanks again!

number of electron 12.0000004 magnetization 1.9998508 -0.0000000 0.0000000
augmentation part 1.1260084 magnetization 0.4917336 -0.0000000 0.0000000

magnetization (x)

# of ion s p d tot
------------------------------------------
1 0.009 0.806 0.000 0.816
2 0.009 0.806 0.000 0.816
--------------------------------------------------
tot 0.019 1.613 0.000 1.631

[fabien_tran1]

As indicated at RWIGS, the value of RWIGS is ignored if LORBIT>=10. With for instance LORBIT=5 and RWIGS=1.5, the atomic magnetic moment is 1.

[liang_zhang]

Thanks for the reply! I have a futher question for the following description.

In the Wiki https://www.vasp.at/wiki/index.php/I_CONSTRAINED_M, it said: "E_p is the contribution to the total energy arising from the penalty functional. Under M_int VASP lists the integrated magnetic moment at each atomic site. The column labeled MW_int shows the result of the integration of magnetization density which has been smoothed towards the boundary of the sphere. It is actually the smoothed integrated moment which enters in the penalty terms (the smoothing ensures that the total local potential remains continuous at the sphere boundary). One should look at the latter numbers to check whether enough of the magnetization density around each atomic site is contained within the integration sphere and increase RWIGS accordingly. What exactly constitutes "enough" in this context is hard to say. It is best to set RWIGS in such a manner that the integration spheres do not overlap and are otherwise as large as possible."

1. Does the "latter numbers" refer to the number under the label M_int?
2. How do I determine if the integration spheres do not overlap? It would be helpful if you can give me some idea.

[fabien_tran1]

1. It is indeed not so clear what "latter numbers" refers to. Anyway, the important is to choose values for RWIGS that are reasonable from the physical point of view.
2. It seems that no warning is issued if spheres overlap. Thus, one needs to compare manually the distance between two atoms and the sum of their radii RWIGS.

[liang_zhang]

The O-O distance is about 1.25 angstrom. If RWIGS = 1.5, the sum of RWIGS is 3.0, I am worried whether the value of 1.5 is too big? Does this indicate that a value of 1.5 would cause overlap?

In fact, the default value of RWIGS in POTCAR is 0.820, the sum is 1.64 which is still larger than the distance of O-O?

[fabien_tran1]

RWIGS has to be smaller than half of the O-O distance to avoid overlap of the spheres.

[liang_zhang]

Thanks, but the problem is the atomic magnetic moment of O would not be 1 if RWIGS is smaller than half of the O-O distance. As you pointed out in your previous reply, "With for instance LORBIT=5 and RWIGS=1.5, the atomic magnetic moment is 1.", it seems that the RWIGS has to be 1.5 such that the atomic magnetic moments is 1. Can you please clarify more on this? Thanks.

[fabien_tran1]

Yes, the magnetic moment inside the sphere depends on the size of the sphere.

[liang_zhang]

Sorry for repeatedly asking. Does this means RWIGS has to be set to 1.5 to get the correct atomic integrated magnetic moment even though the integration spheres will overlap?

[fabien_tran1]

Yes. Try yourself.