| 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;     | 
| 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 | 
  | 
 | 
| 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 | 
  | 
 | 
| 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; | 
| 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 | 
  | 
 |