| 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" | 
| 171 | 
  | 
  int i, j, k; | 
| 172 | 
  | 
  int exI, exJ, exK, exL, slI, slJ; | 
| 173 | 
  | 
  int tempI, tempJ, tempK, tempL; | 
| 174 | 
< | 
  int molI; | 
| 175 | 
< | 
  int stampID, atomOffset, rbOffset; | 
| 174 | 
> | 
  int molI, globalID; | 
| 175 | 
> | 
  int stampID, atomOffset, rbOffset, groupOffset; | 
| 176 | 
  | 
  molInit molInfo; | 
| 177 | 
  | 
  DirectionalAtom* dAtom; | 
| 178 | 
  | 
  RigidBody* myRB; | 
| 186 | 
  | 
  RigidBodyStamp* currentRigidBody; | 
| 187 | 
  | 
  CutoffGroupStamp* currentCutoffGroup; | 
| 188 | 
  | 
  CutoffGroup* myCutoffGroup; | 
| 189 | 
< | 
   | 
| 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; | 
| 194 | 
  | 
  torsion_set* theTorsions; | 
| 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(); | 
| 215 | 
  | 
 | 
| 216 | 
  | 
  for (k = 0; k < nInfo; k++){ | 
| 217 | 
  | 
    the_ff->setSimInfo(&(info[k])); | 
| 218 | 
+ | 
 | 
| 219 | 
+ | 
#ifdef IS_MPI | 
| 220 | 
+ | 
    info[k].globalGroupMembership = new int[mpiSim->getNAtomsGlobal()]; | 
| 221 | 
+ | 
    for (i = 0; i < mpiSim->getNAtomsGlobal(); i++)  | 
| 222 | 
+ | 
      info[k].globalGroupMembership[i] = 0; | 
| 223 | 
+ | 
#else | 
| 224 | 
+ | 
    info[k].globalGroupMembership = new int[info[k].n_atoms]; | 
| 225 | 
+ | 
    for (i = 0; i < info[k].n_atoms; i++)  | 
| 226 | 
+ | 
      info[k].globalGroupMembership[i] = 0; | 
| 227 | 
+ | 
#endif | 
| 228 | 
  | 
 | 
| 229 | 
  | 
    atomOffset = 0; | 
| 230 | 
+ | 
    groupOffset = 0; | 
| 231 | 
  | 
 | 
| 232 | 
  | 
    for (i = 0; i < info[k].n_mol; i++){ | 
| 233 | 
  | 
      stampID = info[k].molecules[i].getStampID(); | 
| 238 | 
  | 
      molInfo.nBends = comp_stamps[stampID]->getNBends(); | 
| 239 | 
  | 
      molInfo.nTorsions = comp_stamps[stampID]->getNTorsions(); | 
| 240 | 
  | 
      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies(); | 
| 241 | 
< | 
      molInfo.nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups(); | 
| 241 | 
> | 
 | 
| 242 | 
> | 
      nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups(); | 
| 243 | 
  | 
       | 
| 244 | 
  | 
      molInfo.myAtoms = &(info[k].atoms[atomOffset]); | 
| 245 | 
  | 
 | 
| 246 | 
  | 
      if (molInfo.nBonds > 0)  | 
| 247 | 
< | 
        molInfo.myBonds = new (Bond *) [molInfo.nBonds]; | 
| 247 | 
> | 
        molInfo.myBonds = new Bond*[molInfo.nBonds]; | 
| 248 | 
  | 
      else  | 
| 249 | 
  | 
        molInfo.myBonds = NULL; | 
| 250 | 
  | 
 | 
| 251 | 
  | 
      if (molInfo.nBends > 0)  | 
| 252 | 
< | 
        molInfo.myBends = new (Bend *) [molInfo.nBends]; | 
| 252 | 
> | 
        molInfo.myBends = new Bend*[molInfo.nBends]; | 
| 253 | 
  | 
      else  | 
| 254 | 
  | 
        molInfo.myBends = NULL; | 
| 255 | 
  | 
 | 
| 256 | 
  | 
      if (molInfo.nTorsions > 0)  | 
| 257 | 
< | 
        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions]; | 
| 257 | 
> | 
        molInfo.myTorsions = new Torsion *[molInfo.nTorsions]; | 
| 258 | 
  | 
      else  | 
| 259 | 
  | 
        molInfo.myTorsions = NULL; | 
| 260 | 
  | 
 | 
| 292 | 
  | 
 | 
| 293 | 
  | 
        molInfo.myAtoms[j]->setType(currentAtom->getType()); | 
| 294 | 
  | 
#ifdef IS_MPI | 
| 274 | 
– | 
 | 
| 295 | 
  | 
        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]); | 
