--- trunk/OOPSE/libmdtools/SimSetup.cpp 2004/06/01 18:42:58 1214 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2004/06/11 14:14:10 1261 @@ -11,8 +11,9 @@ #include "simError.h" #include "RigidBody.hpp" #include "OOPSEMinimizer.hpp" -//#include "ConstraintElement.hpp" -//#include "ConstraintPair.hpp" +#include "ConstraintElement.hpp" +#include "ConstraintPair.hpp" +#include "ConstraintManager.hpp" #ifdef IS_MPI #include "mpiBASS.h" @@ -157,6 +158,10 @@ void SimSetup::createSim(void){ initFortran(); + //creat constraint manager + for(int i = 0; i < nInfo; i++) + info[i].consMan = new ConstraintManager(&info[i]); + if (globals->haveMinimizer()) // make minimizer makeMinimizer(); @@ -199,12 +204,13 @@ void SimSetup::makeMolecules(void){ 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; + 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(); @@ -516,9 +522,9 @@ void SimSetup::makeMolecules(void){ myCutoffGroup = new CutoffGroup(); #ifdef IS_MPI - myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]); + myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); #else - myCutoffGroup->setGlobalIndex(j + groupOffset); + myCutoffGroup->setGlobalIndex(groupOffset); #endif for (int cg = 0; cg < nMembers; cg++) { @@ -531,16 +537,12 @@ void SimSetup::makeMolecules(void){ #ifdef IS_MPI globalID = info[k].atoms[tempI]->getGlobalIndex(); - info[k].globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset]; + info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; #else globalID = info[k].atoms[tempI]->getIndex(); - info[k].globalGroupMembership[globalID] = j + groupOffset; -#endif - - - - myCutoffGroup->addAtom(info[k].atoms[tempI]); - + info[k].globalGroupMembership[globalID] = groupOffset; +#endif + myCutoffGroup->addAtom(info[k].atoms[tempI]); cutoffAtomSet.insert(tempI); } @@ -549,27 +551,28 @@ void SimSetup::makeMolecules(void){ }//end for (j = 0; j < molInfo.nCutoffGroups; j++) - //creat a cutoff group for every atom in current molecule which does not belong to cutoffgroup defined at mdl file + // 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[j + groupOffset]); + myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); globalID = info[k].atoms[atomOffset + j]->getGlobalIndex(); - info[k].globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset]; + info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; #else - myCutoffGroup->setGlobalIndex(j + groupOffset); + myCutoffGroup->setGlobalIndex(groupOffset); globalID = info[k].atoms[atomOffset + j]->getIndex(); - info[k].globalGroupMembership[globalID] = j+groupOffset; + info[k].globalGroupMembership[globalID] = groupOffset; #endif molInfo.myCutoffGroups.push_back(myCutoffGroup); groupOffset++; - } - + } } // After this is all set up, scan through the atoms to @@ -602,37 +605,50 @@ void SimSetup::makeMolecules(void){ 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); - /* - //creat ConstraintPair. - molInfo.myConstraintPair.clear(); + molInfo.myConstraintPairs.clear(); for (j = 0; j < molInfo.nBonds; j++){ - //if both atoms are in the same rigid body, just skip it - currentBond = comp_stamps[stampID]->getBond(j); - if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){ + //if bond is constrained bond, add it into constraint pair + if(molInfo.myBonds[j]->is_constrained()){ - 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]); + //if both atoms are in the same rigid body, just skip it + currentBond = comp_stamps[stampID]->getBond(j); + + if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){ - 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]); + 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]); - consPair = new DistanceConstraintPair(consElement1, consElement2); - molInfo.myConstraintPairs.push_back(consPair); - } + 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 HingeConstraintPair + //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++){ @@ -649,14 +665,6 @@ void SimSetup::makeMolecules(void){ } } -*/ - // 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); @@ -908,8 +916,12 @@ void SimSetup::gatherInfo(void){ painCave.isFatal = 1; simError(); } - - // get the ensemble + if (globals->haveForceFieldVariant()) { + strcpy(forcefield_variant, globals->getForceFieldVariant()); + has_forcefield_variant = 1; + } + + // get the ensemble strcpy(ensemble, globals->getEnsemble()); @@ -1506,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: @@ -1514,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: @@ -1528,6 +1543,7 @@ void SimSetup::createFF(void){ simError(); } + #ifdef IS_MPI strcpy(checkPointMsg, "ForceField creation successful"); MPIcheckPoint();