ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/SimInfo.cpp
(Generate patch)

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 326 by tim, Sun Feb 13 20:05:42 2005 UTC vs.
Revision 523 by chrisfen, Thu May 5 14:47:35 2005 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 65 | Line 65 | namespace oopse {
65  
66   namespace oopse {
67  
68 < SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
69 <                                ForceField* ff, Globals* simParams) :
70 <                                forceField_(ff), simParams_(simParams),
71 <                                ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
72 <                                nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
73 <                                nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
74 <                                nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
75 <                                nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
76 <                                sman_(NULL), fortranInitialized_(false), selectMan_(NULL) {
68 >  SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
69 >                   ForceField* ff, Globals* simParams) :
70 >    stamps_(stamps), forceField_(ff), simParams_(simParams),
71 >    ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
72 >    nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
73 >    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
74 >    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
75 >    nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
76 >    sman_(NULL), fortranInitialized_(false) {
77  
78              
79 <    std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
80 <    MoleculeStamp* molStamp;
81 <    int nMolWithSameStamp;
82 <    int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
83 <    int nGroups = 0;          //total cutoff groups defined in meta-data file
84 <    CutoffGroupStamp* cgStamp;    
85 <    RigidBodyStamp* rbStamp;
86 <    int nRigidAtoms = 0;
79 >      std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
80 >      MoleculeStamp* molStamp;
81 >      int nMolWithSameStamp;
82 >      int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
83 >      int nGroups = 0;          //total cutoff groups defined in meta-data file
84 >      CutoffGroupStamp* cgStamp;    
85 >      RigidBodyStamp* rbStamp;
86 >      int nRigidAtoms = 0;
87      
88 <    for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
88 >      for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
89          molStamp = i->first;
90          nMolWithSameStamp = i->second;
91          
# Line 100 | Line 100 | SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*,
100          int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
101          
102          for (int j=0; j < nCutoffGroupsInStamp; j++) {
103 <            cgStamp = molStamp->getCutoffGroup(j);
104 <            nAtomsInGroups += cgStamp->getNMembers();
103 >          cgStamp = molStamp->getCutoffGroup(j);
104 >          nAtomsInGroups += cgStamp->getNMembers();
105          }
106  
107          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
# Line 112 | Line 112 | SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*,
112          int nRigidBodiesInStamp = molStamp->getNRigidBodies();
113          
114          for (int j=0; j < nRigidBodiesInStamp; j++) {
115 <            rbStamp = molStamp->getRigidBody(j);
116 <            nAtomsInRigidBodies += rbStamp->getNMembers();
115 >          rbStamp = molStamp->getRigidBody(j);
116 >          nAtomsInRigidBodies += rbStamp->getNMembers();
117          }
118  
119          nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
120          nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
121          
122 <    }
122 >      }
123  
124 <    //every free atom (atom does not belong to cutoff groups) is a cutoff group
125 <    //therefore the total number of cutoff groups in the system is equal to
126 <    //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
127 <    //file plus the number of cutoff groups defined in meta-data file
128 <    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
124 >      //every free atom (atom does not belong to cutoff groups) is a cutoff group
125 >      //therefore the total number of cutoff groups in the system is equal to
126 >      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
127 >      //file plus the number of cutoff groups defined in meta-data file
128 >      nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
129  
130 <    //every free atom (atom does not belong to rigid bodies) is an integrable object
131 <    //therefore the total number of  integrable objects in the system is equal to
132 <    //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
133 <    //file plus the number of  rigid bodies defined in meta-data file
134 <    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
130 >      //every free atom (atom does not belong to rigid bodies) is an integrable object
131 >      //therefore the total number of  integrable objects in the system is equal to
132 >      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
133 >      //file plus the number of  rigid bodies defined in meta-data file
134 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
135  
136 <    nGlobalMols_ = molStampIds_.size();
136 >      nGlobalMols_ = molStampIds_.size();
137  
138   #ifdef IS_MPI    
139 <    molToProcMap_.resize(nGlobalMols_);
139 >      molToProcMap_.resize(nGlobalMols_);
140   #endif
141  
142 <    selectMan_ = new SelectionManager(this);
143 <    selectMan_->selectAll();
144 < }
142 >    }
143  
144 < SimInfo::~SimInfo() {
145 <    //MemoryUtils::deleteVectorOfPointer(molecules_);
146 <
147 <    MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
148 <    
144 >  SimInfo::~SimInfo() {
145 >    std::map<int, Molecule*>::iterator i;
146 >    for (i = molecules_.begin(); i != molecules_.end(); ++i) {
147 >      delete i->second;
148 >    }
149 >    molecules_.clear();
150 >      
151 >    delete stamps_;
152      delete sman_;
153      delete simParams_;
154      delete forceField_;
155 <    delete selectMan_;
155 < }
155 >  }
156  
157 < int SimInfo::getNGlobalConstraints() {
157 >  int SimInfo::getNGlobalConstraints() {
158      int nGlobalConstraints;
159   #ifdef IS_MPI
160      MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
# Line 163 | Line 163 | int SimInfo::getNGlobalConstraints() {
163      nGlobalConstraints =  nConstraints_;
164   #endif
165      return nGlobalConstraints;
166 < }
166 >  }
167  
168 < bool SimInfo::addMolecule(Molecule* mol) {
168 >  bool SimInfo::addMolecule(Molecule* mol) {
169      MoleculeIterator i;
170  
171      i = molecules_.find(mol->getGlobalIndex());
172      if (i == molecules_.end() ) {
173  
174 <        molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
174 >      molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
175          
176 <        nAtoms_ += mol->getNAtoms();
177 <        nBonds_ += mol->getNBonds();
178 <        nBends_ += mol->getNBends();
179 <        nTorsions_ += mol->getNTorsions();
180 <        nRigidBodies_ += mol->getNRigidBodies();
181 <        nIntegrableObjects_ += mol->getNIntegrableObjects();
182 <        nCutoffGroups_ += mol->getNCutoffGroups();
183 <        nConstraints_ += mol->getNConstraintPairs();
184 <
185 <        addExcludePairs(mol);
176 >      nAtoms_ += mol->getNAtoms();
177 >      nBonds_ += mol->getNBonds();
178 >      nBends_ += mol->getNBends();
179 >      nTorsions_ += mol->getNTorsions();
180 >      nRigidBodies_ += mol->getNRigidBodies();
181 >      nIntegrableObjects_ += mol->getNIntegrableObjects();
182 >      nCutoffGroups_ += mol->getNCutoffGroups();
183 >      nConstraints_ += mol->getNConstraintPairs();
184 >
185 >      addExcludePairs(mol);
186          
187 <        return true;
187 >      return true;
188      } else {
189 <        return false;
189 >      return false;
190      }
191 < }
191 >  }
192  
193 < bool SimInfo::removeMolecule(Molecule* mol) {
193 >  bool SimInfo::removeMolecule(Molecule* mol) {
194      MoleculeIterator i;
195      i = molecules_.find(mol->getGlobalIndex());
196  
197      if (i != molecules_.end() ) {
198  
199 <        assert(mol == i->second);
199 >      assert(mol == i->second);
200          
201 <        nAtoms_ -= mol->getNAtoms();
202 <        nBonds_ -= mol->getNBonds();
203 <        nBends_ -= mol->getNBends();
204 <        nTorsions_ -= mol->getNTorsions();
205 <        nRigidBodies_ -= mol->getNRigidBodies();
206 <        nIntegrableObjects_ -= mol->getNIntegrableObjects();
207 <        nCutoffGroups_ -= mol->getNCutoffGroups();
208 <        nConstraints_ -= mol->getNConstraintPairs();
201 >      nAtoms_ -= mol->getNAtoms();
202 >      nBonds_ -= mol->getNBonds();
203 >      nBends_ -= mol->getNBends();
204 >      nTorsions_ -= mol->getNTorsions();
205 >      nRigidBodies_ -= mol->getNRigidBodies();
206 >      nIntegrableObjects_ -= mol->getNIntegrableObjects();
207 >      nCutoffGroups_ -= mol->getNCutoffGroups();
208 >      nConstraints_ -= mol->getNConstraintPairs();
209  
210 <        removeExcludePairs(mol);
211 <        molecules_.erase(mol->getGlobalIndex());
210 >      removeExcludePairs(mol);
211 >      molecules_.erase(mol->getGlobalIndex());
212  
213 <        delete mol;
213 >      delete mol;
214          
215 <        return true;
215 >      return true;
216      } else {
217 <        return false;
217 >      return false;
218      }
219  
220  
221 < }    
221 >  }    
222  
223          
224 < Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
224 >  Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
225      i = molecules_.begin();
226      return i == molecules_.end() ? NULL : i->second;
227 < }    
227 >  }    
228  
229 < Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
229 >  Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
230      ++i;
231      return i == molecules_.end() ? NULL : i->second;    
232 < }
232 >  }
233  
234  
235 < void SimInfo::calcNdf() {
235 >  void SimInfo::calcNdf() {
236      int ndf_local;
237      MoleculeIterator i;
238      std::vector<StuntDouble*>::iterator j;
# Line 242 | Line 242 | void SimInfo::calcNdf() {
242      ndf_local = 0;
243      
244      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
245 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
246 <               integrableObject = mol->nextIntegrableObject(j)) {
245 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
246 >           integrableObject = mol->nextIntegrableObject(j)) {
247  
248 <            ndf_local += 3;
248 >        ndf_local += 3;
249  
250 <            if (integrableObject->isDirectional()) {
251 <                if (integrableObject->isLinear()) {
252 <                    ndf_local += 2;
253 <                } else {
254 <                    ndf_local += 3;
255 <                }
256 <            }
250 >        if (integrableObject->isDirectional()) {
251 >          if (integrableObject->isLinear()) {
252 >            ndf_local += 2;
253 >          } else {
254 >            ndf_local += 3;
255 >          }
256 >        }
257              
258 <        }//end for (integrableObject)
258 >      }//end for (integrableObject)
259      }// end for (mol)
260      
261      // n_constraints is local, so subtract them on each processor
# Line 271 | Line 271 | void SimInfo::calcNdf() {
271      // entire system:
272      ndf_ = ndf_ - 3 - nZconstraint_;
273  
274 < }
274 >  }
275  
276 < void SimInfo::calcNdfRaw() {
276 >  void SimInfo::calcNdfRaw() {
277      int ndfRaw_local;
278  
279      MoleculeIterator i;
# Line 285 | Line 285 | void SimInfo::calcNdfRaw() {
285      ndfRaw_local = 0;
286      
287      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
288 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
289 <               integrableObject = mol->nextIntegrableObject(j)) {
288 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
289 >           integrableObject = mol->nextIntegrableObject(j)) {
290  
291 <            ndfRaw_local += 3;
291 >        ndfRaw_local += 3;
292  
293 <            if (integrableObject->isDirectional()) {
294 <                if (integrableObject->isLinear()) {
295 <                    ndfRaw_local += 2;
296 <                } else {
297 <                    ndfRaw_local += 3;
298 <                }
299 <            }
293 >        if (integrableObject->isDirectional()) {
294 >          if (integrableObject->isLinear()) {
295 >            ndfRaw_local += 2;
296 >          } else {
297 >            ndfRaw_local += 3;
298 >          }
299 >        }
300              
301 <        }
301 >      }
302      }
303      
304   #ifdef IS_MPI
# Line 306 | Line 306 | void SimInfo::calcNdfRaw() {
306   #else
307      ndfRaw_ = ndfRaw_local;
308   #endif
309 < }
309 >  }
310  
311 < void SimInfo::calcNdfTrans() {
311 >  void SimInfo::calcNdfTrans() {
312      int ndfTrans_local;
313  
314      ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
# Line 322 | Line 322 | void SimInfo::calcNdfTrans() {
322  
323      ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
324  
325 < }
325 >  }
326  
327 < void SimInfo::addExcludePairs(Molecule* mol) {
327 >  void SimInfo::addExcludePairs(Molecule* mol) {
328      std::vector<Bond*>::iterator bondIter;
329      std::vector<Bend*>::iterator bendIter;
330      std::vector<Torsion*>::iterator torsionIter;
# Line 337 | Line 337 | void SimInfo::addExcludePairs(Molecule* mol) {
337      int d;
338      
339      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
340 <        a = bond->getAtomA()->getGlobalIndex();
341 <        b = bond->getAtomB()->getGlobalIndex();        
342 <        exclude_.addPair(a, b);
340 >      a = bond->getAtomA()->getGlobalIndex();
341 >      b = bond->getAtomB()->getGlobalIndex();        
342 >      exclude_.addPair(a, b);
343      }
344  
345      for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
346 <        a = bend->getAtomA()->getGlobalIndex();
347 <        b = bend->getAtomB()->getGlobalIndex();        
348 <        c = bend->getAtomC()->getGlobalIndex();
346 >      a = bend->getAtomA()->getGlobalIndex();
347 >      b = bend->getAtomB()->getGlobalIndex();        
348 >      c = bend->getAtomC()->getGlobalIndex();
349  
350 <        exclude_.addPair(a, b);
351 <        exclude_.addPair(a, c);
352 <        exclude_.addPair(b, c);        
350 >      exclude_.addPair(a, b);
351 >      exclude_.addPair(a, c);
352 >      exclude_.addPair(b, c);        
353      }
354  
355      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
356 <        a = torsion->getAtomA()->getGlobalIndex();
357 <        b = torsion->getAtomB()->getGlobalIndex();        
358 <        c = torsion->getAtomC()->getGlobalIndex();        
359 <        d = torsion->getAtomD()->getGlobalIndex();        
356 >      a = torsion->getAtomA()->getGlobalIndex();
357 >      b = torsion->getAtomB()->getGlobalIndex();        
358 >      c = torsion->getAtomC()->getGlobalIndex();        
359 >      d = torsion->getAtomD()->getGlobalIndex();        
360  
361 <        exclude_.addPair(a, b);
362 <        exclude_.addPair(a, c);
363 <        exclude_.addPair(a, d);
364 <        exclude_.addPair(b, c);
365 <        exclude_.addPair(b, d);
366 <        exclude_.addPair(c, d);        
361 >      exclude_.addPair(a, b);
362 >      exclude_.addPair(a, c);
363 >      exclude_.addPair(a, d);
364 >      exclude_.addPair(b, c);
365 >      exclude_.addPair(b, d);
366 >      exclude_.addPair(c, d);        
367      }
368  
369 <    
370 < }
369 >    Molecule::RigidBodyIterator rbIter;
370 >    RigidBody* rb;
371 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
372 >      std::vector<Atom*> atoms = rb->getAtoms();
373 >      for (int i = 0; i < atoms.size() -1 ; ++i) {
374 >        for (int j = i + 1; j < atoms.size(); ++j) {
375 >          a = atoms[i]->getGlobalIndex();
376 >          b = atoms[j]->getGlobalIndex();
377 >          exclude_.addPair(a, b);
378 >        }
379 >      }
380 >    }        
381  
382 < void SimInfo::removeExcludePairs(Molecule* mol) {
382 >  }
383 >
384 >  void SimInfo::removeExcludePairs(Molecule* mol) {
385      std::vector<Bond*>::iterator bondIter;
386      std::vector<Bend*>::iterator bendIter;
387      std::vector<Torsion*>::iterator torsionIter;
# Line 382 | Line 394 | void SimInfo::removeExcludePairs(Molecule* mol) {
394      int d;
395      
396      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
397 <        a = bond->getAtomA()->getGlobalIndex();
398 <        b = bond->getAtomB()->getGlobalIndex();        
399 <        exclude_.removePair(a, b);
397 >      a = bond->getAtomA()->getGlobalIndex();
398 >      b = bond->getAtomB()->getGlobalIndex();        
399 >      exclude_.removePair(a, b);
400      }
401  
402      for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
403 <        a = bend->getAtomA()->getGlobalIndex();
404 <        b = bend->getAtomB()->getGlobalIndex();        
405 <        c = bend->getAtomC()->getGlobalIndex();
403 >      a = bend->getAtomA()->getGlobalIndex();
404 >      b = bend->getAtomB()->getGlobalIndex();        
405 >      c = bend->getAtomC()->getGlobalIndex();
406  
407 <        exclude_.removePair(a, b);
408 <        exclude_.removePair(a, c);
409 <        exclude_.removePair(b, c);        
407 >      exclude_.removePair(a, b);
408 >      exclude_.removePair(a, c);
409 >      exclude_.removePair(b, c);        
410      }
411  
412      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
413 <        a = torsion->getAtomA()->getGlobalIndex();
414 <        b = torsion->getAtomB()->getGlobalIndex();        
415 <        c = torsion->getAtomC()->getGlobalIndex();        
416 <        d = torsion->getAtomD()->getGlobalIndex();        
413 >      a = torsion->getAtomA()->getGlobalIndex();
414 >      b = torsion->getAtomB()->getGlobalIndex();        
415 >      c = torsion->getAtomC()->getGlobalIndex();        
416 >      d = torsion->getAtomD()->getGlobalIndex();        
417  
418 <        exclude_.removePair(a, b);
419 <        exclude_.removePair(a, c);
420 <        exclude_.removePair(a, d);
421 <        exclude_.removePair(b, c);
422 <        exclude_.removePair(b, d);
423 <        exclude_.removePair(c, d);        
418 >      exclude_.removePair(a, b);
419 >      exclude_.removePair(a, c);
420 >      exclude_.removePair(a, d);
421 >      exclude_.removePair(b, c);
422 >      exclude_.removePair(b, d);
423 >      exclude_.removePair(c, d);        
424      }
425  
426 < }
426 >    Molecule::RigidBodyIterator rbIter;
427 >    RigidBody* rb;
428 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
429 >      std::vector<Atom*> atoms = rb->getAtoms();
430 >      for (int i = 0; i < atoms.size() -1 ; ++i) {
431 >        for (int j = i + 1; j < atoms.size(); ++j) {
432 >          a = atoms[i]->getGlobalIndex();
433 >          b = atoms[j]->getGlobalIndex();
434 >          exclude_.removePair(a, b);
435 >        }
436 >      }
437 >    }        
438  
439 +  }
440  
441 < void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
441 >
442 >  void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
443      int curStampId;
444  
445      //index from 0
# Line 422 | Line 447 | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp
447  
448      moleculeStamps_.push_back(molStamp);
449      molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
450 < }
450 >  }
451  
452 < void SimInfo::update() {
452 >  void SimInfo::update() {
453  
454      setupSimType();
455  
# Line 439 | Line 464 | void SimInfo::update() {
464      int isError = 0;
465      initFortranFF( &fInfo_.SIM_uses_RF , &isError );
466      if(isError){
467 <        sprintf( painCave.errMsg,
468 <         "ForceField error: There was an error initializing the forceField in fortran.\n" );
469 <        painCave.isFatal = 1;
470 <        simError();
467 >      sprintf( painCave.errMsg,
468 >               "ForceField error: There was an error initializing the forceField in fortran.\n" );
469 >      painCave.isFatal = 1;
470 >      simError();
471      }
472    
473      
# Line 453 | Line 478 | void SimInfo::update() {
478      calcNdfTrans();
479  
480      fortranInitialized_ = true;
481 < }
481 >  }
482  
483 < std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
483 >  std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
484      SimInfo::MoleculeIterator mi;
485      Molecule* mol;
486      Molecule::AtomIterator ai;
# Line 464 | Line 489 | std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
489  
490      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
491  
492 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
493 <            atomTypes.insert(atom->getAtomType());
494 <        }
492 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
493 >        atomTypes.insert(atom->getAtomType());
494 >      }
495          
496      }
497  
498      return atomTypes;        
499 < }
499 >  }
500  
501 < void SimInfo::setupSimType() {
501 >  void SimInfo::setupSimType() {
502      std::set<AtomType*>::iterator i;
503      std::set<AtomType*> atomTypes;
504      atomTypes = getUniqueAtomTypes();
# Line 486 | Line 511 | void SimInfo::setupSimType() {
511      int useDipole = 0;
512      int useGayBerne = 0;
513      int useSticky = 0;
514 +    int useStickyPower = 0;
515      int useShape = 0;
516      int useFLARB = 0; //it is not in AtomType yet
517      int useDirectionalAtom = 0;    
# Line 496 | Line 522 | void SimInfo::setupSimType() {
522  
523      //loop over all of the atom types
524      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
525 <        useLennardJones |= (*i)->isLennardJones();
526 <        useElectrostatic |= (*i)->isElectrostatic();
527 <        useEAM |= (*i)->isEAM();
528 <        useCharge |= (*i)->isCharge();
529 <        useDirectional |= (*i)->isDirectional();
530 <        useDipole |= (*i)->isDipole();
531 <        useGayBerne |= (*i)->isGayBerne();
532 <        useSticky |= (*i)->isSticky();
533 <        useShape |= (*i)->isShape();
525 >      useLennardJones |= (*i)->isLennardJones();
526 >      useElectrostatic |= (*i)->isElectrostatic();
527 >      useEAM |= (*i)->isEAM();
528 >      useCharge |= (*i)->isCharge();
529 >      useDirectional |= (*i)->isDirectional();
530 >      useDipole |= (*i)->isDipole();
531 >      useGayBerne |= (*i)->isGayBerne();
532 >      useSticky |= (*i)->isSticky();
533 >      useStickyPower |= (*i)->isStickyPower();
534 >      useShape |= (*i)->isShape();
535      }
536  
537 <    if (useSticky || useDipole || useGayBerne || useShape) {
538 <        useDirectionalAtom = 1;
537 >    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
538 >      useDirectionalAtom = 1;
539      }
540  
541      if (useCharge || useDipole) {
542 <        useElectrostatics = 1;
542 >      useElectrostatics = 1;
543      }
544  
545   #ifdef IS_MPI    
# Line 539 | Line 566 | void SimInfo::setupSimType() {
566      temp = useSticky;
567      MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
568  
569 +    temp = useStickyPower;
570 +    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
571 +    
572      temp = useGayBerne;
573      MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
574  
# Line 563 | Line 593 | void SimInfo::setupSimType() {
593      fInfo_.SIM_uses_Charges = useCharge;
594      fInfo_.SIM_uses_Dipoles = useDipole;
595      fInfo_.SIM_uses_Sticky = useSticky;
596 +    fInfo_.SIM_uses_StickyPower = useStickyPower;
597      fInfo_.SIM_uses_GayBerne = useGayBerne;
598      fInfo_.SIM_uses_EAM = useEAM;
599      fInfo_.SIM_uses_Shapes = useShape;
# Line 571 | Line 602 | void SimInfo::setupSimType() {
602  
603      if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
604  
605 <        if (simParams_->haveDielectric()) {
606 <            fInfo_.dielect = simParams_->getDielectric();
607 <        } else {
608 <            sprintf(painCave.errMsg,
609 <                    "SimSetup Error: No Dielectric constant was set.\n"
610 <                    "\tYou are trying to use Reaction Field without"
611 <                    "\tsetting a dielectric constant!\n");
612 <            painCave.isFatal = 1;
613 <            simError();
614 <        }
605 >      if (simParams_->haveDielectric()) {
606 >        fInfo_.dielect = simParams_->getDielectric();
607 >      } else {
608 >        sprintf(painCave.errMsg,
609 >                "SimSetup Error: No Dielectric constant was set.\n"
610 >                "\tYou are trying to use Reaction Field without"
611 >                "\tsetting a dielectric constant!\n");
612 >        painCave.isFatal = 1;
613 >        simError();
614 >      }
615          
616      } else {
617 <        fInfo_.dielect = 0.0;
617 >      fInfo_.dielect = 0.0;
618      }
619  
620 < }
620 >  }
621  
622 < void SimInfo::setupFortranSim() {
622 >  void SimInfo::setupFortranSim() {
623      int isError;
624      int nExclude;
625      std::vector<int> fortranGlobalGroupMembership;
# Line 598 | Line 629 | void SimInfo::setupFortranSim() {
629  
630      //globalGroupMembership_ is filled by SimCreator    
631      for (int i = 0; i < nGlobalAtoms_; i++) {
632 <        fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
632 >      fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
633      }
634  
635      //calculate mass ratio of cutoff group
# Line 615 | Line 646 | void SimInfo::setupFortranSim() {
646      mfact.reserve(getNCutoffGroups());
647      
648      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
649 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
649 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
650  
651 <            totalMass = cg->getMass();
652 <            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
653 <                        mfact.push_back(atom->getMass()/totalMass);
654 <            }
651 >        totalMass = cg->getMass();
652 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
653 >          mfact.push_back(atom->getMass()/totalMass);
654 >        }
655  
656 <        }      
656 >      }      
657      }
658  
659      //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
# Line 632 | Line 663 | void SimInfo::setupFortranSim() {
663      identArray.reserve(getNAtoms());
664      
665      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
666 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
667 <            identArray.push_back(atom->getIdent());
668 <        }
666 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
667 >        identArray.push_back(atom->getIdent());
668 >      }
669      }    
670  
671      //fill molMembershipArray
672      //molMembershipArray is filled by SimCreator    
673      std::vector<int> molMembershipArray(nGlobalAtoms_);
674      for (int i = 0; i < nGlobalAtoms_; i++) {
675 <        molMembershipArray[i] = globalMolMembership_[i] + 1;
675 >      molMembershipArray[i] = globalMolMembership_[i] + 1;
676      }
677      
678      //setup fortran simulation
648 –    //gloalExcludes and molMembershipArray should go away (They are never used)
649 –    //why the hell fortran need to know molecule?
650 –    //OOPSE = Object-Obfuscated Parallel Simulation Engine
679      int nGlobalExcludes = 0;
680      int* globalExcludes = NULL;
681      int* excludeList = exclude_.getExcludeList();
682      setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
683 <                  &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
684 <                  &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
683 >                   &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
684 >                   &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
685  
686      if( isError ){
687  
688 <        sprintf( painCave.errMsg,
689 <                 "There was an error setting the simulation information in fortran.\n" );
690 <        painCave.isFatal = 1;
691 <        painCave.severity = OOPSE_ERROR;
692 <        simError();
688 >      sprintf( painCave.errMsg,
689 >               "There was an error setting the simulation information in fortran.\n" );
690 >      painCave.isFatal = 1;
691 >      painCave.severity = OOPSE_ERROR;
692 >      simError();
693      }
694  
695   #ifdef IS_MPI
696      sprintf( checkPointMsg,
697 <       "succesfully sent the simulation information to fortran.\n");
697 >             "succesfully sent the simulation information to fortran.\n");
698      MPIcheckPoint();
699   #endif // is_mpi
700 < }
700 >  }
701  
702  
703   #ifdef IS_MPI
704 < void SimInfo::setupFortranParallel() {
704 >  void SimInfo::setupFortranParallel() {
705      
706      //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
707      std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
# Line 689 | Line 717 | void SimInfo::setupFortranParallel() {
717  
718      for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
719  
720 <        //local index(index in DataStorge) of atom is important
721 <        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
722 <            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
723 <        }
720 >      //local index(index in DataStorge) of atom is important
721 >      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
722 >        localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
723 >      }
724  
725 <        //local index of cutoff group is trivial, it only depends on the order of travesing
726 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
727 <            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
728 <        }        
725 >      //local index of cutoff group is trivial, it only depends on the order of travesing
726 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
727 >        localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
728 >      }        
729          
730      }
731  
# Line 717 | Line 745 | void SimInfo::setupFortranParallel() {
745                      &localToGlobalCutoffGroupIndex[0], &isError);
746  
747      if (isError) {
748 <        sprintf(painCave.errMsg,
749 <                "mpiRefresh errror: fortran didn't like something we gave it.\n");
750 <        painCave.isFatal = 1;
751 <        simError();
748 >      sprintf(painCave.errMsg,
749 >              "mpiRefresh errror: fortran didn't like something we gave it.\n");
750 >      painCave.isFatal = 1;
751 >      simError();
752      }
753  
754      sprintf(checkPointMsg, " mpiRefresh successful.\n");
755      MPIcheckPoint();
756  
757  
758 < }
758 >  }
759  
760   #endif
761  
762 < double SimInfo::calcMaxCutoffRadius() {
762 >  double SimInfo::calcMaxCutoffRadius() {
763  
764  
765      std::set<AtomType*> atomTypes;
# Line 743 | Line 771 | double SimInfo::calcMaxCutoffRadius() {
771  
772      //query the max cutoff radius among these atom types
773      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
774 <        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
774 >      cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
775      }
776  
777      double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
# Line 752 | Line 780 | double SimInfo::calcMaxCutoffRadius() {
780   #endif
781  
782      return maxCutoffRadius;
783 < }
783 >  }
784  
785 < void SimInfo::getCutoff(double& rcut, double& rsw) {
785 >  void SimInfo::getCutoff(double& rcut, double& rsw) {
786      
787      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
788          
789 <        if (!simParams_->haveRcut()){
790 <            sprintf(painCave.errMsg,
789 >      if (!simParams_->haveRcut()){
790 >        sprintf(painCave.errMsg,
791                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
792                  "\tOOPSE will use a default value of 15.0 angstroms"
793                  "\tfor the cutoffRadius.\n");
794 <            painCave.isFatal = 0;
795 <            simError();
796 <            rcut_ = 15.0;
797 <        } else{
798 <            rcut_ = simParams_->getRcut();
799 <        }
794 >        painCave.isFatal = 0;
795 >        simError();
796 >        rcut = 15.0;
797 >      } else{
798 >        rcut = simParams_->getRcut();
799 >      }
800  
801 <        if (!simParams_->haveRsw()){
802 <            sprintf(painCave.errMsg,
801 >      if (!simParams_->haveRsw()){
802 >        sprintf(painCave.errMsg,
803                  "SimCreator Warning: No value was set for switchingRadius.\n"
804                  "\tOOPSE will use a default value of\n"
805                  "\t0.95 * cutoffRadius for the switchingRadius\n");
806 <            painCave.isFatal = 0;
807 <            simError();
808 <            rsw_ = 0.95 * rcut_;
809 <        } else{
810 <            rsw_ = simParams_->getRsw();
811 <        }
806 >        painCave.isFatal = 0;
807 >        simError();
808 >        rsw = 0.95 * rcut;
809 >      } else{
810 >        rsw = simParams_->getRsw();
811 >      }
812  
813      } else {
814 <        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
815 <        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
814 >      // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
815 >      //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
816          
817 <        if (simParams_->haveRcut()) {
818 <            rcut_ = simParams_->getRcut();
819 <        } else {
820 <            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
821 <            rcut_ = calcMaxCutoffRadius();
822 <        }
817 >      if (simParams_->haveRcut()) {
818 >        rcut = simParams_->getRcut();
819 >      } else {
820 >        //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
821 >        rcut = calcMaxCutoffRadius();
822 >      }
823  
824 <        if (simParams_->haveRsw()) {
825 <            rsw_  = simParams_->getRsw();
826 <        } else {
827 <            rsw_ = rcut_;
828 <        }
824 >      if (simParams_->haveRsw()) {
825 >        rsw  = simParams_->getRsw();
826 >      } else {
827 >        rsw = rcut;
828 >      }
829      
830      }
831 < }
831 >  }
832  
833 < void SimInfo::setupCutoff() {
833 >  void SimInfo::setupCutoff() {
834      getCutoff(rcut_, rsw_);    
835      double rnblist = rcut_ + 1; // skin of neighbor list
836  
837      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
838      notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
839 < }
839 >  }
840  
841 < void SimInfo::addProperty(GenericData* genData) {
841 >  void SimInfo::addProperty(GenericData* genData) {
842      properties_.addProperty(genData);  
843 < }
843 >  }
844  
845 < void SimInfo::removeProperty(const std::string& propName) {
845 >  void SimInfo::removeProperty(const std::string& propName) {
846      properties_.removeProperty(propName);  
847 < }
847 >  }
848  
849 < void SimInfo::clearProperties() {
849 >  void SimInfo::clearProperties() {
850      properties_.clearProperties();
851 < }
851 >  }
852  
853 < std::vector<std::string> SimInfo::getPropertyNames() {
853 >  std::vector<std::string> SimInfo::getPropertyNames() {
854      return properties_.getPropertyNames();  
855 < }
855 >  }
856        
857 < std::vector<GenericData*> SimInfo::getProperties() {
857 >  std::vector<GenericData*> SimInfo::getProperties() {
858      return properties_.getProperties();
859 < }
859 >  }
860  
861 < GenericData* SimInfo::getPropertyByName(const std::string& propName) {
861 >  GenericData* SimInfo::getPropertyByName(const std::string& propName) {
862      return properties_.getPropertyByName(propName);
863 < }
863 >  }
864  
865 < void SimInfo::setSnapshotManager(SnapshotManager* sman) {
865 >  void SimInfo::setSnapshotManager(SnapshotManager* sman) {
866 >    if (sman_ == sman) {
867 >      return;
868 >    }    
869 >    delete sman_;
870      sman_ = sman;
871  
872      Molecule* mol;
# Line 846 | Line 878 | void SimInfo::setSnapshotManager(SnapshotManager* sman
878  
879      for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
880          
881 <        for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
882 <            atom->setSnapshotManager(sman_);
883 <        }
881 >      for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
882 >        atom->setSnapshotManager(sman_);
883 >      }
884          
885 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
886 <            rb->setSnapshotManager(sman_);
887 <        }
885 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
886 >        rb->setSnapshotManager(sman_);
887 >      }
888      }    
889      
890 < }
890 >  }
891  
892 < Vector3d SimInfo::getComVel(){
892 >  Vector3d SimInfo::getComVel(){
893      SimInfo::MoleculeIterator i;
894      Molecule* mol;
895  
# Line 866 | Line 898 | Vector3d SimInfo::getComVel(){
898      
899  
900      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
901 <        double mass = mol->getMass();
902 <        totalMass += mass;
903 <        comVel += mass * mol->getComVel();
901 >      double mass = mol->getMass();
902 >      totalMass += mass;
903 >      comVel += mass * mol->getComVel();
904      }  
905  
906   #ifdef IS_MPI
# Line 881 | Line 913 | Vector3d SimInfo::getComVel(){
913      comVel /= totalMass;
914  
915      return comVel;
916 < }
916 >  }
917  
918 < Vector3d SimInfo::getCom(){
918 >  Vector3d SimInfo::getCom(){
919      SimInfo::MoleculeIterator i;
920      Molecule* mol;
921  
# Line 891 | Line 923 | Vector3d SimInfo::getCom(){
923      double totalMass = 0.0;
924      
925      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
926 <        double mass = mol->getMass();
927 <        totalMass += mass;
928 <        com += mass * mol->getCom();
926 >      double mass = mol->getMass();
927 >      totalMass += mass;
928 >      com += mass * mol->getCom();
929      }  
930  
931   #ifdef IS_MPI
# Line 907 | Line 939 | Vector3d SimInfo::getCom(){
939  
940      return com;
941  
942 < }        
942 >  }        
943  
944 < std::ostream& operator <<(std::ostream& o, SimInfo& info) {
944 >  std::ostream& operator <<(std::ostream& o, SimInfo& info) {
945  
946      return o;
947 < }
947 >  }
948  
949   }//end namespace oopse
950  

Diff Legend

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