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

Comparing:
trunk/src/brains/SimInfo.cpp (property svn:keywords), Revision 1129 by chrisfen, Fri Apr 20 18:15:48 2007 UTC vs.
branches/development/src/brains/SimInfo.cpp (property svn:keywords), Revision 1528 by gezelter, Fri Dec 17 20:11:05 2010 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines