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 490 by tim, Wed Apr 13 18:41:17 2005 UTC vs.
Revision 598 by chrisfen, Thu Sep 15 00:14: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 52 | Line 52
52   #include "brains/SimInfo.hpp"
53   #include "math/Vector3.hpp"
54   #include "primitives/Molecule.hpp"
55 + #include "UseTheForce/fCutoffPolicy.h"
56 + #include "UseTheForce/fCoulombicCorrection.h"
57   #include "UseTheForce/doForces_interface.h"
58   #include "UseTheForce/notifyCutoffs_interface.h"
59   #include "utils/MemoryUtils.hpp"
# Line 65 | Line 67 | namespace oopse {
67  
68   namespace oopse {
69  
70 < SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
71 <                                ForceField* ff, Globals* simParams) :
72 <                                stamps_(stamps), forceField_(ff), simParams_(simParams),
73 <                                ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
74 <                                nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
75 <                                nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
76 <                                nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
77 <                                nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
78 <                                sman_(NULL), fortranInitialized_(false) {
70 >  SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
71 >                   ForceField* ff, Globals* simParams) :
72 >    stamps_(stamps), forceField_(ff), simParams_(simParams),
73 >    ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
74 >    nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
75 >    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
76 >    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
77 >    nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
78 >    sman_(NULL), fortranInitialized_(false) {
79  
80              
81 <    std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
82 <    MoleculeStamp* molStamp;
83 <    int nMolWithSameStamp;
84 <    int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
85 <    int nGroups = 0;          //total cutoff groups defined in meta-data file
86 <    CutoffGroupStamp* cgStamp;    
87 <    RigidBodyStamp* rbStamp;
88 <    int nRigidAtoms = 0;
81 >      std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
82 >      MoleculeStamp* molStamp;
83 >      int nMolWithSameStamp;
84 >      int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
85 >      int nGroups = 0;          //total cutoff groups defined in meta-data file
86 >      CutoffGroupStamp* cgStamp;    
87 >      RigidBodyStamp* rbStamp;
88 >      int nRigidAtoms = 0;
89      
90 <    for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
90 >      for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
91          molStamp = i->first;
92          nMolWithSameStamp = i->second;
93          
# Line 100 | Line 102 | SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::
102          int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
103          
104          for (int j=0; j < nCutoffGroupsInStamp; j++) {
105 <            cgStamp = molStamp->getCutoffGroup(j);
106 <            nAtomsInGroups += cgStamp->getNMembers();
105 >          cgStamp = molStamp->getCutoffGroup(j);
106 >          nAtomsInGroups += cgStamp->getNMembers();
107          }
108  
109          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
# Line 112 | Line 114 | SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::
114          int nRigidBodiesInStamp = molStamp->getNRigidBodies();
115          
116          for (int j=0; j < nRigidBodiesInStamp; j++) {
117 <            rbStamp = molStamp->getRigidBody(j);
118 <            nAtomsInRigidBodies += rbStamp->getNMembers();
117 >          rbStamp = molStamp->getRigidBody(j);
118 >          nAtomsInRigidBodies += rbStamp->getNMembers();
119          }
120  
121          nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
122          nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
123          
124 <    }
124 >      }
125  
126 <    //every free atom (atom does not belong to cutoff groups) is a cutoff group
127 <    //therefore the total number of cutoff groups in the system is equal to
128 <    //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
129 <    //file plus the number of cutoff groups defined in meta-data file
130 <    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
126 >      //every free atom (atom does not belong to cutoff groups) is a cutoff group
127 >      //therefore the total number of cutoff groups in the system is equal to
128 >      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
129 >      //file plus the number of cutoff groups defined in meta-data file
130 >      nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
131  
132 <    //every free atom (atom does not belong to rigid bodies) is an integrable object
133 <    //therefore the total number of  integrable objects in the system is equal to
134 <    //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
135 <    //file plus the number of  rigid bodies defined in meta-data file
136 <    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
132 >      //every free atom (atom does not belong to rigid bodies) is an integrable object
133 >      //therefore the total number of  integrable objects in the system is equal to
134 >      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
135 >      //file plus the number of  rigid bodies defined in meta-data file
136 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
137  
138 <    nGlobalMols_ = molStampIds_.size();
138 >      nGlobalMols_ = molStampIds_.size();
139  
140   #ifdef IS_MPI    
141 <    molToProcMap_.resize(nGlobalMols_);
141 >      molToProcMap_.resize(nGlobalMols_);
142   #endif
143  
144 < }
144 >    }
145  
146 < SimInfo::~SimInfo() {
146 >  SimInfo::~SimInfo() {
147      std::map<int, Molecule*>::iterator i;
148      for (i = molecules_.begin(); i != molecules_.end(); ++i) {
149 <        delete i->second;
149 >      delete i->second;
150      }
151      molecules_.clear();
152        
# Line 152 | Line 154 | SimInfo::~SimInfo() {
154      delete sman_;
155      delete simParams_;
156      delete forceField_;
157 < }
157 >  }
158  
159 < int SimInfo::getNGlobalConstraints() {
159 >  int SimInfo::getNGlobalConstraints() {
160      int nGlobalConstraints;
161   #ifdef IS_MPI
162      MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
# Line 163 | Line 165 | int SimInfo::getNGlobalConstraints() {
165      nGlobalConstraints =  nConstraints_;
166   #endif
167      return nGlobalConstraints;
168 < }
168 >  }
169  
170 < bool SimInfo::addMolecule(Molecule* mol) {
170 >  bool SimInfo::addMolecule(Molecule* mol) {
171      MoleculeIterator i;
172  
173      i = molecules_.find(mol->getGlobalIndex());
174      if (i == molecules_.end() ) {
175  
176 <        molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
176 >      molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
177          
178 <        nAtoms_ += mol->getNAtoms();
179 <        nBonds_ += mol->getNBonds();
180 <        nBends_ += mol->getNBends();
181 <        nTorsions_ += mol->getNTorsions();
182 <        nRigidBodies_ += mol->getNRigidBodies();
183 <        nIntegrableObjects_ += mol->getNIntegrableObjects();
184 <        nCutoffGroups_ += mol->getNCutoffGroups();
185 <        nConstraints_ += mol->getNConstraintPairs();
178 >      nAtoms_ += mol->getNAtoms();
179 >      nBonds_ += mol->getNBonds();
180 >      nBends_ += mol->getNBends();
181 >      nTorsions_ += mol->getNTorsions();
182 >      nRigidBodies_ += mol->getNRigidBodies();
183 >      nIntegrableObjects_ += mol->getNIntegrableObjects();
184 >      nCutoffGroups_ += mol->getNCutoffGroups();
185 >      nConstraints_ += mol->getNConstraintPairs();
186  
187 <        addExcludePairs(mol);
187 >      addExcludePairs(mol);
188          
189 <        return true;
189 >      return true;
190      } else {
191 <        return false;
191 >      return false;
192      }
193 < }
193 >  }
194  
195 < bool SimInfo::removeMolecule(Molecule* mol) {
195 >  bool SimInfo::removeMolecule(Molecule* mol) {
196      MoleculeIterator i;
197      i = molecules_.find(mol->getGlobalIndex());
198  
199      if (i != molecules_.end() ) {
200  
201 <        assert(mol == i->second);
201 >      assert(mol == i->second);
202          
203 <        nAtoms_ -= mol->getNAtoms();
204 <        nBonds_ -= mol->getNBonds();
205 <        nBends_ -= mol->getNBends();
206 <        nTorsions_ -= mol->getNTorsions();
207 <        nRigidBodies_ -= mol->getNRigidBodies();
208 <        nIntegrableObjects_ -= mol->getNIntegrableObjects();
209 <        nCutoffGroups_ -= mol->getNCutoffGroups();
210 <        nConstraints_ -= mol->getNConstraintPairs();
211 <
212 <        removeExcludePairs(mol);
213 <        molecules_.erase(mol->getGlobalIndex());
214 <
215 <        delete mol;
203 >      nAtoms_ -= mol->getNAtoms();
204 >      nBonds_ -= mol->getNBonds();
205 >      nBends_ -= mol->getNBends();
206 >      nTorsions_ -= mol->getNTorsions();
207 >      nRigidBodies_ -= mol->getNRigidBodies();
208 >      nIntegrableObjects_ -= mol->getNIntegrableObjects();
209 >      nCutoffGroups_ -= mol->getNCutoffGroups();
210 >      nConstraints_ -= mol->getNConstraintPairs();
211 >
212 >      removeExcludePairs(mol);
213 >      molecules_.erase(mol->getGlobalIndex());
214 >
215 >      delete mol;
216          
217 <        return true;
217 >      return true;
218      } else {
219 <        return false;
219 >      return false;
220      }
221  
222  
223 < }    
223 >  }    
224  
225          
226 < Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
226 >  Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
227      i = molecules_.begin();
228      return i == molecules_.end() ? NULL : i->second;
229 < }    
229 >  }    
230  
231 < Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
231 >  Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
232      ++i;
233      return i == molecules_.end() ? NULL : i->second;    
234 < }
234 >  }
235  
236  
237 < void SimInfo::calcNdf() {
237 >  void SimInfo::calcNdf() {
238      int ndf_local;
239      MoleculeIterator i;
240      std::vector<StuntDouble*>::iterator j;
# Line 242 | Line 244 | void SimInfo::calcNdf() {
244      ndf_local = 0;
245      
246      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
247 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
248 <               integrableObject = mol->nextIntegrableObject(j)) {
247 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
248 >           integrableObject = mol->nextIntegrableObject(j)) {
249  
250 <            ndf_local += 3;
250 >        ndf_local += 3;
251  
252 <            if (integrableObject->isDirectional()) {
253 <                if (integrableObject->isLinear()) {
254 <                    ndf_local += 2;
255 <                } else {
256 <                    ndf_local += 3;
257 <                }
258 <            }
252 >        if (integrableObject->isDirectional()) {
253 >          if (integrableObject->isLinear()) {
254 >            ndf_local += 2;
255 >          } else {
256 >            ndf_local += 3;
257 >          }
258 >        }
259              
260 <        }//end for (integrableObject)
260 >      }//end for (integrableObject)
261      }// end for (mol)
262      
263      // n_constraints is local, so subtract them on each processor
# Line 271 | Line 273 | void SimInfo::calcNdf() {
273      // entire system:
274      ndf_ = ndf_ - 3 - nZconstraint_;
275  
276 < }
276 >  }
277  
278 < void SimInfo::calcNdfRaw() {
278 >  void SimInfo::calcNdfRaw() {
279      int ndfRaw_local;
280  
281      MoleculeIterator i;
# Line 285 | Line 287 | void SimInfo::calcNdfRaw() {
287      ndfRaw_local = 0;
288      
289      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
290 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
291 <               integrableObject = mol->nextIntegrableObject(j)) {
290 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
291 >           integrableObject = mol->nextIntegrableObject(j)) {
292  
293 <            ndfRaw_local += 3;
293 >        ndfRaw_local += 3;
294  
295 <            if (integrableObject->isDirectional()) {
296 <                if (integrableObject->isLinear()) {
297 <                    ndfRaw_local += 2;
298 <                } else {
299 <                    ndfRaw_local += 3;
300 <                }
301 <            }
295 >        if (integrableObject->isDirectional()) {
296 >          if (integrableObject->isLinear()) {
297 >            ndfRaw_local += 2;
298 >          } else {
299 >            ndfRaw_local += 3;
300 >          }
301 >        }
302              
303 <        }
303 >      }
304      }
305      
306   #ifdef IS_MPI
# Line 306 | Line 308 | void SimInfo::calcNdfRaw() {
308   #else
309      ndfRaw_ = ndfRaw_local;
310   #endif
311 < }
311 >  }
312  
313 < void SimInfo::calcNdfTrans() {
313 >  void SimInfo::calcNdfTrans() {
314      int ndfTrans_local;
315  
316      ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
# Line 322 | Line 324 | void SimInfo::calcNdfTrans() {
324  
325      ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
326  
327 < }
327 >  }
328  
329 < void SimInfo::addExcludePairs(Molecule* mol) {
329 >  void SimInfo::addExcludePairs(Molecule* mol) {
330      std::vector<Bond*>::iterator bondIter;
331      std::vector<Bend*>::iterator bendIter;
332      std::vector<Torsion*>::iterator torsionIter;
# Line 337 | Line 339 | void SimInfo::addExcludePairs(Molecule* mol) {
339      int d;
340      
341      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
342 <        a = bond->getAtomA()->getGlobalIndex();
343 <        b = bond->getAtomB()->getGlobalIndex();        
344 <        exclude_.addPair(a, b);
342 >      a = bond->getAtomA()->getGlobalIndex();
343 >      b = bond->getAtomB()->getGlobalIndex();        
344 >      exclude_.addPair(a, b);
345      }
346  
347      for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
348 <        a = bend->getAtomA()->getGlobalIndex();
349 <        b = bend->getAtomB()->getGlobalIndex();        
350 <        c = bend->getAtomC()->getGlobalIndex();
348 >      a = bend->getAtomA()->getGlobalIndex();
349 >      b = bend->getAtomB()->getGlobalIndex();        
350 >      c = bend->getAtomC()->getGlobalIndex();
351  
352 <        exclude_.addPair(a, b);
353 <        exclude_.addPair(a, c);
354 <        exclude_.addPair(b, c);        
352 >      exclude_.addPair(a, b);
353 >      exclude_.addPair(a, c);
354 >      exclude_.addPair(b, c);        
355      }
356  
357      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
358 <        a = torsion->getAtomA()->getGlobalIndex();
359 <        b = torsion->getAtomB()->getGlobalIndex();        
360 <        c = torsion->getAtomC()->getGlobalIndex();        
361 <        d = torsion->getAtomD()->getGlobalIndex();        
358 >      a = torsion->getAtomA()->getGlobalIndex();
359 >      b = torsion->getAtomB()->getGlobalIndex();        
360 >      c = torsion->getAtomC()->getGlobalIndex();        
361 >      d = torsion->getAtomD()->getGlobalIndex();        
362  
363 <        exclude_.addPair(a, b);
364 <        exclude_.addPair(a, c);
365 <        exclude_.addPair(a, d);
366 <        exclude_.addPair(b, c);
367 <        exclude_.addPair(b, d);
368 <        exclude_.addPair(c, d);        
363 >      exclude_.addPair(a, b);
364 >      exclude_.addPair(a, c);
365 >      exclude_.addPair(a, d);
366 >      exclude_.addPair(b, c);
367 >      exclude_.addPair(b, d);
368 >      exclude_.addPair(c, d);        
369      }
370  
371      Molecule::RigidBodyIterator rbIter;
372      RigidBody* rb;
373      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
374 <        std::vector<Atom*> atoms = rb->getAtoms();
375 <        for (int i = 0; i < atoms.size() -1 ; ++i) {
376 <            for (int j = i + 1; j < atoms.size(); ++j) {
377 <                a = atoms[i]->getGlobalIndex();
378 <                b = atoms[j]->getGlobalIndex();
379 <                exclude_.addPair(a, b);
380 <            }
381 <        }
374 >      std::vector<Atom*> atoms = rb->getAtoms();
375 >      for (int i = 0; i < atoms.size() -1 ; ++i) {
376 >        for (int j = i + 1; j < atoms.size(); ++j) {
377 >          a = atoms[i]->getGlobalIndex();
378 >          b = atoms[j]->getGlobalIndex();
379 >          exclude_.addPair(a, b);
380 >        }
381 >      }
382      }        
383  
384 < }
384 >  }
385  
386 < void SimInfo::removeExcludePairs(Molecule* mol) {
386 >  void SimInfo::removeExcludePairs(Molecule* mol) {
387      std::vector<Bond*>::iterator bondIter;
388      std::vector<Bend*>::iterator bendIter;
389      std::vector<Torsion*>::iterator torsionIter;
# Line 394 | Line 396 | void SimInfo::removeExcludePairs(Molecule* mol) {
396      int d;
397      
398      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
399 <        a = bond->getAtomA()->getGlobalIndex();
400 <        b = bond->getAtomB()->getGlobalIndex();        
401 <        exclude_.removePair(a, b);
399 >      a = bond->getAtomA()->getGlobalIndex();
400 >      b = bond->getAtomB()->getGlobalIndex();        
401 >      exclude_.removePair(a, b);
402      }
403  
404      for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
405 <        a = bend->getAtomA()->getGlobalIndex();
406 <        b = bend->getAtomB()->getGlobalIndex();        
407 <        c = bend->getAtomC()->getGlobalIndex();
405 >      a = bend->getAtomA()->getGlobalIndex();
406 >      b = bend->getAtomB()->getGlobalIndex();        
407 >      c = bend->getAtomC()->getGlobalIndex();
408  
409 <        exclude_.removePair(a, b);
410 <        exclude_.removePair(a, c);
411 <        exclude_.removePair(b, c);        
409 >      exclude_.removePair(a, b);
410 >      exclude_.removePair(a, c);
411 >      exclude_.removePair(b, c);        
412      }
413  
414      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
415 <        a = torsion->getAtomA()->getGlobalIndex();
416 <        b = torsion->getAtomB()->getGlobalIndex();        
417 <        c = torsion->getAtomC()->getGlobalIndex();        
418 <        d = torsion->getAtomD()->getGlobalIndex();        
415 >      a = torsion->getAtomA()->getGlobalIndex();
416 >      b = torsion->getAtomB()->getGlobalIndex();        
417 >      c = torsion->getAtomC()->getGlobalIndex();        
418 >      d = torsion->getAtomD()->getGlobalIndex();        
419  
420 <        exclude_.removePair(a, b);
421 <        exclude_.removePair(a, c);
422 <        exclude_.removePair(a, d);
423 <        exclude_.removePair(b, c);
424 <        exclude_.removePair(b, d);
425 <        exclude_.removePair(c, d);        
420 >      exclude_.removePair(a, b);
421 >      exclude_.removePair(a, c);
422 >      exclude_.removePair(a, d);
423 >      exclude_.removePair(b, c);
424 >      exclude_.removePair(b, d);
425 >      exclude_.removePair(c, d);        
426      }
427  
428      Molecule::RigidBodyIterator rbIter;
429      RigidBody* rb;
430      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
431 <        std::vector<Atom*> atoms = rb->getAtoms();
432 <        for (int i = 0; i < atoms.size() -1 ; ++i) {
433 <            for (int j = i + 1; j < atoms.size(); ++j) {
434 <                a = atoms[i]->getGlobalIndex();
435 <                b = atoms[j]->getGlobalIndex();
436 <                exclude_.removePair(a, b);
437 <            }
438 <        }
431 >      std::vector<Atom*> atoms = rb->getAtoms();
432 >      for (int i = 0; i < atoms.size() -1 ; ++i) {
433 >        for (int j = i + 1; j < atoms.size(); ++j) {
434 >          a = atoms[i]->getGlobalIndex();
435 >          b = atoms[j]->getGlobalIndex();
436 >          exclude_.removePair(a, b);
437 >        }
438 >      }
439      }        
440  
441 < }
441 >  }
442  
443  
444 < void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
444 >  void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
445      int curStampId;
446  
447      //index from 0
# Line 447 | Line 449 | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp
449  
450      moleculeStamps_.push_back(molStamp);
451      molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
452 < }
452 >  }
453  
454 < void SimInfo::update() {
454 >  void SimInfo::update() {
455  
456      setupSimType();
457  
# Line 462 | Line 464 | void SimInfo::update() {
464      //setup fortran force field
465      /** @deprecate */    
466      int isError = 0;
467 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
467 >    
468 >    setupCoulombicCorrection( isError );
469 >
470      if(isError){
471 <        sprintf( painCave.errMsg,
472 <         "ForceField error: There was an error initializing the forceField in fortran.\n" );
473 <        painCave.isFatal = 1;
474 <        simError();
471 >      sprintf( painCave.errMsg,
472 >               "ForceField error: There was an error initializing the forceField in fortran.\n" );
473 >      painCave.isFatal = 1;
474 >      simError();
475      }
476    
477      
# Line 478 | Line 482 | void SimInfo::update() {
482      calcNdfTrans();
483  
484      fortranInitialized_ = true;
485 < }
485 >  }
486  
487 < std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
487 >  std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
488      SimInfo::MoleculeIterator mi;
489      Molecule* mol;
490      Molecule::AtomIterator ai;
# Line 489 | Line 493 | std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
493  
494      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
495  
496 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
497 <            atomTypes.insert(atom->getAtomType());
498 <        }
496 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
497 >        atomTypes.insert(atom->getAtomType());
498 >      }
499          
500      }
501  
502      return atomTypes;        
503 < }
503 >  }
504  
505 < void SimInfo::setupSimType() {
505 >  void SimInfo::setupSimType() {
506      std::set<AtomType*>::iterator i;
507      std::set<AtomType*> atomTypes;
508      atomTypes = getUniqueAtomTypes();
# Line 511 | Line 515 | void SimInfo::setupSimType() {
515      int useDipole = 0;
516      int useGayBerne = 0;
517      int useSticky = 0;
518 +    int useStickyPower = 0;
519      int useShape = 0;
520      int useFLARB = 0; //it is not in AtomType yet
521      int useDirectionalAtom = 0;    
# Line 521 | Line 526 | void SimInfo::setupSimType() {
526  
527      //loop over all of the atom types
528      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
529 <        useLennardJones |= (*i)->isLennardJones();
530 <        useElectrostatic |= (*i)->isElectrostatic();
531 <        useEAM |= (*i)->isEAM();
532 <        useCharge |= (*i)->isCharge();
533 <        useDirectional |= (*i)->isDirectional();
534 <        useDipole |= (*i)->isDipole();
535 <        useGayBerne |= (*i)->isGayBerne();
536 <        useSticky |= (*i)->isSticky();
537 <        useShape |= (*i)->isShape();
529 >      useLennardJones |= (*i)->isLennardJones();
530 >      useElectrostatic |= (*i)->isElectrostatic();
531 >      useEAM |= (*i)->isEAM();
532 >      useCharge |= (*i)->isCharge();
533 >      useDirectional |= (*i)->isDirectional();
534 >      useDipole |= (*i)->isDipole();
535 >      useGayBerne |= (*i)->isGayBerne();
536 >      useSticky |= (*i)->isSticky();
537 >      useStickyPower |= (*i)->isStickyPower();
538 >      useShape |= (*i)->isShape();
539      }
540  
541 <    if (useSticky || useDipole || useGayBerne || useShape) {
542 <        useDirectionalAtom = 1;
541 >    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
542 >      useDirectionalAtom = 1;
543      }
544  
545      if (useCharge || useDipole) {
546 <        useElectrostatics = 1;
546 >      useElectrostatics = 1;
547      }
548  
549   #ifdef IS_MPI    
# Line 564 | Line 570 | void SimInfo::setupSimType() {
570      temp = useSticky;
571      MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
572  
573 +    temp = useStickyPower;
574 +    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
575 +    
576      temp = useGayBerne;
577      MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
578  
# Line 578 | Line 587 | void SimInfo::setupSimType() {
587  
588      temp = useRF;
589      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
590 +
591 +    temp = useUW;
592 +    MPI_Allreduce(&temp, &useUW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
593 +
594 +    temp = useDW;
595 +    MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
596      
597   #endif
598  
# Line 588 | Line 603 | void SimInfo::setupSimType() {
603      fInfo_.SIM_uses_Charges = useCharge;
604      fInfo_.SIM_uses_Dipoles = useDipole;
605      fInfo_.SIM_uses_Sticky = useSticky;
606 +    fInfo_.SIM_uses_StickyPower = useStickyPower;
607      fInfo_.SIM_uses_GayBerne = useGayBerne;
608      fInfo_.SIM_uses_EAM = useEAM;
609      fInfo_.SIM_uses_Shapes = useShape;
# Line 596 | Line 612 | void SimInfo::setupSimType() {
612  
613      if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
614  
615 <        if (simParams_->haveDielectric()) {
616 <            fInfo_.dielect = simParams_->getDielectric();
617 <        } else {
618 <            sprintf(painCave.errMsg,
619 <                    "SimSetup Error: No Dielectric constant was set.\n"
620 <                    "\tYou are trying to use Reaction Field without"
621 <                    "\tsetting a dielectric constant!\n");
622 <            painCave.isFatal = 1;
623 <            simError();
624 <        }
615 >      if (simParams_->haveDielectric()) {
616 >        fInfo_.dielect = simParams_->getDielectric();
617 >      } else {
618 >        sprintf(painCave.errMsg,
619 >                "SimSetup Error: No Dielectric constant was set.\n"
620 >                "\tYou are trying to use Reaction Field without"
621 >                "\tsetting a dielectric constant!\n");
622 >        painCave.isFatal = 1;
623 >        simError();
624 >      }
625          
626      } else {
627 <        fInfo_.dielect = 0.0;
627 >      fInfo_.dielect = 0.0;
628      }
629  
630 < }
630 >  }
631  
632 < void SimInfo::setupFortranSim() {
632 >  void SimInfo::setupFortranSim() {
633      int isError;
634      int nExclude;
635      std::vector<int> fortranGlobalGroupMembership;
# Line 623 | Line 639 | void SimInfo::setupFortranSim() {
639  
640      //globalGroupMembership_ is filled by SimCreator    
641      for (int i = 0; i < nGlobalAtoms_; i++) {
642 <        fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
642 >      fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
643      }
644  
645      //calculate mass ratio of cutoff group
# Line 640 | Line 656 | void SimInfo::setupFortranSim() {
656      mfact.reserve(getNCutoffGroups());
657      
658      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
659 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
659 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
660  
661 <            totalMass = cg->getMass();
662 <            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
663 <                        mfact.push_back(atom->getMass()/totalMass);
664 <            }
661 >        totalMass = cg->getMass();
662 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
663 >          mfact.push_back(atom->getMass()/totalMass);
664 >        }
665  
666 <        }      
666 >      }      
667      }
668  
669      //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
# Line 657 | Line 673 | void SimInfo::setupFortranSim() {
673      identArray.reserve(getNAtoms());
674      
675      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
676 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
677 <            identArray.push_back(atom->getIdent());
678 <        }
676 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
677 >        identArray.push_back(atom->getIdent());
678 >      }
679      }    
680  
681      //fill molMembershipArray
682      //molMembershipArray is filled by SimCreator    
683      std::vector<int> molMembershipArray(nGlobalAtoms_);
684      for (int i = 0; i < nGlobalAtoms_; i++) {
685 <        molMembershipArray[i] = globalMolMembership_[i] + 1;
685 >      molMembershipArray[i] = globalMolMembership_[i] + 1;
686      }
687      
688      //setup fortran simulation
# Line 674 | Line 690 | void SimInfo::setupFortranSim() {
690      int* globalExcludes = NULL;
691      int* excludeList = exclude_.getExcludeList();
692      setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
693 <                  &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
694 <                  &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
693 >                   &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
694 >                   &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
695  
696      if( isError ){
697  
698 <        sprintf( painCave.errMsg,
699 <                 "There was an error setting the simulation information in fortran.\n" );
700 <        painCave.isFatal = 1;
701 <        painCave.severity = OOPSE_ERROR;
702 <        simError();
698 >      sprintf( painCave.errMsg,
699 >               "There was an error setting the simulation information in fortran.\n" );
700 >      painCave.isFatal = 1;
701 >      painCave.severity = OOPSE_ERROR;
702 >      simError();
703      }
704  
705   #ifdef IS_MPI
706      sprintf( checkPointMsg,
707 <       "succesfully sent the simulation information to fortran.\n");
707 >             "succesfully sent the simulation information to fortran.\n");
708      MPIcheckPoint();
709   #endif // is_mpi
710 < }
710 >  }
711  
712  
713   #ifdef IS_MPI
714 < void SimInfo::setupFortranParallel() {
714 >  void SimInfo::setupFortranParallel() {
715      
716      //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
717      std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
# Line 711 | Line 727 | void SimInfo::setupFortranParallel() {
727  
728      for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
729  
730 <        //local index(index in DataStorge) of atom is important
731 <        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
732 <            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
733 <        }
730 >      //local index(index in DataStorge) of atom is important
731 >      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
732 >        localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
733 >      }
734  
735 <        //local index of cutoff group is trivial, it only depends on the order of travesing
736 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
737 <            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
738 <        }        
735 >      //local index of cutoff group is trivial, it only depends on the order of travesing
736 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
737 >        localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
738 >      }        
739          
740      }
741  
# Line 739 | Line 755 | void SimInfo::setupFortranParallel() {
755                      &localToGlobalCutoffGroupIndex[0], &isError);
756  
757      if (isError) {
758 <        sprintf(painCave.errMsg,
759 <                "mpiRefresh errror: fortran didn't like something we gave it.\n");
760 <        painCave.isFatal = 1;
761 <        simError();
758 >      sprintf(painCave.errMsg,
759 >              "mpiRefresh errror: fortran didn't like something we gave it.\n");
760 >      painCave.isFatal = 1;
761 >      simError();
762      }
763  
764      sprintf(checkPointMsg, " mpiRefresh successful.\n");
765      MPIcheckPoint();
766  
767  
768 < }
768 >  }
769  
770   #endif
771  
772 < double SimInfo::calcMaxCutoffRadius() {
772 >  double SimInfo::calcMaxCutoffRadius() {
773  
774  
775      std::set<AtomType*> atomTypes;
# Line 765 | Line 781 | double SimInfo::calcMaxCutoffRadius() {
781  
782      //query the max cutoff radius among these atom types
783      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
784 <        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
784 >      cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
785      }
786  
787      double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
# Line 774 | Line 790 | double SimInfo::calcMaxCutoffRadius() {
790   #endif
791  
792      return maxCutoffRadius;
793 < }
793 >  }
794  
795 < void SimInfo::getCutoff(double& rcut, double& rsw) {
795 >  void SimInfo::getCutoff(double& rcut, double& rsw) {
796      
797      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
798          
799 <        if (!simParams_->haveRcut()){
800 <            sprintf(painCave.errMsg,
799 >      if (!simParams_->haveRcut()){
800 >        sprintf(painCave.errMsg,
801                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
802                  "\tOOPSE will use a default value of 15.0 angstroms"
803                  "\tfor the cutoffRadius.\n");
804 <            painCave.isFatal = 0;
805 <            simError();
806 <            rcut = 15.0;
807 <        } else{
808 <            rcut = simParams_->getRcut();
809 <        }
804 >        painCave.isFatal = 0;
805 >        simError();
806 >        rcut = 15.0;
807 >      } else{
808 >        rcut = simParams_->getRcut();
809 >      }
810  
811 <        if (!simParams_->haveRsw()){
812 <            sprintf(painCave.errMsg,
811 >      if (!simParams_->haveRsw()){
812 >        sprintf(painCave.errMsg,
813                  "SimCreator Warning: No value was set for switchingRadius.\n"
814                  "\tOOPSE will use a default value of\n"
815                  "\t0.95 * cutoffRadius for the switchingRadius\n");
816 <            painCave.isFatal = 0;
817 <            simError();
818 <            rsw = 0.95 * rcut;
819 <        } else{
820 <            rsw = simParams_->getRsw();
821 <        }
816 >        painCave.isFatal = 0;
817 >        simError();
818 >        rsw = 0.95 * rcut;
819 >      } else{
820 >        rsw = simParams_->getRsw();
821 >      }
822  
823      } else {
824 <        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
825 <        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
824 >      // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
825 >      //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
826          
827 <        if (simParams_->haveRcut()) {
828 <            rcut = simParams_->getRcut();
829 <        } else {
830 <            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
831 <            rcut = calcMaxCutoffRadius();
832 <        }
827 >      if (simParams_->haveRcut()) {
828 >        rcut = simParams_->getRcut();
829 >      } else {
830 >        //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
831 >        rcut = calcMaxCutoffRadius();
832 >      }
833  
834 <        if (simParams_->haveRsw()) {
835 <            rsw  = simParams_->getRsw();
836 <        } else {
837 <            rsw = rcut;
838 <        }
834 >      if (simParams_->haveRsw()) {
835 >        rsw  = simParams_->getRsw();
836 >      } else {
837 >        rsw = rcut;
838 >      }
839      
840      }
841 < }
841 >  }
842  
843 < void SimInfo::setupCutoff() {
843 >  void SimInfo::setupCutoff() {    
844      getCutoff(rcut_, rsw_);    
845      double rnblist = rcut_ + 1; // skin of neighbor list
846  
847      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
848 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
849 < }
848 >    
849 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
850 >    if (simParams_->haveCutoffPolicy()) {
851 >      std::string myPolicy = simParams_->getCutoffPolicy();
852 >      if (myPolicy == "MIX") {
853 >        cp = MIX_CUTOFF_POLICY;
854 >      } else {
855 >        if (myPolicy == "MAX") {
856 >          cp = MAX_CUTOFF_POLICY;
857 >        } else {
858 >          if (myPolicy == "TRADITIONAL") {            
859 >            cp = TRADITIONAL_CUTOFF_POLICY;
860 >          } else {
861 >            // throw error        
862 >            sprintf( painCave.errMsg,
863 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
864 >            painCave.isFatal = 1;
865 >            simError();
866 >          }    
867 >        }          
868 >      }
869 >    }
870 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
871 >  }
872  
873 < void SimInfo::addProperty(GenericData* genData) {
873 >  void SimInfo::setupCoulombicCorrection( int isError ) {    
874 >    
875 >    int errorOut;
876 >    int cc =  NONE;
877 >    double alphaVal;
878 >
879 >    errorOut = isError;
880 >
881 >    if (simParams_->haveCoulombicCorrection()) {
882 >      std::string myCorrection = simParams_->getCoulombicCorrection();
883 >      if (myCorrection == "NONE") {
884 >        cc = NONE;
885 >      } else {
886 >        if (myCorrection == "UNDAMPED_WOLF") {
887 >          cc = UNDAMPED_WOLF;
888 >        } else {
889 >          if (myCorrection == "WOLF") {            
890 >            cc = WOLF;
891 >            if (!simParams_->haveDampingAlpha()) {
892 >              //throw error
893 >              sprintf( painCave.errMsg,
894 >                       "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Wolf Coulombic Correction.", simParams_->getDampingAlpha());
895 >              painCave.isFatal = 0;
896 >              simError();
897 >            }
898 >            alphaVal = simParams_->getDampingAlpha();
899 >          } else {
900 >            if (myCorrection == "REACTION_FIELD") {
901 >              cc = REACTION_FIELD;
902 >            } else {
903 >              // throw error        
904 >              sprintf( painCave.errMsg,
905 >                       "SimInfo error: Unknown coulombicCorrection. (Input file specified %s .)\n\tcoulombicCorrection must be one of: \"none\", \"undamped_wolf\", \"wolf\", or \"reaction_field\".", myCorrection.c_str() );
906 >              painCave.isFatal = 1;
907 >              simError();
908 >            }    
909 >          }          
910 >        }
911 >      }
912 >    }
913 >    initFortranFF( &fInfo_.SIM_uses_RF, &cc, &alphaVal, &errorOut );
914 >  }
915 >
916 >  void SimInfo::addProperty(GenericData* genData) {
917      properties_.addProperty(genData);  
918 < }
918 >  }
919  
920 < void SimInfo::removeProperty(const std::string& propName) {
920 >  void SimInfo::removeProperty(const std::string& propName) {
921      properties_.removeProperty(propName);  
922 < }
922 >  }
923  
924 < void SimInfo::clearProperties() {
924 >  void SimInfo::clearProperties() {
925      properties_.clearProperties();
926 < }
926 >  }
927  
928 < std::vector<std::string> SimInfo::getPropertyNames() {
928 >  std::vector<std::string> SimInfo::getPropertyNames() {
929      return properties_.getPropertyNames();  
930 < }
930 >  }
931        
932 < std::vector<GenericData*> SimInfo::getProperties() {
932 >  std::vector<GenericData*> SimInfo::getProperties() {
933      return properties_.getProperties();
934 < }
934 >  }
935  
936 < GenericData* SimInfo::getPropertyByName(const std::string& propName) {
936 >  GenericData* SimInfo::getPropertyByName(const std::string& propName) {
937      return properties_.getPropertyByName(propName);
938 < }
938 >  }
939  
940 < void SimInfo::setSnapshotManager(SnapshotManager* sman) {
940 >  void SimInfo::setSnapshotManager(SnapshotManager* sman) {
941      if (sman_ == sman) {
942 <        return;
942 >      return;
943      }    
944      delete sman_;
945      sman_ = sman;
# Line 872 | Line 953 | void SimInfo::setSnapshotManager(SnapshotManager* sman
953  
954      for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
955          
956 <        for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
957 <            atom->setSnapshotManager(sman_);
958 <        }
956 >      for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
957 >        atom->setSnapshotManager(sman_);
958 >      }
959          
960 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
961 <            rb->setSnapshotManager(sman_);
962 <        }
960 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
961 >        rb->setSnapshotManager(sman_);
962 >      }
963      }    
964      
965 < }
965 >  }
966  
967 < Vector3d SimInfo::getComVel(){
967 >  Vector3d SimInfo::getComVel(){
968      SimInfo::MoleculeIterator i;
969      Molecule* mol;
970  
# Line 892 | Line 973 | Vector3d SimInfo::getComVel(){
973      
974  
975      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
976 <        double mass = mol->getMass();
977 <        totalMass += mass;
978 <        comVel += mass * mol->getComVel();
976 >      double mass = mol->getMass();
977 >      totalMass += mass;
978 >      comVel += mass * mol->getComVel();
979      }  
980  
981   #ifdef IS_MPI
# Line 907 | Line 988 | Vector3d SimInfo::getComVel(){
988      comVel /= totalMass;
989  
990      return comVel;
991 < }
991 >  }
992  
993 < Vector3d SimInfo::getCom(){
993 >  Vector3d SimInfo::getCom(){
994      SimInfo::MoleculeIterator i;
995      Molecule* mol;
996  
# Line 917 | Line 998 | Vector3d SimInfo::getCom(){
998      double totalMass = 0.0;
999      
1000      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1001 <        double mass = mol->getMass();
1002 <        totalMass += mass;
1003 <        com += mass * mol->getCom();
1001 >      double mass = mol->getMass();
1002 >      totalMass += mass;
1003 >      com += mass * mol->getCom();
1004      }  
1005  
1006   #ifdef IS_MPI
# Line 933 | Line 1014 | Vector3d SimInfo::getCom(){
1014  
1015      return com;
1016  
1017 < }        
1017 >  }        
1018  
1019 < std::ostream& operator <<(std::ostream& o, SimInfo& info) {
1019 >  std::ostream& operator <<(std::ostream& o, SimInfo& info) {
1020  
1021      return o;
1022 < }
1022 >  }
1023 >  
1024 >  
1025 >   /*
1026 >   Returns center of mass and center of mass velocity in one function call.
1027 >   */
1028 >  
1029 >   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1030 >      SimInfo::MoleculeIterator i;
1031 >      Molecule* mol;
1032 >      
1033 >    
1034 >      double totalMass = 0.0;
1035 >    
1036  
1037 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1038 +         double mass = mol->getMass();
1039 +         totalMass += mass;
1040 +         com += mass * mol->getCom();
1041 +         comVel += mass * mol->getComVel();          
1042 +      }  
1043 +      
1044 + #ifdef IS_MPI
1045 +      double tmpMass = totalMass;
1046 +      Vector3d tmpCom(com);  
1047 +      Vector3d tmpComVel(comVel);
1048 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1049 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1050 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1051 + #endif
1052 +      
1053 +      com /= totalMass;
1054 +      comVel /= totalMass;
1055 +   }        
1056 +  
1057 +   /*
1058 +   Return intertia tensor for entire system and angular momentum Vector.
1059 +
1060 +
1061 +       [  Ixx -Ixy  -Ixz ]
1062 +  J =| -Iyx  Iyy  -Iyz |
1063 +       [ -Izx -Iyz   Izz ]
1064 +    */
1065 +
1066 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1067 +      
1068 +
1069 +      double xx = 0.0;
1070 +      double yy = 0.0;
1071 +      double zz = 0.0;
1072 +      double xy = 0.0;
1073 +      double xz = 0.0;
1074 +      double yz = 0.0;
1075 +      Vector3d com(0.0);
1076 +      Vector3d comVel(0.0);
1077 +      
1078 +      getComAll(com, comVel);
1079 +      
1080 +      SimInfo::MoleculeIterator i;
1081 +      Molecule* mol;
1082 +      
1083 +      Vector3d thisq(0.0);
1084 +      Vector3d thisv(0.0);
1085 +
1086 +      double thisMass = 0.0;
1087 +    
1088 +      
1089 +      
1090 +  
1091 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1092 +        
1093 +         thisq = mol->getCom()-com;
1094 +         thisv = mol->getComVel()-comVel;
1095 +         thisMass = mol->getMass();
1096 +         // Compute moment of intertia coefficients.
1097 +         xx += thisq[0]*thisq[0]*thisMass;
1098 +         yy += thisq[1]*thisq[1]*thisMass;
1099 +         zz += thisq[2]*thisq[2]*thisMass;
1100 +        
1101 +         // compute products of intertia
1102 +         xy += thisq[0]*thisq[1]*thisMass;
1103 +         xz += thisq[0]*thisq[2]*thisMass;
1104 +         yz += thisq[1]*thisq[2]*thisMass;
1105 +            
1106 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1107 +            
1108 +      }  
1109 +      
1110 +      
1111 +      inertiaTensor(0,0) = yy + zz;
1112 +      inertiaTensor(0,1) = -xy;
1113 +      inertiaTensor(0,2) = -xz;
1114 +      inertiaTensor(1,0) = -xy;
1115 +      inertiaTensor(1,1) = xx + zz;
1116 +      inertiaTensor(1,2) = -yz;
1117 +      inertiaTensor(2,0) = -xz;
1118 +      inertiaTensor(2,1) = -yz;
1119 +      inertiaTensor(2,2) = xx + yy;
1120 +      
1121 + #ifdef IS_MPI
1122 +      Mat3x3d tmpI(inertiaTensor);
1123 +      Vector3d tmpAngMom;
1124 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1125 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1126 + #endif
1127 +              
1128 +      return;
1129 +   }
1130 +
1131 +   //Returns the angular momentum of the system
1132 +   Vector3d SimInfo::getAngularMomentum(){
1133 +      
1134 +      Vector3d com(0.0);
1135 +      Vector3d comVel(0.0);
1136 +      Vector3d angularMomentum(0.0);
1137 +      
1138 +      getComAll(com,comVel);
1139 +      
1140 +      SimInfo::MoleculeIterator i;
1141 +      Molecule* mol;
1142 +      
1143 +      Vector3d thisr(0.0);
1144 +      Vector3d thisp(0.0);
1145 +      
1146 +      double thisMass;
1147 +      
1148 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1149 +        thisMass = mol->getMass();
1150 +        thisr = mol->getCom()-com;
1151 +        thisp = (mol->getComVel()-comVel)*thisMass;
1152 +        
1153 +        angularMomentum += cross( thisr, thisp );
1154 +        
1155 +      }  
1156 +      
1157 + #ifdef IS_MPI
1158 +      Vector3d tmpAngMom;
1159 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1160 + #endif
1161 +      
1162 +      return angularMomentum;
1163 +   }
1164 +  
1165 +  
1166   }//end namespace oopse
1167  

Diff Legend

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