Coulomb Operator

The Coulomb interaction is given by $\frac{1}{r_i - r_j}$. In order to make this operator calculable we expand this function on spherical harmonics. See the documentation for more information on how and why we do this.

After such an expansion the angular part can be solved analytical. The radial part results in integrals that depend on the radial functions and the multipole moment of the expansion on spherical harmonics. For the interaction between two $d$ electrons we have a monopole ($k=0$), quadrupole ($k=2$) and $k=4$ interaction.

Within a single shell these multi-pole dependent integrals are called the Slater integrals ($F[k]$). The function NewOperator(“U”, NF, IndexUp, IndexDn, {F[0], F[2], F[4]}) takes these Slater integrals as an input and you provide them as a list of allowed integrals.

See the documentation of the standard operators for more information.

The example in the tutorials does not calculate the multiplets (Eigen-states of the Coulomb operator) for a real system. Without radial wave-functions we take dummy values for the Slater integrals. It turns out that within crystal field theory $F[0]$ is unimportant as the number of $d$ electrons is constant. $F[2]$ is around 10 for the transition metals and $F[4]$ is roughly 6.25 for the transition metals. The values taken here roughly (but indeed only very roughly) follow these values. Good enough for a simple example, but lacking material specificity.

Best wishes, Maurits

Thanks a lot for your reply

Cheers Riccardo

Still referring to the program that calculates the multiplet terms, why if I put in conf.in l=1 and Nelec=6 instead of giving me the 1S_0 term only, the program gives the following error: Eigenstates_of_Coulomb.Quanty 141: Error operation “.1()” not defined on wavefunction

Thanks in advance

Riccardo. This happens because for l=1 and Nelec=6 is a filled shell and there is only one Eigenvector (which is generate by the function EigenSystem). Therefore, the variable psiEigen is a “wavefunction”. For other cases, where there are more then one Eigenvector, the variable psiEigen is a table of Eigenvectors (each table entry is a eigenvector).

To circumvent this you could make the variable psiEigen a table by adding this piece of code:

psiEigen = Eigensystem(Ham, Restrictions, Nstates);
-- Add these lines:
if (type(psiEigen)=="userdata") then
  psi_temporary = psiEigen
  psiEigen = {}
  psiEigen[1] = psi_temporary
end

I picked up this tutorial again because there were things I didn't understand. In analyzing it I came across this equation L = floor(0.5 * (sqrt(abs(L2) * 4 +1) - 1) + 0.5) Where does it come from?