ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/primitives/Molecule.cpp
Revision: 1665
Committed: Tue Nov 22 20:38:56 2011 UTC (13 years, 5 months ago) by gezelter
File size: 10211 byte(s)
Log Message:
updated copyright notices

File Contents

# User Rev Content
1 gezelter 507 /*
2 gezelter 246 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     *
4     * The University of Notre Dame grants you ("Licensee") a
5     * non-exclusive, royalty free, license to use, modify and
6     * redistribute this software in source and binary code form, provided
7     * that the following conditions are met:
8     *
9 gezelter 1390 * 1. Redistributions of source code must retain the above copyright
10 gezelter 246 * notice, this list of conditions and the following disclaimer.
11     *
12 gezelter 1390 * 2. Redistributions in binary form must reproduce the above copyright
13 gezelter 246 * notice, this list of conditions and the following disclaimer in the
14     * documentation and/or other materials provided with the
15     * distribution.
16     *
17     * This software is provided "AS IS," without a warranty of any
18     * kind. All express or implied conditions, representations and
19     * warranties, including any implied warranty of merchantability,
20     * fitness for a particular purpose or non-infringement, are hereby
21     * excluded. The University of Notre Dame and its licensors shall not
22     * be liable for any damages suffered by licensee as a result of
23     * using, modifying or distributing the software or its
24     * derivatives. In no event will the University of Notre Dame or its
25     * licensors be liable for any lost revenue, profit or data, or for
26     * direct, indirect, special, consequential, incidental or punitive
27     * damages, however caused and regardless of the theory of liability,
28     * arising out of the use of or inability to use software, even if the
29     * University of Notre Dame has been advised of the possibility of
30     * such damages.
31 gezelter 1390 *
32     * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33     * research, please cite the appropriate papers when you publish your
34     * work. Good starting points are:
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 gezelter 1665 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40     * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 gezelter 246 */
42    
43     /**
44     * @file Molecule.cpp
45     * @author tlin
46     * @date 10/28/2004
47     * @version 1.0
48     */
49 gezelter 2
50 gezelter 246 #include <algorithm>
51     #include <set>
52 gezelter 2
53 tim 3 #include "primitives/Molecule.hpp"
54 gezelter 246 #include "utils/MemoryUtils.hpp"
55 tim 3 #include "utils/simError.h"
56 gezelter 2
57 gezelter 1390 namespace OpenMD {
58 gezelter 507 Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
59 gezelter 246 : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
60 gezelter 1211 }
61    
62 gezelter 507 Molecule::~Molecule() {
63 gezelter 1211
64 tim 398 MemoryUtils::deletePointers(atoms_);
65     MemoryUtils::deletePointers(bonds_);
66     MemoryUtils::deletePointers(bends_);
67     MemoryUtils::deletePointers(torsions_);
68 gezelter 1277 MemoryUtils::deletePointers(inversions_);
69 tim 398 MemoryUtils::deletePointers(rigidBodies_);
70     MemoryUtils::deletePointers(cutoffGroups_);
71     MemoryUtils::deletePointers(constraintPairs_);
72     MemoryUtils::deletePointers(constraintElems_);
73 gezelter 1211 // integrableObjects_ don't own the objects
74 gezelter 246 integrableObjects_.clear();
75    
76 gezelter 507 }
77 gezelter 1211
78 gezelter 507 void Molecule::addAtom(Atom* atom) {
79 gezelter 246 if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
80 gezelter 507 atoms_.push_back(atom);
81 gezelter 246 }
82 gezelter 507 }
83 gezelter 1211
84 gezelter 507 void Molecule::addBond(Bond* bond) {
85 gezelter 246 if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
86 gezelter 507 bonds_.push_back(bond);
87 gezelter 246 }
88 gezelter 507 }
89 gezelter 1211
90 gezelter 507 void Molecule::addBend(Bend* bend) {
91 gezelter 246 if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
92 gezelter 507 bends_.push_back(bend);
93 gezelter 246 }
94 gezelter 507 }
95 gezelter 1211
96 gezelter 507 void Molecule::addTorsion(Torsion* torsion) {
97 gezelter 1211 if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
98     torsions_.end()) {
99 gezelter 507 torsions_.push_back(torsion);
100 gezelter 246 }
101 gezelter 507 }
102 gezelter 1277
103     void Molecule::addInversion(Inversion* inversion) {
104     if (std::find(inversions_.begin(), inversions_.end(), inversion) ==
105     inversions_.end()) {
106     inversions_.push_back(inversion);
107     }
108     }
109 gezelter 1211
110 gezelter 507 void Molecule::addRigidBody(RigidBody *rb) {
111 gezelter 1211 if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
112     rigidBodies_.end()) {
113 gezelter 507 rigidBodies_.push_back(rb);
114 gezelter 246 }
115 gezelter 507 }
116 gezelter 1211
117 gezelter 507 void Molecule::addCutoffGroup(CutoffGroup* cp) {
118 gezelter 1211 if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
119     cutoffGroups_.end()) {
120 gezelter 507 cutoffGroups_.push_back(cp);
121 gezelter 1211 }
122 gezelter 507 }
123 gezelter 1211
124 gezelter 507 void Molecule::addConstraintPair(ConstraintPair* cp) {
125 gezelter 1211 if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
126     constraintPairs_.end()) {
127 gezelter 507 constraintPairs_.push_back(cp);
128 gezelter 1211 }
129 gezelter 507 }
130 gezelter 1211
131 gezelter 507 void Molecule::addConstraintElem(ConstraintElem* cp) {
132 gezelter 1211 if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
133     constraintElems_.end()) {
134 gezelter 507 constraintElems_.push_back(cp);
135 gezelter 246 }
136 gezelter 507 }
137 gezelter 1211
138 gezelter 507 void Molecule::complete() {
139 gezelter 246
140     std::set<Atom*> rigidAtoms;
141     RigidBody* rb;
142     std::vector<RigidBody*>::iterator rbIter;
143 gezelter 2
144 gezelter 246
145     for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
146 gezelter 507 rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
147 gezelter 246 }
148 gezelter 1211
149 gezelter 246 Atom* atom;
150     AtomIterator ai;
151     for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
152 gezelter 1211
153 gezelter 507 if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
154 gezelter 1211
155     // If an atom does not belong to a rigid body, it is an
156     // integrable object
157    
158 gezelter 507 integrableObjects_.push_back(*ai);
159     }
160 gezelter 246 }
161 gezelter 1211
162 gezelter 246 //find all free atoms (which do not belong to rigid bodies)
163 gezelter 1211 // performs the "difference" operation from set theory, the output
164     // range contains a copy of every element that is contained in
165     // [allAtoms.begin(), allAtoms.end()) and not contained in
166     // [rigidAtoms.begin(), rigidAtoms.end()).
167 gezelter 246 //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
168     // std::back_inserter(integrableObjects_));
169 gezelter 2
170 gezelter 246 //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
171     // //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
172     // sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
173     //
174     // painCave.isFatal = 1;
175     // simError();
176     //}
177 tim 372 for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
178     integrableObjects_.push_back(rb);
179     }
180     //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
181 gezelter 507 }
182 gezelter 2
183 tim 963 RealType Molecule::getMass() {
184 gezelter 246 StuntDouble* sd;
185     std::vector<StuntDouble*>::iterator i;
186 tim 963 RealType mass = 0.0;
187 gezelter 1211
188     for (sd = beginIntegrableObject(i); sd != NULL; sd =
189     nextIntegrableObject(i)){
190 gezelter 507 mass += sd->getMass();
191 gezelter 246 }
192 gezelter 1211
193     return mass;
194 gezelter 507 }
195 gezelter 2
196 gezelter 507 Vector3d Molecule::getCom() {
197 gezelter 246 StuntDouble* sd;
198     std::vector<StuntDouble*>::iterator i;
199     Vector3d com;
200 tim 963 RealType totalMass = 0;
201     RealType mass;
202 gezelter 246
203 gezelter 1211 for (sd = beginIntegrableObject(i); sd != NULL; sd =
204     nextIntegrableObject(i)){
205 gezelter 507 mass = sd->getMass();
206     totalMass += mass;
207     com += sd->getPos() * mass;
208 gezelter 2 }
209 gezelter 1211
210 gezelter 246 com /= totalMass;
211 gezelter 2
212 gezelter 246 return com;
213 gezelter 507 }
214 gezelter 2
215 gezelter 507 void Molecule::moveCom(const Vector3d& delta) {
216 gezelter 246 StuntDouble* sd;
217     std::vector<StuntDouble*>::iterator i;
218    
219 gezelter 1211 for (sd = beginIntegrableObject(i); sd != NULL; sd =
220     nextIntegrableObject(i)){
221 gezelter 507 sd->setPos(sd->getPos() + delta);
222 gezelter 1211 }
223 gezelter 507 }
224 gezelter 2
225 gezelter 507 Vector3d Molecule::getComVel() {
226 gezelter 246 StuntDouble* sd;
227     std::vector<StuntDouble*>::iterator i;
228     Vector3d velCom;
229 tim 963 RealType totalMass = 0;
230     RealType mass;
231 gezelter 246
232 gezelter 1211 for (sd = beginIntegrableObject(i); sd != NULL; sd =
233     nextIntegrableObject(i)){
234 gezelter 507 mass = sd->getMass();
235     totalMass += mass;
236     velCom += sd->getVel() * mass;
237 gezelter 246 }
238 gezelter 1211
239 gezelter 246 velCom /= totalMass;
240 gezelter 1211
241 gezelter 246 return velCom;
242 gezelter 507 }
243 gezelter 2
244 tim 963 RealType Molecule::getPotential() {
245 gezelter 2
246 gezelter 246 Bond* bond;
247     Bend* bend;
248     Torsion* torsion;
249 gezelter 1277 Inversion* inversion;
250 gezelter 246 Molecule::BondIterator bondIter;;
251     Molecule::BendIterator bendIter;
252     Molecule::TorsionIterator torsionIter;
253 gezelter 1277 Molecule::InversionIterator inversionIter;
254 gezelter 2
255 tim 963 RealType potential = 0.0;
256 gezelter 2
257 gezelter 246 for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
258 gezelter 507 potential += bond->getPotential();
259 gezelter 246 }
260 gezelter 2
261 gezelter 246 for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
262 gezelter 507 potential += bend->getPotential();
263 gezelter 2 }
264    
265 gezelter 1211 for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
266     nextTorsion(torsionIter)) {
267 gezelter 507 potential += torsion->getPotential();
268 gezelter 2 }
269 gezelter 1211
270 gezelter 1277 for (inversion = beginInversion(inversionIter); torsion != NULL;
271     inversion = nextInversion(inversionIter)) {
272     potential += inversion->getPotential();
273     }
274    
275 gezelter 246 return potential;
276 gezelter 1211
277 gezelter 507 }
278 gezelter 1211
279 cli2 1360 void Molecule::addProperty(GenericData* genData) {
280     properties_.addProperty(genData);
281     }
282    
283     void Molecule::removeProperty(const std::string& propName) {
284     properties_.removeProperty(propName);
285     }
286    
287     void Molecule::clearProperties() {
288     properties_.clearProperties();
289     }
290    
291     std::vector<std::string> Molecule::getPropertyNames() {
292     return properties_.getPropertyNames();
293     }
294    
295     std::vector<GenericData*> Molecule::getProperties() {
296     return properties_.getProperties();
297     }
298    
299     GenericData* Molecule::getPropertyByName(const std::string& propName) {
300     return properties_.getPropertyByName(propName);
301     }
302    
303    
304    
305    
306 gezelter 507 std::ostream& operator <<(std::ostream& o, Molecule& mol) {
307 gezelter 246 o << std::endl;
308     o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
309     o << mol.getNAtoms() << " atoms" << std::endl;
310     o << mol.getNBonds() << " bonds" << std::endl;
311     o << mol.getNBends() << " bends" << std::endl;
312     o << mol.getNTorsions() << " torsions" << std::endl;
313 gezelter 1277 o << mol.getNInversions() << " inversions" << std::endl;
314 gezelter 246 o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
315     o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
316     o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
317     o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
318     return o;
319 gezelter 507 }
320 gezelter 1211
321 gezelter 1390 }//end namespace OpenMD

Properties

Name Value
svn:keywords Author Id Revision Date