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 616 by mmeineke, Tue Jul 15 18:52:16 2003 UTC vs.
Revision 660 by tim, Thu Jul 31 19:59:34 2003 UTC

# Line 1 | Line 1
1 + #include <algorithm>
2   #include <cstdlib>
3   #include <iostream>
4   #include <cmath>
5 + #include <string>
6  
7   #include "SimSetup.hpp"
8   #include "parse_me.h"
# Line 14 | Line 16
16  
17   // some defines for ensemble and Forcefield  cases
18  
19 < #define NVE_ENS   0
20 < #define NVT_ENS   1
21 < #define NPTi_ENS  2
22 < #define NPTf_ENS  3
23 < #define NPTim_ENS 4
24 < #define NPTfm_ENS 5
19 > #define NVE_ENS        0
20 > #define NVT_ENS        1
21 > #define NPTi_ENS       2
22 > #define NPTf_ENS       3
23 > #define NPTim_ENS      4
24 > #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
31  
24 –
32   #define FF_DUFF 0
33   #define FF_LJ   1
34 + #define FF_EAM  2
35  
36 + using namespace std;
37  
38   SimSetup::SimSetup(){
39 +  
40 +  isInfoArray = 0;
41 +  nInfo = 1;
42 +  
43    stamps = new MakeStamps();
44    globals = new Globals();
45    
46 +  
47   #ifdef IS_MPI
48    strcpy( checkPointMsg, "SimSetup creation successful" );
49    MPIcheckPoint();
# Line 41 | Line 55 | SimSetup::~SimSetup(){
55    delete globals;
56   }
57  
58 + void SimSetup::setSimInfo( SimInfo* the_info, int theNinfo ) {
59 +    info = the_info;
60 +    nInfo = theNinfo;
61 +    isInfoArray = 1;
62 +  }
63 +
64 +
65   void SimSetup::parseFile( char* fileName ){
66  
67   #ifdef IS_MPI
# Line 398 | Line 419 | void SimSetup::initFromBass( void ){
419      have_extra =1;
420  
421      n_cells = (int)temp3 - 1;
422 <    cellx = info->boxLx / temp3;
423 <    celly = info->boxLy / temp3;
424 <    cellz = info->boxLz / temp3;
422 >    cellx = info->boxL[0] / temp3;
423 >    celly = info->boxL[1] / temp3;
424 >    cellz = info->boxL[2] / temp3;
425      n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
426      temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
427      n_per_extra = (int)ceil( temp1 );
# Line 415 | Line 436 | void SimSetup::initFromBass( void ){
436    }
437    else{
438      n_cells = (int)temp3;
439 <    cellx = info->boxLx / temp3;
440 <    celly = info->boxLy / temp3;
441 <    cellz = info->boxLz / temp3;
439 >    cellx = info->boxL[0] / temp3;
440 >    celly = info->boxL[1] / temp3;
441 >    cellz = info->boxL[2] / temp3;
442    }
443  
444    current_mol = 0;
# Line 586 | Line 607 | void SimSetup::gatherInfo( void ){
607  
608    if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
609    else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
610 +  else if( !strcasecmp( force_field, "EAM" )) ffCase = FF_EAM;
611    else{
612      sprintf( painCave.errMsg,
613               "SimSetup Error. Unrecognized force field -> %s\n",
# Line 605 | Line 627 | void SimSetup::gatherInfo( void ){
627    else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
628    else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
629    else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS;
630 +
631 +  else if( !strcasecmp( ensemble, "NVEZCONS")) ensembleCase = NVEZCONS_ENS;
632 +  else if( !strcasecmp( ensemble, "NVTZCONS"))  ensembleCase = NVTZCONS_ENS;
633 +  else if( !strcasecmp( ensemble, "NPTiZCONS") || !strcasecmp( ensemble, "NPT"))
634 +    ensembleCase = NPTiZCONS_ENS;
635 +  else if( !strcasecmp( ensemble, "NPTfCONS"))  ensembleCase = NPTfZCONS_ENS;
636 +  else if( !strcasecmp( ensemble, "NPTimZCONS"))  ensembleCase = NPTimZCONS_ENS;
637 +  else if( !strcasecmp( ensemble, "NPTfmCONS"))  ensembleCase = NPTfmZCONS_ENS;
638 +  
639    else{
640      sprintf( painCave.errMsg,
641               "SimSetup Warning. Unrecognized Ensemble -> %s, "
# Line 762 | Line 793 | void SimSetup::finalInfoCheck( void ){
793    MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
794   #endif //is_mpi
795  
796 +  double theEcr, theEst;
797  
798    if (globals->getUseRF() ) {
799      info->useReactionField = 1;
# Line 774 | Line 806 | void SimSetup::finalInfoCheck( void ){
806        painCave.isFatal = 0;
807        simError();
808        double smallest;
809 <      smallest = info->boxLx;
810 <      if (info->boxLy <= smallest) smallest = info->boxLy;
811 <      if (info->boxLz <= smallest) smallest = info->boxLz;
812 <      info->ecr = 0.5 * smallest;
809 >      smallest = info->boxL[0];
810 >      if (info->boxL[1] <= smallest) smallest = info->boxL[1];
811 >      if (info->boxL[2] <= smallest) smallest = info->boxL[2];
812 >      theEcr = 0.5 * smallest;
813      } else {
814 <      info->ecr        = globals->getECR();
814 >      theEcr = globals->getECR();
815      }
816  
817      if( !globals->haveEST() ){
# Line 789 | Line 821 | void SimSetup::finalInfoCheck( void ){
821                 );
822        painCave.isFatal = 0;
823        simError();
824 <      info->est = 0.05 * info->ecr;
824 >      theEst = 0.05 * theEcr;
825      } else {
826 <      info->est        = globals->getEST();
826 >      theEst= globals->getEST();
827      }
828 +
829 +    info->setEcr( theEcr, theEst );
830      
831      if(!globals->haveDielectric() ){
832        sprintf( painCave.errMsg,
# Line 808 | Line 842 | void SimSetup::finalInfoCheck( void ){
842      if (usesDipoles) {
843        
844        if( !globals->haveECR() ){
845 <        sprintf( painCave.errMsg,
846 <                 "SimSetup Warning: using default value of 1/2 the smallest "
847 <                 "box length for the electrostaticCutoffRadius.\n"
848 <                 "I hope you have a very fast processor!\n");
849 <        painCave.isFatal = 0;
850 <        simError();
851 <        double smallest;
852 <        smallest = info->boxLx;
853 <        if (info->boxLy <= smallest) smallest = info->boxLy;
854 <        if (info->boxLz <= smallest) smallest = info->boxLz;
855 <        info->ecr = 0.5 * smallest;
845 >        sprintf( painCave.errMsg,
846 >                 "SimSetup Warning: using default value of 1/2 the smallest "
847 >                 "box length for the electrostaticCutoffRadius.\n"
848 >                 "I hope you have a very fast processor!\n");
849 >        painCave.isFatal = 0;
850 >        simError();
851 >        double smallest;
852 >        smallest = info->boxL[0];
853 >        if (info->boxL[1] <= smallest) smallest = info->boxL[1];
854 >        if (info->boxL[2] <= smallest) smallest = info->boxL[2];
855 >        theEcr = 0.5 * smallest;
856        } else {
857 <        info->ecr        = globals->getECR();
857 >        theEcr = globals->getECR();
858        }
859        
860        if( !globals->haveEST() ){
861 <        sprintf( painCave.errMsg,
862 <                 "SimSetup Warning: using default value of 5%% of the "
863 <                 "electrostaticCutoffRadius for the "
864 <                 "electrostaticSkinThickness\n"
865 <                 );
866 <        painCave.isFatal = 0;
867 <        simError();
868 <        info->est = 0.05 * info->ecr;
861 >        sprintf( painCave.errMsg,
862 >                 "SimSetup Warning: using default value of 0.05 * the "
863 >                 "electrostaticCutoffRadius for the "
864 >                 "electrostaticSkinThickness\n"
865 >                 );
866 >        painCave.isFatal = 0;
867 >        simError();
868 >        theEst = 0.05 * theEcr;
869        } else {
870 <        info->est        = globals->getEST();
870 >        theEst= globals->getEST();
871        }
872 +
873 +      info->setEcr( theEcr, theEst );
874      }
875    }  
876  
# Line 857 | Line 893 | void SimSetup::initSystemCoords( void ){
893   #ifdef IS_MPI
894       }else fileInit = new InitializeFromFile( NULL );
895   #endif
896 <   fileInit->read_xyz( info ); // default velocities on
896 >   fileInit->readInit( info ); // default velocities on
897  
898     delete fileInit;
899   }
# Line 1027 | Line 1063 | void SimSetup::createFF( void ){
1063      the_ff = new LJFF();
1064      break;
1065  
1066 +  case FF_EAM:
1067 +    the_ff = new EAM_FF();
1068 +    break;
1069 +
1070    default:
1071      sprintf( painCave.errMsg,
1072               "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1276 | Line 1316 | void SimSetup::makeIntegrator( void ){
1316  
1317   void SimSetup::makeIntegrator( void ){
1318  
1319 <  NVT*  myNVT = NULL;
1320 <  NPTi* myNPTi = NULL;
1321 <  NPTf* myNPTf = NULL;
1322 <  NPTim* myNPTim = NULL;
1323 <  NPTfm* myNPTfm = NULL;
1324 <
1319 >  NVT<RealIntegrator>*  myNVT = NULL;
1320 >  NPTi<RealIntegrator>* myNPTi = NULL;
1321 >  NPTf<RealIntegrator>* myNPTf = NULL;
1322 >  NPTim<RealIntegrator>* myNPTim = NULL;
1323 >  NPTfm<RealIntegrator>* myNPTfm = NULL;
1324 >  ZConstraint<NVE<RealIntegrator> >* myNVEZCons = NULL;
1325 >  ZConstraint<NVT<RealIntegrator> >* myNVTZCons = NULL;
1326 >  ZConstraint<NPTi<RealIntegrator> >* myNPTiZCons = NULL;
1327 >  ZConstraint<NPTf<RealIntegrator> >* myNPTfZCons = NULL;
1328 >  ZConstraint<NPTim<RealIntegrator> >* myNPTimZCons = NULL;
1329 >  ZConstraint<NPTfm<RealIntegrator> >* myNPTfmZCons = NULL;
1330 >      
1331 >  cerr << "setting integrator" <<endl;    
1332 >  
1333    switch( ensembleCase ){
1334  
1335    case NVE_ENS:
1336 <    new NVE( info, the_ff );
1336 >    new NVE<RealIntegrator>( info, the_ff );
1337      break;
1338  
1339    case NVT_ENS:
1340 <    myNVT = new NVT( info, the_ff );
1340 >    myNVT = new NVT<RealIntegrator>( info, the_ff );
1341      myNVT->setTargetTemp(globals->getTargetTemp());
1342  
1343      if (globals->haveTauThermostat())
# Line 1305 | Line 1353 | void SimSetup::makeIntegrator( void ){
1353      break;
1354  
1355    case NPTi_ENS:
1356 <    myNPTi = new NPTi( info, the_ff );
1356 >    myNPTi = new NPTi<RealIntegrator>( info, the_ff );
1357      myNPTi->setTargetTemp( globals->getTargetTemp() );
1358  
1359      if (globals->haveTargetPressure())
# Line 1340 | Line 1388 | void SimSetup::makeIntegrator( void ){
1388      break;
1389  
1390    case NPTf_ENS:
1391 <    myNPTf = new NPTf( info, the_ff );
1391 >    myNPTf = new NPTf<RealIntegrator>( info, the_ff );
1392      myNPTf->setTargetTemp( globals->getTargetTemp());
1393  
1394      if (globals->haveTargetPressure())
# Line 1375 | Line 1423 | void SimSetup::makeIntegrator( void ){
1423      break;
1424      
1425    case NPTim_ENS:
1426 <    myNPTim = new NPTim( info, the_ff );
1426 >    myNPTim = new NPTim<RealIntegrator>( info, the_ff );
1427      myNPTim->setTargetTemp( globals->getTargetTemp());
1428  
1429      if (globals->haveTargetPressure())
# Line 1410 | Line 1458 | void SimSetup::makeIntegrator( void ){
1458      break;
1459  
1460    case NPTfm_ENS:
1461 <    myNPTfm = new NPTfm( info, the_ff );
1461 >    myNPTfm = new NPTfm<RealIntegrator>( info, the_ff );
1462      myNPTfm->setTargetTemp( globals->getTargetTemp());
1463  
1464      if (globals->haveTargetPressure())
# Line 1443 | Line 1491 | void SimSetup::makeIntegrator( void ){
1491        simError();
1492      }
1493      break;
1494 +    
1495 +  case NVEZCONS_ENS:
1496  
1497 +
1498 +    //setup index of z-constraint molecules, z-constraint sampel time
1499 +    //and z-constraint force output name. These parameter should be known
1500 +    //before constructing the z-constraint integrator
1501 +    setupZConstraint();
1502 +      
1503 +    myNVEZCons = new ZConstraint<NVE<RealIntegrator> >( info, the_ff );
1504 +        
1505 +    break;
1506 +    
1507 +    
1508 +  case NVTZCONS_ENS:
1509 +  
1510 +    setupZConstraint();
1511 +    
1512 +    myNVTZCons = new ZConstraint<NVT<RealIntegrator> >( info, the_ff );
1513 +    myNVTZCons->setTargetTemp(globals->getTargetTemp());
1514 +
1515 +    if (globals->haveTauThermostat())
1516 +      myNVTZCons->setTauThermostat(globals->getTauThermostat());
1517 +
1518 +    else {
1519 +      sprintf( painCave.errMsg,
1520 +               "SimSetup error: If you use the NVT\n"
1521 +               "    ensemble, you must set tauThermostat.\n");
1522 +      painCave.isFatal = 1;
1523 +      simError();
1524 +    }    
1525 +    break;    
1526 +    
1527 +  case NPTiZCONS_ENS:
1528 +  
1529 +    setupZConstraint();
1530 +    
1531 +    myNPTiZCons = new ZConstraint<NPTi<RealIntegrator> >( info, the_ff );
1532 +    myNPTiZCons->setTargetTemp( globals->getTargetTemp() );
1533 +
1534 +    if (globals->haveTargetPressure())
1535 +      myNPTiZCons->setTargetPressure(globals->getTargetPressure());
1536 +    else {
1537 +      sprintf( painCave.errMsg,
1538 +               "SimSetup error: If you use a constant pressure\n"
1539 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1540 +      painCave.isFatal = 1;
1541 +      simError();
1542 +    }
1543 +    
1544 +    if( globals->haveTauThermostat() )
1545 +      myNPTiZCons->setTauThermostat( globals->getTauThermostat() );
1546 +    else{
1547 +      sprintf( painCave.errMsg,
1548 +               "SimSetup error: If you use an NPT\n"
1549 +               "    ensemble, you must set tauThermostat.\n");
1550 +      painCave.isFatal = 1;
1551 +      simError();
1552 +    }
1553 +
1554 +    if( globals->haveTauBarostat() )
1555 +      myNPTiZCons->setTauBarostat( globals->getTauBarostat() );
1556 +    else{
1557 +      sprintf( painCave.errMsg,
1558 +               "SimSetup error: If you use an NPT\n"
1559 +               "    ensemble, you must set tauBarostat.\n");
1560 +      painCave.isFatal = 1;
1561 +      simError();
1562 +    }  
1563 +    
1564 +    break;
1565 +    
1566 +  case NPTfZCONS_ENS:
1567 +  
1568 +    setupZConstraint();
1569 +  
1570 +    myNPTfZCons = new ZConstraint<NPTf<RealIntegrator> >( info, the_ff );
1571 +    myNPTfZCons->setTargetTemp( globals->getTargetTemp());
1572 +
1573 +    if (globals->haveTargetPressure())
1574 +      myNPTfZCons->setTargetPressure(globals->getTargetPressure());
1575 +    else {
1576 +      sprintf( painCave.errMsg,
1577 +               "SimSetup error: If you use a constant pressure\n"
1578 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1579 +      painCave.isFatal = 1;
1580 +      simError();
1581 +    }    
1582 +
1583 +    if( globals->haveTauThermostat() )
1584 +      myNPTfZCons->setTauThermostat( globals->getTauThermostat() );
1585 +    else{
1586 +      sprintf( painCave.errMsg,
1587 +               "SimSetup error: If you use an NPT\n"
1588 +               "    ensemble, you must set tauThermostat.\n");
1589 +      painCave.isFatal = 1;
1590 +      simError();
1591 +    }
1592 +
1593 +    if( globals->haveTauBarostat() )
1594 +      myNPTfZCons->setTauBarostat( globals->getTauBarostat() );
1595 +    else{
1596 +      sprintf( painCave.errMsg,
1597 +               "SimSetup error: If you use an NPT\n"
1598 +               "    ensemble, you must set tauBarostat.\n");
1599 +      painCave.isFatal = 1;
1600 +      simError();
1601 +    }  
1602 +    
1603 +    break;  
1604 +      
1605 +  case NPTimZCONS_ENS:
1606 +  
1607 +    setupZConstraint();
1608 +  
1609 +    myNPTimZCons = new ZConstraint<NPTim<RealIntegrator> >( info, the_ff );
1610 +    myNPTimZCons->setTargetTemp( globals->getTargetTemp());
1611 +
1612 +    if (globals->haveTargetPressure())
1613 +      myNPTimZCons->setTargetPressure(globals->getTargetPressure());
1614 +    else {
1615 +      sprintf( painCave.errMsg,
1616 +               "SimSetup error: If you use a constant pressure\n"
1617 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1618 +      painCave.isFatal = 1;
1619 +      simError();
1620 +    }
1621 +    
1622 +    if( globals->haveTauThermostat() )
1623 +      myNPTimZCons->setTauThermostat( globals->getTauThermostat() );
1624 +    else{
1625 +      sprintf( painCave.errMsg,
1626 +               "SimSetup error: If you use an NPT\n"
1627 +               "    ensemble, you must set tauThermostat.\n");
1628 +      painCave.isFatal = 1;
1629 +      simError();
1630 +    }
1631 +
1632 +    if( globals->haveTauBarostat() )
1633 +      myNPTimZCons->setTauBarostat( globals->getTauBarostat() );
1634 +    else{
1635 +      sprintf( painCave.errMsg,
1636 +               "SimSetup error: If you use an NPT\n"
1637 +               "    ensemble, you must set tauBarostat.\n");
1638 +      painCave.isFatal = 1;
1639 +      simError();
1640 +    }  
1641 +    
1642 +    break;
1643 +    
1644 +  case NPTfmZCONS_ENS:
1645 +  
1646 +    setupZConstraint();
1647 +    
1648 +    myNPTfmZCons = new ZConstraint<NPTfm<RealIntegrator> >( info, the_ff );
1649 +    myNPTfmZCons->setTargetTemp( globals->getTargetTemp());
1650 +
1651 +    if (globals->haveTargetPressure())
1652 +      myNPTfmZCons->setTargetPressure(globals->getTargetPressure());
1653 +    else {
1654 +      sprintf( painCave.errMsg,
1655 +               "SimSetup error: If you use a constant pressure\n"
1656 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1657 +      painCave.isFatal = 1;
1658 +      simError();
1659 +    }
1660 +    
1661 +    if( globals->haveTauThermostat() )
1662 +      myNPTfmZCons->setTauThermostat( globals->getTauThermostat() );
1663 +    else{
1664 +      sprintf( painCave.errMsg,
1665 +               "SimSetup error: If you use an NPT\n"
1666 +               "    ensemble, you must set tauThermostat.\n");
1667 +      painCave.isFatal = 1;
1668 +      simError();
1669 +    }
1670 +
1671 +    if( globals->haveTauBarostat() )
1672 +      myNPTfmZCons->setTauBarostat( globals->getTauBarostat() );
1673 +    else{
1674 +      sprintf( painCave.errMsg,
1675 +               "SimSetup error: If you use an NPT\n"
1676 +               "    ensemble, you must set tauBarostat.\n");
1677 +      painCave.isFatal = 1;
1678 +      simError();
1679 +    }    
1680 +    break;      
1681 +      
1682 +  
1683 +    
1684    default:
1685      sprintf( painCave.errMsg,
1686               "SimSetup Error. Unrecognized ensemble in case statement.\n");
# Line 1479 | Line 1716 | void SimSetup::initFortran( void ){
1716   #endif // is_mpi
1717  
1718   }
1719 +
1720 + void SimSetup::setupZConstraint()
1721 + {
1722 +  if(globals->haveZConsTime()){  
1723 +
1724 +  //add sample time of z-constraint  into SimInfo's property list                    
1725 +  DoubleData* zconsTimeProp = new DoubleData();
1726 +  zconsTimeProp->setID("zconstime");
1727 +  zconsTimeProp->setData(globals->getZConsTime());
1728 +  info->addProperty(zconsTimeProp);
1729 +  }
1730 +  else{
1731 +    sprintf( painCave.errMsg,
1732 +             "ZConstraint error: If you use an ZConstraint\n"
1733 +             " , you must set sample time.\n");
1734 +    painCave.isFatal = 1;
1735 +    simError();      
1736 +  }
1737 +      
1738 +  if(globals->haveIndexOfAllZConsMols()){
1739 +
1740 +        //add index of z-constraint molecules into SimInfo's property list
1741 +        vector<int> tempIndex = globals->getIndexOfAllZConsMols();
1742 +        
1743 +        //sort the index
1744 +        sort(tempIndex.begin(), tempIndex.end());
1745 +        
1746 +        IndexData* zconsIndex = new IndexData();
1747 +        zconsIndex->setID("zconsindex");
1748 +        zconsIndex->setIndexData(tempIndex);
1749 +        info->addProperty(zconsIndex);
1750 +  }
1751 +  else{
1752 +    sprintf( painCave.errMsg,
1753 +             "SimSetup error: If you use an ZConstraint\n"
1754 +             " , you must set index of z-constraint molecules.\n");
1755 +    painCave.isFatal = 1;
1756 +    simError();    
1757 +      
1758 +  }
1759 +
1760 +  //Determine the name of ouput file and add it into SimInfo's property list
1761 +  //Be careful, do not use inFileName, since it is a pointer which
1762 +  //point to a string at master node, and slave nodes do not contain that string
1763 +    
1764 +  string zconsOutput(info->finalName);
1765 +            
1766 +  zconsOutput = zconsOutput.substr(0, zconsOutput.rfind(".")) + ".fz";
1767 +                
1768 +  StringData* zconsFilename = new StringData();
1769 +  zconsFilename->setID("zconsfilename");
1770 +  zconsFilename->setData(zconsOutput);
1771 +
1772 +  info->addProperty(zconsFilename);      
1773 +
1774 + }

Diff Legend

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