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 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 726 by chrisfen, Fri Nov 11 15:22:11 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;
129 –
130 –      //every free atom (atom does not belong to rigid bodies) is an integrable object
131 –      //therefore the total number of  integrable objects in the system is equal to
132 –      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
133 –      //file plus the number of  rigid bodies defined in meta-data file
134 –      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
137  
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 511 | Line 524 | namespace oopse {
524      int useDipole = 0;
525      int useGayBerne = 0;
526      int useSticky = 0;
527 +    int useStickyPower = 0;
528      int useShape = 0;
529      int useFLARB = 0; //it is not in AtomType yet
530      int useDirectionalAtom = 0;    
531      int useElectrostatics = 0;
532      //usePBC and useRF are from simParams
533 <    int usePBC = simParams_->getPBC();
534 <    int useRF = simParams_->getUseRF();
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) {
557        useLennardJones |= (*i)->isLennardJones();
# Line 529 | Line 562 | namespace oopse {
562        useDipole |= (*i)->isDipole();
563        useGayBerne |= (*i)->isGayBerne();
564        useSticky |= (*i)->isSticky();
565 +      useStickyPower |= (*i)->isStickyPower();
566        useShape |= (*i)->isShape();
567      }
568  
569 <    if (useSticky || useDipole || useGayBerne || useShape) {
569 >    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
570        useDirectionalAtom = 1;
571      }
572  
# Line 564 | Line 598 | namespace oopse {
598      temp = useSticky;
599      MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
600  
601 +    temp = useStickyPower;
602 +    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
603 +    
604      temp = useGayBerne;
605      MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
606  
# Line 578 | Line 615 | namespace oopse {
615  
616      temp = useRF;
617      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
618 <    
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 588 | Line 628 | namespace oopse {
628      fInfo_.SIM_uses_Charges = useCharge;
629      fInfo_.SIM_uses_Dipoles = useDipole;
630      fInfo_.SIM_uses_Sticky = useSticky;
631 +    fInfo_.SIM_uses_StickyPower = useStickyPower;
632      fInfo_.SIM_uses_GayBerne = useGayBerne;
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 605 | Line 647 | namespace oopse {
647                  "\tsetting a dielectric constant!\n");
648          painCave.isFatal = 1;
649          simError();
650 <      }
609 <        
610 <    } else {
611 <      fInfo_.dielect = 0.0;
650 >      }      
651      }
652  
653    }
# Line 644 | 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 780 | 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 789 | 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 824 | Line 867 | namespace oopse {
867      }
868    }
869  
870 <  void SimInfo::setupCutoff() {
870 >  void SimInfo::setupCutoff() {    
871      getCutoff(rcut_, rsw_);    
872      double rnblist = rcut_ + 1; // skin of neighbor list
873  
874      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
875 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
875 >    
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 {
883 >        if (myPolicy == "MAX") {
884 >          cp = MAX_CUTOFF_POLICY;
885 >        } else {
886 >          if (myPolicy == "TRADITIONAL") {            
887 >            cp = TRADITIONAL_CUTOFF_POLICY;
888 >          } else {
889 >            // throw error        
890 >            sprintf( painCave.errMsg,
891 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
892 >            painCave.isFatal = 1;
893 >            simError();
894 >          }    
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_, &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 +
917 +    errorOut = isError;
918 +    alphaVal = simParams_->getDampingAlpha();
919 +    dielectric = simParams_->getDielectric();
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 == "SWITCHING_FUNCTION") {
928 +          esm = SWITCHING_FUNCTION;
929 +        } else {
930 +          if (myMethod == "SHIFTED_POTENTIAL") {
931 +            esm = SHIFTED_POTENTIAL;
932 +          } else {
933 +            if (myMethod == "SHIFTED_FORCE") {            
934 +              esm = SHIFTED_FORCE;
935 +            } else {
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 +    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) {
1011      properties_.addProperty(genData);  
1012    }
# Line 939 | Line 1114 | namespace oopse {
1114  
1115      return o;
1116    }
1117 +  
1118 +  
1119 +   /*
1120 +   Returns center of mass and center of mass velocity in one function call.
1121 +   */
1122 +  
1123 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1124 +      SimInfo::MoleculeIterator i;
1125 +      Molecule* mol;
1126 +      
1127 +    
1128 +      double totalMass = 0.0;
1129 +    
1130  
1131 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1132 +         double mass = mol->getMass();
1133 +         totalMass += mass;
1134 +         com += mass * mol->getCom();
1135 +         comVel += mass * mol->getComVel();          
1136 +      }  
1137 +      
1138 + #ifdef IS_MPI
1139 +      double tmpMass = totalMass;
1140 +      Vector3d tmpCom(com);  
1141 +      Vector3d tmpComVel(comVel);
1142 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1143 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1144 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1145 + #endif
1146 +      
1147 +      com /= totalMass;
1148 +      comVel /= totalMass;
1149 +   }        
1150 +  
1151 +   /*
1152 +   Return intertia tensor for entire system and angular momentum Vector.
1153 +
1154 +
1155 +       [  Ixx -Ixy  -Ixz ]
1156 +  J =| -Iyx  Iyy  -Iyz |
1157 +       [ -Izx -Iyz   Izz ]
1158 +    */
1159 +
1160 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1161 +      
1162 +
1163 +      double xx = 0.0;
1164 +      double yy = 0.0;
1165 +      double zz = 0.0;
1166 +      double xy = 0.0;
1167 +      double xz = 0.0;
1168 +      double yz = 0.0;
1169 +      Vector3d com(0.0);
1170 +      Vector3d comVel(0.0);
1171 +      
1172 +      getComAll(com, comVel);
1173 +      
1174 +      SimInfo::MoleculeIterator i;
1175 +      Molecule* mol;
1176 +      
1177 +      Vector3d thisq(0.0);
1178 +      Vector3d thisv(0.0);
1179 +
1180 +      double thisMass = 0.0;
1181 +    
1182 +      
1183 +      
1184 +  
1185 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1186 +        
1187 +         thisq = mol->getCom()-com;
1188 +         thisv = mol->getComVel()-comVel;
1189 +         thisMass = mol->getMass();
1190 +         // Compute moment of intertia coefficients.
1191 +         xx += thisq[0]*thisq[0]*thisMass;
1192 +         yy += thisq[1]*thisq[1]*thisMass;
1193 +         zz += thisq[2]*thisq[2]*thisMass;
1194 +        
1195 +         // compute products of intertia
1196 +         xy += thisq[0]*thisq[1]*thisMass;
1197 +         xz += thisq[0]*thisq[2]*thisMass;
1198 +         yz += thisq[1]*thisq[2]*thisMass;
1199 +            
1200 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1201 +            
1202 +      }  
1203 +      
1204 +      
1205 +      inertiaTensor(0,0) = yy + zz;
1206 +      inertiaTensor(0,1) = -xy;
1207 +      inertiaTensor(0,2) = -xz;
1208 +      inertiaTensor(1,0) = -xy;
1209 +      inertiaTensor(1,1) = xx + zz;
1210 +      inertiaTensor(1,2) = -yz;
1211 +      inertiaTensor(2,0) = -xz;
1212 +      inertiaTensor(2,1) = -yz;
1213 +      inertiaTensor(2,2) = xx + yy;
1214 +      
1215 + #ifdef IS_MPI
1216 +      Mat3x3d tmpI(inertiaTensor);
1217 +      Vector3d tmpAngMom;
1218 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1219 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1220 + #endif
1221 +              
1222 +      return;
1223 +   }
1224 +
1225 +   //Returns the angular momentum of the system
1226 +   Vector3d SimInfo::getAngularMomentum(){
1227 +      
1228 +      Vector3d com(0.0);
1229 +      Vector3d comVel(0.0);
1230 +      Vector3d angularMomentum(0.0);
1231 +      
1232 +      getComAll(com,comVel);
1233 +      
1234 +      SimInfo::MoleculeIterator i;
1235 +      Molecule* mol;
1236 +      
1237 +      Vector3d thisr(0.0);
1238 +      Vector3d thisp(0.0);
1239 +      
1240 +      double thisMass;
1241 +      
1242 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1243 +        thisMass = mol->getMass();
1244 +        thisr = mol->getCom()-com;
1245 +        thisp = (mol->getComVel()-comVel)*thisMass;
1246 +        
1247 +        angularMomentum += cross( thisr, thisp );
1248 +        
1249 +      }  
1250 +      
1251 + #ifdef IS_MPI
1252 +      Vector3d tmpAngMom;
1253 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1254 + #endif
1255 +      
1256 +      return angularMomentum;
1257 +   }
1258 +  
1259 +  
1260   }//end namespace oopse
1261  

Diff Legend

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