--- trunk/OOPSE/libmdtools/SimSetup.cpp 2004/05/12 16:38:45 1167 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2004/08/23 15:11:36 1452 @@ -10,8 +10,10 @@ #include "Integrator.hpp" #include "simError.h" #include "RigidBody.hpp" -//#include "ConjugateMinimizer.hpp" #include "OOPSEMinimizer.hpp" +#include "ConstraintElement.hpp" +#include "ConstraintPair.hpp" +#include "ConstraintManager.hpp" #ifdef IS_MPI #include "mpiBASS.h" @@ -156,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(); @@ -170,8 +176,8 @@ void SimSetup::makeMolecules(void){ int i, j, k; int exI, exJ, exK, exL, slI, slJ; int tempI, tempJ, tempK, tempL; - int molI; - int stampID, atomOffset, rbOffset; + int molI, globalID; + int stampID, atomOffset, rbOffset, groupOffset; molInit molInfo; DirectionalAtom* dAtom; RigidBody* myRB; @@ -198,6 +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; + double bondLength2; //init the forceField paramters the_ff->readParams(); @@ -209,7 +222,18 @@ void SimSetup::makeMolecules(void){ 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; + groupOffset = 0; for (i = 0; i < info[k].n_mol; i++){ stampID = info[k].molecules[i].getStampID(); @@ -226,17 +250,17 @@ void SimSetup::makeMolecules(void){ molInfo.myAtoms = &(info[k].atoms[atomOffset]); if (molInfo.nBonds > 0) - molInfo.myBonds = new (Bond *) [molInfo.nBonds]; + molInfo.myBonds = new Bond*[molInfo.nBonds]; else molInfo.myBonds = NULL; if (molInfo.nBends > 0) - molInfo.myBends = new (Bend *) [molInfo.nBends]; + molInfo.myBends = new Bend*[molInfo.nBends]; else molInfo.myBends = NULL; if (molInfo.nTorsions > 0) - molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions]; + molInfo.myTorsions = new Torsion *[molInfo.nTorsions]; else molInfo.myTorsions = NULL; @@ -274,9 +298,7 @@ void SimSetup::makeMolecules(void){ molInfo.myAtoms[j]->setType(currentAtom->getType()); #ifdef IS_MPI - molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]); - #endif // is_mpi } @@ -499,6 +521,12 @@ void SimSetup::makeMolecules(void){ 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 @@ -506,30 +534,47 @@ void SimSetup::makeMolecules(void){ // tempI is atom numbering on local processor tempI = molI + atomOffset; - - myCutoffGroup->addAtom(info[k].atoms[tempI]); - + +#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); - }//end for (j = 0; j < molInfo.nCutoffGroups; j++) + groupOffset++; - //creat a cutoff group for every atom in current molecule which does not belong to cutoffgroup defined at mdl file - + }//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]); - molInfo.myCutoffGroups.push_back(myCutoffGroup); - } +#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: @@ -560,8 +605,7 @@ 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); @@ -570,14 +614,92 @@ void SimSetup::makeMolecules(void){ the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions, 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 @@ -794,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")){ @@ -948,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 @@ -1258,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 @@ -1310,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: @@ -1318,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: @@ -1332,6 +1543,7 @@ void SimSetup::createFF(void){ simError(); } + #ifdef IS_MPI strcpy(checkPointMsg, "ForceField creation successful"); MPIcheckPoint(); @@ -1398,21 +1610,33 @@ 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; @@ -1425,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; } } @@ -1436,13 +1660,16 @@ 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; + int local_rigid, local_groups; vector globalMolIndex; + int ncutgroups, atomsingroups, ngroupsinstamp; + CutoffGroupStamp* cg; mpiSim = new mpiSimulation(info); mpiSim->divideLabor(); globalAtomIndex = mpiSim->getGlobalAtomIndex(); + globalGroupIndex = mpiSim->getGlobalGroupIndex(); //globalMolIndex = mpiSim->getGlobalMolIndex(); // set up the local variables @@ -1457,6 +1684,7 @@ void SimSetup::mpiMolDivide(void){ local_bends = 0; local_torsions = 0; local_rigid = 0; + local_groups = 0; globalAtomCounter = 0; for (i = 0; i < n_components; i++){ @@ -1467,6 +1695,17 @@ void SimSetup::mpiMolDivide(void){ 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++){ @@ -1479,9 +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" @@ -1491,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; @@ -1527,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); @@ -1572,11 +1820,11 @@ void SimSetup::makeIntegrator(void){ void SimSetup::makeIntegrator(void){ int k; - NVE* myNVE = NULL; - NVT* myNVT = NULL; - NPTi >* myNPTi = NULL; - NPTf >* myNPTf = NULL; - NPTxyz >* myNPTxyz = NULL; + NVE >* myNVE = NULL; + NVT >* myNVT = NULL; + NPTi > >* myNPTi = NULL; + NPTf > >* myNPTf = NULL; + NPTxyz > >* myNPTxyz = NULL; for (k = 0; k < nInfo; k++){ switch (ensembleCase){ @@ -1586,7 +1834,15 @@ void SimSetup::makeIntegrator(void){ myNVE = new ZConstraint >(&(info[k]), the_ff); } else{ - myNVE = new NVE(&(info[k]), the_ff); + if (globals->haveQuaternion()){ + if (globals->getUseQuaternion()) + info->the_integrator = new NVE >(&(info[k]), the_ff); + } + else + info->the_integrator = new NVE(&(info[k]), the_ff); + break; + + //myNVE = new NVE(&(info[k]), the_ff); } info->the_integrator = myNVE;