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

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 292 by tim, Fri Feb 4 22:44:15 2005 UTC vs.
Revision 507 by gezelter, Fri Apr 15 22:04:00 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 56 | Line 56
56   #include "UseTheForce/notifyCutoffs_interface.h"
57   #include "utils/MemoryUtils.hpp"
58   #include "utils/simError.h"
59 + #include "selection/SelectionManager.hpp"
60  
61   #ifdef IS_MPI
62   #include "UseTheForce/mpiComponentPlan.h"
# Line 64 | 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 99 | 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 111 | 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(nGlobalAtoms_ + nGlobalRigidBodies_);
142 <    selectMan_->selectAll();
143 < }
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_;
154 < }
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 162 | 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 241 | 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 270 | 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 284 | 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 305 | 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 321 | 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 336 | 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 381 | 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 421 | 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 438 | 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 452 | 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 463 | 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 495 | Line 521 | void SimInfo::setupSimType() {
521  
522      //loop over all of the atom types
523      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
524 <        useLennardJones |= (*i)->isLennardJones();
525 <        useElectrostatic |= (*i)->isElectrostatic();
526 <        useEAM |= (*i)->isEAM();
527 <        useCharge |= (*i)->isCharge();
528 <        useDirectional |= (*i)->isDirectional();
529 <        useDipole |= (*i)->isDipole();
530 <        useGayBerne |= (*i)->isGayBerne();
531 <        useSticky |= (*i)->isSticky();
532 <        useShape |= (*i)->isShape();
524 >      useLennardJones |= (*i)->isLennardJones();
525 >      useElectrostatic |= (*i)->isElectrostatic();
526 >      useEAM |= (*i)->isEAM();
527 >      useCharge |= (*i)->isCharge();
528 >      useDirectional |= (*i)->isDirectional();
529 >      useDipole |= (*i)->isDipole();
530 >      useGayBerne |= (*i)->isGayBerne();
531 >      useSticky |= (*i)->isSticky();
532 >      useShape |= (*i)->isShape();
533      }
534  
535      if (useSticky || useDipole || useGayBerne || useShape) {
536 <        useDirectionalAtom = 1;
536 >      useDirectionalAtom = 1;
537      }
538  
539      if (useCharge || useDipole) {
540 <        useElectrostatics = 1;
540 >      useElectrostatics = 1;
541      }
542  
543   #ifdef IS_MPI    
# Line 570 | Line 596 | void SimInfo::setupSimType() {
596  
597      if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
598  
599 <        if (simParams_->haveDielectric()) {
600 <            fInfo_.dielect = simParams_->getDielectric();
601 <        } else {
602 <            sprintf(painCave.errMsg,
603 <                    "SimSetup Error: No Dielectric constant was set.\n"
604 <                    "\tYou are trying to use Reaction Field without"
605 <                    "\tsetting a dielectric constant!\n");
606 <            painCave.isFatal = 1;
607 <            simError();
608 <        }
599 >      if (simParams_->haveDielectric()) {
600 >        fInfo_.dielect = simParams_->getDielectric();
601 >      } else {
602 >        sprintf(painCave.errMsg,
603 >                "SimSetup Error: No Dielectric constant was set.\n"
604 >                "\tYou are trying to use Reaction Field without"
605 >                "\tsetting a dielectric constant!\n");
606 >        painCave.isFatal = 1;
607 >        simError();
608 >      }
609          
610      } else {
611 <        fInfo_.dielect = 0.0;
611 >      fInfo_.dielect = 0.0;
612      }
613  
614 < }
614 >  }
615  
616 < void SimInfo::setupFortranSim() {
616 >  void SimInfo::setupFortranSim() {
617      int isError;
618      int nExclude;
619      std::vector<int> fortranGlobalGroupMembership;
# Line 597 | Line 623 | void SimInfo::setupFortranSim() {
623  
624      //globalGroupMembership_ is filled by SimCreator    
625      for (int i = 0; i < nGlobalAtoms_; i++) {
626 <        fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
626 >      fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
627      }
628  
629      //calculate mass ratio of cutoff group
# Line 614 | Line 640 | void SimInfo::setupFortranSim() {
640      mfact.reserve(getNCutoffGroups());
641      
642      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
643 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
643 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
644  
645 <            totalMass = cg->getMass();
646 <            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
647 <                        mfact.push_back(atom->getMass()/totalMass);
648 <            }
645 >        totalMass = cg->getMass();
646 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
647 >          mfact.push_back(atom->getMass()/totalMass);
648 >        }
649  
650 <        }      
650 >      }      
651      }
652  
653      //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
# Line 631 | Line 657 | void SimInfo::setupFortranSim() {
657      identArray.reserve(getNAtoms());
658      
659      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
660 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
661 <            identArray.push_back(atom->getIdent());
662 <        }
660 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
661 >        identArray.push_back(atom->getIdent());
662 >      }
663      }    
664  
665      //fill molMembershipArray
666      //molMembershipArray is filled by SimCreator    
667      std::vector<int> molMembershipArray(nGlobalAtoms_);
668      for (int i = 0; i < nGlobalAtoms_; i++) {
669 <        molMembershipArray[i] = globalMolMembership_[i] + 1;
669 >      molMembershipArray[i] = globalMolMembership_[i] + 1;
670      }
671      
672      //setup fortran simulation
647 –    //gloalExcludes and molMembershipArray should go away (They are never used)
648 –    //why the hell fortran need to know molecule?
649 –    //OOPSE = Object-Obfuscated Parallel Simulation Engine
673      int nGlobalExcludes = 0;
674      int* globalExcludes = NULL;
675      int* excludeList = exclude_.getExcludeList();
676      setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
677 <                  &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
678 <                  &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
677 >                   &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
678 >                   &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
679  
680      if( isError ){
681  
682 <        sprintf( painCave.errMsg,
683 <                 "There was an error setting the simulation information in fortran.\n" );
684 <        painCave.isFatal = 1;
685 <        painCave.severity = OOPSE_ERROR;
686 <        simError();
682 >      sprintf( painCave.errMsg,
683 >               "There was an error setting the simulation information in fortran.\n" );
684 >      painCave.isFatal = 1;
685 >      painCave.severity = OOPSE_ERROR;
686 >      simError();
687      }
688  
689   #ifdef IS_MPI
690      sprintf( checkPointMsg,
691 <       "succesfully sent the simulation information to fortran.\n");
691 >             "succesfully sent the simulation information to fortran.\n");
692      MPIcheckPoint();
693   #endif // is_mpi
694 < }
694 >  }
695  
696  
697   #ifdef IS_MPI
698 < void SimInfo::setupFortranParallel() {
698 >  void SimInfo::setupFortranParallel() {
699      
700      //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
701      std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
# Line 688 | Line 711 | void SimInfo::setupFortranParallel() {
711  
712      for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
713  
714 <        //local index(index in DataStorge) of atom is important
715 <        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
716 <            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
717 <        }
714 >      //local index(index in DataStorge) of atom is important
715 >      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
716 >        localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
717 >      }
718  
719 <        //local index of cutoff group is trivial, it only depends on the order of travesing
720 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
721 <            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
722 <        }        
719 >      //local index of cutoff group is trivial, it only depends on the order of travesing
720 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
721 >        localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
722 >      }        
723          
724      }
725  
# Line 716 | Line 739 | void SimInfo::setupFortranParallel() {
739                      &localToGlobalCutoffGroupIndex[0], &isError);
740  
741      if (isError) {
742 <        sprintf(painCave.errMsg,
743 <                "mpiRefresh errror: fortran didn't like something we gave it.\n");
744 <        painCave.isFatal = 1;
745 <        simError();
742 >      sprintf(painCave.errMsg,
743 >              "mpiRefresh errror: fortran didn't like something we gave it.\n");
744 >      painCave.isFatal = 1;
745 >      simError();
746      }
747  
748      sprintf(checkPointMsg, " mpiRefresh successful.\n");
749      MPIcheckPoint();
750  
751  
752 < }
752 >  }
753  
754   #endif
755  
756 < double SimInfo::calcMaxCutoffRadius() {
756 >  double SimInfo::calcMaxCutoffRadius() {
757  
758  
759      std::set<AtomType*> atomTypes;
# Line 742 | Line 765 | double SimInfo::calcMaxCutoffRadius() {
765  
766      //query the max cutoff radius among these atom types
767      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
768 <        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
768 >      cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
769      }
770  
771      double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
# Line 751 | Line 774 | double SimInfo::calcMaxCutoffRadius() {
774   #endif
775  
776      return maxCutoffRadius;
777 < }
777 >  }
778  
779 < void SimInfo::setupCutoff() {
757 <    double rcut_;  //cutoff radius
758 <    double rsw_; //switching radius
779 >  void SimInfo::getCutoff(double& rcut, double& rsw) {
780      
781      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
782          
783 <        if (!simParams_->haveRcut()){
784 <            sprintf(painCave.errMsg,
783 >      if (!simParams_->haveRcut()){
784 >        sprintf(painCave.errMsg,
785                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
786                  "\tOOPSE will use a default value of 15.0 angstroms"
787                  "\tfor the cutoffRadius.\n");
788 <            painCave.isFatal = 0;
789 <            simError();
790 <            rcut_ = 15.0;
791 <        } else{
792 <            rcut_ = simParams_->getRcut();
793 <        }
788 >        painCave.isFatal = 0;
789 >        simError();
790 >        rcut = 15.0;
791 >      } else{
792 >        rcut = simParams_->getRcut();
793 >      }
794  
795 <        if (!simParams_->haveRsw()){
796 <            sprintf(painCave.errMsg,
795 >      if (!simParams_->haveRsw()){
796 >        sprintf(painCave.errMsg,
797                  "SimCreator Warning: No value was set for switchingRadius.\n"
798                  "\tOOPSE will use a default value of\n"
799                  "\t0.95 * cutoffRadius for the switchingRadius\n");
800 <            painCave.isFatal = 0;
801 <            simError();
802 <            rsw_ = 0.95 * rcut_;
803 <        } else{
804 <            rsw_ = simParams_->getRsw();
805 <        }
800 >        painCave.isFatal = 0;
801 >        simError();
802 >        rsw = 0.95 * rcut;
803 >      } else{
804 >        rsw = simParams_->getRsw();
805 >      }
806  
807      } else {
808 <        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
809 <        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
808 >      // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
809 >      //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
810          
811 <        if (simParams_->haveRcut()) {
812 <            rcut_ = simParams_->getRcut();
813 <        } else {
814 <            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
815 <            rcut_ = calcMaxCutoffRadius();
816 <        }
811 >      if (simParams_->haveRcut()) {
812 >        rcut = simParams_->getRcut();
813 >      } else {
814 >        //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
815 >        rcut = calcMaxCutoffRadius();
816 >      }
817  
818 <        if (simParams_->haveRsw()) {
819 <            rsw_  = simParams_->getRsw();
820 <        } else {
821 <            rsw_ = rcut_;
822 <        }
818 >      if (simParams_->haveRsw()) {
819 >        rsw  = simParams_->getRsw();
820 >      } else {
821 >        rsw = rcut;
822 >      }
823      
824      }
825 <        
825 >  }
826 >
827 >  void SimInfo::setupCutoff() {
828 >    getCutoff(rcut_, rsw_);    
829      double rnblist = rcut_ + 1; // skin of neighbor list
830  
831      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
832      notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
833 < }
833 >  }
834  
835 < void SimInfo::addProperty(GenericData* genData) {
835 >  void SimInfo::addProperty(GenericData* genData) {
836      properties_.addProperty(genData);  
837 < }
837 >  }
838  
839 < void SimInfo::removeProperty(const std::string& propName) {
839 >  void SimInfo::removeProperty(const std::string& propName) {
840      properties_.removeProperty(propName);  
841 < }
841 >  }
842  
843 < void SimInfo::clearProperties() {
843 >  void SimInfo::clearProperties() {
844      properties_.clearProperties();
845 < }
845 >  }
846  
847 < std::vector<std::string> SimInfo::getPropertyNames() {
847 >  std::vector<std::string> SimInfo::getPropertyNames() {
848      return properties_.getPropertyNames();  
849 < }
849 >  }
850        
851 < std::vector<GenericData*> SimInfo::getProperties() {
851 >  std::vector<GenericData*> SimInfo::getProperties() {
852      return properties_.getProperties();
853 < }
853 >  }
854  
855 < GenericData* SimInfo::getPropertyByName(const std::string& propName) {
855 >  GenericData* SimInfo::getPropertyByName(const std::string& propName) {
856      return properties_.getPropertyByName(propName);
857 < }
857 >  }
858  
859 < void SimInfo::setSnapshotManager(SnapshotManager* sman) {
859 >  void SimInfo::setSnapshotManager(SnapshotManager* sman) {
860 >    if (sman_ == sman) {
861 >      return;
862 >    }    
863 >    delete sman_;
864      sman_ = sman;
865  
866      Molecule* mol;
# Line 844 | Line 872 | void SimInfo::setSnapshotManager(SnapshotManager* sman
872  
873      for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
874          
875 <        for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
876 <            atom->setSnapshotManager(sman_);
877 <        }
875 >      for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
876 >        atom->setSnapshotManager(sman_);
877 >      }
878          
879 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
880 <            rb->setSnapshotManager(sman_);
881 <        }
879 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
880 >        rb->setSnapshotManager(sman_);
881 >      }
882      }    
883      
884 < }
884 >  }
885  
886 < Vector3d SimInfo::getComVel(){
886 >  Vector3d SimInfo::getComVel(){
887      SimInfo::MoleculeIterator i;
888      Molecule* mol;
889  
# Line 864 | Line 892 | Vector3d SimInfo::getComVel(){
892      
893  
894      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
895 <        double mass = mol->getMass();
896 <        totalMass += mass;
897 <        comVel += mass * mol->getComVel();
895 >      double mass = mol->getMass();
896 >      totalMass += mass;
897 >      comVel += mass * mol->getComVel();
898      }  
899  
900   #ifdef IS_MPI
# Line 879 | Line 907 | Vector3d SimInfo::getComVel(){
907      comVel /= totalMass;
908  
909      return comVel;
910 < }
910 >  }
911  
912 < Vector3d SimInfo::getCom(){
912 >  Vector3d SimInfo::getCom(){
913      SimInfo::MoleculeIterator i;
914      Molecule* mol;
915  
# Line 889 | Line 917 | Vector3d SimInfo::getCom(){
917      double totalMass = 0.0;
918      
919      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
920 <        double mass = mol->getMass();
921 <        totalMass += mass;
922 <        com += mass * mol->getCom();
920 >      double mass = mol->getMass();
921 >      totalMass += mass;
922 >      com += mass * mol->getCom();
923      }  
924  
925   #ifdef IS_MPI
# Line 905 | Line 933 | Vector3d SimInfo::getCom(){
933  
934      return com;
935  
936 < }        
936 >  }        
937  
938 < std::ostream& operator <<(std::ostream& o, SimInfo& info) {
938 >  std::ostream& operator <<(std::ostream& o, SimInfo& info) {
939  
940      return o;
941 < }
941 >  }
942  
943   }//end namespace oopse
944  

Diff Legend

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