# | 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/doForces_interface.h" | |
57 | #include "UseTheForce/notifyCutoffs_interface.h" | |
58 | #include "utils/MemoryUtils.hpp" | |
# | Line 462 | Line 463 | namespace oopse { | |
463 | //setup fortran force field | |
464 | /** @deprecate */ | |
465 | int isError = 0; | |
466 | < | initFortranFF( &fInfo_.SIM_uses_RF , &isError ); |
466 | > | initFortranFF( &fInfo_.SIM_uses_RF, &fInfo_.SIM_uses_UW, |
467 | > | &fInfo_.SIM_uses_DW, &isError ); |
468 | if(isError){ | |
469 | sprintf( painCave.errMsg, | |
470 | "ForceField error: There was an error initializing the forceField in fortran.\n" ); | |
# | Line 519 | Line 521 | namespace oopse { | |
521 | //usePBC and useRF are from simParams | |
522 | int usePBC = simParams_->getPBC(); | |
523 | int useRF = simParams_->getUseRF(); | |
524 | + | int useUW = simParams_->getUseUndampedWolf(); |
525 | + | int useDW = simParams_->getUseDampedWolf(); |
526 | ||
527 | //loop over all of the atom types | |
528 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
# | Line 583 | Line 587 | namespace oopse { | |
587 | ||
588 | temp = useRF; | |
589 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
590 | + | |
591 | + | temp = useUW; |
592 | + | MPI_Allreduce(&temp, &useUW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
593 | + | |
594 | + | temp = useDW; |
595 | + | MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
596 | ||
597 | #endif | |
598 | ||
# | Line 599 | Line 609 | namespace oopse { | |
609 | fInfo_.SIM_uses_Shapes = useShape; | |
610 | fInfo_.SIM_uses_FLARB = useFLARB; | |
611 | fInfo_.SIM_uses_RF = useRF; | |
612 | + | fInfo_.SIM_uses_UW = useUW; |
613 | + | fInfo_.SIM_uses_DW = useDW; |
614 | ||
615 | if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { | |
616 | ||
# | Line 830 | Line 842 | namespace oopse { | |
842 | } | |
843 | } | |
844 | ||
845 | < | void SimInfo::setupCutoff() { |
845 | > | void SimInfo::setupCutoff() { |
846 | getCutoff(rcut_, rsw_); | |
847 | double rnblist = rcut_ + 1; // skin of neighbor list | |
848 | ||
849 | //Pass these cutoff radius etc. to fortran. This function should be called once and only once | |
850 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist); |
850 | > | |
851 | > | int cp = TRADITIONAL_CUTOFF_POLICY; |
852 | > | if (simParams_->haveCutoffPolicy()) { |
853 | > | std::string myPolicy = simParams_->getCutoffPolicy(); |
854 | > | if (myPolicy == "MIX") { |
855 | > | cp = MIX_CUTOFF_POLICY; |
856 | > | } else { |
857 | > | if (myPolicy == "MAX") { |
858 | > | cp = MAX_CUTOFF_POLICY; |
859 | > | } else { |
860 | > | if (myPolicy == "TRADITIONAL") { |
861 | > | cp = TRADITIONAL_CUTOFF_POLICY; |
862 | > | } else { |
863 | > | // throw error |
864 | > | sprintf( painCave.errMsg, |
865 | > | "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); |
866 | > | painCave.isFatal = 1; |
867 | > | simError(); |
868 | > | } |
869 | > | } |
870 | > | } |
871 | > | } |
872 | > | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); |
873 | } | |
874 | ||
875 | void SimInfo::addProperty(GenericData* genData) { | |
# | Line 945 | Line 979 | namespace oopse { | |
979 | ||
980 | return o; | |
981 | } | |
982 | + | |
983 | + | |
984 | + | /* |
985 | + | Returns center of mass and center of mass velocity in one function call. |
986 | + | */ |
987 | + | |
988 | + | void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ |
989 | + | SimInfo::MoleculeIterator i; |
990 | + | Molecule* mol; |
991 | + | |
992 | + | |
993 | + | double totalMass = 0.0; |
994 | + | |
995 | + | |
996 | + | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
997 | + | double mass = mol->getMass(); |
998 | + | totalMass += mass; |
999 | + | com += mass * mol->getCom(); |
1000 | + | comVel += mass * mol->getComVel(); |
1001 | + | } |
1002 | + | |
1003 | + | #ifdef IS_MPI |
1004 | + | double tmpMass = totalMass; |
1005 | + | Vector3d tmpCom(com); |
1006 | + | Vector3d tmpComVel(comVel); |
1007 | + | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1008 | + | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1009 | + | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1010 | + | #endif |
1011 | + | |
1012 | + | com /= totalMass; |
1013 | + | comVel /= totalMass; |
1014 | + | } |
1015 | + | |
1016 | + | /* |
1017 | + | Return intertia tensor for entire system and angular momentum Vector. |
1018 | + | |
1019 | + | |
1020 | + | [ Ixx -Ixy -Ixz ] |
1021 | + | J =| -Iyx Iyy -Iyz | |
1022 | + | [ -Izx -Iyz Izz ] |
1023 | + | */ |
1024 | + | |
1025 | + | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ |
1026 | + | |
1027 | + | |
1028 | + | double xx = 0.0; |
1029 | + | double yy = 0.0; |
1030 | + | double zz = 0.0; |
1031 | + | double xy = 0.0; |
1032 | + | double xz = 0.0; |
1033 | + | double yz = 0.0; |
1034 | + | Vector3d com(0.0); |
1035 | + | Vector3d comVel(0.0); |
1036 | + | |
1037 | + | getComAll(com, comVel); |
1038 | + | |
1039 | + | SimInfo::MoleculeIterator i; |
1040 | + | Molecule* mol; |
1041 | + | |
1042 | + | Vector3d thisq(0.0); |
1043 | + | Vector3d thisv(0.0); |
1044 | ||
1045 | + | double thisMass = 0.0; |
1046 | + | |
1047 | + | |
1048 | + | |
1049 | + | |
1050 | + | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1051 | + | |
1052 | + | thisq = mol->getCom()-com; |
1053 | + | thisv = mol->getComVel()-comVel; |
1054 | + | thisMass = mol->getMass(); |
1055 | + | // Compute moment of intertia coefficients. |
1056 | + | xx += thisq[0]*thisq[0]*thisMass; |
1057 | + | yy += thisq[1]*thisq[1]*thisMass; |
1058 | + | zz += thisq[2]*thisq[2]*thisMass; |
1059 | + | |
1060 | + | // compute products of intertia |
1061 | + | xy += thisq[0]*thisq[1]*thisMass; |
1062 | + | xz += thisq[0]*thisq[2]*thisMass; |
1063 | + | yz += thisq[1]*thisq[2]*thisMass; |
1064 | + | |
1065 | + | angularMomentum += cross( thisq, thisv ) * thisMass; |
1066 | + | |
1067 | + | } |
1068 | + | |
1069 | + | |
1070 | + | inertiaTensor(0,0) = yy + zz; |
1071 | + | inertiaTensor(0,1) = -xy; |
1072 | + | inertiaTensor(0,2) = -xz; |
1073 | + | inertiaTensor(1,0) = -xy; |
1074 | + | inertiaTensor(1,1) = xx + zz; |
1075 | + | inertiaTensor(1,2) = -yz; |
1076 | + | inertiaTensor(2,0) = -xz; |
1077 | + | inertiaTensor(2,1) = -yz; |
1078 | + | inertiaTensor(2,2) = xx + yy; |
1079 | + | |
1080 | + | #ifdef IS_MPI |
1081 | + | Mat3x3d tmpI(inertiaTensor); |
1082 | + | Vector3d tmpAngMom; |
1083 | + | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1084 | + | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1085 | + | #endif |
1086 | + | |
1087 | + | return; |
1088 | + | } |
1089 | + | |
1090 | + | //Returns the angular momentum of the system |
1091 | + | Vector3d SimInfo::getAngularMomentum(){ |
1092 | + | |
1093 | + | Vector3d com(0.0); |
1094 | + | Vector3d comVel(0.0); |
1095 | + | Vector3d angularMomentum(0.0); |
1096 | + | |
1097 | + | getComAll(com,comVel); |
1098 | + | |
1099 | + | SimInfo::MoleculeIterator i; |
1100 | + | Molecule* mol; |
1101 | + | |
1102 | + | Vector3d thisr(0.0); |
1103 | + | Vector3d thisp(0.0); |
1104 | + | |
1105 | + | double thisMass; |
1106 | + | |
1107 | + | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1108 | + | thisMass = mol->getMass(); |
1109 | + | thisr = mol->getCom()-com; |
1110 | + | thisp = (mol->getComVel()-comVel)*thisMass; |
1111 | + | |
1112 | + | angularMomentum += cross( thisr, thisp ); |
1113 | + | |
1114 | + | } |
1115 | + | |
1116 | + | #ifdef IS_MPI |
1117 | + | Vector3d tmpAngMom; |
1118 | + | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1119 | + | #endif |
1120 | + | |
1121 | + | return angularMomentum; |
1122 | + | } |
1123 | + | |
1124 | + | |
1125 | }//end namespace oopse | |
1126 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |