CrawfordGroup · amjames · Apr 5, 2017 · Mar 31, 2017 · Mar 31, 2017 · Mar 31, 2017
diff --git a/Project#03/figures/compound-index-restrictions2.png b/Project#03/figures/compound-index-restrictions2.png
diff --git a/Project#03/figures/eri.png b/Project#03/figures/eri.png
diff --git a/Project#04/README.md b/Project#04/README.md
@@ -1,77 +1,57 @@
 # Project #4: The second-order Moller-Plesset perturbation (MP2) energy
 Unlike the Hartree-Fock energy, correlation energies like the MP2 energy are usually expressed in terms of MO-basis quantities (integrals, MO energies). 
-The most expensive part of the calculation is the transformation of the two-electron integrals from the AO to the MO basis representation.  
-The purpose of this project is to understand this transformation and fundamental techniques for its efficient implementation. 
-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]].
+The most expensive part of the calculation is the transformation of the two-electron integrals from the AO to the MO basis representation. 
+The purpose of this project is to understand this transformation and fundamental techniques for its efficient implementation.
+The theoretical background and a concise set of instructions for this project may be found [here](./project4-instructions.pdf).
 
 ## Step #1: Read the Two-Electron Integrals
 The Mulliken-ordered integrals are defined as:
 
-```
-EQUATION
-(\mu \nu|\lambda \sigma) \equiv \int \phi_\mu({\mathbf r_1})
-\phi_\nu({\mathbf r_1}) r_{12}^{-1} \phi_\lambda({\mathbf r_2})
-\phi_\sigma({\mathbf r_2}) d{\mathbf r_1} d{\mathbf r_2}
-```
+<img src="./figures/eri.png" height="45">
 
-Concise instructions for this step can be found in 
-[Project #3](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2303)
+Concise instructions for this step can be found in [Project #3](../Project%2303).
 
 ## Step #2: Obtain the MO Coefficients and MO Energies 
-Use the values you computed in the Hartree-Fock program of 
-[Project #3](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2303)
+Use the values you computed in the Hartree-Fock program of [Project #3](../Project%2303).
 ## Step #3: Transform the Two-Electron Integrals to the MO Basis 
 ### The Noddy Algorithm 
 
 The most straightforward expression of the AO/MO integral transformation is
 
-```
-EQUATION
-(pq|rs) = \sum_\mu \sum_\nu \sum_\lambda \sum_\sigma C_\mu^p C_\nu^q
-(\mu \nu|\lambda \sigma) C_\lambda^r C_\sigma^s
-```
+<img src="./figures/noddy-transform.png" height="50">
 
-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.
+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.
 
-  * Hint 1: Noddy transformation code
+  * [Hint 1](./hints/hint1.md): Noddy transformation code
 ### The Smarter Algorithm
 
 Notice that none of the *C* coefficients in the above expression have any indices in common.  Thus, the summation could be rearranged such that:
-```
-EQUATION 
-(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].
-```
+
+<img src="./figures/smart-transform.png" height="60">
 
 This means that each summation within brackets could be carried out separately, 
-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.
+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.
 
 Be careful about the permutational symmetry of the integrals and intermediates at each step: 
 while the AO-basis integrals have eight-fold permutational symmetry, the partially transformed integrals have much less symmetry.
 
 After you have the noddy algorithm working and timed, modify it to use this smarter algorithm and time the test cases again.  Enjoy!
 
-  * Hint 2: Smarter transformation code
+  * [Hint 2](./hints/hint2.md): Smarter transformation code
 
 ## Step #4: Compute the MP2 Energy
 
-```
-EQUATION
-E_{\rm MP2} = \sum_{ij} \sum_{ab} \frac{(ia|jb) \left[ 2 (ia|jb) -
-(ib|ja) \right]}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b},
-```
+<img src="./figures/mp2-energy.png" height="60">
 
 where *i* and *j* denote doubly-occupied orbitals and *a* and *b* unoccupied orbitals, and the denominator involves the MO energies.
 
-  * Hint 3: Occupied versus unoccupied orbitals
+  * [Hint 3](./hints/hint3.md): Occupied versus unoccupied orbitals
 
 ## Test Cases
-The input structures, integrals, etc. for these examples are found in the [input directory](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2304/input).
+The input structures, integrals, etc. for these examples are found in the [input directory](./input).
 
-* STO-3G Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2304/output/h2o/STO-3G/output.txt)
-* DZ Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2304/output/h2o/DZ/output.txt)
-* DZP Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2304/output/h2o/DZP/output.txt)
-* STO-3G Methane | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2304/output/ch4/STO-3G/output.txt)
+* STO-3G Water | [output](./output/h2o/STO-3G/output.txt)
+* DZ Water | [output](./output/h2o/DZ/output.txt)
+* DZP Water | [output](./output/h2o/DZP/output.txt)
+* STO-3G Methane | [output](./output/ch4/STO-3G/output.txt)
 
diff --git a/Project#04/figures/eri.png b/Project#04/figures/eri.png
diff --git a/Project#04/figures/mp2-energy.png b/Project#04/figures/mp2-energy.png
diff --git a/Project#04/figures/noddy-transform.png b/Project#04/figures/noddy-transform.png
diff --git a/Project#04/figures/smart-transform.png b/Project#04/figures/smart-transform.png
diff --git a/Project#04/hints/hint1.md b/Project#04/hints/hint1.md
@@ -0,0 +1,37 @@
+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).
+
+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.
+
+```c++
+  #define INDEX(i,j) ((i>j) ? (((i)*((i)+1)/2)+(j)) : (((j)*((j)+1)/2)+(i)))
+  ...
+  double *TEI_AO, *TEI_MO;
+  int i, j, k, l, ijkl;
+  int p, q, r, s, pq, rs, pqrs;
+  ...
+
+  for(i=0,ijkl=0; i < nao; i++) {
+    for(j=0; j <= i; j++) {
+      for(k=0; k <= i; k++) {
+        for(l=0; l <= (i==k ? j : k); l++,ijkl++) {
+
+          for(p=0; p < nao; p++) {
+            for(q=0; q < nao; q++) {
+              pq = INDEX(p,q);
+              for(r=0; r < nao; r++) {
+                for(s=0; s < nao; s++) {
+                  rs = INDEX(r,s);
+                  pqrs = INDEX(pq,rs);
+
+                  TEI_MO[ijkl] += C[p][i] * C[q][j] * C[r][k] * C[s][l] * TEI_AO[pqrs];
+
+                }
+              }
+            }
+          }
+
+        }
+      }
+    }
+  }
+```
diff --git a/Project#04/hints/hint2.md b/Project#04/hints/hint2.md
@@ -0,0 +1,59 @@
+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.
+
+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).
+
+```c++
+  #define INDEX(i,j) ((i>j) ? (((i)*((i)+1)/2)+(j)) : (((j)*((j)+1)/2)+(i)))
+  ...
+  double **X, **Y, **TMP, *TEI;
+  int i, j, k, l, ij, kl, ijkl, klij;
+
+  ...
+  X = init_matrix(nao, nao);
+  Y = init_matrix(nao, nao);
+
+  TMP = init_matrix((nao*(nao+1)/2),(nao*(nao+1)/2));
+  for(i=0,ij=0; i < nao; i++)
+    for(j=0; j <= i; j++,ij++) {
+      for(k=0,kl=0; k < nao; k++)
+        for(l=0; l <= k; l++,kl++) {
+          ijkl = INDEX(ij,kl);
+          X[k][l] = X[l][k] = TEI[ijkl];
+        }
+      zero_matrix(Y, nao, nao);
+      mmult(C, 1, X, 0, Y, nao, nao, nao);
+      zero_matrix(X, nao, nao);
+      mmult(Y, 0, C, 0, X, nao, nao, nao);
+      for(k=0, kl=0; k < nao; k++)
+        for(l=0; l <= k; l++, kl++)
+          TMP[kl][ij] = X[k][l];
+    }
+
+  zero_array(TEI, (nao*(nao+1)/2)*((nao*(nao+1)/2)+1)/2);
+
+  for(k=0,kl=0; k < nao; k++)
+    for(l=0; l <= k; l++,kl++) {
+      zero_matrix(X, nao, nao);
+      zero_matrix(Y, nao, nao);
+      for(i=0,ij=0; i < nao; i++)
+        for(j=0; j <=i; j++,ij++)
+          X[i][j] = X[j][i] = TMP[kl][ij];
+      zero_matrix(Y, nao, nao);
+      mmult(C, 1, X, 0, Y, nao, nao, nao);
+      zero_matrix(X, nao, nao);
+      mmult(Y, 0, C, 0, X, nao, nao, nao);
+      for(i=0, ij=0; i < nao; i++)
+        for(j=0; j <= i; j++,ij++) {
+          klij = INDEX(kl,ij);
+          TEI[klij] = X[i][j];
+        }
+    }
+
+  ...
+
+  free_matrix(X, nao);
+  free_matrix(Y, nao);
+  free_matrix(TMP, (nao*(nao+1)/2));
+
+```
+
diff --git a/Project#04/hints/hint3.md b/Project#04/hints/hint3.md
@@ -0,0 +1,20 @@
+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*:
+
+```c++
+  Emp2 = 0.0;
+  for(i=0; i < ndocc; i++) {
+    for(a=ndocc; a < nao; a++) {
+      ia = INDEX(i,a);
+      for(j=0; j < ndocc; j++) {
+        ja = INDEX(j,a);
+        for(b=ndocc; b < nao; b++) {
+          jb = INDEX(j,b);
+          ib = INDEX(i,b);
+          iajb = INDEX(ia,jb);
+          ibja = INDEX(ib,ja);
+          Emp2 += TEI[iajb] * (2 * TEI[iajb] - TEI[ibja])/(eps[i] + eps[j] - eps[a] - eps[b]);
+        }
+      }
+    }
+  }
+```
diff --git a/Project#04/project4-instructions.pdf b/Project#04/project4-instructions.pdf
diff --git a/Project#05/README.md b/Project#05/README.md
@@ -3,75 +3,50 @@ The coupled cluster model provides a higher level of accuracy beyond the MP2 app
 ## Step #1: Preparing the Spin-Orbital Basis Integrals 
 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> 
 amplitudes as well as Fock-matrix elements and antisymmetrized, Dirac-notation two-electron integrals, all given in the molecular spin-orbital basis 
-(as opposed to the spatial-orbital basis used in the earlier [MP2 Project](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2304).
+(as opposed to the spatial-orbital basis used in the earlier [MP2 Project](../Project%2304).
 Thus, the transformation of the AO-basis integrals into the spatial-MO basis must also include their translation into the spin-orbital basis:
 
-```
-EQUATION
-\begin{eqnarray*}
-\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) \\
-& \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) \\
-& \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)
-\end{eqnarry*}
-```
+<img src="./figures/spin-orbital-eri.png" height="160">
 
-
-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.
+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.
 
 Don't forget that you must also create the spin-orbital Fock matrix:
 
-```
-EQUATION
-f_{pq} = h_{pq} + \sum_m^{\rm occ} \langle pm || qm \rangle
-```
+<img src="./figures/spin-orbital-fock.png" height="60">
 
 Suggestion: For simplicity, store the two-electron integrals in a four-dimensional array.  This will greatly facilitate debugging of the complicated CCSD equations.
 
-  * Hint: Sample spatial- to spin-orbital translation code.
+  * [Hint 1](./hints/hint1.md): Sample spatial- to spin-orbital translation code.
 
 ## Step #2: Build the Initial-Guess Cluster Amplitudes 
 For Hartree-Fock reference determinants, the most common initial guess for the cluster amplitudes are the Moller-Plesset first-order perturbed wave function:
-```
-EQUATION
-t_i^a = 0
-```
-
-```
-EQUATION
-t_{ij}^{ab} = \frac{\langle ij || ab \rangle}{\epsilon_i + \epsilon_j - \epsilon_a - \epsilon_b}
 
-```
+<img src="./figures/init-t-amps.png" height="120">
 
 Note that if you have constructed the Fock matrix, two-electron integrals, and initial-guess amplitudes correctly at this point, 
-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] (https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2304):
+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):
 
-```
-EQUATION
-E_{\rm MP2} = \frac{1}{4} \sum_{ijab} \langle ij||ab\rangle t_{ij}^{ab}
-```
+<img src="./figures/mp2-energy.png" height="60">
 
 ## Step #3: Calculate the CC Intermediates 
-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).
+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;).
 
 ## Step #4: Compute the Updated Cluster Amplitudes 
 Use Eqs. 1 and 2 from Stanton's paper to compute the updated T<sub>1</sub> and T<sub>2</sub> cluster amplitudes.
 
 ## Step #5: Check for Convergence and Iterate 
 Calculate the current CC correlation energy:
 
-```
-EQUATION
-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
-```
-Compare energies and cluster amplitudes (using RMS differences) between iterations to check for convergence to some specified cutoff.  
+<img src="./figures/cc-correlation-energy.png" height="60">
+
+Compare energies and cluster amplitudes (using RMS differences) between iterations to check for convergence to some specified cutoff.
 If convergence is reached, you're done; if not, return to Step #3 and continue.
 
 ## Test Cases 
 The input structures, integrals, etc. for these examples are found in the 
-[input directory](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2305/input).
+[input directory](./input).
 
-* STO-3G Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2305/output/h2o/STO-3G/output.txt)
-* DZ Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2305/output/h2o/DZ/output.txt)
-* DZP Water | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2305/output/h2o/DZP/output.txt)
-* STO-3G Methane | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2305/output/ch4/STO-3G/output.txt)
+* STO-3G Water | [output](./output/h2o/STO-3G/output.txt)
+* DZ Water | [output](./output/h2o/DZ/output.txt)
+* DZP Water | [output](./output/h2o/DZP/output.txt)
+* STO-3G Methane | [output](./output/ch4/STO-3G/output.txt)
diff --git a/Project#05/figures/cc-correlation-energy.png b/Project#05/figures/cc-correlation-energy.png
diff --git a/Project#05/figures/init-t-amps.png b/Project#05/figures/init-t-amps.png
diff --git a/Project#05/figures/mp2-energy.png b/Project#05/figures/mp2-energy.png
diff --git a/Project#05/figures/spin-orbital-eri.png b/Project#05/figures/spin-orbital-eri.png
diff --git a/Project#05/figures/spin-orbital-fock.png b/Project#05/figures/spin-orbital-fock.png
diff --git a/Project#05/hints/hint1.md b/Project#05/hints/hint1.md
@@ -0,0 +1,19 @@
+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).
+
+```c++
+  /* Translate integrals to spin-orbital basis */
+  for(p=0; p < nmo; p++)
+    for(q=0; q < nmo; q++)
+      for(r=0; r < nmo; r++)
+        for(s=0; s < nmo; s++) {
+          pr = INDEX(p/2,r/2);
+          qs = INDEX(q/2,s/2);
+          prqs = INDEX(pr,qs);
+          value1 = TEI[prqs] * (p%2 == r%2) * (q%2 == s%2);
+          ps = INDEX(p/2,s/2);
+          qr = INDEX(q/2,r/2);
+          psqr = INDEX(ps,qr);
+          value2 = TEI[psqr] * (p%2 == s%2) * (q%2 == r%2);
+          ints[p][q][r][s] = value1 - value2;
+        }
+```