--- trunk/src/applications/hydrodynamics/Hydro.cpp 2006/02/23 23:16:43 892 +++ 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 @@ -47,22 +48,26 @@ #include "applications/hydrodynamics/HydrodynamicsModel.hpp" #include "applications/hydrodynamics/HydrodynamicsModelCreator.hpp" #include "applications/hydrodynamics/HydrodynamicsModelFactory.hpp" +#include "applications/hydrodynamics/AnalyticalModel.hpp" #include "applications/hydrodynamics/BeadModel.hpp" #include "applications/hydrodynamics/RoughShell.hpp" +#include "applications/hydrodynamics/ShapeBuilder.hpp" #include "brains/Register.hpp" #include "brains/SimCreator.hpp" #include "brains/SimInfo.hpp" +#include "utils/StringUtils.hpp" +#include "utils/simError.h" +#include "utils/MemoryUtils.hpp" +using namespace OpenMD; -using namespace oopse; - -/** Register different hydrodynamics models */ +struct SDShape{ + StuntDouble* sd; + Shape* shape; +}; void registerHydrodynamicsModels(); +void writeHydroProps(std::ostream& os); -bool calcHydrodynamicsProp(const std::string& modelType, StuntDouble* sd, const DynamicProperty& param, std::ostream& os, const std::string& prefix); - int main(int argc, char* argv[]){ - //register force fields - registerForceFields(); registerHydrodynamicsModels(); gengetopt_args_info args_info; @@ -83,92 +88,116 @@ 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 = "hydro"; + prefix = getPrefix(dumpFileName); } std::string outputFilename = prefix + ".diff"; - - DynamicProperty param; - param.insert(DynamicProperty::value_type("Viscosity", args_info.viscosity_arg)); - param.insert(DynamicProperty::value_type("Temperature", args_info.temperature_arg)); - - if (args_info.sigma_given) { - param.insert(DynamicProperty::value_type("Sigma", args_info.sigma_arg)); - } - - + //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; - - - 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()) { - uniqueStuntDoubles.insert(std::map::value_type(integrableObject->getType(), integrableObject)); - integrableObject->setPos(V3Zero); - integrableObject->setA(identMat); - if (integrableObject->isRigidBody()) { - RigidBody* rb = static_cast(integrableObject); - rb->updateAtoms(); - } - } - } + + Globals* simParams = info->getSimParams(); + RealType temperature(0.0); + RealType viscosity(0.0); + + if (simParams->haveViscosity()) { + viscosity = simParams->getViscosity(); + } else { + sprintf(painCave.errMsg, "viscosity must be set\n"); + painCave.isFatal = 1; + simError(); } - std::map::iterator iter; - std::ofstream outputDiff(outputFilename.c_str()); - for (iter = uniqueStuntDoubles.begin(); iter != uniqueStuntDoubles.end(); ++iter) { - calcHydrodynamicsProp(args_info.model_arg, iter->second, param, outputDiff, prefix); + if (simParams->haveTargetTemp()) { + temperature = simParams->getTargetTemp(); + } else { + 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 (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::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); + } + + 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; + } + + + //MemoryUtils::deletePointers(shapes); delete info; } void registerHydrodynamicsModels() { - HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder("RoughShell")); - HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder("BeadModel")); - + HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder("RoughShell")); + HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder("BeadModel")); + HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder("AnalyticalModel")); } - -bool calcHydrodynamicsProp(const std::string& modelType, StuntDouble* sd, const DynamicProperty& param, std::ostream& os, const std::string& prefix) { - HydrodynamicsModel* hydroModel = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(modelType, sd, param); - bool ret = false; - if (hydroModel == NULL) { - std::cout << "Integrator Factory can not create " << modelType <calcHydrodyanmicsProps()) { - ret = true; - hydroModel->writeDiffCenterAndDiffTensor(os); - - std::ofstream ofs; - std::stringstream outputBeads; - outputBeads << prefix << "_" << sd->getType() << ".xyz"; - ofs.open(outputBeads.str().c_str()); - hydroModel->writeBeads(ofs); - ofs.close(); - } - - delete hydroModel; - - return ret; -}