--- trunk/OOPSE/libmdtools/SimSetup.cpp 2003/07/15 18:52:16 616 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2003/07/31 15:35:07 658 @@ -1,6 +1,8 @@ +#include #include #include #include +#include #include "SimSetup.hpp" #include "parse_me.h" @@ -20,16 +22,24 @@ #define NPTf_ENS 3 #define NPTim_ENS 4 #define NPTfm_ENS 5 +#define NVEZCONS_ENS 6 #define FF_DUFF 0 #define FF_LJ 1 +#define FF_EAM 2 +using namespace std; SimSetup::SimSetup(){ + + isInfoArray = 0; + nInfo = 1; + stamps = new MakeStamps(); globals = new Globals(); + #ifdef IS_MPI strcpy( checkPointMsg, "SimSetup creation successful" ); MPIcheckPoint(); @@ -41,6 +51,13 @@ SimSetup::~SimSetup(){ delete globals; } +void SimSetup::setSimInfo( SimInfo* the_info, int theNinfo ) { + info = the_info; + nInfo = theNinfo; + isInfoArray = 1; + } + + void SimSetup::parseFile( char* fileName ){ #ifdef IS_MPI @@ -398,9 +415,9 @@ void SimSetup::initFromBass( void ){ have_extra =1; n_cells = (int)temp3 - 1; - cellx = info->boxLx / temp3; - celly = info->boxLy / temp3; - cellz = info->boxLz / temp3; + cellx = info->boxL[0] / temp3; + celly = info->boxL[1] / temp3; + cellz = info->boxL[2] / temp3; n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells ); temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) ); n_per_extra = (int)ceil( temp1 ); @@ -415,9 +432,9 @@ void SimSetup::initFromBass( void ){ } else{ n_cells = (int)temp3; - cellx = info->boxLx / temp3; - celly = info->boxLy / temp3; - cellz = info->boxLz / temp3; + cellx = info->boxL[0] / temp3; + celly = info->boxL[1] / temp3; + cellz = info->boxL[2] / temp3; } current_mol = 0; @@ -586,6 +603,7 @@ void SimSetup::gatherInfo( void ){ if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF; else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ; + else if( !strcasecmp( force_field, "EAM" )) ffCase = FF_EAM; else{ sprintf( painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n", @@ -605,6 +623,7 @@ void SimSetup::gatherInfo( void ){ else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS; else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS; else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS; + else if( !strcasecmp( ensemble, "NVEZCONS")) ensembleCase = NVEZCONS_ENS; else{ sprintf( painCave.errMsg, "SimSetup Warning. Unrecognized Ensemble -> %s, " @@ -762,6 +781,7 @@ void SimSetup::finalInfoCheck( void ){ MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD ); #endif //is_mpi + double theEcr, theEst; if (globals->getUseRF() ) { info->useReactionField = 1; @@ -774,12 +794,12 @@ void SimSetup::finalInfoCheck( void ){ painCave.isFatal = 0; simError(); double smallest; - smallest = info->boxLx; - if (info->boxLy <= smallest) smallest = info->boxLy; - if (info->boxLz <= smallest) smallest = info->boxLz; - info->ecr = 0.5 * smallest; + smallest = info->boxL[0]; + if (info->boxL[1] <= smallest) smallest = info->boxL[1]; + if (info->boxL[2] <= smallest) smallest = info->boxL[2]; + theEcr = 0.5 * smallest; } else { - info->ecr = globals->getECR(); + theEcr = globals->getECR(); } if( !globals->haveEST() ){ @@ -789,10 +809,12 @@ void SimSetup::finalInfoCheck( void ){ ); painCave.isFatal = 0; simError(); - info->est = 0.05 * info->ecr; + theEst = 0.05 * theEcr; } else { - info->est = globals->getEST(); + theEst= globals->getEST(); } + + info->setEcr( theEcr, theEst ); if(!globals->haveDielectric() ){ sprintf( painCave.errMsg, @@ -808,33 +830,35 @@ void SimSetup::finalInfoCheck( void ){ if (usesDipoles) { if( !globals->haveECR() ){ - sprintf( painCave.errMsg, - "SimSetup Warning: using default value of 1/2 the smallest " - "box length for the electrostaticCutoffRadius.\n" - "I hope you have a very fast processor!\n"); - painCave.isFatal = 0; - simError(); - double smallest; - smallest = info->boxLx; - if (info->boxLy <= smallest) smallest = info->boxLy; - if (info->boxLz <= smallest) smallest = info->boxLz; - info->ecr = 0.5 * smallest; + sprintf( painCave.errMsg, + "SimSetup Warning: using default value of 1/2 the smallest " + "box length for the electrostaticCutoffRadius.\n" + "I hope you have a very fast processor!\n"); + painCave.isFatal = 0; + simError(); + double smallest; + smallest = info->boxL[0]; + if (info->boxL[1] <= smallest) smallest = info->boxL[1]; + if (info->boxL[2] <= smallest) smallest = info->boxL[2]; + theEcr = 0.5 * smallest; } else { - info->ecr = globals->getECR(); + theEcr = globals->getECR(); } if( !globals->haveEST() ){ - sprintf( painCave.errMsg, - "SimSetup Warning: using default value of 5%% of the " - "electrostaticCutoffRadius for the " - "electrostaticSkinThickness\n" - ); - painCave.isFatal = 0; - simError(); - info->est = 0.05 * info->ecr; + sprintf( painCave.errMsg, + "SimSetup Warning: using default value of 0.05 * the " + "electrostaticCutoffRadius for the " + "electrostaticSkinThickness\n" + ); + painCave.isFatal = 0; + simError(); + theEst = 0.05 * theEcr; } else { - info->est = globals->getEST(); + theEst= globals->getEST(); } + + info->setEcr( theEcr, theEst ); } } @@ -857,7 +881,7 @@ void SimSetup::initSystemCoords( void ){ #ifdef IS_MPI }else fileInit = new InitializeFromFile( NULL ); #endif - fileInit->read_xyz( info ); // default velocities on + fileInit->readInit( info ); // default velocities on delete fileInit; } @@ -1027,6 +1051,10 @@ void SimSetup::createFF( void ){ the_ff = new LJFF(); break; + case FF_EAM: + the_ff = new EAM_FF(); + break; + default: sprintf( painCave.errMsg, "SimSetup Error. Unrecognized force field in case statement.\n"); @@ -1276,20 +1304,23 @@ void SimSetup::makeIntegrator( void ){ void SimSetup::makeIntegrator( void ){ - NVT* myNVT = NULL; - NPTi* myNPTi = NULL; - NPTf* myNPTf = NULL; - NPTim* myNPTim = NULL; - NPTfm* myNPTfm = NULL; - + NVT* myNVT = NULL; + NPTi* myNPTi = NULL; + NPTf* myNPTf = NULL; + NPTim* myNPTim = NULL; + NPTfm* myNPTfm = NULL; + ZConstraint >* myNVEZCons = NULL; + + cerr << "setting integrator" <( info, the_ff ); break; case NVT_ENS: - myNVT = new NVT( info, the_ff ); + myNVT = new NVT( info, the_ff ); myNVT->setTargetTemp(globals->getTargetTemp()); if (globals->haveTauThermostat()) @@ -1305,7 +1336,7 @@ void SimSetup::makeIntegrator( void ){ break; case NPTi_ENS: - myNPTi = new NPTi( info, the_ff ); + myNPTi = new NPTi( info, the_ff ); myNPTi->setTargetTemp( globals->getTargetTemp() ); if (globals->haveTargetPressure()) @@ -1340,7 +1371,7 @@ void SimSetup::makeIntegrator( void ){ break; case NPTf_ENS: - myNPTf = new NPTf( info, the_ff ); + myNPTf = new NPTf( info, the_ff ); myNPTf->setTargetTemp( globals->getTargetTemp()); if (globals->haveTargetPressure()) @@ -1375,7 +1406,7 @@ void SimSetup::makeIntegrator( void ){ break; case NPTim_ENS: - myNPTim = new NPTim( info, the_ff ); + myNPTim = new NPTim( info, the_ff ); myNPTim->setTargetTemp( globals->getTargetTemp()); if (globals->haveTargetPressure()) @@ -1410,7 +1441,7 @@ void SimSetup::makeIntegrator( void ){ break; case NPTfm_ENS: - myNPTfm = new NPTfm( info, the_ff ); + myNPTfm = new NPTfm( info, the_ff ); myNPTfm->setTargetTemp( globals->getTargetTemp()); if (globals->haveTargetPressure()) @@ -1443,7 +1474,65 @@ void SimSetup::makeIntegrator( void ){ simError(); } break; + + case NVEZCONS_ENS: + { + if(globals->haveZConsTime()){ + + //add sample time of z-constraint into SimInfo's property list + DoubleData* zconsTimeProp = new DoubleData(); + zconsTimeProp->setID("zconstime"); + zconsTimeProp->setData(globals->getZConsTime()); + info->addProperty(zconsTimeProp); + } + else{ + sprintf( painCave.errMsg, + "ZConstraint error: If you use an ZConstraint\n" + " , you must set sample time.\n"); + painCave.isFatal = 1; + simError(); + } + + if(globals->haveIndexOfAllZConsMols()){ + + //add index of z-constraint molecules into SimInfo's property list + vector tempIndex = globals->getIndexOfAllZConsMols(); + sort(tempIndex.begin(), tempIndex.end()); + + IndexData* zconsIndex = new IndexData(); + zconsIndex->setID("zconsindex"); + zconsIndex->setIndexData(tempIndex); + info->addProperty(zconsIndex); + } + else{ + sprintf( painCave.errMsg, + "SimSetup error: If you use an ZConstraint\n" + " , you must set index of z-constraint molecules.\n"); + painCave.isFatal = 1; + simError(); + + } + + //Determine the name of ouput file and add it into SimInfo's property list + //Be careful, do not use inFileName, since it is a pointer which + //point to a string at master node, and slave nodes do not contain that string + + string zconsOutput(info->finalName); + + zconsOutput = zconsOutput.substr(0, zconsOutput.rfind(".")) + ".fz"; + + StringData* zconsFilename = new StringData(); + zconsFilename->setID("zconsfilename"); + zconsFilename->setData(zconsOutput); + + info->addProperty(zconsFilename); + + myNVEZCons = new ZConstraint >( info, the_ff ); + + break; + } + default: sprintf( painCave.errMsg, "SimSetup Error. Unrecognized ensemble in case statement.\n");