--- trunk/src/applications/hydrodynamics/Hydro.cpp 2006/03/17 23:20:35 906 +++ trunk/src/applications/hydrodynamics/Hydro.cpp 2015/03/07 21:41:51 2071 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,6 +28,16 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [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, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include @@ -57,19 +58,17 @@ #include "utils/StringUtils.hpp" #include "utils/simError.h" #include "utils/MemoryUtils.hpp" -using namespace oopse; +using namespace OpenMD; struct SDShape{ - StuntDouble* sd; - Shape* shape; + StuntDouble* sd; + Shape* shape; }; void registerHydrodynamicsModels(); -void calcHydrodynamicsProp(HydrodynamicsModel* model, Shape* shape,double viscosity, double temperature, std::ostream& os, const std::string& prefix); +void writeHydroProps(std::ostream& os); int main(int argc, char* argv[]){ - //register force fields - registerForceFields(); - + registerHydrodynamicsModels(); gengetopt_args_info args_info; std::string dumpFileName; @@ -89,37 +88,34 @@ int main(int argc, char* argv[]){ exit(1); } - mdFileName = dumpFileName; - mdFileName = mdFileName.substr(0, mdFileName.rfind(".")) + ".md"; - if (args_info.output_given){ prefix = args_info.output_arg; } else { - prefix = getPrefix(mdFileName); + prefix = getPrefix(dumpFileName); } std::string outputFilename = prefix + ".diff"; //parse md file and set up the system SimCreator creator; - SimInfo* info = creator.createSim(mdFileName, true); + SimInfo* info = creator.createSim(dumpFileName, true); SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::IntegrableObjectIterator ii; - StuntDouble* integrableObject; + StuntDouble* sd; Mat3x3d identMat; identMat(0,0) = 1.0; identMat(1,1) = 1.0; identMat(2,2) = 1.0; Globals* simParams = info->getSimParams(); - double temperature; - double viscosity; + RealType temperature(0.0); + RealType viscosity(0.0); if (simParams->haveViscosity()) { viscosity = simParams->getViscosity(); } else { - sprintf(painCave.errMsg, "target temperature must be set\n"); + sprintf(painCave.errMsg, "viscosity must be set\n"); painCave.isFatal = 1; simError(); } @@ -127,49 +123,71 @@ int main(int argc, char* argv[]){ if (simParams->haveTargetTemp()) { temperature = simParams->getTargetTemp(); } else { - sprintf(painCave.errMsg, "viscosity must be set\n"); + sprintf(painCave.errMsg, "target temperature must be set\n"); painCave.isFatal = 1; simError(); } std::map 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()) { - - SDShape tmp; - tmp.shape = ShapeBuilder::createShape(integrableObject); - tmp.sd = integrableObject; - uniqueStuntDoubles.insert(std::map::value_type(integrableObject->getType(), tmp)); - integrableObject->setPos(V3Zero); - integrableObject->setA(identMat); - if (integrableObject->isRigidBody()) { - RigidBody* rb = static_cast(integrableObject); - rb->updateAtoms(); - } - } - } + 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(sd); + rb->updateAtoms(); + } + + SDShape tmp; + tmp.shape = ShapeBuilder::createShape(sd); + tmp.sd = sd; + uniqueStuntDoubles.insert(std::map::value_type(sd->getType(), tmp)); + + } + } } - - - std::ofstream outputDiff(outputFilename.c_str()); + + std::map::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); - } - calcHydrodynamicsProp(model, shape, viscosity, temperature, outputDiff, prefix); - delete model; + 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(); + + 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; } @@ -179,19 +197,7 @@ void registerHydrodynamicsModels() { } void registerHydrodynamicsModels() { - HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder("RoughShell")); - HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder("BeadModel")); - HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder("AnalyticalModel")); - + HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder("RoughShell")); + HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder("BeadModel")); + HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder("AnalyticalModel")); } -void calcHydrodynamicsProp(HydrodynamicsModel* model, Shape* shape,double viscosity, double temperature, std::ostream& os, const std::string& prefix) { - - shape->calcHydroProps(model, viscosity, temperature); - model->writeHydroProps(os); - std::ofstream ofs; - std::stringstream outputBeads; - outputBeads << prefix << "_" << model->getStuntDoubleName() << ".xyz"; - ofs.open(outputBeads.str().c_str()); - model->writeBeads(ofs); - ofs.close(); -}