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 682 by chrisfen, Tue Oct 18 15:01:42 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 80 | Line 83 | namespace oopse {
83        MoleculeStamp* molStamp;
84        int nMolWithSameStamp;
85        int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
86 <      int nGroups = 0;          //total cutoff groups defined in meta-data file
86 >      int nGroups = 0;      //total cutoff groups defined in meta-data file
87        CutoffGroupStamp* cgStamp;    
88        RigidBodyStamp* rbStamp;
89        int nRigidAtoms = 0;
# Line 105 | Line 108 | namespace oopse {
108          }
109  
110          nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
111 +
112          nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
113  
114          //calculate atoms in rigid bodies
# Line 121 | Line 125 | namespace oopse {
125          
126        }
127  
128 <      //every free atom (atom does not belong to cutoff groups) is a cutoff group
129 <      //therefore the total number of cutoff groups in the system is equal to
130 <      //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
131 <      //file plus the number of cutoff groups defined in meta-data file
128 >      //every free atom (atom does not belong to cutoff groups) is a cutoff
129 >      //group therefore the total number of cutoff groups in the system is
130 >      //equal to the total number of atoms minus number of atoms belong to
131 >      //cutoff group defined in meta-data file plus the number of cutoff
132 >      //groups defined in meta-data file
133        nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
134  
135 <      //every free atom (atom does not belong to rigid bodies) is an integrable object
136 <      //therefore the total number of  integrable objects in the system is equal to
137 <      //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
138 <      //file plus the number of  rigid bodies defined in meta-data file
139 <      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
140 <
135 >      //every free atom (atom does not belong to rigid bodies) is an
136 >      //integrable object therefore the total number of integrable objects
137 >      //in the system is equal to the total number of atoms minus number of
138 >      //atoms belong to rigid body defined in meta-data file plus the number
139 >      //of rigid bodies defined in meta-data file
140 >      nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
141 >                                                + nGlobalRigidBodies_;
142 >  
143        nGlobalMols_ = molStampIds_.size();
144  
145   #ifdef IS_MPI    
# Line 462 | Line 469 | namespace oopse {
469      //setup fortran force field
470      /** @deprecate */    
471      int isError = 0;
472 <    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
472 >    
473 >    setupElectrostaticSummationMethod( isError );
474 >
475      if(isError){
476        sprintf( painCave.errMsg,
477                 "ForceField error: There was an error initializing the forceField in fortran.\n" );
# Line 511 | Line 520 | namespace oopse {
520      int useDipole = 0;
521      int useGayBerne = 0;
522      int useSticky = 0;
523 +    int useStickyPower = 0;
524      int useShape = 0;
525      int useFLARB = 0; //it is not in AtomType yet
526      int useDirectionalAtom = 0;    
527      int useElectrostatics = 0;
528      //usePBC and useRF are from simParams
529 <    int usePBC = simParams_->getPBC();
530 <    int useRF = simParams_->getUseRF();
529 >    int usePBC = simParams_->getUsePeriodicBoundaryConditions();
530 >    int useRF;
531 >    std::string myMethod;
532  
533 +    // set the useRF logical
534 +    useRF = 0;
535 +    if (simParams_->haveElectrostaticSummationMethod()) {
536 +        myMethod = simParams_->getElectrostaticSummationMethod();
537 +        if (myMethod == "REACTION_FIELD")
538 +             useRF = 1;
539 +    }
540 +
541      //loop over all of the atom types
542      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
543        useLennardJones |= (*i)->isLennardJones();
# Line 529 | Line 548 | namespace oopse {
548        useDipole |= (*i)->isDipole();
549        useGayBerne |= (*i)->isGayBerne();
550        useSticky |= (*i)->isSticky();
551 +      useStickyPower |= (*i)->isStickyPower();
552        useShape |= (*i)->isShape();
553      }
554  
555 <    if (useSticky || useDipole || useGayBerne || useShape) {
555 >    if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
556        useDirectionalAtom = 1;
557      }
558  
# Line 564 | Line 584 | namespace oopse {
584      temp = useSticky;
585      MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
586  
587 +    temp = useStickyPower;
588 +    MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
589 +    
590      temp = useGayBerne;
591      MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
592  
# Line 578 | Line 601 | namespace oopse {
601  
602      temp = useRF;
603      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
604 <    
604 >
605   #endif
606  
607      fInfo_.SIM_uses_PBC = usePBC;    
# Line 588 | Line 611 | namespace oopse {
611      fInfo_.SIM_uses_Charges = useCharge;
612      fInfo_.SIM_uses_Dipoles = useDipole;
613      fInfo_.SIM_uses_Sticky = useSticky;
614 +    fInfo_.SIM_uses_StickyPower = useStickyPower;
615      fInfo_.SIM_uses_GayBerne = useGayBerne;
616      fInfo_.SIM_uses_EAM = useEAM;
617      fInfo_.SIM_uses_Shapes = useShape;
618      fInfo_.SIM_uses_FLARB = useFLARB;
619      fInfo_.SIM_uses_RF = useRF;
620  
621 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
621 >    if( fInfo_.SIM_uses_Dipoles && myMethod == "REACTION_FIELD") {
622  
623        if (simParams_->haveDielectric()) {
624          fInfo_.dielect = simParams_->getDielectric();
# Line 644 | Line 668 | namespace oopse {
668  
669          totalMass = cg->getMass();
670          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
671 <          mfact.push_back(atom->getMass()/totalMass);
671 >          // Check for massless groups - set mfact to 1 if true
672 >          if (totalMass != 0)
673 >            mfact.push_back(atom->getMass()/totalMass);
674 >          else
675 >            mfact.push_back( 1.0 );
676          }
677  
678        }      
# Line 780 | Line 808 | namespace oopse {
808      
809      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
810          
811 <      if (!simParams_->haveRcut()){
811 >      if (!simParams_->haveCutoffRadius()){
812          sprintf(painCave.errMsg,
813                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
814                  "\tOOPSE will use a default value of 15.0 angstroms"
# Line 789 | Line 817 | namespace oopse {
817          simError();
818          rcut = 15.0;
819        } else{
820 <        rcut = simParams_->getRcut();
820 >        rcut = simParams_->getCutoffRadius();
821        }
822  
823 <      if (!simParams_->haveRsw()){
823 >      if (!simParams_->haveSwitchingRadius()){
824          sprintf(painCave.errMsg,
825                  "SimCreator Warning: No value was set for switchingRadius.\n"
826                  "\tOOPSE will use a default value of\n"
# Line 801 | Line 829 | namespace oopse {
829          simError();
830          rsw = 0.95 * rcut;
831        } else{
832 <        rsw = simParams_->getRsw();
832 >        rsw = simParams_->getSwitchingRadius();
833        }
834  
835      } else {
836        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
837        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
838          
839 <      if (simParams_->haveRcut()) {
840 <        rcut = simParams_->getRcut();
839 >      if (simParams_->haveCutoffRadius()) {
840 >        rcut = simParams_->getCutoffRadius();
841        } else {
842          //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
843          rcut = calcMaxCutoffRadius();
844        }
845  
846 <      if (simParams_->haveRsw()) {
847 <        rsw  = simParams_->getRsw();
846 >      if (simParams_->haveSwitchingRadius()) {
847 >        rsw  = simParams_->getSwitchingRadius();
848        } else {
849          rsw = rcut;
850        }
# Line 824 | Line 852 | namespace oopse {
852      }
853    }
854  
855 <  void SimInfo::setupCutoff() {
855 >  void SimInfo::setupCutoff() {    
856      getCutoff(rcut_, rsw_);    
857      double rnblist = rcut_ + 1; // skin of neighbor list
858  
859      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
860 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
860 >    
861 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
862 >    if (simParams_->haveCutoffPolicy()) {
863 >      std::string myPolicy = simParams_->getCutoffPolicy();
864 >      toUpper(myPolicy);
865 >      if (myPolicy == "MIX") {
866 >        cp = MIX_CUTOFF_POLICY;
867 >      } else {
868 >        if (myPolicy == "MAX") {
869 >          cp = MAX_CUTOFF_POLICY;
870 >        } else {
871 >          if (myPolicy == "TRADITIONAL") {            
872 >            cp = TRADITIONAL_CUTOFF_POLICY;
873 >          } else {
874 >            // throw error        
875 >            sprintf( painCave.errMsg,
876 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
877 >            painCave.isFatal = 1;
878 >            simError();
879 >          }    
880 >        }          
881 >      }
882 >    }
883 >
884 >
885 >    if (simParams_->haveSkinThickness()) {
886 >      double skinThickness = simParams_->getSkinThickness();
887 >    }
888 >
889 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
890 >    // also send cutoff notification to electrostatics
891 >    setElectrostaticCutoffRadius(&rcut_, &rsw_);
892    }
893  
894 +  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
895 +    
896 +    int errorOut;
897 +    int esm =  NONE;
898 +    double alphaVal;
899 +    double dielectric;
900 +
901 +    errorOut = isError;
902 +    alphaVal = simParams_->getDampingAlpha();
903 +    dielectric = simParams_->getDielectric();
904 +
905 +    if (simParams_->haveElectrostaticSummationMethod()) {
906 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
907 +      toUpper(myMethod);
908 +      if (myMethod == "NONE") {
909 +        esm = NONE;
910 +      } else {
911 +        if (myMethod == "UNDAMPED_WOLF") {
912 +          esm = UNDAMPED_WOLF;
913 +        } else {
914 +          if (myMethod == "DAMPED_WOLF") {            
915 +            esm = DAMPED_WOLF;
916 +            if (!simParams_->haveDampingAlpha()) {
917 +              //throw error
918 +              sprintf( painCave.errMsg,
919 +                       "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);
920 +              painCave.isFatal = 0;
921 +              simError();
922 +            }
923 +          } else {
924 +            if (myMethod == "REACTION_FIELD") {
925 +              esm = REACTION_FIELD;
926 +            } else {
927 +              // throw error        
928 +              sprintf( painCave.errMsg,
929 +                       "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() );
930 +              painCave.isFatal = 1;
931 +              simError();
932 +            }    
933 +          }          
934 +        }
935 +      }
936 +    }
937 +    // let's pass some summation method variables to fortran
938 +    setElectrostaticSummationMethod( &esm );
939 +    setDampedWolfAlpha( &alphaVal );
940 +    setReactionFieldDielectric( &dielectric );
941 +    initFortranFF( &esm, &errorOut );
942 +  }
943 +
944    void SimInfo::addProperty(GenericData* genData) {
945      properties_.addProperty(genData);  
946    }
# Line 939 | Line 1048 | namespace oopse {
1048  
1049      return o;
1050    }
1051 +  
1052 +  
1053 +   /*
1054 +   Returns center of mass and center of mass velocity in one function call.
1055 +   */
1056 +  
1057 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1058 +      SimInfo::MoleculeIterator i;
1059 +      Molecule* mol;
1060 +      
1061 +    
1062 +      double totalMass = 0.0;
1063 +    
1064  
1065 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1066 +         double mass = mol->getMass();
1067 +         totalMass += mass;
1068 +         com += mass * mol->getCom();
1069 +         comVel += mass * mol->getComVel();          
1070 +      }  
1071 +      
1072 + #ifdef IS_MPI
1073 +      double tmpMass = totalMass;
1074 +      Vector3d tmpCom(com);  
1075 +      Vector3d tmpComVel(comVel);
1076 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1077 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1078 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1079 + #endif
1080 +      
1081 +      com /= totalMass;
1082 +      comVel /= totalMass;
1083 +   }        
1084 +  
1085 +   /*
1086 +   Return intertia tensor for entire system and angular momentum Vector.
1087 +
1088 +
1089 +       [  Ixx -Ixy  -Ixz ]
1090 +  J =| -Iyx  Iyy  -Iyz |
1091 +       [ -Izx -Iyz   Izz ]
1092 +    */
1093 +
1094 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1095 +      
1096 +
1097 +      double xx = 0.0;
1098 +      double yy = 0.0;
1099 +      double zz = 0.0;
1100 +      double xy = 0.0;
1101 +      double xz = 0.0;
1102 +      double yz = 0.0;
1103 +      Vector3d com(0.0);
1104 +      Vector3d comVel(0.0);
1105 +      
1106 +      getComAll(com, comVel);
1107 +      
1108 +      SimInfo::MoleculeIterator i;
1109 +      Molecule* mol;
1110 +      
1111 +      Vector3d thisq(0.0);
1112 +      Vector3d thisv(0.0);
1113 +
1114 +      double thisMass = 0.0;
1115 +    
1116 +      
1117 +      
1118 +  
1119 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1120 +        
1121 +         thisq = mol->getCom()-com;
1122 +         thisv = mol->getComVel()-comVel;
1123 +         thisMass = mol->getMass();
1124 +         // Compute moment of intertia coefficients.
1125 +         xx += thisq[0]*thisq[0]*thisMass;
1126 +         yy += thisq[1]*thisq[1]*thisMass;
1127 +         zz += thisq[2]*thisq[2]*thisMass;
1128 +        
1129 +         // compute products of intertia
1130 +         xy += thisq[0]*thisq[1]*thisMass;
1131 +         xz += thisq[0]*thisq[2]*thisMass;
1132 +         yz += thisq[1]*thisq[2]*thisMass;
1133 +            
1134 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1135 +            
1136 +      }  
1137 +      
1138 +      
1139 +      inertiaTensor(0,0) = yy + zz;
1140 +      inertiaTensor(0,1) = -xy;
1141 +      inertiaTensor(0,2) = -xz;
1142 +      inertiaTensor(1,0) = -xy;
1143 +      inertiaTensor(1,1) = xx + zz;
1144 +      inertiaTensor(1,2) = -yz;
1145 +      inertiaTensor(2,0) = -xz;
1146 +      inertiaTensor(2,1) = -yz;
1147 +      inertiaTensor(2,2) = xx + yy;
1148 +      
1149 + #ifdef IS_MPI
1150 +      Mat3x3d tmpI(inertiaTensor);
1151 +      Vector3d tmpAngMom;
1152 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1153 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1154 + #endif
1155 +              
1156 +      return;
1157 +   }
1158 +
1159 +   //Returns the angular momentum of the system
1160 +   Vector3d SimInfo::getAngularMomentum(){
1161 +      
1162 +      Vector3d com(0.0);
1163 +      Vector3d comVel(0.0);
1164 +      Vector3d angularMomentum(0.0);
1165 +      
1166 +      getComAll(com,comVel);
1167 +      
1168 +      SimInfo::MoleculeIterator i;
1169 +      Molecule* mol;
1170 +      
1171 +      Vector3d thisr(0.0);
1172 +      Vector3d thisp(0.0);
1173 +      
1174 +      double thisMass;
1175 +      
1176 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1177 +        thisMass = mol->getMass();
1178 +        thisr = mol->getCom()-com;
1179 +        thisp = (mol->getComVel()-comVel)*thisMass;
1180 +        
1181 +        angularMomentum += cross( thisr, thisp );
1182 +        
1183 +      }  
1184 +      
1185 + #ifdef IS_MPI
1186 +      Vector3d tmpAngMom;
1187 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1188 + #endif
1189 +      
1190 +      return angularMomentum;
1191 +   }
1192 +  
1193 +  
1194   }//end namespace oopse
1195  

Diff Legend

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