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
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
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:
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:
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.