test/brains/Molecule.cpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */

/**
 * @file Molecule.cpp
 * @author    tlin
 * @date  10/28/2004
 * @version 1.0
 */ 

#include "primitives/Molecule.hpp"
#include <algorithm>

namespace oopse {

Molecule::~Molecule() {

    deleteVectorOfPointer(atoms_);
    deleteVectorOfPointer(bonds_);
    deleteVectorOfPointer(bends_);
    deleteVectorOfPointer(torsions_);
    deleteVectorOfPointer(rigidbodies_);
    deleteVectorOfPointer(cutoffGroups_);

    integrableObjects_.clear();
    
}

void Molecule::addAtom(Atom* atom) {
    if (atoms_.find(atom) == atoms_.end()) {
        atoms_.push_back(atom);
    }
}

void Molecule::addBond(Bond* bond) {
    if (bonds_.find(bond) == bonds_.end()) {
        bonds_.push_back(bond);
    }
}

void Molecule::addBend(Bend* bend) {
    if (bends_.find(bend) == bends_.end()) {
        bends_.push_back(bend);
    }
}

void Molecule::addTorsion(Torsion* torsion) {
    if (torsions_.find(torsion) == torsions_.end()) {
        torsions_.push_back(torsion);
    }
}

void Molecule::addRigidBody(RigidBody *rb) {
    if (rigidBodies_.find(bond) == bonds_.end()) {
        rigidBodies_.push_back(rb);
    }
}

void Molecule::addCutoffGroup(CutoffGroup* cp) {
    if (cutoffGroups_.find(bond) == bonds_.end()) {
        cutoffGroups_.push_back(cp);
    }

}

void Molecule::complete() {
    
    std::set<Atom*> allAtoms;
    allAtoms.insert(atoms_.begin(), atoms_.end());

    std::set<Atom*> rigidAtoms;
    RigidBody* rb;
    std::vector<RigidBody*> rbIter;

    
    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
        rigidAtoms.insert(rb->beginAtomIter(), rb->endAtomIter());
    }

    //find all free atoms (which do not belong to rigid bodies)  
    //performs the "difference" operation from set theory,  the output range contains a copy of every
    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in 
    //[rigidAtoms.begin(), rigidAtoms.end()).
    std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
                            std::back_inserter(integrableObjects_));

    if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
        //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
    }

    integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
}

Atom* Molecule::beginAtom(std::vector<Atom*>::iterator& i) {
    i = atoms_.begin();
    return (i == atoms_.end()) ? NULL : *i;
}

Atom* Molecule::nextAtom(std::vector<Atom*>::iterator& i) {
    ++i;
    return (i == atoms_.end()) ? NULL : *i;    
}

Bond* Molecule::beginBond(std::vector<Bond*>::iterator& i) {
    i = bonds_.begin();
    return (i == bonds_.end()) ? NULL : *i;
}

Bond* Molecule::nextBond(std::vector<Bond*>::iterator& i) {
    ++i;
    return (i == bonds_.end()) ? NULL : *i;    

}


Bend* Molecule::beginBend(std::vector<Bend*>::iterator& i) {
    i = bends_.begin();
    return (i == bends_.end()) ? NULL : *i;
}

Bend* Molecule::nextBend(std::vector<Bend*>::iterator& i) {
    ++i;
    return (i == bends_.end()) ? NULL : *i;    
}

Torsion* Molecule::beginTorsion(std::vector<Torsion*>::iterator& i) {
    i = torsions_.begin();
    return (i == torsions_.end()) ? NULL : *i;
}

Torsion* Molecule::nextTorsion(std::vector<Torsion*>::iterator& i) {
    ++i;
    return (i == torsions_.end()) ? NULL : *i;    
}    

RigidBody* Molecule::beginRigidBody(std::vector<RigidBody*>::iterator& i) {
    i = rigidBodies_.begin();
    return (i == rigidBodies_.end()) ? NULL : *i;
}

