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 557 by chuckv, Tue May 31 22:31:54 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 64 | 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 80 | Line 96 | namespace oopse {
96        MoleculeStamp* molStamp;
97        int nMolWithSameStamp;
98        int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
99 <      int nGroups = 0;          //total cutoff groups defined in meta-data file
99 >      int nGroups = 0;      //total cutoff groups defined in meta-data file
100        CutoffGroupStamp* cgStamp;    
101        RigidBodyStamp* rbStamp;
102        int nRigidAtoms = 0;
# Line 105 | Line 121 | namespace oopse {
121          }
122  
123          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
124 +
125          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
126  
127          //calculate atoms in rigid bodies
# Line 121 | Line 138 | namespace oopse {
138          
139        }
140  
141 <      //every free atom (atom does not belong to cutoff groups) is a cutoff group
142 <      //therefore the total number of cutoff groups in the system is equal to
143 <      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
144 <      //file plus the number of cutoff groups defined in meta-data file
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        nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
147  
148 <      //every free atom (atom does not belong to rigid bodies) is an integrable object
149 <      //therefore the total number of  integrable objects in the system is equal to
150 <      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
151 <      //file plus the number of  rigid bodies defined in meta-data file
152 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
153 <
148 >      //every free atom (atom does not belong to rigid bodies) is an
149 >      //integrable object therefore the total number of integrable objects
150 >      //in the system is equal to the total number of atoms minus number of
151 >      //atoms belong to rigid body defined in meta-data file plus the number
152 >      //of rigid bodies defined in meta-data file
153 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
154 >                                                + nGlobalRigidBodies_;
155 >  
156        nGlobalMols_ = molStampIds_.size();
157  
158   #ifdef IS_MPI    
# Line 335 | 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 346 | 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 357 | 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  
369 –    Molecule::RigidBodyIterator rbIter;
370 –    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 392 | 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 404 | 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 415 | 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  
426 –    Molecule::RigidBodyIterator rbIter;
427 –    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 462 | Line 593 | namespace oopse {
593      //setup fortran force field
594      /** @deprecate */    
595      int isError = 0;
596 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
596 >    
597 >    setupElectrostaticSummationMethod( isError );
598 >    setupSwitchingFunction();
599 >
600      if(isError){
601        sprintf( painCave.errMsg,
602                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 506 | 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 517 | Line 652 | namespace oopse {
652      int useDirectionalAtom = 0;    
653      int useElectrostatics = 0;
654      //usePBC and useRF are from simParams
655 <    int usePBC = simParams_->getPBC();
656 <    int useRF = simParams_->getUseRF();
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()) {
666 >      std::string myMethod = simParams_->getElectrostaticSummationMethod();
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  
677      //loop over all of the atom types
678      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
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 575 | 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 583 | Line 741 | namespace oopse {
741  
742      temp = useRF;
743      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
744 <    
744 >
745 >    temp = useSF;
746 >    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
747 >
748   #endif
749  
750      fInfo_.SIM_uses_PBC = usePBC;    
# Line 596 | 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( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
767 <
766 >    if( myMethod == "REACTION_FIELD") {
767 >      
768        if (simParams_->haveDielectric()) {
769          fInfo_.dielect = simParams_->getDielectric();
770        } else {
# Line 611 | Line 774 | namespace oopse {
774                  "\tsetting a dielectric constant!\n");
775          painCave.isFatal = 1;
776          simError();
777 <      }
615 <        
616 <    } else {
617 <      fInfo_.dielect = 0.0;
777 >      }      
778      }
779  
780    }
# Line 650 | Line 810 | namespace oopse {
810  
811          totalMass = cg->getMass();
812          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
813 <          mfact.push_back(atom->getMass()/totalMass);
813 >          // Check for massless groups - set mfact to 1 if true
814 >          if (totalMass != 0)
815 >            mfact.push_back(atom->getMass()/totalMass);
816 >          else
817 >            mfact.push_back( 1.0 );
818          }
819  
820        }      
# Line 786 | Line 950 | namespace oopse {
950      
951      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
952          
953 <      if (!simParams_->haveRcut()){
953 >      if (!simParams_->haveCutoffRadius()){
954          sprintf(painCave.errMsg,
955                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
956                  "\tOOPSE will use a default value of 15.0 angstroms"
# Line 795 | Line 959 | namespace oopse {
959          simError();
960          rcut = 15.0;
961        } else{
962 <        rcut = simParams_->getRcut();
962 >        rcut = simParams_->getCutoffRadius();
963        }
964  
965 <      if (!simParams_->haveRsw()){
965 >      if (!simParams_->haveSwitchingRadius()){
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_->getRsw();
974 >        rsw = simParams_->getSwitchingRadius();
975        }
976  
977      } else {
978        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
979        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
980          
981 <      if (simParams_->haveRcut()) {
982 <        rcut = simParams_->getRcut();
981 >      if (simParams_->haveCutoffRadius()) {
982 >        rcut = simParams_->getCutoffRadius();
983        } else {
984          //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
985          rcut = calcMaxCutoffRadius();
986        }
987  
988 <      if (simParams_->haveRsw()) {
989 <        rsw  = simParams_->getRsw();
988 >      if (simParams_->haveSwitchingRadius()) {
989 >        rsw  = simParams_->getSwitchingRadius();
990        } else {
991          rsw = rcut;
992        }
# Line 830 | Line 994 | namespace oopse {
994      }
995    }
996  
997 <  void SimInfo::setupCutoff() {
997 >  void SimInfo::setupCutoff() {    
998      getCutoff(rcut_, rsw_);    
999      double rnblist = rcut_ + 1; // skin of neighbor list
1000  
1001      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
1002 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
1002 >    
1003 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
1004 >    if (simParams_->haveCutoffPolicy()) {
1005 >      std::string myPolicy = simParams_->getCutoffPolicy();
1006 >      toUpper(myPolicy);
1007 >      if (myPolicy == "MIX") {
1008 >        cp = MIX_CUTOFF_POLICY;
1009 >      } else {
1010 >        if (myPolicy == "MAX") {
1011 >          cp = MAX_CUTOFF_POLICY;
1012 >        } else {
1013 >          if (myPolicy == "TRADITIONAL") {            
1014 >            cp = TRADITIONAL_CUTOFF_POLICY;
1015 >          } else {
1016 >            // throw error        
1017 >            sprintf( painCave.errMsg,
1018 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
1019 >            painCave.isFatal = 1;
1020 >            simError();
1021 >          }    
1022 >        }          
1023 >      }
1024 >    }
1025 >
1026 >
1027 >    if (simParams_->haveSkinThickness()) {
1028 >      double skinThickness = simParams_->getSkinThickness();
1029 >    }
1030 >
1031 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
1032 >    // also send cutoff notification to electrostatics
1033 >    setElectrostaticCutoffRadius(&rcut_, &rsw_);
1034 >  }
1035 >
1036 >  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
1037 >    
1038 >    int errorOut;
1039 >    int esm =  NONE;
1040 >    int sm = UNDAMPED;
1041 >    double alphaVal;
1042 >    double dielectric;
1043 >
1044 >    errorOut = isError;
1045 >    alphaVal = simParams_->getDampingAlpha();
1046 >    dielectric = simParams_->getDielectric();
1047 >
1048 >    if (simParams_->haveElectrostaticSummationMethod()) {
1049 >      std::string myMethod = simParams_->getElectrostaticSummationMethod();
1050 >      toUpper(myMethod);
1051 >      if (myMethod == "NONE") {
1052 >        esm = NONE;
1053 >      } else {
1054 >        if (myMethod == "SWITCHING_FUNCTION") {
1055 >          esm = SWITCHING_FUNCTION;
1056 >        } else {
1057 >          if (myMethod == "SHIFTED_POTENTIAL") {
1058 >            esm = SHIFTED_POTENTIAL;
1059 >          } else {
1060 >            if (myMethod == "SHIFTED_FORCE") {            
1061 >              esm = SHIFTED_FORCE;
1062 >            } else {
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 >    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) {

Diff Legend

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