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 580 by chrisfen, Tue Aug 30 18:23:50 2005 UTC

# Line 462 | Line 462 | namespace oopse {
462      //setup fortran force field
463      /** @deprecate */    
464      int isError = 0;
465 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
465 >    initFortranFF( &fInfo_.SIM_uses_RF, &fInfo_.SIM_uses_UW,
466 >                   &fInfo_.SIM_uses_DW, &isError );
467      if(isError){
468        sprintf( painCave.errMsg,
469                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 511 | Line 512 | namespace oopse {
512      int useDipole = 0;
513      int useGayBerne = 0;
514      int useSticky = 0;
515 +    int useStickyPower = 0;
516      int useShape = 0;
517      int useFLARB = 0; //it is not in AtomType yet
518      int useDirectionalAtom = 0;    
# Line 518 | Line 520 | namespace oopse {
520      //usePBC and useRF are from simParams
521      int usePBC = simParams_->getPBC();
522      int useRF = simParams_->getUseRF();
523 +    int useUW = simParams_->getUseUndampedWolf();
524 +    int useDW = simParams_->getUseDampedWolf();
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 578 | Line 586 | namespace oopse {
586  
587      temp = useRF;
588      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
589 +
590 +    temp = useUW;
591 +    MPI_Allreduce(&temp, &useUW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
592 +
593 +    temp = useDW;
594 +    MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
595      
596   #endif
597  
# Line 588 | Line 602 | namespace oopse {
602      fInfo_.SIM_uses_Charges = useCharge;
603      fInfo_.SIM_uses_Dipoles = useDipole;
604      fInfo_.SIM_uses_Sticky = useSticky;
605 +    fInfo_.SIM_uses_StickyPower = useStickyPower;
606      fInfo_.SIM_uses_GayBerne = useGayBerne;
607      fInfo_.SIM_uses_EAM = useEAM;
608      fInfo_.SIM_uses_Shapes = useShape;
609      fInfo_.SIM_uses_FLARB = useFLARB;
610      fInfo_.SIM_uses_RF = useRF;
611 +    fInfo_.SIM_uses_UW = useUW;
612 +    fInfo_.SIM_uses_DW = useDW;
613  
614      if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
615  
# Line 939 | Line 956 | namespace oopse {
956  
957      return o;
958    }
959 +  
960 +  
961 +   /*
962 +   Returns center of mass and center of mass velocity in one function call.
963 +   */
964 +  
965 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
966 +      SimInfo::MoleculeIterator i;
967 +      Molecule* mol;
968 +      
969 +    
970 +      double totalMass = 0.0;
971 +    
972  
973 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
974 +         double mass = mol->getMass();
975 +         totalMass += mass;
976 +         com += mass * mol->getCom();
977 +         comVel += mass * mol->getComVel();          
978 +      }  
979 +      
980 + #ifdef IS_MPI
981 +      double tmpMass = totalMass;
982 +      Vector3d tmpCom(com);  
983 +      Vector3d tmpComVel(comVel);
984 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
985 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
986 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
987 + #endif
988 +      
989 +      com /= totalMass;
990 +      comVel /= totalMass;
991 +   }        
992 +  
993 +   /*
994 +   Return intertia tensor for entire system and angular momentum Vector.
995 +
996 +
997 +       [  Ixx -Ixy  -Ixz ]
998 +  J =| -Iyx  Iyy  -Iyz |
999 +       [ -Izx -Iyz   Izz ]
1000 +    */
1001 +
1002 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1003 +      
1004 +
1005 +      double xx = 0.0;
1006 +      double yy = 0.0;
1007 +      double zz = 0.0;
1008 +      double xy = 0.0;
1009 +      double xz = 0.0;
1010 +      double yz = 0.0;
1011 +      Vector3d com(0.0);
1012 +      Vector3d comVel(0.0);
1013 +      
1014 +      getComAll(com, comVel);
1015 +      
1016 +      SimInfo::MoleculeIterator i;
1017 +      Molecule* mol;
1018 +      
1019 +      Vector3d thisq(0.0);
1020 +      Vector3d thisv(0.0);
1021 +
1022 +      double thisMass = 0.0;
1023 +    
1024 +      
1025 +      
1026 +  
1027 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1028 +        
1029 +         thisq = mol->getCom()-com;
1030 +         thisv = mol->getComVel()-comVel;
1031 +         thisMass = mol->getMass();
1032 +         // Compute moment of intertia coefficients.
1033 +         xx += thisq[0]*thisq[0]*thisMass;
1034 +         yy += thisq[1]*thisq[1]*thisMass;
1035 +         zz += thisq[2]*thisq[2]*thisMass;
1036 +        
1037 +         // compute products of intertia
1038 +         xy += thisq[0]*thisq[1]*thisMass;
1039 +         xz += thisq[0]*thisq[2]*thisMass;
1040 +         yz += thisq[1]*thisq[2]*thisMass;
1041 +            
1042 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1043 +            
1044 +      }  
1045 +      
1046 +      
1047 +      inertiaTensor(0,0) = yy + zz;
1048 +      inertiaTensor(0,1) = -xy;
1049 +      inertiaTensor(0,2) = -xz;
1050 +      inertiaTensor(1,0) = -xy;
1051 +      inertiaTensor(1,1) = xx + zz;
1052 +      inertiaTensor(1,2) = -yz;
1053 +      inertiaTensor(2,0) = -xz;
1054 +      inertiaTensor(2,1) = -yz;
1055 +      inertiaTensor(2,2) = xx + yy;
1056 +      
1057 + #ifdef IS_MPI
1058 +      Mat3x3d tmpI(inertiaTensor);
1059 +      Vector3d tmpAngMom;
1060 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1061 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1062 + #endif
1063 +              
1064 +      return;
1065 +   }
1066 +
1067 +   //Returns the angular momentum of the system
1068 +   Vector3d SimInfo::getAngularMomentum(){
1069 +      
1070 +      Vector3d com(0.0);
1071 +      Vector3d comVel(0.0);
1072 +      Vector3d angularMomentum(0.0);
1073 +      
1074 +      getComAll(com,comVel);
1075 +      
1076 +      SimInfo::MoleculeIterator i;
1077 +      Molecule* mol;
1078 +      
1079 +      Vector3d thisr(0.0);
1080 +      Vector3d thisp(0.0);
1081 +      
1082 +      double thisMass;
1083 +      
1084 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1085 +        thisMass = mol->getMass();
1086 +        thisr = mol->getCom()-com;
1087 +        thisp = (mol->getComVel()-comVel)*thisMass;
1088 +        
1089 +        angularMomentum += cross( thisr, thisp );
1090 +        
1091 +      }  
1092 +      
1093 + #ifdef IS_MPI
1094 +      Vector3d tmpAngMom;
1095 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1096 + #endif
1097 +      
1098 +      return angularMomentum;
1099 +   }
1100 +  
1101 +  
1102   }//end namespace oopse
1103  

Diff Legend

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