[Wien] Constrained DFT+U on Fe-doped BaTiO3: Large discrepancies between setup density matrix, IPDOS, and output in .scfdm
Peter Blaha
peter.blaha at tuwien.ac.at
Fri Feb 20 15:53:50 CET 2026
i) Yes, WIEN2k will always use its symmetry operation. If your
site-symmetry is such that d-xz+yz are degenerate, you should not break
symmetry in your dmat/vorb files.
ii) If you want to go to lower symmetry, you should displace some
atoms/or make them inequivalent during the init_lapw phase. This should
produces symmetries, which allow at least in principle a breaking of
this degeneracy. AFTER init_lapw, you can move the atoms back to the
original position.
iii) Your constraint dmat is NOT dyz but refers to an occupation of
an complex orbital of l=2, m=-1.
The "real" spherical harmonics are linear combinations of the complex
ones. Checkout some text books.
iv) You cannot force a value of "1" for Ti/Fe-d. Note, these are charges
just inside the sphere and the 3d orbital of Ti is too localized to be
fully inside the sphere. Checkout case.outputst to see how much of the
atomic d-orbital is within the sphere.
iv) Depending on the value of U, you may NOT be able to force a
particular occupation. It could be that you need a huge U value to do
this, but in principle it looks "ok".
up-dmat has 1.32 dxz+yz electrons (see above), "nothing" else.
dn-dmat: has the smallest d-xz+yz (0.42 per single orbital), because
this has "0 + 1" in your dmat, while the others have all full "1" and
are all much higher.
v) I don't quite understand your outputs.
Eg:
case.scfdmdn (Target was dxy, dz2, dxz, dx2-y2):
irtest 1 8 2.07000000000000
Density matrix UPUP block, real part. L= 2
0.94475 0.00000 0.00000 0.00000 0.01576
0.00000 0.48222 0.00000 0.00000 0.00000
0.00000 0.00000 0.90325 0.00000 0.00000
0.00000 0.00000 0.00000 0.48222 0.00000
0.01576 0.00000 0.00000 0.00000 0.94475
DN d: 3.70166, dz2: 0.72795, dxy: 0.82023, dx2y2: 1.30642, dxz+dyz:
0.84706 (from where are these numbers ??)
Not even the trace is correct. The :QTLxxx from lapw2 should be
"identical" to your dmat.
PS: In general I'm skeptical about the method of "multiplet" tests. This
would work only with very large U values, because you have to overcome
the "crystal field" and this is meaningless.
In my experience it is best to start with a fully converged PBE
calculation. While these occupations are often too "spherical", they
usually point into the right direction (i.e. a particular orbital is a
bit less occupied than another). This effect gets then "stronger" with U.
On 2/20/26 12:25 PM, Li, Zhiyuan wrote:
> (Apologies to the list for the previous email sent without a subject due
> to a misclick. Please refer to this message for the complete details.)
>
> Dear WIEN2k community and Prof. Blaha,
>
>
> I am currently performing constrained DFT+U calculations on a 2x2x2
> BaTiO3 supercell with an Fe dopant and an oxygen vacancy using WIEN2k
> 24.1. My ultimate goal is to systematically investigate all possible d-
> orbital electronic multiplets (following the methodology in Dorado et
> al., PRB 79, 235125 (2009) and previous discussions on this mailing list
> by Bin Shao in 2015 http://zeus.theochem.tuwien.ac.at/pipermail/
> wien/2015-July/023196.html <http://zeus.theochem.tuwien.ac.at/pipermail/
> wien/2015-July/023196.html> and Jayangani in 2019 http://
> zeus.theochem.tuwien.ac.at/pipermail/wien/2019-July/029531.html <http://
> zeus.theochem.tuwien.ac.at/pipermail/wien/2019-July/029531.html>).
>
> *The Core Issue:*
>
> I am observing a significant deviation (> 0.3 e-) between my explicitly
> constructed initial density matrix (|.dmatup/dn|), the IPDOS up to
> EF(Fe), and the final output density matrix obtained via |lapwdm|, even
> though the calculation converges successfully and i haven't released the
> constraint.
>
> To test the constraint, I intentionally chose a non-standard
> configuration to ensure I can apply constraints beyond estimated egt2g
> spliting:
>
> Target UP: Only dyz occupied (1.0)
>
> Target DN: dxy, dz2, dxz, dx2-y2 occupied (1.0 each)
>
> Here is a figure of the summary table of the discrepancies I observed:
> (with all smaller than 0.1 components ignored)
>
> PS: I can only set d,dxy,dz2,dx2-y2,dxz+yz for this system in the qtl.
>
> (I also considered the local rotation matrix, but since they share the
> same rotation and the tetragonal Ti site is not severely distorted, I
> believe it cannot fully explain this large deviation.)
> *
> My Questions:*
> large deviation while still constrained: I observe this huge discrepancy
> even before releasing the constraint. For example, my single dyz target
> was split into two ~0.85 occupations in .scfdmup. Does this imply that
> the system's symmetry is still actively forcing a degeneracy between xz
> and yz?
>
> Normal behavior vs. Setup error: Is it normal to see such a large
> deviation between the setup .dmat and the output .scfdm while the
> orbital potential is still being actively applied, or does this indicate
> an error in my initialization/setup?
>
> Physical/Algorithmic interpretation: If my setup is perfectly fine, how
> should I understand what is happening here? Is the imposed U potential
> simply not strong enough to overcome the hybridization and SCF mixing,
> thereby "washing out" my imposed constraint during the cycles?
>
> Alternative methods for a "hard lock": Are there other techniques to
> more rigidly lock the desired electronic multiplet? For instance, should
> I artificially increase the $U$ value during the initial SCF cycles?
>
> *ALL Technical Details:*
>
> 1. setup:
> init_lapw -b -sp
> prepare .indm & .inorb & .dmatup/dn
> x orb -up -p
> x orb -dn -p
> runsp_lapw -p -orbc -i 500 -ec 0.0001 -cc 0.001
>
> 2. .indm:
>
> -12. Emin cutoff energy
> 1 number of atoms for which density matrix is
> calculated
> 8 1 2 index of 1st atom, number of L's, L1
> 0 0 r-index, (l,s)index
>
> 3. .inorb
> 1 1 0 nmod, natorb, ipr
> PRATT 1.0 BROYD/PRATT, mixing
> 8 1 2 iatom nlorb, lorb
> 1 nsic 0..AMF, 1..SIC, 2..HFM
> 0.37 0.07 U J (Ry)
>
> 4. initial .dmatup(target dyz)
>
> 8 atom density matrix
> 2 0.000000 0.000000 0.000000 L, Lx,Ly,Lz in global orthogonal system
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 1.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00
>
> 5. initial .dmatdn(target dxy, dz2, dxz, dx2-y2)
> 8 atom density matrix
> 2 0.000000 0.000000 0.000000 L, Lx,Ly,Lz in global orthogonal system
>
> 1.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 1.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 1.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 0.0000000000000E+00 0.0000000000000E+00 0.0000000000000E+00
> 0.0000000000000E+00
>
> 1.0000000000000E+00 0.0000000000000E+00
>
> 6. Post-processing for IPDOS up to EF:
> x lapw2 -qtl -up -p
> x lapw2 -qtl -dn -p
> configure_int_lapw -b total 8 d,dz2,dxy,dx2y2,dxz+dyz
> x tetra -up -p
> x tetra -dn -p
>
> UP (Target was only dyz = 1.0): d: 1.49108, dz2: 0.03198, dxy: 0.04066,
> dx2y2: 0.09713, dxz+dyz: 1.32129
>
> DN (Target was dxy, dz2, dxz, dx2-y2 = 1.0 each): d: 3.70166, dz2:
> 0.72795, dxy: 0.82023, dx2y2: 1.30642, dxz+dyz: 0.84706
>
> 7. Post-processing for Output Density Matrix (lapwdm):
> I used a |.indmc| file same to my |.indm| and ran:
> x lapwdm -up -p
> x lapwdm -dn -p
>
> case.scfdmup (Target was dyz):
> irtest 1 8 2.07000000000000
> Density matrix UPUP block, real part. L= 2
> 0.08807 0.00000 0.00000 0.00000 0.04017
> 0.00000 0.85159 0.00000 0.00000 0.00000
> 0.00000 0.00000 0.04101 0.00000 0.00000
> 0.00000 0.00000 0.00000 0.85159 0.00000
> 0.04017 0.00000 0.00000 0.00000 0.08807
> ...
> :TRA008: TRACE of UPUP MATRIX= 1.92035 0.00000
>
> case.scfdmdn (Target was dxy, dz2, dxz, dx2-y2):
> irtest 1 8 2.07000000000000
> Density matrix UPUP block, real part. L= 2
> 0.94475 0.00000 0.00000 0.00000 0.01576
> 0.00000 0.48222 0.00000 0.00000 0.00000
> 0.00000 0.00000 0.90325 0.00000 0.00000
> 0.00000 0.00000 0.00000 0.48222 0.00000
> 0.01576 0.00000 0.00000 0.00000 0.94475
> ...
> :TRA008: TRACE of UPUP MATRIX= 3.75719 0.00000
>
> Thank you very much for your attention and time!
>
> Best regards,
> Zhiyuan Li
>
>
> _______________________________________________
> Wien mailing list
> Wien at zeus.theochem.tuwien.ac.at
> http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
> SEARCH the MAILING-LIST at: http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
--
Peter Blaha,
Inst. f. Materials Chemistry, TU Vienna, A-1060 Vienna
Email:peter.blaha at tuwien.ac.at <mailto:peter.blaha at tuwien.ac.at>
WIEN2k:http://www.wien2k.at <http://www.wien2k.at/>
More information about the Wien
mailing list