--- trunk/OOPSE/libmdtools/SimSetup.cpp 2003/07/15 18:52:16 616 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2003/08/01 16:18:13 661 @@ -1,6 +1,8 @@ +#include #include #include #include +#include #include "SimSetup.hpp" #include "parse_me.h" @@ -14,22 +16,34 @@ // some defines for ensemble and Forcefield cases -#define NVE_ENS 0 -#define NVT_ENS 1 -#define NPTi_ENS 2 -#define NPTf_ENS 3 -#define NPTim_ENS 4 -#define NPTfm_ENS 5 +#define NVE_ENS 0 +#define NVT_ENS 1 +#define NPTi_ENS 2 +#define NPTf_ENS 3 +#define NPTim_ENS 4 +#define NPTfm_ENS 5 +#define NVEZCONS_ENS 6 +#define NVTZCONS_ENS 7 +#define NPTiZCONS_ENS 8 +#define NPTfZCONS_ENS 9 +#define NPTimZCONS_ENS 10 +#define NPTfmZCONS_ENS 11 - #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(); @@ -40,6 +54,13 @@ SimSetup::~SimSetup(){ delete stamps; delete globals; } + +void SimSetup::setSimInfo( SimInfo* the_info, int theNinfo ) { + info = the_info; + nInfo = theNinfo; + isInfoArray = 1; + } + void SimSetup::parseFile( char* fileName ){ @@ -398,9 +419,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 +436,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 +607,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 +627,15 @@ 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 if( !strcasecmp( ensemble, "NVTZCONS")) ensembleCase = NVTZCONS_ENS; + else if( !strcasecmp( ensemble, "NPTiZCONS") || !strcasecmp( ensemble, "NPTZCONS")) + ensembleCase = NPTiZCONS_ENS; + else if( !strcasecmp( ensemble, "NPTfZCONS")) ensembleCase = NPTfZCONS_ENS; + else if( !strcasecmp( ensemble, "NPTimZCONS")) ensembleCase = NPTimZCONS_ENS; + else if( !strcasecmp( ensemble, "NPTfmZCONS")) ensembleCase = NPTfmZCONS_ENS; + else{ sprintf( painCave.errMsg, "SimSetup Warning. Unrecognized Ensemble -> %s, " @@ -762,6 +793,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 +806,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 +821,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 +842,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; - } else { - info->est = globals->getEST(); + 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 { + theEst= globals->getEST(); } + + info->setEcr( theEcr, theEst ); } } @@ -857,7 +893,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 +1063,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 +1316,26 @@ 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; + ZConstraint >* myNVTZCons = NULL; + ZConstraint >* myNPTiZCons = NULL; + ZConstraint >* myNPTfZCons = NULL; + ZConstraint >* myNPTimZCons = NULL; + ZConstraint >* myNPTfmZCons = NULL; + switch( ensembleCase ){ case NVE_ENS: - new NVE( info, the_ff ); + new NVE( 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 +1351,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 +1386,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 +1421,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 +1456,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 +1489,196 @@ void SimSetup::makeIntegrator( void ){ simError(); } break; + + case NVEZCONS_ENS: + + //setup index of z-constraint molecules, z-constraint sampel time + //and z-constraint force output name. These parameter should be known + //before constructing the z-constraint integrator + setupZConstraint(); + + myNVEZCons = new ZConstraint >( info, the_ff ); + + break; + + + case NVTZCONS_ENS: + + setupZConstraint(); + + myNVTZCons = new ZConstraint >( info, the_ff ); + myNVTZCons->setTargetTemp(globals->getTargetTemp()); + + if (globals->haveTauThermostat()) + myNVTZCons->setTauThermostat(globals->getTauThermostat()); + + else { + sprintf( painCave.errMsg, + "SimSetup error: If you use the NVT\n" + " ensemble, you must set tauThermostat.\n"); + painCave.isFatal = 1; + simError(); + } + break; + + case NPTiZCONS_ENS: + + setupZConstraint(); + + myNPTiZCons = new ZConstraint >( info, the_ff ); + myNPTiZCons->setTargetTemp( globals->getTargetTemp() ); + + if (globals->haveTargetPressure()) + myNPTiZCons->setTargetPressure(globals->getTargetPressure()); + else { + sprintf( painCave.errMsg, + "SimSetup error: If you use a constant pressure\n" + " ensemble, you must set targetPressure in the BASS file.\n"); + painCave.isFatal = 1; + simError(); + } + + if( globals->haveTauThermostat() ) + myNPTiZCons->setTauThermostat( globals->getTauThermostat() ); + else{ + sprintf( painCave.errMsg, + "SimSetup error: If you use an NPT\n" + " ensemble, you must set tauThermostat.\n"); + painCave.isFatal = 1; + simError(); + } + + if( globals->haveTauBarostat() ) + myNPTiZCons->setTauBarostat( globals->getTauBarostat() ); + else{ + sprintf( painCave.errMsg, + "SimSetup error: If you use an NPT\n" + " ensemble, you must set tauBarostat.\n"); + painCave.isFatal = 1; + simError(); + } + + break; + + case NPTfZCONS_ENS: + + setupZConstraint(); + + myNPTfZCons = new ZConstraint >( info, the_ff ); + myNPTfZCons->setTargetTemp( globals->getTargetTemp()); + + if (globals->haveTargetPressure()) + myNPTfZCons->setTargetPressure(globals->getTargetPressure()); + else { + sprintf( painCave.errMsg, + "SimSetup error: If you use a constant pressure\n" + " ensemble, you must set targetPressure in the BASS file.\n"); + painCave.isFatal = 1; + simError(); + } + + if( globals->haveTauThermostat() ) + myNPTfZCons->setTauThermostat( globals->getTauThermostat() ); + else{ + sprintf( painCave.errMsg, + "SimSetup error: If you use an NPT\n" + " ensemble, you must set tauThermostat.\n"); + painCave.isFatal = 1; + simError(); + } + + if( globals->haveTauBarostat() ) + myNPTfZCons->setTauBarostat( globals->getTauBarostat() ); + else{ + sprintf( painCave.errMsg, + "SimSetup error: If you use an NPT\n" + " ensemble, you must set tauBarostat.\n"); + painCave.isFatal = 1; + simError(); + } + + break; + + case NPTimZCONS_ENS: + + setupZConstraint(); + + myNPTimZCons = new ZConstraint >( info, the_ff ); + myNPTimZCons->setTargetTemp( globals->getTargetTemp()); + + if (globals->haveTargetPressure()) + myNPTimZCons->setTargetPressure(globals->getTargetPressure()); + else { + sprintf( painCave.errMsg, + "SimSetup error: If you use a constant pressure\n" + " ensemble, you must set targetPressure in the BASS file.\n"); + painCave.isFatal = 1; + simError(); + } + + if( globals->haveTauThermostat() ) + myNPTimZCons->setTauThermostat( globals->getTauThermostat() ); + else{ + sprintf( painCave.errMsg, + "SimSetup error: If you use an NPT\n" + " ensemble, you must set tauThermostat.\n"); + painCave.isFatal = 1; + simError(); + } + + if( globals->haveTauBarostat() ) + myNPTimZCons->setTauBarostat( globals->getTauBarostat() ); + else{ + sprintf( painCave.errMsg, + "SimSetup error: If you use an NPT\n" + " ensemble, you must set tauBarostat.\n"); + painCave.isFatal = 1; + simError(); + } + + break; + + case NPTfmZCONS_ENS: + + setupZConstraint(); + + myNPTfmZCons = new ZConstraint >( info, the_ff ); + myNPTfmZCons->setTargetTemp( globals->getTargetTemp()); + + if (globals->haveTargetPressure()) + myNPTfmZCons->setTargetPressure(globals->getTargetPressure()); + else { + sprintf( painCave.errMsg, + "SimSetup error: If you use a constant pressure\n" + " ensemble, you must set targetPressure in the BASS file.\n"); + painCave.isFatal = 1; + simError(); + } + + if( globals->haveTauThermostat() ) + myNPTfmZCons->setTauThermostat( globals->getTauThermostat() ); + else{ + sprintf( painCave.errMsg, + "SimSetup error: If you use an NPT\n" + " ensemble, you must set tauThermostat.\n"); + painCave.isFatal = 1; + simError(); + } + + if( globals->haveTauBarostat() ) + myNPTfmZCons->setTauBarostat( globals->getTauBarostat() ); + else{ + sprintf( painCave.errMsg, + "SimSetup error: If you use an NPT\n" + " ensemble, you must set tauBarostat.\n"); + painCave.isFatal = 1; + simError(); + } + break; + + + default: sprintf( painCave.errMsg, "SimSetup Error. Unrecognized ensemble in case statement.\n"); @@ -1479,3 +1714,59 @@ void SimSetup::initFortran( void ){ #endif // is_mpi } + +void SimSetup::setupZConstraint() +{ + 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 the index + 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); + +}