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 629 by chuckv, Mon Sep 26 15:58:17 2005 UTC vs.
Revision 720 by chrisfen, Tue Nov 8 13:32:06 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/doForces_interface.h"
59   #include "UseTheForce/DarkSide/electrostatic_interface.h"
60   #include "UseTheForce/notifyCutoffs_interface.h"
# Line 83 | Line 84 | namespace oopse {
84        MoleculeStamp* molStamp;
85        int nMolWithSameStamp;
86        int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
87 <      int nGroups = 0;          //total cutoff groups defined in meta-data file
87 >      int nGroups = 0;      //total cutoff groups defined in meta-data file
88        CutoffGroupStamp* cgStamp;    
89        RigidBodyStamp* rbStamp;
90        int nRigidAtoms = 0;
# Line 108 | Line 109 | namespace oopse {
109          }
110  
111          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
112 +
113          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
114  
115          //calculate atoms in rigid bodies
# Line 124 | Line 126 | namespace oopse {
126          
127        }
128  
129 <      //every free atom (atom does not belong to cutoff groups) is a cutoff group
130 <      //therefore the total number of cutoff groups in the system is equal to
131 <      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
132 <      //file plus the number of cutoff groups defined in meta-data file
129 >      //every free atom (atom does not belong to cutoff groups) is a cutoff
130 >      //group therefore the total number of cutoff groups in the system is
131 >      //equal to the total number of atoms minus number of atoms belong to
132 >      //cutoff group defined in meta-data file plus the number of cutoff
133 >      //groups defined in meta-data file
134        nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
135  
136 <      //every free atom (atom does not belong to rigid bodies) is an integrable object
137 <      //therefore the total number of  integrable objects in the system is equal to
138 <      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
139 <      //file plus the number of  rigid bodies defined in meta-data file
140 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
141 <
136 >      //every free atom (atom does not belong to rigid bodies) is an
137 >      //integrable object therefore the total number of integrable objects
138 >      //in the system is equal to the total number of atoms minus number of
139 >      //atoms belong to rigid body defined in meta-data file plus the number
140 >      //of rigid bodies defined in meta-data file
141 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
142 >                                                + nGlobalRigidBodies_;
143 >  
144        nGlobalMols_ = molStampIds_.size();
145  
146   #ifdef IS_MPI    
# Line 522 | Line 527 | namespace oopse {
527      int useDirectionalAtom = 0;    
528      int useElectrostatics = 0;
529      //usePBC and useRF are from simParams
530 <    int usePBC = simParams_->getPBC();
530 >    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
531      int useRF;
532 +    int useSF;
533 +    std::string myMethod;
534  
535      // set the useRF logical
536 <    std::string myMethod = simParams_->getElectrostaticSummationMethod();
537 <    if (myMethod == "REACTION_FIELD")
531 <      useRF = 1;
532 <    else
533 <      useRF = 0;
536 >    useRF = 0;
537 >    useSF = 0;
538  
539 +
540 +    if (simParams_->haveElectrostaticSummationMethod()) {
541 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
542 +      toUpper(myMethod);
543 +      if (myMethod == "REACTION_FIELD") {
544 +        useRF=1;
545 +      } else {
546 +        if (myMethod == "SHIFTED_FORCE") {
547 +          useSF = 1;
548 +        }
549 +      }
550 +    }
551 +
552      //loop over all of the atom types
553      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
554        useLennardJones |= (*i)->isLennardJones();
# Line 596 | Line 613 | namespace oopse {
613      temp = useRF;
614      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
615  
616 +    temp = useSF;
617 +    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
618 +
619   #endif
620  
621      fInfo_.SIM_uses_PBC = usePBC;    
# Line 611 | Line 631 | namespace oopse {
631      fInfo_.SIM_uses_Shapes = useShape;
632      fInfo_.SIM_uses_FLARB = useFLARB;
633      fInfo_.SIM_uses_RF = useRF;
634 +    fInfo_.SIM_uses_SF = useSF;
635  
636 <    if( fInfo_.SIM_uses_Dipoles && myMethod == "REACTION_FIELD") {
637 <
636 >    if( myMethod == "REACTION_FIELD") {
637 >      
638        if (simParams_->haveDielectric()) {
639          fInfo_.dielect = simParams_->getDielectric();
640        } else {
# Line 623 | Line 644 | namespace oopse {
644                  "\tsetting a dielectric constant!\n");
645          painCave.isFatal = 1;
646          simError();
647 <      }
627 <        
628 <    } else {
629 <      fInfo_.dielect = 0.0;
647 >      }      
648      }
649  
650    }
# Line 662 | Line 680 | namespace oopse {
680  
681          totalMass = cg->getMass();
682          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
683 <          mfact.push_back(atom->getMass()/totalMass);
683 >          // Check for massless groups - set mfact to 1 if true
684 >          if (totalMass != 0)
685 >            mfact.push_back(atom->getMass()/totalMass);
686 >          else
687 >            mfact.push_back( 1.0 );
688          }
689  
690        }      
# Line 798 | Line 820 | namespace oopse {
820      
821      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
822          
823 <      if (!simParams_->haveRcut()){
823 >      if (!simParams_->haveCutoffRadius()){
824          sprintf(painCave.errMsg,
825                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
826                  "\tOOPSE will use a default value of 15.0 angstroms"
# Line 807 | Line 829 | namespace oopse {
829          simError();
830          rcut = 15.0;
831        } else{
832 <        rcut = simParams_->getRcut();
832 >        rcut = simParams_->getCutoffRadius();
833        }
834  
835 <      if (!simParams_->haveRsw()){
835 >      if (!simParams_->haveSwitchingRadius()){
836          sprintf(painCave.errMsg,
837                  "SimCreator Warning: No value was set for switchingRadius.\n"
838                  "\tOOPSE will use a default value of\n"
839 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
839 >                "\t0.85 * cutoffRadius for the switchingRadius\n");
840          painCave.isFatal = 0;
841          simError();
842 <        rsw = 0.95 * rcut;
842 >        rsw = 0.85 * rcut;
843        } else{
844 <        rsw = simParams_->getRsw();
844 >        rsw = simParams_->getSwitchingRadius();
845        }
846  
847      } else {
848        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
849        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
850          
851 <      if (simParams_->haveRcut()) {
852 <        rcut = simParams_->getRcut();
851 >      if (simParams_->haveCutoffRadius()) {
852 >        rcut = simParams_->getCutoffRadius();
853        } else {
854          //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
855          rcut = calcMaxCutoffRadius();
856        }
857  
858 <      if (simParams_->haveRsw()) {
859 <        rsw  = simParams_->getRsw();
858 >      if (simParams_->haveSwitchingRadius()) {
859 >        rsw  = simParams_->getSwitchingRadius();
860        } else {
861          rsw = rcut;
862        }
# Line 851 | Line 873 | namespace oopse {
873      int cp =  TRADITIONAL_CUTOFF_POLICY;
874      if (simParams_->haveCutoffPolicy()) {
875        std::string myPolicy = simParams_->getCutoffPolicy();
876 +      toUpper(myPolicy);
877        if (myPolicy == "MIX") {
878          cp = MIX_CUTOFF_POLICY;
879        } else {
# Line 877 | Line 900 | namespace oopse {
900  
901      notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
902      // also send cutoff notification to electrostatics
903 <    setElectrostaticCutoffRadius(&rcut_);
903 >    setElectrostaticCutoffRadius(&rcut_, &rsw_);
904    }
905  
906    void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
907      
908      int errorOut;
909      int esm =  NONE;
910 +    int sm = UNDAMPED;
911      double alphaVal;
912      double dielectric;
913  
# Line 893 | Line 917 | namespace oopse {
917  
918      if (simParams_->haveElectrostaticSummationMethod()) {
919        std::string myMethod = simParams_->getElectrostaticSummationMethod();
920 +      toUpper(myMethod);
921        if (myMethod == "NONE") {
922          esm = NONE;
923        } else {
924 <        if (myMethod == "UNDAMPED_WOLF") {
925 <          esm = UNDAMPED_WOLF;
924 >        if (myMethod == "SWITCHING_FUNCTION") {
925 >          esm = SWITCHING_FUNCTION;
926          } else {
927 <          if (myMethod == "DAMPED_WOLF") {            
928 <            esm = DAMPED_WOLF;
929 <            if (!simParams_->haveDampingAlpha()) {
930 <              //throw error
931 <              sprintf( painCave.errMsg,
907 <                       "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);
908 <              painCave.isFatal = 0;
909 <              simError();
910 <            }
911 <          } else {
912 <            if (myMethod == "REACTION_FIELD") {
913 <              esm = REACTION_FIELD;
927 >          if (myMethod == "SHIFTED_POTENTIAL") {
928 >            esm = SHIFTED_POTENTIAL;
929 >          } else {
930 >            if (myMethod == "SHIFTED_FORCE") {            
931 >              esm = SHIFTED_FORCE;
932              } else {
933 <              // throw error        
934 <              sprintf( painCave.errMsg,
935 <                       "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() );
936 <              painCave.isFatal = 1;
937 <              simError();
938 <            }    
939 <          }          
933 >              if (myMethod == "REACTION_FIELD") {            
934 >                esm = REACTION_FIELD;
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\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() );
939 >                painCave.isFatal = 1;
940 >                simError();
941 >              }    
942 >            }          
943 >          }
944          }
945        }
946      }
947 +    
948 +    if (simParams_->haveElectrostaticScreeningMethod()) {
949 +      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
950 +      toUpper(myScreen);
951 +      if (myScreen == "UNDAMPED") {
952 +        sm = UNDAMPED;
953 +      } else {
954 +        if (myScreen == "DAMPED") {
955 +          sm = DAMPED;
956 +          if (!simParams_->haveDampingAlpha()) {
957 +            //throw error
958 +            sprintf( painCave.errMsg,
959 +                     "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal);
960 +            painCave.isFatal = 0;
961 +            simError();
962 +          }
963 +        } else {
964 +          // throw error        
965 +          sprintf( painCave.errMsg,
966 +                   "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() );
967 +          painCave.isFatal = 1;
968 +          simError();
969 +        }
970 +      }
971 +    }
972 +    
973      // let's pass some summation method variables to fortran
974      setElectrostaticSummationMethod( &esm );
975 <    setDampedWolfAlpha( &alphaVal );
975 >    setScreeningMethod( &sm );
976 >    setDampingAlpha( &alphaVal );
977      setReactionFieldDielectric( &dielectric );
978      initFortranFF( &esm, &errorOut );
979    }

Diff Legend

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