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 2053 by gezelter, Mon Jan 12 16:04:43 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/ElementsTable.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 = ayb[ayb.size()-1]->getName();
200	>
201	>	if (bn.compare("O")==0 || bn.compare("N")==0 || bn.compare("F")==0)
202	>	hBondAcceptors_.push_back( atom );
203
204	<	Vector3d Molecule::getCom() {
204	>	}
205	>
206	>	// find electronegative atoms that are either bonded to
207	>	// hydrogens or are present in the same rigid bodies:
208	>
209	>	for (bond = beginBond(bi); bond != NULL; bond = nextBond(bi)) {
210	>	Atom* atom1 = bond->getAtomA();
211	>	Atom* atom2 = bond->getAtomB();
212	>	AtomType* at1 = atom1->getAtomType();
213	>	AtomType* at2 = atom1->getAtomType();
214	>	// get the chain of base types for this atom type:
215	>	std::vector<AtomType*> ayb1 = at1->allYourBase();
216	>	std::vector<AtomType*> ayb2 = at2->allYourBase();
217	>	// use the last type in the chain of base types for the name:
218	>	std::string bn1 = ayb1[ayb1.size()-1]->getName();
219	>	std::string bn2 = ayb2[ayb2.size()-1]->getName();
220	>
221	>	if (bn1.compare("H")==0) {
222	>	if (bn2.compare("O")==0 || bn2.compare("N")==0
223	>	|| bn2.compare("F")==0) {
224	>	HBondDonor* donor = new HBondDonor();
225	>	donor->donorAtom = atom2;
226	>	donor->donatedHydrogen = atom1;
227	>	hBondDonors_.push_back( donor );
228	>	}
229	>	}
230	>	if (bn2.compare("H")==0) {
231	>	if (bn1.compare("O")==0 || bn1.compare("N")==0
232	>	|| bn1.compare("F")==0) {
233	>	HBondDonor* donor = new HBondDonor();
234	>	donor->donorAtom = atom1;
235	>	donor->donatedHydrogen = atom2;
236	>	hBondDonors_.push_back( donor );
237	>	}
238	>	}
239	>	}
240	>
241	>	for (rb = beginRigidBody(rbIter); rb != NULL;
242	>	rb = nextRigidBody(rbIter)) {
243	>	for(atom1 = rb->beginAtom(ai); atom1 != NULL;
244	>	atom1 = rb->nextAtom(ai)) {
245	>	AtomType* at1 = atom1->getAtomType();
246	>	// get the chain of base types for this atom type:
247	>	std::vector<AtomType*> ayb1 = at1->allYourBase();
248	>	// use the last type in the chain of base types for the name:
249	>	std::string bn1 = ayb1[ayb1.size()-1]->getName();
250	>
251	>	if (bn1.compare("O")==0 || bn1.compare("N")==0
252	>	|| bn1.compare("F")==0) {
253	>	for(atom2 = rb->beginAtom(aj); atom2 != NULL;
254	>	atom2 = rb->nextAtom(aj)) {
255	>	AtomType* at2 = atom2->getAtomType();
256	>	// get the chain of base types for this atom type:
257	>	std::vector<AtomType*> ayb2 = at2->allYourBase();
258	>	// use the last type in the chain of base types for the name:
259	>	std::string bn2 = ayb2[ayb2.size()-1]->getName();
260	>	if (bn2.compare("H")==0) {
261	>	HBondDonor* donor = new HBondDonor();
262	>	donor->donorAtom = atom1;
263	>	donor->donatedHydrogen = atom2;
264	>	hBondDonors_.push_back( donor );
265	>	}
266	>	}
267	>	}
268	>	}
269	>	}
270	>	}
271	>
272	>	RealType Molecule::getMass() {
273		StuntDouble* sd;
274		std::vector<StuntDouble*>::iterator i;
275	<	Vector3d com;
186	<	double totalMass = 0;
187	<	double mass;
275	>	RealType mass = 0.0;
276
277	<	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
278	<	mass = sd->getMass();
279	<	totalMass += mass;
192	<	com += sd->getPos() * mass;
277	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
278	>	nextIntegrableObject(i)){
279	>	mass += sd->getMass();
280		}
281	+
282	+	return mass;
283	+	}
284
285	+	Vector3d Molecule::getCom() {
286	+	StuntDouble* sd;
287	+	std::vector<StuntDouble*>::iterator i;
288	+	Vector3d com;
289	+	RealType totalMass = 0;
290	+	RealType mass;
291	+
292	+	for (sd = beginIntegrableObject(i); sd != NULL; sd =
293	+	nextIntegrableObject(i)){
294	+	mass = sd->getMass();
295	+	totalMass += mass;
296	+	com += sd->getPos() * mass;
297	+	}
298	+
299		com /= totalMass;
300
301		return com;
302	<	}
303	<
304	<	void Molecule::moveCom(const Vector3d& delta) {
302	>	}
303	>
304	>	Vector3d Molecule::getCom(int snapshotNo) {
305		StuntDouble* sd;
306		std::vector<StuntDouble*>::iterator i;
307	+	Vector3d com;
308	+	RealType totalMass = 0;
309	+	RealType mass;
310
311	<	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
312	<	sd->setPos(sd->getPos() + delta);
311	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
312	>	nextIntegrableObject(i)){
313	>	mass = sd->getMass();
314	>	totalMass += mass;
315	>	com += sd->getPos(snapshotNo) * mass;
316		}
317	+
318	+	com /= totalMass;
319
320	<	}
320	>	return com;
321	>	}
322
323	<	Vector3d Molecule::getComVel() {
323	>	void Molecule::moveCom(const Vector3d& delta) {
324		StuntDouble* sd;
325		std::vector<StuntDouble*>::iterator i;
326	+
327	+	for (sd = beginIntegrableObject(i); sd != NULL; sd =
328	+	nextIntegrableObject(i)){
329	+	sd->setPos(sd->getPos() + delta);
330	+	}
331	+	}
332	+
333	+	Vector3d Molecule::getComVel() {
334	+	StuntDouble* sd;
335	+	std::vector<StuntDouble*>::iterator i;
336		Vector3d velCom;
337	<	double totalMass = 0;
338	<	double mass;
337	>	RealType totalMass = 0;
338	>	RealType mass;
339
340	<	for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
341	<	mass = sd->getMass();
342	<	totalMass += mass;
343	<	velCom += sd->getVel() * mass;
340	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
341	>	nextIntegrableObject(i)){
342	>	mass = sd->getMass();
343	>	totalMass += mass;
344	>	velCom += sd->getVel() * mass;
345		}
346	<
346	>
347		velCom /= totalMass;
348	<
348	>
349		return velCom;
350	<	}
350	>	}
351
352	<	double Molecule::getPotential() {
352	>	RealType Molecule::getPotential() {
353
354		Bond* bond;
355		Bend* bend;
356		Torsion* torsion;
357	+	Inversion* inversion;
358		Molecule::BondIterator bondIter;;
359		Molecule::BendIterator  bendIter;
360		Molecule::TorsionIterator  torsionIter;
361	+	Molecule::InversionIterator  inversionIter;
362
363	<	double potential = 0.0;
363	>	RealType potential = 0.0;
364
365		for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
366	<	potential += bond->getPotential();
366	>	potential += bond->getPotential();
367		}
368
369		for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
370	<	potential += bend->getPotential();
370	>	potential += bend->getPotential();
371		}
372
373	<	for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
374	<	potential += torsion->getPotential();
373	>	for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
374	>	nextTorsion(torsionIter)) {
375	>	potential += torsion->getPotential();
376		}
377	<
377	>
378	>	for (inversion = beginInversion(inversionIter); torsion != NULL;
379	>	inversion =  nextInversion(inversionIter)) {
380	>	potential += inversion->getPotential();
381	>	}
382	>
383		return potential;
384	+
385	+	}
386	+
387	+	void Molecule::addProperty(GenericData* genData) {
388	+	properties_.addProperty(genData);
389	+	}
390
391	<	}
391	>	void Molecule::removeProperty(const std::string& propName) {
392	>	properties_.removeProperty(propName);
393	>	}
394
395	<	std::ostream& operator <<(std::ostream& o, Molecule& mol) {
395	>	void Molecule::clearProperties() {
396	>	properties_.clearProperties();
397	>	}
398	>
399	>	std::vector<std::string> Molecule::getPropertyNames() {
400	>	return properties_.getPropertyNames();
401	>	}
402	>
403	>	std::vector<GenericData*> Molecule::getProperties() {
404	>	return properties_.getProperties();
405	>	}
406	>
407	>	GenericData* Molecule::getPropertyByName(const std::string& propName) {
408	>	return properties_.getPropertyByName(propName);
409	>	}
410	>
411	>	std::ostream& operator <<(std::ostream& o, Molecule& mol) {
412		o << std::endl;
413		o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
414		o << mol.getNAtoms() << " atoms" << std::endl;
415		o << mol.getNBonds() << " bonds" << std::endl;
416		o << mol.getNBends() << " bends" << std::endl;
417		o << mol.getNTorsions() << " torsions" << std::endl;
418	+	o << mol.getNInversions() << " inversions" << std::endl;
419		o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
420	<	o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
421	<	o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
422	<	o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
420	>	o << mol.getNIntegrableObjects() << " integrable objects" << std::endl;
421	>	o << mol.getNCutoffGroups() << " cutoff groups" << std::endl;
422	>	o << mol.getNConstraintPairs() << " constraint pairs" << std::endl;
423	>	o << mol.getNFluctuatingCharges() << " fluctuating charges" << std::endl;
424		return o;
425	<	}
426	<
427	<	}//end namespace oopse
425	>	}
426	>
427	>	}//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 2053 by gezelter, Mon Jan 12 16:04:43 2015 UTC

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

Diff Legend

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