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 1203 by gezelter, Thu May 27 18:59:17 2004 UTC vs.
Revision 1211 by tim, Tue Jun 1 15:57:30 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 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 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 1500 | 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 1521 | 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 1531 | 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 1545 | Line 1629 | void SimSetup::mpiMolDivide(void){
1629  
1630    info[0].n_atoms = mpiSim->getNAtomsLocal();  
1631    
1548 –
1632    if (local_atoms != info[0].n_atoms){
1633      sprintf(painCave.errMsg,
1634              "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
# Line 1555 | 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;

Diff Legend

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