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 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 598 by chrisfen, Thu Sep 15 00:14:35 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/fCoulombicCorrection.h"
57   #include "UseTheForce/doForces_interface.h"
58   #include "UseTheForce/notifyCutoffs_interface.h"
59   #include "utils/MemoryUtils.hpp"
# Line 462 | Line 464 | namespace oopse {
464      //setup fortran force field
465      /** @deprecate */    
466      int isError = 0;
467 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
467 >    
468 >    setupCoulombicCorrection( isError );
469 >
470      if(isError){
471        sprintf( painCave.errMsg,
472                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 511 | Line 515 | namespace oopse {
515      int useDipole = 0;
516      int useGayBerne = 0;
517      int useSticky = 0;
518 +    int useStickyPower = 0;
519      int useShape = 0;
520      int useFLARB = 0; //it is not in AtomType yet
521      int useDirectionalAtom = 0;    
# Line 529 | Line 534 | namespace oopse {
534        useDipole |= (*i)->isDipole();
535        useGayBerne |= (*i)->isGayBerne();
536        useSticky |= (*i)->isSticky();
537 +      useStickyPower |= (*i)->isStickyPower();
538        useShape |= (*i)->isShape();
539      }
540  
541 <    if (useSticky || useDipole || useGayBerne || useShape) {
541 >    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
542        useDirectionalAtom = 1;
543      }
544  
# Line 564 | Line 570 | namespace oopse {
570      temp = useSticky;
571      MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
572  
573 +    temp = useStickyPower;
574 +    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
575 +    
576      temp = useGayBerne;
577      MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
578  
# Line 578 | Line 587 | namespace oopse {
587  
588      temp = useRF;
589      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
590 +
591 +    temp = useUW;
592 +    MPI_Allreduce(&temp, &useUW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
593 +
594 +    temp = useDW;
595 +    MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
596      
597   #endif
598  
# Line 588 | Line 603 | namespace oopse {
603      fInfo_.SIM_uses_Charges = useCharge;
604      fInfo_.SIM_uses_Dipoles = useDipole;
605      fInfo_.SIM_uses_Sticky = useSticky;
606 +    fInfo_.SIM_uses_StickyPower = useStickyPower;
607      fInfo_.SIM_uses_GayBerne = useGayBerne;
608      fInfo_.SIM_uses_EAM = useEAM;
609      fInfo_.SIM_uses_Shapes = useShape;
# Line 824 | Line 840 | namespace oopse {
840      }
841    }
842  
843 <  void SimInfo::setupCutoff() {
843 >  void SimInfo::setupCutoff() {    
844      getCutoff(rcut_, rsw_);    
845      double rnblist = rcut_ + 1; // skin of neighbor list
846  
847      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
848 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
848 >    
849 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
850 >    if (simParams_->haveCutoffPolicy()) {
851 >      std::string myPolicy = simParams_->getCutoffPolicy();
852 >      if (myPolicy == "MIX") {
853 >        cp = MIX_CUTOFF_POLICY;
854 >      } else {
855 >        if (myPolicy == "MAX") {
856 >          cp = MAX_CUTOFF_POLICY;
857 >        } else {
858 >          if (myPolicy == "TRADITIONAL") {            
859 >            cp = TRADITIONAL_CUTOFF_POLICY;
860 >          } else {
861 >            // throw error        
862 >            sprintf( painCave.errMsg,
863 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
864 >            painCave.isFatal = 1;
865 >            simError();
866 >          }    
867 >        }          
868 >      }
869 >    }
870 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
871    }
872  
873 +  void SimInfo::setupCoulombicCorrection( int isError ) {    
874 +    
875 +    int errorOut;
876 +    int cc =  NONE;
877 +    double alphaVal;
878 +
879 +    errorOut = isError;
880 +
881 +    if (simParams_->haveCoulombicCorrection()) {
882 +      std::string myCorrection = simParams_->getCoulombicCorrection();
883 +      if (myCorrection == "NONE") {
884 +        cc = NONE;
885 +      } else {
886 +        if (myCorrection == "UNDAMPED_WOLF") {
887 +          cc = UNDAMPED_WOLF;
888 +        } else {
889 +          if (myCorrection == "WOLF") {            
890 +            cc = WOLF;
891 +            if (!simParams_->haveDampingAlpha()) {
892 +              //throw error
893 +              sprintf( painCave.errMsg,
894 +                       "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Wolf Coulombic Correction.", simParams_->getDampingAlpha());
895 +              painCave.isFatal = 0;
896 +              simError();
897 +            }
898 +            alphaVal = simParams_->getDampingAlpha();
899 +          } else {
900 +            if (myCorrection == "REACTION_FIELD") {
901 +              cc = REACTION_FIELD;
902 +            } else {
903 +              // throw error        
904 +              sprintf( painCave.errMsg,
905 +                       "SimInfo error: Unknown coulombicCorrection. (Input file specified %s .)\n\tcoulombicCorrection must be one of: \"none\", \"undamped_wolf\", \"wolf\", or \"reaction_field\".", myCorrection.c_str() );
906 +              painCave.isFatal = 1;
907 +              simError();
908 +            }    
909 +          }          
910 +        }
911 +      }
912 +    }
913 +    initFortranFF( &fInfo_.SIM_uses_RF, &cc, &alphaVal, &errorOut );
914 +  }
915 +
916    void SimInfo::addProperty(GenericData* genData) {
917      properties_.addProperty(genData);  
918    }
# Line 939 | Line 1020 | namespace oopse {
1020  
1021      return o;
1022    }
1023 +  
1024 +  
1025 +   /*
1026 +   Returns center of mass and center of mass velocity in one function call.
1027 +   */
1028 +  
1029 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1030 +      SimInfo::MoleculeIterator i;
1031 +      Molecule* mol;
1032 +      
1033 +    
1034 +      double totalMass = 0.0;
1035 +    
1036  
1037 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1038 +         double mass = mol->getMass();
1039 +         totalMass += mass;
1040 +         com += mass * mol->getCom();
1041 +         comVel += mass * mol->getComVel();          
1042 +      }  
1043 +      
1044 + #ifdef IS_MPI
1045 +      double tmpMass = totalMass;
1046 +      Vector3d tmpCom(com);  
1047 +      Vector3d tmpComVel(comVel);
1048 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1049 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1050 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1051 + #endif
1052 +      
1053 +      com /= totalMass;
1054 +      comVel /= totalMass;
1055 +   }        
1056 +  
1057 +   /*
1058 +   Return intertia tensor for entire system and angular momentum Vector.
1059 +
1060 +
1061 +       [  Ixx -Ixy  -Ixz ]
1062 +  J =| -Iyx  Iyy  -Iyz |
1063 +       [ -Izx -Iyz   Izz ]
1064 +    */
1065 +
1066 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1067 +      
1068 +
1069 +      double xx = 0.0;
1070 +      double yy = 0.0;
1071 +      double zz = 0.0;
1072 +      double xy = 0.0;
1073 +      double xz = 0.0;
1074 +      double yz = 0.0;
1075 +      Vector3d com(0.0);
1076 +      Vector3d comVel(0.0);
1077 +      
1078 +      getComAll(com, comVel);
1079 +      
1080 +      SimInfo::MoleculeIterator i;
1081 +      Molecule* mol;
1082 +      
1083 +      Vector3d thisq(0.0);
1084 +      Vector3d thisv(0.0);
1085 +
1086 +      double thisMass = 0.0;
1087 +    
1088 +      
1089 +      
1090 +  
1091 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1092 +        
1093 +         thisq = mol->getCom()-com;
1094 +         thisv = mol->getComVel()-comVel;
1095 +         thisMass = mol->getMass();
1096 +         // Compute moment of intertia coefficients.
1097 +         xx += thisq[0]*thisq[0]*thisMass;
1098 +         yy += thisq[1]*thisq[1]*thisMass;
1099 +         zz += thisq[2]*thisq[2]*thisMass;
1100 +        
1101 +         // compute products of intertia
1102 +         xy += thisq[0]*thisq[1]*thisMass;
1103 +         xz += thisq[0]*thisq[2]*thisMass;
1104 +         yz += thisq[1]*thisq[2]*thisMass;
1105 +            
1106 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1107 +            
1108 +      }  
1109 +      
1110 +      
1111 +      inertiaTensor(0,0) = yy + zz;
1112 +      inertiaTensor(0,1) = -xy;
1113 +      inertiaTensor(0,2) = -xz;
1114 +      inertiaTensor(1,0) = -xy;
1115 +      inertiaTensor(1,1) = xx + zz;
1116 +      inertiaTensor(1,2) = -yz;
1117 +      inertiaTensor(2,0) = -xz;
1118 +      inertiaTensor(2,1) = -yz;
1119 +      inertiaTensor(2,2) = xx + yy;
1120 +      
1121 + #ifdef IS_MPI
1122 +      Mat3x3d tmpI(inertiaTensor);
1123 +      Vector3d tmpAngMom;
1124 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1125 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1126 + #endif
1127 +              
1128 +      return;
1129 +   }
1130 +
1131 +   //Returns the angular momentum of the system
1132 +   Vector3d SimInfo::getAngularMomentum(){
1133 +      
1134 +      Vector3d com(0.0);
1135 +      Vector3d comVel(0.0);
1136 +      Vector3d angularMomentum(0.0);
1137 +      
1138 +      getComAll(com,comVel);
1139 +      
1140 +      SimInfo::MoleculeIterator i;
1141 +      Molecule* mol;
1142 +      
1143 +      Vector3d thisr(0.0);
1144 +      Vector3d thisp(0.0);
1145 +      
1146 +      double thisMass;
1147 +      
1148 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1149 +        thisMass = mol->getMass();
1150 +        thisr = mol->getCom()-com;
1151 +        thisp = (mol->getComVel()-comVel)*thisMass;
1152 +        
1153 +        angularMomentum += cross( thisr, thisp );
1154 +        
1155 +      }  
1156 +      
1157 + #ifdef IS_MPI
1158 +      Vector3d tmpAngMom;
1159 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1160 + #endif
1161 +      
1162 +      return angularMomentum;
1163 +   }
1164 +  
1165 +  
1166   }//end namespace oopse
1167  

Diff Legend

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