ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimSetup.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 1167 by tim, Wed May 12 16:38:45 2004 UTC vs.
Revision 1214 by gezelter, Tue Jun 1 18:42:58 2004 UTC

# 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 170 | Line 171 | void SimSetup::makeMolecules(void){
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;
# 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 208 | Line 215 | void SimSetup::makeMolecules(void){
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();
# Line 226 | Line 244 | void SimSetup::makeMolecules(void){
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  
# Line 274 | Line 292 | void SimSetup::makeMolecules(void){
292  
293          molInfo.myAtoms[j]->setType(currentAtom->getType());
294   #ifdef IS_MPI
277 –
295          molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
279 –
296   #endif // is_mpi
297        }
298  
# Line 498 | Line 514 | void SimSetup::makeMolecules(void){
514          nMembers = currentCutoffGroup->getNMembers();
515  
516          myCutoffGroup = new CutoffGroup();
517 +        
518 + #ifdef IS_MPI
519 +        myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
520 + #else
521 +        myCutoffGroup->setGlobalIndex(j + groupOffset);
522 + #endif
523          
524          for (int cg = 0; cg < nMembers; cg++) {
525  
# Line 506 | Line 528 | void SimSetup::makeMolecules(void){
528  
529            // tempI is atom numbering on local processor
530            tempI = molI + atomOffset;
531 +
532 + #ifdef IS_MPI
533 +          globalID = info[k].atoms[tempI]->getGlobalIndex();
534 +          info[k].globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset];
535 + #else
536 +          globalID = info[k].atoms[tempI]->getIndex();
537 +          info[k].globalGroupMembership[globalID] = j + groupOffset;
538 + #endif
539            
510 –          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
540  
541 +          
542 +          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
543 +          
544            cutoffAtomSet.insert(tempI);
545          }
546 <
546 >        
547          molInfo.myCutoffGroups.push_back(myCutoffGroup);
548 <      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
548 >        groupOffset++;
549  
550 +      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
551 +      
552        //creat a cutoff group for every atom  in current molecule which does not belong to cutoffgroup defined at mdl file
553 <
553 >      
554        for(j = 0; j < molInfo.nAtoms; j++){
555 <
555 >        
556          if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
557            myCutoffGroup = new CutoffGroup();
558            myCutoffGroup->addAtom(molInfo.myAtoms[j]);
559 +
560 + #ifdef IS_MPI
561 +          myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
562 +          globalID = info[k].atoms[atomOffset + j]->getGlobalIndex();
563 +          info[k].globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset];
564 + #else
565 +          myCutoffGroup->setGlobalIndex(j + groupOffset);
566 +          globalID = info[k].atoms[atomOffset + j]->getIndex();
567 +          info[k].globalGroupMembership[globalID] = j+groupOffset;
568 + #endif
569            molInfo.myCutoffGroups.push_back(myCutoffGroup);
570 +          groupOffset++;
571          }
572            
573        }
529 –
530 –              
531 –
574  
575        // After this is all set up, scan through the atoms to
576        // see if they can be added to the integrableObjects:
# Line 560 | Line 602 | void SimSetup::makeMolecules(void){
602          info[k].integrableObjects.push_back(mySD);      
603          molInfo.myIntegrableObjects.push_back(mySD);
604        }
605 <    
605 >
606 >
607 >    /*
608 >
609 >      //creat ConstraintPair.
610 >      molInfo.myConstraintPair.clear();
611        
612 +      for (j = 0; j < molInfo.nBonds; j++){
613 +
614 +        //if both atoms are in the same rigid body, just skip it
615 +        currentBond = comp_stamps[stampID]->getBond(j);
616 +        if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){
617 +
618 +          tempI = currentBond->getA() + atomOffset;
619 +          if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex))
620 +            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
621 +          else
622 +             consElement1 = new ConstraintAtom(info[k].atoms[tempI]);      
623 +
624 +          tempJ =  currentBond->getB() + atomOffset;
625 +          if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex))
626 +            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
627 +          else
628 +             consElement2 = new ConstraintAtom(info[k].atoms[tempJ]);    
629 +
630 +          consPair = new DistanceConstraintPair(consElement1, consElement2);
631 +          molInfo.myConstraintPairs.push_back(consPair);
632 +        }
633 +      }  
634 +      
635 +      //loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair
636 +      for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){
637 +        for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){
638 +          
639 +          jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2);
640 +
641 +          for(size_t m = 0; m < jointAtoms.size(); m++){          
642 +            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first);
643 +            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second);
644 +
645 +            consPair = new JointConstraintPair(consElement1, consElement2);  
646 +            molInfo.myConstraintPairs.push_back(consPair);            
647 +          }
648 +
649 +        }
650 +      }
651 +      
652 + */      
653        // send the arrays off to the forceField for init.
654        
655        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
# Line 571 | Line 659 | void SimSetup::makeMolecules(void){
659                                   theTorsions);
660  
661        info[k].molecules[i].initialize(molInfo);
662 <
663 <
662 >      
663 >      
664        atomOffset += molInfo.nAtoms;
665        delete[] theBonds;
666        delete[] theBends;
667        delete[] theTorsions;
668 <    }    
668 >    }
669 >
670 >
671 >
672 > #ifdef IS_MPI    
673 >    // Since the globalGroupMembership has been zero filled and we've only
674 >    // poked values into the atoms we know, we can do an Allreduce
675 >    // to get the full globalGroupMembership array (We think).
676 >    // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
677 >    // docs said we could.
678 >
679 >    int* ggMjunk = new int[mpiSim->getNAtomsGlobal()];    
680 >
681 >    MPI_Allreduce(info[k].globalGroupMembership,
682 >                  ggMjunk,
683 >                  mpiSim->getNAtomsGlobal(),
684 >                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
685 >
686 >    for (i = 0; i < mpiSim->getNAtomsGlobal(); i++)
687 >      info[k].globalGroupMembership[i] = ggMjunk[i];
688 >
689 >    delete[] ggMjunk;
690 >    
691 > #endif
692 >
693 >
694 >
695    }
696  
697   #ifdef IS_MPI
# Line 948 | Line 1062 | void SimSetup::gatherInfo(void){
1062  
1063      info[i].useInitXSstate = globals->getUseInitXSstate();
1064      info[i].orthoTolerance = globals->getOrthoBoxTolerance();
1065 <    
1065 >
1066 >    // check for thermodynamic integration
1067 >    if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) {
1068 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1069 >        info[i].useSolidThermInt = globals->getUseSolidThermInt();
1070 >        info[i].thermIntLambda = globals->getThermIntLambda();
1071 >        info[i].thermIntK = globals->getThermIntK();
1072 >        
1073 >        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
1074 >        info[i].restraint = myRestraint;
1075 >      }
1076 >      else {
1077 >        sprintf(painCave.errMsg,
1078 >                "SimSetup Error:\n"
1079 >                "\tKeyword useSolidThermInt was set to 'true' but\n"
1080 >                "\tthermodynamicIntegrationLambda (and/or\n"
1081 >                "\tthermodynamicIntegrationK) was not specified.\n"
1082 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1083 >        painCave.isFatal = 1;
1084 >        simError();    
1085 >      }
1086 >    }
1087 >    else if(globals->getUseLiquidThermInt()) {
1088 >      if (globals->getUseSolidThermInt()) {
1089 >        sprintf( painCave.errMsg,
1090 >                 "SimSetup Warning: It appears that you have both solid and\n"
1091 >                 "\tliquid thermodynamic integration activated in your .bass\n"
1092 >                 "\tfile. To avoid confusion, specify only one technique in\n"
1093 >                 "\tyour .bass file. Liquid-state thermodynamic integration\n"
1094 >                 "\twill be assumed for the current simulation. If this is not\n"
1095 >                 "\twhat you desire, set useSolidThermInt to 'true' and\n"
1096 >                 "\tuseLiquidThermInt to 'false' in your .bass file.\n");
1097 >        painCave.isFatal = 0;
1098 >        simError();
1099 >      }
1100 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1101 >        info[i].useLiquidThermInt = globals->getUseLiquidThermInt();
1102 >        info[i].thermIntLambda = globals->getThermIntLambda();
1103 >        info[i].thermIntK = globals->getThermIntK();
1104 >      }
1105 >      else {
1106 >        sprintf(painCave.errMsg,
1107 >                "SimSetup Error:\n"
1108 >                "\tKeyword useLiquidThermInt was set to 'true' but\n"
1109 >                "\tthermodynamicIntegrationLambda (and/or\n"
1110 >                "\tthermodynamicIntegrationK) was not specified.\n"
1111 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1112 >        painCave.isFatal = 1;
1113 >        simError();    
1114 >      }
1115 >    }
1116 >    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
1117 >        sprintf(painCave.errMsg,
1118 >                "SimSetup Warning: If you want to use Thermodynamic\n"
1119 >                "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n"
1120 >                "\t'true' in your .bass file.  These keywords are set to\n"
1121 >                "\t'false' by default, so your lambda and/or k values are\n"
1122 >                "\tbeing ignored.\n");
1123 >        painCave.isFatal = 0;
1124 >        simError();  
1125 >    }
1126    }
1127    
1128    //setup seed for random number generator
# Line 1261 | Line 1435 | void SimSetup::makeOutNames(void){
1435          }
1436        }
1437  
1438 +      strcpy(info[k].rawPotName, inFileName);
1439 +      nameLength = strlen(info[k].rawPotName);
1440 +      endTest = &(info[k].rawPotName[nameLength - 5]);
1441 +      if (!strcmp(endTest, ".bass")){
1442 +        strcpy(endTest, ".raw");
1443 +      }
1444 +      else if (!strcmp(endTest, ".BASS")){
1445 +        strcpy(endTest, ".raw");
1446 +      }
1447 +      else{
1448 +        endTest = &(info[k].rawPotName[nameLength - 4]);
1449 +        if (!strcmp(endTest, ".bss")){
1450 +          strcpy(endTest, ".raw");
1451 +        }
1452 +        else if (!strcmp(endTest, ".mdl")){
1453 +          strcpy(endTest, ".raw");
1454 +        }
1455 +        else{
1456 +          strcat(info[k].rawPotName, ".raw");
1457 +        }
1458 +      }
1459 +
1460   #ifdef IS_MPI
1461  
1462      }
# Line 1398 | Line 1594 | void SimSetup::calcSysValues(void){
1594   }
1595  
1596   void SimSetup::calcSysValues(void){
1597 <  int i;
1597 >  int i, j;
1598 >  int ncutgroups, atomsingroups, ngroupsinstamp;
1599  
1600    int* molMembershipArray;
1601 +  CutoffGroupStamp* cg;
1602  
1603    tot_atoms = 0;
1604    tot_bonds = 0;
1605    tot_bends = 0;
1606    tot_torsions = 0;
1607    tot_rigid = 0;
1608 +  tot_groups = 0;
1609    for (i = 0; i < n_components; i++){
1610      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1611      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1612      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1613      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1614      tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1615 +
1616 +    ncutgroups = comp_stamps[i]->getNCutoffGroups();
1617 +    atomsingroups = 0;
1618 +    for (j=0; j < ncutgroups; j++) {
1619 +      cg = comp_stamps[i]->getCutoffGroup(j);
1620 +      atomsingroups += cg->getNMembers();
1621 +    }
1622 +    ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups;
1623 +    tot_groups += components_nmol[i] * ngroupsinstamp;    
1624    }
1625    
1626    tot_SRI = tot_bonds + tot_bends + tot_torsions;
# Line 1425 | Line 1633 | void SimSetup::calcSysValues(void){
1633      info[i].n_torsions = tot_torsions;
1634      info[i].n_SRI = tot_SRI;
1635      info[i].n_mol = tot_nmol;
1636 <
1636 >    info[i].ngroup = tot_groups;
1637      info[i].molMembershipArray = molMembershipArray;
1638    }
1639   }
# Line 1436 | Line 1644 | void SimSetup::mpiMolDivide(void){
1644    int i, j, k;
1645    int localMol, allMol;
1646    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1647 <  int local_rigid;
1647 >  int local_rigid, local_groups;
1648    vector<int> globalMolIndex;
1649 +  int ncutgroups, atomsingroups, ngroupsinstamp;
1650 +  CutoffGroupStamp* cg;
1651  
1652    mpiSim = new mpiSimulation(info);
1653  
1654    mpiSim->divideLabor();
1655    globalAtomIndex = mpiSim->getGlobalAtomIndex();
1656 +  globalGroupIndex = mpiSim->getGlobalGroupIndex();
1657    //globalMolIndex = mpiSim->getGlobalMolIndex();
1658  
1659    // set up the local variables
# Line 1457 | Line 1668 | void SimSetup::mpiMolDivide(void){
1668    local_bends = 0;
1669    local_torsions = 0;
1670    local_rigid = 0;
1671 +  local_groups = 0;
1672    globalAtomCounter = 0;
1673  
1674    for (i = 0; i < n_components; i++){
# Line 1467 | Line 1679 | void SimSetup::mpiMolDivide(void){
1679          local_bends += comp_stamps[i]->getNBends();
1680          local_torsions += comp_stamps[i]->getNTorsions();
1681          local_rigid += comp_stamps[i]->getNRigidBodies();
1682 +
1683 +        ncutgroups = comp_stamps[i]->getNCutoffGroups();
1684 +        atomsingroups = 0;
1685 +        for (k=0; k < ncutgroups; k++) {
1686 +          cg = comp_stamps[i]->getCutoffGroup(k);
1687 +          atomsingroups += cg->getNMembers();
1688 +        }
1689 +        ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups +
1690 +          ncutgroups;
1691 +        local_groups += ngroupsinstamp;    
1692 +
1693          localMol++;
1694        }      
1695        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
# Line 1479 | Line 1702 | void SimSetup::mpiMolDivide(void){
1702    }
1703    local_SRI = local_bonds + local_bends + local_torsions;
1704  
1705 <  info[0].n_atoms = mpiSim->getMyNlocal();  
1705 >  info[0].n_atoms = mpiSim->getNAtomsLocal();  
1706    
1484 –
1707    if (local_atoms != info[0].n_atoms){
1708      sprintf(painCave.errMsg,
1709              "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
# Line 1491 | Line 1713 | void SimSetup::mpiMolDivide(void){
1713      simError();
1714    }
1715  
1716 +  info[0].ngroup = mpiSim->getNGroupsLocal();  
1717 +  if (local_groups != info[0].ngroup){
1718 +    sprintf(painCave.errMsg,
1719 +            "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n"
1720 +            "\tlocalGroups (%d) are not equal.\n",
1721 +            info[0].ngroup, local_groups);
1722 +    painCave.isFatal = 1;
1723 +    simError();
1724 +  }
1725 +  
1726    info[0].n_bonds = local_bonds;
1727    info[0].n_bends = local_bends;
1728    info[0].n_torsions = local_torsions;
# Line 1527 | Line 1759 | void SimSetup::makeSysArrays(void){
1759  
1760  
1761      molIndex = 0;
1762 <    for (i = 0; i < mpiSim->getTotNmol(); i++){
1762 >    for (i = 0; i < mpiSim->getNMolGlobal(); i++){
1763        if (mol2proc[i] == worldRank){
1764          the_molecules[molIndex].setStampID(molCompType[i]);
1765          the_molecules[molIndex].setMyIndex(molIndex);

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines