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 557 by chuckv, Tue May 31 22:31:54 2005 UTC

# Line 511 | Line 511 | namespace oopse {
511      int useDipole = 0;
512      int useGayBerne = 0;
513      int useSticky = 0;
514 +    int useStickyPower = 0;
515      int useShape = 0;
516      int useFLARB = 0; //it is not in AtomType yet
517      int useDirectionalAtom = 0;    
# Line 529 | Line 530 | namespace oopse {
530        useDipole |= (*i)->isDipole();
531        useGayBerne |= (*i)->isGayBerne();
532        useSticky |= (*i)->isSticky();
533 +      useStickyPower |= (*i)->isStickyPower();
534        useShape |= (*i)->isShape();
535      }
536  
537 <    if (useSticky || useDipole || useGayBerne || useShape) {
537 >    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
538        useDirectionalAtom = 1;
539      }
540  
# Line 564 | Line 566 | namespace oopse {
566      temp = useSticky;
567      MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
568  
569 +    temp = useStickyPower;
570 +    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
571 +    
572      temp = useGayBerne;
573      MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
574  
# 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;
# Line 939 | Line 945 | namespace oopse {
945  
946      return o;
947    }
948 +  
949 +  
950 +   /*
951 +   Returns center of mass and center of mass velocity in one function call.
952 +   */
953 +  
954 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
955 +      SimInfo::MoleculeIterator i;
956 +      Molecule* mol;
957 +      
958 +    
959 +      double totalMass = 0.0;
960 +    
961 +
962 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
963 +         double mass = mol->getMass();
964 +         totalMass += mass;
965 +         com += mass * mol->getCom();
966 +         comVel += mass * mol->getComVel();          
967 +      }  
968 +      
969 + #ifdef IS_MPI
970 +      double tmpMass = totalMass;
971 +      Vector3d tmpCom(com);  
972 +      Vector3d tmpComVel(comVel);
973 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
974 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
975 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
976 + #endif
977 +      
978 +      com /= totalMass;
979 +      comVel /= totalMass;
980 +   }        
981 +  
982 +   /*
983 +   Return intertia tensor for entire system and angular momentum Vector.
984 +
985 +
986 +       [  Ixx -Ixy  -Ixz ]
987 +  J =| -Iyx  Iyy  -Iyz |
988 +       [ -Izx -Iyz   Izz ]
989 +    */
990 +
991 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
992 +      
993 +
994 +      double xx = 0.0;
995 +      double yy = 0.0;
996 +      double zz = 0.0;
997 +      double xy = 0.0;
998 +      double xz = 0.0;
999 +      double yz = 0.0;
1000 +      Vector3d com(0.0);
1001 +      Vector3d comVel(0.0);
1002 +      
1003 +      getComAll(com, comVel);
1004 +      
1005 +      SimInfo::MoleculeIterator i;
1006 +      Molecule* mol;
1007 +      
1008 +      Vector3d thisq(0.0);
1009 +      Vector3d thisv(0.0);
1010 +
1011 +      double thisMass = 0.0;
1012 +    
1013 +      
1014 +      
1015 +  
1016 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1017 +        
1018 +         thisq = mol->getCom()-com;
1019 +         thisv = mol->getComVel()-comVel;
1020 +         thisMass = mol->getMass();
1021 +         // Compute moment of intertia coefficients.
1022 +         xx += thisq[0]*thisq[0]*thisMass;
1023 +         yy += thisq[1]*thisq[1]*thisMass;
1024 +         zz += thisq[2]*thisq[2]*thisMass;
1025 +        
1026 +         // compute products of intertia
1027 +         xy += thisq[0]*thisq[1]*thisMass;
1028 +         xz += thisq[0]*thisq[2]*thisMass;
1029 +         yz += thisq[1]*thisq[2]*thisMass;
1030 +            
1031 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1032 +            
1033 +      }  
1034 +      
1035 +      
1036 +      inertiaTensor(0,0) = yy + zz;
1037 +      inertiaTensor(0,1) = -xy;
1038 +      inertiaTensor(0,2) = -xz;
1039 +      inertiaTensor(1,0) = -xy;
1040 +      inertiaTensor(1,1) = xx + zz;
1041 +      inertiaTensor(1,2) = -yz;
1042 +      inertiaTensor(2,0) = -xz;
1043 +      inertiaTensor(2,1) = -yz;
1044 +      inertiaTensor(2,2) = xx + yy;
1045 +      
1046 + #ifdef IS_MPI
1047 +      Mat3x3d tmpI(inertiaTensor);
1048 +      Vector3d tmpAngMom;
1049 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1050 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1051 + #endif
1052 +              
1053 +      return;
1054 +   }
1055  
1056 +   //Returns the angular momentum of the system
1057 +   Vector3d SimInfo::getAngularMomentum(){
1058 +      
1059 +      Vector3d com(0.0);
1060 +      Vector3d comVel(0.0);
1061 +      Vector3d angularMomentum(0.0);
1062 +      
1063 +      getComAll(com,comVel);
1064 +      
1065 +      SimInfo::MoleculeIterator i;
1066 +      Molecule* mol;
1067 +      
1068 +      Vector3d thisr(0.0);
1069 +      Vector3d thisp(0.0);
1070 +      
1071 +      double thisMass;
1072 +      
1073 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1074 +        thisMass = mol->getMass();
1075 +        thisr = mol->getCom()-com;
1076 +        thisp = (mol->getComVel()-comVel)*thisMass;
1077 +        
1078 +        angularMomentum += cross( thisr, thisp );
1079 +        
1080 +      }  
1081 +      
1082 + #ifdef IS_MPI
1083 +      Vector3d tmpAngMom;
1084 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1085 + #endif
1086 +      
1087 +      return angularMomentum;
1088 +   }
1089 +  
1090 +  
1091   }//end namespace oopse
1092  

Diff Legend

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