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 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 629 by chuckv, Mon Sep 26 15:58:17 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 511 | Line 516 | namespace oopse {
516      int useDipole = 0;
517      int useGayBerne = 0;
518      int useSticky = 0;
519 +    int useStickyPower = 0;
520      int useShape = 0;
521      int useFLARB = 0; //it is not in AtomType yet
522      int useDirectionalAtom = 0;    
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 529 | Line 542 | namespace oopse {
542        useDipole |= (*i)->isDipole();
543        useGayBerne |= (*i)->isGayBerne();
544        useSticky |= (*i)->isSticky();
545 +      useStickyPower |= (*i)->isStickyPower();
546        useShape |= (*i)->isShape();
547      }
548  
549 <    if (useSticky || useDipole || useGayBerne || useShape) {
549 >    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
550        useDirectionalAtom = 1;
551      }
552  
# Line 564 | Line 578 | namespace oopse {
578      temp = useSticky;
579      MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
580  
581 +    temp = useStickyPower;
582 +    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
583 +    
584      temp = useGayBerne;
585      MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
586  
# Line 578 | 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 588 | Line 605 | namespace oopse {
605      fInfo_.SIM_uses_Charges = useCharge;
606      fInfo_.SIM_uses_Dipoles = useDipole;
607      fInfo_.SIM_uses_Sticky = useSticky;
608 +    fInfo_.SIM_uses_StickyPower = useStickyPower;
609      fInfo_.SIM_uses_GayBerne = useGayBerne;
610      fInfo_.SIM_uses_EAM = useEAM;
611      fInfo_.SIM_uses_Shapes = useShape;
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 824 | 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 >
873 >
874 >    if (simParams_->haveSkinThickness()) {
875 >      double skinThickness = simParams_->getSkinThickness();
876 >    }
877 >
878 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
879 >    // also send cutoff notification to electrostatics
880 >    setElectrostaticCutoffRadius(&rcut_);
881 >  }
882 >
883 >  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
884 >    
885 >    int errorOut;
886 >    int esm =  NONE;
887 >    double alphaVal;
888 >    double dielectric;
889 >
890 >    errorOut = isError;
891 >    alphaVal = simParams_->getDampingAlpha();
892 >    dielectric = simParams_->getDielectric();
893 >
894 >    if (simParams_->haveElectrostaticSummationMethod()) {
895 >      std::string myMethod = simParams_->getElectrostaticSummationMethod();
896 >      if (myMethod == "NONE") {
897 >        esm = NONE;
898 >      } else {
899 >        if (myMethod == "UNDAMPED_WOLF") {
900 >          esm = UNDAMPED_WOLF;
901 >        } else {
902 >          if (myMethod == "DAMPED_WOLF") {            
903 >            esm = DAMPED_WOLF;
904 >            if (!simParams_->haveDampingAlpha()) {
905 >              //throw error
906 >              sprintf( painCave.errMsg,
907 >                       "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);
908 >              painCave.isFatal = 0;
909 >              simError();
910 >            }
911 >          } else {
912 >            if (myMethod == "REACTION_FIELD") {
913 >              esm = REACTION_FIELD;
914 >            } else {
915 >              // throw error        
916 >              sprintf( painCave.errMsg,
917 >                       "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() );
918 >              painCave.isFatal = 1;
919 >              simError();
920 >            }    
921 >          }          
922 >        }
923 >      }
924 >    }
925 >    // let's pass some summation method variables to fortran
926 >    setElectrostaticSummationMethod( &esm );
927 >    setDampedWolfAlpha( &alphaVal );
928 >    setReactionFieldDielectric( &dielectric );
929 >    initFortranFF( &esm, &errorOut );
930    }
931  
932    void SimInfo::addProperty(GenericData* genData) {
# Line 939 | Line 1036 | namespace oopse {
1036  
1037      return o;
1038    }
1039 +  
1040 +  
1041 +   /*
1042 +   Returns center of mass and center of mass velocity in one function call.
1043 +   */
1044 +  
1045 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1046 +      SimInfo::MoleculeIterator i;
1047 +      Molecule* mol;
1048 +      
1049 +    
1050 +      double totalMass = 0.0;
1051 +    
1052  
1053 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1054 +         double mass = mol->getMass();
1055 +         totalMass += mass;
1056 +         com += mass * mol->getCom();
1057 +         comVel += mass * mol->getComVel();          
1058 +      }  
1059 +      
1060 + #ifdef IS_MPI
1061 +      double tmpMass = totalMass;
1062 +      Vector3d tmpCom(com);  
1063 +      Vector3d tmpComVel(comVel);
1064 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1065 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1066 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1067 + #endif
1068 +      
1069 +      com /= totalMass;
1070 +      comVel /= totalMass;
1071 +   }        
1072 +  
1073 +   /*
1074 +   Return intertia tensor for entire system and angular momentum Vector.
1075 +
1076 +
1077 +       [  Ixx -Ixy  -Ixz ]
1078 +  J =| -Iyx  Iyy  -Iyz |
1079 +       [ -Izx -Iyz   Izz ]
1080 +    */
1081 +
1082 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1083 +      
1084 +
1085 +      double xx = 0.0;
1086 +      double yy = 0.0;
1087 +      double zz = 0.0;
1088 +      double xy = 0.0;
1089 +      double xz = 0.0;
1090 +      double yz = 0.0;
1091 +      Vector3d com(0.0);
1092 +      Vector3d comVel(0.0);
1093 +      
1094 +      getComAll(com, comVel);
1095 +      
1096 +      SimInfo::MoleculeIterator i;
1097 +      Molecule* mol;
1098 +      
1099 +      Vector3d thisq(0.0);
1100 +      Vector3d thisv(0.0);
1101 +
1102 +      double thisMass = 0.0;
1103 +    
1104 +      
1105 +      
1106 +  
1107 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1108 +        
1109 +         thisq = mol->getCom()-com;
1110 +         thisv = mol->getComVel()-comVel;
1111 +         thisMass = mol->getMass();
1112 +         // Compute moment of intertia coefficients.
1113 +         xx += thisq[0]*thisq[0]*thisMass;
1114 +         yy += thisq[1]*thisq[1]*thisMass;
1115 +         zz += thisq[2]*thisq[2]*thisMass;
1116 +        
1117 +         // compute products of intertia
1118 +         xy += thisq[0]*thisq[1]*thisMass;
1119 +         xz += thisq[0]*thisq[2]*thisMass;
1120 +         yz += thisq[1]*thisq[2]*thisMass;
1121 +            
1122 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1123 +            
1124 +      }  
1125 +      
1126 +      
1127 +      inertiaTensor(0,0) = yy + zz;
1128 +      inertiaTensor(0,1) = -xy;
1129 +      inertiaTensor(0,2) = -xz;
1130 +      inertiaTensor(1,0) = -xy;
1131 +      inertiaTensor(1,1) = xx + zz;
1132 +      inertiaTensor(1,2) = -yz;
1133 +      inertiaTensor(2,0) = -xz;
1134 +      inertiaTensor(2,1) = -yz;
1135 +      inertiaTensor(2,2) = xx + yy;
1136 +      
1137 + #ifdef IS_MPI
1138 +      Mat3x3d tmpI(inertiaTensor);
1139 +      Vector3d tmpAngMom;
1140 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1141 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1142 + #endif
1143 +              
1144 +      return;
1145 +   }
1146 +
1147 +   //Returns the angular momentum of the system
1148 +   Vector3d SimInfo::getAngularMomentum(){
1149 +      
1150 +      Vector3d com(0.0);
1151 +      Vector3d comVel(0.0);
1152 +      Vector3d angularMomentum(0.0);
1153 +      
1154 +      getComAll(com,comVel);
1155 +      
1156 +      SimInfo::MoleculeIterator i;
1157 +      Molecule* mol;
1158 +      
1159 +      Vector3d thisr(0.0);
1160 +      Vector3d thisp(0.0);
1161 +      
1162 +      double thisMass;
1163 +      
1164 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1165 +        thisMass = mol->getMass();
1166 +        thisr = mol->getCom()-com;
1167 +        thisp = (mol->getComVel()-comVel)*thisMass;
1168 +        
1169 +        angularMomentum += cross( thisr, thisp );
1170 +        
1171 +      }  
1172 +      
1173 + #ifdef IS_MPI
1174 +      Vector3d tmpAngMom;
1175 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1176 + #endif
1177 +      
1178 +      return angularMomentum;
1179 +   }
1180 +  
1181 +  
1182   }//end namespace oopse
1183  

Diff Legend

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