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 691 by chrisfen, Wed Oct 19 19:24:40 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"
55   #include "primitives/Molecule.hpp"
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 67 | 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 342 | 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 353 | 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 364 | 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  
376 –    Molecule::RigidBodyIterator rbIter;
377 –    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 399 | 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();
501        b = bond->getAtomB()->getGlobalIndex();        
# Line 411 | 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 422 | 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  
433 –    Molecule::RigidBodyIterator rbIter;
434 –    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 471 | Line 595 | namespace oopse {
595      int isError = 0;
596      
597      setupElectrostaticSummationMethod( isError );
598 +    setupSwitchingFunction();
599  
600      if(isError){
601        sprintf( painCave.errMsg,
# Line 515 | 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 528 | Line 654 | namespace oopse {
654      //usePBC and useRF are from simParams
655      int usePBC = simParams_->getUsePeriodicBoundaryConditions();
656      int useRF;
657 +    int useSF;
658      std::string myMethod;
659  
660      // set the useRF logical
661      useRF = 0;
662 +    useSF = 0;
663  
664  
665      if (simParams_->haveElectrostaticSummationMethod()) {
# Line 539 | Line 667 | namespace oopse {
667        toUpper(myMethod);
668        if (myMethod == "REACTION_FIELD") {
669          useRF=1;
670 +      } else {
671 +        if (myMethod == "SHIFTED_FORCE") {
672 +          useSF = 1;
673 +        }
674        }
675      }
676  
# Line 547 | 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 597 | 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 605 | Line 741 | namespace oopse {
741  
742      temp = useRF;
743      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
744 +
745 +    temp = useSF;
746 +    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
747  
748   #endif
749  
# Line 618 | 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;
764 +    fInfo_.SIM_uses_SF = useSF;
765  
766      if( myMethod == "REACTION_FIELD") {
767        
# Line 635 | Line 776 | namespace oopse {
776          simError();
777        }      
778      }
779 +
780    }
781  
782    void SimInfo::setupFortranSim() {
# Line 824 | Line 966 | namespace oopse {
966          sprintf(painCave.errMsg,
967                  "SimCreator Warning: No value was set for switchingRadius.\n"
968                  "\tOOPSE will use a default value of\n"
969 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
969 >                "\t0.85 * cutoffRadius for the switchingRadius\n");
970          painCave.isFatal = 0;
971          simError();
972 <        rsw = 0.95 * rcut;
972 >        rsw = 0.85 * rcut;
973        } else{
974          rsw = simParams_->getSwitchingRadius();
975        }
# Line 895 | Line 1037 | namespace oopse {
1037      
1038      int errorOut;
1039      int esm =  NONE;
1040 +    int sm = UNDAMPED;
1041      double alphaVal;
1042      double dielectric;
1043  
# Line 908 | Line 1051 | namespace oopse {
1051        if (myMethod == "NONE") {
1052          esm = NONE;
1053        } else {
1054 <        if (myMethod == "UNDAMPED_WOLF") {
1055 <          esm = UNDAMPED_WOLF;
1054 >        if (myMethod == "SWITCHING_FUNCTION") {
1055 >          esm = SWITCHING_FUNCTION;
1056          } else {
1057 <          if (myMethod == "DAMPED_WOLF") {            
1058 <            esm = DAMPED_WOLF;
1059 <            if (!simParams_->haveDampingAlpha()) {
1060 <              //throw error
1061 <              sprintf( painCave.errMsg,
919 <                       "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", alphaVal);
920 <              painCave.isFatal = 0;
921 <              simError();
922 <            }
923 <          } else {
924 <            if (myMethod == "REACTION_FIELD") {      
925 <              esm = REACTION_FIELD;
1057 >          if (myMethod == "SHIFTED_POTENTIAL") {
1058 >            esm = SHIFTED_POTENTIAL;
1059 >          } else {
1060 >            if (myMethod == "SHIFTED_FORCE") {            
1061 >              esm = SHIFTED_FORCE;
1062              } else {
1063 <              // throw error        
1064 <              sprintf( painCave.errMsg,
1065 <                       "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"undamped_wolf\", \"damped_wolf\", or \"reaction_field\".", myMethod.c_str() );
1066 <              painCave.isFatal = 1;
1067 <              simError();
1068 <            }    
1069 <          }          
1063 >              if (myMethod == "REACTION_FIELD") {            
1064 >                esm = REACTION_FIELD;
1065 >              } else {
1066 >                // throw error        
1067 >                sprintf( painCave.errMsg,
1068 >                         "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() );
1069 >                painCave.isFatal = 1;
1070 >                simError();
1071 >              }    
1072 >            }          
1073 >          }
1074          }
1075        }
1076      }
1077 +    
1078 +    if (simParams_->haveElectrostaticScreeningMethod()) {
1079 +      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
1080 +      toUpper(myScreen);
1081 +      if (myScreen == "UNDAMPED") {
1082 +        sm = UNDAMPED;
1083 +      } else {
1084 +        if (myScreen == "DAMPED") {
1085 +          sm = DAMPED;
1086 +          if (!simParams_->haveDampingAlpha()) {
1087 +            //throw error
1088 +            sprintf( painCave.errMsg,
1089 +                     "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal);
1090 +            painCave.isFatal = 0;
1091 +            simError();
1092 +          }
1093 +        } else {
1094 +          // throw error        
1095 +          sprintf( painCave.errMsg,
1096 +                   "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() );
1097 +          painCave.isFatal = 1;
1098 +          simError();
1099 +        }
1100 +      }
1101 +    }
1102 +    
1103      // let's pass some summation method variables to fortran
1104      setElectrostaticSummationMethod( &esm );
1105 <    setDampedWolfAlpha( &alphaVal );
1105 >    setScreeningMethod( &sm );
1106 >    setDampingAlpha( &alphaVal );
1107      setReactionFieldDielectric( &dielectric );
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