OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
MoLocator.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 "utils/MoLocator.hpp"
46
47#include <cmath>
48#include <cstdlib>
49#include <iostream>
50
51#include "types/AtomType.hpp"
52#include "utils/simError.h"
53
54namespace OpenMD {
55 MoLocator::MoLocator(MoleculeStamp* theStamp, ForceField* theFF) {
56 myStamp = theStamp;
57 myFF = theFF;
58 nIntegrableObjects = myStamp->getNIntegrable();
59 calcRef();
60 }
61
62 void MoLocator::placeMol(const Vector3d& offset, const Vector3d& ort,
63 Molecule* mol) {
64 Vector3d newCoor;
65 Vector3d curRefCoor;
66 RotMat3x3d rotMat = latVec2RotMat(ort);
67
68 if (mol->getNIntegrableObjects() != nIntegrableObjects) {
69 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
70 "MoLocator::placeMol error.\n"
71 "\tThe number of integrable objects of MoleculeStamp is not\n"
72 "\tthe same as that of Molecule\n");
73 painCave.isFatal = 1;
74 simError();
75 }
76
77 Molecule::IntegrableObjectIterator ii;
78 StuntDouble* sd;
79 int i;
80 for (sd = mol->beginIntegrableObject(ii), i = 0; sd != NULL;
81 sd = mol->nextIntegrableObject(ii), ++i) {
82 newCoor = rotMat * refCoords[i];
83 newCoor += offset;
84
85 sd->setPos(newCoor);
86 sd->setVel(V3Zero);
87
88 if (sd->isDirectional()) {
89 sd->setA(rotMat * sd->getA());
90 sd->setJ(V3Zero);
91 }
92 }
93 }
94
95 void MoLocator::calcRef(void) {
96 AtomStamp* currAtomStamp;
97 RigidBodyStamp* rbStamp;
98 std::vector<RealType> mass;
99 Vector3d coor;
100 Vector3d refMolCom;
101 RealType totMassInRb;
102 RealType currAtomMass;
103 RealType molMass;
104
105 std::size_t nAtoms = myStamp->getNAtoms();
106 std::size_t nRigidBodies = myStamp->getNRigidBodies();
107
108 for (std::size_t i = 0; i < nAtoms; i++) {
109 currAtomStamp = myStamp->getAtomStamp(i);
110
111 if (!currAtomStamp->havePosition()) {
112 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
113 "MoLocator::calcRef error.\n"
114 "\tComponent %s, atom %s does not have a position specified.\n"
115 "\tThis means MoLocator cannot initalize it's position.\n",
116 myStamp->getName().c_str(), currAtomStamp->getType().c_str());
117
118 painCave.isFatal = 1;
119 simError();
120 }
121
122 // if atom belongs to rigidbody, just skip it
123 if (myStamp->isAtomInRigidBody(i)) continue;
124 // get mass and the reference coordinate
125 else {
126 currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);
127 mass.push_back(currAtomMass);
128 coor.x() = currAtomStamp->getPosX();
129 coor.y() = currAtomStamp->getPosY();
130 coor.z() = currAtomStamp->getPosZ();
131 refCoords.push_back(coor);
132 }
133 }
134
135 for (std::size_t i = 0; i < nRigidBodies; i++) {
136 rbStamp = myStamp->getRigidBodyStamp(i);
137 std::size_t nAtomsInRb = rbStamp->getNMembers();
138
139 coor.x() = 0.0;
140 coor.y() = 0.0;
141 coor.z() = 0.0;
142 totMassInRb = 0.0;
143
144 for (std::size_t j = 0; j < nAtomsInRb; j++) {
145 currAtomStamp = myStamp->getAtomStamp(rbStamp->getMemberAt(j));
146 currAtomMass = getAtomMass(currAtomStamp->getType(), myFF);
147 totMassInRb += currAtomMass;
148
149 coor.x() += currAtomStamp->getPosX() * currAtomMass;
150 coor.y() += currAtomStamp->getPosY() * currAtomMass;
151 coor.z() += currAtomStamp->getPosZ() * currAtomMass;
152 }
153
154 mass.push_back(totMassInRb);
155 coor /= totMassInRb;
156 refCoords.push_back(coor);
157 }
158
159 // calculate the reference center of mass
160 molMass = 0;
161 refMolCom.x() = 0;
162 refMolCom.y() = 0;
163 refMolCom.z() = 0;
164
165 for (std::size_t i = 0; i < nIntegrableObjects; i++) {
166 refMolCom += refCoords[i] * mass[i];
167 molMass += mass[i];
168 }
169
170 refMolCom /= molMass;
171
172 // move the reference center of mass to (0,0,0) and adjust the
173 // reference coordinate of the integrabel objects
174 for (std::size_t i = 0; i < nIntegrableObjects; i++)
175 refCoords[i] -= refMolCom;
176 }
177
178 RealType MoLocator::getAtomMass(const std::string& at, ForceField* myFF) {
179 RealType mass;
180 AtomType* atomType = myFF->getAtomType(at);
181 if (atomType != NULL) {
182 mass = atomType->getMass();
183 } else {
184 mass = 0.0;
185 std::cerr << "Can not find AtomType: " << at << std::endl;
186 }
187 return mass;
188 }
189
190 RealType MoLocator::getMolMass(MoleculeStamp* molStamp, ForceField* myFF) {
191 unsigned int nAtoms;
192 RealType totMass = 0;
193 nAtoms = molStamp->getNAtoms();
194
195 for (std::size_t i = 0; i < nAtoms; i++) {
196 AtomStamp* currAtomStamp = molStamp->getAtomStamp(i);
197 totMass += getAtomMass(currAtomStamp->getType(), myFF);
198 }
199 return totMass;
200 }
201
202 RotMat3x3d MoLocator::latVec2RotMat(const Vector3d& lv) {
203 RealType theta = acos(lv[2]);
204 RealType phi = atan2(lv[1], lv[0]);
205 RealType psi = 0;
206
207 return RotMat3x3d(phi, theta, psi);
208 }
209} // namespace OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.