# | Line 10 | Line 10 | |
---|---|---|
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" | |
# | Line 198 | Line 199 | void SimSetup::makeMolecules(void){ | |
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(); | |
# | Line 210 | Line 217 | void SimSetup::makeMolecules(void){ | |
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(); | |
# | Line 226 | Line 234 | void SimSetup::makeMolecules(void){ | |
234 | molInfo.myAtoms = &(info[k].atoms[atomOffset]); | |
235 | ||
236 | if (molInfo.nBonds > 0) | |
237 | < | molInfo.myBonds = new (Bond *) [molInfo.nBonds]; |
237 | > | molInfo.myBonds = new Bond*[molInfo.nBonds]; |
238 | else | |
239 | molInfo.myBonds = NULL; | |
240 | ||
241 | if (molInfo.nBends > 0) | |
242 | < | molInfo.myBends = new (Bend *) [molInfo.nBends]; |
242 | > | molInfo.myBends = new Bend*[molInfo.nBends]; |
243 | else | |
244 | molInfo.myBends = NULL; | |
245 | ||
246 | if (molInfo.nTorsions > 0) | |
247 | < | molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions]; |
247 | > | molInfo.myTorsions = new Torsion *[molInfo.nTorsions]; |
248 | else | |
249 | molInfo.myTorsions = NULL; | |
250 | ||
# | Line 498 | Line 506 | void SimSetup::makeMolecules(void){ | |
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 | ||
# | Line 506 | Line 515 | void SimSetup::makeMolecules(void){ | |
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 | |
# | Line 522 | Line 541 | void SimSetup::makeMolecules(void){ | |
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 | ||
# | Line 560 | Line 584 | void SimSetup::makeMolecules(void){ | |
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); | |
# | Line 948 | Line 1018 | void SimSetup::gatherInfo(void){ | |
1018 | ||
1019 | info[i].useInitXSstate = globals->getUseInitXSstate(); | |
1020 | info[i].orthoTolerance = globals->getOrthoBoxTolerance(); | |
1021 | < | |
1021 | > | |
1022 | > | // check for thermodynamic integration |
1023 | > | if (globals->getUseThermInt()) { |
1024 | > | if (globals->haveThermIntLambda() && globals->haveThermIntK()) { |
1025 | > | info[i].useThermInt = globals->getUseThermInt(); |
1026 | > | info[i].thermIntLambda = globals->getThermIntLambda(); |
1027 | > | info[i].thermIntK = globals->getThermIntK(); |
1028 | > | |
1029 | > | Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK); |
1030 | > | info[i].restraint = myRestraint; |
1031 | > | } |
1032 | > | else { |
1033 | > | sprintf(painCave.errMsg, |
1034 | > | "SimSetup Error:\n" |
1035 | > | "\tKeyword useThermInt was set to 'true' but\n" |
1036 | > | "\tthermodynamicIntegrationLambda (and/or\n" |
1037 | > | "\tthermodynamicIntegrationK) was not specified.\n" |
1038 | > | "\tPlease provide a lambda value and k value in your .bass file.\n"); |
1039 | > | painCave.isFatal = 1; |
1040 | > | simError(); |
1041 | > | } |
1042 | > | } |
1043 | > | else if(globals->haveThermIntLambda() || globals->haveThermIntK()){ |
1044 | > | sprintf(painCave.errMsg, |
1045 | > | "SimSetup Warning: If you want to use Thermodynamic\n" |
1046 | > | "\tIntegration, set useThermInt to 'true' in your .bass file.\n" |
1047 | > | "\tThe useThermInt keyword is 'false' by default, so your\n" |
1048 | > | "\tlambda and/or k values are being ignored.\n"); |
1049 | > | painCave.isFatal = 0; |
1050 | > | simError(); |
1051 | > | } |
1052 | } | |
1053 | ||
1054 | //setup seed for random number generator | |
# | Line 1261 | Line 1361 | void SimSetup::makeOutNames(void){ | |
1361 | } | |
1362 | } | |
1363 | ||
1364 | + | strcpy(info[k].rawPotName, inFileName); |
1365 | + | nameLength = strlen(info[k].rawPotName); |
1366 | + | endTest = &(info[k].rawPotName[nameLength - 5]); |
1367 | + | if (!strcmp(endTest, ".bass")){ |
1368 | + | strcpy(endTest, ".raw"); |
1369 | + | } |
1370 | + | else if (!strcmp(endTest, ".BASS")){ |
1371 | + | strcpy(endTest, ".raw"); |
1372 | + | } |
1373 | + | else{ |
1374 | + | endTest = &(info[k].rawPotName[nameLength - 4]); |
1375 | + | if (!strcmp(endTest, ".bss")){ |
1376 | + | strcpy(endTest, ".raw"); |
1377 | + | } |
1378 | + | else if (!strcmp(endTest, ".mdl")){ |
1379 | + | strcpy(endTest, ".raw"); |
1380 | + | } |
1381 | + | else{ |
1382 | + | strcat(info[k].rawPotName, ".raw"); |
1383 | + | } |
1384 | + | } |
1385 | + | |
1386 | #ifdef IS_MPI | |
1387 | ||
1388 | } | |
# | Line 1398 | Line 1520 | void SimSetup::calcSysValues(void){ | |
1520 | } | |
1521 | ||
1522 | void SimSetup::calcSysValues(void){ | |
1523 | < | int i; |
1523 | > | int i, j; |
1524 | > | int ncutgroups, atomsingroups, ngroupsinstamp; |
1525 | ||
1526 | int* molMembershipArray; | |
1527 | + | CutoffGroupStamp* cg; |
1528 | ||
1529 | tot_atoms = 0; | |
1530 | tot_bonds = 0; | |
1531 | tot_bends = 0; | |
1532 | tot_torsions = 0; | |
1533 | tot_rigid = 0; | |
1534 | + | tot_groups = 0; |
1535 | for (i = 0; i < n_components; i++){ | |
1536 | tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); | |
1537 | tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); | |
1538 | tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); | |
1539 | tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); | |
1540 | tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies(); | |
1541 | + | |
1542 | + | ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1543 | + | atomsingroups = 0; |
1544 | + | for (j=0; j < ncutgroups; j++) { |
1545 | + | cg = comp_stamps[i]->getCutoffGroup(j); |
1546 | + | atomsingroups += cg->getNMembers(); |
1547 | + | } |
1548 | + | ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups; |
1549 | + | tot_groups += components_nmol[i] * ngroupsinstamp; |
1550 | } | |
1551 | ||
1552 | tot_SRI = tot_bonds + tot_bends + tot_torsions; | |
# | Line 1425 | Line 1559 | void SimSetup::calcSysValues(void){ | |
1559 | info[i].n_torsions = tot_torsions; | |
1560 | info[i].n_SRI = tot_SRI; | |
1561 | info[i].n_mol = tot_nmol; | |
1562 | < | |
1562 | > | info[i].ngroup = tot_groups; |
1563 | info[i].molMembershipArray = molMembershipArray; | |
1564 | } | |
1565 | } | |
# | Line 1436 | Line 1570 | void SimSetup::mpiMolDivide(void){ | |
1570 | int i, j, k; | |
1571 | int localMol, allMol; | |
1572 | int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; | |
1573 | < | int local_rigid; |
1573 | > | int local_rigid, local_groups; |
1574 | vector<int> globalMolIndex; | |
1575 | + | int ncutgroups, atomsingroups, ngroupsinstamp; |
1576 | + | CutoffGroupStamp* cg; |
1577 | ||
1578 | mpiSim = new mpiSimulation(info); | |
1579 | ||
# | Line 1457 | Line 1593 | void SimSetup::mpiMolDivide(void){ | |
1593 | local_bends = 0; | |
1594 | local_torsions = 0; | |
1595 | local_rigid = 0; | |
1596 | + | local_groups = 0; |
1597 | globalAtomCounter = 0; | |
1598 | ||
1599 | for (i = 0; i < n_components; i++){ | |
# | Line 1467 | Line 1604 | void SimSetup::mpiMolDivide(void){ | |
1604 | local_bends += comp_stamps[i]->getNBends(); | |
1605 | local_torsions += comp_stamps[i]->getNTorsions(); | |
1606 | local_rigid += comp_stamps[i]->getNRigidBodies(); | |
1607 | + | |
1608 | + | ncutgroups = comp_stamps[i]->getNCutoffGroups(); |
1609 | + | atomsingroups = 0; |
1610 | + | for (k=0; k < ncutgroups; k++) { |
1611 | + | cg = comp_stamps[i]->getCutoffGroup(k); |
1612 | + | atomsingroups += cg->getNMembers(); |
1613 | + | } |
1614 | + | ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + |
1615 | + | ncutgroups; |
1616 | + | local_groups += ngroupsinstamp; |
1617 | + | |
1618 | localMol++; | |
1619 | } | |
1620 | for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){ | |
# | Line 1479 | Line 1627 | void SimSetup::mpiMolDivide(void){ | |
1627 | } | |
1628 | local_SRI = local_bonds + local_bends + local_torsions; | |
1629 | ||
1630 | < | info[0].n_atoms = mpiSim->getMyNlocal(); |
1630 | > | info[0].n_atoms = mpiSim->getNAtomsLocal(); |
1631 | ||
1484 | – | |
1632 | if (local_atoms != info[0].n_atoms){ | |
1633 | sprintf(painCave.errMsg, | |
1634 | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" | |
# | Line 1491 | Line 1638 | void SimSetup::mpiMolDivide(void){ | |
1638 | simError(); | |
1639 | } | |
1640 | ||
1641 | + | info[0].ngroup = mpiSim->getNGroupsLocal(); |
1642 | + | if (local_groups != info[0].ngroup){ |
1643 | + | sprintf(painCave.errMsg, |
1644 | + | "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n" |
1645 | + | "\tlocalGroups (%d) are not equal.\n", |
1646 | + | info[0].ngroup, local_groups); |
1647 | + | painCave.isFatal = 1; |
1648 | + | simError(); |
1649 | + | } |
1650 | + | |
1651 | info[0].n_bonds = local_bonds; | |
1652 | info[0].n_bends = local_bends; | |
1653 | info[0].n_torsions = local_torsions; | |
# | Line 1527 | Line 1684 | void SimSetup::makeSysArrays(void){ | |
1684 | ||
1685 | ||
1686 | molIndex = 0; | |
1687 | < | for (i = 0; i < mpiSim->getTotNmol(); i++){ |
1687 | > | for (i = 0; i < mpiSim->getNMolGlobal(); i++){ |
1688 | if (mol2proc[i] == worldRank){ | |
1689 | the_molecules[molIndex].setStampID(molCompType[i]); | |
1690 | the_molecules[molIndex].setMyIndex(molIndex); |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |