--- trunk/OOPSE/libmdtools/SimSetup.cpp 2004/04/14 15:37:41 1108 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2004/05/11 16:00:22 1154 @@ -147,13 +147,6 @@ void SimSetup::createSim(void){ // make the output filenames makeOutNames(); - - if (globals->haveMinimizer()) - // make minimizer - makeMinimizer(); - else - // make the integrator - makeIntegrator(); #ifdef IS_MPI mpiSim->mpiRefresh(); @@ -162,6 +155,14 @@ void SimSetup::createSim(void){ // initialize the Fortran initFortran(); + + if (globals->haveMinimizer()) + // make minimizer + makeMinimizer(); + else + // make the integrator + makeIntegrator(); + } @@ -262,10 +263,10 @@ void SimSetup::makeMolecules(void){ else{ molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration()); + } molInfo.myAtoms[j]->setType(currentAtom->getType()); - #ifdef IS_MPI molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]); @@ -409,6 +410,9 @@ void SimSetup::makeMolecules(void){ info[k].excludes->addPair(exK, exL); } + + molInfo.myRigidBodies.clear(); + for (j = 0; j < molInfo.nRigidBodies; j++){ currentRigidBody = comp_stamps[stampID]->getRigidBody(j); @@ -418,7 +422,7 @@ void SimSetup::makeMolecules(void){ myRB = new RigidBody(); - sprintf(rbName,"%s_RB_%s", molName, j); + sprintf(rbName,"%s_RB_%d", molName, j); myRB->setType(rbName); for (rb1 = 0; rb1 < nMembers; rb1++) { @@ -460,11 +464,11 @@ void SimSetup::makeMolecules(void){ // used for the exclude list: #ifdef IS_MPI - exI = info[k].atoms[tempI]->getGlobalIndex() + 1; - exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; + exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1; + exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1; #else - exI = tempI + 1; - exJ = tempJ + 1; + exI = molInfo.myAtoms[tempI]->getIndex() + 1; + exJ = molInfo.myAtoms[tempJ]->getIndex() + 1; #endif info[k].excludes->addPair(exI, exJ); @@ -479,6 +483,9 @@ void SimSetup::makeMolecules(void){ // After this is all set up, scan through the atoms to // see if they can be added to the integrableObjects: + + molInfo.myIntegrableObjects.clear(); + for (j = 0; j < molInfo.nAtoms; j++){ @@ -531,13 +538,13 @@ void SimSetup::makeMolecules(void){ // clean up the forcefield - if (!globals->haveLJrcut()){ + if (!globals->haveRcut()){ the_ff->calcRcut(); } else { - the_ff->setRcut( globals->getLJrcut() ); + the_ff->setRcut( globals->getRcut() ); } the_ff->cleanMe(); @@ -827,7 +834,7 @@ void SimSetup::gatherInfo(void){ } //check whether sample time, status time, thermal time and reset time are divisble by dt - if (!isDivisible(globals->getSampleTime(), globals->getDt())){ + if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ sprintf(painCave.errMsg, "Sample time is not divisible by dt.\n" "\tThis will result in samples that are not uniformly\n" @@ -837,7 +844,7 @@ void SimSetup::gatherInfo(void){ simError(); } - if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ + if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){ sprintf(painCave.errMsg, "Status time is not divisible by dt.\n" "\tThis will result in status reports that are not uniformly\n" @@ -873,12 +880,10 @@ void SimSetup::gatherInfo(void){ if (globals->haveSampleTime()){ info[i].sampleTime = globals->getSampleTime(); info[i].statusTime = info[i].sampleTime; - info[i].thermalTime = info[i].sampleTime; } else{ info[i].sampleTime = globals->getRunTime(); info[i].statusTime = info[i].sampleTime; - info[i].thermalTime = info[i].sampleTime; } if (globals->haveStatusTime()){ @@ -887,6 +892,8 @@ void SimSetup::gatherInfo(void){ if (globals->haveThermalTime()){ info[i].thermalTime = globals->getThermalTime(); + } else { + info[i].thermalTime = globals->getRunTime(); } info[i].resetIntegrator = 0; @@ -957,6 +964,7 @@ void SimSetup::finalInfoCheck(void){ void SimSetup::finalInfoCheck(void){ int index; int usesDipoles; + int usesCharges; int i; for (i = 0; i < nInfo; i++){ @@ -968,45 +976,49 @@ void SimSetup::finalInfoCheck(void){ usesDipoles = (info[i].atoms[index])->hasDipole(); index++; } - + index = 0; + usesCharges = 0; + while ((index < info[i].n_atoms) && !usesCharges){ + usesCharges= (info[i].atoms[index])->hasCharge(); + index++; + } #ifdef IS_MPI int myUse = usesDipoles; 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->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"); + "\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, @@ -1019,35 +1031,36 @@ void SimSetup::finalInfoCheck(void){ info[i].dielectric = globals->getDielectric(); } else{ - if (usesDipoles){ - if (!globals->haveECR()){ + if (usesDipoles || usesCharges){ + + 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(); + 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].setDefaultRcut(theRcut, theRsw); - info[i].setDefaultEcr(theEcr, theEst); } } } @@ -1352,14 +1365,13 @@ void SimSetup::mpiMolDivide(void){ int localMol, allMol; int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; int local_rigid; - vector globalAtomIndex; vector globalMolIndex; mpiSim = new mpiSimulation(info); mpiSim->divideLabor(); globalAtomIndex = mpiSim->getGlobalAtomIndex(); - globalMolIndex = mpiSim->getGlobalMolIndex(); + //globalMolIndex = mpiSim->getGlobalMolIndex(); // set up the local variables