ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/primitives/Molecule.cpp

Comparing trunk/src/primitives/Molecule.cpp (file contents):
Revision 372 by tim, Mon Feb 21 15:28:56 2005 UTC vs.
Revision 2054 by gezelter, Tue Jan 13 20:14:57 2015 UTC

#	Line 1 | Line 1
1	<	/*
1	>	/*
2		* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3		*
4		* The University of Notre Dame grants you ("Licensee") a
#	Line 6 | Line 6
6		* redistribute this software in source and binary code form, provided
7		* that the following conditions are met:
8		*
9	<	* 1. Acknowledgement of the program authors must be made in any
10	<	*    publication of scientific results based in part on use of the
11	<	*    program.  An acceptable form of acknowledgement is citation of
12	<	*    the article in which the program was described (Matthew
13	<	*    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	<	*    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	<	*    Parallel Simulation Engine for Molecular Dynamics,"
16	<	*    J. Comput. Chem. 26, pp. 252-271 (2005))
17	<	*
18	<	* 2. Redistributions of source code must retain the above copyright
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
12	>	* 2. Redistributions in binary form must reproduce the above copyright
13		*    notice, this list of conditions and the following disclaimer in the
14		*    documentation and/or other materials provided with the
15		*    distribution.
#	Line 37 | Line 28
28		* arising out of the use of or inability to use software, even if the
29		* University of Notre Dame has been advised of the possibility of
30		* such damages.
31	+	*
32	+	* SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33	+	* research, please cite the appropriate papers when you publish your
34	+	* work.  Good starting points are:
35	+	*
36	+	* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	+	* [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	+	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	+	* [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	+	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42
43		/**
#	Line 52 | Line 53
53		#include "primitives/Molecule.hpp"
54		#include "utils/MemoryUtils.hpp"
55		#include "utils/simError.h"
56	+	#include "utils/StringUtils.hpp"
57
58	<	namespace oopse {
59	<	Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
60	<	: stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
58	>	namespace OpenMD {
59	>	Molecule::Molecule(int stampId, int globalIndex, const std::string& molName,
60	>	int region) : stampId_(stampId),
61	>	globalIndex_(globalIndex),
62	>	moleculeName_(molName),
63	>	region_(region),
64	>	constrainTotalCharge_(false) {
65	>	}
66	>
67	>	Molecule::~Molecule() {
68	>
69	>	MemoryUtils::deletePointers(atoms_);
70	>	MemoryUtils::deletePointers(bonds_);
71	>	MemoryUtils::deletePointers(bends_);
72	>	MemoryUtils::deletePointers(torsions_);
73	>	MemoryUtils::deletePointers(inversions_);
74	>	MemoryUtils::deletePointers(rigidBodies_);
75	>	MemoryUtils::deletePointers(cutoffGroups_);
76	>	MemoryUtils::deletePointers(constraintPairs_);
77	>	MemoryUtils::deletePointers(constraintElems_);
78
79	<	}
61	<
62	<	Molecule::~Molecule() {
63	<
64	<	MemoryUtils::deleteVectorOfPointer(atoms_);
65	<	MemoryUtils::deleteVectorOfPointer(bonds_);
66	<	MemoryUtils::deleteVectorOfPointer(bends_);
67	<	MemoryUtils::deleteVectorOfPointer(torsions_);
68	<	MemoryUtils::deleteVectorOfPointer(rigidBodies_);
69	<	MemoryUtils::deleteVectorOfPointer(cutoffGroups_);
70	<	MemoryUtils::deleteVectorOfPointer(constraintPairs_);
71	<	MemoryUtils::deleteVectorOfPointer(constraintElems_);
72	<	//integrableObjects_ don't own the objects
79	>	// integrableObjects_ don't own the objects
80		integrableObjects_.clear();
81	+	fluctuatingCharges_.clear();
82
83	<	}
84	<
85	<	void Molecule::addAtom(Atom* atom) {
83	>	}
84	>
85	>	void Molecule::addAtom(Atom* atom) {
86		if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
87	<	atoms_.push_back(atom);
87	>	atoms_.push_back(atom);
88		}
89	<	}
90	<
91	<	void Molecule::addBond(Bond* bond) {
89	>	}
90	>
91	>	void Molecule::addBond(Bond* bond) {
92		if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
93	<	bonds_.push_back(bond);
93	>	bonds_.push_back(bond);
94		}
95	<	}
96	<
97	<	void Molecule::addBend(Bend* bend) {
95	>	}
96	>
97	>	void Molecule::addBend(Bend* bend) {
98		if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
99	<	bends_.push_back(bend);
99	>	bends_.push_back(bend);
100		}
101	<	}
102	<
103	<	void Molecule::addTorsion(Torsion* torsion) {
104	<	if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
105	<	torsions_.push_back(torsion);
101	>	}
102	>
103	>	void Molecule::addTorsion(Torsion* torsion) {
104	>	if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
105	>	torsions_.end()) {
106	>	torsions_.push_back(torsion);
107		}
108	<	}
108	>	}
109
110	<	void Molecule::addRigidBody(RigidBody *rb) {
111	<	if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
112	<	rigidBodies_.push_back(rb);
110	>	void Molecule::addInversion(Inversion* inversion) {
111	>	if (std::find(inversions_.begin(), inversions_.end(), inversion) ==
112	>	inversions_.end()) {
113	>	inversions_.push_back(inversion);
114		}
115	<	}
116	<
117	<	void Molecule::addCutoffGroup(CutoffGroup* cp) {
118	<	if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
119	<	cutoffGroups_.push_back(cp);
115	>	}
116	>
117	>	void Molecule::addRigidBody(RigidBody *rb) {
118	>	if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
119	>	rigidBodies_.end()) {
120	>	rigidBodies_.push_back(rb);
121		}
122	<
123	<	}
124	<
125	<	void Molecule::addConstraintPair(ConstraintPair* cp) {
126	<	if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) {
127	<	constraintPairs_.push_back(cp);
122	>	}
123	>
124	>	void Molecule::addCutoffGroup(CutoffGroup* cp) {
125	>	if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
126	>	cutoffGroups_.end()) {
127	>	cutoffGroups_.push_back(cp);
128	>	}
129	>	}
130	>
131	>	void Molecule::addConstraintPair(ConstraintPair* cp) {
132	>	if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
133	>	constraintPairs_.end()) {
134	>	constraintPairs_.push_back(cp);
135	>	}
136	>	}
137	>
138	>	void Molecule::addConstraintElem(ConstraintElem* cp) {
139	>	if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
140	>	constraintElems_.end()) {
141	>	constraintElems_.push_back(cp);
142		}
143	<
144	<	}
145	<
121	<	void Molecule::addConstraintElem(ConstraintElem* cp) {
122	<	if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == constraintElems_.end()) {
123	<	constraintElems_.push_back(cp);
124	<	}
125	<
126	<	}
127	<
128	<	void Molecule::complete() {
143	>	}
144	>
145	>	void Molecule::complete() {
146
147		std::set<Atom*> rigidAtoms;
148	+	Atom* atom;
149	+	Atom* atom1;
150	+	Atom* atom2;
151	+	AtomIterator ai, aj;
152		RigidBody* rb;
153	<	std::vector<RigidBody*>::iterator rbIter;
153	>	RigidBodyIterator rbIter;
154	>	Bond* bond;
155	>	BondIterator bi;
156
157	<
157	>	// Get list of all the atoms that are part of rigid bodies
158	>
159		for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
160	<	rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
160	>	rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
161		}
162	+
163	+	// add any atom that wasn't part of a rigid body to the list of integrableObjects
164
139	–	Atom* atom;
140	–	AtomIterator ai;
165		for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
166	<
167	<	if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
144	<	//if an atom does not belong to a rigid body, it is an integrable object
145	<	integrableObjects_.push_back(*ai);
146	<	}
147	<	}
166	>
167	>	if (rigidAtoms.find(atom) == rigidAtoms.end()) {
168
169	<	//find all free atoms (which do not belong to rigid bodies)
170	<	//performs the "difference" operation from set theory,  the output range contains a copy of every
151	<	//element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
152	<	//[rigidAtoms.begin(), rigidAtoms.end()).
153	<	//std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
154	<	//                        std::back_inserter(integrableObjects_));
169	>	// If an atom does not belong to a rigid body, it is an
170	>	// integrable object
171
172	<	//if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
173	<	//    //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
174	<	//    sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
175	<	//
176	<	//    painCave.isFatal = 1;
177	<	//    simError();
162	<	//}
172	>	integrableObjects_.push_back(atom);
173	>	}
174	>	}
175	>
176	>	// then add the rigid bodies themselves to the integrableObjects
177	>
178		for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
179		integrableObjects_.push_back(rb);
180		}
166	–	//integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
167	–	}
181
182	<	double Molecule::getMass() {
183	<	StuntDouble* sd;
171	<	std::vector<StuntDouble*>::iterator i;
172	<	double mass = 0.0;
182	>	// find the atoms that are fluctuating charges and add them to the
183	>	// fluctuatingCharges_ vector
184
185	<	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
186	<	mass += sd->getMass();
185	>	for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
186	>	if ( atom->isFluctuatingCharge() )
187	>	fluctuatingCharges_.push_back( atom );
188		}
189
178	–	return mass;
190
191	<	}
191	>	// find the electronegative atoms and add them to the
192	>	// hBondAcceptors_ vector:
193	>
194	>	for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
195	>	AtomType* at = atom->getAtomType();
196	>	// get the chain of base types for this atom type:
197	>	std::vector<AtomType*> ayb = at->allYourBase();
198	>	// use the last type in the chain of base types for the name:
199	>	std::string bn = UpperCase(ayb[ayb.size()-1]->getName());
200
201	<	Vector3d Molecule::getCom() {
201	>	if (bn.compare("O")==0 || bn.compare("N")==0
202	>	|| bn.compare("F")==0)
203	>	hBondAcceptors_.push_back( atom );
204	>
205	>	}
206	>
207	>	// find electronegative atoms that are either bonded to
208	>	// hydrogens or are present in the same rigid bodies:
209	>
210	>	for (bond = beginBond(bi); bond != NULL; bond = nextBond(bi)) {
211	>	Atom* atom1 = bond->getAtomA();
212	>	Atom* atom2 = bond->getAtomB();
213	>	AtomType* at1 = atom1->getAtomType();
214	>	AtomType* at2 = atom1->getAtomType();
215	>	// get the chain of base types for this atom type:
216	>	std::vector<AtomType*> ayb1 = at1->allYourBase();
217	>	std::vector<AtomType*> ayb2 = at2->allYourBase();
218	>	// use the last type in the chain of base types for the name:
219	>	std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName());
220	>	std::string bn2 = UpperCase(ayb2[ayb2.size()-1]->getName());
221	>
222	>	if (bn1.compare("H")==0) {
223	>	if (bn2.compare("O")==0 || bn2.compare("N")==0
224	>	|| bn2.compare("F")==0) {
225	>	HBondDonor* donor = new HBondDonor();
226	>	donor->donorAtom = atom2;
227	>	donor->donatedHydrogen = atom1;
228	>	hBondDonors_.push_back( donor );
229	>	}
230	>	}
231	>	if (bn2.compare("H")==0) {
232	>	if (bn1.compare("O")==0 || bn1.compare("N")==0
233	>	|| bn1.compare("F")==0) {
234	>	HBondDonor* donor = new HBondDonor();
235	>	donor->donorAtom = atom1;
236	>	donor->donatedHydrogen = atom2;
237	>	hBondDonors_.push_back( donor );
238	>	}
239	>	}
240	>	}
241	>
242	>	for (rb = beginRigidBody(rbIter); rb != NULL;
243	>	rb = nextRigidBody(rbIter)) {
244	>	for(atom1 = rb->beginAtom(ai); atom1 != NULL;
245	>	atom1 = rb->nextAtom(ai)) {
246	>	AtomType* at1 = atom1->getAtomType();
247	>	// get the chain of base types for this atom type:
248	>	std::vector<AtomType*> ayb1 = at1->allYourBase();
249	>	// use the last type in the chain of base types for the name:
250	>	std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName());
251	>
252	>	if (bn1.compare("O")==0 || bn1.compare("N")==0
253	>	|| bn1.compare("F")==0) {
254	>	for(atom2 = rb->beginAtom(aj); atom2 != NULL;
255	>	atom2 = rb->nextAtom(aj)) {
256	>	AtomType* at2 = atom2->getAtomType();
257	>	// get the chain of base types for this atom type:
258	>	std::vector<AtomType*> ayb2 = at2->allYourBase();
259	>	// use the last type in the chain of base types for the name:
260	>	std::string bn2 = UpperCase(ayb2[ayb2.size()-1]->getName());
261	>	if (bn2.compare("H")==0) {
262	>	HBondDonor* donor = new HBondDonor();
263	>	donor->donorAtom = atom1;
264	>	donor->donatedHydrogen = atom2;
265	>	hBondDonors_.push_back( donor );
266	>	}
267	>	}
268	>	}
269	>	}
270	>	}
271	>	}
272	>
273	>	RealType Molecule::getMass() {
274		StuntDouble* sd;
275		std::vector<StuntDouble*>::iterator i;
276	<	Vector3d com;
186	<	double totalMass = 0;
187	<	double mass;
276	>	RealType mass = 0.0;
277
278	<	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
279	<	mass = sd->getMass();
280	<	totalMass += mass;
192	<	com += sd->getPos() * mass;
278	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
279	>	nextIntegrableObject(i)){
280	>	mass += sd->getMass();
281		}
282	+
283	+	return mass;
284	+	}
285
286	+	Vector3d Molecule::getCom() {
287	+	StuntDouble* sd;
288	+	std::vector<StuntDouble*>::iterator i;
289	+	Vector3d com;
290	+	RealType totalMass = 0;
291	+	RealType mass;
292	+
293	+	for (sd = beginIntegrableObject(i); sd != NULL; sd =
294	+	nextIntegrableObject(i)){
295	+	mass = sd->getMass();
296	+	totalMass += mass;
297	+	com += sd->getPos() * mass;
298	+	}
299	+
300		com /= totalMass;
301
302		return com;
303	<	}
304	<
305	<	void Molecule::moveCom(const Vector3d& delta) {
303	>	}
304	>
305	>	Vector3d Molecule::getCom(int snapshotNo) {
306		StuntDouble* sd;
307		std::vector<StuntDouble*>::iterator i;
308	+	Vector3d com;
309	+	RealType totalMass = 0;
310	+	RealType mass;
311
312	<	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
313	<	sd->setPos(sd->getPos() + delta);
312	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
313	>	nextIntegrableObject(i)){
314	>	mass = sd->getMass();
315	>	totalMass += mass;
316	>	com += sd->getPos(snapshotNo) * mass;
317		}
318	+
319	+	com /= totalMass;
320
321	<	}
321	>	return com;
322	>	}
323
324	<	Vector3d Molecule::getComVel() {
324	>	void Molecule::moveCom(const Vector3d& delta) {
325		StuntDouble* sd;
326		std::vector<StuntDouble*>::iterator i;
327	+
328	+	for (sd = beginIntegrableObject(i); sd != NULL; sd =
329	+	nextIntegrableObject(i)){
330	+	sd->setPos(sd->getPos() + delta);
331	+	}
332	+	}
333	+
334	+	Vector3d Molecule::getComVel() {
335	+	StuntDouble* sd;
336	+	std::vector<StuntDouble*>::iterator i;
337		Vector3d velCom;
338	<	double totalMass = 0;
339	<	double mass;
338	>	RealType totalMass = 0;
339	>	RealType mass;
340
341	<	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
342	<	mass = sd->getMass();
343	<	totalMass += mass;
344	<	velCom += sd->getVel() * mass;
341	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
342	>	nextIntegrableObject(i)){
343	>	mass = sd->getMass();
344	>	totalMass += mass;
345	>	velCom += sd->getVel() * mass;
346		}
347	<
347	>
348		velCom /= totalMass;
349	<
349	>
350		return velCom;
351	<	}
351	>	}
352
353	<	double Molecule::getPotential() {
353	>	RealType Molecule::getPotential() {
354
355		Bond* bond;
356		Bend* bend;
357		Torsion* torsion;
358	+	Inversion* inversion;
359		Molecule::BondIterator bondIter;;
360		Molecule::BendIterator  bendIter;
361		Molecule::TorsionIterator  torsionIter;
362	+	Molecule::InversionIterator  inversionIter;
363
364	<	double potential = 0.0;
364	>	RealType potential = 0.0;
365
366		for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
367	<	potential += bond->getPotential();
367	>	potential += bond->getPotential();
368		}
369
370		for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
371	<	potential += bend->getPotential();
371	>	potential += bend->getPotential();
372		}
373
374	<	for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
375	<	potential += torsion->getPotential();
374	>	for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
375	>	nextTorsion(torsionIter)) {
376	>	potential += torsion->getPotential();
377		}
378	<
378	>
379	>	for (inversion = beginInversion(inversionIter); torsion != NULL;
380	>	inversion =  nextInversion(inversionIter)) {
381	>	potential += inversion->getPotential();
382	>	}
383	>
384		return potential;
385	+
386	+	}
387	+
388	+	void Molecule::addProperty(GenericData* genData) {
389	+	properties_.addProperty(genData);
390	+	}
391
392	<	}
392	>	void Molecule::removeProperty(const std::string& propName) {
393	>	properties_.removeProperty(propName);
394	>	}
395
396	<	std::ostream& operator <<(std::ostream& o, Molecule& mol) {
396	>	void Molecule::clearProperties() {
397	>	properties_.clearProperties();
398	>	}
399	>
400	>	std::vector<std::string> Molecule::getPropertyNames() {
401	>	return properties_.getPropertyNames();
402	>	}
403	>
404	>	std::vector<GenericData*> Molecule::getProperties() {
405	>	return properties_.getProperties();
406	>	}
407	>
408	>	GenericData* Molecule::getPropertyByName(const std::string& propName) {
409	>	return properties_.getPropertyByName(propName);
410	>	}
411	>
412	>	std::ostream& operator <<(std::ostream& o, Molecule& mol) {
413		o << std::endl;
414		o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
415		o << mol.getNAtoms() << " atoms" << std::endl;
416		o << mol.getNBonds() << " bonds" << std::endl;
417		o << mol.getNBends() << " bends" << std::endl;
418		o << mol.getNTorsions() << " torsions" << std::endl;
419	+	o << mol.getNInversions() << " inversions" << std::endl;
420		o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
421	<	o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
422	<	o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
423	<	o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
421	>	o << mol.getNIntegrableObjects() << " integrable objects" << std::endl;
422	>	o << mol.getNCutoffGroups() << " cutoff groups" << std::endl;
423	>	o << mol.getNConstraintPairs() << " constraint pairs" << std::endl;
424	>	o << mol.getNFluctuatingCharges() << " fluctuating charges" << std::endl;
425		return o;
426	<	}
427	<
428	<	}//end namespace oopse
426	>	}
427	>
428	>	}//end namespace OpenMD

Comparing trunk/src/primitives/Molecule.cpp (property svn:keywords):
Revision 372 by tim, Mon Feb 21 15:28:56 2005 UTC vs.
Revision 2054 by gezelter, Tue Jan 13 20:14:57 2015 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines