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 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
branches/development/src/brains/SimInfo.cpp (file contents), Revision 1597 by gezelter, Tue Jul 26 15:49:24 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 "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  
17 #ifdef IS_MPI
18 #include "brains/mpiSimulation.hpp"
19 #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 < }
163 <
164 < SimInfo* currentInfo;
165 <
166 < SimInfo::SimInfo(){
167 <
168 <  n_constraints = 0;
169 <  nZconstraints = 0;
170 <  n_oriented = 0;
171 <  n_dipoles = 0;
172 <  ndf = 0;
173 <  ndfRaw = 0;
174 <  nZconstraints = 0;
175 <  the_integrator = NULL;
176 <  setTemp = 0;
177 <  thermalTime = 0.0;
178 <  currentTime = 0.0;
179 <  rCut = 0.0;
180 <  rSw = 0.0;
46 <
47 <  haveRcut = 0;
48 <  haveRsw = 0;
49 <  boxIsInit = 0;
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 <  resetTime = 1e99;
182 >  bool SimInfo::removeMolecule(Molecule* mol) {
183 >    MoleculeIterator i;
184 >    i = molecules_.find(mol->getGlobalIndex());
185  
186 <  orthoRhombic = 0;
54 <  orthoTolerance = 1E-6;
55 <  useInitXSstate = true;
186 >    if (i != molecules_.end() ) {
187  
188 <  usePBC = 0;
189 <  useLJ = 0;
190 <  useSticky = 0;
191 <  useCharges = 0;
192 <  useDipoles = 0;
193 <  useReactionField = 0;
194 <  useGB = 0;
195 <  useEAM = 0;
196 <  useSolidThermInt = 0;
197 <  useLiquidThermInt = 0;
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 <  haveCutoffGroups = false;
200 >      removeInteractionPairs(mol);
201 >      molecules_.erase(mol->getGlobalIndex());
202  
203 <  excludes = Exclude::Instance();
203 >      delete mol;
204 >        
205 >      return true;
206 >    } else {
207 >      return false;
208 >    }
209 >  }    
210  
211 <  myConfiguration = new SimState();
211 >        
212 >  Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
213 >    i = molecules_.begin();
214 >    return i == molecules_.end() ? NULL : i->second;
215 >  }    
216  
217 <  has_minimizer = false;
218 <  the_minimizer =NULL;
217 >  Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
218 >    ++i;
219 >    return i == molecules_.end() ? NULL : i->second;    
220 >  }
221  
77  ngroup = 0;
222  
223 <  wrapMeSimInfo( this );
224 < }
223 >  void SimInfo::calcNdf() {
224 >    int ndf_local;
225 >    MoleculeIterator i;
226 >    vector<StuntDouble*>::iterator j;
227 >    Molecule* mol;
228 >    StuntDouble* integrableObject;
229  
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 < SimInfo::~SimInfo(){
236 >        ndf_local += 3;
237  
238 <  delete myConfiguration;
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 <  map<string, GenericData*>::iterator i;
253 <  
254 <  for(i = properties.begin(); i != properties.end(); i++)
255 <    delete (*i).second;
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 < void SimInfo::setBox(double newBox[3]) {
95 <  
96 <  int i, j;
97 <  double tempMat[3][3];
262 >  }
263  
264 <  for(i=0; i<3; i++)
265 <    for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
266 <
267 <  tempMat[0][0] = newBox[0];
268 <  tempMat[1][1] = newBox[1];
269 <  tempMat[2][2] = newBox[2];
270 <
271 <  setBoxM( tempMat );
107 <
108 < }
109 <
110 < void SimInfo::setBoxM( double theBox[3][3] ){
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 <  int i, j;
274 <  double FortranHmat[9]; // to preserve compatibility with Fortran the
275 <                         // ordering in the array is as follows:
276 <                         // [ 0 3 6 ]
277 <                         // [ 1 4 7 ]
278 <                         // [ 2 5 8 ]
279 <  double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
280 <
281 <  if( !boxIsInit ) boxIsInit = 1;
282 <
283 <  for(i=0; i < 3; i++)
284 <    for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
285 <  
286 <  calcBoxL();
126 <  calcHmatInv();
127 <
128 <  for(i=0; i < 3; i++) {
129 <    for (j=0; j < 3; j++) {
130 <      FortranHmat[3*j + i] = Hmat[i][j];
131 <      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 <  setFortranBoxSize(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;
143 <  for(i=0; i<3; i++)
144 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
145 < }
307 >        ndfRaw_local += 3;
308  
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 < void SimInfo::scaleBox(double scale) {
328 <  double theBox[3][3];
150 <  int i, j;
327 >  void SimInfo::calcNdfTrans() {
328 >    int ndfTrans_local;
329  
330 <  // cerr << "Scaling box by " << scale << "\n";
330 >    ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
331  
154  for(i=0; i<3; i++)
155    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
332  
333 <  setBoxM(theBox);
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 < }
339 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
340 >
341 >  }
342  
343 < void SimInfo::calcHmatInv( void ) {
344 <  
345 <  int oldOrtho;
346 <  int i,j;
347 <  double smallDiag;
348 <  double tol;
349 <  double sanity[3][3];
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 <  invertMat3( Hmat, HmatInv );
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 <  // check to see if Hmat is orthorhombic
367 <  
368 <  oldOrtho = orthoRhombic;
369 <
370 <  smallDiag = fabs(Hmat[0][0]);
371 <  if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
372 <  if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
373 <  tol = smallDiag * orthoTolerance;
374 <
375 <  orthoRhombic = 1;
376 <  
377 <  for (i = 0; i < 3; i++ ) {
378 <    for (j = 0 ; j < 3; j++) {
379 <      if (i != j) {
380 <        if (orthoRhombic) {
381 <          if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
382 <        }        
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 <  }
391 >    }  
392 >          
393 >    for (bond= mol->beginBond(bondIter); bond != NULL;
394 >         bond = mol->nextBond(bondIter)) {
395  
396 <  if( oldOrtho != orthoRhombic ){
396 >      a = bond->getAtomA()->getGlobalIndex();
397 >      b = bond->getAtomB()->getGlobalIndex();  
398      
399 <    if( orthoRhombic ) {
400 <      sprintf( painCave.errMsg,
401 <               "OOPSE is switching from the default Non-Orthorhombic\n"
402 <               "\tto the faster Orthorhombic periodic boundary computations.\n"
403 <               "\tThis is usually a good thing, but if you wan't the\n"
199 <               "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
200 <               "\tvariable ( currently set to %G ) smaller.\n",
201 <               orthoTolerance);
202 <      painCave.severity = OOPSE_INFO;
203 <      simError();
399 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
400 >        oneTwoInteractions_.addPair(a, b);
401 >      } else {
402 >        excludedInteractions_.addPair(a, b);
403 >      }
404      }
205    else {
206      sprintf( painCave.errMsg,
207               "OOPSE is switching from the faster Orthorhombic to the more\n"
208               "\tflexible Non-Orthorhombic periodic boundary computations.\n"
209               "\tThis is usually because the box has deformed under\n"
210               "\tNPTf integration. If you wan't to live on the edge with\n"
211               "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
212               "\tvariable ( currently set to %G ) larger.\n",
213               orthoTolerance);
214      painCave.severity = OOPSE_WARNING;
215      simError();
216    }
217  }
218 }
405  
406 < void SimInfo::calcBoxL( void ){
406 >    for (bend= mol->beginBend(bendIter); bend != NULL;
407 >         bend = mol->nextBend(bendIter)) {
408  
409 <  double dx, dy, dz, dsq;
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 <  // boxVol = Determinant of Hmat
421 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
422 >        oneThreeInteractions_.addPair(a, c);      
423 >      } else {
424 >        excludedInteractions_.addPair(a, c);
425 >      }
426 >    }
427  
428 <  boxVol = matDet3( Hmat );
428 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
429 >         torsion = mol->nextTorsion(torsionIter)) {
430  
431 <  // boxLx
432 <  
433 <  dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
434 <  dsq = dx*dx + dy*dy + dz*dz;
232 <  boxL[0] = sqrt( dsq );
233 <  //maxCutoff = 0.5 * boxL[0];
431 >      a = torsion->getAtomA()->getGlobalIndex();
432 >      b = torsion->getAtomB()->getGlobalIndex();        
433 >      c = torsion->getAtomC()->getGlobalIndex();        
434 >      d = torsion->getAtomD()->getGlobalIndex();      
435  
436 <  // boxLy
437 <  
438 <  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
439 <  dsq = dx*dx + dy*dy + dz*dz;
440 <  boxL[1] = sqrt( dsq );
441 <  //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
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 +      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 <  // boxLz
455 <  
456 <  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
457 <  dsq = dx*dx + dy*dy + dz*dz;
458 <  boxL[2] = sqrt( dsq );
459 <  //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
454 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
455 >        oneFourInteractions_.addPair(a, d);      
456 >      } else {
457 >        excludedInteractions_.addPair(a, d);
458 >      }
459 >    }
460  
461 <  //calculate the max cutoff
462 <  maxCutoff =  calcMaxCutOff();
252 <  
253 <  checkCutOffs();
461 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
462 >         inversion = mol->nextInversion(inversionIter)) {
463  
464 < }
464 >      a = inversion->getAtomA()->getGlobalIndex();
465 >      b = inversion->getAtomB()->getGlobalIndex();        
466 >      c = inversion->getAtomC()->getGlobalIndex();        
467 >      d = inversion->getAtomD()->getGlobalIndex();        
468  
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 < double SimInfo::calcMaxCutOff(){
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 <  double ri[3], rj[3], rk[3];
491 <  double rij[3], rjk[3], rki[3];
492 <  double minDist;
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 <  ri[0] = Hmat[0][0];
265 <  ri[1] = Hmat[1][0];
266 <  ri[2] = Hmat[2][0];
502 >  }
503  
504 <  rj[0] = Hmat[0][1];
505 <  rj[1] = Hmat[1][1];
506 <  rj[2] = Hmat[2][1];
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 <  rk[0] = Hmat[0][2];
520 <  rk[1] = Hmat[1][2];
521 <  rk[2] = Hmat[2][2];
519 >    map<int, set<int> > atomGroups;
520 >    Molecule::RigidBodyIterator rbIter;
521 >    RigidBody* rb;
522 >    Molecule::IntegrableObjectIterator ii;
523 >    StuntDouble* integrableObject;
524      
525 <  crossProduct3(ri, rj, rij);
526 <  distXY = dotProduct3(rk,rij) / norm3(rij);
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 <  crossProduct3(rj,rk, rjk);
547 <  distYZ = dotProduct3(ri,rjk) / norm3(rjk);
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 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
553 >        oneTwoInteractions_.removePair(a, b);
554 >      } else {
555 >        excludedInteractions_.removePair(a, b);
556 >      }
557 >    }
558  
559 <  crossProduct3(rk,ri, rki);
560 <  distZX = dotProduct3(rj,rki) / norm3(rki);
559 >    for (bend= mol->beginBend(bendIter); bend != NULL;
560 >         bend = mol->nextBend(bendIter)) {
561  
562 <  minDist = min(min(distXY, distYZ), distZX);
563 <  return minDist/2;
564 <  
565 < }
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 < void SimInfo::wrapVector( double thePos[3] ){
574 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
575 >        oneThreeInteractions_.removePair(a, c);      
576 >      } else {
577 >        excludedInteractions_.removePair(a, c);
578 >      }
579 >    }
580  
581 <  int i;
582 <  double scaled[3];
581 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
582 >         torsion = mol->nextTorsion(torsionIter)) {
583  
584 <  if( !orthoRhombic ){
585 <    // calc the scaled coordinates.
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 <    matVecMul3(HmatInv, thePos, scaled);
600 <    
601 <    for(i=0; i<3; i++)
602 <      scaled[i] -= roundMe(scaled[i]);
603 <    
604 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
605 <    
306 <    matVecMul3(Hmat, scaled, thePos);
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 <  }
608 <  else{
609 <    // calc the scaled coordinates.
610 <    
611 <    for(i=0; i<3; i++)
612 <      scaled[i] = thePos[i]*HmatInv[i][i];
314 <    
315 <    // wrap the scaled coordinates
316 <    
317 <    for(i=0; i<3; i++)
318 <      scaled[i] -= roundMe(scaled[i]);
319 <    
320 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
321 <    
322 <    for(i=0; i<3; i++)
323 <      thePos[i] = scaled[i]*Hmat[i][i];
324 <  }
325 <    
326 < }
607 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
608 >        oneFourInteractions_.removePair(a, d);      
609 >      } else {
610 >        excludedInteractions_.removePair(a, d);
611 >      }
612 >    }
613  
614 +    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
615 +         inversion = mol->nextInversion(inversionIter)) {
616  
617 < int SimInfo::getNDF(){
618 <  int ndf_local;
617 >      a = inversion->getAtomA()->getGlobalIndex();
618 >      b = inversion->getAtomB()->getGlobalIndex();        
619 >      c = inversion->getAtomC()->getGlobalIndex();        
620 >      d = inversion->getAtomD()->getGlobalIndex();        
621  
622 <  ndf_local = 0;
623 <  
624 <  for(int i = 0; i < integrableObjects.size(); i++){
625 <    ndf_local += 3;
626 <    if (integrableObjects[i]->isDirectional()) {
627 <      if (integrableObjects[i]->isLinear())
628 <        ndf_local += 2;
629 <      else
630 <        ndf_local += 3;
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 >      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      }
342  }
642  
643 <  // n_constraints is local, so subtract them on each processor:
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 >    //index from 0
662 >    curStampId = moleculeStamps_.size();
663  
664 <  ndf_local -= n_constraints;
664 >    moleculeStamps_.push_back(molStamp);
665 >    molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
666 >  }
667  
348 #ifdef IS_MPI
349  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
350 #else
351  ndf = ndf_local;
352 #endif
668  
669 <  // nZconstraints is global, as are the 3 COM translations for the
670 <  // entire system:
671 <
672 <  ndf = ndf - 3 - nZconstraints;
673 <
674 <  return ndf;
675 < }
676 <
677 < int SimInfo::getNDFraw() {
678 <  int ndfRaw_local;
679 <
680 <  // Raw degrees of freedom that we have to set
366 <  ndfRaw_local = 0;
367 <
368 <  for(int i = 0; i < integrableObjects.size(); i++){
369 <    ndfRaw_local += 3;
370 <    if (integrableObjects[i]->isDirectional()) {
371 <       if (integrableObjects[i]->isLinear())
372 <        ndfRaw_local += 2;
373 <      else
374 <        ndfRaw_local += 3;
375 <    }
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;
699 +          atom = mol->nextAtom(ai)) {
700 +        atomTypes.insert(atom->getAtomType());
701 +      }      
702 +    }    
703 +    
704   #ifdef IS_MPI
379  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
380 #else
381  ndfRaw = ndfRaw_local;
382 #endif
705  
706 <  return ndfRaw;
707 < }
706 >    // loop over the found atom types on this processor, and add their
707 >    // numerical idents to a vector:
708 >    
709 >    vector<int> foundTypes;
710 >    set<AtomType*>::iterator i;
711 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i)
712 >      foundTypes.push_back( (*i)->getIdent() );
713  
714 < int SimInfo::getNDFtranslational() {
715 <  int ndfTrans_local;
714 >    // count_local holds the number of found types on this processor
715 >    int count_local = foundTypes.size();
716  
717 <  ndfTrans_local = 3 * integrableObjects.size() - n_constraints;
717 >    int nproc = MPI::COMM_WORLD.Get_size();
718  
719 +    // we need arrays to hold the counts and displacement vectors for
720 +    // all processors
721 +    vector<int> counts(nproc, 0);
722 +    vector<int> disps(nproc, 0);
723  
724 < #ifdef IS_MPI
725 <  MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
726 < #else
727 <  ndfTrans = ndfTrans_local;
728 < #endif
724 >    // fill the counts array
725 >    MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0],
726 >                              1, MPI::INT);
727 >  
728 >    // use the processor counts to compute the displacement array
729 >    disps[0] = 0;    
730 >    int totalCount = counts[0];
731 >    for (int iproc = 1; iproc < nproc; iproc++) {
732 >      disps[iproc] = disps[iproc-1] + counts[iproc-1];
733 >      totalCount += counts[iproc];
734 >    }
735  
736 <  ndfTrans = ndfTrans - 3 - nZconstraints;
736 >    // we need a (possibly redundant) set of all found types:
737 >    vector<int> ftGlobal(totalCount);
738 >    
739 >    // now spray out the foundTypes to all the other processors:    
740 >    MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT,
741 >                               &ftGlobal[0], &counts[0], &disps[0],
742 >                               MPI::INT);
743  
744 <  return ndfTrans;
402 < }
744 >    vector<int>::iterator j;
745  
746 < int SimInfo::getTotIntegrableObjects() {
747 <  int nObjs_local;
748 <  int nObjs;
746 >    // foundIdents is a stl set, so inserting an already found ident
747 >    // will have no effect.
748 >    set<int> foundIdents;
749  
750 <  nObjs_local =  integrableObjects.size();
751 <
752 <
753 < #ifdef IS_MPI
754 <  MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
755 < #else
756 <  nObjs = nObjs_local;
750 >    for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j)
751 >      foundIdents.insert((*j));
752 >    
753 >    // now iterate over the foundIdents and get the actual atom types
754 >    // that correspond to these:
755 >    set<int>::iterator it;
756 >    for (it = foundIdents.begin(); it != foundIdents.end(); ++it)
757 >      atomTypes.insert( forceField_->getAtomType((*it)) );
758 >
759   #endif
760  
761 <
418 <  return nObjs;
419 < }
420 <
421 < void SimInfo::refreshSim(){
422 <
423 <  simtype fInfo;
424 <  int isError;
425 <  int n_global;
426 <  int* excl;
427 <
428 <  fInfo.dielect = 0.0;
429 <
430 <  if( useDipoles ){
431 <    if( useReactionField )fInfo.dielect = dielectric;
761 >    return atomTypes;        
762    }
763  
764 <  fInfo.SIM_uses_PBC = usePBC;
765 <  //fInfo.SIM_uses_LJ = 0;
766 <  fInfo.SIM_uses_LJ = useLJ;
767 <  fInfo.SIM_uses_sticky = useSticky;
768 <  //fInfo.SIM_uses_sticky = 0;
769 <  fInfo.SIM_uses_charges = useCharges;
770 <  fInfo.SIM_uses_dipoles = useDipoles;
771 <  //fInfo.SIM_uses_dipoles = 0;
772 <  fInfo.SIM_uses_RF = useReactionField;
773 <  //fInfo.SIM_uses_RF = 0;
774 <  fInfo.SIM_uses_GB = useGB;
775 <  fInfo.SIM_uses_EAM = useEAM;
776 <
777 <  n_exclude = excludes->getSize();
778 <  excl = excludes->getFortranArray();
779 <  
780 < #ifdef IS_MPI
781 <  n_global = mpiSim->getNAtomsGlobal();
764 >  void SimInfo::setupSimVariables() {
765 >    useAtomicVirial_ = simParams_->getUseAtomicVirial();
766 >    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
767 >    calcBoxDipole_ = false;
768 >    if ( simParams_->haveAccumulateBoxDipole() )
769 >      if ( simParams_->getAccumulateBoxDipole() ) {
770 >        calcBoxDipole_ = true;      
771 >      }
772 >    
773 >    set<AtomType*>::iterator i;
774 >    set<AtomType*> atomTypes;
775 >    atomTypes = getSimulatedAtomTypes();    
776 >    int usesElectrostatic = 0;
777 >    int usesMetallic = 0;
778 >    int usesDirectional = 0;
779 >    //loop over all of the atom types
780 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
781 >      usesElectrostatic |= (*i)->isElectrostatic();
782 >      usesMetallic |= (*i)->isMetal();
783 >      usesDirectional |= (*i)->isDirectional();
784 >    }
785 >    
786 > #ifdef IS_MPI    
787 >    int temp;
788 >    temp = usesDirectional;
789 >    MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
790 >    
791 >    temp = usesMetallic;
792 >    MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
793 >    
794 >    temp = usesElectrostatic;
795 >    MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
796   #else
453  n_global = n_atoms;
454 #endif
455  
456  isError = 0;
457  
458  getFortranGroupArrays(this, FglobalGroupMembership, mfact);
459  //it may not be a good idea to pass the address of first element in vector
460  //since c++ standard does not require vector to be stored continuously in meomory
461  //Most of the compilers will organize the memory of vector continuously
462  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
463                  &nGlobalExcludes, globalExcludes, molMembershipArray,
464                  &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError);
797  
798 <  if( isError ){
798 >    usesDirectionalAtoms_ = usesDirectional;
799 >    usesMetallicAtoms_ = usesMetallic;
800 >    usesElectrostaticAtoms_ = usesElectrostatic;
801 >
802 > #endif
803      
804 <    sprintf( painCave.errMsg,
805 <             "There was an error setting the simulation information in fortran.\n" );
806 <    painCave.isFatal = 1;
471 <    painCave.severity = OOPSE_ERROR;
472 <    simError();
804 >    requiresPrepair_ = usesMetallicAtoms_ ? true : false;
805 >    requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false;
806 >    requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false;    
807    }
474  
475 #ifdef IS_MPI
476  sprintf( checkPointMsg,
477           "succesfully sent the simulation information to fortran.\n");
478  MPIcheckPoint();
479 #endif // is_mpi
480  
481  this->ndf = this->getNDF();
482  this->ndfRaw = this->getNDFraw();
483  this->ndfTrans = this->getNDFtranslational();
484 }
808  
486 void SimInfo::setDefaultRcut( double theRcut ){
487  
488  haveRcut = 1;
489  rCut = theRcut;
490  rList = rCut + 1.0;
491  
492  notifyFortranCutOffs( &rCut, &rSw, &rList );
493 }
809  
810 < void SimInfo::setDefaultRcut( double theRcut, double theRsw ){
810 >  vector<int> SimInfo::getGlobalAtomIndices() {
811 >    SimInfo::MoleculeIterator mi;
812 >    Molecule* mol;
813 >    Molecule::AtomIterator ai;
814 >    Atom* atom;
815  
816 <  rSw = theRsw;
498 <  setDefaultRcut( theRcut );
499 < }
500 <
501 <
502 < void SimInfo::checkCutOffs( void ){
503 <  
504 <  if( boxIsInit ){
816 >    vector<int> GlobalAtomIndices(getNAtoms(), 0);
817      
818 <    //we need to check cutOffs against the box
819 <    
820 <    if( rCut > maxCutoff ){
821 <      sprintf( painCave.errMsg,
822 <               "cutoffRadius is too large for the current periodic box.\n"
823 <               "\tCurrent Value of cutoffRadius = %G at time %G\n "
824 <               "\tThis is larger than half of at least one of the\n"
825 <               "\tperiodic box vectors.  Right now, the Box matrix is:\n"
514 <               "\n"
515 <               "\t[ %G %G %G ]\n"
516 <               "\t[ %G %G %G ]\n"
517 <               "\t[ %G %G %G ]\n",
518 <               rCut, currentTime,
519 <               Hmat[0][0], Hmat[0][1], Hmat[0][2],
520 <               Hmat[1][0], Hmat[1][1], Hmat[1][2],
521 <               Hmat[2][0], Hmat[2][1], Hmat[2][2]);
522 <      painCave.severity = OOPSE_ERROR;
523 <      painCave.isFatal = 1;
524 <      simError();
525 <    }    
526 <  } else {
527 <    // initialize this stuff before using it, OK?
528 <    sprintf( painCave.errMsg,
529 <             "Trying to check cutoffs without a box.\n"
530 <             "\tOOPSE should have better programmers than that.\n" );
531 <    painCave.severity = OOPSE_ERROR;
532 <    painCave.isFatal = 1;
533 <    simError();      
818 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
819 >      
820 >      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
821 >        GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex();
822 >        cerr << "LI = " << atom->getLocalIndex() << "GAI = " << GlobalAtomIndices[atom->getLocalIndex()] << "\n";
823 >      }
824 >    }
825 >    return GlobalAtomIndices;
826    }
535  
536 }
827  
538 void SimInfo::addProperty(GenericData* prop){
828  
829 <  map<string, GenericData*>::iterator result;
830 <  result = properties.find(prop->getID());
831 <  
832 <  //we can't simply use  properties[prop->getID()] = prop,
833 <  //it will cause memory leak if we already contain a propery which has the same name of prop
834 <  
835 <  if(result != properties.end()){
829 >  vector<int> SimInfo::getGlobalGroupIndices() {
830 >    SimInfo::MoleculeIterator mi;
831 >    Molecule* mol;
832 >    Molecule::CutoffGroupIterator ci;
833 >    CutoffGroup* cg;
834 >
835 >    vector<int> GlobalGroupIndices;
836      
837 <    delete (*result).second;
549 <    (*result).second = prop;
837 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
838        
839 +      //local index of cutoff group is trivial, it only depends on the
840 +      //order of travesing
841 +      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
842 +           cg = mol->nextCutoffGroup(ci)) {
843 +        GlobalGroupIndices.push_back(cg->getGlobalIndex());
844 +        cerr << "LI, GGI = " << GlobalGroupIndices.size() << " " << cg->getGlobalIndex() << "\n";
845 +      }        
846 +    }
847 +    return GlobalGroupIndices;
848    }
552  else{
849  
554    properties[prop->getID()] = prop;
850  
851 <  }
851 >  void SimInfo::prepareTopology() {
852 >    int nExclude, nOneTwo, nOneThree, nOneFour;
853 >
854 >    //calculate mass ratio of cutoff group
855 >    SimInfo::MoleculeIterator mi;
856 >    Molecule* mol;
857 >    Molecule::CutoffGroupIterator ci;
858 >    CutoffGroup* cg;
859 >    Molecule::AtomIterator ai;
860 >    Atom* atom;
861 >    RealType totalMass;
862 >
863 >    /**
864 >     * The mass factor is the relative mass of an atom to the total
865 >     * mass of the cutoff group it belongs to.  By default, all atoms
866 >     * are their own cutoff groups, and therefore have mass factors of
867 >     * 1.  We need some special handling for massless atoms, which
868 >     * will be treated as carrying the entire mass of the cutoff
869 >     * group.
870 >     */
871 >    massFactors_.clear();
872 >    massFactors_.resize(getNAtoms(), 1.0);
873      
874 < }
874 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
875 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
876 >           cg = mol->nextCutoffGroup(ci)) {
877  
878 < GenericData* SimInfo::getProperty(const string& propName){
879 <
880 <  map<string, GenericData*>::iterator result;
881 <  
882 <  //string lowerCaseName = ();
883 <  
884 <  result = properties.find(propName);
885 <  
886 <  if(result != properties.end())
887 <    return (*result).second;  
570 <  else  
571 <    return NULL;  
572 < }
878 >        totalMass = cg->getMass();
879 >        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
880 >          // Check for massless groups - set mfact to 1 if true
881 >          if (totalMass != 0)
882 >            massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass;
883 >          else
884 >            massFactors_[atom->getLocalIndex()] = 1.0;
885 >        }
886 >      }      
887 >    }
888  
889 +    // Build the identArray_
890  
891 < void SimInfo::getFortranGroupArrays(SimInfo* info,
892 <                                    vector<int>& FglobalGroupMembership,
893 <                                    vector<double>& mfact){
894 <  
895 <  Molecule* myMols;
896 <  Atom** myAtoms;
897 <  int numAtom;
898 <  double mtot;
899 <  int numMol;
584 <  int numCutoffGroups;
585 <  CutoffGroup* myCutoffGroup;
586 <  vector<CutoffGroup*>::iterator iterCutoff;
587 <  Atom* cutoffAtom;
588 <  vector<Atom*>::iterator iterAtom;
589 <  int atomIndex;
590 <  double totalMass;
591 <  
592 <  mfact.clear();
593 <  FglobalGroupMembership.clear();
594 <  
891 >    identArray_.clear();
892 >    identArray_.reserve(getNAtoms());    
893 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
894 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
895 >        identArray_.push_back(atom->getIdent());
896 >      }
897 >    }    
898 >    
899 >    //scan topology
900  
901 <  // Fix the silly fortran indexing problem
901 >    nExclude = excludedInteractions_.getSize();
902 >    nOneTwo = oneTwoInteractions_.getSize();
903 >    nOneThree = oneThreeInteractions_.getSize();
904 >    nOneFour = oneFourInteractions_.getSize();
905 >
906 >    int* excludeList = excludedInteractions_.getPairList();
907 >    int* oneTwoList = oneTwoInteractions_.getPairList();
908 >    int* oneThreeList = oneThreeInteractions_.getPairList();
909 >    int* oneFourList = oneFourInteractions_.getPairList();
910 >
911 >    topologyDone_ = true;
912 >  }
913 >
914 >  void SimInfo::addProperty(GenericData* genData) {
915 >    properties_.addProperty(genData);  
916 >  }
917 >
918 >  void SimInfo::removeProperty(const string& propName) {
919 >    properties_.removeProperty(propName);  
920 >  }
921 >
922 >  void SimInfo::clearProperties() {
923 >    properties_.clearProperties();
924 >  }
925 >
926 >  vector<string> SimInfo::getPropertyNames() {
927 >    return properties_.getPropertyNames();  
928 >  }
929 >      
930 >  vector<GenericData*> SimInfo::getProperties() {
931 >    return properties_.getProperties();
932 >  }
933 >
934 >  GenericData* SimInfo::getPropertyByName(const string& propName) {
935 >    return properties_.getPropertyByName(propName);
936 >  }
937 >
938 >  void SimInfo::setSnapshotManager(SnapshotManager* sman) {
939 >    if (sman_ == sman) {
940 >      return;
941 >    }    
942 >    delete sman_;
943 >    sman_ = sman;
944 >
945 >    Molecule* mol;
946 >    RigidBody* rb;
947 >    Atom* atom;
948 >    CutoffGroup* cg;
949 >    SimInfo::MoleculeIterator mi;
950 >    Molecule::RigidBodyIterator rbIter;
951 >    Molecule::AtomIterator atomIter;
952 >    Molecule::CutoffGroupIterator cgIter;
953 >
954 >    for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
955 >        
956 >      for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
957 >        atom->setSnapshotManager(sman_);
958 >      }
959 >        
960 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
961 >        rb->setSnapshotManager(sman_);
962 >      }
963 >
964 >      for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) {
965 >        cg->setSnapshotManager(sman_);
966 >      }
967 >    }    
968 >    
969 >  }
970 >
971 >  Vector3d SimInfo::getComVel(){
972 >    SimInfo::MoleculeIterator i;
973 >    Molecule* mol;
974 >
975 >    Vector3d comVel(0.0);
976 >    RealType totalMass = 0.0;
977 >    
978 >
979 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
980 >      RealType mass = mol->getMass();
981 >      totalMass += mass;
982 >      comVel += mass * mol->getComVel();
983 >    }  
984 >
985   #ifdef IS_MPI
986 <  numAtom = mpiSim->getNAtomsGlobal();
987 < #else
988 <  numAtom = n_atoms;
986 >    RealType tmpMass = totalMass;
987 >    Vector3d tmpComVel(comVel);    
988 >    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
989 >    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
990   #endif
602  for (int i = 0; i < numAtom; i++)
603    FglobalGroupMembership.push_back(globalGroupMembership[i] + 1);
604  
991  
992 <  myMols = info->molecules;
607 <  numMol = info->n_mol;
608 <  for(int i  = 0; i < numMol; i++){
609 <    numCutoffGroups = myMols[i].getNCutoffGroups();
610 <    for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff);
611 <        myCutoffGroup != NULL;
612 <        myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){
992 >    comVel /= totalMass;
993  
994 <      totalMass = myCutoffGroup->getMass();
994 >    return comVel;
995 >  }
996 >
997 >  Vector3d SimInfo::getCom(){
998 >    SimInfo::MoleculeIterator i;
999 >    Molecule* mol;
1000 >
1001 >    Vector3d com(0.0);
1002 >    RealType totalMass = 0.0;
1003 >    
1004 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1005 >      RealType mass = mol->getMass();
1006 >      totalMass += mass;
1007 >      com += mass * mol->getCom();
1008 >    }  
1009 >
1010 > #ifdef IS_MPI
1011 >    RealType tmpMass = totalMass;
1012 >    Vector3d tmpCom(com);    
1013 >    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1014 >    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1015 > #endif
1016 >
1017 >    com /= totalMass;
1018 >
1019 >    return com;
1020 >
1021 >  }        
1022 >
1023 >  ostream& operator <<(ostream& o, SimInfo& info) {
1024 >
1025 >    return o;
1026 >  }
1027 >  
1028 >  
1029 >   /*
1030 >   Returns center of mass and center of mass velocity in one function call.
1031 >   */
1032 >  
1033 >   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1034 >      SimInfo::MoleculeIterator i;
1035 >      Molecule* mol;
1036        
1037 <      for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom);
1038 <          cutoffAtom != NULL;
1039 <          cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){
1040 <        mfact.push_back(cutoffAtom->getMass()/totalMass);
1037 >    
1038 >      RealType totalMass = 0.0;
1039 >    
1040 >
1041 >      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1042 >         RealType mass = mol->getMass();
1043 >         totalMass += mass;
1044 >         com += mass * mol->getCom();
1045 >         comVel += mass * mol->getComVel();          
1046        }  
1047 +      
1048 + #ifdef IS_MPI
1049 +      RealType tmpMass = totalMass;
1050 +      Vector3d tmpCom(com);  
1051 +      Vector3d tmpComVel(comVel);
1052 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1053 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1054 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1055 + #endif
1056 +      
1057 +      com /= totalMass;
1058 +      comVel /= totalMass;
1059 +   }        
1060 +  
1061 +   /*
1062 +   Return intertia tensor for entire system and angular momentum Vector.
1063 +
1064 +
1065 +       [  Ixx -Ixy  -Ixz ]
1066 +    J =| -Iyx  Iyy  -Iyz |
1067 +       [ -Izx -Iyz   Izz ]
1068 +    */
1069 +
1070 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1071 +      
1072 +
1073 +      RealType xx = 0.0;
1074 +      RealType yy = 0.0;
1075 +      RealType zz = 0.0;
1076 +      RealType xy = 0.0;
1077 +      RealType xz = 0.0;
1078 +      RealType yz = 0.0;
1079 +      Vector3d com(0.0);
1080 +      Vector3d comVel(0.0);
1081 +      
1082 +      getComAll(com, comVel);
1083 +      
1084 +      SimInfo::MoleculeIterator i;
1085 +      Molecule* mol;
1086 +      
1087 +      Vector3d thisq(0.0);
1088 +      Vector3d thisv(0.0);
1089 +
1090 +      RealType thisMass = 0.0;
1091 +    
1092 +      
1093 +      
1094 +  
1095 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1096 +        
1097 +         thisq = mol->getCom()-com;
1098 +         thisv = mol->getComVel()-comVel;
1099 +         thisMass = mol->getMass();
1100 +         // Compute moment of intertia coefficients.
1101 +         xx += thisq[0]*thisq[0]*thisMass;
1102 +         yy += thisq[1]*thisq[1]*thisMass;
1103 +         zz += thisq[2]*thisq[2]*thisMass;
1104 +        
1105 +         // compute products of intertia
1106 +         xy += thisq[0]*thisq[1]*thisMass;
1107 +         xz += thisq[0]*thisq[2]*thisMass;
1108 +         yz += thisq[1]*thisq[2]*thisMass;
1109 +            
1110 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1111 +            
1112 +      }  
1113 +      
1114 +      
1115 +      inertiaTensor(0,0) = yy + zz;
1116 +      inertiaTensor(0,1) = -xy;
1117 +      inertiaTensor(0,2) = -xz;
1118 +      inertiaTensor(1,0) = -xy;
1119 +      inertiaTensor(1,1) = xx + zz;
1120 +      inertiaTensor(1,2) = -yz;
1121 +      inertiaTensor(2,0) = -xz;
1122 +      inertiaTensor(2,1) = -yz;
1123 +      inertiaTensor(2,2) = xx + yy;
1124 +      
1125 + #ifdef IS_MPI
1126 +      Mat3x3d tmpI(inertiaTensor);
1127 +      Vector3d tmpAngMom;
1128 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1129 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1130 + #endif
1131 +              
1132 +      return;
1133 +   }
1134 +
1135 +   //Returns the angular momentum of the system
1136 +   Vector3d SimInfo::getAngularMomentum(){
1137 +      
1138 +      Vector3d com(0.0);
1139 +      Vector3d comVel(0.0);
1140 +      Vector3d angularMomentum(0.0);
1141 +      
1142 +      getComAll(com,comVel);
1143 +      
1144 +      SimInfo::MoleculeIterator i;
1145 +      Molecule* mol;
1146 +      
1147 +      Vector3d thisr(0.0);
1148 +      Vector3d thisp(0.0);
1149 +      
1150 +      RealType thisMass;
1151 +      
1152 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1153 +        thisMass = mol->getMass();
1154 +        thisr = mol->getCom()-com;
1155 +        thisp = (mol->getComVel()-comVel)*thisMass;
1156 +        
1157 +        angularMomentum += cross( thisr, thisp );
1158 +        
1159 +      }  
1160 +      
1161 + #ifdef IS_MPI
1162 +      Vector3d tmpAngMom;
1163 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1164 + #endif
1165 +      
1166 +      return angularMomentum;
1167 +   }
1168 +  
1169 +  StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
1170 +    return IOIndexToIntegrableObject.at(index);
1171 +  }
1172 +  
1173 +  void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) {
1174 +    IOIndexToIntegrableObject= v;
1175 +  }
1176 +
1177 +  /* Returns the Volume of the simulation based on a ellipsoid with semi-axes
1178 +     based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3
1179 +     where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to
1180 +     V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
1181 +  */
1182 +  void SimInfo::getGyrationalVolume(RealType &volume){
1183 +    Mat3x3d intTensor;
1184 +    RealType det;
1185 +    Vector3d dummyAngMom;
1186 +    RealType sysconstants;
1187 +    RealType geomCnst;
1188 +
1189 +    geomCnst = 3.0/2.0;
1190 +    /* Get the inertial tensor and angular momentum for free*/
1191 +    getInertiaTensor(intTensor,dummyAngMom);
1192 +    
1193 +    det = intTensor.determinant();
1194 +    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1195 +    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det);
1196 +    return;
1197 +  }
1198 +
1199 +  void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){
1200 +    Mat3x3d intTensor;
1201 +    Vector3d dummyAngMom;
1202 +    RealType sysconstants;
1203 +    RealType geomCnst;
1204 +
1205 +    geomCnst = 3.0/2.0;
1206 +    /* Get the inertial tensor and angular momentum for free*/
1207 +    getInertiaTensor(intTensor,dummyAngMom);
1208 +    
1209 +    detI = intTensor.determinant();
1210 +    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1211 +    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI);
1212 +    return;
1213 +  }
1214 + /*
1215 +   void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) {
1216 +      assert( v.size() == nAtoms_ + nRigidBodies_);
1217 +      sdByGlobalIndex_ = v;
1218      }
1219 +
1220 +    StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) {
1221 +      //assert(index < nAtoms_ + nRigidBodies_);
1222 +      return sdByGlobalIndex_.at(index);
1223 +    }  
1224 + */  
1225 +  int SimInfo::getNGlobalConstraints() {
1226 +    int nGlobalConstraints;
1227 + #ifdef IS_MPI
1228 +    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
1229 +                  MPI_COMM_WORLD);    
1230 + #else
1231 +    nGlobalConstraints =  nConstraints_;
1232 + #endif
1233 +    return nGlobalConstraints;
1234    }
1235  
1236 < }
1236 > }//end namespace OpenMD
1237 >

Comparing:
trunk/src/brains/SimInfo.cpp (property svn:keywords), Revision 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
branches/development/src/brains/SimInfo.cpp (property svn:keywords), Revision 1597 by gezelter, Tue Jul 26 15:49:24 2011 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines