--- trunk/OOPSE/libmdtools/SimSetup.cpp 2004/03/16 19:22:56 1091 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2004/04/12 20:32:20 1097 @@ -9,6 +9,7 @@ #include "parse_me.h" #include "Integrator.hpp" #include "simError.h" +#include "RigidBody.hpp" //#include "ConjugateMinimizer.hpp" #include "OOPSEMinimizer.hpp" @@ -165,38 +166,44 @@ void SimSetup::makeMolecules(void){ 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; + 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; bond_pair* theBonds; bend_set* theBends; torsion_set* theTorsions; + set skipList; + + double phi, theta, psi; + //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])); atomOffset = 0; - excludeOffset = 0; + for (i = 0; i < info[k].n_mol; i++){ stampID = info[k].molecules[i].getStampID(); @@ -204,22 +211,34 @@ void SimSetup::makeMolecules(void){ 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.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()); @@ -233,42 +252,15 @@ 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); - - u = sqrt(uSqr); - ux = ux / u; - uy = uy / u; - uz = uz / u; - - dAtom->setSUx(ux); - dAtom->setSUy(uy); - dAtom->setSUz(uz); + 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()); #ifdef IS_MPI @@ -284,28 +276,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; +#ifdef IS_MPI + 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++){ @@ -355,33 +338,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); @@ -390,38 +381,100 @@ void SimSetup::makeMolecules(void){ theTorsions[j].c = currentTorsion->getC() + atomOffset; theTorsions[j].d = currentTorsion->getD() + atomOffset; - exI = theTorsions[j].a; - exJ = theTorsions[j].d; - - // exclude_I must always be the smaller of the pair - if (exI > exJ){ - tempEx = exI; - exI = exJ; - exJ = tempEx; - } + tempI = theTorsions[j].a; + tempJ = theTorsions[j].b; + tempK = theTorsions[j].c; + tempL = theTorsions[j].d; + #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; + 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(); + + 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 = 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); + + } + } + } + + // 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); @@ -429,7 +482,28 @@ void SimSetup::makeMolecules(void){ delete[] theBonds; delete[] theBends; delete[] theTorsions; + } + + // build up the integrableObjects vector: + + for (i = 0; i < info[k].n_atoms; i++) { + +#ifdef IS_MPI + slI = info[k].atoms[i]->getGlobalIndex(); +#else + slI = i; +#endif + + if (skipList.find(slI) == skipList.end()) { + mySD = (StuntDouble *) info[k].atoms[i]; + info[k].integrableObjects.push_back(mySD); + } } + for (i = 0; i < info[k].rigidBodies.size(); i++) { + mySD = (StuntDouble *) info[k].rigidBodies[i]; + info[k].integrableObjects.push_back(mySD); + } + } #ifdef IS_MPI @@ -439,7 +513,15 @@ void SimSetup::makeMolecules(void){ // clean up the forcefield - the_ff->calcRcut(); + if (!globals->haveLJrcut()){ + + the_ff->calcRcut(); + + } else { + + the_ff->setRcut( globals->getLJrcut() ); + } + the_ff->cleanMe(); } @@ -1221,13 +1303,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]; @@ -1249,6 +1333,7 @@ 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; mpiSim = new mpiSimulation(info); @@ -1265,9 +1350,9 @@ void SimSetup::mpiMolDivide(void){ local_bonds = 0; local_bends = 0; local_torsions = 0; + local_rigid = 0; globalAtomIndex = 0; - for (i = 0; i < n_components; i++){ for (j = 0; j < components_nmol[i]; j++){ if (mol2proc[allMol] == worldRank){ @@ -1275,6 +1360,7 @@ 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++){ @@ -1288,6 +1374,7 @@ 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, @@ -1320,9 +1407,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 @@ -1365,32 +1450,14 @@ 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); } @@ -1681,7 +1748,15 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ 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 @@ -1712,7 +1787,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); }