I have been trying to calculate the electronic structure of DyMnO3 with the PBE functionals. My calculations diverge instantly with swings more than 1000 Ha in energy between cycles. This does not seem to be an issue with rgkmax/gmaxvr, since my rgkmax is quite low and setting trimvg to 'true' does not fix the issue. I have not worked on lanthanide compounds before, so I don't know if I'm just missing something simple. The Eu example provided with ELK runs fine, and I did a test on Dysprosium metal which worked fine as well.
DyMnO3 is orthorhombic with the space group 62. I'm using the Pnma configuration with the crystal structure by Mori et al. ( http://dx.doi.org/10.1016/S0167-577X(99)00216-5 ). I'm using ELk 3.1.12 with the original species files Dy.in, Mn.in and O.in with the RMTs reduced to 2.6, 2.0 and 1.5 au, respectively. I tried a ferromagnetic calculation first (which failed in the same way) and went back to non-magnetic to more quickly try out different parameters. I have tried both Adaptive linear mixing and Broyden mixing, but I get the same divergence in both cases.
An example from INFO.OUT. Change in total energy seems to often be 100s or 1000s of Hartrees between cycles, as seen here. The main changes seem to be in the kinetic and Hartree energies.
I would appreciate any help with regards to this issue. Thank you!
Best regards,
Otto Mustonen
+--------------------+
| Loop number : 54 |
+--------------------+
Energies :
Fermi : 0.368344478857
sum of eigenvalues : -32992.8536033
electron kinetic : 59613.6150626
core electron kinetic : 57261.3796128
Coulomb : -111323.646732
Coulomb potential : -90844.7886251
nuclear-nuclear : -4489.11393464
electron-nuclear : -122824.276969
Hartree : 15989.7441720
Madelung : -65901.2524191
xc potential : -1761.68004092
exchange : -1327.97976603
correlation : -22.9480122077
electron entropic : -0.501961889793E-02
total energy : -53060.9644669
+--------------------+
| Loop number : 55 |
+--------------------+
Energies :
Fermi : 0.251133032256
sum of eigenvalues : -33088.1797285
electron kinetic : 58190.5569471
core electron kinetic : 57158.9117055
Coulomb : -112099.896666
Coulomb potential : -89466.7746535
nuclear-nuclear : -4489.11393464
electron-nuclear : -125754.790809
Hartree : 18144.0080779
Madelung : -67366.5093393
xc potential : -1811.96202203
exchange : -1365.11394497
correlation : -23.5930737165
electron entropic : -0.324199331623E-02
total energy : -55298.0499797
RMS change in Kohn-Sham potential (target) : 1.936629289 ( 0.1000000000$
Absolute change in total energy (target) : 1765.686264 ( 0.1000000000$
Here is my elk.in file:
tasks
0
10
xctype
20
swidth
0.001
beta0
0.05
betamax
0.5
reducebf
0.6
nwrite
1
ngridk
8 6 8
rgkmax
6.0
gmaxvr
18
lmaxapw
10
lmaxvr
8
lmaxmat
7
lradstp
2
nempty
8
maxscl
500
epsengy
1.e-4
epspot
1.e-6
wplot
960 160 1
-0.441 0.441
tshift
.false.
sppath
'./'
avec
11.03978000 0.000000000 0.000000000
0.000000000 13.94239900 0.000000000
0.000000000 0.000000000 9.977754000
this looks like the pecularity one runs into trying to calculate compounds including rare-earth elements within the framework of density functional theory (DFT). Partially filled 4f shells are a problem DFT cannot tackle directly because between iterations the 4f occupations changes too fast (the 4f electrons are very localized -> high density of states close to the Fermi level). If you look at the occupations numbers of Dy (or the interstitial) you will notice that they change severely from iteration to iteration. There are two solutions:
(a) DFT+U: force a magnetic solution with a close to integer 4f occupancy
(b) "open core" or "frozen core": fix the 4f core charge to a certain configuration by editing / generating the Dy.in file.
(c) DFT+DMFT ....
All the best,
Marc
p.s. a very small linear mixing is only a fake solution. The charge density changes might get small, the energy difference however won't converge nicely. The reason that the Eu example works is because of the high-spin ground state in combination with spin-orbit coupling
p.p.s. mixed valent behavior or paramagnetic ground states can only be achieved with (b) and (c)
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Thank you for the help! I wasn't able to get DFT+U to converge by myself, but I will try your elk.in and work from there. I was using FLL and not AMF like Kay, maybe that's one reason it didn't work out great.
Using Marc's suggestion to force the Dy 4f occupancy to 4f9, I was able to converge calculations. I moved the 4f states to core, and moved one electron from 4f to 5d1 to keep charge neutrality. This resulted in quite high core leakage and the Dy moment was only around 0.3 muB (spin-polarized core on).
I haven't had much time to work on this problem due to other headaches, but I'm very hopeful about getting DFT+U to work on this system starting from Kay's elk.in (and adding U and J for Mn 3d).
Best regards,
- Otto
Last edit: Otto Mustonen 2016-02-22
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
The reason for my DFT+U calculations failing appears to be having cmagz on. With full spin-orbit coupling the system slowly converges with very significant x and y components on the manganese. The dysprosium moments are FM and mostly in z direction. With cmagz on, on the other hand, there are huge swings in total energy between cycles and no convergence.
It's my understanding that one should only put fields in the z direction in elk.in. The experimental magnetic structure of DyMnO3 is cycloidal AFM. Thus, I would like to model system either in a simple AFM cell, an AFM supercell or FM => spin-spiral calculation. Which of these are applicable when full SOC is on, and how should I try to get an AFM state? Should I have fields in z direction only on the atoms or also in x and y?
One last question, INFO.OUT gives magnetic moments inside the atomic muffin tins. Now that spin-orbit coupling is on, how do I know what is the total moment, what is the spin-only moment and what is the orbital moment?
Best regards,
- Otto Mustonen
Last edit: Otto Mustonen 2016-02-25
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
first to get an AFM state you can specfiy the magnetic moments on Mn specfically like shown in the example section in the folder magnetism/Fe-AFM. Then in general since you are calculating a transition metal oxide you should also put a U and J on the Mn-d states !! This may also give you an insulating state for the FM magnetic arrangement. As a general strategy I would suggest you first to set up the Mn magnetic moments and do some iterations and then set also the f-states accordingly to how they have ordered in the first iterations. This strategy should allow you to get a good converged magnetic sate. I also would suggest to use a larger initial magnetic field as suggested in kays example. In my experience a value of 0.5-1 in atomic units with a reducebf of 0.75-0.85 is a good value for fast convergence with linear mixing. In this way you should be able to obtain an AFM state in the primitive unit-cell.
Now comes the hard part conerning the spiral. In princple you could use for this problem the generalized Bloch Theorem extension of the spin-spiral as implemented in elk (see example magnetism/Fe-spiral). With this approach you could set up a q-vector of your spiral in combination with the AFM magnetic ordering and calculate it without much effort in the primitive unit-cell (you should see in the manual for all tags related to that). However, in your case you have spin-orbit coupling switched on and thus you can not use the spin-spiral extension !!!. Thus either you manage to converge this calculation without the spin-orbit coupling or you have to set up a supercell with the spin spiral ordering.
I do not know the details of your system but I could imagine that in these systems like TbMnO3 a q=(1/3,0,0) give you a supercell of 3x1x1 the AFM ordering and thus the entire calculation becomes a huge thing. If you want to try this way you should use the big supercell in combination with ''iterative diagonalization'' tags, which can speed up the calculation in this case.
To its end I would suggest you to try again the LS off calculation to converge, which I managed for some f-electron systems by playing around with the initial magnetic field on the rare earth elements to ''trick'' the system into a stable state. I hope this information could help you.
Last point, could you show how you have putted the f-electrons in the core in the Dy.in file. I never tried this approach but would be keen to play around with it. 2nd You should alos check in your calculation what is the easy axis since this may also crucial for the spin spiral calculation !!
best regards
Michael
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
I played around with elk last weekend and found two possible solutions to your problem:
f electrons treated as frozen core
f without frozen core
I will post below how to setup the calculations and both converge nicely for me. In the first case, however, I obtained a quite big core leakage of 0.7E-1 on Dy and I am not sure how good tolaratable this still is.
1st case with f as frozen core
first the edited Dy.in with the f in the core and removed local orbital of third order of the f-states
'Dy' : spsymb
'dysprosium' : spname
-66.0000 : spzn
296219.3788 : spmass
0.246183E-06 2.8000 52.4581 700 : rminsp, rmt, rmaxsp, nrmt
21 : nstsp
1 0 1 2.00000 T : nsp, lsp, ksp, occsp, spcore
2 0 1 2.00000 T
2 1 1 2.00000 T
2 1 2 4.00000 T
3 0 1 2.00000 T
3 1 1 2.00000 T
3 1 2 4.00000 T
3 2 2 4.00000 T
3 2 3 6.00000 T
4 0 1 2.00000 T
4 1 1 2.00000 T
4 1 2 4.00000 T
4 2 2 4.00000 T
4 2 3 6.00000 T
4 3 3 6.00000 T
4 3 4 3.00000 T
5 0 1 2.00000 F
5 1 1 2.00000 F
5 1 2 4.00000 F
6 0 1 2.00000 F
7 0 1 1.00000 F
1 : apword
0.1500 0 F : apwe0, apwdm, apwve
0 : nlx
6 : nlorb
0 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
1 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
2 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
3 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
0 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-1.7987 0 T
1 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-0.8849 0 T
for the other in files I used the std. from elk 3.3.17
my in file with comments
tasks
0
mixtype
1
xctype
20
swidth
1E-4
! U on the Mn d-states as typical for TM oxides
! optional one can add a U on the empty Dy f-states which will otherwise end up slighlty above the
! Fermi energy
dft+u
1 1 : dftu, inpdftu
2 2 0.15 0.04 : is, l, U, J
spincore
.true.
spinpol
.true.
spinorb
.false.
ngridk
4 4 4
rgkmax
6.0
gmaxvr
18.0
lmaxapw
10
lmaxvr
8
lmaxmat
7
nempty
8
reducebf
0.8
nrmtscf
2.000000000
avec
11.03978000 0.000000000 0.000000000
0.000000000 13.94239900 0.000000000
0.000000000 0.000000000 9.977754000
!please note I set up and AFM state but FM should work as well
this calculation converges within 45 steps !! I found a moment of 3.5 muB on Mn confirming the Mn3+ state as expected.
2. option with f-states regular species files for all elements and
tasks
0
mixtype
1
xctype
20
swidth
1E-4
! Here I add the U on the f-states
dft+u
1 1 : dftu, inpdftu
1 3 0.5000000 0.00000000
2 2 0.15 0.04 : is, l, U, J
spinpol
.true.
spinorb
.false.
ngridk
4 4 4
rgkmax
6.0
gmaxvr
18.0
lmaxapw
10
lmaxvr
8
lmaxmat
7
nempty
8
reducebf
0.8
nrmtscf
2.000000000
dosmsum
.true.
wplot
960 160 1
-0.8 0.8
avec
11.03978000 0.000000000 0.000000000
0.000000000 13.94239900 0.000000000
0.000000000 0.000000000 9.977754000
! similiar to the first case big initial field on Mn and additional big field on the Dy atoms to polarize the
! f-states. The order of the field is obtained as explained in my previous posts above
!
This calculation converges within 60 steps and reveals again 3.48 muB on Mn confiriming the Mn3+ state and in addition a moment of 5 muB on Dy confiriming Dy 3+ with 9 electrons in the 4f shell
hope this helped if you have further questions please feel free to ask. Finally if somebody has a suggestion for a critical value for the core leakage I would also appreaciate that.
best regard
Michael
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
The Dy 3+ ion gives away the 6s2 and one of the 4f electrons consequenlty we have 9 electrons in the 4f-shell. Considering Paulis and Hunds Rule we have 7 electrons with up spin and 2 with down spin summing up to <S>=51/2 which results in a magnetic moment of m=g<S>=5 muB. (g=2 for spins !!).
Please note that the maximal magnetic moment of an f shell is 7 mub (7 electrons with all parallel spin) like for Gd3+ or Eu2+
I thought that you have obtained the effective magnetic moment for Dy3+.
The main term for Dy3+ is 6H7.5 with S=2.5, L = 5 ===> J=7.5. With gj = 1.33 and using m=gj(sqr(J(J+1))) we get ~10.62 muB and that is why I was wondering about 5muB. But now I see what kind of "moment" you were talking about. :)
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
This calculation is indeed a tricky case and convergence is not straightforward to obtain. From the structure, you send me I failed to converge here anything with the frozen core method. Okay then I had a look at the chemistry, and we have from the periodic table the following configuration for Ce 6s2 4f1 5d1. Reflecting this configuration with the actual settings in elk, I realized that we are missing 5d1 to create an f1 core configuration. Hence I introduced the 5d1 and lowered the occupations of the 4f state by one. The resulting element infiile reads :
'Ce' : spsymb
'cerium' : spname
-58.0000 : spzn
255415.8429 : spmass
0.262613E-06 2.8000 70.2137 700 : rminsp, rmt, rmaxsp, nrmt
20 : nstsp << -- this number was increased by one since we have a orbital more
1 0 1 2.00000 T : nsp, lsp, ksp, occsp, spcore
2 0 1 2.00000 T
2 1 1 2.00000 T
2 1 2 4.00000 T
3 0 1 2.00000 T
3 1 1 2.00000 T
3 1 2 4.00000 T
3 2 2 4.00000 T
3 2 3 6.00000 T
4 0 1 2.00000 T
4 1 1 2.00000 T
4 1 2 4.00000 T
4 2 2 4.00000 T
4 2 3 6.00000 T
4 3 3 1.00000 T << -- Here the 4f1 in core
5 0 1 2.00000 F
5 1 1 2.00000 F
5 1 2 4.00000 F
5 2 2 1.00000 F <<-- Here the added 5d1 orbital
6 0 1 2.00000 F
1 : apword
0.1500 0 F : apwe0, apwdm, apwve
0 : nlx
5 : nlorb
0 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
1 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
2 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
0 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-1.3797 0 T
1 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-0.7361 0 T
Please note the comments marked by <<. With this file, your configuration runs smoothly through. However, here some points to consider:
1.) I noticed from "EVALCORE.OUT", that the 4f1 state is quite high in energy, actually crossing already the valence states listed in EIGVAL.OUT. Consequently, I wonder if one really can separate things here in this way. The compound you consider is, moreover, not an insulator and hence the f electrons might contribute to the bonding. I am not an expert, but you might want to think about this point.
2.) I also checked the regular computation of your compound, which also went through fine. Hence I am wondering what the advantage of having here the 4f1 in the core is if there is any.
best regards
Michael
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
First of all, thank you very much for you help. I really appreciate your willingness to assist.
1) "I had a look at the chemistry, and we have from the periodic table the following configuration for Ce 6s2 4f1 5d1.
Reflecting this configuration with the actual settings in elk, I realized that we are missing 5d1 to create an f1 core configuration." Could you clarify - you mean that the
configuration of Ce is wrong in ELK species files?
2) I want to move Ce-4f electrons away from EF. I tried to do LDA+U, but even if I use U=7-8, some parts of Ce4f modify states around EF. For example, here authors just remove Ce-4f and agreement between theory and experiment is pretty good. https://arxiv.org/pdf/1411.4013.pdf
Thank you
Jammi
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
regarding your question 1) The term ''wrong'' sound more dramatic as it is. If you have a look on the isolated atomic filling of Ce (https://www.ptable.com/#Orbital) you find a single 5d1 electron and a single 4f1 one, whereas the confirguration elk uses to start with is 5d0 and 4f2. That different in occupation is what I was refereing to. In the end it shouldn't matter for a regular calculation but if you want to put something into the core you should be taking care.
regarding two but the big ''peaks'' of the 4f are away in these calculation ?
best regards
Michael
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Dear ELK users and developers
I have been trying to calculate the electronic structure of DyMnO3 with the PBE functionals. My calculations diverge instantly with swings more than 1000 Ha in energy between cycles. This does not seem to be an issue with rgkmax/gmaxvr, since my rgkmax is quite low and setting trimvg to 'true' does not fix the issue. I have not worked on lanthanide compounds before, so I don't know if I'm just missing something simple. The Eu example provided with ELK runs fine, and I did a test on Dysprosium metal which worked fine as well.
DyMnO3 is orthorhombic with the space group 62. I'm using the Pnma configuration with the crystal structure by Mori et al. ( http://dx.doi.org/10.1016/S0167-577X(99)00216-5 ). I'm using ELk 3.1.12 with the original species files Dy.in, Mn.in and O.in with the RMTs reduced to 2.6, 2.0 and 1.5 au, respectively. I tried a ferromagnetic calculation first (which failed in the same way) and went back to non-magnetic to more quickly try out different parameters. I have tried both Adaptive linear mixing and Broyden mixing, but I get the same divergence in both cases.
An example from INFO.OUT. Change in total energy seems to often be 100s or 1000s of Hartrees between cycles, as seen here. The main changes seem to be in the kinetic and Hartree energies.
I would appreciate any help with regards to this issue. Thank you!
Best regards,
Otto Mustonen
+--------------------+
| Loop number : 54 |
+--------------------+
Energies :
Fermi : 0.368344478857
sum of eigenvalues : -32992.8536033
electron kinetic : 59613.6150626
core electron kinetic : 57261.3796128
Coulomb : -111323.646732
Coulomb potential : -90844.7886251
nuclear-nuclear : -4489.11393464
electron-nuclear : -122824.276969
Hartree : 15989.7441720
Madelung : -65901.2524191
xc potential : -1761.68004092
exchange : -1327.97976603
correlation : -22.9480122077
electron entropic : -0.501961889793E-02
total energy : -53060.9644669
+--------------------+
| Loop number : 55 |
+--------------------+
Energies :
Fermi : 0.251133032256
sum of eigenvalues : -33088.1797285
electron kinetic : 58190.5569471
core electron kinetic : 57158.9117055
Coulomb : -112099.896666
Coulomb potential : -89466.7746535
nuclear-nuclear : -4489.11393464
electron-nuclear : -125754.790809
Hartree : 18144.0080779
Madelung : -67366.5093393
xc potential : -1811.96202203
exchange : -1365.11394497
correlation : -23.5930737165
electron entropic : -0.324199331623E-02
total energy : -55298.0499797
RMS change in Kohn-Sham potential (target) : 1.936629289 ( 0.1000000000$
Absolute change in total energy (target) : 1765.686264 ( 0.1000000000$
Here is my elk.in file:
tasks
0
10
xctype
20
swidth
0.001
beta0
0.05
betamax
0.5
reducebf
0.6
nwrite
1
ngridk
8 6 8
rgkmax
6.0
gmaxvr
18
lmaxapw
10
lmaxvr
8
lmaxmat
7
lradstp
2
nempty
8
maxscl
500
epsengy
1.e-4
epspot
1.e-6
wplot
960 160 1
-0.441 0.441
tshift
.false.
sppath
'./'
avec
11.03978000 0.000000000 0.000000000
0.000000000 13.94239900 0.000000000
0.000000000 0.000000000 9.977754000
atoms
3 : nspecies
'Dy.in' : spfname
4 : natoms; atposl, bfcmt below
0.08248000 0.25000000 0.98260000 0.00000000 0.00000000 0.00000000
0.91752000 0.75000000 0.01740000 0.00000000 0.00000000 0.00000000
0.41752000 0.75000000 0.48260000 0.00000000 0.00000000 0.00000000
0.58248000 0.25000000 0.51740000 0.00000000 0.00000000 0.00000000
'Mn.in' : spfname
4 : natoms; atposl, bfcmt below
0.00000000 0.00000000 0.50000000 0.00000000 0.00000000 0.00000000
0.50000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
0.50000000 0.50000000 0.00000000 0.00000000 0.00000000 0.00000000
0.00000000 0.50000000 0.50000000 0.00000000 0.00000000 0.00000000
'O.in' : spfname
12 : natoms; atposl, bfcmt below
0.46300000 0.25000000 0.11400000 0.00000000 0.00000000 0.00000000
0.53700000 0.75000000 0.88600000 0.00000000 0.00000000 0.00000000
0.03700000 0.75000000 0.61400000 0.00000000 0.00000000 0.00000000
0.96300000 0.25000000 0.38600000 0.00000000 0.00000000 0.00000000
0.17300000 0.55340000 0.20100000 0.00000000 0.00000000 0.00000000
0.82700000 0.44660000 0.79900000 0.00000000 0.00000000 0.00000000
0.32700000 0.44660000 0.70100000 0.00000000 0.00000000 0.00000000
0.67300000 0.55340000 0.29900000 0.00000000 0.00000000 0.00000000
0.67300000 0.94660000 0.29900000 0.00000000 0.00000000 0.00000000
0.32700000 0.05340000 0.70100000 0.00000000 0.00000000 0.00000000
0.82700000 0.05340000 0.79900000 0.00000000 0.00000000 0.00000000
0.17300000 0.94660000 0.20100000 0.00000000 0.00000000 0.00000000
Dear Otto,
this looks like the pecularity one runs into trying to calculate compounds including rare-earth elements within the framework of density functional theory (DFT). Partially filled 4f shells are a problem DFT cannot tackle directly because between iterations the 4f occupations changes too fast (the 4f electrons are very localized -> high density of states close to the Fermi level). If you look at the occupations numbers of Dy (or the interstitial) you will notice that they change severely from iteration to iteration. There are two solutions:
(a) DFT+U: force a magnetic solution with a close to integer 4f occupancy
(b) "open core" or "frozen core": fix the 4f core charge to a certain configuration by editing / generating the Dy.in file.
(c) DFT+DMFT ....
All the best,
Marc
p.s. a very small linear mixing is only a fake solution. The charge density changes might get small, the energy difference however won't converge nicely. The reason that the Eu example works is because of the high-spin ground state in combination with spin-orbit coupling
p.p.s. mixed valent behavior or paramagnetic ground states can only be achieved with (b) and (c)
Dear supermarche
Thank you for explaining the issue! I will try at least a) and c), this was very helpful.
Best regards,
- Otto Mustonen
Dear Otto and Marc,
Marc is quite correct about the difficultly of converging some f-electron systems.
I've managed to converge your material using DFT+U. Here is my elk.in:
... note that this is not properly converged, and the U and J parameters were arbitrarily chosen. Note also that spin-orbit coupling is switched on.
Also, you used 'reducebf' but without a magnetic field to break the symmetry or setting 'spinpol=.true.'
This system should yield high orbital moments on the Dy atoms.
Regards,
Kay.
Dear Marc and Kay
Thank you for the help! I wasn't able to get DFT+U to converge by myself, but I will try your elk.in and work from there. I was using FLL and not AMF like Kay, maybe that's one reason it didn't work out great.
Using Marc's suggestion to force the Dy 4f occupancy to 4f9, I was able to converge calculations. I moved the 4f states to core, and moved one electron from 4f to 5d1 to keep charge neutrality. This resulted in quite high core leakage and the Dy moment was only around 0.3 muB (spin-polarized core on).
I haven't had much time to work on this problem due to other headaches, but I'm very hopeful about getting DFT+U to work on this system starting from Kay's elk.in (and adding U and J for Mn 3d).
Best regards,
- Otto
Last edit: Otto Mustonen 2016-02-22
Dear Marc and Kay
The reason for my DFT+U calculations failing appears to be having cmagz on. With full spin-orbit coupling the system slowly converges with very significant x and y components on the manganese. The dysprosium moments are FM and mostly in z direction. With cmagz on, on the other hand, there are huge swings in total energy between cycles and no convergence.
It's my understanding that one should only put fields in the z direction in elk.in. The experimental magnetic structure of DyMnO3 is cycloidal AFM. Thus, I would like to model system either in a simple AFM cell, an AFM supercell or FM => spin-spiral calculation. Which of these are applicable when full SOC is on, and how should I try to get an AFM state? Should I have fields in z direction only on the atoms or also in x and y?
One last question, INFO.OUT gives magnetic moments inside the atomic muffin tins. Now that spin-orbit coupling is on, how do I know what is the total moment, what is the spin-only moment and what is the orbital moment?
Best regards,
- Otto Mustonen
Last edit: Otto Mustonen 2016-02-25
Dear Otto,
first to get an AFM state you can specfiy the magnetic moments on Mn specfically like shown in the example section in the folder magnetism/Fe-AFM. Then in general since you are calculating a transition metal oxide you should also put a U and J on the Mn-d states !! This may also give you an insulating state for the FM magnetic arrangement. As a general strategy I would suggest you first to set up the Mn magnetic moments and do some iterations and then set also the f-states accordingly to how they have ordered in the first iterations. This strategy should allow you to get a good converged magnetic sate. I also would suggest to use a larger initial magnetic field as suggested in kays example. In my experience a value of 0.5-1 in atomic units with a reducebf of 0.75-0.85 is a good value for fast convergence with linear mixing. In this way you should be able to obtain an AFM state in the primitive unit-cell.
Now comes the hard part conerning the spiral. In princple you could use for this problem the generalized Bloch Theorem extension of the spin-spiral as implemented in elk (see example magnetism/Fe-spiral). With this approach you could set up a q-vector of your spiral in combination with the AFM magnetic ordering and calculate it without much effort in the primitive unit-cell (you should see in the manual for all tags related to that). However, in your case you have spin-orbit coupling switched on and thus you can not use the spin-spiral extension !!!. Thus either you manage to converge this calculation without the spin-orbit coupling or you have to set up a supercell with the spin spiral ordering.
I do not know the details of your system but I could imagine that in these systems like TbMnO3 a q=(1/3,0,0) give you a supercell of 3x1x1 the AFM ordering and thus the entire calculation becomes a huge thing. If you want to try this way you should use the big supercell in combination with ''iterative diagonalization'' tags, which can speed up the calculation in this case.
To its end I would suggest you to try again the LS off calculation to converge, which I managed for some f-electron systems by playing around with the initial magnetic field on the rare earth elements to ''trick'' the system into a stable state. I hope this information could help you.
Last point, could you show how you have putted the f-electrons in the core in the Dy.in file. I never tried this approach but would be keen to play around with it. 2nd You should alos check in your calculation what is the easy axis since this may also crucial for the spin spiral calculation !!
best regards
Michael
one more point more directed to kay,
is there a reason why you use
thsift
.false.
in your calculation? Is there a deeper reason ?
Dear Otto,
I played around with elk last weekend and found two possible solutions to your problem:
f electrons treated as frozen core
f without frozen core
I will post below how to setup the calculations and both converge nicely for me. In the first case, however, I obtained a quite big core leakage of 0.7E-1 on Dy and I am not sure how good tolaratable this still is.
1st case with f as frozen core
first the edited Dy.in with the f in the core and removed local orbital of third order of the f-states
'Dy' : spsymb
'dysprosium' : spname
-66.0000 : spzn
296219.3788 : spmass
0.246183E-06 2.8000 52.4581 700 : rminsp, rmt, rmaxsp, nrmt
21 : nstsp
1 0 1 2.00000 T : nsp, lsp, ksp, occsp, spcore
2 0 1 2.00000 T
2 1 1 2.00000 T
2 1 2 4.00000 T
3 0 1 2.00000 T
3 1 1 2.00000 T
3 1 2 4.00000 T
3 2 2 4.00000 T
3 2 3 6.00000 T
4 0 1 2.00000 T
4 1 1 2.00000 T
4 1 2 4.00000 T
4 2 2 4.00000 T
4 2 3 6.00000 T
4 3 3 6.00000 T
4 3 4 3.00000 T
5 0 1 2.00000 F
5 1 1 2.00000 F
5 1 2 4.00000 F
6 0 1 2.00000 F
7 0 1 1.00000 F
1 : apword
0.1500 0 F : apwe0, apwdm, apwve
0 : nlx
6 : nlorb
0 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
1 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
2 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
3 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
0 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-1.7987 0 T
1 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-0.8849 0 T
for the other in files I used the std. from elk 3.3.17
my in file with comments
tasks
0
mixtype
1
xctype
20
swidth
1E-4
! U on the Mn d-states as typical for TM oxides
! optional one can add a U on the empty Dy f-states which will otherwise end up slighlty above the
! Fermi energy
dft+u
1 1 : dftu, inpdftu
2 2 0.15 0.04 : is, l, U, J
spincore
.true.
spinpol
.true.
spinorb
.false.
ngridk
4 4 4
rgkmax
6.0
gmaxvr
18.0
lmaxapw
10
lmaxvr
8
lmaxmat
7
nempty
8
reducebf
0.8
nrmtscf
2.000000000
avec
11.03978000 0.000000000 0.000000000
0.000000000 13.94239900 0.000000000
0.000000000 0.000000000 9.977754000
!please note I set up and AFM state but FM should work as well
atoms
3 : nspecies
'Dy.in' : spfname
4 : natoms; atposl, bfcmt below
0.08248000 0.25000000 0.98260000
0.91752000 0.75000000 0.01740000
0.41752000 0.75000000 0.48260000
0.58248000 0.25000000 0.51740000
Mn.in'
4
0.00000000 0.00000000 0.50000000 0 0 15
0.50000000 0.00000000 0.00000000 0 0 -15
0.50000000 0.50000000 0.00000000 0 0 -15
0.00000000 0.50000000 0.50000000 0 0 15
O.in'
12
0.46300000 0.25000000 0.11400000
0.53700000 0.75000000 0.88600000
0.03700000 0.75000000 0.61400000
0.96300000 0.25000000 0.38600000
0.17300000 0.55340000 0.20100000
0.82700000 0.44660000 0.79900000
0.32700000 0.44660000 0.70100000
0.67300000 0.55340000 0.29900000
0.67300000 0.94660000 0.29900000
0.32700000 0.05340000 0.70100000
0.82700000 0.05340000 0.79900000
0.17300000 0.94660000 0.20100000
this calculation converges within 45 steps !! I found a moment of 3.5 muB on Mn confirming the Mn3+ state as expected.
2. option with f-states regular species files for all elements and
tasks
0
mixtype
1
xctype
20
swidth
1E-4
! Here I add the U on the f-states
dft+u
1 1 : dftu, inpdftu
1 3 0.5000000 0.00000000
2 2 0.15 0.04 : is, l, U, J
spinpol
.true.
spinorb
.false.
ngridk
4 4 4
rgkmax
6.0
gmaxvr
18.0
lmaxapw
10
lmaxvr
8
lmaxmat
7
nempty
8
reducebf
0.8
nrmtscf
2.000000000
dosmsum
.true.
wplot
960 160 1
-0.8 0.8
avec
11.03978000 0.000000000 0.000000000
0.000000000 13.94239900 0.000000000
0.000000000 0.000000000 9.977754000
! similiar to the first case big initial field on Mn and additional big field on the Dy atoms to polarize the
! f-states. The order of the field is obtained as explained in my previous posts above
!
atoms
3 : nspecies
'Dy.in' : spfname
4 : natoms; atposl, bfcmt below
0.08248000 0.25000000 0.98260000 0 0 -15
0.91752000 0.75000000 0.01740000 0 0 -15
0.41752000 0.75000000 0.48260000 0 0 15
0.58248000 0.25000000 0.51740000 0 0 15
'Mn.in'
4
0.00000000 0.00000000 0.50000000 0 0 15
0.50000000 0.00000000 0.00000000 0 0 -15
0.50000000 0.50000000 0.00000000 0 0 -15
0.00000000 0.50000000 0.50000000 0 0 15
'O.in'
12
0.46300000 0.25000000 0.11400000
0.53700000 0.75000000 0.88600000
0.03700000 0.75000000 0.61400000
0.96300000 0.25000000 0.38600000
0.17300000 0.55340000 0.20100000
0.82700000 0.44660000 0.79900000
0.32700000 0.44660000 0.70100000
0.67300000 0.55340000 0.29900000
0.67300000 0.94660000 0.29900000
0.32700000 0.05340000 0.70100000
0.82700000 0.05340000 0.79900000
0.17300000 0.94660000 0.20100000
This calculation converges within 60 steps and reveals again 3.48 muB on Mn confiriming the Mn3+ state and in addition a moment of 5 muB on Dy confiriming Dy 3+ with 9 electrons in the 4f shell
hope this helped if you have further questions please feel free to ask. Finally if somebody has a suggestion for a critical value for the core leakage I would also appreaciate that.
best regard
Michael
Only 5muB for Dy3+ ion seems to be not enough. It should be at least 10.64 muB per atom(ion). Or. I have missed something?
The Dy 3+ ion gives away the 6s2 and one of the 4f electrons consequenlty we have 9 electrons in the 4f-shell. Considering Paulis and Hunds Rule we have 7 electrons with up spin and 2 with down spin summing up to <S>=51/2 which results in a magnetic moment of m=g<S>=5 muB. (g=2 for spins !!).
Please note that the maximal magnetic moment of an f shell is 7 mub (7 electrons with all parallel spin) like for Gd3+ or Eu2+
For more details I always recommend --> http://www.ptable.com/#Orbital and then Dy
Last edit: mfechner 2016-03-03
I thought that you have obtained the effective magnetic moment for Dy3+.
The main term for Dy3+ is 6H7.5 with S=2.5, L = 5 ===> J=7.5. With gj = 1.33 and using m=gj(sqr(J(J+1))) we get ~10.62 muB and that is why I was wondering about 5muB. But now I see what kind of "moment" you were talking about. :)
Hi All,
I know this is an "old method", but how to setup the calculations for Ce 4f electrons as a core frozen? (4f1 as a core state)
Here is my Ce.in - which lo should I remove? When I remove 3 2 the DTOTENERGY and core leaking of Ce fluctuate like a crazy
'Ce' : spsymb
'cerium' : spname
-58.0000 : spzn
255415.8429 : spmass
0.262613E-06 2.8000 70.2137 700 : rminsp, rmt, rmaxsp, nrmt
19 : nstsp
1 0 1 2.00000 T : nsp, lsp, ksp, occsp, spcore
2 0 1 2.00000 T
2 1 1 2.00000 T
2 1 2 4.00000 T
3 0 1 2.00000 T
3 1 1 2.00000 T
3 1 2 4.00000 T
3 2 2 4.00000 T
3 2 3 6.00000 T
4 0 1 2.00000 T
4 1 1 2.00000 T
4 1 2 4.00000 T
4 2 2 4.00000 T
4 2 3 6.00000 T
4 3 3 2.00000 T <--- from F to T...
5 0 1 2.00000 F
5 1 1 2.00000 F
5 1 2 4.00000 F
6 0 1 2.00000 F
1 : apword
0.1500 0 F : apwe0, apwdm, apwve
0 : nlx
6 : nlorb
0 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
1 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
2 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
3 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
0 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-1.3797 0 T
1 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-0.7361 0 T
Thanks
Jammi
Dear Jammi,
This calculation is indeed a tricky case and convergence is not straightforward to obtain. From the structure, you send me I failed to converge here anything with the frozen core method. Okay then I had a look at the chemistry, and we have from the periodic table the following configuration for Ce 6s2 4f1 5d1. Reflecting this configuration with the actual settings in elk, I realized that we are missing 5d1 to create an f1 core configuration. Hence I introduced the 5d1 and lowered the occupations of the 4f state by one. The resulting element infiile reads :
'Ce' : spsymb
'cerium' : spname
-58.0000 : spzn
255415.8429 : spmass
0.262613E-06 2.8000 70.2137 700 : rminsp, rmt, rmaxsp, nrmt
20 : nstsp << -- this number was increased by one since we have a orbital more
1 0 1 2.00000 T : nsp, lsp, ksp, occsp, spcore
2 0 1 2.00000 T
2 1 1 2.00000 T
2 1 2 4.00000 T
3 0 1 2.00000 T
3 1 1 2.00000 T
3 1 2 4.00000 T
3 2 2 4.00000 T
3 2 3 6.00000 T
4 0 1 2.00000 T
4 1 1 2.00000 T
4 1 2 4.00000 T
4 2 2 4.00000 T
4 2 3 6.00000 T
4 3 3 1.00000 T << -- Here the 4f1 in core
5 0 1 2.00000 F
5 1 1 2.00000 F
5 1 2 4.00000 F
5 2 2 1.00000 F <<-- Here the added 5d1 orbital
6 0 1 2.00000 F
1 : apword
0.1500 0 F : apwe0, apwdm, apwve
0 : nlx
5 : nlorb
0 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
1 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
2 2 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
0 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-1.3797 0 T
1 3 : lorbl, lorbord
0.1500 0 F : lorbe0, lorbdm, lorbve
0.1500 1 F
-0.7361 0 T
Please note the comments marked by <<. With this file, your configuration runs smoothly through. However, here some points to consider:
1.) I noticed from "EVALCORE.OUT", that the 4f1 state is quite high in energy, actually crossing already the valence states listed in EIGVAL.OUT. Consequently, I wonder if one really can separate things here in this way. The compound you consider is, moreover, not an insulator and hence the f electrons might contribute to the bonding. I am not an expert, but you might want to think about this point.
2.) I also checked the regular computation of your compound, which also went through fine. Hence I am wondering what the advantage of having here the 4f1 in the core is if there is any.
best regards
Michael
Dear Michael,
Thank you very much for contributing generously.
I modify the Sm.in file according to your method.
But the following warning was received. It seems that ELK directly ignores the electrons in my modified orbit.
I also tried to delete the following contents in the Sm.in file, but the warnings received were same.
Here is my elk.in file.
I look forward to your reply.
Last edit: Long Yu 2018-10-26
Dear Michael,
First of all, thank you very much for you help. I really appreciate your willingness to assist.
1) "I had a look at the chemistry, and we have from the periodic table the following configuration for Ce 6s2 4f1 5d1.
Reflecting this configuration with the actual settings in elk, I realized that we are missing 5d1 to create an f1 core configuration." Could you clarify - you mean that the
configuration of Ce is wrong in ELK species files?
2) I want to move Ce-4f electrons away from EF. I tried to do LDA+U, but even if I use U=7-8, some parts of Ce4f modify states around EF. For example, here authors just remove Ce-4f and agreement between theory and experiment is pretty good.
https://arxiv.org/pdf/1411.4013.pdf
Thank you
Jammi
Dear Jammi,
regarding your question 1) The term ''wrong'' sound more dramatic as it is. If you have a look on the isolated atomic filling of Ce (https://www.ptable.com/#Orbital) you find a single 5d1 electron and a single 4f1 one, whereas the confirguration elk uses to start with is 5d0 and 4f2. That different in occupation is what I was refereing to. In the end it shouldn't matter for a regular calculation but if you want to put something into the core you should be taking care.
regarding two but the big ''peaks'' of the 4f are away in these calculation ?
best regards
Michael