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"
#include "utils/simError.h"

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_);
    MemoryUtils::deleteVectorOfPointer(constraintPairs_);
    MemoryUtils::deleteVectorOfPointer(constraintElems_);
    //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::addConstraintPair(ConstraintPair* cp) {
    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) {
        constraintPairs_.push_back(cp);
    }

}

void Molecule::addConstraintElem(ConstraintElem* cp) {
    if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == constraintElems_.end()) {
        constraintElems_.push_back(cp);
    }

}

void Molecule::complete() {
    
    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());
    }

    Atom* atom;
    AtomIterator ai;
    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
   
        if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
            //if an atom does not belong to a rigid body, it is an integrable object
            integrableObjects_.push_back(*ai);
        }
    }

    //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
    //    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());
}

double Molecule::getMass() {
    StuntDouble* sd;
    std::vector<StuntDouble*>::iterator i;
    double mass = 0.0;

    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
        mass += sd->getMass();
    }

    return mass;

}

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;
}

double Molecule::getPotential() {

    Bond* bond;
    Bend* bend;
    Torsion* torsion;
    Molecule::BondIterator bondIter;;
    Molecule::BendIterator  bendIter;
    Molecule::TorsionIterator  torsionIter;

    double potential = 0.0;

    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;

}

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;
    o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
    return o;
}

}//end namespace oopse
Revision:	1907
Committed:	Thu Jan 6 22:31:07 2005 UTC (20 years, 6 months ago) by tim
File size:	7676 byte(s)
Log Message:	constraint is almost working
#	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	tim	1907	MemoryUtils::deleteVectorOfPointer(constraintPairs_);
55			MemoryUtils::deleteVectorOfPointer(constraintElems_);
56	tim	1695	//integrableObjects_ don't own the objects
57	tim	1692	integrableObjects_.clear();
58
59	gezelter	1490	}
60
61	tim	1692	void Molecule::addAtom(Atom* atom) {
62	tim	1695	if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
63	tim	1692	atoms_.push_back(atom);
64			}
65			}
66	gezelter	1490
67	tim	1692	void Molecule::addBond(Bond* bond) {
68	tim	1695	if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
69	tim	1692	bonds_.push_back(bond);
70			}
71			}
72	gezelter	1490
73	tim	1692	void Molecule::addBend(Bend* bend) {
74	tim	1695	if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
75	tim	1692	bends_.push_back(bend);
76			}
77			}
78	gezelter	1490
79	tim	1692	void Molecule::addTorsion(Torsion* torsion) {
80	tim	1695	if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
81	tim	1692	torsions_.push_back(torsion);
82			}
83			}
84	gezelter	1490
85	tim	1692	void Molecule::addRigidBody(RigidBody *rb) {
86	tim	1695	if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
87	tim	1692	rigidBodies_.push_back(rb);
88			}
89	gezelter	1490	}
90
91	tim	1692	void Molecule::addCutoffGroup(CutoffGroup* cp) {
92	tim	1695	if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
93	tim	1692	cutoffGroups_.push_back(cp);
94			}
95	gezelter	1490
96	tim	1692	}
97	gezelter	1490
98	tim	1901	void Molecule::addConstraintPair(ConstraintPair* cp) {
99			if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) {
100			constraintPairs_.push_back(cp);
101			}
102
103			}
104
105	tim	1907	void Molecule::addConstraintElem(ConstraintElem* cp) {
106			if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == constraintElems_.end()) {
107			constraintElems_.push_back(cp);
108			}
109	tim	1901
110	tim	1907	}
111
112	tim	1692	void Molecule::complete() {
113
114			std::set<Atom*> rigidAtoms;
115			RigidBody* rb;
116	tim	1695	std::vector<RigidBody*>::iterator rbIter;
117	gezelter	1490
118	tim	1692
119			for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
120	tim	1695	rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
121	tim	1692	}
122	gezelter	1490
123	tim	1843	Atom* atom;
124			AtomIterator ai;
125			for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
126	tim	1844
127			if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
128			//if an atom does not belong to a rigid body, it is an integrable object
129	tim	1843	integrableObjects_.push_back(*ai);
130			}
131			}
132
133	tim	1692	//find all free atoms (which do not belong to rigid bodies)
134			//performs the "difference" operation from set theory, the output range contains a copy of every
135			//element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
136			//[rigidAtoms.begin(), rigidAtoms.end()).
137	tim	1843	//std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
138			// std::back_inserter(integrableObjects_));
139	gezelter	1490
140	tim	1843	//if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
141			// //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
142			// sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
143			//
144			// painCave.isFatal = 1;
145			// simError();
146			//}
147	tim	1733
148	tim	1692	integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
149	gezelter	1490	}
150
151	tim	1843	double Molecule::getMass() {
152			StuntDouble* sd;
153			std::vector<StuntDouble*>::iterator i;
154			double mass = 0.0;
155	gezelter	1490
156	tim	1843	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
157			mass += sd->getMass();
158			}
159
160			return mass;
161
162			}
163
164	tim	1692	Vector3d Molecule::getCom() {
165			StuntDouble* sd;
166			std::vector<StuntDouble*>::iterator i;
167			Vector3d com;
168			double totalMass = 0;
169			double mass;
170
171			for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
172			mass = sd->getMass();
173			totalMass += mass;
174			com += sd->getPos() * mass;
175			}
176	gezelter	1490
177	tim	1692	com /= totalMass;
178	gezelter	1490
179	tim	1692	return com;
180			}
181	gezelter	1490
182	tim	1695	void Molecule::moveCom(const Vector3d& delta) {
183	tim	1692	StuntDouble* sd;
184			std::vector<StuntDouble*>::iterator i;
185
186			for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
187	tim	1695	sd->setPos(sd->getPos() + delta);
188	gezelter	1490	}
189
190			}
191
192	tim	1692	Vector3d Molecule::getComVel() {
193			StuntDouble* sd;
194			std::vector<StuntDouble*>::iterator i;
195			Vector3d velCom;
196			double totalMass = 0;
197			double mass;
198
199			for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
200			mass = sd->getMass();
201			totalMass += mass;
202			velCom += sd->getVel() * mass;
203	gezelter	1490	}
204
205	tim	1692	velCom /= totalMass;
206	gezelter	1490
207	tim	1692	return velCom;
208	gezelter	1490	}
209
210	tim	1843	double Molecule::getPotential() {
211
212			Bond* bond;
213			Bend* bend;
214			Torsion* torsion;
215			Molecule::BondIterator bondIter;;
216			Molecule::BendIterator bendIter;
217			Molecule::TorsionIterator torsionIter;
218
219			double potential = 0.0;
220
221			for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
222			potential += bond->getPotential();
223			}
224
225			for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
226			potential += bend->getPotential();
227			}
228
229			for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
230			potential += torsion->getPotential();
231			}
232
233			return potential;
234
235			}
236
237	tim	1695	std::ostream& operator <<(std::ostream& o, Molecule& mol) {
238	tim	1692	o << std::endl;
239			o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
240			o << mol.getNAtoms() << " atoms" << std::endl;
241			o << mol.getNBonds() << " bonds" << std::endl;
242			o << mol.getNBends() << " bends" << std::endl;
243			o << mol.getNTorsions() << " torsions" << std::endl;
244			o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
245			o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
246			o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
247	tim	1901	o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
248	tim	1692	return o;
249			}
250	gezelter	1490
251	tim	1692	}//end namespace oopse