--- trunk/src/applications/hydrodynamics/Hydro.cpp	2012/08/21 13:51:14	1781
+++ trunk/src/applications/hydrodynamics/Hydro.cpp	2012/08/22 02:28:28	1782
@@ -36,7 +36,8 @@
  * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
  * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
  * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
- * [4]  Vardeman & Gezelter, in progress (2009).                        
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
  
 #include <iostream>
@@ -60,15 +61,13 @@ struct SDShape{
 using namespace OpenMD;
 
 struct SDShape{
-    StuntDouble* sd;
-    Shape* shape;
+  StuntDouble* sd;
+  Shape* shape;
 };
 void registerHydrodynamicsModels();
 void writeHydroProps(std::ostream& os);
 
 int main(int argc, char* argv[]){
-  //register force fields
-  registerForceFields();    
   registerHydrodynamicsModels();
   
   gengetopt_args_info args_info;
@@ -103,7 +102,7 @@ int main(int argc, char* argv[]){
   SimInfo::MoleculeIterator mi;
   Molecule* mol;
   Molecule::IntegrableObjectIterator  ii;
-  StuntDouble* integrableObject;
+  StuntDouble* sd;
   Mat3x3d identMat;
   identMat(0,0) = 1.0;
   identMat(1,1) = 1.0;
@@ -131,61 +130,64 @@ int main(int argc, char* argv[]){
  
   std::map<std::string, SDShape> uniqueStuntDoubles;
   
-  for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
-      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
-	   integrableObject = mol->nextIntegrableObject(ii)) {
-          if (uniqueStuntDoubles.find(integrableObject->getType()) ==  uniqueStuntDoubles.end()) {
-
-            integrableObject->setPos(V3Zero);
-            integrableObject->setA(identMat);
-            if (integrableObject->isRigidBody()) {
-                RigidBody* rb = static_cast<RigidBody*>(integrableObject);
-                rb->updateAtoms();
-            }
-
-            SDShape tmp;
-            tmp.shape = ShapeBuilder::createShape(integrableObject);
-            tmp.sd = integrableObject;    
-            uniqueStuntDoubles.insert(std::map<std::string, SDShape>::value_type(integrableObject->getType(), tmp));
-
-          }
-	}
+  for (mol = info->beginMolecule(mi); mol != NULL; 
+       mol = info->nextMolecule(mi)) {
+    
+    for (sd = mol->beginIntegrableObject(ii); sd != NULL;
+         sd = mol->nextIntegrableObject(ii)) {
+      
+      if (uniqueStuntDoubles.find(sd->getType()) ==  uniqueStuntDoubles.end()) {
+        
+        sd->setPos(V3Zero);
+        sd->setA(identMat);
+        if (sd->isRigidBody()) {
+          RigidBody* rb = static_cast<RigidBody*>(sd);
+          rb->updateAtoms();
+        }
+        
+        SDShape tmp;
+        tmp.shape = ShapeBuilder::createShape(sd);
+        tmp.sd = sd;    
+        uniqueStuntDoubles.insert(std::map<std::string, SDShape>::value_type(sd->getType(), tmp));
+        
+      }
+    }
   }
-
-
   
+  
+  
   std::map<std::string, SDShape>::iterator si;
   for (si = uniqueStuntDoubles.begin(); si != uniqueStuntDoubles.end(); ++si) {
-      HydrodynamicsModel* model;
-      Shape* shape = si->second.shape;
-      StuntDouble* sd = si->second.sd;;
-      if (args_info.model_given) {  
-        model = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(args_info.model_arg, sd, info);
-      } else if (shape->hasAnalyticalSolution()) {
-        model = new AnalyticalModel(sd, info);
-      } else {
-        model = new BeadModel(sd, info);
-      }
-
-        model->init();
+    HydrodynamicsModel* model;
+    Shape* shape = si->second.shape;
+    StuntDouble* sd = si->second.sd;;
+    if (args_info.model_given) {  
+      model = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(args_info.model_arg, sd, info);
+    } else if (shape->hasAnalyticalSolution()) {
+      model = new AnalyticalModel(sd, info);
+    } else {
+      model = new BeadModel(sd, info);
+    }
+    
+    model->init();
+    
+    std::ofstream ofs;
+    std::stringstream outputBeads;
+    outputBeads << prefix << "_" << model->getStuntDoubleName() << ".xyz";
+    ofs.open(outputBeads.str().c_str());        
+    model->writeBeads(ofs);
+    ofs.close();
+    
+    //if beads option is turned on, skip the calculation
+    if (!args_info.beads_flag) {
+      model->calcHydroProps(shape, viscosity, temperature);
+      std::ofstream outputDiff;
+      outputDiff.open(outputFilename.c_str());
+      model->writeHydroProps(outputDiff);
+      outputDiff.close();
+    }
         
-        std::ofstream ofs;
-        std::stringstream outputBeads;
-        outputBeads << prefix << "_" << model->getStuntDoubleName() << ".xyz";
-        ofs.open(outputBeads.str().c_str());        
-        model->writeBeads(ofs);
-        ofs.close();
-
-        //if beads option is turned on, skip the calculation
-        if (!args_info.beads_flag) {
-          model->calcHydroProps(shape, viscosity, temperature);
-	  std::ofstream outputDiff;
-	  outputDiff.open(outputFilename.c_str());
-          model->writeHydroProps(outputDiff);
-	  outputDiff.close();
-        }
-        
-        delete model;
+    delete model;
   }
 
 
@@ -195,8 +197,7 @@ void registerHydrodynamicsModels() {
 }
 
 void registerHydrodynamicsModels() {
-    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<RoughShell>("RoughShell"));
-    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<BeadModel>("BeadModel"));
-    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<AnalyticalModel>("AnalyticalModel"));
-
+  HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<RoughShell>("RoughShell"));
+  HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<BeadModel>("BeadModel"));
+  HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<AnalyticalModel>("AnalyticalModel"));
 }