ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/no-template-branch/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 658 by tim, Thu Jul 31 15:35:07 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 20 | Line 22
22   #define NPTf_ENS  3
23   #define NPTim_ENS 4
24   #define NPTfm_ENS 5
25 + #define NVEZCONS_ENS 6
26  
27  
28   #define FF_DUFF 0
29   #define FF_LJ   1
30 + #define FF_EAM  2
31  
32 + using namespace std;
33  
34   SimSetup::SimSetup(){
35 +  
36 +  isInfoArray = 0;
37 +  nInfo = 1;
38 +  
39    stamps = new MakeStamps();
40    globals = new Globals();
41    
42 +  
43   #ifdef IS_MPI
44    strcpy( checkPointMsg, "SimSetup creation successful" );
45    MPIcheckPoint();
# Line 41 | Line 51 | SimSetup::~SimSetup(){
51    delete globals;
52   }
53  
54 + void SimSetup::setSimInfo( SimInfo* the_info, int theNinfo ) {
55 +    info = the_info;
56 +    nInfo = theNinfo;
57 +    isInfoArray = 1;
58 +  }
59 +
60 +
61   void SimSetup::parseFile( char* fileName ){
62  
63   #ifdef IS_MPI
# Line 398 | Line 415 | void SimSetup::initFromBass( void ){
415      have_extra =1;
416  
417      n_cells = (int)temp3 - 1;
418 <    cellx = info->boxLx / temp3;
419 <    celly = info->boxLy / temp3;
420 <    cellz = info->boxLz / temp3;
418 >    cellx = info->boxL[0] / temp3;
419 >    celly = info->boxL[1] / temp3;
420 >    cellz = info->boxL[2] / temp3;
421      n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
422      temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
423      n_per_extra = (int)ceil( temp1 );
# Line 415 | Line 432 | void SimSetup::initFromBass( void ){
432    }
433    else{
434      n_cells = (int)temp3;
435 <    cellx = info->boxLx / temp3;
436 <    celly = info->boxLy / temp3;
437 <    cellz = info->boxLz / temp3;
435 >    cellx = info->boxL[0] / temp3;
436 >    celly = info->boxL[1] / temp3;
437 >    cellz = info->boxL[2] / temp3;
438    }
439  
440    current_mol = 0;
# Line 586 | Line 603 | void SimSetup::gatherInfo( void ){
603  
604    if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
605    else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
606 +  else if( !strcasecmp( force_field, "EAM" )) ffCase = FF_EAM;
607    else{
608      sprintf( painCave.errMsg,
609               "SimSetup Error. Unrecognized force field -> %s\n",
# Line 605 | Line 623 | void SimSetup::gatherInfo( void ){
623    else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
624    else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
625    else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS;
626 +  else if( !strcasecmp( ensemble, "NVEZCONS")) ensembleCase = NVEZCONS_ENS;
627    else{
628      sprintf( painCave.errMsg,
629               "SimSetup Warning. Unrecognized Ensemble -> %s, "
# Line 762 | Line 781 | void SimSetup::finalInfoCheck( void ){
781    MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
782   #endif //is_mpi
783  
784 +  double theEcr, theEst;
785  
786    if (globals->getUseRF() ) {
787      info->useReactionField = 1;
# Line 774 | Line 794 | void SimSetup::finalInfoCheck( void ){
794        painCave.isFatal = 0;
795        simError();
796        double smallest;
797 <      smallest = info->boxLx;
798 <      if (info->boxLy <= smallest) smallest = info->boxLy;
799 <      if (info->boxLz <= smallest) smallest = info->boxLz;
800 <      info->ecr = 0.5 * smallest;
797 >      smallest = info->boxL[0];
798 >      if (info->boxL[1] <= smallest) smallest = info->boxL[1];
799 >      if (info->boxL[2] <= smallest) smallest = info->boxL[2];
800 >      theEcr = 0.5 * smallest;
801      } else {
802 <      info->ecr        = globals->getECR();
802 >      theEcr = globals->getECR();
803      }
804  
805      if( !globals->haveEST() ){
# Line 789 | Line 809 | void SimSetup::finalInfoCheck( void ){
809                 );
810        painCave.isFatal = 0;
811        simError();
812 <      info->est = 0.05 * info->ecr;
812 >      theEst = 0.05 * theEcr;
813      } else {
814 <      info->est        = globals->getEST();
814 >      theEst= globals->getEST();
815      }
816 +
817 +    info->setEcr( theEcr, theEst );
818      
819      if(!globals->haveDielectric() ){
820        sprintf( painCave.errMsg,
# Line 808 | Line 830 | void SimSetup::finalInfoCheck( void ){
830      if (usesDipoles) {
831        
832        if( !globals->haveECR() ){
833 <        sprintf( painCave.errMsg,
834 <                 "SimSetup Warning: using default value of 1/2 the smallest "
835 <                 "box length for the electrostaticCutoffRadius.\n"
836 <                 "I hope you have a very fast processor!\n");
837 <        painCave.isFatal = 0;
838 <        simError();
839 <        double smallest;
840 <        smallest = info->boxLx;
841 <        if (info->boxLy <= smallest) smallest = info->boxLy;
842 <        if (info->boxLz <= smallest) smallest = info->boxLz;
843 <        info->ecr = 0.5 * smallest;
833 >        sprintf( painCave.errMsg,
834 >                 "SimSetup Warning: using default value of 1/2 the smallest "
835 >                 "box length for the electrostaticCutoffRadius.\n"
836 >                 "I hope you have a very fast processor!\n");
837 >        painCave.isFatal = 0;
838 >        simError();
839 >        double smallest;
840 >        smallest = info->boxL[0];
841 >        if (info->boxL[1] <= smallest) smallest = info->boxL[1];
842 >        if (info->boxL[2] <= smallest) smallest = info->boxL[2];
843 >        theEcr = 0.5 * smallest;
844        } else {
845 <        info->ecr        = globals->getECR();
845 >        theEcr = globals->getECR();
846        }
847        
848        if( !globals->haveEST() ){
849 <        sprintf( painCave.errMsg,
850 <                 "SimSetup Warning: using default value of 5%% of the "
851 <                 "electrostaticCutoffRadius for the "
852 <                 "electrostaticSkinThickness\n"
853 <                 );
854 <        painCave.isFatal = 0;
855 <        simError();
856 <        info->est = 0.05 * info->ecr;
849 >        sprintf( painCave.errMsg,
850 >                 "SimSetup Warning: using default value of 0.05 * the "
851 >                 "electrostaticCutoffRadius for the "
852 >                 "electrostaticSkinThickness\n"
853 >                 );
854 >        painCave.isFatal = 0;
855 >        simError();
856 >        theEst = 0.05 * theEcr;
857        } else {
858 <        info->est        = globals->getEST();
858 >        theEst= globals->getEST();
859        }
860 +
861 +      info->setEcr( theEcr, theEst );
862      }
863    }  
864  
# Line 857 | Line 881 | void SimSetup::initSystemCoords( void ){
881   #ifdef IS_MPI
882       }else fileInit = new InitializeFromFile( NULL );
883   #endif
884 <   fileInit->read_xyz( info ); // default velocities on
884 >   fileInit->readInit( info ); // default velocities on
885  
886     delete fileInit;
887   }
# Line 1027 | Line 1051 | void SimSetup::createFF( void ){
1051      the_ff = new LJFF();
1052      break;
1053  
1054 +  case FF_EAM:
1055 +    the_ff = new EAM_FF();
1056 +    break;
1057 +
1058    default:
1059      sprintf( painCave.errMsg,
1060               "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1276 | Line 1304 | void SimSetup::makeIntegrator( void ){
1304  
1305   void SimSetup::makeIntegrator( void ){
1306  
1307 <  NVT*  myNVT = NULL;
1308 <  NPTi* myNPTi = NULL;
1309 <  NPTf* myNPTf = NULL;
1310 <  NPTim* myNPTim = NULL;
1311 <  NPTfm* myNPTfm = NULL;
1312 <
1307 >  NVT<RealIntegrator>*  myNVT = NULL;
1308 >  NPTi<RealIntegrator>* myNPTi = NULL;
1309 >  NPTf<RealIntegrator>* myNPTf = NULL;
1310 >  NPTim<RealIntegrator>* myNPTim = NULL;
1311 >  NPTfm<RealIntegrator>* myNPTfm = NULL;
1312 >  ZConstraint<NVE<RealIntegrator> >* myNVEZCons = NULL;
1313 >      
1314 >  cerr << "setting integrator" <<endl;    
1315 >  
1316    switch( ensembleCase ){
1317  
1318    case NVE_ENS:
1319 <    new NVE( info, the_ff );
1319 >    new NVE<RealIntegrator>( info, the_ff );
1320      break;
1321  
1322    case NVT_ENS:
1323 <    myNVT = new NVT( info, the_ff );
1323 >    myNVT = new NVT<RealIntegrator>( info, the_ff );
1324      myNVT->setTargetTemp(globals->getTargetTemp());
1325  
1326      if (globals->haveTauThermostat())
# Line 1305 | Line 1336 | void SimSetup::makeIntegrator( void ){
1336      break;
1337  
1338    case NPTi_ENS:
1339 <    myNPTi = new NPTi( info, the_ff );
1339 >    myNPTi = new NPTi<RealIntegrator>( info, the_ff );
1340      myNPTi->setTargetTemp( globals->getTargetTemp() );
1341  
1342      if (globals->haveTargetPressure())
# Line 1340 | Line 1371 | void SimSetup::makeIntegrator( void ){
1371      break;
1372  
1373    case NPTf_ENS:
1374 <    myNPTf = new NPTf( info, the_ff );
1374 >    myNPTf = new NPTf<RealIntegrator>( info, the_ff );
1375      myNPTf->setTargetTemp( globals->getTargetTemp());
1376  
1377      if (globals->haveTargetPressure())
# Line 1375 | Line 1406 | void SimSetup::makeIntegrator( void ){
1406      break;
1407      
1408    case NPTim_ENS:
1409 <    myNPTim = new NPTim( info, the_ff );
1409 >    myNPTim = new NPTim<RealIntegrator>( info, the_ff );
1410      myNPTim->setTargetTemp( globals->getTargetTemp());
1411  
1412      if (globals->haveTargetPressure())
# Line 1410 | Line 1441 | void SimSetup::makeIntegrator( void ){
1441      break;
1442  
1443    case NPTfm_ENS:
1444 <    myNPTfm = new NPTfm( info, the_ff );
1444 >    myNPTfm = new NPTfm<RealIntegrator>( info, the_ff );
1445      myNPTfm->setTargetTemp( globals->getTargetTemp());
1446  
1447      if (globals->haveTargetPressure())
# Line 1443 | Line 1474 | void SimSetup::makeIntegrator( void ){
1474        simError();
1475      }
1476      break;
1477 +    
1478 +  case NVEZCONS_ENS:
1479 +    {
1480  
1481 +      if(globals->haveZConsTime()){  
1482 +
1483 +        //add sample time of z-constraint  into SimInfo's property list                    
1484 +        DoubleData* zconsTimeProp = new DoubleData();
1485 +        zconsTimeProp->setID("zconstime");
1486 +        zconsTimeProp->setData(globals->getZConsTime());
1487 +        info->addProperty(zconsTimeProp);
1488 +      }
1489 +      else{
1490 +        sprintf( painCave.errMsg,
1491 +                 "ZConstraint error: If you use an ZConstraint\n"
1492 +                 " , you must set sample time.\n");
1493 +        painCave.isFatal = 1;
1494 +        simError();      
1495 +      }
1496 +      
1497 +      if(globals->haveIndexOfAllZConsMols()){
1498 +
1499 +        //add index of z-constraint molecules into SimInfo's property list
1500 +        vector<int> tempIndex = globals->getIndexOfAllZConsMols();
1501 +        sort(tempIndex.begin(), tempIndex.end());
1502 +        
1503 +        IndexData* zconsIndex = new IndexData();
1504 +        zconsIndex->setID("zconsindex");
1505 +        zconsIndex->setIndexData(tempIndex);
1506 +        info->addProperty(zconsIndex);
1507 +      }
1508 +      else{
1509 +        sprintf( painCave.errMsg,
1510 +                 "SimSetup error: If you use an ZConstraint\n"
1511 +                 " , you must set index of z-constraint molecules.\n");
1512 +        painCave.isFatal = 1;
1513 +        simError();    
1514 +      
1515 +      }
1516 +
1517 +      //Determine the name of ouput file and add it into SimInfo's property list
1518 +      //Be careful, do not use inFileName, since it is a pointer which
1519 +      //point to a string at master node, and slave nodes do not contain that string
1520 +    
1521 +      string zconsOutput(info->finalName);
1522 +            
1523 +      zconsOutput = zconsOutput.substr(0, zconsOutput.rfind(".")) + ".fz";
1524 +                
1525 +      StringData* zconsFilename = new StringData();
1526 +      zconsFilename->setID("zconsfilename");
1527 +      zconsFilename->setData(zconsOutput);
1528 +
1529 +      info->addProperty(zconsFilename);      
1530 +      
1531 +      myNVEZCons = new ZConstraint<NVE<RealIntegrator> >( info, the_ff );
1532 +        
1533 +    break;
1534 +    }
1535 +    
1536    default:
1537      sprintf( painCave.errMsg,
1538               "SimSetup Error. Unrecognized ensemble in case statement.\n");

Diff Legend

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