--- trunk/OOPSE/libmdtools/SimSetup.cpp 2004/04/22 03:29:30 1129 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2004/05/12 16:38:45 1167 @@ -183,6 +183,10 @@ void SimSetup::makeMolecules(void){ BendStamp* currentBend; TorsionStamp* currentTorsion; RigidBodyStamp* currentRigidBody; + CutoffGroupStamp* currentCutoffGroup; + CutoffGroup* myCutoffGroup; + int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file + set cutoffAtomSet; //atoms belong to cutoffgroup defined at mdl file bond_pair* theBonds; bend_set* theBends; @@ -216,6 +220,8 @@ void SimSetup::makeMolecules(void){ molInfo.nBends = comp_stamps[stampID]->getNBends(); molInfo.nTorsions = comp_stamps[stampID]->getNTorsions(); molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies(); + + nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups(); molInfo.myAtoms = &(info[k].atoms[atomOffset]); @@ -480,7 +486,50 @@ void SimSetup::makeMolecules(void){ info[k].rigidBodies.push_back(myRB); } + + //create cutoff group for molecule + + cutoffAtomSet.clear(); + molInfo.myCutoffGroups.clear(); + + for (j = 0; j < nCutoffGroups; j++){ + + currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j); + nMembers = currentCutoffGroup->getNMembers(); + + myCutoffGroup = new CutoffGroup(); + + for (int cg = 0; cg < nMembers; cg++) { + + // molI is atom numbering inside this molecule + molI = currentCutoffGroup->getMember(cg); + + // tempI is atom numbering on local processor + tempI = molI + atomOffset; + + myCutoffGroup->addAtom(info[k].atoms[tempI]); + + cutoffAtomSet.insert(tempI); + } + + molInfo.myCutoffGroups.push_back(myCutoffGroup); + }//end for (j = 0; j < molInfo.nCutoffGroups; j++) + + //creat a cutoff group for every atom in current molecule which does not belong to cutoffgroup defined at mdl file + + for(j = 0; j < molInfo.nAtoms; j++){ + + if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){ + myCutoffGroup = new CutoffGroup(); + myCutoffGroup->addAtom(molInfo.myAtoms[j]); + molInfo.myCutoffGroups.push_back(myCutoffGroup); + } + + } + + + // After this is all set up, scan through the atoms to // see if they can be added to the integrableObjects: @@ -536,18 +585,6 @@ void SimSetup::makeMolecules(void){ MPIcheckPoint(); #endif // is_mpi - // clean up the forcefield - - if (!globals->haveLJrcut()){ - - the_ff->calcRcut(); - - } else { - - the_ff->setRcut( globals->getLJrcut() ); - } - - the_ff->cleanMe(); } void SimSetup::initFromBass(void){ @@ -987,39 +1024,61 @@ void SimSetup::finalInfoCheck(void){ MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); #endif //is_mpi - double theEcr, theEst; + double theRcut, theRsw; - if (globals->getUseRF()){ - info[i].useReactionField = 1; + if (globals->haveRcut()) { + theRcut = globals->getRcut(); - if (!globals->haveECR()){ + if (globals->haveRsw()) + theRsw = globals->getRsw(); + else + theRsw = theRcut; + + info[i].setDefaultRcut(theRcut, theRsw); + + } else { + + the_ff->calcRcut(); + theRcut = info[i].getRcut(); + + if (globals->haveRsw()) + theRsw = globals->getRsw(); + else + theRsw = theRcut; + + info[i].setDefaultRcut(theRcut, theRsw); + } + + if (globals->getUseRF()){ + info[i].useReactionField = 1; + + if (!globals->haveRcut()){ sprintf(painCave.errMsg, - "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" + "SimSetup Warning: No value was set for the cutoffRadius.\n" "\tOOPSE will use a default value of 15.0 angstroms" - "\tfor the electrostaticCutoffRadius.\n"); + "\tfor the cutoffRadius.\n"); painCave.isFatal = 0; simError(); - theEcr = 15.0; + theRcut = 15.0; } else{ - theEcr = globals->getECR(); + theRcut = globals->getRcut(); } - if (!globals->haveEST()){ + if (!globals->haveRsw()){ sprintf(painCave.errMsg, - "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" + "SimSetup Warning: No value was set for switchingRadius.\n" "\tOOPSE will use a default value of\n" - "\t0.05 * electrostaticCutoffRadius\n" - "\tfor the electrostaticSkinThickness\n"); + "\t0.95 * cutoffRadius for the switchingRadius\n"); painCave.isFatal = 0; simError(); - theEst = 0.05 * theEcr; + theRsw = 0.95 * theRcut; } else{ - theEst = globals->getEST(); + theRsw = globals->getRsw(); } - info[i].setDefaultEcr(theEcr, theEst); + info[i].setDefaultRcut(theRcut, theRsw); if (!globals->haveDielectric()){ sprintf(painCave.errMsg, @@ -1033,34 +1092,35 @@ void SimSetup::finalInfoCheck(void){ } else{ if (usesDipoles || usesCharges){ - if (!globals->haveECR()){ + + if (!globals->haveRcut()){ sprintf(painCave.errMsg, - "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" + "SimSetup Warning: No value was set for the cutoffRadius.\n" "\tOOPSE will use a default value of 15.0 angstroms" - "\tfor the electrostaticCutoffRadius.\n"); - painCave.isFatal = 0; - simError(); - theEcr = 15.0; - } + "\tfor the cutoffRadius.\n"); + painCave.isFatal = 0; + simError(); + theRcut = 15.0; + } else{ - theEcr = globals->getECR(); - } - - if (!globals->haveEST()){ + theRcut = globals->getRcut(); + } + + if (!globals->haveRsw()){ sprintf(painCave.errMsg, - "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" + "SimSetup Warning: No value was set for switchingRadius.\n" "\tOOPSE will use a default value of\n" - "\t0.05 * electrostaticCutoffRadius\n" - "\tfor the electrostaticSkinThickness\n"); + "\t0.95 * cutoffRadius for the switchingRadius\n"); painCave.isFatal = 0; simError(); - theEst = 0.05 * theEcr; + theRsw = 0.95 * theRcut; } else{ - theEst = globals->getEST(); + theRsw = globals->getRsw(); } + + info[i].setDefaultRcut(theRcut, theRsw); - info[i].setDefaultEcr(theEcr, theEst); } } } @@ -1068,6 +1128,9 @@ void SimSetup::finalInfoCheck(void){ strcpy(checkPointMsg, "post processing checks out"); MPIcheckPoint(); #endif // is_mpi + + // clean up the forcefield + the_ff->cleanMe(); } void SimSetup::initSystemCoords(void){ @@ -1282,7 +1345,10 @@ void SimSetup::compList(void){ LinkedMolStamp* headStamp = new LinkedMolStamp(); LinkedMolStamp* currentStamp = NULL; comp_stamps = new MoleculeStamp * [n_components]; + bool haveCutoffGroups; + haveCutoffGroups = false; + // make an array of molecule stamps that match the components used. // also extract the used stamps out into a separate linked list @@ -1317,7 +1383,13 @@ void SimSetup::compList(void){ headStamp->add(currentStamp); comp_stamps[i] = headStamp->match(id); } + + if(comp_stamps[i]->getNCutoffGroups() > 0) + haveCutoffGroups = true; } + + for (i = 0; i < nInfo; i++) + info[i].haveCutoffGroups = haveCutoffGroups; #ifdef IS_MPI strcpy(checkPointMsg, "Component stamps successfully extracted\n"); @@ -1492,7 +1564,7 @@ void SimSetup::makeSysArrays(void){ info[l].atoms = the_atoms; info[l].molecules = the_molecules; info[l].nGlobalExcludes = 0; - + the_ff->setSimInfo(info); } }