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 720 by chrisfen, Tue Nov 8 13:32:06 2005 UTC vs.
Revision 749 by tim, Wed Nov 16 23:10:02 2005 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"
# Line 55 | Line 56
56   #include "UseTheForce/fCutoffPolicy.h"
57   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
58   #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
59 + #include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
60   #include "UseTheForce/doForces_interface.h"
61   #include "UseTheForce/DarkSide/electrostatic_interface.h"
62   #include "UseTheForce/notifyCutoffs_interface.h"
63 + #include "UseTheForce/DarkSide/switcheroo_interface.h"
64   #include "utils/MemoryUtils.hpp"
65   #include "utils/simError.h"
66   #include "selection/SelectionManager.hpp"
# Line 68 | Line 71 | namespace oopse {
71   #endif
72  
73   namespace oopse {
74 +  std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) {
75 +    std::map<int, std::set<int> >::iterator i = container.find(index);
76 +    std::set<int> result;
77 +    if (i != container.end()) {
78 +        result = i->second;
79 +    }
80  
81 +    return result;
82 +  }
83 +  
84    SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
85                     ForceField* ff, Globals* simParams) :
86      stamps_(stamps), forceField_(ff), simParams_(simParams),
# Line 343 | Line 355 | namespace oopse {
355      int b;
356      int c;
357      int d;
358 +
359 +    std::map<int, std::set<int> > atomGroups;
360 +
361 +    Molecule::RigidBodyIterator rbIter;
362 +    RigidBody* rb;
363 +    Molecule::IntegrableObjectIterator ii;
364 +    StuntDouble* integrableObject;
365 +    
366 +    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
367 +           integrableObject = mol->nextIntegrableObject(ii)) {
368 +
369 +      if (integrableObject->isRigidBody()) {
370 +          rb = static_cast<RigidBody*>(integrableObject);
371 +          std::vector<Atom*> atoms = rb->getAtoms();
372 +          std::set<int> rigidAtoms;
373 +          for (int i = 0; i < atoms.size(); ++i) {
374 +            rigidAtoms.insert(atoms[i]->getGlobalIndex());
375 +          }
376 +          for (int i = 0; i < atoms.size(); ++i) {
377 +            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
378 +          }      
379 +      } else {
380 +        std::set<int> oneAtomSet;
381 +        oneAtomSet.insert(integrableObject->getGlobalIndex());
382 +        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
383 +      }
384 +    }  
385 +
386      
387 +    
388      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
389        a = bond->getAtomA()->getGlobalIndex();
390        b = bond->getAtomB()->getGlobalIndex();        
# Line 354 | Line 395 | namespace oopse {
395        a = bend->getAtomA()->getGlobalIndex();
396        b = bend->getAtomB()->getGlobalIndex();        
397        c = bend->getAtomC()->getGlobalIndex();
398 +      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
399 +      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
400 +      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
401  
402 <      exclude_.addPair(a, b);
403 <      exclude_.addPair(a, c);
404 <      exclude_.addPair(b, c);        
402 >      exclude_.addPairs(rigidSetA, rigidSetB);
403 >      exclude_.addPairs(rigidSetA, rigidSetC);
404 >      exclude_.addPairs(rigidSetB, rigidSetC);
405 >      
406 >      //exclude_.addPair(a, b);
407 >      //exclude_.addPair(a, c);
408 >      //exclude_.addPair(b, c);        
409      }
410  
411      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
# Line 365 | Line 413 | namespace oopse {
413        b = torsion->getAtomB()->getGlobalIndex();        
414        c = torsion->getAtomC()->getGlobalIndex();        
415        d = torsion->getAtomD()->getGlobalIndex();        
416 +      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
417 +      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
418 +      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
419 +      std::set<int> rigidSetD = getRigidSet(d, atomGroups);
420  
421 +      exclude_.addPairs(rigidSetA, rigidSetB);
422 +      exclude_.addPairs(rigidSetA, rigidSetC);
423 +      exclude_.addPairs(rigidSetA, rigidSetD);
424 +      exclude_.addPairs(rigidSetB, rigidSetC);
425 +      exclude_.addPairs(rigidSetB, rigidSetD);
426 +      exclude_.addPairs(rigidSetC, rigidSetD);
427 +
428 +      /*
429 +      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
430 +      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
431 +      exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
432 +      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
433 +      exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
434 +      exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
435 +        
436 +      
437        exclude_.addPair(a, b);
438        exclude_.addPair(a, c);
439        exclude_.addPair(a, d);
440        exclude_.addPair(b, c);
441        exclude_.addPair(b, d);
442        exclude_.addPair(c, d);        
443 +      */
444      }
445  
377 –    Molecule::RigidBodyIterator rbIter;
378 –    RigidBody* rb;
446      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
447        std::vector<Atom*> atoms = rb->getAtoms();
448        for (int i = 0; i < atoms.size() -1 ; ++i) {
# Line 400 | Line 467 | namespace oopse {
467      int b;
468      int c;
469      int d;
470 +
471 +    std::map<int, std::set<int> > atomGroups;
472 +
473 +    Molecule::RigidBodyIterator rbIter;
474 +    RigidBody* rb;
475 +    Molecule::IntegrableObjectIterator ii;
476 +    StuntDouble* integrableObject;
477 +    
478 +    for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
479 +           integrableObject = mol->nextIntegrableObject(ii)) {
480 +
481 +      if (integrableObject->isRigidBody()) {
482 +          rb = static_cast<RigidBody*>(integrableObject);
483 +          std::vector<Atom*> atoms = rb->getAtoms();
484 +          std::set<int> rigidAtoms;
485 +          for (int i = 0; i < atoms.size(); ++i) {
486 +            rigidAtoms.insert(atoms[i]->getGlobalIndex());
487 +          }
488 +          for (int i = 0; i < atoms.size(); ++i) {
489 +            atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
490 +          }      
491 +      } else {
492 +        std::set<int> oneAtomSet;
493 +        oneAtomSet.insert(integrableObject->getGlobalIndex());
494 +        atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));        
495 +      }
496 +    }  
497 +
498      
499      for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
500        a = bond->getAtomA()->getGlobalIndex();
# Line 412 | Line 507 | namespace oopse {
507        b = bend->getAtomB()->getGlobalIndex();        
508        c = bend->getAtomC()->getGlobalIndex();
509  
510 <      exclude_.removePair(a, b);
511 <      exclude_.removePair(a, c);
512 <      exclude_.removePair(b, c);        
510 >      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
511 >      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
512 >      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
513 >
514 >      exclude_.removePairs(rigidSetA, rigidSetB);
515 >      exclude_.removePairs(rigidSetA, rigidSetC);
516 >      exclude_.removePairs(rigidSetB, rigidSetC);
517 >      
518 >      //exclude_.removePair(a, b);
519 >      //exclude_.removePair(a, c);
520 >      //exclude_.removePair(b, c);        
521      }
522  
523      for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
# Line 423 | Line 526 | namespace oopse {
526        c = torsion->getAtomC()->getGlobalIndex();        
527        d = torsion->getAtomD()->getGlobalIndex();        
528  
529 +      std::set<int> rigidSetA = getRigidSet(a, atomGroups);
530 +      std::set<int> rigidSetB = getRigidSet(b, atomGroups);
531 +      std::set<int> rigidSetC = getRigidSet(c, atomGroups);
532 +      std::set<int> rigidSetD = getRigidSet(d, atomGroups);
533 +
534 +      exclude_.removePairs(rigidSetA, rigidSetB);
535 +      exclude_.removePairs(rigidSetA, rigidSetC);
536 +      exclude_.removePairs(rigidSetA, rigidSetD);
537 +      exclude_.removePairs(rigidSetB, rigidSetC);
538 +      exclude_.removePairs(rigidSetB, rigidSetD);
539 +      exclude_.removePairs(rigidSetC, rigidSetD);
540 +
541 +      /*
542 +      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end());
543 +      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end());
544 +      exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end());
545 +      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end());
546 +      exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end());
547 +      exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end());
548 +
549 +      
550        exclude_.removePair(a, b);
551        exclude_.removePair(a, c);
552        exclude_.removePair(a, d);
553        exclude_.removePair(b, c);
554        exclude_.removePair(b, d);
555        exclude_.removePair(c, d);        
556 +      */
557      }
558  
434 –    Molecule::RigidBodyIterator rbIter;
435 –    RigidBody* rb;
559      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
560        std::vector<Atom*> atoms = rb->getAtoms();
561        for (int i = 0; i < atoms.size() -1 ; ++i) {
# Line 472 | Line 595 | namespace oopse {
595      int isError = 0;
596      
597      setupElectrostaticSummationMethod( isError );
598 +    setupSwitchingFunction();
599  
600      if(isError){
601        sprintf( painCave.errMsg,
# Line 516 | Line 640 | namespace oopse {
640      int useLennardJones = 0;
641      int useElectrostatic = 0;
642      int useEAM = 0;
643 +    int useSC = 0;
644      int useCharge = 0;
645      int useDirectional = 0;
646      int useDipole = 0;
# Line 554 | Line 679 | namespace oopse {
679        useLennardJones |= (*i)->isLennardJones();
680        useElectrostatic |= (*i)->isElectrostatic();
681        useEAM |= (*i)->isEAM();
682 +      useSC |= (*i)->isSC();
683        useCharge |= (*i)->isCharge();
684        useDirectional |= (*i)->isDirectional();
685        useDipole |= (*i)->isDipole();
# Line 604 | Line 730 | namespace oopse {
730      temp = useEAM;
731      MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
732  
733 +    temp = useSC;
734 +    MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
735 +    
736      temp = useShape;
737      MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
738  
# Line 628 | Line 757 | namespace oopse {
757      fInfo_.SIM_uses_StickyPower = useStickyPower;
758      fInfo_.SIM_uses_GayBerne = useGayBerne;
759      fInfo_.SIM_uses_EAM = useEAM;
760 +    fInfo_.SIM_uses_SC = useSC;
761      fInfo_.SIM_uses_Shapes = useShape;
762      fInfo_.SIM_uses_FLARB = useFLARB;
763      fInfo_.SIM_uses_RF = useRF;
# Line 978 | Line 1108 | namespace oopse {
1108      initFortranFF( &esm, &errorOut );
1109    }
1110  
1111 +  void SimInfo::setupSwitchingFunction() {    
1112 +    int ft = CUBIC;
1113 +
1114 +    if (simParams_->haveSwitchingFunctionType()) {
1115 +      std::string funcType = simParams_->getSwitchingFunctionType();
1116 +      toUpper(funcType);
1117 +      if (funcType == "CUBIC") {
1118 +        ft = CUBIC;
1119 +      } else {
1120 +        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
1121 +          ft = FIFTH_ORDER_POLY;
1122 +        } else {
1123 +          // throw error        
1124 +          sprintf( painCave.errMsg,
1125 +                   "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
1126 +          painCave.isFatal = 1;
1127 +          simError();
1128 +        }          
1129 +      }
1130 +    }
1131 +
1132 +    // send switching function notification to switcheroo
1133 +    setFunctionType(&ft);
1134 +
1135 +  }
1136 +
1137    void SimInfo::addProperty(GenericData* genData) {
1138      properties_.addProperty(genData);  
1139    }

Diff Legend

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