| 276 | 
– | 
 | 
| 296 | 
  | 
#endif // is_mpi | 
| 297 | 
  | 
      }  | 
| 298 | 
  | 
 | 
| 503 | 
  | 
      } | 
| 504 | 
  | 
       | 
| 505 | 
  | 
 | 
| 506 | 
< | 
      //creat cutoff group for molecule | 
| 506 | 
> | 
      //create cutoff group for molecule | 
| 507 | 
> | 
 | 
| 508 | 
> | 
      cutoffAtomSet.clear(); | 
| 509 | 
  | 
      molInfo.myCutoffGroups.clear(); | 
| 510 | 
< | 
      for (j = 0; j < molInfo.nCutoffGroups; j++){ | 
| 510 | 
> | 
       | 
| 511 | 
> | 
      for (j = 0; j < nCutoffGroups; j++){ | 
| 512 | 
  | 
 | 
| 513 | 
  | 
        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j); | 
| 514 | 
  | 
        nMembers = currentCutoffGroup->getNMembers();  | 
| 515 | 
  | 
 | 
| 516 | 
  | 
        myCutoffGroup = new CutoffGroup(); | 
| 517 | 
  | 
         | 
| 518 | 
+ | 
#ifdef IS_MPI | 
| 519 | 
+ | 
        myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]); | 
| 520 | 
+ | 
#else | 
| 521 | 
+ | 
        myCutoffGroup->setGlobalIndex(j + groupOffset); | 
| 522 | 
+ | 
#endif | 
| 523 | 
+ | 
         | 
| 524 | 
  | 
        for (int cg = 0; cg < nMembers; cg++) { | 
| 525 | 
  | 
 | 
| 526 | 
  | 
          // molI is atom numbering inside this molecule | 
| 529 | 
  | 
          // tempI is atom numbering on local processor | 
| 530 | 
  | 
          tempI = molI + atomOffset; | 
| 531 | 
  | 
 | 
| 532 | 
+ | 
#ifdef IS_MPI | 
| 533 | 
+ | 
          globalID = info[k].atoms[tempI]->getGlobalIndex(); | 
| 534 | 
+ | 
          info[k].globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset]; | 
| 535 | 
+ | 
#else  | 
| 536 | 
+ | 
          globalID = info[k].atoms[tempI]->getIndex(); | 
| 537 | 
+ | 
          info[k].globalGroupMembership[globalID] = j + groupOffset; | 
| 538 | 
+ | 
#endif | 
| 539 | 
+ | 
           | 
| 540 | 
+ | 
 | 
| 541 | 
+ | 
           | 
| 542 | 
  | 
          myCutoffGroup->addAtom(info[k].atoms[tempI]);           | 
| 543 | 
+ | 
           | 
| 544 | 
+ | 
          cutoffAtomSet.insert(tempI); | 
| 545 | 
  | 
        } | 
| 546 | 
< | 
 | 
| 546 | 
> | 
         | 
| 547 | 
  | 
        molInfo.myCutoffGroups.push_back(myCutoffGroup); | 
| 548 | 
+ | 
        groupOffset++; | 
| 549 | 
+ | 
 | 
| 550 | 
  | 
      }//end for (j = 0; j < molInfo.nCutoffGroups; j++) | 
| 551 | 
  | 
       | 
| 552 | 
+ | 
      //creat a cutoff group for every atom  in current molecule which does not belong to cutoffgroup defined at mdl file | 
| 553 | 
+ | 
       | 
| 554 | 
+ | 
      for(j = 0; j < molInfo.nAtoms; j++){ | 
| 555 | 
+ | 
         | 
| 556 | 
+ | 
        if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){ | 
| 557 | 
+ | 
          myCutoffGroup = new CutoffGroup(); | 
| 558 | 
+ | 
          myCutoffGroup->addAtom(molInfo.myAtoms[j]); | 
| 559 | 
  | 
 | 
| 560 | 
+ | 
#ifdef IS_MPI | 
| 561 | 
+ | 
          myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]); | 
