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

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

Diff Legend

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