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 1163 by gezelter, Wed May 12 14:30:12 2004 UTC vs.
Revision 1234 by tim, Fri Jun 4 03:15:31 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 + #include "ConstraintManager.hpp"
17  
18   #ifdef IS_MPI
19   #include "mpiBASS.h"
# Line 156 | Line 158 | void SimSetup::createSim(void){
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();
# Line 170 | Line 176 | void SimSetup::makeMolecules(void){
176    int i, j, k;
177    int exI, exJ, exK, exL, slI, slJ;
178    int tempI, tempJ, tempK, tempL;
179 <  int molI;
180 <  int stampID, atomOffset, rbOffset;
179 >  int molI, globalID;
180 >  int stampID, atomOffset, rbOffset, groupOffset;
181    molInit molInfo;
182    DirectionalAtom* dAtom;
183    RigidBody* myRB;
# Line 185 | Line 191 | void SimSetup::makeMolecules(void){
191    RigidBodyStamp* currentRigidBody;
192    CutoffGroupStamp* currentCutoffGroup;
193    CutoffGroup* myCutoffGroup;
194 <  
194 >  int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file
195 >  set<int> cutoffAtomSet; //atoms belong to  cutoffgroup defined at mdl file
196 >
197    bond_pair* theBonds;
198    bend_set* theBends;
199    torsion_set* theTorsions;
# Line 196 | Line 204 | void SimSetup::makeMolecules(void){
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;
213 +  double bondLength2;
214    //init the forceField paramters
215  
216    the_ff->readParams();
# Line 207 | Line 222 | void SimSetup::makeMolecules(void){
222    for (k = 0; k < nInfo; k++){
223      the_ff->setSimInfo(&(info[k]));
224  
225 + #ifdef IS_MPI
226 +    info[k].globalGroupMembership = new int[mpiSim->getNAtomsGlobal()];
227 +    for (i = 0; i < mpiSim->getNAtomsGlobal(); i++)
228 +      info[k].globalGroupMembership[i] = 0;
229 + #else
230 +    info[k].globalGroupMembership = new int[info[k].n_atoms];
231 +    for (i = 0; i < info[k].n_atoms; i++)
232 +      info[k].globalGroupMembership[i] = 0;
233 + #endif
234 +
235      atomOffset = 0;
236 +    groupOffset = 0;
237  
238      for (i = 0; i < info[k].n_mol; i++){
239        stampID = info[k].molecules[i].getStampID();
# Line 218 | Line 244 | void SimSetup::makeMolecules(void){
244        molInfo.nBends = comp_stamps[stampID]->getNBends();
245        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
246        molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
247 <      molInfo.nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
247 >
248 >      nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
249        
250        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
251  
252        if (molInfo.nBonds > 0)
253 <        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
253 >        molInfo.myBonds = new Bond*[molInfo.nBonds];
254        else
255          molInfo.myBonds = NULL;
256  
257        if (molInfo.nBends > 0)
258 <        molInfo.myBends = new (Bend *) [molInfo.nBends];
258 >        molInfo.myBends = new Bend*[molInfo.nBends];
259        else
260          molInfo.myBends = NULL;
261  
262        if (molInfo.nTorsions > 0)
263 <        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
263 >        molInfo.myTorsions = new Torsion *[molInfo.nTorsions];
264        else
265          molInfo.myTorsions = NULL;
266  
# Line 271 | Line 298 | void SimSetup::makeMolecules(void){
298  
299          molInfo.myAtoms[j]->setType(currentAtom->getType());
300   #ifdef IS_MPI
274 –
301          molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
276 –
302   #endif // is_mpi
303        }
304  
# Line 485 | Line 510 | void SimSetup::makeMolecules(void){
510        
511  
512        //create cutoff group for molecule
513 +
514 +      cutoffAtomSet.clear();
515        molInfo.myCutoffGroups.clear();
516 <      for (j = 0; j < molInfo.nCutoffGroups; j++){
516 >      
517 >      for (j = 0; j < nCutoffGroups; j++){
518  
519          currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
520          nMembers = currentCutoffGroup->getNMembers();
521  
522          myCutoffGroup = new CutoffGroup();
523          
524 + #ifdef IS_MPI
525 +        myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]);
526 + #else
527 +        myCutoffGroup->setGlobalIndex(groupOffset);
528 + #endif
529 +        
530          for (int cg = 0; cg < nMembers; cg++) {
531  
532            // molI is atom numbering inside this molecule
# Line 501 | Line 535 | void SimSetup::makeMolecules(void){
535            // tempI is atom numbering on local processor
536            tempI = molI + atomOffset;
537  
538 <          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
538 > #ifdef IS_MPI
539 >          globalID = info[k].atoms[tempI]->getGlobalIndex();
540 >          info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset];
541 > #else
542 >          globalID = info[k].atoms[tempI]->getIndex();
543 >          info[k].globalGroupMembership[globalID] = groupOffset;
544 > #endif                    
545 >          myCutoffGroup->addAtom(info[k].atoms[tempI]);
546 >          cutoffAtomSet.insert(tempI);
547          }
548 <
548 >        
549          molInfo.myCutoffGroups.push_back(myCutoffGroup);
550 +        groupOffset++;
551 +
552        }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
553        
554 <
554 >      
555 >      // create a cutoff group for every atom in current molecule which
556 >      // does not belong to cutoffgroup defined at mdl file
557 >      
558 >      for(j = 0; j < molInfo.nAtoms; j++){
559 >        
560 >        if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
561 >          myCutoffGroup = new CutoffGroup();
562 >          myCutoffGroup->addAtom(molInfo.myAtoms[j]);
563 >          
564 > #ifdef IS_MPI
565 >          myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]);
566 >          globalID = info[k].atoms[atomOffset + j]->getGlobalIndex();
567 >          info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset];
568 > #else
569 >          myCutoffGroup->setGlobalIndex(groupOffset);
570 >          globalID = info[k].atoms[atomOffset + j]->getIndex();
571 >          info[k].globalGroupMembership[globalID] = groupOffset;
572 > #endif
573 >          molInfo.myCutoffGroups.push_back(myCutoffGroup);
574 >          groupOffset++;
575 >        }          
576 >      }
577  
578        // After this is all set up, scan through the atoms to
579        // see if they can be added to the integrableObjects:
# Line 539 | Line 605 | void SimSetup::makeMolecules(void){
605          info[k].integrableObjects.push_back(mySD);      
606          molInfo.myIntegrableObjects.push_back(mySD);
607        }
608 <    
543 <      
608 >        
609        // send the arrays off to the forceField for init.
610        
611        the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
# Line 549 | Line 614 | void SimSetup::makeMolecules(void){
614        the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
615                                   theTorsions);
616  
552 –      info[k].molecules[i].initialize(molInfo);
617  
618 +      //creat ConstraintPair.
619 +      molInfo.myConstraintPairs.clear();
620 +      
621 +      for (j = 0; j < molInfo.nBonds; j++){
622  
623 +        //if bond is constrained bond, add it into constraint pair
624 +        if(molInfo.myBonds[j]->is_constrained()){
625 +
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 +            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 JointConstraintPair
652 +      for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){
653 +        for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){
654 +          
655 +          jointAtoms = comp_stamps[stampID]->getJointAtoms(rb1, rb2);
656 +
657 +          for(size_t m = 0; m < jointAtoms.size(); m++){          
658 +            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[rb1], jointAtoms[m].first);
659 +            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[rb2], jointAtoms[m].second);
660 +
661 +            consPair = new JointConstraintPair(consElement1, consElement2);  
662 +            molInfo.myConstraintPairs.push_back(consPair);            
663 +          }
664 +
665 +        }
666 +      }
667 +      
668 +
669 +      info[k].molecules[i].initialize(molInfo);
670 +      
671 +      
672        atomOffset += molInfo.nAtoms;
673        delete[] theBonds;
674        delete[] theBends;
675        delete[] theTorsions;
676 <    }    
676 >    }
677 >
678 >
679 >
680 > #ifdef IS_MPI    
681 >    // Since the globalGroupMembership has been zero filled and we've only
682 >    // poked values into the atoms we know, we can do an Allreduce
683 >    // to get the full globalGroupMembership array (We think).
684 >    // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
685 >    // docs said we could.
686 >
687 >    int* ggMjunk = new int[mpiSim->getNAtomsGlobal()];    
688 >
689 >    MPI_Allreduce(info[k].globalGroupMembership,
690 >                  ggMjunk,
691 >                  mpiSim->getNAtomsGlobal(),
692 >                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
693 >
694 >    for (i = 0; i < mpiSim->getNAtomsGlobal(); i++)
695 >      info[k].globalGroupMembership[i] = ggMjunk[i];
696 >
697 >    delete[] ggMjunk;
698 >    
699 > #endif
700 >
701 >
702 >
703    }
704  
705   #ifdef IS_MPI
# Line 927 | Line 1070 | void SimSetup::gatherInfo(void){
1070  
1071      info[i].useInitXSstate = globals->getUseInitXSstate();
1072      info[i].orthoTolerance = globals->getOrthoBoxTolerance();
1073 <    
1073 >
1074 >    // check for thermodynamic integration
1075 >    if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) {
1076 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1077 >        info[i].useSolidThermInt = globals->getUseSolidThermInt();
1078 >        info[i].thermIntLambda = globals->getThermIntLambda();
1079 >        info[i].thermIntK = globals->getThermIntK();
1080 >        
1081 >        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
1082 >        info[i].restraint = myRestraint;
1083 >      }
1084 >      else {
1085 >        sprintf(painCave.errMsg,
1086 >                "SimSetup Error:\n"
1087 >                "\tKeyword useSolidThermInt was set to 'true' but\n"
1088 >                "\tthermodynamicIntegrationLambda (and/or\n"
1089 >                "\tthermodynamicIntegrationK) was not specified.\n"
1090 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1091 >        painCave.isFatal = 1;
1092 >        simError();    
1093 >      }
1094 >    }
1095 >    else if(globals->getUseLiquidThermInt()) {
1096 >      if (globals->getUseSolidThermInt()) {
1097 >        sprintf( painCave.errMsg,
1098 >                 "SimSetup Warning: It appears that you have both solid and\n"
1099 >                 "\tliquid thermodynamic integration activated in your .bass\n"
1100 >                 "\tfile. To avoid confusion, specify only one technique in\n"
1101 >                 "\tyour .bass file. Liquid-state thermodynamic integration\n"
1102 >                 "\twill be assumed for the current simulation. If this is not\n"
1103 >                 "\twhat you desire, set useSolidThermInt to 'true' and\n"
1104 >                 "\tuseLiquidThermInt to 'false' in your .bass file.\n");
1105 >        painCave.isFatal = 0;
1106 >        simError();
1107 >      }
1108 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1109 >        info[i].useLiquidThermInt = globals->getUseLiquidThermInt();
1110 >        info[i].thermIntLambda = globals->getThermIntLambda();
1111 >        info[i].thermIntK = globals->getThermIntK();
1112 >      }
1113 >      else {
1114 >        sprintf(painCave.errMsg,
1115 >                "SimSetup Error:\n"
1116 >                "\tKeyword useLiquidThermInt was set to 'true' but\n"
1117 >                "\tthermodynamicIntegrationLambda (and/or\n"
1118 >                "\tthermodynamicIntegrationK) was not specified.\n"
1119 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1120 >        painCave.isFatal = 1;
1121 >        simError();    
1122 >      }
1123 >    }
1124 >    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
1125 >        sprintf(painCave.errMsg,
1126 >                "SimSetup Warning: If you want to use Thermodynamic\n"
1127 >                "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n"
1128 >                "\t'true' in your .bass file.  These keywords are set to\n"
1129 >                "\t'false' by default, so your lambda and/or k values are\n"
1130 >                "\tbeing ignored.\n");
1131 >        painCave.isFatal = 0;
1132 >        simError();  
1133 >    }
1134    }
1135    
1136    //setup seed for random number generator
# Line 1237 | Line 1440 | void SimSetup::makeOutNames(void){
1440          }
1441          else{
1442            strcat(info[k].statusName, ".stat");
1443 +        }
1444 +      }
1445 +
1446 +      strcpy(info[k].rawPotName, inFileName);
1447 +      nameLength = strlen(info[k].rawPotName);
1448 +      endTest = &(info[k].rawPotName[nameLength - 5]);
1449 +      if (!strcmp(endTest, ".bass")){
1450 +        strcpy(endTest, ".raw");
1451 +      }
1452 +      else if (!strcmp(endTest, ".BASS")){
1453 +        strcpy(endTest, ".raw");
1454 +      }
1455 +      else{
1456 +        endTest = &(info[k].rawPotName[nameLength - 4]);
1457 +        if (!strcmp(endTest, ".bss")){
1458 +          strcpy(endTest, ".raw");
1459          }
1460 +        else if (!strcmp(endTest, ".mdl")){
1461 +          strcpy(endTest, ".raw");
1462 +        }
1463 +        else{
1464 +          strcat(info[k].rawPotName, ".raw");
1465 +        }
1466        }
1467  
1468   #ifdef IS_MPI
# Line 1377 | Line 1602 | void SimSetup::calcSysValues(void){
1602   }
1603  
1604   void SimSetup::calcSysValues(void){
1605 <  int i;
1605 >  int i, j;
1606 >  int ncutgroups, atomsingroups, ngroupsinstamp;
1607  
1608    int* molMembershipArray;
1609 +  CutoffGroupStamp* cg;
1610  
1611    tot_atoms = 0;
1612    tot_bonds = 0;
1613    tot_bends = 0;
1614    tot_torsions = 0;
1615    tot_rigid = 0;
1616 +  tot_groups = 0;
1617    for (i = 0; i < n_components; i++){
1618      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1619      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1620      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1621      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1622      tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1623 +
1624 +    ncutgroups = comp_stamps[i]->getNCutoffGroups();
1625 +    atomsingroups = 0;
1626 +    for (j=0; j < ncutgroups; j++) {
1627 +      cg = comp_stamps[i]->getCutoffGroup(j);
1628 +      atomsingroups += cg->getNMembers();
1629 +    }
1630 +    ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups;
1631 +    tot_groups += components_nmol[i] * ngroupsinstamp;    
1632    }
1633    
1634    tot_SRI = tot_bonds + tot_bends + tot_torsions;
# Line 1404 | Line 1641 | void SimSetup::calcSysValues(void){
1641      info[i].n_torsions = tot_torsions;
1642      info[i].n_SRI = tot_SRI;
1643      info[i].n_mol = tot_nmol;
1644 <
1644 >    info[i].ngroup = tot_groups;
1645      info[i].molMembershipArray = molMembershipArray;
1646    }
1647   }
# Line 1415 | Line 1652 | void SimSetup::mpiMolDivide(void){
1652    int i, j, k;
1653    int localMol, allMol;
1654    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1655 <  int local_rigid;
1655 >  int local_rigid, local_groups;
1656    vector<int> globalMolIndex;
1657 +  int ncutgroups, atomsingroups, ngroupsinstamp;
1658 +  CutoffGroupStamp* cg;
1659  
1660    mpiSim = new mpiSimulation(info);
1661  
1662    mpiSim->divideLabor();
1663    globalAtomIndex = mpiSim->getGlobalAtomIndex();
1664 +  globalGroupIndex = mpiSim->getGlobalGroupIndex();
1665    //globalMolIndex = mpiSim->getGlobalMolIndex();
1666  
1667    // set up the local variables
# Line 1436 | Line 1676 | void SimSetup::mpiMolDivide(void){
1676    local_bends = 0;
1677    local_torsions = 0;
1678    local_rigid = 0;
1679 +  local_groups = 0;
1680    globalAtomCounter = 0;
1681  
1682    for (i = 0; i < n_components; i++){
# Line 1446 | Line 1687 | void SimSetup::mpiMolDivide(void){
1687          local_bends += comp_stamps[i]->getNBends();
1688          local_torsions += comp_stamps[i]->getNTorsions();
1689          local_rigid += comp_stamps[i]->getNRigidBodies();
1690 +
1691 +        ncutgroups = comp_stamps[i]->getNCutoffGroups();
1692 +        atomsingroups = 0;
1693 +        for (k=0; k < ncutgroups; k++) {
1694 +          cg = comp_stamps[i]->getCutoffGroup(k);
1695 +          atomsingroups += cg->getNMembers();
1696 +        }
1697 +        ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups +
1698 +          ncutgroups;
1699 +        local_groups += ngroupsinstamp;    
1700 +
1701          localMol++;
1702        }      
1703        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
# Line 1458 | Line 1710 | void SimSetup::mpiMolDivide(void){
1710    }
1711    local_SRI = local_bonds + local_bends + local_torsions;
1712  
1713 <  info[0].n_atoms = mpiSim->getMyNlocal();  
1713 >  info[0].n_atoms = mpiSim->getNAtomsLocal();  
1714    
1463 –
1715    if (local_atoms != info[0].n_atoms){
1716      sprintf(painCave.errMsg,
1717              "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
# Line 1470 | Line 1721 | void SimSetup::mpiMolDivide(void){
1721      simError();
1722    }
1723  
1724 +  info[0].ngroup = mpiSim->getNGroupsLocal();  
1725 +  if (local_groups != info[0].ngroup){
1726 +    sprintf(painCave.errMsg,
1727 +            "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n"
1728 +            "\tlocalGroups (%d) are not equal.\n",
1729 +            info[0].ngroup, local_groups);
1730 +    painCave.isFatal = 1;
1731 +    simError();
1732 +  }
1733 +  
1734    info[0].n_bonds = local_bonds;
1735    info[0].n_bends = local_bends;
1736    info[0].n_torsions = local_torsions;
# Line 1506 | Line 1767 | void SimSetup::makeSysArrays(void){
1767  
1768  
1769      molIndex = 0;
1770 <    for (i = 0; i < mpiSim->getTotNmol(); i++){
1770 >    for (i = 0; i < mpiSim->getNMolGlobal(); i++){
1771        if (mol2proc[i] == worldRank){
1772          the_molecules[molIndex].setStampID(molCompType[i]);
1773          the_molecules[molIndex].setMyIndex(molIndex);

Diff Legend

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