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

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

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

‎Project#04/project4-instructions.pdf

71.1 KB
Binary file not shown.

‎Project#05/README.md

+17-42
Original file line numberDiff line numberDiff line change
@@ -3,75 +3,50 @@ The coupled cluster model provides a higher level of accuracy beyond the MP2 app
33
## Step #1: Preparing the Spin-Orbital Basis Integrals
44
Each term in the equations given in the above paper by Stanton et al. depends on the T<sub>1</sub> or T<sub>2</sub>
55
amplitudes as well as Fock-matrix elements and antisymmetrized, Dirac-notation two-electron integrals, all given in the molecular spin-orbital basis
6-
(as opposed to the spatial-orbital basis used in the earlier [MP2 Project](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2304).
6+
(as opposed to the spatial-orbital basis used in the earlier [MP2 Project](../Project%2304).
77
Thus, the transformation of the AO-basis integrals into the spatial-MO basis must also include their translation into the spin-orbital basis:
88

9-
```
10-
EQUATION
11-
\begin{eqnarray*}
12-
\langle p q | r s \rangle & \equiv & \int d{\mathbf r}_1 d{\mathbf r}_2 d\omega_1 d\omega_1 \phi_p({\mathbf r}_1)\sigma_p(\omega_1) \phi_q({\mathbf r}_2)\sigma_q(\omega_2) \frac{1}{{\mathbf r}_{12}} \phi_r({\mathbf r}_1)\sigma_r(\omega_1) \phi_s({\mathbf r}_2)\sigma_s(\omega_2) \\
13-
& \equiv & \int d{\mathbf r}_1 d{\mathbf r}_2 \phi_p({\mathbf r}_1) \phi_q({\mathbf r}_2) \frac{1}{{\mathbf r}_{12}} \phi_r({\mathbf r}_1) \phi_s({\mathbf r}_2) \int d\omega_1 d\omega_1 \sigma_p(\omega_1) \sigma_q(\omega_2) \sigma_r(\omega_1) \sigma_s(\omega_2) \\
14-
& \equiv & (p r | q s) \int d\omega_1 d\omega_1 \sigma_p(\omega_1) \sigma_q(\omega_2) \sigma_r(\omega_1) \sigma_s(\omega_2)
15-
\end{eqnarry*}
16-
```
9+
<img src="./figures/spin-orbital-eri.png" height="160">
1710

18-
19-
Thus, if you know the ordering of the orbitals (e.g. all occupied orbitals before virtual orbitals, perhaps alternating between alpha and beta spins), it is straightforward to carry out the integration over the spin components (the sigmas) in the above expression. Thus, each spatial-orbital MO-basis two-electron integral translates to 16 possible spin-orbital integrals, only four of which are non-zero.
11+
Thus, if you know the ordering of the orbitals (e.g. all occupied orbitals before virtual orbitals, perhaps alternating between alpha and beta spins), it is straightforward to carry out the integration over the spin components (the &sigma;'s) in the above expression. Thus, each spatial-orbital MO-basis two-electron integral translates to 16 possible spin-orbital integrals, only four of which are non-zero.
2012

2113
Don't forget that you must also create the spin-orbital Fock matrix:
2214

23-
```
24-
EQUATION
25-
f_{pq} = h_{pq} + \sum_m^{\rm occ} \langle pm || qm \rangle
26-
```
15+
<img src="./figures/spin-orbital-fock.png" height="60">
2716

2817
Suggestion: For simplicity, store the two-electron integrals in a four-dimensional array. This will greatly facilitate debugging of the complicated CCSD equations.
2918

30-
* Hint: Sample spatial- to spin-orbital translation code.
19+
* [Hint 1](./hints/hint1.md): Sample spatial- to spin-orbital translation code.
3120

3221
## Step #2: Build the Initial-Guess Cluster Amplitudes
3322
For Hartree-Fock reference determinants, the most common initial guess for the cluster amplitudes are the Moller-Plesset first-order perturbed wave function:
34-
```
35-
EQUATION
36-
t_i^a = 0
37-
```
38-
39-
```
40-
EQUATION
41-
t_{ij}^{ab} = \frac{\langle ij || ab \rangle}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b}
4223

43-
```
24+
<img src="./figures/init-t-amps.png" height="120">
4425

4526
Note that if you have constructed the Fock matrix, two-electron integrals, and initial-guess amplitudes correctly at this point,
46-
you should be able to compute the MP2 correlation energy using the simple spin-orbital expression and get identical results to those from
47-
[Project #4] (https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2304):
27+
you should be able to compute the MP2 correlation energy using the simple spin-orbital expression and get identical results to those from [Project #4](../Project%2304):
4828

49-
```
50-
EQUATION
51-
E_{\rm MP2} = \frac{1}{4} \sum_{ijab} \langle ij||ab\rangle t_{ij}^{ab}
52-
```
29+
<img src="./figures/mp2-energy.png" height="60">
5330

5431
## Step #3: Calculate the CC Intermediates
55-
Use the spin-orbital Eqs. 3-13 from Stanton's paper to build the two-index (F) and four-index (W) intermediates, as well as the effective doubles (labelled with the Greek letter tau).
32+
Use the spin-orbital Eqs. 3-13 from Stanton's paper to build the two-index (F) and four-index (W) intermediates, as well as the effective doubles (labelled with the Greek letter &tau;).
5633

5734
## Step #4: Compute the Updated Cluster Amplitudes
5835
Use Eqs. 1 and 2 from Stanton's paper to compute the updated T<sub>1</sub> and T<sub>2</sub> cluster amplitudes.
5936

6037
## Step #5: Check for Convergence and Iterate
6138
Calculate the current CC correlation energy:
6239

63-
```
64-
EQUATION
65-
E_{\rm CC} = \sum_{ia} f_{ia} t_i^a + \frac{1}{4} \sum_{ijab} \langle ij||ab\rangle t_{ij}^{ab} + \frac{1}{2} \sum_{ijab} \langle ij||ab\rangle t_i^a t_j^b
66-
```
67-
Compare energies and cluster amplitudes (using RMS differences) between iterations to check for convergence to some specified cutoff.
40+
<img src="./figures/cc-correlation-energy.png" height="60">
41+
42+
Compare energies and cluster amplitudes (using RMS differences) between iterations to check for convergence to some specified cutoff.
6843
If convergence is reached, you're done; if not, return to Step #3 and continue.
6944

7045
## Test Cases
7146
The input structures, integrals, etc. for these examples are found in the
72-
[input directory](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2305/input).
47+
[input directory](./input).
7348

74-
* STO-3G Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2305/output/h2o/STO-3G/output.txt)
75-
* DZ Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2305/output/h2o/DZ/output.txt)
76-
* DZP Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2305/output/h2o/DZP/output.txt)
77-
* STO-3G Methane | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2305/output/ch4/STO-3G/output.txt)
49+
* STO-3G Water | [output](./output/h2o/STO-3G/output.txt)
50+
* DZ Water | [output](./output/h2o/DZ/output.txt)
51+
* DZP Water | [output](./output/h2o/DZP/output.txt)
52+
* STO-3G Methane | [output](./output/ch4/STO-3G/output.txt)

‎Project#05/figures/cc-correlation-energy.png

13.5 KB
Loading

‎Project#05/figures/init-t-amps.png

8.82 KB
Loading

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

8.66 KB
Loading

‎Project#05/figures/spin-orbital-eri.png

36.3 KB
Loading

‎Project#05/figures/spin-orbital-fock.png

7.59 KB
Loading

‎Project#05/hints/hint1.md

+19
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
This block of code takes the one-dimensional array of MO-basis spatial-orbital two-electron integrals we've been using thus far and translates it to a four-dimensional array of antisymmetrized integrals over spin-orbitals. Note that alpha and beta spin functions alternate here, with even-numbered orbitals corresponding to alpha, and odd-numbered to beta. The variable `nmo` is the number of spin orbitals, and the `INDEX` function is the same as that [used before](../../Project%2303/hints/hint7-2.md).
2+
3+
```c++
4+
/* Translate integrals to spin-orbital basis */
5+
for(p=0; p < nmo; p++)
6+
for(q=0; q < nmo; q++)
7+
for(r=0; r < nmo; r++)
8+
for(s=0; s < nmo; s++) {
9+
pr = INDEX(p/2,r/2);
10+
qs = INDEX(q/2,s/2);
11+
prqs = INDEX(pr,qs);
12+
value1 = TEI[prqs] * (p%2 == r%2) * (q%2 == s%2);
13+
ps = INDEX(p/2,s/2);
14+
qr = INDEX(q/2,r/2);
15+
psqr = INDEX(ps,qr);
16+
value2 = TEI[psqr] * (p%2 == s%2) * (q%2 == r%2);
17+
ints[p][q][r][s] = value1 - value2;
18+
}
19+
```

0 commit comments

Comments
 (0)