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

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 555 by chuckv, Mon May 30 14:01:52 2005 UTC vs.
Revision 610 by chrisfen, Sun Sep 18 20:45:38 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 462 | Line 465 | namespace oopse {
465      //setup fortran force field
466      /** @deprecate */    
467      int isError = 0;
468 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
468 >    
469 >    setupElectrostaticSummationMethod( isError );
470 >
471      if(isError){
472        sprintf( painCave.errMsg,
473                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 518 | Line 523 | namespace oopse {
523      int useElectrostatics = 0;
524      //usePBC and useRF are from simParams
525      int usePBC = simParams_->getPBC();
521 –    int useRF = simParams_->getUseRF();
526  
527      //loop over all of the atom types
528      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
# Line 581 | Line 585 | namespace oopse {
585      temp = useFLARB;
586      MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
587  
584 –    temp = useRF;
585 –    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
586 –    
588   #endif
589  
590      fInfo_.SIM_uses_PBC = usePBC;    
# Line 598 | Line 599 | namespace oopse {
599      fInfo_.SIM_uses_EAM = useEAM;
600      fInfo_.SIM_uses_Shapes = useShape;
601      fInfo_.SIM_uses_FLARB = useFLARB;
601 –    fInfo_.SIM_uses_RF = useRF;
602  
603      if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
604  
# Line 830 | Line 830 | namespace oopse {
830      }
831    }
832  
833 <  void SimInfo::setupCutoff() {
833 >  void SimInfo::setupCutoff() {    
834      getCutoff(rcut_, rsw_);    
835      double rnblist = rcut_ + 1; // skin of neighbor list
836  
837      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
838 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
838 >    
839 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
840 >    if (simParams_->haveCutoffPolicy()) {
841 >      std::string myPolicy = simParams_->getCutoffPolicy();
842 >      if (myPolicy == "MIX") {
843 >        cp = MIX_CUTOFF_POLICY;
844 >      } else {
845 >        if (myPolicy == "MAX") {
846 >          cp = MAX_CUTOFF_POLICY;
847 >        } else {
848 >          if (myPolicy == "TRADITIONAL") {            
849 >            cp = TRADITIONAL_CUTOFF_POLICY;
850 >          } else {
851 >            // throw error        
852 >            sprintf( painCave.errMsg,
853 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
854 >            painCave.isFatal = 1;
855 >            simError();
856 >          }    
857 >        }          
858 >      }
859 >    }
860 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
861 >    // also send cutoff notification to electrostatics
862 >    setElectrostaticCutoffRadius(&rcut_);
863    }
864  
865 +  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
866 +    
867 +    int errorOut;
868 +    int esm =  NONE;
869 +    double alphaVal;
870 +    double dielectric;
871 +
872 +    errorOut = isError;
873 +    alphaVal = simParams_->getDampingAlpha();
874 +    dielectric = simParams_->getDielectric();
875 +
876 +    if (simParams_->haveElectrostaticSummationMethod()) {
877 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
878 +      if (myMethod == "NONE") {
879 +        esm = NONE;
880 +      } else {
881 +        if (myMethod == "UNDAMPED_WOLF") {
882 +          esm = UNDAMPED_WOLF;
883 +        } else {
884 +          if (myMethod == "DAMPED_WOLF") {            
885 +            esm = DAMPED_WOLF;
886 +            if (!simParams_->haveDampingAlpha()) {
887 +              //throw error
888 +              sprintf( painCave.errMsg,
889 +                       "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);
890 +              painCave.isFatal = 0;
891 +              simError();
892 +            }
893 +          } else {
894 +            if (myMethod == "REACTION_FIELD") {
895 +              esm = REACTION_FIELD;
896 +            } else {
897 +              // throw error        
898 +              sprintf( painCave.errMsg,
899 +                       "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() );
900 +              painCave.isFatal = 1;
901 +              simError();
902 +            }    
903 +          }          
904 +        }
905 +      }
906 +    }
907 +    // let's pass some summation method variables to fortran
908 +    setElectrostaticSummationMethod( &esm );
909 +    setDampedWolfAlpha( &alphaVal );
910 +    setReactionFieldDielectric( &dielectric );
911 +    initFortranFF( &esm, &errorOut );
912 +  }
913 +
914    void SimInfo::addProperty(GenericData* genData) {
915      properties_.addProperty(genData);  
916    }
# Line 981 | Line 1054 | namespace oopse {
1054    
1055     /*
1056     Return intertia tensor for entire system and angular momentum Vector.
1057 +
1058 +
1059 +       [  Ixx -Ixy  -Ixz ]
1060 +  J =| -Iyx  Iyy  -Iyz |
1061 +       [ -Izx -Iyz   Izz ]
1062      */
1063  
1064     void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
# Line 1032 | Line 1110 | namespace oopse {
1110        inertiaTensor(0,1) = -xy;
1111        inertiaTensor(0,2) = -xz;
1112        inertiaTensor(1,0) = -xy;
1113 <      inertiaTensor(2,0) = xx + zz;
1113 >      inertiaTensor(1,1) = xx + zz;
1114        inertiaTensor(1,2) = -yz;
1115        inertiaTensor(2,0) = -xz;
1116        inertiaTensor(2,1) = -yz;
# Line 1060 | Line 1138 | namespace oopse {
1138        SimInfo::MoleculeIterator i;
1139        Molecule* mol;
1140        
1141 <      Vector3d thisq(0.0);
1142 <      Vector3d thisv(0.0);
1141 >      Vector3d thisr(0.0);
1142 >      Vector3d thisp(0.0);
1143        
1144 <      double thisMass = 0.0;
1144 >      double thisMass;
1145        
1146        for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1147 <         thisq = mol->getCom()-com;
1148 <         thisv = mol->getComVel()-comVel;
1149 <         thisMass = mol->getMass();
1072 <         angularMomentum += cross( thisq, thisv ) * thisMass;
1147 >        thisMass = mol->getMass();
1148 >        thisr = mol->getCom()-com;
1149 >        thisp = (mol->getComVel()-comVel)*thisMass;
1150          
1151 +        angularMomentum += cross( thisr, thisp );
1152 +        
1153        }  
1154        
1155   #ifdef IS_MPI

Diff Legend

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