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 1229 by gezelter, Thu Jun 3 20:02:25 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  
170   void SimSetup::makeMolecules(void){
171    int i, j, k;
172 <  int exI, exJ, exK, exL, slI;
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;
# 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 201 | Line 215 | void SimSetup::makeMolecules(void){
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();
234 +      molName = comp_stamps[stampID]->getID();
235  
236        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
237        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
238        molInfo.nBends = comp_stamps[stampID]->getNBends();
239        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
240        molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
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  
# Line 259 | Line 287 | void SimSetup::makeMolecules(void){
287          else{
288  
289            molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
290 +
291          }
292  
293          molInfo.myAtoms[j]->setType(currentAtom->getType());
265 –
294   #ifdef IS_MPI
295 <
268 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
269 <
295 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
296   #endif // is_mpi
297        }
298  
# Line 406 | Line 432 | void SimSetup::makeMolecules(void){
432          info[k].excludes->addPair(exK, exL);
433        }
434  
435 +      
436 +      molInfo.myRigidBodies.clear();
437 +      
438        for (j = 0; j < molInfo.nRigidBodies; j++){
439  
440          currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
# Line 414 | Line 443 | void SimSetup::makeMolecules(void){
443          // Create the Rigid Body:
444  
445          myRB = new RigidBody();
446 +
447 +        sprintf(rbName,"%s_RB_%d", molName, j);
448 +        myRB->setType(rbName);
449          
450          for (rb1 = 0; rb1 < nMembers; rb1++) {
451  
# Line 454 | Line 486 | void SimSetup::makeMolecules(void){
486              // used for the exclude list:
487              
488   #ifdef IS_MPI
489 <            exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
490 <            exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
489 >            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
490 >            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
491   #else
492 <            exI = tempI + 1;
493 <            exJ = tempJ + 1;
492 >            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
493 >            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
494   #endif
495              
496              info[k].excludes->addPair(exI, exJ);
497              
498            }
499          }
500 +
501 +        molInfo.myRigidBodies.push_back(myRB);
502 +        info[k].rigidBodies.push_back(myRB);
503        }
504        
505 +
506 +      //create cutoff group for molecule
507 +
508 +      cutoffAtomSet.clear();
509 +      molInfo.myCutoffGroups.clear();
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[groupOffset]);
520 + #else
521 +        myCutoffGroup->setGlobalIndex(groupOffset);
522 + #endif
523 +        
524 +        for (int cg = 0; cg < nMembers; cg++) {
525 +
526 +          // molI is atom numbering inside this molecule
527 +          molI = currentCutoffGroup->getMember(cg);    
528 +
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[groupOffset];
535 + #else
536 +          globalID = info[k].atoms[tempI]->getIndex();
537 +          info[k].globalGroupMembership[globalID] = groupOffset;
538 + #endif                    
539 +          myCutoffGroup->addAtom(info[k].atoms[tempI]);
540 +          cutoffAtomSet.insert(tempI);
541 +        }
542 +        
543 +        molInfo.myCutoffGroups.push_back(myCutoffGroup);
544 +        groupOffset++;
545 +
546 +      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
547 +      
548 +      
549 +      // create a cutoff group for every atom in current molecule which
550 +      // does not belong to cutoffgroup defined at mdl file
551 +      
552 +      for(j = 0; j < molInfo.nAtoms; j++){
553 +        
554 +        if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
555 +          myCutoffGroup = new CutoffGroup();
556 +          myCutoffGroup->addAtom(molInfo.myAtoms[j]);
557 +          
558 + #ifdef IS_MPI
559 +          myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]);
560 +          globalID = info[k].atoms[atomOffset + j]->getGlobalIndex();
561 +          info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset];
562 + #else
563 +          myCutoffGroup->setGlobalIndex(groupOffset);
564 +          globalID = info[k].atoms[atomOffset + j]->getIndex();
565 +          info[k].globalGroupMembership[globalID] = groupOffset;
566 + #endif
567 +          molInfo.myCutoffGroups.push_back(myCutoffGroup);
568 +          groupOffset++;
569 +        }          
570 +      }
571 +
572 +      // After this is all set up, scan through the atoms to
573 +      // see if they can be added to the integrableObjects:
574 +
575 +      molInfo.myIntegrableObjects.clear();
576 +      
577 +
578 +      for (j = 0; j < molInfo.nAtoms; j++){
579 +
580 + #ifdef IS_MPI
581 +        slJ = molInfo.myAtoms[j]->getGlobalIndex();
582 + #else
583 +        slJ = j+atomOffset;
584 + #endif
585 +
586 +        // if they aren't on the skip list, then they can be integrated
587 +
588 +        if (skipList.find(slJ) == skipList.end()) {
589 +          mySD = (StuntDouble *) molInfo.myAtoms[j];
590 +          info[k].integrableObjects.push_back(mySD);
591 +          molInfo.myIntegrableObjects.push_back(mySD);
592 +        }
593 +      }
594 +
595 +      // all rigid bodies are integrated:
596 +
597 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
598 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
599 +        info[k].integrableObjects.push_back(mySD);      
600 +        molInfo.myIntegrableObjects.push_back(mySD);
601 +      }
602 +
603 +
604 +    /*
605 +
606 +      //creat ConstraintPair.
607 +      molInfo.myConstraintPair.clear();
608 +      
609 +      for (j = 0; j < molInfo.nBonds; j++){
610 +
611 +        //if both atoms are in the same rigid body, just skip it
612 +        currentBond = comp_stamps[stampID]->getBond(j);
613 +        if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){
614 +
615 +          tempI = currentBond->getA() + atomOffset;
616 +          if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex))
617 +            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
618 +          else
619 +             consElement1 = new ConstraintAtom(info[k].atoms[tempI]);      
620 +
621 +          tempJ =  currentBond->getB() + atomOffset;
622 +          if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex))
623 +            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
624 +          else
625 +             consElement2 = new ConstraintAtom(info[k].atoms[tempJ]);    
626 +
627 +          consPair = new DistanceConstraintPair(consElement1, consElement2);
628 +          molInfo.myConstraintPairs.push_back(consPair);
629 +        }
630 +      }  
631 +      
632 +      //loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair
633 +      for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){
634 +        for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){
635 +          
636 +          jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2);
637 +
638 +          for(size_t m = 0; m < jointAtoms.size(); m++){          
639 +            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first);
640 +            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second);
641 +
642 +            consPair = new JointConstraintPair(consElement1, consElement2);  
643 +            molInfo.myConstraintPairs.push_back(consPair);            
644 +          }
645 +
646 +        }
647 +      }
648 +      
649 + */      
650        // send the arrays off to the forceField for init.
651        
652        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
# Line 476 | Line 656 | void SimSetup::makeMolecules(void){
656                                   theTorsions);
657  
658        info[k].molecules[i].initialize(molInfo);
659 <
660 <
659 >      
660 >      
661        atomOffset += molInfo.nAtoms;
662        delete[] theBonds;
663        delete[] theBends;
664        delete[] theTorsions;
665      }
666  
487 –    // build up the integrableObjects vector:
667  
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
668  
669 <      if (skipList.find(slI) == skipList.end()) {
670 <        mySD = (StuntDouble *) info[k].atoms[i];
671 <        info[k].integrableObjects.push_back(mySD);
672 <      }
673 <    }
674 <    for (i = 0; i < info[k].rigidBodies.size(); i++) {
675 <      mySD = (StuntDouble *) info[k].rigidBodies[i];
676 <      info[k].integrableObjects.push_back(mySD);      
677 <    }
669 > #ifdef IS_MPI    
670 >    // Since the globalGroupMembership has been zero filled and we've only
671 >    // poked values into the atoms we know, we can do an Allreduce
672 >    // to get the full globalGroupMembership array (We think).
673 >    // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
674 >    // docs said we could.
675 >
676 >    int* ggMjunk = new int[mpiSim->getNAtomsGlobal()];    
677 >
678 >    MPI_Allreduce(info[k].globalGroupMembership,
679 >                  ggMjunk,
680 >                  mpiSim->getNAtomsGlobal(),
681 >                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
682 >
683 >    for (i = 0; i < mpiSim->getNAtomsGlobal(); i++)
684 >      info[k].globalGroupMembership[i] = ggMjunk[i];
685 >
686 >    delete[] ggMjunk;
687      
688 + #endif
689 +
690 +
691 +
692    }
693  
694   #ifdef IS_MPI
695    sprintf(checkPointMsg, "all molecules initialized succesfully");
696    MPIcheckPoint();
697   #endif // is_mpi
513 –
514 –  // clean up the forcefield
515 –
516 –  if (!globals->haveLJrcut()){
698  
518 –    the_ff->calcRcut();
519 –
520 –  } else {
521 –    
522 –    the_ff->setRcut( globals->getLJrcut() );
523 –  }
524 –
525 –  the_ff->cleanMe();
699   }
700  
701   void SimSetup::initFromBass(void){
# Line 809 | Line 982 | void SimSetup::gatherInfo(void){
982    }
983  
984    //check whether sample time, status time, thermal time and reset time are divisble by dt
985 <  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
985 >  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
986      sprintf(painCave.errMsg,
987              "Sample time is not divisible by dt.\n"
988              "\tThis will result in samples that are not uniformly\n"
# Line 819 | Line 992 | void SimSetup::gatherInfo(void){
992      simError();    
993    }
994  
995 <  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
995 >  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
996      sprintf(painCave.errMsg,
997              "Status time is not divisible by dt.\n"
998              "\tThis will result in status reports that are not uniformly\n"
# Line 855 | Line 1028 | void SimSetup::gatherInfo(void){
1028      if (globals->haveSampleTime()){
1029        info[i].sampleTime = globals->getSampleTime();
1030        info[i].statusTime = info[i].sampleTime;
858 –      info[i].thermalTime = info[i].sampleTime;
1031      }
1032      else{
1033        info[i].sampleTime = globals->getRunTime();
1034        info[i].statusTime = info[i].sampleTime;
863 –      info[i].thermalTime = info[i].sampleTime;
1035      }
1036  
1037      if (globals->haveStatusTime()){
# Line 869 | Line 1040 | void SimSetup::gatherInfo(void){
1040  
1041      if (globals->haveThermalTime()){
1042        info[i].thermalTime = globals->getThermalTime();
1043 +    } else {
1044 +      info[i].thermalTime = globals->getRunTime();
1045      }
1046  
1047      info[i].resetIntegrator = 0;
# Line 886 | Line 1059 | void SimSetup::gatherInfo(void){
1059  
1060      info[i].useInitXSstate = globals->getUseInitXSstate();
1061      info[i].orthoTolerance = globals->getOrthoBoxTolerance();
1062 <    
1062 >
1063 >    // check for thermodynamic integration
1064 >    if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) {
1065 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1066 >        info[i].useSolidThermInt = globals->getUseSolidThermInt();
1067 >        info[i].thermIntLambda = globals->getThermIntLambda();
1068 >        info[i].thermIntK = globals->getThermIntK();
1069 >        
1070 >        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
1071 >        info[i].restraint = myRestraint;
1072 >      }
1073 >      else {
1074 >        sprintf(painCave.errMsg,
1075 >                "SimSetup Error:\n"
1076 >                "\tKeyword useSolidThermInt was set to 'true' but\n"
1077 >                "\tthermodynamicIntegrationLambda (and/or\n"
1078 >                "\tthermodynamicIntegrationK) was not specified.\n"
1079 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1080 >        painCave.isFatal = 1;
1081 >        simError();    
1082 >      }
1083 >    }
1084 >    else if(globals->getUseLiquidThermInt()) {
1085 >      if (globals->getUseSolidThermInt()) {
1086 >        sprintf( painCave.errMsg,
1087 >                 "SimSetup Warning: It appears that you have both solid and\n"
1088 >                 "\tliquid thermodynamic integration activated in your .bass\n"
1089 >                 "\tfile. To avoid confusion, specify only one technique in\n"
1090 >                 "\tyour .bass file. Liquid-state thermodynamic integration\n"
1091 >                 "\twill be assumed for the current simulation. If this is not\n"
1092 >                 "\twhat you desire, set useSolidThermInt to 'true' and\n"
1093 >                 "\tuseLiquidThermInt to 'false' in your .bass file.\n");
1094 >        painCave.isFatal = 0;
1095 >        simError();
1096 >      }
1097 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1098 >        info[i].useLiquidThermInt = globals->getUseLiquidThermInt();
1099 >        info[i].thermIntLambda = globals->getThermIntLambda();
1100 >        info[i].thermIntK = globals->getThermIntK();
1101 >      }
1102 >      else {
1103 >        sprintf(painCave.errMsg,
1104 >                "SimSetup Error:\n"
1105 >                "\tKeyword useLiquidThermInt was set to 'true' but\n"
1106 >                "\tthermodynamicIntegrationLambda (and/or\n"
1107 >                "\tthermodynamicIntegrationK) was not specified.\n"
1108 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1109 >        painCave.isFatal = 1;
1110 >        simError();    
1111 >      }
1112 >    }
1113 >    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
1114 >        sprintf(painCave.errMsg,
1115 >                "SimSetup Warning: If you want to use Thermodynamic\n"
1116 >                "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n"
1117 >                "\t'true' in your .bass file.  These keywords are set to\n"
1118 >                "\t'false' by default, so your lambda and/or k values are\n"
1119 >                "\tbeing ignored.\n");
1120 >        painCave.isFatal = 0;
1121 >        simError();  
1122 >    }
1123    }
1124    
1125    //setup seed for random number generator
# Line 939 | Line 1172 | void SimSetup::finalInfoCheck(void){
1172   void SimSetup::finalInfoCheck(void){
1173    int index;
1174    int usesDipoles;
1175 +  int usesCharges;
1176    int i;
1177  
1178    for (i = 0; i < nInfo; i++){
# Line 950 | Line 1184 | void SimSetup::finalInfoCheck(void){
1184        usesDipoles = (info[i].atoms[index])->hasDipole();
1185        index++;
1186      }
1187 <
1187 >    index = 0;
1188 >    usesCharges = 0;
1189 >    while ((index < info[i].n_atoms) && !usesCharges){
1190 >      usesCharges= (info[i].atoms[index])->hasCharge();
1191 >      index++;
1192 >    }
1193   #ifdef IS_MPI
1194      int myUse = usesDipoles;
1195      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1196   #endif //is_mpi
1197  
1198 <    double theEcr, theEst;
1198 >    double theRcut, theRsw;
1199  
1200 +    if (globals->haveRcut()) {
1201 +      theRcut = globals->getRcut();
1202 +
1203 +      if (globals->haveRsw())
1204 +        theRsw = globals->getRsw();
1205 +      else
1206 +        theRsw = theRcut;
1207 +      
1208 +      info[i].setDefaultRcut(theRcut, theRsw);
1209 +
1210 +    } else {
1211 +      
1212 +      the_ff->calcRcut();
1213 +      theRcut = info[i].getRcut();
1214 +
1215 +      if (globals->haveRsw())
1216 +        theRsw = globals->getRsw();
1217 +      else
1218 +        theRsw = theRcut;
1219 +      
1220 +      info[i].setDefaultRcut(theRcut, theRsw);
1221 +    }
1222 +
1223      if (globals->getUseRF()){
1224        info[i].useReactionField = 1;
1225 <
1226 <      if (!globals->haveECR()){
1225 >      
1226 >      if (!globals->haveRcut()){
1227          sprintf(painCave.errMsg,
1228 <                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1228 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1229                  "\tOOPSE will use a default value of 15.0 angstroms"
1230 <                "\tfor the electrostaticCutoffRadius.\n");
1230 >                "\tfor the cutoffRadius.\n");
1231          painCave.isFatal = 0;
1232          simError();
1233 <        theEcr = 15.0;
1233 >        theRcut = 15.0;
1234        }
1235        else{
1236 <        theEcr = globals->getECR();
1236 >        theRcut = globals->getRcut();
1237        }
1238  
1239 <      if (!globals->haveEST()){
1239 >      if (!globals->haveRsw()){
1240          sprintf(painCave.errMsg,
1241 <                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1241 >                "SimSetup Warning: No value was set for switchingRadius.\n"
1242                  "\tOOPSE will use a default value of\n"
1243 <                "\t0.05 * electrostaticCutoffRadius\n"
982 <                "\tfor the electrostaticSkinThickness\n");
1243 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
1244          painCave.isFatal = 0;
1245          simError();
1246 <        theEst = 0.05 * theEcr;
1246 >        theRsw = 0.95 * theRcut;
1247        }
1248        else{
1249 <        theEst = globals->getEST();
1249 >        theRsw = globals->getRsw();
1250        }
1251  
1252 <      info[i].setDefaultEcr(theEcr, theEst);
1252 >      info[i].setDefaultRcut(theRcut, theRsw);
1253  
1254        if (!globals->haveDielectric()){
1255          sprintf(painCave.errMsg,
# Line 1001 | Line 1262 | void SimSetup::finalInfoCheck(void){
1262        info[i].dielectric = globals->getDielectric();
1263      }
1264      else{
1265 <      if (usesDipoles){
1266 <        if (!globals->haveECR()){
1265 >      if (usesDipoles || usesCharges){
1266 >
1267 >        if (!globals->haveRcut()){
1268            sprintf(painCave.errMsg,
1269 <                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1269 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1270                    "\tOOPSE will use a default value of 15.0 angstroms"
1271 <                  "\tfor the electrostaticCutoffRadius.\n");
1272 <          painCave.isFatal = 0;
1273 <          simError();
1274 <          theEcr = 15.0;
1275 <        }
1276 <        else{
1277 <          theEcr = globals->getECR();
1271 >                  "\tfor the cutoffRadius.\n");
1272 >          painCave.isFatal = 0;
1273 >          simError();
1274 >          theRcut = 15.0;
1275 >      }
1276 >        else{
1277 >          theRcut = globals->getRcut();
1278          }
1279 <        
1280 <        if (!globals->haveEST()){
1279 >        
1280 >        if (!globals->haveRsw()){
1281            sprintf(painCave.errMsg,
1282 <                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1282 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1283                    "\tOOPSE will use a default value of\n"
1284 <                  "\t0.05 * electrostaticCutoffRadius\n"
1023 <                  "\tfor the electrostaticSkinThickness\n");
1284 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1285            painCave.isFatal = 0;
1286            simError();
1287 <          theEst = 0.05 * theEcr;
1287 >          theRsw = 0.95 * theRcut;
1288          }
1289          else{
1290 <          theEst = globals->getEST();
1290 >          theRsw = globals->getRsw();
1291          }
1292 +        
1293 +        info[i].setDefaultRcut(theRcut, theRsw);
1294          
1032 –        info[i].setDefaultEcr(theEcr, theEst);
1295        }
1296      }
1297    }
# Line 1037 | Line 1299 | void SimSetup::finalInfoCheck(void){
1299    strcpy(checkPointMsg, "post processing checks out");
1300    MPIcheckPoint();
1301   #endif // is_mpi
1302 +
1303 +  // clean up the forcefield
1304 +  the_ff->cleanMe();
1305   }
1306    
1307   void SimSetup::initSystemCoords(void){
# Line 1167 | Line 1432 | void SimSetup::makeOutNames(void){
1432          }
1433        }
1434  
1435 +      strcpy(info[k].rawPotName, inFileName);
1436 +      nameLength = strlen(info[k].rawPotName);
1437 +      endTest = &(info[k].rawPotName[nameLength - 5]);
1438 +      if (!strcmp(endTest, ".bass")){
1439 +        strcpy(endTest, ".raw");
1440 +      }
1441 +      else if (!strcmp(endTest, ".BASS")){
1442 +        strcpy(endTest, ".raw");
1443 +      }
1444 +      else{
1445 +        endTest = &(info[k].rawPotName[nameLength - 4]);
1446 +        if (!strcmp(endTest, ".bss")){
1447 +          strcpy(endTest, ".raw");
1448 +        }
1449 +        else if (!strcmp(endTest, ".mdl")){
1450 +          strcpy(endTest, ".raw");
1451 +        }
1452 +        else{
1453 +          strcat(info[k].rawPotName, ".raw");
1454 +        }
1455 +      }
1456 +
1457   #ifdef IS_MPI
1458  
1459      }
# Line 1251 | Line 1538 | void SimSetup::compList(void){
1538    LinkedMolStamp* headStamp = new LinkedMolStamp();
1539    LinkedMolStamp* currentStamp = NULL;
1540    comp_stamps = new MoleculeStamp * [n_components];
1541 +  bool haveCutoffGroups;
1542  
1543 +  haveCutoffGroups = false;
1544 +  
1545    // make an array of molecule stamps that match the components used.
1546    // also extract the used stamps out into a separate linked list
1547  
# Line 1286 | Line 1576 | void SimSetup::compList(void){
1576        headStamp->add(currentStamp);
1577        comp_stamps[i] = headStamp->match(id);
1578      }
1579 +
1580 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1581 +      haveCutoffGroups = true;    
1582    }
1583 +    
1584 +  for (i = 0; i < nInfo; i++)
1585 +    info[i].haveCutoffGroups = haveCutoffGroups;
1586  
1587   #ifdef IS_MPI
1588    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1295 | Line 1591 | void SimSetup::calcSysValues(void){
1591   }
1592  
1593   void SimSetup::calcSysValues(void){
1594 <  int i;
1594 >  int i, j;
1595 >  int ncutgroups, atomsingroups, ngroupsinstamp;
1596  
1597    int* molMembershipArray;
1598 +  CutoffGroupStamp* cg;
1599  
1600    tot_atoms = 0;
1601    tot_bonds = 0;
1602    tot_bends = 0;
1603    tot_torsions = 0;
1604    tot_rigid = 0;
1605 +  tot_groups = 0;
1606    for (i = 0; i < n_components; i++){
1607      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1608      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1609      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1610      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1611      tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1612 +
1613 +    ncutgroups = comp_stamps[i]->getNCutoffGroups();
1614 +    atomsingroups = 0;
1615 +    for (j=0; j < ncutgroups; j++) {
1616 +      cg = comp_stamps[i]->getCutoffGroup(j);
1617 +      atomsingroups += cg->getNMembers();
1618 +    }
1619 +    ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups;
1620 +    tot_groups += components_nmol[i] * ngroupsinstamp;    
1621    }
1622    
1623    tot_SRI = tot_bonds + tot_bends + tot_torsions;
# Line 1322 | Line 1630 | void SimSetup::calcSysValues(void){
1630      info[i].n_torsions = tot_torsions;
1631      info[i].n_SRI = tot_SRI;
1632      info[i].n_mol = tot_nmol;
1633 <
1633 >    info[i].ngroup = tot_groups;
1634      info[i].molMembershipArray = molMembershipArray;
1635    }
1636   }
# Line 1333 | Line 1641 | void SimSetup::mpiMolDivide(void){
1641    int i, j, k;
1642    int localMol, allMol;
1643    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1644 <  int local_rigid;
1644 >  int local_rigid, local_groups;
1645 >  vector<int> globalMolIndex;
1646 >  int ncutgroups, atomsingroups, ngroupsinstamp;
1647 >  CutoffGroupStamp* cg;
1648  
1649    mpiSim = new mpiSimulation(info);
1650  
1651 <  globalIndex = mpiSim->divideLabor();
1651 >  mpiSim->divideLabor();
1652 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1653 >  globalGroupIndex = mpiSim->getGlobalGroupIndex();
1654 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1655  
1656    // set up the local variables
1657  
# Line 1351 | Line 1665 | void SimSetup::mpiMolDivide(void){
1665    local_bends = 0;
1666    local_torsions = 0;
1667    local_rigid = 0;
1668 <  globalAtomIndex = 0;
1668 >  local_groups = 0;
1669 >  globalAtomCounter = 0;
1670  
1671    for (i = 0; i < n_components; i++){
1672      for (j = 0; j < components_nmol[i]; j++){
# Line 1361 | Line 1676 | void SimSetup::mpiMolDivide(void){
1676          local_bends += comp_stamps[i]->getNBends();
1677          local_torsions += comp_stamps[i]->getNTorsions();
1678          local_rigid += comp_stamps[i]->getNRigidBodies();
1679 +
1680 +        ncutgroups = comp_stamps[i]->getNCutoffGroups();
1681 +        atomsingroups = 0;
1682 +        for (k=0; k < ncutgroups; k++) {
1683 +          cg = comp_stamps[i]->getCutoffGroup(k);
1684 +          atomsingroups += cg->getNMembers();
1685 +        }
1686 +        ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups +
1687 +          ncutgroups;
1688 +        local_groups += ngroupsinstamp;    
1689 +
1690          localMol++;
1691        }      
1692        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1693 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1694 <        globalAtomIndex++;
1693 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1694 >        globalAtomCounter++;
1695        }
1696  
1697        allMol++;
# Line 1373 | Line 1699 | void SimSetup::mpiMolDivide(void){
1699    }
1700    local_SRI = local_bonds + local_bends + local_torsions;
1701  
1702 <  info[0].n_atoms = mpiSim->getMyNlocal();  
1702 >  info[0].n_atoms = mpiSim->getNAtomsLocal();  
1703    
1378 –
1704    if (local_atoms != info[0].n_atoms){
1705      sprintf(painCave.errMsg,
1706              "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
# Line 1385 | Line 1710 | void SimSetup::mpiMolDivide(void){
1710      simError();
1711    }
1712  
1713 +  info[0].ngroup = mpiSim->getNGroupsLocal();  
1714 +  if (local_groups != info[0].ngroup){
1715 +    sprintf(painCave.errMsg,
1716 +            "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n"
1717 +            "\tlocalGroups (%d) are not equal.\n",
1718 +            info[0].ngroup, local_groups);
1719 +    painCave.isFatal = 1;
1720 +    simError();
1721 +  }
1722 +  
1723    info[0].n_bonds = local_bonds;
1724    info[0].n_bends = local_bends;
1725    info[0].n_torsions = local_torsions;
# Line 1421 | Line 1756 | void SimSetup::makeSysArrays(void){
1756  
1757  
1758      molIndex = 0;
1759 <    for (i = 0; i < mpiSim->getTotNmol(); i++){
1759 >    for (i = 0; i < mpiSim->getNMolGlobal(); i++){
1760        if (mol2proc[i] == worldRank){
1761          the_molecules[molIndex].setStampID(molCompType[i]);
1762          the_molecules[molIndex].setMyIndex(molIndex);
# Line 1433 | Line 1768 | void SimSetup::makeSysArrays(void){
1768   #else // is_mpi
1769  
1770      molIndex = 0;
1771 <    globalAtomIndex = 0;
1771 >    globalAtomCounter = 0;
1772      for (i = 0; i < n_components; i++){
1773        for (j = 0; j < components_nmol[i]; j++){
1774          the_molecules[molIndex].setStampID(i);
1775          the_molecules[molIndex].setMyIndex(molIndex);
1776          the_molecules[molIndex].setGlobalIndex(molIndex);
1777          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1778 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1779 <          globalAtomIndex++;
1778 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1779 >          globalAtomCounter++;
1780          }
1781          molIndex++;
1782        }
# Line 1458 | Line 1793 | void SimSetup::makeSysArrays(void){
1793      info[l].atoms = the_atoms;
1794      info[l].molecules = the_molecules;
1795      info[l].nGlobalExcludes = 0;
1796 <
1796 >    
1797      the_ff->setSimInfo(info);
1798    }
1799   }

Diff Legend

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