Wrong (missing) ground state

asked by Hebatalla Elnaggar (2016/11/02 06:33)

Operating system: All

Version:

  • Linux Jun 27 2016
  • Windows Jun 27 2016
  • Mac Jun 28 2016

Error: Missing (ground) state. The calculation doesn't find the correct ground state if we calculate the EigenSystem() with small number of Npsi.

Calculation: d6, Ligand field close to high-low spin transition

Output:

Npsi=1

SingleState.out
Start of BlockGroundState. Converge 1 states to an energy with relative variance smaller than  1.490116119384766E-06
 
Start of BlockOperatorPsiSerialRestricted
Outer loop   1, Number of Determinants:       209       210 last variance  6.779798109699903E+00
Start of BlockOperatorPsiSerialRestricted
Start of BlockGroundState. Converge 1 states to an energy with relative variance smaller than  1.490116119384766E-06
 
Start of BlockOperatorPsiSerial
Outer loop   1, Number of Determinants:       134       784 last variance  2.249384186944060E+01
Start of BlockOperatorPsiSerial
Outer loop   2, Number of Determinants:       282       718 last variance  8.959539344379486E+00
Start of BlockOperatorPsiSerial
Outer loop   3, Number of Determinants:       718      1102 last variance  1.973600950704991E+00
  Restart loop 1 with a Krylov basis of 101 and a full basis of 1102
Start of BlockOperatorPsiSerial
Outer loop   4, Number of Determinants:      1102      1218 last variance  7.338916467507772E-02
  Restart loop 1 with a Krylov basis of 101 and a full basis of 1218
Start of BlockOperatorPsiSerial
 
  <E>    <S^2>    <L^2>    <J^2>    <Sz>     <Lz>     <Np>      <Nd>     <NL>
-8.2124   2.9309  14.8408  20.3187  -0.0000  -0.0000   6.0000   6.6093   9.3907

Npsi=5

FiveStates.out
Start of BlockGroundState. Converge 5 states to an energy with relative variance smaller than  1.490116119384766E-06
 
Start of BlockOperatorPsiSerialRestricted
Outer loop   1, Number of Determinants:       209       210 last variance  6.320029384998617E+00
Start of BlockOperatorPsiSerialRestricted
Start of BlockGroundState. Converge 5 states to an energy with relative variance smaller than  1.490116119384766E-06
 
Start of BlockOperatorPsiSerial
Outer loop   1, Number of Determinants:       210      1162 last variance  2.249262221674925E+01
 Restart loop 1 with a Krylov basis of 105 and a full basis of 1162
 Restart loop 2 with a Krylov basis of 120 and a full basis of 1162
Start of BlockOperatorPsiSerial
Outer loop   2, Number of Determinants:      1162      2896 last variance  9.026020001225831E+00
 Restart loop 1 with a Krylov basis of 105 and a full basis of 2896
 Restart loop 2 with a Krylov basis of 120 and a full basis of 2896
 Restart loop 3 with a Krylov basis of 135 and a full basis of 2896
Start of BlockOperatorPsiSerial
Outer loop   3, Number of Determinants:      2896      4387 last variance  1.834784704827484E+00
 Restart loop 1 with a Krylov basis of 105 and a full basis of 4387
 Restart loop 2 with a Krylov basis of 120 and a full basis of 4387
Start of BlockOperatorPsiSerial
Outer loop   4, Number of Determinants:      4386      4845 last variance  7.780047778494747E-02
 Restart loop 1 with a Krylov basis of 105 and a full basis of 4845
Start of BlockOperatorPsiSerial
 
    <E>    <S^2>    <L^2>    <J^2>    <Sz>     <Lz>     <Np>      <Nd>     <NL>
 -9.1161   0.0349  25.5821  25.5850  -0.0000  -0.0000   6.0000   6.7370   9.2630
 -8.2124   2.9309  14.8407  20.3194  -0.0000  -0.0000   6.0000   6.6093   9.3907
 -8.2124   2.9311  14.8400  20.3193   0.0000   0.0000   6.0000   6.6092   9.3908
 -8.2075   2.7520  15.3105  19.1875   0.0000  -0.0000   6.0000   6.6172   9.3828
 -8.2075   2.7519  15.3104  19.1879   0.4583   0.0858   6.0000   6.6172   9.3828

