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 1121 by chuckv, Mon Feb 26 04:45:42 2007 UTC vs.
branches/development/src/brains/SimInfo.cpp (file contents), Revision 1534 by gezelter, Wed Dec 29 21:53:28 2010 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 54 | Line 54
54   #include "math/Vector3.hpp"
55   #include "primitives/Molecule.hpp"
56   #include "primitives/StuntDouble.hpp"
57 #include "UseTheForce/fCutoffPolicy.h"
58 #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
59 #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
60 #include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
57   #include "UseTheForce/doForces_interface.h"
58   #include "UseTheForce/DarkSide/neighborLists_interface.h"
63 #include "UseTheForce/DarkSide/electrostatic_interface.h"
64 #include "UseTheForce/DarkSide/switcheroo_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  
71
66   #ifdef IS_MPI
67   #include "UseTheForce/mpiComponentPlan.h"
68   #include "UseTheForce/DarkSide/simParallel_interface.h"
69   #endif
70  
71 < namespace oopse {
72 <  std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) {
79 <    std::map<int, std::set<int> >::iterator i = container.find(index);
80 <    std::set<int> result;
81 <    if (i != container.end()) {
82 <        result = i->second;
83 <    }
84 <
85 <    return result;
86 <  }
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), calcBoxDipole_(false) {
82 <
83 <      MoleculeStamp* molStamp;
84 <      int nMolWithSameStamp;
85 <      int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
86 <      int nGroups = 0;      //total cutoff groups defined in meta-data file
87 <      CutoffGroupStamp* cgStamp;    
88 <      RigidBodyStamp* rbStamp;
89 <      int nRigidAtoms = 0;
90 <      std::vector<Component*> components = simParams->getComponents();
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 >    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 <      for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
99 <        molStamp = (*i)->getMoleculeStamp();
100 <        nMolWithSameStamp = (*i)->getNMol();
101 <        
102 <        addMoleculeStamp(molStamp, nMolWithSameStamp);
103 <
104 <        //calculate atoms in molecules
105 <        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
106 <
107 <        //calculate atoms in cutoff groups
108 <        int nAtomsInGroups = 0;
109 <        int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
118 <        
119 <        for (int j=0; j < nCutoffGroupsInStamp; j++) {
120 <          cgStamp = molStamp->getCutoffGroupStamp(j);
121 <          nAtomsInGroups += cgStamp->getNMembers();
122 <        }
123 <
124 <        nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
125 <
126 <        nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
127 <
128 <        //calculate atoms in rigid bodies
129 <        int nAtomsInRigidBodies = 0;
130 <        int nRigidBodiesInStamp = molStamp->getNRigidBodies();
131 <        
132 <        for (int j=0; j < nRigidBodiesInStamp; j++) {
133 <          rbStamp = molStamp->getRigidBodyStamp(j);
134 <          nAtomsInRigidBodies += rbStamp->getNMembers();
135 <        }
136 <
137 <        nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
138 <        nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
139 <        
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();
158 <
159 < #ifdef IS_MPI    
160 <      molToProcMap_.resize(nGlobalMols_);
161 < #endif
162 <
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 <
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 >    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
136 >    
137 >    //every free atom (atom does not belong to rigid bodies) is an
138 >    //integrable object therefore the total number of integrable objects
139 >    //in the system is equal to the total number of atoms minus number of
140 >    //atoms belong to rigid body defined in meta-data file plus the number
141 >    //of rigid bodies defined in meta-data file
142 >    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
143 >      + nGlobalRigidBodies_;
144 >    
145 >    nGlobalMols_ = molStampIds_.size();
146 >    molToProcMap_.resize(nGlobalMols_);
147 >  }
148 >  
149    SimInfo::~SimInfo() {
150 <    std::map<int, Molecule*>::iterator i;
150 >    map<int, Molecule*>::iterator i;
151      for (i = molecules_.begin(); i != molecules_.end(); ++i) {
152        delete i->second;
153      }
# Line 174 | Line 158 | namespace oopse {
158      delete forceField_;
159    }
160  
177  int SimInfo::getNGlobalConstraints() {
178    int nGlobalConstraints;
179 #ifdef IS_MPI
180    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
181                  MPI_COMM_WORLD);    
182 #else
183    nGlobalConstraints =  nConstraints_;
184 #endif
185    return nGlobalConstraints;
186  }
161  
162    bool SimInfo::addMolecule(Molecule* mol) {
163      MoleculeIterator i;
164 <
164 >    
165      i = molecules_.find(mol->getGlobalIndex());
166      if (i == molecules_.end() ) {
167 <
168 <      molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
169 <        
167 >      
168 >      molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
169 >      
170        nAtoms_ += mol->getNAtoms();
171        nBonds_ += mol->getNBonds();
172        nBends_ += mol->getNBends();
173        nTorsions_ += mol->getNTorsions();
174 +      nInversions_ += mol->getNInversions();
175        nRigidBodies_ += mol->getNRigidBodies();
176        nIntegrableObjects_ += mol->getNIntegrableObjects();
177        nCutoffGroups_ += mol->getNCutoffGroups();
178        nConstraints_ += mol->getNConstraintPairs();
179 <
180 <      addExcludePairs(mol);
181 <        
179 >      
180 >      addInteractionPairs(mol);
181 >      
182        return true;
183      } else {
184        return false;
185      }
186    }
187 <
187 >  
188    bool SimInfo::removeMolecule(Molecule* mol) {
189      MoleculeIterator i;
190      i = molecules_.find(mol->getGlobalIndex());
# Line 222 | Line 197 | namespace oopse {
197        nBonds_ -= mol->getNBonds();
198        nBends_ -= mol->getNBends();
199        nTorsions_ -= mol->getNTorsions();
200 +      nInversions_ -= mol->getNInversions();
201        nRigidBodies_ -= mol->getNRigidBodies();
202        nIntegrableObjects_ -= mol->getNIntegrableObjects();
203        nCutoffGroups_ -= mol->getNCutoffGroups();
204        nConstraints_ -= mol->getNConstraintPairs();
205  
206 <      removeExcludePairs(mol);
206 >      removeInteractionPairs(mol);
207        molecules_.erase(mol->getGlobalIndex());
208  
209        delete mol;
# Line 236 | Line 212 | namespace oopse {
212      } else {
213        return false;
214      }
239
240
215    }    
216  
217          
# Line 255 | Line 229 | namespace oopse {
229    void SimInfo::calcNdf() {
230      int ndf_local;
231      MoleculeIterator i;
232 <    std::vector<StuntDouble*>::iterator j;
232 >    vector<StuntDouble*>::iterator j;
233      Molecule* mol;
234      StuntDouble* integrableObject;
235  
# Line 306 | Line 280 | namespace oopse {
280      int ndfRaw_local;
281  
282      MoleculeIterator i;
283 <    std::vector<StuntDouble*>::iterator j;
283 >    vector<StuntDouble*>::iterator j;
284      Molecule* mol;
285      StuntDouble* integrableObject;
286  
# Line 353 | Line 327 | namespace oopse {
327  
328    }
329  
330 <  void SimInfo::addExcludePairs(Molecule* mol) {
331 <    std::vector<Bond*>::iterator bondIter;
332 <    std::vector<Bend*>::iterator bendIter;
333 <    std::vector<Torsion*>::iterator torsionIter;
330 >  void SimInfo::addInteractionPairs(Molecule* mol) {
331 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
332 >    vector<Bond*>::iterator bondIter;
333 >    vector<Bend*>::iterator bendIter;
334 >    vector<Torsion*>::iterator torsionIter;
335 >    vector<Inversion*>::iterator inversionIter;
336      Bond* bond;
337      Bend* bend;
338      Torsion* torsion;
339 +    Inversion* inversion;
340      int a;
341      int b;
342      int c;
343      int d;
344  
345 <    std::map<int, std::set<int> > atomGroups;
345 >    // atomGroups can be used to add special interaction maps between
346 >    // groups of atoms that are in two separate rigid bodies.
347 >    // However, most site-site interactions between two rigid bodies
348 >    // are probably not special, just the ones between the physically
349 >    // bonded atoms.  Interactions *within* a single rigid body should
350 >    // always be excluded.  These are done at the bottom of this
351 >    // function.
352  
353 +    map<int, set<int> > atomGroups;
354      Molecule::RigidBodyIterator rbIter;
355      RigidBody* rb;
356      Molecule::IntegrableObjectIterator ii;
357      StuntDouble* integrableObject;
358      
359 <    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
360 <           integrableObject = mol->nextIntegrableObject(ii)) {
361 <
359 >    for (integrableObject = mol->beginIntegrableObject(ii);
360 >         integrableObject != NULL;
361 >         integrableObject = mol->nextIntegrableObject(ii)) {
362 >      
363        if (integrableObject->isRigidBody()) {
364 <          rb = static_cast<RigidBody*>(integrableObject);
365 <          std::vector<Atom*> atoms = rb->getAtoms();
366 <          std::set<int> rigidAtoms;
367 <          for (int i = 0; i < atoms.size(); ++i) {
368 <            rigidAtoms.insert(atoms[i]->getGlobalIndex());
369 <          }
370 <          for (int i = 0; i < atoms.size(); ++i) {
371 <            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
372 <          }      
364 >        rb = static_cast<RigidBody*>(integrableObject);
365 >        vector<Atom*> atoms = rb->getAtoms();
366 >        set<int> rigidAtoms;
367 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
368 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
369 >        }
370 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
371 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
372 >        }      
373        } else {
374 <        std::set<int> oneAtomSet;
374 >        set<int> oneAtomSet;
375          oneAtomSet.insert(integrableObject->getGlobalIndex());
376 <        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
376 >        atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
377        }
378      }  
379 +          
380 +    for (bond= mol->beginBond(bondIter); bond != NULL;
381 +         bond = mol->nextBond(bondIter)) {
382  
395    
396    
397    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
383        a = bond->getAtomA()->getGlobalIndex();
384 <      b = bond->getAtomB()->getGlobalIndex();        
385 <      exclude_.addPair(a, b);
384 >      b = bond->getAtomB()->getGlobalIndex();  
385 >    
386 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
387 >        oneTwoInteractions_.addPair(a, b);
388 >      } else {
389 >        excludedInteractions_.addPair(a, b);
390 >      }
391      }
392  
393 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
393 >    for (bend= mol->beginBend(bendIter); bend != NULL;
394 >         bend = mol->nextBend(bendIter)) {
395 >
396        a = bend->getAtomA()->getGlobalIndex();
397        b = bend->getAtomB()->getGlobalIndex();        
398        c = bend->getAtomC()->getGlobalIndex();
407      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
408      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
409      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
410
411      exclude_.addPairs(rigidSetA, rigidSetB);
412      exclude_.addPairs(rigidSetA, rigidSetC);
413      exclude_.addPairs(rigidSetB, rigidSetC);
399        
400 <      //exclude_.addPair(a, b);
401 <      //exclude_.addPair(a, c);
402 <      //exclude_.addPair(b, c);        
400 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
401 >        oneTwoInteractions_.addPair(a, b);      
402 >        oneTwoInteractions_.addPair(b, c);
403 >      } else {
404 >        excludedInteractions_.addPair(a, b);
405 >        excludedInteractions_.addPair(b, c);
406 >      }
407 >
408 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
409 >        oneThreeInteractions_.addPair(a, c);      
410 >      } else {
411 >        excludedInteractions_.addPair(a, c);
412 >      }
413      }
414  
415 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
415 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
416 >         torsion = mol->nextTorsion(torsionIter)) {
417 >
418        a = torsion->getAtomA()->getGlobalIndex();
419        b = torsion->getAtomB()->getGlobalIndex();        
420        c = torsion->getAtomC()->getGlobalIndex();        
421 <      d = torsion->getAtomD()->getGlobalIndex();        
425 <      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
426 <      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
427 <      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
428 <      std::set<int> rigidSetD = getRigidSet(d, atomGroups);
421 >      d = torsion->getAtomD()->getGlobalIndex();      
422  
423 <      exclude_.addPairs(rigidSetA, rigidSetB);
424 <      exclude_.addPairs(rigidSetA, rigidSetC);
425 <      exclude_.addPairs(rigidSetA, rigidSetD);
426 <      exclude_.addPairs(rigidSetB, rigidSetC);
427 <      exclude_.addPairs(rigidSetB, rigidSetD);
428 <      exclude_.addPairs(rigidSetC, rigidSetD);
423 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
424 >        oneTwoInteractions_.addPair(a, b);      
425 >        oneTwoInteractions_.addPair(b, c);
426 >        oneTwoInteractions_.addPair(c, d);
427 >      } else {
428 >        excludedInteractions_.addPair(a, b);
429 >        excludedInteractions_.addPair(b, c);
430 >        excludedInteractions_.addPair(c, d);
431 >      }
432  
433 <      /*
434 <      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
435 <      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
436 <      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
437 <      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
438 <      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
439 <      exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
440 <        
441 <      
442 <      exclude_.addPair(a, b);
443 <      exclude_.addPair(a, c);
444 <      exclude_.addPair(a, d);
445 <      exclude_.addPair(b, c);
450 <      exclude_.addPair(b, d);
451 <      exclude_.addPair(c, d);        
452 <      */
433 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
434 >        oneThreeInteractions_.addPair(a, c);      
435 >        oneThreeInteractions_.addPair(b, d);      
436 >      } else {
437 >        excludedInteractions_.addPair(a, c);
438 >        excludedInteractions_.addPair(b, d);
439 >      }
440 >
441 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
442 >        oneFourInteractions_.addPair(a, d);      
443 >      } else {
444 >        excludedInteractions_.addPair(a, d);
445 >      }
446      }
447  
448 <    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
449 <      std::vector<Atom*> atoms = rb->getAtoms();
450 <      for (int i = 0; i < atoms.size() -1 ; ++i) {
451 <        for (int j = i + 1; j < atoms.size(); ++j) {
448 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
449 >         inversion = mol->nextInversion(inversionIter)) {
450 >
451 >      a = inversion->getAtomA()->getGlobalIndex();
452 >      b = inversion->getAtomB()->getGlobalIndex();        
453 >      c = inversion->getAtomC()->getGlobalIndex();        
454 >      d = inversion->getAtomD()->getGlobalIndex();        
455 >
456 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
457 >        oneTwoInteractions_.addPair(a, b);      
458 >        oneTwoInteractions_.addPair(a, c);
459 >        oneTwoInteractions_.addPair(a, d);
460 >      } else {
461 >        excludedInteractions_.addPair(a, b);
462 >        excludedInteractions_.addPair(a, c);
463 >        excludedInteractions_.addPair(a, d);
464 >      }
465 >
466 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
467 >        oneThreeInteractions_.addPair(b, c);    
468 >        oneThreeInteractions_.addPair(b, d);    
469 >        oneThreeInteractions_.addPair(c, d);      
470 >      } else {
471 >        excludedInteractions_.addPair(b, c);
472 >        excludedInteractions_.addPair(b, d);
473 >        excludedInteractions_.addPair(c, d);
474 >      }
475 >    }
476 >
477 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
478 >         rb = mol->nextRigidBody(rbIter)) {
479 >      vector<Atom*> atoms = rb->getAtoms();
480 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
481 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
482            a = atoms[i]->getGlobalIndex();
483            b = atoms[j]->getGlobalIndex();
484 <          exclude_.addPair(a, b);
484 >          excludedInteractions_.addPair(a, b);
485          }
486        }
487      }        
488  
489    }
490  
491 <  void SimInfo::removeExcludePairs(Molecule* mol) {
492 <    std::vector<Bond*>::iterator bondIter;
493 <    std::vector<Bend*>::iterator bendIter;
494 <    std::vector<Torsion*>::iterator torsionIter;
491 >  void SimInfo::removeInteractionPairs(Molecule* mol) {
492 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
493 >    vector<Bond*>::iterator bondIter;
494 >    vector<Bend*>::iterator bendIter;
495 >    vector<Torsion*>::iterator torsionIter;
496 >    vector<Inversion*>::iterator inversionIter;
497      Bond* bond;
498      Bend* bend;
499      Torsion* torsion;
500 +    Inversion* inversion;
501      int a;
502      int b;
503      int c;
504      int d;
505  
506 <    std::map<int, std::set<int> > atomGroups;
481 <
506 >    map<int, set<int> > atomGroups;
507      Molecule::RigidBodyIterator rbIter;
508      RigidBody* rb;
509      Molecule::IntegrableObjectIterator ii;
510      StuntDouble* integrableObject;
511      
512 <    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
513 <           integrableObject = mol->nextIntegrableObject(ii)) {
514 <
512 >    for (integrableObject = mol->beginIntegrableObject(ii);
513 >         integrableObject != NULL;
514 >         integrableObject = mol->nextIntegrableObject(ii)) {
515 >      
516        if (integrableObject->isRigidBody()) {
517 <          rb = static_cast<RigidBody*>(integrableObject);
518 <          std::vector<Atom*> atoms = rb->getAtoms();
519 <          std::set<int> rigidAtoms;
520 <          for (int i = 0; i < atoms.size(); ++i) {
521 <            rigidAtoms.insert(atoms[i]->getGlobalIndex());
522 <          }
523 <          for (int i = 0; i < atoms.size(); ++i) {
524 <            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
525 <          }      
517 >        rb = static_cast<RigidBody*>(integrableObject);
518 >        vector<Atom*> atoms = rb->getAtoms();
519 >        set<int> rigidAtoms;
520 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
521 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
522 >        }
523 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
524 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
525 >        }      
526        } else {
527 <        std::set<int> oneAtomSet;
527 >        set<int> oneAtomSet;
528          oneAtomSet.insert(integrableObject->getGlobalIndex());
529 <        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
529 >        atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
530        }
531      }  
532  
533 <    
534 <    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
533 >    for (bond= mol->beginBond(bondIter); bond != NULL;
534 >         bond = mol->nextBond(bondIter)) {
535 >      
536        a = bond->getAtomA()->getGlobalIndex();
537 <      b = bond->getAtomB()->getGlobalIndex();        
538 <      exclude_.removePair(a, b);
537 >      b = bond->getAtomB()->getGlobalIndex();  
538 >    
539 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
540 >        oneTwoInteractions_.removePair(a, b);
541 >      } else {
542 >        excludedInteractions_.removePair(a, b);
543 >      }
544      }
545  
546 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
546 >    for (bend= mol->beginBend(bendIter); bend != NULL;
547 >         bend = mol->nextBend(bendIter)) {
548 >
549        a = bend->getAtomA()->getGlobalIndex();
550        b = bend->getAtomB()->getGlobalIndex();        
551        c = bend->getAtomC()->getGlobalIndex();
518
519      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
520      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
521      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
522
523      exclude_.removePairs(rigidSetA, rigidSetB);
524      exclude_.removePairs(rigidSetA, rigidSetC);
525      exclude_.removePairs(rigidSetB, rigidSetC);
552        
553 <      //exclude_.removePair(a, b);
554 <      //exclude_.removePair(a, c);
555 <      //exclude_.removePair(b, c);        
553 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
554 >        oneTwoInteractions_.removePair(a, b);      
555 >        oneTwoInteractions_.removePair(b, c);
556 >      } else {
557 >        excludedInteractions_.removePair(a, b);
558 >        excludedInteractions_.removePair(b, c);
559 >      }
560 >
561 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
562 >        oneThreeInteractions_.removePair(a, c);      
563 >      } else {
564 >        excludedInteractions_.removePair(a, c);
565 >      }
566      }
567  
568 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
568 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
569 >         torsion = mol->nextTorsion(torsionIter)) {
570 >
571        a = torsion->getAtomA()->getGlobalIndex();
572        b = torsion->getAtomB()->getGlobalIndex();        
573        c = torsion->getAtomC()->getGlobalIndex();        
574 <      d = torsion->getAtomD()->getGlobalIndex();        
574 >      d = torsion->getAtomD()->getGlobalIndex();      
575 >  
576 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
577 >        oneTwoInteractions_.removePair(a, b);      
578 >        oneTwoInteractions_.removePair(b, c);
579 >        oneTwoInteractions_.removePair(c, d);
580 >      } else {
581 >        excludedInteractions_.removePair(a, b);
582 >        excludedInteractions_.removePair(b, c);
583 >        excludedInteractions_.removePair(c, d);
584 >      }
585  
586 <      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
587 <      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
588 <      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
589 <      std::set<int> rigidSetD = getRigidSet(d, atomGroups);
586 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
587 >        oneThreeInteractions_.removePair(a, c);      
588 >        oneThreeInteractions_.removePair(b, d);      
589 >      } else {
590 >        excludedInteractions_.removePair(a, c);
591 >        excludedInteractions_.removePair(b, d);
592 >      }
593  
594 <      exclude_.removePairs(rigidSetA, rigidSetB);
595 <      exclude_.removePairs(rigidSetA, rigidSetC);
596 <      exclude_.removePairs(rigidSetA, rigidSetD);
597 <      exclude_.removePairs(rigidSetB, rigidSetC);
598 <      exclude_.removePairs(rigidSetB, rigidSetD);
599 <      exclude_.removePairs(rigidSetC, rigidSetD);
594 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
595 >        oneFourInteractions_.removePair(a, d);      
596 >      } else {
597 >        excludedInteractions_.removePair(a, d);
598 >      }
599 >    }
600  
601 <      /*
602 <      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
552 <      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
553 <      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
554 <      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
555 <      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
556 <      exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
601 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
602 >         inversion = mol->nextInversion(inversionIter)) {
603  
604 <      
605 <      exclude_.removePair(a, b);
606 <      exclude_.removePair(a, c);
607 <      exclude_.removePair(a, d);
608 <      exclude_.removePair(b, c);
609 <      exclude_.removePair(b, d);
610 <      exclude_.removePair(c, d);        
611 <      */
604 >      a = inversion->getAtomA()->getGlobalIndex();
605 >      b = inversion->getAtomB()->getGlobalIndex();        
606 >      c = inversion->getAtomC()->getGlobalIndex();        
607 >      d = inversion->getAtomD()->getGlobalIndex();        
608 >
609 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
610 >        oneTwoInteractions_.removePair(a, b);      
611 >        oneTwoInteractions_.removePair(a, c);
612 >        oneTwoInteractions_.removePair(a, d);
613 >      } else {
614 >        excludedInteractions_.removePair(a, b);
615 >        excludedInteractions_.removePair(a, c);
616 >        excludedInteractions_.removePair(a, d);
617 >      }
618 >
619 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
620 >        oneThreeInteractions_.removePair(b, c);    
621 >        oneThreeInteractions_.removePair(b, d);    
622 >        oneThreeInteractions_.removePair(c, d);      
623 >      } else {
624 >        excludedInteractions_.removePair(b, c);
625 >        excludedInteractions_.removePair(b, d);
626 >        excludedInteractions_.removePair(c, d);
627 >      }
628      }
629  
630 <    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
631 <      std::vector<Atom*> atoms = rb->getAtoms();
632 <      for (int i = 0; i < atoms.size() -1 ; ++i) {
633 <        for (int j = i + 1; j < atoms.size(); ++j) {
630 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
631 >         rb = mol->nextRigidBody(rbIter)) {
632 >      vector<Atom*> atoms = rb->getAtoms();
633 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
634 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
635            a = atoms[i]->getGlobalIndex();
636            b = atoms[j]->getGlobalIndex();
637 <          exclude_.removePair(a, b);
637 >          excludedInteractions_.removePair(a, b);
638          }
639        }
640      }        
641 <
641 >    
642    }
643 <
644 <
643 >  
644 >  
645    void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
646      int curStampId;
647 <
647 >    
648      //index from 0
649      curStampId = moleculeStamps_.size();
650  
# Line 589 | Line 652 | namespace oopse {
652      molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
653    }
654  
655 +
656 +  /**
657 +   * update
658 +   *
659 +   *  Performs the global checks and variable settings after the objects have been
660 +   *  created.
661 +   *
662 +   */
663    void SimInfo::update() {
664 +    
665 +    setupSimVariables();
666 +    setupCutoffs();
667 +    setupSwitching();
668 +    setupElectrostatics();
669 +    setupNeighborlists();
670  
594    setupSimType();
595
671   #ifdef IS_MPI
672      setupFortranParallel();
673   #endif
599
674      setupFortranSim();
675 +    fortranInitialized_ = true;
676  
602    //setup fortran force field
603    /** @deprecate */    
604    int isError = 0;
605    
606    setupCutoff();
607    
608    setupElectrostaticSummationMethod( isError );
609    setupSwitchingFunction();
610    setupAccumulateBoxDipole();
611
612    if(isError){
613      sprintf( painCave.errMsg,
614               "ForceField error: There was an error initializing the forceField in fortran.\n" );
615      painCave.isFatal = 1;
616      simError();
617    }
618
677      calcNdf();
678      calcNdfRaw();
679      calcNdfTrans();
622
623    fortranInitialized_ = true;
680    }
681 <
682 <  std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
681 >  
682 >  set<AtomType*> SimInfo::getSimulatedAtomTypes() {
683      SimInfo::MoleculeIterator mi;
684      Molecule* mol;
685      Molecule::AtomIterator ai;
686      Atom* atom;
687 <    std::set<AtomType*> atomTypes;
688 <
689 <    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
634 <
687 >    set<AtomType*> atomTypes;
688 >    
689 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {      
690        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
691          atomTypes.insert(atom->getAtomType());
692 <      }
693 <        
692 >      }      
693 >    }    
694 >    return atomTypes;        
695 >  }
696 >
697 >  /**
698 >   * setupCutoffs
699 >   *
700 >   * Sets the values of cutoffRadius and cutoffMethod
701 >   *
702 >   * cutoffRadius : realType
703 >   *  If the cutoffRadius was explicitly set, use that value.
704 >   *  If the cutoffRadius was not explicitly set:
705 >   *      Are there electrostatic atoms?  Use 12.0 Angstroms.
706 >   *      No electrostatic atoms?  Poll the atom types present in the
707 >   *      simulation for suggested cutoff values (e.g. 2.5 * sigma).
708 >   *      Use the maximum suggested value that was found.
709 >   *
710 >   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL)
711 >   *      If cutoffMethod was explicitly set, use that choice.
712 >   *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
713 >   */
714 >  void SimInfo::setupCutoffs() {
715 >    
716 >    if (simParams_->haveCutoffRadius()) {
717 >      cutoffRadius_ = simParams_->getCutoffRadius();
718 >    } else {      
719 >      if (usesElectrostaticAtoms_) {
720 >        sprintf(painCave.errMsg,
721 >                "SimInfo: No value was set for the cutoffRadius.\n"
722 >                "\tOpenMD will use a default value of 12.0 angstroms"
723 >                "\tfor the cutoffRadius.\n");
724 >        painCave.isFatal = 0;
725 >        painCave.severity = OPENMD_INFO;
726 >        simError();
727 >        cutoffRadius_ = 12.0;
728 >      } else {
729 >        RealType thisCut;
730 >        set<AtomType*>::iterator i;
731 >        set<AtomType*> atomTypes;
732 >        atomTypes = getSimulatedAtomTypes();        
733 >        for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
734 >          thisCut = InteractionManager::Instance()->getSuggestedCutoffRadius((*i));
735 >          cutoffRadius_ = max(thisCut, cutoffRadius_);
736 >        }
737 >        sprintf(painCave.errMsg,
738 >                "SimInfo: No value was set for the cutoffRadius.\n"
739 >                "\tOpenMD will use %lf angstroms.\n",
740 >                cutoffRadius_);
741 >        painCave.isFatal = 0;
742 >        painCave.severity = OPENMD_INFO;
743 >        simError();
744 >      }            
745      }
746  
747 <    return atomTypes;        
747 >    map<string, CutoffMethod> stringToCutoffMethod;
748 >    stringToCutoffMethod["HARD"] = HARD;
749 >    stringToCutoffMethod["SWITCHING_FUNCTION"] = SWITCHING_FUNCTION;
750 >    stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;    
751 >    stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
752 >  
753 >    if (simParams_->haveCutoffMethod()) {
754 >      string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
755 >      map<string, CutoffMethod>::iterator i;
756 >      i = stringToCutoffMethod.find(cutMeth);
757 >      if (i == stringToCutoffMethod.end()) {
758 >        sprintf(painCave.errMsg,
759 >                "SimInfo: Could not find chosen cutoffMethod %s\n"
760 >                "\tShould be one of: "
761 >                "HARD, SWITCHING_FUNCTION, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n",
762 >                cutMeth.c_str());
763 >        painCave.isFatal = 1;
764 >        painCave.severity = OPENMD_ERROR;
765 >        simError();
766 >      } else {
767 >        cutoffMethod_ = i->second;
768 >      }
769 >    } else {
770 >      sprintf(painCave.errMsg,
771 >              "SimInfo: No value was set for the cutoffMethod.\n"
772 >              "\tOpenMD will use SHIFTED_FORCE.\n");
773 >        painCave.isFatal = 0;
774 >        painCave.severity = OPENMD_INFO;
775 >        simError();
776 >        cutoffMethod_ = SHIFTED_FORCE;        
777 >    }
778    }
779 <
780 <  void SimInfo::setupSimType() {
781 <    std::set<AtomType*>::iterator i;
782 <    std::set<AtomType*> atomTypes;
783 <    atomTypes = getUniqueAtomTypes();
779 >  
780 >  /**
781 >   * setupSwitching
782 >   *
783 >   * Sets the values of switchingRadius and
784 >   *  If the switchingRadius was explicitly set, use that value (but check it)
785 >   *  If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
786 >   */
787 >  void SimInfo::setupSwitching() {
788      
789 <    int useLennardJones = 0;
790 <    int useElectrostatic = 0;
791 <    int useEAM = 0;
792 <    int useSC = 0;
793 <    int useCharge = 0;
794 <    int useDirectional = 0;
795 <    int useDipole = 0;
796 <    int useGayBerne = 0;
797 <    int useSticky = 0;
798 <    int useStickyPower = 0;
799 <    int useShape = 0;
800 <    int useFLARB = 0; //it is not in AtomType yet
801 <    int useDirectionalAtom = 0;    
802 <    int useElectrostatics = 0;
803 <    //usePBC and useRF are from simParams
804 <    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
805 <    int useRF;
806 <    int useSF;
807 <    int useSP;
808 <    int useBoxDipole;
809 <    std::string myMethod;
789 >    if (simParams_->haveSwitchingRadius()) {
790 >      switchingRadius_ = simParams_->getSwitchingRadius();
791 >      if (switchingRadius_ > cutoffRadius_) {        
792 >        sprintf(painCave.errMsg,
793 >                "SimInfo: switchingRadius (%f) is larger than cutoffRadius(%f)\n",
794 >                switchingRadius_, cutoffRadius_);
795 >        painCave.isFatal = 1;
796 >        painCave.severity = OPENMD_ERROR;
797 >        simError();
798 >      }
799 >    } else {      
800 >      switchingRadius_ = 0.85 * cutoffRadius_;
801 >      sprintf(painCave.errMsg,
802 >              "SimInfo: No value was set for the switchingRadius.\n"
803 >              "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
804 >              "\tswitchingRadius = %f. for this simulation\n", switchingRadius_);
805 >      painCave.isFatal = 0;
806 >      painCave.severity = OPENMD_WARNING;
807 >      simError();
808 >    }          
809 >    
810 >    if (simParams_->haveSwitchingFunctionType()) {
811 >      string funcType = simParams_->getSwitchingFunctionType();
812 >      toUpper(funcType);
813 >      if (funcType == "CUBIC") {
814 >        sft_ = cubic;
815 >      } else {
816 >        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
817 >          sft_ = fifth_order_poly;
818 >        } else {
819 >          // throw error        
820 >          sprintf( painCave.errMsg,
821 >                   "SimInfo : Unknown switchingFunctionType. (Input file specified %s .)\n"
822 >                   "\tswitchingFunctionType must be one of: "
823 >                   "\"cubic\" or \"fifth_order_polynomial\".",
824 >                   funcType.c_str() );
825 >          painCave.isFatal = 1;
826 >          painCave.severity = OPENMD_ERROR;
827 >          simError();
828 >        }          
829 >      }
830 >    }
831 >  }
832  
833 <    // set the useRF logical
834 <    useRF = 0;
835 <    useSF = 0;
836 <    useSP = 0;
833 >  /**
834 >   * setupNeighborlists
835 >   *
836 >   *  If the skinThickness was explicitly set, use that value (but check it)
837 >   *  If the skinThickness was not explicitly set: use 1.0 angstroms
838 >   */
839 >  void SimInfo::setupNeighborlists() {    
840 >    if (simParams_->haveSkinThickness()) {
841 >      skinThickness_ = simParams_->getSkinThickness();
842 >    } else {      
843 >      skinThickness_ = 1.0;
844 >      sprintf(painCave.errMsg,
845 >              "SimInfo: No value was set for the skinThickness.\n"
846 >              "\tOpenMD will use a default value of %f Angstroms\n"
847 >              "\tfor this simulation\n", skinThickness_);
848 >      painCave.severity = OPENMD_INFO;
849 >      painCave.isFatal = 0;
850 >      simError();
851 >    }            
852 >  }
853  
854 <
855 <    if (simParams_->haveElectrostaticSummationMethod()) {
856 <      std::string myMethod = simParams_->getElectrostaticSummationMethod();
857 <      toUpper(myMethod);
858 <      if (myMethod == "REACTION_FIELD"){
859 <        useRF = 1;
860 <      } else if (myMethod == "SHIFTED_FORCE"){
683 <        useSF = 1;
684 <      } else if (myMethod == "SHIFTED_POTENTIAL"){
685 <        useSP = 1;
854 >  void SimInfo::setupSimVariables() {
855 >    useAtomicVirial_ = simParams_->getUseAtomicVirial();
856 >    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
857 >    calcBoxDipole_ = false;
858 >    if ( simParams_->haveAccumulateBoxDipole() )
859 >      if ( simParams_->getAccumulateBoxDipole() ) {
860 >        calcBoxDipole_ = true;      
861        }
687    }
688    
689    if (simParams_->haveAccumulateBoxDipole())
690      if (simParams_->getAccumulateBoxDipole())
691        useBoxDipole = 1;
862  
863 +    set<AtomType*>::iterator i;
864 +    set<AtomType*> atomTypes;
865 +    atomTypes = getSimulatedAtomTypes();    
866 +    int usesElectrostatic = 0;
867 +    int usesMetallic = 0;
868 +    int usesDirectional = 0;
869      //loop over all of the atom types
870      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
871 <      useLennardJones |= (*i)->isLennardJones();
872 <      useElectrostatic |= (*i)->isElectrostatic();
873 <      useEAM |= (*i)->isEAM();
698 <      useSC |= (*i)->isSC();
699 <      useCharge |= (*i)->isCharge();
700 <      useDirectional |= (*i)->isDirectional();
701 <      useDipole |= (*i)->isDipole();
702 <      useGayBerne |= (*i)->isGayBerne();
703 <      useSticky |= (*i)->isSticky();
704 <      useStickyPower |= (*i)->isStickyPower();
705 <      useShape |= (*i)->isShape();
871 >      usesElectrostatic |= (*i)->isElectrostatic();
872 >      usesMetallic |= (*i)->isMetal();
873 >      usesDirectional |= (*i)->isDirectional();
874      }
875  
708    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
709      useDirectionalAtom = 1;
710    }
711
712    if (useCharge || useDipole) {
713      useElectrostatics = 1;
714    }
715
876   #ifdef IS_MPI    
877      int temp;
878 +    temp = usesDirectional;
879 +    MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
880  
881 <    temp = usePBC;
882 <    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
881 >    temp = usesMetallic;
882 >    MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
883  
884 <    temp = useDirectionalAtom;
885 <    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
724 <
725 <    temp = useLennardJones;
726 <    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
727 <
728 <    temp = useElectrostatics;
729 <    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
730 <
731 <    temp = useCharge;
732 <    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
733 <
734 <    temp = useDipole;
735 <    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
736 <
737 <    temp = useSticky;
738 <    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
739 <
740 <    temp = useStickyPower;
741 <    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
742 <    
743 <    temp = useGayBerne;
744 <    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
745 <
746 <    temp = useEAM;
747 <    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
748 <
749 <    temp = useSC;
750 <    MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
751 <    
752 <    temp = useShape;
753 <    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
754 <
755 <    temp = useFLARB;
756 <    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
757 <
758 <    temp = useRF;
759 <    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
760 <
761 <    temp = useSF;
762 <    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
763 <
764 <    temp = useSP;
765 <    MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
766 <
767 <    temp = useBoxDipole;
768 <    MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
769 <
884 >    temp = usesElectrostatic;
885 >    MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
886   #endif
887 <
888 <    fInfo_.SIM_uses_PBC = usePBC;    
889 <    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
890 <    fInfo_.SIM_uses_LennardJones = useLennardJones;
891 <    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
892 <    fInfo_.SIM_uses_Charges = useCharge;
777 <    fInfo_.SIM_uses_Dipoles = useDipole;
778 <    fInfo_.SIM_uses_Sticky = useSticky;
779 <    fInfo_.SIM_uses_StickyPower = useStickyPower;
780 <    fInfo_.SIM_uses_GayBerne = useGayBerne;
781 <    fInfo_.SIM_uses_EAM = useEAM;
782 <    fInfo_.SIM_uses_SC = useSC;
783 <    fInfo_.SIM_uses_Shapes = useShape;
784 <    fInfo_.SIM_uses_FLARB = useFLARB;
785 <    fInfo_.SIM_uses_RF = useRF;
786 <    fInfo_.SIM_uses_SF = useSF;
787 <    fInfo_.SIM_uses_SP = useSP;
788 <    fInfo_.SIM_uses_BoxDipole = useBoxDipole;
887 >    fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_;    
888 >    fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_;
889 >    fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_;
890 >    fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_;
891 >    fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_;
892 >    fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_;
893    }
894  
895    void SimInfo::setupFortranSim() {
896      int isError;
897 <    int nExclude;
898 <    std::vector<int> fortranGlobalGroupMembership;
897 >    int nExclude, nOneTwo, nOneThree, nOneFour;
898 >    vector<int> fortranGlobalGroupMembership;
899      
900 <    nExclude = exclude_.getSize();
900 >    notifyFortranSkinThickness(&skinThickness_);
901 >
902 >    int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0;
903 >    int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0;
904 >    notifyFortranCutoffs(&cutoffRadius_, &switchingRadius_, &ljsp, &ljsf);
905 >
906      isError = 0;
907  
908      //globalGroupMembership_ is filled by SimCreator    
# Line 802 | Line 911 | namespace oopse {
911      }
912  
913      //calculate mass ratio of cutoff group
914 <    std::vector<RealType> mfact;
914 >    vector<RealType> mfact;
915      SimInfo::MoleculeIterator mi;
916      Molecule* mol;
917      Molecule::CutoffGroupIterator ci;
# Line 825 | Line 934 | namespace oopse {
934            else
935              mfact.push_back( 1.0 );
936          }
828
937        }      
938      }
939  
940      //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
941 <    std::vector<int> identArray;
941 >    vector<int> identArray;
942  
943      //to avoid memory reallocation, reserve enough space identArray
944      identArray.reserve(getNAtoms());
# Line 843 | Line 951 | namespace oopse {
951  
952      //fill molMembershipArray
953      //molMembershipArray is filled by SimCreator    
954 <    std::vector<int> molMembershipArray(nGlobalAtoms_);
954 >    vector<int> molMembershipArray(nGlobalAtoms_);
955      for (int i = 0; i < nGlobalAtoms_; i++) {
956        molMembershipArray[i] = globalMolMembership_[i] + 1;
957      }
958      
959      //setup fortran simulation
852    int nGlobalExcludes = 0;
853    int* globalExcludes = NULL;
854    int* excludeList = exclude_.getExcludeList();
855    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
856                   &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
857                   &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
960  
961 <    if( isError ){
961 >    nExclude = excludedInteractions_.getSize();
962 >    nOneTwo = oneTwoInteractions_.getSize();
963 >    nOneThree = oneThreeInteractions_.getSize();
964 >    nOneFour = oneFourInteractions_.getSize();
965 >
966 >    int* excludeList = excludedInteractions_.getPairList();
967 >    int* oneTwoList = oneTwoInteractions_.getPairList();
968 >    int* oneThreeList = oneThreeInteractions_.getPairList();
969 >    int* oneFourList = oneFourInteractions_.getPairList();
970  
971 +    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0],
972 +                   &nExclude, excludeList,
973 +                   &nOneTwo, oneTwoList,
974 +                   &nOneThree, oneThreeList,
975 +                   &nOneFour, oneFourList,
976 +                   &molMembershipArray[0], &mfact[0], &nCutoffGroups_,
977 +                   &fortranGlobalGroupMembership[0], &isError);
978 +    
979 +    if( isError ){
980 +      
981        sprintf( painCave.errMsg,
982                 "There was an error setting the simulation information in fortran.\n" );
983        painCave.isFatal = 1;
984 <      painCave.severity = OOPSE_ERROR;
984 >      painCave.severity = OPENMD_ERROR;
985        simError();
986      }
987 <
988 < #ifdef IS_MPI
987 >    
988 >    
989      sprintf( checkPointMsg,
990               "succesfully sent the simulation information to fortran.\n");
991 <    MPIcheckPoint();
992 < #endif // is_mpi
993 <
991 >    
992 >    errorCheckPoint();
993 >    
994      // Setup number of neighbors in neighbor list if present
995      if (simParams_->haveNeighborListNeighbors()) {
996        int nlistNeighbors = simParams_->getNeighborListNeighbors();
# Line 881 | Line 1001 | namespace oopse {
1001    }
1002  
1003  
884 #ifdef IS_MPI
1004    void SimInfo::setupFortranParallel() {
1005 <    
1005 > #ifdef IS_MPI    
1006      //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
1007 <    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
1008 <    std::vector<int> localToGlobalCutoffGroupIndex;
1007 >    vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
1008 >    vector<int> localToGlobalCutoffGroupIndex;
1009      SimInfo::MoleculeIterator mi;
1010      Molecule::AtomIterator ai;
1011      Molecule::CutoffGroupIterator ci;
# Line 933 | Line 1052 | namespace oopse {
1052      }
1053  
1054      sprintf(checkPointMsg, " mpiRefresh successful.\n");
1055 <    MPIcheckPoint();
1055 >    errorCheckPoint();
1056  
938
939  }
940
1057   #endif
942
943  void SimInfo::setupCutoff() {          
944    
945    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
946
947    // Check the cutoff policy
948    int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default
949
950    std::string myPolicy;
951    if (forceFieldOptions_.haveCutoffPolicy()){
952      myPolicy = forceFieldOptions_.getCutoffPolicy();
953    }else if (simParams_->haveCutoffPolicy()) {
954      myPolicy = simParams_->getCutoffPolicy();
955    }
956
957    if (!myPolicy.empty()){
958      toUpper(myPolicy);
959      if (myPolicy == "MIX") {
960        cp = MIX_CUTOFF_POLICY;
961      } else {
962        if (myPolicy == "MAX") {
963          cp = MAX_CUTOFF_POLICY;
964        } else {
965          if (myPolicy == "TRADITIONAL") {            
966            cp = TRADITIONAL_CUTOFF_POLICY;
967          } else {
968            // throw error        
969            sprintf( painCave.errMsg,
970                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
971            painCave.isFatal = 1;
972            simError();
973          }    
974        }          
975      }
976    }          
977    notifyFortranCutoffPolicy(&cp);
978
979    // Check the Skin Thickness for neighborlists
980    RealType skin;
981    if (simParams_->haveSkinThickness()) {
982      skin = simParams_->getSkinThickness();
983      notifyFortranSkinThickness(&skin);
984    }            
985        
986    // Check if the cutoff was set explicitly:
987    if (simParams_->haveCutoffRadius()) {
988      rcut_ = simParams_->getCutoffRadius();
989      if (simParams_->haveSwitchingRadius()) {
990        rsw_  = simParams_->getSwitchingRadius();
991      } else {
992        if (fInfo_.SIM_uses_Charges |
993            fInfo_.SIM_uses_Dipoles |
994            fInfo_.SIM_uses_RF) {
995          
996          rsw_ = 0.85 * rcut_;
997          sprintf(painCave.errMsg,
998                  "SimCreator Warning: No value was set for the switchingRadius.\n"
999                  "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n"
1000                  "\tswitchingRadius = %f. for this simulation\n", rsw_);
1001        painCave.isFatal = 0;
1002        simError();
1003        } else {
1004          rsw_ = rcut_;
1005          sprintf(painCave.errMsg,
1006                  "SimCreator Warning: No value was set for the switchingRadius.\n"
1007                  "\tOOPSE will use the same value as the cutoffRadius.\n"
1008                  "\tswitchingRadius = %f. for this simulation\n", rsw_);
1009          painCave.isFatal = 0;
1010          simError();
1011        }
1012      }
1013      
1014      notifyFortranCutoffs(&rcut_, &rsw_);
1015      
1016    } else {
1017      
1018      // For electrostatic atoms, we'll assume a large safe value:
1019      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
1020        sprintf(painCave.errMsg,
1021                "SimCreator Warning: No value was set for the cutoffRadius.\n"
1022                "\tOOPSE will use a default value of 15.0 angstroms"
1023                "\tfor the cutoffRadius.\n");
1024        painCave.isFatal = 0;
1025        simError();
1026        rcut_ = 15.0;
1027      
1028        if (simParams_->haveElectrostaticSummationMethod()) {
1029          std::string myMethod = simParams_->getElectrostaticSummationMethod();
1030          toUpper(myMethod);
1031          if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") {
1032            if (simParams_->haveSwitchingRadius()){
1033              sprintf(painCave.errMsg,
1034                      "SimInfo Warning: A value was set for the switchingRadius\n"
1035                      "\teven though the electrostaticSummationMethod was\n"
1036                      "\tset to %s\n", myMethod.c_str());
1037              painCave.isFatal = 1;
1038              simError();            
1039            }
1040          }
1041        }
1042      
1043        if (simParams_->haveSwitchingRadius()){
1044          rsw_ = simParams_->getSwitchingRadius();
1045        } else {        
1046          sprintf(painCave.errMsg,
1047                  "SimCreator Warning: No value was set for switchingRadius.\n"
1048                  "\tOOPSE will use a default value of\n"
1049                  "\t0.85 * cutoffRadius for the switchingRadius\n");
1050          painCave.isFatal = 0;
1051          simError();
1052          rsw_ = 0.85 * rcut_;
1053        }
1054        notifyFortranCutoffs(&rcut_, &rsw_);
1055      } else {
1056        // We didn't set rcut explicitly, and we don't have electrostatic atoms, so
1057        // We'll punt and let fortran figure out the cutoffs later.
1058        
1059        notifyFortranYouAreOnYourOwn();
1060
1061      }
1062    }
1058    }
1059  
1065  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
1066    
1067    int errorOut;
1068    int esm =  NONE;
1069    int sm = UNDAMPED;
1070    RealType alphaVal;
1071    RealType dielectric;
1072    
1073    errorOut = isError;
1060  
1075    if (simParams_->haveElectrostaticSummationMethod()) {
1076      std::string myMethod = simParams_->getElectrostaticSummationMethod();
1077      toUpper(myMethod);
1078      if (myMethod == "NONE") {
1079        esm = NONE;
1080      } else {
1081        if (myMethod == "SWITCHING_FUNCTION") {
1082          esm = SWITCHING_FUNCTION;
1083        } else {
1084          if (myMethod == "SHIFTED_POTENTIAL") {
1085            esm = SHIFTED_POTENTIAL;
1086          } else {
1087            if (myMethod == "SHIFTED_FORCE") {            
1088              esm = SHIFTED_FORCE;
1089            } else {
1090              if (myMethod == "REACTION_FIELD") {
1091                esm = REACTION_FIELD;
1092                dielectric = simParams_->getDielectric();
1093                if (!simParams_->haveDielectric()) {
1094                  // throw warning
1095                  sprintf( painCave.errMsg,
1096                           "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
1097                           "\tA default value of %f will be used for the dielectric.\n", dielectric);
1098                  painCave.isFatal = 0;
1099                  simError();
1100                }
1101              } else {
1102                // throw error        
1103                sprintf( painCave.errMsg,
1104                         "SimInfo error: Unknown electrostaticSummationMethod.\n"
1105                         "\t(Input file specified %s .)\n"
1106                         "\telectrostaticSummationMethod must be one of: \"none\",\n"
1107                         "\t\"shifted_potential\", \"shifted_force\", or \n"
1108                         "\t\"reaction_field\".\n", myMethod.c_str() );
1109                painCave.isFatal = 1;
1110                simError();
1111              }    
1112            }          
1113          }
1114        }
1115      }
1116    }
1117    
1118    if (simParams_->haveElectrostaticScreeningMethod()) {
1119      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
1120      toUpper(myScreen);
1121      if (myScreen == "UNDAMPED") {
1122        sm = UNDAMPED;
1123      } else {
1124        if (myScreen == "DAMPED") {
1125          sm = DAMPED;
1126          if (!simParams_->haveDampingAlpha()) {
1127            // first set a cutoff dependent alpha value
1128            // we assume alpha depends linearly with rcut from 0 to 20.5 ang
1129            alphaVal = 0.5125 - rcut_* 0.025;
1130            // for values rcut > 20.5, alpha is zero
1131            if (alphaVal < 0) alphaVal = 0;
1132
1133            // throw warning
1134            sprintf( painCave.errMsg,
1135                     "SimInfo warning: dampingAlpha was not specified in the input file.\n"
1136                     "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_);
1137            painCave.isFatal = 0;
1138            simError();
1139          } else {
1140            alphaVal = simParams_->getDampingAlpha();
1141          }
1142          
1143        } else {
1144          // throw error        
1145          sprintf( painCave.errMsg,
1146                   "SimInfo error: Unknown electrostaticScreeningMethod.\n"
1147                   "\t(Input file specified %s .)\n"
1148                   "\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
1149                   "or \"damped\".\n", myScreen.c_str() );
1150          painCave.isFatal = 1;
1151          simError();
1152        }
1153      }
1154    }
1155    
1156    // let's pass some summation method variables to fortran
1157    setElectrostaticSummationMethod( &esm );
1158    setFortranElectrostaticMethod( &esm );
1159    setScreeningMethod( &sm );
1160    setDampingAlpha( &alphaVal );
1161    setReactionFieldDielectric( &dielectric );
1162    initFortranFF( &errorOut );
1163  }
1164
1165  void SimInfo::setupSwitchingFunction() {    
1166    int ft = CUBIC;
1167
1168    if (simParams_->haveSwitchingFunctionType()) {
1169      std::string funcType = simParams_->getSwitchingFunctionType();
1170      toUpper(funcType);
1171      if (funcType == "CUBIC") {
1172        ft = CUBIC;
1173      } else {
1174        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
1175          ft = FIFTH_ORDER_POLY;
1176        } else {
1177          // throw error        
1178          sprintf( painCave.errMsg,
1179                   "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
1180          painCave.isFatal = 1;
1181          simError();
1182        }          
1183      }
1184    }
1185
1186    // send switching function notification to switcheroo
1187    setFunctionType(&ft);
1188
1189  }
1190
1061    void SimInfo::setupAccumulateBoxDipole() {    
1062  
1193    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
1194    if ( simParams_->haveAccumulateBoxDipole() )
1195      if ( simParams_->getAccumulateBoxDipole() ) {
1196        setAccumulateBoxDipole();
1197        calcBoxDipole_ = true;
1198      }
1063  
1064    }
1065  
# Line 1203 | Line 1067 | namespace oopse {
1067      properties_.addProperty(genData);  
1068    }
1069  
1070 <  void SimInfo::removeProperty(const std::string& propName) {
1070 >  void SimInfo::removeProperty(const string& propName) {
1071      properties_.removeProperty(propName);  
1072    }
1073  
# Line 1211 | Line 1075 | namespace oopse {
1075      properties_.clearProperties();
1076    }
1077  
1078 <  std::vector<std::string> SimInfo::getPropertyNames() {
1078 >  vector<string> SimInfo::getPropertyNames() {
1079      return properties_.getPropertyNames();  
1080    }
1081        
1082 <  std::vector<GenericData*> SimInfo::getProperties() {
1082 >  vector<GenericData*> SimInfo::getProperties() {
1083      return properties_.getProperties();
1084    }
1085  
1086 <  GenericData* SimInfo::getPropertyByName(const std::string& propName) {
1086 >  GenericData* SimInfo::getPropertyByName(const string& propName) {
1087      return properties_.getPropertyByName(propName);
1088    }
1089  
# Line 1302 | Line 1166 | namespace oopse {
1166  
1167    }        
1168  
1169 <  std::ostream& operator <<(std::ostream& o, SimInfo& info) {
1169 >  ostream& operator <<(ostream& o, SimInfo& info) {
1170  
1171      return o;
1172    }
# Line 1345 | Line 1209 | namespace oopse {
1209  
1210  
1211         [  Ixx -Ixy  -Ixz ]
1212 <  J =| -Iyx  Iyy  -Iyz |
1212 >    J =| -Iyx  Iyy  -Iyz |
1213         [ -Izx -Iyz   Izz ]
1214      */
1215  
# Line 1452 | Line 1316 | namespace oopse {
1316      return IOIndexToIntegrableObject.at(index);
1317    }
1318    
1319 <  void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) {
1319 >  void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) {
1320      IOIndexToIntegrableObject= v;
1321    }
1322  
# Line 1494 | Line 1358 | namespace oopse {
1358      return;
1359    }
1360   /*
1361 <   void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) {
1361 >   void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) {
1362        assert( v.size() == nAtoms_ + nRigidBodies_);
1363        sdByGlobalIndex_ = v;
1364      }
# Line 1504 | Line 1368 | namespace oopse {
1368        return sdByGlobalIndex_.at(index);
1369      }  
1370   */  
1371 < }//end namespace oopse
1371 >  int SimInfo::getNGlobalConstraints() {
1372 >    int nGlobalConstraints;
1373 > #ifdef IS_MPI
1374 >    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
1375 >                  MPI_COMM_WORLD);    
1376 > #else
1377 >    nGlobalConstraints =  nConstraints_;
1378 > #endif
1379 >    return nGlobalConstraints;
1380 >  }
1381  
1382 + }//end namespace OpenMD
1383 +

Comparing:
trunk/src/brains/SimInfo.cpp (property svn:keywords), Revision 1121 by chuckv, Mon Feb 26 04:45:42 2007 UTC vs.
branches/development/src/brains/SimInfo.cpp (property svn:keywords), Revision 1534 by gezelter, Wed Dec 29 21:53:28 2010 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines