ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/primitives/Molecule.cpp
Revision: 1360
Committed: Mon Sep 7 16:31:51 2009 UTC (15 years, 7 months ago) by cli2
Original Path: trunk/src/primitives/Molecule.cpp
File size: 10102 byte(s)
Log Message:
Added new restraint infrastructure
Added MolecularRestraints
Added ObjectRestraints
Added RestraintStamp
Updated thermodynamic integration to use ObjectRestraints
Added Quaternion mathematics for twist swing decompositions
Significantly updated RestWriter and RestReader to use dump-like files
Added selections for x, y, and z coordinates of atoms
Removed monolithic Restraints class
Fixed a few bugs in gradients of Euler angles in DirectionalAtom and RigidBody
Added some rotational capabilities to prinicpalAxisCalculator

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