Script: (adapted from Crispy)

d6Groundstate.Quanty
--------------------------------------------------------------------------------
-- Define the number of electrons, shells, etc.
--------------------------------------------------------------------------------
 
-- Co3+ --
 
NBosons = 0
NFermions = 26
 
NElectrons_2p = 6
NElectrons_3d = 6
NElectrons_Ld = 10
 
IndexDn_2p = {0, 2, 4}
IndexUp_2p = {1, 3, 5}
IndexDn_3d = {6, 8, 10, 12, 14}
IndexUp_3d = {7, 9, 11, 13, 15}
IndexDn_Ld = {16, 18, 20, 22, 24}
IndexUp_Ld = {17, 19, 21, 23, 25}
 
 
--------------------------------------------------------------------------------
-- Define the Coulomb term.
--------------------------------------------------------------------------------
OppF0_3d_3d = NewOperator('U', NFermions, IndexUp_3d, IndexDn_3d, {1, 0, 0})
OppF2_3d_3d = NewOperator('U', NFermions, IndexUp_3d, IndexDn_3d, {0, 1, 0})
OppF4_3d_3d = NewOperator('U', NFermions, IndexUp_3d, IndexDn_3d, {0, 0, 1})
 
OppF0_2p_3d = NewOperator('U', NFermions, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {1, 0}, {0, 0})
OppF2_2p_3d = NewOperator('U', NFermions, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0, 1}, {0, 0})
OppG1_2p_3d = NewOperator('U', NFermions, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0, 0}, {1, 0})
OppG3_2p_3d = NewOperator('U', NFermions, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0, 0}, {0, 1})
 
OppNUp_2p = NewOperator('Number', NFermions, IndexUp_2p, IndexUp_2p, {1, 1, 1})
OppNDn_2p = NewOperator('Number', NFermions, IndexDn_2p, IndexDn_2p, {1, 1, 1})
OppN_2p   = OppNUp_2p + OppNDn_2p
 
OppNUp_3d = NewOperator('Number', NFermions, IndexUp_3d, IndexUp_3d, {1, 1, 1, 1, 1})
OppNDn_3d = NewOperator('Number', NFermions, IndexDn_3d, IndexDn_3d, {1, 1, 1, 1, 1})
OppN_3d   = OppNUp_3d + OppNDn_3d
 
OppNUp_Ld = NewOperator('Number', NFermions, IndexUp_Ld, IndexUp_Ld, {1, 1, 1, 1, 1})
OppNDn_Ld = NewOperator('Number', NFermions, IndexDn_Ld, IndexDn_Ld, {1, 1, 1, 1, 1})
OppN_Ld   = OppNUp_Ld + OppNDn_Ld
 
-- Co3+ --
 
Delta_sc    =          2.0
U_3d_3d_sc  =          5.0
F2_3d_3d_sc =       12.663*0.8
F4_3d_3d_sc =        7.917*0.8
F0_3d_3d_sc = U_3d_3d_sc + 2 / 63 * F2_3d_3d_sc + 2 / 63 * F4_3d_3d_sc
e_3d_sc     = (10 * Delta_sc - NElectrons_3d * (19 + NElectrons_3d) * U_3d_3d_sc / 2) / (10 + NElectrons_3d)
e_Ld_sc     = NElectrons_3d * ((1 + NElectrons_3d) * U_3d_3d_sc / 2 - Delta_sc) / (10 + NElectrons_3d)
 
Delta_ic    =          2.0
U_3d_3d_ic  =          5.0
F2_3d_3d_ic =       13.422*0.8
F4_3d_3d_ic =        8.395*0.8
F0_3d_3d_ic = U_3d_3d_ic + 2 / 63 * F2_3d_3d_ic + 2 / 63 * F4_3d_3d_ic
U_2p_3d_ic  =          7.0
F2_2p_3d_ic =        7.900*0.8
G1_2p_3d_ic =        5.951*0.8
G3_2p_3d_ic =        3.386*0.8
F0_2p_3d_ic = U_2p_3d_ic + 1 / 15 * G1_2p_3d_ic + 3 / 70 * G3_2p_3d_ic
e_2p_ic    = (10 * Delta_ic + (1 + NElectrons_3d) * (NElectrons_3d * U_3d_3d_ic / 2 - (10 + NElectrons_3d) * U_2p_3d_ic)) / (16 + NElectrons_3d)
e_3d_ic    = (10 * Delta_ic - NElectrons_3d * (31 + NElectrons_3d) * U_3d_3d_ic / 2 - 90 * U_2p_3d_ic) / (16 + NElectrons_3d)
e_Ld_ic    = ((1 + NElectrons_3d) * (NElectrons_3d * U_3d_3d_ic / 2 + 6 * U_2p_3d_ic) - (6 + NElectrons_3d) * Delta_ic) / (16 + NElectrons_3d)
 
Delta_fc    =          2.0
U_3d_3d_fc  =          5.0
F2_3d_3d_fc =       12.663*0.8
F4_3d_3d_fc =        7.917*0.8
F0_3d_3d_fc = U_3d_3d_fc + 2 / 63 * F2_3d_3d_fc + 2 / 63 * F4_3d_3d_fc
e_3d_fc     = (10 * Delta_fc - NElectrons_3d * (19 + NElectrons_3d) * U_3d_3d_fc / 2) / (10 + NElectrons_3d)
e_Ld_fc     = NElectrons_3d * ((1 + NElectrons_3d) * U_3d_3d_fc / 2 - Delta_fc) / (10 + NElectrons_3d)
 
H_coulomb_sc = F0_3d_3d_sc * OppF0_3d_3d
            + F2_3d_3d_sc * OppF2_3d_3d
            + F4_3d_3d_sc * OppF4_3d_3d
            + e_3d_sc     * OppN_3d
            + e_Ld_sc     * OppN_Ld
 
H_coulomb_ic = F0_3d_3d_ic * OppF0_3d_3d
            + F2_3d_3d_ic * OppF2_3d_3d
            + F4_3d_3d_ic * OppF4_3d_3d
            + F0_2p_3d_ic * OppF0_2p_3d
            + F2_2p_3d_ic * OppF2_2p_3d
            + G1_2p_3d_ic * OppG1_2p_3d
            + G3_2p_3d_ic * OppG3_2p_3d
            + e_2p_ic     * OppN_2p
            + e_3d_ic     * OppN_3d
            + e_Ld_ic     * OppN_Ld
 
H_coulomb_fc = F0_3d_3d_fc * OppF0_3d_3d
            + F2_3d_3d_fc * OppF2_3d_3d
            + F4_3d_3d_fc * OppF4_3d_3d
            + e_3d_fc     * OppN_3d
            + e_Ld_fc     * OppN_Ld
 
--------------------------------------------------------------------------------
-- Define the spin-orbit coupling term.
--------------------------------------------------------------------------------
Oppldots_3d = NewOperator('ldots', NFermions, IndexUp_3d, IndexDn_3d)
 
Oppldots_2p = NewOperator('ldots', NFermions, IndexUp_2p, IndexDn_2p)
 
zeta_3d_sc =     0.066
 
zeta_3d_ic =     0.083
zeta_2p_ic =     9.74738
 
zeta_3d_fc =     0.066
 
H_soc_sc = zeta_3d_sc * Oppldots_3d
 
H_soc_ic = zeta_3d_ic * Oppldots_3d
        + zeta_2p_ic * Oppldots_2p
 
H_soc_fc = zeta_3d_fc * Oppldots_3d
 
--------------------------------------------------------------------------------
-- Define the ligand field term.
--------------------------------------------------------------------------------
OpptenDq_3d = NewOperator('CF', NFermions, IndexUp_3d, IndexDn_3d, PotentialExpandedOnClm('Oh', 2, {0.6, -0.4}))
OpptenDq_Ld = NewOperator('CF', NFermions, IndexUp_Ld, IndexDn_Ld, PotentialExpandedOnClm('Oh', 2, {0.6, -0.4}))
 
OppVeg_3d = NewOperator('CF', NFermions, IndexUp_Ld, IndexDn_Ld, IndexUp_3d, IndexDn_3d, PotentialExpandedOnClm('Oh', 2, {1, 0}))
OppVeg_Ld = NewOperator('CF', NFermions, IndexUp_3d, IndexDn_3d, IndexUp_Ld, IndexDn_Ld, PotentialExpandedOnClm('Oh', 2, {1, 0}))
OppVeg = OppVeg_3d + OppVeg_Ld
 
OppVt2g_3d = NewOperator('CF', NFermions, IndexUp_Ld, IndexDn_Ld, IndexUp_3d, IndexDn_3d, PotentialExpandedOnClm('Oh', 2, {0, 1}))
OppVt2g_Ld = NewOperator('CF', NFermions, IndexUp_3d, IndexDn_3d, IndexUp_Ld, IndexDn_Ld, PotentialExpandedOnClm('Oh', 2, {0, 1}))
OppVt2g = OppVt2g_3d + OppVt2g_Ld
 
tenDq_3d_sc =          1.2
tenDq_Ld_sc =          0.0
Veg_sc      =          3.0
Vt2g_sc     =          1.5
 
tenDq_3d_ic =          1.2
tenDq_Ld_ic =          0.0
Veg_ic      =          3.0
Vt2g_ic     =          1.5
 
tenDq_3d_fc =          1.2
tenDq_Ld_fc =          0.0
Veg_fc      =          3.0
Vt2g_fc     =          1.5
 
H_lf_sc = tenDq_3d_sc * OpptenDq_3d
       + tenDq_Ld_sc * OpptenDq_Ld
       + Veg_sc      * OppVeg
       + Vt2g_sc     * OppVt2g
 
H_lf_ic = tenDq_3d_ic * OpptenDq_3d
       + tenDq_Ld_ic * OpptenDq_Ld
       + Veg_ic      * OppVeg
       + Vt2g_ic     * OppVt2g
 
H_lf_fc = tenDq_3d_fc * OpptenDq_3d
       + tenDq_Ld_fc * OpptenDq_Ld
       + Veg_fc      * OppVeg
       + Vt2g_fc     * OppVt2g
 
--------------------------------------------------------------------------------
-- Define the magnetic field term.
--------------------------------------------------------------------------------
OppSx_3d    = NewOperator('Sx'   , NFermions, IndexUp_3d, IndexDn_3d)
OppSy_3d    = NewOperator('Sy'   , NFermions, IndexUp_3d, IndexDn_3d)
OppSz_3d    = NewOperator('Sz'   , NFermions, IndexUp_3d, IndexDn_3d)
OppSsqr_3d  = NewOperator('Ssqr' , NFermions, IndexUp_3d, IndexDn_3d)
OppSplus_3d = NewOperator('Splus', NFermions, IndexUp_3d, IndexDn_3d)
OppSmin_3d  = NewOperator('Smin' , NFermions, IndexUp_3d, IndexDn_3d)
 
OppLx_3d    = NewOperator('Lx'   , NFermions, IndexUp_3d, IndexDn_3d)
OppLy_3d    = NewOperator('Ly'   , NFermions, IndexUp_3d, IndexDn_3d)
OppLz_3d    = NewOperator('Lz'   , NFermions, IndexUp_3d, IndexDn_3d)
OppLsqr_3d  = NewOperator('Lsqr' , NFermions, IndexUp_3d, IndexDn_3d)
OppLplus_3d = NewOperator('Lplus', NFermions, IndexUp_3d, IndexDn_3d)
OppLmin_3d  = NewOperator('Lmin' , NFermions, IndexUp_3d, IndexDn_3d)
 
OppJx_3d    = NewOperator('Jx'   , NFermions, IndexUp_3d, IndexDn_3d)
OppJy_3d    = NewOperator('Jy'   , NFermions, IndexUp_3d, IndexDn_3d)
OppJz_3d    = NewOperator('Jz'   , NFermions, IndexUp_3d, IndexDn_3d)
OppJsqr_3d  = NewOperator('Jsqr' , NFermions, IndexUp_3d, IndexDn_3d)
OppJplus_3d = NewOperator('Jplus', NFermions, IndexUp_3d, IndexDn_3d)
OppJmin_3d  = NewOperator('Jmin' , NFermions, IndexUp_3d, IndexDn_3d)
 
OppSx_Ld    = NewOperator('Sx'   , NFermions, IndexUp_Ld, IndexDn_Ld)
OppSy_Ld    = NewOperator('Sy'   , NFermions, IndexUp_Ld, IndexDn_Ld)
OppSz_Ld    = NewOperator('Sz'   , NFermions, IndexUp_Ld, IndexDn_Ld)
OppSsqr_Ld  = NewOperator('Ssqr' , NFermions, IndexUp_Ld, IndexDn_Ld)
OppSplus_Ld = NewOperator('Splus', NFermions, IndexUp_Ld, IndexDn_Ld)
OppSmin_Ld  = NewOperator('Smin' , NFermions, IndexUp_Ld, IndexDn_Ld)
 
OppLx_Ld    = NewOperator('Lx'   , NFermions, IndexUp_Ld, IndexDn_Ld)
OppLy_Ld    = NewOperator('Ly'   , NFermions, IndexUp_Ld, IndexDn_Ld)
OppLz_Ld    = NewOperator('Lz'   , NFermions, IndexUp_Ld, IndexDn_Ld)
OppLsqr_Ld  = NewOperator('Lsqr' , NFermions, IndexUp_Ld, IndexDn_Ld)
OppLplus_Ld = NewOperator('Lplus', NFermions, IndexUp_Ld, IndexDn_Ld)
OppLmin_Ld  = NewOperator('Lmin' , NFermions, IndexUp_Ld, IndexDn_Ld)
 
OppJx_Ld    = NewOperator('Jx'   , NFermions, IndexUp_Ld, IndexDn_Ld)
OppJy_Ld    = NewOperator('Jy'   , NFermions, IndexUp_Ld, IndexDn_Ld)
OppJz_Ld    = NewOperator('Jz'   , NFermions, IndexUp_Ld, IndexDn_Ld)
OppJsqr_Ld  = NewOperator('Jsqr' , NFermions, IndexUp_Ld, IndexDn_Ld)
OppJplus_Ld = NewOperator('Jplus', NFermions, IndexUp_Ld, IndexDn_Ld)
OppJmin_Ld  = NewOperator('Jmin' , NFermions, IndexUp_Ld, IndexDn_Ld)
 
OppSx   = OppSx_3d + OppSx_Ld
OppSy   = OppSy_3d + OppSy_Ld
OppSz   = OppSz_3d + OppSz_Ld
OppSsqr = OppSx * OppSx + OppSy * OppSy + OppSz * OppSz
 
OppLx   = OppLx_3d + OppLx_Ld
OppLy   = OppLy_3d + OppLy_Ld
OppLz   = OppLz_3d + OppLz_Ld
OppLsqr = OppLx * OppLx + OppLy * OppLy + OppLz * OppLz
 
OppJx   = OppJx_3d + OppJx_Ld
OppJy   = OppJy_3d + OppJy_Ld
OppJz   = OppJz_3d + OppJz_Ld
OppJsqr = OppJx * OppJx + OppJy * OppJy + OppJz * OppJz
 
Bx = 0 * EnergyUnits.Tesla.value
By = 0 * EnergyUnits.Tesla.value
Bz = 1e-6 * EnergyUnits.Tesla.value
 
B = Bx * (2 * OppSx + OppLx)
 + By * (2 * OppSy + OppLy)
 + Bz * (2 * OppSz + OppLz)
 
--------------------------------------------------------------------------------
-- Compose the total Hamiltonian.
--------------------------------------------------------------------------------
H_sc = 1 * H_coulomb_sc + 1 * H_soc_sc + 1 * H_lf_sc + B
H_ic = 1 * H_coulomb_ic + 1 * H_soc_ic + 1 * H_lf_ic + B
H_fc = 1 * H_coulomb_fc + 1 * H_soc_fc + 1 * H_lf_fc + B
 
----
--H_sc.Chop()
--H_ic.Chop()
--H_fc.Chop()
 
--------------------------------------------------------------------------------
-- Define the starting restrictions and set the number of initial states.
--------------------------------------------------------------------------------
StartingRestrictions = {NFermions, NBosons, {'111111 0000000000 0000000000', NElectrons_2p, NElectrons_2p},
                                           {'000000 1111111111 0000000000', NElectrons_3d, NElectrons_3d},
                                           {'000000 0000000000 1111111111', NElectrons_Ld, NElectrons_Ld}}
-- StartingRestrictions = {NFermions, NBosons, {"000000 1111111111 0000000000",6,6},{"111111 0000000000 1111111111",16,16}};
 
 
NPsis = 20
 
Psis = Eigensystem(H_sc, StartingRestrictions, NPsis)
if not (type(Psis) == 'table') then
   Psis = {Psis}
end
 
-- Print some useful information about the calculated states.
OppList = {H_sc, OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz, OppN_2p, OppN_3d, OppN_Ld}
 
print('     <E>    <S^2>    <L^2>    <J^2>    <Sz>     <Lz>     <Np>      <Nd>     <NL>');
for key, Psi in pairs(Psis) do
	expectationValues = Psi * OppList * Psi
	for key, expectationValue in pairs(expectationValues) do
		io.write(string.format('%9.4f', expectationValue))
		--io.write(string.format('%9.4f', Complex.Re(expectationValue)))
	end
	io.write('\n')
end

Answers

, 2016/11/02 14:04, 2016/11/02 14:10

Dear Heba,

Yes this is a well known “problem” or maybe I should even say feature of the method. As one uses iterative schemes to optimize the basis and diagonalize the Hamiltonian one can get trapped into local minima. The d6 case is a notorious case that has two low energy states (low-spin and high-spin) that are only very weakly coupled.

As I get this question basically after every talk I give, and as it is really not a problem once you understand the reason behind it let me elaborate a bit on what the code does.

In the line

Psis = Eigensystem(H_sc, StartingRestrictions, NPsis)

you start a Lanczos algorithm to find the lowest Npsis eigen-states of H_sc. Lanczos works on the principle that $H_{sc} |\psi\rangle$ is closer to the ground-state than $|\psi\rangle$ for an arbitrary, random state $\psi$. In order to start the Lanczos cycle we need to define a starting - random - state $\psi$. This sad one should realize that you are not looking at the lowest eigen-state of $H_{sc}$. You want the lowest eigen-state with a filled core shell and in total 16 electrons in the d plus Ligand shell. We thus set restrictions in order to achieve this. As The Hamiltonian commutes with particle number one only needs to set these restrictions on the random starting state.

The following restrictions would suit your needs:

StartingRestrictions = {NFermions, NBosons, {'111111 0000000000 0000000000', 6, 6},
                                            {'000000 1111111111 1111111111', 16, 16}}

Using those restrictions the program will find the correct lowest ground-state independently of the value of the cubic crystal field.

The starting restrictions you used are:

StartingRestrictions = {NFermions, NBosons, {'111111 0000000000 0000000000', NElectrons_2p, NElectrons_2p},
                                           {'000000 1111111111 0000000000', NElectrons_3d, NElectrons_3d},
                                           {'000000 0000000000 1111111111', NElectrons_Ld, NElectrons_Ld}}

You start from a random state that has no covalent bonding, i.e. perfectly 6 electrons on the Co atom and all O 2p orbitals completely filled. For systems where the eigen-state without covalence has the same symmetry as the state with covalence this will speed up the convergence. For systems where the symmetry is different it will either slow the calculation down, or lead to the “wrong” state to be calculated (namely an excited state). For ligand field theory this means that if you know that the system is in a Hunds-rule high-spin ground-state you want to start from the restrictions that leave the ligand orbitals completely filled. For low-spin states you want to start from the restrictions that include covalence, unless the local crystal-field is so large that the ionic part of the Hamiltonian already gives you a low-spin state. If the number of states you want to calculate is large enough that both are included you are again fastest by starting with the more restricted restrictions.

In practise this is not a problem. For a d6 system with a low-spin and high-spin state close by you probably want to calculate at least 16 (15 High-spin + 1 Low-spin) states anyhow.

You could leave a comment if you were logged in.
Print/export