src/primitives/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 <algorithm>
#include <set>

#include "primitives/Molecule.hpp"
#include "utils/MemoryUtils.hpp"

namespace oopse {

Molecule::~Molecule() {

    MemoryUtils::deleteVectorOfPointer(atoms_);
    MemoryUtils::deleteVectorOfPointer(bonds_);
    MemoryUtils::deleteVectorOfPointer(bends_);
    MemoryUtils::deleteVectorOfPointer(torsions_);
    MemoryUtils::deleteVectorOfPointer(rigidBodies_);
    MemoryUtils::deleteVectorOfPointer(cutoffGroups_);

    //integrableObjects_ don't own the objects
    integrableObjects_.clear();
    
}

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

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

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

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

void Molecule::addRigidBody(RigidBody *rb) {
    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
        rigidBodies_.push_back(rb);
    }
}

void Molecule::addCutoffGroup(CutoffGroup* cp) {
    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.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*>::iterator rbIter;

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

    //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*>::iterator rbIter;
    std::vector<Bond*>::iterator bondIter;;
    std::vector<Bend*>::iterator bendIter;
    std::vector<Torsion*>::iterator torsionIter;

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

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

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

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

double Molecule::getPotential() {
    //RigidBody* rb;
    Bond* bond;
    Bend* bend;
    Torsion* torsion;
    //std::vector<RigidBody*> rbIter;
    std::vector<Bond*>::iterator bondIter;;
    std::vector<Bend*>::iterator bendIter;
    std::vector<Torsion*>::iterator 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->get_potential();
    }

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

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