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 1108 by tim, Wed Apr 14 15:37:41 2004 UTC vs.
Revision 1187 by chrisfen, Sat May 22 18:16:18 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  
# 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 215 | Line 220 | void SimSetup::makeMolecules(void){
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 262 | 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());
268 –
276   #ifdef IS_MPI
277  
278          molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
# Line 409 | 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 418 | Line 428 | void SimSetup::makeMolecules(void){
428  
429          myRB = new RigidBody();
430  
431 <        sprintf(rbName,"%s_RB_%s", molName, j);
431 >        sprintf(rbName,"%s_RB_%d", molName, j);
432          myRB->setType(rbName);
433          
434          for (rb1 = 0; rb1 < nMembers; rb1++) {
# Line 460 | 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);
# Line 477 | Line 487 | void SimSetup::makeMolecules(void){
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
# Line 529 | Line 585 | void SimSetup::makeMolecules(void){
585    MPIcheckPoint();
586   #endif // is_mpi
587  
532 –  // clean up the forcefield
533 –
534 –  if (!globals->haveLJrcut()){
535 –
536 –    the_ff->calcRcut();
537 –
538 –  } else {
539 –    
540 –    the_ff->setRcut( globals->getLJrcut() );
541 –  }
542 –
543 –  the_ff->cleanMe();
588   }
589  
590   void SimSetup::initFromBass(void){
# Line 827 | 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 837 | 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 873 | Line 917 | void SimSetup::gatherInfo(void){
917      if (globals->haveSampleTime()){
918        info[i].sampleTime = globals->getSampleTime();
919        info[i].statusTime = info[i].sampleTime;
876 –      info[i].thermalTime = info[i].sampleTime;
920      }
921      else{
922        info[i].sampleTime = globals->getRunTime();
923        info[i].statusTime = info[i].sampleTime;
881 –      info[i].thermalTime = info[i].sampleTime;
924      }
925  
926      if (globals->haveStatusTime()){
# Line 887 | 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 904 | 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->getUseThermInt()) {
954 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
955 >        info[i].useThermInt = globals->getUseThermInt();
956 >        info[i].thermIntLambda = globals->getThermIntLambda();
957 >        info[i].thermIntK = globals->getThermIntK();
958 >        
959 >        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
960 >        info[i].restraint = myRestraint;
961 >      }
962 >      else {
963 >        sprintf(painCave.errMsg,
964 >                "SimSetup Error:\n"
965 >                "\tKeyword useThermInt was set to 'true' but\n"
966 >                "\tthermodynamicIntegrationLambda (and/or\n"
967 >                "\tthermodynamicIntegrationK) was not specified.\n"
968 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
969 >        painCave.isFatal = 1;
970 >        simError();    
971 >      }
972 >    }
973 >    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
974 >        sprintf(painCave.errMsg,
975 >                "SimSetup Warning: If you want to use Thermodynamic\n"
976 >                "\tIntegration, set useThermInt to 'true' in your .bass file.\n"
977 >                "\tThe useThermInt keyword is 'false' by default, so your\n"
978 >                "\tlambda and/or k values are being ignored.\n");
979 >        painCave.isFatal = 0;
980 >        simError();  
981 >    }
982    }
983    
984    //setup seed for random number generator
# Line 957 | Line 1031 | void SimSetup::finalInfoCheck(void){
1031   void SimSetup::finalInfoCheck(void){
1032    int index;
1033    int usesDipoles;
1034 +  int usesCharges;
1035    int i;
1036  
1037    for (i = 0; i < nInfo; i++){
# Line 968 | Line 1043 | void SimSetup::finalInfoCheck(void){
1043        usesDipoles = (info[i].atoms[index])->hasDipole();
1044        index++;
1045      }
1046 <
1046 >    index = 0;
1047 >    usesCharges = 0;
1048 >    while ((index < info[i].n_atoms) && !usesCharges){
1049 >      usesCharges= (info[i].atoms[index])->hasCharge();
1050 >      index++;
1051 >    }
1052   #ifdef IS_MPI
1053      int myUse = usesDipoles;
1054      MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1055   #endif //is_mpi
1056  
1057 <    double theEcr, theEst;
1057 >    double theRcut, theRsw;
1058 >
1059 >    if (globals->haveRcut()) {
1060 >      theRcut = globals->getRcut();
1061  
1062 +      if (globals->haveRsw())
1063 +        theRsw = globals->getRsw();
1064 +      else
1065 +        theRsw = theRcut;
1066 +      
1067 +      info[i].setDefaultRcut(theRcut, theRsw);
1068 +
1069 +    } else {
1070 +      
1071 +      the_ff->calcRcut();
1072 +      theRcut = info[i].getRcut();
1073 +
1074 +      if (globals->haveRsw())
1075 +        theRsw = globals->getRsw();
1076 +      else
1077 +        theRsw = theRcut;
1078 +      
1079 +      info[i].setDefaultRcut(theRcut, theRsw);
1080 +    }
1081 +
1082      if (globals->getUseRF()){
1083        info[i].useReactionField = 1;
1084 <
1085 <      if (!globals->haveECR()){
1084 >      
1085 >      if (!globals->haveRcut()){
1086          sprintf(painCave.errMsg,
1087 <                "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1087 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1088                  "\tOOPSE will use a default value of 15.0 angstroms"
1089 <                "\tfor the electrostaticCutoffRadius.\n");
1089 >                "\tfor the cutoffRadius.\n");
1090          painCave.isFatal = 0;
1091          simError();
1092 <        theEcr = 15.0;
1092 >        theRcut = 15.0;
1093        }
1094        else{
1095 <        theEcr = globals->getECR();
1095 >        theRcut = globals->getRcut();
1096        }
1097  
1098 <      if (!globals->haveEST()){
1098 >      if (!globals->haveRsw()){
1099          sprintf(painCave.errMsg,
1100 <                "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1100 >                "SimSetup Warning: No value was set for switchingRadius.\n"
1101                  "\tOOPSE will use a default value of\n"
1102 <                "\t0.05 * electrostaticCutoffRadius\n"
1000 <                "\tfor the electrostaticSkinThickness\n");
1102 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
1103          painCave.isFatal = 0;
1104          simError();
1105 <        theEst = 0.05 * theEcr;
1105 >        theRsw = 0.95 * theRcut;
1106        }
1107        else{
1108 <        theEst = globals->getEST();
1108 >        theRsw = globals->getRsw();
1109        }
1110  
1111 <      info[i].setDefaultEcr(theEcr, theEst);
1111 >      info[i].setDefaultRcut(theRcut, theRsw);
1112  
1113        if (!globals->haveDielectric()){
1114          sprintf(painCave.errMsg,
# Line 1019 | Line 1121 | void SimSetup::finalInfoCheck(void){
1121        info[i].dielectric = globals->getDielectric();
1122      }
1123      else{
1124 <      if (usesDipoles){
1125 <        if (!globals->haveECR()){
1124 >      if (usesDipoles || usesCharges){
1125 >
1126 >        if (!globals->haveRcut()){
1127            sprintf(painCave.errMsg,
1128 <                  "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
1128 >                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1129                    "\tOOPSE will use a default value of 15.0 angstroms"
1130 <                  "\tfor the electrostaticCutoffRadius.\n");
1131 <          painCave.isFatal = 0;
1132 <          simError();
1133 <          theEcr = 15.0;
1134 <        }
1130 >                  "\tfor the cutoffRadius.\n");
1131 >          painCave.isFatal = 0;
1132 >          simError();
1133 >          theRcut = 15.0;
1134 >      }
1135          else{
1136 <          theEcr = globals->getECR();
1136 >          theRcut = globals->getRcut();
1137          }
1138 <        
1139 <        if (!globals->haveEST()){
1138 >        
1139 >        if (!globals->haveRsw()){
1140            sprintf(painCave.errMsg,
1141 <                  "SimSetup Warning: No value was set for electrostaticSkinThickness.\n"
1141 >                  "SimSetup Warning: No value was set for switchingRadius.\n"
1142                    "\tOOPSE will use a default value of\n"
1143 <                  "\t0.05 * electrostaticCutoffRadius\n"
1041 <                  "\tfor the electrostaticSkinThickness\n");
1143 >                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1144            painCave.isFatal = 0;
1145            simError();
1146 <          theEst = 0.05 * theEcr;
1146 >          theRsw = 0.95 * theRcut;
1147          }
1148          else{
1149 <          theEst = globals->getEST();
1149 >          theRsw = globals->getRsw();
1150          }
1151 +        
1152 +        info[i].setDefaultRcut(theRcut, theRsw);
1153          
1050 –        info[i].setDefaultEcr(theEcr, theEst);
1154        }
1155      }
1156    }
# Line 1055 | Line 1158 | void SimSetup::finalInfoCheck(void){
1158    strcpy(checkPointMsg, "post processing checks out");
1159    MPIcheckPoint();
1160   #endif // is_mpi
1161 +
1162 +  // clean up the forcefield
1163 +  the_ff->cleanMe();
1164   }
1165    
1166   void SimSetup::initSystemCoords(void){
# Line 1185 | Line 1291 | void SimSetup::makeOutNames(void){
1291          }
1292        }
1293  
1294 +      strcpy(info[k].rawPotName, inFileName);
1295 +      nameLength = strlen(info[k].rawPotName);
1296 +      endTest = &(info[k].rawPotName[nameLength - 5]);
1297 +      if (!strcmp(endTest, ".bass")){
1298 +        strcpy(endTest, ".raw");
1299 +      }
1300 +      else if (!strcmp(endTest, ".BASS")){
1301 +        strcpy(endTest, ".raw");
1302 +      }
1303 +      else{
1304 +        endTest = &(info[k].rawPotName[nameLength - 4]);
1305 +        if (!strcmp(endTest, ".bss")){
1306 +          strcpy(endTest, ".raw");
1307 +        }
1308 +        else if (!strcmp(endTest, ".mdl")){
1309 +          strcpy(endTest, ".raw");
1310 +        }
1311 +        else{
1312 +          strcat(info[k].rawPotName, ".raw");
1313 +        }
1314 +      }
1315 +
1316   #ifdef IS_MPI
1317  
1318      }
# Line 1269 | Line 1397 | void SimSetup::compList(void){
1397    LinkedMolStamp* headStamp = new LinkedMolStamp();
1398    LinkedMolStamp* currentStamp = NULL;
1399    comp_stamps = new MoleculeStamp * [n_components];
1400 +  bool haveCutoffGroups;
1401  
1402 +  haveCutoffGroups = false;
1403 +  
1404    // make an array of molecule stamps that match the components used.
1405    // also extract the used stamps out into a separate linked list
1406  
# Line 1304 | Line 1435 | void SimSetup::compList(void){
1435        headStamp->add(currentStamp);
1436        comp_stamps[i] = headStamp->match(id);
1437      }
1438 +
1439 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1440 +      haveCutoffGroups = true;    
1441    }
1442 +    
1443 +  for (i = 0; i < nInfo; i++)
1444 +    info[i].haveCutoffGroups = haveCutoffGroups;
1445  
1446   #ifdef IS_MPI
1447    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1352 | Line 1489 | void SimSetup::mpiMolDivide(void){
1489    int localMol, allMol;
1490    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1491    int local_rigid;
1355 –  vector<int> globalAtomIndex;
1492    vector<int> globalMolIndex;
1493  
1494    mpiSim = new mpiSimulation(info);
1495  
1496    mpiSim->divideLabor();
1497    globalAtomIndex = mpiSim->getGlobalAtomIndex();
1498 <  globalMolIndex = mpiSim->getGlobalMolIndex();
1498 >  //globalMolIndex = mpiSim->getGlobalMolIndex();
1499  
1500    // set up the local variables
1501  
# Line 1480 | Line 1616 | void SimSetup::makeSysArrays(void){
1616      info[l].atoms = the_atoms;
1617      info[l].molecules = the_molecules;
1618      info[l].nGlobalExcludes = 0;
1619 <
1619 >    
1620      the_ff->setSimInfo(info);
1621    }
1622   }

Diff Legend

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