| 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" | 
| 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; | 
| 111 |  | } | 
| 112 |  |  | 
| 113 |  | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; | 
| 114 | + |  | 
| 115 |  | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; | 
| 116 |  |  | 
| 117 |  | //calculate atoms in rigid bodies | 
| 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 | 
| 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" ); | 
| 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 |  |  | 
| 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 |  |  | 
| 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; | 
| 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 { | 
| 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 |  | } | 
| 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 |  | } | 
| 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" | 
| 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 |  | } | 
| 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 |  | } | 
| 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 |  |  |