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 1180 by chrisfen, Thu May 20 20:24:07 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 +  CutoffGroupStamp* currentCutoffGroup;
187 +  CutoffGroup* myCutoffGroup;
188 +  int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file
189 +  set<int> cutoffAtomSet; //atoms belong to  cutoffgroup defined at mdl file
190  
191    bond_pair* theBonds;
192    bend_set* theBends;
# Line 190 | Line 195 | void SimSetup::makeMolecules(void){
195    set<int> skipList;
196  
197    double phi, theta, psi;
198 +  char* molName;
199 +  char rbName[100];
200  
201    //init the forceField paramters
202  
# Line 206 | Line 213 | void SimSetup::makeMolecules(void){
213  
214      for (i = 0; i < info[k].n_mol; i++){
215        stampID = info[k].molecules[i].getStampID();
216 +      molName = comp_stamps[stampID]->getID();
217  
218        molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
219        molInfo.nBonds = comp_stamps[stampID]->getNBonds();
220        molInfo.nBends = comp_stamps[stampID]->getNBends();
221        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
222        molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
223 +
224 +      nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
225        
226        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
227  
228        if (molInfo.nBonds > 0)
229 <        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
229 >        molInfo.myBonds = new Bond*[molInfo.nBonds];
230        else
231          molInfo.myBonds = NULL;
232  
233        if (molInfo.nBends > 0)
234 <        molInfo.myBends = new (Bend *) [molInfo.nBends];
234 >        molInfo.myBends = new Bend*[molInfo.nBends];
235        else
236          molInfo.myBends = NULL;
237  
238        if (molInfo.nTorsions > 0)
239 <        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
239 >        molInfo.myTorsions = new Torsion *[molInfo.nTorsions];
240        else
241          molInfo.myTorsions = NULL;
242  
# Line 259 | Line 269 | void SimSetup::makeMolecules(void){
269          else{
270  
271            molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
272 +
273          }
274  
275          molInfo.myAtoms[j]->setType(currentAtom->getType());
265 –
276   #ifdef IS_MPI
277  
278 <        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
278 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
279  
280   #endif // is_mpi
281        }
# Line 406 | Line 416 | void SimSetup::makeMolecules(void){
416          info[k].excludes->addPair(exK, exL);
417        }
418  
419 +      
420 +      molInfo.myRigidBodies.clear();
421 +      
422        for (j = 0; j < molInfo.nRigidBodies; j++){
423  
424          currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
# Line 414 | Line 427 | void SimSetup::makeMolecules(void){
427          // Create the Rigid Body:
428  
429          myRB = new RigidBody();
430 +
431 +        sprintf(rbName,"%s_RB_%d", molName, j);
432 +        myRB->setType(rbName);
433          
434          for (rb1 = 0; rb1 < nMembers; rb1++) {
435  
# Line 454 | Line 470 | void SimSetup::makeMolecules(void){
470              // used for the exclude list:
471              
472   #ifdef IS_MPI
473 <            exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
474 <            exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
473 >            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
474 >            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
475   #else
476 <            exI = tempI + 1;
477 <            exJ = tempJ + 1;
476 >            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
477 >            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
478   #endif
479              
480              info[k].excludes->addPair(exI, exJ);
481              
482            }
483          }
484 +
485 +        molInfo.myRigidBodies.push_back(myRB);
486 +        info[k].rigidBodies.push_back(myRB);
487        }
488        
489 +
490 +      //create cutoff group for molecule
491 +
492 +      cutoffAtomSet.clear();
493 +      molInfo.myCutoffGroups.clear();
494 +      
495 +      for (j = 0; j < nCutoffGroups; j++){
496 +
497 +        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
498 +        nMembers = currentCutoffGroup->getNMembers();
499 +
500 +        myCutoffGroup = new CutoffGroup();
501 +        
502 +        for (int cg = 0; cg < nMembers; cg++) {
503 +
504 +          // molI is atom numbering inside this molecule
505 +          molI = currentCutoffGroup->getMember(cg);    
506 +
507 +          // tempI is atom numbering on local processor
508 +          tempI = molI + atomOffset;
509 +          
510 +          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
511 +
512 +          cutoffAtomSet.insert(tempI);
513 +        }
514 +
515 +        molInfo.myCutoffGroups.push_back(myCutoffGroup);
516 +      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
517 +
518 +      //creat a cutoff group for every atom  in current molecule which does not belong to cutoffgroup defined at mdl file
519 +
520 +      for(j = 0; j < molInfo.nAtoms; j++){
521 +
522 +        if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
523 +          myCutoffGroup = new CutoffGroup();
524 +          myCutoffGroup->addAtom(molInfo.myAtoms[j]);
525 +          molInfo.myCutoffGroups.push_back(myCutoffGroup);
526 +        }
527 +          
528 +      }
529 +
530 +              
531 +
532 +
533 +      // After this is all set up, scan through the atoms to
534 +      // see if they can be added to the integrableObjects:
535 +
536 +      molInfo.myIntegrableObjects.clear();
537 +      
538 +
539 +      for (j = 0; j < molInfo.nAtoms; j++){
540 +
541 + #ifdef IS_MPI
542 +        slJ = molInfo.myAtoms[j]->getGlobalIndex();
543 + #else
544 +        slJ = j+atomOffset;
545 + #endif
546 +
547 +        // if they aren't on the skip list, then they can be integrated
548 +
549 +        if (skipList.find(slJ) == skipList.end()) {
550 +          mySD = (StuntDouble *) molInfo.myAtoms[j];
551 +          info[k].integrableObjects.push_back(mySD);
552 +          molInfo.myIntegrableObjects.push_back(mySD);
553 +        }
554 +      }
555 +
556 +      // all rigid bodies are integrated:
557 +
558 +      for (j = 0; j < molInfo.nRigidBodies; j++) {
559 +        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
560 +        info[k].integrableObjects.push_back(mySD);      
561 +        molInfo.myIntegrableObjects.push_back(mySD);
562 +      }
563 +    
564 +      
565        // send the arrays off to the forceField for init.
566        
567        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
# Line 482 | Line 577 | void SimSetup::makeMolecules(void){
577        delete[] theBonds;
578        delete[] theBends;
579        delete[] theTorsions;
580 <    }
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 <    
580 >    }    
581    }
582  
583   #ifdef IS_MPI
584    sprintf(checkPointMsg, "all molecules initialized succesfully");
585    MPIcheckPoint();
586   #endif // is_mpi
513 –
514 –  // clean up the forcefield
515 –
516 –  if (!globals->haveLJrcut()){
517 –
518 –    the_ff->calcRcut();
587  
520 –  } else {
521 –    
522 –    the_ff->setRcut( globals->getLJrcut() );
523 –  }
524 –
525 –  the_ff->cleanMe();
588   }
589  
590   void SimSetup::initFromBass(void){
# Line 809 | Line 871 | void SimSetup::gatherInfo(void){
871    }
872  
873    //check whether sample time, status time, thermal time and reset time are divisble by dt
874 <  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
874 >  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
875      sprintf(painCave.errMsg,
876              "Sample time is not divisible by dt.\n"
877              "\tThis will result in samples that are not uniformly\n"
# Line 819 | Line 881 | void SimSetup::gatherInfo(void){
881      simError();    
882    }
883  
884 <  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
884 >  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
885      sprintf(painCave.errMsg,
886              "Status time is not divisible by dt.\n"
887              "\tThis will result in status reports that are not uniformly\n"
# Line 855 | Line 917 | void SimSetup::gatherInfo(void){
917      if (globals->haveSampleTime()){
918        info[i].sampleTime = globals->getSampleTime();
919        info[i].statusTime = info[i].sampleTime;
858 –      info[i].thermalTime = info[i].sampleTime;
920      }
921      else{
922        info[i].sampleTime = globals->getRunTime();
923        info[i].statusTime = info[i].sampleTime;
863 –      info[i].thermalTime = info[i].sampleTime;
924      }
925  
926      if (globals->haveStatusTime()){
# Line 869 | Line 929 | void SimSetup::gatherInfo(void){
929  
930      if (globals->haveThermalTime()){
931        info[i].thermalTime = globals->getThermalTime();
932 +    } else {
933 +      info[i].thermalTime = globals->getRunTime();
934      }
935  
936      info[i].resetIntegrator = 0;
# Line 886 | Line 948 | void SimSetup::gatherInfo(void){
948  
949      info[i].useInitXSstate = globals->getUseInitXSstate();
950      info[i].orthoTolerance = globals->getOrthoBoxTolerance();
951 <    
951 >
952 >    // check for thermodynamic integration
953 >    if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
954 >      info[i].thermIntLambda = globals->getThermIntLambda();
955 >      info[i].thermIntK = globals->getThermIntK();
956 >      info[i].useThermInt = 1;
957 >      
958 >      Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
959 >      info[i].restraint = myRestraint;
960 >    }
961    }
962    
963    //setup seed for random number generator
# Line 939 | Line 1010 | void SimSetup::finalInfoCheck(void){
1010   void SimSetup::finalInfoCheck(void){
1011    int index;
1012    int usesDipoles;
1013 +  int usesCharges;
1014    int i;
1015  
1016    for (i = 0; i < nInfo; i++){
# Line 950 | Line 1022 | void SimSetup::finalInfoCheck(void){
1022        usesDipoles = (info[i].atoms[index])->hasDipole();
1023        index++;
1024      }
1025 <
1025 >    index = 0;
1026 >    usesCharges = 0;
1027 >    while ((index < info[i].n_atoms) && !usesCharges){
1028 >      usesCharges= (info[i].atoms[index])->hasCharge();
1029 >      index++;
1030 >    }
1031   #ifdef IS_MPI
1032      int myUse = usesDipoles;
1033      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1034   #endif //is_mpi
1035  
1036 <    double theEcr, theEst;
1036 >    double theRcut, theRsw;
1037  
1038 +    if (globals->haveRcut()) {
1039 +      theRcut = globals->getRcut();
1040 +
1041 +      if (globals->haveRsw())
1042 +        theRsw = globals->getRsw();
1043 +      else
1044 +        theRsw = theRcut;
1045 +      
1046 +      info[i].setDefaultRcut(theRcut, theRsw);
1047 +
1048 +    } else {
1049 +      
1050 +      the_ff->calcRcut();
1051 +      theRcut = info[i].getRcut();
1052 +
1053 +      if (globals->haveRsw())
1054 +        theRsw = globals->getRsw();
1055 +      else
1056 +        theRsw = theRcut;
1057 +      
1058 +      info[i].setDefaultRcut(theRcut, theRsw);
1059 +    }
1060 +
1061      if (globals->getUseRF()){
1062        info[i].useReactionField = 1;
1063 <
1064 <      if (!globals->haveECR()){
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");
1068 >                "\tfor the cutoffRadius.\n");
1069          painCave.isFatal = 0;
1070          simError();
1071 <        theEcr = 15.0;
1071 >        theRcut = 15.0;
1072        }
1073        else{
1074 <        theEcr = globals->getECR();
1074 >        theRcut = globals->getRcut();
1075        }
1076  
1077 <      if (!globals->haveEST()){
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"
982 <                "\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].setDefaultEcr(theEcr, theEst);
1090 >      info[i].setDefaultRcut(theRcut, theRsw);
1091  
1092        if (!globals->haveDielectric()){
1093          sprintf(painCave.errMsg,
# Line 1001 | Line 1100 | void SimSetup::finalInfoCheck(void){
1100        info[i].dielectric = globals->getDielectric();
1101      }
1102      else{
1103 <      if (usesDipoles){
1104 <        if (!globals->haveECR()){
1103 >      if (usesDipoles || usesCharges){
1104 >
1105 >        if (!globals->haveRcut()){
1106            sprintf(painCave.errMsg,
1107 <                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1107 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1108                    "\tOOPSE will use a default value of 15.0 angstroms"
1109 <                  "\tfor the electrostaticCutoffRadius.\n");
1110 <          painCave.isFatal = 0;
1111 <          simError();
1112 <          theEcr = 15.0;
1113 <        }
1109 >                  "\tfor the cutoffRadius.\n");
1110 >          painCave.isFatal = 0;
1111 >          simError();
1112 >          theRcut = 15.0;
1113 >      }
1114          else{
1115 <          theEcr = globals->getECR();
1115 >          theRcut = globals->getRcut();
1116          }
1117 <        
1118 <        if (!globals->haveEST()){
1117 >        
1118 >        if (!globals->haveRsw()){
1119            sprintf(painCave.errMsg,
1120 <                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1120 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1121                    "\tOOPSE will use a default value of\n"
1122 <                  "\t0.05 * electrostaticCutoffRadius\n"
1023 <                  "\tfor the electrostaticSkinThickness\n");
1122 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1123            painCave.isFatal = 0;
1124            simError();
1125 <          theEst = 0.05 * theEcr;
1125 >          theRsw = 0.95 * theRcut;
1126          }
1127          else{
1128 <          theEst = globals->getEST();
1128 >          theRsw = globals->getRsw();
1129          }
1130 +        
1131 +        info[i].setDefaultRcut(theRcut, theRsw);
1132          
1032 –        info[i].setDefaultEcr(theEcr, theEst);
1133        }
1134      }
1135    }
# Line 1037 | Line 1137 | void SimSetup::finalInfoCheck(void){
1137    strcpy(checkPointMsg, "post processing checks out");
1138    MPIcheckPoint();
1139   #endif // is_mpi
1140 +
1141 +  // clean up the forcefield
1142 +  the_ff->cleanMe();
1143   }
1144    
1145   void SimSetup::initSystemCoords(void){
# Line 1167 | Line 1270 | void SimSetup::makeOutNames(void){
1270          }
1271        }
1272  
1273 +      strcpy(info[k].rawPotName, inFileName);
1274 +      nameLength = strlen(info[k].rawPotName);
1275 +      endTest = &(info[k].rawPotName[nameLength - 5]);
1276 +      if (!strcmp(endTest, ".bass")){
1277 +        strcpy(endTest, ".raw");
1278 +      }
1279 +      else if (!strcmp(endTest, ".BASS")){
1280 +        strcpy(endTest, ".raw");
1281 +      }
1282 +      else{
1283 +        endTest = &(info[k].rawPotName[nameLength - 4]);
1284 +        if (!strcmp(endTest, ".bss")){
1285 +          strcpy(endTest, ".raw");
1286 +        }
1287 +        else if (!strcmp(endTest, ".mdl")){
1288 +          strcpy(endTest, ".raw");
1289 +        }
1290 +        else{
1291 +          strcat(info[k].rawPotName, ".raw");
1292 +        }
1293 +      }
1294 +
1295   #ifdef IS_MPI
1296  
1297      }
# Line 1251 | Line 1376 | void SimSetup::compList(void){
1376    LinkedMolStamp* headStamp = new LinkedMolStamp();
1377    LinkedMolStamp* currentStamp = NULL;
1378    comp_stamps = new MoleculeStamp * [n_components];
1379 +  bool haveCutoffGroups;
1380  
1381 +  haveCutoffGroups = false;
1382 +  
1383    // make an array of molecule stamps that match the components used.
1384    // also extract the used stamps out into a separate linked list
1385  
# Line 1286 | Line 1414 | void SimSetup::compList(void){
1414        headStamp->add(currentStamp);
1415        comp_stamps[i] = headStamp->match(id);
1416      }
1417 +
1418 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1419 +      haveCutoffGroups = true;    
1420    }
1421 +    
1422 +  for (i = 0; i < nInfo; i++)
1423 +    info[i].haveCutoffGroups = haveCutoffGroups;
1424  
1425   #ifdef IS_MPI
1426    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1334 | Line 1468 | void SimSetup::mpiMolDivide(void){
1468    int localMol, allMol;
1469    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1470    int local_rigid;
1471 +  vector<int> globalMolIndex;
1472  
1473    mpiSim = new mpiSimulation(info);
1474  
1475 <  globalIndex = mpiSim->divideLabor();
1475 >  mpiSim->divideLabor();
1476 >  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1477 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1478  
1479    // set up the local variables
1480  
# Line 1351 | Line 1488 | void SimSetup::mpiMolDivide(void){
1488    local_bends = 0;
1489    local_torsions = 0;
1490    local_rigid = 0;
1491 <  globalAtomIndex = 0;
1491 >  globalAtomCounter = 0;
1492  
1493    for (i = 0; i < n_components; i++){
1494      for (j = 0; j < components_nmol[i]; j++){
# Line 1364 | Line 1501 | void SimSetup::mpiMolDivide(void){
1501          localMol++;
1502        }      
1503        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1504 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1505 <        globalAtomIndex++;
1504 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1505 >        globalAtomCounter++;
1506        }
1507  
1508        allMol++;
# Line 1433 | Line 1570 | void SimSetup::makeSysArrays(void){
1570   #else // is_mpi
1571  
1572      molIndex = 0;
1573 <    globalAtomIndex = 0;
1573 >    globalAtomCounter = 0;
1574      for (i = 0; i < n_components; i++){
1575        for (j = 0; j < components_nmol[i]; j++){
1576          the_molecules[molIndex].setStampID(i);
1577          the_molecules[molIndex].setMyIndex(molIndex);
1578          the_molecules[molIndex].setGlobalIndex(molIndex);
1579          for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1580 <          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1581 <          globalAtomIndex++;
1580 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1581 >          globalAtomCounter++;
1582          }
1583          molIndex++;
1584        }
# Line 1458 | Line 1595 | void SimSetup::makeSysArrays(void){
1595      info[l].atoms = the_atoms;
1596      info[l].molecules = the_molecules;
1597      info[l].nGlobalExcludes = 0;
1598 <
1598 >    
1599      the_ff->setSimInfo(info);
1600    }
1601   }

Diff Legend

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