ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimSetup.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 1108 by tim, Wed Apr 14 15:37:41 2004 UTC vs.
Revision 1163 by gezelter, Wed May 12 14:30:12 2004 UTC

# Line 147 | Line 147 | void SimSetup::createSim(void){
147    // make the output filenames
148  
149    makeOutNames();
150 –
151 –  if (globals->haveMinimizer())
152 –    // make minimizer
153 –    makeMinimizer();
154 –  else
155 –    // make the integrator
156 –    makeIntegrator();
150    
151   #ifdef IS_MPI
152    mpiSim->mpiRefresh();
# Line 162 | Line 155 | void SimSetup::createSim(void){
155    // initialize the Fortran
156  
157    initFortran();
158 +
159 +  if (globals->haveMinimizer())
160 +    // make minimizer
161 +    makeMinimizer();
162 +  else
163 +    // make the integrator
164 +    makeIntegrator();
165 +
166   }
167  
168  
# Line 182 | Line 183 | void SimSetup::makeMolecules(void){
183    BendStamp* currentBend;
184    TorsionStamp* currentTorsion;
185    RigidBodyStamp* currentRigidBody;
186 <
186 >  CutoffGroupStamp* currentCutoffGroup;
187 >  CutoffGroup* myCutoffGroup;
188 >  
189    bond_pair* theBonds;
190    bend_set* theBends;
191    torsion_set* theTorsions;
# Line 215 | Line 218 | void SimSetup::makeMolecules(void){
218        molInfo.nBends = comp_stamps[stampID]->getNBends();
219        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
220        molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
221 +      molInfo.nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
222        
223        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
224  
# Line 262 | Line 266 | void SimSetup::makeMolecules(void){
266          else{
267  
268            molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
269 +
270          }
271  
272          molInfo.myAtoms[j]->setType(currentAtom->getType());
268 –
273   #ifdef IS_MPI
274  
275          molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
# Line 409 | Line 413 | void SimSetup::makeMolecules(void){
413          info[k].excludes->addPair(exK, exL);
414        }
415  
416 +      
417 +      molInfo.myRigidBodies.clear();
418 +      
419        for (j = 0; j < molInfo.nRigidBodies; j++){
420  
421          currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
# Line 418 | Line 425 | void SimSetup::makeMolecules(void){
425  
426          myRB = new RigidBody();
427  
428 <        sprintf(rbName,"%s_RB_%s", molName, j);
428 >        sprintf(rbName,"%s_RB_%d", molName, j);
429          myRB->setType(rbName);
430          
431          for (rb1 = 0; rb1 < nMembers; rb1++) {
# Line 460 | Line 467 | void SimSetup::makeMolecules(void){
467              // used for the exclude list:
468              
469   #ifdef IS_MPI
470 <            exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
471 <            exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
470 >            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
471 >            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
472   #else
473 <            exI = tempI + 1;
474 <            exJ = tempJ + 1;
473 >            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
474 >            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
475   #endif
476              
477              info[k].excludes->addPair(exI, exJ);
# Line 477 | Line 484 | void SimSetup::makeMolecules(void){
484        }
485        
486  
487 +      //create cutoff group for molecule
488 +      molInfo.myCutoffGroups.clear();
489 +      for (j = 0; j < molInfo.nCutoffGroups; j++){
490 +
491 +        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
492 +        nMembers = currentCutoffGroup->getNMembers();
493 +
494 +        myCutoffGroup = new CutoffGroup();
495 +        
496 +        for (int cg = 0; cg < nMembers; cg++) {
497 +
498 +          // molI is atom numbering inside this molecule
499 +          molI = currentCutoffGroup->getMember(cg);    
500 +
501 +          // tempI is atom numbering on local processor
502 +          tempI = molI + atomOffset;
503 +
504 +          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
505 +        }
506 +
507 +        molInfo.myCutoffGroups.push_back(myCutoffGroup);
508 +      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
509 +      
510 +
511 +
512        // After this is all set up, scan through the atoms to
513        // see if they can be added to the integrableObjects:
514  
515 +      molInfo.myIntegrableObjects.clear();
516 +      
517 +
518        for (j = 0; j < molInfo.nAtoms; j++){
519  
520   #ifdef IS_MPI
# Line 528 | Line 563 | void SimSetup::makeMolecules(void){
563    sprintf(checkPointMsg, "all molecules initialized succesfully");
564    MPIcheckPoint();
565   #endif // is_mpi
531 –
532 –  // clean up the forcefield
533 –
534 –  if (!globals->haveLJrcut()){
535 –
536 –    the_ff->calcRcut();
537 –
538 –  } else {
539 –    
540 –    the_ff->setRcut( globals->getLJrcut() );
541 –  }
566  
543 –  the_ff->cleanMe();
567   }
568  
569   void SimSetup::initFromBass(void){
# Line 827 | Line 850 | void SimSetup::gatherInfo(void){
850    }
851  
852    //check whether sample time, status time, thermal time and reset time are divisble by dt
853 <  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
853 >  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
854      sprintf(painCave.errMsg,
855              "Sample time is not divisible by dt.\n"
856              "\tThis will result in samples that are not uniformly\n"
# Line 837 | Line 860 | void SimSetup::gatherInfo(void){
860      simError();    
861    }
862  
863 <  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
863 >  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
864      sprintf(painCave.errMsg,
865              "Status time is not divisible by dt.\n"
866              "\tThis will result in status reports that are not uniformly\n"
# Line 873 | Line 896 | void SimSetup::gatherInfo(void){
896      if (globals->haveSampleTime()){
897        info[i].sampleTime = globals->getSampleTime();
898        info[i].statusTime = info[i].sampleTime;
876 –      info[i].thermalTime = info[i].sampleTime;
899      }
900      else{
901        info[i].sampleTime = globals->getRunTime();
902        info[i].statusTime = info[i].sampleTime;
881 –      info[i].thermalTime = info[i].sampleTime;
903      }
904  
905      if (globals->haveStatusTime()){
# Line 887 | Line 908 | void SimSetup::gatherInfo(void){
908  
909      if (globals->haveThermalTime()){
910        info[i].thermalTime = globals->getThermalTime();
911 +    } else {
912 +      info[i].thermalTime = globals->getRunTime();
913      }
914  
915      info[i].resetIntegrator = 0;
# Line 957 | Line 980 | void SimSetup::finalInfoCheck(void){
980   void SimSetup::finalInfoCheck(void){
981    int index;
982    int usesDipoles;
983 +  int usesCharges;
984    int i;
985  
986    for (i = 0; i < nInfo; i++){
# Line 968 | Line 992 | void SimSetup::finalInfoCheck(void){
992        usesDipoles = (info[i].atoms[index])->hasDipole();
993        index++;
994      }
995 <
995 >    index = 0;
996 >    usesCharges = 0;
997 >    while ((index < info[i].n_atoms) && !usesCharges){
998 >      usesCharges= (info[i].atoms[index])->hasCharge();
999 >      index++;
1000 >    }
1001   #ifdef IS_MPI
1002      int myUse = usesDipoles;
1003      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1004   #endif //is_mpi
1005  
1006 <    double theEcr, theEst;
1006 >    double theRcut, theRsw;
1007  
1008 <    if (globals->getUseRF()){
1009 <      info[i].useReactionField = 1;
1008 >    if (globals->haveRcut()) {
1009 >      theRcut = globals->getRcut();
1010  
1011 <      if (!globals->haveECR()){
1012 <        sprintf(painCave.errMsg,
1013 <                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1011 >      if (globals->haveRsw())
1012 >        theRsw = globals->getRsw();
1013 >      else
1014 >        theRsw = theRcut;
1015 >      
1016 >      info[i].setDefaultRcut(theRcut, theRsw);
1017 >
1018 >    } else {
1019 >      
1020 >      the_ff->calcRcut();
1021 >      theRcut = info[i].getRcut();
1022 >
1023 >      if (globals->haveRsw())
1024 >        theRsw = globals->getRsw();
1025 >      else
1026 >        theRsw = theRcut;
1027 >      
1028 >      info[i].setDefaultRcut(theRcut, theRsw);
1029 >    }
1030 >
1031 >    if (globals->getUseRF()){
1032 >      info[i].useReactionField = 1;
1033 >      
1034 >      if (!globals->haveRcut()){
1035 >        sprintf(painCave.errMsg,
1036 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1037                  "\tOOPSE will use a default value of 15.0 angstroms"
1038 <                "\tfor the electrostaticCutoffRadius.\n");
1038 >                "\tfor the cutoffRadius.\n");
1039          painCave.isFatal = 0;
1040          simError();
1041 <        theEcr = 15.0;
1041 >        theRcut = 15.0;
1042        }
1043        else{
1044 <        theEcr = globals->getECR();
1044 >        theRcut = globals->getRcut();
1045        }
1046  
1047 <      if (!globals->haveEST()){
1047 >      if (!globals->haveRsw()){
1048          sprintf(painCave.errMsg,
1049 <                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1049 >                "SimSetup Warning: No value was set for switchingRadius.\n"
1050                  "\tOOPSE will use a default value of\n"
1051 <                "\t0.05 * electrostaticCutoffRadius\n"
1000 <                "\tfor the electrostaticSkinThickness\n");
1051 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
1052          painCave.isFatal = 0;
1053          simError();
1054 <        theEst = 0.05 * theEcr;
1054 >        theRsw = 0.95 * theRcut;
1055        }
1056        else{
1057 <        theEst = globals->getEST();
1057 >        theRsw = globals->getRsw();
1058        }
1059  
1060 <      info[i].setDefaultEcr(theEcr, theEst);
1060 >      info[i].setDefaultRcut(theRcut, theRsw);
1061  
1062        if (!globals->haveDielectric()){
1063          sprintf(painCave.errMsg,
# Line 1019 | Line 1070 | void SimSetup::finalInfoCheck(void){
1070        info[i].dielectric = globals->getDielectric();
1071      }
1072      else{
1073 <      if (usesDipoles){
1074 <        if (!globals->haveECR()){
1073 >      if (usesDipoles || usesCharges){
1074 >
1075 >        if (!globals->haveRcut()){
1076            sprintf(painCave.errMsg,
1077 <                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1077 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1078                    "\tOOPSE will use a default value of 15.0 angstroms"
1079 <                  "\tfor the electrostaticCutoffRadius.\n");
1080 <          painCave.isFatal = 0;
1081 <          simError();
1082 <          theEcr = 15.0;
1083 <        }
1084 <        else{
1085 <          theEcr = globals->getECR();
1079 >                  "\tfor the cutoffRadius.\n");
1080 >          painCave.isFatal = 0;
1081 >          simError();
1082 >          theRcut = 15.0;
1083 >      }
1084 >        else{
1085 >          theRcut = globals->getRcut();
1086          }
1087 <        
1088 <        if (!globals->haveEST()){
1087 >        
1088 >        if (!globals->haveRsw()){
1089            sprintf(painCave.errMsg,
1090 <                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1090 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1091                    "\tOOPSE will use a default value of\n"
1092 <                  "\t0.05 * electrostaticCutoffRadius\n"
1041 <                  "\tfor the electrostaticSkinThickness\n");
1092 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1093            painCave.isFatal = 0;
1094            simError();
1095 <          theEst = 0.05 * theEcr;
1095 >          theRsw = 0.95 * theRcut;
1096          }
1097          else{
1098 <          theEst = globals->getEST();
1098 >          theRsw = globals->getRsw();
1099          }
1100 +        
1101 +        info[i].setDefaultRcut(theRcut, theRsw);
1102          
1050 –        info[i].setDefaultEcr(theEcr, theEst);
1103        }
1104      }
1105    }
# Line 1055 | Line 1107 | void SimSetup::finalInfoCheck(void){
1107    strcpy(checkPointMsg, "post processing checks out");
1108    MPIcheckPoint();
1109   #endif // is_mpi
1110 +
1111 +  // clean up the forcefield
1112 +  the_ff->cleanMe();
1113   }
1114    
1115   void SimSetup::initSystemCoords(void){
# Line 1269 | Line 1324 | void SimSetup::compList(void){
1324    LinkedMolStamp* headStamp = new LinkedMolStamp();
1325    LinkedMolStamp* currentStamp = NULL;
1326    comp_stamps = new MoleculeStamp * [n_components];
1327 +  bool haveCutoffGroups;
1328  
1329 +  haveCutoffGroups = false;
1330 +  
1331    // make an array of molecule stamps that match the components used.
1332    // also extract the used stamps out into a separate linked list
1333  
# Line 1304 | Line 1362 | void SimSetup::compList(void){
1362        headStamp->add(currentStamp);
1363        comp_stamps[i] = headStamp->match(id);
1364      }
1365 +
1366 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1367 +      haveCutoffGroups = true;    
1368    }
1369 +    
1370 +  for (i = 0; i < nInfo; i++)
1371 +    info[i].haveCutoffGroups = haveCutoffGroups;
1372  
1373   #ifdef IS_MPI
1374    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1352 | Line 1416 | void SimSetup::mpiMolDivide(void){
1416    int localMol, allMol;
1417    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1418    int local_rigid;
1355 –  vector<int> globalAtomIndex;
1419    vector<int> globalMolIndex;
1420  
1421    mpiSim = new mpiSimulation(info);
1422  
1423    mpiSim->divideLabor();
1424    globalAtomIndex = mpiSim->getGlobalAtomIndex();
1425 <  globalMolIndex = mpiSim->getGlobalMolIndex();
1425 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1426  
1427    // set up the local variables
1428  
# Line 1480 | Line 1543 | void SimSetup::makeSysArrays(void){
1543      info[l].atoms = the_atoms;
1544      info[l].molecules = the_molecules;
1545      info[l].nGlobalExcludes = 0;
1546 <
1546 >    
1547      the_ff->setSimInfo(info);
1548    }
1549   }

Diff Legend

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