Full configuration interaction theory
Contents
3. Full configuration interaction theory¶
3.1. Slater determinants as basis states, Repetition¶
The simplest possible choice for many-body wavefunctions are product wavefunctions. That is
because we are really only good at thinking about one particle at a time. Such
product wavefunctions, without correlations, are easy to
work with; for example, if the single-particle states
Similarly, computing matrix elements of operators are relatively easy, because the integrals factorize.
The price we pay is the lack of correlations, which we must build up by using many, many product wavefunctions. (Thus we have a trade-off: compact representation of correlations but difficult integrals versus easy integrals but many states required.)
3.2. Slater determinants as basis states, repetition¶
Because we have fermions, we are required to have antisymmetric wavefunctions, e.g.
etc. This is accomplished formally by using the determinantal formalism
Product wavefunction + antisymmetry = Slater determinant.
3.3. Slater determinants as basis states¶
Properties of the determinant (interchange of any two rows or any two columns yields a change in sign; thus no two rows and no two columns can be the same) lead to the Pauli principle:
No two particles can be at the same place (two columns the same); and
No two particles can be in the same state (two rows the same).
3.4. Slater determinants as basis states¶
As a practical matter, however, Slater determinants beyond
The occupation representation or number representation, using fermion creation and annihilation operators, is compact and efficient. It is also abstract and, at first encounter, not easy to internalize. It is inspired by other operator formalism, such as the ladder operators for the harmonic oscillator or for angular momentum, but unlike those cases, the operators do not have coordinate space representations.
Instead, one can think of fermion creation/annihilation operators as a game of symbols that compactly reproduces what one would do, albeit clumsily, with full coordinate-space Slater determinants.
3.5. Quick repetition of the occupation representation¶
We start with a set of orthonormal single-particle states
To each single-particle state
When acting on the vacuum state
3.6. Quick repetition of the occupation representation¶
But with multiple creation operators we can occupy multiple states:
Now we impose antisymmetry, by having the fermion operators satisfy anticommutation relations:
so that
3.7. Quick repetition of the occupation representation¶
Because of this property, automatically
each index
For some relevant exercises with solutions see chapter 8 of Lecture Notes in Physics, volume 936.
3.8. Full Configuration Interaction Theory¶
We have defined the ansatz for the ground state as
where the index
while a
and a general
3.9. Full Configuration Interaction Theory¶
We can then expand our exact state function for the ground state as
where we have introduced the so-called correlation operator
Since the normalization of
resulting in
3.10. Full Configuration Interaction Theory¶
We rewrite
in a more compact form as
where
and the energy can be written as
3.11. Full Configuration Interaction Theory¶
Normally
is solved by diagonalization setting up the Hamiltonian matrix defined by the basis of all possible Slater determinants. A diagonalization
is equivalent to finding the variational minimum of
where
2 3
< < < ! ! M A T H _ B L O C K
Since the coefficients
3.12. Full Configuration Interaction Theory¶
This leads to
for all sets of
If we then multiply by the corresponding
leading to the identification
3.13. Full Configuration Interaction Theory¶
An alternative way to derive the last equation is to start from
and if this equation is successively projected against all
3.14. Example of a Hamiltonian matrix¶
Suppose, as an example, that we have six fermions below the Fermi level.
This means that we can make at most
$0p-0h$ | $1p-1h$ | $2p-2h$ | $3p-3h$ | $4p-4h$ | $5p-5h$ | $6p-6h$ | |
---|---|---|---|---|---|---|---|
$0p-0h$ | x | x | x | 0 | 0 | 0 | 0 |
$1p-1h$ | x | x | x | x | 0 | 0 | 0 |
$2p-2h$ | x | x | x | x | x | 0 | 0 |
$3p-3h$ | 0 | x | x | x | x | x | 0 |
$4p-4h$ | 0 | 0 | x | x | x | x | x |
$5p-5h$ | 0 | 0 | 0 | x | x | x | x |
$6p-6h$ | 0 | 0 | 0 | 0 | x | x | x |
3.15. Example of a Hamiltonian matrix with a Hartree-Fock basis¶
If we use a Hartree-Fock basis, this corresponds to a particular unitary transformation where matrix elements of the type
$0p-0h$ | $1p-1h$ | $2p-2h$ | $3p-3h$ | $4p-4h$ | $5p-5h$ | $6p-6h$ | |
---|---|---|---|---|---|---|---|
$0p-0h$ | $\tilde{x}$ | 0 | $\tilde{x}$ | 0 | 0 | 0 | 0 |
$1p-1h$ | 0 | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ | 0 | 0 | 0 |
$2p-2h$ | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ | 0 | 0 |
$3p-3h$ | 0 | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ | 0 |
$4p-4h$ | 0 | 0 | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ |
$5p-5h$ | 0 | 0 | 0 | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ |
$6p-6h$ | 0 | 0 | 0 | 0 | $\tilde{x}$ | $\tilde{x}$ | $\tilde{x}$ |
3.16. Shell-model jargon¶
If we do not make any truncations in the possible sets of Slater determinants (many-body states) we can make by distributing
If we make truncations, we have different possibilities
The standard nuclear shell-model. Here we define an effective Hilbert space with respect to a given core. The calculations are normally then performed for all many-body states that can be constructed from the effective Hilbert spaces. This approach requires a properly defined effective Hamiltonian
We can truncate in the number of excitations. For example, we can limit the possible Slater determinants to only
and excitations. This is called a configuration interaction calculation at the level of singles and doubles excitations, or just CISD.We can limit the number of excitations in terms of the excitation energies. If we do not define a core, this defines normally what is called the no-core shell-model approach.
What happens if we have a three-body interaction and a Hartree-Fock basis?
3.17. FCI and the exponential growth¶
Full configuration interaction theory calculations provide in principle, if we can diagonalize numerically, all states of interest. The dimensionality of the problem explodes however quickly.
The total number of Slater determinants which can be built with say
For a model space which comprises the first for major shells only
and multiplying this with the number of proton Slater determinants we end up with approximately with a dimensionality
3.18. Exponential wall¶
This number can be reduced if we look at specific symmetries only. However, the dimensionality explodes quickly!
For Hamiltonian matrices of dimensionalities which are smaller than
, we would use so-called direct methods for diagonalizing the Hamiltonian matrixFor larger dimensionalities iterative eigenvalue solvers like Lanczos’ method are used. The most efficient codes at present can handle matrices of
.
3.19. A non-practical way of solving the eigenvalue problem¶
To see this, we look at the contributions arising from
in Eq. (1), that is we multiply with
If we assume that we have a two-body operator at most, Slater’s rule gives then an equation for the
correlation energy in terms of
or
where the energy
3.20. A non-practical way of solving the eigenvalue problem¶
To see this, we look at the contributions arising from
in Eq. (1), that is we multiply with
3.21. A non-practical way of solving the eigenvalue problem¶
If we assume that we have a two-body operator at most, Slater’s rule gives then an equation for the
correlation energy in terms of
or
where the energy
3.22. Rewriting the FCI equation¶
In our notes on Hartree-Fock calculations,
we have already computed the matrix
3.23. Rewriting the FCI equation¶
Inserting the various matrix elements we can rewrite the previous equation as
This equation determines the correlation energy but not the coefficients
3.24. Rewriting the FCI equation, does not stop here¶
We need more equations. Our next step is to set up
as this equation will allow us to find an expression for the coefficents
3.25. Rewriting the FCI equation, please stop here¶
We see that on the right-hand side we have the energy
3.26. Rewriting the FCI equation, more to add¶
The observant reader will however see that we need an equation for
For
4 4
< < < ! ! M A T H _ B L O C K
and we can isolate the coefficients
3.27. Rewriting the FCI equation, more to add¶
A standard choice for the first iteration is to set
At the end we can rewrite our solution of the Schroedinger equation in terms of
3.28. Summarizing FCI and bringing in approximative methods¶
If we can diagonalize large matrices, FCI is the method of choice since:
It gives all eigenvalues, ground state and excited states
The eigenvectors are obtained directly from the coefficients
which result from the diagonalizationWe can compute easily expectation values of other operators, as well as transition probabilities
Correlations are easy to understand in terms of contributions to a given operator beyond the Hartree-Fock contribution. This is the standard approach in many-body theory.
3.29. Definition of the correlation energy¶
The correlation energy is defined as, with a two-body Hamiltonian,
The coefficients
where the so-called reference energy is the energy we obtain from a Hartree-Fock calculation, that is
3.30. FCI equation and the coefficients¶
However, as we have seen, even for a small case like the four first major shells and a nucleus like oxygen-16, the dimensionality becomes quickly intractable. If we wish to include single-particle states that reflect weakly bound systems, we need a much larger single-particle basis. We need thus approximative methods that sum specific correlations to infinite order.
Popular methods are
Many-body perturbation theory (in essence a Taylor expansion)
Green’s function approaches (matrix inversion)
Similarity group transformation methods (coupled ordinary differential equations)
All these methods start normally with a Hartree-Fock basis as the calculational basis.
3.31. Important ingredients to have in codes¶
Be able to validate and verify the algorithms.
Include concepts like unit testing. Gives the possibility to test and validate several or all parts of the code.
Validation and verification are then included naturally and one can develop a better attitude to what is meant with an ethically sound scientific approach.
3.32. A structured approach to solving problems¶
In the steps that lead to the development of clean code you should think of
How to structure a code in terms of functions (use IDEs or advanced text editors like sublime or atom)
How to make a module
How to read input data flexibly from the command line or files
How to create graphical/web user interfaces
How to write unit tests
How to refactor code in terms of classes (instead of functions only)
How to conduct and automate large-scale numerical experiments
How to write scientific reports in various formats (LaTeX, HTML, doconce)
3.33. Additional benefits¶
Many of the above aspetcs will save you a lot of time when you incrementally extend software over time from simpler to more complicated problems. In particular, you will benefit from many good habits:
New code is added in a modular fashion to a library (modules)
Programs are run through convenient user interfaces
It takes one quick command to let all your code undergo heavy testing
Tedious manual work with running programs is automated,
Your scientific investigations are reproducible, scientific reports with top quality typesetting are produced both for paper and electronic devices. Use version control software like git and repositories like github
3.34. Unit Testing¶
Unit Testing is the practice of testing the smallest testable parts, called units, of an application individually and independently to determine if they behave exactly as expected.
Unit tests (short code fragments) are usually written such that they can be preformed at any time during the development to continually verify the behavior of the code.
In this way, possible bugs will be identified early in the development cycle, making the debugging at later stages much easier.
3.35. Unit Testing, benefits¶
There are many benefits associated with Unit Testing, such as
It increases confidence in changing and maintaining code. Big changes can be made to the code quickly, since the tests will ensure that everything still is working properly.
Since the code needs to be modular to make Unit Testing possible, the code will be easier to reuse. This improves the code design.
Debugging is easier, since when a test fails, only the latest changes need to be debugged.
Different parts of a project can be tested without the need to wait for the other parts to be available.
A unit test can serve as a documentation on the functionality of a unit of the code.
3.36. Simple example of unit test¶
Look up the guide on how to install unit tests for c++ at course webpage. This is the version with classes.
#include <unittest++/UnitTest++.h>
class MyMultiplyClass{
public:
double multiply(double x, double y) {
return x * y;
}
};
TEST(MyMath) {
MyMultiplyClass my;
CHECK_EQUAL(56, my.multiply(7,8));
}
int main()
{
return UnitTest::RunAllTests();
}
3.37. Simple example of unit test¶
And without classes
#include <unittest++/UnitTest++.h>
double multiply(double x, double y) {
return x * y;
}
TEST(MyMath) {
CHECK_EQUAL(56, multiply(7,8));
}
int main()
{
return UnitTest::RunAllTests();
}
For Fortran users, the link at http://sourceforge.net/projects/fortranxunit/ contains a similar software for unit testing. For Python go to https://docs.python.org/2/library/unittest.html.
3.38. Unit tests¶
There are many types of unit test libraries. One which is very popular with C++ programmers is Catch
Catch is header only. All you need to do is drop the file(s) somewhere reachable from your project - either in some central location you can set your header search path to find, or directly into your project tree itself!
This is a particularly good option for other Open-Source projects that want to use Catch for their test suite.
3.39. Examples¶
Computing factorials
inline unsigned int Factorial( unsigned int number ) {
return number > 1 ? Factorial(number-1)*number : 1;
}
3.40. Factorial Example¶
Simple test where we put everything in a single file
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main()
#include "catch.hpp"
inline unsigned int Factorial( unsigned int number ) {
return number > 1 ? Factorial(number-1)*number : 1;
}
TEST_CASE( "Factorials are computed", "[factorial]" ) {
REQUIRE( Factorial(0) == 1 );
REQUIRE( Factorial(1) == 1 );
REQUIRE( Factorial(2) == 2 );
REQUIRE( Factorial(3) == 6 );
REQUIRE( Factorial(10) == 3628800 );
}
This will compile to a complete executable which responds to command line arguments. If you just run it with no arguments it will execute all test cases (in this case there is just one), report any failures, report a summary of how many tests passed and failed and return the number of failed tests.
3.41. What did we do (1)?¶
All we did was
#define
one identifier and
#include
one header and we got everything - even an implementation of main() that will respond to command line arguments. Once you have more than one file with unit tests in you’ll just need to
#include "catch.hpp"
and go. Usually it’s a good idea to have a dedicated implementation file that just has
#define CATCH_CONFIG_MAIN
#include "catch.hpp".
You can also provide your own implementation of main and drive Catch yourself.
3.42. What did we do (2)?¶
We introduce test cases with the
TEST_CASE
macro.
The test name must be unique. You can run sets of tests by specifying a wildcarded test name or a tag expression. All we did was define one identifier and include one header and we got everything.
We write our individual test assertions using the
REQUIRE
macro.
3.43. Unit test summary and testing approach¶
Three levels of tests
Microscopic level: testing small parts of code, use often unit test libraries
Mesoscopic level: testing the integration of various parts of your code
Macroscopic level: testing that the final result is ok
3.44. Coding Recommendations¶
Writing clean and clear code is an art and reflects your understanding of
derivation, verification, and implementation of algorithms
what can go wrong with algorithms
overview of important, known algorithms
how algorithms are used to solve mathematical problems
reproducible science and ethics
algorithmic thinking for gaining deeper insights about scientific problems
Computing is understanding and your understanding is reflected in your abilities to write clear and clean code.
3.45. Summary and recommendations¶
Some simple hints and tips in order to write clean and clear code
Spell out the algorithm and have a top-down approach to the flow of data
Start with coding as close as possible to eventual mathematical expressions
Use meaningful names for variables
Split tasks in simple functions and modules/classes
Functions should return as few as possible variables
Use unit tests and make sure your codes are producing the correct results
Where possible use symbolic coding to autogenerate code and check results
Make a proper timing of your algorithms
Use version control and make your science reproducible
Use IDEs or smart editors with debugging and analysis tools.
Automatize your computations interfacing high-level and compiled languages like C++ and Fortran.
…..
3.46. Building a many-body basis¶
Here we will discuss how we can set up a single-particle basis which we can use in the various parts of our projects, from the simple pairing model to infinite nuclear matter. We will use here the simple pairing model to illustrate in particular how to set up a single-particle basis. We will also use this do discuss standard FCI approaches like:
Standard shell-model basis in one or two major shells
Full CI in a given basis and no truncations
CISD and CISDT approximations
No-core shell model and truncation in excitation energy
3.47. Building a many-body basis¶
An important step in an FCI code is to construct the many-body basis.
While the formalism is independent of the choice of basis, the effectiveness of a calculation will certainly be basis dependent.
Furthermore there are common conventions useful to know.
First, the single-particle basis has angular momentum as a good quantum number. You can imagine the single-particle wavefunctions being generated by a one-body Hamiltonian, for example a harmonic oscillator. Modifications include harmonic oscillator plus spin-orbit splitting, or self-consistent mean-field potentials, or the Woods-Saxon potential which mocks up the self-consistent mean-field. For nuclei, the harmonic oscillator, modified by spin-orbit splitting, provides a useful language for describing single-particle states.
3.48. Building a many-body basis¶
Each single-particle state is labeled by the following quantum numbers:
Orbital angular momentum
Intrinsic spin
= 1/2 for protons and neutronsAngular momentum
-component (or )Some labeling of the radial wavefunction, typically
the number of nodes in the radial wavefunction, but in the case of harmonic oscillator one can also use the principal quantum number , where the harmonic oscillator energy is .
In this format one labels states by
3.49. Building a many-body basis¶
In practice the single-particle space has to be severely truncated. This truncation is typically based upon the single-particle energies, which is the effective energy from a mean-field potential.
Sometimes we freeze the core and only consider a valence space. For example, one
may assume a frozen
Another example is a frozen
Sometimes we refer to nuclei by the valence space where their last nucleons go.
So, for example, we call
3.50. Building a many-body basis¶
There are different kinds of truncations.
For example, one can start with `filled’ orbits (almost always the lowest), and then allow one, two, three… particles excited out of those filled orbits. These are called 1p-1h, 2p-2h, 3p-3h excitations.
Alternately, one can state a maximal orbit and allow all possible configurations with particles occupying states up to that maximum. This is called full configuration.
Finally, for particular use in nuclear physics, there is the energy truncation, also called the
or truncation.
3.51. Building a many-body basis¶
Here one works in a harmonic oscillator basis, with each major oscillator shell assigned a principal quantum number
Excited state are labeled relative to the lowest configuration by the number of harmonic oscillator quanta.
This truncation is useful because if one includes all configuration up to
some
In almost all cases, the many-body Hamiltonian is rotationally invariant. This means
it commutes with the operators
Therefore we can choose to construct a many-body basis which has fixed
Alternately, one can construct a many-body basis which has fixed
3.52. Building a many-body basis¶
The Hamiltonian matrix will have smaller dimensions (a factor of 10 or more) in the
The quantum number
This is not true of
3.53. Building a many-body basis¶
The upshot is that
It is easy to construct a Slater determinant with good total
;It is trivial to calculate
for each Slater determinant;So it is easy to construct an
-scheme basis with fixed total .
Note that the individual
3.54. Building a many-body basis¶
Example: two
Index | $n$ | $l$ | $j$ | $m_j$ |
---|---|---|---|---|
1 | 0 | 0 | 1/2 | -1/2 |
2 | 0 | 0 | 1/2 | 1/2 |
3 | 1 | 0 | 1/2 | -1/2 |
4 | 1 | 0 | 1/2 | 1/2 |
3.55. Building a many-body basis¶
There are
Occupied | $M$ |
---|---|
1,2 | 0 |
1,3 | -1 |
1,4 | 0 |
2,3 | 0 |
2,4 | 1 |
3,4 | 0 |
3.56. Building a many-body basis¶
As another example, consider using only single particle states from the
Index | $n$ | $l$ | $j$ | $m_j$ |
---|---|---|---|---|
1 | 0 | 2 | 5/2 | -5/2 |
2 | 0 | 2 | 5/2 | -3/2 |
3 | 0 | 2 | 5/2 | -1/2 |
4 | 0 | 2 | 5/2 | 1/2 |
5 | 0 | 2 | 5/2 | 3/2 |
6 | 0 | 2 | 5/2 | 5/2 |
3.57. Building a many-body basis¶
There are
Occupied | $M$ | Occupied | $M$ | Occupied | $M$ |
---|---|---|---|---|---|
1,2 | -4 | 2,3 | -2 | 3,5 | 1 |
1,3 | -3 | 2,4 | -1 | 3,6 | 2 |
1,4 | -2 | 2,5 | 0 | 4,5 | 2 |
1,5 | -1 | 2,6 | 1 | 4,6 | 3 |
1,6 | 0 | 3,4 | 0 | 5,6 | 4 |
3.58. Shell-model project¶
The first step is to construct the
The steps could be:
Read in a user-supplied file of single-particle states (examples can be given) or just code these internally;
Ask for the total
of the system and the number of particles ;Construct all the
-particle states with given . You will validate the code by comparing both the number of states and specific states.
3.59. Shell-model project¶
The format of a possible input file could be
Index | $n$ | $l$ | $2j$ | $2m_j$ |
---|---|---|---|---|
1 | 1 | 0 | 1 | -1 |
2 | 1 | 0 | 1 | 1 |
3 | 0 | 2 | 3 | -3 |
4 | 0 | 2 | 3 | -1 |
5 | 0 | 2 | 3 | 1 |
6 | 0 | 2 | 3 | 3 |
7 | 0 | 2 | 5 | -5 |
8 | 0 | 2 | 5 | -3 |
9 | 0 | 2 | 5 | -1 |
10 | 0 | 2 | 5 | 1 |
11 | 0 | 2 | 5 | 3 |
12 | 0 | 2 | 5 | 5 |
3.60. Shell-model project¶
To read in the single-particle states you need to:
Open the file
Read the number of single-particle states (in the above example, 12); allocate memory; all you need is a single array storing
for each state, labeled by the index.Read in the quantum numbers and store
(and anything else you happen to want).
3.61. Shell-model project¶
The next step is to read in the number of particles
You should probably write an error trap to make sure
3.62. Shell-model project¶
The final step is to generate the set of
Hence we can
store the Slater determinant basis states as
3.63. Shell-model project¶
We can construct an occupation representation of Slater determinants by the odometer
method. Consider
(also written as )
Now increase the last occupancy recursively:
Then start over with
and again increase the rightmost digit
3.64. Shell-model project¶
When we restrict ourselves to an
Alternately, and not too difficult, is to run the odometer routine twice: each time, as
as a Slater determinant is calculated, compute
3.65. Shell-model project¶
Some example solutions: Let’s begin with a simple case, the
Index | $n$ | $l$ | $j$ | $m_j$ |
---|---|---|---|---|
1 | 0 | 2 | 5/2 | -5/2 |
2 | 0 | 2 | 5/2 | -3/2 |
3 | 0 | 2 | 5/2 | -1/2 |
4 | 0 | 2 | 5/2 | 1/2 |
5 | 0 | 2 | 5/2 | 3/2 |
6 | 0 | 2 | 5/2 | 5/2 |
, , , , , , , , , , , , , , , , , , ,
Of these, there are only 3 states with
3.66. Shell-model project¶
You should try by hand to show that in this same single-particle space, that for
To test your code, confirm the above.
Also,
for the
3.67. Shell-model project¶
For our project, we will only consider the pairing model.
A simple space is the
Index | $n$ | $l$ | $s$ | $m_s$ |
---|---|---|---|---|
1 | 0 | 0 | 1/2 | -1/2 |
2 | 0 | 0 | 1/2 | 1/2 |
3 | 1 | 0 | 1/2 | -1/2 |
4 | 1 | 0 | 1/2 | 1/2 |
3.68. Shell-model project¶
Another, slightly more challenging space is the
Index | $n$ | $l$ | $s$ | $m_s$ |
---|---|---|---|---|
1 | 0 | 0 | 1/2 | -1/2 |
2 | 0 | 0 | 1/2 | 1/2 |
3 | 1 | 0 | 1/2 | -1/2 |
4 | 1 | 0 | 1/2 | 1/2 |
5 | 2 | 0 | 1/2 | -1/2 |
6 | 2 | 0 | 1/2 | 1/2 |
7 | 3 | 0 | 1/2 | -1/2 |
8 | 3 | 0 | 1/2 | 1/2 |
3.69. Shell-model project¶
In the shell-model context we can interpret this as 4
3.70. Shell-model project¶
For application in the pairing model we can go further and consider only states with
no “broken pairs,” that is, if
3.71. Shell-model project¶
Hints for coding.
Write small modules (routines/functions) ; avoid big functions that do everything. (But not too small.)
Use Unit tests! Write lots of error traps, even for things that are `obvious.’
Document as you go along. The Unit tests serve as documentation. For each function write a header that includes:
a. Main purpose of function and/or unit test
b. names and brief explanation of input variables, if any
c. names and brief explanation of output variables, if any
d. functions called by this function
e. called by which functions
3.72. Shell-model project¶
Hints for coding
Unit tests will save time. Use also IDEs for debugging. If you insist on brute force debugging, print out intermediate values. It’s almost impossible to debug a code by looking at it - the code will almost always win a `staring contest.’
Validate code with SIMPLE CASES. Validate early and often. Unit tests!!
The number one mistake is using a too complex a system to test. For example , if you are computing particles in a potential in a box, try removing the potential - you should get particles in a box. And start with one particle, then two, then three… Don’t start with eight particles.
3.73. Shell-model project¶
Our recommended occupation representation, e.g.
In state-of-the-art shell-model codes, one generally uses bit
representation, i.e.
This is much more compact, but more intricate to code with considerable more overhead. There exist bit-manipulation functions. We will discuss these in more detail at the beginning of the third week.
3.74. Example case: pairing Hamiltonian¶
We consider a space with
The Hamiltonian we consider is
where
and
This problem can be solved using what is called the quasi-spin formalism to obtain the exact results. Thereafter we will try again using the explicit Slater determinant formalism.
3.75. Example case: pairing Hamiltonian¶
One can show (and this is part of the project) that
Now define
Finally you can show
This means the operators
So we rewrite the Hamiltonian to make this explicit:
3.76. Example case: pairing Hamiltonian¶
Because of the SU(2) algebra, we know that the eigenvalues of
But because
We deduce the maximal
Following Racah we introduce the notation
This also works for
3.77. Example case: pairing Hamiltonian¶
Let’s take a specific example:
Now let’s work this out explicitly. The single particle degrees of freedom are defined as
Index | $k$ | $m$ |
---|---|---|
1 | 1 | -1/2 |
2 | -1 | 1/2 |
3 | 2 | -1/2 |
4 | -2 | 1/2 |
5 | 3 | -1/2 |
6 | -3 | 1/2 |
3.78. Example case: pairing Hamiltonian¶
In this basis, the operator
From this we can determine that
so those states all have eigenvalue 0.
3.79. Example case: pairing Hamiltonian¶
Now for further example,
so
The second term vanishes because state 3 is occupied twice, and reordering the last term we get
without picking up a phase.
3.80. Example case: pairing Hamiltonian¶
Continuing in this fashion, with the previous ordering of the many-body states
(
This is useful for our project. One can by hand confirm
that there are 3 eigenvalues
3.81. Example case: pairing Hamiltonian¶
Another example
Using the
Index | $n$ | $l$ | $s$ | $m_s$ |
---|---|---|---|---|
1 | 0 | 0 | 1/2 | -1/2 |
2 | 0 | 0 | 1/2 | 1/2 |
3 | 1 | 0 | 1/2 | -1/2 |
4 | 1 | 0 | 1/2 | 1/2 |
5 | 2 | 0 | 1/2 | -1/2 |
6 | 2 | 0 | 1/2 | 1/2 |
7 | 3 | 0 | 1/2 | -1/2 |
8 | 3 | 0 | 1/2 | 1/2 |
3.82. Example case: pairing Hamiltonian¶
Now we take the following Hamiltonian
where
and
We can write down the
(You should check by hand that this is correct.)
For
3.83. Building a Hamiltonian matrix¶
The goal is to compute the matrix elements of the Hamiltonian, specifically matrix elements between many-body states (Slater determinants) of two-body operators
In particular we will need to compute
where
3.84. Building a Hamiltonian matrix¶
Note: there are other, more efficient ways to do this than the method we describe, but you will be able to produce a working code quickly.
As we coded in the first step,
a Slater determinant
Furthermore, for the two-body matrix elements
What follows here is a more general, but still brute force, approach.
3.85. Building a Hamiltonian matrix¶
Write a function that:
Has as input the single-particle indices
for the two-body operator and the index for the ket Slater determinant;Returns the index
of the unique (if any) Slater determinant such that
as well as the phase
This is equivalent to computing
3.86. Building a Hamiltonian matrix, first step¶
The first step can take as input an initial Slater determinant
(whose position in the list of basis Slater determinants is
It will return another final Slater determinant if the single-particle states
If
3.87. Building a Hamiltonian matrix, second step¶
The second step will take the final Slater determinant
from the first step (if not empty),
and then order by pairwise permutations (i.e., if the Slater determinant is
3.88. Building a Hamiltonian matrix¶
It will also output a phase. If any two single-particle occupancies are repeated,
the phase is
0. Otherwise it is +1 for an even permutation and -1 for an odd permutation to
bring the final
Slater determinant into ascending order,
3.89. Building a Hamiltonian matrix¶
Example: Suppose in the
3.90. Building a Hamiltonian matrix¶
Example: Suppose in the
Lastly, the final step takes the ordered final Slater determinant and
we search through the basis list to
determine its index in the many-body basis, that is,
3.91. Building a Hamiltonian matrix¶
The Hamiltonian is then stored as an
3.92. Building a Hamiltonian matrix¶
Initialize
Set up an outer loop over
Loop over
For each
, loop over and fetch and the single-particle indicesIf
skip. Otherwise, apply to the Slater determinant labeled by .Find, if any, the label
of the resulting Slater determinant and the phase (which is 0, +1, -1).If phase
, then update as . The sum is important because multiple operators might contribute to the same matrix element.Continue loop over
Continue loop over
.End the outer loop over
.
You should force the resulting matrix
3.93. Building a Hamiltonian matrix¶
You will also need to include the single-particle energies. This is easy: they only
contribute to diagonal matrix elements, that is,
Simply find the occupied single-particle states
3.94. Hamiltonian matrix without the bit representation¶
Consider the many-body state
Using the Slater-Condon rules the matrix elements of any one-body
(
3.95. Hamiltonian matrix without the bit representation, one and two-body operators¶
For two determinants which differ only by the substitution of single-particle states
For two determinants which differ by two single-particle states
All other matrix elements involving determinants with more than two substitutions are zero.
3.96. Strategies for setting up an algorithm¶
An efficient implementation of these rules requires
to find the number of single-particle state substitutions between two determinants
to find which single-particle states are involved in the substitution
to compute the phase factor if a reordering of the single-particle states has occured
We can solve this problem using our odometric approach or alternatively using a bit representation as discussed below and in more detail in
We recommend in particular the article by Simen Kvaal. It contains nice general classes for creation and annihilation operators as well as the calculation of the phase (see below).
3.97. Computing expectation values and transitions in the shell-model¶
When we diagonalize the Hamiltonian matrix, the eigenvectors are the coefficients
With these eigenvectors we can compute say the transition likelyhood of a one-body operator as
Writing the one-body operator in second quantization as
we have
3.98. Computing expectation values and transitions in the shell-model and spectroscopic factors¶
The terms we need to evalute then are just the elements
which can be rewritten in terms of spectroscopic factors by inserting a complete set of Slater determinats as
where
3.99. Operators in second quantization¶
In the build-up of a shell-model or FCI code that is meant to tackle large dimensionalities
we need to deal with the action of the Hamiltonian
The time consuming part stems from the action of the Hamiltonian on the above determinant,
A practically useful way to implement this action is to encode a Slater determinant as a bit pattern.
3.100. Operators in second quantization¶
Assume that we have at our disposal
A Slater determinant can then be coded as an integer of
The unoccupied single-particle states have bit value
which translates into the decimal number
We can thus encode a Slater determinant as a bit pattern.
3.101. Operators in second quantization¶
With
The total number of bit patterns is
3.102. Operators in second quantization¶
We assume again that we have at our disposal
in a more compact way as
The action of a creation operator is thus
which becomes
3.103. Operators in second quantization¶
Similarly
which becomes
This gives a simple recipe:
If one of the bits
is and we act with a creation operator on this bit, we return a null vectorIf
, we set it to and return a sign factor , where is the number of bits set before bit .
3.104. Operators in second quantization¶
Consider the action of
What is the simplest way to obtain the phase when we act with one annihilation(creation) operator on the given Slater determinant representation?
3.105. Operators in second quantization¶
We have an SD representation
in a more compact way as
The action of
which becomes
3.106. Operators in second quantization¶
The action
can be obtained by subtracting the logical sum (AND operation) of
from
This operation gives
Similarly, we can form
3.107. Operators in second quantization¶
It is trickier however to get the phase
Let
be a word that represents the 1-bit to be removed and all others set to zero.
In the previous example
Define
as the similar word that represents the bit to be added, that is in our case
Compute then
, which here becomes
Perform then the logical AND operation of
with the word containing
which results in
3.108. Bit counting¶
We include here a python program which may aid in this direction. It uses bit manipulation functions from http://wiki.python.org/moin/BitManipulation.
import math
"""
A simple Python class for Slater determinant manipulation
Bit-manipulation stolen from:
http://wiki.python.org/moin/BitManipulation
"""
# bitCount() counts the number of bits set (not an optimal function)
def bitCount(int_type):
""" Count bits set in integer """
count = 0
while(int_type):
int_type &= int_type - 1
count += 1
return(count)
# testBit() returns a nonzero result, 2**offset, if the bit at 'offset' is one.
def testBit(int_type, offset):
mask = 1 << offset
return(int_type & mask) >> offset
# setBit() returns an integer with the bit at 'offset' set to 1.
def setBit(int_type, offset):
mask = 1 << offset
return(int_type | mask)
# clearBit() returns an integer with the bit at 'offset' cleared.
def clearBit(int_type, offset):
mask = ~(1 << offset)
return(int_type & mask)
# toggleBit() returns an integer with the bit at 'offset' inverted, 0 -> 1 and 1 -> 0.
def toggleBit(int_type, offset):
mask = 1 << offset
return(int_type ^ mask)
# binary string made from number
def bin0(s):
return str(s) if s<=1 else bin0(s>>1) + str(s&1)
def bin(s, L = 0):
ss = bin0(s)
if L > 0:
return '0'*(L-len(ss)) + ss
else:
return ss
class Slater:
""" Class for Slater determinants """
def __init__(self):
self.word = int(0)
def create(self, j):
print "c^+_" + str(j) + " |" + bin(self.word) + "> = ",
# Assume bit j is set, then we return zero.
s = 0
# Check if bit j is set.
isset = testBit(self.word, j)
if isset == 0:
bits = bitCount(self.word & ((1<<j)-1))
s = pow(-1, bits)
self.word = setBit(self.word, j)
print str(s) + " x |" + bin(self.word) + ">"
return s
def annihilate(self, j):
print "c_" + str(j) + " |" + bin(self.word) + "> = ",
# Assume bit j is not set, then we return zero.
s = 0
# Check if bit j is set.
isset = testBit(self.word, j)
if isset == 1:
bits = bitCount(self.word & ((1<<j)-1))
s = pow(-1, bits)
self.word = clearBit(self.word, j)
print str(s) + " x |" + bin(self.word) + ">"
return s
# Do some testing:
phi = Slater()
phi.create(0)
phi.create(1)
phi.create(2)
phi.create(3)
print
s = phi.annihilate(2)
s = phi.create(7)
s = phi.annihilate(0)
s = phi.create(200)
Input In [1]
print "c^+_" + str(j) + " |" + bin(self.word) + "> = ",
^
SyntaxError: invalid syntax
3.109. Eigenvalue problems, basic definitions¶
Let us consider the matrix
where
The above right eigenvector problem is equivalent to a set of
3.110. Eigenvalue problems, basic definitions¶
The eigenvalue problem can be rewritten as
with
which in turn means that the determinant is a polynomial
of degree
3.111. Eigenvalue problems, basic definitions¶
The eigenvalues of a matrix
or
The set of these roots is called the spectrum and is denoted as
and if we define the trace of
then
3.112. Abel-Ruffini Impossibility Theorem¶
The Abel-Ruffini theorem (also known as Abel’s impossibility theorem) states that there is no general solution in radicals to polynomial equations of degree five or higher.
The content of this theorem is frequently misunderstood. It does not assert that higher-degree polynomial equations are unsolvable. In fact, if the polynomial has real or complex coefficients, and we allow complex solutions, then every polynomial equation has solutions; this is the fundamental theorem of algebra. Although these solutions cannot always be computed exactly with radicals, they can be computed to any desired degree of accuracy using numerical methods such as the Newton-Raphson method or Laguerre method, and in this way they are no different from solutions to polynomial equations of the second, third, or fourth degrees.
The theorem only concerns the form that such a solution must take. The content of the theorem is
that the solution of a higher-degree equation cannot in all cases be expressed in terms of the polynomial coefficients with a finite number of operations of addition, subtraction, multiplication, division and root extraction. Some polynomials of arbitrary degree, of which the simplest nontrivial example is the monomial equation
3.113. Abel-Ruffini Impossibility Theorem¶
The Abel-Ruffini theorem says that there are some fifth-degree equations whose solution cannot be so expressed.
The equation
Today, in the modern algebraic context, we say that second, third and fourth degree polynomial
equations can always be solved by radicals because the symmetric groups
3.114. Eigenvalue problems, basic definitions¶
In the present discussion we assume that our matrix is real and symmetric, that is
If
and for
3.115. Eigenvalue problems, basic definitions¶
To obtain the eigenvalues of
We say that a matrix
The importance of a similarity transformation lies in the fact that the resulting matrix has the same eigenvalues, but the eigenvectors are in general different.
3.116. Eigenvalue problems, basic definitions¶
To prove this we
start with the eigenvalue problem and a similarity transformed matrix
We multiply the first equation on the left by
which is the same as
The variable
3.117. Eigenvalue problems, basic definitions¶
The basic philosophy is to
Either apply subsequent similarity transformations (direct method) so that
Or apply subsequent similarity transformations so that
becomes tridiagonal (Householder) or upper/lower triangular (the QR method to be discussed later).Thereafter, techniques for obtaining eigenvalues from tridiagonal matrices can be used.
Or use so-called power methods
Or use iterative methods (Krylov, Lanczos, Arnoldi). These methods are popular for huge matrix problems.
3.118. Discussion of methods for eigenvalues¶
The general overview.
One speaks normally of two main approaches to solving the eigenvalue problem.
The first is the formal method, involving determinants and the characteristic polynomial. This proves how many eigenvalues there are, and is the way most of you learned about how to solve the eigenvalue problem, but for matrices of dimensions greater than 2 or 3, it is rather impractical.
The other general approach is to use similarity or unitary tranformations to reduce a matrix to diagonal form. This is normally done in two steps: first reduce to for example a tridiagonal form, and then to diagonal form. The main algorithms we will discuss in detail, Jacobi’s and Householder’s (so-called direct method) and Lanczos algorithms (an iterative method), follow this methodology.
3.119. Eigenvalues methods¶
Direct or non-iterative methods require for matrices of dimensionality
Year | $n$ | |
---|---|---|
1950 | $n=20$ | (Wilkinson) |
1965 | $n=200$ | (Forsythe et al.) |
1980 | $n=2000$ | Linpack |
1995 | $n=20000$ | Lapack |
This decade | $n\sim 10^5$ | Lapack |
Sloppily speaking, when
3.120. Discussion of methods for eigenvalues¶
If the matrix to diagonalize is large and sparse, direct methods simply become impractical,
also because
many of the direct methods tend to destroy sparsity. As a result large dense matrices may arise during the diagonalization procedure. The idea behind iterative methods is to project the
Matrix | $\mathbf{A}\mathbf{x}=\mathbf{b}$ | $\mathbf{A}\mathbf{x}=\lambda\mathbf{x}$ |
---|---|---|
$\mathbf{A}=\mathbf{A}^*$ | Conjugate gradient | Lanczos |
$\mathbf{A}\ne \mathbf{A}^*$ | GMRES etc | Arnoldi |
3.121. Eigenvalues and Lanczos’ method¶
Basic features with a real symmetric matrix (and normally huge
Lanczos’ algorithm generates a sequence of real tridiagonal matrices
of dimension with , with the property that the extremal eigenvalues of are progressively better estimates of ’ extremal eigenvalues.* The method converges to the extremal eigenvalues.The similarity transformation is
with the first vector
We are going to solve iteratively
with the first vector
3.123. Eigenvalues and Lanczos’ method, tridiagonal and orthogonal matrices¶
Using the fact that
we can rewrite
as
3.125. Eigenvalues and Lanczos’ method, defining the Lanczos’ vectors¶
We have thus
with
and these vectors are called Lanczos vectors.
3.126. Eigenvalues and Lanczos’ method, basic steps¶
We have thus
with
If
is non-zero, then
with