OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
RoughShell.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 "hydrodynamics/RoughShell.hpp"
46
47#include "brains/SimInfo.hpp"
48#include "hydrodynamics/ShapeBuilder.hpp"
49
50namespace OpenMD {
51
52 struct BeadLattice {
53 Vector3d origin;
54 RealType radius;
55 bool interior;
56 };
57
58 struct ExteriorFunctor {
59 bool operator()(const BeadLattice& bead) { return !bead.interior; }
60 };
61
62 struct InteriorFunctor {
63 bool operator()(const BeadLattice& bead) { return bead.interior; }
64 };
65
66 std::size_t RoughShell::assignElements() {
67 std::pair<Vector3d, Vector3d> boxBoundary = shape_->getBoundingBox();
68 RealType firstMin = std::min(
69 std::min(floor(boxBoundary.first[0]), floor(boxBoundary.first[1])),
70 floor(boxBoundary.first[2]));
71 RealType secondMax = std::max(
72 std::max(ceil(boxBoundary.second[0]), ceil(boxBoundary.second[1])),
73 ceil(boxBoundary.second[2]));
74 // std::max is a binary function, i.e., it accepts two values at a time
75 // floor() and ceil() functions return an integer number
76 int len = secondMax - firstMin;
77
78 // number of lattices after shifting positions to integer values (to be
79 // symmetric under reflection)
80
81 int numLattices = ceil(len / sigma_) + 2;
82
83 // +2: one extra lattice point to the left (-1) and one extra
84 // lattice point to the right (+1)
85
86 Grid3D<BeadLattice> grid(numLattices, numLattices, numLattices);
87
88 // fill beads
89 for (int i = 0; i < numLattices; ++i) {
90 for (int j = 0; j < numLattices; ++j) {
91 for (int k = 0; k < numLattices; ++k) {
92 BeadLattice& currentBead = grid(i, j, k);
93 currentBead.origin =
94 Vector3d((i - 1) * sigma_ + floor(boxBoundary.first[0]),
95 (j - 1) * sigma_ + floor(boxBoundary.first[1]),
96 (k - 1) * sigma_ + floor(boxBoundary.first[2]));
97 currentBead.radius = sigma_ / 2.0;
98 currentBead.interior = shape_->isInterior(grid(i, j, k).origin);
99 // returns True if bead is inside the given shape
100 }
101 }
102 }
103
104 // remove embedded beads
105 for (int i = 0; i < numLattices; ++i) {
106 for (int j = 0; j < numLattices; ++j) {
107 for (int k = 0; k < numLattices; ++k) {
108 std::vector<BeadLattice> neighborCells =
109 grid.getAllNeighbors(i, j, k);
110 // if one of its neighbor cells (beads) is exterior, current cell
111 // (bead) is on the surface; loop over beads' center (lattice point)
112
113 if (grid(i, j, k).interior) {
114 bool allNeighborsAreInterior = true;
115 for (std::vector<BeadLattice>::iterator l = neighborCells.begin();
116 l != neighborCells.end(); ++l) {
117 if (!l->interior) {
118 allNeighborsAreInterior = false;
119 break;
120 }
121 }
122
123 if (allNeighborsAreInterior)
124 continue; // if allNeighborsAreInterior == true, skip
125 // the remaining code below, i.e., bead is
126 // not on the surface
127
128 HydrodynamicsElement surfaceBead;
129 // if allNeighborsAreInterior == false, bead is on the
130 // surface; grab this lattice point
131 surfaceBead.name = "H";
132 surfaceBead.pos = grid(i, j, k).origin;
133 // loop over the i, j, k positions (lattice point) of the
134 // grid (loop is above)
135 surfaceBead.radius = grid(i, j, k).radius;
136 elements_.push_back(surfaceBead);
137 }
138 }
139 }
140 }
141 return elements_.size();
142 }
143} // namespace OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.