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 523 by chrisfen, Thu May 5 14:47:35 2005 UTC vs.
Revision 611 by chrisfen, Mon Sep 19 23:21:46 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();
526 <    int useRF = simParams_->getUseRF();
526 >    int useRF;
527  
528 +    // set the useRF logical
529 +    std::string myMethod = simParams_->getElectrostaticSummationMethod();
530 +    if (myMethod == "REACTION_FIELD")
531 +      useRF = 1;
532 +    else
533 +      useRF = 0;
534 +
535      //loop over all of the atom types
536      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
537        useLennardJones |= (*i)->isLennardJones();
# Line 583 | Line 595 | namespace oopse {
595  
596      temp = useRF;
597      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
598 <    
598 >
599   #endif
600  
601      fInfo_.SIM_uses_PBC = usePBC;    
# Line 600 | Line 612 | namespace oopse {
612      fInfo_.SIM_uses_FLARB = useFLARB;
613      fInfo_.SIM_uses_RF = useRF;
614  
615 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
615 >    if( fInfo_.SIM_uses_Dipoles && myMethod == "REACTION_FIELD") {
616  
617        if (simParams_->haveDielectric()) {
618          fInfo_.dielect = simParams_->getDielectric();
# 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 >    // also send cutoff notification to electrostatics
874 >    setElectrostaticCutoffRadius(&rcut_);
875    }
876  
877 +  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
878 +    
879 +    int errorOut;
880 +    int esm =  NONE;
881 +    double alphaVal;
882 +    double dielectric;
883 +
884 +    errorOut = isError;
885 +    alphaVal = simParams_->getDampingAlpha();
886 +    dielectric = simParams_->getDielectric();
887 +
888 +    if (simParams_->haveElectrostaticSummationMethod()) {
889 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
890 +      if (myMethod == "NONE") {
891 +        esm = NONE;
892 +      } else {
893 +        if (myMethod == "UNDAMPED_WOLF") {
894 +          esm = UNDAMPED_WOLF;
895 +        } else {
896 +          if (myMethod == "DAMPED_WOLF") {            
897 +            esm = DAMPED_WOLF;
898 +            if (!simParams_->haveDampingAlpha()) {
899 +              //throw error
900 +              sprintf( painCave.errMsg,
901 +                       "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);
902 +              painCave.isFatal = 0;
903 +              simError();
904 +            }
905 +          } else {
906 +            if (myMethod == "REACTION_FIELD") {
907 +              esm = REACTION_FIELD;
908 +            } else {
909 +              // throw error        
910 +              sprintf( painCave.errMsg,
911 +                       "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() );
912 +              painCave.isFatal = 1;
913 +              simError();
914 +            }    
915 +          }          
916 +        }
917 +      }
918 +    }
919 +    // let's pass some summation method variables to fortran
920 +    setElectrostaticSummationMethod( &esm );
921 +    setDampedWolfAlpha( &alphaVal );
922 +    setReactionFieldDielectric( &dielectric );
923 +    initFortranFF( &esm, &errorOut );
924 +  }
925 +
926    void SimInfo::addProperty(GenericData* genData) {
927      properties_.addProperty(genData);  
928    }
# Line 945 | Line 1030 | namespace oopse {
1030  
1031      return o;
1032    }
1033 +  
1034 +  
1035 +   /*
1036 +   Returns center of mass and center of mass velocity in one function call.
1037 +   */
1038 +  
1039 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1040 +      SimInfo::MoleculeIterator i;
1041 +      Molecule* mol;
1042 +      
1043 +    
1044 +      double totalMass = 0.0;
1045 +    
1046  
1047 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1048 +         double mass = mol->getMass();
1049 +         totalMass += mass;
1050 +         com += mass * mol->getCom();
1051 +         comVel += mass * mol->getComVel();          
1052 +      }  
1053 +      
1054 + #ifdef IS_MPI
1055 +      double tmpMass = totalMass;
1056 +      Vector3d tmpCom(com);  
1057 +      Vector3d tmpComVel(comVel);
1058 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1059 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1060 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1061 + #endif
1062 +      
1063 +      com /= totalMass;
1064 +      comVel /= totalMass;
1065 +   }        
1066 +  
1067 +   /*
1068 +   Return intertia tensor for entire system and angular momentum Vector.
1069 +
1070 +
1071 +       [  Ixx -Ixy  -Ixz ]
1072 +  J =| -Iyx  Iyy  -Iyz |
1073 +       [ -Izx -Iyz   Izz ]
1074 +    */
1075 +
1076 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1077 +      
1078 +
1079 +      double xx = 0.0;
1080 +      double yy = 0.0;
1081 +      double zz = 0.0;
1082 +      double xy = 0.0;
1083 +      double xz = 0.0;
1084 +      double yz = 0.0;
1085 +      Vector3d com(0.0);
1086 +      Vector3d comVel(0.0);
1087 +      
1088 +      getComAll(com, comVel);
1089 +      
1090 +      SimInfo::MoleculeIterator i;
1091 +      Molecule* mol;
1092 +      
1093 +      Vector3d thisq(0.0);
1094 +      Vector3d thisv(0.0);
1095 +
1096 +      double thisMass = 0.0;
1097 +    
1098 +      
1099 +      
1100 +  
1101 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1102 +        
1103 +         thisq = mol->getCom()-com;
1104 +         thisv = mol->getComVel()-comVel;
1105 +         thisMass = mol->getMass();
1106 +         // Compute moment of intertia coefficients.
1107 +         xx += thisq[0]*thisq[0]*thisMass;
1108 +         yy += thisq[1]*thisq[1]*thisMass;
1109 +         zz += thisq[2]*thisq[2]*thisMass;
1110 +        
1111 +         // compute products of intertia
1112 +         xy += thisq[0]*thisq[1]*thisMass;
1113 +         xz += thisq[0]*thisq[2]*thisMass;
1114 +         yz += thisq[1]*thisq[2]*thisMass;
1115 +            
1116 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1117 +            
1118 +      }  
1119 +      
1120 +      
1121 +      inertiaTensor(0,0) = yy + zz;
1122 +      inertiaTensor(0,1) = -xy;
1123 +      inertiaTensor(0,2) = -xz;
1124 +      inertiaTensor(1,0) = -xy;
1125 +      inertiaTensor(1,1) = xx + zz;
1126 +      inertiaTensor(1,2) = -yz;
1127 +      inertiaTensor(2,0) = -xz;
1128 +      inertiaTensor(2,1) = -yz;
1129 +      inertiaTensor(2,2) = xx + yy;
1130 +      
1131 + #ifdef IS_MPI
1132 +      Mat3x3d tmpI(inertiaTensor);
1133 +      Vector3d tmpAngMom;
1134 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1135 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1136 + #endif
1137 +              
1138 +      return;
1139 +   }
1140 +
1141 +   //Returns the angular momentum of the system
1142 +   Vector3d SimInfo::getAngularMomentum(){
1143 +      
1144 +      Vector3d com(0.0);
1145 +      Vector3d comVel(0.0);
1146 +      Vector3d angularMomentum(0.0);
1147 +      
1148 +      getComAll(com,comVel);
1149 +      
1150 +      SimInfo::MoleculeIterator i;
1151 +      Molecule* mol;
1152 +      
1153 +      Vector3d thisr(0.0);
1154 +      Vector3d thisp(0.0);
1155 +      
1156 +      double thisMass;
1157 +      
1158 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1159 +        thisMass = mol->getMass();
1160 +        thisr = mol->getCom()-com;
1161 +        thisp = (mol->getComVel()-comVel)*thisMass;
1162 +        
1163 +        angularMomentum += cross( thisr, thisp );
1164 +        
1165 +      }  
1166 +      
1167 + #ifdef IS_MPI
1168 +      Vector3d tmpAngMom;
1169 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1170 + #endif
1171 +      
1172 +      return angularMomentum;
1173 +   }
1174 +  
1175 +  
1176   }//end namespace oopse
1177  

Diff Legend

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