| 562 | 
+ | 
          globalID = info[k].atoms[atomOffset + j]->getGlobalIndex(); | 
| 563 | 
+ | 
          info[k].globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset];  | 
| 564 | 
+ | 
#else | 
| 565 | 
+ | 
          myCutoffGroup->setGlobalIndex(j + groupOffset); | 
| 566 | 
+ | 
          globalID = info[k].atoms[atomOffset + j]->getIndex(); | 
| 567 | 
+ | 
          info[k].globalGroupMembership[globalID] = j+groupOffset; | 
| 568 | 
+ | 
#endif | 
| 569 | 
+ | 
          molInfo.myCutoffGroups.push_back(myCutoffGroup); | 
| 570 | 
+ | 
          groupOffset++; | 
| 571 | 
+ | 
        } | 
| 572 | 
+ | 
           | 
| 573 | 
+ | 
      } | 
| 574 | 
  | 
 | 
| 575 | 
  | 
      // After this is all set up, scan through the atoms to  | 
| 576 | 
  | 
      // see if they can be added to the integrableObjects: | 
| 602 | 
  | 
        info[k].integrableObjects.push_back(mySD);       | 
| 603 | 
  | 
        molInfo.myIntegrableObjects.push_back(mySD); | 
| 604 | 
  | 
      } | 
| 605 | 
< | 
     | 
| 605 | 
> | 
 | 
| 606 | 
> | 
 | 
| 607 | 
> | 
    /* | 
| 608 | 
> | 
 | 
| 609 | 
> | 
      //creat ConstraintPair. | 
| 610 | 
> | 
      molInfo.myConstraintPair.clear(); | 
| 611 | 
  | 
       | 
| 612 | 
+ | 
      for (j = 0; j < molInfo.nBonds; j++){ | 
| 613 | 
+ | 
 | 
| 614 | 
+ | 
        //if both atoms are in the same rigid body, just skip it | 
| 615 | 
+ | 
        currentBond = comp_stamps[stampID]->getBond(j); | 
| 616 | 
+ | 
        if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){ | 
| 617 | 
+ | 
 | 
| 618 | 
+ | 
          tempI = currentBond->getA() + atomOffset; | 
| 619 | 
+ | 
          if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex)) | 
| 620 | 
+ | 
            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); | 
| 621 | 
+ | 
          else | 
| 622 | 
+ | 
             consElement1 = new ConstraintAtom(info[k].atoms[tempI]);        | 
| 623 | 
+ | 
 | 
| 624 | 
+ | 
          tempJ =  currentBond->getB() + atomOffset; | 
| 625 | 
+ | 
          if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex)) | 
| 626 | 
+ | 
            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); | 
| 627 | 
+ | 
          else | 
| 628 | 
+ | 
             consElement2 = new ConstraintAtom(info[k].atoms[tempJ]);      | 
| 629 | 
+ | 
 | 
| 630 | 
+ | 
          consPair = new DistanceConstraintPair(consElement1, consElement2); | 
| 631 | 
+ | 
          molInfo.myConstraintPairs.push_back(consPair); | 
| 632 | 
+ | 
        } | 
| 633 | 
+ | 
      }   | 
| 634 | 
+ | 
       | 
| 635 | 
+ | 
      //loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair | 
| 636 | 
+ | 
      for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){ | 
| 637 | 
+ | 
        for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){ | 
| 638 | 
+ | 
           | 
| 639 | 
+ | 
          jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2); | 
| 640 | 
+ | 
 | 
| 641 | 
+ | 
          for(size_t m = 0; m < jointAtoms.size(); m++){           | 
| 642 | 
+ | 
            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first); | 
| 643 | 
+ | 
            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second); | 
| 644 | 
+ | 
 | 
| 645 | 
+ | 
            consPair = new JointConstraintPair(consElement1, consElement2);   | 
| 646 | 
+ | 
            molInfo.myConstraintPairs.push_back(consPair);             | 
| 647 | 
+ | 
          } | 
| 648 | 
+ | 
 | 
| 649 | 
+ | 
        } | 
| 650 | 
+ | 
      } | 
| 651 | 
+ | 
       | 
| 652 | 
+ | 
*/       | 
| 653 | 
  | 
      // send the arrays off to the forceField for init. | 
