| 10 |
|
#include "Integrator.hpp" |
| 11 |
|
#include "simError.h" |
| 12 |
|
#include "RigidBody.hpp" |
| 13 |
– |
//#include "ConjugateMinimizer.hpp" |
| 13 |
|
#include "OOPSEMinimizer.hpp" |
| 14 |
+ |
//#include "ConstraintElement.hpp" |
| 15 |
+ |
//#include "ConstraintPair.hpp" |
| 16 |
|
|
| 17 |
|
#ifdef IS_MPI |
| 18 |
|
#include "mpiBASS.h" |
| 199 |
|
char* molName; |
| 200 |
|
char rbName[100]; |
| 201 |
|
|
| 202 |
+ |
//ConstraintPair* consPair; //constraint pair |
| 203 |
+ |
//ConstraintElement* consElement1; //first element of constraint pair |
| 204 |
+ |
//ConstraintElement* consElement2; //second element of constraint pair |
| 205 |
+ |
//int whichRigidBody; |
| 206 |
+ |
//int consAtomIndex; //index of constraint atom in rigid body's atom array |
| 207 |
+ |
//vector<pair<int, int> > jointAtoms; |
| 208 |
|
//init the forceField paramters |
| 209 |
|
|
| 210 |
|
the_ff->readParams(); |
| 217 |
|
the_ff->setSimInfo(&(info[k])); |
| 218 |
|
|
| 219 |
|
atomOffset = 0; |
| 220 |
+ |
groupOffset = 0; |
| 221 |
|
|
| 222 |
|
for (i = 0; i < info[k].n_mol; i++){ |
| 223 |
|
stampID = info[k].molecules[i].getStampID(); |
| 506 |
|
nMembers = currentCutoffGroup->getNMembers(); |
| 507 |
|
|
| 508 |
|
myCutoffGroup = new CutoffGroup(); |
| 509 |
+ |
myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]); |
| 510 |
|
|
| 511 |
|
for (int cg = 0; cg < nMembers; cg++) { |
| 512 |
|
|
| 515 |
|
|
| 516 |
|
// tempI is atom numbering on local processor |
| 517 |
|
tempI = molI + atomOffset; |
| 518 |
< |
|
| 518 |
> |
|
| 519 |
> |
#ifdef IS_MPI |
| 520 |
> |
globalID = info[k].atoms[tempI]->getGlobalIndex() |
| 521 |
> |
#else |
| 522 |
> |
globalID = info[k].atoms[tempI]->getIndex(); |
| 523 |
> |
#endif |
| 524 |
> |
|
| 525 |
> |
globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset]; |
| 526 |
> |
|
| 527 |
|
myCutoffGroup->addAtom(info[k].atoms[tempI]); |
| 528 |
|
|
| 529 |
|
cutoffAtomSet.insert(tempI); |
| 530 |
|
} |
| 531 |
< |
|
| 531 |
> |
|
| 532 |
|
molInfo.myCutoffGroups.push_back(myCutoffGroup); |
| 533 |
+ |
groupOffset++; |
| 534 |
+ |
|
| 535 |
|
}//end for (j = 0; j < molInfo.nCutoffGroups; j++) |
| 536 |
|
|
| 537 |
|
//creat a cutoff group for every atom in current molecule which does not belong to cutoffgroup defined at mdl file |
| 541 |
|
if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){ |
| 542 |
|
myCutoffGroup = new CutoffGroup(); |
| 543 |
|
myCutoffGroup->addAtom(molInfo.myAtoms[j]); |
| 544 |
+ |
myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]); |
| 545 |
+ |
#ifdef IS_MPI |
| 546 |
+ |
globalID = info[k].atoms[atomOffset + j]->getGlobalIndex() |
| 547 |
+ |
#else |
| 548 |
+ |
globalID = info[k].atoms[atomOffset + j]->getIndex(); |
| 549 |
+ |
#endif |
| 550 |
+ |
globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset]; |
| 551 |
|
molInfo.myCutoffGroups.push_back(myCutoffGroup); |
| 552 |
+ |
groupOffset++; |
| 553 |
|
} |
| 554 |
|
|
| 555 |
|
} |
| 556 |
|
|
| 530 |
– |
|
| 531 |
– |
|
| 532 |
– |
|
| 557 |
|
// After this is all set up, scan through the atoms to |
| 558 |
|
// see if they can be added to the integrableObjects: |
| 559 |
|
|
| 584 |
|
info[k].integrableObjects.push_back(mySD); |
| 585 |
|
molInfo.myIntegrableObjects.push_back(mySD); |
| 586 |
|
} |
| 587 |
< |
|
| 587 |
> |
|
| 588 |
> |
|
| 589 |
> |
/* |
| 590 |
> |
|
| 591 |
> |
//creat ConstraintPair. |
| 592 |
> |
molInfo.myConstraintPair.clear(); |
| 593 |
|
|
| 594 |
+ |
for (j = 0; j < molInfo.nBonds; j++){ |
| 595 |
+ |
|
| 596 |
+ |
//if both atoms are in the same rigid body, just skip it |
| 597 |
+ |
currentBond = comp_stamps[stampID]->getBond(j); |
| 598 |
+ |
if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){ |
| 599 |
+ |
|
| 600 |
+ |
tempI = currentBond->getA() + atomOffset; |
| 601 |
+ |
if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex)) |
| 602 |
+ |
consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
| 603 |
+ |
else |
| 604 |
+ |
consElement1 = new ConstraintAtom(info[k].atoms[tempI]); |
| 605 |
+ |
|
| 606 |
+ |
tempJ = currentBond->getB() + atomOffset; |
| 607 |
+ |
if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex)) |
| 608 |
+ |
consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
| 609 |
+ |
else |
| 610 |
+ |
consElement2 = new ConstraintAtom(info[k].atoms[tempJ]); |
| 611 |
+ |
|
| 612 |
+ |
consPair = new DistanceConstraintPair(consElement1, consElement2); |
| 613 |
+ |
molInfo.myConstraintPairs.push_back(consPair); |
| 614 |
+ |
} |
| 615 |
+ |
} |
| 616 |
+ |
|
| 617 |
+ |
//loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair |
| 618 |
+ |
for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){ |
| 619 |
+ |
for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){ |
| 620 |
+ |
|
| 621 |
+ |
jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2); |
| 622 |
+ |
|
| 623 |
+ |
for(size_t m = 0; m < jointAtoms.size(); m++){ |
| 624 |
+ |
consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first); |
| 625 |
+ |
consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second); |
| 626 |
+ |
|
| 627 |
+ |
consPair = new JointConstraintPair(consElement1, consElement2); |
| 628 |
+ |
molInfo.myConstraintPairs.push_back(consPair); |
| 629 |
+ |
} |
| 630 |
+ |
|
| 631 |
+ |
} |
| 632 |
+ |
} |
| 633 |
+ |
|
| 634 |
+ |
*/ |
| 635 |
|
// send the arrays off to the forceField for init. |
| 636 |
|
|
| 637 |
|
the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); |
