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 1211 by tim, Tue Jun 1 15:57:30 2004 UTC vs.
Revision 1234 by tim, Fri Jun 4 03:15:31 2004 UTC

# Line 11 | Line 11
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"
# Line 157 | 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 171 | 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 199 | 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;
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 216 | 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  
# Line 282 | Line 298 | void SimSetup::makeMolecules(void){
298  
299          molInfo.myAtoms[j]->setType(currentAtom->getType());
300   #ifdef IS_MPI
285 –
301          molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
287 –
302   #endif // is_mpi
303        }
304  
# Line 506 | Line 520 | void SimSetup::makeMolecules(void){
520          nMembers = currentCutoffGroup->getNMembers();
521  
522          myCutoffGroup = new CutoffGroup();
509 –        myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
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 517 | Line 536 | void SimSetup::makeMolecules(void){
536            tempI = molI + atomOffset;
537  
538   #ifdef IS_MPI
539 <          globalID = info[k].atoms[tempI]->getGlobalIndex()
539 >          globalID = info[k].atoms[tempI]->getGlobalIndex();
540 >          info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset];
541   #else
542            globalID = info[k].atoms[tempI]->getIndex();
543 < #endif
544 <
545 <          globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset];
526 <
527 <          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
528 <
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 <      //creat a cutoff group for every atom  in current molecule which does not belong to cutoffgroup defined at mdl file
555 <
553 >      
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 <
559 >        
560          if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
561            myCutoffGroup = new CutoffGroup();
562            myCutoffGroup->addAtom(molInfo.myAtoms[j]);
563 <          myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
563 >          
564   #ifdef IS_MPI
565 <          globalID = info[k].atoms[atomOffset + j]->getGlobalIndex()
566 < #else
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
550 –          globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset];
573            molInfo.myCutoffGroups.push_back(myCutoffGroup);
574            groupOffset++;
575 <        }
554 <          
575 >        }          
576        }
577  
578        // After this is all set up, scan through the atoms to
# Line 584 | Line 605 | void SimSetup::makeMolecules(void){
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  
589 –    /*
590 –
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);
598 <        if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){
623 >        //if bond is constrained bond, add it into constraint pair
624 >        if(molInfo.myBonds[j]->is_constrained()){
625  
626 <          tempI = currentBond->getA() + atomOffset;
627 <          if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex))
628 <            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
629 <          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]);    
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            
# Line 631 | Line 665 | void SimSetup::makeMolecules(void){
665          }
666        }
667        
634 – */      
635 –      // send the arrays off to the forceField for init.
636 –      
637 –      the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
638 –      the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
639 –      the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
640 –      the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
641 –                                 theTorsions);
668  
669        info[k].molecules[i].initialize(molInfo);
670 <
671 <
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 1020 | Line 1072 | void SimSetup::gatherInfo(void){
1072      info[i].orthoTolerance = globals->getOrthoBoxTolerance();
1073  
1074      // check for thermodynamic integration
1075 <    if (globals->getUseThermInt()) {
1075 >    if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) {
1076        if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1077 <        info[i].useThermInt = globals->getUseThermInt();
1077 >        info[i].useSolidThermInt = globals->getUseSolidThermInt();
1078          info[i].thermIntLambda = globals->getThermIntLambda();
1079          info[i].thermIntK = globals->getThermIntK();
1080          
# Line 1032 | Line 1084 | void SimSetup::gatherInfo(void){
1084        else {
1085          sprintf(painCave.errMsg,
1086                  "SimSetup Error:\n"
1087 <                "\tKeyword useThermInt was set to 'true' but\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");
# Line 1040 | Line 1092 | void SimSetup::gatherInfo(void){
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 useThermInt to 'true' in your .bass file.\n"
1128 <                "\tThe useThermInt keyword is 'false' by default, so your\n"
1129 <                "\tlambda and/or k values are being ignored.\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      }
# Line 1579 | Line 1661 | void SimSetup::mpiMolDivide(void){
1661  
1662    mpiSim->divideLabor();
1663    globalAtomIndex = mpiSim->getGlobalAtomIndex();
1664 +  globalGroupIndex = mpiSim->getGlobalGroupIndex();
1665    //globalMolIndex = mpiSim->getGlobalMolIndex();
1666  
1667    // set up the local variables

Diff Legend

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