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(int stampId, int globalIndex, const std::string& molName) 
    : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {

}

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_.end()));

    if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
        //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
        sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");

        painCave.isFatal = 1;
        simError();        
    }

    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:	1734
Committed:	Fri Nov 12 07:05:43 2004 UTC (19 years, 8 months ago) by tim
File size:	9538 byte(s)
Log Message:	MoleculeCreator forget to create CutoffGroups for free atoms
#	Content
1	/*
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	* @file Molecule.cpp
28	* @author tlin
29	* @date 10/28/2004
30	* @version 1.0
31	*/
32
33	#include <algorithm>
34	#include <set>
35
36	#include "primitives/Molecule.hpp"
37	#include "utils/MemoryUtils.hpp"
38
39	namespace oopse {
40	Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
41	: stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
42
43	}
44
45	Molecule::~Molecule() {
46
47	MemoryUtils::deleteVectorOfPointer(atoms_);
48	MemoryUtils::deleteVectorOfPointer(bonds_);
49	MemoryUtils::deleteVectorOfPointer(bends_);
50	MemoryUtils::deleteVectorOfPointer(torsions_);
51	MemoryUtils::deleteVectorOfPointer(rigidBodies_);
52	MemoryUtils::deleteVectorOfPointer(cutoffGroups_);
53
54	//integrableObjects_ don't own the objects
55	integrableObjects_.clear();
56
57	}
58
59	void Molecule::addAtom(Atom* atom) {
60	if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
61	atoms_.push_back(atom);
62	}
63	}
64
65	void Molecule::addBond(Bond* bond) {
66	if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
67	bonds_.push_back(bond);
68	}
69	}
70
71	void Molecule::addBend(Bend* bend) {
72	if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
73	bends_.push_back(bend);
74	}
75	}
76
77	void Molecule::addTorsion(Torsion* torsion) {
78	if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
79	torsions_.push_back(torsion);
80	}
81	}
82
83	void Molecule::addRigidBody(RigidBody *rb) {
84	if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
85	rigidBodies_.push_back(rb);
86	}
87	}
88
89	void Molecule::addCutoffGroup(CutoffGroup* cp) {
90	if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
91	cutoffGroups_.push_back(cp);
92	}
93
94	}
95
96	void Molecule::complete() {
97
98	std::set<Atom*> allAtoms;
99	allAtoms.insert(atoms_.begin(), atoms_.end());
100
101	std::set<Atom*> rigidAtoms;
102	RigidBody* rb;
103	std::vector<RigidBody*>::iterator rbIter;
104
105
106	for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
107	rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
108	}
109
110	//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_.end()));
116
117	if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
118	//Some atoms in rigidAtoms are not in allAtoms, something must be wrong
119	sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
120
121	painCave.isFatal = 1;
122	simError();
123	}
124
125	integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
126	}
127
128	Atom* Molecule::beginAtom(std::vector<Atom*>::iterator& i) {
129	i = atoms_.begin();
130	return (i == atoms_.end()) ? NULL : *i;
131	}
132
133	Atom* Molecule::nextAtom(std::vector<Atom*>::iterator& i) {
134	++i;
135	return (i == atoms_.end()) ? NULL : *i;
136	}
137
138	Bond* Molecule::beginBond(std::vector<Bond*>::iterator& i) {
139	i = bonds_.begin();
140	return (i == bonds_.end()) ? NULL : *i;
141	}
142
143	Bond* Molecule::nextBond(std::vector<Bond*>::iterator& i) {
144	++i;
145	return (i == bonds_.end()) ? NULL : *i;
146
147	}
148
149
150	Bend* Molecule::beginBend(std::vector<Bend*>::iterator& i) {
151	i = bends_.begin();
152	return (i == bends_.end()) ? NULL : *i;
153	}
154
155	Bend* Molecule::nextBend(std::vector<Bend*>::iterator& i) {
156	++i;
157	return (i == bends_.end()) ? NULL : *i;
158	}
159
160	Torsion* Molecule::beginTorsion(std::vector<Torsion*>::iterator& i) {
161	i = torsions_.begin();
162	return (i == torsions_.end()) ? NULL : *i;
163	}
164
165	Torsion* Molecule::nextTorsion(std::vector<Torsion*>::iterator& i) {
166	++i;
167	return (i == torsions_.end()) ? NULL : *i;
168	}
169
170	RigidBody* Molecule::beginRigidBody(std::vector<RigidBody*>::iterator& i) {
171	i = rigidBodies_.begin();
172	return (i == rigidBodies_.end()) ? NULL : *i;
173	}
174
175	RigidBody* Molecule::nextRigidBody(std::vector<RigidBody*>::iterator& i) {
176	++i;
177	return (i == rigidBodies_.end()) ? NULL : *i;
178	}
179
180	StuntDouble* Molecule::beginIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
181	i = integrableObjects_.begin();
182	return (i == integrableObjects_.end()) ? NULL : *i;
183	}
184
185	StuntDouble* Molecule::nextIntegrableObject(std::vector<StuntDouble*>::iterator& i) {
186	++i;
187	return (i == integrableObjects_.end()) ? NULL : *i;
188	}
189
190	CutoffGroup* Molecule::beginCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {
191	i = cutoffGroups_.begin();
192	return (i == cutoffGroups_.end()) ? NULL : *i;
193	}
194
195	CutoffGroup* Molecule::nextCutoffGroup(std::vector<CutoffGroup*>::iterator& i) {
196	++i;
197	return (i == cutoffGroups_.end()) ? NULL : *i;
198	}
199
200	void Molecule::calcForces() {
201	RigidBody* rb;
202	Bond* bond;
203	Bend* bend;
204	Torsion* torsion;
205	std::vector<RigidBody*>::iterator rbIter;
206	std::vector<Bond*>::iterator bondIter;;
207	std::vector<Bend*>::iterator bendIter;
208	std::vector<Torsion*>::iterator torsionIter;
209
210	for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
211	rb->updateAtoms();
212	}
213
214	for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
215	bond->calc_forces();
216	}
217
218	for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
219	bend->calc_forces();
220	}
221
222	for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
223	torsion->calc_forces();
224	}
225
226	}
227
228	double Molecule::getPotential() {
229	//RigidBody* rb;
230	Bond* bond;
231	Bend* bend;
232	Torsion* torsion;
233	//std::vector<RigidBody*> rbIter;
234	std::vector<Bond*>::iterator bondIter;;
235	std::vector<Bend*>::iterator bendIter;
236	std::vector<Torsion*>::iterator torsionIter;
237
238	double potential = 0;
239
240	//for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
241	// rb->updateAtoms();
242	//}
243
244	for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
245	potential += bond->get_potential();
246	}
247
248	for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
249	potential += bend->get_potential();
250	}
251
252	for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
253	potential += torsion->get_potential();
254	}
255
256	return potential;
257	}
258
259	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
272	com /= totalMass;
273
274	return com;
275	}
276
277	void Molecule::moveCom(const Vector3d& delta) {
278	StuntDouble* sd;
279	std::vector<StuntDouble*>::iterator i;
280
281	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
282	sd->setPos(sd->getPos() + delta);
283	}
284
285	}
286
287	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	}
299
300	velCom /= totalMass;
301
302	return velCom;
303	}
304
305	std::ostream& operator <<(std::ostream& o, Molecule& mol) {
306	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
316	return o;
317	}
318
319	}//end namespace oopse