First set of hints from Project#1

Author: kcpearce

Project#01/hints/hint2-2.md

@@ -6,11 +6,9 @@ To build the distance matrix, we need a loop for each index:
 
   for(int i=0; i < mol.natom; i++) {
     for(int j=0; j < mol.natom; j++) {
-        R[i][j] = sqrt(
-            (mol.geom[i][0]-mol.geom[j][0])*(mol.geom[i][0]-mol.geom[j][0])
-          + (mol.geom[i][1]-mol.geom[j][1])*(mol.geom[i][1]-mol.geom[j][1])
-          + (mol.geom[i][2]-mol.geom[j][2])*(mol.geom[i][2]-mol.geom[j][2])
-                      );
+        R[i][j] = sqrt( (mol.geom[i][0]-mol.geom[j][0])*(mol.geom[i][0]-mol.geom[j][0])
+                      + (mol.geom[i][1]-mol.geom[j][1])*(mol.geom[i][1]-mol.geom[j][1])
+                      + (mol.geom[i][2]-mol.geom[j][2])*(mol.geom[i][2]-mol.geom[j][2]) );
     }
   }
 ```

Project#01/hints/hint2-3.md

@@ -0,0 +1,7 @@
+To print the interatomic distance matrix, we could just run through its unique values:
+```c++
+  for(int i=0; i < mol.natom; i++)
+    for(int j=0; j < i; j++)
+      printf("%d %d %8.5f\n", i, j, R[i][j]);
+```
+Note the conditional on the index `j` which keeps the code form printing redundant information.

Project#01/hints/hint2-4.md

@@ -0,0 +1,11 @@
+Now that we've written some code to compute bond distances, we can extend the Molecule class a bit more to define our `bond()` function:
+```c++
+double Molecule::bond(int a, int b)
+{
+  return sqrt( (geom[a][0]-geom[b][0])*(geom[a][0]-geom[b][0])
+             + (geom[a][1]-geom[b][1])*(geom[a][1]-geom[b][1])
+             + (geom[a][2]-geom[b][2])*(geom[a][2]-geom[b][2]) );
+}
+```
+
+Note that since `bond()` is a member function of the Molecule class, we can access the geom array with just `geom` rather than `mol.geom`

Project#01/hints/hint3-1.md

@@ -0,0 +1,18 @@
+Each unit vector points from one atom to another, hence each Cartesian component should be treated as a matrix, much like the interatomic distance matrix.  For now, let's store the unit vectors in memory, but later we'll write a function to recompute them as needed:
+```c++
+  double **ex = new double* [mol.natom];
+  double **ey = new double* [mol.natom];
+  double **ez = new double* [mol.natom];
+  for(int i=0; i < mol.natom; i++) {
+    ex[i] = new double[mol.natom];
+    ey[i] = new double[mol.natom];
+    ez[i] = new double[mol.natom];
+  }
+```
+And don't forget to delete[] them at the end:
+```c++
+  for(int i=0; i < mol.natom; i++) {
+    delete[] ex[i]; delete[] ey[i]; delete[] ez[i];
+  }
+  delete[] ex; delete[] ey; delete[] ez;
+```

Project#01/hints/hint3-2.md

@@ -0,0 +1,10 @@
+Be careful about the diagonal elements of the unit vector matrices because of the division by the bond distance between atoms i and j.  You should restrict the loops to skip those elements involving a divide-by-zero condition!
+```c++
+  for(i=0; i < mol.natom; i++) {
+    for(j=0; j < i; j++) {
+      ex[i][j] = ex[j][i] = -(x[i] - x[j])/R[i][j];
+      ey[i][j] = ey[j][i] = -(y[i] - y[j])/R[i][j];
+      ez[i][j] = ez[j][i] = -(z[i] - z[j])/R[i][j];
+    }
+  }
+```

Project#01/hints/hint3-3.md

@@ -0,0 +1,10 @@
+If you choose to store the bond angles for later use (not absolutely necessary, as we'll see), you need a three-dimensional array:
+```c++
+  double ***phi = new double** [mol.natom];
+  for(int i=0; i < mol.natom; i++) {
+    phi[i] = new double* [mol.natom];
+    for(int j=0; j < mol.natom; j++) {
+      phi[i][j] = new double[mol.natom];
+    }
+  }
+```

Project#01/hints/hint3-4.md

@@ -0,0 +1,22 @@
+If you print the angle between every possible combination of three atoms, you'll get lots of zeroes and angles between atoms that are far apart.   To print mainly interesting angles, use if-else blocks to enforce these restrictions:
+```c++
+if(i!=j && i!=k && j!=k) { }  /* Skip coincidences */
+```
+
+and
+```c++
+if(i < j && j < k && R[i][j] < 4.0 && R[j][k] < 4.0) { } /* Skip atoms far apart (specifically with bond distances > 4.0 bohr) */
+```
+
+Alternatively, you can limit the loop structure and just filter out the bond distances:
+```c++
+for(int i=0; i < mol.natom; i++) {
+  for(int j=0; j < i; j++) {
+    for(int k=0; k < j; k++) {
+      if(R[i][j] < 4.0 && R[j][k] < 4.0) {
+         ...
+      }
+    }
+  }
+}
+```

Project#01/hints/hint3-5.md

@@ -0,0 +1 @@
+To compute the final bond angle, you must use the `acos()` function (arccos), which, like the `sqrt()` function, is part of the [C math library](http://en.wikipedia.org/wiki/Math.h).

Project#01/hints/step2-solution.md

@@ -0,0 +1,63 @@
+Here we use the Molecule class we've been working on so far, and thanks to the new "bond()" function, we don't need to worry about storing the distances in a matrix:
+
+```c++
+#include "molecule.h"
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include <cstdio>
+#include <cmath>
+
+using namespace std;
+ 
+int main()
+{
+  Molecule mol("geom.dat", 0);
+ 
+  cout << "Number of atoms: " << mol.natom << endl;
+  cout << "Input Cartesian coordinates:\n";
+  mol.print_geom();
+ 
+  cout << "Interatomic distances (bohr):\n";
+  for(int i=0; i < mol.natom; i++)
+    for(int j=0; j < i; j++)
+      printf("%d %d %8.5f\n", i, j, mol.bond(i,j));
+      
+  return 0;
+}
+```
+
+The output from the above program for the acetaldehyde test case is:
+```
+Number of atoms: 7
+Input Cartesian coordinates:
+6       0.000000000000       0.000000000000       0.000000000000
+6       0.000000000000       0.000000000000       2.845112131228
+8       1.899115961744       0.000000000000       4.139062527233
+1      -1.894048308506       0.000000000000       3.747688672216
+1       1.942500819960       0.000000000000      -0.701145981971
+1      -1.007295466862      -1.669971842687      -0.705916966833
+1      -1.007295466862       1.669971842687      -0.705916966833
+Interatomic distances (bohr):
+1 0  2.84511
+2 0  4.55395
+2 1  2.29803
+3 0  4.19912
+3 1  2.09811
+3 2  3.81330
+4 0  2.06517
+4 1  4.04342
+4 2  4.84040
+4 3  5.87463
+5 0  2.07407
+5 1  4.05133
+5 2  5.89151
+5 3  4.83836
+5 4  3.38971
+6 0  2.07407
+6 1  4.05133
+6 2  5.89151
+6 3  4.83836
+6 4  3.38971
+6 5  3.33994
+```

Project#01/hints/step3-solution.md

@@ -0,0 +1,128 @@
+Here again we've extended the Molecule class to compute the bond angles without storing them.  This also required the addition of a new function to compute the Cartesian unit vectors on the fly, whose declaration we add to the member functions of molecule.h:
+```c++
+#include <string>
+
+using namespace std;
+
+class Molecule
+{
+  public:
+    int natom;
+    int charge;
+    int *zvals;
+    double **geom;
+    string point_group;
+
+    void print_geom();
+    void rotate(double phi);
+    void translate(double x, double y, double z);
+    double bond(int atom1, int atom2);
+    double angle(int atom1, int atom2, int atom3);
+    double torsion(int atom1, int atom2, int atom3, int atom4);
+    double unit(int cart, int atom1, int atom2);
+
+    Molecule(const char *filename, int q);
+    ~Molecule();
+};
+```
+
+Now the two new functions added to molecule.cc:
+```c++
+// Returns the value of the unit vector between atoms a and b
+// in the cart direction (cart=0=x, cart=1=y, cart=2=z)
+double Molecule::unit(int cart, int a, int b)
+{
+  return -(geom[a][cart]-geom[b][cart])/bond(a,b);
+}
+
+// Returns the angle between atoms a, b, and c in radians
+double Molecule::angle(int a, int b, int c)
+{
+  return acos(unit(0,b,a) * unit(0,b,c) + unit(1,b,a) * unit(1,b,c) + unit(2,b,a) * unit(2,b,c));
+}
+```
+
+And finally, the code that makes use of the above: 
+
+```c++
+#include "molecule.h"
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include <cstdio>
+#include <cmath>
+
+using namespace std;
+ 
+int main()
+{
+  Molecule mol("geom.dat", 0);
+ 
+  cout << "Number of atoms: " << mol.natom << endl;
+  cout << "Input Cartesian coordinates:\n";
+  mol.print_geom();
+  
+  cout << "Interatomic distances (bohr):\n";
+  for(int i=0; i < mol.natom; i++)
+    for(int j=0; j < i; j++)
+      printf("%d %d %8.5f\n", i, j, mol.bond(i,j));
+
+  cout << "\nBond angles:\n";
+  for(int i=0; i < mol.natom; i++) {
+    for(int j=0; j < i; j++) {
+      for(int k=0; k < j; k++) {
+        if(mol.bond(i,j) < 4.0 && mol.bond(j,k) < 4.0)
+          printf("%2d-%2d-%2d %10.6f\n", i, j, k, mol.angle(i,j,k)*(180.0/acos(-1.0)));
+      }
+    }
+  }
+
+  return 0;
+}
+```
+Note that the value of &pi; is obtained with machine precision using `acos(-1.0)`.
+
+The above code produces the following output for the acetaldehyde test case:
+
+```
+Number of atoms: 7
+Input Cartesian coordinates:
+6       0.000000000000       0.000000000000       0.000000000000
+6       0.000000000000       0.000000000000       2.845112131228
+8       1.899115961744       0.000000000000       4.139062527233
+1      -1.894048308506       0.000000000000       3.747688672216
+1       1.942500819960       0.000000000000      -0.701145981971
+1      -1.007295466862      -1.669971842687      -0.705916966833
+1      -1.007295466862       1.669971842687      -0.705916966833
+Interatomic distances (bohr):
+1 0  2.84511
+2 0  4.55395
+2 1  2.29803
+3 0  4.19912
+3 1  2.09811
+3 2  3.81330
+4 0  2.06517
+4 1  4.04342
+4 2  4.84040
+4 3  5.87463
+5 0  2.07407
+5 1  4.05133
+5 2  5.89151
+5 3  4.83836
+5 4  3.38971
+6 0  2.07407
+6 1  4.05133
+6 2  5.89151
+6 3  4.83836
+6 4  3.38971
+6 5  3.33994
+
+Bond angles:
+ 2- 1- 0 124.268308
+ 3- 1- 0 115.479341
+ 3- 2- 1  28.377448
+ 5- 4- 0  35.109529
+ 6- 4- 0  35.109529
+ 6- 5- 0  36.373677
+ 6- 5- 4  60.484476
+```