OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
OtherVisitor.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 "visitors/OtherVisitor.hpp"
46
47#include <memory>
48
49#include "brains/SimInfo.hpp"
50#include "brains/Thermo.hpp"
54#include "selection/SelectionManager.hpp"
55
56namespace OpenMD {
57
58 void WrappingVisitor::visit(Atom* atom) { internalVisit(atom); }
59
60 void WrappingVisitor::visit(DirectionalAtom* datom) { internalVisit(datom); }
61
62 void WrappingVisitor::visit(RigidBody* rb) { internalVisit(rb); }
63
64 void WrappingVisitor::internalVisit(StuntDouble* sd) {
65 std::shared_ptr<GenericData> data;
66 std::shared_ptr<AtomData> atomData;
67 std::shared_ptr<AtomInfo> atomInfo;
68 std::vector<std::shared_ptr<AtomInfo>>::iterator i;
69
70 data = sd->getPropertyByName("ATOMDATA");
71
72 if (data != nullptr) {
73 atomData = std::dynamic_pointer_cast<AtomData>(data);
74
75 if (atomData == nullptr) return;
76 } else
77 return;
78
79 Snapshot* currSnapshot = info->getSnapshotManager()->getCurrentSnapshot();
80
81 for (atomInfo = atomData->beginAtomInfo(i); atomInfo;
82 atomInfo = atomData->nextAtomInfo(i)) {
83 Vector3d newPos = atomInfo->pos - origin_;
84 currSnapshot->wrapVector(newPos);
85 atomInfo->pos = newPos;
86 }
87 }
88
89 void WrappingVisitor::update() {
90 if (useCom_) {
91 Thermo thermo(info);
92 origin_ = thermo.getCom();
93 }
94 }
95
96 const std::string WrappingVisitor::toString() {
97 char buffer[65535];
98 std::string result;
99
100 snprintf(
101 buffer, 65535,
102 "------------------------------------------------------------------\n");
103 result += buffer;
104
105 snprintf(buffer, 65535, "Visitor name: %s\n", visitorName.c_str());
106 result += buffer;
107
108 snprintf(buffer, 65535,
109 "Visitor Description: wrapping atoms back to periodic box\n");
110 result += buffer;
111
112 snprintf(
113 buffer, 65535,
114 "------------------------------------------------------------------\n");
115 result += buffer;
116
117 return result;
118 }
119
120 //----------------------------------------------------------------------------//
121
122 ReplicateVisitor::ReplicateVisitor(SimInfo* info, Vector3i opt) :
123 BaseVisitor(), replicateOpt(opt) {
124 this->info = info;
125 visitorName = "ReplicateVisitor";
126
127 // generate the replicate directions
128 for (int i = 0; i <= replicateOpt[0]; i++) {
129 for (int j = 0; j <= replicateOpt[1]; j++) {
130 for (int k = 0; k <= replicateOpt[2]; k++) {
131 // skip original frame
132 if (i == 0 && j == 0 && k == 0) {
133 continue;
134 } else {
135 dir.push_back(Vector3d((RealType)i, (RealType)j, (RealType)k));
136 }
137 }
138 }
139 }
140 }
141
142 void ReplicateVisitor::visit(Atom* atom) { internalVisit(atom); }
143
144 void ReplicateVisitor::visit(DirectionalAtom* datom) { internalVisit(datom); }
145
146 void ReplicateVisitor::visit(RigidBody* rb) { internalVisit(rb); }
147
148 void ReplicateVisitor::internalVisit(StuntDouble* sd) {
149 std::shared_ptr<GenericData> data;
150 std::shared_ptr<AtomData> atomData;
151
152 // if there is not atom data, just skip it
153 data = sd->getPropertyByName("ATOMDATA");
154
155 if (data != nullptr) {
156 atomData = std::dynamic_pointer_cast<AtomData>(data);
157
158 if (atomData == nullptr) { return; }
159 } else {
160 return;
161 }
162
163 Snapshot* currSnapshot = info->getSnapshotManager()->getCurrentSnapshot();
164 Mat3x3d box = currSnapshot->getHmat();
165
166 std::vector<std::shared_ptr<AtomInfo>> atomInfoList = atomData->getData();
167
168 replicate(atomInfoList, atomData, box);
169 }
170
171 void ReplicateVisitor::replicate(
172 std::vector<std::shared_ptr<AtomInfo>>& infoList,
173 std::shared_ptr<AtomData> data, const Mat3x3d& box) {
174 std::shared_ptr<AtomInfo> newAtomInfo;
175 std::vector<Vector3d>::iterator dirIter;
176 std::vector<std::shared_ptr<AtomInfo>>::iterator i;
177
178 for (dirIter = dir.begin(); dirIter != dir.end(); ++dirIter) {
179 for (i = infoList.begin(); i != infoList.end(); ++i) {
180 newAtomInfo = std::make_shared<AtomInfo>();
181 *newAtomInfo = *(*i);
182 newAtomInfo->pos += box * (*dirIter);
183 data->addAtomInfo(newAtomInfo);
184 }
185 }
186 }
187
188 const std::string ReplicateVisitor::toString() {
189 char buffer[65535];
190 std::string result;
191
192 snprintf(
193 buffer, 65535,
194 "--------------------------------------------------------------\n");
195 result += buffer;
196
197 snprintf(buffer, 65535, "Visitor name: %s\n", visitorName.c_str());
198 result += buffer;
199
200 snprintf(
201 buffer, 65535,
202 "Visitor Description: replicate the atoms in different direction\n");
203 result += buffer;
204
205 // print the replicate direction
206 snprintf(buffer, 65535, "repeatX = %d:\n", replicateOpt[0]);
207 result += buffer;
208
209 snprintf(buffer, 65535, "repeatY = %d:\n", replicateOpt[1]);
210 result += buffer;
211
212 snprintf(buffer, 65535, "repeatZ = %d:\n", replicateOpt[2]);
213 result += buffer;
214
215 snprintf(
216 buffer, 65535,
217 "--------------------------------------------------------------\n");
218 result += buffer;
219
220 return result;
221 }
222
223 //------------------------------------------------------------------------//
224
225 XYZVisitor::XYZVisitor(SimInfo* info) :
226 BaseVisitor(), seleMan(info), evaluator(info), doPositions_(true),
227 doVelocities_(false), doForces_(false), doVectors_(false),
228 doCharges_(false), doElectricFields_(false), doGlobalIDs_(false) {
229 this->info = info;
230 visitorName = "XYZVisitor";
231
232 evaluator.loadScriptString("select all");
233
234 if (!evaluator.isDynamic()) {
235 seleMan.setSelectionSet(evaluator.evaluate());
236 }
237 }
238
239 XYZVisitor::XYZVisitor(SimInfo* info, const std::string& script) :
240 BaseVisitor(), seleMan(info), evaluator(info), doPositions_(true),
241 doVelocities_(false), doForces_(false), doVectors_(false),
242 doCharges_(false), doElectricFields_(false), doGlobalIDs_(false) {
243 this->info = info;
244 visitorName = "XYZVisitor";
245
246 evaluator.loadScriptString(script);
247
248 if (!evaluator.isDynamic()) {
249 seleMan.setSelectionSet(evaluator.evaluate());
250 }
251 }
252
253 void XYZVisitor::visit(Atom* atom) {
254 if (isSelected(atom)) internalVisit(atom);
255 }
256
257 void XYZVisitor::visit(DirectionalAtom* datom) {
258 if (isSelected(datom)) internalVisit(datom);
259 }
260
261 void XYZVisitor::visit(RigidBody* rb) {
262 if (isSelected(rb)) internalVisit(rb);
263 }
264
265 void XYZVisitor::update() {
266 // if dynamic, we need to re-evaluate the selection
267 if (evaluator.isDynamic()) {
268 seleMan.setSelectionSet(evaluator.evaluate());
269 }
270 }
271
272 void XYZVisitor::internalVisit(StuntDouble* sd) {
273 std::shared_ptr<GenericData> data;
274 std::shared_ptr<AtomData> atomData;
275 std::shared_ptr<AtomInfo> atomInfo;
276 std::vector<std::shared_ptr<AtomInfo>>::iterator i;
277 char buffer[1024];
278
279 // if there is not atom data, just skip it
280 data = sd->getPropertyByName("ATOMDATA");
281
282 if (data != nullptr) {
283 atomData = std::dynamic_pointer_cast<AtomData>(data);
284
285 if (atomData == nullptr) return;
286 } else
287 return;
288
289 for (atomInfo = atomData->beginAtomInfo(i); atomInfo;
290 atomInfo = atomData->nextAtomInfo(i)) {
291 std::string line;
292 snprintf(buffer, 1024, "%s", atomInfo->atomTypeName.c_str());
293 line += buffer;
294
295 if (doPositions_) {
296 snprintf(buffer, 1024, "%15.8f%15.8f%15.8f", atomInfo->pos[0],
297 atomInfo->pos[1], atomInfo->pos[2]);
298 line += buffer;
299 }
300 if (doCharges_ && atomInfo->hasCharge) {
301 snprintf(buffer, 1024, "%15.8f", atomInfo->charge);
302 line += buffer;
303 }
304 if (doVectors_ && atomInfo->hasVector) {
305 snprintf(buffer, 1024, "%15.8f%15.8f%15.8f", atomInfo->vec[0],
306 atomInfo->vec[1], atomInfo->vec[2]);
307 line += buffer;
308 }
309 if (doVelocities_ && atomInfo->hasVelocity) {
310 snprintf(buffer, 1024, "%15.8f%15.8f%15.8f", atomInfo->vel[0],
311 atomInfo->vel[1], atomInfo->vel[2]);
312 line += buffer;
313 }
314 if (doForces_ && atomInfo->hasForce) {
315 snprintf(buffer, 1024, "%15.8f%15.8f%15.8f", atomInfo->frc[0],
316 atomInfo->frc[1], atomInfo->frc[2]);
317 line += buffer;
318 }
319 if (doElectricFields_ && atomInfo->hasElectricField) {
320 snprintf(buffer, 1024, "%15.8f%15.8f%15.8f", atomInfo->eField[0],
321 atomInfo->eField[1], atomInfo->eField[2]);
322 line += buffer;
323 }
324 if (doGlobalIDs_) {
325 snprintf(buffer, 1024, "%10d", atomInfo->globalID);
326 line += buffer;
327 }
328
329 frame.push_back(line);
330 }
331 }
332
333 bool XYZVisitor::isSelected(StuntDouble* sd) {
334 return seleMan.isSelected(sd);
335 }
336
337 void XYZVisitor::writeFrame(std::ostream& outStream) {
338 std::vector<std::string>::iterator i;
339 char buffer[1024];
340
341 if (frame.empty())
342 std::cerr << "Current Frame does not contain any atoms" << std::endl;
343
344 // total number of atoms
345 outStream << frame.size() << std::endl;
346
347 // write comment line
348 Snapshot* currSnapshot = info->getSnapshotManager()->getCurrentSnapshot();
349 Mat3x3d box = currSnapshot->getHmat();
350
351 snprintf(buffer, 1024,
352 "%15.8f;%15.8f%15.8f%15.8f;%15.8f%15.8f%15.8f;%15.8f%15.8f%15.8f",
353 currSnapshot->getTime(), box(0, 0), box(0, 1), box(0, 2),
354 box(1, 0), box(1, 1), box(1, 2), box(2, 0), box(2, 1), box(2, 2));
355
356 outStream << buffer << std::endl;
357
358 for (i = frame.begin(); i != frame.end(); ++i)
359 outStream << *i << std::endl;
360 }
361
362 std::string XYZVisitor::trimmedName(const std::string& atomTypeName) {
363 return atomTypeName.substr(0, atomTypeName.find('-'));
364 }
365
366 const std::string XYZVisitor::toString() {
367 char buffer[65535];
368 std::string result;
369
370 snprintf(
371 buffer, 65535,
372 "------------------------------------------------------------------\n");
373 result += buffer;
374
375 snprintf(buffer, 65535, "Visitor name: %s\n", visitorName.c_str());
376 result += buffer;
377
378 snprintf(
379 buffer, 65535,
380 "Visitor Description: assemble the atom data and output xyz file\n");
381 result += buffer;
382
383 snprintf(
384 buffer, 65535,
385 "------------------------------------------------------------------\n");
386 result += buffer;
387
388 return result;
389 }
390
391 //----------------------------------------------------------------------------//
392
393 void PrepareVisitor::internalVisit(Atom* atom) {
394 std::shared_ptr<GenericData> data;
395
396 // if visited property is existed, remove it
397 data = atom->getPropertyByName("VISITED");
398
399 if (data != nullptr) { atom->removeProperty("VISITED"); }
400
401 // remove atomdata
402 data = atom->getPropertyByName("ATOMDATA");
403
404 if (data != nullptr) {
405 std::shared_ptr<AtomData> atomData =
406 std::dynamic_pointer_cast<AtomData>(data);
407
408 if (atomData != NULL) atom->removeProperty("ATOMDATA");
409 }
410 }
411
412 void PrepareVisitor::internalVisit(RigidBody* rb) {
413 std::shared_ptr<GenericData> data;
414 std::shared_ptr<AtomData> atomData;
415 std::vector<Atom*> myAtoms;
416 std::vector<Atom*>::iterator atomIter;
417
418 // if visited property is existed, remove it
419 data = rb->getPropertyByName("VISITED");
420
421 if (data != nullptr) { rb->removeProperty("VISITED"); }
422
423 // remove atomdata
424 data = rb->getPropertyByName("ATOMDATA");
425
426 if (data != nullptr) {
427 atomData = std::dynamic_pointer_cast<AtomData>(data);
428
429 if (atomData != NULL) rb->removeProperty("ATOMDATA");
430 }
431
432 myAtoms = rb->getAtoms();
433
434 for (atomIter = myAtoms.begin(); atomIter != myAtoms.end(); ++atomIter)
435 internalVisit(*atomIter);
436 }
437
438 const std::string PrepareVisitor::toString() {
439 char buffer[65535];
440 std::string result;
441
442 snprintf(
443 buffer, 65535,
444 "------------------------------------------------------------------\n");
445 result += buffer;
446
447 snprintf(buffer, 65535, "Visitor name: %s", visitorName.c_str());
448 result += buffer;
449
450 snprintf(buffer, 65535,
451 "Visitor Description: prepare for operation of other vistors\n");
452 result += buffer;
453
454 snprintf(
455 buffer, 65535,
456 "------------------------------------------------------------------\n");
457 result += buffer;
458
459 return result;
460 }
461
462 //----------------------------------------------------------------------------//
463
464 WaterTypeVisitor::WaterTypeVisitor() {
465 visitorName = "WaterTypeVisitor";
466 waterTypeList.insert("TIP3P_RB_0");
467 waterTypeList.insert("TIP4P_RB_0");
468 waterTypeList.insert("TIP4P-Ice_RB_0");
469 waterTypeList.insert("TIP4P-Ew_RB_0");
470 waterTypeList.insert("TIP4P-2005_RB_0");
471 waterTypeList.insert("TIP5P_RB_0");
472 waterTypeList.insert("TIP5P-E_RB_0");
473 waterTypeList.insert("SPCE_RB_0");
474 waterTypeList.insert("SPC_RB_0");
475 waterTypeList.insert("SPC-HW_RB_0");
476 waterTypeList.insert("NE6_RB_0");
477 }
478
479 void WaterTypeVisitor::visit(RigidBody* rb) {
480 std::string rbName;
481 std::vector<Atom*> myAtoms;
482 std::vector<Atom*>::iterator atomIter;
483 std::vector<std::shared_ptr<AtomInfo>>::iterator i;
484 std::shared_ptr<AtomData> atomData;
485
486 rbName = rb->getType();
487
488 if (waterTypeList.find(rbName) != waterTypeList.end()) {
489 myAtoms = rb->getAtoms();
490
491 for (atomIter = myAtoms.begin(); atomIter != myAtoms.end(); ++atomIter) {
492 std::shared_ptr<GenericData> data =
493 (*atomIter)->getPropertyByName("ATOMDATA");
494
495 if (data != nullptr) {
496 atomData = std::dynamic_pointer_cast<AtomData>(data);
497
498 if (atomData == nullptr) continue;
499 } else
500 continue;
501
502 for (std::shared_ptr<AtomInfo> atomInfo = atomData->beginAtomInfo(i);
503 atomInfo; atomInfo = atomData->nextAtomInfo(i)) {
504 atomInfo->atomTypeName = trimmedName(atomInfo->atomTypeName);
505 }
506 }
507 }
508 }
509
510 std::string WaterTypeVisitor::trimmedName(const std::string& atomTypeName) {
511 return atomTypeName.substr(0, atomTypeName.find('_'));
512 }
513
514 const std::string WaterTypeVisitor::toString() {
515 char buffer[65535];
516 std::string result;
517
518 snprintf(
519 buffer, 65535,
520 "------------------------------------------------------------------\n");
521 result += buffer;
522
523 snprintf(buffer, 65535, "Visitor name: %s\n", visitorName.c_str());
524 result += buffer;
525
526 snprintf(buffer, 65535,
527 "Visitor Description: Replace the atom type in water model\n");
528 result += buffer;
529
530 snprintf(
531 buffer, 65535,
532 "------------------------------------------------------------------\n");
533 result += buffer;
534
535 return result;
536 }
537
538} // namespace OpenMD
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
Definition SimInfo.hpp:248
Snapshot * getCurrentSnapshot()
Returns the pointer of current snapshot.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.