ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/primitives/Molecule.cpp
Revision: 1901
Committed: Tue Jan 4 22:18:36 2005 UTC (19 years, 6 months ago) by tim
File size: 9724 byte(s)
Log Message:
constraints in progress

File Contents

# User Rev Content
1 tim 1692 /*
2     * Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3     *
4     * Contact: oopse@oopse.org
5     *
6     * This program is free software; you can redistribute it and/or
7     * modify it under the terms of the GNU Lesser General Public License
8     * as published by the Free Software Foundation; either version 2.1
9     * of the License, or (at your option) any later version.
10     * All we ask is that proper credit is given for our work, which includes
11     * - but is not limited to - adding the above copyright notice to the beginning
12     * of your source code files, and to any copyright notice that you may distribute
13     * with programs based on this work.
14     *
15     * This program is distributed in the hope that it will be useful,
16     * but WITHOUT ANY WARRANTY; without even the implied warranty of
17     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18     * GNU Lesser General Public License for more details.
19     *
20     * You should have received a copy of the GNU Lesser General Public License
21     * along with this program; if not, write to the Free Software
22     * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23     *
24     */
25 gezelter 1490
26 tim 1692 /**
27     * @file Molecule.cpp
28     * @author tlin
29     * @date 10/28/2004
30     * @version 1.0
31     */
32 gezelter 1490
33 tim 1692 #include <algorithm>
34 tim 1695 #include <set>
35 gezelter 1490
36 tim 1695 #include "primitives/Molecule.hpp"
37     #include "utils/MemoryUtils.hpp"
38 tim 1782 #include "utils/simError.h"
39 tim 1695
40 tim 1692 namespace oopse {
41 tim 1733 Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
42     : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
43 gezelter 1490
44 tim 1733 }
45    
46 tim 1692 Molecule::~Molecule() {
47 gezelter 1490
48 tim 1695 MemoryUtils::deleteVectorOfPointer(atoms_);
49     MemoryUtils::deleteVectorOfPointer(bonds_);
50     MemoryUtils::deleteVectorOfPointer(bends_);
51     MemoryUtils::deleteVectorOfPointer(torsions_);
52     MemoryUtils::deleteVectorOfPointer(rigidBodies_);
53     MemoryUtils::deleteVectorOfPointer(cutoffGroups_);
54 gezelter 1490
55 tim 1695 //integrableObjects_ don't own the objects
56 tim 1692 integrableObjects_.clear();
57    
58 gezelter 1490 }
59    
60 tim 1692 void Molecule::addAtom(Atom* atom) {
61 tim 1695 if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
62 tim 1692 atoms_.push_back(atom);
63     }
64     }
65 gezelter 1490
66 tim 1692 void Molecule::addBond(Bond* bond) {
67 tim 1695 if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
68 tim 1692 bonds_.push_back(bond);
69     }
70     }
71 gezelter 1490
72 tim 1692 void Molecule::addBend(Bend* bend) {
73 tim 1695 if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
74 tim 1692 bends_.push_back(bend);
75     }
76     }
77 gezelter 1490
78 tim 1692 void Molecule::addTorsion(Torsion* torsion) {
79 tim 1695 if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
80 tim 1692 torsions_.push_back(torsion);
81     }
82     }
83 gezelter 1490
84 tim 1692 void Molecule::addRigidBody(RigidBody *rb) {
85 tim 1695 if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
86 tim 1692 rigidBodies_.push_back(rb);
87     }
88 gezelter 1490 }
89    
90 tim 1692 void Molecule::addCutoffGroup(CutoffGroup* cp) {
91 tim 1695 if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
92 tim 1692 cutoffGroups_.push_back(cp);
93     }
94 gezelter 1490
95 tim 1692 }
96 gezelter 1490
97 tim 1901 void Molecule::addConstraintPair(ConstraintPair* cp) {
98     if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) {
99     constraintPairs_.push_back(cp);
100     }
101    
102     }
103    
104    
105 tim 1692 void Molecule::complete() {
106    
107     std::set<Atom*> rigidAtoms;
108     RigidBody* rb;
109 tim 1695 std::vector<RigidBody*>::iterator rbIter;
110 gezelter 1490
111 tim 1692
112     for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
113 tim 1695 rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
114 tim 1692 }
115 gezelter 1490
116 tim 1843 Atom* atom;
117     AtomIterator ai;
118     for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
119 tim 1844
120     if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
121     //if an atom does not belong to a rigid body, it is an integrable object
122 tim 1843 integrableObjects_.push_back(*ai);
123     }
124     }
125    
126 tim 1692 //find all free atoms (which do not belong to rigid bodies)
127     //performs the "difference" operation from set theory, the output range contains a copy of every
128     //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
129     //[rigidAtoms.begin(), rigidAtoms.end()).
130 tim 1843 //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
131     // std::back_inserter(integrableObjects_));
132 gezelter 1490
133 tim 1843 //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
134     // //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
135     // sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
136     //
137     // painCave.isFatal = 1;
138     // simError();
139     //}
140 tim 1733
141 tim 1692 integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
142 gezelter 1490 }
143    
144 tim 1692 Atom* Molecule::beginAtom(std::vector<Atom*>::iterator& i) {
145     i = atoms_.begin();
146     return (i == atoms_.end()) ? NULL : *i;
147     }
148 gezelter 1490
149 tim 1692 Atom* Molecule::nextAtom(std::vector<Atom*>::iterator& i) {
150     ++i;
151     return (i == atoms_.end()) ? NULL : *i;
152     }
153 gezelter 1490
154 tim 1692 Bond* Molecule::beginBond(std::vector<Bond*>::iterator& i) {
155     i = bonds_.begin();
156     return (i == bonds_.end()) ? NULL : *i;
157     }
158 gezelter 1490
159 tim 1692 Bond* Molecule::nextBond(std::vector<Bond*>::iterator& i) {
160     ++i;
161     return (i == bonds_.end()) ? NULL : *i;
162 gezelter 1490
163 tim 1692 }
164 gezelter 1490
165    
166 tim 1692 Bend* Molecule::beginBend(std::vector<Bend*>::iterator& i) {
167     i = bends_.begin();
168     return (i == bends_.end()) ? NULL : *i;
169 gezelter 1490 }
170    
171 tim 1692 Bend* Molecule::nextBend(std::vector<Bend*>::iterator& i) {
172     ++i;
173     return (i == bends_.end()) ? NULL : *i;
174     }
175 gezelter 1490
176 tim 1692 Torsion* Molecule::beginTorsion(std::vector<Torsion*>::iterator& i) {
177     i = torsions_.begin();
178     return (i == torsions_.end()) ? NULL : *i;
179     }
180 gezelter 1490
181 tim 1692 Torsion* Molecule::nextTorsion(std::vector<Torsion*>::iterator& i) {
182     ++i;
183     return (i == torsions_.end()) ? NULL : *i;
184     }
185 gezelter 1490
186 tim 1692 RigidBody* Molecule::beginRigidBody(std::vector<RigidBody*>::iterator& i) {
187     i = rigidBodies_.begin();
188     return (i == rigidBodies_.end()) ? NULL : *i;
189     }
190 gezelter 1490
191 tim 1692 RigidBody* Molecule::nextRigidBody(std::vector<RigidBody*>::iterator& i) {
192     ++i;
193     return (i == rigidBodies_.end()) ? NULL : *i;
194     }
195 gezelter 1490
196 tim 1692 StuntDouble* Molecule::beginIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
197     i = integrableObjects_.begin();
198     return (i == integrableObjects_.end()) ? NULL : *i;
199 gezelter 1490 }
200    
201 tim 1692 StuntDouble* Molecule::nextIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
202     ++i;
203     return (i == integrableObjects_.end()) ? NULL : *i;
204     }
205 gezelter 1490
206 tim 1692 CutoffGroup* Molecule::beginCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {
207     i = cutoffGroups_.begin();
208     return (i == cutoffGroups_.end()) ? NULL : *i;
209 gezelter 1490 }
210    
211 tim 1692 CutoffGroup* Molecule::nextCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {
212     ++i;
213     return (i == cutoffGroups_.end()) ? NULL : *i;
214     }
215 gezelter 1490
216 tim 1901 ConstraintPair* Molecule::beginConstraintPair(std::vector<ConstraintPair*>::iterator& i) {
217     i = constraintPairs_.begin();
218     return (i == constraintPairs_.end()) ? NULL : *i;
219     }
220    
221     ConstraintPair* Molecule::nextConstraintPair(std::vector<ConstraintPair*>::iterator& i) {
222     ++i;
223     return (i == constraintPairs_.end()) ? NULL : *i;
224     }
225    
226 tim 1843 double Molecule::getMass() {
227     StuntDouble* sd;
228     std::vector<StuntDouble*>::iterator i;
229     double mass = 0.0;
230 gezelter 1490
231 tim 1843 for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
232     mass += sd->getMass();
233     }
234    
235     return mass;
236    
237     }
238    
239 tim 1692 Vector3d Molecule::getCom() {
240     StuntDouble* sd;
241     std::vector<StuntDouble*>::iterator i;
242     Vector3d com;
243     double totalMass = 0;
244     double mass;
245    
246     for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
247     mass = sd->getMass();
248     totalMass += mass;
249     com += sd->getPos() * mass;
250     }
251 gezelter 1490
252 tim 1692 com /= totalMass;
253 gezelter 1490
254 tim 1692 return com;
255     }
256 gezelter 1490
257 tim 1695 void Molecule::moveCom(const Vector3d& delta) {
258 tim 1692 StuntDouble* sd;
259     std::vector<StuntDouble*>::iterator i;
260    
261     for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
262 tim 1695 sd->setPos(sd->getPos() + delta);
263 gezelter 1490 }
264    
265     }
266    
267 tim 1692 Vector3d Molecule::getComVel() {
268     StuntDouble* sd;
269     std::vector<StuntDouble*>::iterator i;
270     Vector3d velCom;
271     double totalMass = 0;
272     double mass;
273    
274     for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
275     mass = sd->getMass();
276     totalMass += mass;
277     velCom += sd->getVel() * mass;
278 gezelter 1490 }
279    
280 tim 1692 velCom /= totalMass;
281 gezelter 1490
282 tim 1692 return velCom;
283 gezelter 1490 }
284    
285 tim 1843 double Molecule::getPotential() {
286    
287     Bond* bond;
288     Bend* bend;
289     Torsion* torsion;
290     Molecule::BondIterator bondIter;;
291     Molecule::BendIterator bendIter;
292     Molecule::TorsionIterator torsionIter;
293    
294     double potential = 0.0;
295    
296     for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
297     potential += bond->getPotential();
298     }
299    
300     for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
301     potential += bend->getPotential();
302     }
303    
304     for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
305     potential += torsion->getPotential();
306     }
307    
308     return potential;
309    
310     }
311    
312 tim 1695 std::ostream& operator <<(std::ostream& o, Molecule& mol) {
313 tim 1692 o << std::endl;
314     o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
315     o << mol.getNAtoms() << " atoms" << std::endl;
316     o << mol.getNBonds() << " bonds" << std::endl;
317     o << mol.getNBends() << " bends" << std::endl;
318     o << mol.getNTorsions() << " torsions" << std::endl;
319     o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
320     o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
321     o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
322 tim 1901 o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
323 tim 1692 return o;
324     }
325 gezelter 1490
326 tim 1692 }//end namespace oopse