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 1163 by gezelter, Wed May 12 14:30:12 2004 UTC vs.
Revision 1187 by chrisfen, Sat May 22 18:16:18 2004 UTC

# Line 185 | Line 185 | void SimSetup::makeMolecules(void){
185    RigidBodyStamp* currentRigidBody;
186    CutoffGroupStamp* currentCutoffGroup;
187    CutoffGroup* myCutoffGroup;
188 <  
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;
193    torsion_set* theTorsions;
# Line 218 | 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 <      molInfo.nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
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 485 | Line 488 | void SimSetup::makeMolecules(void){
488        
489  
490        //create cutoff group for molecule
491 +
492 +      cutoffAtomSet.clear();
493        molInfo.myCutoffGroups.clear();
494 <      for (j = 0; j < molInfo.nCutoffGroups; j++){
494 >      
495 >      for (j = 0; j < nCutoffGroups; j++){
496  
497          currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
498          nMembers = currentCutoffGroup->getNMembers();
# Line 500 | Line 506 | void SimSetup::makeMolecules(void){
506  
507            // tempI is atom numbering on local processor
508            tempI = molI + atomOffset;
509 <
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++)
509      
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:
# Line 927 | 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 1240 | 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      }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines