Project 1 hint/solutions complete

Author: kcpearce

Project#01/README.md

@@ -1,10 +1,10 @@
 # TODO:
 - [ ] Add embedded equations
-- [ ] Add hints and solutions for each step
+- [x] Add hints and solutions for each step
 - [x] Add input/output links for test cases
 
 # Molecular-Geometry-Analysis
-The purpose of this project is to introduce you to fundamental C-language (or C++) programming techniques in the context of a scientific problem, viz. the calculation of the internal coordinates (bond lengths, bond angles, dihedral angles), moments of inertia, and rotational constants of a polyatomic molecule. A concise set of instructions for this project may be found [here](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2301/project1-instructions.pdf).
+The purpose of this project is to introduce you to fundamental C-language (or C++) programming techniques in the context of a scientific problem, viz. the calculation of the internal coordinates (bond lengths, bond angles, dihedral angles), moments of inertia, and rotational constants of a polyatomic molecule. A concise set of instructions for this project may be found [here](./project1-instructions.pdf).
 
 We thank Dr. Yukio Yamaguchi of the University of Georgia for the original version of this project.
 
@@ -20,14 +20,14 @@ The input to the program is the set of Cartesian coordinates of the atoms (in bo
     1 -1.007295466862    -1.669971842687    -0.705916966833
     1 -1.007295466862     1.669971842687    -0.705916966833
 
-The first line above is the number of atoms (an integer), while the remaining lines contain the z-values and x-, y-, and z-coordinates of each atom (one integer followed by three double-precision floating-point numbers). This [input file](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2301/input/acetaldehyde.dat) ("acetaldehyde.dat") along with a few other test cases can be found in this repository in the [input directory](https://github.com/CrawfordGroup/ProgrammingProjects/tree/master/Project%2301/input).
+The first line above is the number of atoms (an integer), while the remaining lines contain the z-values and x-, y-, and z-coordinates of each atom (one integer followed by three double-precision floating-point numbers). This [input file](./input/acetaldehyde.dat) ("acetaldehyde.dat") along with a few other test cases can be found in this repository in the [input directory](./input).
 
 After downloading the file to your computer (to a file called “geom.dat”, for example), you must open the file, read the data from each line into appropriate variables, and finally close the file.
 
-- Hint #1: Opening and closing the file stream
-- Hint #2: Reading the number of atoms
-- Hint #3: Storing the z-values and the coordinates
-- Solution
+- [Hint #1](./hints/hint1-1.md): Opening and closing the file stream
+- [Hint #2](./hints/hint1-2.md): Reading the number of atoms
+- [Hint #3](./hints/hint1-3.md): Storing the z-values and the coordinates
+- [Solution](./hints/step1-solution.md)
 
 ## Step 2: Bond Lengths
 Calculate the interatomic distances using the expression:
@@ -38,11 +38,11 @@ EQUATION
 
 where x, y, and z are Cartesian coordinates and i and j denote atomic indices.
 
-- Hint 1: Memory allocation
-- Hint 2: Loop structure
-- Hint 3: Printing the results
-- Hint 4: Extending the Molecule class
-- Solution
+- [Hint 1](./hints/hint2-1.md): Memory allocation
+- [Hint 2](./hints/hint2-2.md): Loop structure
+- [Hint 3](./hints/hint2-3.md): Printing the results
+- [Hint 4](./hints/hint2-4.md): Extending the Molecule class
+- [Solution](./hints/step2-solution.md)
 
 ## Step 3: Bond Angles
 Calculate all possible bond angles. For example, the angle, &phi;<sub>ijk</sub>, between atoms i-j-k, where j is the central atom is given by:
@@ -57,12 +57,12 @@ where the e&#8407;<sub>ij</sub> are unit vectors between the atoms, e.g.,
 EQUATION
 ```
 
-- Hint 1: Memory allocation for the unit vectors
-- Hint 2: Avoiding a divide-by-zero
-- Hint 3: Memory allocation for the bond angles
-- Hint 4: Smart printing
-- Hint 5: Trigonometric functions
-- Solution
+- [Hint 1](./hints/hint3-1.md): Memory allocation for the unit vectors
+- [Hint 2](./hints/hint3-2.md): Avoiding a divide-by-zero
+- [Hint 3](./hints/hint3-3.md): Memory allocation for the bond angles
+- [Hint 4](./hints/hint3-4.md): Smart printing
+- [Hint 5](./hints/hint3-5.md): Trigonometric functions
+- [Solution](./hints/step3-solution.md)
 
 ## Step 4: Out-of-Plane Angles
 Calculate all possible out-of-plane angles. For example, the angle &theta;<sub>ijkl</sub> for atom i out of the plane containing atoms j-k-l (with k as the central atom, connected to i) is given by:
@@ -71,11 +71,11 @@ Calculate all possible out-of-plane angles. For example, the angle &theta;<sub>i
 EQUATION
 ```
 
-- Hint 1: Memory allocation?
-- Hint 2: Cross products
-- Hint 3: Numerical precision
-- Hint 4: Smarter printing
-- Solution
+- [Hint 1](./hints/hint4-1.md): Memory allocation?
+- [Hint 2](./hints/hint4-2.md): Cross products
+- [Hint 3](./hints/hint4-3.md): Numerical precision
+- [Hint 4](./hints/hint4-4.md): Smarter printing
+- [Solution](./hints/step4-solution.md)
 
 ## Step 5: Torsion/Dihedral Angles
 Calculate all possible torsional angles. For example, the torsional angle &tau;<sub>ijkl</sub> for the atom connectivity i-j-k-l is given by:
@@ -86,11 +86,11 @@ EQUATION
 
 Can you also determine the sign of the torsional angle?
 
-- Hint 1: Memory allocation?
-- Hint 2: Numerical precision
-- Hint 3: Smart printing
-- Hint 4: Sign
-- Solution
+- [Hint 1](./hints/hint5-1.md): Memory allocation?
+- [Hint 2](./hints/hint5-2.md): Numerical precision
+- [Hint 3](./hints/hint5-3.md): Smart printing
+- [Hint 4](./hints/hint5-4.md): Sign
+- [Solution](./hints/step5-solution.md)
 
 ## Step 6: Center-of-Mass Translation
 Find the center of mass of the molecule:
@@ -103,9 +103,9 @@ where m<sub>i</sub> is the mass of atom i and the summation runs over all atoms
 
 Translate the input coordinates of the molecule to the center-of-mass.
 
-- Hint 1: Atomic masses
-- Hint 2: Translating between atomic number and atomic mass
-- Solution
+- [Hint 1](./hints/hint6-1.md): Atomic masses
+- [Hint 2](./hints/hint6-2.md): Translating between atomic number and atomic mass
+- [Solution](./hints/step6-solution.md)
 
 ## Step 7: Principal Moments of Inertia
 Calculate elements of the [moment of inertia tensor](http://en.wikipedia.org/wiki/Moment_of_inertia_tensor).
@@ -129,9 +129,9 @@ Report the moments of inertia in amu bohr<sup>2</sup>, amu &#8491;<sup>2</sup>,
 
 Based on the relative values of the principal moments, determine the [molecular rotor type](http://en.wikipedia.org/wiki/Rotational_spectroscopy): linear, oblate, prolate, asymmetric.
 
-- Hint 1: Diagonalization of a 3×3 matrix
-- Hint 2: Physical constants
-- Solution
+- [Hint 1](./hints/hint7-1.md): Diagonalization of a 3×3 matrix
+- [Hint 2](./hints/hint7-2.md): Physical constants
+- [Solution](./hints/step7-solution.md)
 
 ## Step 8: Rotational Constants
 Compute the rotational constants in cm<sup>-1</sup> and MHz:
@@ -140,13 +140,13 @@ Compute the rotational constants in cm<sup>-1</sup> and MHz:
 EQUATION
 ```
 
-- Solution
+- [Solution](./hints/step8-solution.md)
 
 
 ## Test Cases
-- Acetaldehyde: [input coordinates](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2301/input/acetaldehyde.dat) | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2301/output/acetaldehyde_out.txt)
-- Benzene: [input coordinates](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2301/input/benzene.dat) | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2301/output/benzene_out.txt)
-- Allene: [input coordinates](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2301/input/allene.dat) | [output](https://github.com/CrawfordGroup/ProgrammingProjects/blob/master/Project%2301/output/allene_out.txt)
+- Acetaldehyde: [input coordinates](./input/acetaldehyde.dat) | [output](./output/acetaldehyde_out.txt)
+- Benzene: [input coordinates](./input/benzene.dat) | [output](./output/benzene_out.txt)
+- Allene: [input coordinates](./input/allene.dat) | [output](./output/allene_out.txt)
 
 ## References
 E.B. Wilson, J.C. Decius, and P.C. Cross, __Molecular Vibrations__, McGraw-Hill, 1955.

Project#01/eigen.tar.gz

@@ -0,0 +1 @@
+The sign of a torsional/dihedral angle among atoms **_i-j-k-l_** is positive (negative) if the vector along **_k-l_** lies to the right (left) of the plane formed by **_i-j-k_** when the plane is viewed along the **_j-k_** vector.

Project#01/hints/hint5-4.md

@@ -0,0 +1 @@
+An excellent source for atomic masses and other physical constants is the [National Institute of Standard and Technology (NIST) website](http://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=&ascii=html&isotype=some).

Project#01/hints/hint6-1.md

@@ -0,0 +1,10 @@
+An elegant way to translate between the atomic number (z-value) of a given atom and its mass is to prepare a static array of the masses of the most abundant isotope of each element.  I suggest preparing a header file containing a [global array](https://github.com/CrawfordGroup/ProgrammingProjects/wiki/Variable-Scope-and-Reference-Types#global-variables), e.g.:
+```c++
+  double masses[] = {
+  0.0000000, 
+  1.007825, 
+  4.002603,  
+  6.015123,
+  ...};
+```
+Note that, for example, `masses[1] = 1.007825`, which is the correct atomic mass for a hydrogen atom.

Project#01/hints/hint6-2.md

@@ -0,0 +1,51 @@
+Here are two approaches for the [diagonalization](http://en.wikipedia.org/wiki/Diagonalizable_matrix) of the moment of inertia tensor:
+
+## Secular Determinant
+Since the moment of inertia tensor is only a 3x3 matrix, a brute-force approach via the secular determinant is feasible:
+
+```latex
+\left|\begin{array}{ccc}
+(I_{11}-\lambda) & I_{12} & I_{13} \\
+I_{21} & (I_{22}-\lambda) & I_{23} \\
+I_{31} & I_{32} & (I_{33}-\lambda) \\
+\end{array}\right| = 0
+```
+
+This leads to a cubic equation in λ, which one can solve directly.  Have fun with that.
+
+## More General Algorithms
+"Canned" algorithms are definitely the way to go for general matrix diagonalization.  Most such algorithms are based on a two-step procedure:
+  - Reduction of the matrix to a tridiagonal form using the [Householder](http://en.wikipedia.org/wiki/Householder's_method) or Givens approaches.
+  - Diagonalization of the tridiagonal structure, either by solving its secular determinant or by other methods, e.g. [QR or QL decompositions](http://en.wikipedia.org/wiki/QR_decomposition).
+
+A convenient canned library for a wide range of linear algebraic operations is the [Eigen package](http://eigen.tuxfamily.org).  This is a template-only library that provides a very clean interface for manipulating a large number of matrix types.  You can either download and install the library from  the [official website](http://eigen.tuxfamily.org) or just grab the [gzipped tarfile from here](../eigen.tar.gz).  Unpack the library in the same directory as your source code, and you're ready to get started.
+
+To use the library to diagonalize your moment inertia tensor, follow these steps:
+
+- Add the following lines to your main source file below the inclusion of other headers:
+```c++
+#include "Eigen/Dense"
+#include "Eigen/Eigenvalues"
+#include "Eigen/Core"
+
+typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Matrix;
+```
+This code defines a new type called a `Matrix` that may be dynamically allocated and contains only doubles.
+- Allocate your moment of inertia tensor by a line of code like:
+```c++
+Matrix I(3,3);
+```
+- Access or assign individual elements of the Matrix using parenthetical notation rather than brackets, e.g.:
+```c++
+I(0,0) = 2.0;
+```
+- The Eigen package makes it easy to examine your matrix using cout:
+```c++
+cout << I << endl;
+```
+- After you have built the moment of inertia tensor, you may compute its eigenvalues and eigenvectors as follows:
+```c++
+  Eigen::SelfAdjointEigenSolver<Matrix> solver(I);
+  Matrix evecs = solver.eigenvectors();
+  Matrix evals = solver.eigenvalues();
+```

Project#01/hints/hint7-1.md

@@ -0,0 +1 @@
+Lots of useful and precise physical constants are available at the [National Institute of Standards and Technology website](http://physics.nist.gov/cuu/Constants/index.html?/codata86.html).

Project#01/hints/hint7-2.md

@@ -0,0 +1,177 @@
+Now we've added a new member function to the Molecule class:
+```c++
+// Computes the angle between planes a-b-c and b-c-d
+double Molecule::torsion(int a, int b, int c, int d)
+{
+  double eabc_x = (unit(1,b,a)*unit(2,b,c) - unit(2,b,a)*unit(1,b,c));
+  double eabc_y = (unit(2,b,a)*unit(0,b,c) - unit(0,b,a)*unit(2,b,c));
+  double eabc_z = (unit(0,b,a)*unit(1,b,c) - unit(1,b,a)*unit(0,b,c));
+
+  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 = eabc_x * ebcd_x;
+  double eyy = eabc_y * ebcd_y;
+  double ezz = eabc_z * ebcd_z;
+
+  double tau = (exx + eyy + ezz)/(sin(angle(a,b,c)) * sin(angle(b,c,d)));
+
+  if(tau < -1.0) tau = acos(-1.0);
+  else if(tau > 1.0) tau = acos(1.0);
+  else tau = acos(tau);
+
+  // Compute the sign of the torsion 
+  double cross_x = eabc_y * ebcd_z - eabc_z * ebcd_y;
+  double cross_y = eabc_z * ebcd_x - eabc_x * ebcd_z;
+  double cross_z = eabc_x * ebcd_y - eabc_y * ebcd_x;
+  double norm = cross_x*cross_x + cross_y*cross_y + cross_z*cross_z;
+  cross_x /= norm;
+  cross_y /= norm;
+  cross_z /= norm;
+  double sign = 1.0;
+  double dot = cross_x*unit(0,b,c)+cross_y*unit(1,b,c)+cross_z*unit(2,b,c);
+  if(dot < 0.0) sign = -1.0;
+
+  return tau*sign;
+}
+```
+
+And we use the new function in the code as follows:
+
+```c++
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include <cstdio>
+#include <cmath>
+#include "molecule.h"
+
+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)));
+        }
+      }
+    }
+  }
+
+  cout << "\nTorsional angles:\n\n";
+  for(int i=0; i < mol.natom; i++) {
+    for(int j=0; j < i; j++) {
+      for(int k=0; k < j; k++) {
+        for(int l=0; l < k; l++) {
+          if(mol.bond(i,j) < 4.0 && mol.bond(j,k) < 4.0 && mol.bond(k,l) < 4.0)
+            printf("%2d-%2d-%2d-%2d %10.6f\n", i, j, k, l, mol.torsion(i,j,k,l)*(180.0/acos(-1.0)));
+        }
+      }
+    }
+  }
+
+  return 0;
+}
+```
+
+The above code produces the following output when applied to 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:
+ 0- 1- 2 124.268308
+ 0- 1- 3 115.479341
+ 0- 4- 5  35.109529
+ 0- 4- 6  35.109529
+ 0- 5- 6  36.373677
+ 1- 2- 3  28.377448
+ 4- 5- 6  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
+
+Torsional angles:
+
+ 3- 2- 1- 0 180.000000
+ 6- 5- 4- 0  36.366799
+```