Skip to content

Commit 9826cc6

Browse files
amjamesamjames
authored
Merge pull request CrawfordGroup#19 from kcpearce/add-eqns
Add Equations to Projects 4-14
2 parents e569e65 + 12b5a00 commit 9826cc6

File tree

85 files changed

+1205
-553
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

85 files changed

+1205
-553
lines changed

‎Project#03/figures/compound-index-restrictions2.png‎

15 Bytes
Loading

‎Project#03/figures/eri.png‎

62 Bytes
Loading

‎Project#04/README.md‎

Lines changed: 20 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -1,77 +1,57 @@
11
# Project #4: The second-order Moller-Plesset perturbation (MP2) energy
22
Unlike the Hartree-Fock energy, correlation energies like the MP2 energy are usually expressed in terms of MO-basis quantities (integrals, MO energies).
3-
The most expensive part of the calculation is the transformation of the two-electron integrals from the AO to the MO basis representation.
4-
The purpose of this project is to understand this transformation and fundamental techniques for its efficient implementation.
5-
The theoretical background and a concise set of instructions for this project may be found [[http://sirius.chem.vt.edu/~crawdad/programming/mp2.pdf|here]].
3+
The most expensive part of the calculation is the transformation of the two-electron integrals from the AO to the MO basis representation.
4+
The purpose of this project is to understand this transformation and fundamental techniques for its efficient implementation.
5+
The theoretical background and a concise set of instructions for this project may be found [here](./project4-instructions.pdf).
66

77
## Step #1: Read the Two-Electron Integrals
88
The Mulliken-ordered integrals are defined as:
99

10-
```
11-
EQUATION
12-
(\mu \nu|\lambda \sigma) \equiv \int \phi_\mu({\mathbf r_1})
13-
\phi_\nu({\mathbf r_1}) r_{12}^{-1} \phi_\lambda({\mathbf r_2})
14-
\phi_\sigma({\mathbf r_2}) d{\mathbf r_1} d{\mathbf r_2}
15-
```
10+
<img src="./figures/eri.png" height="45">
1611

17-
Concise instructions for this step can be found in
18-
[Project #3](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2303)
12+
Concise instructions for this step can be found in [Project #3](../Project%2303).
1913

2014
## Step #2: Obtain the MO Coefficients and MO Energies
21-
Use the values you computed in the Hartree-Fock program of
22-
[Project #3](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2303)
15+
Use the values you computed in the Hartree-Fock program of [Project #3](../Project%2303).
2316
## Step #3: Transform the Two-Electron Integrals to the MO Basis
2417
### The Noddy Algorithm
2518

2619
The most straightforward expression of the AO/MO integral transformation is
2720

28-
```
29-
EQUATION
30-
(pq|rs) = \sum_\mu \sum_\nu \sum_\lambda \sum_\sigma C_\mu^p C_\nu^q
31-
(\mu \nu|\lambda \sigma) C_\lambda^r C_\sigma^s
32-
```
21+
<img src="./figures/noddy-transform.png" height="50">
3322

34-
This approach is easy to implement (hence the word "[noddy](http://www.hackerslang.com/noddy.html)" above), but is expensive due to its N<sup>8</sup> computational order.
35-
Nevertheless, you should start with this algorithm to get your code working, and run timings (use the UNIX "time" command)
36-
for the test cases below to get an idea of its computational cost.
23+
This approach is easy to implement (hence the word "[noddy](http://www.hackerslang.com/noddy.html)" above), but is expensive due to its N<sup>8</sup> computational order. Nevertheless, you should start with this algorithm to get your code working, and run timings (use the UNIX "time" command) for the test cases below to get an idea of its computational cost.
3724

38-
* Hint 1: Noddy transformation code
25+
* [Hint 1](./hints/hint1.md): Noddy transformation code
3926
### The Smarter Algorithm
4027

4128
Notice that none of the *C* coefficients in the above expression have any indices in common. Thus, the summation could be rearranged such that:
42-
```
43-
EQUATION
44-
(pq|rs) = \sum_\mu C_\mu^p \left[ \sum_\nu C_\nu^q \left[ \sum_\lambda C_\lambda^r \left[ \sum_\sigma C_\sigma^s (\mu \nu|\lambda \sigma) \right] \right] \right].
45-
```
29+
30+
<img src="./figures/smart-transform.png" height="60">
4631

4732
This means that each summation within brackets could be carried out separately,
48-
starting from the innermost summation over <html>&#963;</html>, if we store the results at each step.
49-
This reduces the N<sup>8</sup> algorithm above to four N<sup>5</sup> steps.
33+
starting from the innermost summation over <html>&#963;</html>, if we store the results at each step. This reduces the N<sup>8</sup> algorithm above to four N<sup>5</sup> steps.
5034

5135
Be careful about the permutational symmetry of the integrals and intermediates at each step:
5236
while the AO-basis integrals have eight-fold permutational symmetry, the partially transformed integrals have much less symmetry.
5337

5438
After you have the noddy algorithm working and timed, modify it to use this smarter algorithm and time the test cases again. Enjoy!
5539

56-
* Hint 2: Smarter transformation code
40+
* [Hint 2](./hints/hint2.md): Smarter transformation code
5741

5842
## Step #4: Compute the MP2 Energy
5943

60-
```
61-
EQUATION
62-
E_{\rm MP2} = \sum_{ij} \sum_{ab} \frac{(ia|jb) \left[ 2 (ia|jb) -
63-
(ib|ja) \right]}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b},
64-
```
44+
<img src="./figures/mp2-energy.png" height="60">
6545

6646
where *i* and *j* denote doubly-occupied orbitals and *a* and *b* unoccupied orbitals, and the denominator involves the MO energies.
6747

68-
* Hint 3: Occupied versus unoccupied orbitals
48+
* [Hint 3](./hints/hint3.md): Occupied versus unoccupied orbitals
6949

7050
## Test Cases
71-
The input structures, integrals, etc. for these examples are found in the [input directory](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2304/input).
51+
The input structures, integrals, etc. for these examples are found in the [input directory](./input).
7252

73-
* STO-3G Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2304/output/h2o/STO-3G/output.txt)
74-
* DZ Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2304/output/h2o/DZ/output.txt)
75-
* DZP Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2304/output/h2o/DZP/output.txt)
76-
* STO-3G Methane | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2304/output/ch4/STO-3G/output.txt)
53+
* STO-3G Water | [output](./output/h2o/STO-3G/output.txt)
54+
* DZ Water | [output](./output/h2o/DZ/output.txt)
55+
* DZP Water | [output](./output/h2o/DZP/output.txt)
56+
* STO-3G Methane | [output](./output/ch4/STO-3G/output.txt)
7757

‎Project#04/figures/eri.png‎

10.6 KB
Loading

‎Project#04/figures/mp2-energy.png‎

11.8 KB
Loading

‎Project#04/figures/noddy-transform.png‎

10.4 KB
Loading

‎Project#04/figures/smart-transform.png‎

12.9 KB
Loading

‎Project#04/hints/hint1.md‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
Here's a code block that demonstrates how to carry out the AO to MO two-electron integral transformation using a single N<sup>8</sup> step. Note that the original AO- and transformed MO-basis integrals are stored in a one-dimensional array taking full advantage of their eight-fold permutational symmetry, as described in [Project #3](../../Project%2303).
2+
3+
Note how the loop structure in the MO-indices (*i*, *j*, *k*, and *l*) automatically takes into account the permutational symmetry of the fully transformed integrals.
4+
5+
```c++
6+
#define INDEX(i,j) ((i>j) ? (((i)*((i)+1)/2)+(j)) : (((j)*((j)+1)/2)+(i)))
7+
...
8+
double *TEI_AO, *TEI_MO;
9+
int i, j, k, l, ijkl;
10+
int p, q, r, s, pq, rs, pqrs;
11+
...
12+
13+
for(i=0,ijkl=0; i < nao; i++) {
14+
for(j=0; j <= i; j++) {
15+
for(k=0; k <= i; k++) {
16+
for(l=0; l <= (i==k ? j : k); l++,ijkl++) {
17+
18+
for(p=0; p < nao; p++) {
19+
for(q=0; q < nao; q++) {
20+
pq = INDEX(p,q);
21+
for(r=0; r < nao; r++) {
22+
for(s=0; s < nao; s++) {
23+
rs = INDEX(r,s);
24+
pqrs = INDEX(pq,rs);
25+
26+
TEI_MO[ijkl] += C[p][i] * C[q][j] * C[r][k] * C[s][l] * TEI_AO[pqrs];
27+
28+
}
29+
}
30+
}
31+
}
32+
33+
}
34+
}
35+
}
36+
}
37+
```

‎Project#04/hints/hint2.md‎

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
Here's a code block that demonstrates how to carry out the AO to MO two-electron integral transformation using four N<sup>5</sup> steps. Note that the original AO- and transformed MO-basis integrals are stored in a one-dimensional array taking advantage of their eight-fold permutational symmetry, as described in [Project #3](../../Project%2303). However, the half-transformed integrals, which lack bra-ket permutational symmetry, are stored in a two-dimensional array.
2+
3+
A number of convenient functions [`mmult()`, `init_matrix()`, etc.] used in this code block can be found in the [diag.cc](http://sirius.chem.vt.edu/~crawdad/programming/diag.cc) discussed in [Project #1](../../Project%2301).
4+
5+
```c++
6+
#define INDEX(i,j) ((i>j) ? (((i)*((i)+1)/2)+(j)) : (((j)*((j)+1)/2)+(i)))
7+
...
8+
double **X, **Y, **TMP, *TEI;
9+
int i, j, k, l, ij, kl, ijkl, klij;
10+
11+
...
12+
X = init_matrix(nao, nao);
13+
Y = init_matrix(nao, nao);
14+
15+
TMP = init_matrix((nao*(nao+1)/2),(nao*(nao+1)/2));
16+
for(i=0,ij=0; i < nao; i++)
17+
for(j=0; j <= i; j++,ij++) {
18+
for(k=0,kl=0; k < nao; k++)
19+
for(l=0; l <= k; l++,kl++) {
20+
ijkl = INDEX(ij,kl);
21+
X[k][l] = X[l][k] = TEI[ijkl];
22+
}
23+
zero_matrix(Y, nao, nao);
24+
mmult(C, 1, X, 0, Y, nao, nao, nao);
25+
zero_matrix(X, nao, nao);
26+
mmult(Y, 0, C, 0, X, nao, nao, nao);
27+
for(k=0, kl=0; k < nao; k++)
28+
for(l=0; l <= k; l++, kl++)
29+
TMP[kl][ij] = X[k][l];
30+
}
31+
32+
zero_array(TEI, (nao*(nao+1)/2)*((nao*(nao+1)/2)+1)/2);
33+
34+
for(k=0,kl=0; k < nao; k++)
35+
for(l=0; l <= k; l++,kl++) {
36+
zero_matrix(X, nao, nao);
37+
zero_matrix(Y, nao, nao);
38+
for(i=0,ij=0; i < nao; i++)
39+
for(j=0; j <=i; j++,ij++)
40+
X[i][j] = X[j][i] = TMP[kl][ij];
41+
zero_matrix(Y, nao, nao);
42+
mmult(C, 1, X, 0, Y, nao, nao, nao);
43+
zero_matrix(X, nao, nao);
44+
mmult(Y, 0, C, 0, X, nao, nao, nao);
45+
for(i=0, ij=0; i < nao; i++)
46+
for(j=0; j <= i; j++,ij++) {
47+
klij = INDEX(kl,ij);
48+
TEI[klij] = X[i][j];
49+
}
50+
}
51+
52+
...
53+
54+
free_matrix(X, nao);
55+
free_matrix(Y, nao);
56+
free_matrix(TMP, (nao*(nao+1)/2));
57+
58+
```
59+

‎Project#04/hints/hint3.md‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
The following code block illustrates a simple-minded MP2 energy calculation with the two-electron integrals stored in a one-dimensional array. The molecular orbitals are assumed to be ordered with all doubly-occupied orbitals first, followed by the unoccupied/virtual orbitals. Notice the difference in limits of the loops over *i* and *a*:
2+
3+
```c++
4+
Emp2 = 0.0;
5+
for(i=0; i < ndocc; i++) {
6+
for(a=ndocc; a < nao; a++) {
7+
ia = INDEX(i,a);
8+
for(j=0; j < ndocc; j++) {
9+
ja = INDEX(j,a);
10+
for(b=ndocc; b < nao; b++) {
11+
jb = INDEX(j,b);
12+
ib = INDEX(i,b);
13+
iajb = INDEX(ia,jb);
14+
ibja = INDEX(ib,ja);
15+
Emp2 += TEI[iajb] * (2 * TEI[iajb] - TEI[ibja])/(eps[i] + eps[j] - eps[a] - eps[b]);
16+
}
17+
}
18+
}
19+
}
20+
```

0 commit comments

Comments
 (0)