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 331 by tim, Sun Feb 13 21:18:27 2005 UTC vs.
Revision 1938 by gezelter, Thu Oct 31 15:32:17 2013 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 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>
56  
57   #include "brains/SimInfo.hpp"
58   #include "math/Vector3.hpp"
59   #include "primitives/Molecule.hpp"
60 < #include "UseTheForce/doForces_interface.h"
56 < #include "UseTheForce/notifyCutoffs_interface.h"
60 > #include "primitives/StuntDouble.hpp"
61   #include "utils/MemoryUtils.hpp"
62   #include "utils/simError.h"
63   #include "selection/SelectionManager.hpp"
64 + #include "io/ForceFieldOptions.hpp"
65 + #include "brains/ForceField.hpp"
66 + #include "nonbonded/SwitchingFunction.hpp"
67  
68 < #ifdef IS_MPI
69 < #include "UseTheForce/mpiComponentPlan.h"
70 < #include "UseTheForce/DarkSide/simParallel_interface.h"
71 < #endif
72 <
73 < namespace oopse {
74 <
75 < SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
76 <                                ForceField* ff, Globals* simParams) :
77 <                                forceField_(ff), simParams_(simParams),
78 <                                ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
79 <                                nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
80 <                                nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
74 <                                nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
75 <                                nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
76 <                                sman_(NULL), fortranInitialized_(false), selectMan_(NULL) {
77 <
78 <            
79 <    std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
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),
74 >    nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
75 >    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nGlobalFluctuatingCharges_(0),
76 >    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nInversions_(0),
77 >    nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0),
78 >    nConstraints_(0), nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false),
79 >    calcBoxDipole_(false), useAtomicVirial_(true) {    
80 >    
81      MoleculeStamp* molStamp;
82      int nMolWithSameStamp;
83      int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
84 <    int nGroups = 0;          //total cutoff groups defined in meta-data file
84 >    int nGroups = 0;       //total cutoff groups defined in meta-data file
85      CutoffGroupStamp* cgStamp;    
86      RigidBodyStamp* rbStamp;
87      int nRigidAtoms = 0;
88      
89 <    for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
90 <        molStamp = i->first;
91 <        nMolWithSameStamp = i->second;
92 <        
93 <        addMoleculeStamp(molStamp, nMolWithSameStamp);
89 >    vector<Component*> components = simParams->getComponents();
90 >    
91 >    for (vector<Component*>::iterator i = components.begin();
92 >         i !=components.end(); ++i) {
93 >      molStamp = (*i)->getMoleculeStamp();
94 >      if ( (*i)->haveRegion() ) {        
95 >        molStamp->setRegion( (*i)->getRegion() );
96 >      } else {
97 >        // set the region to a disallowed value:
98 >        molStamp->setRegion( -1 );
99 >      }
100  
101 <        //calculate atoms in molecules
102 <        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
103 <
104 <
105 <        //calculate atoms in cutoff groups
106 <        int nAtomsInGroups = 0;
107 <        int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
108 <        
109 <        for (int j=0; j < nCutoffGroupsInStamp; j++) {
110 <            cgStamp = molStamp->getCutoffGroup(j);
111 <            nAtomsInGroups += cgStamp->getNMembers();
112 <        }
113 <
114 <        nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
115 <        nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
116 <
117 <        //calculate atoms in rigid bodies
118 <        int nAtomsInRigidBodies = 0;
119 <        int nRigidBodiesInStamp = molStamp->getNRigidBodies();
120 <        
121 <        for (int j=0; j < nRigidBodiesInStamp; j++) {
122 <            rbStamp = molStamp->getRigidBody(j);
123 <            nAtomsInRigidBodies += rbStamp->getNMembers();
124 <        }
125 <
126 <        nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
127 <        nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
128 <        
101 >      nMolWithSameStamp = (*i)->getNMol();
102 >      
103 >      addMoleculeStamp(molStamp, nMolWithSameStamp);
104 >      
105 >      //calculate atoms in molecules
106 >      nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
107 >      
108 >      //calculate atoms in cutoff groups
109 >      int nAtomsInGroups = 0;
110 >      int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
111 >      
112 >      for (int j=0; j < nCutoffGroupsInStamp; j++) {
113 >        cgStamp = molStamp->getCutoffGroupStamp(j);
114 >        nAtomsInGroups += cgStamp->getNMembers();
115 >      }
116 >      
117 >      nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
118 >      
119 >      nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
120 >      
121 >      //calculate atoms in rigid bodies
122 >      int nAtomsInRigidBodies = 0;
123 >      int nRigidBodiesInStamp = molStamp->getNRigidBodies();
124 >      
125 >      for (int j=0; j < nRigidBodiesInStamp; j++) {
126 >        rbStamp = molStamp->getRigidBodyStamp(j);
127 >        nAtomsInRigidBodies += rbStamp->getNMembers();
128 >      }
129 >      
130 >      nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
131 >      nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
132 >      
133      }
134 +    
135 +    //every free atom (atom does not belong to cutoff groups) is a cutoff
136 +    //group therefore the total number of cutoff groups in the system is
137 +    //equal to the total number of atoms minus number of atoms belong to
138 +    //cutoff group defined in meta-data file plus the number of cutoff
139 +    //groups defined in meta-data file
140  
124    //every free atom (atom does not belong to cutoff groups) is a cutoff group
125    //therefore the total number of cutoff groups in the system is equal to
126    //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
127    //file plus the number of cutoff groups defined in meta-data file
141      nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
142 <
143 <    //every free atom (atom does not belong to rigid bodies) is an integrable object
144 <    //therefore the total number of  integrable objects in the system is equal to
145 <    //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
146 <    //file plus the number of  rigid bodies defined in meta-data file
147 <    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
148 <
142 >    
143 >    //every free atom (atom does not belong to rigid bodies) is an
144 >    //integrable object therefore the total number of integrable objects
145 >    //in the system is equal to the total number of atoms minus number of
146 >    //atoms belong to rigid body defined in meta-data file plus the number
147 >    //of rigid bodies defined in meta-data file
148 >    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
149 >      + nGlobalRigidBodies_;
150 >    
151      nGlobalMols_ = molStampIds_.size();
137
138 #ifdef IS_MPI    
152      molToProcMap_.resize(nGlobalMols_);
153 < #endif
154 <
155 <    selectMan_ = new SelectionManager(this);
156 <    selectMan_->selectAll();
157 < }
158 <
159 < SimInfo::~SimInfo() {
160 <    //MemoryUtils::deleteVectorOfPointer(molecules_);
161 <
149 <    MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
150 <    
153 >  }
154 >  
155 >  SimInfo::~SimInfo() {
156 >    map<int, Molecule*>::iterator i;
157 >    for (i = molecules_.begin(); i != molecules_.end(); ++i) {
158 >      delete i->second;
159 >    }
160 >    molecules_.clear();
161 >      
162      delete sman_;
163      delete simParams_;
164      delete forceField_;
165 <    delete selectMan_;
155 < }
165 >  }
166  
157 int SimInfo::getNGlobalConstraints() {
158    int nGlobalConstraints;
159 #ifdef IS_MPI
160    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
161                  MPI_COMM_WORLD);    
162 #else
163    nGlobalConstraints =  nConstraints_;
164 #endif
165    return nGlobalConstraints;
166 }
167  
168 < bool SimInfo::addMolecule(Molecule* mol) {
168 >  bool SimInfo::addMolecule(Molecule* mol) {
169      MoleculeIterator i;
170 <
170 >    
171      i = molecules_.find(mol->getGlobalIndex());
172      if (i == molecules_.end() ) {
173 <
174 <        molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
175 <        
176 <        nAtoms_ += mol->getNAtoms();
177 <        nBonds_ += mol->getNBonds();
178 <        nBends_ += mol->getNBends();
179 <        nTorsions_ += mol->getNTorsions();
180 <        nRigidBodies_ += mol->getNRigidBodies();
181 <        nIntegrableObjects_ += mol->getNIntegrableObjects();
182 <        nCutoffGroups_ += mol->getNCutoffGroups();
183 <        nConstraints_ += mol->getNConstraintPairs();
184 <
185 <        addExcludePairs(mol);
186 <        
187 <        return true;
173 >      
174 >      molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
175 >      
176 >      nAtoms_ += mol->getNAtoms();
177 >      nBonds_ += mol->getNBonds();
178 >      nBends_ += mol->getNBends();
179 >      nTorsions_ += mol->getNTorsions();
180 >      nInversions_ += mol->getNInversions();
181 >      nRigidBodies_ += mol->getNRigidBodies();
182 >      nIntegrableObjects_ += mol->getNIntegrableObjects();
183 >      nCutoffGroups_ += mol->getNCutoffGroups();
184 >      nConstraints_ += mol->getNConstraintPairs();
185 >      
186 >      addInteractionPairs(mol);
187 >      
188 >      return true;
189      } else {
190 <        return false;
190 >      return false;
191      }
192 < }
193 <
194 < bool SimInfo::removeMolecule(Molecule* mol) {
192 >  }
193 >  
194 >  bool SimInfo::removeMolecule(Molecule* mol) {
195      MoleculeIterator i;
196      i = molecules_.find(mol->getGlobalIndex());
197  
198      if (i != molecules_.end() ) {
199  
200 <        assert(mol == i->second);
200 >      assert(mol == i->second);
201          
202 <        nAtoms_ -= mol->getNAtoms();
203 <        nBonds_ -= mol->getNBonds();
204 <        nBends_ -= mol->getNBends();
205 <        nTorsions_ -= mol->getNTorsions();
206 <        nRigidBodies_ -= mol->getNRigidBodies();
207 <        nIntegrableObjects_ -= mol->getNIntegrableObjects();
208 <        nCutoffGroups_ -= mol->getNCutoffGroups();
209 <        nConstraints_ -= mol->getNConstraintPairs();
202 >      nAtoms_ -= mol->getNAtoms();
203 >      nBonds_ -= mol->getNBonds();
204 >      nBends_ -= mol->getNBends();
205 >      nTorsions_ -= mol->getNTorsions();
206 >      nInversions_ -= mol->getNInversions();
207 >      nRigidBodies_ -= mol->getNRigidBodies();
208 >      nIntegrableObjects_ -= mol->getNIntegrableObjects();
209 >      nCutoffGroups_ -= mol->getNCutoffGroups();
210 >      nConstraints_ -= mol->getNConstraintPairs();
211  
212 <        removeExcludePairs(mol);
213 <        molecules_.erase(mol->getGlobalIndex());
212 >      removeInteractionPairs(mol);
213 >      molecules_.erase(mol->getGlobalIndex());
214  
215 <        delete mol;
215 >      delete mol;
216          
217 <        return true;
217 >      return true;
218      } else {
219 <        return false;
219 >      return false;
220      }
221 +  }    
222  
220
221 }    
222
223          
224 < Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
224 >  Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
225      i = molecules_.begin();
226      return i == molecules_.end() ? NULL : i->second;
227 < }    
227 >  }    
228  
229 < Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
229 >  Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
230      ++i;
231      return i == molecules_.end() ? NULL : i->second;    
232 < }
232 >  }
233  
234  
235 < void SimInfo::calcNdf() {
236 <    int ndf_local;
235 >  void SimInfo::calcNdf() {
236 >    int ndf_local, nfq_local;
237      MoleculeIterator i;
238 <    std::vector<StuntDouble*>::iterator j;
238 >    vector<StuntDouble*>::iterator j;
239 >    vector<Atom*>::iterator k;
240 >
241      Molecule* mol;
242 <    StuntDouble* integrableObject;
242 >    StuntDouble* sd;
243 >    Atom* atom;
244  
245      ndf_local = 0;
246 +    nfq_local = 0;
247      
248      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
245        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
246               integrableObject = mol->nextIntegrableObject(j)) {
249  
250 <            ndf_local += 3;
250 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
251 >           sd = mol->nextIntegrableObject(j)) {
252  
253 <            if (integrableObject->isDirectional()) {
254 <                if (integrableObject->isLinear()) {
255 <                    ndf_local += 2;
256 <                } else {
257 <                    ndf_local += 3;
258 <                }
259 <            }
260 <            
261 <        }//end for (integrableObject)
262 <    }// end for (mol)
253 >        ndf_local += 3;
254 >
255 >        if (sd->isDirectional()) {
256 >          if (sd->isLinear()) {
257 >            ndf_local += 2;
258 >          } else {
259 >            ndf_local += 3;
260 >          }
261 >        }
262 >      }
263 >
264 >      for (atom = mol->beginFluctuatingCharge(k); atom != NULL;
265 >           atom = mol->nextFluctuatingCharge(k)) {
266 >        if (atom->isFluctuatingCharge()) {
267 >          nfq_local++;
268 >        }
269 >      }
270 >    }
271      
272 +    ndfLocal_ = ndf_local;
273 +
274      // n_constraints is local, so subtract them on each processor
275      ndf_local -= nConstraints_;
276  
277   #ifdef IS_MPI
278 <    MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
278 >    MPI::COMM_WORLD.Allreduce(&ndf_local, &ndf_, 1, MPI::INT,MPI::SUM);
279 >    MPI::COMM_WORLD.Allreduce(&nfq_local, &nGlobalFluctuatingCharges_, 1,
280 >                              MPI::INT, MPI::SUM);
281   #else
282      ndf_ = ndf_local;
283 +    nGlobalFluctuatingCharges_ = nfq_local;
284   #endif
285  
286      // nZconstraints_ is global, as are the 3 COM translations for the
287      // entire system:
288      ndf_ = ndf_ - 3 - nZconstraint_;
289  
290 < }
290 >  }
291  
292 < void SimInfo::calcNdfRaw() {
292 >  int SimInfo::getFdf() {
293 > #ifdef IS_MPI
294 >    MPI::COMM_WORLD.Allreduce(&fdf_local, &fdf_, 1, MPI::INT, MPI::SUM);
295 > #else
296 >    fdf_ = fdf_local;
297 > #endif
298 >    return fdf_;
299 >  }
300 >  
301 >  unsigned int SimInfo::getNLocalCutoffGroups(){
302 >    int nLocalCutoffAtoms = 0;
303 >    Molecule* mol;
304 >    MoleculeIterator mi;
305 >    CutoffGroup* cg;
306 >    Molecule::CutoffGroupIterator ci;
307 >    
308 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
309 >      
310 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
311 >           cg = mol->nextCutoffGroup(ci)) {
312 >        nLocalCutoffAtoms += cg->getNumAtom();
313 >        
314 >      }        
315 >    }
316 >    
317 >    return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_;
318 >  }
319 >    
320 >  void SimInfo::calcNdfRaw() {
321      int ndfRaw_local;
322  
323      MoleculeIterator i;
324 <    std::vector<StuntDouble*>::iterator j;
324 >    vector<StuntDouble*>::iterator j;
325      Molecule* mol;
326 <    StuntDouble* integrableObject;
326 >    StuntDouble* sd;
327  
328      // Raw degrees of freedom that we have to set
329      ndfRaw_local = 0;
330      
331      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
288        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
289               integrableObject = mol->nextIntegrableObject(j)) {
332  
333 <            ndfRaw_local += 3;
333 >      for (sd = mol->beginIntegrableObject(j); sd != NULL;
334 >           sd = mol->nextIntegrableObject(j)) {
335  
336 <            if (integrableObject->isDirectional()) {
337 <                if (integrableObject->isLinear()) {
338 <                    ndfRaw_local += 2;
339 <                } else {
340 <                    ndfRaw_local += 3;
341 <                }
342 <            }
336 >        ndfRaw_local += 3;
337 >
338 >        if (sd->isDirectional()) {
339 >          if (sd->isLinear()) {
340 >            ndfRaw_local += 2;
341 >          } else {
342 >            ndfRaw_local += 3;
343 >          }
344 >        }
345              
346 <        }
346 >      }
347      }
348      
349   #ifdef IS_MPI
350 <    MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
350 >    MPI::COMM_WORLD.Allreduce(&ndfRaw_local, &ndfRaw_, 1, MPI::INT, MPI::SUM);
351   #else
352      ndfRaw_ = ndfRaw_local;
353   #endif
354 < }
354 >  }
355  
356 < void SimInfo::calcNdfTrans() {
356 >  void SimInfo::calcNdfTrans() {
357      int ndfTrans_local;
358  
359      ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
360  
361  
362   #ifdef IS_MPI
363 <    MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
363 >    MPI::COMM_WORLD.Allreduce(&ndfTrans_local, &ndfTrans_, 1,
364 >                              MPI::INT, MPI::SUM);
365   #else
366      ndfTrans_ = ndfTrans_local;
367   #endif
368  
369      ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
370  
371 < }
371 >  }
372  
373 < void SimInfo::addExcludePairs(Molecule* mol) {
374 <    std::vector<Bond*>::iterator bondIter;
375 <    std::vector<Bend*>::iterator bendIter;
376 <    std::vector<Torsion*>::iterator torsionIter;
373 >  void SimInfo::addInteractionPairs(Molecule* mol) {
374 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
375 >    vector<Atom*>::iterator atomIter;
376 >    vector<Bond*>::iterator bondIter;
377 >    vector<Bend*>::iterator bendIter;
378 >    vector<Torsion*>::iterator torsionIter;
379 >    vector<Inversion*>::iterator inversionIter;
380 >    Atom* atom;
381      Bond* bond;
382      Bend* bend;
383      Torsion* torsion;
384 +    Inversion* inversion;
385      int a;
386      int b;
387      int c;
388      int d;
389 +
390 +    // atomGroups can be used to add special interaction maps between
391 +    // groups of atoms that are in two separate rigid bodies.
392 +    // However, most site-site interactions between two rigid bodies
393 +    // are probably not special, just the ones between the physically
394 +    // bonded atoms.  Interactions *within* a single rigid body should
395 +    // always be excluded.  These are done at the bottom of this
396 +    // function.
397 +
398 +    map<int, set<int> > atomGroups;
399 +    Molecule::RigidBodyIterator rbIter;
400 +    RigidBody* rb;
401 +    Molecule::IntegrableObjectIterator ii;
402 +    StuntDouble* sd;
403      
404 <    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
405 <        a = bond->getAtomA()->getGlobalIndex();
406 <        b = bond->getAtomB()->getGlobalIndex();        
407 <        exclude_.addPair(a, b);
404 >    for (sd = mol->beginIntegrableObject(ii); sd != NULL;
405 >         sd = mol->nextIntegrableObject(ii)) {
406 >      
407 >      if (sd->isRigidBody()) {
408 >        rb = static_cast<RigidBody*>(sd);
409 >        vector<Atom*> atoms = rb->getAtoms();
410 >        set<int> rigidAtoms;
411 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
412 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
413 >        }
414 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
415 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
416 >        }      
417 >      } else {
418 >        set<int> oneAtomSet;
419 >        oneAtomSet.insert(sd->getGlobalIndex());
420 >        atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet));        
421 >      }
422 >    }  
423 >
424 >          
425 >    for (bond= mol->beginBond(bondIter); bond != NULL;
426 >         bond = mol->nextBond(bondIter)) {
427 >
428 >      a = bond->getAtomA()->getGlobalIndex();
429 >      b = bond->getAtomB()->getGlobalIndex();  
430 >
431 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
432 >        oneTwoInteractions_.addPair(a, b);
433 >      } else {
434 >        excludedInteractions_.addPair(a, b);
435 >      }
436      }
437  
438 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
439 <        a = bend->getAtomA()->getGlobalIndex();
347 <        b = bend->getAtomB()->getGlobalIndex();        
348 <        c = bend->getAtomC()->getGlobalIndex();
438 >    for (bend= mol->beginBend(bendIter); bend != NULL;
439 >         bend = mol->nextBend(bendIter)) {
440  
441 <        exclude_.addPair(a, b);
442 <        exclude_.addPair(a, c);
443 <        exclude_.addPair(b, c);        
441 >      a = bend->getAtomA()->getGlobalIndex();
442 >      b = bend->getAtomB()->getGlobalIndex();        
443 >      c = bend->getAtomC()->getGlobalIndex();
444 >      
445 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
446 >        oneTwoInteractions_.addPair(a, b);      
447 >        oneTwoInteractions_.addPair(b, c);
448 >      } else {
449 >        excludedInteractions_.addPair(a, b);
450 >        excludedInteractions_.addPair(b, c);
451 >      }
452 >
453 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
454 >        oneThreeInteractions_.addPair(a, c);      
455 >      } else {
456 >        excludedInteractions_.addPair(a, c);
457 >      }
458      }
459  
460 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
461 <        a = torsion->getAtomA()->getGlobalIndex();
357 <        b = torsion->getAtomB()->getGlobalIndex();        
358 <        c = torsion->getAtomC()->getGlobalIndex();        
359 <        d = torsion->getAtomD()->getGlobalIndex();        
460 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
461 >         torsion = mol->nextTorsion(torsionIter)) {
462  
463 <        exclude_.addPair(a, b);
464 <        exclude_.addPair(a, c);
465 <        exclude_.addPair(a, d);
466 <        exclude_.addPair(b, c);
467 <        exclude_.addPair(b, d);
468 <        exclude_.addPair(c, d);        
463 >      a = torsion->getAtomA()->getGlobalIndex();
464 >      b = torsion->getAtomB()->getGlobalIndex();        
465 >      c = torsion->getAtomC()->getGlobalIndex();        
466 >      d = torsion->getAtomD()->getGlobalIndex();      
467 >
468 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
469 >        oneTwoInteractions_.addPair(a, b);      
470 >        oneTwoInteractions_.addPair(b, c);
471 >        oneTwoInteractions_.addPair(c, d);
472 >      } else {
473 >        excludedInteractions_.addPair(a, b);
474 >        excludedInteractions_.addPair(b, c);
475 >        excludedInteractions_.addPair(c, d);
476 >      }
477 >
478 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
479 >        oneThreeInteractions_.addPair(a, c);      
480 >        oneThreeInteractions_.addPair(b, d);      
481 >      } else {
482 >        excludedInteractions_.addPair(a, c);
483 >        excludedInteractions_.addPair(b, d);
484 >      }
485 >
486 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
487 >        oneFourInteractions_.addPair(a, d);      
488 >      } else {
489 >        excludedInteractions_.addPair(a, d);
490 >      }
491      }
492  
493 <    
494 < }
493 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
494 >         inversion = mol->nextInversion(inversionIter)) {
495  
496 < void SimInfo::removeExcludePairs(Molecule* mol) {
497 <    std::vector<Bond*>::iterator bondIter;
498 <    std::vector<Bend*>::iterator bendIter;
499 <    std::vector<Torsion*>::iterator torsionIter;
496 >      a = inversion->getAtomA()->getGlobalIndex();
497 >      b = inversion->getAtomB()->getGlobalIndex();        
498 >      c = inversion->getAtomC()->getGlobalIndex();        
499 >      d = inversion->getAtomD()->getGlobalIndex();        
500 >
501 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
502 >        oneTwoInteractions_.addPair(a, b);      
503 >        oneTwoInteractions_.addPair(a, c);
504 >        oneTwoInteractions_.addPair(a, d);
505 >      } else {
506 >        excludedInteractions_.addPair(a, b);
507 >        excludedInteractions_.addPair(a, c);
508 >        excludedInteractions_.addPair(a, d);
509 >      }
510 >
511 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
512 >        oneThreeInteractions_.addPair(b, c);    
513 >        oneThreeInteractions_.addPair(b, d);    
514 >        oneThreeInteractions_.addPair(c, d);      
515 >      } else {
516 >        excludedInteractions_.addPair(b, c);
517 >        excludedInteractions_.addPair(b, d);
518 >        excludedInteractions_.addPair(c, d);
519 >      }
520 >    }
521 >
522 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
523 >         rb = mol->nextRigidBody(rbIter)) {
524 >      vector<Atom*> atoms = rb->getAtoms();
525 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
526 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
527 >          a = atoms[i]->getGlobalIndex();
528 >          b = atoms[j]->getGlobalIndex();
529 >          excludedInteractions_.addPair(a, b);
530 >        }
531 >      }
532 >    }        
533 >
534 >  }
535 >
536 >  void SimInfo::removeInteractionPairs(Molecule* mol) {
537 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
538 >    vector<Bond*>::iterator bondIter;
539 >    vector<Bend*>::iterator bendIter;
540 >    vector<Torsion*>::iterator torsionIter;
541 >    vector<Inversion*>::iterator inversionIter;
542      Bond* bond;
543      Bend* bend;
544      Torsion* torsion;
545 +    Inversion* inversion;
546      int a;
547      int b;
548      int c;
549      int d;
550 +
551 +    map<int, set<int> > atomGroups;
552 +    Molecule::RigidBodyIterator rbIter;
553 +    RigidBody* rb;
554 +    Molecule::IntegrableObjectIterator ii;
555 +    StuntDouble* sd;
556      
557 <    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
558 <        a = bond->getAtomA()->getGlobalIndex();
559 <        b = bond->getAtomB()->getGlobalIndex();        
560 <        exclude_.removePair(a, b);
557 >    for (sd = mol->beginIntegrableObject(ii); sd != NULL;
558 >         sd = mol->nextIntegrableObject(ii)) {
559 >      
560 >      if (sd->isRigidBody()) {
561 >        rb = static_cast<RigidBody*>(sd);
562 >        vector<Atom*> atoms = rb->getAtoms();
563 >        set<int> rigidAtoms;
564 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
565 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
566 >        }
567 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
568 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
569 >        }      
570 >      } else {
571 >        set<int> oneAtomSet;
572 >        oneAtomSet.insert(sd->getGlobalIndex());
573 >        atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet));        
574 >      }
575 >    }  
576 >
577 >    for (bond= mol->beginBond(bondIter); bond != NULL;
578 >         bond = mol->nextBond(bondIter)) {
579 >      
580 >      a = bond->getAtomA()->getGlobalIndex();
581 >      b = bond->getAtomB()->getGlobalIndex();  
582 >    
583 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
584 >        oneTwoInteractions_.removePair(a, b);
585 >      } else {
586 >        excludedInteractions_.removePair(a, b);
587 >      }
588      }
589  
590 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
591 <        a = bend->getAtomA()->getGlobalIndex();
392 <        b = bend->getAtomB()->getGlobalIndex();        
393 <        c = bend->getAtomC()->getGlobalIndex();
590 >    for (bend= mol->beginBend(bendIter); bend != NULL;
591 >         bend = mol->nextBend(bendIter)) {
592  
593 <        exclude_.removePair(a, b);
594 <        exclude_.removePair(a, c);
595 <        exclude_.removePair(b, c);        
593 >      a = bend->getAtomA()->getGlobalIndex();
594 >      b = bend->getAtomB()->getGlobalIndex();        
595 >      c = bend->getAtomC()->getGlobalIndex();
596 >      
597 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
598 >        oneTwoInteractions_.removePair(a, b);      
599 >        oneTwoInteractions_.removePair(b, c);
600 >      } else {
601 >        excludedInteractions_.removePair(a, b);
602 >        excludedInteractions_.removePair(b, c);
603 >      }
604 >
605 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
606 >        oneThreeInteractions_.removePair(a, c);      
607 >      } else {
608 >        excludedInteractions_.removePair(a, c);
609 >      }
610      }
611  
612 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
613 <        a = torsion->getAtomA()->getGlobalIndex();
402 <        b = torsion->getAtomB()->getGlobalIndex();        
403 <        c = torsion->getAtomC()->getGlobalIndex();        
404 <        d = torsion->getAtomD()->getGlobalIndex();        
612 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
613 >         torsion = mol->nextTorsion(torsionIter)) {
614  
615 <        exclude_.removePair(a, b);
616 <        exclude_.removePair(a, c);
617 <        exclude_.removePair(a, d);
618 <        exclude_.removePair(b, c);
619 <        exclude_.removePair(b, d);
620 <        exclude_.removePair(c, d);        
615 >      a = torsion->getAtomA()->getGlobalIndex();
616 >      b = torsion->getAtomB()->getGlobalIndex();        
617 >      c = torsion->getAtomC()->getGlobalIndex();        
618 >      d = torsion->getAtomD()->getGlobalIndex();      
619 >  
620 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
621 >        oneTwoInteractions_.removePair(a, b);      
622 >        oneTwoInteractions_.removePair(b, c);
623 >        oneTwoInteractions_.removePair(c, d);
624 >      } else {
625 >        excludedInteractions_.removePair(a, b);
626 >        excludedInteractions_.removePair(b, c);
627 >        excludedInteractions_.removePair(c, d);
628 >      }
629 >
630 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
631 >        oneThreeInteractions_.removePair(a, c);      
632 >        oneThreeInteractions_.removePair(b, d);      
633 >      } else {
634 >        excludedInteractions_.removePair(a, c);
635 >        excludedInteractions_.removePair(b, d);
636 >      }
637 >
638 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
639 >        oneFourInteractions_.removePair(a, d);      
640 >      } else {
641 >        excludedInteractions_.removePair(a, d);
642 >      }
643      }
644  
645 < }
645 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
646 >         inversion = mol->nextInversion(inversionIter)) {
647  
648 +      a = inversion->getAtomA()->getGlobalIndex();
649 +      b = inversion->getAtomB()->getGlobalIndex();        
650 +      c = inversion->getAtomC()->getGlobalIndex();        
651 +      d = inversion->getAtomD()->getGlobalIndex();        
652  
653 < void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
654 <    int curStampId;
653 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
654 >        oneTwoInteractions_.removePair(a, b);      
655 >        oneTwoInteractions_.removePair(a, c);
656 >        oneTwoInteractions_.removePair(a, d);
657 >      } else {
658 >        excludedInteractions_.removePair(a, b);
659 >        excludedInteractions_.removePair(a, c);
660 >        excludedInteractions_.removePair(a, d);
661 >      }
662  
663 +      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
664 +        oneThreeInteractions_.removePair(b, c);    
665 +        oneThreeInteractions_.removePair(b, d);    
666 +        oneThreeInteractions_.removePair(c, d);      
667 +      } else {
668 +        excludedInteractions_.removePair(b, c);
669 +        excludedInteractions_.removePair(b, d);
670 +        excludedInteractions_.removePair(c, d);
671 +      }
672 +    }
673 +
674 +    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
675 +         rb = mol->nextRigidBody(rbIter)) {
676 +      vector<Atom*> atoms = rb->getAtoms();
677 +      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
678 +        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
679 +          a = atoms[i]->getGlobalIndex();
680 +          b = atoms[j]->getGlobalIndex();
681 +          excludedInteractions_.removePair(a, b);
682 +        }
683 +      }
684 +    }        
685 +    
686 +  }
687 +  
688 +  
689 +  void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
690 +    int curStampId;
691 +    
692      //index from 0
693      curStampId = moleculeStamps_.size();
694  
695      moleculeStamps_.push_back(molStamp);
696      molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
697 < }
697 >  }
698  
427 void SimInfo::update() {
699  
700 <    setupSimType();
701 <
702 < #ifdef IS_MPI
703 <    setupFortranParallel();
704 < #endif
705 <
706 <    setupFortranSim();
707 <
708 <    //setup fortran force field
438 <    /** @deprecate */    
439 <    int isError = 0;
440 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
441 <    if(isError){
442 <        sprintf( painCave.errMsg,
443 <         "ForceField error: There was an error initializing the forceField in fortran.\n" );
444 <        painCave.isFatal = 1;
445 <        simError();
446 <    }
447 <  
448 <    
449 <    setupCutoff();
450 <
700 >  /**
701 >   * update
702 >   *
703 >   *  Performs the global checks and variable settings after the
704 >   *  objects have been created.
705 >   *
706 >   */
707 >  void SimInfo::update() {  
708 >    setupSimVariables();
709      calcNdf();
710      calcNdfRaw();
711      calcNdfTrans();
712 <
713 <    fortranInitialized_ = true;
714 < }
715 <
716 < std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
712 >  }
713 >  
714 >  /**
715 >   * getSimulatedAtomTypes
716 >   *
717 >   * Returns an STL set of AtomType* that are actually present in this
718 >   * simulation.  Must query all processors to assemble this information.
719 >   *
720 >   */
721 >  set<AtomType*> SimInfo::getSimulatedAtomTypes() {
722      SimInfo::MoleculeIterator mi;
723      Molecule* mol;
724      Molecule::AtomIterator ai;
725      Atom* atom;
726 <    std::set<AtomType*> atomTypes;
727 <
726 >    set<AtomType*> atomTypes;
727 >    
728      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
729 +      for(atom = mol->beginAtom(ai); atom != NULL;
730 +          atom = mol->nextAtom(ai)) {
731 +        atomTypes.insert(atom->getAtomType());
732 +      }      
733 +    }    
734 +    
735 + #ifdef IS_MPI
736  
737 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
738 <            atomTypes.insert(atom->getAtomType());
739 <        }
740 <        
741 <    }
737 >    // loop over the found atom types on this processor, and add their
738 >    // numerical idents to a vector:
739 >    
740 >    vector<int> foundTypes;
741 >    set<AtomType*>::iterator i;
742 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i)
743 >      foundTypes.push_back( (*i)->getIdent() );
744  
745 <    return atomTypes;        
746 < }
745 >    // count_local holds the number of found types on this processor
746 >    int count_local = foundTypes.size();
747  
748 < void SimInfo::setupSimType() {
477 <    std::set<AtomType*>::iterator i;
478 <    std::set<AtomType*> atomTypes;
479 <    atomTypes = getUniqueAtomTypes();
480 <    
481 <    int useLennardJones = 0;
482 <    int useElectrostatic = 0;
483 <    int useEAM = 0;
484 <    int useCharge = 0;
485 <    int useDirectional = 0;
486 <    int useDipole = 0;
487 <    int useGayBerne = 0;
488 <    int useSticky = 0;
489 <    int useShape = 0;
490 <    int useFLARB = 0; //it is not in AtomType yet
491 <    int useDirectionalAtom = 0;    
492 <    int useElectrostatics = 0;
493 <    //usePBC and useRF are from simParams
494 <    int usePBC = simParams_->getPBC();
495 <    int useRF = simParams_->getUseRF();
748 >    int nproc = MPI::COMM_WORLD.Get_size();
749  
750 <    //loop over all of the atom types
751 <    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
752 <        useLennardJones |= (*i)->isLennardJones();
753 <        useElectrostatic |= (*i)->isElectrostatic();
501 <        useEAM |= (*i)->isEAM();
502 <        useCharge |= (*i)->isCharge();
503 <        useDirectional |= (*i)->isDirectional();
504 <        useDipole |= (*i)->isDipole();
505 <        useGayBerne |= (*i)->isGayBerne();
506 <        useSticky |= (*i)->isSticky();
507 <        useShape |= (*i)->isShape();
508 <    }
750 >    // we need arrays to hold the counts and displacement vectors for
751 >    // all processors
752 >    vector<int> counts(nproc, 0);
753 >    vector<int> disps(nproc, 0);
754  
755 <    if (useSticky || useDipole || useGayBerne || useShape) {
756 <        useDirectionalAtom = 1;
755 >    // fill the counts array
756 >    MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0],
757 >                              1, MPI::INT);
758 >  
759 >    // use the processor counts to compute the displacement array
760 >    disps[0] = 0;    
761 >    int totalCount = counts[0];
762 >    for (int iproc = 1; iproc < nproc; iproc++) {
763 >      disps[iproc] = disps[iproc-1] + counts[iproc-1];
764 >      totalCount += counts[iproc];
765      }
766  
767 <    if (useCharge || useDipole) {
768 <        useElectrostatics = 1;
769 <    }
767 >    // we need a (possibly redundant) set of all found types:
768 >    vector<int> ftGlobal(totalCount);
769 >    
770 >    // now spray out the foundTypes to all the other processors:    
771 >    MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT,
772 >                               &ftGlobal[0], &counts[0], &disps[0],
773 >                               MPI::INT);
774  
775 < #ifdef IS_MPI    
519 <    int temp;
775 >    vector<int>::iterator j;
776  
777 <    temp = usePBC;
778 <    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
777 >    // foundIdents is a stl set, so inserting an already found ident
778 >    // will have no effect.
779 >    set<int> foundIdents;
780  
781 <    temp = useDirectionalAtom;
782 <    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
781 >    for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j)
782 >      foundIdents.insert((*j));
783 >    
784 >    // now iterate over the foundIdents and get the actual atom types
785 >    // that correspond to these:
786 >    set<int>::iterator it;
787 >    for (it = foundIdents.begin(); it != foundIdents.end(); ++it)
788 >      atomTypes.insert( forceField_->getAtomType((*it)) );
789 >
790 > #endif
791  
792 <    temp = useLennardJones;
793 <    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
792 >    return atomTypes;        
793 >  }
794  
530    temp = useElectrostatics;
531    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
795  
796 <    temp = useCharge;
797 <    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
796 >  int getGlobalCountOfType(AtomType* atype) {
797 >    /*
798 >    set<AtomType*> atypes = getSimulatedAtomTypes();
799 >    map<AtomType*, int> counts_;
800  
801 <    temp = useDipole;
802 <    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
801 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
802 >      for(atom = mol->beginAtom(ai); atom != NULL;
803 >          atom = mol->nextAtom(ai)) {
804 >        atom->getAtomType();
805 >      }      
806 >    }    
807 >    */
808 >    return 0;
809 >  }
810  
811 <    temp = useSticky;
812 <    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
811 >  void SimInfo::setupSimVariables() {
812 >    useAtomicVirial_ = simParams_->getUseAtomicVirial();
813 >    // we only call setAccumulateBoxDipole if the accumulateBoxDipole
814 >    // parameter is true
815 >    calcBoxDipole_ = false;
816 >    if ( simParams_->haveAccumulateBoxDipole() )
817 >      if ( simParams_->getAccumulateBoxDipole() ) {
818 >        calcBoxDipole_ = true;      
819 >      }
820 >    
821 >    set<AtomType*>::iterator i;
822 >    set<AtomType*> atomTypes;
823 >    atomTypes = getSimulatedAtomTypes();    
824 >    bool usesElectrostatic = false;
825 >    bool usesMetallic = false;
826 >    bool usesDirectional = false;
827 >    bool usesFluctuatingCharges =  false;
828 >    //loop over all of the atom types
829 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
830 >      usesElectrostatic |= (*i)->isElectrostatic();
831 >      usesMetallic |= (*i)->isMetal();
832 >      usesDirectional |= (*i)->isDirectional();
833 >      usesFluctuatingCharges |= (*i)->isFluctuatingCharge();
834 >    }
835  
836 <    temp = useGayBerne;
837 <    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
836 > #ifdef IS_MPI
837 >    bool temp;
838 >    temp = usesDirectional;
839 >    MPI::COMM_WORLD.Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI::BOOL,
840 >                              MPI::LOR);
841 >        
842 >    temp = usesMetallic;
843 >    MPI::COMM_WORLD.Allreduce(&temp, &usesMetallicAtoms_, 1, MPI::BOOL,
844 >                              MPI::LOR);
845 >    
846 >    temp = usesElectrostatic;
847 >    MPI::COMM_WORLD.Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI::BOOL,
848 >                              MPI::LOR);
849  
850 <    temp = useEAM;
851 <    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
850 >    temp = usesFluctuatingCharges;
851 >    MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL,
852 >                              MPI::LOR);
853 > #else
854  
855 <    temp = useShape;
856 <    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
855 >    usesDirectionalAtoms_ = usesDirectional;
856 >    usesMetallicAtoms_ = usesMetallic;
857 >    usesElectrostaticAtoms_ = usesElectrostatic;
858 >    usesFluctuatingCharges_ = usesFluctuatingCharges;
859  
551    temp = useFLARB;
552    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
553
554    temp = useRF;
555    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
556    
860   #endif
861 +    
862 +    requiresPrepair_ = usesMetallicAtoms_ ? true : false;
863 +    requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false;
864 +    requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false;    
865 +  }
866  
559    fInfo_.SIM_uses_PBC = usePBC;    
560    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
561    fInfo_.SIM_uses_LennardJones = useLennardJones;
562    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
563    fInfo_.SIM_uses_Charges = useCharge;
564    fInfo_.SIM_uses_Dipoles = useDipole;
565    fInfo_.SIM_uses_Sticky = useSticky;
566    fInfo_.SIM_uses_GayBerne = useGayBerne;
567    fInfo_.SIM_uses_EAM = useEAM;
568    fInfo_.SIM_uses_Shapes = useShape;
569    fInfo_.SIM_uses_FLARB = useFLARB;
570    fInfo_.SIM_uses_RF = useRF;
867  
868 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
868 >  vector<int> SimInfo::getGlobalAtomIndices() {
869 >    SimInfo::MoleculeIterator mi;
870 >    Molecule* mol;
871 >    Molecule::AtomIterator ai;
872 >    Atom* atom;
873  
874 <        if (simParams_->haveDielectric()) {
875 <            fInfo_.dielect = simParams_->getDielectric();
876 <        } else {
877 <            sprintf(painCave.errMsg,
878 <                    "SimSetup Error: No Dielectric constant was set.\n"
879 <                    "\tYou are trying to use Reaction Field without"
880 <                    "\tsetting a dielectric constant!\n");
581 <            painCave.isFatal = 1;
582 <            simError();
583 <        }
584 <        
585 <    } else {
586 <        fInfo_.dielect = 0.0;
874 >    vector<int> GlobalAtomIndices(getNAtoms(), 0);
875 >    
876 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
877 >      
878 >      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
879 >        GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex();
880 >      }
881      }
882 +    return GlobalAtomIndices;
883 +  }
884  
589 }
885  
886 < void SimInfo::setupFortranSim() {
887 <    int isError;
888 <    int nExclude;
889 <    std::vector<int> fortranGlobalGroupMembership;
890 <    
596 <    nExclude = exclude_.getSize();
597 <    isError = 0;
886 >  vector<int> SimInfo::getGlobalGroupIndices() {
887 >    SimInfo::MoleculeIterator mi;
888 >    Molecule* mol;
889 >    Molecule::CutoffGroupIterator ci;
890 >    CutoffGroup* cg;
891  
892 <    //globalGroupMembership_ is filled by SimCreator    
893 <    for (int i = 0; i < nGlobalAtoms_; i++) {
894 <        fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
892 >    vector<int> GlobalGroupIndices;
893 >    
894 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
895 >      
896 >      //local index of cutoff group is trivial, it only depends on the
897 >      //order of travesing
898 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
899 >           cg = mol->nextCutoffGroup(ci)) {
900 >        GlobalGroupIndices.push_back(cg->getGlobalIndex());
901 >      }        
902      }
903 +    return GlobalGroupIndices;
904 +  }
905  
906 +
907 +  void SimInfo::prepareTopology() {
908 +
909      //calculate mass ratio of cutoff group
605    std::vector<double> mfact;
910      SimInfo::MoleculeIterator mi;
911      Molecule* mol;
912      Molecule::CutoffGroupIterator ci;
913      CutoffGroup* cg;
914      Molecule::AtomIterator ai;
915      Atom* atom;
916 <    double totalMass;
916 >    RealType totalMass;
917  
918 <    //to avoid memory reallocation, reserve enough space for mfact
919 <    mfact.reserve(getNCutoffGroups());
918 >    /**
919 >     * The mass factor is the relative mass of an atom to the total
920 >     * mass of the cutoff group it belongs to.  By default, all atoms
921 >     * are their own cutoff groups, and therefore have mass factors of
922 >     * 1.  We need some special handling for massless atoms, which
923 >     * will be treated as carrying the entire mass of the cutoff
924 >     * group.
925 >     */
926 >    massFactors_.clear();
927 >    massFactors_.resize(getNAtoms(), 1.0);
928      
929      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
930 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
930 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
931 >           cg = mol->nextCutoffGroup(ci)) {
932  
933 <            totalMass = cg->getMass();
934 <            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
935 <                        mfact.push_back(atom->getMass()/totalMass);
936 <            }
937 <
938 <        }      
933 >        totalMass = cg->getMass();
934 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
935 >          // Check for massless groups - set mfact to 1 if true
936 >          if (totalMass != 0)
937 >            massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass;
938 >          else
939 >            massFactors_[atom->getLocalIndex()] = 1.0;
940 >        }
941 >      }      
942      }
943  
944 <    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
629 <    std::vector<int> identArray;
944 >    // Build the identArray_ and regions_
945  
946 <    //to avoid memory reallocation, reserve enough space identArray
947 <    identArray.reserve(getNAtoms());
948 <    
949 <    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
950 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
951 <            identArray.push_back(atom->getIdent());
952 <        }
946 >    identArray_.clear();
947 >    identArray_.reserve(getNAtoms());  
948 >    regions_.clear();
949 >    regions_.reserve(getNAtoms());
950 >
951 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {      
952 >      int reg = mol->getRegion();      
953 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
954 >        identArray_.push_back(atom->getIdent());
955 >        regions_.push_back(reg);
956 >      }
957      }    
958 +      
959 +    topologyDone_ = true;
960 +  }
961  
962 <    //fill molMembershipArray
641 <    //molMembershipArray is filled by SimCreator    
642 <    std::vector<int> molMembershipArray(nGlobalAtoms_);
643 <    for (int i = 0; i < nGlobalAtoms_; i++) {
644 <        molMembershipArray[i] = globalMolMembership_[i] + 1;
645 <    }
646 <    
647 <    //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
651 <    int nGlobalExcludes = 0;
652 <    int* globalExcludes = NULL;
653 <    int* excludeList = exclude_.getExcludeList();
654 <    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
655 <                  &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
656 <                  &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
657 <
658 <    if( isError ){
659 <
660 <        sprintf( painCave.errMsg,
661 <                 "There was an error setting the simulation information in fortran.\n" );
662 <        painCave.isFatal = 1;
663 <        painCave.severity = OOPSE_ERROR;
664 <        simError();
665 <    }
666 <
667 < #ifdef IS_MPI
668 <    sprintf( checkPointMsg,
669 <       "succesfully sent the simulation information to fortran.\n");
670 <    MPIcheckPoint();
671 < #endif // is_mpi
672 < }
673 <
674 <
675 < #ifdef IS_MPI
676 < void SimInfo::setupFortranParallel() {
677 <    
678 <    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
679 <    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
680 <    std::vector<int> localToGlobalCutoffGroupIndex;
681 <    SimInfo::MoleculeIterator mi;
682 <    Molecule::AtomIterator ai;
683 <    Molecule::CutoffGroupIterator ci;
684 <    Molecule* mol;
685 <    Atom* atom;
686 <    CutoffGroup* cg;
687 <    mpiSimData parallelData;
688 <    int isError;
689 <
690 <    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
691 <
692 <        //local index(index in DataStorge) of atom is important
693 <        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
694 <            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
695 <        }
696 <
697 <        //local index of cutoff group is trivial, it only depends on the order of travesing
698 <        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
699 <            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
700 <        }        
701 <        
702 <    }
703 <
704 <    //fill up mpiSimData struct
705 <    parallelData.nMolGlobal = getNGlobalMolecules();
706 <    parallelData.nMolLocal = getNMolecules();
707 <    parallelData.nAtomsGlobal = getNGlobalAtoms();
708 <    parallelData.nAtomsLocal = getNAtoms();
709 <    parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
710 <    parallelData.nGroupsLocal = getNCutoffGroups();
711 <    parallelData.myNode = worldRank;
712 <    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors));
713 <
714 <    //pass mpiSimData struct and index arrays to fortran
715 <    setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
716 <                    &localToGlobalAtomIndex[0],  &(parallelData.nGroupsLocal),
717 <                    &localToGlobalCutoffGroupIndex[0], &isError);
718 <
719 <    if (isError) {
720 <        sprintf(painCave.errMsg,
721 <                "mpiRefresh errror: fortran didn't like something we gave it.\n");
722 <        painCave.isFatal = 1;
723 <        simError();
724 <    }
725 <
726 <    sprintf(checkPointMsg, " mpiRefresh successful.\n");
727 <    MPIcheckPoint();
728 <
729 <
730 < }
731 <
732 < #endif
733 <
734 < double SimInfo::calcMaxCutoffRadius() {
735 <
736 <
737 <    std::set<AtomType*> atomTypes;
738 <    std::set<AtomType*>::iterator i;
739 <    std::vector<double> cutoffRadius;
740 <
741 <    //get the unique atom types
742 <    atomTypes = getUniqueAtomTypes();
743 <
744 <    //query the max cutoff radius among these atom types
745 <    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
746 <        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
747 <    }
748 <
749 <    double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
750 < #ifdef IS_MPI
751 <    //pick the max cutoff radius among the processors
752 < #endif
753 <
754 <    return maxCutoffRadius;
755 < }
756 <
757 < void SimInfo::getCutoff(double& rcut, double& rsw) {
758 <    
759 <    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
760 <        
761 <        if (!simParams_->haveRcut()){
762 <            sprintf(painCave.errMsg,
763 <                "SimCreator Warning: No value was set for the cutoffRadius.\n"
764 <                "\tOOPSE will use a default value of 15.0 angstroms"
765 <                "\tfor the cutoffRadius.\n");
766 <            painCave.isFatal = 0;
767 <            simError();
768 <            rcut = 15.0;
769 <        } else{
770 <            rcut = simParams_->getRcut();
771 <        }
772 <
773 <        if (!simParams_->haveRsw()){
774 <            sprintf(painCave.errMsg,
775 <                "SimCreator Warning: No value was set for switchingRadius.\n"
776 <                "\tOOPSE will use a default value of\n"
777 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
778 <            painCave.isFatal = 0;
779 <            simError();
780 <            rsw = 0.95 * rcut;
781 <        } else{
782 <            rsw = simParams_->getRsw();
783 <        }
784 <
785 <    } else {
786 <        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
787 <        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
788 <        
789 <        if (simParams_->haveRcut()) {
790 <            rcut = simParams_->getRcut();
791 <        } else {
792 <            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
793 <            rcut = calcMaxCutoffRadius();
794 <        }
795 <
796 <        if (simParams_->haveRsw()) {
797 <            rsw  = simParams_->getRsw();
798 <        } else {
799 <            rsw = rcut;
800 <        }
801 <    
802 <    }
803 < }
804 <
805 < void SimInfo::setupCutoff() {
806 <    getCutoff(rcut_, rsw_);    
807 <    double rnblist = rcut_ + 1; // skin of neighbor list
808 <
809 <    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
810 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
811 < }
812 <
813 < void SimInfo::addProperty(GenericData* genData) {
962 >  void SimInfo::addProperty(GenericData* genData) {
963      properties_.addProperty(genData);  
964 < }
964 >  }
965  
966 < void SimInfo::removeProperty(const std::string& propName) {
966 >  void SimInfo::removeProperty(const string& propName) {
967      properties_.removeProperty(propName);  
968 < }
968 >  }
969  
970 < void SimInfo::clearProperties() {
970 >  void SimInfo::clearProperties() {
971      properties_.clearProperties();
972 < }
972 >  }
973  
974 < std::vector<std::string> SimInfo::getPropertyNames() {
974 >  vector<string> SimInfo::getPropertyNames() {
975      return properties_.getPropertyNames();  
976 < }
976 >  }
977        
978 < std::vector<GenericData*> SimInfo::getProperties() {
978 >  vector<GenericData*> SimInfo::getProperties() {
979      return properties_.getProperties();
980 < }
980 >  }
981  
982 < GenericData* SimInfo::getPropertyByName(const std::string& propName) {
982 >  GenericData* SimInfo::getPropertyByName(const string& propName) {
983      return properties_.getPropertyByName(propName);
984 < }
984 >  }
985  
986 < void SimInfo::setSnapshotManager(SnapshotManager* sman) {
987 <    if (sman_ == sman_) {
988 <        return;
989 <    }
841 <    
986 >  void SimInfo::setSnapshotManager(SnapshotManager* sman) {
987 >    if (sman_ == sman) {
988 >      return;
989 >    }    
990      delete sman_;
991      sman_ = sman;
992  
993      Molecule* mol;
994      RigidBody* rb;
995      Atom* atom;
996 +    CutoffGroup* cg;
997      SimInfo::MoleculeIterator mi;
998      Molecule::RigidBodyIterator rbIter;
999 <    Molecule::AtomIterator atomIter;;
999 >    Molecule::AtomIterator atomIter;
1000 >    Molecule::CutoffGroupIterator cgIter;
1001  
1002      for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
1003          
1004 <        for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
1005 <            atom->setSnapshotManager(sman_);
1006 <        }
1004 >      for (atom = mol->beginAtom(atomIter); atom != NULL;
1005 >           atom = mol->nextAtom(atomIter)) {
1006 >        atom->setSnapshotManager(sman_);
1007 >      }
1008          
1009 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
1010 <            rb->setSnapshotManager(sman_);
1011 <        }
1009 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
1010 >           rb = mol->nextRigidBody(rbIter)) {
1011 >        rb->setSnapshotManager(sman_);
1012 >      }
1013 >
1014 >      for (cg = mol->beginCutoffGroup(cgIter); cg != NULL;
1015 >           cg = mol->nextCutoffGroup(cgIter)) {
1016 >        cg->setSnapshotManager(sman_);
1017 >      }
1018      }    
1019      
1020 < }
1020 >  }
1021  
865 Vector3d SimInfo::getComVel(){
866    SimInfo::MoleculeIterator i;
867    Molecule* mol;
1022  
1023 <    Vector3d comVel(0.0);
870 <    double totalMass = 0.0;
871 <    
872 <
873 <    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
874 <        double mass = mol->getMass();
875 <        totalMass += mass;
876 <        comVel += mass * mol->getComVel();
877 <    }  
1023 >  ostream& operator <<(ostream& o, SimInfo& info) {
1024  
1025 < #ifdef IS_MPI
1026 <    double tmpMass = totalMass;
1027 <    Vector3d tmpComVel(comVel);    
1028 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1029 <    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1030 < #endif
1025 >    return o;
1026 >  }
1027 >  
1028 >  
1029 >  StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
1030 >    if (index >= int(IOIndexToIntegrableObject.size())) {
1031 >      sprintf(painCave.errMsg,
1032 >              "SimInfo::getIOIndexToIntegrableObject Error: Integrable Object\n"
1033 >              "\tindex exceeds number of known objects!\n");
1034 >      painCave.isFatal = 1;
1035 >      simError();
1036 >      return NULL;
1037 >    } else
1038 >      return IOIndexToIntegrableObject.at(index);
1039 >  }
1040 >  
1041 >  void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) {
1042 >    IOIndexToIntegrableObject= v;
1043 >  }
1044  
1045 <    comVel /= totalMass;
1046 <
888 <    return comVel;
889 < }
890 <
891 < Vector3d SimInfo::getCom(){
892 <    SimInfo::MoleculeIterator i;
893 <    Molecule* mol;
894 <
895 <    Vector3d com(0.0);
896 <    double totalMass = 0.0;
897 <    
898 <    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
899 <        double mass = mol->getMass();
900 <        totalMass += mass;
901 <        com += mass * mol->getCom();
902 <    }  
903 <
1045 >  int SimInfo::getNGlobalConstraints() {
1046 >    int nGlobalConstraints;
1047   #ifdef IS_MPI
1048 <    double tmpMass = totalMass;
1049 <    Vector3d tmpCom(com);    
1050 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1051 <    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1048 >    MPI::COMM_WORLD.Allreduce(&nConstraints_, &nGlobalConstraints, 1,
1049 >                              MPI::INT, MPI::SUM);
1050 > #else
1051 >    nGlobalConstraints =  nConstraints_;
1052   #endif
1053 +    return nGlobalConstraints;
1054 +  }
1055  
1056 <    com /= totalMass;
1056 > }//end namespace OpenMD
1057  
913    return com;
914
915 }        
916
917 std::ostream& operator <<(std::ostream& o, SimInfo& info) {
918
919    return o;
920 }
921
922 }//end namespace oopse
923

Comparing trunk/src/brains/SimInfo.cpp (property svn:keywords):
Revision 331 by tim, Sun Feb 13 21:18:27 2005 UTC vs.
Revision 1938 by gezelter, Thu Oct 31 15:32:17 2013 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines