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 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC vs.
Revision 1157 by tim, Tue May 11 20:33:41 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  
169   void SimSetup::makeMolecules(void){
170    int i, j, k;
171 <  int exI, exJ, exK, exL, slI;
171 >  int exI, exJ, exK, exL, slI, slJ;
172    int tempI, tempJ, tempK, tempL;
173    int molI;
174    int stampID, atomOffset, rbOffset;
# 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 190 | Line 193 | void SimSetup::makeMolecules(void){
193    set<int> skipList;
194  
195    double phi, theta, psi;
196 +  char* molName;
197 +  char rbName[100];
198  
199    //init the forceField paramters
200  
# Line 206 | Line 211 | void SimSetup::makeMolecules(void){
211  
212      for (i = 0; i < info[k].n_mol; i++){
213        stampID = info[k].molecules[i].getStampID();
214 +      molName = comp_stamps[stampID]->getID();
215  
216        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
217        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
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 259 | 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());
265 –
273   #ifdef IS_MPI
274  
275 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
275 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
276  
277   #endif // is_mpi
278        }
# Line 406 | 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 414 | Line 424 | void SimSetup::makeMolecules(void){
424          // Create the Rigid Body:
425  
426          myRB = new RigidBody();
427 +
428 +        sprintf(rbName,"%s_RB_%d", molName, j);
429 +        myRB->setType(rbName);
430          
431          for (rb1 = 0; rb1 < nMembers; rb1++) {
432  
# Line 454 | 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);
478              
479            }
480          }
481 +
482 +        molInfo.myRigidBodies.push_back(myRB);
483 +        info[k].rigidBodies.push_back(myRB);
484        }
485        
486 +
487 +      //creat 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
521 +        slJ = molInfo.myAtoms[j]->getGlobalIndex();
522 + #else
523 +        slJ = j+atomOffset;
524 + #endif
525 +
526 +        // if they aren't on the skip list, then they can be integrated
527 +
528 +        if (skipList.find(slJ) == skipList.end()) {
529 +          mySD = (StuntDouble *) molInfo.myAtoms[j];
530 +          info[k].integrableObjects.push_back(mySD);
531 +          molInfo.myIntegrableObjects.push_back(mySD);
532 +        }
533 +      }
534 +
535 +      // all rigid bodies are integrated:
536 +
537 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
538 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
539 +        info[k].integrableObjects.push_back(mySD);      
540 +        molInfo.myIntegrableObjects.push_back(mySD);
541 +      }
542 +    
543 +      
544        // send the arrays off to the forceField for init.
545        
546        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
# Line 482 | Line 556 | void SimSetup::makeMolecules(void){
556        delete[] theBonds;
557        delete[] theBends;
558        delete[] theTorsions;
559 <    }
486 <
487 <    // build up the integrableObjects vector:
488 <
489 <    for (i = 0; i < info[k].n_atoms; i++) {
490 <      
491 < #ifdef IS_MPI
492 <      slI = info[k].atoms[i]->getGlobalIndex();
493 < #else
494 <      slI = i;
495 < #endif
496 <
497 <      if (skipList.find(slI) == skipList.end()) {
498 <        mySD = (StuntDouble *) info[k].atoms[i];
499 <        info[k].integrableObjects.push_back(mySD);
500 <      }
501 <    }
502 <    for (i = 0; i < info[k].rigidBodies.size(); i++) {
503 <      mySD = (StuntDouble *) info[k].rigidBodies[i];
504 <      info[k].integrableObjects.push_back(mySD);      
505 <    }
506 <    
559 >    }    
560    }
561  
562   #ifdef IS_MPI
# Line 513 | Line 566 | void SimSetup::makeMolecules(void){
566  
567    // clean up the forcefield
568  
569 <  if (!globals->haveLJrcut()){
569 >  if (!globals->haveRcut()){
570  
571      the_ff->calcRcut();
572  
573    } else {
574      
575 <    the_ff->setRcut( globals->getLJrcut() );
575 >    the_ff->setRcut( globals->getRcut() );
576    }
577  
578    the_ff->cleanMe();
# Line 809 | Line 862 | void SimSetup::gatherInfo(void){
862    }
863  
864    //check whether sample time, status time, thermal time and reset time are divisble by dt
865 <  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
865 >  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
866      sprintf(painCave.errMsg,
867              "Sample time is not divisible by dt.\n"
868              "\tThis will result in samples that are not uniformly\n"
# Line 819 | Line 872 | void SimSetup::gatherInfo(void){
872      simError();    
873    }
874  
875 <  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
875 >  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
876      sprintf(painCave.errMsg,
877              "Status time is not divisible by dt.\n"
878              "\tThis will result in status reports that are not uniformly\n"
# Line 855 | Line 908 | void SimSetup::gatherInfo(void){
908      if (globals->haveSampleTime()){
909        info[i].sampleTime = globals->getSampleTime();
910        info[i].statusTime = info[i].sampleTime;
858 –      info[i].thermalTime = info[i].sampleTime;
911      }
912      else{
913        info[i].sampleTime = globals->getRunTime();
914        info[i].statusTime = info[i].sampleTime;
863 –      info[i].thermalTime = info[i].sampleTime;
915      }
916  
917      if (globals->haveStatusTime()){
# Line 869 | Line 920 | void SimSetup::gatherInfo(void){
920  
921      if (globals->haveThermalTime()){
922        info[i].thermalTime = globals->getThermalTime();
923 +    } else {
924 +      info[i].thermalTime = globals->getRunTime();
925      }
926  
927      info[i].resetIntegrator = 0;
# Line 939 | Line 992 | void SimSetup::finalInfoCheck(void){
992   void SimSetup::finalInfoCheck(void){
993    int index;
994    int usesDipoles;
995 +  int usesCharges;
996    int i;
997  
998    for (i = 0; i < nInfo; i++){
# Line 950 | Line 1004 | void SimSetup::finalInfoCheck(void){
1004        usesDipoles = (info[i].atoms[index])->hasDipole();
1005        index++;
1006      }
1007 <
1007 >    index = 0;
1008 >    usesCharges = 0;
1009 >    while ((index < info[i].n_atoms) && !usesCharges){
1010 >      usesCharges= (info[i].atoms[index])->hasCharge();
1011 >      index++;
1012 >    }
1013   #ifdef IS_MPI
1014      int myUse = usesDipoles;
1015      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1016   #endif //is_mpi
1017  
1018 <    double theEcr, theEst;
1018 >    double theRcut, theRsw;
1019  
1020      if (globals->getUseRF()){
1021        info[i].useReactionField = 1;
1022  
1023 <      if (!globals->haveECR()){
1023 >      if (!globals->haveRcut()){
1024          sprintf(painCave.errMsg,
1025 <                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1025 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1026                  "\tOOPSE will use a default value of 15.0 angstroms"
1027 <                "\tfor the electrostaticCutoffRadius.\n");
1027 >                "\tfor the cutoffRadius.\n");
1028          painCave.isFatal = 0;
1029          simError();
1030 <        theEcr = 15.0;
1030 >        theRcut = 15.0;
1031        }
1032        else{
1033 <        theEcr = globals->getECR();
1033 >        theRcut = globals->getRcut();
1034        }
1035  
1036 <      if (!globals->haveEST()){
1036 >      if (!globals->haveRsw()){
1037          sprintf(painCave.errMsg,
1038 <                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1038 >                "SimSetup Warning: No value was set for switchingRadius.\n"
1039                  "\tOOPSE will use a default value of\n"
1040 <                "\t0.05 * electrostaticCutoffRadius\n"
982 <                "\tfor the electrostaticSkinThickness\n");
1040 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
1041          painCave.isFatal = 0;
1042          simError();
1043 <        theEst = 0.05 * theEcr;
1043 >        theRsw = 0.95 * theRcut;
1044        }
1045        else{
1046 <        theEst = globals->getEST();
1046 >        theRsw = globals->getRsw();
1047        }
1048  
1049 <      info[i].setDefaultEcr(theEcr, theEst);
1049 >      info[i].setDefaultRcut(theRcut, theRsw);
1050  
1051        if (!globals->haveDielectric()){
1052          sprintf(painCave.errMsg,
# Line 1001 | Line 1059 | void SimSetup::finalInfoCheck(void){
1059        info[i].dielectric = globals->getDielectric();
1060      }
1061      else{
1062 <      if (usesDipoles){
1063 <        if (!globals->haveECR()){
1062 >      if (usesDipoles || usesCharges){
1063 >
1064 >        if (!globals->haveRcut()){
1065            sprintf(painCave.errMsg,
1066 <                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1066 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1067                    "\tOOPSE will use a default value of 15.0 angstroms"
1068 <                  "\tfor the electrostaticCutoffRadius.\n");
1069 <          painCave.isFatal = 0;
1070 <          simError();
1071 <          theEcr = 15.0;
1072 <        }
1073 <        else{
1074 <          theEcr = globals->getECR();
1075 <        }
1076 <        
1077 <        if (!globals->haveEST()){
1068 >                  "\tfor the cutoffRadius.\n");
1069 >          painCave.isFatal = 0;
1070 >          simError();
1071 >          theRcut = 15.0;
1072 >      }
1073 >        else{
1074 >          theRcut = globals->getRcut();
1075 >        }
1076 >        
1077 >        if (!globals->haveRsw()){
1078            sprintf(painCave.errMsg,
1079 <                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1079 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1080                    "\tOOPSE will use a default value of\n"
1081 <                  "\t0.05 * electrostaticCutoffRadius\n"
1023 <                  "\tfor the electrostaticSkinThickness\n");
1081 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1082            painCave.isFatal = 0;
1083            simError();
1084 <          theEst = 0.05 * theEcr;
1084 >          theRsw = 0.95 * theRcut;
1085          }
1086          else{
1087 <          theEst = globals->getEST();
1087 >          theRsw = globals->getRsw();
1088          }
1089 +        
1090 +        info[i].setDefaultRcut(theRcut, theRsw);
1091          
1032 –        info[i].setDefaultEcr(theEcr, theEst);
1092        }
1093      }
1094    }
# Line 1251 | Line 1310 | void SimSetup::compList(void){
1310    LinkedMolStamp* headStamp = new LinkedMolStamp();
1311    LinkedMolStamp* currentStamp = NULL;
1312    comp_stamps = new MoleculeStamp * [n_components];
1313 +  bool haveCutoffGroups;
1314  
1315 +  haveCutoffGroups = false;
1316 +  
1317    // make an array of molecule stamps that match the components used.
1318    // also extract the used stamps out into a separate linked list
1319  
# Line 1286 | Line 1348 | void SimSetup::compList(void){
1348        headStamp->add(currentStamp);
1349        comp_stamps[i] = headStamp->match(id);
1350      }
1351 +
1352 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1353 +      haveCutoffGroups = true;    
1354    }
1355 +    
1356 +  for (i = 0; i < nInfo; i++)
1357 +    info[i].haveCutoffGroups = haveCutoffGroups;
1358  
1359   #ifdef IS_MPI
1360    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1334 | Line 1402 | void SimSetup::mpiMolDivide(void){
1402    int localMol, allMol;
1403    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1404    int local_rigid;
1405 +  vector<int> globalMolIndex;
1406  
1407    mpiSim = new mpiSimulation(info);
1408  
1409 <  globalIndex = mpiSim->divideLabor();
1409 >  mpiSim->divideLabor();
1410 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1411 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1412  
1413    // set up the local variables
1414  
# Line 1351 | Line 1422 | void SimSetup::mpiMolDivide(void){
1422    local_bends = 0;
1423    local_torsions = 0;
1424    local_rigid = 0;
1425 <  globalAtomIndex = 0;
1425 >  globalAtomCounter = 0;
1426  
1427    for (i = 0; i < n_components; i++){
1428      for (j = 0; j < components_nmol[i]; j++){
# Line 1364 | Line 1435 | void SimSetup::mpiMolDivide(void){
1435          localMol++;
1436        }      
1437        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1438 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1439 <        globalAtomIndex++;
1438 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1439 >        globalAtomCounter++;
1440        }
1441  
1442        allMol++;
# Line 1433 | Line 1504 | void SimSetup::makeSysArrays(void){
1504   #else // is_mpi
1505  
1506      molIndex = 0;
1507 <    globalAtomIndex = 0;
1507 >    globalAtomCounter = 0;
1508      for (i = 0; i < n_components; i++){
1509        for (j = 0; j < components_nmol[i]; j++){
1510          the_molecules[molIndex].setStampID(i);
1511          the_molecules[molIndex].setMyIndex(molIndex);
1512          the_molecules[molIndex].setGlobalIndex(molIndex);
1513          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1514 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1515 <          globalAtomIndex++;
1514 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1515 >          globalAtomCounter++;
1516          }
1517          molIndex++;
1518        }
# Line 1458 | Line 1529 | void SimSetup::makeSysArrays(void){
1529      info[l].atoms = the_atoms;
1530      info[l].molecules = the_molecules;
1531      info[l].nGlobalExcludes = 0;
1532 <
1532 >    
1533      the_ff->setSimInfo(info);
1534    }
1535   }

Diff Legend

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