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 705 by chrisfen, Tue Nov 1 19:14:27 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 517 | Line 526 | namespace oopse {
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 >    int useDW;
532 >    std::string myMethod;
533 >
534 >    // set the useRF logical
535 >    useRF = 0;
536 >    useDW = 0;
537  
538 +
539 +    if (simParams_->haveElectrostaticSummationMethod()) {
540 +      std::string myMethod = simParams_->getElectrostaticSummationMethod();
541 +      toUpper(myMethod);
542 +      if (myMethod == "REACTION_FIELD") {
543 +        useRF=1;
544 +      } else {
545 +        if (myMethod == "DAMPED_WOLF") {
546 +          useDW = 1;
547 +        }
548 +      }
549 +    }
550 +
551      //loop over all of the atom types
552      for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
553        useLennardJones |= (*i)->isLennardJones();
# Line 583 | Line 611 | namespace oopse {
611  
612      temp = useRF;
613      MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
614 <    
614 >
615 >    temp = useDW;
616 >    MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
617 >
618   #endif
619  
620      fInfo_.SIM_uses_PBC = usePBC;    
# Line 599 | Line 630 | namespace oopse {
630      fInfo_.SIM_uses_Shapes = useShape;
631      fInfo_.SIM_uses_FLARB = useFLARB;
632      fInfo_.SIM_uses_RF = useRF;
633 +    fInfo_.SIM_uses_DampedWolf = useDW;
634  
635 <    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
636 <
635 >    if( myMethod == "REACTION_FIELD") {
636 >      
637        if (simParams_->haveDielectric()) {
638          fInfo_.dielect = simParams_->getDielectric();
639        } else {
# Line 611 | Line 643 | namespace oopse {
643                  "\tsetting a dielectric constant!\n");
644          painCave.isFatal = 1;
645          simError();
646 <      }
615 <        
616 <    } else {
617 <      fInfo_.dielect = 0.0;
646 >      }      
647      }
648  
649    }
# Line 650 | Line 679 | namespace oopse {
679  
680          totalMass = cg->getMass();
681          for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
682 <          mfact.push_back(atom->getMass()/totalMass);
682 >          // Check for massless groups - set mfact to 1 if true
683 >          if (totalMass != 0)
684 >            mfact.push_back(atom->getMass()/totalMass);
685 >          else
686 >            mfact.push_back( 1.0 );
687          }
688  
689        }      
# Line 786 | Line 819 | namespace oopse {
819      
820      if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
821          
822 <      if (!simParams_->haveRcut()){
822 >      if (!simParams_->haveCutoffRadius()){
823          sprintf(painCave.errMsg,
824                  "SimCreator Warning: No value was set for the cutoffRadius.\n"
825                  "\tOOPSE will use a default value of 15.0 angstroms"
# Line 795 | Line 828 | namespace oopse {
828          simError();
829          rcut = 15.0;
830        } else{
831 <        rcut = simParams_->getRcut();
831 >        rcut = simParams_->getCutoffRadius();
832        }
833  
834 <      if (!simParams_->haveRsw()){
834 >      if (!simParams_->haveSwitchingRadius()){
835          sprintf(painCave.errMsg,
836                  "SimCreator Warning: No value was set for switchingRadius.\n"
837                  "\tOOPSE will use a default value of\n"
838 <                "\t0.95 * cutoffRadius for the switchingRadius\n");
838 >                "\t0.85 * cutoffRadius for the switchingRadius\n");
839          painCave.isFatal = 0;
840          simError();
841 <        rsw = 0.95 * rcut;
841 >        rsw = 0.85 * rcut;
842        } else{
843 <        rsw = simParams_->getRsw();
843 >        rsw = simParams_->getSwitchingRadius();
844        }
845  
846      } else {
847        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
848        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
849          
850 <      if (simParams_->haveRcut()) {
851 <        rcut = simParams_->getRcut();
850 >      if (simParams_->haveCutoffRadius()) {
851 >        rcut = simParams_->getCutoffRadius();
852        } else {
853          //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
854          rcut = calcMaxCutoffRadius();
855        }
856  
857 <      if (simParams_->haveRsw()) {
858 <        rsw  = simParams_->getRsw();
857 >      if (simParams_->haveSwitchingRadius()) {
858 >        rsw  = simParams_->getSwitchingRadius();
859        } else {
860          rsw = rcut;
861        }
# Line 830 | Line 863 | namespace oopse {
863      }
864    }
865  
866 <  void SimInfo::setupCutoff() {
866 >  void SimInfo::setupCutoff() {    
867      getCutoff(rcut_, rsw_);    
868      double rnblist = rcut_ + 1; // skin of neighbor list
869  
870      //Pass these cutoff radius etc. to fortran. This function should be called once and only once
871 <    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
871 >    
872 >    int cp =  TRADITIONAL_CUTOFF_POLICY;
873 >    if (simParams_->haveCutoffPolicy()) {
874 >      std::string myPolicy = simParams_->getCutoffPolicy();
875 >      toUpper(myPolicy);
876 >      if (myPolicy == "MIX") {
877 >        cp = MIX_CUTOFF_POLICY;
878 >      } else {
879 >        if (myPolicy == "MAX") {
880 >          cp = MAX_CUTOFF_POLICY;
881 >        } else {
882 >          if (myPolicy == "TRADITIONAL") {            
883 >            cp = TRADITIONAL_CUTOFF_POLICY;
884 >          } else {
885 >            // throw error        
886 >            sprintf( painCave.errMsg,
887 >                     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
888 >            painCave.isFatal = 1;
889 >            simError();
890 >          }    
891 >        }          
892 >      }
893 >    }
894 >
895 >
896 >    if (simParams_->haveSkinThickness()) {
897 >      double skinThickness = simParams_->getSkinThickness();
898 >    }
899 >
900 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
901 >    // also send cutoff notification to electrostatics
902 >    setElectrostaticCutoffRadius(&rcut_, &rsw_);
903 >  }
904 >
905 >  void SimInfo::setupElectrostaticSummationMethod( int isError ) {    
906 >    
907 >    int errorOut;
908 >    int esm =  NONE;
909 >    double alphaVal;
910 >    double dielectric;
911 >
912 >    errorOut = isError;
913 >    alphaVal = simParams_->getDampingAlpha();
914 >    dielectric = simParams_->getDielectric();
915 >
916 >    if (simParams_->haveElectrostaticSummationMethod()) {
917 >      std::string myMethod = simParams_->getElectrostaticSummationMethod();
918 >      toUpper(myMethod);
919 >      if (myMethod == "NONE") {
920 >        esm = NONE;
921 >      } else {
922 >        if (myMethod == "UNDAMPED_WOLF") {
923 >          esm = UNDAMPED_WOLF;
924 >        } else {
925 >          if (myMethod == "DAMPED_WOLF") {            
926 >            esm = DAMPED_WOLF;
927 >            if (!simParams_->haveDampingAlpha()) {
928 >              //throw error
929 >              sprintf( painCave.errMsg,
930 >                       "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);
931 >              painCave.isFatal = 0;
932 >              simError();
933 >            }
934 >          } else {
935 >            if (myMethod == "REACTION_FIELD") {      
936 >              esm = REACTION_FIELD;
937 >            } else {
938 >              // throw error        
939 >              sprintf( painCave.errMsg,
940 >                       "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() );
941 >              painCave.isFatal = 1;
942 >              simError();
943 >            }    
944 >          }          
945 >        }
946 >      }
947 >    }
948 >    // let's pass some summation method variables to fortran
949 >    setElectrostaticSummationMethod( &esm );
950 >    setDampedWolfAlpha( &alphaVal );
951 >    setReactionFieldDielectric( &dielectric );
952 >    initFortranFF( &esm, &errorOut );
953    }
954  
955    void SimInfo::addProperty(GenericData* genData) {
# Line 945 | Line 1059 | namespace oopse {
1059  
1060      return o;
1061    }
1062 +  
1063 +  
1064 +   /*
1065 +   Returns center of mass and center of mass velocity in one function call.
1066 +   */
1067 +  
1068 +   void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1069 +      SimInfo::MoleculeIterator i;
1070 +      Molecule* mol;
1071 +      
1072 +    
1073 +      double totalMass = 0.0;
1074 +    
1075  
1076 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1077 +         double mass = mol->getMass();
1078 +         totalMass += mass;
1079 +         com += mass * mol->getCom();
1080 +         comVel += mass * mol->getComVel();          
1081 +      }  
1082 +      
1083 + #ifdef IS_MPI
1084 +      double tmpMass = totalMass;
1085 +      Vector3d tmpCom(com);  
1086 +      Vector3d tmpComVel(comVel);
1087 +      MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1088 +      MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1089 +      MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1090 + #endif
1091 +      
1092 +      com /= totalMass;
1093 +      comVel /= totalMass;
1094 +   }        
1095 +  
1096 +   /*
1097 +   Return intertia tensor for entire system and angular momentum Vector.
1098 +
1099 +
1100 +       [  Ixx -Ixy  -Ixz ]
1101 +  J =| -Iyx  Iyy  -Iyz |
1102 +       [ -Izx -Iyz   Izz ]
1103 +    */
1104 +
1105 +   void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1106 +      
1107 +
1108 +      double xx = 0.0;
1109 +      double yy = 0.0;
1110 +      double zz = 0.0;
1111 +      double xy = 0.0;
1112 +      double xz = 0.0;
1113 +      double yz = 0.0;
1114 +      Vector3d com(0.0);
1115 +      Vector3d comVel(0.0);
1116 +      
1117 +      getComAll(com, comVel);
1118 +      
1119 +      SimInfo::MoleculeIterator i;
1120 +      Molecule* mol;
1121 +      
1122 +      Vector3d thisq(0.0);
1123 +      Vector3d thisv(0.0);
1124 +
1125 +      double thisMass = 0.0;
1126 +    
1127 +      
1128 +      
1129 +  
1130 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1131 +        
1132 +         thisq = mol->getCom()-com;
1133 +         thisv = mol->getComVel()-comVel;
1134 +         thisMass = mol->getMass();
1135 +         // Compute moment of intertia coefficients.
1136 +         xx += thisq[0]*thisq[0]*thisMass;
1137 +         yy += thisq[1]*thisq[1]*thisMass;
1138 +         zz += thisq[2]*thisq[2]*thisMass;
1139 +        
1140 +         // compute products of intertia
1141 +         xy += thisq[0]*thisq[1]*thisMass;
1142 +         xz += thisq[0]*thisq[2]*thisMass;
1143 +         yz += thisq[1]*thisq[2]*thisMass;
1144 +            
1145 +         angularMomentum += cross( thisq, thisv ) * thisMass;
1146 +            
1147 +      }  
1148 +      
1149 +      
1150 +      inertiaTensor(0,0) = yy + zz;
1151 +      inertiaTensor(0,1) = -xy;
1152 +      inertiaTensor(0,2) = -xz;
1153 +      inertiaTensor(1,0) = -xy;
1154 +      inertiaTensor(1,1) = xx + zz;
1155 +      inertiaTensor(1,2) = -yz;
1156 +      inertiaTensor(2,0) = -xz;
1157 +      inertiaTensor(2,1) = -yz;
1158 +      inertiaTensor(2,2) = xx + yy;
1159 +      
1160 + #ifdef IS_MPI
1161 +      Mat3x3d tmpI(inertiaTensor);
1162 +      Vector3d tmpAngMom;
1163 +      MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1164 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1165 + #endif
1166 +              
1167 +      return;
1168 +   }
1169 +
1170 +   //Returns the angular momentum of the system
1171 +   Vector3d SimInfo::getAngularMomentum(){
1172 +      
1173 +      Vector3d com(0.0);
1174 +      Vector3d comVel(0.0);
1175 +      Vector3d angularMomentum(0.0);
1176 +      
1177 +      getComAll(com,comVel);
1178 +      
1179 +      SimInfo::MoleculeIterator i;
1180 +      Molecule* mol;
1181 +      
1182 +      Vector3d thisr(0.0);
1183 +      Vector3d thisp(0.0);
1184 +      
1185 +      double thisMass;
1186 +      
1187 +      for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {        
1188 +        thisMass = mol->getMass();
1189 +        thisr = mol->getCom()-com;
1190 +        thisp = (mol->getComVel()-comVel)*thisMass;
1191 +        
1192 +        angularMomentum += cross( thisr, thisp );
1193 +        
1194 +      }  
1195 +      
1196 + #ifdef IS_MPI
1197 +      Vector3d tmpAngMom;
1198 +      MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1199 + #endif
1200 +      
1201 +      return angularMomentum;
1202 +   }
1203 +  
1204 +  
1205   }//end namespace oopse
1206  

Diff Legend

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