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 586 by gezelter, Wed Sep 7 20:46:46 2005 UTC vs.
Revision 701 by chrisfen, Wed Oct 26 23:32:25 2005 UTC

# Line 53 | Line 53
53   #include "math/Vector3.hpp"
54   #include "primitives/Molecule.hpp"
55   #include "UseTheForce/fCutoffPolicy.h"
56 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
57   #include "UseTheForce/doForces_interface.h"
58 + #include "UseTheForce/DarkSide/electrostatic_interface.h"
59   #include "UseTheForce/notifyCutoffs_interface.h"
60   #include "utils/MemoryUtils.hpp"
61   #include "utils/simError.h"
# Line 81 | Line 83 | namespace oopse {
83        MoleculeStamp* molStamp;
84        int nMolWithSameStamp;
85        int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
86 <      int nGroups = 0;          //total cutoff groups defined in meta-data file
86 >      int nGroups = 0;      //total cutoff groups defined in meta-data file
87        CutoffGroupStamp* cgStamp;    
88        RigidBodyStamp* rbStamp;
89        int nRigidAtoms = 0;
# Line 106 | Line 108 | namespace oopse {
108          }
109  
110          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
111 +
112          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
113  
114          //calculate atoms in rigid bodies
# Line 122 | Line 125 | namespace oopse {
125          
126        }
127  
128 <      //every free atom (atom does not belong to cutoff groups) is a cutoff group
129 <      //therefore the total number of cutoff groups in the system is equal to
130 <      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
131 <      //file plus the number of cutoff groups defined in meta-data file
128 >      //every free atom (atom does not belong to cutoff groups) is a cutoff
129 >      //group therefore the total number of cutoff groups in the system is
130 >      //equal to the total number of atoms minus number of atoms belong to
131 >      //cutoff group defined in meta-data file plus the number of cutoff
132 >      //groups defined in meta-data file
133        nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
130 –
131 –      //every free atom (atom does not belong to rigid bodies) is an integrable object
132 –      //therefore the total number of  integrable objects in the system is equal to
133 –      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
134 –      //file plus the number of  rigid bodies defined in meta-data file
135 –      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
134  
135 +      //every free atom (atom does not belong to rigid bodies) is an
136 +      //integrable object therefore the total number of integrable objects
137 +      //in the system is equal to the total number of atoms minus number of
138 +      //atoms belong to rigid body defined in meta-data file plus the number
139 +      //of rigid bodies defined in meta-data file
140 +      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
141 +                                                + nGlobalRigidBodies_;
142 +  
143        nGlobalMols_ = molStampIds_.size();
144  
145   #ifdef IS_MPI    
# Line 463 | Line 469 | namespace oopse {
469      //setup fortran force field
470      /** @deprecate */    
471      int isError = 0;
472 <    initFortranFF( &fInfo_.SIM_uses_RF, &fInfo_.SIM_uses_UW,
473 <                   &fInfo_.SIM_uses_DW, &isError );
472 >    
473 >    setupElectrostaticSummationMethod( isError );
474 >
475      if(isError){
476        sprintf( painCave.errMsg,
477                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 519 | Line 526 | namespace oopse {
526      int useDirectionalAtom = 0;    
527      int useElectrostatics = 0;
528      //usePBC and useRF are from simParams
529 <    int usePBC = simParams_->getPBC();
530 <    int useRF = simParams_->getUseRF();
531 <    int useUW = simParams_->getUseUndampedWolf();
532 <    int useDW = simParams_->getUseDampedWolf();
529 >    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
530 >    int useRF;
531 >    std::string myMethod;
532 >
533 >    // set the useRF logical
534 >    useRF = 0;
535  
536 +
537 +    if (simParams_->haveElectrostaticSummationMethod()) {
538 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
539 +      toUpper(myMethod);
540 +      if (myMethod == "REACTION_FIELD") {
541 +        useRF=1;
542 +      }
543 +    }
544 +
545      //loop over all of the atom types
546      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
547        useLennardJones |= (*i)->isLennardJones();
# Line 588 | Line 606 | namespace oopse {
606      temp = useRF;
607      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
608  
591 –    temp = useUW;
592 –    MPI_Allreduce(&temp, &useUW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
593 –
594 –    temp = useDW;
595 –    MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
596 –    
609   #endif
610  
611      fInfo_.SIM_uses_PBC = usePBC;    
# Line 609 | Line 621 | namespace oopse {
621      fInfo_.SIM_uses_Shapes = useShape;
622      fInfo_.SIM_uses_FLARB = useFLARB;
623      fInfo_.SIM_uses_RF = useRF;
612 –    fInfo_.SIM_uses_UW = useUW;
613 –    fInfo_.SIM_uses_DW = useDW;
624  
625 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
626 <
625 >    if( myMethod == "REACTION_FIELD") {
626 >      
627        if (simParams_->haveDielectric()) {
628          fInfo_.dielect = simParams_->getDielectric();
629        } else {
# Line 623 | Line 633 | namespace oopse {
633                  "\tsetting a dielectric constant!\n");
634          painCave.isFatal = 1;
635          simError();
636 <      }
627 <        
628 <    } else {
629 <      fInfo_.dielect = 0.0;
636 >      }      
637      }
631 –
638    }
639  
640    void SimInfo::setupFortranSim() {
# Line 662 | Line 668 | namespace oopse {
668  
669          totalMass = cg->getMass();
670          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
671 <          mfact.push_back(atom->getMass()/totalMass);
671 >          // Check for massless groups - set mfact to 1 if true
672 >          if (totalMass != 0)
673 >            mfact.push_back(atom->getMass()/totalMass);
674 >          else
675 >            mfact.push_back( 1.0 );
676          }
677  
678        }      
# Line 798 | Line 808 | namespace oopse {
808      
809      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
810          
811 <      if (!simParams_->haveRcut()){
811 >      if (!simParams_->haveCutoffRadius()){
812          sprintf(painCave.errMsg,
813                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
814                  "\tOOPSE will use a default value of 15.0 angstroms"
# Line 807 | Line 817 | namespace oopse {
817          simError();
818          rcut = 15.0;
819        } else{
820 <        rcut = simParams_->getRcut();
820 >        rcut = simParams_->getCutoffRadius();
821        }
822  
823 <      if (!simParams_->haveRsw()){
823 >      if (!simParams_->haveSwitchingRadius()){
824          sprintf(painCave.errMsg,
825                  "SimCreator Warning: No value was set for switchingRadius.\n"
826                  "\tOOPSE will use a default value of\n"
827 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
827 >                "\t0.85 * cutoffRadius for the switchingRadius\n");
828          painCave.isFatal = 0;
829          simError();
830 <        rsw = 0.95 * rcut;
830 >        rsw = 0.85 * rcut;
831        } else{
832 <        rsw = simParams_->getRsw();
832 >        rsw = simParams_->getSwitchingRadius();
833        }
834  
835      } else {
836        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
837        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
838          
839 <      if (simParams_->haveRcut()) {
840 <        rcut = simParams_->getRcut();
839 >      if (simParams_->haveCutoffRadius()) {
840 >        rcut = simParams_->getCutoffRadius();
841        } else {
842          //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
843          rcut = calcMaxCutoffRadius();
844        }
845  
846 <      if (simParams_->haveRsw()) {
847 <        rsw  = simParams_->getRsw();
846 >      if (simParams_->haveSwitchingRadius()) {
847 >        rsw  = simParams_->getSwitchingRadius();
848        } else {
849          rsw = rcut;
850        }
# Line 851 | Line 861 | namespace oopse {
861      int cp =  TRADITIONAL_CUTOFF_POLICY;
862      if (simParams_->haveCutoffPolicy()) {
863        std::string myPolicy = simParams_->getCutoffPolicy();
864 +      toUpper(myPolicy);
865        if (myPolicy == "MIX") {
866          cp = MIX_CUTOFF_POLICY;
867        } else {
# Line 869 | Line 880 | namespace oopse {
880          }          
881        }
882      }
883 +
884 +
885 +    if (simParams_->haveSkinThickness()) {
886 +      double skinThickness = simParams_->getSkinThickness();
887 +    }
888 +
889      notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
890 +    // also send cutoff notification to electrostatics
891 +    setElectrostaticCutoffRadius(&rcut_, &rsw_);
892 +  }
893 +
894 +  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
895 +    
896 +    int errorOut;
897 +    int esm =  NONE;
898 +    double alphaVal;
899 +    double dielectric;
900 +
901 +    errorOut = isError;
902 +    alphaVal = simParams_->getDampingAlpha();
903 +    dielectric = simParams_->getDielectric();
904 +
905 +    if (simParams_->haveElectrostaticSummationMethod()) {
906 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
907 +      toUpper(myMethod);
908 +      if (myMethod == "NONE") {
909 +        esm = NONE;
910 +      } else {
911 +        if (myMethod == "UNDAMPED_WOLF") {
912 +          esm = UNDAMPED_WOLF;
913 +        } else {
914 +          if (myMethod == "DAMPED_WOLF") {            
915 +            esm = DAMPED_WOLF;
916 +            if (!simParams_->haveDampingAlpha()) {
917 +              //throw error
918 +              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;
926 +            } else {
927 +              // throw error        
928 +              sprintf( painCave.errMsg,
929 +                       "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() );
930 +              painCave.isFatal = 1;
931 +              simError();
932 +            }    
933 +          }          
934 +        }
935 +      }
936 +    }
937 +    // let's pass some summation method variables to fortran
938 +    setElectrostaticSummationMethod( &esm );
939 +    setDampedWolfAlpha( &alphaVal );
940 +    setReactionFieldDielectric( &dielectric );
941 +    initFortranFF( &esm, &errorOut );
942    }
943  
944    void SimInfo::addProperty(GenericData* genData) {

Diff Legend

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