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 720 by chrisfen, Tue Nov 8 13:32:06 2005 UTC vs.
branches/development/src/brains/SimInfo.cpp (file contents), Revision 1547 by gezelter, Mon Apr 11 18:44:16 2011 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42   /**
# Line 48 | Line 48
48  
49   #include <algorithm>
50   #include <set>
51 + #include <map>
52  
53   #include "brains/SimInfo.hpp"
54   #include "math/Vector3.hpp"
55   #include "primitives/Molecule.hpp"
56 < #include "UseTheForce/fCutoffPolicy.h"
57 < #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
57 < #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
56 > #include "primitives/StuntDouble.hpp"
57 > #include "UseTheForce/DarkSide/neighborLists_interface.h"
58   #include "UseTheForce/doForces_interface.h"
59 #include "UseTheForce/DarkSide/electrostatic_interface.h"
60 #include "UseTheForce/notifyCutoffs_interface.h"
59   #include "utils/MemoryUtils.hpp"
60   #include "utils/simError.h"
61   #include "selection/SelectionManager.hpp"
62 + #include "io/ForceFieldOptions.hpp"
63 + #include "UseTheForce/ForceField.hpp"
64 + #include "nonbonded/SwitchingFunction.hpp"
65  
66   #ifdef IS_MPI
67   #include "UseTheForce/mpiComponentPlan.h"
68   #include "UseTheForce/DarkSide/simParallel_interface.h"
69   #endif
70  
71 < namespace oopse {
72 <
73 <  SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
74 <                   ForceField* ff, Globals* simParams) :
75 <    stamps_(stamps), forceField_(ff), simParams_(simParams),
76 <    ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
71 > using namespace std;
72 > namespace OpenMD {
73 >  
74 >  SimInfo::SimInfo(ForceField* ff, Globals* simParams) :
75 >    forceField_(ff), simParams_(simParams),
76 >    ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
77      nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
78      nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
79 <    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
80 <    nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
81 <    sman_(NULL), fortranInitialized_(false) {
82 <
82 <            
83 <      std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
84 <      MoleculeStamp* molStamp;
85 <      int nMolWithSameStamp;
86 <      int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
87 <      int nGroups = 0;      //total cutoff groups defined in meta-data file
88 <      CutoffGroupStamp* cgStamp;    
89 <      RigidBodyStamp* rbStamp;
90 <      int nRigidAtoms = 0;
79 >    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nInversions_(0),
80 >    nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0),
81 >    nConstraints_(0), sman_(NULL), fortranInitialized_(false),
82 >    calcBoxDipole_(false), useAtomicVirial_(true) {    
83      
84 <      for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
85 <        molStamp = i->first;
86 <        nMolWithSameStamp = i->second;
87 <        
88 <        addMoleculeStamp(molStamp, nMolWithSameStamp);
89 <
90 <        //calculate atoms in molecules
91 <        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
92 <
93 <
94 <        //calculate atoms in cutoff groups
95 <        int nAtomsInGroups = 0;
96 <        int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
97 <        
98 <        for (int j=0; j < nCutoffGroupsInStamp; j++) {
99 <          cgStamp = molStamp->getCutoffGroup(j);
100 <          nAtomsInGroups += cgStamp->getNMembers();
101 <        }
102 <
103 <        nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
104 <
105 <        nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
106 <
107 <        //calculate atoms in rigid bodies
108 <        int nAtomsInRigidBodies = 0;
109 <        int nRigidBodiesInStamp = molStamp->getNRigidBodies();
118 <        
119 <        for (int j=0; j < nRigidBodiesInStamp; j++) {
120 <          rbStamp = molStamp->getRigidBody(j);
121 <          nAtomsInRigidBodies += rbStamp->getNMembers();
122 <        }
123 <
124 <        nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
125 <        nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
126 <        
84 >    MoleculeStamp* molStamp;
85 >    int nMolWithSameStamp;
86 >    int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
87 >    int nGroups = 0;       //total cutoff groups defined in meta-data file
88 >    CutoffGroupStamp* cgStamp;    
89 >    RigidBodyStamp* rbStamp;
90 >    int nRigidAtoms = 0;
91 >    
92 >    vector<Component*> components = simParams->getComponents();
93 >    
94 >    for (vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
95 >      molStamp = (*i)->getMoleculeStamp();
96 >      nMolWithSameStamp = (*i)->getNMol();
97 >      
98 >      addMoleculeStamp(molStamp, nMolWithSameStamp);
99 >      
100 >      //calculate atoms in molecules
101 >      nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
102 >      
103 >      //calculate atoms in cutoff groups
104 >      int nAtomsInGroups = 0;
105 >      int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
106 >      
107 >      for (int j=0; j < nCutoffGroupsInStamp; j++) {
108 >        cgStamp = molStamp->getCutoffGroupStamp(j);
109 >        nAtomsInGroups += cgStamp->getNMembers();
110        }
111 <
112 <      //every free atom (atom does not belong to cutoff groups) is a cutoff
113 <      //group therefore the total number of cutoff groups in the system is
114 <      //equal to the total number of atoms minus number of atoms belong to
115 <      //cutoff group defined in meta-data file plus the number of cutoff
116 <      //groups defined in meta-data file
117 <      nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
118 <
119 <      //every free atom (atom does not belong to rigid bodies) is an
120 <      //integrable object therefore the total number of integrable objects
121 <      //in the system is equal to the total number of atoms minus number of
122 <      //atoms belong to rigid body defined in meta-data file plus the number
123 <      //of rigid bodies defined in meta-data file
124 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
125 <                                                + nGlobalRigidBodies_;
126 <  
127 <      nGlobalMols_ = molStampIds_.size();
145 <
146 < #ifdef IS_MPI    
147 <      molToProcMap_.resize(nGlobalMols_);
148 < #endif
149 <
111 >      
112 >      nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
113 >      
114 >      nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
115 >      
116 >      //calculate atoms in rigid bodies
117 >      int nAtomsInRigidBodies = 0;
118 >      int nRigidBodiesInStamp = molStamp->getNRigidBodies();
119 >      
120 >      for (int j=0; j < nRigidBodiesInStamp; j++) {
121 >        rbStamp = molStamp->getRigidBodyStamp(j);
122 >        nAtomsInRigidBodies += rbStamp->getNMembers();
123 >      }
124 >      
125 >      nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
126 >      nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
127 >      
128      }
129 +    
130 +    //every free atom (atom does not belong to cutoff groups) is a cutoff
131 +    //group therefore the total number of cutoff groups in the system is
132 +    //equal to the total number of atoms minus number of atoms belong to
133 +    //cutoff group defined in meta-data file plus the number of cutoff
134 +    //groups defined in meta-data file
135 +    std::cerr << "nGA = " << nGlobalAtoms_ << "\n";
136 +    std::cerr << "nCA = " << nCutoffAtoms << "\n";
137 +    std::cerr << "nG = " << nGroups << "\n";
138  
139 +    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
140 +
141 +    std::cerr << "nGCG = " << nGlobalCutoffGroups_ << "\n";
142 +    
143 +    //every free atom (atom does not belong to rigid bodies) is an
144 +    //integrable object therefore the total number of integrable objects
145 +    //in the system is equal to the total number of atoms minus number of
146 +    //atoms belong to rigid body defined in meta-data file plus the number
147 +    //of rigid bodies defined in meta-data file
148 +    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
149 +      + nGlobalRigidBodies_;
150 +    
151 +    nGlobalMols_ = molStampIds_.size();
152 +    molToProcMap_.resize(nGlobalMols_);
153 +  }
154 +  
155    SimInfo::~SimInfo() {
156 <    std::map<int, Molecule*>::iterator i;
156 >    map<int, Molecule*>::iterator i;
157      for (i = molecules_.begin(); i != molecules_.end(); ++i) {
158        delete i->second;
159      }
160      molecules_.clear();
161        
159    delete stamps_;
162      delete sman_;
163      delete simParams_;
164      delete forceField_;
165    }
166  
165  int SimInfo::getNGlobalConstraints() {
166    int nGlobalConstraints;
167 #ifdef IS_MPI
168    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
169                  MPI_COMM_WORLD);    
170 #else
171    nGlobalConstraints =  nConstraints_;
172 #endif
173    return nGlobalConstraints;
174  }
167  
168    bool SimInfo::addMolecule(Molecule* mol) {
169      MoleculeIterator i;
170 <
170 >    
171      i = molecules_.find(mol->getGlobalIndex());
172      if (i == molecules_.end() ) {
173 <
174 <      molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
175 <        
173 >      
174 >      molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
175 >      
176        nAtoms_ += mol->getNAtoms();
177        nBonds_ += mol->getNBonds();
178        nBends_ += mol->getNBends();
179        nTorsions_ += mol->getNTorsions();
180 +      nInversions_ += mol->getNInversions();
181        nRigidBodies_ += mol->getNRigidBodies();
182        nIntegrableObjects_ += mol->getNIntegrableObjects();
183        nCutoffGroups_ += mol->getNCutoffGroups();
184        nConstraints_ += mol->getNConstraintPairs();
185 <
186 <      addExcludePairs(mol);
187 <        
185 >      
186 >      addInteractionPairs(mol);
187 >      
188        return true;
189      } else {
190        return false;
191      }
192    }
193 <
193 >  
194    bool SimInfo::removeMolecule(Molecule* mol) {
195      MoleculeIterator i;
196      i = molecules_.find(mol->getGlobalIndex());
# Line 210 | Line 203 | namespace oopse {
203        nBonds_ -= mol->getNBonds();
204        nBends_ -= mol->getNBends();
205        nTorsions_ -= mol->getNTorsions();
206 +      nInversions_ -= mol->getNInversions();
207        nRigidBodies_ -= mol->getNRigidBodies();
208        nIntegrableObjects_ -= mol->getNIntegrableObjects();
209        nCutoffGroups_ -= mol->getNCutoffGroups();
210        nConstraints_ -= mol->getNConstraintPairs();
211  
212 <      removeExcludePairs(mol);
212 >      removeInteractionPairs(mol);
213        molecules_.erase(mol->getGlobalIndex());
214  
215        delete mol;
# Line 224 | Line 218 | namespace oopse {
218      } else {
219        return false;
220      }
227
228
221    }    
222  
223          
# Line 243 | Line 235 | namespace oopse {
235    void SimInfo::calcNdf() {
236      int ndf_local;
237      MoleculeIterator i;
238 <    std::vector<StuntDouble*>::iterator j;
238 >    vector<StuntDouble*>::iterator j;
239      Molecule* mol;
240      StuntDouble* integrableObject;
241  
# Line 263 | Line 255 | namespace oopse {
255            }
256          }
257              
258 <      }//end for (integrableObject)
259 <    }// end for (mol)
258 >      }
259 >    }
260      
261      // n_constraints is local, so subtract them on each processor
262      ndf_local -= nConstraints_;
# Line 281 | Line 273 | namespace oopse {
273  
274    }
275  
276 +  int SimInfo::getFdf() {
277 + #ifdef IS_MPI
278 +    MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
279 + #else
280 +    fdf_ = fdf_local;
281 + #endif
282 +    return fdf_;
283 +  }
284 +    
285    void SimInfo::calcNdfRaw() {
286      int ndfRaw_local;
287  
288      MoleculeIterator i;
289 <    std::vector<StuntDouble*>::iterator j;
289 >    vector<StuntDouble*>::iterator j;
290      Molecule* mol;
291      StuntDouble* integrableObject;
292  
# Line 332 | Line 333 | namespace oopse {
333  
334    }
335  
336 <  void SimInfo::addExcludePairs(Molecule* mol) {
337 <    std::vector<Bond*>::iterator bondIter;
338 <    std::vector<Bend*>::iterator bendIter;
339 <    std::vector<Torsion*>::iterator torsionIter;
336 >  void SimInfo::addInteractionPairs(Molecule* mol) {
337 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
338 >    vector<Bond*>::iterator bondIter;
339 >    vector<Bend*>::iterator bendIter;
340 >    vector<Torsion*>::iterator torsionIter;
341 >    vector<Inversion*>::iterator inversionIter;
342      Bond* bond;
343      Bend* bend;
344      Torsion* torsion;
345 +    Inversion* inversion;
346      int a;
347      int b;
348      int c;
349      int d;
350 +
351 +    // atomGroups can be used to add special interaction maps between
352 +    // groups of atoms that are in two separate rigid bodies.
353 +    // However, most site-site interactions between two rigid bodies
354 +    // are probably not special, just the ones between the physically
355 +    // bonded atoms.  Interactions *within* a single rigid body should
356 +    // always be excluded.  These are done at the bottom of this
357 +    // function.
358 +
359 +    map<int, set<int> > atomGroups;
360 +    Molecule::RigidBodyIterator rbIter;
361 +    RigidBody* rb;
362 +    Molecule::IntegrableObjectIterator ii;
363 +    StuntDouble* integrableObject;
364      
365 <    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
365 >    for (integrableObject = mol->beginIntegrableObject(ii);
366 >         integrableObject != NULL;
367 >         integrableObject = mol->nextIntegrableObject(ii)) {
368 >      
369 >      if (integrableObject->isRigidBody()) {
370 >        rb = static_cast<RigidBody*>(integrableObject);
371 >        vector<Atom*> atoms = rb->getAtoms();
372 >        set<int> rigidAtoms;
373 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
374 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
375 >        }
376 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
377 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
378 >        }      
379 >      } else {
380 >        set<int> oneAtomSet;
381 >        oneAtomSet.insert(integrableObject->getGlobalIndex());
382 >        atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
383 >      }
384 >    }  
385 >          
386 >    for (bond= mol->beginBond(bondIter); bond != NULL;
387 >         bond = mol->nextBond(bondIter)) {
388 >
389        a = bond->getAtomA()->getGlobalIndex();
390 <      b = bond->getAtomB()->getGlobalIndex();        
391 <      exclude_.addPair(a, b);
390 >      b = bond->getAtomB()->getGlobalIndex();  
391 >    
392 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
393 >        oneTwoInteractions_.addPair(a, b);
394 >      } else {
395 >        excludedInteractions_.addPair(a, b);
396 >      }
397      }
398  
399 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
399 >    for (bend= mol->beginBend(bendIter); bend != NULL;
400 >         bend = mol->nextBend(bendIter)) {
401 >
402        a = bend->getAtomA()->getGlobalIndex();
403        b = bend->getAtomB()->getGlobalIndex();        
404        c = bend->getAtomC()->getGlobalIndex();
405 +      
406 +      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
407 +        oneTwoInteractions_.addPair(a, b);      
408 +        oneTwoInteractions_.addPair(b, c);
409 +      } else {
410 +        excludedInteractions_.addPair(a, b);
411 +        excludedInteractions_.addPair(b, c);
412 +      }
413  
414 <      exclude_.addPair(a, b);
415 <      exclude_.addPair(a, c);
416 <      exclude_.addPair(b, c);        
414 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
415 >        oneThreeInteractions_.addPair(a, c);      
416 >      } else {
417 >        excludedInteractions_.addPair(a, c);
418 >      }
419      }
420  
421 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
421 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
422 >         torsion = mol->nextTorsion(torsionIter)) {
423 >
424        a = torsion->getAtomA()->getGlobalIndex();
425        b = torsion->getAtomB()->getGlobalIndex();        
426        c = torsion->getAtomC()->getGlobalIndex();        
427 <      d = torsion->getAtomD()->getGlobalIndex();        
427 >      d = torsion->getAtomD()->getGlobalIndex();      
428 >
429 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
430 >        oneTwoInteractions_.addPair(a, b);      
431 >        oneTwoInteractions_.addPair(b, c);
432 >        oneTwoInteractions_.addPair(c, d);
433 >      } else {
434 >        excludedInteractions_.addPair(a, b);
435 >        excludedInteractions_.addPair(b, c);
436 >        excludedInteractions_.addPair(c, d);
437 >      }
438  
439 <      exclude_.addPair(a, b);
440 <      exclude_.addPair(a, c);
441 <      exclude_.addPair(a, d);
442 <      exclude_.addPair(b, c);
443 <      exclude_.addPair(b, d);
444 <      exclude_.addPair(c, d);        
439 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
440 >        oneThreeInteractions_.addPair(a, c);      
441 >        oneThreeInteractions_.addPair(b, d);      
442 >      } else {
443 >        excludedInteractions_.addPair(a, c);
444 >        excludedInteractions_.addPair(b, d);
445 >      }
446 >
447 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
448 >        oneFourInteractions_.addPair(a, d);      
449 >      } else {
450 >        excludedInteractions_.addPair(a, d);
451 >      }
452      }
453  
454 <    Molecule::RigidBodyIterator rbIter;
455 <    RigidBody* rb;
456 <    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
457 <      std::vector<Atom*> atoms = rb->getAtoms();
458 <      for (int i = 0; i < atoms.size() -1 ; ++i) {
459 <        for (int j = i + 1; j < atoms.size(); ++j) {
454 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
455 >         inversion = mol->nextInversion(inversionIter)) {
456 >
457 >      a = inversion->getAtomA()->getGlobalIndex();
458 >      b = inversion->getAtomB()->getGlobalIndex();        
459 >      c = inversion->getAtomC()->getGlobalIndex();        
460 >      d = inversion->getAtomD()->getGlobalIndex();        
461 >
462 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
463 >        oneTwoInteractions_.addPair(a, b);      
464 >        oneTwoInteractions_.addPair(a, c);
465 >        oneTwoInteractions_.addPair(a, d);
466 >      } else {
467 >        excludedInteractions_.addPair(a, b);
468 >        excludedInteractions_.addPair(a, c);
469 >        excludedInteractions_.addPair(a, d);
470 >      }
471 >
472 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
473 >        oneThreeInteractions_.addPair(b, c);    
474 >        oneThreeInteractions_.addPair(b, d);    
475 >        oneThreeInteractions_.addPair(c, d);      
476 >      } else {
477 >        excludedInteractions_.addPair(b, c);
478 >        excludedInteractions_.addPair(b, d);
479 >        excludedInteractions_.addPair(c, d);
480 >      }
481 >    }
482 >
483 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
484 >         rb = mol->nextRigidBody(rbIter)) {
485 >      vector<Atom*> atoms = rb->getAtoms();
486 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
487 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
488            a = atoms[i]->getGlobalIndex();
489            b = atoms[j]->getGlobalIndex();
490 <          exclude_.addPair(a, b);
490 >          excludedInteractions_.addPair(a, b);
491          }
492        }
493      }        
494  
495    }
496  
497 <  void SimInfo::removeExcludePairs(Molecule* mol) {
498 <    std::vector<Bond*>::iterator bondIter;
499 <    std::vector<Bend*>::iterator bendIter;
500 <    std::vector<Torsion*>::iterator torsionIter;
497 >  void SimInfo::removeInteractionPairs(Molecule* mol) {
498 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
499 >    vector<Bond*>::iterator bondIter;
500 >    vector<Bend*>::iterator bendIter;
501 >    vector<Torsion*>::iterator torsionIter;
502 >    vector<Inversion*>::iterator inversionIter;
503      Bond* bond;
504      Bend* bend;
505      Torsion* torsion;
506 +    Inversion* inversion;
507      int a;
508      int b;
509      int c;
510      int d;
511 +
512 +    map<int, set<int> > atomGroups;
513 +    Molecule::RigidBodyIterator rbIter;
514 +    RigidBody* rb;
515 +    Molecule::IntegrableObjectIterator ii;
516 +    StuntDouble* integrableObject;
517      
518 <    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
518 >    for (integrableObject = mol->beginIntegrableObject(ii);
519 >         integrableObject != NULL;
520 >         integrableObject = mol->nextIntegrableObject(ii)) {
521 >      
522 >      if (integrableObject->isRigidBody()) {
523 >        rb = static_cast<RigidBody*>(integrableObject);
524 >        vector<Atom*> atoms = rb->getAtoms();
525 >        set<int> rigidAtoms;
526 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
527 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
528 >        }
529 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
530 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
531 >        }      
532 >      } else {
533 >        set<int> oneAtomSet;
534 >        oneAtomSet.insert(integrableObject->getGlobalIndex());
535 >        atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
536 >      }
537 >    }  
538 >
539 >    for (bond= mol->beginBond(bondIter); bond != NULL;
540 >         bond = mol->nextBond(bondIter)) {
541 >      
542        a = bond->getAtomA()->getGlobalIndex();
543 <      b = bond->getAtomB()->getGlobalIndex();        
544 <      exclude_.removePair(a, b);
543 >      b = bond->getAtomB()->getGlobalIndex();  
544 >    
545 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
546 >        oneTwoInteractions_.removePair(a, b);
547 >      } else {
548 >        excludedInteractions_.removePair(a, b);
549 >      }
550      }
551  
552 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
552 >    for (bend= mol->beginBend(bendIter); bend != NULL;
553 >         bend = mol->nextBend(bendIter)) {
554 >
555        a = bend->getAtomA()->getGlobalIndex();
556        b = bend->getAtomB()->getGlobalIndex();        
557        c = bend->getAtomC()->getGlobalIndex();
558 +      
559 +      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
560 +        oneTwoInteractions_.removePair(a, b);      
561 +        oneTwoInteractions_.removePair(b, c);
562 +      } else {
563 +        excludedInteractions_.removePair(a, b);
564 +        excludedInteractions_.removePair(b, c);
565 +      }
566  
567 <      exclude_.removePair(a, b);
568 <      exclude_.removePair(a, c);
569 <      exclude_.removePair(b, c);        
567 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
568 >        oneThreeInteractions_.removePair(a, c);      
569 >      } else {
570 >        excludedInteractions_.removePair(a, c);
571 >      }
572      }
573  
574 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
574 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
575 >         torsion = mol->nextTorsion(torsionIter)) {
576 >
577        a = torsion->getAtomA()->getGlobalIndex();
578        b = torsion->getAtomB()->getGlobalIndex();        
579        c = torsion->getAtomC()->getGlobalIndex();        
580 <      d = torsion->getAtomD()->getGlobalIndex();        
580 >      d = torsion->getAtomD()->getGlobalIndex();      
581 >  
582 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
583 >        oneTwoInteractions_.removePair(a, b);      
584 >        oneTwoInteractions_.removePair(b, c);
585 >        oneTwoInteractions_.removePair(c, d);
586 >      } else {
587 >        excludedInteractions_.removePair(a, b);
588 >        excludedInteractions_.removePair(b, c);
589 >        excludedInteractions_.removePair(c, d);
590 >      }
591  
592 <      exclude_.removePair(a, b);
593 <      exclude_.removePair(a, c);
594 <      exclude_.removePair(a, d);
595 <      exclude_.removePair(b, c);
596 <      exclude_.removePair(b, d);
597 <      exclude_.removePair(c, d);        
592 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
593 >        oneThreeInteractions_.removePair(a, c);      
594 >        oneThreeInteractions_.removePair(b, d);      
595 >      } else {
596 >        excludedInteractions_.removePair(a, c);
597 >        excludedInteractions_.removePair(b, d);
598 >      }
599 >
600 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
601 >        oneFourInteractions_.removePair(a, d);      
602 >      } else {
603 >        excludedInteractions_.removePair(a, d);
604 >      }
605      }
606  
607 <    Molecule::RigidBodyIterator rbIter;
608 <    RigidBody* rb;
609 <    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
610 <      std::vector<Atom*> atoms = rb->getAtoms();
611 <      for (int i = 0; i < atoms.size() -1 ; ++i) {
612 <        for (int j = i + 1; j < atoms.size(); ++j) {
607 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
608 >         inversion = mol->nextInversion(inversionIter)) {
609 >
610 >      a = inversion->getAtomA()->getGlobalIndex();
611 >      b = inversion->getAtomB()->getGlobalIndex();        
612 >      c = inversion->getAtomC()->getGlobalIndex();        
613 >      d = inversion->getAtomD()->getGlobalIndex();        
614 >
615 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
616 >        oneTwoInteractions_.removePair(a, b);      
617 >        oneTwoInteractions_.removePair(a, c);
618 >        oneTwoInteractions_.removePair(a, d);
619 >      } else {
620 >        excludedInteractions_.removePair(a, b);
621 >        excludedInteractions_.removePair(a, c);
622 >        excludedInteractions_.removePair(a, d);
623 >      }
624 >
625 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
626 >        oneThreeInteractions_.removePair(b, c);    
627 >        oneThreeInteractions_.removePair(b, d);    
628 >        oneThreeInteractions_.removePair(c, d);      
629 >      } else {
630 >        excludedInteractions_.removePair(b, c);
631 >        excludedInteractions_.removePair(b, d);
632 >        excludedInteractions_.removePair(c, d);
633 >      }
634 >    }
635 >
636 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
637 >         rb = mol->nextRigidBody(rbIter)) {
638 >      vector<Atom*> atoms = rb->getAtoms();
639 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
640 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
641            a = atoms[i]->getGlobalIndex();
642            b = atoms[j]->getGlobalIndex();
643 <          exclude_.removePair(a, b);
643 >          excludedInteractions_.removePair(a, b);
644          }
645        }
646      }        
647 <
647 >    
648    }
649 <
650 <
649 >  
650 >  
651    void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
652      int curStampId;
653 <
653 >    
654      //index from 0
655      curStampId = moleculeStamps_.size();
656  
# Line 457 | Line 658 | namespace oopse {
658      molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
659    }
660  
460  void SimInfo::update() {
661  
662 <    setupSimType();
663 <
664 < #ifdef IS_MPI
665 <    setupFortranParallel();
666 < #endif
667 <
668 <    setupFortranSim();
669 <
670 <    //setup fortran force field
471 <    /** @deprecate */    
472 <    int isError = 0;
473 <    
474 <    setupElectrostaticSummationMethod( isError );
475 <
476 <    if(isError){
477 <      sprintf( painCave.errMsg,
478 <               "ForceField error: There was an error initializing the forceField in fortran.\n" );
479 <      painCave.isFatal = 1;
480 <      simError();
481 <    }
482 <  
483 <    
484 <    setupCutoff();
485 <
662 >  /**
663 >   * update
664 >   *
665 >   *  Performs the global checks and variable settings after the
666 >   *  objects have been created.
667 >   *
668 >   */
669 >  void SimInfo::update() {  
670 >    setupSimVariables();
671      calcNdf();
672      calcNdfRaw();
673      calcNdfTrans();
489
490    fortranInitialized_ = true;
674    }
675 <
676 <  std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
675 >  
676 >  /**
677 >   * getSimulatedAtomTypes
678 >   *
679 >   * Returns an STL set of AtomType* that are actually present in this
680 >   * simulation.  Must query all processors to assemble this information.
681 >   *
682 >   */
683 >  set<AtomType*> SimInfo::getSimulatedAtomTypes() {
684      SimInfo::MoleculeIterator mi;
685      Molecule* mol;
686      Molecule::AtomIterator ai;
687      Atom* atom;
688 <    std::set<AtomType*> atomTypes;
689 <
690 <    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
501 <
688 >    set<AtomType*> atomTypes;
689 >    
690 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {      
691        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
692          atomTypes.insert(atom->getAtomType());
693 <      }
694 <        
506 <    }
693 >      }      
694 >    }    
695  
696 <    return atomTypes;        
697 <  }
696 > #ifdef IS_MPI
697 >
698 >    // loop over the found atom types on this processor, and add their
699 >    // numerical idents to a vector:
700 >
701 >    vector<int> foundTypes;
702 >    set<AtomType*>::iterator i;
703 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i)
704 >      foundTypes.push_back( (*i)->getIdent() );
705 >
706 >    // count_local holds the number of found types on this processor
707 >    int count_local = foundTypes.size();
708 >
709 >    // count holds the total number of found types on all processors
710 >    // (some will be redundant with the ones found locally):
711 >    int count;
712 >    MPI::COMM_WORLD.Allreduce(&count_local, &count, 1, MPI::INT, MPI::SUM);
713 >
714 >    // create a vector to hold the globally found types, and resize it:
715 >    vector<int> ftGlobal;
716 >    ftGlobal.resize(count);
717 >    vector<int> counts;
718  
719 <  void SimInfo::setupSimType() {
720 <    std::set<AtomType*>::iterator i;
721 <    std::set<AtomType*> atomTypes;
722 <    atomTypes = getUniqueAtomTypes();
719 >    int nproc = MPI::COMM_WORLD.Get_size();
720 >    counts.resize(nproc);
721 >    vector<int> disps;
722 >    disps.resize(nproc);
723 >
724 >    // now spray out the foundTypes to all the other processors:
725      
726 <    int useLennardJones = 0;
727 <    int useElectrostatic = 0;
518 <    int useEAM = 0;
519 <    int useCharge = 0;
520 <    int useDirectional = 0;
521 <    int useDipole = 0;
522 <    int useGayBerne = 0;
523 <    int useSticky = 0;
524 <    int useStickyPower = 0;
525 <    int useShape = 0;
526 <    int useFLARB = 0; //it is not in AtomType yet
527 <    int useDirectionalAtom = 0;    
528 <    int useElectrostatics = 0;
529 <    //usePBC and useRF are from simParams
530 <    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
531 <    int useRF;
532 <    int useSF;
533 <    std::string myMethod;
726 >    MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT,
727 >                               &ftGlobal[0], &counts[0], &disps[0], MPI::INT);
728  
729 <    // set the useRF logical
730 <    useRF = 0;
731 <    useSF = 0;
729 >    // foundIdents is a stl set, so inserting an already found ident
730 >    // will have no effect.
731 >    set<int> foundIdents;
732 >    vector<int>::iterator j;
733 >    for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j)
734 >      foundIdents.insert((*j));
735 >    
736 >    // now iterate over the foundIdents and get the actual atom types
737 >    // that correspond to these:
738 >    set<int>::iterator it;
739 >    for (it = foundIdents.begin(); it != foundIdents.end(); ++it)
740 >      atomTypes.insert( forceField_->getAtomType((*it)) );
741 >
742 > #endif
743 >    
744 >    return atomTypes;        
745 >  }
746  
747 <
748 <    if (simParams_->haveElectrostaticSummationMethod()) {
749 <      std::string myMethod = simParams_->getElectrostaticSummationMethod();
750 <      toUpper(myMethod);
751 <      if (myMethod == "REACTION_FIELD") {
752 <        useRF=1;
753 <      } else {
546 <        if (myMethod == "SHIFTED_FORCE") {
547 <          useSF = 1;
548 <        }
747 >  void SimInfo::setupSimVariables() {
748 >    useAtomicVirial_ = simParams_->getUseAtomicVirial();
749 >    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
750 >    calcBoxDipole_ = false;
751 >    if ( simParams_->haveAccumulateBoxDipole() )
752 >      if ( simParams_->getAccumulateBoxDipole() ) {
753 >        calcBoxDipole_ = true;      
754        }
550    }
755  
756 +    set<AtomType*>::iterator i;
757 +    set<AtomType*> atomTypes;
758 +    atomTypes = getSimulatedAtomTypes();    
759 +    int usesElectrostatic = 0;
760 +    int usesMetallic = 0;
761 +    int usesDirectional = 0;
762      //loop over all of the atom types
763      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
764 <      useLennardJones |= (*i)->isLennardJones();
765 <      useElectrostatic |= (*i)->isElectrostatic();
766 <      useEAM |= (*i)->isEAM();
557 <      useCharge |= (*i)->isCharge();
558 <      useDirectional |= (*i)->isDirectional();
559 <      useDipole |= (*i)->isDipole();
560 <      useGayBerne |= (*i)->isGayBerne();
561 <      useSticky |= (*i)->isSticky();
562 <      useStickyPower |= (*i)->isStickyPower();
563 <      useShape |= (*i)->isShape();
764 >      usesElectrostatic |= (*i)->isElectrostatic();
765 >      usesMetallic |= (*i)->isMetal();
766 >      usesDirectional |= (*i)->isDirectional();
767      }
768  
566    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
567      useDirectionalAtom = 1;
568    }
569
570    if (useCharge || useDipole) {
571      useElectrostatics = 1;
572    }
573
769   #ifdef IS_MPI    
770      int temp;
771 +    temp = usesDirectional;
772 +    MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
773  
774 <    temp = usePBC;
775 <    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
774 >    temp = usesMetallic;
775 >    MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
776  
777 <    temp = useDirectionalAtom;
778 <    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
582 <
583 <    temp = useLennardJones;
584 <    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
585 <
586 <    temp = useElectrostatics;
587 <    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
588 <
589 <    temp = useCharge;
590 <    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
591 <
592 <    temp = useDipole;
593 <    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
594 <
595 <    temp = useSticky;
596 <    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
597 <
598 <    temp = useStickyPower;
599 <    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
600 <    
601 <    temp = useGayBerne;
602 <    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
603 <
604 <    temp = useEAM;
605 <    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
606 <
607 <    temp = useShape;
608 <    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
609 <
610 <    temp = useFLARB;
611 <    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
612 <
613 <    temp = useRF;
614 <    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
615 <
616 <    temp = useSF;
617 <    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
618 <
777 >    temp = usesElectrostatic;
778 >    MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
779   #endif
780 <
781 <    fInfo_.SIM_uses_PBC = usePBC;    
782 <    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
783 <    fInfo_.SIM_uses_LennardJones = useLennardJones;
784 <    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
785 <    fInfo_.SIM_uses_Charges = useCharge;
626 <    fInfo_.SIM_uses_Dipoles = useDipole;
627 <    fInfo_.SIM_uses_Sticky = useSticky;
628 <    fInfo_.SIM_uses_StickyPower = useStickyPower;
629 <    fInfo_.SIM_uses_GayBerne = useGayBerne;
630 <    fInfo_.SIM_uses_EAM = useEAM;
631 <    fInfo_.SIM_uses_Shapes = useShape;
632 <    fInfo_.SIM_uses_FLARB = useFLARB;
633 <    fInfo_.SIM_uses_RF = useRF;
634 <    fInfo_.SIM_uses_SF = useSF;
635 <
636 <    if( myMethod == "REACTION_FIELD") {
637 <      
638 <      if (simParams_->haveDielectric()) {
639 <        fInfo_.dielect = simParams_->getDielectric();
640 <      } else {
641 <        sprintf(painCave.errMsg,
642 <                "SimSetup Error: No Dielectric constant was set.\n"
643 <                "\tYou are trying to use Reaction Field without"
644 <                "\tsetting a dielectric constant!\n");
645 <        painCave.isFatal = 1;
646 <        simError();
647 <      }      
648 <    }
649 <
780 >    fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_;    
781 >    fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_;
782 >    fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_;
783 >    fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_;
784 >    fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_;
785 >    fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_;
786    }
787  
788 <  void SimInfo::setupFortranSim() {
788 >  void SimInfo::setupFortran() {
789      int isError;
790 <    int nExclude;
791 <    std::vector<int> fortranGlobalGroupMembership;
790 >    int nExclude, nOneTwo, nOneThree, nOneFour;
791 >    vector<int> fortranGlobalGroupMembership;
792      
657    nExclude = exclude_.getSize();
793      isError = 0;
794  
795      //globalGroupMembership_ is filled by SimCreator    
# Line 663 | Line 798 | namespace oopse {
798      }
799  
800      //calculate mass ratio of cutoff group
801 <    std::vector<double> mfact;
801 >    vector<RealType> mfact;
802      SimInfo::MoleculeIterator mi;
803      Molecule* mol;
804      Molecule::CutoffGroupIterator ci;
805      CutoffGroup* cg;
806      Molecule::AtomIterator ai;
807      Atom* atom;
808 <    double totalMass;
808 >    RealType totalMass;
809  
810      //to avoid memory reallocation, reserve enough space for mfact
811      mfact.reserve(getNCutoffGroups());
# Line 686 | Line 821 | namespace oopse {
821            else
822              mfact.push_back( 1.0 );
823          }
689
824        }      
825      }
826  
827 <    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
694 <    std::vector<int> identArray;
827 >    // Build the identArray_
828  
829 <    //to avoid memory reallocation, reserve enough space identArray
830 <    identArray.reserve(getNAtoms());
698 <    
829 >    identArray_.clear();
830 >    identArray_.reserve(getNAtoms());    
831      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
832        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
833 <        identArray.push_back(atom->getIdent());
833 >        identArray_.push_back(atom->getIdent());
834        }
835      }    
836  
837      //fill molMembershipArray
838      //molMembershipArray is filled by SimCreator    
839 <    std::vector<int> molMembershipArray(nGlobalAtoms_);
839 >    vector<int> molMembershipArray(nGlobalAtoms_);
840      for (int i = 0; i < nGlobalAtoms_; i++) {
841        molMembershipArray[i] = globalMolMembership_[i] + 1;
842      }
843      
844      //setup fortran simulation
713    int nGlobalExcludes = 0;
714    int* globalExcludes = NULL;
715    int* excludeList = exclude_.getExcludeList();
716    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
717                   &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
718                   &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
845  
846 <    if( isError ){
846 >    nExclude = excludedInteractions_.getSize();
847 >    nOneTwo = oneTwoInteractions_.getSize();
848 >    nOneThree = oneThreeInteractions_.getSize();
849 >    nOneFour = oneFourInteractions_.getSize();
850  
851 +    int* excludeList = excludedInteractions_.getPairList();
852 +    int* oneTwoList = oneTwoInteractions_.getPairList();
853 +    int* oneThreeList = oneThreeInteractions_.getPairList();
854 +    int* oneFourList = oneFourInteractions_.getPairList();
855 +
856 +    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray_[0],
857 +                   &nExclude, excludeList,
858 +                   &nOneTwo, oneTwoList,
859 +                   &nOneThree, oneThreeList,
860 +                   &nOneFour, oneFourList,
861 +                   &molMembershipArray[0], &mfact[0], &nCutoffGroups_,
862 +                   &fortranGlobalGroupMembership[0], &isError);
863 +    
864 +    if( isError ){
865 +      
866        sprintf( painCave.errMsg,
867                 "There was an error setting the simulation information in fortran.\n" );
868        painCave.isFatal = 1;
869 <      painCave.severity = OOPSE_ERROR;
869 >      painCave.severity = OPENMD_ERROR;
870        simError();
871      }
872 <
873 < #ifdef IS_MPI
872 >    
873 >    
874      sprintf( checkPointMsg,
875               "succesfully sent the simulation information to fortran.\n");
732    MPIcheckPoint();
733 #endif // is_mpi
734  }
735
736
737 #ifdef IS_MPI
738  void SimInfo::setupFortranParallel() {
876      
877 <    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
878 <    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
879 <    std::vector<int> localToGlobalCutoffGroupIndex;
880 <    SimInfo::MoleculeIterator mi;
881 <    Molecule::AtomIterator ai;
882 <    Molecule::CutoffGroupIterator ci;
883 <    Molecule* mol;
884 <    Atom* atom;
885 <    CutoffGroup* cg;
877 >    errorCheckPoint();
878 >    
879 >    // Setup number of neighbors in neighbor list if present
880 >    if (simParams_->haveNeighborListNeighbors()) {
881 >      int nlistNeighbors = simParams_->getNeighborListNeighbors();
882 >      setNeighbors(&nlistNeighbors);
883 >    }
884 >  
885 > #ifdef IS_MPI    
886 >    //SimInfo is responsible for creating localToGlobalAtomIndex and
887 >    //localToGlobalGroupIndex
888 >    vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
889 >    vector<int> localToGlobalCutoffGroupIndex;
890      mpiSimData parallelData;
750    int isError;
891  
892      for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
893  
# Line 786 | Line 926 | namespace oopse {
926      }
927  
928      sprintf(checkPointMsg, " mpiRefresh successful.\n");
929 <    MPIcheckPoint();
790 <
791 <
792 <  }
793 <
929 >    errorCheckPoint();
930   #endif
931  
932 <  double SimInfo::calcMaxCutoffRadius() {
933 <
934 <
935 <    std::set<AtomType*> atomTypes;
936 <    std::set<AtomType*>::iterator i;
937 <    std::vector<double> cutoffRadius;
802 <
803 <    //get the unique atom types
804 <    atomTypes = getUniqueAtomTypes();
805 <
806 <    //query the max cutoff radius among these atom types
807 <    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
808 <      cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
932 >    initFortranFF(&isError);
933 >    if (isError) {
934 >      sprintf(painCave.errMsg,
935 >              "initFortranFF errror: fortran didn't like something we gave it.\n");
936 >      painCave.isFatal = 1;
937 >      simError();
938      }
939 <
811 <    double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
812 < #ifdef IS_MPI
813 <    //pick the max cutoff radius among the processors
814 < #endif
815 <
816 <    return maxCutoffRadius;
817 <  }
818 <
819 <  void SimInfo::getCutoff(double& rcut, double& rsw) {
820 <    
821 <    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
822 <        
823 <      if (!simParams_->haveCutoffRadius()){
824 <        sprintf(painCave.errMsg,
825 <                "SimCreator Warning: No value was set for the cutoffRadius.\n"
826 <                "\tOOPSE will use a default value of 15.0 angstroms"
827 <                "\tfor the cutoffRadius.\n");
828 <        painCave.isFatal = 0;
829 <        simError();
830 <        rcut = 15.0;
831 <      } else{
832 <        rcut = simParams_->getCutoffRadius();
833 <      }
834 <
835 <      if (!simParams_->haveSwitchingRadius()){
836 <        sprintf(painCave.errMsg,
837 <                "SimCreator Warning: No value was set for switchingRadius.\n"
838 <                "\tOOPSE will use a default value of\n"
839 <                "\t0.85 * cutoffRadius for the switchingRadius\n");
840 <        painCave.isFatal = 0;
841 <        simError();
842 <        rsw = 0.85 * rcut;
843 <      } else{
844 <        rsw = simParams_->getSwitchingRadius();
845 <      }
846 <
847 <    } else {
848 <      // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
849 <      //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
850 <        
851 <      if (simParams_->haveCutoffRadius()) {
852 <        rcut = simParams_->getCutoffRadius();
853 <      } else {
854 <        //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
855 <        rcut = calcMaxCutoffRadius();
856 <      }
857 <
858 <      if (simParams_->haveSwitchingRadius()) {
859 <        rsw  = simParams_->getSwitchingRadius();
860 <      } else {
861 <        rsw = rcut;
862 <      }
863 <    
864 <    }
939 >    fortranInitialized_ = true;
940    }
941  
867  void SimInfo::setupCutoff() {    
868    getCutoff(rcut_, rsw_);    
869    double rnblist = rcut_ + 1; // skin of neighbor list
870
871    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
872    
873    int cp =  TRADITIONAL_CUTOFF_POLICY;
874    if (simParams_->haveCutoffPolicy()) {
875      std::string myPolicy = simParams_->getCutoffPolicy();
876      toUpper(myPolicy);
877      if (myPolicy == "MIX") {
878        cp = MIX_CUTOFF_POLICY;
879      } else {
880        if (myPolicy == "MAX") {
881          cp = MAX_CUTOFF_POLICY;
882        } else {
883          if (myPolicy == "TRADITIONAL") {            
884            cp = TRADITIONAL_CUTOFF_POLICY;
885          } else {
886            // throw error        
887            sprintf( painCave.errMsg,
888                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
889            painCave.isFatal = 1;
890            simError();
891          }    
892        }          
893      }
894    }
895
896
897    if (simParams_->haveSkinThickness()) {
898      double skinThickness = simParams_->getSkinThickness();
899    }
900
901    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
902    // also send cutoff notification to electrostatics
903    setElectrostaticCutoffRadius(&rcut_, &rsw_);
904  }
905
906  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
907    
908    int errorOut;
909    int esm =  NONE;
910    int sm = UNDAMPED;
911    double alphaVal;
912    double dielectric;
913
914    errorOut = isError;
915    alphaVal = simParams_->getDampingAlpha();
916    dielectric = simParams_->getDielectric();
917
918    if (simParams_->haveElectrostaticSummationMethod()) {
919      std::string myMethod = simParams_->getElectrostaticSummationMethod();
920      toUpper(myMethod);
921      if (myMethod == "NONE") {
922        esm = NONE;
923      } else {
924        if (myMethod == "SWITCHING_FUNCTION") {
925          esm = SWITCHING_FUNCTION;
926        } else {
927          if (myMethod == "SHIFTED_POTENTIAL") {
928            esm = SHIFTED_POTENTIAL;
929          } else {
930            if (myMethod == "SHIFTED_FORCE") {            
931              esm = SHIFTED_FORCE;
932            } else {
933              if (myMethod == "REACTION_FIELD") {            
934                esm = REACTION_FIELD;
935              } else {
936                // throw error        
937                sprintf( painCave.errMsg,
938                         "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() );
939                painCave.isFatal = 1;
940                simError();
941              }    
942            }          
943          }
944        }
945      }
946    }
947    
948    if (simParams_->haveElectrostaticScreeningMethod()) {
949      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
950      toUpper(myScreen);
951      if (myScreen == "UNDAMPED") {
952        sm = UNDAMPED;
953      } else {
954        if (myScreen == "DAMPED") {
955          sm = DAMPED;
956          if (!simParams_->haveDampingAlpha()) {
957            //throw error
958            sprintf( painCave.errMsg,
959                     "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal);
960            painCave.isFatal = 0;
961            simError();
962          }
963        } else {
964          // throw error        
965          sprintf( painCave.errMsg,
966                   "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() );
967          painCave.isFatal = 1;
968          simError();
969        }
970      }
971    }
972    
973    // let's pass some summation method variables to fortran
974    setElectrostaticSummationMethod( &esm );
975    setScreeningMethod( &sm );
976    setDampingAlpha( &alphaVal );
977    setReactionFieldDielectric( &dielectric );
978    initFortranFF( &esm, &errorOut );
979  }
980
942    void SimInfo::addProperty(GenericData* genData) {
943      properties_.addProperty(genData);  
944    }
945  
946 <  void SimInfo::removeProperty(const std::string& propName) {
946 >  void SimInfo::removeProperty(const string& propName) {
947      properties_.removeProperty(propName);  
948    }
949  
# Line 990 | Line 951 | namespace oopse {
951      properties_.clearProperties();
952    }
953  
954 <  std::vector<std::string> SimInfo::getPropertyNames() {
954 >  vector<string> SimInfo::getPropertyNames() {
955      return properties_.getPropertyNames();  
956    }
957        
958 <  std::vector<GenericData*> SimInfo::getProperties() {
958 >  vector<GenericData*> SimInfo::getProperties() {
959      return properties_.getProperties();
960    }
961  
962 <  GenericData* SimInfo::getPropertyByName(const std::string& propName) {
962 >  GenericData* SimInfo::getPropertyByName(const string& propName) {
963      return properties_.getPropertyByName(propName);
964    }
965  
# Line 1012 | Line 973 | namespace oopse {
973      Molecule* mol;
974      RigidBody* rb;
975      Atom* atom;
976 +    CutoffGroup* cg;
977      SimInfo::MoleculeIterator mi;
978      Molecule::RigidBodyIterator rbIter;
979 <    Molecule::AtomIterator atomIter;;
979 >    Molecule::AtomIterator atomIter;
980 >    Molecule::CutoffGroupIterator cgIter;
981  
982      for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
983          
# Line 1025 | Line 988 | namespace oopse {
988        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
989          rb->setSnapshotManager(sman_);
990        }
991 +
992 +      for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) {
993 +        cg->setSnapshotManager(sman_);
994 +      }
995      }    
996      
997    }
# Line 1034 | Line 1001 | namespace oopse {
1001      Molecule* mol;
1002  
1003      Vector3d comVel(0.0);
1004 <    double totalMass = 0.0;
1004 >    RealType totalMass = 0.0;
1005      
1006  
1007      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1008 <      double mass = mol->getMass();
1008 >      RealType mass = mol->getMass();
1009        totalMass += mass;
1010        comVel += mass * mol->getComVel();
1011      }  
1012  
1013   #ifdef IS_MPI
1014 <    double tmpMass = totalMass;
1014 >    RealType tmpMass = totalMass;
1015      Vector3d tmpComVel(comVel);    
1016 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1017 <    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1016 >    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1017 >    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1018   #endif
1019  
1020      comVel /= totalMass;
# Line 1060 | Line 1027 | namespace oopse {
1027      Molecule* mol;
1028  
1029      Vector3d com(0.0);
1030 <    double totalMass = 0.0;
1030 >    RealType totalMass = 0.0;
1031      
1032      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1033 <      double mass = mol->getMass();
1033 >      RealType mass = mol->getMass();
1034        totalMass += mass;
1035        com += mass * mol->getCom();
1036      }  
1037  
1038   #ifdef IS_MPI
1039 <    double tmpMass = totalMass;
1039 >    RealType tmpMass = totalMass;
1040      Vector3d tmpCom(com);    
1041 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1042 <    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
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   #endif
1044  
1045      com /= totalMass;
# Line 1081 | Line 1048 | namespace oopse {
1048  
1049    }        
1050  
1051 <  std::ostream& operator <<(std::ostream& o, SimInfo& info) {
1051 >  ostream& operator <<(ostream& o, SimInfo& info) {
1052  
1053      return o;
1054    }
# Line 1096 | Line 1063 | namespace oopse {
1063        Molecule* mol;
1064        
1065      
1066 <      double totalMass = 0.0;
1066 >      RealType totalMass = 0.0;
1067      
1068  
1069        for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1070 <         double mass = mol->getMass();
1070 >         RealType mass = mol->getMass();
1071           totalMass += mass;
1072           com += mass * mol->getCom();
1073           comVel += mass * mol->getComVel();          
1074        }  
1075        
1076   #ifdef IS_MPI
1077 <      double tmpMass = totalMass;
1077 >      RealType tmpMass = totalMass;
1078        Vector3d tmpCom(com);  
1079        Vector3d tmpComVel(comVel);
1080 <      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1081 <      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1082 <      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1080 >      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1081 >      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1082 >      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1083   #endif
1084        
1085        com /= totalMass;
# Line 1124 | Line 1091 | namespace oopse {
1091  
1092  
1093         [  Ixx -Ixy  -Ixz ]
1094 <  J =| -Iyx  Iyy  -Iyz |
1094 >    J =| -Iyx  Iyy  -Iyz |
1095         [ -Izx -Iyz   Izz ]
1096      */
1097  
1098     void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1099        
1100  
1101 <      double xx = 0.0;
1102 <      double yy = 0.0;
1103 <      double zz = 0.0;
1104 <      double xy = 0.0;
1105 <      double xz = 0.0;
1106 <      double yz = 0.0;
1101 >      RealType xx = 0.0;
1102 >      RealType yy = 0.0;
1103 >      RealType zz = 0.0;
1104 >      RealType xy = 0.0;
1105 >      RealType xz = 0.0;
1106 >      RealType yz = 0.0;
1107        Vector3d com(0.0);
1108        Vector3d comVel(0.0);
1109        
# Line 1148 | Line 1115 | namespace oopse {
1115        Vector3d thisq(0.0);
1116        Vector3d thisv(0.0);
1117  
1118 <      double thisMass = 0.0;
1118 >      RealType thisMass = 0.0;
1119      
1120        
1121        
# Line 1186 | Line 1153 | namespace oopse {
1153   #ifdef IS_MPI
1154        Mat3x3d tmpI(inertiaTensor);
1155        Vector3d tmpAngMom;
1156 <      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1157 <      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1156 >      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1157 >      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1158   #endif
1159                
1160        return;
# Line 1208 | Line 1175 | namespace oopse {
1175        Vector3d thisr(0.0);
1176        Vector3d thisp(0.0);
1177        
1178 <      double thisMass;
1178 >      RealType thisMass;
1179        
1180        for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1181          thisMass = mol->getMass();
# Line 1221 | Line 1188 | namespace oopse {
1188        
1189   #ifdef IS_MPI
1190        Vector3d tmpAngMom;
1191 <      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1191 >      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1192   #endif
1193        
1194        return angularMomentum;
1195     }
1196    
1197 <  
1198 < }//end namespace oopse
1197 >  StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
1198 >    return IOIndexToIntegrableObject.at(index);
1199 >  }
1200 >  
1201 >  void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) {
1202 >    IOIndexToIntegrableObject= v;
1203 >  }
1204  
1205 +  /* Returns the Volume of the simulation based on a ellipsoid with semi-axes
1206 +     based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3
1207 +     where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to
1208 +     V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
1209 +  */
1210 +  void SimInfo::getGyrationalVolume(RealType &volume){
1211 +    Mat3x3d intTensor;
1212 +    RealType det;
1213 +    Vector3d dummyAngMom;
1214 +    RealType sysconstants;
1215 +    RealType geomCnst;
1216 +
1217 +    geomCnst = 3.0/2.0;
1218 +    /* Get the inertial tensor and angular momentum for free*/
1219 +    getInertiaTensor(intTensor,dummyAngMom);
1220 +    
1221 +    det = intTensor.determinant();
1222 +    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1223 +    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det);
1224 +    return;
1225 +  }
1226 +
1227 +  void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){
1228 +    Mat3x3d intTensor;
1229 +    Vector3d dummyAngMom;
1230 +    RealType sysconstants;
1231 +    RealType geomCnst;
1232 +
1233 +    geomCnst = 3.0/2.0;
1234 +    /* Get the inertial tensor and angular momentum for free*/
1235 +    getInertiaTensor(intTensor,dummyAngMom);
1236 +    
1237 +    detI = intTensor.determinant();
1238 +    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1239 +    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI);
1240 +    return;
1241 +  }
1242 + /*
1243 +   void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) {
1244 +      assert( v.size() == nAtoms_ + nRigidBodies_);
1245 +      sdByGlobalIndex_ = v;
1246 +    }
1247 +
1248 +    StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) {
1249 +      //assert(index < nAtoms_ + nRigidBodies_);
1250 +      return sdByGlobalIndex_.at(index);
1251 +    }  
1252 + */  
1253 +  int SimInfo::getNGlobalConstraints() {
1254 +    int nGlobalConstraints;
1255 + #ifdef IS_MPI
1256 +    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
1257 +                  MPI_COMM_WORLD);    
1258 + #else
1259 +    nGlobalConstraints =  nConstraints_;
1260 + #endif
1261 +    return nGlobalConstraints;
1262 +  }
1263 +
1264 + }//end namespace OpenMD
1265 +

Comparing:
trunk/src/brains/SimInfo.cpp (property svn:keywords), Revision 720 by chrisfen, Tue Nov 8 13:32:06 2005 UTC vs.
branches/development/src/brains/SimInfo.cpp (property svn:keywords), Revision 1547 by gezelter, Mon Apr 11 18:44:16 2011 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines