--- trunk/OOPSE/libmdtools/SimSetup.cpp 2004/04/12 20:32:20 1097 +++ 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,12 +155,20 @@ void SimSetup::createSim(void){ // initialize the Fortran initFortran(); + + if (globals->haveMinimizer()) + // make minimizer + makeMinimizer(); + else + // make the integrator + makeIntegrator(); + } void SimSetup::makeMolecules(void){ int i, j, k; - int exI, exJ, exK, exL, slI; + int exI, exJ, exK, exL, slI, slJ; int tempI, tempJ, tempK, tempL; int molI; int stampID, atomOffset, rbOffset; @@ -190,6 +191,8 @@ void SimSetup::makeMolecules(void){ set skipList; double phi, theta, psi; + char* molName; + char rbName[100]; //init the forceField paramters @@ -206,6 +209,7 @@ void SimSetup::makeMolecules(void){ for (i = 0; i < info[k].n_mol; i++){ stampID = info[k].molecules[i].getStampID(); + molName = comp_stamps[stampID]->getID(); molInfo.nAtoms = comp_stamps[stampID]->getNAtoms(); molInfo.nBonds = comp_stamps[stampID]->getNBonds(); @@ -259,13 +263,13 @@ 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(globalIndex[j + atomOffset]); + molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]); #endif // is_mpi } @@ -406,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); @@ -414,6 +421,9 @@ void SimSetup::makeMolecules(void){ // Create the Rigid Body: myRB = new RigidBody(); + + sprintf(rbName,"%s_RB_%d", molName, j); + myRB->setType(rbName); for (rb1 = 0; rb1 < nMembers; rb1++) { @@ -454,18 +464,54 @@ 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); } + } + + molInfo.myRigidBodies.push_back(myRB); + info[k].rigidBodies.push_back(myRB); + } + + + // 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++){ + +#ifdef IS_MPI + slJ = molInfo.myAtoms[j]->getGlobalIndex(); +#else + slJ = j+atomOffset; +#endif + + // if they aren't on the skip list, then they can be integrated + + if (skipList.find(slJ) == skipList.end()) { + mySD = (StuntDouble *) molInfo.myAtoms[j]; + info[k].integrableObjects.push_back(mySD); + molInfo.myIntegrableObjects.push_back(mySD); } } + + // all rigid bodies are integrated: + + for (j = 0; j < molInfo.nRigidBodies; j++) { + mySD = (StuntDouble *) molInfo.myRigidBodies[j]; + info[k].integrableObjects.push_back(mySD); + molInfo.myIntegrableObjects.push_back(mySD); + } + // send the arrays off to the forceField for init. @@ -482,28 +528,7 @@ void SimSetup::makeMolecules(void){ delete[] theBonds; delete[] theBends; delete[] theTorsions; - } - - // build up the integrableObjects vector: - - for (i = 0; i < info[k].n_atoms; i++) { - -#ifdef IS_MPI - slI = info[k].atoms[i]->getGlobalIndex(); -#else - slI = i; -#endif - - if (skipList.find(slI) == skipList.end()) { - mySD = (StuntDouble *) info[k].atoms[i]; - info[k].integrableObjects.push_back(mySD); - } - } - for (i = 0; i < info[k].rigidBodies.size(); i++) { - mySD = (StuntDouble *) info[k].rigidBodies[i]; - info[k].integrableObjects.push_back(mySD); - } - + } } #ifdef IS_MPI @@ -513,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(); @@ -809,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" @@ -819,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" @@ -855,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()){ @@ -869,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; @@ -939,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++){ @@ -950,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, @@ -1001,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; - } - else{ - theEcr = globals->getECR(); - } - - if (!globals->haveEST()){ + "\tfor the cutoffRadius.\n"); + painCave.isFatal = 0; + simError(); + theRcut = 15.0; + } + else{ + 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); } } } @@ -1334,10 +1365,13 @@ void SimSetup::mpiMolDivide(void){ int localMol, allMol; int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; int local_rigid; + vector globalMolIndex; mpiSim = new mpiSimulation(info); - globalIndex = mpiSim->divideLabor(); + mpiSim->divideLabor(); + globalAtomIndex = mpiSim->getGlobalAtomIndex(); + //globalMolIndex = mpiSim->getGlobalMolIndex(); // set up the local variables @@ -1351,7 +1385,7 @@ void SimSetup::mpiMolDivide(void){ local_bends = 0; local_torsions = 0; local_rigid = 0; - globalAtomIndex = 0; + globalAtomCounter = 0; for (i = 0; i < n_components; i++){ for (j = 0; j < components_nmol[i]; j++){ @@ -1364,8 +1398,8 @@ void SimSetup::mpiMolDivide(void){ localMol++; } for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ - info[0].molMembershipArray[globalAtomIndex] = allMol; - globalAtomIndex++; + info[0].molMembershipArray[globalAtomCounter] = allMol; + globalAtomCounter++; } allMol++; @@ -1433,15 +1467,15 @@ void SimSetup::makeSysArrays(void){ #else // is_mpi molIndex = 0; - globalAtomIndex = 0; + globalAtomCounter = 0; for (i = 0; i < n_components; i++){ for (j = 0; j < components_nmol[i]; j++){ the_molecules[molIndex].setStampID(i); the_molecules[molIndex].setMyIndex(molIndex); the_molecules[molIndex].setGlobalIndex(molIndex); for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ - info[l].molMembershipArray[globalAtomIndex] = molIndex; - globalAtomIndex++; + info[l].molMembershipArray[globalAtomCounter] = molIndex; + globalAtomCounter++; } molIndex++; }