ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/brains/SimInfo.cpp
(Generate patch)

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 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, 234107 (2008).          
39 + * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 + * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   /**
# Line 46 | Line 47
47   * @version 1.0
48   */
49  
50 + #ifdef IS_MPI
51 + #include <mpi.h>
52 + #endif
53   #include <algorithm>
54   #include <set>
55 + #include <map>
56  
57   #include "brains/SimInfo.hpp"
58   #include "math/Vector3.hpp"
59   #include "primitives/Molecule.hpp"
60 < #include "UseTheForce/doForces_interface.h"
56 < #include "UseTheForce/notifyCutoffs_interface.h"
60 > #include "primitives/StuntDouble.hpp"
61   #include "utils/MemoryUtils.hpp"
62   #include "utils/simError.h"
63   #include "selection/SelectionManager.hpp"
64 + #include "io/ForceFieldOptions.hpp"
65 + #include "brains/ForceField.hpp"
66 + #include "nonbonded/SwitchingFunction.hpp"
67  
68 < #ifdef IS_MPI
69 < #include "UseTheForce/mpiComponentPlan.h"
70 < #include "UseTheForce/DarkSide/simParallel_interface.h"
71 < #endif
72 <
73 < namespace oopse {
67 <
68 <  SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
69 <                   ForceField* ff, Globals* simParams) :
70 <    stamps_(stamps), forceField_(ff), simParams_(simParams),
71 <    ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
68 > using namespace std;
69 > namespace OpenMD {
70 >  
71 >  SimInfo::SimInfo(ForceField* ff, Globals* simParams) :
72 >    forceField_(ff), simParams_(simParams),
73 >    ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
74      nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
75 <    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
76 <    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
77 <    nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
78 <    sman_(NULL), fortranInitialized_(false) {
79 <
80 <            
81 <      std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
80 <      MoleculeStamp* molStamp;
81 <      int nMolWithSameStamp;
82 <      int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
83 <      int nGroups = 0;          //total cutoff groups defined in meta-data file
84 <      CutoffGroupStamp* cgStamp;    
85 <      RigidBodyStamp* rbStamp;
86 <      int nRigidAtoms = 0;
75 >    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
76 >    nGlobalFluctuatingCharges_(0), nGlobalBonds_(0), nGlobalBends_(0),
77 >    nGlobalTorsions_(0), nGlobalInversions_(0), nAtoms_(0), nBonds_(0),
78 >    nBends_(0), nTorsions_(0), nInversions_(0), nRigidBodies_(0),
79 >    nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0),
80 >    nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false),
81 >    calcBoxDipole_(false), useAtomicVirial_(true) {    
82      
83 <      for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
84 <        molStamp = i->first;
85 <        nMolWithSameStamp = i->second;
86 <        
87 <        addMoleculeStamp(molStamp, nMolWithSameStamp);
88 <
89 <        //calculate atoms in molecules
90 <        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
91 <
92 <
93 <        //calculate atoms in cutoff groups
94 <        int nAtomsInGroups = 0;
95 <        int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
96 <        
97 <        for (int j=0; j < nCutoffGroupsInStamp; j++) {
98 <          cgStamp = molStamp->getCutoffGroup(j);
99 <          nAtomsInGroups += cgStamp->getNMembers();
100 <        }
106 <
107 <        nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
108 <        nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
109 <
110 <        //calculate atoms in rigid bodies
111 <        int nAtomsInRigidBodies = 0;
112 <        int nRigidBodiesInStamp = molStamp->getNRigidBodies();
113 <        
114 <        for (int j=0; j < nRigidBodiesInStamp; j++) {
115 <          rbStamp = molStamp->getRigidBody(j);
116 <          nAtomsInRigidBodies += rbStamp->getNMembers();
117 <        }
118 <
119 <        nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
120 <        nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
121 <        
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 >    
91 >    vector<Component*> components = simParams->getComponents();
92 >    
93 >    for (vector<Component*>::iterator i = components.begin();
94 >         i !=components.end(); ++i) {
95 >      molStamp = (*i)->getMoleculeStamp();
96 >      if ( (*i)->haveRegion() ) {        
97 >        molStamp->setRegion( (*i)->getRegion() );
98 >      } else {
99 >        // set the region to a disallowed value:
100 >        molStamp->setRegion( -1 );
101        }
102  
103 <      //every free atom (atom does not belong to cutoff groups) is a cutoff group
104 <      //therefore the total number of cutoff groups in the system is equal to
105 <      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
106 <      //file plus the number of cutoff groups defined in meta-data file
107 <      nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
108 <
109 <      //every free atom (atom does not belong to rigid bodies) is an integrable object
110 <      //therefore the total number of  integrable objects in the system is equal to
111 <      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
112 <      //file plus the number of  rigid bodies defined in meta-data file
113 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
114 <
115 <      nGlobalMols_ = molStampIds_.size();
116 <
117 < #ifdef IS_MPI    
118 <      molToProcMap_.resize(nGlobalMols_);
119 < #endif
120 <
103 >      nMolWithSameStamp = (*i)->getNMol();
104 >      
105 >      addMoleculeStamp(molStamp, nMolWithSameStamp);
106 >      
107 >      //calculate atoms in molecules
108 >      nGlobalAtoms_ += molStamp->getNAtoms() * nMolWithSameStamp;
109 >      nGlobalBonds_ += molStamp->getNBonds() * nMolWithSameStamp;
110 >      nGlobalBends_ += molStamp->getNBends() * nMolWithSameStamp;
111 >      nGlobalTorsions_ += molStamp->getNTorsions() * nMolWithSameStamp;
112 >      nGlobalInversions_ += molStamp->getNInversions() * nMolWithSameStamp;
113 >      
114 >      //calculate atoms in cutoff groups
115 >      int nAtomsInGroups = 0;
116 >      int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
117 >      
118 >      for (int j=0; j < nCutoffGroupsInStamp; j++) {
119 >        cgStamp = molStamp->getCutoffGroupStamp(j);
120 >        nAtomsInGroups += cgStamp->getNMembers();
121 >      }
122 >      
123 >      nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
124 >      
125 >      nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
126 >      
127 >      //calculate atoms in rigid bodies
128 >      int nAtomsInRigidBodies = 0;
129 >      int nRigidBodiesInStamp = molStamp->getNRigidBodies();
130 >      
131 >      for (int j=0; j < nRigidBodiesInStamp; j++) {
132 >        rbStamp = molStamp->getRigidBodyStamp(j);
133 >        nAtomsInRigidBodies += rbStamp->getNMembers();
134 >      }
135 >      
136 >      nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
137 >      nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
138 >      
139      }
140 +    
141 +    //every free atom (atom does not belong to cutoff groups) is a cutoff
142 +    //group therefore the total number of cutoff groups in the system is
143 +    //equal to the total number of atoms minus number of atoms belong to
144 +    //cutoff group defined in meta-data file plus the number of cutoff
145 +    //groups defined in meta-data file
146  
147 +    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
148 +    
149 +    //every free atom (atom does not belong to rigid bodies) is an
150 +    //integrable object therefore the total number of integrable objects
151 +    //in the system is equal to the total number of atoms minus number of
152 +    //atoms belong to rigid body defined in meta-data file plus the number
153 +    //of rigid bodies defined in meta-data file
154 +    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
155 +      + nGlobalRigidBodies_;
156 +    
157 +    nGlobalMols_ = molStampIds_.size();
158 +    molToProcMap_.resize(nGlobalMols_);
159 +  }
160 +  
161    SimInfo::~SimInfo() {
162 <    std::map<int, Molecule*>::iterator i;
162 >    map<int, Molecule*>::iterator i;
163      for (i = molecules_.begin(); i != molecules_.end(); ++i) {
164        delete i->second;
165      }
166      molecules_.clear();
167        
151    delete stamps_;
168      delete sman_;
169      delete simParams_;
170      delete forceField_;
171    }
172  
157  int SimInfo::getNGlobalConstraints() {
158    int nGlobalConstraints;
159 #ifdef IS_MPI
160    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
161                  MPI_COMM_WORLD);    
162 #else
163    nGlobalConstraints =  nConstraints_;
164 #endif
165    return nGlobalConstraints;
166  }
173  
174    bool SimInfo::addMolecule(Molecule* mol) {
175      MoleculeIterator i;
176 <
176 >    
177      i = molecules_.find(mol->getGlobalIndex());
178      if (i == molecules_.end() ) {
179 <
180 <      molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
181 <        
179 >      
180 >      molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
181 >      
182        nAtoms_ += mol->getNAtoms();
183        nBonds_ += mol->getNBonds();
184        nBends_ += mol->getNBends();
185        nTorsions_ += mol->getNTorsions();
186 +      nInversions_ += mol->getNInversions();
187        nRigidBodies_ += mol->getNRigidBodies();
188        nIntegrableObjects_ += mol->getNIntegrableObjects();
189        nCutoffGroups_ += mol->getNCutoffGroups();
190        nConstraints_ += mol->getNConstraintPairs();
191 <
192 <      addExcludePairs(mol);
193 <        
191 >      
192 >      addInteractionPairs(mol);
193 >      
194        return true;
195      } else {
196        return false;
197      }
198    }
199 <
199 >  
200    bool SimInfo::removeMolecule(Molecule* mol) {
201      MoleculeIterator i;
202      i = molecules_.find(mol->getGlobalIndex());
# Line 202 | Line 209 | namespace oopse {
209        nBonds_ -= mol->getNBonds();
210        nBends_ -= mol->getNBends();
211        nTorsions_ -= mol->getNTorsions();
212 +      nInversions_ -= mol->getNInversions();
213        nRigidBodies_ -= mol->getNRigidBodies();
214        nIntegrableObjects_ -= mol->getNIntegrableObjects();
215        nCutoffGroups_ -= mol->getNCutoffGroups();
216        nConstraints_ -= mol->getNConstraintPairs();
217  
218 <      removeExcludePairs(mol);
218 >      removeInteractionPairs(mol);
219        molecules_.erase(mol->getGlobalIndex());
220  
221        delete mol;
# Line 216 | Line 224 | namespace oopse {
224      } else {
225        return false;
226      }
219
220
227    }    
228  
229          
# Line 233 | Line 239 | namespace oopse {
239  
240  
241    void SimInfo::calcNdf() {
242 <    int ndf_local;
242 >    int ndf_local, nfq_local;
243      MoleculeIterator i;
244 <    std::vector<StuntDouble*>::iterator j;
244 >    vector<StuntDouble*>::iterator j;
245 >    vector<Atom*>::iterator k;
246 >
247      Molecule* mol;
248 <    StuntDouble* integrableObject;
248 >    StuntDouble* sd;
249 >    Atom* atom;
250  
251      ndf_local = 0;
252 +    nfq_local = 0;
253      
254      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
245      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
246           integrableObject = mol->nextIntegrableObject(j)) {
255  
256 +      for (sd = mol->beginIntegrableObject(j); sd != NULL;
257 +           sd = mol->nextIntegrableObject(j)) {
258 +
259          ndf_local += 3;
260  
261 <        if (integrableObject->isDirectional()) {
262 <          if (integrableObject->isLinear()) {
261 >        if (sd->isDirectional()) {
262 >          if (sd->isLinear()) {
263              ndf_local += 2;
264            } else {
265              ndf_local += 3;
266            }
267          }
268 <            
269 <      }//end for (integrableObject)
270 <    }// end for (mol)
268 >      }
269 >
270 >      for (atom = mol->beginFluctuatingCharge(k); atom != NULL;
271 >           atom = mol->nextFluctuatingCharge(k)) {
272 >        if (atom->isFluctuatingCharge()) {
273 >          nfq_local++;
274 >        }
275 >      }
276 >    }
277      
278 +    ndfLocal_ = ndf_local;
279 +
280      // n_constraints is local, so subtract them on each processor
281      ndf_local -= nConstraints_;
282  
283   #ifdef IS_MPI
284 <    MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
284 >    MPI_Allreduce(&ndf_local, &ndf_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
285 >    MPI_Allreduce(&nfq_local, &nGlobalFluctuatingCharges_, 1,
286 >      MPI_INT, MPI_SUM, MPI_COMM_WORLD);
287 >    // MPI::COMM_WORLD.Allreduce(&ndf_local, &ndf_, 1, MPI::INT,MPI::SUM);
288 >    // MPI::COMM_WORLD.Allreduce(&nfq_local, &nGlobalFluctuatingCharges_, 1,
289 >    //                           MPI::INT, MPI::SUM);
290   #else
291      ndf_ = ndf_local;
292 +    nGlobalFluctuatingCharges_ = nfq_local;
293   #endif
294  
295      // nZconstraints_ is global, as are the 3 COM translations for the
# Line 273 | Line 298 | namespace oopse {
298  
299    }
300  
301 +  int SimInfo::getFdf() {
302 + #ifdef IS_MPI
303 +    MPI_Allreduce(&fdf_local, &fdf_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
304 +    // MPI::COMM_WORLD.Allreduce(&fdf_local, &fdf_, 1, MPI::INT, MPI::SUM);
305 + #else
306 +    fdf_ = fdf_local;
307 + #endif
308 +    return fdf_;
309 +  }
310 +  
311 +  unsigned int SimInfo::getNLocalCutoffGroups(){
312 +    int nLocalCutoffAtoms = 0;
313 +    Molecule* mol;
314 +    MoleculeIterator mi;
315 +    CutoffGroup* cg;
316 +    Molecule::CutoffGroupIterator ci;
317 +    
318 +    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
319 +      
320 +      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
321 +           cg = mol->nextCutoffGroup(ci)) {
322 +        nLocalCutoffAtoms += cg->getNumAtom();
323 +        
324 +      }        
325 +    }
326 +    
327 +    return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_;
328 +  }
329 +    
330    void SimInfo::calcNdfRaw() {
331      int ndfRaw_local;
332  
333      MoleculeIterator i;
334 <    std::vector<StuntDouble*>::iterator j;
334 >    vector<StuntDouble*>::iterator j;
335      Molecule* mol;
336 <    StuntDouble* integrableObject;
336 >    StuntDouble* sd;
337  
338      // Raw degrees of freedom that we have to set
339      ndfRaw_local = 0;
340      
341      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
288      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
289           integrableObject = mol->nextIntegrableObject(j)) {
342  
343 +      for (sd = mol->beginIntegrableObject(j); sd != NULL;
344 +           sd = mol->nextIntegrableObject(j)) {
345 +
346          ndfRaw_local += 3;
347  
348 <        if (integrableObject->isDirectional()) {
349 <          if (integrableObject->isLinear()) {
348 >        if (sd->isDirectional()) {
349 >          if (sd->isLinear()) {
350              ndfRaw_local += 2;
351            } else {
352              ndfRaw_local += 3;
# Line 302 | Line 357 | namespace oopse {
357      }
358      
359   #ifdef IS_MPI
360 <    MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
360 >    MPI_Allreduce(&ndfRaw_local, &ndfRaw_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
361 >    // MPI::COMM_WORLD.Allreduce(&ndfRaw_local, &ndfRaw_, 1, MPI::INT, MPI::SUM);
362   #else
363      ndfRaw_ = ndfRaw_local;
364   #endif
# Line 315 | Line 371 | namespace oopse {
371  
372  
373   #ifdef IS_MPI
374 <    MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
374 >    MPI_Allreduce(&ndfTrans_local, &ndfTrans_, 1,
375 >      MPI_INT, MPI_SUM, MPI_COMM_WORLD);
376 >    // MPI::COMM_WORLD.Allreduce(&ndfTrans_local, &ndfTrans_, 1,
377 >    //                           MPI::INT, MPI::SUM);
378   #else
379      ndfTrans_ = ndfTrans_local;
380   #endif
# Line 324 | Line 383 | namespace oopse {
383  
384    }
385  
386 <  void SimInfo::addExcludePairs(Molecule* mol) {
387 <    std::vector<Bond*>::iterator bondIter;
388 <    std::vector<Bend*>::iterator bendIter;
389 <    std::vector<Torsion*>::iterator torsionIter;
386 >  void SimInfo::addInteractionPairs(Molecule* mol) {
387 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
388 >    vector<Bond*>::iterator bondIter;
389 >    vector<Bend*>::iterator bendIter;
390 >    vector<Torsion*>::iterator torsionIter;
391 >    vector<Inversion*>::iterator inversionIter;
392      Bond* bond;
393      Bend* bend;
394      Torsion* torsion;
395 +    Inversion* inversion;
396      int a;
397      int b;
398      int c;
399      int d;
400 +
401 +    // atomGroups can be used to add special interaction maps between
402 +    // groups of atoms that are in two separate rigid bodies.
403 +    // However, most site-site interactions between two rigid bodies
404 +    // are probably not special, just the ones between the physically
405 +    // bonded atoms.  Interactions *within* a single rigid body should
406 +    // always be excluded.  These are done at the bottom of this
407 +    // function.
408 +
409 +    map<int, set<int> > atomGroups;
410 +    Molecule::RigidBodyIterator rbIter;
411 +    RigidBody* rb;
412 +    Molecule::IntegrableObjectIterator ii;
413 +    StuntDouble* sd;
414      
415 <    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
415 >    for (sd = mol->beginIntegrableObject(ii); sd != NULL;
416 >         sd = mol->nextIntegrableObject(ii)) {
417 >      
418 >      if (sd->isRigidBody()) {
419 >        rb = static_cast<RigidBody*>(sd);
420 >        vector<Atom*> atoms = rb->getAtoms();
421 >        set<int> rigidAtoms;
422 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
423 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
424 >        }
425 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
426 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
427 >        }      
428 >      } else {
429 >        set<int> oneAtomSet;
430 >        oneAtomSet.insert(sd->getGlobalIndex());
431 >        atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet));        
432 >      }
433 >    }  
434 >
435 >          
436 >    for (bond= mol->beginBond(bondIter); bond != NULL;
437 >         bond = mol->nextBond(bondIter)) {
438 >
439        a = bond->getAtomA()->getGlobalIndex();
440 <      b = bond->getAtomB()->getGlobalIndex();        
441 <      exclude_.addPair(a, b);
440 >      b = bond->getAtomB()->getGlobalIndex();  
441 >
442 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
443 >        oneTwoInteractions_.addPair(a, b);
444 >      } else {
445 >        excludedInteractions_.addPair(a, b);
446 >      }
447      }
448  
449 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
449 >    for (bend= mol->beginBend(bendIter); bend != NULL;
450 >         bend = mol->nextBend(bendIter)) {
451 >
452        a = bend->getAtomA()->getGlobalIndex();
453        b = bend->getAtomB()->getGlobalIndex();        
454        c = bend->getAtomC()->getGlobalIndex();
455 +      
456 +      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
457 +        oneTwoInteractions_.addPair(a, b);      
458 +        oneTwoInteractions_.addPair(b, c);
459 +      } else {
460 +        excludedInteractions_.addPair(a, b);
461 +        excludedInteractions_.addPair(b, c);
462 +      }
463  
464 <      exclude_.addPair(a, b);
465 <      exclude_.addPair(a, c);
466 <      exclude_.addPair(b, c);        
464 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
465 >        oneThreeInteractions_.addPair(a, c);      
466 >      } else {
467 >        excludedInteractions_.addPair(a, c);
468 >      }
469      }
470  
471 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
471 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
472 >         torsion = mol->nextTorsion(torsionIter)) {
473 >
474        a = torsion->getAtomA()->getGlobalIndex();
475        b = torsion->getAtomB()->getGlobalIndex();        
476        c = torsion->getAtomC()->getGlobalIndex();        
477 <      d = torsion->getAtomD()->getGlobalIndex();        
477 >      d = torsion->getAtomD()->getGlobalIndex();      
478  
479 <      exclude_.addPair(a, b);
480 <      exclude_.addPair(a, c);
481 <      exclude_.addPair(a, d);
482 <      exclude_.addPair(b, c);
483 <      exclude_.addPair(b, d);
484 <      exclude_.addPair(c, d);        
479 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
480 >        oneTwoInteractions_.addPair(a, b);      
481 >        oneTwoInteractions_.addPair(b, c);
482 >        oneTwoInteractions_.addPair(c, d);
483 >      } else {
484 >        excludedInteractions_.addPair(a, b);
485 >        excludedInteractions_.addPair(b, c);
486 >        excludedInteractions_.addPair(c, d);
487 >      }
488 >
489 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
490 >        oneThreeInteractions_.addPair(a, c);      
491 >        oneThreeInteractions_.addPair(b, d);      
492 >      } else {
493 >        excludedInteractions_.addPair(a, c);
494 >        excludedInteractions_.addPair(b, d);
495 >      }
496 >
497 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
498 >        oneFourInteractions_.addPair(a, d);      
499 >      } else {
500 >        excludedInteractions_.addPair(a, d);
501 >      }
502      }
503  
504 <    Molecule::RigidBodyIterator rbIter;
505 <    RigidBody* rb;
506 <    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
507 <      std::vector<Atom*> atoms = rb->getAtoms();
508 <      for (int i = 0; i < atoms.size() -1 ; ++i) {
509 <        for (int j = i + 1; j < atoms.size(); ++j) {
504 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
505 >         inversion = mol->nextInversion(inversionIter)) {
506 >
507 >      a = inversion->getAtomA()->getGlobalIndex();
508 >      b = inversion->getAtomB()->getGlobalIndex();        
509 >      c = inversion->getAtomC()->getGlobalIndex();        
510 >      d = inversion->getAtomD()->getGlobalIndex();        
511 >
512 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
513 >        oneTwoInteractions_.addPair(a, b);      
514 >        oneTwoInteractions_.addPair(a, c);
515 >        oneTwoInteractions_.addPair(a, d);
516 >      } else {
517 >        excludedInteractions_.addPair(a, b);
518 >        excludedInteractions_.addPair(a, c);
519 >        excludedInteractions_.addPair(a, d);
520 >      }
521 >
522 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
523 >        oneThreeInteractions_.addPair(b, c);    
524 >        oneThreeInteractions_.addPair(b, d);    
525 >        oneThreeInteractions_.addPair(c, d);      
526 >      } else {
527 >        excludedInteractions_.addPair(b, c);
528 >        excludedInteractions_.addPair(b, d);
529 >        excludedInteractions_.addPair(c, d);
530 >      }
531 >    }
532 >
533 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
534 >         rb = mol->nextRigidBody(rbIter)) {
535 >      vector<Atom*> atoms = rb->getAtoms();
536 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
537 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
538            a = atoms[i]->getGlobalIndex();
539            b = atoms[j]->getGlobalIndex();
540 <          exclude_.addPair(a, b);
540 >          excludedInteractions_.addPair(a, b);
541          }
542        }
543      }        
544  
545    }
546  
547 <  void SimInfo::removeExcludePairs(Molecule* mol) {
548 <    std::vector<Bond*>::iterator bondIter;
549 <    std::vector<Bend*>::iterator bendIter;
550 <    std::vector<Torsion*>::iterator torsionIter;
547 >  void SimInfo::removeInteractionPairs(Molecule* mol) {
548 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
549 >    vector<Bond*>::iterator bondIter;
550 >    vector<Bend*>::iterator bendIter;
551 >    vector<Torsion*>::iterator torsionIter;
552 >    vector<Inversion*>::iterator inversionIter;
553      Bond* bond;
554      Bend* bend;
555      Torsion* torsion;
556 +    Inversion* inversion;
557      int a;
558      int b;
559      int c;
560      int d;
561 +
562 +    map<int, set<int> > atomGroups;
563 +    Molecule::RigidBodyIterator rbIter;
564 +    RigidBody* rb;
565 +    Molecule::IntegrableObjectIterator ii;
566 +    StuntDouble* sd;
567      
568 <    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
568 >    for (sd = mol->beginIntegrableObject(ii); sd != NULL;
569 >         sd = mol->nextIntegrableObject(ii)) {
570 >      
571 >      if (sd->isRigidBody()) {
572 >        rb = static_cast<RigidBody*>(sd);
573 >        vector<Atom*> atoms = rb->getAtoms();
574 >        set<int> rigidAtoms;
575 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
576 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
577 >        }
578 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
579 >          atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
580 >        }      
581 >      } else {
582 >        set<int> oneAtomSet;
583 >        oneAtomSet.insert(sd->getGlobalIndex());
584 >        atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet));        
585 >      }
586 >    }  
587 >
588 >    for (bond= mol->beginBond(bondIter); bond != NULL;
589 >         bond = mol->nextBond(bondIter)) {
590 >      
591        a = bond->getAtomA()->getGlobalIndex();
592 <      b = bond->getAtomB()->getGlobalIndex();        
593 <      exclude_.removePair(a, b);
592 >      b = bond->getAtomB()->getGlobalIndex();  
593 >    
594 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
595 >        oneTwoInteractions_.removePair(a, b);
596 >      } else {
597 >        excludedInteractions_.removePair(a, b);
598 >      }
599      }
600  
601 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
601 >    for (bend= mol->beginBend(bendIter); bend != NULL;
602 >         bend = mol->nextBend(bendIter)) {
603 >
604        a = bend->getAtomA()->getGlobalIndex();
605        b = bend->getAtomB()->getGlobalIndex();        
606        c = bend->getAtomC()->getGlobalIndex();
607 +      
608 +      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
609 +        oneTwoInteractions_.removePair(a, b);      
610 +        oneTwoInteractions_.removePair(b, c);
611 +      } else {
612 +        excludedInteractions_.removePair(a, b);
613 +        excludedInteractions_.removePair(b, c);
614 +      }
615  
616 <      exclude_.removePair(a, b);
617 <      exclude_.removePair(a, c);
618 <      exclude_.removePair(b, c);        
616 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
617 >        oneThreeInteractions_.removePair(a, c);      
618 >      } else {
619 >        excludedInteractions_.removePair(a, c);
620 >      }
621      }
622  
623 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
623 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
624 >         torsion = mol->nextTorsion(torsionIter)) {
625 >
626        a = torsion->getAtomA()->getGlobalIndex();
627        b = torsion->getAtomB()->getGlobalIndex();        
628        c = torsion->getAtomC()->getGlobalIndex();        
629 <      d = torsion->getAtomD()->getGlobalIndex();        
629 >      d = torsion->getAtomD()->getGlobalIndex();      
630 >  
631 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
632 >        oneTwoInteractions_.removePair(a, b);      
633 >        oneTwoInteractions_.removePair(b, c);
634 >        oneTwoInteractions_.removePair(c, d);
635 >      } else {
636 >        excludedInteractions_.removePair(a, b);
637 >        excludedInteractions_.removePair(b, c);
638 >        excludedInteractions_.removePair(c, d);
639 >      }
640  
641 <      exclude_.removePair(a, b);
642 <      exclude_.removePair(a, c);
643 <      exclude_.removePair(a, d);
644 <      exclude_.removePair(b, c);
645 <      exclude_.removePair(b, d);
646 <      exclude_.removePair(c, d);        
641 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
642 >        oneThreeInteractions_.removePair(a, c);      
643 >        oneThreeInteractions_.removePair(b, d);      
644 >      } else {
645 >        excludedInteractions_.removePair(a, c);
646 >        excludedInteractions_.removePair(b, d);
647 >      }
648 >
649 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
650 >        oneFourInteractions_.removePair(a, d);      
651 >      } else {
652 >        excludedInteractions_.removePair(a, d);
653 >      }
654      }
655  
656 <    Molecule::RigidBodyIterator rbIter;
657 <    RigidBody* rb;
658 <    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
659 <      std::vector<Atom*> atoms = rb->getAtoms();
660 <      for (int i = 0; i < atoms.size() -1 ; ++i) {
661 <        for (int j = i + 1; j < atoms.size(); ++j) {
656 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
657 >         inversion = mol->nextInversion(inversionIter)) {
658 >
659 >      a = inversion->getAtomA()->getGlobalIndex();
660 >      b = inversion->getAtomB()->getGlobalIndex();        
661 >      c = inversion->getAtomC()->getGlobalIndex();        
662 >      d = inversion->getAtomD()->getGlobalIndex();        
663 >
664 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
665 >        oneTwoInteractions_.removePair(a, b);      
666 >        oneTwoInteractions_.removePair(a, c);
667 >        oneTwoInteractions_.removePair(a, d);
668 >      } else {
669 >        excludedInteractions_.removePair(a, b);
670 >        excludedInteractions_.removePair(a, c);
671 >        excludedInteractions_.removePair(a, d);
672 >      }
673 >
674 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
675 >        oneThreeInteractions_.removePair(b, c);    
676 >        oneThreeInteractions_.removePair(b, d);    
677 >        oneThreeInteractions_.removePair(c, d);      
678 >      } else {
679 >        excludedInteractions_.removePair(b, c);
680 >        excludedInteractions_.removePair(b, d);
681 >        excludedInteractions_.removePair(c, d);
682 >      }
683 >    }
684 >
685 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
686 >         rb = mol->nextRigidBody(rbIter)) {
687 >      vector<Atom*> atoms = rb->getAtoms();
688 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
689 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
690            a = atoms[i]->getGlobalIndex();
691            b = atoms[j]->getGlobalIndex();
692 <          exclude_.removePair(a, b);
692 >          excludedInteractions_.removePair(a, b);
693          }
694        }
695      }        
696 <
696 >    
697    }
698 <
699 <
698 >  
699 >  
700    void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
701      int curStampId;
702 <
702 >    
703      //index from 0
704      curStampId = moleculeStamps_.size();
705  
# Line 449 | Line 707 | namespace oopse {
707      molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
708    }
709  
452  void SimInfo::update() {
710  
711 <    setupSimType();
712 <
713 < #ifdef IS_MPI
714 <    setupFortranParallel();
715 < #endif
716 <
717 <    setupFortranSim();
718 <
719 <    //setup fortran force field
463 <    /** @deprecate */    
464 <    int isError = 0;
465 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
466 <    if(isError){
467 <      sprintf( painCave.errMsg,
468 <               "ForceField error: There was an error initializing the forceField in fortran.\n" );
469 <      painCave.isFatal = 1;
470 <      simError();
471 <    }
472 <  
473 <    
474 <    setupCutoff();
475 <
711 >  /**
712 >   * update
713 >   *
714 >   *  Performs the global checks and variable settings after the
715 >   *  objects have been created.
716 >   *
717 >   */
718 >  void SimInfo::update() {  
719 >    setupSimVariables();
720      calcNdf();
721      calcNdfRaw();
722      calcNdfTrans();
479
480    fortranInitialized_ = true;
723    }
724 <
725 <  std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
724 >  
725 >  /**
726 >   * getSimulatedAtomTypes
727 >   *
728 >   * Returns an STL set of AtomType* that are actually present in this
729 >   * simulation.  Must query all processors to assemble this information.
730 >   *
731 >   */
732 >  set<AtomType*> SimInfo::getSimulatedAtomTypes() {
733      SimInfo::MoleculeIterator mi;
734      Molecule* mol;
735      Molecule::AtomIterator ai;
736      Atom* atom;
737 <    std::set<AtomType*> atomTypes;
738 <
737 >    set<AtomType*> atomTypes;
738 >    
739      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
740 <
741 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
740 >      for(atom = mol->beginAtom(ai); atom != NULL;
741 >          atom = mol->nextAtom(ai)) {
742          atomTypes.insert(atom->getAtomType());
743 <      }
744 <        
743 >      }      
744 >    }    
745 >    
746 > #ifdef IS_MPI
747 >
748 >    // loop over the found atom types on this processor, and add their
749 >    // numerical idents to a vector:
750 >    
751 >    vector<int> foundTypes;
752 >    set<AtomType*>::iterator i;
753 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i)
754 >      foundTypes.push_back( (*i)->getIdent() );
755 >
756 >    // count_local holds the number of found types on this processor
757 >    int count_local = foundTypes.size();
758 >
759 >    int nproc;
760 >    MPI_Comm_size( MPI_COMM_WORLD, &nproc);
761 >    // int nproc = MPI::COMM_WORLD.Get_size();
762 >
763 >    // we need arrays to hold the counts and displacement vectors for
764 >    // all processors
765 >    vector<int> counts(nproc, 0);
766 >    vector<int> disps(nproc, 0);
767 >
768 >    // fill the counts array
769 >    MPI_Allgather(&count_local, 1, MPI_INT, &counts[0],
770 >                  1, MPI_INT, MPI_COMM_WORLD);
771 >    // MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0],
772 >    //                           1, MPI::INT);
773 >  
774 >    // use the processor counts to compute the displacement array
775 >    disps[0] = 0;    
776 >    int totalCount = counts[0];
777 >    for (int iproc = 1; iproc < nproc; iproc++) {
778 >      disps[iproc] = disps[iproc-1] + counts[iproc-1];
779 >      totalCount += counts[iproc];
780      }
781  
782 +    // we need a (possibly redundant) set of all found types:
783 +    vector<int> ftGlobal(totalCount);
784 +    
785 +    // now spray out the foundTypes to all the other processors:    
786 +    MPI_Allgatherv(&foundTypes[0], count_local, MPI_INT,
787 +                   &ftGlobal[0], &counts[0], &disps[0],
788 +                   MPI_INT, MPI_COMM_WORLD);
789 +    // MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT,
790 +    //                            &ftGlobal[0], &counts[0], &disps[0],
791 +    //                            MPI::INT);
792 +
793 +    vector<int>::iterator j;
794 +
795 +    // foundIdents is a stl set, so inserting an already found ident
796 +    // will have no effect.
797 +    set<int> foundIdents;
798 +
799 +    for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j)
800 +      foundIdents.insert((*j));
801 +    
802 +    // now iterate over the foundIdents and get the actual atom types
803 +    // that correspond to these:
804 +    set<int>::iterator it;
805 +    for (it = foundIdents.begin(); it != foundIdents.end(); ++it)
806 +      atomTypes.insert( forceField_->getAtomType((*it)) );
807 +
808 + #endif
809 +
810      return atomTypes;        
811    }
812  
501  void SimInfo::setupSimType() {
502    std::set<AtomType*>::iterator i;
503    std::set<AtomType*> atomTypes;
504    atomTypes = getUniqueAtomTypes();
505    
506    int useLennardJones = 0;
507    int useElectrostatic = 0;
508    int useEAM = 0;
509    int useCharge = 0;
510    int useDirectional = 0;
511    int useDipole = 0;
512    int useGayBerne = 0;
513    int useSticky = 0;
514    int useShape = 0;
515    int useFLARB = 0; //it is not in AtomType yet
516    int useDirectionalAtom = 0;    
517    int useElectrostatics = 0;
518    //usePBC and useRF are from simParams
519    int usePBC = simParams_->getPBC();
520    int useRF = simParams_->getUseRF();
813  
814 +  int getGlobalCountOfType(AtomType* atype) {
815 +    /*
816 +    set<AtomType*> atypes = getSimulatedAtomTypes();
817 +    map<AtomType*, int> counts_;
818 +
819 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
820 +      for(atom = mol->beginAtom(ai); atom != NULL;
821 +          atom = mol->nextAtom(ai)) {
822 +        atom->getAtomType();
823 +      }      
824 +    }    
825 +    */
826 +    return 0;
827 +  }
828 +
829 +  void SimInfo::setupSimVariables() {
830 +    useAtomicVirial_ = simParams_->getUseAtomicVirial();
831 +    // we only call setAccumulateBoxDipole if the accumulateBoxDipole
832 +    // parameter is true
833 +    calcBoxDipole_ = false;
834 +    if ( simParams_->haveAccumulateBoxDipole() )
835 +      if ( simParams_->getAccumulateBoxDipole() ) {
836 +        calcBoxDipole_ = true;      
837 +      }
838 +    
839 +    set<AtomType*>::iterator i;
840 +    set<AtomType*> atomTypes;
841 +    atomTypes = getSimulatedAtomTypes();    
842 +    bool usesElectrostatic = false;
843 +    bool usesMetallic = false;
844 +    bool usesDirectional = false;
845 +    bool usesFluctuatingCharges =  false;
846      //loop over all of the atom types
847      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
848 <      useLennardJones |= (*i)->isLennardJones();
849 <      useElectrostatic |= (*i)->isElectrostatic();
850 <      useEAM |= (*i)->isEAM();
851 <      useCharge |= (*i)->isCharge();
528 <      useDirectional |= (*i)->isDirectional();
529 <      useDipole |= (*i)->isDipole();
530 <      useGayBerne |= (*i)->isGayBerne();
531 <      useSticky |= (*i)->isSticky();
532 <      useShape |= (*i)->isShape();
848 >      usesElectrostatic |= (*i)->isElectrostatic();
849 >      usesMetallic |= (*i)->isMetal();
850 >      usesDirectional |= (*i)->isDirectional();
851 >      usesFluctuatingCharges |= (*i)->isFluctuatingCharge();
852      }
853  
854 <    if (useSticky || useDipole || useGayBerne || useShape) {
536 <      useDirectionalAtom = 1;
537 <    }
538 <
539 <    if (useCharge || useDipole) {
540 <      useElectrostatics = 1;
541 <    }
542 <
543 < #ifdef IS_MPI    
854 > #ifdef IS_MPI
855      int temp;
856  
857 <    temp = usePBC;
858 <    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
857 >    temp = usesDirectional;
858 >    MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD);
859 >    usesDirectionalAtoms_ = (temp == 0) ? false : true;
860  
861 <    temp = useDirectionalAtom;
862 <    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
861 >    // MPI::COMM_WORLD.Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI::BOOL,
862 >    //                           MPI::LOR);
863 >    
864 >    temp = usesMetallic;
865 >    MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD);
866 >    usesMetallicAtoms_ = (temp == 0) ? false : true;
867  
868 <    temp = useLennardJones;
869 <    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
868 >    // MPI::COMM_WORLD.Allreduce(&temp, &usesMetallicAtoms_, 1, MPI::BOOL,
869 >    //                           MPI::LOR);
870 >    
871 >    temp = usesElectrostatic;
872 >    MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD);
873 >    usesElectrostaticAtoms_ = (temp == 0) ? false : true;
874  
875 <    temp = useElectrostatics;
876 <    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
875 >    // MPI::COMM_WORLD.Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI::BOOL,
876 >    //                           MPI::LOR);
877  
878 <    temp = useCharge;
879 <    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
878 >    temp = usesFluctuatingCharges;
879 >    MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD);
880 >    usesFluctuatingCharges_ = (temp == 0) ? false : true;
881  
882 <    temp = useDipole;
883 <    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
882 >    // MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL,
883 >    //                           MPI::LOR);
884  
885 <    temp = useSticky;
565 <    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
885 > #else
886  
887 <    temp = useGayBerne;
888 <    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
887 >    usesDirectionalAtoms_ = usesDirectional;
888 >    usesMetallicAtoms_ = usesMetallic;
889 >    usesElectrostaticAtoms_ = usesElectrostatic;
890 >    usesFluctuatingCharges_ = usesFluctuatingCharges;
891  
892 <    temp = useEAM;
893 <    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
892 > #endif
893 >    
894 >    requiresPrepair_ = usesMetallicAtoms_ ? true : false;
895 >    requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false;
896 >    requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false;    
897 >  }
898  
573    temp = useShape;
574    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
899  
900 <    temp = useFLARB;
901 <    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
900 >  vector<int> SimInfo::getGlobalAtomIndices() {
901 >    SimInfo::MoleculeIterator mi;
902 >    Molecule* mol;
903 >    Molecule::AtomIterator ai;
904 >    Atom* atom;
905  
906 <    temp = useRF;
580 <    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
906 >    vector<int> GlobalAtomIndices(getNAtoms(), 0);
907      
908 < #endif
908 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
909 >      
910 >      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
911 >        GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex();
912 >      }
913 >    }
914 >    return GlobalAtomIndices;
915 >  }
916  
584    fInfo_.SIM_uses_PBC = usePBC;    
585    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
586    fInfo_.SIM_uses_LennardJones = useLennardJones;
587    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
588    fInfo_.SIM_uses_Charges = useCharge;
589    fInfo_.SIM_uses_Dipoles = useDipole;
590    fInfo_.SIM_uses_Sticky = useSticky;
591    fInfo_.SIM_uses_GayBerne = useGayBerne;
592    fInfo_.SIM_uses_EAM = useEAM;
593    fInfo_.SIM_uses_Shapes = useShape;
594    fInfo_.SIM_uses_FLARB = useFLARB;
595    fInfo_.SIM_uses_RF = useRF;
917  
918 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
918 >  vector<int> SimInfo::getGlobalGroupIndices() {
919 >    SimInfo::MoleculeIterator mi;
920 >    Molecule* mol;
921 >    Molecule::CutoffGroupIterator ci;
922 >    CutoffGroup* cg;
923  
924 <      if (simParams_->haveDielectric()) {
925 <        fInfo_.dielect = simParams_->getDielectric();
926 <      } else {
927 <        sprintf(painCave.errMsg,
928 <                "SimSetup Error: No Dielectric constant was set.\n"
929 <                "\tYou are trying to use Reaction Field without"
930 <                "\tsetting a dielectric constant!\n");
931 <        painCave.isFatal = 1;
932 <        simError();
933 <      }
609 <        
610 <    } else {
611 <      fInfo_.dielect = 0.0;
924 >    vector<int> GlobalGroupIndices;
925 >    
926 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
927 >      
928 >      //local index of cutoff group is trivial, it only depends on the
929 >      //order of travesing
930 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
931 >           cg = mol->nextCutoffGroup(ci)) {
932 >        GlobalGroupIndices.push_back(cg->getGlobalIndex());
933 >      }        
934      }
935 <
935 >    return GlobalGroupIndices;
936    }
937  
616  void SimInfo::setupFortranSim() {
617    int isError;
618    int nExclude;
619    std::vector<int> fortranGlobalGroupMembership;
620    
621    nExclude = exclude_.getSize();
622    isError = 0;
938  
939 <    //globalGroupMembership_ is filled by SimCreator    
625 <    for (int i = 0; i < nGlobalAtoms_; i++) {
626 <      fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
627 <    }
939 >  void SimInfo::prepareTopology() {
940  
941      //calculate mass ratio of cutoff group
630    std::vector<double> mfact;
942      SimInfo::MoleculeIterator mi;
943      Molecule* mol;
944      Molecule::CutoffGroupIterator ci;
945      CutoffGroup* cg;
946      Molecule::AtomIterator ai;
947      Atom* atom;
948 <    double totalMass;
948 >    RealType totalMass;
949  
950 <    //to avoid memory reallocation, reserve enough space for mfact
951 <    mfact.reserve(getNCutoffGroups());
950 >    /**
951 >     * The mass factor is the relative mass of an atom to the total
952 >     * mass of the cutoff group it belongs to.  By default, all atoms
953 >     * are their own cutoff groups, and therefore have mass factors of
954 >     * 1.  We need some special handling for massless atoms, which
955 >     * will be treated as carrying the entire mass of the cutoff
956 >     * group.
957 >     */
958 >    massFactors_.clear();
959 >    massFactors_.resize(getNAtoms(), 1.0);
960      
961      for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
962 <      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
962 >      for (cg = mol->beginCutoffGroup(ci); cg != NULL;
963 >           cg = mol->nextCutoffGroup(ci)) {
964  
965          totalMass = cg->getMass();
966          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
967 <          mfact.push_back(atom->getMass()/totalMass);
967 >          // Check for massless groups - set mfact to 1 if true
968 >          if (totalMass != 0)
969 >            massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass;
970 >          else
971 >            massFactors_[atom->getLocalIndex()] = 1.0;
972          }
649
973        }      
974      }
975  
976 <    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
654 <    std::vector<int> identArray;
976 >    // Build the identArray_ and regions_
977  
978 <    //to avoid memory reallocation, reserve enough space identArray
979 <    identArray.reserve(getNAtoms());
980 <    
981 <    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
978 >    identArray_.clear();
979 >    identArray_.reserve(getNAtoms());  
980 >    regions_.clear();
981 >    regions_.reserve(getNAtoms());
982 >
983 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {      
984 >      int reg = mol->getRegion();      
985        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
986 <        identArray.push_back(atom->getIdent());
986 >        identArray_.push_back(atom->getIdent());
987 >        regions_.push_back(reg);
988        }
989      }    
990 <
991 <    //fill molMembershipArray
666 <    //molMembershipArray is filled by SimCreator    
667 <    std::vector<int> molMembershipArray(nGlobalAtoms_);
668 <    for (int i = 0; i < nGlobalAtoms_; i++) {
669 <      molMembershipArray[i] = globalMolMembership_[i] + 1;
670 <    }
671 <    
672 <    //setup fortran simulation
673 <    int nGlobalExcludes = 0;
674 <    int* globalExcludes = NULL;
675 <    int* excludeList = exclude_.getExcludeList();
676 <    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
677 <                   &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
678 <                   &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
679 <
680 <    if( isError ){
681 <
682 <      sprintf( painCave.errMsg,
683 <               "There was an error setting the simulation information in fortran.\n" );
684 <      painCave.isFatal = 1;
685 <      painCave.severity = OOPSE_ERROR;
686 <      simError();
687 <    }
688 <
689 < #ifdef IS_MPI
690 <    sprintf( checkPointMsg,
691 <             "succesfully sent the simulation information to fortran.\n");
692 <    MPIcheckPoint();
693 < #endif // is_mpi
990 >      
991 >    topologyDone_ = true;
992    }
993  
696
697 #ifdef IS_MPI
698  void SimInfo::setupFortranParallel() {
699    
700    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
701    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
702    std::vector<int> localToGlobalCutoffGroupIndex;
703    SimInfo::MoleculeIterator mi;
704    Molecule::AtomIterator ai;
705    Molecule::CutoffGroupIterator ci;
706    Molecule* mol;
707    Atom* atom;
708    CutoffGroup* cg;
709    mpiSimData parallelData;
710    int isError;
711
712    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
713
714      //local index(index in DataStorge) of atom is important
715      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
716        localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
717      }
718
719      //local index of cutoff group is trivial, it only depends on the order of travesing
720      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
721        localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
722      }        
723        
724    }
725
726    //fill up mpiSimData struct
727    parallelData.nMolGlobal = getNGlobalMolecules();
728    parallelData.nMolLocal = getNMolecules();
729    parallelData.nAtomsGlobal = getNGlobalAtoms();
730    parallelData.nAtomsLocal = getNAtoms();
731    parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
732    parallelData.nGroupsLocal = getNCutoffGroups();
733    parallelData.myNode = worldRank;
734    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors));
735
736    //pass mpiSimData struct and index arrays to fortran
737    setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
738                    &localToGlobalAtomIndex[0],  &(parallelData.nGroupsLocal),
739                    &localToGlobalCutoffGroupIndex[0], &isError);
740
741    if (isError) {
742      sprintf(painCave.errMsg,
743              "mpiRefresh errror: fortran didn't like something we gave it.\n");
744      painCave.isFatal = 1;
745      simError();
746    }
747
748    sprintf(checkPointMsg, " mpiRefresh successful.\n");
749    MPIcheckPoint();
750
751
752  }
753
754 #endif
755
756  double SimInfo::calcMaxCutoffRadius() {
757
758
759    std::set<AtomType*> atomTypes;
760    std::set<AtomType*>::iterator i;
761    std::vector<double> cutoffRadius;
762
763    //get the unique atom types
764    atomTypes = getUniqueAtomTypes();
765
766    //query the max cutoff radius among these atom types
767    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
768      cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
769    }
770
771    double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
772 #ifdef IS_MPI
773    //pick the max cutoff radius among the processors
774 #endif
775
776    return maxCutoffRadius;
777  }
778
779  void SimInfo::getCutoff(double& rcut, double& rsw) {
780    
781    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
782        
783      if (!simParams_->haveRcut()){
784        sprintf(painCave.errMsg,
785                "SimCreator Warning: No value was set for the cutoffRadius.\n"
786                "\tOOPSE will use a default value of 15.0 angstroms"
787                "\tfor the cutoffRadius.\n");
788        painCave.isFatal = 0;
789        simError();
790        rcut = 15.0;
791      } else{
792        rcut = simParams_->getRcut();
793      }
794
795      if (!simParams_->haveRsw()){
796        sprintf(painCave.errMsg,
797                "SimCreator Warning: No value was set for switchingRadius.\n"
798                "\tOOPSE will use a default value of\n"
799                "\t0.95 * cutoffRadius for the switchingRadius\n");
800        painCave.isFatal = 0;
801        simError();
802        rsw = 0.95 * rcut;
803      } else{
804        rsw = simParams_->getRsw();
805      }
806
807    } else {
808      // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
809      //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
810        
811      if (simParams_->haveRcut()) {
812        rcut = simParams_->getRcut();
813      } else {
814        //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
815        rcut = calcMaxCutoffRadius();
816      }
817
818      if (simParams_->haveRsw()) {
819        rsw  = simParams_->getRsw();
820      } else {
821        rsw = rcut;
822      }
823    
824    }
825  }
826
827  void SimInfo::setupCutoff() {
828    getCutoff(rcut_, rsw_);    
829    double rnblist = rcut_ + 1; // skin of neighbor list
830
831    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
832    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
833  }
834
994    void SimInfo::addProperty(GenericData* genData) {
995      properties_.addProperty(genData);  
996    }
997  
998 <  void SimInfo::removeProperty(const std::string& propName) {
998 >  void SimInfo::removeProperty(const string& propName) {
999      properties_.removeProperty(propName);  
1000    }
1001  
# Line 844 | Line 1003 | namespace oopse {
1003      properties_.clearProperties();
1004    }
1005  
1006 <  std::vector<std::string> SimInfo::getPropertyNames() {
1006 >  vector<string> SimInfo::getPropertyNames() {
1007      return properties_.getPropertyNames();  
1008    }
1009        
1010 <  std::vector<GenericData*> SimInfo::getProperties() {
1010 >  vector<GenericData*> SimInfo::getProperties() {
1011      return properties_.getProperties();
1012    }
1013  
1014 <  GenericData* SimInfo::getPropertyByName(const std::string& propName) {
1014 >  GenericData* SimInfo::getPropertyByName(const string& propName) {
1015      return properties_.getPropertyByName(propName);
1016    }
1017  
# Line 863 | Line 1022 | namespace oopse {
1022      delete sman_;
1023      sman_ = sman;
1024  
866    Molecule* mol;
867    RigidBody* rb;
868    Atom* atom;
1025      SimInfo::MoleculeIterator mi;
1026 +    Molecule::AtomIterator ai;
1027      Molecule::RigidBodyIterator rbIter;
1028 <    Molecule::AtomIterator atomIter;;
1028 >    Molecule::CutoffGroupIterator cgIter;
1029 >    Molecule::BondIterator bondIter;
1030 >    Molecule::BendIterator bendIter;
1031 >    Molecule::TorsionIterator torsionIter;
1032 >    Molecule::InversionIterator inversionIter;
1033  
1034 +    Molecule* mol;
1035 +    Atom* atom;
1036 +    RigidBody* rb;
1037 +    CutoffGroup* cg;
1038 +    Bond* bond;
1039 +    Bend* bend;
1040 +    Torsion* torsion;
1041 +    Inversion* inversion;    
1042 +
1043      for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
1044          
1045 <      for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
1045 >      for (atom = mol->beginAtom(ai); atom != NULL;
1046 >           atom = mol->nextAtom(ai)) {
1047          atom->setSnapshotManager(sman_);
1048 <      }
1049 <        
1050 <      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
1048 >      }        
1049 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
1050 >           rb = mol->nextRigidBody(rbIter)) {
1051          rb->setSnapshotManager(sman_);
1052        }
1053 <    }    
1054 <    
1053 >      for (cg = mol->beginCutoffGroup(cgIter); cg != NULL;
1054 >           cg = mol->nextCutoffGroup(cgIter)) {
1055 >        cg->setSnapshotManager(sman_);
1056 >      }
1057 >      for (bond = mol->beginBond(bondIter); bond != NULL;
1058 >           bond = mol->nextBond(bondIter)) {
1059 >        bond->setSnapshotManager(sman_);
1060 >      }
1061 >      for (bend = mol->beginBend(bendIter); bend != NULL;
1062 >           bend = mol->nextBend(bendIter)) {
1063 >        bend->setSnapshotManager(sman_);
1064 >      }
1065 >      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
1066 >           torsion = mol->nextTorsion(torsionIter)) {
1067 >        torsion->setSnapshotManager(sman_);
1068 >      }
1069 >      for (inversion = mol->beginInversion(inversionIter); inversion != NULL;
1070 >           inversion = mol->nextInversion(inversionIter)) {
1071 >        inversion->setSnapshotManager(sman_);
1072 >      }
1073 >    }
1074    }
1075  
886  Vector3d SimInfo::getComVel(){
887    SimInfo::MoleculeIterator i;
888    Molecule* mol;
1076  
1077 <    Vector3d comVel(0.0);
891 <    double totalMass = 0.0;
892 <    
893 <
894 <    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
895 <      double mass = mol->getMass();
896 <      totalMass += mass;
897 <      comVel += mass * mol->getComVel();
898 <    }  
1077 >  ostream& operator <<(ostream& o, SimInfo& info) {
1078  
1079 < #ifdef IS_MPI
901 <    double tmpMass = totalMass;
902 <    Vector3d tmpComVel(comVel);    
903 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
904 <    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
905 < #endif
906 <
907 <    comVel /= totalMass;
908 <
909 <    return comVel;
1079 >    return o;
1080    }
1081 +  
1082 +  
1083 +  StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
1084 +    if (index >= int(IOIndexToIntegrableObject.size())) {
1085 +      sprintf(painCave.errMsg,
1086 +              "SimInfo::getIOIndexToIntegrableObject Error: Integrable Object\n"
1087 +              "\tindex exceeds number of known objects!\n");
1088 +      painCave.isFatal = 1;
1089 +      simError();
1090 +      return NULL;
1091 +    } else
1092 +      return IOIndexToIntegrableObject.at(index);
1093 +  }
1094 +  
1095 +  void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) {
1096 +    IOIndexToIntegrableObject= v;
1097 +  }
1098  
1099 <  Vector3d SimInfo::getCom(){
1100 <    SimInfo::MoleculeIterator i;
914 <    Molecule* mol;
915 <
916 <    Vector3d com(0.0);
917 <    double totalMass = 0.0;
918 <    
919 <    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
920 <      double mass = mol->getMass();
921 <      totalMass += mass;
922 <      com += mass * mol->getCom();
923 <    }  
924 <
1099 >  int SimInfo::getNGlobalConstraints() {
1100 >    int nGlobalConstraints;
1101   #ifdef IS_MPI
1102 <    double tmpMass = totalMass;
1103 <    Vector3d tmpCom(com);    
1104 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1105 <    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1102 >    MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1,  
1103 >                              MPI_INT, MPI_SUM, MPI_COMM_WORLD);
1104 >    // MPI::COMM_WORLD.Allreduce(&nConstraints_, &nGlobalConstraints, 1,
1105 >    //                           MPI::INT, MPI::SUM);
1106 > #else
1107 >    nGlobalConstraints =  nConstraints_;
1108   #endif
1109 <
932 <    com /= totalMass;
933 <
934 <    return com;
935 <
936 <  }        
937 <
938 <  std::ostream& operator <<(std::ostream& o, SimInfo& info) {
939 <
940 <    return o;
1109 >    return nGlobalConstraints;
1110    }
1111  
1112 < }//end namespace oopse
1112 > }//end namespace OpenMD
1113  

Comparing trunk/src/brains/SimInfo.cpp (property svn:keywords):
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines