asked by Hebatalla Elnaggar (2016/11/02 06:33)
Operating system: All
Version:
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
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
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)
--------------------------------------------------------------------------------
-- 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