CCD for the electron gas

Log, Autumn 2014, Audun Skau Hansen

In this project we performed some simple Coupled-Cluster Doubles (CCD) calculations for the electron gas. We used the basis and interactions described in Gustav Baardsens PhD-thesis, while the main part of the code was written by Ole Nordli as a part of his master thesis.

The first objective of the project was to reproduce some results from the literature, while the overall goal was to gradually extend the scope to the heterogenous neutron gas with interactions.

Thursday, september 18th.

We wrote a class in C++ which sets up the basis and Fock-matrix for the electron gas. This class was successfully integrated into the prewritten code from Ole, and produced some nice results both for the Hartree-Fock energy and the CCD-calculations very quickly. (The results was in the same order of magnitude as those found in the literature) As we still lacked some of the insight into the theoretical foundation, we was not able to initialize the basis with the exact same parameters as found in the literature.

Pending work: * To alter the parameter $L$ ($L^3$ is the total volume of the box) so that it corresponds to the given $r_s$ from the literature.

To look into how to implement periodic boundary conditions.*

Monday, October 6th.

I have written a coupled cluster extension (class) for the restricted Hartree-Fock-Roothaan solver that I and Gøran created for the project in FYS4411. Despite extensive debugging, benchmarking and comparisons, I am unable to reproduce results from the literature.

I have beeen trying to reproduce the results from a RHF-CCSD calculation by NWChem, as described in the following document:

http://institute.loni.org/NWChem2012/documents/tce-session.pdf.

I am supposed to find

CCSD iterations

Iter Residuum Correlation Cpu Wall
1 0.0891232379551 -0.0358672469179 0.4 0.4
2 0.0317596201320 -0.0454068882657 0.4 0.4
3 0.0126828916023 -0.0483870059027 0.4 0.4
4 0.0053832778844 -0.0494370597647 0.4 0.4
5 0.0023954452285 -0.0498391184890 0.4 0.4
6 0.0011108272683 -0.0500021724029 0.4 0.4
7 0.0005330437725 -0.0500711904756 0.4 0.4
8 0.0002625570400 -0.0501014381364 0.4 0.4
9 0.0001317707641 -0.0501150974135 0.4 0.4
10 0.0000669953658 -0.0501214303300 0.4 0.4
11 0.0000343654836 -0.0501244348663 0.4 0.4
12 0.0000177357555 -0.0501258887096 0.4 0.4
13 0.0000091923943 -0.0501266039080 0.4 0.4
14 0.0000047789221 -0.0501269605251 0.4 0.4
15 0.0000024899940 -0.0501271402835 0.4 0.4
16 0.0000012995442 -0.0501272316751 0.4 0.4
17 0.0000006791083 -0.0501272784536 0.4 0.4
18 0.0000003552400 -0.0501273025229 0.4 0.4
19 0.0000001859739 -0.0501273149581 0.4 0.4
20 0.0000000974237 -0.0501273214031 0.4 0.4

My own RHS-CCD calculation (the result is comparable for CCSD) finds instead a value of -0.0875531 a.u., overestimating the energy by a factor of approximately 1.7. The table below shows the iterations and lists the separate contributions to the energy.

Iteration CCDEnergy 0th Cont. Lin.Cont. Quad.cont.
1 -0.0358976 0 0 0
2 -0.0562349 -0.0358976 -0.0203373 0
3 -0.0682 -0.0358976 -0.0325562 0.000253831
4 -0.0753594 -0.0358976 -0.040085 0.000623145
5 -0.0797398 -0.0358976 -0.044754 0.000911875
6 -0.082475 -0.0358976 -0.0476843 0.00110689
7 -0.0842106 -0.0358976 -0.0495466 0.00123363
8 -0.0853254 -0.0358976 -0.0507435 0.00131574
9 -0.0860489 -0.0358976 -0.0515206 0.00136931
10 -0.0865231 -0.0358976 -0.0520301 0.00140458
11 -0.0868371 -0.0358976 -0.0523675 0.00142802
12 -0.0870474 -0.0358976 -0.0525935 0.00144376
13 -0.0871899 -0.0358976 -0.0527467 0.00145447
14 -0.0872878 -0.0358976 -0.052852 0.00146186
15 -0.087356 -0.0358976 -0.0529254 0.00146703
16 -0.0874043 -0.0358976 -0.0529774 0.00147072
17 -0.0874391 -0.0358976 -0.0530149 0.00147339
18 -0.0874645 -0.0358976 -0.0530423 0.00147536
19 -0.0874834 -0.0358976 -0.0530627 0.00147684
20 -0.0874977 -0.0358976 -0.053078 0.00147797
21 -0.0875086 -0.0358976 -0.0530898 0.00147885
22 -0.0875171 -0.0358976 -0.053099 0.00147954
23 -0.0875237 -0.0358976 -0.0531062 0.00148009
24 -0.087529 -0.0358976 -0.0531119 0.00148053
25 -0.0875333 -0.0358976 -0.0531165 0.00148089
26 -0.0875367 -0.0358976 -0.0531202 0.00148118
27 -0.0875395 -0.0358976 -0.0531233 0.00148142
28 -0.0875418 -0.0358976 -0.0531258 0.00148162
29 -0.0875436 -0.0358976 -0.0531278 0.00148178
30 -0.0875452 -0.0358976 -0.0531295 0.00148192
31 -0.0875465 -0.0358976 -0.0531309 0.00148204
32 -0.0875475 -0.0358976 -0.0531321 0.00148213
33 -0.0875484 -0.0358976 -0.053133 0.00148221
34 -0.0875492 -0.0358976 -0.0531338 0.00148228
35 -0.0875498 -0.0358976 -0.0531345 0.00148234
36 -0.0875503 -0.0358976 -0.0531351 0.00148239
37 -0.0875508 -0.0358976 -0.0531356 0.00148243
38 -0.0875511 -0.0358976 -0.053136 0.00148246
39 -0.0875515 -0.0358976 -0.0531363 0.00148249
40 -0.0875517 -0.0358976 -0.0531366 0.00148251
41 -0.0875519 -0.0358976 -0.0531368 0.00148253
42 -0.0875521 -0.0358976 -0.0531371 0.00148255
43 -0.0875523 -0.0358976 -0.0531372 0.00148256
44 -0.0875524 -0.0358976 -0.0531374 0.00148258
45 -0.0875525 -0.0358976 -0.0531375 0.00148259
46 -0.0875526 -0.0358976 -0.0531376 0.0014826
47 -0.0875527 -0.0358976 -0.0531377 0.0014826
48 -0.0875528 -0.0358976 -0.0531378 0.00148261
49 -0.0875528 -0.0358976 -0.0531378 0.00148262
50 -0.0875529 -0.0358976 -0.0531379 0.00148262
51 -0.0875529 -0.0358976 -0.0531379 0.00148262
52 -0.087553 -0.0358976 -0.053138 0.00148263
53 -0.087553 -0.0358976 -0.053138 0.00148263
54 -0.087553 -0.0358976 -0.053138 0.00148263
55 -0.087553 -0.0358976 -0.053138 0.00148263
56 -0.0875531 -0.0358976 -0.0531381 0.00148264
57 -0.0875531 -0.0358976 -0.0531381 0.00148264
58 -0.0875531 -0.0358976 -0.0531381 0.00148264
59 -0.0875531 -0.0358976 -0.0531381 0.00148264

As already noted, I am unable to find any errors in the implementation. There is, however, a number of aspects in the method that remains unclear to me. Some of these are:

(1) Shavitt-Bartlett derives the CCD equations, ending in the T2 amplitude equation (9.126) on page 288. In this equation they restrict the summations to $$ i>j, a>b $$ It is unclear to me if this is due to the symmetries in the matrix elementes $\langle ab || ij \rangle$ so that one may manually set the others for each iteration, or if it means that the T2 amplitude tensor will contain a lot of zeros. In my code I have interpreted this in the former manner.

(2) I have manually ensured that my coefficient matrix obtained from the RHF-procedure contains the same elements at the same locations as Ole Tobias Norlis elements. This means that the CC-iterations have a very similar initial condition as the ones described in the literature, as may be seen from the initial energy in the table above. However, upon initialization, Ole sets up his elements like

$\langle pr||qs\rangle = \langle pq|rs\rangle - \langle ps|rq\rangle$

(permuting the r and q index). I adopted this in my code to reproduce his basis, but I did not compensate for this in the CC iterations since I found no reason to do so.

Conclusion

The error in the code needs to be identified and fixed. As the initial energy is very comparable to the ones found in the literature, the error is most likely related to the iterative process itself.

Wednesday, October 8th

I have made the output from the code more specific as to how much each term in the CCD T2 amplitudes contributes. Following Shavitt-Bartlett, this means breaking the full correlation energy down to the following parts:

  • MBPT(2)
  • Linear contributions (functions of t)
    • L2a
    • L2b
    • L2c
  • Quadratic contributions (functions of tt)
    • Qa
    • Qb
    • Qc
    • Qd

My aim for today is to compare each constituent contribution to similar results from Ole Tobias Norlis. I have prepared calculations for both H2O and H2 using a STO-3g basis through a RHF procedure prior to the CCD calculation.

Breakthrough, 03.28 pm.

With help from Ole I finally identified the error in the code. It was actually a bug occuring in the for-loops for the linear L2b contribution.

I now have a working CCD and possibly a CCSD (needs more testing) implementation for my RHF solver.

Thursday, October 9th

With the new progress it is possible to perform benchmark calculations on certain molecules and atoms. I therefore plan to perform the following experiments today (if time permits it):

  • H2O CCSD
  • H2 CCSD
  • Electron gas CCSD (CCD)

I observed some minor discrepancies in the H2O CCSD calculation performed yesterday, so I will probably spend some time reading through the singles contribution code today. I suspect there is a very minor indexing error in one contribution.

Conclusion, 04.21 pm.

I performed some of the tests mentioned above. The results can be found in the notebook "Benchmarking.ipynb". The results was satisfactory compared to other similar software.

I also inspected all the contributions to the CCSD code and found no errors. I now feel a bit more confident that I have a working CCD and CCSD code. After a meeting with Morten, we agreed on some steps ahead:

  • Implement an UHF extension in HFSolve
  • Begin optimization of the code
  • Do more benchmarking with Oles code

In [ ]: