This is an old revision of the document!
Density matrix plots
asked by Hebatalla Elnaggar (2019/01/31 20:45)
Hi Maurits,
I am trying to plot the GS density matrix for a d5 HS cubic system however I think something goes wrong. Here I compare four calculations (the problem occurs in the fourth case ):
1- Jex parallel to the X axis with no spin-orbit coupling. The first 3 ground-states are: # <E> <S^2> <L^2> <J^2> <Sx> <Lx> <Np> <Nd> <NL>
1 -5.6576 8.7500 -0.0000 8.7500 -2.5000 0.0000 6.0000 5.0000 10.0000 2 -5.5676 8.7500 -0.0000 8.7500 -1.5000 0.0000 6.0000 5.0000 10.0000 3 -5.4776 8.7500 -0.0000 8.7500 -0.5000 0.0000 6.0000 5.0000 10.0000
plotting the 1st GS I get a spherical state fully red as expected projecting along the x-axis
2- Jex parallel to the X axis with 100% spin-orbit coupling. The first 3 ground-states are:
# <E> <S^2> <L^2> <J^2> <Sx> <Lx> <Np> <Nd> <NL> 1 -5.6621 8.7439 0.0030 8.7505 -2.4988 -0.0012 6.0000 5.0000 10.0000 2 -5.5722 8.7437 0.0031 8.7505 -1.4993 -0.0007 6.0000 5.0000 10.0000 3 -5.4822 8.7436 0.0032 8.7506 -0.4998 -0.0002 6.0000 5.0000 10.0000
plotting the 1st GS I get an almost spherical state fully red as expected projecting along the x-axis
3- Jex aligned 30 degrees from the Y axis (rotation about the Z-axis) with no spin-orbit coupling. The first 3 ground-states are:
# <E> <S^2> <L^2> <J^2> <S||> <L||> <Np> <Nd> <NL> 1 -5.6576 8.7500 -0.0000 8.7500 -2.5000 0.0000 6.0000 5.0000 10.0000 2 -5.5676 8.7500 -0.0000 8.7500 -1.5000 0.0000 6.0000 5.0000 10.0000 3 -5.4776 8.7500 -0.0000 8.7500 -0.5000 0.0000 6.0000 5.0000 10.0000
plotting the 1st GS I get an almost spherical state nearly blue as expected for projecting along the || axis
4- Jex aligned 30 degrees from the Y axis (rotation about the Z-axis) with 100% spin-orbit coupling. The first 3 ground-states are: # <E> <S^2> <L^2> <J^2> <S||> <L||> <Np> <Nd> <NL>
1 -5.6621 8.7439 0.0030 8.7505 -2.4988 -0.0012 6.0000 5.0000 10.0000 2 -5.5722 8.7437 0.0031 8.7505 -1.4993 -0.0007 6.0000 5.0000 10.0000 3 -5.4822 8.7436 0.0032 8.7506 -0.4998 -0.0002 6.0000 5.0000 10.0000
plotting the 1st GS I get a very strange non-spherical state with a mixed spin. I do not understand why the plot fails here. The GS is almost identical to the case of 2 but the resulting density matrix is very different.
Below is the script I used:
function TableToMathematica(t)
Chop(t)
local ret = "{ "
for k,v in pairs(t) do
if k~=1 then
ret = ret.." , "
end
if (type(v) == "table") then
ret = ret..TableToMathematica(v)
else
ret = ret..string.format("+ %18.15f ",Complex.Re(v))
ret = ret..string.format("+ I %18.15f ",Complex.Im(v))
end
end
ret = ret.." }"
return ret
end
NBosons = 0 NFermions = 26
NElectrons_2p = 6 NElectrons_3d = 5 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
OppConfNd={} for i=0,10 do
OppConfNd[i] = NewOperator("Identity", NFermions)
OppConfNd[i].Restrictions = {NFermions,NBosons,{"000000 1111111111 0000000000",i,i}}
end
-- Fe3+ --
Delta_sc = 1.5*1 U_3d_3d_sc = 6.5*1 F2_3d_3d_sc = 10.965*0.74 F4_3d_3d_sc = 7.5351*0.74 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) / (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 = 2)
+ 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.059*1 zeta_3d_ic = 0.075*1 zeta_2p_ic = 8.199 zeta_3d_fc = zeta_3d_sc H_soc_sc = zeta_3d_sc * Oppldots_3d H_soc_ic = zeta_3d_ic * Oppldots_3d
+ zeta_2p_ic * Oppldots_2pH_soc_fc = zeta_3d_fc * Oppldots_3d
-- Define the ligand field term. Akm = 14 OpptenDq_3d = NewOperator(“CF”, NFermions, IndexUp_3d, IndexDn_3d, Akm) OpptenDq_Ld = NewOperator(“CF”, NFermions, IndexUp_Ld, IndexDn_Ld, Akm) Akm = 2_0_-7 OppDs_3d = NewOperator(“CF”, NFermions, IndexUp_3d, IndexDn_3d, Akm) OppDs_Ld = NewOperator(“CF”, NFermions, IndexUp_Ld, IndexDn_Ld, Akm) Akm = 4_0_-21 OppDt_3d = NewOperator(“CF”, NFermions, IndexUp_3d, IndexDn_3d, Akm) OppDt_Ld = NewOperator(“CF”, NFermions, IndexUp_Ld, IndexDn_Ld, Akm) Akm = PotentialExpandedOnClm('D4h', 2, {0,0,0,1}) OppVe = NewOperator(“CF”, NFermions, IndexUp_3d, IndexDn_3d, IndexUp_Ld, IndexDn_Ld,Akm) + NewOperator(“CF”, NFermions, IndexUp_Ld, IndexDn_Ld, IndexUp_3d, IndexDn_3d, Akm) Akm = PotentialExpandedOnClm('D4h', 2, {0,0,1,0}) OppVb2 = NewOperator(“CF”, NFermions, IndexUp_3d, IndexDn_3d, IndexUp_Ld, IndexDn_Ld,Akm) + NewOperator(“CF”, NFermions, IndexUp_Ld, IndexDn_Ld, IndexUp_3d, IndexDn_3d, Akm) Akm = PotentialExpandedOnClm('D4h', 2, {1,0,0,0}) OppVa1 = NewOperator(“CF”, NFermions, IndexUp_3d, IndexDn_3d, IndexUp_Ld, IndexDn_Ld,Akm) + NewOperator(“CF”, NFermions, IndexUp_Ld, IndexDn_Ld, IndexUp_3d, IndexDn_3d, Akm) Akm = PotentialExpandedOnClm('D4h', 2, {0,1,0,0}) OppVb1 = NewOperator(“CF”, NFermions, IndexUp_3d, IndexDn_3d, IndexUp_Ld, IndexDn_Ld,Akm) + NewOperator(“CF”, NFermions, IndexUp_Ld, IndexDn_Ld, IndexUp_3d, IndexDn_3d, Akm) Akm = PotentialExpandedOnClm('D4h', 2, {0,0,0,1}) OppNe = NewOperator(“CF”, NFermions, IndexUp_3d, IndexDn_3d, Akm) Akm = PotentialExpandedOnClm('D4h', 2, {0,0,1,0}) OppNb2 = NewOperator(“CF”, NFermions, IndexUp_3d, IndexDn_3d,Akm) Akm = PotentialExpandedOnClm('D4h', 2, {1,0,0,0}) OppNa1 = NewOperator(“CF”, NFermions, IndexUp_3d, IndexDn_3d,Akm) Akm = PotentialExpandedOnClm('D4h', 2, {0,1,0,0}) OppNb1 = NewOperator(“CF”, NFermions, IndexUp_3d, IndexDn_3d, Akm) Ds_3d = 0.0 Dt_3d = Ds_3d*0.15 tenDq_3d_sc = -0.6 + Dt_3d*10*7/12 tenDq_Ld_sc = tenDq_3d_sc/3--0.0 Veg_sc = 3.2*0 Vt2g_sc = 1.705 *0 Ve = Vt2g_sc Vb2 = Vt2g_sc Va1 = Veg_sc /1000 Vb1 = Veg_sc tenDq_3d_ic = tenDq_3d_sc tenDq_Ld_ic = tenDq_Ld_sc Veg_ic = Veg_sc Vt2g_ic = Vt2g_sc tenDq_3d_fc = tenDq_3d_sc tenDq_Ld_fc = tenDq_Ld_sc Veg_fc = Veg_sc Vt2g_fc = Vt2g_sc H_lf_sc = tenDq_3d_sc * OpptenDq_3d
+ tenDq_Ld_sc * OpptenDq_Ld
+ Ds_3d * OppDs_3d
+ Dt_3d * OppDt_3d
+ Ve * OppVe
+ Vb2 * OppVb2
+ Va1 * OppVa1
+ Vb1 * OppVb1
-- + Veg_sc * OppVeg
-- + Vt2g_sc * OppVt2g
H_lf_ic = tenDq_3d_ic * OpptenDq_3d + tenDq_Ld_ic * OpptenDq_Ld
+ Ds_3d * OppDs_3d
+ Dt_3d * OppDt_3d
+ Ve * OppVe
+ Vb2 * OppVb2
+ Va1 * OppVa1
+ Vb1 * OppVb1
-- + Veg_ic * OppVeg
-- + Vt2g_ic * OppVt2g
H_lf_fc = tenDq_3d_fc * OpptenDq_3d + tenDq_Ld_fc * OpptenDq_Ld
+ Ds_3d * OppDs_3d
+ Dt_3d * OppDt_3d
+ Ve * OppVe
+ Vb2 * OppVb2
+ Va1 * OppVa1
+ Vb1 * OppVb1
-- + 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 Jvec= {1,0,0} J0 = 1e-3 * 90 B0 = 1e-4 Bx = Jvec[1]*B0 By = Jvec[2]*B0 Bz = Jvec[3]*B0 Jx = Jvec[1]*J0 Jy = Jvec[2]*J0 Jz = Jvec[3]*J0 Jex = Jx * OppSx
+ Jy * OppSy + Jz * OppSzB = Bx * ( OppLx)
+ By * ( OppLy) + Bz * ( OppLz)
-- Compose the total Hamiltonian.
H_sc = 1 * H_coulomb_sc + 1 * H_soc_sc + 1 * H_lf_sc + B + Jex H_ic = 1 * H_coulomb_ic + 1 * H_soc_ic + 1 * H_lf_ic + B + Jex H_fc = 1 * H_coulomb_fc + 1 * H_soc_fc + 1 * H_lf_fc + B + Jex 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}}
NPsis = 6
Restrictions = {NFermions, NBosons, {“000000 0000000000 1111111111”,10,10}}
Psis = Eigensystem(H_sc, StartingRestrictions, NPsis,restrictions_restrictions)
if not (type(Psis) == 'table') then Psis = {Psis}
end
-- Plotting mathematicaInput = [[
Needs["Quanty`PlotTools`"];
rho=%s;
pl = Table[ Rasterize[ DensityMatrixPlot[ rho[ [i] ],QuantizationAxes->"x", PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}] ], {i, 1, Length[rho]}];
For[i = 1, i <= Length[pl], i++,
Export[",." <> ToString[i] <> ".png", pl[ [i] ] ];
];
Quit[];
]]
-- Plotting density plots:
rhoList1 = DensityMatrix(Psis, {6,7,8,9,10,11,12,13,14,15})
rhoListMathematicaForm1 = TableToMathematica(rhoList1)
file = io.open("./Densitymatrix2.nb", "w")
file:write( mathematicaInput:format(rhoListMathematicaForm1 ) )
file:close()
--os.execute("/Applications/Mathematica.app/Contents/MacOS/MathKernel -run '<<".."Densitymatrix1.nb'")
print('Finished the density matrix')
-- Print some useful information about the calculated states.
file = io.open(“Expect2.txt”, “w”);
OppSpar= Jvec[3]*OppSz_3d+Jvec[2]*OppSy_3d+Jvec[1]*OppSx_3d
OppLpar= Jvec[3]*OppLz_3d+Jvec[2]*OppLy_3d+Jvec[1]*OppLx_3d
OppList = {H_sc, OppSsqr, OppLsqr, OppJsqr, OppSpar, OppLpar, OppN_2p, OppN_3d, OppN_Ld}
ConfNds={OppConfNd[6], OppConfNd[7], OppConfNd[8], OppConfNd[9], OppConfNd[10]}
Psitemp={}
print(' # <E> <S^2> <L^2> <J^2> <S||> <L||> <Np> <Nd> <NL>');
file:write(' # <E> <S^2> <L^2> <J^2> <Sz> <Lz> <Np> <Nd> <NL>');
file:write('\n')
for key, Psi in pairs(Psis) doexpectationValues = Psi * OppList * Psi
file:write(string.format('%3d', key))
for key, expectationValue in pairs(expectationValues) do
io.write(string.format('%9.4f', Complex.Re(expectationValue)))
file:write(string.format('%9.4f', Complex.Re(expectationValue)))
--io.write(string.format('%9.4f', Complex.Re(expectationValue)))
end
-- for k=6,10 do
-- Psitemp = OppConfNd[k] * Psi
-- ConfNdValue = Psitemp * OppConfNd[k] * Psitemp
-- file:write(string.format('%9.4f', ConfNdValue))
-- end
io.write('\n')
file:write('\n')
end
file:close()
OppList = {H_sc, OppNa1, OppNb1, OppNe, OppNb2}
print(' # <E> <Na1> <Nb1> <Ne> <Nb2> ');
for key, Psi in pairs(Psis) doexpectationValues = Psi * OppList * Psi
for key, expectationValue in pairs(expectationValues) do
io.write(string.format('%9.4f', Complex.Re(expectationValue)))
end
io.write('\n')
end
os.exit()
</WRAP>
~~DISCUSSION|Answers~~
Answers
Dear Heba,
Your function TableToMathematica(t) can not handle complex numbers with negative imaginary part. In that case 1 - 2 I becomes 1 + I - 2. Below you find a correct version that should solve all your problems.
Maurits
function TableToMathematica(t) Chop(t) local ret = "{ " for k,v in pairs(t) do if k~=1 then ret = ret.." , " end if (type(v) == "table") then ret = ret..TableToMathematica(v) else if( Complex.Re(v) < 0) then ret = ret..string.format("- %18.15f ",Abs(Complex.Re(v))) else ret = ret..string.format("+ %18.15f ",Abs(Complex.Re(v))) end if( Complex.Im(v) < 0) then ret = ret..string.format("- I %18.15f ",Abs(Complex.Im(v))) else ret = ret..string.format("+ I %18.15f ",Abs(Complex.Im(v))) end end end ret = ret.." }" return ret endDear Maurits,
I see. That indeed solved the problem.
Thanks!