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 1154 by gezelter, Tue May 11 16:00:22 2004 UTC vs.
Revision 1180 by chrisfen, Thu May 20 20:24:07 2004 UTC

# Line 183 | 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 216 | 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 479 | Line 485 | void SimSetup::makeMolecules(void){
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  
# Line 535 | Line 584 | void SimSetup::makeMolecules(void){
584    sprintf(checkPointMsg, "all molecules initialized succesfully");
585    MPIcheckPoint();
586   #endif // is_mpi
538
539  // clean up the forcefield
540
541  if (!globals->haveRcut()){
542
543    the_ff->calcRcut();
544
545  } else {
546    
547    the_ff->setRcut( globals->getRcut() );
548  }
587  
550  the_ff->cleanMe();
588   }
589  
590   void SimSetup::initFromBass(void){
# Line 911 | 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 989 | Line 1035 | void SimSetup::finalInfoCheck(void){
1035  
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 <
1063 >      
1064        if (!globals->haveRcut()){
1065          sprintf(painCave.errMsg,
1066                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
# Line 1068 | 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 1195 | Line 1267 | void SimSetup::makeOutNames(void){
1267          }
1268          else{
1269            strcat(info[k].statusName, ".stat");
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
# Line 1282 | 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 1317 | 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 1492 | 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