ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/primitives/Molecule.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/primitives/Molecule.cpp (file contents):
Revision 1492 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 1930 by gezelter, Wed Jan 12 22:41:40 2005 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. 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
19 > *    notice, this list of conditions and the following disclaimer.
20 > *
21 > * 3. Redistributions in binary form must reproduce the above copyright
22 > *    notice, this list of conditions and the following disclaimer in the
23 > *    documentation and/or other materials provided with the
24 > *    distribution.
25 > *
26 > * This software is provided "AS IS," without a warranty of any
27 > * kind. All express or implied conditions, representations and
28 > * warranties, including any implied warranty of merchantability,
29 > * fitness for a particular purpose or non-infringement, are hereby
30 > * excluded.  The University of Notre Dame and its licensors shall not
31 > * be liable for any damages suffered by licensee as a result of
32 > * using, modifying or distributing the software or its
33 > * derivatives. In no event will the University of Notre Dame or its
34 > * licensors be liable for any lost revenue, profit or data, or for
35 > * direct, indirect, special, consequential, incidental or punitive
36 > * damages, however caused and regardless of the theory of liability,
37 > * arising out of the use of or inability to use software, even if the
38 > * University of Notre Dame has been advised of the possibility of
39 > * such damages.
40 > */
41 >
42 > /**
43 > * @file Molecule.cpp
44 > * @author    tlin
45 > * @date  10/28/2004
46 > * @version 1.0
47 > */
48  
49 + #include <algorithm>
50 + #include <set>
51  
52   #include "primitives/Molecule.hpp"
53 + #include "utils/MemoryUtils.hpp"
54   #include "utils/simError.h"
55  
56 + namespace oopse {
57 + Molecule::Molecule(int stampId, int globalIndex, const std::string& molName)
58 +    : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) {
59  
60 + }
61  
62 < Molecule::Molecule( void ){
62 > Molecule::~Molecule() {
63  
64 <  myAtoms = NULL;
65 <  myBonds = NULL;
66 <  myBends = NULL;
67 <  myTorsions = NULL;
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
73 >    integrableObjects_.clear();
74 >    
75   }
76  
77 < Molecule::~Molecule( void ){
78 <  int i;
79 <  CutoffGroup* cg;
80 <  vector<CutoffGroup*>::iterator iter;
81 <  
22 <  if( myAtoms != NULL ){
23 <    for(i=0; i<nAtoms; i++) if(myAtoms[i] != NULL ) delete myAtoms[i];
24 <    delete[] myAtoms;
25 <  }
77 > void Molecule::addAtom(Atom* atom) {
78 >    if (std::find(atoms_.begin(), atoms_.end(), atom) == atoms_.end()) {
79 >        atoms_.push_back(atom);
80 >    }
81 > }
82  
83 <  if( myBonds != NULL ){
84 <    for(i=0; i<nBonds; i++) if(myBonds[i] != NULL ) delete myBonds[i];
85 <    delete[] myBonds;
86 <  }
83 > void Molecule::addBond(Bond* bond) {
84 >    if (std::find(bonds_.begin(), bonds_.end(), bond) == bonds_.end()) {
85 >        bonds_.push_back(bond);
86 >    }
87 > }
88  
89 <  if( myBends != NULL ){
90 <    for(i=0; i<nBends; i++) if(myBends[i] != NULL ) delete myBends[i];
91 <    delete[] myBends;
92 <  }
89 > void Molecule::addBend(Bend* bend) {
90 >    if (std::find(bends_.begin(), bends_.end(), bend) == bends_.end()) {
91 >        bends_.push_back(bend);
92 >    }
93 > }
94  
95 <  if( myTorsions != NULL ){
96 <    for(i=0; i<nTorsions; i++) if(myTorsions[i] != NULL ) delete myTorsions[i];
97 <    delete[] myTorsions;
98 <  }
95 > void Molecule::addTorsion(Torsion* torsion) {
96 >    if (std::find(torsions_.begin(), torsions_.end(), torsion) == torsions_.end()) {
97 >        torsions_.push_back(torsion);
98 >    }
99 > }
100  
101 <  for(cg = beginCutoffGroup(iter);  cg != NULL; cg = nextCutoffGroup(iter))
102 <    delete cg;
103 <  myCutoffGroups.clear();
104 <  
101 > void Molecule::addRigidBody(RigidBody *rb) {
102 >    if (std::find(rigidBodies_.begin(), rigidBodies_.end(), rb) == rigidBodies_.end()) {
103 >        rigidBodies_.push_back(rb);
104 >    }
105   }
106  
107 + void Molecule::addCutoffGroup(CutoffGroup* cp) {
108 +    if (std::find(cutoffGroups_.begin(), cutoffGroups_.end(), cp) == cutoffGroups_.end()) {
109 +        cutoffGroups_.push_back(cp);
110 +    }
111  
112 < void Molecule::initialize( molInit &theInit ){
112 > }
113  
114 <  CutoffGroup* curCutoffGroup;
115 <  vector<CutoffGroup*>::iterator iterCutoff;
116 <  Atom* cutoffAtom;
117 <  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;
114 > void Molecule::addConstraintPair(ConstraintPair* cp) {
115 >    if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) {
116 >        constraintPairs_.push_back(cp);
117 >    }
118  
119 <  myAtoms = theInit.myAtoms;
66 <  myBonds = theInit.myBonds;
67 <  myBends = theInit.myBends;
68 <  myTorsions = theInit.myTorsions;
69 <  myRigidBodies = theInit.myRigidBodies;
119 > }
120  
121 <  myIntegrableObjects = theInit.myIntegrableObjects;
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 <  for (int i = 0; i < myRigidBodies.size(); i++)
74 <      myRigidBodies[i]->calcRefCoords();
126 > }
127  
128 <  myCutoffGroups = theInit.myCutoffGroups;
129 <  nCutoffGroups = myCutoffGroups.size();
128 > void Molecule::complete() {
129 >    
130 >    std::set<Atom*> rigidAtoms;
131 >    RigidBody* rb;
132 >    std::vector<RigidBody*>::iterator rbIter;
133  
134 < }
134 >    
135 >    for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) {
136 >        rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter());
137 >    }
138  
139 < void Molecule::calcForces( void ){
140 <  
141 <  int i;
142 <  double com[3];
139 >    Atom* atom;
140 >    AtomIterator ai;
141 >    for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) {
142 >  
143 >        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 >    }
148  
149 <  for(i=0; i<myRigidBodies.size(); i++) {
150 <    myRigidBodies[i]->updateAtoms();
151 <  }
149 >    //find all free atoms (which do not belong to rigid bodies)  
150 >    //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_));
155  
156 <  for(i=0; i<nBonds; i++){
157 <    myBonds[i]->calc_forces();
158 <  }
156 >    //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) {
157 >    //    //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
158 >    //    sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule");
159 >    //
160 >    //    painCave.isFatal = 1;
161 >    //    simError();        
162 >    //}
163  
164 <  for(i=0; i<nBends; i++){
165 <    myBends[i]->calc_forces();
96 <  }
164 >    integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end());
165 > }
166  
167 <  for(i=0; i<nTorsions; i++){
168 <    myTorsions[i]->calc_forces();
169 <  }
167 > double Molecule::getMass() {
168 >    StuntDouble* sd;
169 >    std::vector<StuntDouble*>::iterator i;
170 >    double mass = 0.0;
171  
172 <  // Rigid Body forces and torques are done after the fortran force loop
172 >    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
173 >        mass += sd->getMass();
174 >    }
175  
176 +    return mass;
177 +
178   }
179  
180 + Vector3d Molecule::getCom() {
181 +    StuntDouble* sd;
182 +    std::vector<StuntDouble*>::iterator i;
183 +    Vector3d com;
184 +    double totalMass = 0;
185 +    double mass;
186 +    
187 +    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
188 +        mass = sd->getMass();
189 +        totalMass += mass;
190 +        com += sd->getPos() * mass;    
191 +    }
192  
193 < double Molecule::getPotential( void ){
108 <  
109 <  int i;
110 <  double myPot = 0.0;
193 >    com /= totalMass;
194  
195 <  for(i=0; i<myRigidBodies.size(); i++) {
113 <    myRigidBodies[i]->updateAtoms();
114 <  }
115 <  
116 <  for(i=0; i<nBonds; i++){
117 <    myPot += myBonds[i]->get_potential();
118 <  }
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;
195 >    return com;
196   }
197  
198 < void Molecule::printMe( void ){
199 <  
200 <  int i;
198 > void Molecule::moveCom(const Vector3d& delta) {
199 >    StuntDouble* sd;
200 >    std::vector<StuntDouble*>::iterator i;
201 >    
202 >    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
203 >        sd->setPos(sd->getPos() + delta);
204 >    }
205  
135  for(i=0; i<nBonds; i++){
136    myBonds[i]->printMe();
137  }
138
139  for(i=0; i<nBends; i++){
140    myBends[i]->printMe();
141  }
142
143  for(i=0; i<nTorsions; i++){
144    myTorsions[i]->printMe();
145  }
146
206   }
207  
208 < void Molecule::moveCOM(double delta[3]){
209 <  double aPos[3];
210 <  int i, j;
211 <
212 <  for(i=0; i<myIntegrableObjects.size(); i++) {
213 <    if(myIntegrableObjects[i] != NULL ) {
214 <      
215 <      myIntegrableObjects[i]->getPos( aPos );
216 <      
217 <      for (j=0; j< 3; j++)
218 <        aPos[j] += delta[j];
160 <
161 <      myIntegrableObjects[i]->setPos( aPos );
208 > Vector3d Molecule::getComVel() {
209 >    StuntDouble* sd;
210 >    std::vector<StuntDouble*>::iterator i;
211 >    Vector3d velCom;
212 >    double totalMass = 0;
213 >    double mass;
214 >    
215 >    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
216 >        mass = sd->getMass();
217 >        totalMass += mass;
218 >        velCom += sd->getVel() * mass;    
219      }
163  }
220  
221 <  for(i=0; i<myRigidBodies.size(); i++) {
221 >    velCom /= totalMass;
222  
223 <      myRigidBodies[i]->getPos( aPos );
168 <
169 <      for (j=0; j< 3; j++)
170 <        aPos[j] += delta[j];
171 <      
172 <      myRigidBodies[i]->setPos( aPos );
173 <    }
223 >    return velCom;
224   }
225  
226 < void Molecule::atoms2rigidBodies( void ) {
177 <  int i;
178 <  for (i = 0; i < myRigidBodies.size(); i++) {
179 <    myRigidBodies[i]->calcForcesAndTorques();  
180 <  }
181 < }
226 > double Molecule::getPotential() {
227  
228 < void Molecule::getCOM( double COM[3] ) {
228 >    Bond* bond;
229 >    Bend* bend;
230 >    Torsion* torsion;
231 >    Molecule::BondIterator bondIter;;
232 >    Molecule::BendIterator  bendIter;
233 >    Molecule::TorsionIterator  torsionIter;
234  
235 <  double mass, mtot;
186 <  double aPos[3];
187 <  int i, j;
235 >    double potential = 0.0;
236  
237 <  for (j=0; j<3; j++)
238 <    COM[j] = 0.0;
237 >    for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) {
238 >        potential += bond->getPotential();
239 >    }
240  
241 <  mtot   = 0.0;
242 <
194 <  for (i=0; i < myIntegrableObjects.size(); i++) {
195 <    if (myIntegrableObjects[i] != NULL) {
196 <
197 <      mass = myIntegrableObjects[i]->getMass();
198 <      mtot   += mass;
199 <      
200 <      myIntegrableObjects[i]->getPos( aPos );
201 <
202 <      for( j = 0; j < 3; j++)
203 <        COM[j] += aPos[j] * mass;
204 <
241 >    for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) {
242 >        potential += bend->getPotential();
243      }
206  }
244  
245 <  for (j = 0; j < 3; j++)
246 <    COM[j] /= mtot;
210 < }
211 <
212 < double Molecule::getCOMvel( double COMvel[3] ) {
213 <
214 <  double mass, mtot;
215 <  double aVel[3];
216 <  int i, j;
217 <
218 <
219 <  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 <
245 >    for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) {
246 >        potential += torsion->getPotential();
247      }
236  }
248  
249 <  for (j=0; j<3; j++)
239 <    COMvel[j] /= mtot;
240 <
241 <  return mtot;
249 >    return potential;
250  
251   }
252  
253 < double Molecule::getTotalMass()
254 < {
255 <
256 <  double totalMass;
257 <  
258 <  totalMass = 0;
259 <  for(int i =0; i < myIntegrableObjects.size(); i++){
260 <    totalMass += myIntegrableObjects[i]->getMass();
261 <  }
262 <
263 <  return totalMass;
253 > std::ostream& operator <<(std::ostream& o, Molecule& mol) {
254 >    o << std::endl;
255 >    o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl;
256 >    o << mol.getNAtoms() << " atoms" << std::endl;
257 >    o << mol.getNBonds() << " bonds" << std::endl;
258 >    o << mol.getNBends() << " bends" << std::endl;
259 >    o << mol.getNTorsions() << " torsions" << std::endl;
260 >    o << mol.getNRigidBodies() << " rigid bodies" << std::endl;
261 >    o << mol.getNIntegrableObjects() << "integrable objects" << std::endl;
262 >    o << mol.getNCutoffGroups() << "cutoff groups" << std::endl;
263 >    o << mol.getNConstraintPairs() << "constraint pairs" << std::endl;
264 >    return o;
265   }
266 +
267 + }//end namespace oopse

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines