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 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 604 by chrisfen, Fri Sep 16 19:00:12 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/doForces_interface.h"
58   #include "UseTheForce/notifyCutoffs_interface.h"
59   #include "utils/MemoryUtils.hpp"
# Line 462 | Line 464 | namespace oopse {
464      //setup fortran force field
465      /** @deprecate */    
466      int isError = 0;
467 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
467 >    
468 >    setupElectrostaticSummationMethod( isError );
469 >
470      if(isError){
471        sprintf( painCave.errMsg,
472                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 511 | Line 515 | namespace oopse {
515      int useDipole = 0;
516      int useGayBerne = 0;
517      int useSticky = 0;
518 +    int useStickyPower = 0;
519      int useShape = 0;
520      int useFLARB = 0; //it is not in AtomType yet
521      int useDirectionalAtom = 0;    
522      int useElectrostatics = 0;
523      //usePBC and useRF are from simParams
524      int usePBC = simParams_->getPBC();
520 –    int useRF = simParams_->getUseRF();
525  
526      //loop over all of the atom types
527      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
# Line 529 | Line 533 | namespace oopse {
533        useDipole |= (*i)->isDipole();
534        useGayBerne |= (*i)->isGayBerne();
535        useSticky |= (*i)->isSticky();
536 +      useStickyPower |= (*i)->isStickyPower();
537        useShape |= (*i)->isShape();
538      }
539  
540 <    if (useSticky || useDipole || useGayBerne || useShape) {
540 >    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
541        useDirectionalAtom = 1;
542      }
543  
# Line 564 | Line 569 | namespace oopse {
569      temp = useSticky;
570      MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
571  
572 +    temp = useStickyPower;
573 +    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
574 +    
575      temp = useGayBerne;
576      MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
577  
# Line 576 | Line 584 | namespace oopse {
584      temp = useFLARB;
585      MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
586  
579 –    temp = useRF;
580 –    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
581 –    
587   #endif
588  
589      fInfo_.SIM_uses_PBC = usePBC;    
# Line 588 | Line 593 | namespace oopse {
593      fInfo_.SIM_uses_Charges = useCharge;
594      fInfo_.SIM_uses_Dipoles = useDipole;
595      fInfo_.SIM_uses_Sticky = useSticky;
596 +    fInfo_.SIM_uses_StickyPower = useStickyPower;
597      fInfo_.SIM_uses_GayBerne = useGayBerne;
598      fInfo_.SIM_uses_EAM = useEAM;
599      fInfo_.SIM_uses_Shapes = useShape;
600      fInfo_.SIM_uses_FLARB = useFLARB;
595 –    fInfo_.SIM_uses_RF = useRF;
601  
602      if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
603  
# Line 824 | Line 829 | namespace oopse {
829      }
830    }
831  
832 <  void SimInfo::setupCutoff() {
832 >  void SimInfo::setupCutoff() {    
833      getCutoff(rcut_, rsw_);    
834      double rnblist = rcut_ + 1; // skin of neighbor list
835  
836      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
837 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
837 >    
838 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
839 >    if (simParams_->haveCutoffPolicy()) {
840 >      std::string myPolicy = simParams_->getCutoffPolicy();
841 >      if (myPolicy == "MIX") {
842 >        cp = MIX_CUTOFF_POLICY;
843 >      } else {
844 >        if (myPolicy == "MAX") {
845 >          cp = MAX_CUTOFF_POLICY;
846 >        } else {
847 >          if (myPolicy == "TRADITIONAL") {            
848 >            cp = TRADITIONAL_CUTOFF_POLICY;
849 >          } else {
850 >            // throw error        
851 >            sprintf( painCave.errMsg,
852 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
853 >            painCave.isFatal = 1;
854 >            simError();
855 >          }    
856 >        }          
857 >      }
858 >    }
859 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
860    }
861  
862 +  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
863 +    
864 +    int errorOut;
865 +    int esm =  NONE;
866 +    double alphaVal;
867 +
868 +    errorOut = isError;
869 +
870 +    if (simParams_->haveElectrostaticSummationMethod()) {
871 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
872 +      if (myMethod == "NONE") {
873 +        esm = NONE;
874 +      } else {
875 +        if (myMethod == "UNDAMPED_WOLF") {
876 +          esm = UNDAMPED_WOLF;
877 +        } else {
878 +          if (myMethod == "DAMPED_WOLF") {            
879 +            esm = DAMPED_WOLF;
880 +            if (!simParams_->haveDampingAlpha()) {
881 +              //throw error
882 +              sprintf( painCave.errMsg,
883 +                       "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.", simParams_->getDampingAlpha());
884 +              painCave.isFatal = 0;
885 +              simError();
886 +            }
887 +            alphaVal = simParams_->getDampingAlpha();
888 +          } else {
889 +            if (myMethod == "REACTION_FIELD") {
890 +              esm = REACTION_FIELD;
891 +            } else {
892 +              // throw error        
893 +              sprintf( painCave.errMsg,
894 +                       "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() );
895 +              painCave.isFatal = 1;
896 +              simError();
897 +            }    
898 +          }          
899 +        }
900 +      }
901 +    }
902 +    initFortranFF( &fInfo_.SIM_uses_RF, &esm, &alphaVal, &errorOut );
903 +  }
904 +
905    void SimInfo::addProperty(GenericData* genData) {
906      properties_.addProperty(genData);  
907    }
# Line 939 | Line 1009 | namespace oopse {
1009  
1010      return o;
1011    }
1012 +  
1013 +  
1014 +   /*
1015 +   Returns center of mass and center of mass velocity in one function call.
1016 +   */
1017 +  
1018 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1019 +      SimInfo::MoleculeIterator i;
1020 +      Molecule* mol;
1021 +      
1022 +    
1023 +      double totalMass = 0.0;
1024 +    
1025 +
1026 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1027 +         double mass = mol->getMass();
1028 +         totalMass += mass;
1029 +         com += mass * mol->getCom();
1030 +         comVel += mass * mol->getComVel();          
1031 +      }  
1032 +      
1033 + #ifdef IS_MPI
1034 +      double tmpMass = totalMass;
1035 +      Vector3d tmpCom(com);  
1036 +      Vector3d tmpComVel(comVel);
1037 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1038 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1039 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1040 + #endif
1041 +      
1042 +      com /= totalMass;
1043 +      comVel /= totalMass;
1044 +   }        
1045 +  
1046 +   /*
1047 +   Return intertia tensor for entire system and angular momentum Vector.
1048  
1049 +
1050 +       [  Ixx -Ixy  -Ixz ]
1051 +  J =| -Iyx  Iyy  -Iyz |
1052 +       [ -Izx -Iyz   Izz ]
1053 +    */
1054 +
1055 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1056 +      
1057 +
1058 +      double xx = 0.0;
1059 +      double yy = 0.0;
1060 +      double zz = 0.0;
1061 +      double xy = 0.0;
1062 +      double xz = 0.0;
1063 +      double yz = 0.0;
1064 +      Vector3d com(0.0);
1065 +      Vector3d comVel(0.0);
1066 +      
1067 +      getComAll(com, comVel);
1068 +      
1069 +      SimInfo::MoleculeIterator i;
1070 +      Molecule* mol;
1071 +      
1072 +      Vector3d thisq(0.0);
1073 +      Vector3d thisv(0.0);
1074 +
1075 +      double thisMass = 0.0;
1076 +    
1077 +      
1078 +      
1079 +  
1080 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1081 +        
1082 +         thisq = mol->getCom()-com;
1083 +         thisv = mol->getComVel()-comVel;
1084 +         thisMass = mol->getMass();
1085 +         // Compute moment of intertia coefficients.
1086 +         xx += thisq[0]*thisq[0]*thisMass;
1087 +         yy += thisq[1]*thisq[1]*thisMass;
1088 +         zz += thisq[2]*thisq[2]*thisMass;
1089 +        
1090 +         // compute products of intertia
1091 +         xy += thisq[0]*thisq[1]*thisMass;
1092 +         xz += thisq[0]*thisq[2]*thisMass;
1093 +         yz += thisq[1]*thisq[2]*thisMass;
1094 +            
1095 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1096 +            
1097 +      }  
1098 +      
1099 +      
1100 +      inertiaTensor(0,0) = yy + zz;
1101 +      inertiaTensor(0,1) = -xy;
1102 +      inertiaTensor(0,2) = -xz;
1103 +      inertiaTensor(1,0) = -xy;
1104 +      inertiaTensor(1,1) = xx + zz;
1105 +      inertiaTensor(1,2) = -yz;
1106 +      inertiaTensor(2,0) = -xz;
1107 +      inertiaTensor(2,1) = -yz;
1108 +      inertiaTensor(2,2) = xx + yy;
1109 +      
1110 + #ifdef IS_MPI
1111 +      Mat3x3d tmpI(inertiaTensor);
1112 +      Vector3d tmpAngMom;
1113 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1114 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1115 + #endif
1116 +              
1117 +      return;
1118 +   }
1119 +
1120 +   //Returns the angular momentum of the system
1121 +   Vector3d SimInfo::getAngularMomentum(){
1122 +      
1123 +      Vector3d com(0.0);
1124 +      Vector3d comVel(0.0);
1125 +      Vector3d angularMomentum(0.0);
1126 +      
1127 +      getComAll(com,comVel);
1128 +      
1129 +      SimInfo::MoleculeIterator i;
1130 +      Molecule* mol;
1131 +      
1132 +      Vector3d thisr(0.0);
1133 +      Vector3d thisp(0.0);
1134 +      
1135 +      double thisMass;
1136 +      
1137 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1138 +        thisMass = mol->getMass();
1139 +        thisr = mol->getCom()-com;
1140 +        thisp = (mol->getComVel()-comVel)*thisMass;
1141 +        
1142 +        angularMomentum += cross( thisr, thisp );
1143 +        
1144 +      }  
1145 +      
1146 + #ifdef IS_MPI
1147 +      Vector3d tmpAngMom;
1148 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1149 + #endif
1150 +      
1151 +      return angularMomentum;
1152 +   }
1153 +  
1154 +  
1155   }//end namespace oopse
1156  

Diff Legend

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