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 610 by chrisfen, Sun Sep 18 20:45:38 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 83 | 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 108 | Line 121 | namespace oopse {
121          }
122  
123          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
124 +
125          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
126  
127          //calculate atoms in rigid bodies
# Line 124 | 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 338 | 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 349 | 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 360 | 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  
372 –    Molecule::RigidBodyIterator rbIter;
373 –    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 395 | 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 407 | 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 418 | 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  
429 –    Molecule::RigidBodyIterator rbIter;
430 –    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 467 | Line 595 | namespace oopse {
595      int isError = 0;
596      
597      setupElectrostaticSummationMethod( isError );
598 +    setupSwitchingFunction();
599  
600      if(isError){
601        sprintf( painCave.errMsg,
# Line 511 | 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 522 | Line 652 | namespace oopse {
652      int useDirectionalAtom = 0;    
653      int useElectrostatics = 0;
654      //usePBC and useRF are from simParams
655 <    int usePBC = simParams_->getPBC();
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 579 | 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  
739      temp = useFLARB;
740      MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
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  
750      fInfo_.SIM_uses_PBC = usePBC;    
# Line 597 | 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 839 | Line 1003 | namespace oopse {
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 {
# Line 857 | Line 1022 | namespace oopse {
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_);
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  
# Line 875 | Line 1047 | namespace oopse {
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 == "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,
889 <                       "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);
890 <              painCave.isFatal = 0;
891 <              simError();
892 <            }
893 <          } else {
894 <            if (myMethod == "REACTION_FIELD") {
895 <              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