ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/brains/SimInfo.cpp
(Generate patch)

Comparing trunk/OOPSE-2.0/src/brains/SimInfo.cpp (file contents):
Revision 2252 by chuckv, Mon May 30 14:01:52 2005 UTC vs.
Revision 2344 by chrisfen, Tue Oct 4 19:34:03 2005 UTC

# 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/doForces_interface.h"
58 + #include "UseTheForce/DarkSide/electrostatic_interface.h"
59   #include "UseTheForce/notifyCutoffs_interface.h"
60   #include "utils/MemoryUtils.hpp"
61   #include "utils/simError.h"
# Line 80 | Line 83 | namespace oopse {
83        MoleculeStamp* molStamp;
84        int nMolWithSameStamp;
85        int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
86 <      int nGroups = 0;          //total cutoff groups defined in meta-data file
86 >      int nGroups = 0;      //total cutoff groups defined in meta-data file
87        CutoffGroupStamp* cgStamp;    
88        RigidBodyStamp* rbStamp;
89        int nRigidAtoms = 0;
# Line 105 | Line 108 | namespace oopse {
108          }
109  
110          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
111 +
112          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
113  
114          //calculate atoms in rigid bodies
# Line 121 | Line 125 | namespace oopse {
125          
126        }
127  
128 <      //every free atom (atom does not belong to cutoff groups) is a cutoff group
129 <      //therefore the total number of cutoff groups in the system is equal to
130 <      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
131 <      //file plus the number of cutoff groups defined in meta-data file
128 >      //every free atom (atom does not belong to cutoff groups) is a cutoff
129 >      //group therefore the total number of cutoff groups in the system is
130 >      //equal to the total number of atoms minus number of atoms belong to
131 >      //cutoff group defined in meta-data file plus the number of cutoff
132 >      //groups defined in meta-data file
133        nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
134  
135 <      //every free atom (atom does not belong to rigid bodies) is an integrable object
136 <      //therefore the total number of  integrable objects in the system is equal to
137 <      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
138 <      //file plus the number of  rigid bodies defined in meta-data file
139 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
140 <
135 >      //every free atom (atom does not belong to rigid bodies) is an
136 >      //integrable object therefore the total number of integrable objects
137 >      //in the system is equal to the total number of atoms minus number of
138 >      //atoms belong to rigid body defined in meta-data file plus the number
139 >      //of rigid bodies defined in meta-data file
140 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
141 >                                                + nGlobalRigidBodies_;
142 >  
143        nGlobalMols_ = molStampIds_.size();
144  
145   #ifdef IS_MPI    
# Line 462 | Line 469 | namespace oopse {
469      //setup fortran force field
470      /** @deprecate */    
471      int isError = 0;
472 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
472 >    
473 >    setupElectrostaticSummationMethod( isError );
474 >
475      if(isError){
476        sprintf( painCave.errMsg,
477                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 518 | Line 527 | namespace oopse {
527      int useElectrostatics = 0;
528      //usePBC and useRF are from simParams
529      int usePBC = simParams_->getPBC();
530 <    int useRF = simParams_->getUseRF();
530 >    int useRF;
531  
532 +    // set the useRF logical
533 +    std::string myMethod = simParams_->getElectrostaticSummationMethod();
534 +    if (myMethod == "REACTION_FIELD")
535 +      useRF = 1;
536 +    else
537 +      useRF = 0;
538 +
539      //loop over all of the atom types
540      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
541        useLennardJones |= (*i)->isLennardJones();
# Line 583 | Line 599 | namespace oopse {
599  
600      temp = useRF;
601      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
602 <    
602 >
603   #endif
604  
605      fInfo_.SIM_uses_PBC = usePBC;    
# Line 600 | Line 616 | namespace oopse {
616      fInfo_.SIM_uses_FLARB = useFLARB;
617      fInfo_.SIM_uses_RF = useRF;
618  
619 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
619 >    if( fInfo_.SIM_uses_Dipoles && myMethod == "REACTION_FIELD") {
620  
621        if (simParams_->haveDielectric()) {
622          fInfo_.dielect = simParams_->getDielectric();
# Line 650 | Line 666 | namespace oopse {
666  
667          totalMass = cg->getMass();
668          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
669 <          mfact.push_back(atom->getMass()/totalMass);
669 >          // Check for massless groups - set mfact to 1 if true
670 >          if (totalMass != 0)
671 >            mfact.push_back(atom->getMass()/totalMass);
672 >          else
673 >            mfact.push_back( 1.0 );
674          }
675  
676        }      
# Line 830 | Line 850 | namespace oopse {
850      }
851    }
852  
853 <  void SimInfo::setupCutoff() {
853 >  void SimInfo::setupCutoff() {    
854      getCutoff(rcut_, rsw_);    
855      double rnblist = rcut_ + 1; // skin of neighbor list
856  
857      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
858 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
858 >    
859 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
860 >    if (simParams_->haveCutoffPolicy()) {
861 >      std::string myPolicy = simParams_->getCutoffPolicy();
862 >      if (myPolicy == "MIX") {
863 >        cp = MIX_CUTOFF_POLICY;
864 >      } else {
865 >        if (myPolicy == "MAX") {
866 >          cp = MAX_CUTOFF_POLICY;
867 >        } else {
868 >          if (myPolicy == "TRADITIONAL") {            
869 >            cp = TRADITIONAL_CUTOFF_POLICY;
870 >          } else {
871 >            // throw error        
872 >            sprintf( painCave.errMsg,
873 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
874 >            painCave.isFatal = 1;
875 >            simError();
876 >          }    
877 >        }          
878 >      }
879 >    }
880 >
881 >
882 >    if (simParams_->haveSkinThickness()) {
883 >      double skinThickness = simParams_->getSkinThickness();
884 >    }
885 >
886 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
887 >    // also send cutoff notification to electrostatics
888 >    setElectrostaticCutoffRadius(&rcut_);
889    }
890  
891 +  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
892 +    
893 +    int errorOut;
894 +    int esm =  NONE;
895 +    double alphaVal;
896 +    double dielectric;
897 +
898 +    errorOut = isError;
899 +    alphaVal = simParams_->getDampingAlpha();
900 +    dielectric = simParams_->getDielectric();
901 +
902 +    if (simParams_->haveElectrostaticSummationMethod()) {
903 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
904 +      if (myMethod == "NONE") {
905 +        esm = NONE;
906 +      } else {
907 +        if (myMethod == "UNDAMPED_WOLF") {
908 +          esm = UNDAMPED_WOLF;
909 +        } else {
910 +          if (myMethod == "DAMPED_WOLF") {            
911 +            esm = DAMPED_WOLF;
912 +            if (!simParams_->haveDampingAlpha()) {
913 +              //throw error
914 +              sprintf( painCave.errMsg,
915 +                       "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", alphaVal);
916 +              painCave.isFatal = 0;
917 +              simError();
918 +            }
919 +          } else {
920 +            if (myMethod == "REACTION_FIELD") {
921 +              esm = REACTION_FIELD;
922 +            } else {
923 +              // throw error        
924 +              sprintf( painCave.errMsg,
925 +                       "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"undamped_wolf\", \"damped_wolf\", or \"reaction_field\".", myMethod.c_str() );
926 +              painCave.isFatal = 1;
927 +              simError();
928 +            }    
929 +          }          
930 +        }
931 +      }
932 +    }
933 +    // let's pass some summation method variables to fortran
934 +    setElectrostaticSummationMethod( &esm );
935 +    setDampedWolfAlpha( &alphaVal );
936 +    setReactionFieldDielectric( &dielectric );
937 +    initFortranFF( &esm, &errorOut );
938 +  }
939 +
940    void SimInfo::addProperty(GenericData* genData) {
941      properties_.addProperty(genData);  
942    }
# Line 981 | Line 1080 | namespace oopse {
1080    
1081     /*
1082     Return intertia tensor for entire system and angular momentum Vector.
1083 +
1084 +
1085 +       [  Ixx -Ixy  -Ixz ]
1086 +  J =| -Iyx  Iyy  -Iyz |
1087 +       [ -Izx -Iyz   Izz ]
1088      */
1089  
1090     void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
# Line 1032 | Line 1136 | namespace oopse {
1136        inertiaTensor(0,1) = -xy;
1137        inertiaTensor(0,2) = -xz;
1138        inertiaTensor(1,0) = -xy;
1139 <      inertiaTensor(2,0) = xx + zz;
1139 >      inertiaTensor(1,1) = xx + zz;
1140        inertiaTensor(1,2) = -yz;
1141        inertiaTensor(2,0) = -xz;
1142        inertiaTensor(2,1) = -yz;
# Line 1060 | Line 1164 | namespace oopse {
1164        SimInfo::MoleculeIterator i;
1165        Molecule* mol;
1166        
1167 <      Vector3d thisq(0.0);
1168 <      Vector3d thisv(0.0);
1167 >      Vector3d thisr(0.0);
1168 >      Vector3d thisp(0.0);
1169        
1170 <      double thisMass = 0.0;
1170 >      double thisMass;
1171        
1172        for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1173 <         thisq = mol->getCom()-com;
1174 <         thisv = mol->getComVel()-comVel;
1175 <         thisMass = mol->getMass();
1072 <         angularMomentum += cross( thisq, thisv ) * thisMass;
1173 >        thisMass = mol->getMass();
1174 >        thisr = mol->getCom()-com;
1175 >        thisp = (mol->getComVel()-comVel)*thisMass;
1176          
1177 +        angularMomentum += cross( thisr, thisp );
1178 +        
1179        }  
1180        
1181   #ifdef IS_MPI

Diff Legend

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