fastlib/la/matrix.h

@@ -646,7 +646,7 @@ class Matrix {
    *
    * @param col the column number
    * @return an array where the i'th element is the i'th row of that
-   *         par ticular column
+   *         particular column
    */
   double *GetColumnPtr(index_t col) {
     DEBUG_BOUNDS(col, n_cols_);

fastlib/u/jim/md/AtomTree.h

@@ -1,14 +1,10 @@
 /**
  * KD-tree stucture for molecular dynamics simulation. 
  * Each node stores a bounding box, centroid, and number of
- * atoms. Velocity for dynamics problems can be stored in 
- * a separate matrix. 
- *
- * Eventually, we will have to add stats to permit multiple
- * types of atoms in a simulation.
+ * atoms. Leaf nodes also store the velocity of the corresponding 
+ * atom. Atoms are assumed to be homogeneous, as this is the
+ * largely the case for applications of the LJ potential.
  * 
- * J. Waters 
- * Begun 11-13-2007 
  */
 
 #include "fastlib/fastlib.h"
@@ -17,14 +13,14 @@
 #include "tree/bounds.h"
 
 
-// We need to track total number of atoms and centroid for each node.
 struct AtomStat {
   double mass;
   Vector centroid;
-  Vector velocity; // At present, only the velocities at leaf nodes
-                    // are used and calculated properly.
- 
-  // Basic Initialization
+  Vector velocity; 
+
+  /**
+   * Default Initialization
+   */
   void Init(){
     centroid.Init(3);
     velocity.Init(3);
@@ -33,7 +29,9 @@ struct AtomStat {
     mass = 0;
   }
 
-  // Leaf node initialization
+  /**
+   * Init funciton for leaf node. Each leaf corresponds to a single atom.
+   */
   void Init(const Matrix& dataset, int start, int count){
     centroid.Init(3);
     centroid.SetZero();
@@ -50,7 +48,11 @@ struct AtomStat {
   }
 
 
-  // Non-leaf node initialization
+  /**
+   * Init function to build node from two children, tracking mass and
+   * centroid of each node. Since the updating of velocities is done
+   * as a single tree search, non-leaves do not need to store velocity.
+   */
   void Init(const Matrix& dataset, int start, int count, 
 	    const AtomStat &left_stat, const AtomStat &right_stat){
     

fastlib/u/jim/md/LennardJones.h

@@ -1,10 +1,15 @@
 /**
+ * @file LennardJones.h
+ *
  * Molecular Dynamics via Lennard-Jones potential
  *
  * Accelerations are calculated via single-tree search. 
  * Equations of motion are integrated by a leapfrogging method.
+ *
+ * @see LennardJones_main.cc
  */
 
+
 #include "fastlib/fastlib.h"
 #include "fastlib/fastlib_int.h"
 #include "AtomTree.h"
@@ -34,6 +39,8 @@ private:
    * Compute the change in the velocity of center_1, according to Lennard-Jones
    * potential between centers 1 & 2. This force will be scaled when it is returned
    * to UpdateVelocityRecursion, if center_2 corresponds to more than one atom.
+   * Note also that this computes the acceleration integrated over the time step,
+   * equal to the change in velocity, rather than the instantaneous acceleration.
    */
   void Acceleration_(Vector &center_1, Vector &center_2, Vector &delta_v_){
     int i;
@@ -43,13 +50,13 @@ private:
     double dist_sq_ = la::DistanceSqEuclidean(center_1, center_2);
     double r_scaled_ = sig*sig / dist_sq_;
     r_scaled_ = r_scaled_*r_scaled_*r_scaled_;
-    double force_mag_ = time_step*24*eps*r_scaled_*(1 - 2*r_scaled_) / (dist_sq_*mass);
+    double force_mag_ = time_step*24*eps*r_scaled_*(2*r_scaled_ - 1) / (dist_sq_*mass);
     la::Scale(force_mag_, &delta_v_);         
   }
 
   /**
    * Update centroids and bounding boxes. Note that we may develop
-   * intersections between bounding boxes.
+   * intersections between bounding boxes as the simulation progresses.
    */
   void UpdatePositionsRecursion_(AtomTree *current_node){
     int i;     
@@ -63,7 +70,8 @@ private:
 	current_node->stat().centroid[i] = 
 	  current_node->left()->stat().mass*current_node->left()->stat().centroid[i] + 
 	  current_node->right()->stat().mass*current_node->right()->stat().centroid[i];
-	current_node->stat().centroid[i] = current_node->stat().centroid[i] / current_node->stat().mass;
+	current_node->stat().centroid[i] = 
+	  current_node->stat().centroid[i] / current_node->stat().mass;
       } 
     } else {  // Base Case
       for (i = 0; i < 3; i++){
@@ -75,8 +83,12 @@ private:
     }      
   }
 
-  void UpdateVelocityRecursion_(AtomTree* vel_query_, AtomTree* vel_ref_){
-    //    int i;
+  /**
+   * Compute the effect of the reference node on the velocity of
+   * the query node. The effect of distant atoms is approximated
+   * from the centroid of the reference node atoms.
+   */
+  void UpdateVelocityRecursion_(AtomTree* vel_query_, AtomTree* vel_ref_){   
     if (unlikely(vel_ref_->count() == 1)){
       UpdateVelocityBase_(vel_query_, vel_ref_);
     } else {
@@ -95,6 +107,10 @@ private:
     }
   } //UpdateVelocityRecursion
 
+
+  /**
+   * Base case calculates pairwise interactions between nearby atoms.
+   */
   void UpdateVelocityBase_(AtomTree* vel_query_, AtomTree* vel_ref_){
     if (likely(vel_query_->begin() != vel_ref_->begin())){
       Vector delta_v_;
@@ -102,7 +118,7 @@ private:
       Acceleration_(vel_query_->stat().centroid, vel_ref_->stat().centroid, delta_v_);
       la::AddTo(delta_v_, &vel_query_->stat().velocity);
     }
-  }
+  } // UpdateVelocityBase
 
 
   ////////////////////////////// Constructors ///////////////////////////////////////
@@ -138,6 +154,10 @@ public:
   } //Init
 
 
+  /**
+   * Naive implementation computes all pairwise interactions, and can be used to
+   * validate approximations made by tree implementation.
+   */
   void InitNaive(const Matrix& atoms_in, double eps_in, double sig_in, double mass_in){
 
     atoms_.Copy(atoms_in);
@@ -149,7 +169,7 @@ public:
     mass = mass_in;
     velocities_.Init(3, n_atoms_);
     velocities_.SetZero();
-  }
+  } // InitNaive
 
 
   void UpdatePositions(double time_step_in){
@@ -219,7 +239,7 @@ public:
       }
     }
 
-  }
+  } // WritePositions
 
 }; // class LennardJones
 

