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 1116 by tim, Thu Apr 15 22:15:21 2004 UTC vs.
Revision 1174 by gezelter, Wed May 12 20:54:10 2004 UTC

# Line 183 | Line 183 | void SimSetup::makeMolecules(void){
183    BendStamp* currentBend;
184    TorsionStamp* currentTorsion;
185    RigidBodyStamp* currentRigidBody;
186 +  CutoffGroupStamp* currentCutoffGroup;
187 +  CutoffGroup* myCutoffGroup;
188 +  int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file
189 +  set<int> cutoffAtomSet; //atoms belong to  cutoffgroup defined at mdl file
190  
191    bond_pair* theBonds;
192    bend_set* theBends;
# Line 216 | Line 220 | void SimSetup::makeMolecules(void){
220        molInfo.nBends = comp_stamps[stampID]->getNBends();
221        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
222        molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
223 +
224 +      nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
225        
226        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
227  
228        if (molInfo.nBonds > 0)
229 <        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
229 >        molInfo.myBonds = new Bond*[molInfo.nBonds];
230        else
231          molInfo.myBonds = NULL;
232  
233        if (molInfo.nBends > 0)
234 <        molInfo.myBends = new (Bend *) [molInfo.nBends];
234 >        molInfo.myBends = new Bend*[molInfo.nBends];
235        else
236          molInfo.myBends = NULL;
237  
238        if (molInfo.nTorsions > 0)
239 <        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
239 >        molInfo.myTorsions = new Torsion *[molInfo.nTorsions];
240        else
241          molInfo.myTorsions = NULL;
242  
# Line 480 | Line 486 | void SimSetup::makeMolecules(void){
486          info[k].rigidBodies.push_back(myRB);
487        }
488        
489 +
490 +      //create cutoff group for molecule
491 +
492 +      cutoffAtomSet.clear();
493 +      molInfo.myCutoffGroups.clear();
494 +      
495 +      for (j = 0; j < nCutoffGroups; j++){
496 +
497 +        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
498 +        nMembers = currentCutoffGroup->getNMembers();
499 +
500 +        myCutoffGroup = new CutoffGroup();
501 +        
502 +        for (int cg = 0; cg < nMembers; cg++) {
503  
504 +          // molI is atom numbering inside this molecule
505 +          molI = currentCutoffGroup->getMember(cg);    
506 +
507 +          // tempI is atom numbering on local processor
508 +          tempI = molI + atomOffset;
509 +          
510 +          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
511 +
512 +          cutoffAtomSet.insert(tempI);
513 +        }
514 +
515 +        molInfo.myCutoffGroups.push_back(myCutoffGroup);
516 +      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
517 +
518 +      //creat a cutoff group for every atom  in current molecule which does not belong to cutoffgroup defined at mdl file
519 +
520 +      for(j = 0; j < molInfo.nAtoms; j++){
521 +
522 +        if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
523 +          myCutoffGroup = new CutoffGroup();
524 +          myCutoffGroup->addAtom(molInfo.myAtoms[j]);
525 +          molInfo.myCutoffGroups.push_back(myCutoffGroup);
526 +        }
527 +          
528 +      }
529 +
530 +              
531 +
532 +
533        // After this is all set up, scan through the atoms to
534        // see if they can be added to the integrableObjects:
535  
# Line 536 | Line 585 | void SimSetup::makeMolecules(void){
585    MPIcheckPoint();
586   #endif // is_mpi
587  
539 –  // clean up the forcefield
540 –
541 –  if (!globals->haveLJrcut()){
542 –
543 –    the_ff->calcRcut();
544 –
545 –  } else {
546 –    
547 –    the_ff->setRcut( globals->getLJrcut() );
548 –  }
549 –
550 –  the_ff->cleanMe();
588   }
589  
590   void SimSetup::initFromBass(void){
# Line 834 | Line 871 | void SimSetup::gatherInfo(void){
871    }
872  
873    //check whether sample time, status time, thermal time and reset time are divisble by dt
874 <  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
874 >  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
875      sprintf(painCave.errMsg,
876              "Sample time is not divisible by dt.\n"
877              "\tThis will result in samples that are not uniformly\n"
# Line 844 | Line 881 | void SimSetup::gatherInfo(void){
881      simError();    
882    }
883  
884 <  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
884 >  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
885      sprintf(painCave.errMsg,
886              "Status time is not divisible by dt.\n"
887              "\tThis will result in status reports that are not uniformly\n"
# Line 880 | Line 917 | void SimSetup::gatherInfo(void){
917      if (globals->haveSampleTime()){
918        info[i].sampleTime = globals->getSampleTime();
919        info[i].statusTime = info[i].sampleTime;
883 –      info[i].thermalTime = info[i].sampleTime;
920      }
921      else{
922        info[i].sampleTime = globals->getRunTime();
923        info[i].statusTime = info[i].sampleTime;
888 –      info[i].thermalTime = info[i].sampleTime;
924      }
925  
926      if (globals->haveStatusTime()){
# Line 894 | Line 929 | void SimSetup::gatherInfo(void){
929  
930      if (globals->haveThermalTime()){
931        info[i].thermalTime = globals->getThermalTime();
932 +    } else {
933 +      info[i].thermalTime = globals->getRunTime();
934      }
935  
936      info[i].resetIntegrator = 0;
# Line 987 | Line 1024 | void SimSetup::finalInfoCheck(void){
1024      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1025   #endif //is_mpi
1026  
1027 <    double theEcr, theEst;
1027 >    double theRcut, theRsw;
1028  
1029 +    if (globals->haveRcut()) {
1030 +      theRcut = globals->getRcut();
1031 +
1032 +      if (globals->haveRsw())
1033 +        theRsw = globals->getRsw();
1034 +      else
1035 +        theRsw = theRcut;
1036 +      
1037 +      info[i].setDefaultRcut(theRcut, theRsw);
1038 +
1039 +    } else {
1040 +      
1041 +      the_ff->calcRcut();
1042 +      theRcut = info[i].getRcut();
1043 +
1044 +      if (globals->haveRsw())
1045 +        theRsw = globals->getRsw();
1046 +      else
1047 +        theRsw = theRcut;
1048 +      
1049 +      info[i].setDefaultRcut(theRcut, theRsw);
1050 +    }
1051 +
1052      if (globals->getUseRF()){
1053        info[i].useReactionField = 1;
1054 <
1055 <      if (!globals->haveECR()){
1054 >      
1055 >      if (!globals->haveRcut()){
1056          sprintf(painCave.errMsg,
1057 <                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1057 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1058                  "\tOOPSE will use a default value of 15.0 angstroms"
1059 <                "\tfor the electrostaticCutoffRadius.\n");
1059 >                "\tfor the cutoffRadius.\n");
1060          painCave.isFatal = 0;
1061          simError();
1062 <        theEcr = 15.0;
1062 >        theRcut = 15.0;
1063        }
1064        else{
1065 <        theEcr = globals->getECR();
1065 >        theRcut = globals->getRcut();
1066        }
1067  
1068 <      if (!globals->haveEST()){
1068 >      if (!globals->haveRsw()){
1069          sprintf(painCave.errMsg,
1070 <                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1070 >                "SimSetup Warning: No value was set for switchingRadius.\n"
1071                  "\tOOPSE will use a default value of\n"
1072 <                "\t0.05 * electrostaticCutoffRadius\n"
1013 <                "\tfor the electrostaticSkinThickness\n");
1072 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
1073          painCave.isFatal = 0;
1074          simError();
1075 <        theEst = 0.05 * theEcr;
1075 >        theRsw = 0.95 * theRcut;
1076        }
1077        else{
1078 <        theEst = globals->getEST();
1078 >        theRsw = globals->getRsw();
1079        }
1080  
1081 <      info[i].setDefaultEcr(theEcr, theEst);
1081 >      info[i].setDefaultRcut(theRcut, theRsw);
1082  
1083        if (!globals->haveDielectric()){
1084          sprintf(painCave.errMsg,
# Line 1033 | Line 1092 | void SimSetup::finalInfoCheck(void){
1092      }
1093      else{
1094        if (usesDipoles || usesCharges){
1095 <        if (!globals->haveECR()){
1095 >
1096 >        if (!globals->haveRcut()){
1097            sprintf(painCave.errMsg,
1098 <                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1098 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1099                    "\tOOPSE will use a default value of 15.0 angstroms"
1100 <                  "\tfor the electrostaticCutoffRadius.\n");
1101 <          painCave.isFatal = 0;
1102 <          simError();
1103 <          theEcr = 15.0;
1104 <        }
1100 >                  "\tfor the cutoffRadius.\n");
1101 >          painCave.isFatal = 0;
1102 >          simError();
1103 >          theRcut = 15.0;
1104 >      }
1105          else{
1106 <          theEcr = globals->getECR();
1107 <        }
1108 <        
1109 <        if (!globals->haveEST()){
1106 >          theRcut = globals->getRcut();
1107 >        }
1108 >        
1109 >        if (!globals->haveRsw()){
1110            sprintf(painCave.errMsg,
1111 <                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1111 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1112                    "\tOOPSE will use a default value of\n"
1113 <                  "\t0.05 * electrostaticCutoffRadius\n"
1054 <                  "\tfor the electrostaticSkinThickness\n");
1113 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1114            painCave.isFatal = 0;
1115            simError();
1116 <          theEst = 0.05 * theEcr;
1116 >          theRsw = 0.95 * theRcut;
1117          }
1118          else{
1119 <          theEst = globals->getEST();
1119 >          theRsw = globals->getRsw();
1120          }
1121 +        
1122 +        info[i].setDefaultRcut(theRcut, theRsw);
1123          
1063 –        info[i].setDefaultEcr(theEcr, theEst);
1124        }
1125      }
1126    }
# Line 1068 | Line 1128 | void SimSetup::finalInfoCheck(void){
1128    strcpy(checkPointMsg, "post processing checks out");
1129    MPIcheckPoint();
1130   #endif // is_mpi
1131 +
1132 +  // clean up the forcefield
1133 +  the_ff->cleanMe();
1134   }
1135    
1136   void SimSetup::initSystemCoords(void){
# Line 1282 | Line 1345 | void SimSetup::compList(void){
1345    LinkedMolStamp* headStamp = new LinkedMolStamp();
1346    LinkedMolStamp* currentStamp = NULL;
1347    comp_stamps = new MoleculeStamp * [n_components];
1348 +  bool haveCutoffGroups;
1349  
1350 +  haveCutoffGroups = false;
1351 +  
1352    // make an array of molecule stamps that match the components used.
1353    // also extract the used stamps out into a separate linked list
1354  
# Line 1317 | Line 1383 | void SimSetup::compList(void){
1383        headStamp->add(currentStamp);
1384        comp_stamps[i] = headStamp->match(id);
1385      }
1386 +
1387 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1388 +      haveCutoffGroups = true;    
1389    }
1390 +    
1391 +  for (i = 0; i < nInfo; i++)
1392 +    info[i].haveCutoffGroups = haveCutoffGroups;
1393  
1394   #ifdef IS_MPI
1395    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1365 | Line 1437 | void SimSetup::mpiMolDivide(void){
1437    int localMol, allMol;
1438    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1439    int local_rigid;
1368 –  vector<int> globalAtomIndex;
1440    vector<int> globalMolIndex;
1441  
1442    mpiSim = new mpiSimulation(info);
1443  
1444    mpiSim->divideLabor();
1445    globalAtomIndex = mpiSim->getGlobalAtomIndex();
1446 <  globalMolIndex = mpiSim->getGlobalMolIndex();
1446 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1447  
1448    // set up the local variables
1449  
# Line 1493 | Line 1564 | void SimSetup::makeSysArrays(void){
1564      info[l].atoms = the_atoms;
1565      info[l].molecules = the_molecules;
1566      info[l].nGlobalExcludes = 0;
1567 <
1567 >    
1568      the_ff->setSimInfo(info);
1569    }
1570   }

Diff Legend

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