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 555 by chuckv, Mon May 30 14:01:52 2005 UTC vs.
Revision 734 by chuckv, Tue Nov 15 16:05:38 2005 UTC

# Line 52 | Line 52
52   #include "brains/SimInfo.hpp"
53   #include "math/Vector3.hpp"
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 80 | 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 105 | Line 111 | namespace oopse {
111          }
112  
113          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
114 +
115          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
116  
117          //calculate atoms in rigid bodies
# Line 121 | 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 462 | Line 472 | namespace oopse {
472      //setup fortran force field
473      /** @deprecate */    
474      int isError = 0;
475 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
475 >    
476 >    setupElectrostaticSummationMethod( isError );
477 >    setupSwitchingFunction();
478 >
479      if(isError){
480        sprintf( painCave.errMsg,
481                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 506 | Line 519 | namespace oopse {
519      int useLennardJones = 0;
520      int useElectrostatic = 0;
521      int useEAM = 0;
522 +    int useSC = 0;
523      int useCharge = 0;
524      int useDirectional = 0;
525      int useDipole = 0;
# Line 517 | Line 531 | namespace oopse {
531      int useDirectionalAtom = 0;    
532      int useElectrostatics = 0;
533      //usePBC and useRF are from simParams
534 <    int usePBC = simParams_->getPBC();
535 <    int useRF = simParams_->getUseRF();
534 >    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
535 >    int useRF;
536 >    int useSF;
537 >    std::string myMethod;
538  
539 +    // set the useRF logical
540 +    useRF = 0;
541 +    useSF = 0;
542 +
543 +
544 +    if (simParams_->haveElectrostaticSummationMethod()) {
545 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
546 +      toUpper(myMethod);
547 +      if (myMethod == "REACTION_FIELD") {
548 +        useRF=1;
549 +      } else {
550 +        if (myMethod == "SHIFTED_FORCE") {
551 +          useSF = 1;
552 +        }
553 +      }
554 +    }
555 +
556      //loop over all of the atom types
557      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
558        useLennardJones |= (*i)->isLennardJones();
559        useElectrostatic |= (*i)->isElectrostatic();
560        useEAM |= (*i)->isEAM();
561 +      useSC |= (*i)->isSC();
562        useCharge |= (*i)->isCharge();
563        useDirectional |= (*i)->isDirectional();
564        useDipole |= (*i)->isDipole();
# Line 575 | Line 609 | namespace oopse {
609      temp = useEAM;
610      MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
611  
612 +    temp = useSC;
613 +    MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
614 +    
615      temp = useShape;
616      MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
617  
# Line 583 | Line 620 | namespace oopse {
620  
621      temp = useRF;
622      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
623 <    
623 >
624 >    temp = useSF;
625 >    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
626 >
627   #endif
628  
629      fInfo_.SIM_uses_PBC = usePBC;    
# Line 596 | Line 636 | namespace oopse {
636      fInfo_.SIM_uses_StickyPower = useStickyPower;
637      fInfo_.SIM_uses_GayBerne = useGayBerne;
638      fInfo_.SIM_uses_EAM = useEAM;
639 +    fInfo_.SIM_uses_SC = useSC;
640      fInfo_.SIM_uses_Shapes = useShape;
641      fInfo_.SIM_uses_FLARB = useFLARB;
642      fInfo_.SIM_uses_RF = useRF;
643 +    fInfo_.SIM_uses_SF = useSF;
644  
645 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
646 <
645 >    if( myMethod == "REACTION_FIELD") {
646 >      
647        if (simParams_->haveDielectric()) {
648          fInfo_.dielect = simParams_->getDielectric();
649        } else {
# Line 611 | Line 653 | namespace oopse {
653                  "\tsetting a dielectric constant!\n");
654          painCave.isFatal = 1;
655          simError();
656 <      }
615 <        
616 <    } else {
617 <      fInfo_.dielect = 0.0;
656 >      }      
657      }
658  
659    }
# Line 650 | Line 689 | namespace oopse {
689  
690          totalMass = cg->getMass();
691          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
692 <          mfact.push_back(atom->getMass()/totalMass);
692 >          // Check for massless groups - set mfact to 1 if true
693 >          if (totalMass != 0)
694 >            mfact.push_back(atom->getMass()/totalMass);
695 >          else
696 >            mfact.push_back( 1.0 );
697          }
698  
699        }      
# Line 786 | Line 829 | namespace oopse {
829      
830      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
831          
832 <      if (!simParams_->haveRcut()){
832 >      if (!simParams_->haveCutoffRadius()){
833          sprintf(painCave.errMsg,
834                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
835                  "\tOOPSE will use a default value of 15.0 angstroms"
# Line 795 | Line 838 | namespace oopse {
838          simError();
839          rcut = 15.0;
840        } else{
841 <        rcut = simParams_->getRcut();
841 >        rcut = simParams_->getCutoffRadius();
842        }
843  
844 <      if (!simParams_->haveRsw()){
844 >      if (!simParams_->haveSwitchingRadius()){
845          sprintf(painCave.errMsg,
846                  "SimCreator Warning: No value was set for switchingRadius.\n"
847                  "\tOOPSE will use a default value of\n"
848 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
848 >                "\t0.85 * cutoffRadius for the switchingRadius\n");
849          painCave.isFatal = 0;
850          simError();
851 <        rsw = 0.95 * rcut;
851 >        rsw = 0.85 * rcut;
852        } else{
853 <        rsw = simParams_->getRsw();
853 >        rsw = simParams_->getSwitchingRadius();
854        }
855  
856      } else {
857        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
858        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
859          
860 <      if (simParams_->haveRcut()) {
861 <        rcut = simParams_->getRcut();
860 >      if (simParams_->haveCutoffRadius()) {
861 >        rcut = simParams_->getCutoffRadius();
862        } else {
863          //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
864          rcut = calcMaxCutoffRadius();
865        }
866  
867 <      if (simParams_->haveRsw()) {
868 <        rsw  = simParams_->getRsw();
867 >      if (simParams_->haveSwitchingRadius()) {
868 >        rsw  = simParams_->getSwitchingRadius();
869        } else {
870          rsw = rcut;
871        }
# Line 830 | Line 873 | namespace oopse {
873      }
874    }
875  
876 <  void SimInfo::setupCutoff() {
876 >  void SimInfo::setupCutoff() {    
877      getCutoff(rcut_, rsw_);    
878      double rnblist = rcut_ + 1; // skin of neighbor list
879  
880      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
881 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
881 >    
882 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
883 >    if (simParams_->haveCutoffPolicy()) {
884 >      std::string myPolicy = simParams_->getCutoffPolicy();
885 >      toUpper(myPolicy);
886 >      if (myPolicy == "MIX") {
887 >        cp = MIX_CUTOFF_POLICY;
888 >      } else {
889 >        if (myPolicy == "MAX") {
890 >          cp = MAX_CUTOFF_POLICY;
891 >        } else {
892 >          if (myPolicy == "TRADITIONAL") {            
893 >            cp = TRADITIONAL_CUTOFF_POLICY;
894 >          } else {
895 >            // throw error        
896 >            sprintf( painCave.errMsg,
897 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
898 >            painCave.isFatal = 1;
899 >            simError();
900 >          }    
901 >        }          
902 >      }
903 >    }
904 >
905 >
906 >    if (simParams_->haveSkinThickness()) {
907 >      double skinThickness = simParams_->getSkinThickness();
908 >    }
909 >
910 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
911 >    // also send cutoff notification to electrostatics
912 >    setElectrostaticCutoffRadius(&rcut_, &rsw_);
913    }
914  
915 +  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
916 +    
917 +    int errorOut;
918 +    int esm =  NONE;
919 +    int sm = UNDAMPED;
920 +    double alphaVal;
921 +    double dielectric;
922 +
923 +    errorOut = isError;
924 +    alphaVal = simParams_->getDampingAlpha();
925 +    dielectric = simParams_->getDielectric();
926 +
927 +    if (simParams_->haveElectrostaticSummationMethod()) {
928 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
929 +      toUpper(myMethod);
930 +      if (myMethod == "NONE") {
931 +        esm = NONE;
932 +      } else {
933 +        if (myMethod == "SWITCHING_FUNCTION") {
934 +          esm = SWITCHING_FUNCTION;
935 +        } else {
936 +          if (myMethod == "SHIFTED_POTENTIAL") {
937 +            esm = SHIFTED_POTENTIAL;
938 +          } else {
939 +            if (myMethod == "SHIFTED_FORCE") {            
940 +              esm = SHIFTED_FORCE;
941 +            } else {
942 +              if (myMethod == "REACTION_FIELD") {            
943 +                esm = REACTION_FIELD;
944 +              } else {
945 +                // throw error        
946 +                sprintf( painCave.errMsg,
947 +                         "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() );
948 +                painCave.isFatal = 1;
949 +                simError();
950 +              }    
951 +            }          
952 +          }
953 +        }
954 +      }
955 +    }
956 +    
957 +    if (simParams_->haveElectrostaticScreeningMethod()) {
958 +      std::string myScreen = simParams_->getElectrostaticScreeningMethod();
959 +      toUpper(myScreen);
960 +      if (myScreen == "UNDAMPED") {
961 +        sm = UNDAMPED;
962 +      } else {
963 +        if (myScreen == "DAMPED") {
964 +          sm = DAMPED;
965 +          if (!simParams_->haveDampingAlpha()) {
966 +            //throw error
967 +            sprintf( painCave.errMsg,
968 +                     "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal);
969 +            painCave.isFatal = 0;
970 +            simError();
971 +          }
972 +        } else {
973 +          // throw error        
974 +          sprintf( painCave.errMsg,
975 +                   "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() );
976 +          painCave.isFatal = 1;
977 +          simError();
978 +        }
979 +      }
980 +    }
981 +    
982 +    // let's pass some summation method variables to fortran
983 +    setElectrostaticSummationMethod( &esm );
984 +    setScreeningMethod( &sm );
985 +    setDampingAlpha( &alphaVal );
986 +    setReactionFieldDielectric( &dielectric );
987 +    initFortranFF( &esm, &errorOut );
988 +  }
989 +
990 +  void SimInfo::setupSwitchingFunction() {    
991 +    int ft = CUBIC;
992 +
993 +    if (simParams_->haveSwitchingFunctionType()) {
994 +      std::string funcType = simParams_->getSwitchingFunctionType();
995 +      toUpper(funcType);
996 +      if (funcType == "CUBIC") {
997 +        ft = CUBIC;
998 +      } else {
999 +        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
1000 +          ft = FIFTH_ORDER_POLY;
1001 +        } else {
1002 +          // throw error        
1003 +          sprintf( painCave.errMsg,
1004 +                   "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
1005 +          painCave.isFatal = 1;
1006 +          simError();
1007 +        }          
1008 +      }
1009 +    }
1010 +
1011 +    // send switching function notification to switcheroo
1012 +    setFunctionType(&ft);
1013 +
1014 +  }
1015 +
1016    void SimInfo::addProperty(GenericData* genData) {
1017      properties_.addProperty(genData);  
1018    }
# Line 981 | Line 1156 | namespace oopse {
1156    
1157     /*
1158     Return intertia tensor for entire system and angular momentum Vector.
1159 +
1160 +
1161 +       [  Ixx -Ixy  -Ixz ]
1162 +  J =| -Iyx  Iyy  -Iyz |
1163 +       [ -Izx -Iyz   Izz ]
1164      */
1165  
1166     void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
# Line 1032 | Line 1212 | namespace oopse {
1212        inertiaTensor(0,1) = -xy;
1213        inertiaTensor(0,2) = -xz;
1214        inertiaTensor(1,0) = -xy;
1215 <      inertiaTensor(2,0) = xx + zz;
1215 >      inertiaTensor(1,1) = xx + zz;
1216        inertiaTensor(1,2) = -yz;
1217        inertiaTensor(2,0) = -xz;
1218        inertiaTensor(2,1) = -yz;
# Line 1060 | Line 1240 | namespace oopse {
1240        SimInfo::MoleculeIterator i;
1241        Molecule* mol;
1242        
1243 <      Vector3d thisq(0.0);
1244 <      Vector3d thisv(0.0);
1243 >      Vector3d thisr(0.0);
1244 >      Vector3d thisp(0.0);
1245        
1246 <      double thisMass = 0.0;
1246 >      double thisMass;
1247        
1248        for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1249 <         thisq = mol->getCom()-com;
1250 <         thisv = mol->getComVel()-comVel;
1251 <         thisMass = mol->getMass();
1072 <         angularMomentum += cross( thisq, thisv ) * thisMass;
1249 >        thisMass = mol->getMass();
1250 >        thisr = mol->getCom()-com;
1251 >        thisp = (mol->getComVel()-comVel)*thisMass;
1252          
1253 +        angularMomentum += cross( thisr, thisp );
1254 +        
1255        }  
1256        
1257   #ifdef IS_MPI

Diff Legend

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