fastlib/u/jim/md/LennardJones_main.cc

@@ -1,5 +1,11 @@
 /**
  * @file LennardJones_main.cc
+ * 
+ * This program creates an instance of the LennardJones problem
+ * class, and updates the velocities using a leapfrogging scheme
+ * until a specified end time is reached. 
+ *
+ * @see LennardJones.h
  */
 
 
@@ -7,7 +13,6 @@
 
 int main(int argc, char *argv[])
 {
-  // Parse Inputs
   fx_init(argc, argv);
   const char* fp;
   fp = fx_param_str(NULL, "data", "default.txt");
@@ -37,15 +42,16 @@ int main(int argc, char *argv[])
   // Read Atom Matrix
   data::Load(fp, &atom_mat);
 
+
+  /**
+   * Begin simulation, and run to end time.
+   */
   simulation.Init(atom_mat, eps, sig, mass, cutoff);
   simulation.UpdateVelocities(time_step/2);
   
   naive_test.InitNaive(atom_mat, eps, sig, mass);
   naive_test.UpdateVelocitiesNaive(time_step/2);
   
-  
-  
-
   while (time < stop_time){
     simulation.UpdateVelocities(time_step);
     simulation.UpdatePositions(time_step);
@@ -57,7 +63,10 @@ int main(int argc, char *argv[])
   
   }
 
-simulation.WritePositions(tree);
+  /**
+   * Record final positions according to both methods.
+   */
+  simulation.WritePositions(tree);
   naive_test.WritePositions(naive);
 
   fx_done();