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 586 by gezelter, Wed Sep 7 20:46:46 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/doForces_interface.h"
57   #include "UseTheForce/notifyCutoffs_interface.h"
58   #include "utils/MemoryUtils.hpp"
# Line 462 | Line 463 | namespace oopse {
463      //setup fortran force field
464      /** @deprecate */    
465      int isError = 0;
466 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
466 >    initFortranFF( &fInfo_.SIM_uses_RF, &fInfo_.SIM_uses_UW,
467 >                   &fInfo_.SIM_uses_DW, &isError );
468      if(isError){
469        sprintf( painCave.errMsg,
470                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 511 | Line 513 | namespace oopse {
513      int useDipole = 0;
514      int useGayBerne = 0;
515      int useSticky = 0;
516 +    int useStickyPower = 0;
517      int useShape = 0;
518      int useFLARB = 0; //it is not in AtomType yet
519      int useDirectionalAtom = 0;    
# Line 518 | Line 521 | namespace oopse {
521      //usePBC and useRF are from simParams
522      int usePBC = simParams_->getPBC();
523      int useRF = simParams_->getUseRF();
524 +    int useUW = simParams_->getUseUndampedWolf();
525 +    int useDW = simParams_->getUseDampedWolf();
526  
527      //loop over all of the atom types
528      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
# Line 529 | Line 534 | namespace oopse {
534        useDipole |= (*i)->isDipole();
535        useGayBerne |= (*i)->isGayBerne();
536        useSticky |= (*i)->isSticky();
537 +      useStickyPower |= (*i)->isStickyPower();
538        useShape |= (*i)->isShape();
539      }
540  
541 <    if (useSticky || useDipole || useGayBerne || useShape) {
541 >    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
542        useDirectionalAtom = 1;
543      }
544  
# Line 564 | Line 570 | namespace oopse {
570      temp = useSticky;
571      MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
572  
573 +    temp = useStickyPower;
574 +    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
575 +    
576      temp = useGayBerne;
577      MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
578  
# Line 578 | Line 587 | namespace oopse {
587  
588      temp = useRF;
589      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
590 +
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      
597   #endif
598  
# Line 588 | Line 603 | namespace oopse {
603      fInfo_.SIM_uses_Charges = useCharge;
604      fInfo_.SIM_uses_Dipoles = useDipole;
605      fInfo_.SIM_uses_Sticky = useSticky;
606 +    fInfo_.SIM_uses_StickyPower = useStickyPower;
607      fInfo_.SIM_uses_GayBerne = useGayBerne;
608      fInfo_.SIM_uses_EAM = useEAM;
609      fInfo_.SIM_uses_Shapes = useShape;
610      fInfo_.SIM_uses_FLARB = useFLARB;
611      fInfo_.SIM_uses_RF = useRF;
612 +    fInfo_.SIM_uses_UW = useUW;
613 +    fInfo_.SIM_uses_DW = useDW;
614  
615      if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
616  
# Line 824 | Line 842 | namespace oopse {
842      }
843    }
844  
845 <  void SimInfo::setupCutoff() {
845 >  void SimInfo::setupCutoff() {    
846      getCutoff(rcut_, rsw_);    
847      double rnblist = rcut_ + 1; // skin of neighbor list
848  
849      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
850 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
850 >    
851 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
852 >    if (simParams_->haveCutoffPolicy()) {
853 >      std::string myPolicy = simParams_->getCutoffPolicy();
854 >      if (myPolicy == "MIX") {
855 >        cp = MIX_CUTOFF_POLICY;
856 >      } else {
857 >        if (myPolicy == "MAX") {
858 >          cp = MAX_CUTOFF_POLICY;
859 >        } else {
860 >          if (myPolicy == "TRADITIONAL") {            
861 >            cp = TRADITIONAL_CUTOFF_POLICY;
862 >          } else {
863 >            // throw error        
864 >            sprintf( painCave.errMsg,
865 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
866 >            painCave.isFatal = 1;
867 >            simError();
868 >          }    
869 >        }          
870 >      }
871 >    }
872 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
873    }
874  
875    void SimInfo::addProperty(GenericData* genData) {
# Line 939 | Line 979 | namespace oopse {
979  
980      return o;
981    }
982 +  
983 +  
984 +   /*
985 +   Returns center of mass and center of mass velocity in one function call.
986 +   */
987 +  
988 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
989 +      SimInfo::MoleculeIterator i;
990 +      Molecule* mol;
991 +      
992 +    
993 +      double totalMass = 0.0;
994 +    
995 +
996 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
997 +         double mass = mol->getMass();
998 +         totalMass += mass;
999 +         com += mass * mol->getCom();
1000 +         comVel += mass * mol->getComVel();          
1001 +      }  
1002 +      
1003 + #ifdef IS_MPI
1004 +      double tmpMass = totalMass;
1005 +      Vector3d tmpCom(com);  
1006 +      Vector3d tmpComVel(comVel);
1007 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1008 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1009 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1010 + #endif
1011 +      
1012 +      com /= totalMass;
1013 +      comVel /= totalMass;
1014 +   }        
1015 +  
1016 +   /*
1017 +   Return intertia tensor for entire system and angular momentum Vector.
1018 +
1019 +
1020 +       [  Ixx -Ixy  -Ixz ]
1021 +  J =| -Iyx  Iyy  -Iyz |
1022 +       [ -Izx -Iyz   Izz ]
1023 +    */
1024 +
1025 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1026 +      
1027 +
1028 +      double xx = 0.0;
1029 +      double yy = 0.0;
1030 +      double zz = 0.0;
1031 +      double xy = 0.0;
1032 +      double xz = 0.0;
1033 +      double yz = 0.0;
1034 +      Vector3d com(0.0);
1035 +      Vector3d comVel(0.0);
1036 +      
1037 +      getComAll(com, comVel);
1038 +      
1039 +      SimInfo::MoleculeIterator i;
1040 +      Molecule* mol;
1041 +      
1042 +      Vector3d thisq(0.0);
1043 +      Vector3d thisv(0.0);
1044  
1045 +      double thisMass = 0.0;
1046 +    
1047 +      
1048 +      
1049 +  
1050 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1051 +        
1052 +         thisq = mol->getCom()-com;
1053 +         thisv = mol->getComVel()-comVel;
1054 +         thisMass = mol->getMass();
1055 +         // Compute moment of intertia coefficients.
1056 +         xx += thisq[0]*thisq[0]*thisMass;
1057 +         yy += thisq[1]*thisq[1]*thisMass;
1058 +         zz += thisq[2]*thisq[2]*thisMass;
1059 +        
1060 +         // compute products of intertia
1061 +         xy += thisq[0]*thisq[1]*thisMass;
1062 +         xz += thisq[0]*thisq[2]*thisMass;
1063 +         yz += thisq[1]*thisq[2]*thisMass;
1064 +            
1065 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1066 +            
1067 +      }  
1068 +      
1069 +      
1070 +      inertiaTensor(0,0) = yy + zz;
1071 +      inertiaTensor(0,1) = -xy;
1072 +      inertiaTensor(0,2) = -xz;
1073 +      inertiaTensor(1,0) = -xy;
1074 +      inertiaTensor(1,1) = xx + zz;
1075 +      inertiaTensor(1,2) = -yz;
1076 +      inertiaTensor(2,0) = -xz;
1077 +      inertiaTensor(2,1) = -yz;
1078 +      inertiaTensor(2,2) = xx + yy;
1079 +      
1080 + #ifdef IS_MPI
1081 +      Mat3x3d tmpI(inertiaTensor);
1082 +      Vector3d tmpAngMom;
1083 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1084 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1085 + #endif
1086 +              
1087 +      return;
1088 +   }
1089 +
1090 +   //Returns the angular momentum of the system
1091 +   Vector3d SimInfo::getAngularMomentum(){
1092 +      
1093 +      Vector3d com(0.0);
1094 +      Vector3d comVel(0.0);
1095 +      Vector3d angularMomentum(0.0);
1096 +      
1097 +      getComAll(com,comVel);
1098 +      
1099 +      SimInfo::MoleculeIterator i;
1100 +      Molecule* mol;
1101 +      
1102 +      Vector3d thisr(0.0);
1103 +      Vector3d thisp(0.0);
1104 +      
1105 +      double thisMass;
1106 +      
1107 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1108 +        thisMass = mol->getMass();
1109 +        thisr = mol->getCom()-com;
1110 +        thisp = (mol->getComVel()-comVel)*thisMass;
1111 +        
1112 +        angularMomentum += cross( thisr, thisp );
1113 +        
1114 +      }  
1115 +      
1116 + #ifdef IS_MPI
1117 +      Vector3d tmpAngMom;
1118 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1119 + #endif
1120 +      
1121 +      return angularMomentum;
1122 +   }
1123 +  
1124 +  
1125   }//end namespace oopse
1126  

Diff Legend

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