--- trunk/OOPSE/libmdtools/SimSetup.cpp 2003/09/02 14:30:12 738 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2004/05/12 14:30:12 1163 @@ -1,15 +1,17 @@ #include -#include +#include #include -#include +#include #include #include - #include "SimSetup.hpp" #include "ReadWrite.hpp" #include "parse_me.h" #include "Integrator.hpp" #include "simError.h" +#include "RigidBody.hpp" +//#include "ConjugateMinimizer.hpp" +#include "OOPSEMinimizer.hpp" #ifdef IS_MPI #include "mpiBASS.h" @@ -22,16 +24,43 @@ #define NVT_ENS 1 #define NPTi_ENS 2 #define NPTf_ENS 3 -#define NPTim_ENS 4 -#define NPTfm_ENS 5 +#define NPTxyz_ENS 4 -#define FF_DUFF 0 -#define FF_LJ 1 -#define FF_EAM 2 +#define FF_DUFF 0 +#define FF_LJ 1 +#define FF_EAM 2 +#define FF_H2O 3 + using namespace std; +/** + * Check whether dividend is divisble by divisor or not + */ +bool isDivisible(double dividend, double divisor){ + double tolerance = 0.000001; + double quotient; + double diff; + int intQuotient; + + quotient = dividend / divisor; + + if (quotient < 0) + quotient = -quotient; + + intQuotient = int (quotient + tolerance); + + diff = fabs(fabs(dividend) - intQuotient * fabs(divisor)); + + if (diff <= tolerance) + return true; + else + return false; +} + SimSetup::SimSetup(){ + + initSuspend = false; isInfoArray = 0; nInfo = 1; @@ -54,6 +83,7 @@ void SimSetup::setSimInfo(SimInfo* the_info, int theNi info = the_info; nInfo = theNinfo; isInfoArray = 1; + initSuspend = true; } @@ -92,7 +122,6 @@ void SimSetup::createSim(void){ #endif // is_mpi void SimSetup::createSim(void){ - int i, j, k, globalAtomIndex; // gather all of the information from the Bass file @@ -108,18 +137,17 @@ void SimSetup::createSim(void){ // initialize the system coordinates - if (!isInfoArray){ + if ( !initSuspend ){ initSystemCoords(); + + if( !(globals->getUseInitTime()) ) + info[0].currentTime = 0.0; } // make the output filenames makeOutNames(); - - // make the integrator - - makeIntegrator(); - + #ifdef IS_MPI mpiSim->mpiRefresh(); #endif @@ -127,93 +155,124 @@ void SimSetup::createSim(void){ // initialize the Fortran initFortran(); + + if (globals->haveMinimizer()) + // make minimizer + makeMinimizer(); + else + // make the integrator + makeIntegrator(); + } void SimSetup::makeMolecules(void){ - int k, l; - int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; + int i, j, k; + int exI, exJ, exK, exL, slI, slJ; + int tempI, tempJ, tempK, tempL; + int molI; + int stampID, atomOffset, rbOffset; molInit molInfo; DirectionalAtom* dAtom; + RigidBody* myRB; + StuntDouble* mySD; LinkedAssign* extras; LinkedAssign* current_extra; AtomStamp* currentAtom; BondStamp* currentBond; BendStamp* currentBend; TorsionStamp* currentTorsion; - + RigidBodyStamp* currentRigidBody; + CutoffGroupStamp* currentCutoffGroup; + CutoffGroup* myCutoffGroup; + bond_pair* theBonds; bend_set* theBends; torsion_set* theTorsions; + set skipList; + double phi, theta, psi; + char* molName; + char rbName[100]; + //init the forceField paramters the_ff->readParams(); - // init the atoms - double ux, uy, uz, u, uSqr; + int nMembers, nNew, rb1, rb2; for (k = 0; k < nInfo; k++){ the_ff->setSimInfo(&(info[k])); atomOffset = 0; - excludeOffset = 0; + 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(); molInfo.nBends = comp_stamps[stampID]->getNBends(); molInfo.nTorsions = comp_stamps[stampID]->getNTorsions(); - molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions; - + molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies(); + molInfo.nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups(); + molInfo.myAtoms = &(info[k].atoms[atomOffset]); - molInfo.myExcludes = &(info[k].excludes[excludeOffset]); - molInfo.myBonds = new Bond * [molInfo.nBonds]; - molInfo.myBends = new Bend * [molInfo.nBends]; - molInfo.myTorsions = new Torsion * [molInfo.nTorsions]; + if (molInfo.nBonds > 0) + molInfo.myBonds = new (Bond *) [molInfo.nBonds]; + else + molInfo.myBonds = NULL; + + if (molInfo.nBends > 0) + molInfo.myBends = new (Bend *) [molInfo.nBends]; + else + molInfo.myBends = NULL; + + if (molInfo.nTorsions > 0) + molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions]; + else + molInfo.myTorsions = NULL; + theBonds = new bond_pair[molInfo.nBonds]; theBends = new bend_set[molInfo.nBends]; theTorsions = new torsion_set[molInfo.nTorsions]; - + // make the Atoms for (j = 0; j < molInfo.nAtoms; j++){ currentAtom = comp_stamps[stampID]->getAtom(j); + if (currentAtom->haveOrientation()){ dAtom = new DirectionalAtom((j + atomOffset), info[k].getConfiguration()); info[k].n_oriented++; molInfo.myAtoms[j] = dAtom; - ux = currentAtom->getOrntX(); - uy = currentAtom->getOrntY(); - uz = currentAtom->getOrntZ(); + // Directional Atoms have standard unit vectors which are oriented + // in space using the three Euler angles. We assume the standard + // unit vector was originally along the z axis below. - uSqr = (ux * ux) + (uy * uy) + (uz * uz); - - u = sqrt(uSqr); - ux = ux / u; - uy = uy / u; - uz = uz / u; - - dAtom->setSUx(ux); - dAtom->setSUy(uy); - dAtom->setSUz(uz); + phi = currentAtom->getEulerPhi() * M_PI / 180.0; + theta = currentAtom->getEulerTheta() * M_PI / 180.0; + psi = currentAtom->getEulerPsi()* M_PI / 180.0; + + dAtom->setUnitFrameFromEuler(phi, theta, psi); + } else{ - molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset), - info[k].getConfiguration()); + + molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration()); + } - molInfo.myAtoms[j]->setType(currentAtom->getType()); + 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 } @@ -224,28 +283,19 @@ void SimSetup::makeMolecules(void){ theBonds[j].a = currentBond->getA() + atomOffset; theBonds[j].b = currentBond->getB() + atomOffset; - exI = theBonds[j].a; - exJ = theBonds[j].b; + tempI = theBonds[j].a; + tempJ = theBonds[j].b; - // exclude_I must always be the smaller of the pair - if (exI > exJ){ - tempEx = exI; - exI = exJ; - exJ = tempEx; - } #ifdef IS_MPI - tempEx = exI; - exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; - tempEx = exJ; - exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; + exI = info[k].atoms[tempI]->getGlobalIndex() + 1; + exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; +#else + exI = tempI + 1; + exJ = tempJ + 1; +#endif - info[k].excludes[j + excludeOffset]->setPair(exI, exJ); -#else // isn't MPI - - info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); -#endif //is_mpi + info[k].excludes->addPair(exI, exJ); } - excludeOffset += molInfo.nBonds; //make the bends for (j = 0; j < molInfo.nBends; j++){ @@ -295,33 +345,41 @@ void SimSetup::makeMolecules(void){ } } - if (!theBends[j].isGhost){ - exI = theBends[j].a; - exJ = theBends[j].c; - } - else{ - exI = theBends[j].a; - exJ = theBends[j].b; - } - - // exclude_I must always be the smaller of the pair - if (exI > exJ){ - tempEx = exI; - exI = exJ; - exJ = tempEx; - } + if (theBends[j].isGhost) { + + tempI = theBends[j].a; + tempJ = theBends[j].b; + #ifdef IS_MPI - tempEx = exI; - exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; - tempEx = exJ; - exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; + exI = info[k].atoms[tempI]->getGlobalIndex() + 1; + exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; +#else + exI = tempI + 1; + exJ = tempJ + 1; +#endif + info[k].excludes->addPair(exI, exJ); - info[k].excludes[j + excludeOffset]->setPair(exI, exJ); -#else // isn't MPI - info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); -#endif //is_mpi + } else { + + tempI = theBends[j].a; + tempJ = theBends[j].b; + tempK = theBends[j].c; + +#ifdef IS_MPI + exI = info[k].atoms[tempI]->getGlobalIndex() + 1; + exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; + exK = info[k].atoms[tempK]->getGlobalIndex() + 1; +#else + exI = tempI + 1; + exJ = tempJ + 1; + exK = tempK + 1; +#endif + + info[k].excludes->addPair(exI, exK); + info[k].excludes->addPair(exI, exJ); + info[k].excludes->addPair(exJ, exK); + } } - excludeOffset += molInfo.nBends; for (j = 0; j < molInfo.nTorsions; j++){ currentTorsion = comp_stamps[stampID]->getTorsion(j); @@ -330,38 +388,167 @@ void SimSetup::makeMolecules(void){ theTorsions[j].c = currentTorsion->getC() + atomOffset; theTorsions[j].d = currentTorsion->getD() + atomOffset; - exI = theTorsions[j].a; - exJ = theTorsions[j].d; + tempI = theTorsions[j].a; + tempJ = theTorsions[j].b; + tempK = theTorsions[j].c; + tempL = theTorsions[j].d; - // exclude_I must always be the smaller of the pair - if (exI > exJ){ - tempEx = exI; - exI = exJ; - exJ = tempEx; - } #ifdef IS_MPI - tempEx = exI; - exI = info[k].atoms[tempEx]->getGlobalIndex() + 1; - tempEx = exJ; - exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1; - - info[k].excludes[j + excludeOffset]->setPair(exI, exJ); -#else // isn't MPI - info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1)); -#endif //is_mpi - } - excludeOffset += molInfo.nTorsions; + exI = info[k].atoms[tempI]->getGlobalIndex() + 1; + exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1; + exK = info[k].atoms[tempK]->getGlobalIndex() + 1; + exL = info[k].atoms[tempL]->getGlobalIndex() + 1; +#else + exI = tempI + 1; + exJ = tempJ + 1; + exK = tempK + 1; + exL = tempL + 1; +#endif + info[k].excludes->addPair(exI, exJ); + info[k].excludes->addPair(exI, exK); + info[k].excludes->addPair(exI, exL); + info[k].excludes->addPair(exJ, exK); + info[k].excludes->addPair(exJ, exL); + info[k].excludes->addPair(exK, exL); + } - // send the arrays off to the forceField for init. + + molInfo.myRigidBodies.clear(); + + for (j = 0; j < molInfo.nRigidBodies; j++){ + currentRigidBody = comp_stamps[stampID]->getRigidBody(j); + nMembers = currentRigidBody->getNMembers(); + + // Create the Rigid Body: + + myRB = new RigidBody(); + + sprintf(rbName,"%s_RB_%d", molName, j); + myRB->setType(rbName); + + for (rb1 = 0; rb1 < nMembers; rb1++) { + + // molI is atom numbering inside this molecule + molI = currentRigidBody->getMember(rb1); + + // tempI is atom numbering on local processor + tempI = molI + atomOffset; + + // currentAtom is the AtomStamp (which we need for + // rigid body reference positions) + currentAtom = comp_stamps[stampID]->getAtom(molI); + + // When we add to the rigid body, add the atom itself and + // the stamp info: + + myRB->addAtom(info[k].atoms[tempI], currentAtom); + + // Add this atom to the Skip List for the integrators +#ifdef IS_MPI + slI = info[k].atoms[tempI]->getGlobalIndex(); +#else + slI = tempI; +#endif + skipList.insert(slI); + + } + + for(rb1 = 0; rb1 < nMembers - 1; rb1++) { + for(rb2 = rb1+1; rb2 < nMembers; rb2++) { + + tempI = currentRigidBody->getMember(rb1); + tempJ = currentRigidBody->getMember(rb2); + + // Some explanation is required here. + // Fortran indexing starts at 1, while c indexing starts at 0 + // Also, in parallel computations, the GlobalIndex is + // used for the exclude list: + +#ifdef IS_MPI + exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1; + exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1; +#else + 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); + } + + + //create cutoff group for molecule + molInfo.myCutoffGroups.clear(); + for (j = 0; j < molInfo.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]); + } + + molInfo.myCutoffGroups.push_back(myCutoffGroup); + }//end for (j = 0; j < molInfo.nCutoffGroups; j++) + + + + // 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. + the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds); the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends); the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions, theTorsions); - info[k].molecules[i].initialize(molInfo); @@ -369,7 +556,7 @@ void SimSetup::makeMolecules(void){ delete[] theBonds; delete[] theBends; delete[] theTorsions; - } + } } #ifdef IS_MPI @@ -377,10 +564,6 @@ void SimSetup::makeMolecules(void){ MPIcheckPoint(); #endif // is_mpi - // clean up the forcefield - - the_ff->calcRcut(); - the_ff->cleanMe(); } void SimSetup::initFromBass(void){ @@ -553,7 +736,7 @@ void SimSetup::gatherInfo(void){ void SimSetup::gatherInfo(void){ - int i, j, k; + int i; ensembleCase = -1; ffCase = -1; @@ -581,6 +764,9 @@ void SimSetup::gatherInfo(void){ else if (!strcasecmp(force_field, "EAM")){ ffCase = FF_EAM; } + else if (!strcasecmp(force_field, "WATER")){ + ffCase = FF_H2O; + } else{ sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n", force_field); @@ -604,16 +790,13 @@ void SimSetup::gatherInfo(void){ else if (!strcasecmp(ensemble, "NPTf")){ ensembleCase = NPTf_ENS; } - else if (!strcasecmp(ensemble, "NPTim")){ - ensembleCase = NPTim_ENS; + else if (!strcasecmp(ensemble, "NPTxyz")){ + ensembleCase = NPTxyz_ENS; } - else if (!strcasecmp(ensemble, "NPTfm")){ - ensembleCase = NPTfm_ENS; - } else{ sprintf(painCave.errMsg, - "SimSetup Warning. Unrecognized Ensemble -> %s, " - "reverting to NVE for this simulation.\n", + "SimSetup Warning. Unrecognized Ensemble -> %s \n" + "\treverting to NVE for this simulation.\n", ensemble); painCave.isFatal = 0; simError(); @@ -645,8 +828,8 @@ void SimSetup::gatherInfo(void){ if (!the_components[i]->haveNMol()){ // we have a problem sprintf(painCave.errMsg, - "SimSetup Error. No global NMol or component NMol" - " given. Cannot calculate the number of atoms.\n"); + "SimSetup Error. No global NMol or component NMol given.\n" + "\tCannot calculate the number of atoms.\n"); painCave.isFatal = 1; simError(); } @@ -664,20 +847,59 @@ void SimSetup::gatherInfo(void){ " Please give nMol in the components.\n"); painCave.isFatal = 1; simError(); + } + + //check whether sample time, status time, thermal time and reset time are divisble by dt + 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" + "\tdistributed in time. If this is a problem, change\n" + "\tyour sampleTime variable.\n"); + painCave.isFatal = 0; + simError(); + } + + 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" + "\tdistributed in time. If this is a problem, change \n" + "\tyour statusTime variable.\n"); + painCave.isFatal = 0; + simError(); } + if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){ + sprintf(painCave.errMsg, + "Thermal time is not divisible by dt.\n" + "\tThis will result in thermalizations that are not uniformly\n" + "\tdistributed in time. If this is a problem, change \n" + "\tyour thermalTime variable.\n"); + painCave.isFatal = 0; + simError(); + } + + if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){ + sprintf(painCave.errMsg, + "Reset time is not divisible by dt.\n" + "\tThis will result in integrator resets that are not uniformly\n" + "\tdistributed in time. If this is a problem, change\n" + "\tyour resetTime variable.\n"); + painCave.isFatal = 0; + simError(); + } + // set the status, sample, and thermal kick times for (i = 0; i < nInfo; i++){ 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()){ @@ -686,61 +908,28 @@ void SimSetup::gatherInfo(void){ if (globals->haveThermalTime()){ info[i].thermalTime = globals->getThermalTime(); + } else { + info[i].thermalTime = globals->getRunTime(); } - // check for the temperature set flag + info[i].resetIntegrator = 0; + if( globals->haveResetTime() ){ + info[i].resetTime = globals->getResetTime(); + info[i].resetIntegrator = 1; + } + // check for the temperature set flag + if (globals->haveTempSet()) info[i].setTemp = globals->getTempSet(); - // get some of the tricky things that may still be in the globals + // check for the extended State init - double boxVector[3]; - if (globals->haveBox()){ - boxVector[0] = globals->getBox(); - boxVector[1] = globals->getBox(); - boxVector[2] = globals->getBox(); - - info[i].setBox(boxVector); - } - else if (globals->haveDensity()){ - double vol; - vol = (double) tot_nmol / globals->getDensity(); - boxVector[0] = pow(vol, (1.0 / 3.0)); - boxVector[1] = boxVector[0]; - boxVector[2] = boxVector[0]; - - info[i].setBox(boxVector); - } - else{ - if (!globals->haveBoxX()){ - sprintf(painCave.errMsg, - "SimSetup error, no periodic BoxX size given.\n"); - painCave.isFatal = 1; - simError(); - } - boxVector[0] = globals->getBoxX(); - - if (!globals->haveBoxY()){ - sprintf(painCave.errMsg, - "SimSetup error, no periodic BoxY size given.\n"); - painCave.isFatal = 1; - simError(); - } - boxVector[1] = globals->getBoxY(); - - if (!globals->haveBoxZ()){ - sprintf(painCave.errMsg, - "SimSetup error, no periodic BoxZ size given.\n"); - painCave.isFatal = 1; - simError(); - } - boxVector[2] = globals->getBoxZ(); - - info[i].setBox(boxVector); - } + info[i].useInitXSstate = globals->getUseInitXSstate(); + info[i].orthoTolerance = globals->getOrthoBoxTolerance(); + } - + //setup seed for random number generator int seedValue; @@ -780,9 +969,9 @@ void SimSetup::gatherInfo(void){ for (int i = 0; i < nInfo; i++){ info[i].setSeed(seedValue); } - + #ifdef IS_MPI - strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n"); + strcpy(checkPointMsg, "Successfully gathered all information from Bass\n"); MPIcheckPoint(); #endif // is_mpi } @@ -791,6 +980,7 @@ void SimSetup::finalInfoCheck(void){ void SimSetup::finalInfoCheck(void){ int index; int usesDipoles; + int usesCharges; int i; for (i = 0; i < nInfo; i++){ @@ -802,104 +992,126 @@ 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->haveRcut()) { + theRcut = globals->getRcut(); + + 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->haveECR()){ + + if (!globals->haveRcut()){ sprintf(painCave.errMsg, - "SimSetup Warning: using default value of 1/2 the smallest " - "box length for the electrostaticCutoffRadius.\n" - "I hope you have a very fast processor!\n"); + "SimSetup Warning: No value was set for the cutoffRadius.\n" + "\tOOPSE will use a default value of 15.0 angstroms" + "\tfor the cutoffRadius.\n"); painCave.isFatal = 0; simError(); - double smallest; - smallest = info[i].boxL[0]; - if (info[i].boxL[1] <= smallest) - smallest = info[i].boxL[1]; - if (info[i].boxL[2] <= smallest) - smallest = info[i].boxL[2]; - theEcr = 0.5 * smallest; + theRcut = 15.0; } else{ - theEcr = globals->getECR(); + theRcut = globals->getRcut(); } - if (!globals->haveEST()){ + if (!globals->haveRsw()){ sprintf(painCave.errMsg, - "SimSetup Warning: using default value of 0.05 * the " - "electrostaticCutoffRadius for the electrostaticSkinThickness\n"); + "SimSetup Warning: No value was set for switchingRadius.\n" + "\tOOPSE will use a default value of\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].setEcr(theEcr, theEst); + info[i].setDefaultRcut(theRcut, theRsw); if (!globals->haveDielectric()){ sprintf(painCave.errMsg, - "SimSetup Error: You are trying to use Reaction Field without" - "setting a dielectric constant!\n"); + "SimSetup Error: No Dielectric constant was set.\n" + "\tYou are trying to use Reaction Field without" + "\tsetting a dielectric constant!\n"); painCave.isFatal = 1; simError(); } info[i].dielectric = globals->getDielectric(); } else{ - if (usesDipoles){ - if (!globals->haveECR()){ + if (usesDipoles || usesCharges){ + + if (!globals->haveRcut()){ sprintf(painCave.errMsg, - "SimSetup Warning: using default value of 1/2 the smallest " - "box length for the electrostaticCutoffRadius.\n" - "I hope you have a very fast processor!\n"); + "SimSetup Warning: No value was set for the cutoffRadius.\n" + "\tOOPSE will use a default value of 15.0 angstroms" + "\tfor the cutoffRadius.\n"); painCave.isFatal = 0; simError(); - double smallest; - smallest = info[i].boxL[0]; - if (info[i].boxL[1] <= smallest) - smallest = info[i].boxL[1]; - if (info[i].boxL[2] <= smallest) - smallest = info[i].boxL[2]; - theEcr = 0.5 * smallest; - } + theRcut = 15.0; + } else{ - theEcr = globals->getECR(); + theRcut = globals->getRcut(); } - - if (!globals->haveEST()){ + + if (!globals->haveRsw()){ sprintf(painCave.errMsg, - "SimSetup Warning: using default value of 0.05 * the " - "electrostaticCutoffRadius for the " - "electrostaticSkinThickness\n"); + "SimSetup Warning: No value was set for switchingRadius.\n" + "\tOOPSE will use a default value of\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].setEcr(theEcr, theEst); + + info[i].setDefaultRcut(theRcut, theRsw); + } } } - #ifdef IS_MPI strcpy(checkPointMsg, "post processing checks out"); MPIcheckPoint(); #endif // is_mpi -} + // clean up the forcefield + the_ff->cleanMe(); +} + void SimSetup::initSystemCoords(void){ int i; @@ -916,7 +1128,6 @@ void SimSetup::initSystemCoords(void){ if (worldRank == 0){ #endif //is_mpi inName = globals->getInitialConfig(); - double* tempDouble = new double[1000000]; fileInit = new InitializeFromFile(inName); #ifdef IS_MPI } @@ -928,21 +1139,14 @@ void SimSetup::initSystemCoords(void){ delete fileInit; } else{ -#ifdef IS_MPI - + // no init from bass - - sprintf(painCave.errMsg, - "Cannot intialize a parallel simulation without an initial configuration file.\n"); - painCave.isFatal; + + sprintf(painCave.errMsg, + "Cannot intialize a simulation without an initial configuration file.\n"); + painCave.isFatal = 1;; simError(); - -#else - - initFromBass(); - - -#endif + } #ifdef IS_MPI @@ -1096,6 +1300,10 @@ void SimSetup::createFF(void){ the_ff = new EAM_FF(); break; + case FF_H2O: + the_ff = new WATER(); + break; + default: sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field in case statement.\n"); @@ -1116,7 +1324,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 @@ -1151,7 +1362,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"); @@ -1160,7 +1377,7 @@ void SimSetup::calcSysValues(void){ } void SimSetup::calcSysValues(void){ - int i, j, k; + int i; int* molMembershipArray; @@ -1168,13 +1385,15 @@ void SimSetup::calcSysValues(void){ tot_bonds = 0; tot_bends = 0; tot_torsions = 0; + tot_rigid = 0; for (i = 0; i < n_components; i++){ tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); + tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies(); } - + tot_SRI = tot_bonds + tot_bends + tot_torsions; molMembershipArray = new int[tot_atoms]; @@ -1196,10 +1415,14 @@ void SimSetup::mpiMolDivide(void){ int i, j, k; 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 @@ -1212,9 +1435,9 @@ void SimSetup::mpiMolDivide(void){ local_bonds = 0; local_bends = 0; local_torsions = 0; - globalAtomIndex = 0; + local_rigid = 0; + globalAtomCounter = 0; - for (i = 0; i < n_components; i++){ for (j = 0; j < components_nmol[i]; j++){ if (mol2proc[allMol] == worldRank){ @@ -1222,11 +1445,12 @@ void SimSetup::mpiMolDivide(void){ local_bonds += comp_stamps[i]->getNBonds(); local_bends += comp_stamps[i]->getNBends(); local_torsions += comp_stamps[i]->getNTorsions(); + local_rigid += comp_stamps[i]->getNRigidBodies(); localMol++; } for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ - info[0].molMembershipArray[globalAtomIndex] = allMol; - globalAtomIndex++; + info[0].molMembershipArray[globalAtomCounter] = allMol; + globalAtomCounter++; } allMol++; @@ -1235,11 +1459,12 @@ void SimSetup::mpiMolDivide(void){ local_SRI = local_bonds + local_bends + local_torsions; info[0].n_atoms = mpiSim->getMyNlocal(); + if (local_atoms != info[0].n_atoms){ sprintf(painCave.errMsg, - "SimSetup error: mpiSim's localAtom (%d) and SimSetup's" - " localAtom (%d) are not equal.\n", + "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" + "\tlocalAtom (%d) are not equal.\n", info[0].n_atoms, local_atoms); painCave.isFatal = 1; simError(); @@ -1259,13 +1484,15 @@ void SimSetup::makeSysArrays(void){ void SimSetup::makeSysArrays(void){ - int i, j, k, l; + +#ifndef IS_MPI + int k, j; +#endif // is_mpi + int i, l; Atom** the_atoms; Molecule* the_molecules; - Exclude** the_excludes; - for (l = 0; l < nInfo; l++){ // create the atom and short range interaction arrays @@ -1291,15 +1518,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++; } @@ -1308,33 +1535,15 @@ void SimSetup::makeSysArrays(void){ #endif // is_mpi - - if (info[l].n_SRI){ - Exclude::createArray(info[l].n_SRI); - the_excludes = new Exclude * [info[l].n_SRI]; - for (int ex = 0; ex < info[l].n_SRI; ex++){ - the_excludes[ex] = new Exclude(ex); - } - info[l].globalExcludes = new int; - info[l].n_exclude = info[l].n_SRI; - } - else{ - Exclude::createArray(1); - the_excludes = new Exclude * ; - the_excludes[0] = new Exclude(0); - the_excludes[0]->setPair(0, 0); - info[l].globalExcludes = new int; - info[l].globalExcludes[0] = 0; - info[l].n_exclude = 0; - } - + info[l].globalExcludes = new int; + info[l].globalExcludes[0] = 0; + // set the arrays into the SimInfo object info[l].atoms = the_atoms; info[l].molecules = the_molecules; info[l].nGlobalExcludes = 0; - info[l].excludes = the_excludes; - + the_ff->setSimInfo(info); } } @@ -1342,21 +1551,24 @@ void SimSetup::makeIntegrator(void){ void SimSetup::makeIntegrator(void){ int k; + NVE* myNVE = NULL; NVT* myNVT = NULL; - NPTi* myNPTi = NULL; - NPTf* myNPTf = NULL; - NPTim* myNPTim = NULL; - NPTfm* myNPTfm = NULL; + NPTi >* myNPTi = NULL; + NPTf >* myNPTf = NULL; + NPTxyz >* myNPTxyz = NULL; for (k = 0; k < nInfo; k++){ switch (ensembleCase){ case NVE_ENS: if (globals->haveZconstraints()){ setupZConstraint(info[k]); - new ZConstraint >(&(info[k]), the_ff); + myNVE = new ZConstraint >(&(info[k]), the_ff); } - else - new NVE(&(info[k]), the_ff); + else{ + myNVE = new NVE(&(info[k]), the_ff); + } + + info->the_integrator = myNVE; break; case NVT_ENS: @@ -1374,19 +1586,21 @@ void SimSetup::makeIntegrator(void){ else{ sprintf(painCave.errMsg, "SimSetup error: If you use the NVT\n" - " ensemble, you must set tauThermostat.\n"); + "\tensemble, you must set tauThermostat.\n"); painCave.isFatal = 1; simError(); } + + info->the_integrator = myNVT; break; case NPTi_ENS: if (globals->haveZconstraints()){ setupZConstraint(info[k]); - myNPTi = new ZConstraint >(&(info[k]), the_ff); + myNPTi = new ZConstraint > >(&(info[k]), the_ff); } else - myNPTi = new NPTi(&(info[k]), the_ff); + myNPTi = new NPTi >(&(info[k]), the_ff); myNPTi->setTargetTemp(globals->getTargetTemp()); @@ -1395,7 +1609,7 @@ void SimSetup::makeIntegrator(void){ else{ sprintf(painCave.errMsg, "SimSetup error: If you use a constant pressure\n" - " ensemble, you must set targetPressure in the BASS file.\n"); + "\tensemble, you must set targetPressure in the BASS file.\n"); painCave.isFatal = 1; simError(); } @@ -1405,7 +1619,7 @@ void SimSetup::makeIntegrator(void){ else{ sprintf(painCave.errMsg, "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauThermostat.\n"); + "\tensemble, you must set tauThermostat.\n"); painCave.isFatal = 1; simError(); } @@ -1415,19 +1629,21 @@ void SimSetup::makeIntegrator(void){ else{ sprintf(painCave.errMsg, "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauBarostat.\n"); + "\tensemble, you must set tauBarostat.\n"); painCave.isFatal = 1; simError(); } + + info->the_integrator = myNPTi; break; case NPTf_ENS: if (globals->haveZconstraints()){ setupZConstraint(info[k]); - myNPTf = new ZConstraint >(&(info[k]), the_ff); + myNPTf = new ZConstraint > >(&(info[k]), the_ff); } else - myNPTf = new NPTf(&(info[k]), the_ff); + myNPTf = new NPTf >(&(info[k]), the_ff); myNPTf->setTargetTemp(globals->getTargetTemp()); @@ -1436,112 +1652,77 @@ void SimSetup::makeIntegrator(void){ else{ sprintf(painCave.errMsg, "SimSetup error: If you use a constant pressure\n" - " ensemble, you must set targetPressure in the BASS file.\n"); + "\tensemble, you must set targetPressure in the BASS file.\n"); painCave.isFatal = 1; simError(); } if (globals->haveTauThermostat()) myNPTf->setTauThermostat(globals->getTauThermostat()); + else{ sprintf(painCave.errMsg, "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauThermostat.\n"); + "\tensemble, you must set tauThermostat.\n"); painCave.isFatal = 1; simError(); } if (globals->haveTauBarostat()) myNPTf->setTauBarostat(globals->getTauBarostat()); - else{ - sprintf(painCave.errMsg, - "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauBarostat.\n"); - painCave.isFatal = 1; - simError(); - } - break; - case NPTim_ENS: - if (globals->haveZconstraints()){ - setupZConstraint(info[k]); - myNPTim = new ZConstraint >(&(info[k]), the_ff); - } - else - myNPTim = new NPTim(&(info[k]), the_ff); - - myNPTim->setTargetTemp(globals->getTargetTemp()); - - if (globals->haveTargetPressure()) - myNPTim->setTargetPressure(globals->getTargetPressure()); else{ sprintf(painCave.errMsg, - "SimSetup error: If you use a constant pressure\n" - " ensemble, you must set targetPressure in the BASS file.\n"); - painCave.isFatal = 1; - simError(); - } - - if (globals->haveTauThermostat()) - myNPTim->setTauThermostat(globals->getTauThermostat()); - else{ - sprintf(painCave.errMsg, "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauThermostat.\n"); + "\tensemble, you must set tauBarostat.\n"); painCave.isFatal = 1; simError(); } - if (globals->haveTauBarostat()) - myNPTim->setTauBarostat(globals->getTauBarostat()); - else{ - sprintf(painCave.errMsg, - "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauBarostat.\n"); - painCave.isFatal = 1; - simError(); - } + info->the_integrator = myNPTf; break; - case NPTfm_ENS: + case NPTxyz_ENS: if (globals->haveZconstraints()){ setupZConstraint(info[k]); - myNPTfm = new ZConstraint >(&(info[k]), the_ff); + myNPTxyz = new ZConstraint > >(&(info[k]), the_ff); } else - myNPTfm = new NPTfm(&(info[k]), the_ff); + myNPTxyz = new NPTxyz >(&(info[k]), the_ff); - myNPTfm->setTargetTemp(globals->getTargetTemp()); + myNPTxyz->setTargetTemp(globals->getTargetTemp()); if (globals->haveTargetPressure()) - myNPTfm->setTargetPressure(globals->getTargetPressure()); + myNPTxyz->setTargetPressure(globals->getTargetPressure()); else{ sprintf(painCave.errMsg, "SimSetup error: If you use a constant pressure\n" - " ensemble, you must set targetPressure in the BASS file.\n"); + "\tensemble, you must set targetPressure in the BASS file.\n"); painCave.isFatal = 1; simError(); - } + } if (globals->haveTauThermostat()) - myNPTfm->setTauThermostat(globals->getTauThermostat()); + myNPTxyz->setTauThermostat(globals->getTauThermostat()); else{ sprintf(painCave.errMsg, "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauThermostat.\n"); + "\tensemble, you must set tauThermostat.\n"); painCave.isFatal = 1; simError(); } if (globals->haveTauBarostat()) - myNPTfm->setTauBarostat(globals->getTauBarostat()); + myNPTxyz->setTauBarostat(globals->getTauBarostat()); else{ sprintf(painCave.errMsg, "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauBarostat.\n"); + "\tensemble, you must set tauBarostat.\n"); painCave.isFatal = 1; simError(); } + + info->the_integrator = myNPTxyz; break; default: @@ -1589,8 +1770,8 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ } else{ sprintf(painCave.errMsg, - "ZConstraint error: If you use an ZConstraint\n" - " , you must set sample time.\n"); + "ZConstraint error: If you use a ZConstraint,\n" + "\tyou must set zconsTime.\n"); painCave.isFatal = 1; simError(); } @@ -1605,8 +1786,9 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ else{ double defaultZConsTol = 0.01; sprintf(painCave.errMsg, - "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n" - " , default value %f is used.\n", + "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n" + "\tOOPSE will use a default value of %f.\n" + "\tTo set the tolerance, use the zconsTol variable.\n", defaultZConsTol); painCave.isFatal = 0; simError(); @@ -1624,8 +1806,9 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ } else{ sprintf(painCave.errMsg, - "ZConstraint Warning: User does not set force Subtraction policy, " - "PolicyByMass is used\n"); + "ZConstraint Warning: No force subtraction policy was set.\n" + "\tOOPSE will use PolicyByMass.\n" + "\tTo set the policy, use the zconsForcePolicy variable.\n"); painCave.isFatal = 0; simError(); zconsForcePolicy->setData("BYMASS"); @@ -1633,6 +1816,33 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ theInfo.addProperty(zconsForcePolicy); + //set zcons gap + DoubleData* zconsGap = new DoubleData(); + zconsGap->setID(ZCONSGAP_ID); + + if (globals->haveZConsGap()){ + zconsGap->setData(globals->getZconsGap()); + theInfo.addProperty(zconsGap); + } + + //set zcons fixtime + DoubleData* zconsFixtime = new DoubleData(); + zconsFixtime->setID(ZCONSFIXTIME_ID); + + if (globals->haveZConsFixTime()){ + zconsFixtime->setData(globals->getZconsFixtime()); + theInfo.addProperty(zconsFixtime); + } + + //set zconsUsingSMD + IntData* zconsUsingSMD = new IntData(); + zconsUsingSMD->setID(ZCONSUSINGSMD_ID); + + if (globals->haveZConsUsingSMD()){ + zconsUsingSMD->setData(globals->getZconsUsingSMD()); + theInfo.addProperty(zconsUsingSMD); + } + //Determine the name of ouput file and add it into SimInfo's property list //Be careful, do not use inFileName, since it is a pointer which //point to a string at master node, and slave nodes do not contain that string @@ -1662,14 +1872,15 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ tempParaItem.zPos = zconStamp[i]->getZpos(); tempParaItem.zconsIndex = zconStamp[i]->getMolIndex(); tempParaItem.kRatio = zconStamp[i]->getKratio(); - + tempParaItem.havingCantVel = zconStamp[i]->haveCantVel(); + tempParaItem.cantVel = zconStamp[i]->getCantVel(); zconsParaData->addItem(tempParaItem); } //check the uniqueness of index if(!zconsParaData->isIndexUnique()){ sprintf(painCave.errMsg, - "ZConstraint Error: molIndex is not unique\n"); + "ZConstraint Error: molIndex is not unique!\n"); painCave.isFatal = 1; simError(); } @@ -1680,3 +1891,73 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ //push data into siminfo, therefore, we can retrieve later theInfo.addProperty(zconsParaData); } + +void SimSetup::makeMinimizer(){ + + OOPSEMinimizer* myOOPSEMinimizer; + MinimizerParameterSet* param; + char minimizerName[100]; + + for (int i = 0; i < nInfo; i++){ + + //prepare parameter set for minimizer + param = new MinimizerParameterSet(); + param->setDefaultParameter(); + + if (globals->haveMinimizer()){ + param->setFTol(globals->getMinFTol()); + } + + if (globals->haveMinGTol()){ + param->setGTol(globals->getMinGTol()); + } + + if (globals->haveMinMaxIter()){ + param->setMaxIteration(globals->getMinMaxIter()); + } + + if (globals->haveMinWriteFrq()){ + param->setMaxIteration(globals->getMinMaxIter()); + } + + if (globals->haveMinWriteFrq()){ + param->setWriteFrq(globals->getMinWriteFrq()); + } + + if (globals->haveMinStepSize()){ + param->setStepSize(globals->getMinStepSize()); + } + + if (globals->haveMinLSMaxIter()){ + param->setLineSearchMaxIteration(globals->getMinLSMaxIter()); + } + + if (globals->haveMinLSTol()){ + param->setLineSearchTol(globals->getMinLSTol()); + } + + strcpy(minimizerName, globals->getMinimizer()); + + if (!strcasecmp(minimizerName, "CG")){ + myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); + } + else if (!strcasecmp(minimizerName, "SD")){ + //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param); + myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param); + } + else{ + sprintf(painCave.errMsg, + "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n"); + painCave.isFatal = 0; + simError(); + + myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); + } + info[i].the_integrator = myOOPSEMinimizer; + + //store the minimizer into simInfo + info[i].the_minimizer = myOOPSEMinimizer; + info[i].has_minimizer = true; + } + +}