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 670 by mmeineke, Thu Aug 7 21:47:18 2003 UTC vs.
Revision 699 by tim, Fri Aug 15 19:24:13 2003 UTC

# Line 5 | Line 5
5   #include <string>
6  
7   #include "SimSetup.hpp"
8 + #include "ReadWrite.hpp"
9   #include "parse_me.h"
10   #include "Integrator.hpp"
11   #include "simError.h"
# Line 22 | Line 23
23   #define NPTf_ENS       3
24   #define NPTim_ENS      4
25   #define NPTfm_ENS      5
25 – #define NVEZCONS_ENS   6
26 – #define NVTZCONS_ENS   7
27 – #define NPTiZCONS_ENS  8
28 – #define NPTfZCONS_ENS  9
29 – #define NPTimZCONS_ENS 10
30 – #define NPTfmZCONS_ENS 11
26  
27   #define FF_DUFF 0
28   #define FF_LJ   1
# Line 102 | Line 97 | void SimSetup::createSim(void){
97    int i, j, k, globalAtomIndex;
98    
99    // gather all of the information from the Bass file
100 <  
100 >
101    gatherInfo();
102  
103    // creation of complex system objects
# Line 110 | Line 105 | void SimSetup::createSim(void){
105    sysObjectsCreation();
106  
107    // check on the post processing info
108 <  
108 >
109    finalInfoCheck();
110  
111    // initialize the system coordinates
# Line 639 | Line 634 | void SimSetup::gatherInfo( void ){
634    else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
635    else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
636    else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS;
642 –
643 –  else if( !strcasecmp( ensemble, "NVEZCONS")) ensembleCase = NVEZCONS_ENS;
644 –  else if( !strcasecmp( ensemble, "NVTZCONS"))  ensembleCase = NVTZCONS_ENS;
645 –  else if( !strcasecmp( ensemble, "NPTiZCONS") || !strcasecmp( ensemble, "NPTZCONS"))
646 –    ensembleCase = NPTiZCONS_ENS;
647 –  else if( !strcasecmp( ensemble, "NPTfZCONS"))  ensembleCase = NPTfZCONS_ENS;
648 –  else if( !strcasecmp( ensemble, "NPTimZCONS"))  ensembleCase = NPTimZCONS_ENS;
649 –  else if( !strcasecmp( ensemble, "NPTfmZCONS"))  ensembleCase = NPTfmZCONS_ENS;
650 –  
637    else{
638      sprintf( painCave.errMsg,
639               "SimSetup Warning. Unrecognized Ensemble -> %s, "
# Line 903 | Line 889 | void SimSetup::initSystemCoords( void ){
889   void SimSetup::initSystemCoords( void ){
890    int i;
891    
892 <  std::cerr << "Setting atom Coords\n";
892 >  char* inName;
893  
894 +
895    (info[0].getConfiguration())->createArrays( info[0].n_atoms );
896    
897    for(i=0; i<info[0].n_atoms; i++) info[0].atoms[i]->setCoords();
# Line 915 | Line 902 | void SimSetup::initSystemCoords( void ){
902   #ifdef IS_MPI // is_mpi
903      if( worldRank == 0 ){
904   #endif //is_mpi
905 <      fileInit = new InitializeFromFile( globals->getInitialConfig() );
905 >      inName = globals->getInitialConfig();
906 >      double* tempDouble = new double[1000000];
907 >      fileInit = new InitializeFromFile( inName );
908   #ifdef IS_MPI
909      }else fileInit = new InitializeFromFile( NULL );
910   #endif
# Line 1048 | Line 1037 | void SimSetup::sysObjectsCreation( void ){
1037    int i,k;
1038    
1039    // create the forceField
1040 <  
1040 >
1041    createFF();
1042  
1043    // extract componentList
# Line 1066 | Line 1055 | void SimSetup::sysObjectsCreation( void ){
1055   #endif //is_mpi
1056    
1057    // create the atom and SRI arrays. Also initialize Molecule Stamp ID's
1058 <  
1058 >
1059    makeSysArrays();
1060  
1061    // make and initialize the molecules (all but atomic coordinates)
1062 <  
1062 >
1063    makeMolecules();
1064    
1065    for(k=0; k<nInfo; k++){
# Line 1367 | Line 1356 | void SimSetup::makeIntegrator( void ){
1356    NPTf<RealIntegrator>* myNPTf = NULL;
1357    NPTim<RealIntegrator>* myNPTim = NULL;
1358    NPTfm<RealIntegrator>* myNPTfm = NULL;
1370 –  ZConstraint<NVE<RealIntegrator> >* myNVEZCons = NULL;
1371 –  ZConstraint<NVT<RealIntegrator> >* myNVTZCons = NULL;
1372 –  ZConstraint<NPTi<RealIntegrator> >* myNPTiZCons = NULL;
1373 –  ZConstraint<NPTf<RealIntegrator> >* myNPTfZCons = NULL;
1374 –  ZConstraint<NPTim<RealIntegrator> >* myNPTimZCons = NULL;
1375 –  ZConstraint<NPTfm<RealIntegrator> >* myNPTfmZCons = NULL;
1359          
1360    for(k=0; k<nInfo; k++){
1361      
1362      switch( ensembleCase ){
1363        
1364      case NVE_ENS:
1365 <      new NVE<RealIntegrator>( &(info[k]), the_ff );
1365 >        if (globals->haveZconstraints()){
1366 >         setupZConstraint(info[k]);
1367 >           new ZConstraint<NVE<RealIntegrator> >( &(info[k]), the_ff );
1368 >        }
1369 >
1370 >        else
1371 >        new NVE<RealIntegrator>( &(info[k]), the_ff );
1372        break;
1373        
1374      case NVT_ENS:
1375 <      myNVT = new NVT<RealIntegrator>( &(info[k]), the_ff );
1375 >        if (globals->haveZconstraints()){
1376 >         setupZConstraint(info[k]);
1377 >           myNVT = new ZConstraint<NVT<RealIntegrator> >( &(info[k]), the_ff );
1378 >        }
1379 >        else
1380 >        myNVT = new NVT<RealIntegrator>( &(info[k]), the_ff );
1381 >
1382        myNVT->setTargetTemp(globals->getTargetTemp());
1383        
1384        if (globals->haveTauThermostat())
# Line 1399 | Line 1394 | void SimSetup::makeIntegrator( void ){
1394        break;
1395        
1396      case NPTi_ENS:
1397 <      myNPTi = new NPTi<RealIntegrator>( &(info[k]), the_ff );
1398 <      myNPTi->setTargetTemp( globals->getTargetTemp() );
1397 >        if (globals->haveZconstraints()){
1398 >         setupZConstraint(info[k]);
1399 >           myNPTi = new ZConstraint<NPTi<RealIntegrator> >( &(info[k]), the_ff );
1400 >        }
1401 >        else
1402 >        myNPTi = new NPTi<RealIntegrator>( &(info[k]), the_ff );
1403 >
1404 >        myNPTi->setTargetTemp( globals->getTargetTemp() );
1405        
1406        if (globals->haveTargetPressure())
1407          myNPTi->setTargetPressure(globals->getTargetPressure());
# Line 1434 | Line 1435 | void SimSetup::makeIntegrator( void ){
1435        break;
1436        
1437      case NPTf_ENS:
1438 <      myNPTf = new NPTf<RealIntegrator>( &(info[k]), the_ff );
1438 >        if (globals->haveZconstraints()){
1439 >         setupZConstraint(info[k]);
1440 >           myNPTf = new ZConstraint<NPTf<RealIntegrator> >( &(info[k]), the_ff );
1441 >        }
1442 >        else
1443 >        myNPTf = new NPTf<RealIntegrator>( &(info[k]), the_ff );
1444 >
1445        myNPTf->setTargetTemp( globals->getTargetTemp());
1446        
1447        if (globals->haveTargetPressure())
# Line 1469 | Line 1476 | void SimSetup::makeIntegrator( void ){
1476        break;
1477        
1478      case NPTim_ENS:
1479 <      myNPTim = new NPTim<RealIntegrator>( &(info[k]), the_ff );
1480 <      myNPTim->setTargetTemp( globals->getTargetTemp());
1479 >        if (globals->haveZconstraints()){
1480 >         setupZConstraint(info[k]);
1481 >           myNPTim = new ZConstraint<NPTim<RealIntegrator> >( &(info[k]), the_ff );
1482 >        }
1483 >        else
1484 >        myNPTim = new NPTim<RealIntegrator>( &(info[k]), the_ff );
1485 >
1486 >        myNPTim->setTargetTemp( globals->getTargetTemp());
1487        
1488        if (globals->haveTargetPressure())
1489          myNPTim->setTargetPressure(globals->getTargetPressure());
# Line 1504 | Line 1517 | void SimSetup::makeIntegrator( void ){
1517        break;
1518        
1519      case NPTfm_ENS:
1520 <      myNPTfm = new NPTfm<RealIntegrator>( &(info[k]), the_ff );
1521 <      myNPTfm->setTargetTemp( globals->getTargetTemp());
1520 >        if (globals->haveZconstraints()){
1521 >         setupZConstraint(info[k]);
1522 >           myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >( &(info[k]), the_ff );
1523 >        }
1524 >        else
1525 >        myNPTfm = new NPTfm<RealIntegrator>( &(info[k]), the_ff );
1526 >
1527 >        myNPTfm->setTargetTemp( globals->getTargetTemp());
1528        
1529        if (globals->haveTargetPressure())
1530          myNPTfm->setTargetPressure(globals->getTargetPressure());
# Line 1538 | Line 1557 | void SimSetup::makeIntegrator( void ){
1557        }
1558        break;
1559        
1541 –    case NVEZCONS_ENS:
1542 –      
1543 –      
1544 –      //setup index of z-constraint molecules, z-constraint sampel time
1545 –      //and z-constraint force output name. These parameter should be known
1546 –      //before constructing the z-constraint integrator
1547 –      setupZConstraint();
1548 –      
1549 –      myNVEZCons = new ZConstraint<NVE<RealIntegrator> >( &(info[k]), the_ff );
1550 –      
1551 –      break;
1552 –      
1553 –      
1554 –    case NVTZCONS_ENS:
1555 –      
1556 –      setupZConstraint();
1557 –      
1558 –      myNVTZCons = new ZConstraint<NVT<RealIntegrator> >( &(info[k]), the_ff );
1559 –      myNVTZCons->setTargetTemp(globals->getTargetTemp());
1560 –      
1561 –      if (globals->haveTauThermostat())
1562 –        myNVTZCons->setTauThermostat(globals->getTauThermostat());
1563 –      
1564 –      else {
1565 –        sprintf( painCave.errMsg,
1566 –                 "SimSetup error: If you use the NVT\n"
1567 –                 "    ensemble, you must set tauThermostat.\n");
1568 –        painCave.isFatal = 1;
1569 –        simError();
1570 –      }    
1571 –      break;    
1572 –      
1573 –    case NPTiZCONS_ENS:
1574 –      
1575 –      setupZConstraint();
1576 –      
1577 –      myNPTiZCons = new ZConstraint<NPTi<RealIntegrator> >( &(info[k]), the_ff );
1578 –      myNPTiZCons->setTargetTemp( globals->getTargetTemp() );
1579 –      
1580 –      if (globals->haveTargetPressure())
1581 –        myNPTiZCons->setTargetPressure(globals->getTargetPressure());
1582 –      else {
1583 –        sprintf( painCave.errMsg,
1584 –                 "SimSetup error: If you use a constant pressure\n"
1585 –                 "    ensemble, you must set targetPressure in the BASS file.\n");
1586 –        painCave.isFatal = 1;
1587 –        simError();
1588 –      }
1589 –      
1590 –      if( globals->haveTauThermostat() )
1591 –        myNPTiZCons->setTauThermostat( globals->getTauThermostat() );
1592 –      else{
1593 –        sprintf( painCave.errMsg,
1594 –                 "SimSetup error: If you use an NPT\n"
1595 –                 "    ensemble, you must set tauThermostat.\n");
1596 –        painCave.isFatal = 1;
1597 –        simError();
1598 –      }
1599 –      
1600 –      if( globals->haveTauBarostat() )
1601 –        myNPTiZCons->setTauBarostat( globals->getTauBarostat() );
1602 –      else{
1603 –        sprintf( painCave.errMsg,
1604 –                 "SimSetup error: If you use an NPT\n"
1605 –                 "    ensemble, you must set tauBarostat.\n");
1606 –        painCave.isFatal = 1;
1607 –        simError();
1608 –      }  
1609 –      
1610 –      break;
1611 –      
1612 –    case NPTfZCONS_ENS:
1613 –      
1614 –      setupZConstraint();
1615 –      
1616 –      myNPTfZCons = new ZConstraint<NPTf<RealIntegrator> >( &(info[k]), the_ff );
1617 –      myNPTfZCons->setTargetTemp( globals->getTargetTemp());
1618 –      
1619 –      if (globals->haveTargetPressure())
1620 –        myNPTfZCons->setTargetPressure(globals->getTargetPressure());
1621 –      else {
1622 –        sprintf( painCave.errMsg,
1623 –                 "SimSetup error: If you use a constant pressure\n"
1624 –                 "    ensemble, you must set targetPressure in the BASS file.\n");
1625 –        painCave.isFatal = 1;
1626 –        simError();
1627 –      }    
1628 –      
1629 –      if( globals->haveTauThermostat() )
1630 –        myNPTfZCons->setTauThermostat( globals->getTauThermostat() );
1631 –      else{
1632 –        sprintf( painCave.errMsg,
1633 –                 "SimSetup error: If you use an NPT\n"
1634 –                 "    ensemble, you must set tauThermostat.\n");
1635 –        painCave.isFatal = 1;
1636 –        simError();
1637 –      }
1638 –      
1639 –      if( globals->haveTauBarostat() )
1640 –        myNPTfZCons->setTauBarostat( globals->getTauBarostat() );
1641 –      else{
1642 –        sprintf( painCave.errMsg,
1643 –                 "SimSetup error: If you use an NPT\n"
1644 –                 "    ensemble, you must set tauBarostat.\n");
1645 –        painCave.isFatal = 1;
1646 –        simError();
1647 –      }  
1648 –      
1649 –      break;  
1650 –      
1651 –    case NPTimZCONS_ENS:
1652 –      
1653 –      setupZConstraint();
1654 –      
1655 –      myNPTimZCons = new ZConstraint<NPTim<RealIntegrator> >( &(info[k]), the_ff );
1656 –      myNPTimZCons->setTargetTemp( globals->getTargetTemp());
1657 –      
1658 –      if (globals->haveTargetPressure())
1659 –        myNPTimZCons->setTargetPressure(globals->getTargetPressure());
1660 –      else {
1661 –        sprintf( painCave.errMsg,
1662 –                 "SimSetup error: If you use a constant pressure\n"
1663 –                 "    ensemble, you must set targetPressure in the BASS file.\n");
1664 –        painCave.isFatal = 1;
1665 –        simError();
1666 –      }
1667 –      
1668 –      if( globals->haveTauThermostat() )
1669 –        myNPTimZCons->setTauThermostat( globals->getTauThermostat() );
1670 –      else{
1671 –        sprintf( painCave.errMsg,
1672 –                 "SimSetup error: If you use an NPT\n"
1673 –                 "    ensemble, you must set tauThermostat.\n");
1674 –        painCave.isFatal = 1;
1675 –        simError();
1676 –      }
1677 –      
1678 –      if( globals->haveTauBarostat() )
1679 –        myNPTimZCons->setTauBarostat( globals->getTauBarostat() );
1680 –      else{
1681 –        sprintf( painCave.errMsg,
1682 –                 "SimSetup error: If you use an NPT\n"
1683 –                 "    ensemble, you must set tauBarostat.\n");
1684 –        painCave.isFatal = 1;
1685 –        simError();
1686 –      }  
1687 –      
1688 –      break;
1689 –      
1690 –    case NPTfmZCONS_ENS:
1691 –      
1692 –      setupZConstraint();
1693 –      
1694 –      myNPTfmZCons = new ZConstraint<NPTfm<RealIntegrator> >( &(info[k]), the_ff );
1695 –      myNPTfmZCons->setTargetTemp( globals->getTargetTemp());
1696 –      
1697 –      if (globals->haveTargetPressure())
1698 –        myNPTfmZCons->setTargetPressure(globals->getTargetPressure());
1699 –      else {
1700 –        sprintf( painCave.errMsg,
1701 –                 "SimSetup error: If you use a constant pressure\n"
1702 –                 "    ensemble, you must set targetPressure in the BASS file.\n");
1703 –        painCave.isFatal = 1;
1704 –        simError();
1705 –      }
1706 –      
1707 –      if( globals->haveTauThermostat() )
1708 –        myNPTfmZCons->setTauThermostat( globals->getTauThermostat() );
1709 –      else{
1710 –        sprintf( painCave.errMsg,
1711 –                 "SimSetup error: If you use an NPT\n"
1712 –                 "    ensemble, you must set tauThermostat.\n");
1713 –        painCave.isFatal = 1;
1714 –        simError();
1715 –      }
1716 –      
1717 –      if( globals->haveTauBarostat() )
1718 –        myNPTfmZCons->setTauBarostat( globals->getTauBarostat() );
1719 –      else{
1720 –        sprintf( painCave.errMsg,
1721 –                 "SimSetup error: If you use an NPT\n"
1722 –                 "    ensemble, you must set tauBarostat.\n");
1723 –        painCave.isFatal = 1;
1724 –        simError();
1725 –      }    
1726 –      break;      
1727 –      
1728 –      
1729 –      
1560      default:
1561        sprintf( painCave.errMsg,
1562                 "SimSetup Error. Unrecognized ensemble in case statement.\n");
# Line 1763 | Line 1593 | void SimSetup::setupZConstraint()
1593  
1594   }
1595  
1596 < void SimSetup::setupZConstraint()
1596 > void SimSetup::setupZConstraint(SimInfo& theInfo)
1597   {
1598 <  int k;
1599 <
1600 <  for(k=0; k<nInfo; k++){
1601 <    
1772 <    if(globals->haveZConsTime()){  
1598 >    int nZConstraints;
1599 >    ZconStamp** zconStamp;
1600 >        
1601 >    if(globals->haveZconstraintTime()){  
1602        
1603        //add sample time of z-constraint  into SimInfo's property list                    
1604        DoubleData* zconsTimeProp = new DoubleData();
1605 <      zconsTimeProp->setID("zconstime");
1606 <      zconsTimeProp->setData(globals->getZConsTime());
1607 <      info[k].addProperty(zconsTimeProp);
1605 >      zconsTimeProp->setID(ZCONSTIME_ID);
1606 >      zconsTimeProp->setData(globals->getZconsTime());
1607 >      theInfo.addProperty(zconsTimeProp);
1608      }
1609      else{
1610        sprintf( painCave.errMsg,
# Line 1784 | Line 1613 | void SimSetup::setupZConstraint()
1613        painCave.isFatal = 1;
1614        simError();      
1615      }
1616 <    
1617 <    if(globals->haveIndexOfAllZConsMols()){
1618 <
1619 <      //add index of z-constraint molecules into SimInfo's property list
1620 <      vector<int> tempIndex = globals->getIndexOfAllZConsMols();
1621 <      
1622 <      //sort the index
1794 <      sort(tempIndex.begin(), tempIndex.end());
1795 <      
1796 <      IndexData* zconsIndex = new IndexData();
1797 <      zconsIndex->setID("zconsindex");
1798 <      zconsIndex->setIndexData(tempIndex);
1799 <      info[k].addProperty(zconsIndex);
1616 >
1617 >    //push zconsTol into siminfo, if user does not specify
1618 >    //value for zconsTol, a default value will be used
1619 >    DoubleData* zconsTol = new DoubleData();
1620 >    zconsTol->setID(ZCONSTOL_ID);
1621 >    if(globals->haveZconsTol()){
1622 >      zconsTol->setData(globals->getZconsTol());
1623      }
1624 <    else{
1624 >         else{
1625 >                double defaultZConsTol = 0.01;
1626        sprintf( painCave.errMsg,
1627 <               "SimSetup error: If you use an ZConstraint\n"
1628 <               " , you must set index of z-constraint molecules.\n");
1629 <      painCave.isFatal = 1;
1630 <      simError();    
1631 <      
1632 <    }
1633 <    
1627 >               "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
1628 >               " , default value %f is used.\n", defaultZConsTol);
1629 >      painCave.isFatal = 0;
1630 >      simError();      
1631 >
1632 >      zconsTol->setData(defaultZConsTol);
1633 >         }
1634 >    theInfo.addProperty(zconsTol);
1635 >
1636 >    //set Force Substraction Policy
1637 >    StringData* zconsForcePolicy =  new StringData();
1638 >    zconsForcePolicy->setID(ZCONSFORCEPOLICY_ID);
1639 >                
1640 >         if(globals->haveZconsForcePolicy()){
1641 >      zconsForcePolicy->setData(globals->getZconsForcePolicy());
1642 >         }      
1643 >         else{
1644 >       sprintf( painCave.errMsg,
1645 >               "ZConstraint Warning: User does not set force substraction policy, "
1646 >               "average force substraction policy is used\n");
1647 >       painCave.isFatal = 0;
1648 >       simError();
1649 >                 zconsForcePolicy->setData("BYNUMBER");
1650 >         }
1651 >        
1652 >         theInfo.addProperty(zconsForcePolicy);
1653 >        
1654      //Determine the name of ouput file and add it into SimInfo's property list
1655      //Be careful, do not use inFileName, since it is a pointer which
1656      //point to a string at master node, and slave nodes do not contain that string
1657      
1658 <    string zconsOutput(info[k].finalName);
1658 >    string zconsOutput(theInfo.finalName);
1659      
1660      zconsOutput = zconsOutput.substr(0, zconsOutput.rfind(".")) + ".fz";
1661      
1662      StringData* zconsFilename = new StringData();
1663 <    zconsFilename->setID("zconsfilename");
1663 >    zconsFilename->setID(ZCONSFILENAME_ID);
1664      zconsFilename->setData(zconsOutput);
1665      
1666 <    info[k].addProperty(zconsFilename);      
1667 <  }
1666 >    theInfo.addProperty(zconsFilename);
1667 >                
1668 >    //setup index, pos and other parameters of z-constraint molecules
1669 >    nZConstraints = globals->getNzConstraints();
1670 >    theInfo.nZconstraints = nZConstraints;
1671 >        
1672 >    zconStamp = globals->getZconStamp();
1673 >    ZConsParaItem tempParaItem;
1674 >
1675 >    ZConsParaData* zconsParaData = new ZConsParaData();
1676 >    zconsParaData->setID(ZCONSPARADATA_ID);
1677 >  
1678 >    for(int i = 0; i < nZConstraints; i++){
1679 >    tempParaItem.havingZPos = zconStamp[i]->haveZpos();
1680 >    tempParaItem.zPos = zconStamp[i]->getZpos();
1681 >    tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1682 >    tempParaItem.kRatio = zconStamp[i]->getKratio();
1683 >
1684 >    zconsParaData->addItem(tempParaItem);
1685 >    }
1686 >
1687 >    //sort the parameters by index of molecules
1688 >    zconsParaData->sortByIndex();
1689 >        
1690 >    //push data into siminfo, therefore, we can retrieve later
1691 >    theInfo.addProperty(zconsParaData);
1692 >      
1693   }

Diff Legend

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