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 764 by gezelter, Mon Nov 21 22:59:21 2005 UTC vs.
Revision 1129 by chrisfen, Fri Apr 20 18:15:48 2007 UTC

# Line 53 | Line 53
53   #include "brains/SimInfo.hpp"
54   #include "math/Vector3.hpp"
55   #include "primitives/Molecule.hpp"
56 + #include "primitives/StuntDouble.hpp"
57   #include "UseTheForce/fCutoffPolicy.h"
58   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
59   #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
60   #include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
61   #include "UseTheForce/doForces_interface.h"
62 + #include "UseTheForce/DarkSide/neighborLists_interface.h"
63   #include "UseTheForce/DarkSide/electrostatic_interface.h"
64   #include "UseTheForce/DarkSide/switcheroo_interface.h"
65   #include "utils/MemoryUtils.hpp"
66   #include "utils/simError.h"
67   #include "selection/SelectionManager.hpp"
68 + #include "io/ForceFieldOptions.hpp"
69 + #include "UseTheForce/ForceField.hpp"
70  
71 +
72   #ifdef IS_MPI
73   #include "UseTheForce/mpiComponentPlan.h"
74   #include "UseTheForce/DarkSide/simParallel_interface.h"
# Line 80 | Line 85 | namespace oopse {
85      return result;
86    }
87    
88 <  SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
89 <                   ForceField* ff, Globals* simParams) :
90 <    stamps_(stamps), forceField_(ff), simParams_(simParams),
86 <    ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
88 >  SimInfo::SimInfo(ForceField* ff, Globals* simParams) :
89 >    forceField_(ff), simParams_(simParams),
90 >    ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
91      nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
92      nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
93      nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
94      nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
95 <    sman_(NULL), fortranInitialized_(false) {
95 >    sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false),
96 >    useAtomicVirial_(true) {
97  
93 –            
94 –      std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
98        MoleculeStamp* molStamp;
99        int nMolWithSameStamp;
100        int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
# Line 99 | Line 102 | namespace oopse {
102        CutoffGroupStamp* cgStamp;    
103        RigidBodyStamp* rbStamp;
104        int nRigidAtoms = 0;
105 <    
106 <      for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
107 <        molStamp = i->first;
108 <        nMolWithSameStamp = i->second;
105 >      std::vector<Component*> components = simParams->getComponents();
106 >      
107 >      for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
108 >        molStamp = (*i)->getMoleculeStamp();
109 >        nMolWithSameStamp = (*i)->getNMol();
110          
111          addMoleculeStamp(molStamp, nMolWithSameStamp);
112  
113          //calculate atoms in molecules
114          nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
115  
112 –
116          //calculate atoms in cutoff groups
117          int nAtomsInGroups = 0;
118          int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
119          
120          for (int j=0; j < nCutoffGroupsInStamp; j++) {
121 <          cgStamp = molStamp->getCutoffGroup(j);
121 >          cgStamp = molStamp->getCutoffGroupStamp(j);
122            nAtomsInGroups += cgStamp->getNMembers();
123          }
124  
# Line 128 | Line 131 | namespace oopse {
131          int nRigidBodiesInStamp = molStamp->getNRigidBodies();
132          
133          for (int j=0; j < nRigidBodiesInStamp; j++) {
134 <          rbStamp = molStamp->getRigidBody(j);
134 >          rbStamp = molStamp->getRigidBodyStamp(j);
135            nAtomsInRigidBodies += rbStamp->getNMembers();
136          }
137  
# Line 167 | Line 170 | namespace oopse {
170      }
171      molecules_.clear();
172        
170 –    delete stamps_;
173      delete sman_;
174      delete simParams_;
175      delete forceField_;
# Line 274 | Line 276 | namespace oopse {
276            }
277          }
278              
279 <      }//end for (integrableObject)
280 <    }// end for (mol)
279 >      }
280 >    }
281      
282      // n_constraints is local, so subtract them on each processor
283      ndf_local -= nConstraints_;
# Line 292 | Line 294 | namespace oopse {
294  
295    }
296  
297 +  int SimInfo::getFdf() {
298 + #ifdef IS_MPI
299 +    MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
300 + #else
301 +    fdf_ = fdf_local;
302 + #endif
303 +    return fdf_;
304 +  }
305 +    
306    void SimInfo::calcNdfRaw() {
307      int ndfRaw_local;
308  
# Line 593 | Line 604 | namespace oopse {
604      /** @deprecate */    
605      int isError = 0;
606      
607 +    setupCutoff();
608 +    
609      setupElectrostaticSummationMethod( isError );
610      setupSwitchingFunction();
611 +    setupAccumulateBoxDipole();
612  
613      if(isError){
614        sprintf( painCave.errMsg,
# Line 602 | Line 616 | namespace oopse {
616        painCave.isFatal = 1;
617        simError();
618      }
605 –  
606 –    
607 –    setupCutoff();
619  
620      calcNdf();
621      calcNdfRaw();
# Line 654 | Line 665 | namespace oopse {
665      int usePBC = simParams_->getUsePeriodicBoundaryConditions();
666      int useRF;
667      int useSF;
668 +    int useSP;
669 +    int useBoxDipole;
670 +
671      std::string myMethod;
672  
673      // set the useRF logical
674      useRF = 0;
675      useSF = 0;
676 +    useSP = 0;
677  
678  
679      if (simParams_->haveElectrostaticSummationMethod()) {
680        std::string myMethod = simParams_->getElectrostaticSummationMethod();
681        toUpper(myMethod);
682 <      if (myMethod == "REACTION_FIELD") {
683 <        useRF=1;
684 <      } else {
685 <        if (myMethod == "SHIFTED_FORCE") {
686 <          useSF = 1;
687 <        }
682 >      if (myMethod == "REACTION_FIELD"){
683 >        useRF = 1;
684 >      } else if (myMethod == "SHIFTED_FORCE"){
685 >        useSF = 1;
686 >      } else if (myMethod == "SHIFTED_POTENTIAL"){
687 >        useSP = 1;
688        }
689      }
690 +    
691 +    if (simParams_->haveAccumulateBoxDipole())
692 +      if (simParams_->getAccumulateBoxDipole())
693 +        useBoxDipole = 1;
694  
695 +    useAtomicVirial_ = simParams_->getUseAtomicVirial();
696 +
697      //loop over all of the atom types
698      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
699        useLennardJones |= (*i)->isLennardJones();
# Line 742 | Line 763 | namespace oopse {
763      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
764  
765      temp = useSF;
766 <    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
766 >    MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
767  
768 +    temp = useSP;
769 +    MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
770 +
771 +    temp = useBoxDipole;
772 +    MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
773 +
774 +    temp = useAtomicVirial_;
775 +    MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
776 +
777   #endif
778  
779      fInfo_.SIM_uses_PBC = usePBC;    
# Line 761 | Line 791 | namespace oopse {
791      fInfo_.SIM_uses_FLARB = useFLARB;
792      fInfo_.SIM_uses_RF = useRF;
793      fInfo_.SIM_uses_SF = useSF;
794 <
795 <    if( myMethod == "REACTION_FIELD") {
796 <      
767 <      if (simParams_->haveDielectric()) {
768 <        fInfo_.dielect = simParams_->getDielectric();
769 <      } else {
770 <        sprintf(painCave.errMsg,
771 <                "SimSetup Error: No Dielectric constant was set.\n"
772 <                "\tYou are trying to use Reaction Field without"
773 <                "\tsetting a dielectric constant!\n");
774 <        painCave.isFatal = 1;
775 <        simError();
776 <      }      
777 <    }
778 <
794 >    fInfo_.SIM_uses_SP = useSP;
795 >    fInfo_.SIM_uses_BoxDipole = useBoxDipole;
796 >    fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_;
797    }
798  
799    void SimInfo::setupFortranSim() {
# Line 792 | Line 810 | namespace oopse {
810      }
811  
812      //calculate mass ratio of cutoff group
813 <    std::vector<double> mfact;
813 >    std::vector<RealType> mfact;
814      SimInfo::MoleculeIterator mi;
815      Molecule* mol;
816      Molecule::CutoffGroupIterator ci;
817      CutoffGroup* cg;
818      Molecule::AtomIterator ai;
819      Atom* atom;
820 <    double totalMass;
820 >    RealType totalMass;
821  
822      //to avoid memory reallocation, reserve enough space for mfact
823      mfact.reserve(getNCutoffGroups());
# Line 860 | Line 878 | namespace oopse {
878               "succesfully sent the simulation information to fortran.\n");
879      MPIcheckPoint();
880   #endif // is_mpi
881 +
882 +    // Setup number of neighbors in neighbor list if present
883 +    if (simParams_->haveNeighborListNeighbors()) {
884 +      int nlistNeighbors = simParams_->getNeighborListNeighbors();
885 +      setNeighbors(&nlistNeighbors);
886 +    }
887 +  
888 +
889    }
890  
891  
# Line 924 | Line 950 | namespace oopse {
950  
951    void SimInfo::setupCutoff() {          
952      
953 +    ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
954 +
955      // Check the cutoff policy
956 <    int cp =  TRADITIONAL_CUTOFF_POLICY;
957 <    if (simParams_->haveCutoffPolicy()) {
958 <      std::string myPolicy = simParams_->getCutoffPolicy();
956 >    int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default
957 >
958 >    // Set LJ shifting bools to false
959 >    ljsp_ = false;
960 >    ljsf_ = false;
961 >
962 >    std::string myPolicy;
963 >    if (forceFieldOptions_.haveCutoffPolicy()){
964 >      myPolicy = forceFieldOptions_.getCutoffPolicy();
965 >    }else if (simParams_->haveCutoffPolicy()) {
966 >      myPolicy = simParams_->getCutoffPolicy();
967 >    }
968 >
969 >    if (!myPolicy.empty()){
970        toUpper(myPolicy);
971        if (myPolicy == "MIX") {
972          cp = MIX_CUTOFF_POLICY;
# Line 950 | Line 989 | namespace oopse {
989      notifyFortranCutoffPolicy(&cp);
990  
991      // Check the Skin Thickness for neighborlists
992 <    double skin;
992 >    RealType skin;
993      if (simParams_->haveSkinThickness()) {
994        skin = simParams_->getSkinThickness();
995        notifyFortranSkinThickness(&skin);
# Line 962 | Line 1001 | namespace oopse {
1001        if (simParams_->haveSwitchingRadius()) {
1002          rsw_  = simParams_->getSwitchingRadius();
1003        } else {
1004 <        rsw_ = rcut_;
1004 >        if (fInfo_.SIM_uses_Charges |
1005 >            fInfo_.SIM_uses_Dipoles |
1006 >            fInfo_.SIM_uses_RF) {
1007 >          
1008 >          rsw_ = 0.85 * rcut_;
1009 >          sprintf(painCave.errMsg,
1010 >                  "SimCreator Warning: No value was set for the switchingRadius.\n"
1011 >                  "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n"
1012 >                  "\tswitchingRadius = %f. for this simulation\n", rsw_);
1013 >        painCave.isFatal = 0;
1014 >        simError();
1015 >        } else {
1016 >          rsw_ = rcut_;
1017 >          sprintf(painCave.errMsg,
1018 >                  "SimCreator Warning: No value was set for the switchingRadius.\n"
1019 >                  "\tOOPSE will use the same value as the cutoffRadius.\n"
1020 >                  "\tswitchingRadius = %f. for this simulation\n", rsw_);
1021 >          painCave.isFatal = 0;
1022 >          simError();
1023 >        }
1024        }
1025 <      notifyFortranCutoffs(&rcut_, &rsw_);
1025 >
1026 >      if (simParams_->haveElectrostaticSummationMethod()) {
1027 >        std::string myMethod = simParams_->getElectrostaticSummationMethod();
1028 >        toUpper(myMethod);
1029 >        
1030 >        if (myMethod == "SHIFTED_POTENTIAL") {
1031 >          ljsp_ = true;
1032 >        } else if (myMethod == "SHIFTED_FORCE") {
1033 >          ljsf_ = true;
1034 >        }
1035 >      }
1036 >      notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
1037        
1038      } else {
1039        
# Line 981 | Line 1050 | namespace oopse {
1050          if (simParams_->haveElectrostaticSummationMethod()) {
1051            std::string myMethod = simParams_->getElectrostaticSummationMethod();
1052            toUpper(myMethod);
1053 <          if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") {
1053 >      
1054 >      // For the time being, we're tethering the LJ shifted behavior to the
1055 >      // electrostaticSummationMethod keyword options
1056 >          if (myMethod == "SHIFTED_POTENTIAL") {
1057 >            ljsp_ = true;
1058 >          } else if (myMethod == "SHIFTED_FORCE") {
1059 >            ljsf_ = true;
1060 >          }
1061 >          if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") {
1062              if (simParams_->haveSwitchingRadius()){
1063                sprintf(painCave.errMsg,
1064                        "SimInfo Warning: A value was set for the switchingRadius\n"
# Line 1004 | Line 1081 | namespace oopse {
1081            simError();
1082            rsw_ = 0.85 * rcut_;
1083          }
1084 <        notifyFortranCutoffs(&rcut_, &rsw_);
1084 >
1085 >        notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
1086 >
1087        } else {
1088          // We didn't set rcut explicitly, and we don't have electrostatic atoms, so
1089          // We'll punt and let fortran figure out the cutoffs later.
# Line 1020 | Line 1099 | namespace oopse {
1099      int errorOut;
1100      int esm =  NONE;
1101      int sm = UNDAMPED;
1102 <    double alphaVal;
1103 <    double dielectric;
1104 <
1102 >    RealType alphaVal;
1103 >    RealType dielectric;
1104 >    
1105      errorOut = isError;
1027 –    alphaVal = simParams_->getDampingAlpha();
1028 –    dielectric = simParams_->getDielectric();
1106  
1107      if (simParams_->haveElectrostaticSummationMethod()) {
1108        std::string myMethod = simParams_->getElectrostaticSummationMethod();
# Line 1042 | Line 1119 | namespace oopse {
1119              if (myMethod == "SHIFTED_FORCE") {            
1120                esm = SHIFTED_FORCE;
1121              } else {
1122 <              if (myMethod == "REACTION_FIELD") {            
1122 >              if (myMethod == "REACTION_FIELD") {
1123                  esm = REACTION_FIELD;
1124 +                dielectric = simParams_->getDielectric();
1125 +                if (!simParams_->haveDielectric()) {
1126 +                  // throw warning
1127 +                  sprintf( painCave.errMsg,
1128 +                           "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
1129 +                           "\tA default value of %f will be used for the dielectric.\n", dielectric);
1130 +                  painCave.isFatal = 0;
1131 +                  simError();
1132 +                }
1133                } else {
1134                  // throw error        
1135                  sprintf( painCave.errMsg,
# Line 1070 | Line 1156 | namespace oopse {
1156          if (myScreen == "DAMPED") {
1157            sm = DAMPED;
1158            if (!simParams_->haveDampingAlpha()) {
1159 <            //throw error
1159 >            // first set a cutoff dependent alpha value
1160 >            // we assume alpha depends linearly with rcut from 0 to 20.5 ang
1161 >            alphaVal = 0.5125 - rcut_* 0.025;
1162 >            // for values rcut > 20.5, alpha is zero
1163 >            if (alphaVal < 0) alphaVal = 0;
1164 >
1165 >            // throw warning
1166              sprintf( painCave.errMsg,
1167                       "SimInfo warning: dampingAlpha was not specified in the input file.\n"
1168 <                     "\tA default value of %f (1/ang) will be used.\n", alphaVal);
1168 >                     "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_);
1169              painCave.isFatal = 0;
1170              simError();
1171 +          } else {
1172 +            alphaVal = simParams_->getDampingAlpha();
1173            }
1174 +          
1175          } else {
1176            // throw error        
1177            sprintf( painCave.errMsg,
# Line 1092 | Line 1187 | namespace oopse {
1187      
1188      // let's pass some summation method variables to fortran
1189      setElectrostaticSummationMethod( &esm );
1190 <    notifyFortranElectrostaticMethod( &esm );
1190 >    setFortranElectrostaticMethod( &esm );
1191      setScreeningMethod( &sm );
1192      setDampingAlpha( &alphaVal );
1193      setReactionFieldDielectric( &dielectric );
# Line 1125 | Line 1220 | namespace oopse {
1220  
1221    }
1222  
1223 +  void SimInfo::setupAccumulateBoxDipole() {    
1224 +
1225 +    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
1226 +    if ( simParams_->haveAccumulateBoxDipole() )
1227 +      if ( simParams_->getAccumulateBoxDipole() ) {
1228 +        setAccumulateBoxDipole();
1229 +        calcBoxDipole_ = true;
1230 +      }
1231 +
1232 +  }
1233 +
1234    void SimInfo::addProperty(GenericData* genData) {
1235      properties_.addProperty(genData);  
1236    }
# Line 1181 | Line 1287 | namespace oopse {
1287      Molecule* mol;
1288  
1289      Vector3d comVel(0.0);
1290 <    double totalMass = 0.0;
1290 >    RealType totalMass = 0.0;
1291      
1292  
1293      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1294 <      double mass = mol->getMass();
1294 >      RealType mass = mol->getMass();
1295        totalMass += mass;
1296        comVel += mass * mol->getComVel();
1297      }  
1298  
1299   #ifdef IS_MPI
1300 <    double tmpMass = totalMass;
1300 >    RealType tmpMass = totalMass;
1301      Vector3d tmpComVel(comVel);    
1302 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1303 <    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1302 >    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1303 >    MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1304   #endif
1305  
1306      comVel /= totalMass;
# Line 1207 | Line 1313 | namespace oopse {
1313      Molecule* mol;
1314  
1315      Vector3d com(0.0);
1316 <    double totalMass = 0.0;
1316 >    RealType totalMass = 0.0;
1317      
1318      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1319 <      double mass = mol->getMass();
1319 >      RealType mass = mol->getMass();
1320        totalMass += mass;
1321        com += mass * mol->getCom();
1322      }  
1323  
1324   #ifdef IS_MPI
1325 <    double tmpMass = totalMass;
1325 >    RealType tmpMass = totalMass;
1326      Vector3d tmpCom(com);    
1327 <    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1328 <    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1327 >    MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1328 >    MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1329   #endif
1330  
1331      com /= totalMass;
# Line 1243 | Line 1349 | namespace oopse {
1349        Molecule* mol;
1350        
1351      
1352 <      double totalMass = 0.0;
1352 >      RealType totalMass = 0.0;
1353      
1354  
1355        for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1356 <         double mass = mol->getMass();
1356 >         RealType mass = mol->getMass();
1357           totalMass += mass;
1358           com += mass * mol->getCom();
1359           comVel += mass * mol->getComVel();          
1360        }  
1361        
1362   #ifdef IS_MPI
1363 <      double tmpMass = totalMass;
1363 >      RealType tmpMass = totalMass;
1364        Vector3d tmpCom(com);  
1365        Vector3d tmpComVel(comVel);
1366 <      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1367 <      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1368 <      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1366 >      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1367 >      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1368 >      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1369   #endif
1370        
1371        com /= totalMass;
# Line 1278 | Line 1384 | namespace oopse {
1384     void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1385        
1386  
1387 <      double xx = 0.0;
1388 <      double yy = 0.0;
1389 <      double zz = 0.0;
1390 <      double xy = 0.0;
1391 <      double xz = 0.0;
1392 <      double yz = 0.0;
1387 >      RealType xx = 0.0;
1388 >      RealType yy = 0.0;
1389 >      RealType zz = 0.0;
1390 >      RealType xy = 0.0;
1391 >      RealType xz = 0.0;
1392 >      RealType yz = 0.0;
1393        Vector3d com(0.0);
1394        Vector3d comVel(0.0);
1395        
# Line 1295 | Line 1401 | namespace oopse {
1401        Vector3d thisq(0.0);
1402        Vector3d thisv(0.0);
1403  
1404 <      double thisMass = 0.0;
1404 >      RealType thisMass = 0.0;
1405      
1406        
1407        
# Line 1333 | Line 1439 | namespace oopse {
1439   #ifdef IS_MPI
1440        Mat3x3d tmpI(inertiaTensor);
1441        Vector3d tmpAngMom;
1442 <      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1443 <      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1442 >      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1443 >      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1444   #endif
1445                
1446        return;
# Line 1355 | Line 1461 | namespace oopse {
1461        Vector3d thisr(0.0);
1462        Vector3d thisp(0.0);
1463        
1464 <      double thisMass;
1464 >      RealType thisMass;
1465        
1466        for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1467          thisMass = mol->getMass();
# Line 1368 | Line 1474 | namespace oopse {
1474        
1475   #ifdef IS_MPI
1476        Vector3d tmpAngMom;
1477 <      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1477 >      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1478   #endif
1479        
1480        return angularMomentum;
1481     }
1482    
1483 <  
1483 >  StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
1484 >    return IOIndexToIntegrableObject.at(index);
1485 >  }
1486 >  
1487 >  void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) {
1488 >    IOIndexToIntegrableObject= v;
1489 >  }
1490 >
1491 >  /* Returns the Volume of the simulation based on a ellipsoid with semi-axes
1492 >     based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3
1493 >     where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to
1494 >     V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
1495 >  */
1496 >  void SimInfo::getGyrationalVolume(RealType &volume){
1497 >    Mat3x3d intTensor;
1498 >    RealType det;
1499 >    Vector3d dummyAngMom;
1500 >    RealType sysconstants;
1501 >    RealType geomCnst;
1502 >
1503 >    geomCnst = 3.0/2.0;
1504 >    /* Get the inertial tensor and angular momentum for free*/
1505 >    getInertiaTensor(intTensor,dummyAngMom);
1506 >    
1507 >    det = intTensor.determinant();
1508 >    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1509 >    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det);
1510 >    return;
1511 >  }
1512 >
1513 >  void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){
1514 >    Mat3x3d intTensor;
1515 >    Vector3d dummyAngMom;
1516 >    RealType sysconstants;
1517 >    RealType geomCnst;
1518 >
1519 >    geomCnst = 3.0/2.0;
1520 >    /* Get the inertial tensor and angular momentum for free*/
1521 >    getInertiaTensor(intTensor,dummyAngMom);
1522 >    
1523 >    detI = intTensor.determinant();
1524 >    sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1525 >    volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI);
1526 >    return;
1527 >  }
1528 > /*
1529 >   void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) {
1530 >      assert( v.size() == nAtoms_ + nRigidBodies_);
1531 >      sdByGlobalIndex_ = v;
1532 >    }
1533 >
1534 >    StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) {
1535 >      //assert(index < nAtoms_ + nRigidBodies_);
1536 >      return sdByGlobalIndex_.at(index);
1537 >    }  
1538 > */  
1539   }//end namespace oopse
1540  

Diff Legend

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