| 654 | 
  | 
       | 
| 655 | 
  | 
      the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); | 
| 659 | 
  | 
                                 theTorsions); | 
| 660 | 
  | 
 | 
| 661 | 
  | 
      info[k].molecules[i].initialize(molInfo); | 
| 662 | 
< | 
 | 
| 663 | 
< | 
 | 
| 662 | 
> | 
       | 
| 663 | 
> | 
       | 
| 664 | 
  | 
      atomOffset += molInfo.nAtoms; | 
| 665 | 
  | 
      delete[] theBonds; | 
| 666 | 
  | 
      delete[] theBends; | 
| 667 | 
  | 
      delete[] theTorsions; | 
| 668 | 
< | 
    }     | 
| 560 | 
< | 
  } | 
| 668 | 
> | 
    } | 
| 669 | 
  | 
 | 
| 562 | 
– | 
#ifdef IS_MPI | 
| 563 | 
– | 
  sprintf(checkPointMsg, "all molecules initialized succesfully"); | 
| 564 | 
– | 
  MPIcheckPoint(); | 
| 565 | 
– | 
#endif // is_mpi | 
| 670 | 
  | 
 | 
| 567 | 
– | 
  // clean up the forcefield | 
| 671 | 
  | 
 | 
| 672 | 
< | 
  if (!globals->haveRcut()){ | 
| 672 | 
> | 
#ifdef IS_MPI     | 
| 673 | 
> | 
    // Since the globalGroupMembership has been zero filled and we've only | 
| 674 | 
> | 
    // poked values into the atoms we know, we can do an Allreduce | 
| 675 | 
> | 
    // to get the full globalGroupMembership array (We think). | 
| 676 | 
> | 
    // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 | 
| 677 | 
> | 
    // docs said we could. | 
| 678 | 
  | 
 | 
| 679 | 
< | 
    the_ff->calcRcut(); | 
| 679 | 
> | 
    int* ggMjunk = new int[mpiSim->getNAtomsGlobal()];     | 
| 680 | 
  | 
 | 
| 681 | 
< | 
  } else { | 
| 681 | 
> | 
    MPI_Allreduce(info[k].globalGroupMembership, | 
| 682 | 
> | 
                  ggMjunk, | 
| 683 | 
> | 
                  mpiSim->getNAtomsGlobal(), | 
| 684 | 
> | 
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 685 | 
> | 
 | 
| 686 | 
> | 
    for (i = 0; i < mpiSim->getNAtomsGlobal(); i++)  | 
| 687 | 
> | 
      info[k].globalGroupMembership[i] = ggMjunk[i]; | 
| 688 | 
> | 
 | 
| 689 | 
> | 
    delete[] ggMjunk; | 
| 690 | 
  | 
     | 
| 691 | 
< | 
    the_ff->setRcut( globals->getRcut() ); | 
| 691 | 
> | 
#endif | 
| 692 | 
> | 
 | 
| 693 | 
> | 
 | 
| 694 | 
> | 
 | 
| 695 | 
  | 
  } | 
| 696 | 
  | 
 | 
| 697 | 
< | 
  the_ff->cleanMe(); | 
| 697 | 
> | 
#ifdef IS_MPI | 
| 698 | 
> | 
  sprintf(checkPointMsg, "all molecules initialized succesfully"); | 
| 699 | 
> | 
  MPIcheckPoint(); | 
| 700 | 
> | 
#endif // is_mpi | 
| 701 | 
> | 
 | 
| 702 | 
  | 
} | 
| 703 | 
  | 
 | 
| 704 | 
  | 
void SimSetup::initFromBass(void){ | 
| 1062 | 
  | 
 | 
| 1063 | 
  | 
    info[i].useInitXSstate = globals->getUseInitXSstate(); | 
| 1064 | 
  | 
    info[i].orthoTolerance = globals->getOrthoBoxTolerance(); | 
| 1065 | 
< | 
     | 
| 1065 | 
> | 
 | 
| 1066 | 
> | 
    // check for thermodynamic integration | 
| 1067 | 
> | 
    if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) { | 
| 1068 | 
> | 
      if (globals->haveThermIntLambda() && globals->haveThermIntK()) { | 
| 1069 | 
> | 
        info[i].useSolidThermInt = globals->getUseSolidThermInt(); | 
| 1070 | 
> | 
        info[i].thermIntLambda = globals->getThermIntLambda(); | 
| 1071 | 
> | 
        info[i].thermIntK = globals->getThermIntK(); | 
| 1072 | 
> | 
         | 
| 1073 | 
> | 
        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK); | 
| 1074 | 
> | 
        info[i].restraint = myRestraint; | 
| 1075 | 
> | 
      } | 
| 1076 | 
> | 
      else { | 
| 1077 | 
> | 
        sprintf(painCave.errMsg, | 
| 1078 | 
> | 
                "SimSetup Error:\n" | 
| 1079 | 
> | 
                "\tKeyword useSolidThermInt was set to 'true' but\n" | 
| 1080 | 
> | 
                "\tthermodynamicIntegrationLambda (and/or\n" | 
| 1081 | 
> | 
                "\tthermodynamicIntegrationK) was not specified.\n" | 
| 1082 | 
> | 
                "\tPlease provide a lambda value and k value in your .bass file.\n"); | 
| 1083 | 
> | 
        painCave.isFatal = 1; | 
| 1084 | 
> | 
        simError();     | 
| 1085 | 
> | 
      } | 
| 1086 | 
> | 
    } | 
| 1087 | 
> | 
    else if(globals->getUseLiquidThermInt()) { | 
| 1088 | 
> | 
      if (globals->getUseSolidThermInt()) { | 
| 1089 | 
> | 
        sprintf( painCave.errMsg, | 
| 1090 | 
> | 
                 "SimSetup Warning: It appears that you have both solid and\n" | 
| 1091 | 
> | 
                 "\tliquid thermodynamic integration activated in your .bass\n" | 
| 1092 | 
> | 
                 "\tfile. To avoid confusion, specify only one technique in\n" | 
| 1093 | 
> | 
                 "\tyour .bass file. Liquid-state thermodynamic integration\n" | 
| 1094 | 
> | 
                 "\twill be assumed for the current simulation. If this is not\n" | 
| 1095 | 
> | 
                 "\twhat you desire, set useSolidThermInt to 'true' and\n" | 
| 1096 | 
> | 
                 "\tuseLiquidThermInt to 'false' in your .bass file.\n"); | 
| 1097 | 
> | 
        painCave.isFatal = 0; | 
| 1098 | 
> | 
        simError(); | 
| 1099 | 
> | 
      } | 
| 1100 | 
> | 
      if (globals->haveThermIntLambda() && globals->haveThermIntK()) { | 
| 1101 | 
> | 
        info[i].useLiquidThermInt = globals->getUseLiquidThermInt(); | 
| 1102 | 
> | 
        info[i].thermIntLambda = globals->getThermIntLambda(); | 
| 1103 | 
> | 
        info[i].thermIntK = globals->getThermIntK(); | 
| 1104 | 
> | 
      } | 
| 1105 | 
> | 
      else { | 
| 1106 | 
> | 
        sprintf(painCave.errMsg, | 
| 1107 | 
> | 
                "SimSetup Error:\n" | 
| 1108 | 
> | 
                "\tKeyword useLiquidThermInt was set to 'true' but\n" | 
| 1109 | 
> | 
                "\tthermodynamicIntegrationLambda (and/or\n" | 
| 1110 | 
> | 
                "\tthermodynamicIntegrationK) was not specified.\n" | 
| 1111 | 
> | 
                "\tPlease provide a lambda value and k value in your .bass file.\n"); | 
| 1112 | 
> | 
        painCave.isFatal = 1; | 
| 1113 | 
> | 
        simError();     | 
| 1114 | 
> | 
      } | 
| 1115 | 
> | 
    } | 
| 1116 | 
> | 
    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){ | 
| 1117 | 
> | 
        sprintf(painCave.errMsg, | 
| 1118 | 
> | 
                "SimSetup Warning: If you want to use Thermodynamic\n" | 
| 1119 | 
> | 
                "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n" | 
| 1120 | 
> | 
                "\t'true' in your .bass file.  These keywords are set to\n" | 
| 1121 | 
> | 
                "\t'false' by default, so your lambda and/or k values are\n" | 
| 1122 | 
> | 
                "\tbeing ignored.\n"); | 
| 1123 | 
> | 
        painCave.isFatal = 0; | 
| 1124 | 
> | 
        simError();    | 
| 1125 | 
> | 
    } | 
| 1126 | 
  | 
  } | 
| 1127 | 
  | 
   | 
| 1128 | 
  | 
  //setup seed for random number generator | 
| 1199 | 
  | 
#endif //is_mpi | 
| 1200 | 
  | 
 | 
| 1201 | 
  | 
    double theRcut, theRsw; | 
| 1202 | 
+ | 
 | 
| 1203 | 
+ | 
    if (globals->haveRcut()) { | 
| 1204 | 
+ | 
      theRcut = globals->getRcut(); | 
| 1205 | 
+ | 
 | 
| 1206 | 
+ | 
      if (globals->haveRsw())  | 
| 1207 | 
+ | 
        theRsw = globals->getRsw(); | 
| 1208 | 
+ | 
      else  | 
| 1209 | 
+ | 
        theRsw = theRcut; | 
| 1210 | 
+ | 
       | 
| 1211 | 
+ | 
      info[i].setDefaultRcut(theRcut, theRsw); | 
| 1212 | 
+ | 
 | 
| 1213 | 
+ | 
    } else { | 
| 1214 | 
+ | 
       | 
| 1215 | 
+ | 
      the_ff->calcRcut(); | 
| 1216 | 
+ | 
      theRcut = info[i].getRcut(); | 
| 1217 | 
+ | 
 | 
| 1218 | 
+ | 
      if (globals->haveRsw())  | 
| 1219 | 
+ | 
        theRsw = globals->getRsw(); | 
| 1220 | 
+ | 
      else  | 
| 1221 | 
+ | 
        theRsw = theRcut; | 
| 1222 | 
+ | 
       | 
| 1223 | 
+ | 
      info[i].setDefaultRcut(theRcut, theRsw); | 
| 1224 | 
+ | 
    } | 
| 1225 | 
  | 
 | 
| 1226 | 
  | 
    if (globals->getUseRF()){ | 
| 1227 | 
  | 
      info[i].useReactionField = 1; | 
| 1228 | 
< | 
 | 
| 1228 | 
> | 
       | 
| 1229 | 
  | 
      if (!globals->haveRcut()){ | 
| 1230 | 
  | 
        sprintf(painCave.errMsg, | 
| 1231 | 
  | 
                "SimSetup Warning: No value was set for the cutoffRadius.\n" | 
| 1302 | 
  | 
  strcpy(checkPointMsg, "post processing checks out"); | 
| 1303 | 
  | 
  MPIcheckPoint(); | 
| 1304 | 
  | 
#endif // is_mpi | 
| 1305 | 
+ | 
 | 
| 1306 | 
+ | 
  // clean up the forcefield | 
| 1307 | 
+ | 
  the_ff->cleanMe(); | 
| 1308 | 
  | 
} | 
| 1309 | 
  | 
   | 
| 1310 | 
  | 
void SimSetup::initSystemCoords(void){ | 
| 1432 | 
  | 
        } | 
| 1433 | 
  | 
        else{ | 
| 1434 | 
  | 
          strcat(info[k].statusName, ".stat"); | 
| 1435 | 
+ | 
        } | 
| 1436 | 
+ | 
      } | 
| 1437 | 
+ | 
 | 
| 1438 | 
+ | 
      strcpy(info[k].rawPotName, inFileName); | 
| 1439 | 
+ | 
      nameLength = strlen(info[k].rawPotName); | 
| 1440 | 
+ | 
      endTest = &(info[k].rawPotName[nameLength - 5]); | 
| 1441 | 
+ | 
      if (!strcmp(endTest, ".bass")){ | 
| 1442 | 
+ | 
        strcpy(endTest, ".raw"); | 
| 1443 | 
+ | 
      } | 
| 1444 | 
+ | 
      else if (!strcmp(endTest, ".BASS")){ | 
| 1445 | 
+ | 
        strcpy(endTest, ".raw"); | 
| 1446 | 
+ | 
      } | 
| 1447 | 
+ | 
      else{ | 
| 1448 | 
+ | 
        endTest = &(info[k].rawPotName[nameLength - 4]); | 
| 1449 | 
+ | 
        if (!strcmp(endTest, ".bss")){ | 
| 1450 | 
+ | 
          strcpy(endTest, ".raw"); | 
| 1451 | 
  | 
        } | 
| 1452 | 
+ | 
        else if (!strcmp(endTest, ".mdl")){ | 
| 1453 | 
+ | 
          strcpy(endTest, ".raw"); | 
| 1454 | 
+ | 
        } | 
| 1455 | 
+ | 
        else{ | 
| 1456 | 
+ | 
          strcat(info[k].rawPotName, ".raw"); | 
| 1457 | 
+ | 
        } | 
| 1458 | 
  | 
      } | 
| 1459 | 
  | 
 | 
| 1460 | 
  | 
#ifdef IS_MPI | 
| 1594 | 
  | 
} | 
| 1595 | 
  | 
 | 
| 1596 | 
  | 
void SimSetup::calcSysValues(void){ | 
| 1597 | 
< | 
  int i; | 
| 1597 | 
> | 
  int i, j; | 
| 1598 | 
> | 
  int ncutgroups, atomsingroups, ngroupsinstamp; | 
| 1599 | 
  | 
 | 
| 1600 | 
  | 
  int* molMembershipArray; | 
| 1601 | 
+ | 
  CutoffGroupStamp* cg; | 
| 1602 | 
  | 
 | 
| 1603 | 
  | 
  tot_atoms = 0; | 
| 1604 | 
  | 
  tot_bonds = 0; | 
| 1605 | 
  | 
  tot_bends = 0; | 
| 1606 | 
  | 
  tot_torsions = 0; | 
| 1607 | 
  | 
  tot_rigid = 0; | 
| 1608 | 
+ | 
  tot_groups = 0; | 
| 1609 | 
  | 
  for (i = 0; i < n_components; i++){ | 
| 1610 | 
  | 
    tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); | 
| 1611 | 
  | 
    tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); | 
| 1612 | 
  | 
    tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); | 
| 1613 | 
  | 
    tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); | 
| 1614 | 
  | 
    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies(); | 
| 1615 | 
+ | 
 | 
| 1616 | 
+ | 
    ncutgroups = comp_stamps[i]->getNCutoffGroups(); | 
| 1617 | 
+ | 
    atomsingroups = 0; | 
| 1618 | 
+ | 
    for (j=0; j < ncutgroups; j++) { | 
| 1619 | 
+ | 
      cg = comp_stamps[i]->getCutoffGroup(j); | 
| 1620 | 
+ | 
      atomsingroups += cg->getNMembers(); | 
| 1621 | 
+ | 
    } | 
| 1622 | 
+ | 
    ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups; | 
| 1623 | 
+ | 
    tot_groups += components_nmol[i] * ngroupsinstamp;     | 
| 1624 | 
  | 
  } | 
| 1625 | 
  | 
   | 
| 1626 | 
  | 
  tot_SRI = tot_bonds + tot_bends + tot_torsions; | 
| 1633 | 
  | 
    info[i].n_torsions = tot_torsions; | 
| 1634 | 
  | 
    info[i].n_SRI = tot_SRI; | 
| 1635 | 
  | 
    info[i].n_mol = tot_nmol; | 
| 1636 | 
< | 
 | 
| 1636 | 
> | 
    info[i].ngroup = tot_groups; | 
| 1637 | 
  | 
    info[i].molMembershipArray = molMembershipArray; | 
| 1638 | 
  | 
  } | 
| 1639 | 
  | 
} | 
| 1644 | 
  | 
  int i, j, k; | 
| 1645 | 
  | 
  int localMol, allMol; | 
| 1646 | 
  | 
  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; | 
| 1647 | 
< | 
  int local_rigid; | 
| 1647 | 
> | 
  int local_rigid, local_groups; | 
| 1648 | 
  | 
  vector<int> globalMolIndex; | 
| 1649 | 
+ | 
  int ncutgroups, atomsingroups, ngroupsinstamp; | 
| 1650 | 
+ | 
  CutoffGroupStamp* cg; | 
| 1651 | 
  | 
 | 
| 1652 | 
  | 
  mpiSim = new mpiSimulation(info); | 
| 1653 | 
  | 
 | 
| 1654 | 
  | 
  mpiSim->divideLabor(); | 
| 1655 | 
  | 
  globalAtomIndex = mpiSim->getGlobalAtomIndex(); | 
| 1656 | 
+ | 
  globalGroupIndex = mpiSim->getGlobalGroupIndex(); | 
