Now we are almost done. We miss the implementation for the normal ordering. Here it is: Term_NormalOrder.C
1 Sum<Term, int> Term::normalOrder() const
2 { return normalOrder(false); }
3
4 Sum<Term, int> Term::normalOrder_fullyContractedOnly() const
5 { return normalOrder(true); }
6
7
8 Sum<Term, int> Term::normalOrder(bool fullyContractedOnly) const
9 {
10 Sum<Term, int> sum;
11 for ( unsigned int i=0 ; i+1<_opProd.size() ; ++i ) // iterate over Product<SQOperator> but the last
12 {
13 // check if two consecutive operators need reordering
14 if ( _opProd[i].gender()==SQOperator::Annihilator && _opProd[i+1].gender()==SQOperator::Creator )
15 { // yes
16
17 // handle 1st term
18 Product<SQOperator> p(_opProd); // copy Product<SQOperator>
19 swap(p[i], p[i+1]); // swap operators
20 // check if we need downward recursion
21 if ( !fullyContractedOnly || p[_opProd.size()-1].gender()==SQOperator::Creator )
22 sum -= Term(p, _kProd).normalOrder(fullyContractedOnly);
23
24 // handle 2nd term
25 Product<SQOperator> q(p); // copy Product<SQOperator>
26 q.erase(find(q.begin(), q.end(), p[i])); // remove
27 q.erase(find(q.begin(), q.end(), p[i+1])); // operators
28
29 Product<Kronecker> d(_kProd);
30 d *= Kronecker(p[i].orb(), p[i+1].orb()); // and add kronecker
31 // check if we need downward recursion
32 if ( !fullyContractedOnly || q.size()==0 || q[q.size()-1].gender()==SQOperator::Creator )
33 sum += Term(q, d).normalOrder(fullyContractedOnly);
34 return sum;
35 }
36 }
37
38 sum += *this; // add current term to the intermediate sum
39 return sum;
40 }
Let's test it: NormalOrder_Test.C
1 #include "Term.H"
2
3 using namespace std;
4
5 int main()
6 {
7 Product<SQOperator> p;
8
9 p *= SQOperator('i');
10 p *= SQOperator('j');
11 p *= SQOperator('A');
12 p *= SQOperator('B');
13 // p *= SQOperator('k');
14 // p *= SQOperator('l');
15 // p *= SQOperator('C');
16 // p *= SQOperator('D');
17 // p *= SQOperator('m');
18 // p *= SQOperator('E');
19 // p *= SQOperator('n');
20 // p *= SQOperator('F');
21
22
23 Term term(p);
24
25 cout << term << endl;
26 cout << "--------------------------" << endl;
27 cout << term.normalOrder() << endl;
28 cout << "--------------------------" << endl;
29 cout << term.normalOrder_fullyContractedOnly() << endl;
30
31
32 return 0;
33 }
The output is (download
tarball
into a fresh directory, extract, compile and link with g++ *.C, run a.out)
i*j*A*B -------------------------- A*B*i*j + d(ai)*B*j - d(ai)*d(bj) - d(aj)*B*i + d(aj)*d(bi) - d(bi)*A*j + d(bj)*A*i -------------------------- -d(ai)*d(bj) + d(aj)*d(bi)
It works. Now this is really useful and no longer trivial at all. Play around with it.
Uncomment the additional operators in NormalOrder_Test.C and see what happens.
Switch on compilation with optimization (g++ -O3 *.C) to make it faster.
However, normal ordering by direct application of the anticommutation relations is pretty inefficient (exponential cost). Remember: If your program is too slow, do not blame the language. Search for a better algorithm. In our case this would be an implementation of Wick's theorem or better diagrams/efficient algebraic approaches. Nevertheless, it was not our goal to write the best normal ordering program ever but a tutorial on object orientation. We focused on the application of the object oriented paradigm to a common problem of theoretical physics. Therefore this tutorial could not be exhaustive. We missed makefiles, dependency checks, header guards and much more.