# | 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/DarkSide/fElectrostaticSummationMethod.h" |
57 | + | #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" |
58 | + | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" |
59 | #include "UseTheForce/doForces_interface.h" | |
60 | + | #include "UseTheForce/DarkSide/electrostatic_interface.h" |
61 | #include "UseTheForce/notifyCutoffs_interface.h" | |
62 | + | #include "UseTheForce/DarkSide/switcheroo_interface.h" |
63 | #include "utils/MemoryUtils.hpp" | |
64 | #include "utils/simError.h" | |
65 | #include "selection/SelectionManager.hpp" | |
# | Line 80 | Line 86 | namespace oopse { | |
86 | MoleculeStamp* molStamp; | |
87 | int nMolWithSameStamp; | |
88 | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | |
89 | < | int nGroups = 0; //total cutoff groups defined in meta-data file |
89 | > | int nGroups = 0; //total cutoff groups defined in meta-data file |
90 | CutoffGroupStamp* cgStamp; | |
91 | RigidBodyStamp* rbStamp; | |
92 | int nRigidAtoms = 0; | |
# | Line 105 | Line 111 | namespace oopse { | |
111 | } | |
112 | ||
113 | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; | |
114 | + | |
115 | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; | |
116 | ||
117 | //calculate atoms in rigid bodies | |
# | Line 121 | Line 128 | namespace oopse { | |
128 | ||
129 | } | |
130 | ||
131 | < | //every free atom (atom does not belong to cutoff groups) is a cutoff group |
132 | < | //therefore the total number of cutoff groups in the system is equal to |
133 | < | //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data |
134 | < | //file plus the number of cutoff groups defined in meta-data file |
131 | > | //every free atom (atom does not belong to cutoff groups) is a cutoff |
132 | > | //group therefore the total number of cutoff groups in the system is |
133 | > | //equal to the total number of atoms minus number of atoms belong to |
134 | > | //cutoff group defined in meta-data file plus the number of cutoff |
135 | > | //groups defined in meta-data file |
136 | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | |
137 | ||
138 | < | //every free atom (atom does not belong to rigid bodies) is an integrable object |
139 | < | //therefore the total number of integrable objects in the system is equal to |
140 | < | //the total number of atoms minus number of atoms belong to rigid body defined in meta-data |
141 | < | //file plus the number of rigid bodies defined in meta-data file |
142 | < | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_; |
143 | < | |
138 | > | //every free atom (atom does not belong to rigid bodies) is an |
139 | > | //integrable object therefore the total number of integrable objects |
140 | > | //in the system is equal to the total number of atoms minus number of |
141 | > | //atoms belong to rigid body defined in meta-data file plus the number |
142 | > | //of rigid bodies defined in meta-data file |
143 | > | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms |
144 | > | + nGlobalRigidBodies_; |
145 | > | |
146 | nGlobalMols_ = molStampIds_.size(); | |
147 | ||
148 | #ifdef IS_MPI | |
# | Line 462 | Line 472 | namespace oopse { | |
472 | //setup fortran force field | |
473 | /** @deprecate */ | |
474 | int isError = 0; | |
475 | < | initFortranFF( &fInfo_.SIM_uses_RF , &isError ); |
475 | > | |
476 | > | setupElectrostaticSummationMethod( isError ); |
477 | > | setupSwitchingFunction(); |
478 | > | |
479 | if(isError){ | |
480 | sprintf( painCave.errMsg, | |
481 | "ForceField error: There was an error initializing the forceField in fortran.\n" ); | |
# | Line 506 | Line 519 | namespace oopse { | |
519 | int useLennardJones = 0; | |
520 | int useElectrostatic = 0; | |
521 | int useEAM = 0; | |
522 | + | int useSC = 0; |
523 | int useCharge = 0; | |
524 | int useDirectional = 0; | |
525 | int useDipole = 0; | |
526 | int useGayBerne = 0; | |
527 | int useSticky = 0; | |
528 | + | int useStickyPower = 0; |
529 | int useShape = 0; | |
530 | int useFLARB = 0; //it is not in AtomType yet | |
531 | int useDirectionalAtom = 0; | |
532 | int useElectrostatics = 0; | |
533 | //usePBC and useRF are from simParams | |
534 | < | int usePBC = simParams_->getPBC(); |
535 | < | int useRF = simParams_->getUseRF(); |
534 | > | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); |
535 | > | int useRF; |
536 | > | int useSF; |
537 | > | std::string myMethod; |
538 | ||
539 | + | // set the useRF logical |
540 | + | useRF = 0; |
541 | + | useSF = 0; |
542 | + | |
543 | + | |
544 | + | if (simParams_->haveElectrostaticSummationMethod()) { |
545 | + | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
546 | + | toUpper(myMethod); |
547 | + | if (myMethod == "REACTION_FIELD") { |
548 | + | useRF=1; |
549 | + | } else { |
550 | + | if (myMethod == "SHIFTED_FORCE") { |
551 | + | useSF = 1; |
552 | + | } |
553 | + | } |
554 | + | } |
555 | + | |
556 | //loop over all of the atom types | |
557 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
558 | useLennardJones |= (*i)->isLennardJones(); | |
559 | useElectrostatic |= (*i)->isElectrostatic(); | |
560 | useEAM |= (*i)->isEAM(); | |
561 | + | useSC |= (*i)->isSC(); |
562 | useCharge |= (*i)->isCharge(); | |
563 | useDirectional |= (*i)->isDirectional(); | |
564 | useDipole |= (*i)->isDipole(); | |
565 | useGayBerne |= (*i)->isGayBerne(); | |
566 | useSticky |= (*i)->isSticky(); | |
567 | + | useStickyPower |= (*i)->isStickyPower(); |
568 | useShape |= (*i)->isShape(); | |
569 | } | |
570 | ||
571 | < | if (useSticky || useDipole || useGayBerne || useShape) { |
571 | > | if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) { |
572 | useDirectionalAtom = 1; | |
573 | } | |
574 | ||
# | Line 564 | Line 600 | namespace oopse { | |
600 | temp = useSticky; | |
601 | MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
602 | ||
603 | + | temp = useStickyPower; |
604 | + | MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
605 | + | |
606 | temp = useGayBerne; | |
607 | MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
608 | ||
609 | temp = useEAM; | |
610 | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
611 | ||
612 | + | temp = useSC; |
613 | + | MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
614 | + | |
615 | temp = useShape; | |
616 | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
617 | ||
# | Line 578 | Line 620 | namespace oopse { | |
620 | ||
621 | temp = useRF; | |
622 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
623 | < | |
623 | > | |
624 | > | temp = useSF; |
625 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
626 | > | |
627 | #endif | |
628 | ||
629 | fInfo_.SIM_uses_PBC = usePBC; | |
# | Line 588 | Line 633 | namespace oopse { | |
633 | fInfo_.SIM_uses_Charges = useCharge; | |
634 | fInfo_.SIM_uses_Dipoles = useDipole; | |
635 | fInfo_.SIM_uses_Sticky = useSticky; | |
636 | + | fInfo_.SIM_uses_StickyPower = useStickyPower; |
637 | fInfo_.SIM_uses_GayBerne = useGayBerne; | |
638 | fInfo_.SIM_uses_EAM = useEAM; | |
639 | + | fInfo_.SIM_uses_SC = useSC; |
640 | fInfo_.SIM_uses_Shapes = useShape; | |
641 | fInfo_.SIM_uses_FLARB = useFLARB; | |
642 | fInfo_.SIM_uses_RF = useRF; | |
643 | + | fInfo_.SIM_uses_SF = useSF; |
644 | ||
645 | < | if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { |
646 | < | |
645 | > | if( myMethod == "REACTION_FIELD") { |
646 | > | |
647 | if (simParams_->haveDielectric()) { | |
648 | fInfo_.dielect = simParams_->getDielectric(); | |
649 | } else { | |
# | Line 605 | Line 653 | namespace oopse { | |
653 | "\tsetting a dielectric constant!\n"); | |
654 | painCave.isFatal = 1; | |
655 | simError(); | |
656 | < | } |
609 | < | |
610 | < | } else { |
611 | < | fInfo_.dielect = 0.0; |
656 | > | } |
657 | } | |
658 | ||
659 | } | |
# | Line 644 | Line 689 | namespace oopse { | |
689 | ||
690 | totalMass = cg->getMass(); | |
691 | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | |
692 | < | mfact.push_back(atom->getMass()/totalMass); |
692 | > | // Check for massless groups - set mfact to 1 if true |
693 | > | if (totalMass != 0) |
694 | > | mfact.push_back(atom->getMass()/totalMass); |
695 | > | else |
696 | > | mfact.push_back( 1.0 ); |
697 | } | |
698 | ||
699 | } | |
# | Line 780 | Line 829 | namespace oopse { | |
829 | ||
830 | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | |
831 | ||
832 | < | if (!simParams_->haveRcut()){ |
832 | > | if (!simParams_->haveCutoffRadius()){ |
833 | sprintf(painCave.errMsg, | |
834 | "SimCreator Warning: No value was set for the cutoffRadius.\n" | |
835 | "\tOOPSE will use a default value of 15.0 angstroms" | |
# | Line 789 | Line 838 | namespace oopse { | |
838 | simError(); | |
839 | rcut = 15.0; | |
840 | } else{ | |
841 | < | rcut = simParams_->getRcut(); |
841 | > | rcut = simParams_->getCutoffRadius(); |
842 | } | |
843 | ||
844 | < | if (!simParams_->haveRsw()){ |
844 | > | if (!simParams_->haveSwitchingRadius()){ |
845 | sprintf(painCave.errMsg, | |
846 | "SimCreator Warning: No value was set for switchingRadius.\n" | |
847 | "\tOOPSE will use a default value of\n" | |
848 | < | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
848 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
849 | painCave.isFatal = 0; | |
850 | simError(); | |
851 | < | rsw = 0.95 * rcut; |
851 | > | rsw = 0.85 * rcut; |
852 | } else{ | |
853 | < | rsw = simParams_->getRsw(); |
853 | > | rsw = simParams_->getSwitchingRadius(); |
854 | } | |
855 | ||
856 | } else { | |
857 | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in | |
858 | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used | |
859 | ||
860 | < | if (simParams_->haveRcut()) { |
861 | < | rcut = simParams_->getRcut(); |
860 | > | if (simParams_->haveCutoffRadius()) { |
861 | > | rcut = simParams_->getCutoffRadius(); |
862 | } else { | |
863 | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system | |
864 | rcut = calcMaxCutoffRadius(); | |
865 | } | |
866 | ||
867 | < | if (simParams_->haveRsw()) { |
868 | < | rsw = simParams_->getRsw(); |
867 | > | if (simParams_->haveSwitchingRadius()) { |
868 | > | rsw = simParams_->getSwitchingRadius(); |
869 | } else { | |
870 | rsw = rcut; | |
871 | } | |
# | Line 824 | Line 873 | namespace oopse { | |
873 | } | |
874 | } | |
875 | ||
876 | < | void SimInfo::setupCutoff() { |
876 | > | void SimInfo::setupCutoff() { |
877 | getCutoff(rcut_, rsw_); | |
878 | double rnblist = rcut_ + 1; // skin of neighbor list | |
879 | ||
880 | //Pass these cutoff radius etc. to fortran. This function should be called once and only once | |
881 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist); |
881 | > | |
882 | > | int cp = TRADITIONAL_CUTOFF_POLICY; |
883 | > | if (simParams_->haveCutoffPolicy()) { |
884 | > | std::string myPolicy = simParams_->getCutoffPolicy(); |
885 | > | toUpper(myPolicy); |
886 | > | if (myPolicy == "MIX") { |
887 | > | cp = MIX_CUTOFF_POLICY; |
888 | > | } else { |
889 | > | if (myPolicy == "MAX") { |
890 | > | cp = MAX_CUTOFF_POLICY; |
891 | > | } else { |
892 | > | if (myPolicy == "TRADITIONAL") { |
893 | > | cp = TRADITIONAL_CUTOFF_POLICY; |
894 | > | } else { |
895 | > | // throw error |
896 | > | sprintf( painCave.errMsg, |
897 | > | "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); |
898 | > | painCave.isFatal = 1; |
899 | > | simError(); |
900 | > | } |
901 | > | } |
902 | > | } |
903 | > | } |
904 | > | |
905 | > | |
906 | > | if (simParams_->haveSkinThickness()) { |
907 | > | double skinThickness = simParams_->getSkinThickness(); |
908 | > | } |
909 | > | |
910 | > | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); |
911 | > | // also send cutoff notification to electrostatics |
912 | > | setElectrostaticCutoffRadius(&rcut_, &rsw_); |
913 | } | |
914 | ||
915 | + | void SimInfo::setupElectrostaticSummationMethod( int isError ) { |
916 | + | |
917 | + | int errorOut; |
918 | + | int esm = NONE; |
919 | + | int sm = UNDAMPED; |
920 | + | double alphaVal; |
921 | + | double dielectric; |
922 | + | |
923 | + | errorOut = isError; |
924 | + | alphaVal = simParams_->getDampingAlpha(); |
925 | + | dielectric = simParams_->getDielectric(); |
926 | + | |
927 | + | if (simParams_->haveElectrostaticSummationMethod()) { |
928 | + | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
929 | + | toUpper(myMethod); |
930 | + | if (myMethod == "NONE") { |
931 | + | esm = NONE; |
932 | + | } else { |
933 | + | if (myMethod == "SWITCHING_FUNCTION") { |
934 | + | esm = SWITCHING_FUNCTION; |
935 | + | } else { |
936 | + | if (myMethod == "SHIFTED_POTENTIAL") { |
937 | + | esm = SHIFTED_POTENTIAL; |
938 | + | } else { |
939 | + | if (myMethod == "SHIFTED_FORCE") { |
940 | + | esm = SHIFTED_FORCE; |
941 | + | } else { |
942 | + | if (myMethod == "REACTION_FIELD") { |
943 | + | esm = REACTION_FIELD; |
944 | + | } else { |
945 | + | // throw error |
946 | + | sprintf( painCave.errMsg, |
947 | + | "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() ); |
948 | + | painCave.isFatal = 1; |
949 | + | simError(); |
950 | + | } |
951 | + | } |
952 | + | } |
953 | + | } |
954 | + | } |
955 | + | } |
956 | + | |
957 | + | if (simParams_->haveElectrostaticScreeningMethod()) { |
958 | + | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); |
959 | + | toUpper(myScreen); |
960 | + | if (myScreen == "UNDAMPED") { |
961 | + | sm = UNDAMPED; |
962 | + | } else { |
963 | + | if (myScreen == "DAMPED") { |
964 | + | sm = DAMPED; |
965 | + | if (!simParams_->haveDampingAlpha()) { |
966 | + | //throw error |
967 | + | sprintf( painCave.errMsg, |
968 | + | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal); |
969 | + | painCave.isFatal = 0; |
970 | + | simError(); |
971 | + | } |
972 | + | } else { |
973 | + | // throw error |
974 | + | sprintf( painCave.errMsg, |
975 | + | "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); |
976 | + | painCave.isFatal = 1; |
977 | + | simError(); |
978 | + | } |
979 | + | } |
980 | + | } |
981 | + | |
982 | + | // let's pass some summation method variables to fortran |
983 | + | setElectrostaticSummationMethod( &esm ); |
984 | + | setScreeningMethod( &sm ); |
985 | + | setDampingAlpha( &alphaVal ); |
986 | + | setReactionFieldDielectric( &dielectric ); |
987 | + | initFortranFF( &esm, &errorOut ); |
988 | + | } |
989 | + | |
990 | + | void SimInfo::setupSwitchingFunction() { |
991 | + | int ft = CUBIC; |
992 | + | |
993 | + | if (simParams_->haveSwitchingFunctionType()) { |
994 | + | std::string funcType = simParams_->getSwitchingFunctionType(); |
995 | + | toUpper(funcType); |
996 | + | if (funcType == "CUBIC") { |
997 | + | ft = CUBIC; |
998 | + | } else { |
999 | + | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
1000 | + | ft = FIFTH_ORDER_POLY; |
1001 | + | } else { |
1002 | + | // throw error |
1003 | + | sprintf( painCave.errMsg, |
1004 | + | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); |
1005 | + | painCave.isFatal = 1; |
1006 | + | simError(); |
1007 | + | } |
1008 | + | } |
1009 | + | } |
1010 | + | |
1011 | + | // send switching function notification to switcheroo |
1012 | + | setFunctionType(&ft); |
1013 | + | |
1014 | + | } |
1015 | + | |
1016 | void SimInfo::addProperty(GenericData* genData) { | |
1017 | properties_.addProperty(genData); | |
1018 | } | |
# | Line 939 | Line 1120 | namespace oopse { | |
1120 | ||
1121 | return o; | |
1122 | } | |
1123 | + | |
1124 | + | |
1125 | + | /* |
1126 | + | Returns center of mass and center of mass velocity in one function call. |
1127 | + | */ |
1128 | + | |
1129 | + | void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ |
1130 | + | SimInfo::MoleculeIterator i; |
1131 | + | Molecule* mol; |
1132 | + | |
1133 | + | |
1134 | + | double totalMass = 0.0; |
1135 | + | |
1136 | ||
1137 | + | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1138 | + | double mass = mol->getMass(); |
1139 | + | totalMass += mass; |
1140 | + | com += mass * mol->getCom(); |
1141 | + | comVel += mass * mol->getComVel(); |
1142 | + | } |
1143 | + | |
1144 | + | #ifdef IS_MPI |
1145 | + | double tmpMass = totalMass; |
1146 | + | Vector3d tmpCom(com); |
1147 | + | Vector3d tmpComVel(comVel); |
1148 | + | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1149 | + | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1150 | + | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1151 | + | #endif |
1152 | + | |
1153 | + | com /= totalMass; |
1154 | + | comVel /= totalMass; |
1155 | + | } |
1156 | + | |
1157 | + | /* |
1158 | + | Return intertia tensor for entire system and angular momentum Vector. |
1159 | + | |
1160 | + | |
1161 | + | [ Ixx -Ixy -Ixz ] |
1162 | + | J =| -Iyx Iyy -Iyz | |
1163 | + | [ -Izx -Iyz Izz ] |
1164 | + | */ |
1165 | + | |
1166 | + | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ |
1167 | + | |
1168 | + | |
1169 | + | double xx = 0.0; |
1170 | + | double yy = 0.0; |
1171 | + | double zz = 0.0; |
1172 | + | double xy = 0.0; |
1173 | + | double xz = 0.0; |
1174 | + | double yz = 0.0; |
1175 | + | Vector3d com(0.0); |
1176 | + | Vector3d comVel(0.0); |
1177 | + | |
1178 | + | getComAll(com, comVel); |
1179 | + | |
1180 | + | SimInfo::MoleculeIterator i; |
1181 | + | Molecule* mol; |
1182 | + | |
1183 | + | Vector3d thisq(0.0); |
1184 | + | Vector3d thisv(0.0); |
1185 | + | |
1186 | + | double thisMass = 0.0; |
1187 | + | |
1188 | + | |
1189 | + | |
1190 | + | |
1191 | + | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1192 | + | |
1193 | + | thisq = mol->getCom()-com; |
1194 | + | thisv = mol->getComVel()-comVel; |
1195 | + | thisMass = mol->getMass(); |
1196 | + | // Compute moment of intertia coefficients. |
1197 | + | xx += thisq[0]*thisq[0]*thisMass; |
1198 | + | yy += thisq[1]*thisq[1]*thisMass; |
1199 | + | zz += thisq[2]*thisq[2]*thisMass; |
1200 | + | |
1201 | + | // compute products of intertia |
1202 | + | xy += thisq[0]*thisq[1]*thisMass; |
1203 | + | xz += thisq[0]*thisq[2]*thisMass; |
1204 | + | yz += thisq[1]*thisq[2]*thisMass; |
1205 | + | |
1206 | + | angularMomentum += cross( thisq, thisv ) * thisMass; |
1207 | + | |
1208 | + | } |
1209 | + | |
1210 | + | |
1211 | + | inertiaTensor(0,0) = yy + zz; |
1212 | + | inertiaTensor(0,1) = -xy; |
1213 | + | inertiaTensor(0,2) = -xz; |
1214 | + | inertiaTensor(1,0) = -xy; |
1215 | + | inertiaTensor(1,1) = xx + zz; |
1216 | + | inertiaTensor(1,2) = -yz; |
1217 | + | inertiaTensor(2,0) = -xz; |
1218 | + | inertiaTensor(2,1) = -yz; |
1219 | + | inertiaTensor(2,2) = xx + yy; |
1220 | + | |
1221 | + | #ifdef IS_MPI |
1222 | + | Mat3x3d tmpI(inertiaTensor); |
1223 | + | Vector3d tmpAngMom; |
1224 | + | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1225 | + | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1226 | + | #endif |
1227 | + | |
1228 | + | return; |
1229 | + | } |
1230 | + | |
1231 | + | //Returns the angular momentum of the system |
1232 | + | Vector3d SimInfo::getAngularMomentum(){ |
1233 | + | |
1234 | + | Vector3d com(0.0); |
1235 | + | Vector3d comVel(0.0); |
1236 | + | Vector3d angularMomentum(0.0); |
1237 | + | |
1238 | + | getComAll(com,comVel); |
1239 | + | |
1240 | + | SimInfo::MoleculeIterator i; |
1241 | + | Molecule* mol; |
1242 | + | |
1243 | + | Vector3d thisr(0.0); |
1244 | + | Vector3d thisp(0.0); |
1245 | + | |
1246 | + | double thisMass; |
1247 | + | |
1248 | + | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1249 | + | thisMass = mol->getMass(); |
1250 | + | thisr = mol->getCom()-com; |
1251 | + | thisp = (mol->getComVel()-comVel)*thisMass; |
1252 | + | |
1253 | + | angularMomentum += cross( thisr, thisp ); |
1254 | + | |
1255 | + | } |
1256 | + | |
1257 | + | #ifdef IS_MPI |
1258 | + | Vector3d tmpAngMom; |
1259 | + | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1260 | + | #endif |
1261 | + | |
1262 | + | return angularMomentum; |
1263 | + | } |
1264 | + | |
1265 | + | |
1266 | }//end namespace oopse | |
1267 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |