Add Project 9 hints and eqns

kcpearce · kcpearce · commit 642c6833cc2a · 2017-04-01T14:11:12.000-04:00
diff --git a/Project#09/README.md b/Project#09/README.md
@@ -3,34 +3,31 @@
 The use of point-group symmetry in quantum chemical calculations can lead to exceptionally efficient programs, 
 both in terms of their memory/disk usage and their computing times.  If a molecule has even one element of symmetry, 
 for example, this can reduce the number of wave function parameters that must be stored by up to a factor of two (the order of the point group) 
-and the computing effort by up to a factor of four (the square of the order of the point group).   
+and the computing effort by up to a factor of four (the square of the order of the point group).
 Most quantum chemical programs can take advantage of Abelian symmetry -- more accurately, 
 D<sub>2h</sub> and its subgroups -- which means that higher symmetry structures can yield up to a factor of 64 reduction in the computational costs.
 
-The purpose of this project is to take advantage of basic point-group symmetry in an SCF calculation.  
+The purpose of this project is to take advantage of basic point-group symmetry in an SCF calculation.
 While the most efficient programs will store only the symmetry-allowed sets of two-electron integrals, 
 for the purposes of this project we'll deal only with the one-electron (two-index) components.
 
 ## Background
 
 The vanishing integral rule of point-group theory may be stated as follows for Abelian groups:
 
-```
-EQUATION
-\langle \phi_i | \hat{A} |\phi_j \rangle = 0 \ \ \ \ {\rm if} \ \ \ \ \Gamma_A \neq \Gamma_{\phi_i} \otimes \Gamma_{\phi_j},
-```
+<img src="./figures/point-group-rule.png" height="20">
 
-where <html>&Gamma;<sub>X</sub></html> denotes the irreducible representation (irrep) of entity X.  
+where <html>&Gamma;<sub>X</sub></html> denotes the irreducible representation (irrep) of entity X.
 That is, unless the direct product of the irreps of the functions or operators in the integrand contains the totally symmetric irrep, the integral must be zero.
 
 The vanishing integral rule of group theory may be readily applied to the various integrals and wave-function parameters in the SCF procedure 
-** *if* ** the AO-basis functions over which the they are evaluated have been properly symmetrized 
+<b><i>if</i></b> the AO-basis functions over which the they are evaluated have been properly symmetrized 
 [i.e. used to form symmetry-adapted linear combinations, normally referred to as symmetry orbitals (SOs)].
 
