Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Next revision		Previous revision
documentation:tutorials:nio_crystal_field:density_matrix_plot [2016/10/08 21:26] – created Maurits W. Haverkortdocumentation:tutorials:nio_crystal_field:density_matrix_plot [2025/11/20 03:29] (current) – external edit 127.0.0.1
Line 1: Line 1:
 +{{indexmenu_n>11}}
 +====== Density matrix plot ======
  
 +###
 +An intuitive way to look at eigenstates is to plot the charge density. For multi-electron wave-functions one can not plot the wave-function (it depends on $3\times n$ coordinates) but the charge density is a well defined simple quantity. We here use the rendering options of \textit{Mathematica} in combination with the package for \textit{Mathematica} (Quanty.nb) to plot the density matrix of the lowest three eigenstates in NiO.
 +###
 +
 +###
 +The input script is:
 +<code Quanty density_matrix_plots_Mathematica.Quanty>
 +-- It is often nice to make plots of the resulting "wave-functions" Now one can not
 +-- simply plot a multi-electron wavefunction \psi(r_1,r_2,r_3, ..., r_n) as it depends
 +-- on 3n coordinates.
 +-- One can plot the resulting charge density.
 +
 +-- This example calculates the density matrix of a state and then calls Mathematica to
 +-- make a density plot of this matrix.
 +
 +-- the script for mathematica is included in the input file
 +-- the input string needs calculated data which will be substituted before calling mathematica
 +
 +Verbosity(0)
 +
 +-- You need to change the absolute path in the script below
 +mathematicaInput = [[
 +Needs["Quanty`PlotTools`"];
 +rho=%s;
 +pl = Table[ Rasterize[ DensityMatrixPlot[IdentityMatrix[10] - rho[ [i] ], PlotRange -> {{-0.4, 0.4}, {-0.4, 0.4}, {-0.4, 0.4}}] ], {i, 1, Length[rho]}];
 +For[i = 1, i <= Length[pl], i++,
 +    Export["/Users/haverkort/Documents/Quanty/Example_and_Testing/History/current/Tutorials/20_NiO_Crystal_Field/Rho" <> ToString[i] <> ".png", pl[ [i] ] ];
 +   ];
 +Quit[];
 +]]
 +
 +NF=10
 +NB=0
 +dIndexDn={0,2,4,6,8}
 +dIndexUp={1,3,5,7,9}
 +
 +OppSx   =NewOperator("Sx"   ,NF, dIndexUp, dIndexDn)
 +OppSy   =NewOperator("Sy"   ,NF, dIndexUp, dIndexDn)
 +OppSz   =NewOperator("Sz"   ,NF, dIndexUp, dIndexDn)
 +OppSsqr =NewOperator("Ssqr" ,NF, dIndexUp, dIndexDn)
 +OppSplus=NewOperator("Splus",NF, dIndexUp, dIndexDn)
 +OppSmin =NewOperator("Smin" ,NF, dIndexUp, dIndexDn)
 +
 +OppLx   =NewOperator("Lx"   ,NF, dIndexUp, dIndexDn)
 +OppLy   =NewOperator("Ly"   ,NF, dIndexUp, dIndexDn)
 +OppLz   =NewOperator("Lz"   ,NF, dIndexUp, dIndexDn)
 +OppLsqr =NewOperator("Lsqr" ,NF, dIndexUp, dIndexDn)
 +OppLplus=NewOperator("Lplus",NF, dIndexUp, dIndexDn)
 +OppLmin =NewOperator("Lmin" ,NF, dIndexUp, dIndexDn)
 +
 +OppJx   =NewOperator("Jx"   ,NF, dIndexUp, dIndexDn)
 +OppJy   =NewOperator("Jy"   ,NF, dIndexUp, dIndexDn)
 +OppJz   =NewOperator("Jz"   ,NF, dIndexUp, dIndexDn)
 +OppJsqr =NewOperator("Jsqr" ,NF, dIndexUp, dIndexDn)
 +OppJplus=NewOperator("Jplus",NF, dIndexUp, dIndexDn)
 +OppJmin =NewOperator("Jmin" ,NF, dIndexUp, dIndexDn)
 +
 +Oppldots=NewOperator("ldots",NF, dIndexUp, dIndexDn)
 +
 +-- define the coulomb operator
 +-- we here define the part depending on F0 seperately from the part depending on F2
 +-- when summing we can put in the numerical values of the slater integrals
 +OppF0 =NewOperator("U", NF, dIndexUp, dIndexDn, {1,0,0})
 +OppF2 =NewOperator("U", NF, dIndexUp, dIndexDn, {0,1,0})
 +OppF4 =NewOperator("U", NF, dIndexUp, dIndexDn, {0,0,1})
 +
 +-- Akm = {{k1,m1,Akm1},{k2,m2,Akm2}, ... }
 +Akm = PotentialExpandedOnClm("Oh", 2, {0.6,-0.4})
 +OpptenDq = NewOperator("CF", NF, dIndexUp, dIndexDn, Akm)
 +
 +Akm = PotentialExpandedOnClm("Oh", 2, {1,0})
 +OppNeg = NewOperator("CF", NF, dIndexUp, dIndexDn, Akm)
 +Akm = PotentialExpandedOnClm("Oh", 2, {0,1})
 +OppNt2g = NewOperator("CF", NF, dIndexUp, dIndexDn, Akm)
 +
 +Hamiltonian =  10.0 * OppF2 + 6.25 * OppF4 + 2.0 * OpptenDq + 0.06 * Oppldots + 0.000001 * (2*OppSz + OppLz)
 +
 +-- we now can create the lowest Npsi eigenstates:
 +Npsi=3
 +-- in order to make sure we have a filling of 2 electrons we need to define some restrictions
 +StartRestrictions = {NF, NB, {"1111111111",8,8}}
 +
 +psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi)
 +oppList={Hamiltonian, OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz, Oppldots, OppF2, OppF4, OppNeg, OppNt2g}
 +
 +print(" <E>    <S^2>  <L^2>  <J^2>  <S_z>  <L_z>  <l.s>  <F[2]> <F[4]> <Neg>  <Nt2g>");
 +for key,psi in pairs(psiList) do
 +  expvalue = psi * oppList * psi
 +  for k,v in pairs(expvalue) do
 +    io.write(string.format("%6.3f ",v))
 +  end;
 +  io.write("\n")
 +end
 +io:flush()
 +
 +function tableToMathematica(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 + (%18.15f) I ",Complex.Re(v),Complex.Im(v))
 +    end
 +  end
 +  ret = ret.." }"
 +  return ret
 +end
 +
 +rhoList = DensityMatrix(psiList, {0,1,2,3,4,5,6,7,8,9})
 +rhoListMathematicaForm = tableToMathematica(rhoList)
 +
 +file = io.open("densitymatrix.nb", "w")
 +file:write( mathematicaInput:format( rhoListMathematicaForm ) )
 +file:close()
 +
 +os.execute("/Applications/Mathematica08.app/Contents/MacOS/MathKernel -run '<<densitymatrix.nb'")
 +os.execute("convert Rho1.png -transparent white temp.png ; convert temp.png Rho1.eps ; rm temp.png")
 +os.execute("convert Rho2.png -transparent white temp.png ; convert temp.png Rho2.eps ; rm temp.png")
 +os.execute("convert Rho3.png -transparent white temp.png ; convert temp.png Rho3.eps ; rm temp.png")
 +</code>
 +###
 +
 +The produced figures are:
 +^ $S_z=-1$ ^ $S_z=0$ ^ $S_z=1$ ^
 +| {{ :documentation:tutorials:nio_crystal_field:rho1.png?nolink&300 |}} | {{ :documentation:tutorials:nio_crystal_field:rho2.png?nolink&300 |}} | {{ :documentation:tutorials:nio_crystal_field:rho3.png?nolink&300 |}} |
 +^ Density matrix of the ground, first and second excited state of Ni in NiO. ^^^
 +
 +
 +###
 +The plots can be seen in the figure above. We show the hole density (identity matrix minus the density matrix) which has lobes towards the O atoms. Ni in NiO has a hole in the $x^2-y^2$ and in the $z^2$ orbital. The color represent the spin down, $S_z=0$ and spin up state of the lowest triplet.
 +###
 +
 +###
 +The output to the screen is:
 +<file Quanty_Output density_matrix_plots_Mathematica.out>
 + <E>    <S^2>  <L^2>  <J^2>  <S_z>  <L_z>  <l.s>  <F[2]> <F[4]> <Neg>  <Nt2g>
 +-18.093  2.000 12.000 14.471 -0.999 -0.119 -0.147 -1.020 -0.878  2.002  5.998 
 +-18.093  2.000 12.000 14.471 -0.000 -0.000 -0.147 -1.020 -0.878  2.002  5.998 
 +-18.093  2.000 12.000 14.471  0.999  0.119 -0.147 -1.020 -0.878  2.002  5.998 
 +Mathematica 8.0 for Mac OS X x86 (64-bit)
 +Copyright 1988-2011 Wolfram Research, Inc.
 +Loaded the Quanty : Plot Tools Package - version 2016.4.13
 +Written by Maurits W. Haverkort
 +</file>
 +###
 +
 +===== Table of contents =====
 +{{indexmenu>.#1|msort}}

You are here: Quanty - a quantum many body script language » Documentation » Tutorials » NiO crystal field » Density matrix plot

Print/export