45#include "hydrodynamics/HydroIO.hpp"
50#include "utils/simError.h"
56 HydroIO::~HydroIO() {}
58 void HydroIO::openWriter(std::ostream& os) {
59 std::string h =
"OpenMD-Hydro";
60#if defined(NLOHMANN_JSON)
61 j_[h] = ordered_json::array();
63#elif defined(RAPID_JSON)
64 osw_ =
new OStreamWrapper(os);
68 w_.SetMaxDecimalPlaces(7);
77 void HydroIO::writeHydroProp(HydroProp* hp, RealType viscosity,
78 RealType temperature, std::ostream& os) {
79 if (!writerOpen_) openWriter(os);
81 std::string h =
"OpenMD-Hydro";
82 std::string name = hp->getName();
83 Vector3d cor = hp->getCenterOfResistance();
84 Mat6x6d Xi = hp->getResistanceTensor();
85 Vector3d cod = hp->getCenterOfDiffusion(temperature);
86 Mat6x6d Xid = hp->getDiffusionTensorAtPos(cod, temperature);
87 Vector3d cop = hp->getCenterOfPitch();
92 hp->pitchAxes(pitchAxes, pitches, pitchScalar);
94#if defined(NLOHMANN_JSON)
98 o[
"viscosity"] = viscosity;
99 o[
"centerOfResistance"] = {cor[0], cor[1], cor[2]};
100 o[
"resistanceTensor"] = json::array();
102 for (
unsigned int i = 0; i < 6; i++) {
103 o[
"resistanceTensor"][i] = {Xi(i, 0), Xi(i, 1), Xi(i, 2),
104 Xi(i, 3), Xi(i, 4), Xi(i, 5)};
107 o[
"temperature"] = temperature;
108 o[
"centerOfDiffusion"] = {cod[0], cod[1], cod[2]};
109 o[
"diffusionTensor"] = json::array();
111 for (
unsigned int i = 0; i < 6; i++) {
112 o[
"diffusionTensor"][i] = {Xid(i, 0), Xid(i, 1), Xid(i, 2),
113 Xid(i, 3), Xid(i, 4), Xid(i, 5)};
116 o[
"pitch"] = pitchScalar;
117 o[
"centerOfPitch"] = {cop[0], cop[1], cop[2]};
118 o[
"pitches"] = {pitches[0], pitches[1], pitches[2]};
120 o[
"pitchAxes"] = json::array();
121 for (
unsigned int i = 0; i < 3; i++) {
122 o[
"pitchAxes"][i] = {pitchAxes(i, 0), pitchAxes(i, 1), pitchAxes(i, 2)};
127#elif defined(RAPID_JSON)
131 w_.String(name.c_str());
134 w_.Double(viscosity);
135 w_.Key(
"centerOfResistance");
137 w_.SetFormatOptions(kFormatSingleLineArray);
139 for (
unsigned i = 0; i < 3; i++)
142 w_.SetFormatOptions(kFormatDefault);
144 w_.Key(
"resistanceTensor");
146 for (
unsigned i = 0; i < 6; i++) {
148 w_.SetFormatOptions(kFormatSingleLineArray);
150 for (
unsigned j = 0; j < 6; j++) {
154 w_.SetFormatOptions(kFormatDefault);
158 w_.Key(
"temperature");
159 w_.Double(temperature);
160 w_.Key(
"centerOfDiffusion");
162 w_.SetFormatOptions(kFormatSingleLineArray);
164 for (
unsigned i = 0; i < 3; i++)
167 w_.SetFormatOptions(kFormatDefault);
169 w_.Key(
"diffusionTensor");
171 for (
unsigned i = 0; i < 6; i++) {
173 w_.SetFormatOptions(kFormatSingleLineArray);
175 for (
unsigned j = 0; j < 6; j++) {
176 w_.Double(Xid(i, j));
179 w_.SetFormatOptions(kFormatDefault);
184 w_.Double(pitchScalar);
185 w_.Key(
"centerOfPitch");
187 w_.SetFormatOptions(kFormatSingleLineArray);
189 for (
unsigned i = 0; i < 3; i++)
192 w_.SetFormatOptions(kFormatDefault);
196 for (
unsigned i = 0; i < 3; i++) {
198 w_.SetFormatOptions(kFormatSingleLineArray);
199 for (
unsigned j = 0; j < 3; j++) {
200 w_.Double(pitchAxes(i, j));
203 w_.SetFormatOptions(kFormatDefault);
211 void HydroIO::closeWriter(std::ostream& os) {
212#if defined(NLOHMANN_JSON)
213 os << j_.dump(2) << std::endl;
214#elif defined(RAPID_JSON)
222 map<string, HydroProp*> HydroIO::parseHydroFile(
const string& f) {
223 map<string, HydroProp*> props;
228 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
229 "HydroIO: Cannot open file: %s\n", f.c_str());
230 painCave.isFatal = 1;
234#if defined(NLOHMANN_JSON)
235 json ij = json::parse(ifs);
237 auto& entries = ij[
"OpenMD-Hydro"];
239 for (
auto& entry : entries) {
240 HydroProp* hp =
new HydroProp();
245 name = entry[
"name"].get<std::string>();
247 for (
unsigned int i = 0; i < 3; i++) {
248 cor[i] = entry[
"centerOfResistance"].get<vector<RealType>>()[i];
251 for (
unsigned int i = 0; i < 6; i++) {
252 for (
unsigned int j = 0; j < 6; j++) {
254 entry[
"resistanceTensor"].get<vector<vector<RealType>>>()[i][j];
259 hp->setCenterOfResistance(cor);
260 hp->setResistanceTensor(Xi);
261 props.insert(map<string, HydroProp*>::value_type(name, hp));
263#elif defined(RAPID_JSON)
267 if (ifs.peek() != EOF) {
268 rapidjson::IStreamWrapper isw(ifs);
270 if (d_.HasParseError()) {
271 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
272 "HydroIO: JSON parse error in file %s\n"
273 "\tError: %zu : %s\n",
274 f.c_str(), d_.GetErrorOffset(),
275 rapidjson::GetParseError_En(d_.GetParseError()));
276 painCave.isFatal = 1;
279 if (!d_.IsObject()) {
280 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
281 "HydroIO: OpenMD-Hydro should be a single object.\n");
282 painCave.isFatal = 1;
287 if (!d_.HasMember(
"OpenMD-Hydro")) {
288 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
289 "HydroIO: File %s does not have a OpenMD-Hydro object.\n",
291 painCave.isFatal = 1;
295 const Value& entries = d_[
"OpenMD-Hydro"];
296 for (
auto& entry : entries.GetArray()) {
297 HydroProp* hp =
new HydroProp();
302 name = entry[
"name"].GetString();
304 for (
unsigned int i = 0; i < 3; i++) {
305 cor[i] = entry[
"centerOfResistance"][i].GetDouble();
308 for (
unsigned int i = 0; i < 6; i++) {
309 for (
unsigned int j = 0; j < 6; j++) {
310 Xi(i, j) = entry[
"resistanceTensor"][i][j].GetDouble();
315 hp->setCenterOfResistance(cor);
316 hp->setResistanceTensor(Xi);
317 props.insert(map<string, HydroProp*>::value_type(name, hp));
323 void HydroIO::interpretHydroProp(HydroProp* hp, RealType viscosity,
324 RealType temperature) {
325 Vector3d ror = hp->getCenterOfResistance();
328 Xi = hp->getResistanceTensor();
334 Xi.getSubMatrix(0, 0, Xirtt);
335 Xi.getSubMatrix(0, 3, Xirrt);
336 Xi.getSubMatrix(3, 0, Xirtr);
337 Xi.getSubMatrix(3, 3, Xirrr);
340 Dr = hp->getDiffusionTensor(temperature);
347 Dr.getSubMatrix(0, 0, Drtt);
348 Dr.getSubMatrix(0, 3, Drrt);
349 Dr.getSubMatrix(3, 0, Drtr);
350 Dr.getSubMatrix(3, 3, Drrr);
353 std::cout <<
"-----------------------------------------\n";
354 std::cout <<
"viscosity = " << viscosity <<
" Poise" << std::endl;
355 std::cout <<
"temperature = " << temperature <<
" K" << std::endl;
356 std::cout <<
"-----------------------------------------\n";
357 std::cout <<
"The centers are based on the elements generated by Hydro "
359 std::cout <<
"which have been placed in an .xyz or .stl file." << std::endl;
360 std::cout <<
"They are not based on the geometry in the .omd file.\n"
362 std::cout <<
"-----------------------------------------\n\n";
363 std::cout <<
"Center of resistance :" << std::endl;
364 std::cout << ror <<
"\n" << std::endl;
365 std::cout <<
"-----------------------------------------\n\n";
366 std::cout <<
"Resistance tensor at center of resistance\n" << std::endl;
367 std::cout <<
"translation [kcal.fs/(mol.A^2)]:" << std::endl;
368 std::cout << Xirtt << std::endl;
369 std::cout <<
"rotation-translation [kcal.fs/(mol.A.radian)]:" << std::endl;
370 std::cout << Xirtr.transpose() << std::endl;
371 std::cout <<
"translation-rotation [kcal.fs/(mol.A.radian)]:" << std::endl;
372 std::cout << Xirtr << std::endl;
373 std::cout <<
"rotation [kcal.fs/(mol.radian^2)]:" << std::endl;
374 std::cout << Xirrr << std::endl;
375 std::cout <<
"-----------------------------------------\n\n";
376 std::cout <<
"Diffusion tensor at center of resistance\n" << std::endl;
377 std::cout <<
"translation (A^2 / fs):" << std::endl;
378 std::cout << Drtt << std::endl;
379 std::cout <<
"rotation-translation (A.radian / fs):" << std::endl;
380 std::cout << Drrt << std::endl;
381 std::cout <<
"translation-rotation (A.radian / fs):" << std::endl;
382 std::cout << Drtr << std::endl;
383 std::cout <<
"rotation (radian^2 / fs):" << std::endl;
384 std::cout << Drrr << std::endl;
385 std::cout <<
"-----------------------------------------\n\n";
390 Vector3d cod = hp->getCenterOfDiffusion(temperature);
391 Mat6x6d Xid = hp->getResistanceTensorAtPos(cod);
398 Xid.getSubMatrix(0, 0, Xidtt);
399 Xid.getSubMatrix(0, 3, Xidrt);
400 Xid.getSubMatrix(3, 0, Xidtr);
401 Xid.getSubMatrix(3, 3, Xidrr);
404 Mat6x6d Dd = hp->getDiffusionTensorAtPos(cod, temperature);
411 Dd.getSubMatrix(0, 0, Ddtt);
412 Dd.getSubMatrix(0, 3, Ddrt);
413 Dd.getSubMatrix(3, 0, Ddtr);
414 Dd.getSubMatrix(3, 3, Ddrr);
416 std::cout <<
"Center of diffusion: " << std::endl;
417 std::cout << cod <<
"\n" << std::endl;
418 std::cout <<
"-----------------------------------------\n\n";
419 std::cout <<
"Diffusion tensor at center of diffusion \n " << std::endl;
420 std::cout <<
"translation (A^2 / fs) :" << std::endl;
421 std::cout << Ddtt << std::endl;
422 std::cout <<
"rotation-translation (A.radian / fs):" << std::endl;
423 std::cout << Ddtr.transpose() << std::endl;
424 std::cout <<
"translation-rotation (A.radian / fs):" << std::endl;
425 std::cout << Ddtr << std::endl;
426 std::cout <<
"rotation (radian^2 / fs):" << std::endl;
427 std::cout << Ddrr << std::endl;
428 std::cout <<
"-----------------------------------------\n\n";
429 std::cout <<
"Resistance tensor at center of diffusion \n " << std::endl;
430 std::cout <<
"translation [kcal.fs/(mol.A^2)]:" << std::endl;
431 std::cout << Xidtt << std::endl;
432 std::cout <<
"rotation-translation [kcal.fs/(mol.A.radian)]:" << std::endl;
433 std::cout << Xidrt << std::endl;
434 std::cout <<
"translation-rotation [kcal.fs/(mol.A.radian)]:" << std::endl;
435 std::cout << Xidtr << std::endl;
436 std::cout <<
"rotation [kcal.fs/(mol.radian^2)]:" << std::endl;
437 std::cout << Xidrr << std::endl;
438 std::cout <<
"-----------------------------------------\n\n";
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.