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 1104 by gezelter, Tue Apr 13 16:26:03 2004 UTC vs.
Revision 1211 by tim, Tue Jun 1 15:57:30 2004 UTC

# Line 10 | Line 10
10   #include "Integrator.hpp"
11   #include "simError.h"
12   #include "RigidBody.hpp"
13 – //#include "ConjugateMinimizer.hpp"
13   #include "OOPSEMinimizer.hpp"
14 + //#include "ConstraintElement.hpp"
15 + //#include "ConstraintPair.hpp"
16  
17   #ifdef IS_MPI
18   #include "mpiBASS.h"
# Line 147 | Line 148 | void SimSetup::createSim(void){
148    // make the output filenames
149  
150    makeOutNames();
150 –
151 –  if (globals->haveMinimizer())
152 –    // make minimizer
153 –    makeMinimizer();
154 –  else
155 –    // make the integrator
156 –    makeIntegrator();
151    
152   #ifdef IS_MPI
153    mpiSim->mpiRefresh();
# Line 162 | Line 156 | void SimSetup::createSim(void){
156    // initialize the Fortran
157  
158    initFortran();
159 +
160 +  if (globals->haveMinimizer())
161 +    // make minimizer
162 +    makeMinimizer();
163 +  else
164 +    // make the integrator
165 +    makeIntegrator();
166 +
167   }
168  
169  
# Line 182 | Line 184 | void SimSetup::makeMolecules(void){
184    BendStamp* currentBend;
185    TorsionStamp* currentTorsion;
186    RigidBodyStamp* currentRigidBody;
187 +  CutoffGroupStamp* currentCutoffGroup;
188 +  CutoffGroup* myCutoffGroup;
189 +  int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file
190 +  set<int> cutoffAtomSet; //atoms belong to  cutoffgroup defined at mdl file
191  
192    bond_pair* theBonds;
193    bend_set* theBends;
# Line 190 | Line 196 | void SimSetup::makeMolecules(void){
196    set<int> skipList;
197  
198    double phi, theta, psi;
199 +  char* molName;
200 +  char rbName[100];
201  
202 +  //ConstraintPair* consPair; //constraint pair
203 +  //ConstraintElement* consElement1;  //first element of constraint pair
204 +  //ConstraintElement* consElement2;  //second element of constraint pair
205 +  //int whichRigidBody;
206 +  //int consAtomIndex;  //index of constraint atom in rigid body's atom array
207 +  //vector<pair<int, int> > jointAtoms;
208    //init the forceField paramters
209  
210    the_ff->readParams();
# Line 203 | Line 217 | void SimSetup::makeMolecules(void){
217      the_ff->setSimInfo(&(info[k]));
218  
219      atomOffset = 0;
220 +    groupOffset = 0;
221  
222      for (i = 0; i < info[k].n_mol; i++){
223        stampID = info[k].molecules[i].getStampID();
224 +      molName = comp_stamps[stampID]->getID();
225  
226        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
227        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
228        molInfo.nBends = comp_stamps[stampID]->getNBends();
229        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
230        molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
231 +
232 +      nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
233        
234        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
235  
236        if (molInfo.nBonds > 0)
237 <        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
237 >        molInfo.myBonds = new Bond*[molInfo.nBonds];
238        else
239          molInfo.myBonds = NULL;
240  
241        if (molInfo.nBends > 0)
242 <        molInfo.myBends = new (Bend *) [molInfo.nBends];
242 >        molInfo.myBends = new Bend*[molInfo.nBends];
243        else
244          molInfo.myBends = NULL;
245  
246        if (molInfo.nTorsions > 0)
247 <        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
247 >        molInfo.myTorsions = new Torsion *[molInfo.nTorsions];
248        else
249          molInfo.myTorsions = NULL;
250  
# Line 259 | Line 277 | void SimSetup::makeMolecules(void){
277          else{
278  
279            molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
280 +
281          }
282  
283          molInfo.myAtoms[j]->setType(currentAtom->getType());
265 –
284   #ifdef IS_MPI
285  
286 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
286 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
287  
288   #endif // is_mpi
289        }
# Line 406 | Line 424 | void SimSetup::makeMolecules(void){
424          info[k].excludes->addPair(exK, exL);
425        }
426  
427 +      
428 +      molInfo.myRigidBodies.clear();
429 +      
430        for (j = 0; j < molInfo.nRigidBodies; j++){
431  
432          currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
# Line 414 | Line 435 | void SimSetup::makeMolecules(void){
435          // Create the Rigid Body:
436  
437          myRB = new RigidBody();
438 +
439 +        sprintf(rbName,"%s_RB_%d", molName, j);
440 +        myRB->setType(rbName);
441          
442          for (rb1 = 0; rb1 < nMembers; rb1++) {
443  
# Line 454 | Line 478 | void SimSetup::makeMolecules(void){
478              // used for the exclude list:
479              
480   #ifdef IS_MPI
481 <            exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
482 <            exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
481 >            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
482 >            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
483   #else
484 <            exI = tempI + 1;
485 <            exJ = tempJ + 1;
484 >            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
485 >            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
486   #endif
487              
488              info[k].excludes->addPair(exI, exJ);
# Line 471 | Line 495 | void SimSetup::makeMolecules(void){
495        }
496        
497  
498 +      //create cutoff group for molecule
499 +
500 +      cutoffAtomSet.clear();
501 +      molInfo.myCutoffGroups.clear();
502 +      
503 +      for (j = 0; j < nCutoffGroups; j++){
504 +
505 +        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
506 +        nMembers = currentCutoffGroup->getNMembers();
507 +
508 +        myCutoffGroup = new CutoffGroup();
509 +        myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
510 +        
511 +        for (int cg = 0; cg < nMembers; cg++) {
512 +
513 +          // molI is atom numbering inside this molecule
514 +          molI = currentCutoffGroup->getMember(cg);    
515 +
516 +          // tempI is atom numbering on local processor
517 +          tempI = molI + atomOffset;
518 +
519 + #ifdef IS_MPI
520 +          globalID = info[k].atoms[tempI]->getGlobalIndex()
521 + #else
522 +          globalID = info[k].atoms[tempI]->getIndex();
523 + #endif
524 +
525 +          globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset];
526 +
527 +          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
528 +
529 +          cutoffAtomSet.insert(tempI);
530 +        }
531 +      
532 +        molInfo.myCutoffGroups.push_back(myCutoffGroup);
533 +        groupOffset++;
534 +
535 +      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
536 +
537 +      //creat a cutoff group for every atom  in current molecule which does not belong to cutoffgroup defined at mdl file
538 +
539 +      for(j = 0; j < molInfo.nAtoms; j++){
540 +
541 +        if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
542 +          myCutoffGroup = new CutoffGroup();
543 +          myCutoffGroup->addAtom(molInfo.myAtoms[j]);
544 +          myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
545 + #ifdef IS_MPI
546 +          globalID = info[k].atoms[atomOffset + j]->getGlobalIndex()
547 + #else
548 +          globalID = info[k].atoms[atomOffset + j]->getIndex();
549 + #endif
550 +          globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset];
551 +          molInfo.myCutoffGroups.push_back(myCutoffGroup);
552 +          groupOffset++;
553 +        }
554 +          
555 +      }
556 +
557        // After this is all set up, scan through the atoms to
558        // see if they can be added to the integrableObjects:
559  
560 +      molInfo.myIntegrableObjects.clear();
561 +      
562 +
563        for (j = 0; j < molInfo.nAtoms; j++){
564  
565   #ifdef IS_MPI
# Line 498 | Line 584 | void SimSetup::makeMolecules(void){
584          info[k].integrableObjects.push_back(mySD);      
585          molInfo.myIntegrableObjects.push_back(mySD);
586        }
587 <    
587 >
588 >
589 >    /*
590 >
591 >      //creat ConstraintPair.
592 >      molInfo.myConstraintPair.clear();
593        
594 +      for (j = 0; j < molInfo.nBonds; j++){
595 +
596 +        //if both atoms are in the same rigid body, just skip it
597 +        currentBond = comp_stamps[stampID]->getBond(j);
598 +        if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){
599 +
600 +          tempI = currentBond->getA() + atomOffset;
601 +          if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex))
602 +            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
603 +          else
604 +             consElement1 = new ConstraintAtom(info[k].atoms[tempI]);      
605 +
606 +          tempJ =  currentBond->getB() + atomOffset;
607 +          if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex))
608 +            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
609 +          else
610 +             consElement2 = new ConstraintAtom(info[k].atoms[tempJ]);    
611 +
612 +          consPair = new DistanceConstraintPair(consElement1, consElement2);
613 +          molInfo.myConstraintPairs.push_back(consPair);
614 +        }
615 +      }  
616 +      
617 +      //loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair
618 +      for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){
619 +        for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){
620 +          
621 +          jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2);
622 +
623 +          for(size_t m = 0; m < jointAtoms.size(); m++){          
624 +            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first);
625 +            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second);
626 +
627 +            consPair = new JointConstraintPair(consElement1, consElement2);  
628 +            molInfo.myConstraintPairs.push_back(consPair);            
629 +          }
630 +
631 +        }
632 +      }
633 +      
634 + */      
635        // send the arrays off to the forceField for init.
636        
637        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
# Line 523 | Line 655 | void SimSetup::makeMolecules(void){
655    MPIcheckPoint();
656   #endif // is_mpi
657  
526 –  // clean up the forcefield
527 –
528 –  if (!globals->haveLJrcut()){
529 –
530 –    the_ff->calcRcut();
531 –
532 –  } else {
533 –    
534 –    the_ff->setRcut( globals->getLJrcut() );
535 –  }
536 –
537 –  the_ff->cleanMe();
658   }
659  
660   void SimSetup::initFromBass(void){
# Line 821 | Line 941 | void SimSetup::gatherInfo(void){
941    }
942  
943    //check whether sample time, status time, thermal time and reset time are divisble by dt
944 <  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
944 >  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
945      sprintf(painCave.errMsg,
946              "Sample time is not divisible by dt.\n"
947              "\tThis will result in samples that are not uniformly\n"
# Line 831 | Line 951 | void SimSetup::gatherInfo(void){
951      simError();    
952    }
953  
954 <  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
954 >  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
955      sprintf(painCave.errMsg,
956              "Status time is not divisible by dt.\n"
957              "\tThis will result in status reports that are not uniformly\n"
# Line 867 | Line 987 | void SimSetup::gatherInfo(void){
987      if (globals->haveSampleTime()){
988        info[i].sampleTime = globals->getSampleTime();
989        info[i].statusTime = info[i].sampleTime;
870 –      info[i].thermalTime = info[i].sampleTime;
990      }
991      else{
992        info[i].sampleTime = globals->getRunTime();
993        info[i].statusTime = info[i].sampleTime;
875 –      info[i].thermalTime = info[i].sampleTime;
994      }
995  
996      if (globals->haveStatusTime()){
# Line 881 | Line 999 | void SimSetup::gatherInfo(void){
999  
1000      if (globals->haveThermalTime()){
1001        info[i].thermalTime = globals->getThermalTime();
1002 +    } else {
1003 +      info[i].thermalTime = globals->getRunTime();
1004      }
1005  
1006      info[i].resetIntegrator = 0;
# Line 898 | Line 1018 | void SimSetup::gatherInfo(void){
1018  
1019      info[i].useInitXSstate = globals->getUseInitXSstate();
1020      info[i].orthoTolerance = globals->getOrthoBoxTolerance();
1021 <    
1021 >
1022 >    // check for thermodynamic integration
1023 >    if (globals->getUseThermInt()) {
1024 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1025 >        info[i].useThermInt = globals->getUseThermInt();
1026 >        info[i].thermIntLambda = globals->getThermIntLambda();
1027 >        info[i].thermIntK = globals->getThermIntK();
1028 >        
1029 >        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
1030 >        info[i].restraint = myRestraint;
1031 >      }
1032 >      else {
1033 >        sprintf(painCave.errMsg,
1034 >                "SimSetup Error:\n"
1035 >                "\tKeyword useThermInt was set to 'true' but\n"
1036 >                "\tthermodynamicIntegrationLambda (and/or\n"
1037 >                "\tthermodynamicIntegrationK) was not specified.\n"
1038 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1039 >        painCave.isFatal = 1;
1040 >        simError();    
1041 >      }
1042 >    }
1043 >    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
1044 >        sprintf(painCave.errMsg,
1045 >                "SimSetup Warning: If you want to use Thermodynamic\n"
1046 >                "\tIntegration, set useThermInt to 'true' in your .bass file.\n"
1047 >                "\tThe useThermInt keyword is 'false' by default, so your\n"
1048 >                "\tlambda and/or k values are being ignored.\n");
1049 >        painCave.isFatal = 0;
1050 >        simError();  
1051 >    }
1052    }
1053    
1054    //setup seed for random number generator
# Line 951 | Line 1101 | void SimSetup::finalInfoCheck(void){
1101   void SimSetup::finalInfoCheck(void){
1102    int index;
1103    int usesDipoles;
1104 +  int usesCharges;
1105    int i;
1106  
1107    for (i = 0; i < nInfo; i++){
# Line 962 | Line 1113 | void SimSetup::finalInfoCheck(void){
1113        usesDipoles = (info[i].atoms[index])->hasDipole();
1114        index++;
1115      }
1116 <
1116 >    index = 0;
1117 >    usesCharges = 0;
1118 >    while ((index < info[i].n_atoms) && !usesCharges){
1119 >      usesCharges= (info[i].atoms[index])->hasCharge();
1120 >      index++;
1121 >    }
1122   #ifdef IS_MPI
1123      int myUse = usesDipoles;
1124      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1125   #endif //is_mpi
1126  
1127 <    double theEcr, theEst;
1127 >    double theRcut, theRsw;
1128  
1129 +    if (globals->haveRcut()) {
1130 +      theRcut = globals->getRcut();
1131 +
1132 +      if (globals->haveRsw())
1133 +        theRsw = globals->getRsw();
1134 +      else
1135 +        theRsw = theRcut;
1136 +      
1137 +      info[i].setDefaultRcut(theRcut, theRsw);
1138 +
1139 +    } else {
1140 +      
1141 +      the_ff->calcRcut();
1142 +      theRcut = info[i].getRcut();
1143 +
1144 +      if (globals->haveRsw())
1145 +        theRsw = globals->getRsw();
1146 +      else
1147 +        theRsw = theRcut;
1148 +      
1149 +      info[i].setDefaultRcut(theRcut, theRsw);
1150 +    }
1151 +
1152      if (globals->getUseRF()){
1153        info[i].useReactionField = 1;
1154 <
1155 <      if (!globals->haveECR()){
1154 >      
1155 >      if (!globals->haveRcut()){
1156          sprintf(painCave.errMsg,
1157 <                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1157 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1158                  "\tOOPSE will use a default value of 15.0 angstroms"
1159 <                "\tfor the electrostaticCutoffRadius.\n");
1159 >                "\tfor the cutoffRadius.\n");
1160          painCave.isFatal = 0;
1161          simError();
1162 <        theEcr = 15.0;
1162 >        theRcut = 15.0;
1163        }
1164        else{
1165 <        theEcr = globals->getECR();
1165 >        theRcut = globals->getRcut();
1166        }
1167  
1168 <      if (!globals->haveEST()){
1168 >      if (!globals->haveRsw()){
1169          sprintf(painCave.errMsg,
1170 <                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1170 >                "SimSetup Warning: No value was set for switchingRadius.\n"
1171                  "\tOOPSE will use a default value of\n"
1172 <                "\t0.05 * electrostaticCutoffRadius\n"
994 <                "\tfor the electrostaticSkinThickness\n");
1172 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
1173          painCave.isFatal = 0;
1174          simError();
1175 <        theEst = 0.05 * theEcr;
1175 >        theRsw = 0.95 * theRcut;
1176        }
1177        else{
1178 <        theEst = globals->getEST();
1178 >        theRsw = globals->getRsw();
1179        }
1180  
1181 <      info[i].setDefaultEcr(theEcr, theEst);
1181 >      info[i].setDefaultRcut(theRcut, theRsw);
1182  
1183        if (!globals->haveDielectric()){
1184          sprintf(painCave.errMsg,
# Line 1013 | Line 1191 | void SimSetup::finalInfoCheck(void){
1191        info[i].dielectric = globals->getDielectric();
1192      }
1193      else{
1194 <      if (usesDipoles){
1195 <        if (!globals->haveECR()){
1194 >      if (usesDipoles || usesCharges){
1195 >
1196 >        if (!globals->haveRcut()){
1197            sprintf(painCave.errMsg,
1198 <                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1198 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1199                    "\tOOPSE will use a default value of 15.0 angstroms"
1200 <                  "\tfor the electrostaticCutoffRadius.\n");
1201 <          painCave.isFatal = 0;
1202 <          simError();
1203 <          theEcr = 15.0;
1204 <        }
1200 >                  "\tfor the cutoffRadius.\n");
1201 >          painCave.isFatal = 0;
1202 >          simError();
1203 >          theRcut = 15.0;
1204 >      }
1205          else{
1206 <          theEcr = globals->getECR();
1206 >          theRcut = globals->getRcut();
1207          }
1208 <        
1209 <        if (!globals->haveEST()){
1208 >        
1209 >        if (!globals->haveRsw()){
1210            sprintf(painCave.errMsg,
1211 <                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1211 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1212                    "\tOOPSE will use a default value of\n"
1213 <                  "\t0.05 * electrostaticCutoffRadius\n"
1035 <                  "\tfor the electrostaticSkinThickness\n");
1213 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1214            painCave.isFatal = 0;
1215            simError();
1216 <          theEst = 0.05 * theEcr;
1216 >          theRsw = 0.95 * theRcut;
1217          }
1218          else{
1219 <          theEst = globals->getEST();
1219 >          theRsw = globals->getRsw();
1220          }
1221 +        
1222 +        info[i].setDefaultRcut(theRcut, theRsw);
1223          
1044 –        info[i].setDefaultEcr(theEcr, theEst);
1224        }
1225      }
1226    }
# Line 1049 | Line 1228 | void SimSetup::finalInfoCheck(void){
1228    strcpy(checkPointMsg, "post processing checks out");
1229    MPIcheckPoint();
1230   #endif // is_mpi
1231 +
1232 +  // clean up the forcefield
1233 +  the_ff->cleanMe();
1234   }
1235    
1236   void SimSetup::initSystemCoords(void){
# Line 1179 | Line 1361 | void SimSetup::makeOutNames(void){
1361          }
1362        }
1363  
1364 +      strcpy(info[k].rawPotName, inFileName);
1365 +      nameLength = strlen(info[k].rawPotName);
1366 +      endTest = &(info[k].rawPotName[nameLength - 5]);
1367 +      if (!strcmp(endTest, ".bass")){
1368 +        strcpy(endTest, ".raw");
1369 +      }
1370 +      else if (!strcmp(endTest, ".BASS")){
1371 +        strcpy(endTest, ".raw");
1372 +      }
1373 +      else{
1374 +        endTest = &(info[k].rawPotName[nameLength - 4]);
1375 +        if (!strcmp(endTest, ".bss")){
1376 +          strcpy(endTest, ".raw");
1377 +        }
1378 +        else if (!strcmp(endTest, ".mdl")){
1379 +          strcpy(endTest, ".raw");
1380 +        }
1381 +        else{
1382 +          strcat(info[k].rawPotName, ".raw");
1383 +        }
1384 +      }
1385 +
1386   #ifdef IS_MPI
1387  
1388      }
# Line 1263 | Line 1467 | void SimSetup::compList(void){
1467    LinkedMolStamp* headStamp = new LinkedMolStamp();
1468    LinkedMolStamp* currentStamp = NULL;
1469    comp_stamps = new MoleculeStamp * [n_components];
1470 +  bool haveCutoffGroups;
1471  
1472 +  haveCutoffGroups = false;
1473 +  
1474    // make an array of molecule stamps that match the components used.
1475    // also extract the used stamps out into a separate linked list
1476  
# Line 1298 | Line 1505 | void SimSetup::compList(void){
1505        headStamp->add(currentStamp);
1506        comp_stamps[i] = headStamp->match(id);
1507      }
1508 +
1509 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1510 +      haveCutoffGroups = true;    
1511    }
1512 +    
1513 +  for (i = 0; i < nInfo; i++)
1514 +    info[i].haveCutoffGroups = haveCutoffGroups;
1515  
1516   #ifdef IS_MPI
1517    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1307 | Line 1520 | void SimSetup::calcSysValues(void){
1520   }
1521  
1522   void SimSetup::calcSysValues(void){
1523 <  int i;
1523 >  int i, j;
1524 >  int ncutgroups, atomsingroups, ngroupsinstamp;
1525  
1526    int* molMembershipArray;
1527 +  CutoffGroupStamp* cg;
1528  
1529    tot_atoms = 0;
1530    tot_bonds = 0;
1531    tot_bends = 0;
1532    tot_torsions = 0;
1533    tot_rigid = 0;
1534 +  tot_groups = 0;
1535    for (i = 0; i < n_components; i++){
1536      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1537      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1538      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1539      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1540      tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1541 +
1542 +    ncutgroups = comp_stamps[i]->getNCutoffGroups();
1543 +    atomsingroups = 0;
1544 +    for (j=0; j < ncutgroups; j++) {
1545 +      cg = comp_stamps[i]->getCutoffGroup(j);
1546 +      atomsingroups += cg->getNMembers();
1547 +    }
1548 +    ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups;
1549 +    tot_groups += components_nmol[i] * ngroupsinstamp;    
1550    }
1551    
1552    tot_SRI = tot_bonds + tot_bends + tot_torsions;
# Line 1334 | Line 1559 | void SimSetup::calcSysValues(void){
1559      info[i].n_torsions = tot_torsions;
1560      info[i].n_SRI = tot_SRI;
1561      info[i].n_mol = tot_nmol;
1562 <
1562 >    info[i].ngroup = tot_groups;
1563      info[i].molMembershipArray = molMembershipArray;
1564    }
1565   }
# Line 1345 | Line 1570 | void SimSetup::mpiMolDivide(void){
1570    int i, j, k;
1571    int localMol, allMol;
1572    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1573 <  int local_rigid;
1573 >  int local_rigid, local_groups;
1574 >  vector<int> globalMolIndex;
1575 >  int ncutgroups, atomsingroups, ngroupsinstamp;
1576 >  CutoffGroupStamp* cg;
1577  
1578    mpiSim = new mpiSimulation(info);
1579  
1580 <  globalIndex = mpiSim->divideLabor();
1580 >  mpiSim->divideLabor();
1581 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1582 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1583  
1584    // set up the local variables
1585  
# Line 1363 | Line 1593 | void SimSetup::mpiMolDivide(void){
1593    local_bends = 0;
1594    local_torsions = 0;
1595    local_rigid = 0;
1596 <  globalAtomIndex = 0;
1596 >  local_groups = 0;
1597 >  globalAtomCounter = 0;
1598  
1599    for (i = 0; i < n_components; i++){
1600      for (j = 0; j < components_nmol[i]; j++){
# Line 1373 | Line 1604 | void SimSetup::mpiMolDivide(void){
1604          local_bends += comp_stamps[i]->getNBends();
1605          local_torsions += comp_stamps[i]->getNTorsions();
1606          local_rigid += comp_stamps[i]->getNRigidBodies();
1607 +
1608 +        ncutgroups = comp_stamps[i]->getNCutoffGroups();
1609 +        atomsingroups = 0;
1610 +        for (k=0; k < ncutgroups; k++) {
1611 +          cg = comp_stamps[i]->getCutoffGroup(k);
1612 +          atomsingroups += cg->getNMembers();
1613 +        }
1614 +        ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups +
1615 +          ncutgroups;
1616 +        local_groups += ngroupsinstamp;    
1617 +
1618          localMol++;
1619        }      
1620        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1621 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1622 <        globalAtomIndex++;
1621 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1622 >        globalAtomCounter++;
1623        }
1624  
1625        allMol++;
# Line 1385 | Line 1627 | void SimSetup::mpiMolDivide(void){
1627    }
1628    local_SRI = local_bonds + local_bends + local_torsions;
1629  
1630 <  info[0].n_atoms = mpiSim->getMyNlocal();  
1630 >  info[0].n_atoms = mpiSim->getNAtomsLocal();  
1631    
1390 –
1632    if (local_atoms != info[0].n_atoms){
1633      sprintf(painCave.errMsg,
1634              "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
# Line 1397 | Line 1638 | void SimSetup::mpiMolDivide(void){
1638      simError();
1639    }
1640  
1641 +  info[0].ngroup = mpiSim->getNGroupsLocal();  
1642 +  if (local_groups != info[0].ngroup){
1643 +    sprintf(painCave.errMsg,
1644 +            "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n"
1645 +            "\tlocalGroups (%d) are not equal.\n",
1646 +            info[0].ngroup, local_groups);
1647 +    painCave.isFatal = 1;
1648 +    simError();
1649 +  }
1650 +  
1651    info[0].n_bonds = local_bonds;
1652    info[0].n_bends = local_bends;
1653    info[0].n_torsions = local_torsions;
# Line 1433 | Line 1684 | void SimSetup::makeSysArrays(void){
1684  
1685  
1686      molIndex = 0;
1687 <    for (i = 0; i < mpiSim->getTotNmol(); i++){
1687 >    for (i = 0; i < mpiSim->getNMolGlobal(); i++){
1688        if (mol2proc[i] == worldRank){
1689          the_molecules[molIndex].setStampID(molCompType[i]);
1690          the_molecules[molIndex].setMyIndex(molIndex);
# Line 1445 | Line 1696 | void SimSetup::makeSysArrays(void){
1696   #else // is_mpi
1697  
1698      molIndex = 0;
1699 <    globalAtomIndex = 0;
1699 >    globalAtomCounter = 0;
1700      for (i = 0; i < n_components; i++){
1701        for (j = 0; j < components_nmol[i]; j++){
1702          the_molecules[molIndex].setStampID(i);
1703          the_molecules[molIndex].setMyIndex(molIndex);
1704          the_molecules[molIndex].setGlobalIndex(molIndex);
1705          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1706 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1707 <          globalAtomIndex++;
1706 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1707 >          globalAtomCounter++;
1708          }
1709          molIndex++;
1710        }
# Line 1470 | Line 1721 | void SimSetup::makeSysArrays(void){
1721      info[l].atoms = the_atoms;
1722      info[l].molecules = the_molecules;
1723      info[l].nGlobalExcludes = 0;
1724 <
1724 >    
1725      the_ff->setSimInfo(info);
1726    }
1727   }

Diff Legend

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