ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-2.0/src/primitives/Molecule.cpp
Revision: 1733
Committed: Fri Nov 12 06:19:04 2004 UTC (19 years, 7 months ago) by tim
File size: 9532 byte(s)
Log Message:
OOPSE = Object-Obfuscated Parallel Simulation Engine

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    
39 tim 1692 namespace oopse {
40 tim 1733 Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
41     : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
42 gezelter 1490
43 tim 1733 }
44    
45 tim 1692 Molecule::~Molecule() {
46 gezelter 1490
47 tim 1695 MemoryUtils::deleteVectorOfPointer(atoms_);
48     MemoryUtils::deleteVectorOfPointer(bonds_);
49     MemoryUtils::deleteVectorOfPointer(bends_);
50     MemoryUtils::deleteVectorOfPointer(torsions_);
51     MemoryUtils::deleteVectorOfPointer(rigidBodies_);
52     MemoryUtils::deleteVectorOfPointer(cutoffGroups_);
53 gezelter 1490
54 tim 1695 //integrableObjects_ don't own the objects
55 tim 1692 integrableObjects_.clear();
56    
57 gezelter 1490 }
58    
59 tim 1692 void Molecule::addAtom(Atom* atom) {
60 tim 1695 if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
61 tim 1692 atoms_.push_back(atom);
62     }
63     }
64 gezelter 1490
65 tim 1692 void Molecule::addBond(Bond* bond) {
66 tim 1695 if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
67 tim 1692 bonds_.push_back(bond);
68     }
69     }
70 gezelter 1490
71 tim 1692 void Molecule::addBend(Bend* bend) {
72 tim 1695 if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
73 tim 1692 bends_.push_back(bend);
74     }
75     }
76 gezelter 1490
77 tim 1692 void Molecule::addTorsion(Torsion* torsion) {
78 tim 1695 if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
79 tim 1692 torsions_.push_back(torsion);
80     }
81     }
82 gezelter 1490
83 tim 1692 void Molecule::addRigidBody(RigidBody *rb) {
84 tim 1695 if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
85 tim 1692 rigidBodies_.push_back(rb);
86     }
87 gezelter 1490 }
88    
89 tim 1692 void Molecule::addCutoffGroup(CutoffGroup* cp) {
90 tim 1695 if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
91 tim 1692 cutoffGroups_.push_back(cp);
92     }
93 gezelter 1490
94 tim 1692 }
95 gezelter 1490
96 tim 1692 void Molecule::complete() {
97    
98     std::set<Atom*> allAtoms;
99     allAtoms.insert(atoms_.begin(), atoms_.end());
100 gezelter 1490
101 tim 1692 std::set<Atom*> rigidAtoms;
102     RigidBody* rb;
103 tim 1695 std::vector<RigidBody*>::iterator rbIter;
104 gezelter 1490
105 tim 1692
106     for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
107 tim 1695 rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
108 tim 1692 }
109 gezelter 1490
110 tim 1692 //find all free atoms (which do not belong to rigid bodies)
111     //performs the "difference" operation from set theory, the output range contains a copy of every
112     //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
113     //[rigidAtoms.begin(), rigidAtoms.end()).
114     std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
115     std::back_inserter(integrableObjects_));
116 gezelter 1490
117 tim 1692 if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
118     //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
119 tim 1733 sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
120    
121     painCave.isFatal = 1;
122     simError();
123 tim 1692 }
124 gezelter 1490
125 tim 1692 integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
126 gezelter 1490 }
127    
128 tim 1692 Atom* Molecule::beginAtom(std::vector<Atom*>::iterator& i) {
129     i = atoms_.begin();
130     return (i == atoms_.end()) ? NULL : *i;
131     }
132 gezelter 1490
133 tim 1692 Atom* Molecule::nextAtom(std::vector<Atom*>::iterator& i) {
134     ++i;
135     return (i == atoms_.end()) ? NULL : *i;
136     }
137 gezelter 1490
138 tim 1692 Bond* Molecule::beginBond(std::vector<Bond*>::iterator& i) {
139     i = bonds_.begin();
140     return (i == bonds_.end()) ? NULL : *i;
141     }
142 gezelter 1490
143 tim 1692 Bond* Molecule::nextBond(std::vector<Bond*>::iterator& i) {
144     ++i;
145     return (i == bonds_.end()) ? NULL : *i;
146 gezelter 1490
147 tim 1692 }
148 gezelter 1490
149    
150 tim 1692 Bend* Molecule::beginBend(std::vector<Bend*>::iterator& i) {
151     i = bends_.begin();
152     return (i == bends_.end()) ? NULL : *i;
153 gezelter 1490 }
154    
155 tim 1692 Bend* Molecule::nextBend(std::vector<Bend*>::iterator& i) {
156     ++i;
157     return (i == bends_.end()) ? NULL : *i;
158     }
159 gezelter 1490
160 tim 1692 Torsion* Molecule::beginTorsion(std::vector<Torsion*>::iterator& i) {
161     i = torsions_.begin();
162     return (i == torsions_.end()) ? NULL : *i;
163     }
164 gezelter 1490
165 tim 1692 Torsion* Molecule::nextTorsion(std::vector<Torsion*>::iterator& i) {
166     ++i;
167     return (i == torsions_.end()) ? NULL : *i;
168     }
169 gezelter 1490
170 tim 1692 RigidBody* Molecule::beginRigidBody(std::vector<RigidBody*>::iterator& i) {
171     i = rigidBodies_.begin();
172     return (i == rigidBodies_.end()) ? NULL : *i;
173     }
174 gezelter 1490
175 tim 1692 RigidBody* Molecule::nextRigidBody(std::vector<RigidBody*>::iterator& i) {
176     ++i;
177     return (i == rigidBodies_.end()) ? NULL : *i;
178     }
179 gezelter 1490
180 tim 1692 StuntDouble* Molecule::beginIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
181     i = integrableObjects_.begin();
182     return (i == integrableObjects_.end()) ? NULL : *i;
183 gezelter 1490 }
184    
185 tim 1692 StuntDouble* Molecule::nextIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
186     ++i;
187     return (i == integrableObjects_.end()) ? NULL : *i;
188     }
189 gezelter 1490
190 tim 1692 CutoffGroup* Molecule::beginCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {
191     i = cutoffGroups_.begin();
192     return (i == cutoffGroups_.end()) ? NULL : *i;
193 gezelter 1490 }
194    
195 tim 1692 CutoffGroup* Molecule::nextCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {
196     ++i;
197     return (i == cutoffGroups_.end()) ? NULL : *i;
198     }
199 gezelter 1490
200 tim 1692 void Molecule::calcForces() {
201     RigidBody* rb;
202     Bond* bond;
203     Bend* bend;
204     Torsion* torsion;
205 tim 1695 std::vector<RigidBody*>::iterator rbIter;
206     std::vector<Bond*>::iterator bondIter;;
207     std::vector<Bend*>::iterator bendIter;
208     std::vector<Torsion*>::iterator torsionIter;
209 gezelter 1490
210 tim 1692 for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
211     rb->updateAtoms();
212 gezelter 1490 }
213    
214 tim 1692 for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
215 tim 1701 bond->calc_forces();
216 tim 1692 }
217 gezelter 1490
218 tim 1692 for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
219 tim 1701 bend->calc_forces();
220 tim 1692 }
221 gezelter 1490
222 tim 1692 for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
223 tim 1701 torsion->calc_forces();
224 gezelter 1490 }
225 tim 1692
226 gezelter 1490 }
227    
228 tim 1692 double Molecule::getPotential() {
229     //RigidBody* rb;
230     Bond* bond;
231     Bend* bend;
232     Torsion* torsion;
233     //std::vector<RigidBody*> rbIter;
234 tim 1695 std::vector<Bond*>::iterator bondIter;;
235     std::vector<Bend*>::iterator bendIter;
236     std::vector<Torsion*>::iterator torsionIter;
237 gezelter 1490
238 tim 1692 double potential = 0;
239    
240     //for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
241     // rb->updateAtoms();
242     //}
243 gezelter 1490
244 tim 1692 for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
245 tim 1701 potential += bond->get_potential();
246 tim 1692 }
247 gezelter 1490
248 tim 1692 for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
249 tim 1701 potential += bend->get_potential();
250 tim 1692 }
251 gezelter 1490
252 tim 1692 for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
253 tim 1701 potential += torsion->get_potential();
254 tim 1692 }
255    
256     return potential;
257     }
258 gezelter 1490
259 tim 1692 Vector3d Molecule::getCom() {
260     StuntDouble* sd;
261     std::vector<StuntDouble*>::iterator i;
262     Vector3d com;
263     double totalMass = 0;
264     double mass;
265    
266     for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
267     mass = sd->getMass();
268     totalMass += mass;
269     com += sd->getPos() * mass;
270     }
271 gezelter 1490
272 tim 1692 com /= totalMass;
273 gezelter 1490
274 tim 1692 return com;
275     }
276 gezelter 1490
277 tim 1695 void Molecule::moveCom(const Vector3d& delta) {
278 tim 1692 StuntDouble* sd;
279     std::vector<StuntDouble*>::iterator i;
280    
281     for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
282 tim 1695 sd->setPos(sd->getPos() + delta);
283 gezelter 1490 }
284    
285     }
286    
287 tim 1692 Vector3d Molecule::getComVel() {
288     StuntDouble* sd;
289     std::vector<StuntDouble*>::iterator i;
290     Vector3d velCom;
291     double totalMass = 0;
292     double mass;
293    
294     for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
295     mass = sd->getMass();
296     totalMass += mass;
297     velCom += sd->getVel() * mass;
298 gezelter 1490 }
299    
300 tim 1692 velCom /= totalMass;
301 gezelter 1490
302 tim 1692 return velCom;
303 gezelter 1490 }
304    
305 tim 1695 std::ostream& operator <<(std::ostream& o, Molecule& mol) {
306 tim 1692 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     o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
313     o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
314     o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
315 gezelter 1490
316 tim 1692 return o;
317     }
318 gezelter 1490
319 tim 1692 }//end namespace oopse