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 555 by chuckv, Mon May 30 14:01:52 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 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
987 +      
988 +
989 +      double xx = 0.0;
990 +      double yy = 0.0;
991 +      double zz = 0.0;
992 +      double xy = 0.0;
993 +      double xz = 0.0;
994 +      double yz = 0.0;
995 +      Vector3d com(0.0);
996 +      Vector3d comVel(0.0);
997 +      
998 +      getComAll(com, comVel);
999 +      
1000 +      SimInfo::MoleculeIterator i;
1001 +      Molecule* mol;
1002 +      
1003 +      Vector3d thisq(0.0);
1004 +      Vector3d thisv(0.0);
1005 +
1006 +      double thisMass = 0.0;
1007 +    
1008 +      
1009 +      
1010 +  
1011 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1012 +        
1013 +         thisq = mol->getCom()-com;
1014 +         thisv = mol->getComVel()-comVel;
1015 +         thisMass = mol->getMass();
1016 +         // Compute moment of intertia coefficients.
1017 +         xx += thisq[0]*thisq[0]*thisMass;
1018 +         yy += thisq[1]*thisq[1]*thisMass;
1019 +         zz += thisq[2]*thisq[2]*thisMass;
1020 +        
1021 +         // compute products of intertia
1022 +         xy += thisq[0]*thisq[1]*thisMass;
1023 +         xz += thisq[0]*thisq[2]*thisMass;
1024 +         yz += thisq[1]*thisq[2]*thisMass;
1025 +            
1026 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1027 +            
1028 +      }  
1029 +      
1030 +      
1031 +      inertiaTensor(0,0) = yy + zz;
1032 +      inertiaTensor(0,1) = -xy;
1033 +      inertiaTensor(0,2) = -xz;
1034 +      inertiaTensor(1,0) = -xy;
1035 +      inertiaTensor(2,0) = xx + zz;
1036 +      inertiaTensor(1,2) = -yz;
1037 +      inertiaTensor(2,0) = -xz;
1038 +      inertiaTensor(2,1) = -yz;
1039 +      inertiaTensor(2,2) = xx + yy;
1040 +      
1041 + #ifdef IS_MPI
1042 +      Mat3x3d tmpI(inertiaTensor);
1043 +      Vector3d tmpAngMom;
1044 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1045 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1046 + #endif
1047 +              
1048 +      return;
1049 +   }
1050  
1051 +   //Returns the angular momentum of the system
1052 +   Vector3d SimInfo::getAngularMomentum(){
1053 +      
1054 +      Vector3d com(0.0);
1055 +      Vector3d comVel(0.0);
1056 +      Vector3d angularMomentum(0.0);
1057 +      
1058 +      getComAll(com,comVel);
1059 +      
1060 +      SimInfo::MoleculeIterator i;
1061 +      Molecule* mol;
1062 +      
1063 +      Vector3d thisq(0.0);
1064 +      Vector3d thisv(0.0);
1065 +      
1066 +      double thisMass = 0.0;
1067 +      
1068 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1069 +         thisq = mol->getCom()-com;
1070 +         thisv = mol->getComVel()-comVel;
1071 +         thisMass = mol->getMass();
1072 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1073 +        
1074 +      }  
1075 +      
1076 + #ifdef IS_MPI
1077 +      Vector3d tmpAngMom;
1078 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1079 + #endif
1080 +      
1081 +      return angularMomentum;
1082 +   }
1083 +  
1084 +  
1085   }//end namespace oopse
1086  

Diff Legend

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