NewOperator(“AtomicU”,indices,RadialWavefunctions1,RadialWavefunctions2,options) calculates the full Coulomb interaction operator between all given atomic orbitals.

• indices : indices of atomic orbitals belonging to the radial wave-functions
• RadialWavefunctions1 : table of InterpolatingFunction objects
• RadialWavefunctions2 : optional - table of InterpolatingFunction objects - if given, the function calculates relativistic Slater integrals where RadialWavefunctions1 are the large parts and RadialWavefunctions2 the small ones
• options : possible Options are {“conserve”, bool} and {“1P-Operator”, bool}

  • “conserve” : if true only occupation conserving Slater integrals are calculated
  • “1P-Operator” : only if “conserve” is true and additional arguments are given, transforms interaction to 1-particle operator (see example)

The radial wave-functions can be read from FPLO with ReadFPLOBasisFunctions(), or calculated with CalculateAtomicBasisFunctions(). For the functions of CalculateAtomicBasisFunctions() the Slater integrals are calculated exactly on the B-spline basis of the atom, without a radial grid (see the second example). All functions must then come from the same call.

The example below shows how to create the full Coulomb interaction between multiple atomic orbitals with relativistic wave-functions. It is recommended to use CreateAtomicIndicesList(orbitallabels,{ {“Kappas”,true} }) for index creation, since it automatically calculates the quantum numbers kappa and stores them in the field ind.kappas which needs to be set when calling NewOperator(“AtomicU”,…) with relativistic wave-functions.

Example.Quanty

-- the Coulomb interaction of the electrons of Ne on its own Dirac-Fock orbitals
G, F, x, Atom = CalculateAtomicBasisFunctions(10)
ind, NF = CreateAtomicIndicesList(Atom.Labels, {{"Kappas", true}})
 
U = NewOperator("AU", NF, ind, G, F)
 
psi = NewWavefunction(NF, 0, { {string.rep("1", NF), 1} })
print(string.format("<psi|U|psi> = %.8f Hartree, Edirect + Eexchange = %.8f Hartree", psi*U*psi, Atom.Edirect + Atom.Eexchange))

Example.out

<psi|U|psi> = 54.08660432 Hartree, Edirect + Eexchange = 54.08660432 Hartree

Note: parts of the functionality of this function will be available in the release of October 2019

NewOperator(“AtomicU”,indices,RadialWavefunctions1,RadialWavefunctions2,options) calculates the full Coulomb interaction operator between all given atomic orbitals.

• indices : indices of atomic orbitals belonging to the radial wave-functions
• RadialWavefunctions1 : table of InterpolatingFunction objects
• RadialWavefunctions2 : optional - table of InterpolatingFunction objects - if given, the function calculates relativistic Slater integrals where RadialWavefunctions1 are the large parts and RadialWavefunctions2 the small ones
• options : possible Options are {“conserve”, bool} and {“1P-Operator”, bool}

  • “conserve” : if true only occupation conserving Slater integrals are calculated
  • “1P-Operator” : only if “conserve” is true and additional arguments are given, transforms interaction to 1-particle operator (see example)

The example below shows how to create the full Coulomb interaction between multiple atomic orbitals with relativistic wave-functions. It is recommended to use CreateAtomicIndicesList(orbitallabels,{ {“Kappas”,true} }) for index creation, since it automatically calculates the quantum numbers kappa and stores them in the field ind.kappas which needs to be set when calling NewOperator(“AtomicU”,…) with relativistic wave-functions.

Example.Quanty

orbitallabels = {"1s1/2","2s1/2","2p1/2","2p3/2"}
ind, NF = CreateAtomicIndicesList(orbitallabels,{ {"Kappas",true} })
G, F = ReadFPLOBasisFunctions(orbitallabels,"+fval.001.1", "+fval.001.2")
 
U = NewOperator("AtomicU",NF,ind,G,F) -- creates full Coulomb interaction
 
-- now occupation of the s-orbitals shall be conserved
 
function Take(list,n1,n2)
  local result={}
  if n2== nil then
    if n1 > 0 then
      for i = 1, n1 do
        result[i]=list[i]
      end
    else
      for i=1,-n1 do
        result[i]=list[#list+n1+i]
      end
    end
   else
     for i = n1, n2 do
       result[i+1-n1] = list[i]
     end
   end
  return result
end
 
n0 = 3
n1 = 4
 
indConserved = Take(ind,1,n0-1)
indConserved.kappas = Take(ind.kappas,1,n0-1)
GConserved = Take(G,1,n0-1)
FConserved = Take(F,1,n0-1)
 
indNonCnsrvd = Take(ind,n0,n1)
indNonCnsrvd.kappas = Take(ind.kappas,n0,n1)
GNonCnsrvd = Take(G,n0,n1)
FNonCnsrvd = Take(F,n0,n1)
 
U0 = NewOperator("AU",NF,indConserved,GConserved,FConserved,{ {"conserve",true} }) -- Coulomb interaction between conserved orbitals only
 
U1 = NewOperator("AU",NF,ind,G,F,indConserved.kappas,indNonCnsrvd.kappas,{ {"conserve",true},{"1P-Operator",true} }) -- Coulomb interaction between conserved orbitals and others transformed to a one-particle operator given that the conserved orbitals are always completely filled.
 
U2 = NewOperator("AU",NF,indNonCnsrvd,GNonCnsrvd,FNonCnsrvd) -- Coulomb interaction between non-conserved orbitals
 
Uconserved = U0 + U1 + U2