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 586 by gezelter, Wed Sep 7 20:46:46 2005 UTC vs.
Revision 1386 by gezelter, Fri Oct 23 18:41:09 2009 UTC

# Line 48 | Line 48
48  
49   #include <algorithm>
50   #include <set>
51 + #include <map>
52  
53   #include "brains/SimInfo.hpp"
54   #include "math/Vector3.hpp"
55   #include "primitives/Molecule.hpp"
56 + #include "primitives/StuntDouble.hpp"
57   #include "UseTheForce/fCutoffPolicy.h"
58 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
59 + #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
60 + #include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
61   #include "UseTheForce/doForces_interface.h"
62 < #include "UseTheForce/notifyCutoffs_interface.h"
62 > #include "UseTheForce/DarkSide/neighborLists_interface.h"
63 > #include "UseTheForce/DarkSide/electrostatic_interface.h"
64 > #include "UseTheForce/DarkSide/switcheroo_interface.h"
65   #include "utils/MemoryUtils.hpp"
66   #include "utils/simError.h"
67   #include "selection/SelectionManager.hpp"
68 + #include "io/ForceFieldOptions.hpp"
69 + #include "UseTheForce/ForceField.hpp"
70  
71 +
72   #ifdef IS_MPI
73   #include "UseTheForce/mpiComponentPlan.h"
74   #include "UseTheForce/DarkSide/simParallel_interface.h"
75   #endif
76  
77   namespace oopse {
78 +  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 <  SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
86 <                   ForceField* ff, Globals* simParams) :
87 <    stamps_(stamps), forceField_(ff), simParams_(simParams),
88 <    ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
85 >    return result;
86 >  }
87 >  
88 >  SimInfo::SimInfo(ForceField* ff, Globals* simParams) :
89 >    forceField_(ff), simParams_(simParams),
90 >    ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
91      nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
92      nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
93 <    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
94 <    nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
95 <    sman_(NULL), fortranInitialized_(false) {
93 >    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nInversions_(0),
94 >    nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0),
95 >    nConstraints_(0), sman_(NULL), fortranInitialized_(false),
96 >    calcBoxDipole_(false), useAtomicVirial_(true) {
97  
98 <            
80 <      std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
98 >
99        MoleculeStamp* molStamp;
100        int nMolWithSameStamp;
101        int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
102 <      int nGroups = 0;          //total cutoff groups defined in meta-data file
102 >      int nGroups = 0;      //total cutoff groups defined in meta-data file
103        CutoffGroupStamp* cgStamp;    
104        RigidBodyStamp* rbStamp;
105        int nRigidAtoms = 0;
106 <    
107 <      for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
108 <        molStamp = i->first;
109 <        nMolWithSameStamp = i->second;
106 >
107 >      std::vector<Component*> components = simParams->getComponents();
108 >      
109 >      for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
110 >        molStamp = (*i)->getMoleculeStamp();
111 >        nMolWithSameStamp = (*i)->getNMol();
112          
113          addMoleculeStamp(molStamp, nMolWithSameStamp);
114  
115          //calculate atoms in molecules
116          nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
117  
98 –
118          //calculate atoms in cutoff groups
119          int nAtomsInGroups = 0;
120          int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
121          
122          for (int j=0; j < nCutoffGroupsInStamp; j++) {
123 <          cgStamp = molStamp->getCutoffGroup(j);
123 >          cgStamp = molStamp->getCutoffGroupStamp(j);
124            nAtomsInGroups += cgStamp->getNMembers();
125          }
126  
127          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
128 +
129          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
130  
131          //calculate atoms in rigid bodies
# Line 113 | Line 133 | namespace oopse {
133          int nRigidBodiesInStamp = molStamp->getNRigidBodies();
134          
135          for (int j=0; j < nRigidBodiesInStamp; j++) {
136 <          rbStamp = molStamp->getRigidBody(j);
136 >          rbStamp = molStamp->getRigidBodyStamp(j);
137            nAtomsInRigidBodies += rbStamp->getNMembers();
138          }
139  
# Line 122 | Line 142 | namespace oopse {
142          
143        }
144  
145 <      //every free atom (atom does not belong to cutoff groups) is a cutoff group
146 <      //therefore the total number of cutoff groups in the system is equal to
147 <      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
148 <      //file plus the number of cutoff groups defined in meta-data file
145 >      //every free atom (atom does not belong to cutoff groups) is a cutoff
146 >      //group therefore the total number of cutoff groups in the system is
147 >      //equal to the total number of atoms minus number of atoms belong to
148 >      //cutoff group defined in meta-data file plus the number of cutoff
149 >      //groups defined in meta-data file
150        nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
151  
152 <      //every free atom (atom does not belong to rigid bodies) is an integrable object
153 <      //therefore the total number of  integrable objects in the system is equal to
154 <      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
155 <      //file plus the number of  rigid bodies defined in meta-data file
156 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
157 <
152 >      //every free atom (atom does not belong to rigid bodies) is an
153 >      //integrable object therefore the total number of integrable objects
154 >      //in the system is equal to the total number of atoms minus number of
155 >      //atoms belong to rigid body defined in meta-data file plus the number
156 >      //of rigid bodies defined in meta-data file
157 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
158 >                                                + nGlobalRigidBodies_;
159 >  
160        nGlobalMols_ = molStampIds_.size();
138 –
139 – #ifdef IS_MPI    
161        molToProcMap_.resize(nGlobalMols_);
141 – #endif
142 –
162      }
163  
164    SimInfo::~SimInfo() {
# Line 149 | Line 168 | namespace oopse {
168      }
169      molecules_.clear();
170        
152 –    delete stamps_;
171      delete sman_;
172      delete simParams_;
173      delete forceField_;
# Line 178 | Line 196 | namespace oopse {
196        nBonds_ += mol->getNBonds();
197        nBends_ += mol->getNBends();
198        nTorsions_ += mol->getNTorsions();
199 +      nInversions_ += mol->getNInversions();
200        nRigidBodies_ += mol->getNRigidBodies();
201        nIntegrableObjects_ += mol->getNIntegrableObjects();
202        nCutoffGroups_ += mol->getNCutoffGroups();
203        nConstraints_ += mol->getNConstraintPairs();
204  
205 <      addExcludePairs(mol);
206 <        
205 >      addInteractionPairs(mol);
206 >  
207        return true;
208      } else {
209        return false;
# Line 203 | Line 222 | namespace oopse {
222        nBonds_ -= mol->getNBonds();
223        nBends_ -= mol->getNBends();
224        nTorsions_ -= mol->getNTorsions();
225 +      nInversions_ -= mol->getNInversions();
226        nRigidBodies_ -= mol->getNRigidBodies();
227        nIntegrableObjects_ -= mol->getNIntegrableObjects();
228        nCutoffGroups_ -= mol->getNCutoffGroups();
229        nConstraints_ -= mol->getNConstraintPairs();
230  
231 <      removeExcludePairs(mol);
231 >      removeInteractionPairs(mol);
232        molecules_.erase(mol->getGlobalIndex());
233  
234        delete mol;
# Line 256 | Line 276 | namespace oopse {
276            }
277          }
278              
279 <      }//end for (integrableObject)
280 <    }// end for (mol)
279 >      }
280 >    }
281      
282      // n_constraints is local, so subtract them on each processor
283      ndf_local -= nConstraints_;
# Line 274 | Line 294 | namespace oopse {
294  
295    }
296  
297 +  int SimInfo::getFdf() {
298 + #ifdef IS_MPI
299 +    MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
300 + #else
301 +    fdf_ = fdf_local;
302 + #endif
303 +    return fdf_;
304 +  }
305 +    
306    void SimInfo::calcNdfRaw() {
307      int ndfRaw_local;
308  
# Line 325 | Line 354 | namespace oopse {
354  
355    }
356  
357 <  void SimInfo::addExcludePairs(Molecule* mol) {
357 >  void SimInfo::addInteractionPairs(Molecule* mol) {
358 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
359      std::vector<Bond*>::iterator bondIter;
360      std::vector<Bend*>::iterator bendIter;
361      std::vector<Torsion*>::iterator torsionIter;
362 +    std::vector<Inversion*>::iterator inversionIter;
363      Bond* bond;
364      Bend* bend;
365      Torsion* torsion;
366 +    Inversion* inversion;
367      int a;
368      int b;
369      int c;
370      int d;
371 +
372 +    // atomGroups can be used to add special interaction maps between
373 +    // groups of atoms that are in two separate rigid bodies.
374 +    // However, most site-site interactions between two rigid bodies
375 +    // are probably not special, just the ones between the physically
376 +    // bonded atoms.  Interactions *within* a single rigid body should
377 +    // always be excluded.  These are done at the bottom of this
378 +    // function.
379 +
380 +    std::map<int, std::set<int> > atomGroups;
381 +    Molecule::RigidBodyIterator rbIter;
382 +    RigidBody* rb;
383 +    Molecule::IntegrableObjectIterator ii;
384 +    StuntDouble* integrableObject;
385      
386 <    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);
386 >    for (integrableObject = mol->beginIntegrableObject(ii);
387 >         integrableObject != NULL;
388 >         integrableObject = mol->nextIntegrableObject(ii)) {
389 >      
390 >      if (integrableObject->isRigidBody()) {
391 >        rb = static_cast<RigidBody*>(integrableObject);
392 >        std::vector<Atom*> atoms = rb->getAtoms();
393 >        std::set<int> rigidAtoms;
394 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
395 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
396 >        }
397 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
398 >          atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
399 >        }      
400 >      } else {
401 >        std::set<int> oneAtomSet;
402 >        oneAtomSet.insert(integrableObject->getGlobalIndex());
403 >        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
404 >      }
405 >    }  
406 >          
407 >    for (bond= mol->beginBond(bondIter); bond != NULL;
408 >         bond = mol->nextBond(bondIter)) {
409 >
410 >      a = bond->getAtomA()->getGlobalIndex();
411 >      b = bond->getAtomB()->getGlobalIndex();  
412 >    
413 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
414 >        oneTwoInteractions_.addPair(a, b);
415 >      } else {
416 >        excludedInteractions_.addPair(a, b);
417 >      }
418      }
419  
420 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
420 >    for (bend= mol->beginBend(bendIter); bend != NULL;
421 >         bend = mol->nextBend(bendIter)) {
422 >
423        a = bend->getAtomA()->getGlobalIndex();
424        b = bend->getAtomB()->getGlobalIndex();        
425        c = bend->getAtomC()->getGlobalIndex();
426 +      
427 +      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
428 +        oneTwoInteractions_.addPair(a, b);      
429 +        oneTwoInteractions_.addPair(b, c);
430 +      } else {
431 +        excludedInteractions_.addPair(a, b);
432 +        excludedInteractions_.addPair(b, c);
433 +      }
434  
435 <      exclude_.addPair(a, b);
436 <      exclude_.addPair(a, c);
437 <      exclude_.addPair(b, c);        
435 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
436 >        oneThreeInteractions_.addPair(a, c);      
437 >      } else {
438 >        excludedInteractions_.addPair(a, c);
439 >      }
440      }
441  
442 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
442 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
443 >         torsion = mol->nextTorsion(torsionIter)) {
444 >
445        a = torsion->getAtomA()->getGlobalIndex();
446        b = torsion->getAtomB()->getGlobalIndex();        
447        c = torsion->getAtomC()->getGlobalIndex();        
448 <      d = torsion->getAtomD()->getGlobalIndex();        
448 >      d = torsion->getAtomD()->getGlobalIndex();      
449  
450 <      exclude_.addPair(a, b);
451 <      exclude_.addPair(a, c);
452 <      exclude_.addPair(a, d);
453 <      exclude_.addPair(b, c);
454 <      exclude_.addPair(b, d);
455 <      exclude_.addPair(c, d);        
450 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
451 >        oneTwoInteractions_.addPair(a, b);      
452 >        oneTwoInteractions_.addPair(b, c);
453 >        oneTwoInteractions_.addPair(c, d);
454 >      } else {
455 >        excludedInteractions_.addPair(a, b);
456 >        excludedInteractions_.addPair(b, c);
457 >        excludedInteractions_.addPair(c, d);
458 >      }
459 >
460 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
461 >        oneThreeInteractions_.addPair(a, c);      
462 >        oneThreeInteractions_.addPair(b, d);      
463 >      } else {
464 >        excludedInteractions_.addPair(a, c);
465 >        excludedInteractions_.addPair(b, d);
466 >      }
467 >
468 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
469 >        oneFourInteractions_.addPair(a, d);      
470 >      } else {
471 >        excludedInteractions_.addPair(a, d);
472 >      }
473      }
474  
475 <    Molecule::RigidBodyIterator rbIter;
476 <    RigidBody* rb;
477 <    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
475 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
476 >         inversion = mol->nextInversion(inversionIter)) {
477 >
478 >      a = inversion->getAtomA()->getGlobalIndex();
479 >      b = inversion->getAtomB()->getGlobalIndex();        
480 >      c = inversion->getAtomC()->getGlobalIndex();        
481 >      d = inversion->getAtomD()->getGlobalIndex();        
482 >
483 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
484 >        oneTwoInteractions_.addPair(a, b);      
485 >        oneTwoInteractions_.addPair(a, c);
486 >        oneTwoInteractions_.addPair(a, d);
487 >      } else {
488 >        excludedInteractions_.addPair(a, b);
489 >        excludedInteractions_.addPair(a, c);
490 >        excludedInteractions_.addPair(a, d);
491 >      }
492 >
493 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
494 >        oneThreeInteractions_.addPair(b, c);    
495 >        oneThreeInteractions_.addPair(b, d);    
496 >        oneThreeInteractions_.addPair(c, d);      
497 >      } else {
498 >        excludedInteractions_.addPair(b, c);
499 >        excludedInteractions_.addPair(b, d);
500 >        excludedInteractions_.addPair(c, d);
501 >      }
502 >    }
503 >
504 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
505 >         rb = mol->nextRigidBody(rbIter)) {
506        std::vector<Atom*> atoms = rb->getAtoms();
507 <      for (int i = 0; i < atoms.size() -1 ; ++i) {
508 <        for (int j = i + 1; j < atoms.size(); ++j) {
507 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
508 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
509            a = atoms[i]->getGlobalIndex();
510            b = atoms[j]->getGlobalIndex();
511 <          exclude_.addPair(a, b);
511 >          excludedInteractions_.addPair(a, b);
512          }
513        }
514      }        
515  
516    }
517  
518 <  void SimInfo::removeExcludePairs(Molecule* mol) {
518 >  void SimInfo::removeInteractionPairs(Molecule* mol) {
519 >    ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
520      std::vector<Bond*>::iterator bondIter;
521      std::vector<Bend*>::iterator bendIter;
522      std::vector<Torsion*>::iterator torsionIter;
523 +    std::vector<Inversion*>::iterator inversionIter;
524      Bond* bond;
525      Bend* bend;
526      Torsion* torsion;
527 +    Inversion* inversion;
528      int a;
529      int b;
530      int c;
531      int d;
532 <    
533 <    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
532 >
533 >    std::map<int, std::set<int> > atomGroups;
534 >    Molecule::RigidBodyIterator rbIter;
535 >    RigidBody* rb;
536 >    Molecule::IntegrableObjectIterator ii;
537 >    StuntDouble* integrableObject;
538 >    
539 >    for (integrableObject = mol->beginIntegrableObject(ii);
540 >         integrableObject != NULL;
541 >         integrableObject = mol->nextIntegrableObject(ii)) {
542 >      
543 >      if (integrableObject->isRigidBody()) {
544 >        rb = static_cast<RigidBody*>(integrableObject);
545 >        std::vector<Atom*> atoms = rb->getAtoms();
546 >        std::set<int> rigidAtoms;
547 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
548 >          rigidAtoms.insert(atoms[i]->getGlobalIndex());
549 >        }
550 >        for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
551 >          atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
552 >        }      
553 >      } else {
554 >        std::set<int> oneAtomSet;
555 >        oneAtomSet.insert(integrableObject->getGlobalIndex());
556 >        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
557 >      }
558 >    }  
559 >
560 >    for (bond= mol->beginBond(bondIter); bond != NULL;
561 >         bond = mol->nextBond(bondIter)) {
562 >      
563        a = bond->getAtomA()->getGlobalIndex();
564 <      b = bond->getAtomB()->getGlobalIndex();        
565 <      exclude_.removePair(a, b);
564 >      b = bond->getAtomB()->getGlobalIndex();  
565 >    
566 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
567 >        oneTwoInteractions_.removePair(a, b);
568 >      } else {
569 >        excludedInteractions_.removePair(a, b);
570 >      }
571      }
572  
573 <    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
573 >    for (bend= mol->beginBend(bendIter); bend != NULL;
574 >         bend = mol->nextBend(bendIter)) {
575 >
576        a = bend->getAtomA()->getGlobalIndex();
577        b = bend->getAtomB()->getGlobalIndex();        
578        c = bend->getAtomC()->getGlobalIndex();
579 +      
580 +      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
581 +        oneTwoInteractions_.removePair(a, b);      
582 +        oneTwoInteractions_.removePair(b, c);
583 +      } else {
584 +        excludedInteractions_.removePair(a, b);
585 +        excludedInteractions_.removePair(b, c);
586 +      }
587  
588 <      exclude_.removePair(a, b);
589 <      exclude_.removePair(a, c);
590 <      exclude_.removePair(b, c);        
588 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
589 >        oneThreeInteractions_.removePair(a, c);      
590 >      } else {
591 >        excludedInteractions_.removePair(a, c);
592 >      }
593      }
594  
595 <    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
595 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
596 >         torsion = mol->nextTorsion(torsionIter)) {
597 >
598        a = torsion->getAtomA()->getGlobalIndex();
599        b = torsion->getAtomB()->getGlobalIndex();        
600        c = torsion->getAtomC()->getGlobalIndex();        
601 <      d = torsion->getAtomD()->getGlobalIndex();        
601 >      d = torsion->getAtomD()->getGlobalIndex();      
602 >  
603 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
604 >        oneTwoInteractions_.removePair(a, b);      
605 >        oneTwoInteractions_.removePair(b, c);
606 >        oneTwoInteractions_.removePair(c, d);
607 >      } else {
608 >        excludedInteractions_.removePair(a, b);
609 >        excludedInteractions_.removePair(b, c);
610 >        excludedInteractions_.removePair(c, d);
611 >      }
612  
613 <      exclude_.removePair(a, b);
614 <      exclude_.removePair(a, c);
615 <      exclude_.removePair(a, d);
616 <      exclude_.removePair(b, c);
617 <      exclude_.removePair(b, d);
618 <      exclude_.removePair(c, d);        
613 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
614 >        oneThreeInteractions_.removePair(a, c);      
615 >        oneThreeInteractions_.removePair(b, d);      
616 >      } else {
617 >        excludedInteractions_.removePair(a, c);
618 >        excludedInteractions_.removePair(b, d);
619 >      }
620 >
621 >      if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
622 >        oneFourInteractions_.removePair(a, d);      
623 >      } else {
624 >        excludedInteractions_.removePair(a, d);
625 >      }
626      }
627  
628 <    Molecule::RigidBodyIterator rbIter;
629 <    RigidBody* rb;
630 <    for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
628 >    for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
629 >         inversion = mol->nextInversion(inversionIter)) {
630 >
631 >      a = inversion->getAtomA()->getGlobalIndex();
632 >      b = inversion->getAtomB()->getGlobalIndex();        
633 >      c = inversion->getAtomC()->getGlobalIndex();        
634 >      d = inversion->getAtomD()->getGlobalIndex();        
635 >
636 >      if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
637 >        oneTwoInteractions_.removePair(a, b);      
638 >        oneTwoInteractions_.removePair(a, c);
639 >        oneTwoInteractions_.removePair(a, d);
640 >      } else {
641 >        excludedInteractions_.removePair(a, b);
642 >        excludedInteractions_.removePair(a, c);
643 >        excludedInteractions_.removePair(a, d);
644 >      }
645 >
646 >      if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
647 >        oneThreeInteractions_.removePair(b, c);    
648 >        oneThreeInteractions_.removePair(b, d);    
649 >        oneThreeInteractions_.removePair(c, d);      
650 >      } else {
651 >        excludedInteractions_.removePair(b, c);
652 >        excludedInteractions_.removePair(b, d);
653 >        excludedInteractions_.removePair(c, d);
654 >      }
655 >    }
656 >
657 >    for (rb = mol->beginRigidBody(rbIter); rb != NULL;
658 >         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) {
660 >      for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
661 >        for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
662            a = atoms[i]->getGlobalIndex();
663            b = atoms[j]->getGlobalIndex();
664 <          exclude_.removePair(a, b);
664 >          excludedInteractions_.removePair(a, b);
665          }
666        }
667      }        
668 <
668 >    
669    }
670 <
671 <
670 >  
671 >  
672    void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
673      int curStampId;
674 <
674 >    
675      //index from 0
676      curStampId = moleculeStamps_.size();
677  
# Line 463 | Line 692 | namespace oopse {
692      //setup fortran force field
693      /** @deprecate */    
694      int isError = 0;
695 <    initFortranFF( &fInfo_.SIM_uses_RF, &fInfo_.SIM_uses_UW,
696 <                   &fInfo_.SIM_uses_DW, &isError );
695 >    
696 >    setupCutoff();
697 >    
698 >    setupElectrostaticSummationMethod( isError );
699 >    setupSwitchingFunction();
700 >    setupAccumulateBoxDipole();
701 >
702      if(isError){
703        sprintf( painCave.errMsg,
704                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
705        painCave.isFatal = 1;
706        simError();
707      }
474 –  
475 –    
476 –    setupCutoff();
708  
709      calcNdf();
710      calcNdfRaw();
# Line 508 | Line 739 | namespace oopse {
739      int useLennardJones = 0;
740      int useElectrostatic = 0;
741      int useEAM = 0;
742 +    int useSC = 0;
743      int useCharge = 0;
744      int useDirectional = 0;
745      int useDipole = 0;
# Line 519 | Line 751 | namespace oopse {
751      int useDirectionalAtom = 0;    
752      int useElectrostatics = 0;
753      //usePBC and useRF are from simParams
754 <    int usePBC = simParams_->getPBC();
755 <    int useRF = simParams_->getUseRF();
756 <    int useUW = simParams_->getUseUndampedWolf();
757 <    int useDW = simParams_->getUseDampedWolf();
754 >    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
755 >    int useRF;
756 >    int useSF;
757 >    int useSP;
758 >    int useBoxDipole;
759 >
760 >    std::string myMethod;
761 >
762 >    // set the useRF logical
763 >    useRF = 0;
764 >    useSF = 0;
765 >    useSP = 0;
766 >    useBoxDipole = 0;
767 >
768 >
769 >    if (simParams_->haveElectrostaticSummationMethod()) {
770 >      std::string myMethod = simParams_->getElectrostaticSummationMethod();
771 >      toUpper(myMethod);
772 >      if (myMethod == "REACTION_FIELD"){
773 >        useRF = 1;
774 >      } else if (myMethod == "SHIFTED_FORCE"){
775 >        useSF = 1;
776 >      } else if (myMethod == "SHIFTED_POTENTIAL"){
777 >        useSP = 1;
778 >      }
779 >    }
780 >    
781 >    if (simParams_->haveAccumulateBoxDipole())
782 >      if (simParams_->getAccumulateBoxDipole())
783 >        useBoxDipole = 1;
784 >
785 >    useAtomicVirial_ = simParams_->getUseAtomicVirial();
786  
787      //loop over all of the atom types
788      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
789        useLennardJones |= (*i)->isLennardJones();
790        useElectrostatic |= (*i)->isElectrostatic();
791        useEAM |= (*i)->isEAM();
792 +      useSC |= (*i)->isSC();
793        useCharge |= (*i)->isCharge();
794        useDirectional |= (*i)->isDirectional();
795        useDipole |= (*i)->isDipole();
# Line 579 | Line 840 | namespace oopse {
840      temp = useEAM;
841      MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
842  
843 +    temp = useSC;
844 +    MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
845 +    
846      temp = useShape;
847      MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
848  
# Line 588 | Line 852 | namespace oopse {
852      temp = useRF;
853      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
854  
855 <    temp = useUW;
856 <    MPI_Allreduce(&temp, &useUW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
855 >    temp = useSF;
856 >    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
857  
858 <    temp = useDW;
859 <    MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
860 <    
858 >    temp = useSP;
859 >    MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
860 >
861 >    temp = useBoxDipole;
862 >    MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
863 >
864 >    temp = useAtomicVirial_;
865 >    MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
866 >
867   #endif
868  
869      fInfo_.SIM_uses_PBC = usePBC;    
# Line 606 | Line 876 | namespace oopse {
876      fInfo_.SIM_uses_StickyPower = useStickyPower;
877      fInfo_.SIM_uses_GayBerne = useGayBerne;
878      fInfo_.SIM_uses_EAM = useEAM;
879 +    fInfo_.SIM_uses_SC = useSC;
880      fInfo_.SIM_uses_Shapes = useShape;
881      fInfo_.SIM_uses_FLARB = useFLARB;
882      fInfo_.SIM_uses_RF = useRF;
883 <    fInfo_.SIM_uses_UW = useUW;
884 <    fInfo_.SIM_uses_DW = useDW;
885 <
886 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
616 <
617 <      if (simParams_->haveDielectric()) {
618 <        fInfo_.dielect = simParams_->getDielectric();
619 <      } else {
620 <        sprintf(painCave.errMsg,
621 <                "SimSetup Error: No Dielectric constant was set.\n"
622 <                "\tYou are trying to use Reaction Field without"
623 <                "\tsetting a dielectric constant!\n");
624 <        painCave.isFatal = 1;
625 <        simError();
626 <      }
627 <        
628 <    } else {
629 <      fInfo_.dielect = 0.0;
630 <    }
631 <
883 >    fInfo_.SIM_uses_SF = useSF;
884 >    fInfo_.SIM_uses_SP = useSP;
885 >    fInfo_.SIM_uses_BoxDipole = useBoxDipole;
886 >    fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_;
887    }
888  
889    void SimInfo::setupFortranSim() {
890      int isError;
891 <    int nExclude;
891 >    int nExclude, nOneTwo, nOneThree, nOneFour;
892      std::vector<int> fortranGlobalGroupMembership;
893      
639 –    nExclude = exclude_.getSize();
894      isError = 0;
895  
896      //globalGroupMembership_ is filled by SimCreator    
# Line 645 | Line 899 | namespace oopse {
899      }
900  
901      //calculate mass ratio of cutoff group
902 <    std::vector<double> mfact;
902 >    std::vector<RealType> mfact;
903      SimInfo::MoleculeIterator mi;
904      Molecule* mol;
905      Molecule::CutoffGroupIterator ci;
906      CutoffGroup* cg;
907      Molecule::AtomIterator ai;
908      Atom* atom;
909 <    double totalMass;
909 >    RealType totalMass;
910  
911      //to avoid memory reallocation, reserve enough space for mfact
912      mfact.reserve(getNCutoffGroups());
# Line 662 | Line 916 | namespace oopse {
916  
917          totalMass = cg->getMass();
918          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
919 <          mfact.push_back(atom->getMass()/totalMass);
919 >          // Check for massless groups - set mfact to 1 if true
920 >          if (totalMass != 0)
921 >            mfact.push_back(atom->getMass()/totalMass);
922 >          else
923 >            mfact.push_back( 1.0 );
924          }
667 –
925        }      
926      }
927  
# Line 688 | Line 945 | namespace oopse {
945      }
946      
947      //setup fortran simulation
948 <    int nGlobalExcludes = 0;
949 <    int* globalExcludes = NULL;
950 <    int* excludeList = exclude_.getExcludeList();
951 <    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
952 <                   &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
696 <                   &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
948 >
949 >    nExclude = excludedInteractions_.getSize();
950 >    nOneTwo = oneTwoInteractions_.getSize();
951 >    nOneThree = oneThreeInteractions_.getSize();
952 >    nOneFour = oneFourInteractions_.getSize();
953  
954 <    if( isError ){
954 >    int* excludeList = excludedInteractions_.getPairList();
955 >    int* oneTwoList = oneTwoInteractions_.getPairList();
956 >    int* oneThreeList = oneThreeInteractions_.getPairList();
957 >    int* oneFourList = oneFourInteractions_.getPairList();
958  
959 +    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0],
960 +                   &nExclude, excludeList,
961 +                   &nOneTwo, oneTwoList,
962 +                   &nOneThree, oneThreeList,
963 +                   &nOneFour, oneFourList,
964 +                   &molMembershipArray[0], &mfact[0], &nCutoffGroups_,
965 +                   &fortranGlobalGroupMembership[0], &isError);
966 +    
967 +    if( isError ){
968 +      
969        sprintf( painCave.errMsg,
970                 "There was an error setting the simulation information in fortran.\n" );
971        painCave.isFatal = 1;
972        painCave.severity = OOPSE_ERROR;
973        simError();
974      }
975 <
976 < #ifdef IS_MPI
975 >    
976 >    
977      sprintf( checkPointMsg,
978               "succesfully sent the simulation information to fortran.\n");
979 <    MPIcheckPoint();
980 < #endif // is_mpi
979 >    
980 >    errorCheckPoint();
981 >    
982 >    // Setup number of neighbors in neighbor list if present
983 >    if (simParams_->haveNeighborListNeighbors()) {
984 >      int nlistNeighbors = simParams_->getNeighborListNeighbors();
985 >      setNeighbors(&nlistNeighbors);
986 >    }
987 >  
988 >
989    }
990  
991  
715 – #ifdef IS_MPI
992    void SimInfo::setupFortranParallel() {
993 <    
993 > #ifdef IS_MPI    
994      //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
995      std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
996      std::vector<int> localToGlobalCutoffGroupIndex;
# Line 764 | Line 1040 | namespace oopse {
1040      }
1041  
1042      sprintf(checkPointMsg, " mpiRefresh successful.\n");
1043 <    MPIcheckPoint();
1043 >    errorCheckPoint();
1044  
1045 <
1045 > #endif
1046    }
1047  
1048 < #endif
1048 >  void SimInfo::setupCutoff() {          
1049 >    
1050 >    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
1051  
1052 <  double SimInfo::calcMaxCutoffRadius() {
1052 >    // Check the cutoff policy
1053 >    int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default
1054  
1055 +    // Set LJ shifting bools to false
1056 +    ljsp_ = 0;
1057 +    ljsf_ = 0;
1058  
1059 <    std::set<AtomType*> atomTypes;
1060 <    std::set<AtomType*>::iterator i;
1061 <    std::vector<double> cutoffRadius;
1062 <
1063 <    //get the unique atom types
782 <    atomTypes = getUniqueAtomTypes();
783 <
784 <    //query the max cutoff radius among these atom types
785 <    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
786 <      cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
1059 >    std::string myPolicy;
1060 >    if (forceFieldOptions_.haveCutoffPolicy()){
1061 >      myPolicy = forceFieldOptions_.getCutoffPolicy();
1062 >    }else if (simParams_->haveCutoffPolicy()) {
1063 >      myPolicy = simParams_->getCutoffPolicy();
1064      }
1065  
1066 <    double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
1067 < #ifdef IS_MPI
1068 <    //pick the max cutoff radius among the processors
1069 < #endif
1066 >    if (!myPolicy.empty()){
1067 >      toUpper(myPolicy);
1068 >      if (myPolicy == "MIX") {
1069 >        cp = MIX_CUTOFF_POLICY;
1070 >      } else {
1071 >        if (myPolicy == "MAX") {
1072 >          cp = MAX_CUTOFF_POLICY;
1073 >        } else {
1074 >          if (myPolicy == "TRADITIONAL") {            
1075 >            cp = TRADITIONAL_CUTOFF_POLICY;
1076 >          } else {
1077 >            // throw error        
1078 >            sprintf( painCave.errMsg,
1079 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
1080 >            painCave.isFatal = 1;
1081 >            simError();
1082 >          }    
1083 >        }          
1084 >      }
1085 >    }          
1086 >    notifyFortranCutoffPolicy(&cp);
1087  
1088 <    return maxCutoffRadius;
1089 <  }
1088 >    // Check the Skin Thickness for neighborlists
1089 >    RealType skin;
1090 >    if (simParams_->haveSkinThickness()) {
1091 >      skin = simParams_->getSkinThickness();
1092 >      notifyFortranSkinThickness(&skin);
1093 >    }            
1094 >        
1095 >    // Check if the cutoff was set explicitly:
1096 >    if (simParams_->haveCutoffRadius()) {
1097 >      rcut_ = simParams_->getCutoffRadius();
1098 >      if (simParams_->haveSwitchingRadius()) {
1099 >        rsw_  = simParams_->getSwitchingRadius();
1100 >      } else {
1101 >        if (fInfo_.SIM_uses_Charges |
1102 >            fInfo_.SIM_uses_Dipoles |
1103 >            fInfo_.SIM_uses_RF) {
1104 >          
1105 >          rsw_ = 0.85 * rcut_;
1106 >          sprintf(painCave.errMsg,
1107 >                  "SimCreator Warning: No value was set for the switchingRadius.\n"
1108 >                  "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n"
1109 >                  "\tswitchingRadius = %f. for this simulation\n", rsw_);
1110 >        painCave.isFatal = 0;
1111 >        simError();
1112 >        } else {
1113 >          rsw_ = rcut_;
1114 >          sprintf(painCave.errMsg,
1115 >                  "SimCreator Warning: No value was set for the switchingRadius.\n"
1116 >                  "\tOOPSE will use the same value as the cutoffRadius.\n"
1117 >                  "\tswitchingRadius = %f. for this simulation\n", rsw_);
1118 >          painCave.isFatal = 0;
1119 >          simError();
1120 >        }
1121 >      }
1122  
1123 <  void SimInfo::getCutoff(double& rcut, double& rsw) {
1124 <    
1125 <    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
1123 >      if (simParams_->haveElectrostaticSummationMethod()) {
1124 >        std::string myMethod = simParams_->getElectrostaticSummationMethod();
1125 >        toUpper(myMethod);
1126          
1127 <      if (!simParams_->haveRcut()){
1128 <        sprintf(painCave.errMsg,
1127 >        if (myMethod == "SHIFTED_POTENTIAL") {
1128 >          ljsp_ = 1;
1129 >        } else if (myMethod == "SHIFTED_FORCE") {
1130 >          ljsf_ = 1;
1131 >        }
1132 >      }
1133 >
1134 >      notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
1135 >      
1136 >    } else {
1137 >      
1138 >      // For electrostatic atoms, we'll assume a large safe value:
1139 >      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
1140 >        sprintf(painCave.errMsg,
1141                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
1142                  "\tOOPSE will use a default value of 15.0 angstroms"
1143                  "\tfor the cutoffRadius.\n");
1144 <        painCave.isFatal = 0;
1144 >        painCave.isFatal = 0;
1145          simError();
1146 <        rcut = 15.0;
1147 <      } else{
1148 <        rcut = simParams_->getRcut();
1149 <      }
1146 >        rcut_ = 15.0;
1147 >      
1148 >        if (simParams_->haveElectrostaticSummationMethod()) {
1149 >          std::string myMethod = simParams_->getElectrostaticSummationMethod();
1150 >          toUpper(myMethod);
1151 >      
1152 >      // For the time being, we're tethering the LJ shifted behavior to the
1153 >      // electrostaticSummationMethod keyword options
1154 >          if (myMethod == "SHIFTED_POTENTIAL") {
1155 >            ljsp_ = 1;
1156 >          } else if (myMethod == "SHIFTED_FORCE") {
1157 >            ljsf_ = 1;
1158 >          }
1159 >          if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") {
1160 >            if (simParams_->haveSwitchingRadius()){
1161 >              sprintf(painCave.errMsg,
1162 >                      "SimInfo Warning: A value was set for the switchingRadius\n"
1163 >                      "\teven though the electrostaticSummationMethod was\n"
1164 >                      "\tset to %s\n", myMethod.c_str());
1165 >              painCave.isFatal = 1;
1166 >              simError();            
1167 >            }
1168 >          }
1169 >        }
1170 >      
1171 >        if (simParams_->haveSwitchingRadius()){
1172 >          rsw_ = simParams_->getSwitchingRadius();
1173 >        } else {        
1174 >          sprintf(painCave.errMsg,
1175 >                  "SimCreator Warning: No value was set for switchingRadius.\n"
1176 >                  "\tOOPSE will use a default value of\n"
1177 >                  "\t0.85 * cutoffRadius for the switchingRadius\n");
1178 >          painCave.isFatal = 0;
1179 >          simError();
1180 >          rsw_ = 0.85 * rcut_;
1181 >        }
1182  
1183 <      if (!simParams_->haveRsw()){
814 <        sprintf(painCave.errMsg,
815 <                "SimCreator Warning: No value was set for switchingRadius.\n"
816 <                "\tOOPSE will use a default value of\n"
817 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
818 <        painCave.isFatal = 0;
819 <        simError();
820 <        rsw = 0.95 * rcut;
821 <      } else{
822 <        rsw = simParams_->getRsw();
823 <      }
1183 >        notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
1184  
825 –    } else {
826 –      // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
827 –      //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
828 –        
829 –      if (simParams_->haveRcut()) {
830 –        rcut = simParams_->getRcut();
1185        } else {
1186 <        //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
1187 <        rcut = calcMaxCutoffRadius();
1186 >        // We didn't set rcut explicitly, and we don't have electrostatic atoms, so
1187 >        // We'll punt and let fortran figure out the cutoffs later.
1188 >        
1189 >        notifyFortranYouAreOnYourOwn();
1190 >
1191        }
1192 +    }
1193 +  }
1194  
1195 <      if (simParams_->haveRsw()) {
1196 <        rsw  = simParams_->getRsw();
1195 >  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
1196 >    
1197 >    int errorOut;
1198 >    int esm =  NONE;
1199 >    int sm = UNDAMPED;
1200 >    RealType alphaVal;
1201 >    RealType dielectric;
1202 >    
1203 >    errorOut = isError;
1204 >
1205 >    if (simParams_->haveElectrostaticSummationMethod()) {
1206 >      std::string myMethod = simParams_->getElectrostaticSummationMethod();
1207 >      toUpper(myMethod);
1208 >      if (myMethod == "NONE") {
1209 >        esm = NONE;
1210        } else {
1211 <        rsw = rcut;
1211 >        if (myMethod == "SWITCHING_FUNCTION") {
1212 >          esm = SWITCHING_FUNCTION;
1213 >        } else {
1214 >          if (myMethod == "SHIFTED_POTENTIAL") {
1215 >            esm = SHIFTED_POTENTIAL;
1216 >          } else {
1217 >            if (myMethod == "SHIFTED_FORCE") {            
1218 >              esm = SHIFTED_FORCE;
1219 >            } else {
1220 >              if (myMethod == "REACTION_FIELD") {
1221 >                esm = REACTION_FIELD;
1222 >                dielectric = simParams_->getDielectric();
1223 >                if (!simParams_->haveDielectric()) {
1224 >                  // throw warning
1225 >                  sprintf( painCave.errMsg,
1226 >                           "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
1227 >                           "\tA default value of %f will be used for the dielectric.\n", dielectric);
1228 >                  painCave.isFatal = 0;
1229 >                  simError();
1230 >                }
1231 >              } else {
1232 >                // throw error        
1233 >                sprintf( painCave.errMsg,
1234 >                         "SimInfo error: Unknown electrostaticSummationMethod.\n"
1235 >                         "\t(Input file specified %s .)\n"
1236 >                         "\telectrostaticSummationMethod must be one of: \"none\",\n"
1237 >                         "\t\"shifted_potential\", \"shifted_force\", or \n"
1238 >                         "\t\"reaction_field\".\n", myMethod.c_str() );
1239 >                painCave.isFatal = 1;
1240 >                simError();
1241 >              }    
1242 >            }          
1243 >          }
1244 >        }
1245        }
1246 +    }
1247      
1248 +    if (simParams_->haveElectrostaticScreeningMethod()) {
1249 +      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
1250 +      toUpper(myScreen);
1251 +      if (myScreen == "UNDAMPED") {
1252 +        sm = UNDAMPED;
1253 +      } else {
1254 +        if (myScreen == "DAMPED") {
1255 +          sm = DAMPED;
1256 +          if (!simParams_->haveDampingAlpha()) {
1257 +            // first set a cutoff dependent alpha value
1258 +            // we assume alpha depends linearly with rcut from 0 to 20.5 ang
1259 +            alphaVal = 0.5125 - rcut_* 0.025;
1260 +            // for values rcut > 20.5, alpha is zero
1261 +            if (alphaVal < 0) alphaVal = 0;
1262 +
1263 +            // throw warning
1264 +            sprintf( painCave.errMsg,
1265 +                     "SimInfo warning: dampingAlpha was not specified in the input file.\n"
1266 +                     "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_);
1267 +            painCave.isFatal = 0;
1268 +            simError();
1269 +          } else {
1270 +            alphaVal = simParams_->getDampingAlpha();
1271 +          }
1272 +          
1273 +        } else {
1274 +          // throw error        
1275 +          sprintf( painCave.errMsg,
1276 +                   "SimInfo error: Unknown electrostaticScreeningMethod.\n"
1277 +                   "\t(Input file specified %s .)\n"
1278 +                   "\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
1279 +                   "or \"damped\".\n", myScreen.c_str() );
1280 +          painCave.isFatal = 1;
1281 +          simError();
1282 +        }
1283 +      }
1284      }
1285 +    
1286 +    // let's pass some summation method variables to fortran
1287 +    setElectrostaticSummationMethod( &esm );
1288 +    setFortranElectrostaticMethod( &esm );
1289 +    setScreeningMethod( &sm );
1290 +    setDampingAlpha( &alphaVal );
1291 +    setReactionFieldDielectric( &dielectric );
1292 +    initFortranFF( &errorOut );
1293    }
1294  
1295 <  void SimInfo::setupCutoff() {    
1296 <    getCutoff(rcut_, rsw_);    
847 <    double rnblist = rcut_ + 1; // skin of neighbor list
1295 >  void SimInfo::setupSwitchingFunction() {    
1296 >    int ft = CUBIC;
1297  
1298 <    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
1299 <    
1300 <    int cp =  TRADITIONAL_CUTOFF_POLICY;
1301 <    if (simParams_->haveCutoffPolicy()) {
1302 <      std::string myPolicy = simParams_->getCutoffPolicy();
854 <      if (myPolicy == "MIX") {
855 <        cp = MIX_CUTOFF_POLICY;
1298 >    if (simParams_->haveSwitchingFunctionType()) {
1299 >      std::string funcType = simParams_->getSwitchingFunctionType();
1300 >      toUpper(funcType);
1301 >      if (funcType == "CUBIC") {
1302 >        ft = CUBIC;
1303        } else {
1304 <        if (myPolicy == "MAX") {
1305 <          cp = MAX_CUTOFF_POLICY;
1306 <        } else {
1307 <          if (myPolicy == "TRADITIONAL") {            
1308 <            cp = TRADITIONAL_CUTOFF_POLICY;
1309 <          } else {
1310 <            // throw error        
1311 <            sprintf( painCave.errMsg,
865 <                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
866 <            painCave.isFatal = 1;
867 <            simError();
868 <          }    
1304 >        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
1305 >          ft = FIFTH_ORDER_POLY;
1306 >        } else {
1307 >          // throw error        
1308 >          sprintf( painCave.errMsg,
1309 >                   "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
1310 >          painCave.isFatal = 1;
1311 >          simError();
1312          }          
1313        }
1314      }
1315 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
1315 >
1316 >    // send switching function notification to switcheroo
1317 >    setFunctionType(&ft);
1318 >
1319    }
1320  
1321 +  void SimInfo::setupAccumulateBoxDipole() {    
1322 +
1323 +    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
1324 +    if ( simParams_->haveAccumulateBoxDipole() )
1325 +      if ( simParams_->getAccumulateBoxDipole() ) {
1326 +        setAccumulateBoxDipole();
1327 +        calcBoxDipole_ = true;
1328 +      }
1329 +
1330 +  }
1331 +
1332    void SimInfo::addProperty(GenericData* genData) {
1333      properties_.addProperty(genData);  
1334    }
# Line 928 | Line 1385 | namespace oopse {
1385      Molecule* mol;
1386  
1387      Vector3d comVel(0.0);
1388 <    double totalMass = 0.0;
1388 >    RealType totalMass = 0.0;
1389      
1390  
1391      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1392 <      double mass = mol->getMass();
1392 >      RealType mass = mol->getMass();
1393        totalMass += mass;
1394        comVel += mass * mol->getComVel();
1395      }  
1396  
1397   #ifdef IS_MPI
1398 <    double tmpMass = totalMass;
1398 >    RealType tmpMass = totalMass;
1399      Vector3d tmpComVel(comVel);    
1400 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1401 <    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1400 >    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1401 >    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1402   #endif
1403  
1404      comVel /= totalMass;
# Line 954 | Line 1411 | namespace oopse {
1411      Molecule* mol;
1412  
1413      Vector3d com(0.0);
1414 <    double totalMass = 0.0;
1414 >    RealType totalMass = 0.0;
1415      
1416      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1417 <      double mass = mol->getMass();
1417 >      RealType mass = mol->getMass();
1418        totalMass += mass;
1419        com += mass * mol->getCom();
1420      }  
1421  
1422   #ifdef IS_MPI
1423 <    double tmpMass = totalMass;
1423 >    RealType tmpMass = totalMass;
1424      Vector3d tmpCom(com);    
1425 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1426 <    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1425 >    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1426 >    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1427   #endif
1428  
1429      com /= totalMass;
# Line 990 | Line 1447 | namespace oopse {
1447        Molecule* mol;
1448        
1449      
1450 <      double totalMass = 0.0;
1450 >      RealType totalMass = 0.0;
1451      
1452  
1453        for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1454 <         double mass = mol->getMass();
1454 >         RealType mass = mol->getMass();
1455           totalMass += mass;
1456           com += mass * mol->getCom();
1457           comVel += mass * mol->getComVel();          
1458        }  
1459        
1460   #ifdef IS_MPI
1461 <      double tmpMass = totalMass;
1461 >      RealType tmpMass = totalMass;
1462        Vector3d tmpCom(com);  
1463        Vector3d tmpComVel(comVel);
1464 <      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1465 <      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1466 <      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1464 >      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1465 >      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1466 >      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1467   #endif
1468        
1469        com /= totalMass;
# Line 1025 | Line 1482 | namespace oopse {
1482     void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1483        
1484  
1485 <      double xx = 0.0;
1486 <      double yy = 0.0;
1487 <      double zz = 0.0;
1488 <      double xy = 0.0;
1489 <      double xz = 0.0;
1490 <      double yz = 0.0;
1485 >      RealType xx = 0.0;
1486 >      RealType yy = 0.0;
1487 >      RealType zz = 0.0;
1488 >      RealType xy = 0.0;
1489 >      RealType xz = 0.0;
1490 >      RealType yz = 0.0;
1491        Vector3d com(0.0);
1492        Vector3d comVel(0.0);
1493        
# Line 1042 | Line 1499 | namespace oopse {
1499        Vector3d thisq(0.0);
1500        Vector3d thisv(0.0);
1501  
1502 <      double thisMass = 0.0;
1502 >      RealType thisMass = 0.0;
1503      
1504        
1505        
# Line 1080 | Line 1537 | namespace oopse {
1537   #ifdef IS_MPI
1538        Mat3x3d tmpI(inertiaTensor);
1539        Vector3d tmpAngMom;
1540 <      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1541 <      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1540 >      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1541 >      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1542   #endif
1543                
1544        return;
# Line 1102 | Line 1559 | namespace oopse {
1559        Vector3d thisr(0.0);
1560        Vector3d thisp(0.0);
1561        
1562 <      double thisMass;
1562 >      RealType thisMass;
1563        
1564        for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1565          thisMass = mol->getMass();
# Line 1115 | Line 1572 | namespace oopse {
1572        
1573   #ifdef IS_MPI
1574        Vector3d tmpAngMom;
1575 <      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1575 >      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1576   #endif
1577        
1578        return angularMomentum;
1579     }
1580    
1581 <  
1581 >  StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
1582 >    return IOIndexToIntegrableObject.at(index);
1583 >  }
1584 >  
1585 >  void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) {
1586 >    IOIndexToIntegrableObject= v;
1587 >  }
1588 >
1589 >  /* Returns the Volume of the simulation based on a ellipsoid with semi-axes
1590 >     based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3
1591 >     where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to
1592 >     V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
1593 >  */
1594 >  void SimInfo::getGyrationalVolume(RealType &volume){
1595 >    Mat3x3d intTensor;
1596 >    RealType det;
1597 >    Vector3d dummyAngMom;
1598 >    RealType sysconstants;
1599 >    RealType geomCnst;
1600 >
1601 >    geomCnst = 3.0/2.0;
1602 >    /* Get the inertial tensor and angular momentum for free*/
1603 >    getInertiaTensor(intTensor,dummyAngMom);
1604 >    
1605 >    det = intTensor.determinant();
1606 >    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1607 >    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det);
1608 >    return;
1609 >  }
1610 >
1611 >  void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){
1612 >    Mat3x3d intTensor;
1613 >    Vector3d dummyAngMom;
1614 >    RealType sysconstants;
1615 >    RealType geomCnst;
1616 >
1617 >    geomCnst = 3.0/2.0;
1618 >    /* Get the inertial tensor and angular momentum for free*/
1619 >    getInertiaTensor(intTensor,dummyAngMom);
1620 >    
1621 >    detI = intTensor.determinant();
1622 >    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1623 >    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI);
1624 >    return;
1625 >  }
1626 > /*
1627 >   void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) {
1628 >      assert( v.size() == nAtoms_ + nRigidBodies_);
1629 >      sdByGlobalIndex_ = v;
1630 >    }
1631 >
1632 >    StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) {
1633 >      //assert(index < nAtoms_ + nRigidBodies_);
1634 >      return sdByGlobalIndex_.at(index);
1635 >    }  
1636 > */  
1637   }//end namespace oopse
1638  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines