--- trunk/OOPSE/libmdtools/SimSetup.cpp 2004/02/06 21:37:59 1035 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2004/06/11 14:14:10 1261 @@ -9,7 +9,11 @@ #include "parse_me.h" #include "Integrator.hpp" #include "simError.h" -#include "ConjugateMinimizer.hpp" +#include "RigidBody.hpp" +#include "OOPSEMinimizer.hpp" +#include "ConstraintElement.hpp" +#include "ConstraintPair.hpp" +#include "ConstraintManager.hpp" #ifdef IS_MPI #include "mpiBASS.h" @@ -28,7 +32,7 @@ #define FF_DUFF 0 #define FF_LJ 1 #define FF_EAM 2 -#define FF_H2O 3 +#define FF_H2O 3 using namespace std; @@ -145,13 +149,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(); @@ -160,65 +157,122 @@ void SimSetup::createSim(void){ // initialize the Fortran initFortran(); + + //creat constraint manager + for(int i = 0; i < nInfo; i++) + info[i].consMan = new ConstraintManager(&info[i]); + + if (globals->haveMinimizer()) + // make minimizer + makeMinimizer(); + else + // make the integrator + makeIntegrator(); + } void SimSetup::makeMolecules(void){ - int k; - 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, globalID; + int stampID, atomOffset, rbOffset, groupOffset; 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; + 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; torsion_set* theTorsions; + set skipList; + + double phi, theta, psi; + char* molName; + char rbName[100]; + + ConstraintPair* consPair; //constraint pair + ConstraintElement* consElement1; //first element of constraint pair + ConstraintElement* consElement2; //second element of constraint pair + int whichRigidBody; + int consAtomIndex; //index of constraint atom in rigid body's atom array + vector > jointAtoms; + double bondLength2; //init the forceField paramters the_ff->readParams(); - // init the atoms - double phi, theta, psi; - double sux, suy, suz; - double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz; - double ux, uy, uz, u, uSqr; + int nMembers, nNew, rb1, rb2; for (k = 0; k < nInfo; k++){ the_ff->setSimInfo(&(info[k])); +#ifdef IS_MPI + info[k].globalGroupMembership = new int[mpiSim->getNAtomsGlobal()]; + for (i = 0; i < mpiSim->getNAtomsGlobal(); i++) + info[k].globalGroupMembership[i] = 0; +#else + info[k].globalGroupMembership = new int[info[k].n_atoms]; + for (i = 0; i < info[k].n_atoms; i++) + info[k].globalGroupMembership[i] = 0; +#endif + atomOffset = 0; - excludeOffset = 0; + groupOffset = 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(); + 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()); @@ -232,48 +286,19 @@ void SimSetup::makeMolecules(void){ phi = currentAtom->getEulerPhi() * M_PI / 180.0; theta = currentAtom->getEulerTheta() * M_PI / 180.0; psi = currentAtom->getEulerPsi()* M_PI / 180.0; - - Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi)); - Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi)); - Axz = sin(theta) * sin(psi); - - Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi)); - Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi)); - Ayz = sin(theta) * cos(psi); - - Azx = sin(phi) * sin(theta); - Azy = -cos(phi) * sin(theta); - Azz = cos(theta); - sux = 0.0; - suy = 0.0; - suz = 1.0; - - ux = (Axx * sux) + (Ayx * suy) + (Azx * suz); - uy = (Axy * sux) + (Ayy * suy) + (Azy * suz); - uz = (Axz * sux) + (Ayz * suy) + (Azz * suz); - - uSqr = (ux * ux) + (uy * uy) + (uz * uz); + dAtom->setUnitFrameFromEuler(phi, theta, psi); + + } + else{ - u = sqrt(uSqr); - ux = ux / u; - uy = uy / u; - uz = uz / u; + molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration()); - dAtom->setSUx(ux); - dAtom->setSUy(uy); - dAtom->setSUz(uz); } - else{ - molInfo.myAtoms[j] = new GeneralAtom((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 } @@ -283,28 +308,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++){ @@ -354,33 +370,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); @@ -389,31 +413,201 @@ 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; + 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[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); + 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); } - excludeOffset += molInfo.nTorsions; + + molInfo.myRigidBodies.clear(); + + for (j = 0; j < molInfo.nRigidBodies; j++){ - // send the arrays off to the forceField for init. + 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 + + cutoffAtomSet.clear(); + molInfo.myCutoffGroups.clear(); + + for (j = 0; j < nCutoffGroups; j++){ + + currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j); + nMembers = currentCutoffGroup->getNMembers(); + + myCutoffGroup = new CutoffGroup(); + +#ifdef IS_MPI + myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); +#else + myCutoffGroup->setGlobalIndex(groupOffset); +#endif + + 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; + +#ifdef IS_MPI + globalID = info[k].atoms[tempI]->getGlobalIndex(); + info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; +#else + globalID = info[k].atoms[tempI]->getIndex(); + info[k].globalGroupMembership[globalID] = groupOffset; +#endif + myCutoffGroup->addAtom(info[k].atoms[tempI]); + cutoffAtomSet.insert(tempI); + } + + molInfo.myCutoffGroups.push_back(myCutoffGroup); + groupOffset++; + + }//end for (j = 0; j < molInfo.nCutoffGroups; j++) + + + // create 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]); + +#ifdef IS_MPI + myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); + globalID = info[k].atoms[atomOffset + j]->getGlobalIndex(); + info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; +#else + myCutoffGroup->setGlobalIndex(groupOffset); + globalID = info[k].atoms[atomOffset + j]->getIndex(); + info[k].globalGroupMembership[globalID] = groupOffset; +#endif + molInfo.myCutoffGroups.push_back(myCutoffGroup); + groupOffset++; + } + } + + // 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); @@ -421,14 +615,91 @@ void SimSetup::makeMolecules(void){ theTorsions); - info[k].molecules[i].initialize(molInfo); + //creat ConstraintPair. + molInfo.myConstraintPairs.clear(); + + for (j = 0; j < molInfo.nBonds; j++){ + //if bond is constrained bond, add it into constraint pair + if(molInfo.myBonds[j]->is_constrained()){ + //if both atoms are in the same rigid body, just skip it + currentBond = comp_stamps[stampID]->getBond(j); + + if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){ + + tempI = currentBond->getA() + atomOffset; + if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex)) + consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); + else + consElement1 = new ConstraintAtom(info[k].atoms[tempI]); + + tempJ = currentBond->getB() + atomOffset; + if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex)) + consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); + else + consElement2 = new ConstraintAtom(info[k].atoms[tempJ]); + + bondLength2 = molInfo.myBonds[j]->get_constraint()->get_dsqr(); + consPair = new DistanceConstraintPair(consElement1, consElement2, bondLength2); + + molInfo.myConstraintPairs.push_back(consPair); + } + }//end if(molInfo.myBonds[j]->is_constrained()) + } + + //loop over rigid bodies, if two rigid bodies share same joint, creat a JointConstraintPair + for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){ + for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){ + + jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2); + + for(size_t m = 0; m < jointAtoms.size(); m++){ + consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first); + consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second); + + consPair = new JointConstraintPair(consElement1, consElement2); + molInfo.myConstraintPairs.push_back(consPair); + } + + } + } + + + info[k].molecules[i].initialize(molInfo); + + atomOffset += molInfo.nAtoms; delete[] theBonds; delete[] theBends; delete[] theTorsions; } + + + +#ifdef IS_MPI + // Since the globalGroupMembership has been zero filled and we've only + // poked values into the atoms we know, we can do an Allreduce + // to get the full globalGroupMembership array (We think). + // This would be prettier if we could use MPI_IN_PLACE like the MPI-2 + // docs said we could. + + int* ggMjunk = new int[mpiSim->getNAtomsGlobal()]; + + MPI_Allreduce(info[k].globalGroupMembership, + ggMjunk, + mpiSim->getNAtomsGlobal(), + MPI_INT, MPI_SUM, MPI_COMM_WORLD); + + for (i = 0; i < mpiSim->getNAtomsGlobal(); i++) + info[k].globalGroupMembership[i] = ggMjunk[i]; + + delete[] ggMjunk; + +#endif + + + } #ifdef IS_MPI @@ -436,10 +707,6 @@ void SimSetup::makeMolecules(void){ MPIcheckPoint(); #endif // is_mpi - // clean up the forcefield - - the_ff->calcRcut(); - the_ff->cleanMe(); } void SimSetup::initFromBass(void){ @@ -649,9 +916,13 @@ void SimSetup::gatherInfo(void){ painCave.isFatal = 1; simError(); } + if (globals->haveForceFieldVariant()) { + strcpy(forcefield_variant, globals->getForceFieldVariant()); + has_forcefield_variant = 1; + } + + // get the ensemble - // get the ensemble - strcpy(ensemble, globals->getEnsemble()); if (!strcasecmp(ensemble, "NVE")){ @@ -726,7 +997,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" @@ -736,7 +1007,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" @@ -772,12 +1043,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()){ @@ -786,6 +1055,8 @@ void SimSetup::gatherInfo(void){ if (globals->haveThermalTime()){ info[i].thermalTime = globals->getThermalTime(); + } else { + info[i].thermalTime = globals->getRunTime(); } info[i].resetIntegrator = 0; @@ -803,7 +1074,67 @@ void SimSetup::gatherInfo(void){ info[i].useInitXSstate = globals->getUseInitXSstate(); info[i].orthoTolerance = globals->getOrthoBoxTolerance(); - + + // check for thermodynamic integration + if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) { + if (globals->haveThermIntLambda() && globals->haveThermIntK()) { + info[i].useSolidThermInt = globals->getUseSolidThermInt(); + info[i].thermIntLambda = globals->getThermIntLambda(); + info[i].thermIntK = globals->getThermIntK(); + + Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK); + info[i].restraint = myRestraint; + } + else { + sprintf(painCave.errMsg, + "SimSetup Error:\n" + "\tKeyword useSolidThermInt was set to 'true' but\n" + "\tthermodynamicIntegrationLambda (and/or\n" + "\tthermodynamicIntegrationK) was not specified.\n" + "\tPlease provide a lambda value and k value in your .bass file.\n"); + painCave.isFatal = 1; + simError(); + } + } + else if(globals->getUseLiquidThermInt()) { + if (globals->getUseSolidThermInt()) { + sprintf( painCave.errMsg, + "SimSetup Warning: It appears that you have both solid and\n" + "\tliquid thermodynamic integration activated in your .bass\n" + "\tfile. To avoid confusion, specify only one technique in\n" + "\tyour .bass file. Liquid-state thermodynamic integration\n" + "\twill be assumed for the current simulation. If this is not\n" + "\twhat you desire, set useSolidThermInt to 'true' and\n" + "\tuseLiquidThermInt to 'false' in your .bass file.\n"); + painCave.isFatal = 0; + simError(); + } + if (globals->haveThermIntLambda() && globals->haveThermIntK()) { + info[i].useLiquidThermInt = globals->getUseLiquidThermInt(); + info[i].thermIntLambda = globals->getThermIntLambda(); + info[i].thermIntK = globals->getThermIntK(); + } + else { + sprintf(painCave.errMsg, + "SimSetup Error:\n" + "\tKeyword useLiquidThermInt was set to 'true' but\n" + "\tthermodynamicIntegrationLambda (and/or\n" + "\tthermodynamicIntegrationK) was not specified.\n" + "\tPlease provide a lambda value and k value in your .bass file.\n"); + painCave.isFatal = 1; + simError(); + } + } + else if(globals->haveThermIntLambda() || globals->haveThermIntK()){ + sprintf(painCave.errMsg, + "SimSetup Warning: If you want to use Thermodynamic\n" + "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n" + "\t'true' in your .bass file. These keywords are set to\n" + "\t'false' by default, so your lambda and/or k values are\n" + "\tbeing ignored.\n"); + painCave.isFatal = 0; + simError(); + } } //setup seed for random number generator @@ -856,6 +1187,7 @@ void SimSetup::finalInfoCheck(void){ void SimSetup::finalInfoCheck(void){ int index; int usesDipoles; + int usesCharges; int i; for (i = 0; i < nInfo; i++){ @@ -867,45 +1199,72 @@ 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: 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, @@ -918,35 +1277,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); } } } @@ -954,6 +1314,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){ @@ -1081,7 +1444,29 @@ void SimSetup::makeOutNames(void){ } else{ strcat(info[k].statusName, ".stat"); + } + } + + strcpy(info[k].rawPotName, inFileName); + nameLength = strlen(info[k].rawPotName); + endTest = &(info[k].rawPotName[nameLength - 5]); + if (!strcmp(endTest, ".bass")){ + strcpy(endTest, ".raw"); + } + else if (!strcmp(endTest, ".BASS")){ + strcpy(endTest, ".raw"); + } + else{ + endTest = &(info[k].rawPotName[nameLength - 4]); + if (!strcmp(endTest, ".bss")){ + strcpy(endTest, ".raw"); } + else if (!strcmp(endTest, ".mdl")){ + strcpy(endTest, ".raw"); + } + else{ + strcat(info[k].rawPotName, ".raw"); + } } #ifdef IS_MPI @@ -1133,7 +1518,7 @@ void SimSetup::createFF(void){ void SimSetup::createFF(void){ switch (ffCase){ case FF_DUFF: - the_ff = new DUFF(); + the_ff = new DUFF(); break; case FF_LJ: @@ -1141,7 +1526,10 @@ void SimSetup::createFF(void){ break; case FF_EAM: - the_ff = new EAM_FF(); + if (has_forcefield_variant) + the_ff = new EAM_FF(forcefield_variant); + else + the_ff = new EAM_FF(); break; case FF_H2O: @@ -1155,6 +1543,7 @@ void SimSetup::createFF(void){ simError(); } + #ifdef IS_MPI strcpy(checkPointMsg, "ForceField creation successful"); MPIcheckPoint(); @@ -1168,7 +1557,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 @@ -1203,7 +1595,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"); @@ -1212,21 +1610,35 @@ void SimSetup::calcSysValues(void){ } void SimSetup::calcSysValues(void){ - int i; + int i, j; + int ncutgroups, atomsingroups, ngroupsinstamp; int* molMembershipArray; + CutoffGroupStamp* cg; tot_atoms = 0; tot_bonds = 0; tot_bends = 0; tot_torsions = 0; + tot_rigid = 0; + tot_groups = 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(); + ncutgroups = comp_stamps[i]->getNCutoffGroups(); + atomsingroups = 0; + for (j=0; j < ncutgroups; j++) { + cg = comp_stamps[i]->getCutoffGroup(j); + atomsingroups += cg->getNMembers(); + } + ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups; + tot_groups += components_nmol[i] * ngroupsinstamp; + } + tot_SRI = tot_bonds + tot_bends + tot_torsions; molMembershipArray = new int[tot_atoms]; @@ -1237,7 +1649,7 @@ void SimSetup::calcSysValues(void){ info[i].n_torsions = tot_torsions; info[i].n_SRI = tot_SRI; info[i].n_mol = tot_nmol; - + info[i].ngroup = tot_groups; info[i].molMembershipArray = molMembershipArray; } } @@ -1248,10 +1660,17 @@ 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, local_groups; + vector globalMolIndex; + int ncutgroups, atomsingroups, ngroupsinstamp; + CutoffGroupStamp* cg; mpiSim = new mpiSimulation(info); - globalIndex = mpiSim->divideLabor(); + mpiSim->divideLabor(); + globalAtomIndex = mpiSim->getGlobalAtomIndex(); + globalGroupIndex = mpiSim->getGlobalGroupIndex(); + //globalMolIndex = mpiSim->getGlobalMolIndex(); // set up the local variables @@ -1264,9 +1683,10 @@ void SimSetup::mpiMolDivide(void){ local_bonds = 0; local_bends = 0; local_torsions = 0; - globalAtomIndex = 0; + local_rigid = 0; + local_groups = 0; + globalAtomCounter = 0; - for (i = 0; i < n_components; i++){ for (j = 0; j < components_nmol[i]; j++){ if (mol2proc[allMol] == worldRank){ @@ -1274,11 +1694,23 @@ 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(); + + ncutgroups = comp_stamps[i]->getNCutoffGroups(); + atomsingroups = 0; + for (k=0; k < ncutgroups; k++) { + cg = comp_stamps[i]->getCutoffGroup(k); + atomsingroups += cg->getNMembers(); + } + ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + + ncutgroups; + local_groups += ngroupsinstamp; + localMol++; } for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ - info[0].molMembershipArray[globalAtomIndex] = allMol; - globalAtomIndex++; + info[0].molMembershipArray[globalAtomCounter] = allMol; + globalAtomCounter++; } allMol++; @@ -1286,8 +1718,8 @@ void SimSetup::mpiMolDivide(void){ } local_SRI = local_bonds + local_bends + local_torsions; - info[0].n_atoms = mpiSim->getMyNlocal(); - + info[0].n_atoms = mpiSim->getNAtomsLocal(); + if (local_atoms != info[0].n_atoms){ sprintf(painCave.errMsg, "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" @@ -1297,6 +1729,16 @@ void SimSetup::mpiMolDivide(void){ simError(); } + info[0].ngroup = mpiSim->getNGroupsLocal(); + if (local_groups != info[0].ngroup){ + sprintf(painCave.errMsg, + "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n" + "\tlocalGroups (%d) are not equal.\n", + info[0].ngroup, local_groups); + painCave.isFatal = 1; + simError(); + } + info[0].n_bonds = local_bonds; info[0].n_bends = local_bends; info[0].n_torsions = local_torsions; @@ -1319,9 +1761,7 @@ void SimSetup::makeSysArrays(void){ Atom** the_atoms; Molecule* the_molecules; - Exclude** the_excludes; - for (l = 0; l < nInfo; l++){ // create the atom and short range interaction arrays @@ -1335,7 +1775,7 @@ void SimSetup::makeSysArrays(void){ molIndex = 0; - for (i = 0; i < mpiSim->getTotNmol(); i++){ + for (i = 0; i < mpiSim->getNMolGlobal(); i++){ if (mol2proc[i] == worldRank){ the_molecules[molIndex].setStampID(molCompType[i]); the_molecules[molIndex].setMyIndex(molIndex); @@ -1347,15 +1787,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++; } @@ -1364,33 +1804,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); } } @@ -1663,6 +2085,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 @@ -1692,7 +2141,8 @@ 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); } @@ -1713,33 +2163,12 @@ void SimSetup::makeMinimizer(){ void SimSetup::makeMinimizer(){ - OOPSEMinimizerBase* myOOPSEMinimizerBase; - ObjFunctor1 * objFunc; - OutputFunctor* outputFunc; - ConcreteNLModel1* nlp; + OOPSEMinimizer* myOOPSEMinimizer; MinimizerParameterSet* param; - ConjugateMinimizerBase* minimizer; - int dim; + char minimizerName[100]; for (int i = 0; i < nInfo; i++){ - //creat - myOOPSEMinimizerBase = new OOPSEMinimizerBase(&(info[i]), the_ff); - - info[i].the_integrator = myOOPSEMinimizerBase; - //creat the object functor; - objFunc = (ObjFunctor1*) new ClassMemObjFunctor1 - (myOOPSEMinimizerBase, &OOPSEMinimizerBase::calcGradient); - - //creat output functor; - outputFunc = new ClassMemOutputFunctor - (myOOPSEMinimizerBase, &OOPSEMinimizerBase::output); - - //creat nonlinear model - dim = myOOPSEMinimizerBase->getDim(); - nlp = new ConcreteNLModel1(dim, objFunc); - - nlp->setX(myOOPSEMinimizerBase->getCoor()); - + //prepare parameter set for minimizer param = new MinimizerParameterSet(); param->setDefaultParameter(); @@ -1764,8 +2193,8 @@ void SimSetup::makeMinimizer(){ param->setWriteFrq(globals->getMinWriteFrq()); } - if (globals->haveMinResetFrq()){ - param->setResetFrq(globals->getMinResetFrq()); + if (globals->haveMinStepSize()){ + param->setStepSize(globals->getMinStepSize()); } if (globals->haveMinLSMaxIter()){ @@ -1775,14 +2204,28 @@ void SimSetup::makeMinimizer(){ if (globals->haveMinLSTol()){ param->setLineSearchTol(globals->getMinLSTol()); } - - //creat the minimizer - minimizer = new PRCGMinimizer(nlp, param); - minimizer->setLineSearchStrategy(nlp, GoldenSection); - minimizer->setOutputFunctor(outputFunc); + 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 = minimizer; + info[i].the_minimizer = myOOPSEMinimizer; info[i].has_minimizer = true; }