| 1657 | 
  | 
  //globalMolIndex = mpiSim->getGlobalMolIndex(); | 
| 1658 | 
  | 
 | 
| 1659 | 
  | 
  // set up the local variables  | 
| 1668 | 
  | 
  local_bends = 0; | 
| 1669 | 
  | 
  local_torsions = 0; | 
| 1670 | 
  | 
  local_rigid = 0; | 
| 1671 | 
+ | 
  local_groups = 0; | 
| 1672 | 
  | 
  globalAtomCounter = 0; | 
| 1673 | 
  | 
 | 
| 1674 | 
  | 
  for (i = 0; i < n_components; i++){ | 
| 1679 | 
  | 
        local_bends += comp_stamps[i]->getNBends(); | 
| 1680 | 
  | 
        local_torsions += comp_stamps[i]->getNTorsions(); | 
| 1681 | 
  | 
        local_rigid += comp_stamps[i]->getNRigidBodies(); | 
| 1682 | 
+ | 
 | 
| 1683 | 
+ | 
        ncutgroups = comp_stamps[i]->getNCutoffGroups(); | 
| 1684 | 
+ | 
        atomsingroups = 0; | 
| 1685 | 
+ | 
        for (k=0; k < ncutgroups; k++) { | 
| 1686 | 
+ | 
          cg = comp_stamps[i]->getCutoffGroup(k); | 
| 1687 | 
+ | 
          atomsingroups += cg->getNMembers(); | 
| 1688 | 
+ | 
        } | 
| 1689 | 
+ | 
        ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups +  | 
| 1690 | 
+ | 
          ncutgroups; | 
| 1691 | 
+ | 
        local_groups += ngroupsinstamp;     | 
| 1692 | 
+ | 
 | 
| 1693 | 
  | 
        localMol++; | 
| 1694 | 
  | 
      }       | 
| 1695 | 
  | 
      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ | 
| 1702 | 
  | 
  } | 
| 1703 | 
  | 
  local_SRI = local_bonds + local_bends + local_torsions; | 
| 1704 | 
  | 
 | 
| 1705 | 
< | 
  info[0].n_atoms = mpiSim->getMyNlocal();   | 
| 1705 | 
> | 
  info[0].n_atoms = mpiSim->getNAtomsLocal();   | 
| 1706 | 
  | 
   | 
| 1449 | 
– | 
 | 
| 1707 | 
  | 
  if (local_atoms != info[0].n_atoms){ | 
| 1708 | 
  | 
    sprintf(painCave.errMsg, | 
| 1709 | 
  | 
            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" | 
| 1713 | 
  | 
    simError(); | 
| 1714 | 
  | 
  } | 
| 1715 | 
  | 
 | 
| 1716 | 
+ | 
  info[0].ngroup = mpiSim->getNGroupsLocal();    | 
| 1717 | 
+ | 
  if (local_groups != info[0].ngroup){ | 
| 1718 | 
+ | 
    sprintf(painCave.errMsg, | 
| 1719 | 
+ | 
            "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n" | 
| 1720 | 
+ | 
            "\tlocalGroups (%d) are not equal.\n", | 
| 1721 | 
+ | 
            info[0].ngroup, local_groups); | 
| 1722 | 
+ | 
    painCave.isFatal = 1; | 
| 1723 | 
+ | 
    simError(); | 
| 1724 | 
+ | 
  } | 
| 1725 | 
+ | 
   | 
| 1726 | 
  | 
  info[0].n_bonds = local_bonds; | 
| 1727 | 
  | 
  info[0].n_bends = local_bends; | 
| 1728 | 
  | 
  info[0].n_torsions = local_torsions; | 
| 1759 | 
  | 
 | 
| 1760 | 
  | 
 | 
| 1761 | 
  | 
    molIndex = 0; | 
| 1762 | 
< | 
    for (i = 0; i < mpiSim->getTotNmol(); i++){ | 
| 1762 | 
> | 
    for (i = 0; i < mpiSim->getNMolGlobal(); i++){ | 
| 1763 | 
  | 
      if (mol2proc[i] == worldRank){ | 
| 1764 | 
  | 
        the_molecules[molIndex].setStampID(molCompType[i]); | 
| 1765 | 
  | 
        the_molecules[molIndex].setMyIndex(molIndex); |