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

Comparing:
trunk/src/brains/SimInfo.cpp (property svn:keywords), Revision 143 by chrisfen, Fri Oct 22 22:54:01 2004 UTC vs.
branches/development/src/brains/SimInfo.cpp (property svn:keywords), Revision 1586 by gezelter, Tue Jun 21 06:34:35 2011 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines