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 1154 by gezelter, Tue May 11 16:00:22 2004 UTC vs.
Revision 1212 by chrisfen, Tue Jun 1 17:15:43 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 183 | Line 184 | void SimSetup::makeMolecules(void){
184    BendStamp* currentBend;
185    TorsionStamp* currentTorsion;
186    RigidBodyStamp* currentRigidBody;
187 +  CutoffGroupStamp* currentCutoffGroup;
188 +  CutoffGroup* myCutoffGroup;
189 +  int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file
190 +  set<int> cutoffAtomSet; //atoms belong to  cutoffgroup defined at mdl file
191  
192    bond_pair* theBonds;
193    bend_set* theBends;
# Line 194 | 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 206 | 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 216 | Line 228 | void SimSetup::makeMolecules(void){
228        molInfo.nBends = comp_stamps[stampID]->getNBends();
229        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
230        molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
231 +
232 +      nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
233        
234        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
235  
236        if (molInfo.nBonds > 0)
237 <        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
237 >        molInfo.myBonds = new Bond*[molInfo.nBonds];
238        else
239          molInfo.myBonds = NULL;
240  
241        if (molInfo.nBends > 0)
242 <        molInfo.myBends = new (Bend *) [molInfo.nBends];
242 >        molInfo.myBends = new Bend*[molInfo.nBends];
243        else
244          molInfo.myBends = NULL;
245  
246        if (molInfo.nTorsions > 0)
247 <        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
247 >        molInfo.myTorsions = new Torsion *[molInfo.nTorsions];
248        else
249          molInfo.myTorsions = NULL;
250  
# Line 480 | Line 494 | void SimSetup::makeMolecules(void){
494          info[k].rigidBodies.push_back(myRB);
495        }
496        
497 +
498 +      //create cutoff group for molecule
499 +
500 +      cutoffAtomSet.clear();
501 +      molInfo.myCutoffGroups.clear();
502 +      
503 +      for (j = 0; j < nCutoffGroups; j++){
504 +
505 +        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
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 +
513 +          // molI is atom numbering inside this molecule
514 +          molI = currentCutoffGroup->getMember(cg);    
515 +
516 +          // tempI is atom numbering on local processor
517 +          tempI = molI + atomOffset;
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 +      
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
538 +
539 +      for(j = 0; j < molInfo.nAtoms; j++){
540 +
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  
557        // After this is all set up, scan through the atoms to
558        // see if they can be added to the integrableObjects:
# Line 511 | 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 536 | Line 655 | void SimSetup::makeMolecules(void){
655    MPIcheckPoint();
656   #endif // is_mpi
657  
539 –  // clean up the forcefield
540 –
541 –  if (!globals->haveRcut()){
542 –
543 –    the_ff->calcRcut();
544 –
545 –  } else {
546 –    
547 –    the_ff->setRcut( globals->getRcut() );
548 –  }
549 –
550 –  the_ff->cleanMe();
658   }
659  
660   void SimSetup::initFromBass(void){
# Line 911 | Line 1018 | void SimSetup::gatherInfo(void){
1018  
1019      info[i].useInitXSstate = globals->getUseInitXSstate();
1020      info[i].orthoTolerance = globals->getOrthoBoxTolerance();
1021 <    
1021 >
1022 >    // check for thermodynamic integration
1023 >    if (globals->getUseSolidThermInt() && !globals->getUseLiquidThermInt()) {
1024 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1025 >        info[i].useSolidThermInt = globals->getUseSolidThermInt();
1026 >        info[i].thermIntLambda = globals->getThermIntLambda();
1027 >        info[i].thermIntK = globals->getThermIntK();
1028 >        
1029 >        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
1030 >        info[i].restraint = myRestraint;
1031 >      }
1032 >      else {
1033 >        sprintf(painCave.errMsg,
1034 >                "SimSetup Error:\n"
1035 >                "\tKeyword useSolidThermInt was set to 'true' but\n"
1036 >                "\tthermodynamicIntegrationLambda (and/or\n"
1037 >                "\tthermodynamicIntegrationK) was not specified.\n"
1038 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1039 >        painCave.isFatal = 1;
1040 >        simError();    
1041 >      }
1042 >    }
1043 >    else if(globals->getUseLiquidThermInt()) {
1044 >      if (globals->getUseSolidThermInt()) {
1045 >        sprintf( painCave.errMsg,
1046 >                 "SimSetup Warning: It appears that you have both solid and\n"
1047 >                 "\tliquid thermodynamic integration activated in your .bass\n"
1048 >                 "\tfile. To avoid confusion, specify only one technique in\n"
1049 >                 "\tyour .bass file. Liquid-state thermodynamic integration\n"
1050 >                 "\twill be assumed for the current simulation. If this is not\n"
1051 >                 "\twhat you desire, set useSolidThermInt to 'true' and\n"
1052 >                 "\tuseLiquidThermInt to 'false' in your .bass file.\n");
1053 >        painCave.isFatal = 0;
1054 >        simError();
1055 >      }
1056 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
1057 >        info[i].useLiquidThermInt = globals->getUseLiquidThermInt();
1058 >        info[i].thermIntLambda = globals->getThermIntLambda();
1059 >        info[i].thermIntK = globals->getThermIntK();
1060 >      }
1061 >      else {
1062 >        sprintf(painCave.errMsg,
1063 >                "SimSetup Error:\n"
1064 >                "\tKeyword useLiquidThermInt was set to 'true' but\n"
1065 >                "\tthermodynamicIntegrationLambda (and/or\n"
1066 >                "\tthermodynamicIntegrationK) was not specified.\n"
1067 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
1068 >        painCave.isFatal = 1;
1069 >        simError();    
1070 >      }
1071 >    }
1072 >    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
1073 >        sprintf(painCave.errMsg,
1074 >                "SimSetup Warning: If you want to use Thermodynamic\n"
1075 >                "\tIntegration, set useSolidThermInt or useLiquidThermInt to\n"
1076 >                "\t'true' in your .bass file.  These keywords are set to\n"
1077 >                "\t'false' by default, so your lambda and/or k values are\n"
1078 >                "\tbeing ignored.\n");
1079 >        painCave.isFatal = 0;
1080 >        simError();  
1081 >    }
1082    }
1083    
1084    //setup seed for random number generator
# Line 988 | Line 1155 | void SimSetup::finalInfoCheck(void){
1155   #endif //is_mpi
1156  
1157      double theRcut, theRsw;
1158 +
1159 +    if (globals->haveRcut()) {
1160 +      theRcut = globals->getRcut();
1161 +
1162 +      if (globals->haveRsw())
1163 +        theRsw = globals->getRsw();
1164 +      else
1165 +        theRsw = theRcut;
1166 +      
1167 +      info[i].setDefaultRcut(theRcut, theRsw);
1168 +
1169 +    } else {
1170 +      
1171 +      the_ff->calcRcut();
1172 +      theRcut = info[i].getRcut();
1173  
1174 +      if (globals->haveRsw())
1175 +        theRsw = globals->getRsw();
1176 +      else
1177 +        theRsw = theRcut;
1178 +      
1179 +      info[i].setDefaultRcut(theRcut, theRsw);
1180 +    }
1181 +
1182      if (globals->getUseRF()){
1183        info[i].useReactionField = 1;
1184 <
1184 >      
1185        if (!globals->haveRcut()){
1186          sprintf(painCave.errMsg,
1187                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
# Line 1068 | Line 1258 | void SimSetup::finalInfoCheck(void){
1258    strcpy(checkPointMsg, "post processing checks out");
1259    MPIcheckPoint();
1260   #endif // is_mpi
1261 +
1262 +  // clean up the forcefield
1263 +  the_ff->cleanMe();
1264   }
1265    
1266   void SimSetup::initSystemCoords(void){
# Line 1195 | Line 1388 | void SimSetup::makeOutNames(void){
1388          }
1389          else{
1390            strcat(info[k].statusName, ".stat");
1391 +        }
1392 +      }
1393 +
1394 +      strcpy(info[k].rawPotName, inFileName);
1395 +      nameLength = strlen(info[k].rawPotName);
1396 +      endTest = &(info[k].rawPotName[nameLength - 5]);
1397 +      if (!strcmp(endTest, ".bass")){
1398 +        strcpy(endTest, ".raw");
1399 +      }
1400 +      else if (!strcmp(endTest, ".BASS")){
1401 +        strcpy(endTest, ".raw");
1402 +      }
1403 +      else{
1404 +        endTest = &(info[k].rawPotName[nameLength - 4]);
1405 +        if (!strcmp(endTest, ".bss")){
1406 +          strcpy(endTest, ".raw");
1407          }
1408 +        else if (!strcmp(endTest, ".mdl")){
1409 +          strcpy(endTest, ".raw");
1410 +        }
1411 +        else{
1412 +          strcat(info[k].rawPotName, ".raw");
1413 +        }
1414        }
1415  
1416   #ifdef IS_MPI
# Line 1282 | Line 1497 | void SimSetup::compList(void){
1497    LinkedMolStamp* headStamp = new LinkedMolStamp();
1498    LinkedMolStamp* currentStamp = NULL;
1499    comp_stamps = new MoleculeStamp * [n_components];
1500 +  bool haveCutoffGroups;
1501  
1502 +  haveCutoffGroups = false;
1503 +  
1504    // make an array of molecule stamps that match the components used.
1505    // also extract the used stamps out into a separate linked list
1506  
# Line 1317 | Line 1535 | void SimSetup::compList(void){
1535        headStamp->add(currentStamp);
1536        comp_stamps[i] = headStamp->match(id);
1537      }
1538 +
1539 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1540 +      haveCutoffGroups = true;    
1541    }
1542 +    
1543 +  for (i = 0; i < nInfo; i++)
1544 +    info[i].haveCutoffGroups = haveCutoffGroups;
1545  
1546   #ifdef IS_MPI
1547    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
# Line 1326 | Line 1550 | void SimSetup::calcSysValues(void){
1550   }
1551  
1552   void SimSetup::calcSysValues(void){
1553 <  int i;
1553 >  int i, j;
1554 >  int ncutgroups, atomsingroups, ngroupsinstamp;
1555  
1556    int* molMembershipArray;
1557 +  CutoffGroupStamp* cg;
1558  
1559    tot_atoms = 0;
1560    tot_bonds = 0;
1561    tot_bends = 0;
1562    tot_torsions = 0;
1563    tot_rigid = 0;
1564 +  tot_groups = 0;
1565    for (i = 0; i < n_components; i++){
1566      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1567      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1568      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1569      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1570      tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1571 +
1572 +    ncutgroups = comp_stamps[i]->getNCutoffGroups();
1573 +    atomsingroups = 0;
1574 +    for (j=0; j < ncutgroups; j++) {
1575 +      cg = comp_stamps[i]->getCutoffGroup(j);
1576 +      atomsingroups += cg->getNMembers();
1577 +    }
1578 +    ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups + ncutgroups;
1579 +    tot_groups += components_nmol[i] * ngroupsinstamp;    
1580    }
1581    
1582    tot_SRI = tot_bonds + tot_bends + tot_torsions;
# Line 1353 | Line 1589 | void SimSetup::calcSysValues(void){
1589      info[i].n_torsions = tot_torsions;
1590      info[i].n_SRI = tot_SRI;
1591      info[i].n_mol = tot_nmol;
1592 <
1592 >    info[i].ngroup = tot_groups;
1593      info[i].molMembershipArray = molMembershipArray;
1594    }
1595   }
# Line 1364 | Line 1600 | void SimSetup::mpiMolDivide(void){
1600    int i, j, k;
1601    int localMol, allMol;
1602    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1603 <  int local_rigid;
1603 >  int local_rigid, local_groups;
1604    vector<int> globalMolIndex;
1605 +  int ncutgroups, atomsingroups, ngroupsinstamp;
1606 +  CutoffGroupStamp* cg;
1607  
1608    mpiSim = new mpiSimulation(info);
1609  
# Line 1385 | Line 1623 | void SimSetup::mpiMolDivide(void){
1623    local_bends = 0;
1624    local_torsions = 0;
1625    local_rigid = 0;
1626 +  local_groups = 0;
1627    globalAtomCounter = 0;
1628  
1629    for (i = 0; i < n_components; i++){
# Line 1395 | Line 1634 | void SimSetup::mpiMolDivide(void){
1634          local_bends += comp_stamps[i]->getNBends();
1635          local_torsions += comp_stamps[i]->getNTorsions();
1636          local_rigid += comp_stamps[i]->getNRigidBodies();
1637 +
1638 +        ncutgroups = comp_stamps[i]->getNCutoffGroups();
1639 +        atomsingroups = 0;
1640 +        for (k=0; k < ncutgroups; k++) {
1641 +          cg = comp_stamps[i]->getCutoffGroup(k);
1642 +          atomsingroups += cg->getNMembers();
1643 +        }
1644 +        ngroupsinstamp = comp_stamps[i]->getNAtoms() - atomsingroups +
1645 +          ncutgroups;
1646 +        local_groups += ngroupsinstamp;    
1647 +
1648          localMol++;
1649        }      
1650        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
# Line 1407 | Line 1657 | void SimSetup::mpiMolDivide(void){
1657    }
1658    local_SRI = local_bonds + local_bends + local_torsions;
1659  
1660 <  info[0].n_atoms = mpiSim->getMyNlocal();  
1660 >  info[0].n_atoms = mpiSim->getNAtomsLocal();  
1661    
1412 –
1662    if (local_atoms != info[0].n_atoms){
1663      sprintf(painCave.errMsg,
1664              "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
# Line 1419 | Line 1668 | void SimSetup::mpiMolDivide(void){
1668      simError();
1669    }
1670  
1671 +  info[0].ngroup = mpiSim->getNGroupsLocal();  
1672 +  if (local_groups != info[0].ngroup){
1673 +    sprintf(painCave.errMsg,
1674 +            "SimSetup error: mpiSim's localGroups (%d) and SimSetup's\n"
1675 +            "\tlocalGroups (%d) are not equal.\n",
1676 +            info[0].ngroup, local_groups);
1677 +    painCave.isFatal = 1;
1678 +    simError();
1679 +  }
1680 +  
1681    info[0].n_bonds = local_bonds;
1682    info[0].n_bends = local_bends;
1683    info[0].n_torsions = local_torsions;
# Line 1455 | Line 1714 | void SimSetup::makeSysArrays(void){
1714  
1715  
1716      molIndex = 0;
1717 <    for (i = 0; i < mpiSim->getTotNmol(); i++){
1717 >    for (i = 0; i < mpiSim->getNMolGlobal(); i++){
1718        if (mol2proc[i] == worldRank){
1719          the_molecules[molIndex].setStampID(molCompType[i]);
1720          the_molecules[molIndex].setMyIndex(molIndex);
# Line 1492 | Line 1751 | void SimSetup::makeSysArrays(void){
1751      info[l].atoms = the_atoms;
1752      info[l].molecules = the_molecules;
1753      info[l].nGlobalExcludes = 0;
1754 <
1754 >    
1755      the_ff->setSimInfo(info);
1756    }
1757   }

Diff Legend

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