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 1214 by gezelter, Tue Jun 1 18:42:58 2004 UTC vs.
Revision 1452 by tim, Mon Aug 23 15:11:36 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 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 516 | Line 522 | void SimSetup::makeMolecules(void){
522          myCutoffGroup = new CutoffGroup();
523          
524   #ifdef IS_MPI
525 <        myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
525 >        myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]);
526   #else
527 <        myCutoffGroup->setGlobalIndex(j + groupOffset);
527 >        myCutoffGroup->setGlobalIndex(groupOffset);
528   #endif
529          
530          for (int cg = 0; cg < nMembers; cg++) {
# Line 531 | Line 537 | void SimSetup::makeMolecules(void){
537  
538   #ifdef IS_MPI
539            globalID = info[k].atoms[tempI]->getGlobalIndex();
540 <          info[k].globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset];
540 >          info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset];
541   #else
542            globalID = info[k].atoms[tempI]->getIndex();
543 <          info[k].globalGroupMembership[globalID] = j + groupOffset;
544 < #endif
545 <          
540 <
541 <          
542 <          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
543 <          
543 >          info[k].globalGroupMembership[globalID] = groupOffset;
544 > #endif                    
545 >          myCutoffGroup->addAtom(info[k].atoms[tempI]);
546            cutoffAtomSet.insert(tempI);
547          }
548          
# Line 549 | Line 551 | void SimSetup::makeMolecules(void){
551  
552        }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
553        
552 –      //creat a cutoff group for every atom  in current molecule which does not belong to cutoffgroup defined at mdl file
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 <
563 >          
564   #ifdef IS_MPI
565 <          myCutoffGroup->setGlobalIndex(globalGroupIndex[j + groupOffset]);
565 >          myCutoffGroup->setGlobalIndex(globalGroupIndex[groupOffset]);
566            globalID = info[k].atoms[atomOffset + j]->getGlobalIndex();
567 <          info[k].globalGroupMembership[globalID] = globalGroupIndex[j+groupOffset];
567 >          info[k].globalGroupMembership[globalID] = globalGroupIndex[groupOffset];
568   #else
569 <          myCutoffGroup->setGlobalIndex(j + groupOffset);
569 >          myCutoffGroup->setGlobalIndex(groupOffset);
570            globalID = info[k].atoms[atomOffset + j]->getIndex();
571 <          info[k].globalGroupMembership[globalID] = j+groupOffset;
571 >          info[k].globalGroupMembership[globalID] = groupOffset;
572   #endif
573            molInfo.myCutoffGroups.push_back(myCutoffGroup);
574            groupOffset++;
575 <        }
572 <          
575 >        }          
576        }
577  
578        // After this is all set up, scan through the atoms to
# Line 602 | 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  
607 –    /*
608 –
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);
616 <        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
622 <             consElement1 = new ConstraintAtom(info[k].atoms[tempI]);      
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 <          tempJ =  currentBond->getB() + atomOffset;
632 <          if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex))
633 <            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
634 <          else
635 <             consElement2 = new ConstraintAtom(info[k].atoms[tempJ]);    
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 <          consPair = new DistanceConstraintPair(consElement1, consElement2);
638 <          molInfo.myConstraintPairs.push_back(consPair);
639 <        }
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 649 | Line 665 | void SimSetup::makeMolecules(void){
665          }
666        }
667        
652 – */      
653 –      // send the arrays off to the forceField for init.
654 –      
655 –      the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
656 –      the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
657 –      the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
658 –      the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
659 –                                 theTorsions);
668  
669        info[k].molecules[i].initialize(molInfo);
670        
# Line 908 | Line 916 | void SimSetup::gatherInfo(void){
916           painCave.isFatal = 1;
917           simError();
918    }
919 <
920 <    // get the ensemble
919 >  if (globals->haveForceFieldVariant()) {
920 >    strcpy(forcefield_variant, globals->getForceFieldVariant());
921 >    has_forcefield_variant = 1;
922 >  }
923 >  
924 >  // get the ensemble
925  
926    strcpy(ensemble, globals->getEnsemble());
927  
# Line 1506 | Line 1518 | void SimSetup::createFF(void){
1518   void SimSetup::createFF(void){
1519    switch (ffCase){
1520      case FF_DUFF:
1521 <      the_ff = new DUFF();
1521 >        the_ff = new DUFF();
1522        break;
1523  
1524      case FF_LJ:
# Line 1514 | Line 1526 | void SimSetup::createFF(void){
1526        break;
1527  
1528      case FF_EAM:
1529 <      the_ff = new EAM_FF();
1529 >      if (has_forcefield_variant)
1530 >        the_ff = new EAM_FF(forcefield_variant);
1531 >      else
1532 >        the_ff = new EAM_FF();
1533        break;
1534  
1535      case FF_H2O:
# Line 1528 | Line 1543 | void SimSetup::createFF(void){
1543        simError();
1544    }
1545  
1546 +
1547   #ifdef IS_MPI
1548    strcpy(checkPointMsg, "ForceField creation successful");
1549    MPIcheckPoint();
# Line 1804 | Line 1820 | void SimSetup::makeIntegrator(void){
1820   void SimSetup::makeIntegrator(void){
1821    int k;
1822  
1823 <  NVE<RealIntegrator>* myNVE = NULL;
1824 <  NVT<RealIntegrator>* myNVT = NULL;
1825 <  NPTi<NPT<RealIntegrator> >* myNPTi = NULL;
1826 <  NPTf<NPT<RealIntegrator> >* myNPTf = NULL;
1827 <  NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL;
1823 >  NVE<Integrator<BaseIntegrator> >* myNVE = NULL;
1824 >  NVT<Integrator<BaseIntegrator> >* myNVT = NULL;
1825 >  NPTi<NPT<Integrator<BaseIntegrator> > >* myNPTi = NULL;
1826 >  NPTf<NPT<Integrator<BaseIntegrator> > >* myNPTf = NULL;
1827 >  NPTxyz<NPT<Integrator<BaseIntegrator> > >* myNPTxyz = NULL;
1828    
1829    for (k = 0; k < nInfo; k++){
1830      switch (ensembleCase){
# Line 1818 | Line 1834 | void SimSetup::makeIntegrator(void){
1834            myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff);
1835          }
1836          else{
1837 <          myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff);
1837 >          if (globals->haveQuaternion()){
1838 >            if (globals->getUseQuaternion())
1839 >              info->the_integrator = new NVE<SQSIntegrator<RealIntegrator> >(&(info[k]), the_ff);
1840 >          }
1841 >          else
1842 >            info->the_integrator = new NVE<RealIntegrator>(&(info[k]), the_ff);
1843 >          break;
1844 >
1845 >          //myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff);
1846          }
1847          
1848          info->the_integrator = myNVE;

Diff Legend

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