ViewVC Help

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

root/OpenMD/branches/development/src/primitives/Molecule.cpp

Comparing:
trunk/src/primitives/Molecule.cpp (file contents), Revision 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
branches/development/src/primitives/Molecule.cpp (file contents), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

#	Line 1 | Line 1
1	<	#include <stdlib.h>
1	>	/*
2	>	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	>	*
4	>	* The University of Notre Dame grants you ("Licensee") a
5	>	* non-exclusive, royalty free, license to use, modify and
6	>	* redistribute this software in source and binary code form, provided
7	>	* that the following conditions are met:
8	>	*
9	>	* 1. Redistributions of source code must retain the above copyright
10	>	*    notice, this list of conditions and the following disclaimer.
11	>	*
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.
16	>	*
17	>	* This software is provided "AS IS," without a warranty of any
18	>	* kind. All express or implied conditions, representations and
19	>	* warranties, including any implied warranty of merchantability,
20	>	* fitness for a particular purpose or non-infringement, are hereby
21	>	* excluded.  The University of Notre Dame and its licensors shall not
22	>	* be liable for any damages suffered by licensee as a result of
23	>	* using, modifying or distributing the software or its
24	>	* derivatives. In no event will the University of Notre Dame or its
25	>	* licensors be liable for any lost revenue, profit or data, or for
26	>	* direct, indirect, special, consequential, incidental or punitive
27	>	* damages, however caused and regardless of the theory of liability,
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, 24107 (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	>	/**
44	>	* @file Molecule.cpp
45	>	* @author    tlin
46	>	* @date  10/28/2004
47	>	* @version 1.0
48	>	*/
49
50	+	#include <algorithm>
51	+	#include <set>
52
53		#include "primitives/Molecule.hpp"
54	+	#include "utils/MemoryUtils.hpp"
55		#include "utils/simError.h"
56
57	<
58	<
59	<	Molecule::Molecule( void ){
60	<
11	<	myAtoms = NULL;
12	<	myBonds = NULL;
13	<	myBends = NULL;
14	<	myTorsions = NULL;
15	<	}
16	<
17	<	Molecule::~Molecule( void ){
18	<	int i;
19	<	CutoffGroup* cg;
20	<	vector<CutoffGroup*>::iterator iter;
57	>	namespace OpenMD {
58	>	Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
59	>	: stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
60	>	}
61
62	<	if( myAtoms != NULL ){
63	<	for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
64	<	delete[] myAtoms;
62	>	Molecule::~Molecule() {
63	>
64	>	MemoryUtils::deletePointers(atoms_);
65	>	MemoryUtils::deletePointers(bonds_);
66	>	MemoryUtils::deletePointers(bends_);
67	>	MemoryUtils::deletePointers(torsions_);
68	>	MemoryUtils::deletePointers(inversions_);
69	>	MemoryUtils::deletePointers(rigidBodies_);
70	>	MemoryUtils::deletePointers(cutoffGroups_);
71	>	MemoryUtils::deletePointers(constraintPairs_);
72	>	MemoryUtils::deletePointers(constraintElems_);
73	>	// integrableObjects_ don't own the objects
74	>	integrableObjects_.clear();
75	>
76		}
77	<
78	<	if( myBonds != NULL ){
79	<	for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i];
80	<	delete[] myBonds;
77	>
78	>	void Molecule::addAtom(Atom* atom) {
79	>	if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
80	>	atoms_.push_back(atom);
81	>	}
82		}
83	<
84	<	if( myBends != NULL ){
85	<	for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i];
86	<	delete[] myBends;
83	>
84	>	void Molecule::addBond(Bond* bond) {
85	>	if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
86	>	bonds_.push_back(bond);
87	>	}
88		}
89	<
90	<	if( myTorsions != NULL ){
91	<	for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i];
92	<	delete[] myTorsions;
89	>
90	>	void Molecule::addBend(Bend* bend) {
91	>	if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
92	>	bends_.push_back(bend);
93	>	}
94		}
41	–
42	–	for(cg = beginCutoffGroup(iter);  cg != NULL; cg = nextCutoffGroup(iter))
43	–	delete cg;
44	–	myCutoffGroups.clear();
95
96	<	}
97	<
98	<
99	<	void Molecule::initialize( molInit &theInit ){
100	<
51	<	CutoffGroup* curCutoffGroup;
52	<	vector<CutoffGroup*>::iterator iterCutoff;
53	<	Atom* cutoffAtom;
54	<	vector<Atom*>::iterator iterAtom;
55	<	int atomIndex;
56	<
57	<	nAtoms = theInit.nAtoms;
58	<	nMembers = nAtoms;
59	<	nBonds = theInit.nBonds;
60	<	nBends = theInit.nBends;
61	<	nTorsions = theInit.nTorsions;
62	<	nRigidBodies = theInit.nRigidBodies;
63	<	nOriented = theInit.nOriented;
64	<
65	<	myAtoms = theInit.myAtoms;
66	<	myBonds = theInit.myBonds;
67	<	myBends = theInit.myBends;
68	<	myTorsions = theInit.myTorsions;
69	<	myRigidBodies = theInit.myRigidBodies;
70	<
71	<	myIntegrableObjects = theInit.myIntegrableObjects;
72	<
73	<	for (int i = 0; i < myRigidBodies.size(); i++)
74	<	myRigidBodies[i]->calcRefCoords();
75	<
76	<	myCutoffGroups = theInit.myCutoffGroups;
77	<	nCutoffGroups = myCutoffGroups.size();
78	<
79	<	}
80	<
81	<	void Molecule::calcForces( void ){
82	<
83	<	int i;
84	<	double com[3];
85	<
86	<	for(i=0; i<myRigidBodies.size(); i++) {
87	<	myRigidBodies[i]->updateAtoms();
96	>	void Molecule::addTorsion(Torsion* torsion) {
97	>	if (std::find(torsions_.begin(), torsions_.end(), torsion) ==
98	>	torsions_.end()) {
99	>	torsions_.push_back(torsion);
100	>	}
101		}
102
103	<	for(i=0; i<nBonds; i++){
104	<	myBonds[i]->calc_forces();
103	>	void Molecule::addInversion(Inversion* inversion) {
104	>	if (std::find(inversions_.begin(), inversions_.end(), inversion) ==
105	>	inversions_.end()) {
106	>	inversions_.push_back(inversion);
107	>	}
108		}
109	<
110	<	for(i=0; i<nBends; i++){
111	<	myBends[i]->calc_forces();
109	>
110	>	void Molecule::addRigidBody(RigidBody *rb) {
111	>	if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) ==
112	>	rigidBodies_.end()) {
113	>	rigidBodies_.push_back(rb);
114	>	}
115		}
116	<
117	<	for(i=0; i<nTorsions; i++){
118	<	myTorsions[i]->calc_forces();
116	>
117	>	void Molecule::addCutoffGroup(CutoffGroup* cp) {
118	>	if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) ==
119	>	cutoffGroups_.end()) {
120	>	cutoffGroups_.push_back(cp);
121	>	}
122		}
101	–
102	–	// Rigid Body forces and torques are done after the fortran force loop
103	–
104	–	}
105	–
106	–
107	–	double Molecule::getPotential( void ){
123
124	<	int i;
125	<	double myPot = 0.0;
126	<
127	<	for(i=0; i<myRigidBodies.size(); i++) {
128	<	myRigidBodies[i]->updateAtoms();
124	>	void Molecule::addConstraintPair(ConstraintPair* cp) {
125	>	if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) ==
126	>	constraintPairs_.end()) {
127	>	constraintPairs_.push_back(cp);
128	>	}
129		}
130	<
131	<	for(i=0; i<nBonds; i++){
132	<	myPot += myBonds[i]->get_potential();
130	>
131	>	void Molecule::addConstraintElem(ConstraintElem* cp) {
132	>	if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) ==
133	>	constraintElems_.end()) {
134	>	constraintElems_.push_back(cp);
135	>	}
136		}
119	–
120	–	for(i=0; i<nBends; i++){
121	–	myPot += myBends[i]->get_potential();
122	–	}
123	–
124	–	for(i=0; i<nTorsions; i++){
125	–	myPot += myTorsions[i]->get_potential();
126	–	}
127	–
128	–	return myPot;
129	–	}
130	–
131	–	void Molecule::printMe( void ){
137
138	<	int i;
138	>	void Molecule::complete() {
139	>
140	>	std::set<Atom*> rigidAtoms;
141	>	RigidBody* rb;
142	>	std::vector<RigidBody*>::iterator rbIter;
143
144	<	for(i=0; i<nBonds; i++){
145	<	myBonds[i]->printMe();
146	<	}
144	>
145	>	for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
146	>	rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
147	>	}
148	>
149	>	Atom* atom;
150	>	AtomIterator ai;
151	>	for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
152	>
153	>	if (rigidAtoms.find(*ai) == rigidAtoms.end()) {
154
155	<	for(i=0; i<nBends; i++){
156	<	myBends[i]->printMe();
141	<	}
155	>	// If an atom does not belong to a rigid body, it is an
156	>	// integrable object
157
158	<	for(i=0; i<nTorsions; i++){
159	<	myTorsions[i]->printMe();
158	>	integrableObjects_.push_back(*ai);
159	>	}
160	>	}
161	>
162	>	//find all free atoms (which do not belong to rigid bodies)
163	>	// performs the "difference" operation from set theory, the output
164	>	// range contains a copy of every element that is contained in
165	>	// [allAtoms.begin(), allAtoms.end()) and not contained in
166	>	// [rigidAtoms.begin(), rigidAtoms.end()).
167	>	//std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(),
168	>	//                        std::back_inserter(integrableObjects_));
169	>
170	>	//if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
171	>	//    //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
172	>	//    sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
173	>	//
174	>	//    painCave.isFatal = 1;
175	>	//    simError();
176	>	//}
177	>	for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
178	>	integrableObjects_.push_back(rb);
179	>	}
180	>	//integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
181		}
182
183	<	}
184	<
185	<	void Molecule::moveCOM(double delta[3]){
186	<	double aPos[3];
187	<	int i, j;
188	<
189	<	for(i=0; i<myIntegrableObjects.size(); i++) {
190	<	if(myIntegrableObjects[i] != NULL ) {
155	<
156	<	myIntegrableObjects[i]->getPos( aPos );
157	<
158	<	for (j=0; j< 3; j++)
159	<	aPos[j] += delta[j];
160	<
161	<	myIntegrableObjects[i]->setPos( aPos );
183	>	RealType Molecule::getMass() {
184	>	StuntDouble* sd;
185	>	std::vector<StuntDouble*>::iterator i;
186	>	RealType mass = 0.0;
187	>
188	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
189	>	nextIntegrableObject(i)){
190	>	mass += sd->getMass();
191		}
192	+
193	+	return mass;
194		}
195
196	<	for(i=0; i<myRigidBodies.size(); i++) {
197	<
198	<	myRigidBodies[i]->getPos( aPos );
199	<
200	<	for (j=0; j< 3; j++)
201	<	aPos[j] += delta[j];
202	<
203	<	myRigidBodies[i]->setPos( aPos );
196	>	Vector3d Molecule::getCom() {
197	>	StuntDouble* sd;
198	>	std::vector<StuntDouble*>::iterator i;
199	>	Vector3d com;
200	>	RealType totalMass = 0;
201	>	RealType mass;
202	>
203	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
204	>	nextIntegrableObject(i)){
205	>	mass = sd->getMass();
206	>	totalMass += mass;
207	>	com += sd->getPos() * mass;
208		}
209	<	}
209	>
210	>	com /= totalMass;
211
212	<	void Molecule::atoms2rigidBodies( void ) {
177	<	int i;
178	<	for (i = 0; i < myRigidBodies.size(); i++) {
179	<	myRigidBodies[i]->calcForcesAndTorques();
212	>	return com;
213		}
181	–	}
214
215	<	void Molecule::getCOM( double COM[3] ) {
215	>	void Molecule::moveCom(const Vector3d& delta) {
216	>	StuntDouble* sd;
217	>	std::vector<StuntDouble*>::iterator i;
218	>
219	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
220	>	nextIntegrableObject(i)){
221	>	sd->setPos(sd->getPos() + delta);
222	>	}
223	>	}
224
225	<	double mass, mtot;
226	<	double aPos[3];
227	<	int i, j;
225	>	Vector3d Molecule::getComVel() {
226	>	StuntDouble* sd;
227	>	std::vector<StuntDouble*>::iterator i;
228	>	Vector3d velCom;
229	>	RealType totalMass = 0;
230	>	RealType mass;
231	>
232	>	for (sd = beginIntegrableObject(i); sd != NULL; sd =
233	>	nextIntegrableObject(i)){
234	>	mass = sd->getMass();
235	>	totalMass += mass;
236	>	velCom += sd->getVel() * mass;
237	>	}
238	>
239	>	velCom /= totalMass;
240	>
241	>	return velCom;
242	>	}
243
244	<	for (j=0; j<3; j++)
190	<	COM[j] = 0.0;
244	>	RealType Molecule::getPotential() {
245
246	<	mtot   = 0.0;
246	>	Bond* bond;
247	>	Bend* bend;
248	>	Torsion* torsion;
249	>	Inversion* inversion;
250	>	Molecule::BondIterator bondIter;;
251	>	Molecule::BendIterator  bendIter;
252	>	Molecule::TorsionIterator  torsionIter;
253	>	Molecule::InversionIterator  inversionIter;
254
255	<	for (i=0; i < myIntegrableObjects.size(); i++) {
195	<	if (myIntegrableObjects[i] != NULL) {
255	>	RealType potential = 0.0;
256
257	<	mass = myIntegrableObjects[i]->getMass();
258	<	mtot   += mass;
259	<
200	<	myIntegrableObjects[i]->getPos( aPos );
257	>	for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
258	>	potential += bond->getPotential();
259	>	}
260
261	<	for( j = 0; j < 3; j++)
262	<	COM[j] += aPos[j] * mass;
261	>	for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
262	>	potential += bend->getPotential();
263	>	}
264
265	+	for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion =
266	+	nextTorsion(torsionIter)) {
267	+	potential += torsion->getPotential();
268		}
269	+
270	+	for (inversion = beginInversion(inversionIter); torsion != NULL;
271	+	inversion =  nextInversion(inversionIter)) {
272	+	potential += inversion->getPotential();
273	+	}
274	+
275	+	return potential;
276	+
277		}
278	+
279	+	void Molecule::addProperty(GenericData* genData) {
280	+	properties_.addProperty(genData);
281	+	}
282
283	<	for (j = 0; j < 3; j++)
284	<	COM[j] /= mtot;
285	<	}
283	>	void Molecule::removeProperty(const std::string& propName) {
284	>	properties_.removeProperty(propName);
285	>	}
286
287	<	double Molecule::getCOMvel( double COMvel[3] ) {
287	>	void Molecule::clearProperties() {
288	>	properties_.clearProperties();
289	>	}
290
291	<	double mass, mtot;
292	<	double aVel[3];
293	<	int i, j;
291	>	std::vector<std::string> Molecule::getPropertyNames() {
292	>	return properties_.getPropertyNames();
293	>	}
294	>
295	>	std::vector<GenericData*> Molecule::getProperties() {
296	>	return properties_.getProperties();
297	>	}
298
299	<
300	<	for (j=0; j<3; j++)
220	<	COMvel[j] = 0.0;
221	<
222	<	mtot   = 0.0;
223	<
224	<	for (i=0; i < myIntegrableObjects.size(); i++) {
225	<	if (myIntegrableObjects[i] != NULL) {
226	<
227	<	mass = myIntegrableObjects[i]->getMass();
228	<	mtot   += mass;
229	<
230	<	myIntegrableObjects[i]->getVel(aVel);
231	<
232	<	for (j=0; j<3; j++)
233	<	COMvel[j] += aVel[j]*mass;
234	<
235	<	}
299	>	GenericData* Molecule::getPropertyByName(const std::string& propName) {
300	>	return properties_.getPropertyByName(propName);
301		}
302
238	–	for (j=0; j<3; j++)
239	–	COMvel[j] /= mtot;
240	–
241	–	return mtot;
303
243	–	}
304
245	–	double Molecule::getTotalMass()
246	–	{
305
306	<	double totalMass;
307	<
308	<	totalMass = 0;
309	<	for(int i =0; i < myIntegrableObjects.size(); i++){
310	<	totalMass += myIntegrableObjects[i]->getMass();
306	>	std::ostream& operator <<(std::ostream& o, Molecule& mol) {
307	>	o << std::endl;
308	>	o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
309	>	o << mol.getNAtoms() << " atoms" << std::endl;
310	>	o << mol.getNBonds() << " bonds" << std::endl;
311	>	o << mol.getNBends() << " bends" << std::endl;
312	>	o << mol.getNTorsions() << " torsions" << std::endl;
313	>	o << mol.getNInversions() << " inversions" << std::endl;
314	>	o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
315	>	o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
316	>	o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
317	>	o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
318	>	return o;
319		}
320	<
321	<	return totalMass;
256	<	}
320	>
321	>	}//end namespace OpenMD

Comparing:
trunk/src/primitives/Molecule.cpp (property svn:keywords), Revision 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
branches/development/src/primitives/Molecule.cpp (property svn:keywords), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

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

Diff Legend

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