OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
Hydro.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include <fstream>
46#include <iostream>
47#include <string>
48
49#include "HydroCmd.hpp"
50#include "brains/Register.hpp"
51#include "brains/SimCreator.hpp"
52#include "brains/SimInfo.hpp"
53#include "hydrodynamics/AnalyticalModel.hpp"
54#include "hydrodynamics/AtomicBeadModel.hpp"
55#include "hydrodynamics/BoundaryElementModel.hpp"
56#include "hydrodynamics/CompositeShape.hpp"
57#include "hydrodynamics/HydroIO.hpp"
58#include "hydrodynamics/HydrodynamicsModel.hpp"
59#include "hydrodynamics/HydrodynamicsModelCreator.hpp"
61#include "hydrodynamics/Mesh.hpp"
62#include "hydrodynamics/RoughShell.hpp"
63#include "hydrodynamics/ShapeBuilder.hpp"
64#include "hydrodynamics/Sphere.hpp"
65#include "io/MSMSFormat.hpp"
66#include "io/XYZFormat.hpp"
67#include "stl_reader.h"
69#include "utils/MemoryUtils.hpp"
71#include "utils/StringUtils.hpp"
72#include "utils/Trim.hpp"
73#include "utils/simError.h"
74
75using namespace OpenMD;
76
77void registerHydrodynamicsModels();
78void writeHydroProps(std::ostream& os);
79
80int main(int argc, char* argv[]) {
81 registerHydrodynamicsModels();
82
83 gengetopt_args_info args_info;
84 // parse the command line option
85 if (cmdline_parser(argc, argv, &args_info) != 0) { exit(1); }
86
87 std::string inFileName;
88 std::string modelName;
89 std::string modelSpecified;
90 std::string modelRequired;
91
92 bool hasInput = false;
93 bool hasModel = false;
94
95 Shape* shape;
96
97 // Check to make sure the model has been given
98 if (args_info.model_given) {
99 switch (args_info.model_arg) {
100 case model_arg_RoughShell:
101 modelSpecified = "RoughShell";
102 break;
103 case model_arg_BoundaryElement:
104 modelSpecified = "BoundaryElementModel";
105 break;
106 case model_arg_AtomicBead:
107 default:
108 modelSpecified = "AtomicBeadModel";
109 break;
110 }
111 hasModel = true;
112 }
113
114 // figure out what kind of input file we have
115
116 if (args_info.input_given) {
117 inFileName = args_info.input_arg;
118 modelRequired = modelSpecified;
119 hasInput = true;
120 } else if (args_info.stl_given) {
121 inFileName = args_info.stl_arg;
122 modelRequired = "BoundaryElementModel";
123 hasInput = true;
124 } else if (args_info.msms_given) {
125 inFileName = args_info.msms_arg;
126 modelRequired = "BoundaryElementModel";
127 hasInput = true;
128 } else if (args_info.xyz_given) {
129 inFileName = args_info.xyz_arg;
130 modelRequired = "AtomicBeadModel";
131 hasInput = true;
132 }
133
134 if (!hasInput) {
135 strcpy(painCave.errMsg, "No input file name was specified.\n");
136 painCave.isFatal = 1;
137 simError();
138 }
139
140 if (hasModel) {
141 if (modelSpecified.compare(modelRequired) != 0) {
142 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
143 "Specified model (%s) does not match model required for input "
144 "type (%s).\n",
145 modelSpecified.c_str(), modelRequired.c_str());
146 painCave.isFatal = 1;
147 simError();
148 } else {
149 modelName = modelSpecified;
150 }
151 } else {
152 modelName = modelRequired;
153 hasModel = true;
154 }
155
156 std::string prefix;
157 if (args_info.output_given) {
158 prefix = args_info.output_arg;
159 } else {
160 prefix = getPrefix(inFileName);
161 }
162
163 std::string outputFilename = prefix + ".hydro";
164
165 RealType temperature = args_info.temperature_arg;
166 RealType viscosity = args_info.viscosity_arg;
167 RealType beadSize = args_info.beadSize_arg;
168
169 // read the files and create the shapes:
170
171 std::map<std::string, Shape*> uniqueShapes;
172
173 if (args_info.stl_given) {
174 try {
175 stl_reader::StlMesh<RealType, unsigned int> mesh(inFileName.c_str());
176 for (size_t isolid = 0; isolid < mesh.num_solids(); ++isolid) {
177 shape = new Mesh();
178 std::cout << "solid " << isolid << std::endl;
179 for (size_t itri = mesh.solid_tris_begin(isolid);
180 itri < mesh.solid_tris_end(isolid); ++itri) {
181 const RealType* c0 = mesh.tri_corner_coords(itri, 0);
182 const RealType* c1 = mesh.tri_corner_coords(itri, 1);
183 const RealType* c2 = mesh.tri_corner_coords(itri, 2);
184 Vector3d vert0(c0[0], c0[1], c0[2]);
185 Vector3d vert1(c1[0], c1[1], c1[2]);
186 Vector3d vert2(c2[0], c2[1], c2[2]);
187 dynamic_cast<Mesh*>(shape)->add(vert0, vert1, vert2);
188 }
189 std::string solidName;
190 if (mesh.num_solids() > 1) {
191 solidName = prefix + "_" + std::to_string(isolid);
192 } else {
193 solidName = prefix;
194 }
195
196 shape->setName(solidName);
197 uniqueShapes.insert(
198 std::map<std::string, Shape*>::value_type(solidName, shape));
199 }
200 } catch (std::exception& e) { std::cout << e.what() << std::endl; }
201
202 } else if (args_info.msms_given) {
203 MSMSFormat* msms = new MSMSFormat(inFileName.c_str());
204 shape = msms->ReadShape();
205 shape->setName(prefix);
206 uniqueShapes.insert(
207 std::map<std::string, Shape*>::value_type(prefix, shape));
208 } else if (args_info.xyz_given) {
209 ifstream in(inFileName);
210 if (!in) {
211 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
212 "Could not open XYZ file\n");
213 painCave.isFatal = 1;
214 simError();
215 }
216
217 XYZFormat* xyz = new XYZFormat();
218 xyz->ReadMolecule(in);
219
220 shape = new CompositeShape();
221 shape->setName(xyz->title_);
222
223 Shape* currShape = NULL;
224
225 size_t natoms = xyz->mol_.size();
226 for (size_t iatom = 0; iatom < natoms; ++iatom) {
227 Vector3d pos = xyz->mol_[iatom]->pos;
228 int anum = xyz->mol_[iatom]->atomicNum;
229 std::string atype = xyz->mol_[iatom]->type;
230 currShape = new Sphere(pos, etab.GetVdwRad(anum));
231 currShape->setName(atype);
232 if (currShape != NULL) {
233 dynamic_cast<CompositeShape*>(shape)->addShape(currShape);
234 }
235 }
236 uniqueShapes.insert(
237 std::map<std::string, Shape*>::value_type(xyz->title_, shape));
238 } else {
239 // parse md file and set up the system
240 SimCreator creator;
241 SimInfo* info = creator.createSim(inFileName, true);
242
243 SimInfo::MoleculeIterator mi;
244 Molecule* mol;
245 Molecule::IntegrableObjectIterator ii;
246 StuntDouble* sd;
247
248 Globals* simParams = info->getSimParams();
249
250 if (simParams->haveViscosity()) {
251 viscosity = simParams->getViscosity();
252 } else {
253 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
254 "viscosity must be set\n");
255 painCave.isFatal = 1;
256 simError();
257 }
258
259 if (simParams->haveTargetTemp()) {
260 temperature = simParams->getTargetTemp();
261 } else {
262 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
263 "target temperature must be set\n");
264 painCave.isFatal = 1;
265 simError();
266 }
267
268 for (mol = info->beginMolecule(mi); mol != NULL;
269 mol = info->nextMolecule(mi)) {
270 for (sd = mol->beginIntegrableObject(ii); sd != NULL;
271 sd = mol->nextIntegrableObject(ii)) {
272 if (uniqueShapes.find(sd->getType()) == uniqueShapes.end()) {
273 sd->setPos(V3Zero);
274 if (sd->isRigidBody()) {
276 RigidBody* rb = static_cast<RigidBody*>(sd);
277 rb->updateAtoms();
278 }
279
280 Shape* tmp = ShapeBuilder::createShape(sd);
281 uniqueShapes.insert(
282 std::map<std::string, Shape*>::value_type(sd->getType(), tmp));
283 }
284 }
285 }
286 delete info;
287 }
288
289 HydrodynamicsModel* model =
290 HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(
291 modelName);
292
293 if (model == NULL) {
294 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
295 "Could not create HydrodynamicsModel\n");
296 painCave.isFatal = 1;
297 simError();
298 } else {
299 model->init();
300
301 if (modelName.compare("RoughShell") == 0) {
302 dynamic_cast<RoughShell*>(model)->setSigma(beadSize);
303 }
304
305 std::ofstream outputHydro;
306 outputHydro.open(outputFilename.c_str());
307 HydroIO* hio = new HydroIO();
308 hio->openWriter(outputHydro);
309
310 std::map<std::string, Shape*>::iterator si;
311 for (si = uniqueShapes.begin(); si != uniqueShapes.end(); ++si) {
312 shape = si->second;
313 model->setShape(shape);
314
315 // write out the elements making up the shape:
316 std::ofstream ofs;
317 std::stringstream elementFile;
318 elementFile << prefix << "_" << shape->getName();
319 if (modelName.compare("BoundaryElementModel") == 0) {
320 elementFile << ".stl";
321 } else {
322 elementFile << ".xyz";
323 }
324 ofs.open(elementFile.str().c_str());
325 model->writeElements(ofs);
326 ofs.close();
327
328 // if elements option is turned on, skip the calculation
329 if (!args_info.elements_flag) {
330 HydroProp* hp = model->calcHydroProps(viscosity);
331 hio->writeHydroProp(hp, viscosity, temperature, outputHydro);
332 hio->interpretHydroProp(hp, viscosity, temperature);
333 }
334 }
335
336 hio->closeWriter(outputHydro);
337 outputHydro.close();
338
339 delete model;
340 }
341}
342
343void registerHydrodynamicsModels() {
344 HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(
345 new HydrodynamicsModelBuilder<RoughShell>("RoughShell"));
346 HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(
347 new HydrodynamicsModelBuilder<AtomicBeadModel>("AtomicBeadModel"));
348 HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(
349 new HydrodynamicsModelBuilder<AnalyticalModel>("AnalyticalModel"));
350 HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(
351 new HydrodynamicsModelBuilder<BoundaryElementModel>(
352 "BoundaryElementModel"));
353}
This basic Periodic Table class was originally taken from the data.h file in OpenBabel.
The header file for the command line option parser generated by GNU Gengetopt version 2....
Combine composite pattern and visitor pattern.
RealType GetVdwRad(int atomicnum)
Container for information about the hydrodynamic behavior of objects interacting with surroundings.
Definition HydroProp.hpp:74
static HydrodynamicsModelFactory * getInstance()
Returns an instance of HydrodynamicsModel factory.
void updateAtoms()
update the positions of atoms belong to this rigidbody
The only responsibility of SimCreator is to parse the meta-data file and create a SimInfo instance ba...
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Definition SimInfo.hpp:93
Molecule * beginMolecule(MoleculeIterator &i)
Returns the first molecule in this SimInfo and intialize the iterator.
Definition SimInfo.cpp:240
Molecule * nextMolecule(MoleculeIterator &i)
Returns the next avaliable Molecule based on the iterator.
Definition SimInfo.cpp:245
"Don't move, or you're dead! Stand up! Captain, we've got them!"
virtual void setA(const RotMat3x3d &a)
Sets the current rotation matrix of this stuntDouble.
virtual std::string getType()=0
Returns the name of this stuntDouble.
void setPos(const Vector3d &pos)
Sets the current position of this stuntDouble.
bool isRigidBody()
Tests if this stuntDouble is a rigid body.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
std::string getPrefix(const std::string &str)
Where the command line options are stored.