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 709 by chrisfen, Wed Nov 2 20:36:15 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/DarkSide/fScreeningMethod.h"
58   #include "UseTheForce/doForces_interface.h"
59 + #include "UseTheForce/DarkSide/electrostatic_interface.h"
60   #include "UseTheForce/notifyCutoffs_interface.h"
61   #include "utils/MemoryUtils.hpp"
62   #include "utils/simError.h"
# Line 80 | Line 84 | namespace oopse {
84        MoleculeStamp* molStamp;
85        int nMolWithSameStamp;
86        int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
87 <      int nGroups = 0;          //total cutoff groups defined in meta-data file
87 >      int nGroups = 0;      //total cutoff groups defined in meta-data file
88        CutoffGroupStamp* cgStamp;    
89        RigidBodyStamp* rbStamp;
90        int nRigidAtoms = 0;
# Line 105 | Line 109 | namespace oopse {
109          }
110  
111          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
112 +
113          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
114  
115          //calculate atoms in rigid bodies
# Line 121 | Line 126 | namespace oopse {
126          
127        }
128  
129 <      //every free atom (atom does not belong to cutoff groups) is a cutoff group
130 <      //therefore the total number of cutoff groups in the system is equal to
131 <      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
132 <      //file plus the number of cutoff groups defined in meta-data file
129 >      //every free atom (atom does not belong to cutoff groups) is a cutoff
130 >      //group therefore the total number of cutoff groups in the system is
131 >      //equal to the total number of atoms minus number of atoms belong to
132 >      //cutoff group defined in meta-data file plus the number of cutoff
133 >      //groups defined in meta-data file
134        nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
135  
136 <      //every free atom (atom does not belong to rigid bodies) is an integrable object
137 <      //therefore the total number of  integrable objects in the system is equal to
138 <      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
139 <      //file plus the number of  rigid bodies defined in meta-data file
140 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
141 <
136 >      //every free atom (atom does not belong to rigid bodies) is an
137 >      //integrable object therefore the total number of integrable objects
138 >      //in the system is equal to the total number of atoms minus number of
139 >      //atoms belong to rigid body defined in meta-data file plus the number
140 >      //of rigid bodies defined in meta-data file
141 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
142 >                                                + nGlobalRigidBodies_;
143 >  
144        nGlobalMols_ = molStampIds_.size();
145  
146   #ifdef IS_MPI    
# Line 462 | Line 470 | namespace oopse {
470      //setup fortran force field
471      /** @deprecate */    
472      int isError = 0;
473 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
473 >    
474 >    setupElectrostaticSummationMethod( isError );
475 >
476      if(isError){
477        sprintf( painCave.errMsg,
478                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 517 | Line 527 | namespace oopse {
527      int useDirectionalAtom = 0;    
528      int useElectrostatics = 0;
529      //usePBC and useRF are from simParams
530 <    int usePBC = simParams_->getPBC();
531 <    int useRF = simParams_->getUseRF();
530 >    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
531 >    int useRF;
532 >    int useDW;
533 >    std::string myMethod;
534  
535 +    // set the useRF logical
536 +    useRF = 0;
537 +    useDW = 0;
538 +
539 +
540 +    if (simParams_->haveElectrostaticSummationMethod()) {
541 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
542 +      toUpper(myMethod);
543 +      if (myMethod == "REACTION_FIELD") {
544 +        useRF=1;
545 +      } else {
546 +        if (myMethod == "SHIFTED_POTENTIAL") {
547 +          useDW = 1;
548 +        }
549 +      }
550 +    }
551 +
552      //loop over all of the atom types
553      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
554        useLennardJones |= (*i)->isLennardJones();
# Line 583 | Line 612 | namespace oopse {
612  
613      temp = useRF;
614      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
615 <    
615 >
616 >    temp = useDW;
617 >    MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
618 >
619   #endif
620  
621      fInfo_.SIM_uses_PBC = usePBC;    
# Line 599 | Line 631 | namespace oopse {
631      fInfo_.SIM_uses_Shapes = useShape;
632      fInfo_.SIM_uses_FLARB = useFLARB;
633      fInfo_.SIM_uses_RF = useRF;
634 +    fInfo_.SIM_uses_DampedWolf = useDW;
635  
636 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
637 <
636 >    if( myMethod == "REACTION_FIELD") {
637 >      
638        if (simParams_->haveDielectric()) {
639          fInfo_.dielect = simParams_->getDielectric();
640        } else {
# Line 611 | Line 644 | namespace oopse {
644                  "\tsetting a dielectric constant!\n");
645          painCave.isFatal = 1;
646          simError();
647 <      }
615 <        
616 <    } else {
617 <      fInfo_.dielect = 0.0;
647 >      }      
648      }
649  
650    }
# Line 650 | Line 680 | namespace oopse {
680  
681          totalMass = cg->getMass();
682          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
683 <          mfact.push_back(atom->getMass()/totalMass);
683 >          // Check for massless groups - set mfact to 1 if true
684 >          if (totalMass != 0)
685 >            mfact.push_back(atom->getMass()/totalMass);
686 >          else
687 >            mfact.push_back( 1.0 );
688          }
689  
690        }      
# Line 786 | Line 820 | namespace oopse {
820      
821      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
822          
823 <      if (!simParams_->haveRcut()){
823 >      if (!simParams_->haveCutoffRadius()){
824          sprintf(painCave.errMsg,
825                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
826                  "\tOOPSE will use a default value of 15.0 angstroms"
# Line 795 | Line 829 | namespace oopse {
829          simError();
830          rcut = 15.0;
831        } else{
832 <        rcut = simParams_->getRcut();
832 >        rcut = simParams_->getCutoffRadius();
833        }
834  
835 <      if (!simParams_->haveRsw()){
835 >      if (!simParams_->haveSwitchingRadius()){
836          sprintf(painCave.errMsg,
837                  "SimCreator Warning: No value was set for switchingRadius.\n"
838                  "\tOOPSE will use a default value of\n"
839 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
839 >                "\t0.85 * cutoffRadius for the switchingRadius\n");
840          painCave.isFatal = 0;
841          simError();
842 <        rsw = 0.95 * rcut;
842 >        rsw = 0.85 * rcut;
843        } else{
844 <        rsw = simParams_->getRsw();
844 >        rsw = simParams_->getSwitchingRadius();
845        }
846  
847      } else {
848        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
849        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
850          
851 <      if (simParams_->haveRcut()) {
852 <        rcut = simParams_->getRcut();
851 >      if (simParams_->haveCutoffRadius()) {
852 >        rcut = simParams_->getCutoffRadius();
853        } else {
854          //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
855          rcut = calcMaxCutoffRadius();
856        }
857  
858 <      if (simParams_->haveRsw()) {
859 <        rsw  = simParams_->getRsw();
858 >      if (simParams_->haveSwitchingRadius()) {
859 >        rsw  = simParams_->getSwitchingRadius();
860        } else {
861          rsw = rcut;
862        }
# Line 830 | Line 864 | namespace oopse {
864      }
865    }
866  
867 <  void SimInfo::setupCutoff() {
867 >  void SimInfo::setupCutoff() {    
868      getCutoff(rcut_, rsw_);    
869      double rnblist = rcut_ + 1; // skin of neighbor list
870  
871      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
872 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
872 >    
873 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
874 >    if (simParams_->haveCutoffPolicy()) {
875 >      std::string myPolicy = simParams_->getCutoffPolicy();
876 >      toUpper(myPolicy);
877 >      if (myPolicy == "MIX") {
878 >        cp = MIX_CUTOFF_POLICY;
879 >      } else {
880 >        if (myPolicy == "MAX") {
881 >          cp = MAX_CUTOFF_POLICY;
882 >        } else {
883 >          if (myPolicy == "TRADITIONAL") {            
884 >            cp = TRADITIONAL_CUTOFF_POLICY;
885 >          } else {
886 >            // throw error        
887 >            sprintf( painCave.errMsg,
888 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
889 >            painCave.isFatal = 1;
890 >            simError();
891 >          }    
892 >        }          
893 >      }
894 >    }
895 >
896 >
897 >    if (simParams_->haveSkinThickness()) {
898 >      double skinThickness = simParams_->getSkinThickness();
899 >    }
900 >
901 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
902 >    // also send cutoff notification to electrostatics
903 >    setElectrostaticCutoffRadius(&rcut_, &rsw_);
904    }
905  
906 +  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
907 +    
908 +    int errorOut;
909 +    int esm =  NONE;
910 +    int sm = UNDAMPED;
911 +    double alphaVal;
912 +    double dielectric;
913 +
914 +    errorOut = isError;
915 +    alphaVal = simParams_->getDampingAlpha();
916 +    dielectric = simParams_->getDielectric();
917 +
918 +    if (simParams_->haveElectrostaticSummationMethod()) {
919 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
920 +      toUpper(myMethod);
921 +      if (myMethod == "NONE") {
922 +        esm = NONE;
923 +      } else {
924 +        if (myMethod == "SWITCHING_FUNCTION") {
925 +          esm = SWITCHING_FUNCTION;
926 +        } else {
927 +          if (myMethod == "SHIFTED_POTENTIAL") {
928 +            esm = SHIFTED_POTENTIAL;
929 +          } else {
930 +            if (myMethod == "SHIFTED_FORCE") {            
931 +              esm = SHIFTED_FORCE;
932 +            } else {
933 +              if (myMethod == "REACTION_FIELD") {            
934 +                esm = REACTION_FIELD;
935 +              } else {
936 +                // throw error        
937 +                sprintf( painCave.errMsg,
938 +                         "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() );
939 +                painCave.isFatal = 1;
940 +                simError();
941 +              }    
942 +            }          
943 +          }
944 +        }
945 +      }
946 +    }
947 +    
948 +    if (simParams_->haveScreeningMethod()) {
949 +      std::string myScreen = simParams_->getScreeningMethod();
950 +      toUpper(myScreen);
951 +      if (myScreen == "UNDAMPED") {
952 +        sm = UNDAMPED;
953 +      } else {
954 +        if (myScreen == "DAMPED") {
955 +          sm = DAMPED;
956 +          if (!simParams_->haveDampingAlpha()) {
957 +            //throw error
958 +            sprintf( painCave.errMsg,
959 +                     "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal);
960 +            painCave.isFatal = 0;
961 +            simError();
962 +          } else {
963 +            // throw error        
964 +            sprintf( painCave.errMsg,
965 +                     "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() );
966 +            painCave.isFatal = 1;
967 +            simError();
968 +          }
969 +        }
970 +      }
971 +    }
972 +    // let's pass some summation method variables to fortran
973 +    setElectrostaticSummationMethod( &esm );
974 +    setScreeningMethod( &sm );
975 +    setDampingAlpha( &alphaVal );
976 +    setReactionFieldDielectric( &dielectric );
977 +    initFortranFF( &esm, &errorOut );
978 +  }
979 +
980    void SimInfo::addProperty(GenericData* genData) {
981      properties_.addProperty(genData);  
982    }
# Line 945 | Line 1084 | namespace oopse {
1084  
1085      return o;
1086    }
1087 +  
1088 +  
1089 +   /*
1090 +   Returns center of mass and center of mass velocity in one function call.
1091 +   */
1092 +  
1093 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1094 +      SimInfo::MoleculeIterator i;
1095 +      Molecule* mol;
1096 +      
1097 +    
1098 +      double totalMass = 0.0;
1099 +    
1100  
1101 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1102 +         double mass = mol->getMass();
1103 +         totalMass += mass;
1104 +         com += mass * mol->getCom();
1105 +         comVel += mass * mol->getComVel();          
1106 +      }  
1107 +      
1108 + #ifdef IS_MPI
1109 +      double tmpMass = totalMass;
1110 +      Vector3d tmpCom(com);  
1111 +      Vector3d tmpComVel(comVel);
1112 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1113 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1114 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1115 + #endif
1116 +      
1117 +      com /= totalMass;
1118 +      comVel /= totalMass;
1119 +   }        
1120 +  
1121 +   /*
1122 +   Return intertia tensor for entire system and angular momentum Vector.
1123 +
1124 +
1125 +       [  Ixx -Ixy  -Ixz ]
1126 +  J =| -Iyx  Iyy  -Iyz |
1127 +       [ -Izx -Iyz   Izz ]
1128 +    */
1129 +
1130 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1131 +      
1132 +
1133 +      double xx = 0.0;
1134 +      double yy = 0.0;
1135 +      double zz = 0.0;
1136 +      double xy = 0.0;
1137 +      double xz = 0.0;
1138 +      double yz = 0.0;
1139 +      Vector3d com(0.0);
1140 +      Vector3d comVel(0.0);
1141 +      
1142 +      getComAll(com, comVel);
1143 +      
1144 +      SimInfo::MoleculeIterator i;
1145 +      Molecule* mol;
1146 +      
1147 +      Vector3d thisq(0.0);
1148 +      Vector3d thisv(0.0);
1149 +
1150 +      double thisMass = 0.0;
1151 +    
1152 +      
1153 +      
1154 +  
1155 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1156 +        
1157 +         thisq = mol->getCom()-com;
1158 +         thisv = mol->getComVel()-comVel;
1159 +         thisMass = mol->getMass();
1160 +         // Compute moment of intertia coefficients.
1161 +         xx += thisq[0]*thisq[0]*thisMass;
1162 +         yy += thisq[1]*thisq[1]*thisMass;
1163 +         zz += thisq[2]*thisq[2]*thisMass;
1164 +        
1165 +         // compute products of intertia
1166 +         xy += thisq[0]*thisq[1]*thisMass;
1167 +         xz += thisq[0]*thisq[2]*thisMass;
1168 +         yz += thisq[1]*thisq[2]*thisMass;
1169 +            
1170 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1171 +            
1172 +      }  
1173 +      
1174 +      
1175 +      inertiaTensor(0,0) = yy + zz;
1176 +      inertiaTensor(0,1) = -xy;
1177 +      inertiaTensor(0,2) = -xz;
1178 +      inertiaTensor(1,0) = -xy;
1179 +      inertiaTensor(1,1) = xx + zz;
1180 +      inertiaTensor(1,2) = -yz;
1181 +      inertiaTensor(2,0) = -xz;
1182 +      inertiaTensor(2,1) = -yz;
1183 +      inertiaTensor(2,2) = xx + yy;
1184 +      
1185 + #ifdef IS_MPI
1186 +      Mat3x3d tmpI(inertiaTensor);
1187 +      Vector3d tmpAngMom;
1188 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1189 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1190 + #endif
1191 +              
1192 +      return;
1193 +   }
1194 +
1195 +   //Returns the angular momentum of the system
1196 +   Vector3d SimInfo::getAngularMomentum(){
1197 +      
1198 +      Vector3d com(0.0);
1199 +      Vector3d comVel(0.0);
1200 +      Vector3d angularMomentum(0.0);
1201 +      
1202 +      getComAll(com,comVel);
1203 +      
1204 +      SimInfo::MoleculeIterator i;
1205 +      Molecule* mol;
1206 +      
1207 +      Vector3d thisr(0.0);
1208 +      Vector3d thisp(0.0);
1209 +      
1210 +      double thisMass;
1211 +      
1212 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1213 +        thisMass = mol->getMass();
1214 +        thisr = mol->getCom()-com;
1215 +        thisp = (mol->getComVel()-comVel)*thisMass;
1216 +        
1217 +        angularMomentum += cross( thisr, thisp );
1218 +        
1219 +      }  
1220 +      
1221 + #ifdef IS_MPI
1222 +      Vector3d tmpAngMom;
1223 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1224 + #endif
1225 +      
1226 +      return angularMomentum;
1227 +   }
1228 +  
1229 +  
1230   }//end namespace oopse
1231  

Diff Legend

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