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 316 by tim, Fri Feb 11 22:41:02 2005 UTC vs.
Revision 734 by chuckv, Tue Nov 15 16:05:38 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/DarkSide/fElectrostaticSummationMethod.h"
57 + #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
58 + #include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
59   #include "UseTheForce/doForces_interface.h"
60 + #include "UseTheForce/DarkSide/electrostatic_interface.h"
61   #include "UseTheForce/notifyCutoffs_interface.h"
62 + #include "UseTheForce/DarkSide/switcheroo_interface.h"
63   #include "utils/MemoryUtils.hpp"
64   #include "utils/simError.h"
65   #include "selection/SelectionManager.hpp"
# Line 65 | Line 71 | namespace oopse {
71  
72   namespace oopse {
73  
74 < SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
75 <                                ForceField* ff, Globals* simParams) :
76 <                                forceField_(ff), simParams_(simParams),
77 <                                ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
78 <                                nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
79 <                                nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
80 <                                nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
81 <                                nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
82 <                                sman_(NULL), fortranInitialized_(false), selectMan_(NULL) {
74 >  SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
75 >                   ForceField* ff, Globals* simParams) :
76 >    stamps_(stamps), forceField_(ff), simParams_(simParams),
77 >    ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
78 >    nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
79 >    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
80 >    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
81 >    nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
82 >    sman_(NULL), fortranInitialized_(false) {
83  
84              
85 <    std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
86 <    MoleculeStamp* molStamp;
87 <    int nMolWithSameStamp;
88 <    int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
89 <    int nGroups = 0;          //total cutoff groups defined in meta-data file
90 <    CutoffGroupStamp* cgStamp;    
91 <    RigidBodyStamp* rbStamp;
92 <    int nRigidAtoms = 0;
85 >      std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
86 >      MoleculeStamp* molStamp;
87 >      int nMolWithSameStamp;
88 >      int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
89 >      int nGroups = 0;      //total cutoff groups defined in meta-data file
90 >      CutoffGroupStamp* cgStamp;    
91 >      RigidBodyStamp* rbStamp;
92 >      int nRigidAtoms = 0;
93      
94 <    for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
94 >      for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
95          molStamp = i->first;
96          nMolWithSameStamp = i->second;
97          
# Line 100 | Line 106 | SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*,
106          int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
107          
108          for (int j=0; j < nCutoffGroupsInStamp; j++) {
109 <            cgStamp = molStamp->getCutoffGroup(j);
110 <            nAtomsInGroups += cgStamp->getNMembers();
109 >          cgStamp = molStamp->getCutoffGroup(j);
110 >          nAtomsInGroups += cgStamp->getNMembers();
111          }
112  
113          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
114 +
115          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
116  
117          //calculate atoms in rigid bodies
# Line 112 | Line 119 | SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*,
119          int nRigidBodiesInStamp = molStamp->getNRigidBodies();
120          
121          for (int j=0; j < nRigidBodiesInStamp; j++) {
122 <            rbStamp = molStamp->getRigidBody(j);
123 <            nAtomsInRigidBodies += rbStamp->getNMembers();
122 >          rbStamp = molStamp->getRigidBody(j);
123 >          nAtomsInRigidBodies += rbStamp->getNMembers();
124          }
125  
126          nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
127          nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
128          
129 <    }
129 >      }
130  
131 <    //every free atom (atom does not belong to cutoff groups) is a cutoff group
132 <    //therefore the total number of cutoff groups in the system is equal to
133 <    //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
134 <    //file plus the number of cutoff groups defined in meta-data file
135 <    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
131 >      //every free atom (atom does not belong to cutoff groups) is a cutoff
132 >      //group therefore the total number of cutoff groups in the system is
133 >      //equal to the total number of atoms minus number of atoms belong to
134 >      //cutoff group defined in meta-data file plus the number of cutoff
135 >      //groups defined in meta-data file
136 >      nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
137  
138 <    //every free atom (atom does not belong to rigid bodies) is an integrable object
139 <    //therefore the total number of  integrable objects in the system is equal to
140 <    //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
141 <    //file plus the number of  rigid bodies defined in meta-data file
142 <    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
138 >      //every free atom (atom does not belong to rigid bodies) is an
139 >      //integrable object therefore the total number of integrable objects
140 >      //in the system is equal to the total number of atoms minus number of
141 >      //atoms belong to rigid body defined in meta-data file plus the number
142 >      //of rigid bodies defined in meta-data file
143 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
144 >                                                + nGlobalRigidBodies_;
145 >  
146 >      nGlobalMols_ = molStampIds_.size();
147  
136 –    nGlobalMols_ = molStampIds_.size();
137 –
148   #ifdef IS_MPI    
149 <    molToProcMap_.resize(nGlobalMols_);
149 >      molToProcMap_.resize(nGlobalMols_);
150   #endif
151  
152 <    selectMan_ = new SelectionManager(this);
143 <    selectMan_->selectAll();
144 < }
152 >    }
153  
154 < SimInfo::~SimInfo() {
155 <    //MemoryUtils::deleteVectorOfPointer(molecules_);
156 <
157 <    MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
158 <    
154 >  SimInfo::~SimInfo() {
155 >    std::map<int, Molecule*>::iterator i;
156 >    for (i = molecules_.begin(); i != molecules_.end(); ++i) {
157 >      delete i->second;
158 >    }
159 >    molecules_.clear();
160 >      
161 >    delete stamps_;
162      delete sman_;
163      delete simParams_;
164      delete forceField_;
165 <    delete selectMan_;
155 < }
165 >  }
166  
167 < int SimInfo::getNGlobalConstraints() {
167 >  int SimInfo::getNGlobalConstraints() {
168      int nGlobalConstraints;
169   #ifdef IS_MPI
170      MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
# Line 163 | Line 173 | int SimInfo::getNGlobalConstraints() {
173      nGlobalConstraints =  nConstraints_;
174   #endif
175      return nGlobalConstraints;
176 < }
176 >  }
177  
178 < bool SimInfo::addMolecule(Molecule* mol) {
178 >  bool SimInfo::addMolecule(Molecule* mol) {
179      MoleculeIterator i;
180  
181      i = molecules_.find(mol->getGlobalIndex());
182      if (i == molecules_.end() ) {
183  
184 <        molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
184 >      molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
185          
186 <        nAtoms_ += mol->getNAtoms();
187 <        nBonds_ += mol->getNBonds();
188 <        nBends_ += mol->getNBends();
189 <        nTorsions_ += mol->getNTorsions();
190 <        nRigidBodies_ += mol->getNRigidBodies();
191 <        nIntegrableObjects_ += mol->getNIntegrableObjects();
192 <        nCutoffGroups_ += mol->getNCutoffGroups();
193 <        nConstraints_ += mol->getNConstraintPairs();
186 >      nAtoms_ += mol->getNAtoms();
187 >      nBonds_ += mol->getNBonds();
188 >      nBends_ += mol->getNBends();
189 >      nTorsions_ += mol->getNTorsions();
190 >      nRigidBodies_ += mol->getNRigidBodies();
191 >      nIntegrableObjects_ += mol->getNIntegrableObjects();
192 >      nCutoffGroups_ += mol->getNCutoffGroups();
193 >      nConstraints_ += mol->getNConstraintPairs();
194  
195 <        addExcludePairs(mol);
195 >      addExcludePairs(mol);
196          
197 <        return true;
197 >      return true;
198      } else {
199 <        return false;
199 >      return false;
200      }
201 < }
201 >  }
202  
203 < bool SimInfo::removeMolecule(Molecule* mol) {
203 >  bool SimInfo::removeMolecule(Molecule* mol) {
204      MoleculeIterator i;
205      i = molecules_.find(mol->getGlobalIndex());
206  
207      if (i != molecules_.end() ) {
208  
209 <        assert(mol == i->second);
209 >      assert(mol == i->second);
210          
211 <        nAtoms_ -= mol->getNAtoms();
212 <        nBonds_ -= mol->getNBonds();
213 <        nBends_ -= mol->getNBends();
214 <        nTorsions_ -= mol->getNTorsions();
215 <        nRigidBodies_ -= mol->getNRigidBodies();
216 <        nIntegrableObjects_ -= mol->getNIntegrableObjects();
217 <        nCutoffGroups_ -= mol->getNCutoffGroups();
218 <        nConstraints_ -= mol->getNConstraintPairs();
209 <
210 <        removeExcludePairs(mol);
211 <        molecules_.erase(mol->getGlobalIndex());
211 >      nAtoms_ -= mol->getNAtoms();
212 >      nBonds_ -= mol->getNBonds();
213 >      nBends_ -= mol->getNBends();
214 >      nTorsions_ -= mol->getNTorsions();
215 >      nRigidBodies_ -= mol->getNRigidBodies();
216 >      nIntegrableObjects_ -= mol->getNIntegrableObjects();
217 >      nCutoffGroups_ -= mol->getNCutoffGroups();
218 >      nConstraints_ -= mol->getNConstraintPairs();
219  
220 <        delete mol;
220 >      removeExcludePairs(mol);
221 >      molecules_.erase(mol->getGlobalIndex());
222 >
223 >      delete mol;
224          
225 <        return true;
225 >      return true;
226      } else {
227 <        return false;
227 >      return false;
228      }
229  
230  
231 < }    
231 >  }    
232  
233          
234 < Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
234 >  Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
235      i = molecules_.begin();
236      return i == molecules_.end() ? NULL : i->second;
237 < }    
237 >  }    
238  
239 < Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
239 >  Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
240      ++i;
241      return i == molecules_.end() ? NULL : i->second;    
242 < }
242 >  }
243  
244  
245 < void SimInfo::calcNdf() {
245 >  void SimInfo::calcNdf() {
246      int ndf_local;
247      MoleculeIterator i;
248      std::vector<StuntDouble*>::iterator j;
# Line 242 | Line 252 | void SimInfo::calcNdf() {
252      ndf_local = 0;
253      
254      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
255 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
256 <               integrableObject = mol->nextIntegrableObject(j)) {
255 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
256 >           integrableObject = mol->nextIntegrableObject(j)) {
257  
258 <            ndf_local += 3;
258 >        ndf_local += 3;
259  
260 <            if (integrableObject->isDirectional()) {
261 <                if (integrableObject->isLinear()) {
262 <                    ndf_local += 2;
263 <                } else {
264 <                    ndf_local += 3;
265 <                }
266 <            }
260 >        if (integrableObject->isDirectional()) {
261 >          if (integrableObject->isLinear()) {
262 >            ndf_local += 2;
263 >          } else {
264 >            ndf_local += 3;
265 >          }
266 >        }
267              
268 <        }//end for (integrableObject)
268 >      }//end for (integrableObject)
269      }// end for (mol)
270      
271      // n_constraints is local, so subtract them on each processor
# Line 271 | Line 281 | void SimInfo::calcNdf() {
281      // entire system:
282      ndf_ = ndf_ - 3 - nZconstraint_;
283  
284 < }
284 >  }
285  
286 < void SimInfo::calcNdfRaw() {
286 >  void SimInfo::calcNdfRaw() {
287      int ndfRaw_local;
288  
289      MoleculeIterator i;
# Line 285 | Line 295 | void SimInfo::calcNdfRaw() {
295      ndfRaw_local = 0;
296      
297      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
298 <        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
299 <               integrableObject = mol->nextIntegrableObject(j)) {
298 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
299 >           integrableObject = mol->nextIntegrableObject(j)) {
300  
301 <            ndfRaw_local += 3;
301 >        ndfRaw_local += 3;
302  
303 <            if (integrableObject->isDirectional()) {
304 <                if (integrableObject->isLinear()) {
305 <                    ndfRaw_local += 2;
306 <                } else {
307 <                    ndfRaw_local += 3;
308 <                }
309 <            }
303 >        if (integrableObject->isDirectional()) {
304 >          if (integrableObject->isLinear()) {
305 >            ndfRaw_local += 2;
306 >          } else {
307 >            ndfRaw_local += 3;
308 >          }
309 >        }
310              
311 <        }
311 >      }
312      }
313      
314   #ifdef IS_MPI
# Line 306 | Line 316 | void SimInfo::calcNdfRaw() {
316   #else
317      ndfRaw_ = ndfRaw_local;
318   #endif
319 < }
319 >  }
320  
321 < void SimInfo::calcNdfTrans() {
321 >  void SimInfo::calcNdfTrans() {
322      int ndfTrans_local;
323  
324      ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
# Line 322 | Line 332 | void SimInfo::calcNdfTrans() {
332  
333      ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
334  
335 < }
335 >  }
336  
337 < void SimInfo::addExcludePairs(Molecule* mol) {
337 >  void SimInfo::addExcludePairs(Molecule* mol) {
338      std::vector<Bond*>::iterator bondIter;
339      std::vector<Bend*>::iterator bendIter;
340      std::vector<Torsion*>::iterator torsionIter;
# Line 337 | Line 347 | void SimInfo::addExcludePairs(Molecule* mol) {
347      int d;
348      
349      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
350 <        a = bond->getAtomA()->getGlobalIndex();
351 <        b = bond->getAtomB()->getGlobalIndex();        
352 <        exclude_.addPair(a, b);
350 >      a = bond->getAtomA()->getGlobalIndex();
351 >      b = bond->getAtomB()->getGlobalIndex();        
352 >      exclude_.addPair(a, b);
353      }
354  
355      for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
356 <        a = bend->getAtomA()->getGlobalIndex();
357 <        b = bend->getAtomB()->getGlobalIndex();        
358 <        c = bend->getAtomC()->getGlobalIndex();
356 >      a = bend->getAtomA()->getGlobalIndex();
357 >      b = bend->getAtomB()->getGlobalIndex();        
358 >      c = bend->getAtomC()->getGlobalIndex();
359  
360 <        exclude_.addPair(a, b);
361 <        exclude_.addPair(a, c);
362 <        exclude_.addPair(b, c);        
360 >      exclude_.addPair(a, b);
361 >      exclude_.addPair(a, c);
362 >      exclude_.addPair(b, c);        
363      }
364  
365      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
366 <        a = torsion->getAtomA()->getGlobalIndex();
367 <        b = torsion->getAtomB()->getGlobalIndex();        
368 <        c = torsion->getAtomC()->getGlobalIndex();        
369 <        d = torsion->getAtomD()->getGlobalIndex();        
366 >      a = torsion->getAtomA()->getGlobalIndex();
367 >      b = torsion->getAtomB()->getGlobalIndex();        
368 >      c = torsion->getAtomC()->getGlobalIndex();        
369 >      d = torsion->getAtomD()->getGlobalIndex();        
370  
371 <        exclude_.addPair(a, b);
372 <        exclude_.addPair(a, c);
373 <        exclude_.addPair(a, d);
374 <        exclude_.addPair(b, c);
375 <        exclude_.addPair(b, d);
376 <        exclude_.addPair(c, d);        
371 >      exclude_.addPair(a, b);
372 >      exclude_.addPair(a, c);
373 >      exclude_.addPair(a, d);
374 >      exclude_.addPair(b, c);
375 >      exclude_.addPair(b, d);
376 >      exclude_.addPair(c, d);        
377      }
378  
379 <    
380 < }
379 >    Molecule::RigidBodyIterator rbIter;
380 >    RigidBody* rb;
381 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
382 >      std::vector<Atom*> atoms = rb->getAtoms();
383 >      for (int i = 0; i < atoms.size() -1 ; ++i) {
384 >        for (int j = i + 1; j < atoms.size(); ++j) {
385 >          a = atoms[i]->getGlobalIndex();
386 >          b = atoms[j]->getGlobalIndex();
387 >          exclude_.addPair(a, b);
388 >        }
389 >      }
390 >    }        
391  
392 < void SimInfo::removeExcludePairs(Molecule* mol) {
392 >  }
393 >
394 >  void SimInfo::removeExcludePairs(Molecule* mol) {
395      std::vector<Bond*>::iterator bondIter;
396      std::vector<Bend*>::iterator bendIter;
397      std::vector<Torsion*>::iterator torsionIter;
# Line 382 | Line 404 | void SimInfo::removeExcludePairs(Molecule* mol) {
404      int d;
405      
406      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
407 <        a = bond->getAtomA()->getGlobalIndex();
408 <        b = bond->getAtomB()->getGlobalIndex();        
409 <        exclude_.removePair(a, b);
407 >      a = bond->getAtomA()->getGlobalIndex();
408 >      b = bond->getAtomB()->getGlobalIndex();        
409 >      exclude_.removePair(a, b);
410      }
411  
412      for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
413 <        a = bend->getAtomA()->getGlobalIndex();
414 <        b = bend->getAtomB()->getGlobalIndex();        
415 <        c = bend->getAtomC()->getGlobalIndex();
413 >      a = bend->getAtomA()->getGlobalIndex();
414 >      b = bend->getAtomB()->getGlobalIndex();        
415 >      c = bend->getAtomC()->getGlobalIndex();
416  
417 <        exclude_.removePair(a, b);
418 <        exclude_.removePair(a, c);
419 <        exclude_.removePair(b, c);        
417 >      exclude_.removePair(a, b);
418 >      exclude_.removePair(a, c);
419 >      exclude_.removePair(b, c);        
420      }
421  
422      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
423 <        a = torsion->getAtomA()->getGlobalIndex();
424 <        b = torsion->getAtomB()->getGlobalIndex();        
425 <        c = torsion->getAtomC()->getGlobalIndex();        
426 <        d = torsion->getAtomD()->getGlobalIndex();        
423 >      a = torsion->getAtomA()->getGlobalIndex();
424 >      b = torsion->getAtomB()->getGlobalIndex();        
425 >      c = torsion->getAtomC()->getGlobalIndex();        
426 >      d = torsion->getAtomD()->getGlobalIndex();        
427  
428 <        exclude_.removePair(a, b);
429 <        exclude_.removePair(a, c);
430 <        exclude_.removePair(a, d);
431 <        exclude_.removePair(b, c);
432 <        exclude_.removePair(b, d);
433 <        exclude_.removePair(c, d);        
428 >      exclude_.removePair(a, b);
429 >      exclude_.removePair(a, c);
430 >      exclude_.removePair(a, d);
431 >      exclude_.removePair(b, c);
432 >      exclude_.removePair(b, d);
433 >      exclude_.removePair(c, d);        
434      }
435  
436 < }
436 >    Molecule::RigidBodyIterator rbIter;
437 >    RigidBody* rb;
438 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
439 >      std::vector<Atom*> atoms = rb->getAtoms();
440 >      for (int i = 0; i < atoms.size() -1 ; ++i) {
441 >        for (int j = i + 1; j < atoms.size(); ++j) {
442 >          a = atoms[i]->getGlobalIndex();
443 >          b = atoms[j]->getGlobalIndex();
444 >          exclude_.removePair(a, b);
445 >        }
446 >      }
447 >    }        
448  
449 +  }
450  
451 < void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
451 >
452 >  void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
453      int curStampId;
454  
455      //index from 0
# Line 422 | Line 457 | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp
457  
458      moleculeStamps_.push_back(molStamp);
459      molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
460 < }
460 >  }
461  
462 < void SimInfo::update() {
462 >  void SimInfo::update() {
463  
464      setupSimType();
465  
# Line 437 | Line 472 | void SimInfo::update() {
472      //setup fortran force field
473      /** @deprecate */    
474      int isError = 0;
475 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
475 >    
476 >    setupElectrostaticSummationMethod( isError );
477 >    setupSwitchingFunction();
478 >
479      if(isError){
480 <        sprintf( painCave.errMsg,
481 <         "ForceField error: There was an error initializing the forceField in fortran.\n" );
482 <        painCave.isFatal = 1;
483 <        simError();
480 >      sprintf( painCave.errMsg,
481 >               "ForceField error: There was an error initializing the forceField in fortran.\n" );
482 >      painCave.isFatal = 1;
483 >      simError();
484      }
485    
486      
# Line 453 | Line 491 | void SimInfo::update() {
491      calcNdfTrans();
492  
493      fortranInitialized_ = true;
494 < }
494 >  }
495  
496 < std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
496 >  std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
497      SimInfo::MoleculeIterator mi;
498      Molecule* mol;
499      Molecule::AtomIterator ai;
# Line 464 | Line 502 | std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
502  
503      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
504  
505 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
506 <            atomTypes.insert(atom->getAtomType());
507 <        }
505 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
506 >        atomTypes.insert(atom->getAtomType());
507 >      }
508          
509      }
510  
511      return atomTypes;        
512 < }
512 >  }
513  
514 < void SimInfo::setupSimType() {
514 >  void SimInfo::setupSimType() {
515      std::set<AtomType*>::iterator i;
516      std::set<AtomType*> atomTypes;
517      atomTypes = getUniqueAtomTypes();
# Line 481 | Line 519 | void SimInfo::setupSimType() {
519      int useLennardJones = 0;
520      int useElectrostatic = 0;
521      int useEAM = 0;
522 +    int useSC = 0;
523      int useCharge = 0;
524      int useDirectional = 0;
525      int useDipole = 0;
526      int useGayBerne = 0;
527      int useSticky = 0;
528 +    int useStickyPower = 0;
529      int useShape = 0;
530      int useFLARB = 0; //it is not in AtomType yet
531      int useDirectionalAtom = 0;    
532      int useElectrostatics = 0;
533      //usePBC and useRF are from simParams
534 <    int usePBC = simParams_->getPBC();
535 <    int useRF = simParams_->getUseRF();
534 >    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
535 >    int useRF;
536 >    int useSF;
537 >    std::string myMethod;
538  
539 +    // set the useRF logical
540 +    useRF = 0;
541 +    useSF = 0;
542 +
543 +
544 +    if (simParams_->haveElectrostaticSummationMethod()) {
545 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
546 +      toUpper(myMethod);
547 +      if (myMethod == "REACTION_FIELD") {
548 +        useRF=1;
549 +      } else {
550 +        if (myMethod == "SHIFTED_FORCE") {
551 +          useSF = 1;
552 +        }
553 +      }
554 +    }
555 +
556      //loop over all of the atom types
557      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
558 <        useLennardJones |= (*i)->isLennardJones();
559 <        useElectrostatic |= (*i)->isElectrostatic();
560 <        useEAM |= (*i)->isEAM();
561 <        useCharge |= (*i)->isCharge();
562 <        useDirectional |= (*i)->isDirectional();
563 <        useDipole |= (*i)->isDipole();
564 <        useGayBerne |= (*i)->isGayBerne();
565 <        useSticky |= (*i)->isSticky();
566 <        useShape |= (*i)->isShape();
558 >      useLennardJones |= (*i)->isLennardJones();
559 >      useElectrostatic |= (*i)->isElectrostatic();
560 >      useEAM |= (*i)->isEAM();
561 >      useSC |= (*i)->isSC();
562 >      useCharge |= (*i)->isCharge();
563 >      useDirectional |= (*i)->isDirectional();
564 >      useDipole |= (*i)->isDipole();
565 >      useGayBerne |= (*i)->isGayBerne();
566 >      useSticky |= (*i)->isSticky();
567 >      useStickyPower |= (*i)->isStickyPower();
568 >      useShape |= (*i)->isShape();
569      }
570  
571 <    if (useSticky || useDipole || useGayBerne || useShape) {
572 <        useDirectionalAtom = 1;
571 >    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
572 >      useDirectionalAtom = 1;
573      }
574  
575      if (useCharge || useDipole) {
576 <        useElectrostatics = 1;
576 >      useElectrostatics = 1;
577      }
578  
579   #ifdef IS_MPI    
# Line 539 | Line 600 | void SimInfo::setupSimType() {
600      temp = useSticky;
601      MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
602  
603 +    temp = useStickyPower;
604 +    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
605 +    
606      temp = useGayBerne;
607      MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
608  
609      temp = useEAM;
610      MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
611  
612 +    temp = useSC;
613 +    MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
614 +    
615      temp = useShape;
616      MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
617  
# Line 553 | Line 620 | void SimInfo::setupSimType() {
620  
621      temp = useRF;
622      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
623 <    
623 >
624 >    temp = useSF;
625 >    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
626 >
627   #endif
628  
629      fInfo_.SIM_uses_PBC = usePBC;    
# Line 563 | Line 633 | void SimInfo::setupSimType() {
633      fInfo_.SIM_uses_Charges = useCharge;
634      fInfo_.SIM_uses_Dipoles = useDipole;
635      fInfo_.SIM_uses_Sticky = useSticky;
636 +    fInfo_.SIM_uses_StickyPower = useStickyPower;
637      fInfo_.SIM_uses_GayBerne = useGayBerne;
638      fInfo_.SIM_uses_EAM = useEAM;
639 +    fInfo_.SIM_uses_SC = useSC;
640      fInfo_.SIM_uses_Shapes = useShape;
641      fInfo_.SIM_uses_FLARB = useFLARB;
642      fInfo_.SIM_uses_RF = useRF;
643 +    fInfo_.SIM_uses_SF = useSF;
644  
645 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
646 <
647 <        if (simParams_->haveDielectric()) {
648 <            fInfo_.dielect = simParams_->getDielectric();
649 <        } else {
650 <            sprintf(painCave.errMsg,
651 <                    "SimSetup Error: No Dielectric constant was set.\n"
652 <                    "\tYou are trying to use Reaction Field without"
653 <                    "\tsetting a dielectric constant!\n");
654 <            painCave.isFatal = 1;
655 <            simError();
656 <        }
584 <        
585 <    } else {
586 <        fInfo_.dielect = 0.0;
645 >    if( myMethod == "REACTION_FIELD") {
646 >      
647 >      if (simParams_->haveDielectric()) {
648 >        fInfo_.dielect = simParams_->getDielectric();
649 >      } else {
650 >        sprintf(painCave.errMsg,
651 >                "SimSetup Error: No Dielectric constant was set.\n"
652 >                "\tYou are trying to use Reaction Field without"
653 >                "\tsetting a dielectric constant!\n");
654 >        painCave.isFatal = 1;
655 >        simError();
656 >      }      
657      }
658  
659 < }
659 >  }
660  
661 < void SimInfo::setupFortranSim() {
661 >  void SimInfo::setupFortranSim() {
662      int isError;
663      int nExclude;
664      std::vector<int> fortranGlobalGroupMembership;
# Line 598 | Line 668 | void SimInfo::setupFortranSim() {
668  
669      //globalGroupMembership_ is filled by SimCreator    
670      for (int i = 0; i < nGlobalAtoms_; i++) {
671 <        fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
671 >      fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
672      }
673  
674      //calculate mass ratio of cutoff group
# Line 615 | Line 685 | void SimInfo::setupFortranSim() {
685      mfact.reserve(getNCutoffGroups());
686      
687      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
688 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
688 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
689  
690 <            totalMass = cg->getMass();
691 <            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
692 <                        mfact.push_back(atom->getMass()/totalMass);
693 <            }
690 >        totalMass = cg->getMass();
691 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
692 >          // Check for massless groups - set mfact to 1 if true
693 >          if (totalMass != 0)
694 >            mfact.push_back(atom->getMass()/totalMass);
695 >          else
696 >            mfact.push_back( 1.0 );
697 >        }
698  
699 <        }      
699 >      }      
700      }
701  
702      //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
# Line 632 | Line 706 | void SimInfo::setupFortranSim() {
706      identArray.reserve(getNAtoms());
707      
708      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
709 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
710 <            identArray.push_back(atom->getIdent());
711 <        }
709 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
710 >        identArray.push_back(atom->getIdent());
711 >      }
712      }    
713  
714      //fill molMembershipArray
715      //molMembershipArray is filled by SimCreator    
716      std::vector<int> molMembershipArray(nGlobalAtoms_);
717      for (int i = 0; i < nGlobalAtoms_; i++) {
718 <        molMembershipArray[i] = globalMolMembership_[i] + 1;
718 >      molMembershipArray[i] = globalMolMembership_[i] + 1;
719      }
720      
721      //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
722      int nGlobalExcludes = 0;
723      int* globalExcludes = NULL;
724      int* excludeList = exclude_.getExcludeList();
725      setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
726 <                  &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
727 <                  &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
726 >                   &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
727 >                   &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
728  
729      if( isError ){
730  
731 <        sprintf( painCave.errMsg,
732 <                 "There was an error setting the simulation information in fortran.\n" );
733 <        painCave.isFatal = 1;
734 <        painCave.severity = OOPSE_ERROR;
735 <        simError();
731 >      sprintf( painCave.errMsg,
732 >               "There was an error setting the simulation information in fortran.\n" );
733 >      painCave.isFatal = 1;
734 >      painCave.severity = OOPSE_ERROR;
735 >      simError();
736      }
737  
738   #ifdef IS_MPI
739      sprintf( checkPointMsg,
740 <       "succesfully sent the simulation information to fortran.\n");
740 >             "succesfully sent the simulation information to fortran.\n");
741      MPIcheckPoint();
742   #endif // is_mpi
743 < }
743 >  }
744  
745  
746   #ifdef IS_MPI
747 < void SimInfo::setupFortranParallel() {
747 >  void SimInfo::setupFortranParallel() {
748      
749      //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
750      std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
# Line 689 | Line 760 | void SimInfo::setupFortranParallel() {
760  
761      for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
762  
763 <        //local index(index in DataStorge) of atom is important
764 <        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
765 <            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
766 <        }
763 >      //local index(index in DataStorge) of atom is important
764 >      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
765 >        localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
766 >      }
767  
768 <        //local index of cutoff group is trivial, it only depends on the order of travesing
769 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
770 <            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
771 <        }        
768 >      //local index of cutoff group is trivial, it only depends on the order of travesing
769 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
770 >        localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
771 >      }        
772          
773      }
774  
# Line 717 | Line 788 | void SimInfo::setupFortranParallel() {
788                      &localToGlobalCutoffGroupIndex[0], &isError);
789  
790      if (isError) {
791 <        sprintf(painCave.errMsg,
792 <                "mpiRefresh errror: fortran didn't like something we gave it.\n");
793 <        painCave.isFatal = 1;
794 <        simError();
791 >      sprintf(painCave.errMsg,
792 >              "mpiRefresh errror: fortran didn't like something we gave it.\n");
793 >      painCave.isFatal = 1;
794 >      simError();
795      }
796  
797      sprintf(checkPointMsg, " mpiRefresh successful.\n");
798      MPIcheckPoint();
799  
800  
801 < }
801 >  }
802  
803   #endif
804  
805 < double SimInfo::calcMaxCutoffRadius() {
805 >  double SimInfo::calcMaxCutoffRadius() {
806  
807  
808      std::set<AtomType*> atomTypes;
# Line 743 | Line 814 | double SimInfo::calcMaxCutoffRadius() {
814  
815      //query the max cutoff radius among these atom types
816      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
817 <        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
817 >      cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
818      }
819  
820      double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
# Line 752 | Line 823 | double SimInfo::calcMaxCutoffRadius() {
823   #endif
824  
825      return maxCutoffRadius;
826 < }
826 >  }
827  
828 < void SimInfo::setupCutoff() {
758 <    double rcut_;  //cutoff radius
759 <    double rsw_; //switching radius
828 >  void SimInfo::getCutoff(double& rcut, double& rsw) {
829      
830      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
831          
832 <        if (!simParams_->haveRcut()){
833 <            sprintf(painCave.errMsg,
832 >      if (!simParams_->haveCutoffRadius()){
833 >        sprintf(painCave.errMsg,
834                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
835                  "\tOOPSE will use a default value of 15.0 angstroms"
836                  "\tfor the cutoffRadius.\n");
837 <            painCave.isFatal = 0;
838 <            simError();
839 <            rcut_ = 15.0;
840 <        } else{
841 <            rcut_ = simParams_->getRcut();
842 <        }
837 >        painCave.isFatal = 0;
838 >        simError();
839 >        rcut = 15.0;
840 >      } else{
841 >        rcut = simParams_->getCutoffRadius();
842 >      }
843  
844 <        if (!simParams_->haveRsw()){
845 <            sprintf(painCave.errMsg,
844 >      if (!simParams_->haveSwitchingRadius()){
845 >        sprintf(painCave.errMsg,
846                  "SimCreator Warning: No value was set for switchingRadius.\n"
847                  "\tOOPSE will use a default value of\n"
848 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
849 <            painCave.isFatal = 0;
850 <            simError();
851 <            rsw_ = 0.95 * rcut_;
852 <        } else{
853 <            rsw_ = simParams_->getRsw();
854 <        }
848 >                "\t0.85 * cutoffRadius for the switchingRadius\n");
849 >        painCave.isFatal = 0;
850 >        simError();
851 >        rsw = 0.85 * rcut;
852 >      } else{
853 >        rsw = simParams_->getSwitchingRadius();
854 >      }
855  
856      } else {
857 <        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
858 <        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
857 >      // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
858 >      //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
859          
860 <        if (simParams_->haveRcut()) {
861 <            rcut_ = simParams_->getRcut();
862 <        } else {
863 <            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
864 <            rcut_ = calcMaxCutoffRadius();
865 <        }
860 >      if (simParams_->haveCutoffRadius()) {
861 >        rcut = simParams_->getCutoffRadius();
862 >      } else {
863 >        //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
864 >        rcut = calcMaxCutoffRadius();
865 >      }
866  
867 <        if (simParams_->haveRsw()) {
868 <            rsw_  = simParams_->getRsw();
869 <        } else {
870 <            rsw_ = rcut_;
871 <        }
867 >      if (simParams_->haveSwitchingRadius()) {
868 >        rsw  = simParams_->getSwitchingRadius();
869 >      } else {
870 >        rsw = rcut;
871 >      }
872      
873      }
874 <        
874 >  }
875 >
876 >  void SimInfo::setupCutoff() {    
877 >    getCutoff(rcut_, rsw_);    
878      double rnblist = rcut_ + 1; // skin of neighbor list
879  
880      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
881 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
882 < }
881 >    
882 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
883 >    if (simParams_->haveCutoffPolicy()) {
884 >      std::string myPolicy = simParams_->getCutoffPolicy();
885 >      toUpper(myPolicy);
886 >      if (myPolicy == "MIX") {
887 >        cp = MIX_CUTOFF_POLICY;
888 >      } else {
889 >        if (myPolicy == "MAX") {
890 >          cp = MAX_CUTOFF_POLICY;
891 >        } else {
892 >          if (myPolicy == "TRADITIONAL") {            
893 >            cp = TRADITIONAL_CUTOFF_POLICY;
894 >          } else {
895 >            // throw error        
896 >            sprintf( painCave.errMsg,
897 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
898 >            painCave.isFatal = 1;
899 >            simError();
900 >          }    
901 >        }          
902 >      }
903 >    }
904  
812 – void SimInfo::addProperty(GenericData* genData) {
813 –    properties_.addProperty(genData);  
814 – }
905  
906 < void SimInfo::removeProperty(const std::string& propName) {
907 <    properties_.removeProperty(propName);  
908 < }
906 >    if (simParams_->haveSkinThickness()) {
907 >      double skinThickness = simParams_->getSkinThickness();
908 >    }
909  
910 < void SimInfo::clearProperties() {
911 <    properties_.clearProperties();
912 < }
910 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
911 >    // also send cutoff notification to electrostatics
912 >    setElectrostaticCutoffRadius(&rcut_, &rsw_);
913 >  }
914  
915 < std::vector<std::string> SimInfo::getPropertyNames() {
916 <    return properties_.getPropertyNames();  
917 < }
918 <      
919 < std::vector<GenericData*> SimInfo::getProperties() {
920 <    return properties_.getProperties();
921 < }
915 >  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
916 >    
917 >    int errorOut;
918 >    int esm =  NONE;
919 >    int sm = UNDAMPED;
920 >    double alphaVal;
921 >    double dielectric;
922  
923 < GenericData* SimInfo::getPropertyByName(const std::string& propName) {
924 <    return properties_.getPropertyByName(propName);
925 < }
923 >    errorOut = isError;
924 >    alphaVal = simParams_->getDampingAlpha();
925 >    dielectric = simParams_->getDielectric();
926  
927 < void SimInfo::setSnapshotManager(SnapshotManager* sman) {
927 >    if (simParams_->haveElectrostaticSummationMethod()) {
928 >      std::string myMethod = simParams_->getElectrostaticSummationMethod();
929 >      toUpper(myMethod);
930 >      if (myMethod == "NONE") {
931 >        esm = NONE;
932 >      } else {
933 >        if (myMethod == "SWITCHING_FUNCTION") {
934 >          esm = SWITCHING_FUNCTION;
935 >        } else {
936 >          if (myMethod == "SHIFTED_POTENTIAL") {
937 >            esm = SHIFTED_POTENTIAL;
938 >          } else {
939 >            if (myMethod == "SHIFTED_FORCE") {            
940 >              esm = SHIFTED_FORCE;
941 >            } else {
942 >              if (myMethod == "REACTION_FIELD") {            
943 >                esm = REACTION_FIELD;
944 >              } else {
945 >                // throw error        
946 >                sprintf( painCave.errMsg,
947 >                         "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() );
948 >                painCave.isFatal = 1;
949 >                simError();
950 >              }    
951 >            }          
952 >          }
953 >        }
954 >      }
955 >    }
956 >    
957 >    if (simParams_->haveElectrostaticScreeningMethod()) {
958 >      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
959 >      toUpper(myScreen);
960 >      if (myScreen == "UNDAMPED") {
961 >        sm = UNDAMPED;
962 >      } else {
963 >        if (myScreen == "DAMPED") {
964 >          sm = DAMPED;
965 >          if (!simParams_->haveDampingAlpha()) {
966 >            //throw error
967 >            sprintf( painCave.errMsg,
968 >                     "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal);
969 >            painCave.isFatal = 0;
970 >            simError();
971 >          }
972 >        } else {
973 >          // throw error        
974 >          sprintf( painCave.errMsg,
975 >                   "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() );
976 >          painCave.isFatal = 1;
977 >          simError();
978 >        }
979 >      }
980 >    }
981 >    
982 >    // let's pass some summation method variables to fortran
983 >    setElectrostaticSummationMethod( &esm );
984 >    setScreeningMethod( &sm );
985 >    setDampingAlpha( &alphaVal );
986 >    setReactionFieldDielectric( &dielectric );
987 >    initFortranFF( &esm, &errorOut );
988 >  }
989 >
990 >  void SimInfo::setupSwitchingFunction() {    
991 >    int ft = CUBIC;
992 >
993 >    if (simParams_->haveSwitchingFunctionType()) {
994 >      std::string funcType = simParams_->getSwitchingFunctionType();
995 >      toUpper(funcType);
996 >      if (funcType == "CUBIC") {
997 >        ft = CUBIC;
998 >      } else {
999 >        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
1000 >          ft = FIFTH_ORDER_POLY;
1001 >        } else {
1002 >          // throw error        
1003 >          sprintf( painCave.errMsg,
1004 >                   "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
1005 >          painCave.isFatal = 1;
1006 >          simError();
1007 >        }          
1008 >      }
1009 >    }
1010 >
1011 >    // send switching function notification to switcheroo
1012 >    setFunctionType(&ft);
1013 >
1014 >  }
1015 >
1016 >  void SimInfo::addProperty(GenericData* genData) {
1017 >    properties_.addProperty(genData);  
1018 >  }
1019 >
1020 >  void SimInfo::removeProperty(const std::string& propName) {
1021 >    properties_.removeProperty(propName);  
1022 >  }
1023 >
1024 >  void SimInfo::clearProperties() {
1025 >    properties_.clearProperties();
1026 >  }
1027 >
1028 >  std::vector<std::string> SimInfo::getPropertyNames() {
1029 >    return properties_.getPropertyNames();  
1030 >  }
1031 >      
1032 >  std::vector<GenericData*> SimInfo::getProperties() {
1033 >    return properties_.getProperties();
1034 >  }
1035 >
1036 >  GenericData* SimInfo::getPropertyByName(const std::string& propName) {
1037 >    return properties_.getPropertyByName(propName);
1038 >  }
1039 >
1040 >  void SimInfo::setSnapshotManager(SnapshotManager* sman) {
1041 >    if (sman_ == sman) {
1042 >      return;
1043 >    }    
1044 >    delete sman_;
1045      sman_ = sman;
1046  
1047      Molecule* mol;
# Line 845 | Line 1053 | void SimInfo::setSnapshotManager(SnapshotManager* sman
1053  
1054      for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
1055          
1056 <        for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
1057 <            atom->setSnapshotManager(sman_);
1058 <        }
1056 >      for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
1057 >        atom->setSnapshotManager(sman_);
1058 >      }
1059          
1060 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
1061 <            rb->setSnapshotManager(sman_);
1062 <        }
1060 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
1061 >        rb->setSnapshotManager(sman_);
1062 >      }
1063      }    
1064      
1065 < }
1065 >  }
1066  
1067 < Vector3d SimInfo::getComVel(){
1067 >  Vector3d SimInfo::getComVel(){
1068      SimInfo::MoleculeIterator i;
1069      Molecule* mol;
1070  
# Line 865 | Line 1073 | Vector3d SimInfo::getComVel(){
1073      
1074  
1075      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1076 <        double mass = mol->getMass();
1077 <        totalMass += mass;
1078 <        comVel += mass * mol->getComVel();
1076 >      double mass = mol->getMass();
1077 >      totalMass += mass;
1078 >      comVel += mass * mol->getComVel();
1079      }  
1080  
1081   #ifdef IS_MPI
# Line 880 | Line 1088 | Vector3d SimInfo::getComVel(){
1088      comVel /= totalMass;
1089  
1090      return comVel;
1091 < }
1091 >  }
1092  
1093 < Vector3d SimInfo::getCom(){
1093 >  Vector3d SimInfo::getCom(){
1094      SimInfo::MoleculeIterator i;
1095      Molecule* mol;
1096  
# Line 890 | Line 1098 | Vector3d SimInfo::getCom(){
1098      double totalMass = 0.0;
1099      
1100      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1101 <        double mass = mol->getMass();
1102 <        totalMass += mass;
1103 <        com += mass * mol->getCom();
1101 >      double mass = mol->getMass();
1102 >      totalMass += mass;
1103 >      com += mass * mol->getCom();
1104      }  
1105  
1106   #ifdef IS_MPI
# Line 906 | Line 1114 | Vector3d SimInfo::getCom(){
1114  
1115      return com;
1116  
1117 < }        
1117 >  }        
1118  
1119 < std::ostream& operator <<(std::ostream& o, SimInfo& info) {
1119 >  std::ostream& operator <<(std::ostream& o, SimInfo& info) {
1120  
1121      return o;
1122 < }
1122 >  }
1123 >  
1124 >  
1125 >   /*
1126 >   Returns center of mass and center of mass velocity in one function call.
1127 >   */
1128 >  
1129 >   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1130 >      SimInfo::MoleculeIterator i;
1131 >      Molecule* mol;
1132 >      
1133 >    
1134 >      double totalMass = 0.0;
1135 >    
1136  
1137 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1138 +         double mass = mol->getMass();
1139 +         totalMass += mass;
1140 +         com += mass * mol->getCom();
1141 +         comVel += mass * mol->getComVel();          
1142 +      }  
1143 +      
1144 + #ifdef IS_MPI
1145 +      double tmpMass = totalMass;
1146 +      Vector3d tmpCom(com);  
1147 +      Vector3d tmpComVel(comVel);
1148 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1149 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1150 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1151 + #endif
1152 +      
1153 +      com /= totalMass;
1154 +      comVel /= totalMass;
1155 +   }        
1156 +  
1157 +   /*
1158 +   Return intertia tensor for entire system and angular momentum Vector.
1159 +
1160 +
1161 +       [  Ixx -Ixy  -Ixz ]
1162 +  J =| -Iyx  Iyy  -Iyz |
1163 +       [ -Izx -Iyz   Izz ]
1164 +    */
1165 +
1166 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1167 +      
1168 +
1169 +      double xx = 0.0;
1170 +      double yy = 0.0;
1171 +      double zz = 0.0;
1172 +      double xy = 0.0;
1173 +      double xz = 0.0;
1174 +      double yz = 0.0;
1175 +      Vector3d com(0.0);
1176 +      Vector3d comVel(0.0);
1177 +      
1178 +      getComAll(com, comVel);
1179 +      
1180 +      SimInfo::MoleculeIterator i;
1181 +      Molecule* mol;
1182 +      
1183 +      Vector3d thisq(0.0);
1184 +      Vector3d thisv(0.0);
1185 +
1186 +      double thisMass = 0.0;
1187 +    
1188 +      
1189 +      
1190 +  
1191 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1192 +        
1193 +         thisq = mol->getCom()-com;
1194 +         thisv = mol->getComVel()-comVel;
1195 +         thisMass = mol->getMass();
1196 +         // Compute moment of intertia coefficients.
1197 +         xx += thisq[0]*thisq[0]*thisMass;
1198 +         yy += thisq[1]*thisq[1]*thisMass;
1199 +         zz += thisq[2]*thisq[2]*thisMass;
1200 +        
1201 +         // compute products of intertia
1202 +         xy += thisq[0]*thisq[1]*thisMass;
1203 +         xz += thisq[0]*thisq[2]*thisMass;
1204 +         yz += thisq[1]*thisq[2]*thisMass;
1205 +            
1206 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1207 +            
1208 +      }  
1209 +      
1210 +      
1211 +      inertiaTensor(0,0) = yy + zz;
1212 +      inertiaTensor(0,1) = -xy;
1213 +      inertiaTensor(0,2) = -xz;
1214 +      inertiaTensor(1,0) = -xy;
1215 +      inertiaTensor(1,1) = xx + zz;
1216 +      inertiaTensor(1,2) = -yz;
1217 +      inertiaTensor(2,0) = -xz;
1218 +      inertiaTensor(2,1) = -yz;
1219 +      inertiaTensor(2,2) = xx + yy;
1220 +      
1221 + #ifdef IS_MPI
1222 +      Mat3x3d tmpI(inertiaTensor);
1223 +      Vector3d tmpAngMom;
1224 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1225 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1226 + #endif
1227 +              
1228 +      return;
1229 +   }
1230 +
1231 +   //Returns the angular momentum of the system
1232 +   Vector3d SimInfo::getAngularMomentum(){
1233 +      
1234 +      Vector3d com(0.0);
1235 +      Vector3d comVel(0.0);
1236 +      Vector3d angularMomentum(0.0);
1237 +      
1238 +      getComAll(com,comVel);
1239 +      
1240 +      SimInfo::MoleculeIterator i;
1241 +      Molecule* mol;
1242 +      
1243 +      Vector3d thisr(0.0);
1244 +      Vector3d thisp(0.0);
1245 +      
1246 +      double thisMass;
1247 +      
1248 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1249 +        thisMass = mol->getMass();
1250 +        thisr = mol->getCom()-com;
1251 +        thisp = (mol->getComVel()-comVel)*thisMass;
1252 +        
1253 +        angularMomentum += cross( thisr, thisp );
1254 +        
1255 +      }  
1256 +      
1257 + #ifdef IS_MPI
1258 +      Vector3d tmpAngMom;
1259 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1260 + #endif
1261 +      
1262 +      return angularMomentum;
1263 +   }
1264 +  
1265 +  
1266   }//end namespace oopse
1267  

Diff Legend

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