35 |
|
* |
36 |
|
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
37 |
|
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
38 |
< |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). |
39 |
< |
* [4] Vardeman & Gezelter, in progress (2009). |
38 |
> |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
39 |
> |
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 |
> |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 |
|
*/ |
42 |
|
|
43 |
|
#include <iostream> |
61 |
|
using namespace OpenMD; |
62 |
|
|
63 |
|
struct SDShape{ |
64 |
< |
StuntDouble* sd; |
65 |
< |
Shape* shape; |
64 |
> |
StuntDouble* sd; |
65 |
> |
Shape* shape; |
66 |
|
}; |
67 |
|
void registerHydrodynamicsModels(); |
68 |
|
void writeHydroProps(std::ostream& os); |
69 |
|
|
70 |
|
int main(int argc, char* argv[]){ |
70 |
– |
//register force fields |
71 |
– |
registerForceFields(); |
71 |
|
registerHydrodynamicsModels(); |
72 |
|
|
73 |
|
gengetopt_args_info args_info; |
102 |
|
SimInfo::MoleculeIterator mi; |
103 |
|
Molecule* mol; |
104 |
|
Molecule::IntegrableObjectIterator ii; |
105 |
< |
StuntDouble* integrableObject; |
105 |
> |
StuntDouble* sd; |
106 |
|
Mat3x3d identMat; |
107 |
|
identMat(0,0) = 1.0; |
108 |
|
identMat(1,1) = 1.0; |
130 |
|
|
131 |
|
std::map<std::string, SDShape> uniqueStuntDoubles; |
132 |
|
|
133 |
< |
for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { |
134 |
< |
for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; |
135 |
< |
integrableObject = mol->nextIntegrableObject(ii)) { |
136 |
< |
if (uniqueStuntDoubles.find(integrableObject->getType()) == uniqueStuntDoubles.end()) { |
137 |
< |
|
138 |
< |
integrableObject->setPos(V3Zero); |
139 |
< |
integrableObject->setA(identMat); |
140 |
< |
if (integrableObject->isRigidBody()) { |
141 |
< |
RigidBody* rb = static_cast<RigidBody*>(integrableObject); |
142 |
< |
rb->updateAtoms(); |
143 |
< |
} |
144 |
< |
|
145 |
< |
SDShape tmp; |
146 |
< |
tmp.shape = ShapeBuilder::createShape(integrableObject); |
147 |
< |
tmp.sd = integrableObject; |
148 |
< |
uniqueStuntDoubles.insert(std::map<std::string, SDShape>::value_type(integrableObject->getType(), tmp)); |
149 |
< |
|
150 |
< |
} |
151 |
< |
} |
133 |
> |
for (mol = info->beginMolecule(mi); mol != NULL; |
134 |
> |
mol = info->nextMolecule(mi)) { |
135 |
> |
|
136 |
> |
for (sd = mol->beginIntegrableObject(ii); sd != NULL; |
137 |
> |
sd = mol->nextIntegrableObject(ii)) { |
138 |
> |
|
139 |
> |
if (uniqueStuntDoubles.find(sd->getType()) == uniqueStuntDoubles.end()) { |
140 |
> |
|
141 |
> |
sd->setPos(V3Zero); |
142 |
> |
sd->setA(identMat); |
143 |
> |
if (sd->isRigidBody()) { |
144 |
> |
RigidBody* rb = static_cast<RigidBody*>(sd); |
145 |
> |
rb->updateAtoms(); |
146 |
> |
} |
147 |
> |
|
148 |
> |
SDShape tmp; |
149 |
> |
tmp.shape = ShapeBuilder::createShape(sd); |
150 |
> |
tmp.sd = sd; |
151 |
> |
uniqueStuntDoubles.insert(std::map<std::string, SDShape>::value_type(sd->getType(), tmp)); |
152 |
> |
|
153 |
> |
} |
154 |
> |
} |
155 |
|
} |
154 |
– |
|
155 |
– |
|
156 |
|
|
157 |
+ |
|
158 |
+ |
|
159 |
|
std::map<std::string, SDShape>::iterator si; |
160 |
|
for (si = uniqueStuntDoubles.begin(); si != uniqueStuntDoubles.end(); ++si) { |
161 |
< |
HydrodynamicsModel* model; |
162 |
< |
Shape* shape = si->second.shape; |
163 |
< |
StuntDouble* sd = si->second.sd;; |
164 |
< |
if (args_info.model_given) { |
165 |
< |
model = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(args_info.model_arg, sd, info); |
166 |
< |
} else if (shape->hasAnalyticalSolution()) { |
167 |
< |
model = new AnalyticalModel(sd, info); |
168 |
< |
} else { |
169 |
< |
model = new BeadModel(sd, info); |
170 |
< |
} |
171 |
< |
|
172 |
< |
model->init(); |
161 |
> |
HydrodynamicsModel* model; |
162 |
> |
Shape* shape = si->second.shape; |
163 |
> |
StuntDouble* sd = si->second.sd;; |
164 |
> |
if (args_info.model_given) { |
165 |
> |
model = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(args_info.model_arg, sd, info); |
166 |
> |
} else if (shape->hasAnalyticalSolution()) { |
167 |
> |
model = new AnalyticalModel(sd, info); |
168 |
> |
} else { |
169 |
> |
model = new BeadModel(sd, info); |
170 |
> |
} |
171 |
> |
|
172 |
> |
model->init(); |
173 |
> |
|
174 |
> |
std::ofstream ofs; |
175 |
> |
std::stringstream outputBeads; |
176 |
> |
outputBeads << prefix << "_" << model->getStuntDoubleName() << ".xyz"; |
177 |
> |
ofs.open(outputBeads.str().c_str()); |
178 |
> |
model->writeBeads(ofs); |
179 |
> |
ofs.close(); |
180 |
> |
|
181 |
> |
//if beads option is turned on, skip the calculation |
182 |
> |
if (!args_info.beads_flag) { |
183 |
> |
model->calcHydroProps(shape, viscosity, temperature); |
184 |
> |
std::ofstream outputDiff; |
185 |
> |
outputDiff.open(outputFilename.c_str()); |
186 |
> |
model->writeHydroProps(outputDiff); |
187 |
> |
outputDiff.close(); |
188 |
> |
} |
189 |
|
|
190 |
< |
std::ofstream ofs; |
173 |
< |
std::stringstream outputBeads; |
174 |
< |
outputBeads << prefix << "_" << model->getStuntDoubleName() << ".xyz"; |
175 |
< |
ofs.open(outputBeads.str().c_str()); |
176 |
< |
model->writeBeads(ofs); |
177 |
< |
ofs.close(); |
178 |
< |
|
179 |
< |
//if beads option is turned on, skip the calculation |
180 |
< |
if (!args_info.beads_flag) { |
181 |
< |
model->calcHydroProps(shape, viscosity, temperature); |
182 |
< |
std::ofstream outputDiff; |
183 |
< |
outputDiff.open(outputFilename.c_str()); |
184 |
< |
model->writeHydroProps(outputDiff); |
185 |
< |
outputDiff.close(); |
186 |
< |
} |
187 |
< |
|
188 |
< |
delete model; |
190 |
> |
delete model; |
191 |
|
} |
192 |
|
|
193 |
|
|
197 |
|
} |
198 |
|
|
199 |
|
void registerHydrodynamicsModels() { |
200 |
< |
HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<RoughShell>("RoughShell")); |
201 |
< |
HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<BeadModel>("BeadModel")); |
202 |
< |
HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<AnalyticalModel>("AnalyticalModel")); |
201 |
< |
|
200 |
> |
HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<RoughShell>("RoughShell")); |
201 |
> |
HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<BeadModel>("BeadModel")); |
202 |
> |
HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<AnalyticalModel>("AnalyticalModel")); |
203 |
|
} |