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" |
171 |
|
int i, j, k; |
172 |
|
int exI, exJ, exK, exL, slI, slJ; |
173 |
|
int tempI, tempJ, tempK, tempL; |
174 |
< |
int molI; |
175 |
< |
int stampID, atomOffset, rbOffset; |
174 |
> |
int molI, globalID; |
175 |
> |
int stampID, atomOffset, rbOffset, groupOffset; |
176 |
|
molInit molInfo; |
177 |
|
DirectionalAtom* dAtom; |
178 |
|
RigidBody* myRB; |
186 |
|
RigidBodyStamp* currentRigidBody; |
187 |
|
CutoffGroupStamp* currentCutoffGroup; |
188 |
|
CutoffGroup* myCutoffGroup; |
189 |
< |
|
189 |
> |
int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file |
190 |
> |
set<int> cutoffAtomSet; //atoms belong to cutoffgroup defined at mdl file |
191 |
> |
|
192 |
|
bond_pair* theBonds; |
193 |
|
bend_set* theBends; |
194 |
|
torsion_set* theTorsions; |
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(); |
215 |
|
|
216 |
|
for (k = 0; k < nInfo; k++){ |
217 |
|
the_ff->setSimInfo(&(info[k])); |
218 |
+ |
|
219 |
+ |
#ifdef IS_MPI |
220 |
+ |
info[k].globalGroupMembership = new int[mpiSim->getNAtomsGlobal()]; |
221 |
+ |
for (i = 0; i < mpiSim->getNAtomsGlobal(); i++) |
222 |
+ |
info[k].globalGroupMembership[i] = 0; |
223 |
+ |
#else |
224 |
+ |
info[k].globalGroupMembership = new int[info[k].n_atoms]; |
225 |
+ |
for (i = 0; i < info[k].n_atoms; i++) |
226 |
+ |
info[k].globalGroupMembership[i] = 0; |
227 |
+ |
#endif |
228 |
|
|
229 |
|
atomOffset = 0; |
230 |
+ |
groupOffset = 0; |
231 |
|
|
232 |
|
for (i = 0; i < info[k].n_mol; i++){ |
233 |
|
stampID = info[k].molecules[i].getStampID(); |
238 |
|
molInfo.nBends = comp_stamps[stampID]->getNBends(); |
239 |
|
molInfo.nTorsions = comp_stamps[stampID]->getNTorsions(); |
240 |
|
molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies(); |
241 |
< |
molInfo.nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups(); |
241 |
> |
|
242 |
> |
nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups(); |
243 |
|
|
244 |
|
molInfo.myAtoms = &(info[k].atoms[atomOffset]); |
245 |
|
|
246 |
|
if (molInfo.nBonds > 0) |
247 |
< |
molInfo.myBonds = new (Bond *) [molInfo.nBonds]; |
247 |
> |
molInfo.myBonds = new Bond*[molInfo.nBonds]; |
248 |
|
else |
249 |
|
molInfo.myBonds = NULL; |
250 |
|
|
251 |
|
if (molInfo.nBends > 0) |
252 |
< |
molInfo.myBends = new (Bend *) [molInfo.nBends]; |
252 |
> |
molInfo.myBends = new Bend*[molInfo.nBends]; |
253 |
|
else |
254 |
|
molInfo.myBends = NULL; |
255 |
|
|
256 |
|
if (molInfo.nTorsions > 0) |
257 |
< |
molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions]; |
257 |
> |
molInfo.myTorsions = new Torsion *[molInfo.nTorsions]; |
258 |
|
else |
259 |
|
molInfo.myTorsions = NULL; |
260 |
|
|
292 |
|
|
293 |
|
molInfo.myAtoms[j]->setType(currentAtom->getType()); |
294 |
|
#ifdef IS_MPI |
274 |
– |
|
295 |
|
molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]); |
276 |
– |
|
296 |
|
#endif // is_mpi |
297 |
|
} |
298 |
|
|
503 |
|
} |
504 |
|
|
505 |
|
|
506 |
< |
//creat cutoff group for molecule |
506 |
> |
//create cutoff group for molecule |
507 |
> |
|
508 |
> |
cutoffAtomSet.clear(); |
509 |
|
molInfo.myCutoffGroups.clear(); |
510 |
< |
for (j = 0; j < molInfo.nCutoffGroups; j++){ |
510 |
> |
|
511 |
> |
for (j = 0; j < nCutoffGroups; j++){ |
512 |
|
|
513 |
|
currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j); |
514 |
|
nMembers = currentCutoffGroup->getNMembers(); |
515 |
|
|
516 |
|
myCutoffGroup = new CutoffGroup(); |
517 |
|
|
518 |
+ |
#ifdef IS_MPI |
519 |
+ |
myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); |
520 |
+ |
#else |
521 |
+ |
myCutoffGroup->setGlobalIndex(groupOffset); |
522 |
+ |
#endif |
523 |
+ |
|
524 |
|
for (int cg = 0; cg < nMembers; cg++) { |
525 |
|
|
526 |
|
// molI is atom numbering inside this molecule |
529 |
|
// tempI is atom numbering on local processor |
530 |
|
tempI = molI + atomOffset; |
531 |
|
|
532 |
< |
myCutoffGroup->addAtom(info[k].atoms[tempI]); |
532 |
> |
#ifdef IS_MPI |
533 |
> |
globalID = info[k].atoms[tempI]->getGlobalIndex(); |
534 |
> |
info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; |
535 |
> |
#else |
536 |
> |
globalID = info[k].atoms[tempI]->getIndex(); |
537 |
> |
info[k].globalGroupMembership[globalID] = groupOffset; |
538 |
> |
#endif |
539 |
> |
myCutoffGroup->addAtom(info[k].atoms[tempI]); |
540 |
> |
cutoffAtomSet.insert(tempI); |
541 |
|
} |
542 |
< |
|
542 |
> |
|
543 |
|
molInfo.myCutoffGroups.push_back(myCutoffGroup); |
544 |
+ |
groupOffset++; |
545 |
+ |
|
546 |
|
}//end for (j = 0; j < molInfo.nCutoffGroups; j++) |
547 |
|
|
548 |
+ |
|
549 |
+ |
// create a cutoff group for every atom in current molecule which |
550 |
+ |
// does not belong to cutoffgroup defined at mdl file |
551 |
+ |
|
552 |
+ |
for(j = 0; j < molInfo.nAtoms; j++){ |
553 |
+ |
|
554 |
+ |
if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){ |
555 |
+ |
myCutoffGroup = new CutoffGroup(); |
556 |
+ |
myCutoffGroup->addAtom(molInfo.myAtoms[j]); |
557 |
+ |
|
558 |
+ |
#ifdef IS_MPI |
559 |
+ |
myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]); |
560 |
+ |
globalID = info[k].atoms[atomOffset + j]->getGlobalIndex(); |
561 |
+ |
info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset]; |
562 |
+ |
#else |
563 |
+ |
myCutoffGroup->setGlobalIndex(groupOffset); |
564 |
+ |
globalID = info[k].atoms[atomOffset + j]->getIndex(); |
565 |
+ |
info[k].globalGroupMembership[globalID] = groupOffset; |
566 |
+ |
#endif |
567 |
+ |
molInfo.myCutoffGroups.push_back(myCutoffGroup); |
568 |
+ |
groupOffset++; |
569 |
+ |
} |
570 |
+ |
} |
571 |
|
|
511 |
– |
|
572 |
|
// After this is all set up, scan through the atoms to |
573 |
|
// see if they can be added to the integrableObjects: |
574 |
|
|
599 |
|
info[k].integrableObjects.push_back(mySD); |
600 |
|
molInfo.myIntegrableObjects.push_back(mySD); |
601 |
|
} |
602 |
< |
|
602 |
> |
|
603 |
> |
|
604 |
> |
/* |
605 |
> |
|
606 |
> |
//creat ConstraintPair. |
607 |
> |
molInfo.myConstraintPair.clear(); |
608 |
|
|
609 |
+ |
for (j = 0; j < molInfo.nBonds; j++){ |
610 |
+ |
|
611 |
+ |
//if both atoms are in the same rigid body, just skip it |
612 |
+ |
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]); |
620 |
+ |
|
621 |
+ |
tempJ = currentBond->getB() + atomOffset; |
622 |
+ |
if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex)) |
623 |
+ |
consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex); |
624 |
+ |
else |
625 |
+ |
consElement2 = new ConstraintAtom(info[k].atoms[tempJ]); |
626 |
+ |
|
627 |
+ |
consPair = new DistanceConstraintPair(consElement1, consElement2); |
628 |
+ |
molInfo.myConstraintPairs.push_back(consPair); |
629 |
+ |
} |
630 |
+ |
} |
631 |
+ |
|
632 |
+ |
//loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair |
633 |
+ |
for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){ |
634 |
+ |
for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){ |
635 |
+ |
|
636 |
+ |
jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2); |
637 |
+ |
|
638 |
+ |
for(size_t m = 0; m < jointAtoms.size(); m++){ |
639 |
+ |
consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first); |
640 |
+ |
consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second); |
641 |
+ |
|
642 |
+ |
consPair = new JointConstraintPair(consElement1, consElement2); |
643 |
+ |
molInfo.myConstraintPairs.push_back(consPair); |
644 |
+ |
} |
645 |
+ |
|
646 |
+ |
} |
647 |
+ |
} |
648 |
+ |
|
649 |
+ |
*/ |
650 |
|
// send the arrays off to the forceField for init. |
651 |
|
|
652 |
|
the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms); |
656 |
|
theTorsions); |
657 |
|
|
658 |
|
info[k].molecules[i].initialize(molInfo); |
659 |
< |
|
660 |
< |
|
659 |
> |
|
660 |
> |
|
661 |
|
atomOffset += molInfo.nAtoms; |
662 |
|
delete[] theBonds; |
663 |
|
delete[] theBends; |
664 |
|
delete[] theTorsions; |
665 |
< |
} |
560 |
< |
} |
665 |
> |
} |
666 |
|
|
562 |
– |
#ifdef IS_MPI |
563 |
– |
sprintf(checkPointMsg, "all molecules initialized succesfully"); |
564 |
– |
MPIcheckPoint(); |
565 |
– |
#endif // is_mpi |
667 |
|
|
567 |
– |
// clean up the forcefield |
668 |
|
|
669 |
< |
if (!globals->haveRcut()){ |
669 |
> |
#ifdef IS_MPI |
670 |
> |
// Since the globalGroupMembership has been zero filled and we've only |
671 |
> |
// poked values into the atoms we know, we can do an Allreduce |
672 |
> |
// to get the full globalGroupMembership array (We think). |
673 |
> |
// This would be prettier if we could use MPI_IN_PLACE like the MPI-2 |
674 |
> |
// docs said we could. |
675 |
|
|
676 |
< |
the_ff->calcRcut(); |
676 |
> |
int* ggMjunk = new int[mpiSim->getNAtomsGlobal()]; |
677 |
|
|
678 |
< |
} else { |
678 |
> |
MPI_Allreduce(info[k].globalGroupMembership, |
679 |
> |
ggMjunk, |
680 |
> |
mpiSim->getNAtomsGlobal(), |
681 |
> |
MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
682 |
> |
|
683 |
> |
for (i = 0; i < mpiSim->getNAtomsGlobal(); i++) |
684 |
> |
info[k].globalGroupMembership[i] = ggMjunk[i]; |
685 |
> |
|
686 |
> |
delete[] ggMjunk; |
687 |
|
|
688 |
< |
the_ff->setRcut( globals->getRcut() ); |
688 |
> |
#endif |
689 |
> |
|
690 |
> |
|
691 |
> |
|
692 |
|
} |
693 |
|
|
694 |
< |
the_ff->cleanMe(); |
694 |
> |
#ifdef IS_MPI |
695 |
> |
sprintf(checkPointMsg, "all molecules initialized succesfully"); |
696 |
> |
MPIcheckPoint(); |
697 |
> |
#endif // is_mpi |
698 |
> |
|
699 |
|
} |
700 |
|
|
701 |
|
void SimSetup::initFromBass(void){ |
1059 |
|
|
1060 |
|
info[i].useInitXSstate = globals->getUseInitXSstate(); |
1061 |
|
info[i].orthoTolerance = globals->getOrthoBoxTolerance(); |
1062 |
< |
|
1062 |
> |
|
1063 |
> |
// check for thermodynamic integration |
1064 |
> |
if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) { |
1065 |
> |
if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
1066 |
> |
info[i].useSolidThermInt = globals->getUseSolidThermInt(); |
1067 |
> |
info[i].thermIntLambda = globals->getThermIntLambda(); |
1068 |
> |
info[i].thermIntK = globals->getThermIntK(); |
1069 |
> |
|
1070 |
> |
Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK); |
1071 |
> |
info[i].restraint = myRestraint; |
1072 |
> |
} |
1073 |
> |
else { |
1074 |
> |
sprintf(painCave.errMsg, |
1075 |
> |
"SimSetup Error:\n" |
1076 |
> |
"\tKeyword useSolidThermInt was set to 'true' but\n" |
1077 |
> |
"\tthermodynamicIntegrationLambda (and/or\n" |
1078 |
> |
"\tthermodynamicIntegrationK) was not specified.\n" |
1079 |
> |
"\tPlease provide a lambda value and k value in your .bass file.\n"); |
1080 |
> |
painCave.isFatal = 1; |
1081 |
> |
simError(); |
1082 |
> |
} |
1083 |
> |
} |
1084 |
> |
else if(globals->getUseLiquidThermInt()) { |
1085 |
> |
if (globals->getUseSolidThermInt()) { |
1086 |
> |
sprintf( painCave.errMsg, |
1087 |
> |
"SimSetup Warning: It appears that you have both solid and\n" |
1088 |
> |
"\tliquid thermodynamic integration activated in your .bass\n" |
1089 |
> |
"\tfile. To avoid confusion, specify only one technique in\n" |
1090 |
> |
"\tyour .bass file. Liquid-state thermodynamic integration\n" |
1091 |
> |
"\twill be assumed for the current simulation. If this is not\n" |
1092 |
> |
"\twhat you desire, set useSolidThermInt to 'true' and\n" |
1093 |
> |
"\tuseLiquidThermInt to 'false' in your .bass file.\n"); |
1094 |
> |
painCave.isFatal = 0; |
1095 |
> |
simError(); |
1096 |
> |
} |
1097 |
> |
if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
1098 |
> |
info[i].useLiquidThermInt = globals->getUseLiquidThermInt(); |
1099 |
> |
info[i].thermIntLambda = globals->getThermIntLambda(); |
1100 |
> |
info[i].thermIntK = globals->getThermIntK(); |
1101 |
> |
} |
1102 |
> |
else { |
1103 |
> |
sprintf(painCave.errMsg, |
1104 |
> |
"SimSetup Error:\n" |
1105 |
> |
"\tKeyword useLiquidThermInt was set to 'true' but\n" |
1106 |
> |
"\tthermodynamicIntegrationLambda (and/or\n" |
1107 |
> |
"\tthermodynamicIntegrationK) was not specified.\n" |
1108 |
> |
"\tPlease provide a lambda value and k value in your .bass file.\n"); |
1109 |
> |
painCave.isFatal = 1; |
1110 |
> |
simError(); |
1111 |
> |
} |
1112 |
> |
} |
1113 |
> |
else if(globals->haveThermIntLambda() || globals->haveThermIntK()){ |
1114 |
> |
sprintf(painCave.errMsg, |
1115 |
> |
"SimSetup Warning: If you want to use Thermodynamic\n" |
1116 |
> |
"\tIntegration, set useSolidThermInt or useLiquidThermInt to\n" |
1117 |
> |
"\t'true' in your .bass file. These keywords are set to\n" |
1118 |
> |
"\t'false' by default, so your lambda and/or k values are\n" |
1119 |
> |
"\tbeing ignored.\n"); |
1120 |
> |
painCave.isFatal = 0; |
1121 |
> |
simError(); |
1122 |
> |
} |
1123 |
|
} |
1124 |
|
|
1125 |
|
//setup seed for random number generator |
1196 |
|
#endif //is_mpi |
1197 |
|
|
1198 |
|
double theRcut, theRsw; |
1199 |
+ |
|
1200 |
+ |
if (globals->haveRcut()) { |
1201 |
+ |
theRcut = globals->getRcut(); |
1202 |
+ |
|
1203 |
+ |
if (globals->haveRsw()) |
1204 |
+ |
theRsw = globals->getRsw(); |
1205 |
+ |
else |
1206 |
+ |
theRsw = theRcut; |
1207 |
+ |
|
1208 |
+ |
info[i].setDefaultRcut(theRcut, theRsw); |
1209 |
+ |
|
1210 |
+ |
} else { |
1211 |
+ |
|
1212 |
+ |
the_ff->calcRcut(); |
1213 |
+ |
theRcut = info[i].getRcut(); |
1214 |
+ |
|
1215 |
+ |
if (globals->haveRsw()) |
1216 |
+ |
theRsw = globals->getRsw(); |
1217 |
+ |
else |
1218 |
+ |
theRsw = theRcut; |
1219 |
+ |
|
1220 |
+ |
info[i].setDefaultRcut(theRcut, theRsw); |
1221 |
+ |
} |
1222 |
|
|
1223 |
|
if (globals->getUseRF()){ |
1224 |
|
info[i].useReactionField = 1; |
1225 |
< |
|
1225 |
> |
|
1226 |
|
if (!globals->haveRcut()){ |
1227 |
|
sprintf(painCave.errMsg, |
1228 |
|
"SimSetup Warning: No value was set for the cutoffRadius.\n" |
1299 |
|
strcpy(checkPointMsg, "post processing checks out"); |
1300 |
|
MPIcheckPoint(); |
1301 |
|
#endif // is_mpi |
1302 |
+ |
|
1303 |
+ |
// clean up the forcefield |
1304 |
+ |
the_ff->cleanMe(); |
1305 |
|
} |
1306 |
|
|
1307 |
|
void SimSetup::initSystemCoords(void){ |
1429 |
|
} |
1430 |
|
else{ |
1431 |
|
strcat(info[k].statusName, ".stat"); |
1432 |
+ |
} |
1433 |
+ |
} |
1434 |
+ |
|
1435 |
+ |
strcpy(info[k].rawPotName, inFileName); |
1436 |
+ |
nameLength = strlen(info[k].rawPotName); |
1437 |
+ |
endTest = &(info[k].rawPotName[nameLength - 5]); |
1438 |
+ |
if (!strcmp(endTest, ".bass")){ |
1439 |
+ |
strcpy(endTest, ".raw"); |
1440 |
+ |
} |
1441 |
+ |
else if (!strcmp(endTest, ".BASS")){ |
1442 |
+ |
strcpy(endTest, ".raw"); |
1443 |
+ |
} |
1444 |
+ |
else{ |
1445 |
+ |
endTest = &(info[k].rawPotName[nameLength - 4]); |
1446 |
+ |
if (!strcmp(endTest, ".bss")){ |
1447 |
+ |
strcpy(endTest, ".raw"); |
1448 |
|
} |
1449 |
+ |
else if (!strcmp(endTest, ".mdl")){ |
1450 |
+ |
strcpy(endTest, ".raw"); |
1451 |
+ |
} |
1452 |
+ |
else{ |
1453 |
+ |
strcat(info[k].rawPotName, ".raw"); |
1454 |
+ |
} |
1455 |
|
} |
1456 |
|
|
1457 |
|
#ifdef IS_MPI |
1591 |
|
} |
1592 |
|
|
1593 |
|
void SimSetup::calcSysValues(void){ |
1594 |
< |
int i; |
1594 |
> |
int i, j; |
1595 |
> |
int ncutgroups, atomsingroups, ngroupsinstamp; |
1596 |
|
|
1597 |
|
int* molMembershipArray; |
1598 |
+ |
CutoffGroupStamp* cg; |
1599 |
|
|
1600 |
|
tot_atoms = 0; |
1601 |
|
tot_bonds = 0; |
1602 |
|
tot_bends = 0; |
1603 |
|
tot_torsions = 0; |
1604 |
|
tot_rigid = 0; |
1605 |
+ |
tot_groups = 0; |
1606 |
|
for (i = 0; i < n_components; i++){ |
1607 |
|
tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); |
1608 |
|
tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); |
1609 |
|
tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); |
1610 |
|
tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); |
1611 |
|
tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies(); |
1612 |
+ |
|
1613 |
+ |
ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1614 |
+ |
atomsingroups = 0; |
1615 |
+ |
for (j=0; j < ncutgroups; j++) { |
1616 |
+ |
cg = comp_stamps[i]->getCutoffGroup(j); |
1617 |
+ |
atomsingroups += cg->getNMembers(); |
1618 |
+ |
} |
1619 |
+ |
ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups; |
1620 |
+ |
tot_groups += components_nmol[i] * ngroupsinstamp; |
1621 |
|
} |
1622 |
|
|
1623 |
|
tot_SRI = tot_bonds + tot_bends + tot_torsions; |
1630 |
|
info[i].n_torsions = tot_torsions; |
1631 |
|
info[i].n_SRI = tot_SRI; |
1632 |
|
info[i].n_mol = tot_nmol; |
1633 |
< |
|
1633 |
> |
info[i].ngroup = tot_groups; |
1634 |
|
info[i].molMembershipArray = molMembershipArray; |
1635 |
|
} |
1636 |
|
} |
1641 |
|
int i, j, k; |
1642 |
|
int localMol, allMol; |
1643 |
|
int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; |
1644 |
< |
int local_rigid; |
1644 |
> |
int local_rigid, local_groups; |
1645 |
|
vector<int> globalMolIndex; |
1646 |
+ |
int ncutgroups, atomsingroups, ngroupsinstamp; |
1647 |
+ |
CutoffGroupStamp* cg; |
1648 |
|
|
1649 |
|
mpiSim = new mpiSimulation(info); |
1650 |
|
|
1651 |
|
mpiSim->divideLabor(); |
1652 |
|
globalAtomIndex = mpiSim->getGlobalAtomIndex(); |
1653 |
+ |
globalGroupIndex = mpiSim->getGlobalGroupIndex(); |
1654 |
|
//globalMolIndex = mpiSim->getGlobalMolIndex(); |
1655 |
|
|
1656 |
|
// set up the local variables |
1665 |
|
local_bends = 0; |
1666 |
|
local_torsions = 0; |
1667 |
|
local_rigid = 0; |
1668 |
+ |
local_groups = 0; |
1669 |
|
globalAtomCounter = 0; |
1670 |
|
|
1671 |
|
for (i = 0; i < n_components; i++){ |
1676 |
|
local_bends += comp_stamps[i]->getNBends(); |
1677 |
|
local_torsions += comp_stamps[i]->getNTorsions(); |
1678 |
|
local_rigid += comp_stamps[i]->getNRigidBodies(); |
1679 |
+ |
|
1680 |
+ |
ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1681 |
+ |
atomsingroups = 0; |
1682 |
+ |
for (k=0; k < ncutgroups; k++) { |
1683 |
+ |
cg = comp_stamps[i]->getCutoffGroup(k); |
1684 |
+ |
atomsingroups += cg->getNMembers(); |
1685 |
+ |
} |
1686 |
+ |
ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + |
1687 |
+ |
ncutgroups; |
1688 |
+ |
local_groups += ngroupsinstamp; |
1689 |
+ |
|
1690 |
|
localMol++; |
1691 |
|
} |
1692 |
|
for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ |
1699 |
|
} |
1700 |
|
local_SRI = local_bonds + local_bends + local_torsions; |
1701 |
|
|
1702 |
< |
info[0].n_atoms = mpiSim->getMyNlocal(); |
1702 |
> |
info[0].n_atoms = mpiSim->getNAtomsLocal(); |
1703 |
|
|
1449 |
– |
|
1704 |
|
if (local_atoms != info[0].n_atoms){ |
1705 |
|
sprintf(painCave.errMsg, |
1706 |
|
"SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" |
1710 |
|
simError(); |
1711 |
|
} |
1712 |
|
|
1713 |
+ |
info[0].ngroup = mpiSim->getNGroupsLocal(); |
1714 |
+ |
if (local_groups != info[0].ngroup){ |
1715 |
+ |
sprintf(painCave.errMsg, |
1716 |
+ |
"SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n" |
1717 |
+ |
"\tlocalGroups (%d) are not equal.\n", |
1718 |
+ |
info[0].ngroup, local_groups); |
1719 |
+ |
painCave.isFatal = 1; |
1720 |
+ |
simError(); |
1721 |
+ |
} |
1722 |
+ |
|
1723 |
|
info[0].n_bonds = local_bonds; |
1724 |
|
info[0].n_bends = local_bends; |
1725 |
|
info[0].n_torsions = local_torsions; |
1756 |
|
|
1757 |
|
|
1758 |
|
molIndex = 0; |
1759 |
< |
for (i = 0; i < mpiSim->getTotNmol(); i++){ |
1759 |
> |
for (i = 0; i < mpiSim->getNMolGlobal(); i++){ |
1760 |
|
if (mol2proc[i] == worldRank){ |
1761 |
|
the_molecules[molIndex].setStampID(molCompType[i]); |
1762 |
|
the_molecules[molIndex].setMyIndex(molIndex); |