-For example, consider the H<sub>2</sub>O test case with the STO-3G basis set as given in Project #3.  The AO 
+For example, consider the H<sub>2</sub>O test case with the STO-3G basis set as given in [Project #3](../Project%2303).  The AO 
 ([contracted Gaussian](http://sirius.chem.vt.edu/wiki/lib/exe/fetch.php?media=basis_sets.pdf)) basis functions 
 on the oxygen atom consist of 1s, 2s, 2p<sub>x</sub>, 2p<sub>y</sub>, and 2p<sub>z</sub> orbitals, 
-which those for each hydrogen atom consist of only the 1s, giving a total of seven AOs.  
+which those for each hydrogen atom consist of only the 1s, giving a total of seven AOs.
 If we take oxygen to be the first atom and ignore the C<sub>2v</sub> point group symmetry of the molecule, 
 the resulting overlap integral matrix looks like:
 
@@ -48,20 +45,20 @@ the resulting overlap integral matrix looks like:
 The first five row/column indices correspond to the basis functions on the oxygen (in the order listed above),
 while rows/columns six and seven correspond to the basis functions on each hydrogen. 
 Thus, the 1-2 element is the overlap integral between the 1s and 2s orbitals on the oxygen, 
-while element 3-7 is the overlap integral between the oxygen 2p<sub>x</sub> orbital and the second hydrogen 1s orbital.  
+while element 3-7 is the overlap integral between the oxygen 2p<sub>x</sub> orbital and the second hydrogen 1s orbital.
 
 Points to note about this matrix:
-  - The 1s and 2s orbitals on oxygen are not orthogonal to one another, but both are orthogonal to all three p-type orbitals.  
+  - The 1s and 2s orbitals on oxygen are not orthogonal to one another, but both are orthogonal to all three p-type orbitals.
   This occurs because, in Cartesian Gaussian basis sets, 
   orthogonality is generally obtained only between functions of different angular momenta on the same atomic center.
   - The 1s orbitals on the hydrogens have the same overlap integral (to within a sign) 
-  with all the orbitals on the oxygen because they are ** *symmetry equivalent* ** .  
+  with all the orbitals on the oxygen because they are <b><i>symmetry equivalent</i></b>.
   However, the 1s orbitals alone do not transform as irreps of the C<sub>2v</sub> point group.
-  - The hydrogen 1s orbitals are orthogonal to the 2p<sub>z</sub> orbital on the oxygen (cf. elements 5-6 and 5-7).  
+  - The hydrogen 1s orbitals are orthogonal to the 2p<sub>z</sub> orbital on the oxygen (cf. elements 5-6 and 5-7).
   This is because the water molecule was arbitrarily oriented in the xy-plane to compute the above integrals; 
   thus the 2p<sub>z</sub> orbital is perpendicular to the molecular plane containing the hydrogen 1s orbitals.
 
-What happens if we turn on symmetry, i.e. take linear combinations of the AO to obtain SOs?  
+What happens if we turn on symmetry, i.e. take linear combinations of the AO to obtain SOs?
 If we orient the molecule oriented in the yz-plane, so that the z-axis corresponds to the C<sub>2</sub> rotation axis, 
 we find that the AOs transform as the following irreps:
 
@@ -75,7 +72,7 @@ we find that the AOs transform as the following irreps:
 |  H<sub>1</sub> 1s + H<sub>2</sub> 1s  |  A<sub>1</sub>  |
 |  H<sub>1</sub> 1s - H<sub>2</sub> 1s  |  B<sub>2</sub>  |
 
-Thus, after symmetrization, there are four A<sub>1</sub> SOs, one B<sub>1</sub> SO, and two B<sub>2</sub> SOs.  
+Thus, after symmetrization, there are four A<sub>1</sub> SOs, one B<sub>1</sub> SO, and two B<sub>2</sub> SOs.
 (The STO-3G basis set contains no A<sub>2</sub> orbitals.  One must include at least d-type functions on the oxygen and/or 
 p-type functions on the hydrogens to obtain A<sub>2</sub> orbitals.)  Note that a normalization factor of 1/sqrt(2) should be included for the H 1s SOs above.
 
@@ -107,13 +104,13 @@ The PSI checkpoint file contains a great deal of useful symmetry information. To
 
 ## Storing One-Electron Integrals with Symmetry
 
-The block-diagonal structure above implies that one need only store and use the sub-matrices of the integrals rather than the full matrix.  
+The block-diagonal structure above implies that one need only store and use the sub-matrices of the integrals rather than the full matrix.
 One way of doing this is to store the one-electron integrals as an array of smaller matrices, 
 the size of which is dictated by the number of SOs per irrep.  For example, here's a code snippet that will extract the one-electron integrals 
-from the appropriate file (as was done in [Project #7](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2307) ), 
+from the appropriate file (as was done in [Project #7](../Project%2307) ), 
 but stores them in an array of matrices:
 
-```cpp
+```c++
   int ntri, nirreps;
   double *scratch;
   double ***S;
@@ -138,9 +135,9 @@ but stores them in an array of matrices:
 ```
 
 Notice how S is now a `double ***` -- i.e. an array of matrices.  This code uses the same `INDEX()` function found in 
-[Project #3](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2303), 
-and the sopi array is the integer array of SOs per irrep obtained using the `chkpt_rd_sopi()` function.  
-The `symm_offset` array is useful for converting between ** *relative* ** orbital indices and ** *absolute* ** orbital indices.  
+[Project #3](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2303),
+and the sopi array is the integer array of SOs per irrep obtained using the `chkpt_rd_sopi()` function.
+The `symm_offset` array is useful for converting between <b><i>relative</i></b> orbital indices and <b><i>absolute</i></b> orbital indices.
 (For example, orbital 0 in the B<sub>2</sub> irrep block corresponds to orbital 5 
 in the full list of orbitals; the first number is the relative index, the second is the absolute index.)
 
@@ -149,28 +146,23 @@ in the full list of orbitals; the first number is the relative index, the second
 If we store the integral, MO coefficient, density, and Fock matrices in the block-diagonal structure described above, 
 we may compute matrix products of these quantities separately for each diagonal block.  For example, the matrix product may be written as:
 
-```
-EQUATION
-{\mathbf C} = {\mathbf A} {\mathbf B} = \bigoplus_h {\mathbf C}_h = \bigoplus_h {\mathbf A}_h {\mathbf B}_h,
-```
+<img src="./figures/matrix-product.png" height="25">
 
-where *h* denotes a particular irrep of the point group and <b><i>C</i></b><sub>h</sub> denotes the *h*-th irrep subblock of the full matrix ** *C* ** .  
+where *h* denotes a particular irrep of the point group and <b><i>C</i></b><sub>h</sub> denotes the *h*-th irrep subblock of the full matrix <b><i>C</i></b>.
 (The notation of the large plus with a circle around it indicates a direct sum.)
 
 ## Occupied Orbitals and the Density Matrix
 
 The calculation of the density matrix requires information about the number of occupied orbitals:
 
-```
-EQUATION
-D_{\mu\nu} = \sum_m^{\rm occ.} \left( {\mathbf C} \right)_\mu^m \left( {\mathbf C} \right)_\nu^m
-```
+<img src="./figures/density-matrix.png" height="25">
+
 In a calculation without symmetry, one need only know the number of electrons (equivalently, the atomic numbers of the atoms and the overall molecular charge) 
 to evaluate this expression.  However, if the density and MO coefficient matrices are stored in symmetry-blocked form, 
-as above, one needs to know the number of occupied MOs in each irrep.  
+as above, one needs to know the number of occupied MOs in each irrep.
 
 For the case of the C<sub>2v</sub> water molecule, which has two <html>&sigma;</html> O-H bonds, two oxygen lone pairs, 
-and one oxygen core orbital, the orbital occupations are not difficult to determine, given the above information about the number of SOs in each irrep.  
+and one oxygen core orbital, the orbital occupations are not difficult to determine, given the above information about the number of SOs in each irrep.
 Since the MOs are constructed as linear combinations of the SOs, and only SOs of the same irrep can combine, we may rationalize the occupied orbitals of 
 H<sub>2</sub>O as follows:
 
@@ -184,9 +176,9 @@ H<sub>2</sub>O as follows:
 
 Thus, for the water molecule, there are 3 A<sub>1</sub>, 1 B<sub>1</sub>, and 1 B<sub>2</sub> occupied orbitals.
 
-Determining the occupations automatically is sometimes a difficult task, especially for molecules with unusual bonding patterns.  
+Determining the occupations automatically is sometimes a difficult task, especially for molecules with unusual bonding patterns.
 Many SCF programs will simply take the core Hamiltonian matrix (often used as the initial guess for the Fock matrix), 
-diagonalize it (in the othogonal AO basis) in each irrep block, and identify the lowest eigenvalues in each irrep as occupied.  
+diagonalize it (in the othogonal AO basis) in each irrep block, and identify the lowest eigenvalues in each irrep as occupied.
 However, this approach often fails for many molecules (especially those with open shells) because of near-degeneracies 
 that may not be adequately described by the simple core Hamiltonian.  Thus, some programs will make use of more 
 sophisticated guesses at the Fock matrix, e.g. Hueckel orbitals or other (semi)empirical forms.
@@ -195,17 +187,14 @@ sophisticated guesses at the Fock matrix, e.g. Hueckel orbitals or other (semi)e
 
 The only portion of the SCF procedure involving two-electron integrals is the Fock-matrix build:
 
-```
-EQUATION
-F_{\mu\nu} = H^{\rm core}_{\mu\nu} + \sum_{\lambda\sigma}^{\rm AO} D^{i-1}_{\lambda\sigma} \left[ 2 (\mu\nu | \lambda\sigma) - (\mu\lambda|\nu\sigma) \right],
-```
+<img src="./figures/fock-build.png" height="30">
 
-where the notation used in [Project #3](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2303) is retained above.  
+where the notation used in [Project #3](../Project%2303) is retained above.
 If the Fock matrix and density matrix are stored in a symmetry-blocked manner similar to the overlap matrix above, 
 then one may limit the corresponding loops only to the symmetry-allowed combinations, e.g. the irreps of <html>&mu;</html> and 
 <html>&nu;</html> must match for F<html><sub>&mu;&nu;</sub></html> to be non-zero and the irreps of <html>&lambda;</html> 
 and <html>&sigma;</html> must match for D<html><sub>&lambda;&sigma;</sub></html> to be non-zero.  
 Note, however, that, because of the symmetries in the two-electron integrals, 
 the irrep of <html>&lambda;</html> and <html>&sigma;</html> do not necessarily have to match that of <html>&mu;</html> and <html>&nu;</html>.
 
-  * Hint: Fock-Build Code
+  * [Hint](./hints/hint1.md): Fock-Build Code
diff --git a/Project#09/hint1.md b/Project#09/hint1.md
@@ -0,0 +1,28 @@
+```c++
+    for(h=0; h < nirreps; h++) {
+      for(i=0; i < sopi[h]; i++) {
+        I = symm_offset[h] + i;
+        for(j=0; j < sopi[h]; j++) {
+          J = symm_offset[h] + j;
+          F[h][i][j] = H[h][i][j];
+
+          for(g=0; g < nirreps; g++) {
+            for(k=0; k < sopi[g]; k++) {
+              K = symm_offset[g] + k;
+              for(l=0; l < sopi[g]; l++) {
+                L = symm_offset[g] + l;
+
+                ij = INDEX(I,J); kl = INDEX(K,L);
+                ijkl = INDEX(ij,kl);
+                ik = INDEX(I,K); jl = INDEX(J,L);
+                ikjl = INDEX(ik,jl);
+
+                F[h][i][j] += D[g][k][l] * (2.0 * TEI[ijkl] - TEI[ikjl]);
+              }
+            }
+
+          }
+        }
+      }
+    }
+```
