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

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 1121 by chuckv, Mon Feb 26 04:45:42 2007 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 + * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 + * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   /**
# Line 46 | Line 47
47   * @version 1.0
48   */
49  
50 + #ifdef IS_MPI
51 + #include <mpi.h>
52 + #endif
53   #include <algorithm>
54   #include <set>
55   #include <map>
# Line 54 | Line 58
58   #include "math/Vector3.hpp"
59   #include "primitives/Molecule.hpp"
60   #include "primitives/StuntDouble.hpp"
57 #include "UseTheForce/fCutoffPolicy.h"
58 #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
59 #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
60 #include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
61 #include "UseTheForce/doForces_interface.h"
62 #include "UseTheForce/DarkSide/neighborLists_interface.h"
63 #include "UseTheForce/DarkSide/electrostatic_interface.h"
64 #include "UseTheForce/DarkSide/switcheroo_interface.h"
61   #include "utils/MemoryUtils.hpp"
62   #include "utils/simError.h"
63   #include "selection/SelectionManager.hpp"
64   #include "io/ForceFieldOptions.hpp"
65 < #include "UseTheForce/ForceField.hpp"
65 > #include "brains/ForceField.hpp"
66 > #include "nonbonded/SwitchingFunction.hpp"
67  
68 <
69 < #ifdef IS_MPI
73 < #include "UseTheForce/mpiComponentPlan.h"
74 < #include "UseTheForce/DarkSide/simParallel_interface.h"
75 < #endif
76 <
77 < namespace oopse {
78 <  std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) {
79 <    std::map<int, std::set<int> >::iterator i = container.find(index);
80 <    std::set<int> result;
81 <    if (i != container.end()) {
82 <        result = i->second;
83 <    }
84 <
85 <    return result;
86 <  }
68 > using namespace std;
69 > namespace OpenMD {
70    
71    SimInfo::SimInfo(ForceField* ff, Globals* simParams) :
72      forceField_(ff), simParams_(simParams),
73 <    ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
73 >    nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0),
74 >    nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0),
75 >    nConstraints_(0), nFluctuatingCharges_(0),    
76      nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
77 <    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
78 <    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
79 <    nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
80 <    sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) {
77 >    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
78 >    nGlobalFluctuatingCharges_(0), nGlobalBonds_(0), nGlobalBends_(0),
79 >    nGlobalTorsions_(0), nGlobalInversions_(0), nGlobalConstraints_(0),
80 >    hasNGlobalConstraints_(false),
81 >    ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
82 >    sman_(NULL), topologyDone_(false), calcBoxDipole_(false),
83 >    calcBoxQuadrupole_(false), useAtomicVirial_(true) {    
84 >    
85 >    MoleculeStamp* molStamp;
86 >    int nMolWithSameStamp;
87 >    int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
88 >    int nGroups = 0;       //total cutoff groups defined in meta-data file
89 >    CutoffGroupStamp* cgStamp;    
90 >    RigidBodyStamp* rbStamp;
91 >    int nRigidAtoms = 0;
92 >    
93 >    vector<Component*> components = simParams->getComponents();
94 >    
95 >    for (vector<Component*>::iterator i = components.begin();
96 >         i !=components.end(); ++i) {
97 >      molStamp = (*i)->getMoleculeStamp();
98 >      if ( (*i)->haveRegion() ) {        
99 >        molStamp->setRegion( (*i)->getRegion() );
100 >      } else {
101 >        // set the region to a disallowed value:
102 >        molStamp->setRegion( -1 );
103 >      }
104  
105 <      MoleculeStamp* molStamp;
98 <      int nMolWithSameStamp;
99 <      int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
100 <      int nGroups = 0;      //total cutoff groups defined in meta-data file
101 <      CutoffGroupStamp* cgStamp;    
102 <      RigidBodyStamp* rbStamp;
103 <      int nRigidAtoms = 0;
104 <      std::vector<Component*> components = simParams->getComponents();
105 >      nMolWithSameStamp = (*i)->getNMol();
106        
107 <      for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
108 <        molStamp = (*i)->getMoleculeStamp();
109 <        nMolWithSameStamp = (*i)->getNMol();
110 <        
111 <        addMoleculeStamp(molStamp, nMolWithSameStamp);
112 <
113 <        //calculate atoms in molecules
114 <        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
115 <
116 <        //calculate atoms in cutoff groups
117 <        int nAtomsInGroups = 0;
118 <        int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
119 <        
120 <        for (int j=0; j < nCutoffGroupsInStamp; j++) {
121 <          cgStamp = molStamp->getCutoffGroupStamp(j);
122 <          nAtomsInGroups += cgStamp->getNMembers();
122 <        }
123 <
124 <        nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
125 <
126 <        nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
127 <
128 <        //calculate atoms in rigid bodies
129 <        int nAtomsInRigidBodies = 0;
130 <        int nRigidBodiesInStamp = molStamp->getNRigidBodies();
131 <        
132 <        for (int j=0; j < nRigidBodiesInStamp; j++) {
133 <          rbStamp = molStamp->getRigidBodyStamp(j);
134 <          nAtomsInRigidBodies += rbStamp->getNMembers();
135 <        }
136 <
137 <        nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
138 <        nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
139 <        
107 >      addMoleculeStamp(molStamp, nMolWithSameStamp);
108 >      
109 >      //calculate atoms in molecules
110 >      nGlobalAtoms_ += molStamp->getNAtoms() * nMolWithSameStamp;
111 >      nGlobalBonds_ += molStamp->getNBonds() * nMolWithSameStamp;
112 >      nGlobalBends_ += molStamp->getNBends() * nMolWithSameStamp;
113 >      nGlobalTorsions_ += molStamp->getNTorsions() * nMolWithSameStamp;
114 >      nGlobalInversions_ += molStamp->getNInversions() * nMolWithSameStamp;
115 >      
116 >      //calculate atoms in cutoff groups
117 >      int nAtomsInGroups = 0;
118 >      int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
119 >      
120 >      for (int j=0; j < nCutoffGroupsInStamp; j++) {
121 >        cgStamp = molStamp->getCutoffGroupStamp(j);
122 >        nAtomsInGroups += cgStamp->getNMembers();
123        }
124 <
125 <      //every free atom (atom does not belong to cutoff groups) is a cutoff
126 <      //group therefore the total number of cutoff groups in the system is
127 <      //equal to the total number of atoms minus number of atoms belong to
128 <      //cutoff group defined in meta-data file plus the number of cutoff
129 <      //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
133 <      //integrable object therefore the total number of integrable objects
134 <      //in the system is equal to the total number of atoms minus number of
135 <      //atoms belong to rigid body defined in meta-data file plus the number
136 <      //of rigid bodies defined in meta-data file
137 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
138 <                                                + nGlobalRigidBodies_;
139 <  
140 <      nGlobalMols_ = molStampIds_.size();
158 <
159 < #ifdef IS_MPI    
160 <      molToProcMap_.resize(nGlobalMols_);
161 < #endif
162 <
124 >      
125 >      nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
126 >      
127 >      nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
128 >      
129 >      //calculate atoms in rigid bodies
130 >      int nAtomsInRigidBodies = 0;
131 >      int nRigidBodiesInStamp = molStamp->getNRigidBodies();
132 >      
133 >      for (int j=0; j < nRigidBodiesInStamp; j++) {
134 >        rbStamp = molStamp->getRigidBodyStamp(j);
135 >        nAtomsInRigidBodies += rbStamp->getNMembers();
136 >      }
137 >      
138 >      nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
139 >      nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
140 >      
141      }
142 +    
143 +    //every free atom (atom does not belong to cutoff groups) is a cutoff
144 +    //group therefore the total number of cutoff groups in the system is
145 +    //equal to the total number of atoms minus number of atoms belong to
146 +    //cutoff group defined in meta-data file plus the number of cutoff
147 +    //groups defined in meta-data file
148  
149 +    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
150 +    
151 +    //every free atom (atom does not belong to rigid bodies) is an
152 +    //integrable object therefore the total number of integrable objects
153 +    //in the system is equal to the total number of atoms minus number of
154 +    //atoms belong to rigid body defined in meta-data file plus the number
155 +    //of rigid bodies defined in meta-data file
156 +    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
157 +      + nGlobalRigidBodies_;
158 +    
159 +    nGlobalMols_ = molStampIds_.size();
160 +    molToProcMap_.resize(nGlobalMols_);
161 +  }
162 +  
163    SimInfo::~SimInfo() {
164 <    std::map<int, Molecule*>::iterator i;
164 >    map<int, Molecule*>::iterator i;
165      for (i = molecules_.begin(); i != molecules_.end(); ++i) {
166        delete i->second;
167      }
# Line 174 | Line 172 | namespace oopse {
172      delete forceField_;
173    }
174  
177  int SimInfo::getNGlobalConstraints() {
178    int nGlobalConstraints;
179 #ifdef IS_MPI
180    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
181                  MPI_COMM_WORLD);    
182 #else
183    nGlobalConstraints =  nConstraints_;
184 #endif
185    return nGlobalConstraints;
186  }
175  
176    bool SimInfo::addMolecule(Molecule* mol) {
177      MoleculeIterator i;
178 <
178 >    
179      i = molecules_.find(mol->getGlobalIndex());
180      if (i == molecules_.end() ) {
181 <
182 <      molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
183 <        
181 >      
182 >      molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
183 >      
184        nAtoms_ += mol->getNAtoms();
185        nBonds_ += mol->getNBonds();
186        nBends_ += mol->getNBends();
187        nTorsions_ += mol->getNTorsions();
188 +      nInversions_ += mol->getNInversions();
189        nRigidBodies_ += mol->getNRigidBodies();
190        nIntegrableObjects_ += mol->getNIntegrableObjects();
191        nCutoffGroups_ += mol->getNCutoffGroups();
192        nConstraints_ += mol->getNConstraintPairs();
193 <
194 <      addExcludePairs(mol);
195 <        
193 >      
194 >      addInteractionPairs(mol);
195 >      
196        return true;
197      } else {
198        return false;
199      }
200    }
201 <
201 >  
202    bool SimInfo::removeMolecule(Molecule* mol) {
203      MoleculeIterator i;
204      i = molecules_.find(mol->getGlobalIndex());
# Line 222 | Line 211 | namespace oopse {
211        nBonds_ -= mol->getNBonds();
212        nBends_ -= mol->getNBends();
213        nTorsions_ -= mol->getNTorsions();
214 +      nInversions_ -= mol->getNInversions();
215        nRigidBodies_ -= mol->getNRigidBodies();
216        nIntegrableObjects_ -= mol->getNIntegrableObjects();
217        nCutoffGroups_ -= mol->getNCutoffGroups();
218        nConstraints_ -= mol->getNConstraintPairs();
219  
220 <      removeExcludePairs(mol);
220 >      removeInteractionPairs(mol);
221        molecules_.erase(mol->getGlobalIndex());
222  
223        delete mol;
# Line 236 | Line 226 | namespace oopse {
226      } else {
227        return false;
228      }
239
240
229    }    
230  
231          
# Line 253 | Line 241 | namespace oopse {
241  
242  
243    void SimInfo::calcNdf() {
244 <    int ndf_local;
244 >    int ndf_local, nfq_local;
245      MoleculeIterator i;
246 <    std::vector<StuntDouble*>::iterator j;
246 >    vector<StuntDouble*>::iterator j;
247 >    vector<Atom*>::iterator k;
248 >
249      Molecule* mol;
250 <    StuntDouble* integrableObject;
250 >    StuntDouble* sd;
251 >    Atom* atom;
252  
253      ndf_local = 0;
254 +    nfq_local = 0;
255      
256      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
265      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
266           integrableObject = mol->nextIntegrableObject(j)) {
257  
258 +      for (sd = mol->beginIntegrableObject(j); sd != NULL;
259 +           sd = mol->nextIntegrableObject(j)) {
260 +
261          ndf_local += 3;
262  
263 <        if (integrableObject->isDirectional()) {
264 <          if (integrableObject->isLinear()) {
263 >        if (sd->isDirectional()) {
264 >          if (sd->isLinear()) {
265              ndf_local += 2;
266            } else {
267              ndf_local += 3;
268            }
269          }
277            
270        }
271 +
272 +      for (atom = mol->beginFluctuatingCharge(k); atom != NULL;
273 +           atom = mol->nextFluctuatingCharge(k)) {
274 +        if (atom->isFluctuatingCharge()) {
275 +          nfq_local++;
276 +        }
277 +      }
278      }
279      
280 +    ndfLocal_ = ndf_local;
281 +
282      // n_constraints is local, so subtract them on each processor
283      ndf_local -= nConstraints_;
284  
285   #ifdef IS_MPI
286 <    MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
286 >    MPI_Allreduce(&ndf_local, &ndf_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
287 >    MPI_Allreduce(&nfq_local, &nGlobalFluctuatingCharges_, 1,
288 >      MPI_INT, MPI_SUM, MPI_COMM_WORLD);
289   #else
290      ndf_ = ndf_local;
291 +    nGlobalFluctuatingCharges_ = nfq_local;
292   #endif
293  
294      // nZconstraints_ is global, as are the 3 COM translations for the
# Line 295 | Line 299 | namespace oopse {
299  
300    int SimInfo::getFdf() {
301   #ifdef IS_MPI
302 <    MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
302 >    MPI_Allreduce(&fdf_local, &fdf_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
303   #else
304      fdf_ = fdf_local;
305   #endif
306      return fdf_;
307    }
308 +  
309 +  unsigned int SimInfo::getNLocalCutoffGroups(){
310 +    int nLocalCutoffAtoms = 0;
311 +    Molecule* mol;
312 +    MoleculeIterator mi;
313 +    CutoffGroup* cg;
314 +    Molecule::CutoffGroupIterator ci;
315      
316 +    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
317 +      
318 +      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
319 +           cg = mol->nextCutoffGroup(ci)) {
320 +        nLocalCutoffAtoms += cg->getNumAtom();
321 +        
322 +      }        
323 +    }
324 +    
325 +    return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_;
326 +  }
327 +    
328    void SimInfo::calcNdfRaw() {
329      int ndfRaw_local;
330  
331      MoleculeIterator i;
332 <    std::vector<StuntDouble*>::iterator j;
332 >    vector<StuntDouble*>::iterator j;
333      Molecule* mol;
334 <    StuntDouble* integrableObject;
334 >    StuntDouble* sd;
335  
336      // Raw degrees of freedom that we have to set
337      ndfRaw_local = 0;
338      
339      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
317      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
318           integrableObject = mol->nextIntegrableObject(j)) {
340  
341 +      for (sd = mol->beginIntegrableObject(j); sd != NULL;
342 +           sd = mol->nextIntegrableObject(j)) {
343 +
344          ndfRaw_local += 3;
345  
346 <        if (integrableObject->isDirectional()) {
347 <          if (integrableObject->isLinear()) {
346 >        if (sd->isDirectional()) {
347 >          if (sd->isLinear()) {
348              ndfRaw_local += 2;
349            } else {
350              ndfRaw_local += 3;
# Line 331 | Line 355 | namespace oopse {
355      }
356      
357   #ifdef IS_MPI
358 <    MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
358 >    MPI_Allreduce(&ndfRaw_local, &ndfRaw_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
359   #else
360      ndfRaw_ = ndfRaw_local;
361   #endif
# Line 342 | Line 366 | namespace oopse {
366  
367      ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
368  
345
369   #ifdef IS_MPI
370 <    MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
370 >    MPI_Allreduce(&ndfTrans_local, &ndfTrans_, 1, MPI_INT, MPI_SUM,
371 >                  MPI_COMM_WORLD);
372   #else
373      ndfTrans_ = ndfTrans_local;
374   #endif
375  
376      ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
353
377    }
378  
379 <  void SimInfo::addExcludePairs(Molecule* mol) {
380 <    std::vector<Bond*>::iterator bondIter;
381 <    std::vector<Bend*>::iterator bendIter;
382 <    std::vector<Torsion*>::iterator torsionIter;
379 >  void SimInfo::addInteractionPairs(Molecule* mol) {
380 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
381 >    vector<Bond*>::iterator bondIter;
382 >    vector<Bend*>::iterator bendIter;
383 >    vector<Torsion*>::iterator torsionIter;
384 >    vector<Inversion*>::iterator inversionIter;
385      Bond* bond;
386      Bend* bend;
387      Torsion* torsion;
388 +    Inversion* inversion;
389      int a;
390      int b;
391      int c;
392      int d;
393  
394 <    std::map<int, std::set<int> > atomGroups;
394 >    // atomGroups can be used to add special interaction maps between
395 >    // groups of atoms that are in two separate rigid bodies.
396 >    // However, most site-site interactions between two rigid bodies
397 >    // are probably not special, just the ones between the physically
398 >    // bonded atoms.  Interactions *within* a single rigid body should
399 >    // always be excluded.  These are done at the bottom of this
400 >    // function.
401  
402 +    map<int, set<int> > atomGroups;
403      Molecule::RigidBodyIterator rbIter;
404      RigidBody* rb;
405      Molecule::IntegrableObjectIterator ii;
406 <    StuntDouble* integrableObject;
406 >    StuntDouble* sd;
407      
408 <    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
409 <           integrableObject = mol->nextIntegrableObject(ii)) {
410 <
411 <      if (integrableObject->isRigidBody()) {
412 <          rb = static_cast<RigidBody*>(integrableObject);
413 <          std::vector<Atom*> atoms = rb->getAtoms();
414 <          std::set<int> rigidAtoms;
415 <          for (int i = 0; i < atoms.size(); ++i) {
416 <            rigidAtoms.insert(atoms[i]->getGlobalIndex());
417 <          }
418 <          for (int i = 0; i < atoms.size(); ++i) {
419 <            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
420 <          }      
408 >    for (sd = mol->beginIntegrableObject(ii); sd != NULL;
409 >         sd = mol->nextIntegrableObject(ii)) {
410 >      
411 >      if (sd->isRigidBody()) {
412 >        rb = static_cast<RigidBody*>(sd);
413 >        vector<Atom*> atoms = rb->getAtoms();
414 >        set<int> rigidAtoms;
415 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
416 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
417 >        }
418 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
419 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
420 >        }      
421        } else {
422 <        std::set<int> oneAtomSet;
423 <        oneAtomSet.insert(integrableObject->getGlobalIndex());
424 <        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
422 >        set<int> oneAtomSet;
423 >        oneAtomSet.insert(sd->getGlobalIndex());
424 >        atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet));        
425        }
426      }  
427  
428 <    
429 <    
430 <    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
428 >          
429 >    for (bond= mol->beginBond(bondIter); bond != NULL;
430 >         bond = mol->nextBond(bondIter)) {
431 >
432        a = bond->getAtomA()->getGlobalIndex();
433 <      b = bond->getAtomB()->getGlobalIndex();        
434 <      exclude_.addPair(a, b);
433 >      b = bond->getAtomB()->getGlobalIndex();  
434 >
435 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
436 >        oneTwoInteractions_.addPair(a, b);
437 >      } else {
438 >        excludedInteractions_.addPair(a, b);
439 >      }
440      }
441  
442 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
442 >    for (bend= mol->beginBend(bendIter); bend != NULL;
443 >         bend = mol->nextBend(bendIter)) {
444 >
445        a = bend->getAtomA()->getGlobalIndex();
446        b = bend->getAtomB()->getGlobalIndex();        
447        c = bend->getAtomC()->getGlobalIndex();
407      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
408      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
409      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
410
411      exclude_.addPairs(rigidSetA, rigidSetB);
412      exclude_.addPairs(rigidSetA, rigidSetC);
413      exclude_.addPairs(rigidSetB, rigidSetC);
448        
449 <      //exclude_.addPair(a, b);
450 <      //exclude_.addPair(a, c);
451 <      //exclude_.addPair(b, c);        
449 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
450 >        oneTwoInteractions_.addPair(a, b);      
451 >        oneTwoInteractions_.addPair(b, c);
452 >      } else {
453 >        excludedInteractions_.addPair(a, b);
454 >        excludedInteractions_.addPair(b, c);
455 >      }
456 >
457 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
458 >        oneThreeInteractions_.addPair(a, c);      
459 >      } else {
460 >        excludedInteractions_.addPair(a, c);
461 >      }
462      }
463  
464 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
464 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
465 >         torsion = mol->nextTorsion(torsionIter)) {
466 >
467        a = torsion->getAtomA()->getGlobalIndex();
468        b = torsion->getAtomB()->getGlobalIndex();        
469        c = torsion->getAtomC()->getGlobalIndex();        
470 <      d = torsion->getAtomD()->getGlobalIndex();        
425 <      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
426 <      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
427 <      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
428 <      std::set<int> rigidSetD = getRigidSet(d, atomGroups);
470 >      d = torsion->getAtomD()->getGlobalIndex();      
471  
472 <      exclude_.addPairs(rigidSetA, rigidSetB);
473 <      exclude_.addPairs(rigidSetA, rigidSetC);
474 <      exclude_.addPairs(rigidSetA, rigidSetD);
475 <      exclude_.addPairs(rigidSetB, rigidSetC);
476 <      exclude_.addPairs(rigidSetB, rigidSetD);
477 <      exclude_.addPairs(rigidSetC, rigidSetD);
472 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
473 >        oneTwoInteractions_.addPair(a, b);      
474 >        oneTwoInteractions_.addPair(b, c);
475 >        oneTwoInteractions_.addPair(c, d);
476 >      } else {
477 >        excludedInteractions_.addPair(a, b);
478 >        excludedInteractions_.addPair(b, c);
479 >        excludedInteractions_.addPair(c, d);
480 >      }
481  
482 <      /*
483 <      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
484 <      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
485 <      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
486 <      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
487 <      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
488 <      exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
489 <        
490 <      
491 <      exclude_.addPair(a, b);
492 <      exclude_.addPair(a, c);
493 <      exclude_.addPair(a, d);
494 <      exclude_.addPair(b, c);
450 <      exclude_.addPair(b, d);
451 <      exclude_.addPair(c, d);        
452 <      */
482 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
483 >        oneThreeInteractions_.addPair(a, c);      
484 >        oneThreeInteractions_.addPair(b, d);      
485 >      } else {
486 >        excludedInteractions_.addPair(a, c);
487 >        excludedInteractions_.addPair(b, d);
488 >      }
489 >
490 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
491 >        oneFourInteractions_.addPair(a, d);      
492 >      } else {
493 >        excludedInteractions_.addPair(a, d);
494 >      }
495      }
496  
497 <    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
498 <      std::vector<Atom*> atoms = rb->getAtoms();
499 <      for (int i = 0; i < atoms.size() -1 ; ++i) {
500 <        for (int j = i + 1; j < atoms.size(); ++j) {
497 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
498 >         inversion = mol->nextInversion(inversionIter)) {
499 >
500 >      a = inversion->getAtomA()->getGlobalIndex();
501 >      b = inversion->getAtomB()->getGlobalIndex();        
502 >      c = inversion->getAtomC()->getGlobalIndex();        
503 >      d = inversion->getAtomD()->getGlobalIndex();        
504 >
505 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
506 >        oneTwoInteractions_.addPair(a, b);      
507 >        oneTwoInteractions_.addPair(a, c);
508 >        oneTwoInteractions_.addPair(a, d);
509 >      } else {
510 >        excludedInteractions_.addPair(a, b);
511 >        excludedInteractions_.addPair(a, c);
512 >        excludedInteractions_.addPair(a, d);
513 >      }
514 >
515 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
516 >        oneThreeInteractions_.addPair(b, c);    
517 >        oneThreeInteractions_.addPair(b, d);    
518 >        oneThreeInteractions_.addPair(c, d);      
519 >      } else {
520 >        excludedInteractions_.addPair(b, c);
521 >        excludedInteractions_.addPair(b, d);
522 >        excludedInteractions_.addPair(c, d);
523 >      }
524 >    }
525 >
526 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
527 >         rb = mol->nextRigidBody(rbIter)) {
528 >      vector<Atom*> atoms = rb->getAtoms();
529 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
530 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
531            a = atoms[i]->getGlobalIndex();
532            b = atoms[j]->getGlobalIndex();
533 <          exclude_.addPair(a, b);
533 >          excludedInteractions_.addPair(a, b);
534          }
535        }
536      }        
537  
538    }
539  
540 <  void SimInfo::removeExcludePairs(Molecule* mol) {
541 <    std::vector<Bond*>::iterator bondIter;
542 <    std::vector<Bend*>::iterator bendIter;
543 <    std::vector<Torsion*>::iterator torsionIter;
540 >  void SimInfo::removeInteractionPairs(Molecule* mol) {
541 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
542 >    vector<Bond*>::iterator bondIter;
543 >    vector<Bend*>::iterator bendIter;
544 >    vector<Torsion*>::iterator torsionIter;
545 >    vector<Inversion*>::iterator inversionIter;
546      Bond* bond;
547      Bend* bend;
548      Torsion* torsion;
549 +    Inversion* inversion;
550      int a;
551      int b;
552      int c;
553      int d;
554  
555 <    std::map<int, std::set<int> > atomGroups;
481 <
555 >    map<int, set<int> > atomGroups;
556      Molecule::RigidBodyIterator rbIter;
557      RigidBody* rb;
558      Molecule::IntegrableObjectIterator ii;
559 <    StuntDouble* integrableObject;
559 >    StuntDouble* sd;
560      
561 <    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
562 <           integrableObject = mol->nextIntegrableObject(ii)) {
563 <
564 <      if (integrableObject->isRigidBody()) {
565 <          rb = static_cast<RigidBody*>(integrableObject);
566 <          std::vector<Atom*> atoms = rb->getAtoms();
567 <          std::set<int> rigidAtoms;
568 <          for (int i = 0; i < atoms.size(); ++i) {
569 <            rigidAtoms.insert(atoms[i]->getGlobalIndex());
570 <          }
571 <          for (int i = 0; i < atoms.size(); ++i) {
572 <            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
573 <          }      
561 >    for (sd = mol->beginIntegrableObject(ii); sd != NULL;
562 >         sd = mol->nextIntegrableObject(ii)) {
563 >      
564 >      if (sd->isRigidBody()) {
565 >        rb = static_cast<RigidBody*>(sd);
566 >        vector<Atom*> atoms = rb->getAtoms();
567 >        set<int> rigidAtoms;
568 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
569 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
570 >        }
571 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
572 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
573 >        }      
574        } else {
575 <        std::set<int> oneAtomSet;
576 <        oneAtomSet.insert(integrableObject->getGlobalIndex());
577 <        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
575 >        set<int> oneAtomSet;
576 >        oneAtomSet.insert(sd->getGlobalIndex());
577 >        atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet));        
578        }
579      }  
580  
581 <    
582 <    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
581 >    for (bond= mol->beginBond(bondIter); bond != NULL;
582 >         bond = mol->nextBond(bondIter)) {
583 >      
584        a = bond->getAtomA()->getGlobalIndex();
585 <      b = bond->getAtomB()->getGlobalIndex();        
586 <      exclude_.removePair(a, b);
585 >      b = bond->getAtomB()->getGlobalIndex();  
586 >    
587 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
588 >        oneTwoInteractions_.removePair(a, b);
589 >      } else {
590 >        excludedInteractions_.removePair(a, b);
591 >      }
592      }
593  
594 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
594 >    for (bend= mol->beginBend(bendIter); bend != NULL;
595 >         bend = mol->nextBend(bendIter)) {
596 >
597        a = bend->getAtomA()->getGlobalIndex();
598        b = bend->getAtomB()->getGlobalIndex();        
599        c = bend->getAtomC()->getGlobalIndex();
518
519      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
520      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
521      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
522
523      exclude_.removePairs(rigidSetA, rigidSetB);
524      exclude_.removePairs(rigidSetA, rigidSetC);
525      exclude_.removePairs(rigidSetB, rigidSetC);
600        
601 <      //exclude_.removePair(a, b);
602 <      //exclude_.removePair(a, c);
603 <      //exclude_.removePair(b, c);        
601 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
602 >        oneTwoInteractions_.removePair(a, b);      
603 >        oneTwoInteractions_.removePair(b, c);
604 >      } else {
605 >        excludedInteractions_.removePair(a, b);
606 >        excludedInteractions_.removePair(b, c);
607 >      }
608 >
609 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
610 >        oneThreeInteractions_.removePair(a, c);      
611 >      } else {
612 >        excludedInteractions_.removePair(a, c);
613 >      }
614      }
615  
616 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
616 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
617 >         torsion = mol->nextTorsion(torsionIter)) {
618 >
619        a = torsion->getAtomA()->getGlobalIndex();
620        b = torsion->getAtomB()->getGlobalIndex();        
621        c = torsion->getAtomC()->getGlobalIndex();        
622 <      d = torsion->getAtomD()->getGlobalIndex();        
622 >      d = torsion->getAtomD()->getGlobalIndex();      
623 >  
624 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
625 >        oneTwoInteractions_.removePair(a, b);      
626 >        oneTwoInteractions_.removePair(b, c);
627 >        oneTwoInteractions_.removePair(c, d);
628 >      } else {
629 >        excludedInteractions_.removePair(a, b);
630 >        excludedInteractions_.removePair(b, c);
631 >        excludedInteractions_.removePair(c, d);
632 >      }
633  
634 <      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
635 <      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
636 <      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
637 <      std::set<int> rigidSetD = getRigidSet(d, atomGroups);
634 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
635 >        oneThreeInteractions_.removePair(a, c);      
636 >        oneThreeInteractions_.removePair(b, d);      
637 >      } else {
638 >        excludedInteractions_.removePair(a, c);
639 >        excludedInteractions_.removePair(b, d);
640 >      }
641  
642 <      exclude_.removePairs(rigidSetA, rigidSetB);
643 <      exclude_.removePairs(rigidSetA, rigidSetC);
644 <      exclude_.removePairs(rigidSetA, rigidSetD);
645 <      exclude_.removePairs(rigidSetB, rigidSetC);
646 <      exclude_.removePairs(rigidSetB, rigidSetD);
647 <      exclude_.removePairs(rigidSetC, rigidSetD);
642 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
643 >        oneFourInteractions_.removePair(a, d);      
644 >      } else {
645 >        excludedInteractions_.removePair(a, d);
646 >      }
647 >    }
648  
649 <      /*
650 <      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
552 <      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
553 <      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
554 <      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
555 <      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
556 <      exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
649 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
650 >         inversion = mol->nextInversion(inversionIter)) {
651  
652 <      
653 <      exclude_.removePair(a, b);
654 <      exclude_.removePair(a, c);
655 <      exclude_.removePair(a, d);
656 <      exclude_.removePair(b, c);
657 <      exclude_.removePair(b, d);
658 <      exclude_.removePair(c, d);        
659 <      */
652 >      a = inversion->getAtomA()->getGlobalIndex();
653 >      b = inversion->getAtomB()->getGlobalIndex();        
654 >      c = inversion->getAtomC()->getGlobalIndex();        
655 >      d = inversion->getAtomD()->getGlobalIndex();        
656 >
657 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
658 >        oneTwoInteractions_.removePair(a, b);      
659 >        oneTwoInteractions_.removePair(a, c);
660 >        oneTwoInteractions_.removePair(a, d);
661 >      } else {
662 >        excludedInteractions_.removePair(a, b);
663 >        excludedInteractions_.removePair(a, c);
664 >        excludedInteractions_.removePair(a, d);
665 >      }
666 >
667 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
668 >        oneThreeInteractions_.removePair(b, c);    
669 >        oneThreeInteractions_.removePair(b, d);    
670 >        oneThreeInteractions_.removePair(c, d);      
671 >      } else {
672 >        excludedInteractions_.removePair(b, c);
673 >        excludedInteractions_.removePair(b, d);
674 >        excludedInteractions_.removePair(c, d);
675 >      }
676      }
677  
678 <    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
679 <      std::vector<Atom*> atoms = rb->getAtoms();
680 <      for (int i = 0; i < atoms.size() -1 ; ++i) {
681 <        for (int j = i + 1; j < atoms.size(); ++j) {
678 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
679 >         rb = mol->nextRigidBody(rbIter)) {
680 >      vector<Atom*> atoms = rb->getAtoms();
681 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
682 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
683            a = atoms[i]->getGlobalIndex();
684            b = atoms[j]->getGlobalIndex();
685 <          exclude_.removePair(a, b);
685 >          excludedInteractions_.removePair(a, b);
686          }
687        }
688      }        
689 <
689 >    
690    }
691 <
692 <
691 >  
692 >  
693    void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
694      int curStampId;
695 <
695 >    
696      //index from 0
697      curStampId = moleculeStamps_.size();
698  
# Line 589 | Line 700 | namespace oopse {
700      molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
701    }
702  
592  void SimInfo::update() {
703  
704 <    setupSimType();
705 <
706 < #ifdef IS_MPI
707 <    setupFortranParallel();
708 < #endif
709 <
710 <    setupFortranSim();
711 <
712 <    //setup fortran force field
713 <    /** @deprecate */    
604 <    int isError = 0;
605 <    
606 <    setupCutoff();
607 <    
608 <    setupElectrostaticSummationMethod( isError );
609 <    setupSwitchingFunction();
610 <    setupAccumulateBoxDipole();
611 <
612 <    if(isError){
613 <      sprintf( painCave.errMsg,
614 <               "ForceField error: There was an error initializing the forceField in fortran.\n" );
615 <      painCave.isFatal = 1;
616 <      simError();
617 <    }
618 <
704 >  /**
705 >   * update
706 >   *
707 >   *  Performs the global checks and variable settings after the
708 >   *  objects have been created.
709 >   *
710 >   */
711 >  void SimInfo::update() {  
712 >    setupSimVariables();
713 >    calcNConstraints();
714      calcNdf();
715      calcNdfRaw();
716      calcNdfTrans();
622
623    fortranInitialized_ = true;
717    }
718 <
719 <  std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
718 >  
719 >  /**
720 >   * getSimulatedAtomTypes
721 >   *
722 >   * Returns an STL set of AtomType* that are actually present in this
723 >   * simulation.  Must query all processors to assemble this information.
724 >   *
725 >   */
726 >  set<AtomType*> SimInfo::getSimulatedAtomTypes() {
727      SimInfo::MoleculeIterator mi;
728      Molecule* mol;
729      Molecule::AtomIterator ai;
730      Atom* atom;
731 <    std::set<AtomType*> atomTypes;
732 <
731 >    set<AtomType*> atomTypes;
732 >    
733      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
734 <
735 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
734 >      for(atom = mol->beginAtom(ai); atom != NULL;
735 >          atom = mol->nextAtom(ai)) {
736          atomTypes.insert(atom->getAtomType());
737 <      }
738 <        
739 <    }
737 >      }      
738 >    }    
739 >    
740 > #ifdef IS_MPI
741  
742 <    return atomTypes;        
743 <  }
643 <
644 <  void SimInfo::setupSimType() {
645 <    std::set<AtomType*>::iterator i;
646 <    std::set<AtomType*> atomTypes;
647 <    atomTypes = getUniqueAtomTypes();
742 >    // loop over the found atom types on this processor, and add their
743 >    // numerical idents to a vector:
744      
745 <    int useLennardJones = 0;
746 <    int useElectrostatic = 0;
747 <    int useEAM = 0;
748 <    int useSC = 0;
653 <    int useCharge = 0;
654 <    int useDirectional = 0;
655 <    int useDipole = 0;
656 <    int useGayBerne = 0;
657 <    int useSticky = 0;
658 <    int useStickyPower = 0;
659 <    int useShape = 0;
660 <    int useFLARB = 0; //it is not in AtomType yet
661 <    int useDirectionalAtom = 0;    
662 <    int useElectrostatics = 0;
663 <    //usePBC and useRF are from simParams
664 <    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
665 <    int useRF;
666 <    int useSF;
667 <    int useSP;
668 <    int useBoxDipole;
669 <    std::string myMethod;
745 >    vector<int> foundTypes;
746 >    set<AtomType*>::iterator i;
747 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i)
748 >      foundTypes.push_back( (*i)->getIdent() );
749  
750 <    // set the useRF logical
751 <    useRF = 0;
673 <    useSF = 0;
674 <    useSP = 0;
750 >    // count_local holds the number of found types on this processor
751 >    int count_local = foundTypes.size();
752  
753 +    int nproc;
754 +    MPI_Comm_size( MPI_COMM_WORLD, &nproc);
755  
756 <    if (simParams_->haveElectrostaticSummationMethod()) {
757 <      std::string myMethod = simParams_->getElectrostaticSummationMethod();
758 <      toUpper(myMethod);
759 <      if (myMethod == "REACTION_FIELD"){
681 <        useRF = 1;
682 <      } else if (myMethod == "SHIFTED_FORCE"){
683 <        useSF = 1;
684 <      } else if (myMethod == "SHIFTED_POTENTIAL"){
685 <        useSP = 1;
686 <      }
687 <    }
688 <    
689 <    if (simParams_->haveAccumulateBoxDipole())
690 <      if (simParams_->getAccumulateBoxDipole())
691 <        useBoxDipole = 1;
756 >    // we need arrays to hold the counts and displacement vectors for
757 >    // all processors
758 >    vector<int> counts(nproc, 0);
759 >    vector<int> disps(nproc, 0);
760  
761 <    //loop over all of the atom types
762 <    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
763 <      useLennardJones |= (*i)->isLennardJones();
764 <      useElectrostatic |= (*i)->isElectrostatic();
765 <      useEAM |= (*i)->isEAM();
766 <      useSC |= (*i)->isSC();
767 <      useCharge |= (*i)->isCharge();
768 <      useDirectional |= (*i)->isDirectional();
769 <      useDipole |= (*i)->isDipole();
770 <      useGayBerne |= (*i)->isGayBerne();
703 <      useSticky |= (*i)->isSticky();
704 <      useStickyPower |= (*i)->isStickyPower();
705 <      useShape |= (*i)->isShape();
761 >    // fill the counts array
762 >    MPI_Allgather(&count_local, 1, MPI_INT, &counts[0],
763 >                  1, MPI_INT, MPI_COMM_WORLD);
764 >  
765 >    // use the processor counts to compute the displacement array
766 >    disps[0] = 0;    
767 >    int totalCount = counts[0];
768 >    for (int iproc = 1; iproc < nproc; iproc++) {
769 >      disps[iproc] = disps[iproc-1] + counts[iproc-1];
770 >      totalCount += counts[iproc];
771      }
772  
773 <    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
774 <      useDirectionalAtom = 1;
775 <    }
773 >    // we need a (possibly redundant) set of all found types:
774 >    vector<int> ftGlobal(totalCount);
775 >    
776 >    // now spray out the foundTypes to all the other processors:    
777 >    MPI_Allgatherv(&foundTypes[0], count_local, MPI_INT,
778 >                   &ftGlobal[0], &counts[0], &disps[0],
779 >                   MPI_INT, MPI_COMM_WORLD);
780  
781 <    if (useCharge || useDipole) {
713 <      useElectrostatics = 1;
714 <    }
781 >    vector<int>::iterator j;
782  
783 < #ifdef IS_MPI    
784 <    int temp;
783 >    // foundIdents is a stl set, so inserting an already found ident
784 >    // will have no effect.
785 >    set<int> foundIdents;
786  
787 <    temp = usePBC;
788 <    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
787 >    for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j)
788 >      foundIdents.insert((*j));
789 >    
790 >    // now iterate over the foundIdents and get the actual atom types
791 >    // that correspond to these:
792 >    set<int>::iterator it;
793 >    for (it = foundIdents.begin(); it != foundIdents.end(); ++it)
794 >      atomTypes.insert( forceField_->getAtomType((*it)) );
795 >
796 > #endif
797  
798 <    temp = useDirectionalAtom;
799 <    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
798 >    return atomTypes;        
799 >  }
800  
725    temp = useLennardJones;
726    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
801  
802 <    temp = useElectrostatics;
803 <    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
802 >  int getGlobalCountOfType(AtomType* atype) {
803 >    /*
804 >    set<AtomType*> atypes = getSimulatedAtomTypes();
805 >    map<AtomType*, int> counts_;
806  
807 <    temp = useCharge;
808 <    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
809 <
810 <    temp = useDipole;
811 <    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
812 <
813 <    temp = useSticky;
814 <    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
807 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
808 >      for(atom = mol->beginAtom(ai); atom != NULL;
809 >          atom = mol->nextAtom(ai)) {
810 >        atom->getAtomType();
811 >      }      
812 >    }    
813 >    */
814 >    return 0;
815 >  }
816  
817 <    temp = useStickyPower;
818 <    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
817 >  void SimInfo::setupSimVariables() {
818 >    useAtomicVirial_ = simParams_->getUseAtomicVirial();
819 >    // we only call setAccumulateBoxDipole if the accumulateBoxDipole
820 >    // parameter is true
821 >    calcBoxDipole_ = false;
822 >    if ( simParams_->haveAccumulateBoxDipole() )
823 >      if ( simParams_->getAccumulateBoxDipole() ) {
824 >        calcBoxDipole_ = true;      
825 >      }
826 >    // we only call setAccumulateBoxQuadrupole if the accumulateBoxQuadrupole
827 >    // parameter is true
828 >    calcBoxQuadrupole_ = false;
829 >    if ( simParams_->haveAccumulateBoxQuadrupole() )
830 >      if ( simParams_->getAccumulateBoxQuadrupole() ) {
831 >        calcBoxQuadrupole_ = true;      
832 >      }
833      
834 <    temp = useGayBerne;
835 <    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
834 >    set<AtomType*>::iterator i;
835 >    set<AtomType*> atomTypes;
836 >    atomTypes = getSimulatedAtomTypes();    
837 >    bool usesElectrostatic = false;
838 >    bool usesMetallic = false;
839 >    bool usesDirectional = false;
840 >    bool usesFluctuatingCharges =  false;
841 >    //loop over all of the atom types
842 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
843 >      usesElectrostatic |= (*i)->isElectrostatic();
844 >      usesMetallic |= (*i)->isMetal();
845 >      usesDirectional |= (*i)->isDirectional();
846 >      usesFluctuatingCharges |= (*i)->isFluctuatingCharge();
847 >    }
848  
849 <    temp = useEAM;
850 <    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
849 > #ifdef IS_MPI
850 >    int temp;
851  
852 <    temp = useSC;
853 <    MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
852 >    temp = usesDirectional;
853 >    MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD);
854 >    usesDirectionalAtoms_ = (temp == 0) ? false : true;
855      
856 <    temp = useShape;
857 <    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
856 >    temp = usesMetallic;
857 >    MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD);
858 >    usesMetallicAtoms_ = (temp == 0) ? false : true;
859  
860 <    temp = useFLARB;
861 <    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
860 >    temp = usesElectrostatic;
861 >    MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD);
862 >    usesElectrostaticAtoms_ = (temp == 0) ? false : true;
863  
864 <    temp = useRF;
865 <    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
864 >    temp = usesFluctuatingCharges;
865 >    MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD);
866 >    usesFluctuatingCharges_ = (temp == 0) ? false : true;
867 > #else
868  
869 <    temp = useSF;
870 <    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
869 >    usesDirectionalAtoms_ = usesDirectional;
870 >    usesMetallicAtoms_ = usesMetallic;
871 >    usesElectrostaticAtoms_ = usesElectrostatic;
872 >    usesFluctuatingCharges_ = usesFluctuatingCharges;
873  
874 <    temp = useSP;
875 <    MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
874 > #endif
875 >    
876 >    requiresPrepair_ = usesMetallicAtoms_ ? true : false;
877 >    requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false;
878 >    requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false;    
879 >  }
880  
767    temp = useBoxDipole;
768    MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
881  
882 < #endif
882 >  vector<int> SimInfo::getGlobalAtomIndices() {
883 >    SimInfo::MoleculeIterator mi;
884 >    Molecule* mol;
885 >    Molecule::AtomIterator ai;
886 >    Atom* atom;
887  
888 <    fInfo_.SIM_uses_PBC = usePBC;    
889 <    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
890 <    fInfo_.SIM_uses_LennardJones = useLennardJones;
891 <    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
892 <    fInfo_.SIM_uses_Charges = useCharge;
893 <    fInfo_.SIM_uses_Dipoles = useDipole;
894 <    fInfo_.SIM_uses_Sticky = useSticky;
895 <    fInfo_.SIM_uses_StickyPower = useStickyPower;
896 <    fInfo_.SIM_uses_GayBerne = useGayBerne;
781 <    fInfo_.SIM_uses_EAM = useEAM;
782 <    fInfo_.SIM_uses_SC = useSC;
783 <    fInfo_.SIM_uses_Shapes = useShape;
784 <    fInfo_.SIM_uses_FLARB = useFLARB;
785 <    fInfo_.SIM_uses_RF = useRF;
786 <    fInfo_.SIM_uses_SF = useSF;
787 <    fInfo_.SIM_uses_SP = useSP;
788 <    fInfo_.SIM_uses_BoxDipole = useBoxDipole;
888 >    vector<int> GlobalAtomIndices(getNAtoms(), 0);
889 >    
890 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
891 >      
892 >      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
893 >        GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex();
894 >      }
895 >    }
896 >    return GlobalAtomIndices;
897    }
898  
791  void SimInfo::setupFortranSim() {
792    int isError;
793    int nExclude;
794    std::vector<int> fortranGlobalGroupMembership;
795    
796    nExclude = exclude_.getSize();
797    isError = 0;
899  
900 <    //globalGroupMembership_ is filled by SimCreator    
901 <    for (int i = 0; i < nGlobalAtoms_; i++) {
902 <      fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
900 >  vector<int> SimInfo::getGlobalGroupIndices() {
901 >    SimInfo::MoleculeIterator mi;
902 >    Molecule* mol;
903 >    Molecule::CutoffGroupIterator ci;
904 >    CutoffGroup* cg;
905 >
906 >    vector<int> GlobalGroupIndices;
907 >    
908 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
909 >      
910 >      //local index of cutoff group is trivial, it only depends on the
911 >      //order of travesing
912 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
913 >           cg = mol->nextCutoffGroup(ci)) {
914 >        GlobalGroupIndices.push_back(cg->getGlobalIndex());
915 >      }        
916      }
917 +    return GlobalGroupIndices;
918 +  }
919  
920 +
921 +  void SimInfo::prepareTopology() {
922 +
923      //calculate mass ratio of cutoff group
805    std::vector<RealType> mfact;
924      SimInfo::MoleculeIterator mi;
925      Molecule* mol;
926      Molecule::CutoffGroupIterator ci;
# Line 811 | Line 929 | namespace oopse {
929      Atom* atom;
930      RealType totalMass;
931  
932 <    //to avoid memory reallocation, reserve enough space for mfact
933 <    mfact.reserve(getNCutoffGroups());
932 >    /**
933 >     * The mass factor is the relative mass of an atom to the total
934 >     * mass of the cutoff group it belongs to.  By default, all atoms
935 >     * are their own cutoff groups, and therefore have mass factors of
936 >     * 1.  We need some special handling for massless atoms, which
937 >     * will be treated as carrying the entire mass of the cutoff
938 >     * group.
939 >     */
940 >    massFactors_.clear();
941 >    massFactors_.resize(getNAtoms(), 1.0);
942      
943      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
944 <      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
944 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
945 >           cg = mol->nextCutoffGroup(ci)) {
946  
947          totalMass = cg->getMass();
948          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
949            // Check for massless groups - set mfact to 1 if true
950 <          if (totalMass != 0)
951 <            mfact.push_back(atom->getMass()/totalMass);
950 >          if (totalMass != 0)
951 >            massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass;
952            else
953 <            mfact.push_back( 1.0 );
953 >            massFactors_[atom->getLocalIndex()] = 1.0;
954          }
828
955        }      
956      }
957  
958 <    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
833 <    std::vector<int> identArray;
834 <
835 <    //to avoid memory reallocation, reserve enough space identArray
836 <    identArray.reserve(getNAtoms());
837 <    
838 <    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
839 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
840 <        identArray.push_back(atom->getIdent());
841 <      }
842 <    }    
843 <
844 <    //fill molMembershipArray
845 <    //molMembershipArray is filled by SimCreator    
846 <    std::vector<int> molMembershipArray(nGlobalAtoms_);
847 <    for (int i = 0; i < nGlobalAtoms_; i++) {
848 <      molMembershipArray[i] = globalMolMembership_[i] + 1;
849 <    }
850 <    
851 <    //setup fortran simulation
852 <    int nGlobalExcludes = 0;
853 <    int* globalExcludes = NULL;
854 <    int* excludeList = exclude_.getExcludeList();
855 <    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
856 <                   &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
857 <                   &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
858 <
859 <    if( isError ){
860 <
861 <      sprintf( painCave.errMsg,
862 <               "There was an error setting the simulation information in fortran.\n" );
863 <      painCave.isFatal = 1;
864 <      painCave.severity = OOPSE_ERROR;
865 <      simError();
866 <    }
867 <
868 < #ifdef IS_MPI
869 <    sprintf( checkPointMsg,
870 <             "succesfully sent the simulation information to fortran.\n");
871 <    MPIcheckPoint();
872 < #endif // is_mpi
873 <
874 <    // Setup number of neighbors in neighbor list if present
875 <    if (simParams_->haveNeighborListNeighbors()) {
876 <      int nlistNeighbors = simParams_->getNeighborListNeighbors();
877 <      setNeighbors(&nlistNeighbors);
878 <    }
879 <  
880 <
881 <  }
882 <
883 <
884 < #ifdef IS_MPI
885 <  void SimInfo::setupFortranParallel() {
886 <    
887 <    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
888 <    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
889 <    std::vector<int> localToGlobalCutoffGroupIndex;
890 <    SimInfo::MoleculeIterator mi;
891 <    Molecule::AtomIterator ai;
892 <    Molecule::CutoffGroupIterator ci;
893 <    Molecule* mol;
894 <    Atom* atom;
895 <    CutoffGroup* cg;
896 <    mpiSimData parallelData;
897 <    int isError;
958 >    // Build the identArray_ and regions_
959  
960 <    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
961 <
962 <      //local index(index in DataStorge) of atom is important
963 <      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
964 <        localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
965 <      }
966 <
967 <      //local index of cutoff group is trivial, it only depends on the order of travesing
968 <      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
969 <        localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
909 <      }        
910 <        
911 <    }
912 <
913 <    //fill up mpiSimData struct
914 <    parallelData.nMolGlobal = getNGlobalMolecules();
915 <    parallelData.nMolLocal = getNMolecules();
916 <    parallelData.nAtomsGlobal = getNGlobalAtoms();
917 <    parallelData.nAtomsLocal = getNAtoms();
918 <    parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
919 <    parallelData.nGroupsLocal = getNCutoffGroups();
920 <    parallelData.myNode = worldRank;
921 <    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors));
922 <
923 <    //pass mpiSimData struct and index arrays to fortran
924 <    setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
925 <                    &localToGlobalAtomIndex[0],  &(parallelData.nGroupsLocal),
926 <                    &localToGlobalCutoffGroupIndex[0], &isError);
927 <
928 <    if (isError) {
929 <      sprintf(painCave.errMsg,
930 <              "mpiRefresh errror: fortran didn't like something we gave it.\n");
931 <      painCave.isFatal = 1;
932 <      simError();
933 <    }
934 <
935 <    sprintf(checkPointMsg, " mpiRefresh successful.\n");
936 <    MPIcheckPoint();
937 <
938 <
939 <  }
940 <
941 < #endif
942 <
943 <  void SimInfo::setupCutoff() {          
944 <    
945 <    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
946 <
947 <    // Check the cutoff policy
948 <    int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default
949 <
950 <    std::string myPolicy;
951 <    if (forceFieldOptions_.haveCutoffPolicy()){
952 <      myPolicy = forceFieldOptions_.getCutoffPolicy();
953 <    }else if (simParams_->haveCutoffPolicy()) {
954 <      myPolicy = simParams_->getCutoffPolicy();
955 <    }
956 <
957 <    if (!myPolicy.empty()){
958 <      toUpper(myPolicy);
959 <      if (myPolicy == "MIX") {
960 <        cp = MIX_CUTOFF_POLICY;
961 <      } else {
962 <        if (myPolicy == "MAX") {
963 <          cp = MAX_CUTOFF_POLICY;
964 <        } else {
965 <          if (myPolicy == "TRADITIONAL") {            
966 <            cp = TRADITIONAL_CUTOFF_POLICY;
967 <          } else {
968 <            // throw error        
969 <            sprintf( painCave.errMsg,
970 <                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
971 <            painCave.isFatal = 1;
972 <            simError();
973 <          }    
974 <        }          
975 <      }
976 <    }          
977 <    notifyFortranCutoffPolicy(&cp);
978 <
979 <    // Check the Skin Thickness for neighborlists
980 <    RealType skin;
981 <    if (simParams_->haveSkinThickness()) {
982 <      skin = simParams_->getSkinThickness();
983 <      notifyFortranSkinThickness(&skin);
984 <    }            
985 <        
986 <    // Check if the cutoff was set explicitly:
987 <    if (simParams_->haveCutoffRadius()) {
988 <      rcut_ = simParams_->getCutoffRadius();
989 <      if (simParams_->haveSwitchingRadius()) {
990 <        rsw_  = simParams_->getSwitchingRadius();
991 <      } else {
992 <        if (fInfo_.SIM_uses_Charges |
993 <            fInfo_.SIM_uses_Dipoles |
994 <            fInfo_.SIM_uses_RF) {
995 <          
996 <          rsw_ = 0.85 * rcut_;
997 <          sprintf(painCave.errMsg,
998 <                  "SimCreator Warning: No value was set for the switchingRadius.\n"
999 <                  "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n"
1000 <                  "\tswitchingRadius = %f. for this simulation\n", rsw_);
1001 <        painCave.isFatal = 0;
1002 <        simError();
1003 <        } else {
1004 <          rsw_ = rcut_;
1005 <          sprintf(painCave.errMsg,
1006 <                  "SimCreator Warning: No value was set for the switchingRadius.\n"
1007 <                  "\tOOPSE will use the same value as the cutoffRadius.\n"
1008 <                  "\tswitchingRadius = %f. for this simulation\n", rsw_);
1009 <          painCave.isFatal = 0;
1010 <          simError();
1011 <        }
1012 <      }
1013 <      
1014 <      notifyFortranCutoffs(&rcut_, &rsw_);
1015 <      
1016 <    } else {
1017 <      
1018 <      // For electrostatic atoms, we'll assume a large safe value:
1019 <      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
1020 <        sprintf(painCave.errMsg,
1021 <                "SimCreator Warning: No value was set for the cutoffRadius.\n"
1022 <                "\tOOPSE will use a default value of 15.0 angstroms"
1023 <                "\tfor the cutoffRadius.\n");
1024 <        painCave.isFatal = 0;
1025 <        simError();
1026 <        rcut_ = 15.0;
1027 <      
1028 <        if (simParams_->haveElectrostaticSummationMethod()) {
1029 <          std::string myMethod = simParams_->getElectrostaticSummationMethod();
1030 <          toUpper(myMethod);
1031 <          if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") {
1032 <            if (simParams_->haveSwitchingRadius()){
1033 <              sprintf(painCave.errMsg,
1034 <                      "SimInfo Warning: A value was set for the switchingRadius\n"
1035 <                      "\teven though the electrostaticSummationMethod was\n"
1036 <                      "\tset to %s\n", myMethod.c_str());
1037 <              painCave.isFatal = 1;
1038 <              simError();            
1039 <            }
1040 <          }
1041 <        }
1042 <      
1043 <        if (simParams_->haveSwitchingRadius()){
1044 <          rsw_ = simParams_->getSwitchingRadius();
1045 <        } else {        
1046 <          sprintf(painCave.errMsg,
1047 <                  "SimCreator Warning: No value was set for switchingRadius.\n"
1048 <                  "\tOOPSE will use a default value of\n"
1049 <                  "\t0.85 * cutoffRadius for the switchingRadius\n");
1050 <          painCave.isFatal = 0;
1051 <          simError();
1052 <          rsw_ = 0.85 * rcut_;
1053 <        }
1054 <        notifyFortranCutoffs(&rcut_, &rsw_);
1055 <      } else {
1056 <        // We didn't set rcut explicitly, and we don't have electrostatic atoms, so
1057 <        // We'll punt and let fortran figure out the cutoffs later.
1058 <        
1059 <        notifyFortranYouAreOnYourOwn();
1060 <
1061 <      }
1062 <    }
1063 <  }
1064 <
1065 <  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
1066 <    
1067 <    int errorOut;
1068 <    int esm =  NONE;
1069 <    int sm = UNDAMPED;
1070 <    RealType alphaVal;
1071 <    RealType dielectric;
1072 <    
1073 <    errorOut = isError;
1074 <
1075 <    if (simParams_->haveElectrostaticSummationMethod()) {
1076 <      std::string myMethod = simParams_->getElectrostaticSummationMethod();
1077 <      toUpper(myMethod);
1078 <      if (myMethod == "NONE") {
1079 <        esm = NONE;
1080 <      } else {
1081 <        if (myMethod == "SWITCHING_FUNCTION") {
1082 <          esm = SWITCHING_FUNCTION;
1083 <        } else {
1084 <          if (myMethod == "SHIFTED_POTENTIAL") {
1085 <            esm = SHIFTED_POTENTIAL;
1086 <          } else {
1087 <            if (myMethod == "SHIFTED_FORCE") {            
1088 <              esm = SHIFTED_FORCE;
1089 <            } else {
1090 <              if (myMethod == "REACTION_FIELD") {
1091 <                esm = REACTION_FIELD;
1092 <                dielectric = simParams_->getDielectric();
1093 <                if (!simParams_->haveDielectric()) {
1094 <                  // throw warning
1095 <                  sprintf( painCave.errMsg,
1096 <                           "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
1097 <                           "\tA default value of %f will be used for the dielectric.\n", dielectric);
1098 <                  painCave.isFatal = 0;
1099 <                  simError();
1100 <                }
1101 <              } else {
1102 <                // throw error        
1103 <                sprintf( painCave.errMsg,
1104 <                         "SimInfo error: Unknown electrostaticSummationMethod.\n"
1105 <                         "\t(Input file specified %s .)\n"
1106 <                         "\telectrostaticSummationMethod must be one of: \"none\",\n"
1107 <                         "\t\"shifted_potential\", \"shifted_force\", or \n"
1108 <                         "\t\"reaction_field\".\n", myMethod.c_str() );
1109 <                painCave.isFatal = 1;
1110 <                simError();
1111 <              }    
1112 <            }          
1113 <          }
1114 <        }
960 >    identArray_.clear();
961 >    identArray_.reserve(getNAtoms());  
962 >    regions_.clear();
963 >    regions_.reserve(getNAtoms());
964 >
965 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {      
966 >      int reg = mol->getRegion();      
967 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
968 >        identArray_.push_back(atom->getIdent());
969 >        regions_.push_back(reg);
970        }
971 <    }
972 <    
973 <    if (simParams_->haveElectrostaticScreeningMethod()) {
1119 <      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
1120 <      toUpper(myScreen);
1121 <      if (myScreen == "UNDAMPED") {
1122 <        sm = UNDAMPED;
1123 <      } else {
1124 <        if (myScreen == "DAMPED") {
1125 <          sm = DAMPED;
1126 <          if (!simParams_->haveDampingAlpha()) {
1127 <            // first set a cutoff dependent alpha value
1128 <            // we assume alpha depends linearly with rcut from 0 to 20.5 ang
1129 <            alphaVal = 0.5125 - rcut_* 0.025;
1130 <            // for values rcut > 20.5, alpha is zero
1131 <            if (alphaVal < 0) alphaVal = 0;
1132 <
1133 <            // throw warning
1134 <            sprintf( painCave.errMsg,
1135 <                     "SimInfo warning: dampingAlpha was not specified in the input file.\n"
1136 <                     "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_);
1137 <            painCave.isFatal = 0;
1138 <            simError();
1139 <          } else {
1140 <            alphaVal = simParams_->getDampingAlpha();
1141 <          }
1142 <          
1143 <        } else {
1144 <          // throw error        
1145 <          sprintf( painCave.errMsg,
1146 <                   "SimInfo error: Unknown electrostaticScreeningMethod.\n"
1147 <                   "\t(Input file specified %s .)\n"
1148 <                   "\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
1149 <                   "or \"damped\".\n", myScreen.c_str() );
1150 <          painCave.isFatal = 1;
1151 <          simError();
1152 <        }
1153 <      }
1154 <    }
1155 <    
1156 <    // let's pass some summation method variables to fortran
1157 <    setElectrostaticSummationMethod( &esm );
1158 <    setFortranElectrostaticMethod( &esm );
1159 <    setScreeningMethod( &sm );
1160 <    setDampingAlpha( &alphaVal );
1161 <    setReactionFieldDielectric( &dielectric );
1162 <    initFortranFF( &errorOut );
971 >    }    
972 >      
973 >    topologyDone_ = true;
974    }
975  
1165  void SimInfo::setupSwitchingFunction() {    
1166    int ft = CUBIC;
1167
1168    if (simParams_->haveSwitchingFunctionType()) {
1169      std::string funcType = simParams_->getSwitchingFunctionType();
1170      toUpper(funcType);
1171      if (funcType == "CUBIC") {
1172        ft = CUBIC;
1173      } else {
1174        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
1175          ft = FIFTH_ORDER_POLY;
1176        } else {
1177          // throw error        
1178          sprintf( painCave.errMsg,
1179                   "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
1180          painCave.isFatal = 1;
1181          simError();
1182        }          
1183      }
1184    }
1185
1186    // send switching function notification to switcheroo
1187    setFunctionType(&ft);
1188
1189  }
1190
1191  void SimInfo::setupAccumulateBoxDipole() {    
1192
1193    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
1194    if ( simParams_->haveAccumulateBoxDipole() )
1195      if ( simParams_->getAccumulateBoxDipole() ) {
1196        setAccumulateBoxDipole();
1197        calcBoxDipole_ = true;
1198      }
1199
1200  }
1201
976    void SimInfo::addProperty(GenericData* genData) {
977      properties_.addProperty(genData);  
978    }
979  
980 <  void SimInfo::removeProperty(const std::string& propName) {
980 >  void SimInfo::removeProperty(const string& propName) {
981      properties_.removeProperty(propName);  
982    }
983  
# Line 1211 | Line 985 | namespace oopse {
985      properties_.clearProperties();
986    }
987  
988 <  std::vector<std::string> SimInfo::getPropertyNames() {
988 >  vector<string> SimInfo::getPropertyNames() {
989      return properties_.getPropertyNames();  
990    }
991        
992 <  std::vector<GenericData*> SimInfo::getProperties() {
992 >  vector<GenericData*> SimInfo::getProperties() {
993      return properties_.getProperties();
994    }
995  
996 <  GenericData* SimInfo::getPropertyByName(const std::string& propName) {
996 >  GenericData* SimInfo::getPropertyByName(const string& propName) {
997      return properties_.getPropertyByName(propName);
998    }
999  
# Line 1230 | Line 1004 | namespace oopse {
1004      delete sman_;
1005      sman_ = sman;
1006  
1233    Molecule* mol;
1234    RigidBody* rb;
1235    Atom* atom;
1007      SimInfo::MoleculeIterator mi;
1008 +    Molecule::AtomIterator ai;
1009      Molecule::RigidBodyIterator rbIter;
1010 <    Molecule::AtomIterator atomIter;;
1010 >    Molecule::CutoffGroupIterator cgIter;
1011 >    Molecule::BondIterator bondIter;
1012 >    Molecule::BendIterator bendIter;
1013 >    Molecule::TorsionIterator torsionIter;
1014 >    Molecule::InversionIterator inversionIter;
1015  
1016 +    Molecule* mol;
1017 +    Atom* atom;
1018 +    RigidBody* rb;
1019 +    CutoffGroup* cg;
1020 +    Bond* bond;
1021 +    Bend* bend;
1022 +    Torsion* torsion;
1023 +    Inversion* inversion;    
1024 +
1025      for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
1026          
1027 <      for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
1027 >      for (atom = mol->beginAtom(ai); atom != NULL;
1028 >           atom = mol->nextAtom(ai)) {
1029          atom->setSnapshotManager(sman_);
1030 <      }
1031 <        
1032 <      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
1030 >      }        
1031 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
1032 >           rb = mol->nextRigidBody(rbIter)) {
1033          rb->setSnapshotManager(sman_);
1034        }
1035 <    }    
1036 <    
1035 >      for (cg = mol->beginCutoffGroup(cgIter); cg != NULL;
1036 >           cg = mol->nextCutoffGroup(cgIter)) {
1037 >        cg->setSnapshotManager(sman_);
1038 >      }
1039 >      for (bond = mol->beginBond(bondIter); bond != NULL;
1040 >           bond = mol->nextBond(bondIter)) {
1041 >        bond->setSnapshotManager(sman_);
1042 >      }
1043 >      for (bend = mol->beginBend(bendIter); bend != NULL;
1044 >           bend = mol->nextBend(bendIter)) {
1045 >        bend->setSnapshotManager(sman_);
1046 >      }
1047 >      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
1048 >           torsion = mol->nextTorsion(torsionIter)) {
1049 >        torsion->setSnapshotManager(sman_);
1050 >      }
1051 >      for (inversion = mol->beginInversion(inversionIter); inversion != NULL;
1052 >           inversion = mol->nextInversion(inversionIter)) {
1053 >        inversion->setSnapshotManager(sman_);
1054 >      }
1055 >    }
1056    }
1057  
1253  Vector3d SimInfo::getComVel(){
1254    SimInfo::MoleculeIterator i;
1255    Molecule* mol;
1058  
1059 <    Vector3d comVel(0.0);
1258 <    RealType totalMass = 0.0;
1259 <    
1260 <
1261 <    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1262 <      RealType mass = mol->getMass();
1263 <      totalMass += mass;
1264 <      comVel += mass * mol->getComVel();
1265 <    }  
1266 <
1267 < #ifdef IS_MPI
1268 <    RealType tmpMass = totalMass;
1269 <    Vector3d tmpComVel(comVel);    
1270 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1271 <    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1272 < #endif
1273 <
1274 <    comVel /= totalMass;
1275 <
1276 <    return comVel;
1277 <  }
1278 <
1279 <  Vector3d SimInfo::getCom(){
1280 <    SimInfo::MoleculeIterator i;
1281 <    Molecule* mol;
1282 <
1283 <    Vector3d com(0.0);
1284 <    RealType totalMass = 0.0;
1285 <    
1286 <    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1287 <      RealType mass = mol->getMass();
1288 <      totalMass += mass;
1289 <      com += mass * mol->getCom();
1290 <    }  
1291 <
1292 < #ifdef IS_MPI
1293 <    RealType tmpMass = totalMass;
1294 <    Vector3d tmpCom(com);    
1295 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1296 <    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1297 < #endif
1298 <
1299 <    com /= totalMass;
1300 <
1301 <    return com;
1059 >  ostream& operator <<(ostream& o, SimInfo& info) {
1060  
1303  }        
1304
1305  std::ostream& operator <<(std::ostream& o, SimInfo& info) {
1306
1061      return o;
1062    }
1063    
1064 <  
1311 <   /*
1312 <   Returns center of mass and center of mass velocity in one function call.
1313 <   */
1314 <  
1315 <   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1316 <      SimInfo::MoleculeIterator i;
1317 <      Molecule* mol;
1318 <      
1319 <    
1320 <      RealType totalMass = 0.0;
1321 <    
1322 <
1323 <      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1324 <         RealType mass = mol->getMass();
1325 <         totalMass += mass;
1326 <         com += mass * mol->getCom();
1327 <         comVel += mass * mol->getComVel();          
1328 <      }  
1329 <      
1330 < #ifdef IS_MPI
1331 <      RealType tmpMass = totalMass;
1332 <      Vector3d tmpCom(com);  
1333 <      Vector3d tmpComVel(comVel);
1334 <      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1335 <      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1336 <      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1337 < #endif
1338 <      
1339 <      com /= totalMass;
1340 <      comVel /= totalMass;
1341 <   }        
1342 <  
1343 <   /*
1344 <   Return intertia tensor for entire system and angular momentum Vector.
1345 <
1346 <
1347 <       [  Ixx -Ixy  -Ixz ]
1348 <  J =| -Iyx  Iyy  -Iyz |
1349 <       [ -Izx -Iyz   Izz ]
1350 <    */
1351 <
1352 <   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1353 <      
1354 <
1355 <      RealType xx = 0.0;
1356 <      RealType yy = 0.0;
1357 <      RealType zz = 0.0;
1358 <      RealType xy = 0.0;
1359 <      RealType xz = 0.0;
1360 <      RealType yz = 0.0;
1361 <      Vector3d com(0.0);
1362 <      Vector3d comVel(0.0);
1363 <      
1364 <      getComAll(com, comVel);
1365 <      
1366 <      SimInfo::MoleculeIterator i;
1367 <      Molecule* mol;
1368 <      
1369 <      Vector3d thisq(0.0);
1370 <      Vector3d thisv(0.0);
1371 <
1372 <      RealType thisMass = 0.0;
1373 <    
1374 <      
1375 <      
1376 <  
1377 <      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1378 <        
1379 <         thisq = mol->getCom()-com;
1380 <         thisv = mol->getComVel()-comVel;
1381 <         thisMass = mol->getMass();
1382 <         // Compute moment of intertia coefficients.
1383 <         xx += thisq[0]*thisq[0]*thisMass;
1384 <         yy += thisq[1]*thisq[1]*thisMass;
1385 <         zz += thisq[2]*thisq[2]*thisMass;
1386 <        
1387 <         // compute products of intertia
1388 <         xy += thisq[0]*thisq[1]*thisMass;
1389 <         xz += thisq[0]*thisq[2]*thisMass;
1390 <         yz += thisq[1]*thisq[2]*thisMass;
1391 <            
1392 <         angularMomentum += cross( thisq, thisv ) * thisMass;
1393 <            
1394 <      }  
1395 <      
1396 <      
1397 <      inertiaTensor(0,0) = yy + zz;
1398 <      inertiaTensor(0,1) = -xy;
1399 <      inertiaTensor(0,2) = -xz;
1400 <      inertiaTensor(1,0) = -xy;
1401 <      inertiaTensor(1,1) = xx + zz;
1402 <      inertiaTensor(1,2) = -yz;
1403 <      inertiaTensor(2,0) = -xz;
1404 <      inertiaTensor(2,1) = -yz;
1405 <      inertiaTensor(2,2) = xx + yy;
1406 <      
1407 < #ifdef IS_MPI
1408 <      Mat3x3d tmpI(inertiaTensor);
1409 <      Vector3d tmpAngMom;
1410 <      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1411 <      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1412 < #endif
1413 <              
1414 <      return;
1415 <   }
1416 <
1417 <   //Returns the angular momentum of the system
1418 <   Vector3d SimInfo::getAngularMomentum(){
1419 <      
1420 <      Vector3d com(0.0);
1421 <      Vector3d comVel(0.0);
1422 <      Vector3d angularMomentum(0.0);
1423 <      
1424 <      getComAll(com,comVel);
1425 <      
1426 <      SimInfo::MoleculeIterator i;
1427 <      Molecule* mol;
1428 <      
1429 <      Vector3d thisr(0.0);
1430 <      Vector3d thisp(0.0);
1431 <      
1432 <      RealType thisMass;
1433 <      
1434 <      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1435 <        thisMass = mol->getMass();
1436 <        thisr = mol->getCom()-com;
1437 <        thisp = (mol->getComVel()-comVel)*thisMass;
1438 <        
1439 <        angularMomentum += cross( thisr, thisp );
1440 <        
1441 <      }  
1442 <      
1443 < #ifdef IS_MPI
1444 <      Vector3d tmpAngMom;
1445 <      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1446 < #endif
1447 <      
1448 <      return angularMomentum;
1449 <   }
1450 <  
1064 >  
1065    StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
1066 <    return IOIndexToIntegrableObject.at(index);
1066 >    if (index >= int(IOIndexToIntegrableObject.size())) {
1067 >      sprintf(painCave.errMsg,
1068 >              "SimInfo::getIOIndexToIntegrableObject Error: Integrable Object\n"
1069 >              "\tindex exceeds number of known objects!\n");
1070 >      painCave.isFatal = 1;
1071 >      simError();
1072 >      return NULL;
1073 >    } else
1074 >      return IOIndexToIntegrableObject.at(index);
1075    }
1076    
1077 <  void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) {
1077 >  void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) {
1078      IOIndexToIntegrableObject= v;
1079    }
1080  
1081 <  /* Returns the Volume of the simulation based on a ellipsoid with semi-axes
1082 <     based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3
1083 <     where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to
1084 <     V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
1085 <  */
1086 <  void SimInfo::getGyrationalVolume(RealType &volume){
1087 <    Mat3x3d intTensor;
1466 <    RealType det;
1467 <    Vector3d dummyAngMom;
1468 <    RealType sysconstants;
1469 <    RealType geomCnst;
1470 <
1471 <    geomCnst = 3.0/2.0;
1472 <    /* Get the inertial tensor and angular momentum for free*/
1473 <    getInertiaTensor(intTensor,dummyAngMom);
1474 <    
1475 <    det = intTensor.determinant();
1476 <    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1477 <    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det);
1478 <    return;
1081 >  void SimInfo::calcNConstraints() {
1082 > #ifdef IS_MPI
1083 >    MPI_Allreduce(&nConstraints_, &nGlobalConstraints_, 1,  
1084 >                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1085 > #else
1086 >    nGlobalConstraints_ =  nConstraints_;
1087 > #endif
1088    }
1089  
1090 <  void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){
1482 <    Mat3x3d intTensor;
1483 <    Vector3d dummyAngMom;
1484 <    RealType sysconstants;
1485 <    RealType geomCnst;
1090 > }//end namespace OpenMD
1091  
1487    geomCnst = 3.0/2.0;
1488    /* Get the inertial tensor and angular momentum for free*/
1489    getInertiaTensor(intTensor,dummyAngMom);
1490    
1491    detI = intTensor.determinant();
1492    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1493    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI);
1494    return;
1495  }
1496 /*
1497   void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) {
1498      assert( v.size() == nAtoms_ + nRigidBodies_);
1499      sdByGlobalIndex_ = v;
1500    }
1501
1502    StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) {
1503      //assert(index < nAtoms_ + nRigidBodies_);
1504      return sdByGlobalIndex_.at(index);
1505    }  
1506 */  
1507 }//end namespace oopse
1508

Comparing trunk/src/brains/SimInfo.cpp (property svn:keywords):
Revision 1121 by chuckv, Mon Feb 26 04:45:42 2007 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines