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

Comparing:
trunk/src/brains/SimInfo.cpp (file contents), Revision 124 by chuckv, Wed Oct 20 20:46:20 2004 UTC vs.
branches/development/src/brains/SimInfo.cpp (file contents), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

# Line 1 | Line 1
1 < #include <stdlib.h>
2 < #include <string.h>
3 < #include <math.h>
1 > /*
2 > * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 > *
4 > * The University of Notre Dame grants you ("Licensee") a
5 > * non-exclusive, royalty free, license to use, modify and
6 > * redistribute this software in source and binary code form, provided
7 > * that the following conditions are met:
8 > *
9 > * 1. Redistributions of source code must retain the above copyright
10 > *    notice, this list of conditions and the following disclaimer.
11 > *
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.
16 > *
17 > * This software is provided "AS IS," without a warranty of any
18 > * kind. All express or implied conditions, representations and
19 > * warranties, including any implied warranty of merchantability,
20 > * fitness for a particular purpose or non-infringement, are hereby
21 > * excluded.  The University of Notre Dame and its licensors shall not
22 > * be liable for any damages suffered by licensee as a result of
23 > * using, modifying or distributing the software or its
24 > * derivatives. In no event will the University of Notre Dame or its
25 > * licensors be liable for any lost revenue, profit or data, or for
26 > * direct, indirect, special, consequential, incidental or punitive
27 > * damages, however caused and regardless of the theory of liability,
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, 24107 (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 > /**
44 > * @file SimInfo.cpp
45 > * @author    tlin
46 > * @date  11/02/2004
47 > * @version 1.0
48 > */
49  
50 < #include <iostream>
51 < using namespace std;
50 > #include <algorithm>
51 > #include <set>
52 > #include <map>
53  
54   #include "brains/SimInfo.hpp"
55 < #define __C
56 < #include "brains/fSimulation.h"
55 > #include "math/Vector3.hpp"
56 > #include "primitives/Molecule.hpp"
57 > #include "primitives/StuntDouble.hpp"
58 > #include "utils/MemoryUtils.hpp"
59   #include "utils/simError.h"
60 < #include "UseTheForce/DarkSide/simulation_interface.h"
61 < #include "UseTheForce/notifyCutoffs_interface.h"
62 <
63 < //#include "UseTheForce/fortranWrappers.hpp"
16 <
17 < #include "math/MatVec3.h"
18 <
60 > #include "selection/SelectionManager.hpp"
61 > #include "io/ForceFieldOptions.hpp"
62 > #include "UseTheForce/ForceField.hpp"
63 > #include "nonbonded/SwitchingFunction.hpp"
64   #ifdef IS_MPI
65 < #include "brains/mpiSimulation.hpp"
65 > #include <mpi.h>
66   #endif
67  
68 < inline double roundMe( double x ){
69 <  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
70 < }
71 <          
72 < inline double min( double a, double b ){
73 <  return (a < b ) ? a : b;
74 < }
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),
76 >    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nInversions_(0),
77 >    nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0),
78 >    nConstraints_(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
85 >    CutoffGroupStamp* cgStamp;    
86 >    RigidBodyStamp* rbStamp;
87 >    int nRigidAtoms = 0;
88 >    
89 >    vector<Component*> components = simParams->getComponents();
90 >    
91 >    for (vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
92 >      molStamp = (*i)->getMoleculeStamp();
93 >      nMolWithSameStamp = (*i)->getNMol();
94 >      
95 >      addMoleculeStamp(molStamp, nMolWithSameStamp);
96 >      
97 >      //calculate atoms in molecules
98 >      nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
99 >      
100 >      //calculate atoms in cutoff groups
101 >      int nAtomsInGroups = 0;
102 >      int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
103 >      
104 >      for (int j=0; j < nCutoffGroupsInStamp; j++) {
105 >        cgStamp = molStamp->getCutoffGroupStamp(j);
106 >        nAtomsInGroups += cgStamp->getNMembers();
107 >      }
108 >      
109 >      nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
110 >      
111 >      nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
112 >      
113 >      //calculate atoms in rigid bodies
114 >      int nAtomsInRigidBodies = 0;
115 >      int nRigidBodiesInStamp = molStamp->getNRigidBodies();
116 >      
117 >      for (int j=0; j < nRigidBodiesInStamp; j++) {
118 >        rbStamp = molStamp->getRigidBodyStamp(j);
119 >        nAtomsInRigidBodies += rbStamp->getNMembers();
120 >      }
121 >      
122 >      nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
123 >      nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
124 >      
125 >    }
126 >    
127 >    //every free atom (atom does not belong to cutoff groups) is a cutoff
128 >    //group therefore the total number of cutoff groups in the system is
129 >    //equal to the total number of atoms minus number of atoms belong to
130 >    //cutoff group defined in meta-data file plus the number of cutoff
131 >    //groups defined in meta-data file
132  
133 < SimInfo* currentInfo;
133 >    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
134 >    
135 >    //every free atom (atom does not belong to rigid bodies) is an
136 >    //integrable object therefore the total number of integrable objects
137 >    //in the system is equal to the total number of atoms minus number of
138 >    //atoms belong to rigid body defined in meta-data file plus the number
139 >    //of rigid bodies defined in meta-data file
140 >    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
141 >      + nGlobalRigidBodies_;
142 >    
143 >    nGlobalMols_ = molStampIds_.size();
144 >    molToProcMap_.resize(nGlobalMols_);
145 >  }
146 >  
147 >  SimInfo::~SimInfo() {
148 >    map<int, Molecule*>::iterator i;
149 >    for (i = molecules_.begin(); i != molecules_.end(); ++i) {
150 >      delete i->second;
151 >    }
152 >    molecules_.clear();
153 >      
154 >    delete sman_;
155 >    delete simParams_;
156 >    delete forceField_;
157 >  }
158  
33 SimInfo::SimInfo(){
159  
160 <  n_constraints = 0;
161 <  nZconstraints = 0;
162 <  n_oriented = 0;
163 <  n_dipoles = 0;
164 <  ndf = 0;
165 <  ndfRaw = 0;
166 <  nZconstraints = 0;
167 <  the_integrator = NULL;
168 <  setTemp = 0;
169 <  thermalTime = 0.0;
170 <  currentTime = 0.0;
171 <  rCut = 0.0;
172 <  rSw = 0.0;
173 <
174 <  haveRcut = 0;
175 <  haveRsw = 0;
176 <  boxIsInit = 0;
160 >  bool SimInfo::addMolecule(Molecule* mol) {
161 >    MoleculeIterator i;
162 >    
163 >    i = molecules_.find(mol->getGlobalIndex());
164 >    if (i == molecules_.end() ) {
165 >      
166 >      molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
167 >      
168 >      nAtoms_ += mol->getNAtoms();
169 >      nBonds_ += mol->getNBonds();
170 >      nBends_ += mol->getNBends();
171 >      nTorsions_ += mol->getNTorsions();
172 >      nInversions_ += mol->getNInversions();
173 >      nRigidBodies_ += mol->getNRigidBodies();
174 >      nIntegrableObjects_ += mol->getNIntegrableObjects();
175 >      nCutoffGroups_ += mol->getNCutoffGroups();
176 >      nConstraints_ += mol->getNConstraintPairs();
177 >      
178 >      addInteractionPairs(mol);
179 >      
180 >      return true;
181 >    } else {
182 >      return false;
183 >    }
184 >  }
185    
186 <  resetTime = 1e99;
186 >  bool SimInfo::removeMolecule(Molecule* mol) {
187 >    MoleculeIterator i;
188 >    i = molecules_.find(mol->getGlobalIndex());
189  
190 <  orthoRhombic = 0;
56 <  orthoTolerance = 1E-6;
57 <  useInitXSstate = true;
190 >    if (i != molecules_.end() ) {
191  
192 <  usePBC = 0;
193 <  useLJ = 0;
194 <  useSticky = 0;
195 <  useCharges = 0;
196 <  useDipoles = 0;
197 <  useReactionField = 0;
198 <  useGB = 0;
199 <  useEAM = 0;
200 <  useSolidThermInt = 0;
201 <  useLiquidThermInt = 0;
202 <
70 <  haveCutoffGroups = false;
71 <
72 <  excludes = Exclude::Instance();
73 <
74 <  myConfiguration = new SimState();
75 <
76 <  has_minimizer = false;
77 <  the_minimizer =NULL;
78 <
79 <  ngroup = 0;
80 <
81 < }
192 >      assert(mol == i->second);
193 >        
194 >      nAtoms_ -= mol->getNAtoms();
195 >      nBonds_ -= mol->getNBonds();
196 >      nBends_ -= mol->getNBends();
197 >      nTorsions_ -= mol->getNTorsions();
198 >      nInversions_ -= mol->getNInversions();
199 >      nRigidBodies_ -= mol->getNRigidBodies();
200 >      nIntegrableObjects_ -= mol->getNIntegrableObjects();
201 >      nCutoffGroups_ -= mol->getNCutoffGroups();
202 >      nConstraints_ -= mol->getNConstraintPairs();
203  
204 +      removeInteractionPairs(mol);
205 +      molecules_.erase(mol->getGlobalIndex());
206  
207 < SimInfo::~SimInfo(){
207 >      delete mol;
208 >        
209 >      return true;
210 >    } else {
211 >      return false;
212 >    }
213 >  }    
214  
215 <  delete myConfiguration;
215 >        
216 >  Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
217 >    i = molecules_.begin();
218 >    return i == molecules_.end() ? NULL : i->second;
219 >  }    
220  
221 <  map<string, GenericData*>::iterator i;
222 <  
223 <  for(i = properties.begin(); i != properties.end(); i++)
224 <    delete (*i).second;
221 >  Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
222 >    ++i;
223 >    return i == molecules_.end() ? NULL : i->second;    
224 >  }
225  
93 }
226  
227 < void SimInfo::setBox(double newBox[3]) {
228 <  
229 <  int i, j;
230 <  double tempMat[3][3];
227 >  void SimInfo::calcNdf() {
228 >    int ndf_local;
229 >    MoleculeIterator i;
230 >    vector<StuntDouble*>::iterator j;
231 >    Molecule* mol;
232 >    StuntDouble* integrableObject;
233  
234 <  for(i=0; i<3; i++)
235 <    for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
234 >    ndf_local = 0;
235 >    
236 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
237 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
238 >           integrableObject = mol->nextIntegrableObject(j)) {
239  
240 <  tempMat[0][0] = newBox[0];
104 <  tempMat[1][1] = newBox[1];
105 <  tempMat[2][2] = newBox[2];
240 >        ndf_local += 3;
241  
242 <  setBoxM( tempMat );
242 >        if (integrableObject->isDirectional()) {
243 >          if (integrableObject->isLinear()) {
244 >            ndf_local += 2;
245 >          } else {
246 >            ndf_local += 3;
247 >          }
248 >        }
249 >            
250 >      }
251 >    }
252 >    
253 >    // n_constraints is local, so subtract them on each processor
254 >    ndf_local -= nConstraints_;
255  
256 < }
256 > #ifdef IS_MPI
257 >    MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
258 > #else
259 >    ndf_ = ndf_local;
260 > #endif
261  
262 < void SimInfo::setBoxM( double theBox[3][3] ){
263 <  
264 <  int i, j;
114 <  double FortranHmat[9]; // to preserve compatibility with Fortran the
115 <                         // ordering in the array is as follows:
116 <                         // [ 0 3 6 ]
117 <                         // [ 1 4 7 ]
118 <                         // [ 2 5 8 ]
119 <  double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
262 >    // nZconstraints_ is global, as are the 3 COM translations for the
263 >    // entire system:
264 >    ndf_ = ndf_ - 3 - nZconstraint_;
265  
266 <  if( !boxIsInit ) boxIsInit = 1;
266 >  }
267  
268 <  for(i=0; i < 3; i++)
269 <    for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
268 >  int SimInfo::getFdf() {
269 > #ifdef IS_MPI
270 >    MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
271 > #else
272 >    fdf_ = fdf_local;
273 > #endif
274 >    return fdf_;
275 >  }
276    
277 <  calcBoxL();
278 <  calcHmatInv();
279 <
280 <  for(i=0; i < 3; i++) {
281 <    for (j=0; j < 3; j++) {
282 <      FortranHmat[3*j + i] = Hmat[i][j];
283 <      FortranHmatInv[3*j + i] = HmatInv[i][j];
277 >  unsigned int SimInfo::getNLocalCutoffGroups(){
278 >    int nLocalCutoffAtoms = 0;
279 >    Molecule* mol;
280 >    MoleculeIterator mi;
281 >    CutoffGroup* cg;
282 >    Molecule::CutoffGroupIterator ci;
283 >    
284 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
285 >      
286 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
287 >           cg = mol->nextCutoffGroup(ci)) {
288 >        nLocalCutoffAtoms += cg->getNumAtom();
289 >        
290 >      }        
291      }
292 +    
293 +    return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_;
294    }
295 +    
296 +  void SimInfo::calcNdfRaw() {
297 +    int ndfRaw_local;
298  
299 <  setFortranBox(FortranHmat, FortranHmatInv, &orthoRhombic);
300 <
301 < }
302 <
299 >    MoleculeIterator i;
300 >    vector<StuntDouble*>::iterator j;
301 >    Molecule* mol;
302 >    StuntDouble* integrableObject;
303  
304 < void SimInfo::getBoxM (double theBox[3][3]) {
304 >    // Raw degrees of freedom that we have to set
305 >    ndfRaw_local = 0;
306 >    
307 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
308 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
309 >           integrableObject = mol->nextIntegrableObject(j)) {
310  
311 <  int i, j;
144 <  for(i=0; i<3; i++)
145 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
146 < }
311 >        ndfRaw_local += 3;
312  
313 +        if (integrableObject->isDirectional()) {
314 +          if (integrableObject->isLinear()) {
315 +            ndfRaw_local += 2;
316 +          } else {
317 +            ndfRaw_local += 3;
318 +          }
319 +        }
320 +            
321 +      }
322 +    }
323 +    
324 + #ifdef IS_MPI
325 +    MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
326 + #else
327 +    ndfRaw_ = ndfRaw_local;
328 + #endif
329 +  }
330  
331 < void SimInfo::scaleBox(double scale) {
332 <  double theBox[3][3];
151 <  int i, j;
331 >  void SimInfo::calcNdfTrans() {
332 >    int ndfTrans_local;
333  
334 <  // cerr << "Scaling box by " << scale << "\n";
334 >    ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
335  
155  for(i=0; i<3; i++)
156    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
336  
337 <  setBoxM(theBox);
337 > #ifdef IS_MPI
338 >    MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
339 > #else
340 >    ndfTrans_ = ndfTrans_local;
341 > #endif
342  
343 < }
343 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
344 >
345 >  }
346  
347 < void SimInfo::calcHmatInv( void ) {
348 <  
349 <  int oldOrtho;
350 <  int i,j;
351 <  double smallDiag;
352 <  double tol;
353 <  double sanity[3][3];
347 >  void SimInfo::addInteractionPairs(Molecule* mol) {
348 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
349 >    vector<Bond*>::iterator bondIter;
350 >    vector<Bend*>::iterator bendIter;
351 >    vector<Torsion*>::iterator torsionIter;
352 >    vector<Inversion*>::iterator inversionIter;
353 >    Bond* bond;
354 >    Bend* bend;
355 >    Torsion* torsion;
356 >    Inversion* inversion;
357 >    int a;
358 >    int b;
359 >    int c;
360 >    int d;
361  
362 <  invertMat3( Hmat, HmatInv );
362 >    // atomGroups can be used to add special interaction maps between
363 >    // groups of atoms that are in two separate rigid bodies.
364 >    // However, most site-site interactions between two rigid bodies
365 >    // are probably not special, just the ones between the physically
366 >    // bonded atoms.  Interactions *within* a single rigid body should
367 >    // always be excluded.  These are done at the bottom of this
368 >    // function.
369  
370 <  // check to see if Hmat is orthorhombic
371 <  
372 <  oldOrtho = orthoRhombic;
373 <
374 <  smallDiag = fabs(Hmat[0][0]);
375 <  if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
376 <  if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
377 <  tol = smallDiag * orthoTolerance;
378 <
379 <  orthoRhombic = 1;
380 <  
381 <  for (i = 0; i < 3; i++ ) {
382 <    for (j = 0 ; j < 3; j++) {
383 <      if (i != j) {
384 <        if (orthoRhombic) {
385 <          if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
386 <        }        
370 >    map<int, set<int> > atomGroups;
371 >    Molecule::RigidBodyIterator rbIter;
372 >    RigidBody* rb;
373 >    Molecule::IntegrableObjectIterator ii;
374 >    StuntDouble* integrableObject;
375 >    
376 >    for (integrableObject = mol->beginIntegrableObject(ii);
377 >         integrableObject != NULL;
378 >         integrableObject = mol->nextIntegrableObject(ii)) {
379 >      
380 >      if (integrableObject->isRigidBody()) {
381 >        rb = static_cast<RigidBody*>(integrableObject);
382 >        vector<Atom*> atoms = rb->getAtoms();
383 >        set<int> rigidAtoms;
384 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
385 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
386 >        }
387 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
388 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
389 >        }      
390 >      } else {
391 >        set<int> oneAtomSet;
392 >        oneAtomSet.insert(integrableObject->getGlobalIndex());
393 >        atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
394        }
395 <    }
396 <  }
395 >    }  
396 >          
397 >    for (bond= mol->beginBond(bondIter); bond != NULL;
398 >         bond = mol->nextBond(bondIter)) {
399  
400 <  if( oldOrtho != orthoRhombic ){
400 >      a = bond->getAtomA()->getGlobalIndex();
401 >      b = bond->getAtomB()->getGlobalIndex();  
402      
403 <    if( orthoRhombic ) {
404 <      sprintf( painCave.errMsg,
405 <               "OOPSE is switching from the default Non-Orthorhombic\n"
406 <               "\tto the faster Orthorhombic periodic boundary computations.\n"
407 <               "\tThis is usually a good thing, but if you wan't the\n"
200 <               "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
201 <               "\tvariable ( currently set to %G ) smaller.\n",
202 <               orthoTolerance);
203 <      painCave.severity = OOPSE_INFO;
204 <      simError();
403 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
404 >        oneTwoInteractions_.addPair(a, b);
405 >      } else {
406 >        excludedInteractions_.addPair(a, b);
407 >      }
408      }
206    else {
207      sprintf( painCave.errMsg,
208               "OOPSE is switching from the faster Orthorhombic to the more\n"
209               "\tflexible Non-Orthorhombic periodic boundary computations.\n"
210               "\tThis is usually because the box has deformed under\n"
211               "\tNPTf integration. If you wan't to live on the edge with\n"
212               "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
213               "\tvariable ( currently set to %G ) larger.\n",
214               orthoTolerance);
215      painCave.severity = OOPSE_WARNING;
216      simError();
217    }
218  }
219 }
409  
410 < void SimInfo::calcBoxL( void ){
410 >    for (bend= mol->beginBend(bendIter); bend != NULL;
411 >         bend = mol->nextBend(bendIter)) {
412  
413 <  double dx, dy, dz, dsq;
413 >      a = bend->getAtomA()->getGlobalIndex();
414 >      b = bend->getAtomB()->getGlobalIndex();        
415 >      c = bend->getAtomC()->getGlobalIndex();
416 >      
417 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
418 >        oneTwoInteractions_.addPair(a, b);      
419 >        oneTwoInteractions_.addPair(b, c);
420 >      } else {
421 >        excludedInteractions_.addPair(a, b);
422 >        excludedInteractions_.addPair(b, c);
423 >      }
424  
425 <  // boxVol = Determinant of Hmat
425 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
426 >        oneThreeInteractions_.addPair(a, c);      
427 >      } else {
428 >        excludedInteractions_.addPair(a, c);
429 >      }
430 >    }
431  
432 <  boxVol = matDet3( Hmat );
432 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
433 >         torsion = mol->nextTorsion(torsionIter)) {
434  
435 <  // boxLx
436 <  
437 <  dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
438 <  dsq = dx*dx + dy*dy + dz*dz;
233 <  boxL[0] = sqrt( dsq );
234 <  //maxCutoff = 0.5 * boxL[0];
435 >      a = torsion->getAtomA()->getGlobalIndex();
436 >      b = torsion->getAtomB()->getGlobalIndex();        
437 >      c = torsion->getAtomC()->getGlobalIndex();        
438 >      d = torsion->getAtomD()->getGlobalIndex();      
439  
440 <  // boxLy
441 <  
442 <  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
443 <  dsq = dx*dx + dy*dy + dz*dz;
444 <  boxL[1] = sqrt( dsq );
445 <  //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
440 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
441 >        oneTwoInteractions_.addPair(a, b);      
442 >        oneTwoInteractions_.addPair(b, c);
443 >        oneTwoInteractions_.addPair(c, d);
444 >      } else {
445 >        excludedInteractions_.addPair(a, b);
446 >        excludedInteractions_.addPair(b, c);
447 >        excludedInteractions_.addPair(c, d);
448 >      }
449  
450 +      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
451 +        oneThreeInteractions_.addPair(a, c);      
452 +        oneThreeInteractions_.addPair(b, d);      
453 +      } else {
454 +        excludedInteractions_.addPair(a, c);
455 +        excludedInteractions_.addPair(b, d);
456 +      }
457  
458 <  // boxLz
459 <  
460 <  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
461 <  dsq = dx*dx + dy*dy + dz*dz;
462 <  boxL[2] = sqrt( dsq );
463 <  //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
458 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
459 >        oneFourInteractions_.addPair(a, d);      
460 >      } else {
461 >        excludedInteractions_.addPair(a, d);
462 >      }
463 >    }
464  
465 <  //calculate the max cutoff
466 <  maxCutoff =  calcMaxCutOff();
253 <  
254 <  checkCutOffs();
465 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
466 >         inversion = mol->nextInversion(inversionIter)) {
467  
468 < }
468 >      a = inversion->getAtomA()->getGlobalIndex();
469 >      b = inversion->getAtomB()->getGlobalIndex();        
470 >      c = inversion->getAtomC()->getGlobalIndex();        
471 >      d = inversion->getAtomD()->getGlobalIndex();        
472  
473 +      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
474 +        oneTwoInteractions_.addPair(a, b);      
475 +        oneTwoInteractions_.addPair(a, c);
476 +        oneTwoInteractions_.addPair(a, d);
477 +      } else {
478 +        excludedInteractions_.addPair(a, b);
479 +        excludedInteractions_.addPair(a, c);
480 +        excludedInteractions_.addPair(a, d);
481 +      }
482  
483 < double SimInfo::calcMaxCutOff(){
483 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
484 >        oneThreeInteractions_.addPair(b, c);    
485 >        oneThreeInteractions_.addPair(b, d);    
486 >        oneThreeInteractions_.addPair(c, d);      
487 >      } else {
488 >        excludedInteractions_.addPair(b, c);
489 >        excludedInteractions_.addPair(b, d);
490 >        excludedInteractions_.addPair(c, d);
491 >      }
492 >    }
493  
494 <  double ri[3], rj[3], rk[3];
495 <  double rij[3], rjk[3], rki[3];
496 <  double minDist;
494 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
495 >         rb = mol->nextRigidBody(rbIter)) {
496 >      vector<Atom*> atoms = rb->getAtoms();
497 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
498 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
499 >          a = atoms[i]->getGlobalIndex();
500 >          b = atoms[j]->getGlobalIndex();
501 >          excludedInteractions_.addPair(a, b);
502 >        }
503 >      }
504 >    }        
505  
506 <  ri[0] = Hmat[0][0];
266 <  ri[1] = Hmat[1][0];
267 <  ri[2] = Hmat[2][0];
506 >  }
507  
508 <  rj[0] = Hmat[0][1];
509 <  rj[1] = Hmat[1][1];
510 <  rj[2] = Hmat[2][1];
508 >  void SimInfo::removeInteractionPairs(Molecule* mol) {
509 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
510 >    vector<Bond*>::iterator bondIter;
511 >    vector<Bend*>::iterator bendIter;
512 >    vector<Torsion*>::iterator torsionIter;
513 >    vector<Inversion*>::iterator inversionIter;
514 >    Bond* bond;
515 >    Bend* bend;
516 >    Torsion* torsion;
517 >    Inversion* inversion;
518 >    int a;
519 >    int b;
520 >    int c;
521 >    int d;
522  
523 <  rk[0] = Hmat[0][2];
524 <  rk[1] = Hmat[1][2];
525 <  rk[2] = Hmat[2][2];
523 >    map<int, set<int> > atomGroups;
524 >    Molecule::RigidBodyIterator rbIter;
525 >    RigidBody* rb;
526 >    Molecule::IntegrableObjectIterator ii;
527 >    StuntDouble* integrableObject;
528      
529 <  crossProduct3(ri, rj, rij);
530 <  distXY = dotProduct3(rk,rij) / norm3(rij);
529 >    for (integrableObject = mol->beginIntegrableObject(ii);
530 >         integrableObject != NULL;
531 >         integrableObject = mol->nextIntegrableObject(ii)) {
532 >      
533 >      if (integrableObject->isRigidBody()) {
534 >        rb = static_cast<RigidBody*>(integrableObject);
535 >        vector<Atom*> atoms = rb->getAtoms();
536 >        set<int> rigidAtoms;
537 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
538 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
539 >        }
540 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
541 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
542 >        }      
543 >      } else {
544 >        set<int> oneAtomSet;
545 >        oneAtomSet.insert(integrableObject->getGlobalIndex());
546 >        atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
547 >      }
548 >    }  
549  
550 <  crossProduct3(rj,rk, rjk);
551 <  distYZ = dotProduct3(ri,rjk) / norm3(rjk);
550 >    for (bond= mol->beginBond(bondIter); bond != NULL;
551 >         bond = mol->nextBond(bondIter)) {
552 >      
553 >      a = bond->getAtomA()->getGlobalIndex();
554 >      b = bond->getAtomB()->getGlobalIndex();  
555 >    
556 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
557 >        oneTwoInteractions_.removePair(a, b);
558 >      } else {
559 >        excludedInteractions_.removePair(a, b);
560 >      }
561 >    }
562  
563 <  crossProduct3(rk,ri, rki);
564 <  distZX = dotProduct3(rj,rki) / norm3(rki);
563 >    for (bend= mol->beginBend(bendIter); bend != NULL;
564 >         bend = mol->nextBend(bendIter)) {
565  
566 <  minDist = min(min(distXY, distYZ), distZX);
567 <  return minDist/2;
568 <  
569 < }
566 >      a = bend->getAtomA()->getGlobalIndex();
567 >      b = bend->getAtomB()->getGlobalIndex();        
568 >      c = bend->getAtomC()->getGlobalIndex();
569 >      
570 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
571 >        oneTwoInteractions_.removePair(a, b);      
572 >        oneTwoInteractions_.removePair(b, c);
573 >      } else {
574 >        excludedInteractions_.removePair(a, b);
575 >        excludedInteractions_.removePair(b, c);
576 >      }
577  
578 < void SimInfo::wrapVector( double thePos[3] ){
578 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
579 >        oneThreeInteractions_.removePair(a, c);      
580 >      } else {
581 >        excludedInteractions_.removePair(a, c);
582 >      }
583 >    }
584  
585 <  int i;
586 <  double scaled[3];
585 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
586 >         torsion = mol->nextTorsion(torsionIter)) {
587  
588 <  if( !orthoRhombic ){
589 <    // calc the scaled coordinates.
588 >      a = torsion->getAtomA()->getGlobalIndex();
589 >      b = torsion->getAtomB()->getGlobalIndex();        
590 >      c = torsion->getAtomC()->getGlobalIndex();        
591 >      d = torsion->getAtomD()->getGlobalIndex();      
592    
593 +      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
594 +        oneTwoInteractions_.removePair(a, b);      
595 +        oneTwoInteractions_.removePair(b, c);
596 +        oneTwoInteractions_.removePair(c, d);
597 +      } else {
598 +        excludedInteractions_.removePair(a, b);
599 +        excludedInteractions_.removePair(b, c);
600 +        excludedInteractions_.removePair(c, d);
601 +      }
602  
603 <    matVecMul3(HmatInv, thePos, scaled);
604 <    
605 <    for(i=0; i<3; i++)
606 <      scaled[i] -= roundMe(scaled[i]);
607 <    
608 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
609 <    
307 <    matVecMul3(Hmat, scaled, thePos);
603 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
604 >        oneThreeInteractions_.removePair(a, c);      
605 >        oneThreeInteractions_.removePair(b, d);      
606 >      } else {
607 >        excludedInteractions_.removePair(a, c);
608 >        excludedInteractions_.removePair(b, d);
609 >      }
610  
611 <  }
612 <  else{
613 <    // calc the scaled coordinates.
614 <    
615 <    for(i=0; i<3; i++)
616 <      scaled[i] = thePos[i]*HmatInv[i][i];
315 <    
316 <    // wrap the scaled coordinates
317 <    
318 <    for(i=0; i<3; i++)
319 <      scaled[i] -= roundMe(scaled[i]);
320 <    
321 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
322 <    
323 <    for(i=0; i<3; i++)
324 <      thePos[i] = scaled[i]*Hmat[i][i];
325 <  }
326 <    
327 < }
611 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
612 >        oneFourInteractions_.removePair(a, d);      
613 >      } else {
614 >        excludedInteractions_.removePair(a, d);
615 >      }
616 >    }
617  
618 +    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
619 +         inversion = mol->nextInversion(inversionIter)) {
620  
621 < int SimInfo::getNDF(){
622 <  int ndf_local;
621 >      a = inversion->getAtomA()->getGlobalIndex();
622 >      b = inversion->getAtomB()->getGlobalIndex();        
623 >      c = inversion->getAtomC()->getGlobalIndex();        
624 >      d = inversion->getAtomD()->getGlobalIndex();        
625  
626 <  ndf_local = 0;
627 <  
628 <  for(int i = 0; i < integrableObjects.size(); i++){
629 <    ndf_local += 3;
630 <    if (integrableObjects[i]->isDirectional()) {
631 <      if (integrableObjects[i]->isLinear())
632 <        ndf_local += 2;
633 <      else
634 <        ndf_local += 3;
635 <    }
626 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
627 >        oneTwoInteractions_.removePair(a, b);      
628 >        oneTwoInteractions_.removePair(a, c);
629 >        oneTwoInteractions_.removePair(a, d);
630 >      } else {
631 >        excludedInteractions_.removePair(a, b);
632 >        excludedInteractions_.removePair(a, c);
633 >        excludedInteractions_.removePair(a, d);
634 >      }
635 >
636 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
637 >        oneThreeInteractions_.removePair(b, c);    
638 >        oneThreeInteractions_.removePair(b, d);    
639 >        oneThreeInteractions_.removePair(c, d);      
640 >      } else {
641 >        excludedInteractions_.removePair(b, c);
642 >        excludedInteractions_.removePair(b, d);
643 >        excludedInteractions_.removePair(c, d);
644 >      }
645 >    }
646 >
647 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
648 >         rb = mol->nextRigidBody(rbIter)) {
649 >      vector<Atom*> atoms = rb->getAtoms();
650 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
651 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
652 >          a = atoms[i]->getGlobalIndex();
653 >          b = atoms[j]->getGlobalIndex();
654 >          excludedInteractions_.removePair(a, b);
655 >        }
656 >      }
657 >    }        
658 >    
659    }
660 +  
661 +  
662 +  void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
663 +    int curStampId;
664 +    
665 +    //index from 0
666 +    curStampId = moleculeStamps_.size();
667  
668 <  // n_constraints is local, so subtract them on each processor:
668 >    moleculeStamps_.push_back(molStamp);
669 >    molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
670 >  }
671  
347  ndf_local -= n_constraints;
672  
673 < #ifdef IS_MPI
674 <  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
675 < #else
676 <  ndf = ndf_local;
677 < #endif
678 <
679 <  // nZconstraints is global, as are the 3 COM translations for the
680 <  // entire system:
681 <
682 <  ndf = ndf - 3 - nZconstraints;
683 <
684 <  return ndf;
361 < }
362 <
363 < int SimInfo::getNDFraw() {
364 <  int ndfRaw_local;
365 <
366 <  // Raw degrees of freedom that we have to set
367 <  ndfRaw_local = 0;
368 <
369 <  for(int i = 0; i < integrableObjects.size(); i++){
370 <    ndfRaw_local += 3;
371 <    if (integrableObjects[i]->isDirectional()) {
372 <       if (integrableObjects[i]->isLinear())
373 <        ndfRaw_local += 2;
374 <      else
375 <        ndfRaw_local += 3;
376 <    }
673 >  /**
674 >   * update
675 >   *
676 >   *  Performs the global checks and variable settings after the
677 >   *  objects have been created.
678 >   *
679 >   */
680 >  void SimInfo::update() {  
681 >    setupSimVariables();
682 >    calcNdf();
683 >    calcNdfRaw();
684 >    calcNdfTrans();
685    }
686 +  
687 +  /**
688 +   * getSimulatedAtomTypes
689 +   *
690 +   * Returns an STL set of AtomType* that are actually present in this
691 +   * simulation.  Must query all processors to assemble this information.
692 +   *
693 +   */
694 +  set<AtomType*> SimInfo::getSimulatedAtomTypes() {
695 +    SimInfo::MoleculeIterator mi;
696 +    Molecule* mol;
697 +    Molecule::AtomIterator ai;
698 +    Atom* atom;
699 +    set<AtomType*> atomTypes;
700      
701 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
702 +      for(atom = mol->beginAtom(ai); atom != NULL;
703 +          atom = mol->nextAtom(ai)) {
704 +        atomTypes.insert(atom->getAtomType());
705 +      }      
706 +    }    
707 +    
708   #ifdef IS_MPI
380  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
381 #else
382  ndfRaw = ndfRaw_local;
383 #endif
709  
710 <  return ndfRaw;
711 < }
710 >    // loop over the found atom types on this processor, and add their
711 >    // numerical idents to a vector:
712 >    
713 >    vector<int> foundTypes;
714 >    set<AtomType*>::iterator i;
715 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i)
716 >      foundTypes.push_back( (*i)->getIdent() );
717  
718 < int SimInfo::getNDFtranslational() {
719 <  int ndfTrans_local;
718 >    // count_local holds the number of found types on this processor
719 >    int count_local = foundTypes.size();
720  
721 <  ndfTrans_local = 3 * integrableObjects.size() - n_constraints;
721 >    int nproc = MPI::COMM_WORLD.Get_size();
722  
723 +    // we need arrays to hold the counts and displacement vectors for
724 +    // all processors
725 +    vector<int> counts(nproc, 0);
726 +    vector<int> disps(nproc, 0);
727  
728 < #ifdef IS_MPI
729 <  MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
730 < #else
731 <  ndfTrans = ndfTrans_local;
732 < #endif
728 >    // fill the counts array
729 >    MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0],
730 >                              1, MPI::INT);
731 >  
732 >    // use the processor counts to compute the displacement array
733 >    disps[0] = 0;    
734 >    int totalCount = counts[0];
735 >    for (int iproc = 1; iproc < nproc; iproc++) {
736 >      disps[iproc] = disps[iproc-1] + counts[iproc-1];
737 >      totalCount += counts[iproc];
738 >    }
739  
740 <  ndfTrans = ndfTrans - 3 - nZconstraints;
740 >    // we need a (possibly redundant) set of all found types:
741 >    vector<int> ftGlobal(totalCount);
742 >    
743 >    // now spray out the foundTypes to all the other processors:    
744 >    MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT,
745 >                               &ftGlobal[0], &counts[0], &disps[0],
746 >                               MPI::INT);
747  
748 <  return ndfTrans;
403 < }
748 >    vector<int>::iterator j;
749  
750 < int SimInfo::getTotIntegrableObjects() {
751 <  int nObjs_local;
752 <  int nObjs;
750 >    // foundIdents is a stl set, so inserting an already found ident
751 >    // will have no effect.
752 >    set<int> foundIdents;
753  
754 <  nObjs_local =  integrableObjects.size();
755 <
756 <
757 < #ifdef IS_MPI
758 <  MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
759 < #else
760 <  nObjs = nObjs_local;
754 >    for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j)
755 >      foundIdents.insert((*j));
756 >    
757 >    // now iterate over the foundIdents and get the actual atom types
758 >    // that correspond to these:
759 >    set<int>::iterator it;
760 >    for (it = foundIdents.begin(); it != foundIdents.end(); ++it)
761 >      atomTypes.insert( forceField_->getAtomType((*it)) );
762 >
763   #endif
764  
765 <
419 <  return nObjs;
420 < }
421 <
422 < void SimInfo::refreshSim(){
423 <
424 <  simtype fInfo;
425 <  int isError;
426 <  int n_global;
427 <  int* excl;
428 <
429 <  fInfo.dielect = 0.0;
430 <
431 <  if( useDipoles ){
432 <    if( useReactionField )fInfo.dielect = dielectric;
765 >    return atomTypes;        
766    }
767  
768 <  fInfo.SIM_uses_PBC = usePBC;
769 <  //fInfo.SIM_uses_LJ = 0;
770 <  fInfo.SIM_uses_LJ = useLJ;
771 <  fInfo.SIM_uses_sticky = useSticky;
772 <  //fInfo.SIM_uses_sticky = 0;
773 <  fInfo.SIM_uses_charges = useCharges;
774 <  fInfo.SIM_uses_dipoles = useDipoles;
775 <  //fInfo.SIM_uses_dipoles = 0;
776 <  fInfo.SIM_uses_RF = useReactionField;
777 <  //fInfo.SIM_uses_RF = 0;
778 <  fInfo.SIM_uses_GB = useGB;
779 <  fInfo.SIM_uses_EAM = useEAM;
780 <
781 <  n_exclude = excludes->getSize();
782 <  excl = excludes->getFortranArray();
783 <  
784 < #ifdef IS_MPI
785 <  n_global = mpiSim->getNAtomsGlobal();
768 >  void SimInfo::setupSimVariables() {
769 >    useAtomicVirial_ = simParams_->getUseAtomicVirial();
770 >    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
771 >    calcBoxDipole_ = false;
772 >    if ( simParams_->haveAccumulateBoxDipole() )
773 >      if ( simParams_->getAccumulateBoxDipole() ) {
774 >        calcBoxDipole_ = true;      
775 >      }
776 >    
777 >    set<AtomType*>::iterator i;
778 >    set<AtomType*> atomTypes;
779 >    atomTypes = getSimulatedAtomTypes();    
780 >    int usesElectrostatic = 0;
781 >    int usesMetallic = 0;
782 >    int usesDirectional = 0;
783 >    //loop over all of the atom types
784 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
785 >      usesElectrostatic |= (*i)->isElectrostatic();
786 >      usesMetallic |= (*i)->isMetal();
787 >      usesDirectional |= (*i)->isDirectional();
788 >    }
789 >    
790 > #ifdef IS_MPI    
791 >    int temp;
792 >    temp = usesDirectional;
793 >    MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
794 >    
795 >    temp = usesMetallic;
796 >    MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
797 >    
798 >    temp = usesElectrostatic;
799 >    MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
800   #else
454  n_global = n_atoms;
455 #endif
456  
457  isError = 0;
458  
459  getFortranGroupArrays(this, FglobalGroupMembership, mfact);
460  //it may not be a good idea to pass the address of first element in vector
461  //since c++ standard does not require vector to be stored continuously in meomory
462  //Most of the compilers will organize the memory of vector continuously
463  setFortranSim( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
464                  &nGlobalExcludes, globalExcludes, molMembershipArray,
465                  &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError);
801  
802 <  if( isError ){
802 >    usesDirectionalAtoms_ = usesDirectional;
803 >    usesMetallicAtoms_ = usesMetallic;
804 >    usesElectrostaticAtoms_ = usesElectrostatic;
805 >
806 > #endif
807      
808 <    sprintf( painCave.errMsg,
809 <             "There was an error setting the simulation information in fortran.\n" );
810 <    painCave.isFatal = 1;
472 <    painCave.severity = OOPSE_ERROR;
473 <    simError();
808 >    requiresPrepair_ = usesMetallicAtoms_ ? true : false;
809 >    requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false;
810 >    requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false;    
811    }
475  
476 #ifdef IS_MPI
477  sprintf( checkPointMsg,
478           "succesfully sent the simulation information to fortran.\n");
479  MPIcheckPoint();
480 #endif // is_mpi
481  
482  this->ndf = this->getNDF();
483  this->ndfRaw = this->getNDFraw();
484  this->ndfTrans = this->getNDFtranslational();
485 }
812  
487 void SimInfo::setDefaultRcut( double theRcut ){
488  
489  haveRcut = 1;
490  rCut = theRcut;
491  rList = rCut + 1.0;
492  
493  notifyFortranCutoffs( &rCut, &rSw, &rList );
494 }
813  
814 < void SimInfo::setDefaultRcut( double theRcut, double theRsw ){
814 >  vector<int> SimInfo::getGlobalAtomIndices() {
815 >    SimInfo::MoleculeIterator mi;
816 >    Molecule* mol;
817 >    Molecule::AtomIterator ai;
818 >    Atom* atom;
819  
820 <  rSw = theRsw;
499 <  setDefaultRcut( theRcut );
500 < }
501 <
502 <
503 < void SimInfo::checkCutOffs( void ){
504 <  
505 <  if( boxIsInit ){
820 >    vector<int> GlobalAtomIndices(getNAtoms(), 0);
821      
822 <    //we need to check cutOffs against the box
823 <    
824 <    if( rCut > maxCutoff ){
825 <      sprintf( painCave.errMsg,
826 <               "cutoffRadius is too large for the current periodic box.\n"
827 <               "\tCurrent Value of cutoffRadius = %G at time %G\n "
828 <               "\tThis is larger than half of at least one of the\n"
514 <               "\tperiodic box vectors.  Right now, the Box matrix is:\n"
515 <               "\n"
516 <               "\t[ %G %G %G ]\n"
517 <               "\t[ %G %G %G ]\n"
518 <               "\t[ %G %G %G ]\n",
519 <               rCut, currentTime,
520 <               Hmat[0][0], Hmat[0][1], Hmat[0][2],
521 <               Hmat[1][0], Hmat[1][1], Hmat[1][2],
522 <               Hmat[2][0], Hmat[2][1], Hmat[2][2]);
523 <      painCave.severity = OOPSE_ERROR;
524 <      painCave.isFatal = 1;
525 <      simError();
526 <    }    
527 <  } else {
528 <    // initialize this stuff before using it, OK?
529 <    sprintf( painCave.errMsg,
530 <             "Trying to check cutoffs without a box.\n"
531 <             "\tOOPSE should have better programmers than that.\n" );
532 <    painCave.severity = OOPSE_ERROR;
533 <    painCave.isFatal = 1;
534 <    simError();      
822 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
823 >      
824 >      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
825 >        GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex();
826 >      }
827 >    }
828 >    return GlobalAtomIndices;
829    }
536  
537 }
830  
539 void SimInfo::addProperty(GenericData* prop){
831  
832 <  map<string, GenericData*>::iterator result;
833 <  result = properties.find(prop->getID());
834 <  
835 <  //we can't simply use  properties[prop->getID()] = prop,
836 <  //it will cause memory leak if we already contain a propery which has the same name of prop
837 <  
838 <  if(result != properties.end()){
832 >  vector<int> SimInfo::getGlobalGroupIndices() {
833 >    SimInfo::MoleculeIterator mi;
834 >    Molecule* mol;
835 >    Molecule::CutoffGroupIterator ci;
836 >    CutoffGroup* cg;
837 >
838 >    vector<int> GlobalGroupIndices;
839      
840 <    delete (*result).second;
550 <    (*result).second = prop;
840 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
841        
842 +      //local index of cutoff group is trivial, it only depends on the
843 +      //order of travesing
844 +      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
845 +           cg = mol->nextCutoffGroup(ci)) {
846 +        GlobalGroupIndices.push_back(cg->getGlobalIndex());
847 +      }        
848 +    }
849 +    return GlobalGroupIndices;
850    }
553  else{
851  
555    properties[prop->getID()] = prop;
852  
853 <  }
853 >  void SimInfo::prepareTopology() {
854 >    int nExclude, nOneTwo, nOneThree, nOneFour;
855 >
856 >    //calculate mass ratio of cutoff group
857 >    SimInfo::MoleculeIterator mi;
858 >    Molecule* mol;
859 >    Molecule::CutoffGroupIterator ci;
860 >    CutoffGroup* cg;
861 >    Molecule::AtomIterator ai;
862 >    Atom* atom;
863 >    RealType totalMass;
864 >
865 >    /**
866 >     * The mass factor is the relative mass of an atom to the total
867 >     * mass of the cutoff group it belongs to.  By default, all atoms
868 >     * are their own cutoff groups, and therefore have mass factors of
869 >     * 1.  We need some special handling for massless atoms, which
870 >     * will be treated as carrying the entire mass of the cutoff
871 >     * group.
872 >     */
873 >    massFactors_.clear();
874 >    massFactors_.resize(getNAtoms(), 1.0);
875      
876 < }
876 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
877 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
878 >           cg = mol->nextCutoffGroup(ci)) {
879  
880 < GenericData* SimInfo::getProperty(const string& propName){
881 <
882 <  map<string, GenericData*>::iterator result;
883 <  
884 <  //string lowerCaseName = ();
885 <  
886 <  result = properties.find(propName);
887 <  
888 <  if(result != properties.end())
889 <    return (*result).second;  
571 <  else  
572 <    return NULL;  
573 < }
880 >        totalMass = cg->getMass();
881 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
882 >          // Check for massless groups - set mfact to 1 if true
883 >          if (totalMass != 0)
884 >            massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass;
885 >          else
886 >            massFactors_[atom->getLocalIndex()] = 1.0;
887 >        }
888 >      }      
889 >    }
890  
891 +    // Build the identArray_
892  
893 < void SimInfo::getFortranGroupArrays(SimInfo* info,
894 <                                    vector<int>& FglobalGroupMembership,
895 <                                    vector<double>& mfact){
896 <  
897 <  Molecule* myMols;
898 <  Atom** myAtoms;
899 <  int numAtom;
900 <  double mtot;
901 <  int numMol;
585 <  int numCutoffGroups;
586 <  CutoffGroup* myCutoffGroup;
587 <  vector<CutoffGroup*>::iterator iterCutoff;
588 <  Atom* cutoffAtom;
589 <  vector<Atom*>::iterator iterAtom;
590 <  int atomIndex;
591 <  double totalMass;
592 <  
593 <  mfact.clear();
594 <  FglobalGroupMembership.clear();
595 <  
893 >    identArray_.clear();
894 >    identArray_.reserve(getNAtoms());    
895 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
896 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
897 >        identArray_.push_back(atom->getIdent());
898 >      }
899 >    }    
900 >    
901 >    //scan topology
902  
903 <  // Fix the silly fortran indexing problem
903 >    nExclude = excludedInteractions_.getSize();
904 >    nOneTwo = oneTwoInteractions_.getSize();
905 >    nOneThree = oneThreeInteractions_.getSize();
906 >    nOneFour = oneFourInteractions_.getSize();
907 >
908 >    int* excludeList = excludedInteractions_.getPairList();
909 >    int* oneTwoList = oneTwoInteractions_.getPairList();
910 >    int* oneThreeList = oneThreeInteractions_.getPairList();
911 >    int* oneFourList = oneFourInteractions_.getPairList();
912 >
913 >    topologyDone_ = true;
914 >  }
915 >
916 >  void SimInfo::addProperty(GenericData* genData) {
917 >    properties_.addProperty(genData);  
918 >  }
919 >
920 >  void SimInfo::removeProperty(const string& propName) {
921 >    properties_.removeProperty(propName);  
922 >  }
923 >
924 >  void SimInfo::clearProperties() {
925 >    properties_.clearProperties();
926 >  }
927 >
928 >  vector<string> SimInfo::getPropertyNames() {
929 >    return properties_.getPropertyNames();  
930 >  }
931 >      
932 >  vector<GenericData*> SimInfo::getProperties() {
933 >    return properties_.getProperties();
934 >  }
935 >
936 >  GenericData* SimInfo::getPropertyByName(const string& propName) {
937 >    return properties_.getPropertyByName(propName);
938 >  }
939 >
940 >  void SimInfo::setSnapshotManager(SnapshotManager* sman) {
941 >    if (sman_ == sman) {
942 >      return;
943 >    }    
944 >    delete sman_;
945 >    sman_ = sman;
946 >
947 >    Molecule* mol;
948 >    RigidBody* rb;
949 >    Atom* atom;
950 >    CutoffGroup* cg;
951 >    SimInfo::MoleculeIterator mi;
952 >    Molecule::RigidBodyIterator rbIter;
953 >    Molecule::AtomIterator atomIter;
954 >    Molecule::CutoffGroupIterator cgIter;
955 >
956 >    for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
957 >        
958 >      for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
959 >        atom->setSnapshotManager(sman_);
960 >      }
961 >        
962 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
963 >        rb->setSnapshotManager(sman_);
964 >      }
965 >
966 >      for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) {
967 >        cg->setSnapshotManager(sman_);
968 >      }
969 >    }    
970 >    
971 >  }
972 >
973 >  Vector3d SimInfo::getComVel(){
974 >    SimInfo::MoleculeIterator i;
975 >    Molecule* mol;
976 >
977 >    Vector3d comVel(0.0);
978 >    RealType totalMass = 0.0;
979 >    
980 >
981 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
982 >      RealType mass = mol->getMass();
983 >      totalMass += mass;
984 >      comVel += mass * mol->getComVel();
985 >    }  
986 >
987   #ifdef IS_MPI
988 <  numAtom = mpiSim->getNAtomsGlobal();
989 < #else
990 <  numAtom = n_atoms;
988 >    RealType tmpMass = totalMass;
989 >    Vector3d tmpComVel(comVel);    
990 >    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
991 >    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
992   #endif
603  for (int i = 0; i < numAtom; i++)
604    FglobalGroupMembership.push_back(globalGroupMembership[i] + 1);
605  
993  
994 <  myMols = info->molecules;
608 <  numMol = info->n_mol;
609 <  for(int i  = 0; i < numMol; i++){
610 <    numCutoffGroups = myMols[i].getNCutoffGroups();
611 <    for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff);
612 <        myCutoffGroup != NULL;
613 <        myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){
994 >    comVel /= totalMass;
995  
996 <      totalMass = myCutoffGroup->getMass();
996 >    return comVel;
997 >  }
998 >
999 >  Vector3d SimInfo::getCom(){
1000 >    SimInfo::MoleculeIterator i;
1001 >    Molecule* mol;
1002 >
1003 >    Vector3d com(0.0);
1004 >    RealType totalMass = 0.0;
1005 >    
1006 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1007 >      RealType mass = mol->getMass();
1008 >      totalMass += mass;
1009 >      com += mass * mol->getCom();
1010 >    }  
1011 >
1012 > #ifdef IS_MPI
1013 >    RealType tmpMass = totalMass;
1014 >    Vector3d tmpCom(com);    
1015 >    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1016 >    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1017 > #endif
1018 >
1019 >    com /= totalMass;
1020 >
1021 >    return com;
1022 >
1023 >  }        
1024 >
1025 >  ostream& operator <<(ostream& o, SimInfo& info) {
1026 >
1027 >    return o;
1028 >  }
1029 >  
1030 >  
1031 >   /*
1032 >   Returns center of mass and center of mass velocity in one function call.
1033 >   */
1034 >  
1035 >   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1036 >      SimInfo::MoleculeIterator i;
1037 >      Molecule* mol;
1038        
1039 <      for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom);
1040 <          cutoffAtom != NULL;
1041 <          cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){
1042 <        mfact.push_back(cutoffAtom->getMass()/totalMass);
1039 >    
1040 >      RealType totalMass = 0.0;
1041 >    
1042 >
1043 >      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1044 >         RealType mass = mol->getMass();
1045 >         totalMass += mass;
1046 >         com += mass * mol->getCom();
1047 >         comVel += mass * mol->getComVel();          
1048        }  
1049 +      
1050 + #ifdef IS_MPI
1051 +      RealType tmpMass = totalMass;
1052 +      Vector3d tmpCom(com);  
1053 +      Vector3d tmpComVel(comVel);
1054 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1055 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1056 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1057 + #endif
1058 +      
1059 +      com /= totalMass;
1060 +      comVel /= totalMass;
1061 +   }        
1062 +  
1063 +   /*
1064 +   Return intertia tensor for entire system and angular momentum Vector.
1065 +
1066 +
1067 +       [  Ixx -Ixy  -Ixz ]
1068 +    J =| -Iyx  Iyy  -Iyz |
1069 +       [ -Izx -Iyz   Izz ]
1070 +    */
1071 +
1072 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1073 +      
1074 +
1075 +      RealType xx = 0.0;
1076 +      RealType yy = 0.0;
1077 +      RealType zz = 0.0;
1078 +      RealType xy = 0.0;
1079 +      RealType xz = 0.0;
1080 +      RealType yz = 0.0;
1081 +      Vector3d com(0.0);
1082 +      Vector3d comVel(0.0);
1083 +      
1084 +      getComAll(com, comVel);
1085 +      
1086 +      SimInfo::MoleculeIterator i;
1087 +      Molecule* mol;
1088 +      
1089 +      Vector3d thisq(0.0);
1090 +      Vector3d thisv(0.0);
1091 +
1092 +      RealType thisMass = 0.0;
1093 +    
1094 +      
1095 +      
1096 +  
1097 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1098 +        
1099 +         thisq = mol->getCom()-com;
1100 +         thisv = mol->getComVel()-comVel;
1101 +         thisMass = mol->getMass();
1102 +         // Compute moment of intertia coefficients.
1103 +         xx += thisq[0]*thisq[0]*thisMass;
1104 +         yy += thisq[1]*thisq[1]*thisMass;
1105 +         zz += thisq[2]*thisq[2]*thisMass;
1106 +        
1107 +         // compute products of intertia
1108 +         xy += thisq[0]*thisq[1]*thisMass;
1109 +         xz += thisq[0]*thisq[2]*thisMass;
1110 +         yz += thisq[1]*thisq[2]*thisMass;
1111 +            
1112 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1113 +            
1114 +      }  
1115 +      
1116 +      
1117 +      inertiaTensor(0,0) = yy + zz;
1118 +      inertiaTensor(0,1) = -xy;
1119 +      inertiaTensor(0,2) = -xz;
1120 +      inertiaTensor(1,0) = -xy;
1121 +      inertiaTensor(1,1) = xx + zz;
1122 +      inertiaTensor(1,2) = -yz;
1123 +      inertiaTensor(2,0) = -xz;
1124 +      inertiaTensor(2,1) = -yz;
1125 +      inertiaTensor(2,2) = xx + yy;
1126 +      
1127 + #ifdef IS_MPI
1128 +      Mat3x3d tmpI(inertiaTensor);
1129 +      Vector3d tmpAngMom;
1130 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1131 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1132 + #endif
1133 +              
1134 +      return;
1135 +   }
1136 +
1137 +   //Returns the angular momentum of the system
1138 +   Vector3d SimInfo::getAngularMomentum(){
1139 +      
1140 +      Vector3d com(0.0);
1141 +      Vector3d comVel(0.0);
1142 +      Vector3d angularMomentum(0.0);
1143 +      
1144 +      getComAll(com,comVel);
1145 +      
1146 +      SimInfo::MoleculeIterator i;
1147 +      Molecule* mol;
1148 +      
1149 +      Vector3d thisr(0.0);
1150 +      Vector3d thisp(0.0);
1151 +      
1152 +      RealType thisMass;
1153 +      
1154 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1155 +        thisMass = mol->getMass();
1156 +        thisr = mol->getCom()-com;
1157 +        thisp = (mol->getComVel()-comVel)*thisMass;
1158 +        
1159 +        angularMomentum += cross( thisr, thisp );
1160 +        
1161 +      }  
1162 +      
1163 + #ifdef IS_MPI
1164 +      Vector3d tmpAngMom;
1165 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1166 + #endif
1167 +      
1168 +      return angularMomentum;
1169 +   }
1170 +  
1171 +  StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
1172 +    return IOIndexToIntegrableObject.at(index);
1173 +  }
1174 +  
1175 +  void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) {
1176 +    IOIndexToIntegrableObject= v;
1177 +  }
1178 +
1179 +  /* Returns the Volume of the simulation based on a ellipsoid with semi-axes
1180 +     based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3
1181 +     where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to
1182 +     V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
1183 +  */
1184 +  void SimInfo::getGyrationalVolume(RealType &volume){
1185 +    Mat3x3d intTensor;
1186 +    RealType det;
1187 +    Vector3d dummyAngMom;
1188 +    RealType sysconstants;
1189 +    RealType geomCnst;
1190 +
1191 +    geomCnst = 3.0/2.0;
1192 +    /* Get the inertial tensor and angular momentum for free*/
1193 +    getInertiaTensor(intTensor,dummyAngMom);
1194 +    
1195 +    det = intTensor.determinant();
1196 +    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1197 +    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det);
1198 +    return;
1199 +  }
1200 +
1201 +  void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){
1202 +    Mat3x3d intTensor;
1203 +    Vector3d dummyAngMom;
1204 +    RealType sysconstants;
1205 +    RealType geomCnst;
1206 +
1207 +    geomCnst = 3.0/2.0;
1208 +    /* Get the inertial tensor and angular momentum for free*/
1209 +    getInertiaTensor(intTensor,dummyAngMom);
1210 +    
1211 +    detI = intTensor.determinant();
1212 +    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1213 +    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI);
1214 +    return;
1215 +  }
1216 + /*
1217 +   void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) {
1218 +      assert( v.size() == nAtoms_ + nRigidBodies_);
1219 +      sdByGlobalIndex_ = v;
1220      }
1221 +
1222 +    StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) {
1223 +      //assert(index < nAtoms_ + nRigidBodies_);
1224 +      return sdByGlobalIndex_.at(index);
1225 +    }  
1226 + */  
1227 +  int SimInfo::getNGlobalConstraints() {
1228 +    int nGlobalConstraints;
1229 + #ifdef IS_MPI
1230 +    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
1231 +                  MPI_COMM_WORLD);    
1232 + #else
1233 +    nGlobalConstraints =  nConstraints_;
1234 + #endif
1235 +    return nGlobalConstraints;
1236    }
1237  
1238 < }
1238 > }//end namespace OpenMD
1239 >

Comparing:
trunk/src/brains/SimInfo.cpp (property svn:keywords), Revision 124 by chuckv, Wed Oct 20 20:46:20 2004 UTC vs.
branches/development/src/brains/SimInfo.cpp (property svn:keywords), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines