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 726 by chrisfen, Fri Nov 11 15:22:11 2005 UTC

# Line 54 | Line 54
54   #include "primitives/Molecule.hpp"
55   #include "UseTheForce/fCutoffPolicy.h"
56   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
57 + #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
58 + #include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
59   #include "UseTheForce/doForces_interface.h"
60   #include "UseTheForce/DarkSide/electrostatic_interface.h"
61   #include "UseTheForce/notifyCutoffs_interface.h"
62 + #include "UseTheForce/DarkSide/switcheroo_interface.h"
63   #include "utils/MemoryUtils.hpp"
64   #include "utils/simError.h"
65   #include "selection/SelectionManager.hpp"
# Line 83 | Line 86 | namespace oopse {
86        MoleculeStamp* molStamp;
87        int nMolWithSameStamp;
88        int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
89 <      int nGroups = 0;          //total cutoff groups defined in meta-data file
89 >      int nGroups = 0;      //total cutoff groups defined in meta-data file
90        CutoffGroupStamp* cgStamp;    
91        RigidBodyStamp* rbStamp;
92        int nRigidAtoms = 0;
# Line 108 | Line 111 | namespace oopse {
111          }
112  
113          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
114 +
115          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
116  
117          //calculate atoms in rigid bodies
# Line 124 | Line 128 | namespace oopse {
128          
129        }
130  
131 <      //every free atom (atom does not belong to cutoff groups) is a cutoff group
132 <      //therefore the total number of cutoff groups in the system is equal to
133 <      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
134 <      //file plus the number of cutoff groups defined in meta-data file
131 >      //every free atom (atom does not belong to cutoff groups) is a cutoff
132 >      //group therefore the total number of cutoff groups in the system is
133 >      //equal to the total number of atoms minus number of atoms belong to
134 >      //cutoff group defined in meta-data file plus the number of cutoff
135 >      //groups defined in meta-data file
136        nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
137  
138 <      //every free atom (atom does not belong to rigid bodies) is an integrable object
139 <      //therefore the total number of  integrable objects in the system is equal to
140 <      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
141 <      //file plus the number of  rigid bodies defined in meta-data file
142 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
143 <
138 >      //every free atom (atom does not belong to rigid bodies) is an
139 >      //integrable object therefore the total number of integrable objects
140 >      //in the system is equal to the total number of atoms minus number of
141 >      //atoms belong to rigid body defined in meta-data file plus the number
142 >      //of rigid bodies defined in meta-data file
143 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
144 >                                                + nGlobalRigidBodies_;
145 >  
146        nGlobalMols_ = molStampIds_.size();
147  
148   #ifdef IS_MPI    
# Line 467 | Line 474 | namespace oopse {
474      int isError = 0;
475      
476      setupElectrostaticSummationMethod( isError );
477 +    setupSwitchingFunction();
478  
479      if(isError){
480        sprintf( painCave.errMsg,
# Line 522 | Line 530 | namespace oopse {
530      int useDirectionalAtom = 0;    
531      int useElectrostatics = 0;
532      //usePBC and useRF are from simParams
533 <    int usePBC = simParams_->getPBC();
533 >    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
534 >    int useRF;
535 >    int useSF;
536 >    std::string myMethod;
537 >
538 >    // set the useRF logical
539 >    useRF = 0;
540 >    useSF = 0;
541 >
542 >
543 >    if (simParams_->haveElectrostaticSummationMethod()) {
544 >      std::string myMethod = simParams_->getElectrostaticSummationMethod();
545 >      toUpper(myMethod);
546 >      if (myMethod == "REACTION_FIELD") {
547 >        useRF=1;
548 >      } else {
549 >        if (myMethod == "SHIFTED_FORCE") {
550 >          useSF = 1;
551 >        }
552 >      }
553 >    }
554  
555      //loop over all of the atom types
556      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
# Line 585 | Line 613 | namespace oopse {
613      temp = useFLARB;
614      MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
615  
616 +    temp = useRF;
617 +    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
618 +
619 +    temp = useSF;
620 +    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
621 +
622   #endif
623  
624      fInfo_.SIM_uses_PBC = usePBC;    
# Line 599 | Line 633 | namespace oopse {
633      fInfo_.SIM_uses_EAM = useEAM;
634      fInfo_.SIM_uses_Shapes = useShape;
635      fInfo_.SIM_uses_FLARB = useFLARB;
636 +    fInfo_.SIM_uses_RF = useRF;
637 +    fInfo_.SIM_uses_SF = useSF;
638  
639 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
640 <
639 >    if( myMethod == "REACTION_FIELD") {
640 >      
641        if (simParams_->haveDielectric()) {
642          fInfo_.dielect = simParams_->getDielectric();
643        } else {
# Line 611 | Line 647 | namespace oopse {
647                  "\tsetting a dielectric constant!\n");
648          painCave.isFatal = 1;
649          simError();
650 <      }
615 <        
616 <    } else {
617 <      fInfo_.dielect = 0.0;
650 >      }      
651      }
652  
653    }
# Line 650 | Line 683 | namespace oopse {
683  
684          totalMass = cg->getMass();
685          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
686 <          mfact.push_back(atom->getMass()/totalMass);
686 >          // Check for massless groups - set mfact to 1 if true
687 >          if (totalMass != 0)
688 >            mfact.push_back(atom->getMass()/totalMass);
689 >          else
690 >            mfact.push_back( 1.0 );
691          }
692  
693        }      
# Line 786 | Line 823 | namespace oopse {
823      
824      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
825          
826 <      if (!simParams_->haveRcut()){
826 >      if (!simParams_->haveCutoffRadius()){
827          sprintf(painCave.errMsg,
828                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
829                  "\tOOPSE will use a default value of 15.0 angstroms"
# Line 795 | Line 832 | namespace oopse {
832          simError();
833          rcut = 15.0;
834        } else{
835 <        rcut = simParams_->getRcut();
835 >        rcut = simParams_->getCutoffRadius();
836        }
837  
838 <      if (!simParams_->haveRsw()){
838 >      if (!simParams_->haveSwitchingRadius()){
839          sprintf(painCave.errMsg,
840                  "SimCreator Warning: No value was set for switchingRadius.\n"
841                  "\tOOPSE will use a default value of\n"
842 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
842 >                "\t0.85 * cutoffRadius for the switchingRadius\n");
843          painCave.isFatal = 0;
844          simError();
845 <        rsw = 0.95 * rcut;
845 >        rsw = 0.85 * rcut;
846        } else{
847 <        rsw = simParams_->getRsw();
847 >        rsw = simParams_->getSwitchingRadius();
848        }
849  
850      } else {
851        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
852        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
853          
854 <      if (simParams_->haveRcut()) {
855 <        rcut = simParams_->getRcut();
854 >      if (simParams_->haveCutoffRadius()) {
855 >        rcut = simParams_->getCutoffRadius();
856        } else {
857          //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
858          rcut = calcMaxCutoffRadius();
859        }
860  
861 <      if (simParams_->haveRsw()) {
862 <        rsw  = simParams_->getRsw();
861 >      if (simParams_->haveSwitchingRadius()) {
862 >        rsw  = simParams_->getSwitchingRadius();
863        } else {
864          rsw = rcut;
865        }
# Line 839 | Line 876 | namespace oopse {
876      int cp =  TRADITIONAL_CUTOFF_POLICY;
877      if (simParams_->haveCutoffPolicy()) {
878        std::string myPolicy = simParams_->getCutoffPolicy();
879 +      toUpper(myPolicy);
880        if (myPolicy == "MIX") {
881          cp = MIX_CUTOFF_POLICY;
882        } else {
# Line 857 | Line 895 | namespace oopse {
895          }          
896        }
897      }
898 +
899 +
900 +    if (simParams_->haveSkinThickness()) {
901 +      double skinThickness = simParams_->getSkinThickness();
902 +    }
903 +
904      notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
905      // also send cutoff notification to electrostatics
906 <    setElectrostaticCutoffRadius(&rcut_);
906 >    setElectrostaticCutoffRadius(&rcut_, &rsw_);
907    }
908  
909    void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
910      
911      int errorOut;
912      int esm =  NONE;
913 +    int sm = UNDAMPED;
914      double alphaVal;
915      double dielectric;
916  
# Line 875 | Line 920 | namespace oopse {
920  
921      if (simParams_->haveElectrostaticSummationMethod()) {
922        std::string myMethod = simParams_->getElectrostaticSummationMethod();
923 +      toUpper(myMethod);
924        if (myMethod == "NONE") {
925          esm = NONE;
926        } else {
927 <        if (myMethod == "UNDAMPED_WOLF") {
928 <          esm = UNDAMPED_WOLF;
927 >        if (myMethod == "SWITCHING_FUNCTION") {
928 >          esm = SWITCHING_FUNCTION;
929          } else {
930 <          if (myMethod == "DAMPED_WOLF") {            
931 <            esm = DAMPED_WOLF;
932 <            if (!simParams_->haveDampingAlpha()) {
933 <              //throw error
934 <              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;
930 >          if (myMethod == "SHIFTED_POTENTIAL") {
931 >            esm = SHIFTED_POTENTIAL;
932 >          } else {
933 >            if (myMethod == "SHIFTED_FORCE") {            
934 >              esm = SHIFTED_FORCE;
935              } else {
936 <              // throw error        
937 <              sprintf( painCave.errMsg,
938 <                       "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() );
939 <              painCave.isFatal = 1;
940 <              simError();
941 <            }    
942 <          }          
936 >              if (myMethod == "REACTION_FIELD") {            
937 >                esm = REACTION_FIELD;
938 >              } else {
939 >                // throw error        
940 >                sprintf( painCave.errMsg,
941 >                         "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() );
942 >                painCave.isFatal = 1;
943 >                simError();
944 >              }    
945 >            }          
946 >          }
947 >        }
948 >      }
949 >    }
950 >    
951 >    if (simParams_->haveElectrostaticScreeningMethod()) {
952 >      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
953 >      toUpper(myScreen);
954 >      if (myScreen == "UNDAMPED") {
955 >        sm = UNDAMPED;
956 >      } else {
957 >        if (myScreen == "DAMPED") {
958 >          sm = DAMPED;
959 >          if (!simParams_->haveDampingAlpha()) {
960 >            //throw error
961 >            sprintf( painCave.errMsg,
962 >                     "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal);
963 >            painCave.isFatal = 0;
964 >            simError();
965 >          }
966 >        } else {
967 >          // throw error        
968 >          sprintf( painCave.errMsg,
969 >                   "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() );
970 >          painCave.isFatal = 1;
971 >          simError();
972          }
973        }
974      }
975 +    
976      // let's pass some summation method variables to fortran
977      setElectrostaticSummationMethod( &esm );
978 <    setDampedWolfAlpha( &alphaVal );
978 >    setScreeningMethod( &sm );
979 >    setDampingAlpha( &alphaVal );
980      setReactionFieldDielectric( &dielectric );
981      initFortranFF( &esm, &errorOut );
982 +  }
983 +
984 +  void SimInfo::setupSwitchingFunction() {    
985 +    int ft = CUBIC;
986 +
987 +    if (simParams_->haveSwitchingFunctionType()) {
988 +      std::string funcType = simParams_->getSwitchingFunctionType();
989 +      toUpper(funcType);
990 +      if (funcType == "CUBIC") {
991 +        ft = CUBIC;
992 +      } else {
993 +        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
994 +          ft = FIFTH_ORDER_POLY;
995 +        } else {
996 +          // throw error        
997 +          sprintf( painCave.errMsg,
998 +                   "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
999 +          painCave.isFatal = 1;
1000 +          simError();
1001 +        }          
1002 +      }
1003 +    }
1004 +
1005 +    // send switching function notification to switcheroo
1006 +    setFunctionType(&ft);
1007 +
1008    }
1009  
1010    void SimInfo::addProperty(GenericData* genData) {

Diff Legend

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