Skip to content

Commit 6ad2c8c

Browse files
committed
Add Project 4 Hints
1 parent ccdf932 commit 6ad2c8c

File tree

4 files changed

+119
-3
lines changed

4 files changed

+119
-3
lines changed

Project#04/README.md

+3-3
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ The most straightforward expression of the AO/MO integral transformation is
2222

2323
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.
2424

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

2828
Notice that none of the *C* coefficients in the above expression have any indices in common. Thus, the summation could be rearranged such that:
@@ -37,15 +37,15 @@ while the AO-basis integrals have eight-fold permutational symmetry, the partial
3737

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

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

4242
## Step #4: Compute the MP2 Energy
4343

4444
<img src="./figures/mp2-energy.png" height="60">
4545

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

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

5050
## Test Cases
5151
The input structures, integrals, etc. for these examples are found in the [input directory](./input).

Project#04/hints/hint1.md

+37
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

+59
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

+20
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)