Fermi energy. We demonstrate the function for computing Fermi energy of an arbitrary set of energy levels with an arbitrary number of electrons occupying some of these levels. For the purpose of demonstration, we set the energy levels in the list "e" with the 3 elements. Each element represents the energy level. In our case, they are chosen to be -1.0, -0.5, and -0.4 a.u. (Well, the absolute units of energy do not matter here, as long as the kT is chosen in the same units of energy). We then define some essential parameters: Nel = 2 - the number of electrons occupying these levels (well, it must be larger than 0 and not larger than the number of energy levels), degen = 1 - the degeneracy of each energy level - what is the maximal number of electrons can be placed on each orbital. The value 1 implies we consider only alpha or only beta electrons/orbitals. kT = 0.025 is the parameter controlling the spread of the Fermi distribution. The smaller this parameter is, the steeper the change of populations of the frontier orbitals. In the limit of kT --> 0, one would get integer occupations, but that may not be what one wants (e.g because integer occupations could slow down or destroy the convergence of SCF iterations). Finally, etol = 0.0001 - is the parameter controlling the accuracy of the computed Fermi energy. Once the parameters are set (or you could simply plug the numbers into appropriate positions when calling the function), we are ready to compute the Fermi energy:
Ef = fermi_energy(e, Nel, degen, kT, etol)
Ordering bands. Next, simple, but useful function showcased in this tutorial is the "order_bands". As it is apparent from its name, the function orders the energy levels. In this case, it takes a diagonal matrix (MATRIX class) as an argument and returns a list of 2-element lists. Each element of the outer list is a pair containing the integer index of the (unordered) band and the corresponding energy level. The index of the pair is the index in the ordered list. The ordering is made in:
bnds = order_bands(a)where "a" is the input diagonal matrix, and bnds is the output - list of 2-element lists.
Populate bands Essentially, once we know the ordered bands, we can also compute the occupation numbers of all levels. This is done using the "populate_bands" function:
occ = populate_bands(Nel, degen, kT, etol, pop_opt, bnds)As you can see, the function takes the following argument: Nel, degen, kT, etol - the smae parameters as we used in the "fermi_energy" function this is because, in case we choose Fermi populations, the Fermi energy will be computed inside the "populate_bands" function; bnds - is the main element containing input data - the ordered bands (indexed energy levels); finally, the pop_opt - is the parameter selecting the option for population scheme: 0 - integer population; 1 - Fermi-populations. The output of the function, "occ", is the list of 2-element lists - the structure fully congruent to the "bnds". The only difference is that "occ" contains populations numbers instead of the level energies.
Density matrix from the effective Hamiltonian (Fock) matrix As we go further in this tutorial, we start getting more and more advanced and useful functionality. This example shows the function "Fock_to_P", which takes as input: the Fock matrix (or any other effective Hamiltonian) "F", the matrix of basis states overlaps "S", the familiar Nel, degen, kT, etol, and pop_opt parameters controlling the computations of Fermi energy and occupations of the states. For the purposes of this example, we have taken F diagonal with the energies as in the Fermi_energy example, the matrix S is chosen to be identity. The function is then called as:
res = Fock_to_P(F,S, Nel, degen, kT, etol, pop_opt)The function Fock_to_P essentially hides several steps we have covered above, plus it imlements some additional functionality. Essentially, the computation algorithm hidden in that function is this:
An important note about Libra design: According to the design philosophy of the Libra package, the code exports functionality of different levels of complexity. This means the user in the Python-script side should, in principle, be able to implement any of the complex functions using only the lower-level utilities. This idea is well illustrated on the above example of "Fock_to_P" funciton. As we have shown, each step of the calculations can be done within Python with the exported simpler functions that have been showcased above. Moreover, once the Python-side algorithm is implemented, it is very straightforward to implement the complex function withing C++ environment - using practically same (only with C++ signatures) simpler functions. Then the complex function, like "Fock_to_P" can be exported to Python to be used as a "simple" function withing an even more complex one. This design makes Libra not just a black box for some calculations, but a "methodology discovery" tool - allowing one to tinkering with all its diverse functionality.
Excitations as objects. Well, when we are talking about excitaions in quantum mechanics, we often imagine a process of promoting one of few electrons from one of several orbitals to other ones. In principle, calculations involving electronic excitations can be made very non-verbose and very implicit, meaning no objects would be needed. However, it is still convenient to "materialize" the excitation (along the lines of quasiparticles called excitons) and create them in some way. Specifically, we want to create an object that would represent such excitation (exciton).
To represent an exciton, let us first ask ourselves - what type of information defines the excitation. The answer is almost obvious - the occupation numbers of the given set of orbitals. In the above examples, we already created an variable containing such information - the "occ" list of 2-element lists. We now want to generalize this type of variable to also describe excited states. Apparently, the structure will be the same, only the occupation numbers of different orbitals will be different. For instance, if the structure [[0,1], [1,1], [2,1], [3,0], [4,0]] represents a ground state with first 3 orbitals occupied (integer occupation scheme) and two other orbitals unoccupied, then the structure [[0,1], [1,1], [2,0], [3,1], [4,0]] would represent one of the excitations: the one in which an electron is promoted from the orbital 2 to orbital 3, the 2->3 excitation. Here we simply demonstrate a function that allows one to conveniently generate such excitations starting from a given ground state configuration:
ex1 = excite(1, 2, res[4]) ex2 = excite(0, 2, res[4])In our example, res[4] is the ground state configuration of the 2-electron 3-state system: res[4] = [[0,1], [1,1], [2,0]]. The function "excite" takes 3 arguments: the indices of source and target orbitals, and the initial occupation numbers list. It generates a new list of occupation numbers, according to chosen excitation. In our example, we create two excitations ex1 = [[0,1], [1,0], [2,1]] (1->2) and ex2 = [[0,0], [1,1], [2,1]] (0->2). Once we have created the excitations, ex1 and ex2, we can use them in some calculations.
Excited state density matrix and energy And the most logical way to use the excitations just created in the previous example, would be to use them to compute excited density matrix and the electronic energy of the excited state. These operations can done with two more functions implemented in the "libcalculators" module - "compute_density_matrix" and "energy_elec", respectively. The computations of the excited density matrices are:
P_ex1 = compute_density_matrix(ex1,res[1]) P_ex2 = compute_density_matrix(ex2,res[1])Here, we use the ground-state MOs (given by the eigenvector matrix C - see the discussion above). This is the second argument to the function, res[1]. The first argument is the occupation matrix (MO-basis density matrix), rpresented by ex1 or ex2. The "compute_density_matrix" returns an object of type MATRIX containing the density matrix in AO basis for given excitation.
Next, the computed density matrices can be used to compute electronic energy of the system:
print "energy(GS) = ", energy_elec(res[2],F,F) print "energy(Ex1) = ", energy_elec(P_ex1,F,F) print "energy(Ex2) = ", energy_elec(P_ex2,F,F)The "energy_elec" function takes 3 arguments - the density matrix, the core Hamiltonian, and the Fock matrix. In our example, the core Hamiltonian is same as the Fock matrix (tight-binding approximation). The only difference in all calculations is that we use different density matrices. For the reference ground state calculations, we use already computed ground state matrix res[2]. For other two excitations we use the excited density matrices P_ex1 and P_ex2. Enventually, since in this tutorial we use the tight-binding approximation with density-independent Fock matrix (e.g. like in extended Huckel model), the excitation energies (differences between the total electronic energies of excited and ground state configurations) are simply the differences of the 1-electron orbital energies for the orbitals involved in the excitation.