Check some Project 1 hints

Author: kcpearce

Project#01/hints/hint4-1.md

@@ -0,0 +1 @@
+Note that, unlike the interatomic distance matrix and the three-dimensional bond angle array, the out-of-plane angles are not needed later in the program.  Hence, there is no need to allocate any memory to store these angles. Simply print them one at a time.

Project#01/hints/hint4-2.md

@@ -0,0 +1,10 @@
+A cross product between two vectors is another vector. Therefore, to store the cross product needed to compute a given out-of-plane angle, you need to store the x, y, and z components of the resulting vector, which you then use to compute a dot product with yet another vector. It is convenient to assign separate variables for each part of this calculation:
+```c++
+ejkl_x = (ey[k][j] * ez[k][l] - ez[k][j] * ey[k][l]);
+ejkl_y = (ez[k][j] * ex[k][l] - ex[k][j] * ez[k][l]);
+ejkl_z = (ex[k][j] * ey[k][l] - ey[k][j] * ex[k][l]);
+ 
+exx = ejkl_x * ex[k][i];
+eyy = ejkl_y * ey[k][i];
+ezz = ejkl_z * ez[k][i];
+```

Project#01/hints/hint4-3.md

@@ -0,0 +1,8 @@
+The final out-of-plane angle is computed using the `asin()` function. Note that since sine yields values between -1 and +1, the `asin()` function can only take arguments between -1.0 and +1.0. However, numerical precision in the calculation of the cross- and dot-products earlier in the calculation can yield results slightly outside this domain. Hence, you should test the argument of `asin()` before you call it:
+```c++
+  theta = (exx + eyy + ezz)/sin(phi[j][k][l]);
+ 
+  if(theta < -1.0) theta = asin(-1.0);
+  else if(theta > 1.0) theta = asin(1.0);
+  else theta = asin(theta);
+```

Project#01/hints/hint4-4

@@ -0,0 +1,8 @@
+Just as for the bond angle code, we need to exclude ijkl combinations involving coincidences among the indices as well as distant atom pairs:
+```c++
+  if(i!=j && i!=k && i!=l && j!=k && j!=l && k!=l) { } /* Skip coincidences */
+```
+and
+```c++
+  if(R[i][k] < 4.0 && R[k][j] < 4.0 && R[k][l] < 4.0) { } /* Skip distant atom pairs */
+```

Project#01/hints/hint5-1.md

@@ -0,0 +1 @@
+Just as for the out-of-plane angles, there is no need to store the torsional angles after they are printed, so no special memory allocation procedure is needed for this step.

Project#01/hints/hint5-2.md

@@ -0,0 +1 @@
+Just as for `asin()` in the out-of-plane angle calculation, you should check the argument to `acos()` to make sure it lies between -1.0 and +1.0.

Project#01/hints/hint5-3.md

@@ -0,0 +1 @@
+To print only the unique dihedral angles, it is straightforward in this case to limit the loop structure over i, j, k, and l to keep j < i, k < j, and l < k. Also, one should limit the printing only to atom pairs that are close together.

Project#01/hints/step4-solution.md

@@ -0,0 +1,174 @@
+We've added a new function to the Molecule class for computing the out-of-plane angles.  Here's the new 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 oop(int atom1, int atom2, int atom3, int atom4);
+    double unit(int cart, int atom1, int atom2);
+
+    Molecule(const char *filename, int q);
+    ~Molecule();
+};
+```
+
+Here's the new `oop()` function for molecule.cc:
+
+```c++
+double Molecule::oop(int a, int b, int c, int d)
+{
+  double ebcd_x = (unit(1,c,b) * unit(2,c,d) - unit(2,c,b) * unit(1,c,d));
+  double ebcd_y = (unit(2,c,b) * unit(0,c,d) - unit(0,c,b) * unit(2,c,d));
+  double ebcd_z = (unit(0,c,b) * unit(1,c,d) - unit(1,c,b) * unit(0,c,d));
+  
+  double exx = ebcd_x * unit(0,c,a);
+  double eyy = ebcd_y * unit(1,c,a);
+  double ezz = ebcd_z * unit(2,c,a);
+
+  double theta = (exx + eyy + ezz)/sin(angle(b,c,d));
+
+  if(theta < -1.0) theta = asin(-1.0);
+  else if(theta > 1.0) theta = asin(1.0);
+  else theta = asin(theta);
+
+  return theta;
+}
+```
+
+And finally the new code that makes use of the class:
+
+```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)));
+      }
+    }
+  }
+
+  cout << "\nOut-of-Plane angles:\n";
+  for(int i=0; i < mol.natom; i++) {
+    for(int k=0; k < mol.natom; k++) {
+      for(int j=0; j < mol.natom; j++) {
+        for(int l=0; l < j; l++) {
+          if(i!=j && i!=k && i!=l && j!=k && k!=l && mol.bond(i,k) < 4.0 && mol.bond(k,j) < 4.0 && mol.bond(k,l) < 4.0)
+              printf("%2d-%2d-%2d-%2d %10.6f\n", i, j, k, l, mol.oop(i,j,k,l)*(180.0/acos(-1.0)));
+        }
+      }
+    }
+  }
+
+  return 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
+
+Out-of-plane angles:
+ 0- 3- 1- 2  -0.000000
+ 0- 6- 4- 5  19.939726
+ 0- 6- 5- 4 -19.850523
+ 0- 5- 6- 4  19.850523
+ 1- 5- 0- 4  53.678778
+ 1- 6- 0- 4 -53.678778
+ 1- 6- 0- 5  54.977064
+ 2- 3- 1- 0   0.000000
+ 3- 2- 1- 0  -0.000000
+ 4- 5- 0- 1 -53.651534
+ 4- 6- 0- 1  53.651534
+ 4- 6- 0- 5 -54.869992
+ 4- 6- 5- 0  29.885677
+ 4- 5- 6- 0 -29.885677
+ 5- 4- 0- 1  53.626323
+ 5- 6- 0- 1 -56.277112
+ 5- 6- 0- 4  56.194621
+ 5- 6- 4- 0 -30.558964
+ 5- 4- 6- 0  31.064344
+ 6- 4- 0- 1 -53.626323
+ 6- 5- 0- 1  56.277112
+ 6- 5- 0- 4 -56.194621
+ 6- 5- 4- 0  30.558964
+ 6- 4- 5- 0 -31.064344
+```