RigidBody* Molecule::nextRigidBody(std::vector<RigidBody*>::iterator& i) {
    ++i;
    return (i == rigidBodies_.end()) ? NULL : *i;    
}

StuntDouble* Molecule::beginIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
    i = integrableObjects_.begin();
    return (i == integrableObjects_.end()) ? NULL : *i;
}

StuntDouble* Molecule::nextIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
    ++i;
    return (i == integrableObjects_.end()) ? NULL : *i;    
}    

CutoffGroup* Molecule::beginCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {
    i = cutoffGroups_.begin();
    return (i == cutoffGroups_.end()) ? NULL : *i;
}

CutoffGroup* Molecule::nextCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {            
    ++i;
    return (i == cutoffGroups_.end()) ? NULL : *i;    
} 

void Molecule::calcForces() {
    RigidBody* rb;
    Bond* bond;
    Bend* bend;
    Torsion* torsion;
    std::vector<RigidBody*> rbIter;
    std::vector<Bond*> bondIter;;
    std::vector<Bend*> bendIter;
    std::vector<Torsion*> torsionIter;

    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
        rb->updateAtoms();
    }

    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
        bond->calcForce();
    }

    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
        bend->calcForce();
    }

    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
        torsion->calcForce();
    }
    
}

double Molecule::getPotential() {
    //RigidBody* rb;
    Bond* bond;
    Bend* bend;
    Torsion* torsion;
    //std::vector<RigidBody*> rbIter;
    std::vector<Bond*> bondIter;;
    std::vector<Bend*> bendIter;
    std::vector<Torsion*> torsionIter;

    double potential = 0;
    
    //for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
    //    rb->updateAtoms();
    //}

    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
        potential += bond->getPotential();
    }

    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
        potential += bend->getPotential();
    }

    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
        potential += torsion->getPotential();
    }
    
    return potential;
}

Vector3d Molecule::getCom() {
    StuntDouble* sd;
    std::vector<StuntDouble*>::iterator i;
    Vector3d com;
    double totalMass = 0;
    double mass;
    
    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
        mass = sd->getMass();
        totalMass += mass;
        com += sd->getPos() * mass;    
    }

    com /= totalMass;

    return com;
}

void Molecule::moveCom(const Vetor3d& delta) {
    StuntDouble* sd;
    std::vector<StuntDouble*>::iterator i;
    
    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
        s->setPos(sd->getPos() + delta);
    }

}

Vector3d Molecule::getComVel() {
    StuntDouble* sd;
    std::vector<StuntDouble*>::iterator i;
    Vector3d velCom;
    double totalMass = 0;
    double mass;
    
    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
        mass = sd->getMass();
        totalMass += mass;
        velCom += sd->getVel() * mass;    
    }

    velCom /= totalMass;

    return velCom;
}

std::ostream& operator <<(std::ostream& o, const Molecule& mol) {
    o << std::endl;
    o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
    o << mol.getNAtoms() << " atoms" << std::endl;
    o << mol.getNBonds() << " bonds" << std::endl;
    o << mol.getNBends() << " bends" << std::endl;
    o << mol.getNTorsions() << " torsions" << std::endl;
    o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
    o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
    o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;

    return o;
}

}//end namespace oopse
Revision:	1690
Committed:	Fri Oct 29 22:37:59 2004 UTC (19 years, 8 months ago) by tim
File size:	8755 byte(s)
Log Message:	more work on Molecule class
#	User	Rev	Content
1	tim	1686	/*
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
26			/**
27	tim	1690	* @file Molecule.cpp
28	tim	1686	* @author tlin
29	tim	1690	* @date 10/28/2004
30	tim	1686	* @version 1.0
31			*/
32
33			#include "primitives/Molecule.hpp"
34	tim	1690	#include <algorithm>
35	tim	1686
36			namespace oopse {
37
38			Molecule::~Molecule() {
39
40			deleteVectorOfPointer(atoms_);
41			deleteVectorOfPointer(bonds_);
42			deleteVectorOfPointer(bends_);
43			deleteVectorOfPointer(torsions_);
44			deleteVectorOfPointer(rigidbodies_);
45			deleteVectorOfPointer(cutoffGroups_);
46
47			integrableObjects_.clear();
48
49			}
50
51	tim	1690	void Molecule::addAtom(Atom* atom) {
52			if (atoms_.find(atom) == atoms_.end()) {
53			atoms_.push_back(atom);
54			}
55			}
56	tim	1686
57	tim	1690	void Molecule::addBond(Bond* bond) {
58			if (bonds_.find(bond) == bonds_.end()) {
59			bonds_.push_back(bond);
60			}
61			}
62
63			void Molecule::addBend(Bend* bend) {
64			if (bends_.find(bend) == bends_.end()) {
65			bends_.push_back(bend);
66			}
67			}
68
69			void Molecule::addTorsion(Torsion* torsion) {
70			if (torsions_.find(torsion) == torsions_.end()) {
71			torsions_.push_back(torsion);
72			}
73			}
74
75			void Molecule::addRigidBody(RigidBody *rb) {
76			if (rigidBodies_.find(bond) == bonds_.end()) {
77			rigidBodies_.push_back(rb);
78			}
79			}
80
81			void Molecule::addCutoffGroup(CutoffGroup* cp) {
82			if (cutoffGroups_.find(bond) == bonds_.end()) {
83			cutoffGroups_.push_back(cp);
84			}
85
86			}
87
88			void Molecule::complete() {
89
90			std::set<Atom*> allAtoms;
91			allAtoms.insert(atoms_.begin(), atoms_.end());
92
93			std::set<Atom*> rigidAtoms;
94			RigidBody* rb;
95			std::vector<RigidBody*> rbIter;
96
97
98			for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
99			rigidAtoms.insert(rb->beginAtomIter(), rb->endAtomIter());
100			}
101
102			//find all free atoms (which do not belong to rigid bodies)
103			//performs the "difference" operation from set theory, the output range contains a copy of every
104			//element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
105			//[rigidAtoms.begin(), rigidAtoms.end()).
106			std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
107			std::back_inserter(integrableObjects_));
108
109			if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
110			//Some atoms in rigidAtoms are not in allAtoms, something must be wrong
111			}
112
113			integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
114			}
115
116	tim	1686	Atom* Molecule::beginAtom(std::vector<Atom*>::iterator& i) {
117			i = atoms_.begin();
118			return (i == atoms_.end()) ? NULL : *i;
119			}
120
121			Atom* Molecule::nextAtom(std::vector<Atom*>::iterator& i) {
122			++i;
123			return (i == atoms_.end()) ? NULL : *i;
124			}
125
126			Bond* Molecule::beginBond(std::vector<Bond*>::iterator& i) {
127			i = bonds_.begin();
128			return (i == bonds_.end()) ? NULL : *i;
129			}
130
131			Bond* Molecule::nextBond(std::vector<Bond*>::iterator& i) {
132			++i;
133			return (i == bonds_.end()) ? NULL : *i;
134
135			}
136
137
138			Bend* Molecule::beginBend(std::vector<Bend*>::iterator& i) {
139			i = bends_.begin();
140			return (i == bends_.end()) ? NULL : *i;
141			}
142
143			Bend* Molecule::nextBend(std::vector<Bend*>::iterator& i) {
144			++i;
145			return (i == bends_.end()) ? NULL : *i;
146			}
147
148			Torsion* Molecule::beginTorsion(std::vector<Torsion*>::iterator& i) {
149			i = torsions_.begin();
150			return (i == torsions_.end()) ? NULL : *i;
151			}
152
153			Torsion* Molecule::nextTorsion(std::vector<Torsion*>::iterator& i) {
154			++i;
155			return (i == torsions_.end()) ? NULL : *i;
156			}
157
158			RigidBody* Molecule::beginRigidBody(std::vector<RigidBody*>::iterator& i) {
159			i = rigidBodies_.begin();
160			return (i == rigidBodies_.end()) ? NULL : *i;
161			}
162
163			RigidBody* Molecule::nextRigidBody(std::vector<RigidBody*>::iterator& i) {
164			++i;
165			return (i == rigidBodies_.end()) ? NULL : *i;
166			}
167
168			StuntDouble* Molecule::beginIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
169			i = integrableObjects_.begin();
170			return (i == integrableObjects_.end()) ? NULL : *i;
171			}
172
173			StuntDouble* Molecule::nextIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
174			++i;
175			return (i == integrableObjects_.end()) ? NULL : *i;
176			}
177
178			CutoffGroup* Molecule::beginCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {
179			i = cutoffGroups_.begin();
180			return (i == cutoffGroups_.end()) ? NULL : *i;
181			}
182
183			CutoffGroup* Molecule::nextCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {
184			++i;
185			return (i == cutoffGroups_.end()) ? NULL : *i;
186			}
187
188			void Molecule::calcForces() {
189			RigidBody* rb;
190			Bond* bond;
191			Bend* bend;
192			Torsion* torsion;
193			std::vector<RigidBody*> rbIter;
194			std::vector<Bond*> bondIter;;
195			std::vector<Bend*> bendIter;
196			std::vector<Torsion*> torsionIter;
197
198			for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
199			rb->updateAtoms();
200			}
201
202			for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
203			bond->calcForce();
204			}
205
206			for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
207			bend->calcForce();
208			}
209
210			for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
211			torsion->calcForce();
212			}
213
214			}
215
216	tim	1690	double Molecule::getPotential() {
217			//RigidBody* rb;
218			Bond* bond;
219			Bend* bend;
220			Torsion* torsion;
221			//std::vector<RigidBody*> rbIter;
222			std::vector<Bond*> bondIter;;
223			std::vector<Bend*> bendIter;
224			std::vector<Torsion*> torsionIter;
225
226			double potential = 0;
227
228			//for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
229			// rb->updateAtoms();
230			//}
231
232			for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
233			potential += bond->getPotential();
234			}
235
236			for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
237			potential += bend->getPotential();
238			}
239
240			for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
241			potential += torsion->getPotential();
242			}
243
244			return potential;
245			}
246
247			Vector3d Molecule::getCom() {
248			StuntDouble* sd;
249			std::vector<StuntDouble*>::iterator i;
250			Vector3d com;
251			double totalMass = 0;
252			double mass;
253
254			for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
255			mass = sd->getMass();
256			totalMass += mass;
257			com += sd->getPos() * mass;
258			}
259
260			com /= totalMass;
261
262			return com;
263			}
264
265			void Molecule::moveCom(const Vetor3d& delta) {
266			StuntDouble* sd;
267			std::vector<StuntDouble*>::iterator i;
268
269			for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
270			s->setPos(sd->getPos() + delta);
271			}
272
273			}
274
275			Vector3d Molecule::getComVel() {
276			StuntDouble* sd;
277			std::vector<StuntDouble*>::iterator i;
278			Vector3d velCom;
279			double totalMass = 0;
280			double mass;
281
282			for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
283			mass = sd->getMass();
284			totalMass += mass;
285			velCom += sd->getVel() * mass;
286			}
287
288			velCom /= totalMass;
289
290			return velCom;
291			}
292
293			std::ostream& operator <<(std::ostream& o, const Molecule& mol) {
294			o << std::endl;
295			o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
296			o << mol.getNAtoms() << " atoms" << std::endl;
297			o << mol.getNBonds() << " bonds" << std::endl;
298			o << mol.getNBends() << " bends" << std::endl;
299			o << mol.getNTorsions() << " torsions" << std::endl;
300			o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
301			o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
302			o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
303
304			return o;
305			}
306
307	tim	1686	}//end namespace oopse