| 11 |
|
#include "simError.h" |
| 12 |
|
#include "RigidBody.hpp" |
| 13 |
|
#include "OOPSEMinimizer.hpp" |
| 14 |
< |
//#include "ConstraintElement.hpp" |
| 15 |
< |
//#include "ConstraintPair.hpp" |
| 14 |
> |
#include "ConstraintElement.hpp" |
| 15 |
> |
#include "ConstraintPair.hpp" |
| 16 |
> |
#include "ConstraintManager.hpp" |
| 17 |
|
|
| 18 |
|
#ifdef IS_MPI |
| 19 |
|
#include "mpiBASS.h" |
| 158 |
|
|
| 159 |
|
initFortran(); |
| 160 |
|
|
| 161 |
+ |
//creat constraint manager |
| 162 |
+ |
for(int i = 0; i < nInfo; i++) |
| 163 |
+ |
info[i].consMan = new ConstraintManager(&info[i]); |
| 164 |
+ |
|
| 165 |
|
if (globals->haveMinimizer()) |
| 166 |
|
// make minimizer |
| 167 |
|
makeMinimizer(); |
| 204 |
|
char* molName; |
| 205 |
|
char rbName[100]; |
| 206 |
|
|
| 207 |
< |
//ConstraintPair* consPair; //constraint pair |
| 208 |
< |
//ConstraintElement* consElement1; //first element of constraint pair |
| 209 |
< |
//ConstraintElement* consElement2; //second element of constraint pair |
| 210 |
< |
//int whichRigidBody; |
| 211 |
< |
//int consAtomIndex; //index of constraint atom in rigid body's atom array |
| 212 |
< |
//vector<pair<int, int> > jointAtoms; |
| 207 |
> |
ConstraintPair* consPair; //constraint pair |
| 208 |
> |
ConstraintElement* consElement1; //first element of constraint pair |
| 209 |
> |
ConstraintElement* consElement2; //second element of constraint pair |
| 210 |
> |
int whichRigidBody; |
| 211 |
> |
int consAtomIndex; //index of constraint atom in rigid body's atom array |
| 212 |
> |
vector<pair<int, int> > jointAtoms; |
| 213 |
> |
double bondLength2; |
| 214 |
|
//init the forceField paramters |
| 215 |
|
|
| 216 |
|
the_ff->readParams(); |
| 605 |
|
info[k].integrableObjects.push_back(mySD); |
| 606 |
|
molInfo.myIntegrableObjects.push_back(mySD); |
| 607 |
|
} |
| 608 |
+ |
|
| 609 |
+ |
// send the arrays off to the forceField for init. |
| 610 |
+ |
|
| 611 |
+ |
the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); |
| 612 |
+ |
the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds); |
| 613 |
+ |
the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends); |
| 614 |
+ |
the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions, |
| 615 |
+ |
theTorsions); |
| 616 |
|
|
| 617 |
|
|
| 604 |
– |
/* |
| 605 |
– |
|
| 618 |
|
//creat ConstraintPair. |
| 619 |
< |
molInfo.myConstraintPair.clear(); |
| 619 |
> |
molInfo.myConstraintPairs.clear(); |
| 620 |
|
|
| 621 |
|
for (j = 0; j < molInfo.nBonds; j++){ |
| 622 |
|
|
| 623 |
< |
//if both atoms are in the same rigid body, just skip it |
| 624 |
< |
currentBond = comp_stamps[stampID]->getBond(j); |
| 613 |
< |
if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){ |
| 614 |
< |
|
| 615 |
< |
tempI = currentBond->getA() + atomOffset; |
| 616 |
< |
if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex)) |
| 617 |
< |
consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
| 618 |
< |
else |
| 619 |
< |
consElement1 = new ConstraintAtom(info[k].atoms[tempI]); |
| 623 |
> |
//if bond is constrained bond, add it into constraint pair |
| 624 |
> |
if(molInfo.myBonds[j]->is_constrained()){ |
| 625 |
|
|
| 626 |
< |
tempJ = currentBond->getB() + atomOffset; |
| 627 |
< |
if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex)) |
| 628 |
< |
consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
| 629 |
< |
else |
| 625 |
< |
consElement2 = new ConstraintAtom(info[k].atoms[tempJ]); |
| 626 |
> |
//if both atoms are in the same rigid body, just skip it |
| 627 |
> |
currentBond = comp_stamps[stampID]->getBond(j); |
| 628 |
> |
|
| 629 |
> |
if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){ |
| 630 |
|
|
| 631 |
< |
consPair = new DistanceConstraintPair(consElement1, consElement2); |
| 632 |
< |
molInfo.myConstraintPairs.push_back(consPair); |
| 633 |
< |
} |
| 631 |
> |
tempI = currentBond->getA() + atomOffset; |
| 632 |
> |
if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex)) |
| 633 |
> |
consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
| 634 |
> |
else |
| 635 |
> |
consElement1 = new ConstraintAtom(info[k].atoms[tempI]); |
| 636 |
> |
|
| 637 |
> |
tempJ = currentBond->getB() + atomOffset; |
| 638 |
> |
if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex)) |
| 639 |
> |
consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
| 640 |
> |
else |
| 641 |
> |
consElement2 = new ConstraintAtom(info[k].atoms[tempJ]); |
| 642 |
> |
|
| 643 |
> |
bondLength2 = molInfo.myBonds[j]->get_constraint()->get_dsqr(); |
| 644 |
> |
consPair = new DistanceConstraintPair(consElement1, consElement2, bondLength2); |
| 645 |
> |
|
| 646 |
> |
molInfo.myConstraintPairs.push_back(consPair); |
| 647 |
> |
} |
| 648 |
> |
}//end if(molInfo.myBonds[j]->is_constrained()) |
| 649 |
|
} |
| 650 |
|
|
| 651 |
< |
//loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair |
| 651 |
> |
//loop over rigid bodies, if two rigid bodies share same joint, creat a JointConstraintPair |
| 652 |
|
for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){ |
| 653 |
|
for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){ |
| 654 |
|
|
| 664 |
|
|
| 665 |
|
} |
| 666 |
|
} |
| 648 |
– |
|
| 649 |
– |
*/ |
| 650 |
– |
// send the arrays off to the forceField for init. |
| 667 |
|
|
| 652 |
– |
the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); |
| 653 |
– |
the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds); |
| 654 |
– |
the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends); |
| 655 |
– |
the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions, |
| 656 |
– |
theTorsions); |
| 668 |
|
|
| 669 |
|
info[k].molecules[i].initialize(molInfo); |
| 670 |
|
|
