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 2071 by gezelter, Sat Mar 7 21:41:51 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) :
61	>	globalIndex_(globalIndex),
62	>	stampId_(stampId),
63	>	region_(region),
64	>	moleculeName_(molName),
65	>	constrainTotalCharge_(false) {
66	>	}
67	>
68	>	Molecule::~Molecule() {
69	>
70	>	MemoryUtils::deletePointers(atoms_);
71	>	MemoryUtils::deletePointers(bonds_);
72	>	MemoryUtils::deletePointers(bends_);
73	>	MemoryUtils::deletePointers(torsions_);
74	>	MemoryUtils::deletePointers(inversions_);
75	>	MemoryUtils::deletePointers(rigidBodies_);
76	>	MemoryUtils::deletePointers(cutoffGroups_);
77	>	MemoryUtils::deletePointers(constraintPairs_);
78	>	MemoryUtils::deletePointers(constraintElems_);
79
80	<	}
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
80	>	// integrableObjects_ don't own the objects
81		integrableObjects_.clear();
82	+	fluctuatingCharges_.clear();
83
84	<	}
85	<
86	<	void Molecule::addAtom(Atom* atom) {
84	>	}
85	>
86	>	void Molecule::addAtom(Atom* atom) {
87		if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
88	<	atoms_.push_back(atom);
88	>	atoms_.push_back(atom);
89		}
90	<	}
91	<
92	<	void Molecule::addBond(Bond* bond) {
90	>	}
91	>
92	>	void Molecule::addBond(Bond* bond) {
93		if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
94	<	bonds_.push_back(bond);
94	>	bonds_.push_back(bond);
95		}
96	<	}
97	<
98	<	void Molecule::addBend(Bend* bend) {
96	>	}
97	>
98	>	void Molecule::addBend(Bend* bend) {
99		if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
100	<	bends_.push_back(bend);
100	>	bends_.push_back(bend);
101		}
102	<	}
103	<
104	<	void Molecule::addTorsion(Torsion* torsion) {
105	<	if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
106	<	torsions_.push_back(torsion);
102	>	}
103	>
104	>	void Molecule::addTorsion(Torsion* torsion) {
105	>	if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
106	>	torsions_.end()) {
107	>	torsions_.push_back(torsion);
108		}
109	<	}
109	>	}
110
111	<	void Molecule::addRigidBody(RigidBody *rb) {
112	<	if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
113	<	rigidBodies_.push_back(rb);
111	>	void Molecule::addInversion(Inversion* inversion) {
112	>	if (std::find(inversions_.begin(), inversions_.end(), inversion) ==
113	>	inversions_.end()) {
114	>	inversions_.push_back(inversion);
115		}
116	<	}
117	<
118	<	void Molecule::addCutoffGroup(CutoffGroup* cp) {
119	<	if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
120	<	cutoffGroups_.push_back(cp);
116	>	}
117	>
118	>	void Molecule::addRigidBody(RigidBody *rb) {
119	>	if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
120	>	rigidBodies_.end()) {
121	>	rigidBodies_.push_back(rb);
122		}
123	<
124	<	}
125	<
126	<	void Molecule::addConstraintPair(ConstraintPair* cp) {
127	<	if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) {
128	<	constraintPairs_.push_back(cp);
123	>	}
124	>
125	>	void Molecule::addCutoffGroup(CutoffGroup* cp) {
126	>	if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
127	>	cutoffGroups_.end()) {
128	>	cutoffGroups_.push_back(cp);
129	>	}
130	>	}
131	>
132	>	void Molecule::addConstraintPair(ConstraintPair* cp) {
133	>	if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
134	>	constraintPairs_.end()) {
135	>	constraintPairs_.push_back(cp);
136	>	}
137	>	}
138	>
139	>	void Molecule::addConstraintElem(ConstraintElem* cp) {
140	>	if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
141	>	constraintElems_.end()) {
142	>	constraintElems_.push_back(cp);
143		}
144	<
145	<	}
146	<
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() {
144	>	}
145	>
146	>	void Molecule::complete() {
147
148		std::set<Atom*> rigidAtoms;
149	+	Atom* atom;
150	+	Atom* atom1;
151	+	Atom* atom2;
152	+	AtomIterator ai, aj;
153		RigidBody* rb;
154	<	std::vector<RigidBody*>::iterator rbIter;
154	>	RigidBodyIterator rbIter;
155	>	Bond* bond;
156	>	BondIterator bi;
157
158	<
158	>	// Get list of all the atoms that are part of rigid bodies
159	>
160		for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
161	<	rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
161	>	rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
162		}
163	+
164	+	// add any atom that wasn't part of a rigid body to the list of integrableObjects
165
139	–	Atom* atom;
140	–	AtomIterator ai;
166		for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
167	<
168	<	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	<	}
167	>
168	>	if (rigidAtoms.find(atom) == rigidAtoms.end()) {
169
170	<	//find all free atoms (which do not belong to rigid bodies)
171	<	//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_));
170	>	// If an atom does not belong to a rigid body, it is an
171	>	// integrable object
172
173	<	//if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
174	<	//    //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
175	<	//    sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
176	<	//
177	<	//    painCave.isFatal = 1;
178	<	//    simError();
162	<	//}
173	>	integrableObjects_.push_back(atom);
174	>	}
175	>	}
176	>
177	>	// then add the rigid bodies themselves to the integrableObjects
178	>
179		for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
180		integrableObjects_.push_back(rb);
181		}
166	–	//integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
167	–	}
182
183	<	double Molecule::getMass() {
184	<	StuntDouble* sd;
171	<	std::vector<StuntDouble*>::iterator i;
172	<	double mass = 0.0;
183	>	// find the atoms that are fluctuating charges and add them to the
184	>	// fluctuatingCharges_ vector
185
186	<	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
187	<	mass += sd->getMass();
186	>	for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
187	>	if ( atom->isFluctuatingCharge() )
188	>	fluctuatingCharges_.push_back( atom );
189		}
190
178	–	return mass;
191
192	<	}
192	>	// find the electronegative atoms and add them to the
193	>	// hBondAcceptors_ vector:
194	>
195	>	for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
196	>	AtomType* at = atom->getAtomType();
197	>	// get the chain of base types for this atom type:
198	>	std::vector<AtomType*> ayb = at->allYourBase();
199	>	// use the last type in the chain of base types for the name:
200	>	std::string bn = UpperCase(ayb[ayb.size()-1]->getName());
201
202	<	Vector3d Molecule::getCom() {
202	>	if (bn.compare("O")==0 || bn.compare("N")==0
203	>	|| bn.compare("F")==0)
204	>	hBondAcceptors_.push_back( atom );
205	>
206	>	}
207	>
208	>	// find electronegative atoms that are either bonded to
209	>	// hydrogens or are present in the same rigid bodies:
210	>
211	>	for (bond = beginBond(bi); bond != NULL; bond = nextBond(bi)) {
212	>	Atom* atom1 = bond->getAtomA();
213	>	Atom* atom2 = bond->getAtomB();
214	>	AtomType* at1 = atom1->getAtomType();
215	>	AtomType* at2 = atom1->getAtomType();
216	>	// get the chain of base types for this atom type:
217	>	std::vector<AtomType*> ayb1 = at1->allYourBase();
218	>	std::vector<AtomType*> ayb2 = at2->allYourBase();
219	>	// use the last type in the chain of base types for the name:
220	>	std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName());
221	>	std::string bn2 = UpperCase(ayb2[ayb2.size()-1]->getName());
222	>
223	>	if (bn1.compare("H")==0) {
224	>	if (bn2.compare("O")==0 || bn2.compare("N")==0
225	>	|| bn2.compare("F")==0) {
226	>	HBondDonor* donor = new HBondDonor();
227	>	donor->donorAtom = atom2;
228	>	donor->donatedHydrogen = atom1;
229	>	hBondDonors_.push_back( donor );
230	>	}
231	>	}
232	>	if (bn2.compare("H")==0) {
233	>	if (bn1.compare("O")==0 || bn1.compare("N")==0
234	>	|| bn1.compare("F")==0) {
235	>	HBondDonor* donor = new HBondDonor();
236	>	donor->donorAtom = atom1;
237	>	donor->donatedHydrogen = atom2;
238	>	hBondDonors_.push_back( donor );
239	>	}
240	>	}
241	>	}
242	>
243	>	for (rb = beginRigidBody(rbIter); rb != NULL;
244	>	rb = nextRigidBody(rbIter)) {
245	>	for(atom1 = rb->beginAtom(ai); atom1 != NULL;
246	>	atom1 = rb->nextAtom(ai)) {
247	>	AtomType* at1 = atom1->getAtomType();
248	>	// get the chain of base types for this atom type:
249	>	std::vector<AtomType*> ayb1 = at1->allYourBase();
250	>	// use the last type in the chain of base types for the name:
251	>	std::string bn1 = UpperCase(ayb1[ayb1.size()-1]->getName());
252	>
253	>	if (bn1.compare("O")==0 || bn1.compare("N")==0
254	>	|| bn1.compare("F")==0) {
255	>	for(atom2 = rb->beginAtom(aj); atom2 != NULL;
256	>	atom2 = rb->nextAtom(aj)) {
257	>	AtomType* at2 = atom2->getAtomType();
258	>	// get the chain of base types for this atom type:
259	>	std::vector<AtomType*> ayb2 = at2->allYourBase();
260	>	// use the last type in the chain of base types for the name:
261	>	std::string bn2 = UpperCase(ayb2[ayb2.size()-1]->getName());
262	>	if (bn2.compare("H")==0) {
263	>	HBondDonor* donor = new HBondDonor();
264	>	donor->donorAtom = atom1;
265	>	donor->donatedHydrogen = atom2;
266	>	hBondDonors_.push_back( donor );
267	>	}
268	>	}
269	>	}
270	>	}
271	>	}
272	>	}
273	>
274	>	RealType Molecule::getMass() {
275		StuntDouble* sd;
276		std::vector<StuntDouble*>::iterator i;
277	<	Vector3d com;
186	<	double totalMass = 0;
187	<	double mass;
277	>	RealType mass = 0.0;
278
279	<	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
280	<	mass = sd->getMass();
281	<	totalMass += mass;
192	<	com += sd->getPos() * mass;
279	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
280	>	nextIntegrableObject(i)){
281	>	mass += sd->getMass();
282		}
283	+
284	+	return mass;
285	+	}
286
287	+	Vector3d Molecule::getCom() {
288	+	StuntDouble* sd;
289	+	std::vector<StuntDouble*>::iterator i;
290	+	Vector3d com;
291	+	RealType totalMass = 0;
292	+	RealType mass;
293	+
294	+	for (sd = beginIntegrableObject(i); sd != NULL; sd =
295	+	nextIntegrableObject(i)){
296	+	mass = sd->getMass();
297	+	totalMass += mass;
298	+	com += sd->getPos() * mass;
299	+	}
300	+
301		com /= totalMass;
302
303		return com;
304	<	}
305	<
306	<	void Molecule::moveCom(const Vector3d& delta) {
304	>	}
305	>
306	>	Vector3d Molecule::getCom(int snapshotNo) {
307		StuntDouble* sd;
308		std::vector<StuntDouble*>::iterator i;
309	+	Vector3d com;
310	+	RealType totalMass = 0;
311	+	RealType mass;
312
313	<	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
314	<	sd->setPos(sd->getPos() + delta);
313	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
314	>	nextIntegrableObject(i)){
315	>	mass = sd->getMass();
316	>	totalMass += mass;
317	>	com += sd->getPos(snapshotNo) * mass;
318		}
319	+
320	+	com /= totalMass;
321
322	<	}
322	>	return com;
323	>	}
324
325	<	Vector3d Molecule::getComVel() {
325	>	void Molecule::moveCom(const Vector3d& delta) {
326		StuntDouble* sd;
327		std::vector<StuntDouble*>::iterator i;
328	+
329	+	for (sd = beginIntegrableObject(i); sd != NULL; sd =
330	+	nextIntegrableObject(i)){
331	+	sd->setPos(sd->getPos() + delta);
332	+	}
333	+	}
334	+
335	+	Vector3d Molecule::getComVel() {
336	+	StuntDouble* sd;
337	+	std::vector<StuntDouble*>::iterator i;
338		Vector3d velCom;
339	<	double totalMass = 0;
340	<	double mass;
339	>	RealType totalMass = 0;
340	>	RealType mass;
341
342	<	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
343	<	mass = sd->getMass();
344	<	totalMass += mass;
345	<	velCom += sd->getVel() * mass;
342	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
343	>	nextIntegrableObject(i)){
344	>	mass = sd->getMass();
345	>	totalMass += mass;
346	>	velCom += sd->getVel() * mass;
347		}
348	<
348	>
349		velCom /= totalMass;
350	<
350	>
351		return velCom;
352	<	}
352	>	}
353
354	<	double Molecule::getPotential() {
354	>	RealType Molecule::getPotential() {
355
356		Bond* bond;
357		Bend* bend;
358		Torsion* torsion;
359	+	Inversion* inversion;
360		Molecule::BondIterator bondIter;;
361		Molecule::BendIterator  bendIter;
362		Molecule::TorsionIterator  torsionIter;
363	+	Molecule::InversionIterator  inversionIter;
364
365	<	double potential = 0.0;
365	>	RealType potential = 0.0;
366
367		for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
368	<	potential += bond->getPotential();
368	>	potential += bond->getPotential();
369		}
370
371		for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
372	<	potential += bend->getPotential();
372	>	potential += bend->getPotential();
373		}
374
375	<	for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
376	<	potential += torsion->getPotential();
375	>	for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
376	>	nextTorsion(torsionIter)) {
377	>	potential += torsion->getPotential();
378		}
379	<
379	>
380	>	for (inversion = beginInversion(inversionIter); torsion != NULL;
381	>	inversion =  nextInversion(inversionIter)) {
382	>	potential += inversion->getPotential();
383	>	}
384	>
385		return potential;
386	+
387	+	}
388	+
389	+	void Molecule::addProperty(GenericData* genData) {
390	+	properties_.addProperty(genData);
391	+	}
392
393	<	}
393	>	void Molecule::removeProperty(const std::string& propName) {
394	>	properties_.removeProperty(propName);
395	>	}
396
397	<	std::ostream& operator <<(std::ostream& o, Molecule& mol) {
397	>	void Molecule::clearProperties() {
398	>	properties_.clearProperties();
399	>	}
400	>
401	>	std::vector<std::string> Molecule::getPropertyNames() {
402	>	return properties_.getPropertyNames();
403	>	}
404	>
405	>	std::vector<GenericData*> Molecule::getProperties() {
406	>	return properties_.getProperties();
407	>	}
408	>
409	>	GenericData* Molecule::getPropertyByName(const std::string& propName) {
410	>	return properties_.getPropertyByName(propName);
411	>	}
412	>
413	>	std::ostream& operator <<(std::ostream& o, Molecule& mol) {
414		o << std::endl;
415		o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
416		o << mol.getNAtoms() << " atoms" << std::endl;
417		o << mol.getNBonds() << " bonds" << std::endl;
418		o << mol.getNBends() << " bends" << std::endl;
419		o << mol.getNTorsions() << " torsions" << std::endl;
420	+	o << mol.getNInversions() << " inversions" << std::endl;
421		o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
422	<	o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
423	<	o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
424	<	o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
422	>	o << mol.getNIntegrableObjects() << " integrable objects" << std::endl;
423	>	o << mol.getNCutoffGroups() << " cutoff groups" << std::endl;
424	>	o << mol.getNConstraintPairs() << " constraint pairs" << std::endl;
425	>	o << mol.getNFluctuatingCharges() << " fluctuating charges" << std::endl;
426		return o;
427	<	}
428	<
429	<	}//end namespace oopse
427	>	}
428	>
429	>	}//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 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

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

Diff Legend

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