--- trunk/OOPSE/libmdtools/SimSetup.cpp 2003/07/15 17:57:04 614 +++ 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 @@ -78,8 +95,6 @@ void SimSetup::createSim( void ){ void SimSetup::createSim( void ){ - MakeStamps *the_stamps; - Globals* the_globals; int i, j, k, globalAtomIndex; // gather all of the information from the Bass file @@ -90,26 +105,10 @@ void SimSetup::createSim( void ){ sysObjectsCreation(); - - - // initialize the arrays - - - - makeMolecules(); - info->identArray = new int[info->n_atoms]; - for(i=0; in_atoms; i++){ - info->identArray[i] = the_atoms[i]->getIdent(); - } - - // check on the post processing info finalInfoCheck(); - - - // initialize the system coordinates initSystemCoords(); @@ -119,229 +118,27 @@ void SimSetup::createSim( void ){ makeOutNames(); - - - - - - - - // make the integrator + makeIntegrator(); - NVT* myNVT = NULL; - NPTi* myNPTi = NULL; - NPTf* myNPTf = NULL; - NPTim* myNPTim = NULL; - NPTfm* myNPTfm = NULL; - - switch( ensembleCase ){ - - case NVE_ENS: - new NVE( info, the_ff ); - break; - - case NVT_ENS: - myNVT = new NVT( info, the_ff ); - myNVT->setTargetTemp(the_globals->getTargetTemp()); - - if (the_globals->haveTauThermostat()) - myNVT->setTauThermostat(the_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 NPTi_ENS: - myNPTi = new NPTi( info, the_ff ); - myNPTi->setTargetTemp( the_globals->getTargetTemp() ); - - if (the_globals->haveTargetPressure()) - myNPTi->setTargetPressure(the_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( the_globals->haveTauThermostat() ) - myNPTi->setTauThermostat( the_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( the_globals->haveTauBarostat() ) - myNPTi->setTauBarostat( the_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 NPTf_ENS: - myNPTf = new NPTf( info, the_ff ); - myNPTf->setTargetTemp( the_globals->getTargetTemp()); - - if (the_globals->haveTargetPressure()) - myNPTf->setTargetPressure(the_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( the_globals->haveTauThermostat() ) - myNPTf->setTauThermostat( the_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( the_globals->haveTauBarostat() ) - myNPTf->setTauBarostat( the_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 NPTim_ENS: - myNPTim = new NPTim( info, the_ff ); - myNPTim->setTargetTemp( the_globals->getTargetTemp()); - - if (the_globals->haveTargetPressure()) - myNPTim->setTargetPressure(the_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( the_globals->haveTauThermostat() ) - myNPTim->setTauThermostat( the_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( the_globals->haveTauBarostat() ) - myNPTim->setTauBarostat( the_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 NPTfm_ENS: - myNPTfm = new NPTfm( info, the_ff ); - myNPTfm->setTargetTemp( the_globals->getTargetTemp()); - - if (the_globals->haveTargetPressure()) - myNPTfm->setTargetPressure(the_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( the_globals->haveTauThermostat() ) - myNPTfm->setTauThermostat( the_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( the_globals->haveTauBarostat() ) - myNPTfm->setTauBarostat( the_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"); - painCave.isFatal = 1; - simError(); - } - - #ifdef IS_MPI mpiSim->mpiRefresh(); #endif // initialize the Fortran + initFortran(); - info->refreshSim(); - - if( !strcmp( info->mixingRule, "standard") ){ - the_ff->initForceField( LB_MIXING_RULE ); - } - else if( !strcmp( info->mixingRule, "explicit") ){ - the_ff->initForceField( EXPLICIT_MIXING_RULE ); - } - else{ - sprintf( painCave.errMsg, - "SimSetup Error: unknown mixing rule -> \"%s\"\n", - info->mixingRule ); - painCave.isFatal = 1; - simError(); - } -#ifdef IS_MPI - strcpy( checkPointMsg, - "Successfully intialized the mixingRule for Fortran." ); - MPIcheckPoint(); -#endif // is_mpi } void SimSetup::makeMolecules( void ){ int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; - molInit info; + molInit molInfo; DirectionalAtom* dAtom; LinkedAssign* extras; LinkedAssign* current_extra; @@ -370,32 +167,32 @@ void SimSetup::makeMolecules( void ){ stampID = the_molecules[i].getStampID(); - info.nAtoms = comp_stamps[stampID]->getNAtoms(); - info.nBonds = comp_stamps[stampID]->getNBonds(); - info.nBends = comp_stamps[stampID]->getNBends(); - info.nTorsions = comp_stamps[stampID]->getNTorsions(); - info.nExcludes = info.nBonds + info.nBends + info.nTorsions; + molInfo.nAtoms = comp_stamps[stampID]->getNAtoms(); + molInfo.nBonds = comp_stamps[stampID]->getNBonds(); + molInfo.nBends = comp_stamps[stampID]->getNBends(); + molInfo.nTorsions = comp_stamps[stampID]->getNTorsions(); + molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions; - info.myAtoms = &the_atoms[atomOffset]; - info.myExcludes = &the_excludes[excludeOffset]; - info.myBonds = new Bond*[info.nBonds]; - info.myBends = new Bend*[info.nBends]; - info.myTorsions = new Torsion*[info.nTorsions]; + molInfo.myAtoms = &the_atoms[atomOffset]; + molInfo.myExcludes = &the_excludes[excludeOffset]; + molInfo.myBonds = new Bond*[molInfo.nBonds]; + molInfo.myBends = new Bend*[molInfo.nBends]; + molInfo.myTorsions = new Torsion*[molInfo.nTorsions]; - theBonds = new bond_pair[info.nBonds]; - theBends = new bend_set[info.nBends]; - theTorsions = new torsion_set[info.nTorsions]; + theBonds = new bond_pair[molInfo.nBonds]; + theBends = new bend_set[molInfo.nBends]; + theTorsions = new torsion_set[molInfo.nTorsions]; // make the Atoms - for(j=0; jgetAtom( j ); if( currentAtom->haveOrientation() ){ dAtom = new DirectionalAtom(j + atomOffset); info->n_oriented++; - info.myAtoms[j] = dAtom; + molInfo.myAtoms[j] = dAtom; ux = currentAtom->getOrntX(); uy = currentAtom->getOrntY(); @@ -413,19 +210,19 @@ void SimSetup::makeMolecules( void ){ dAtom->setSUz( uz ); } else{ - info.myAtoms[j] = new GeneralAtom(j + atomOffset); + molInfo.myAtoms[j] = new GeneralAtom(j + atomOffset); } - info.myAtoms[j]->setType( currentAtom->getType() ); + molInfo.myAtoms[j]->setType( currentAtom->getType() ); #ifdef IS_MPI - info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] ); + molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] ); #endif // is_mpi } // make the bonds - for(j=0; jgetBond( j ); theBonds[j].a = currentBond->getA() + atomOffset; @@ -452,10 +249,10 @@ void SimSetup::makeMolecules( void ){ the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); #endif //is_mpi } - excludeOffset += info.nBonds; + excludeOffset += molInfo.nBonds; //make the bends - for(j=0; jgetBend( j ); theBends[j].a = currentBend->getA() + atomOffset; @@ -538,9 +335,9 @@ void SimSetup::makeMolecules( void ){ the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); #endif //is_mpi } - excludeOffset += info.nBends; + excludeOffset += molInfo.nBends; - for(j=0; jgetTorsion( j ); theTorsions[j].a = currentTorsion->getA() + atomOffset; @@ -568,21 +365,21 @@ void SimSetup::makeMolecules( void ){ the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); #endif //is_mpi } - excludeOffset += info.nTorsions; + excludeOffset += molInfo.nTorsions; // send the arrays off to the forceField for init. - the_ff->initializeAtoms( info.nAtoms, info.myAtoms ); - the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds ); - the_ff->initializeBends( info.nBends, info.myBends, theBends ); - the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions ); + the_ff->initializeAtoms( molInfo.nAtoms, molInfo.myAtoms ); + the_ff->initializeBonds( molInfo.nBonds, molInfo.myBonds, theBonds ); + the_ff->initializeBends( molInfo.nBends, molInfo.myBends, theBends ); + the_ff->initializeTorsions( molInfo.nTorsions, molInfo.myTorsions, theTorsions ); - the_molecules[i].initialize( info ); + the_molecules[i].initialize( molInfo ); - atomOffset += info.nAtoms; + atomOffset += molInfo.nAtoms; delete[] theBonds; delete[] theBends; delete[] theTorsions; @@ -618,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 ); @@ -635,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; @@ -784,27 +581,29 @@ void SimSetup::gatherInfo( void ){ void SimSetup::gatherInfo( void ){ + int i,j,k; ensembleCase = -1; ffCase = -1; // get the stamps and globals; - the_stamps = stamps; - the_globals = globals; + stamps = stamps; + globals = globals; // set the easy ones first - info->target_temp = the_globals->getTargetTemp(); - info->dt = the_globals->getDt(); - info->run_time = the_globals->getRunTime(); - n_components = the_globals->getNComponents(); + info->target_temp = globals->getTargetTemp(); + info->dt = globals->getDt(); + info->run_time = globals->getRunTime(); + n_components = globals->getNComponents(); // get the forceField - strcpy( force_field, the_globals->getForceField() ); + strcpy( force_field, globals->getForceField() ); 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", @@ -815,7 +614,7 @@ void SimSetup::gatherInfo( void ){ // get the ensemble - strcpy( ensemble, the_globals->getEnsemble() ); + strcpy( ensemble, globals->getEnsemble() ); if( !strcasecmp( ensemble, "NVE" )) ensembleCase = NVE_ENS; else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS; @@ -824,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, " @@ -838,17 +638,17 @@ void SimSetup::gatherInfo( void ){ // get the mixing rule - strcpy( info->mixingRule, the_globals->getMixingRule() ); - info->usePBC = the_globals->getPBC(); + strcpy( info->mixingRule, globals->getMixingRule() ); + info->usePBC = globals->getPBC(); // get the components and calculate the tot_nMol and indvidual n_mol - the_components = the_globals->getComponents(); + the_components = globals->getComponents(); components_nmol = new int[n_components]; - if( !the_globals->haveNMol() ){ + if( !globals->haveNMol() ){ // we don't have the total number of molecules, so we assume it is // given in each component @@ -881,43 +681,43 @@ void SimSetup::gatherInfo( void ){ // set the status, sample, and thermal kick times - if( the_globals->haveSampleTime() ){ - info->sampleTime = the_globals->getSampleTime(); + if( globals->haveSampleTime() ){ + info->sampleTime = globals->getSampleTime(); info->statusTime = info->sampleTime; info->thermalTime = info->sampleTime; } else{ - info->sampleTime = the_globals->getRunTime(); + info->sampleTime = globals->getRunTime(); info->statusTime = info->sampleTime; info->thermalTime = info->sampleTime; } - if( the_globals->haveStatusTime() ){ - info->statusTime = the_globals->getStatusTime(); + if( globals->haveStatusTime() ){ + info->statusTime = globals->getStatusTime(); } - if( the_globals->haveThermalTime() ){ - info->thermalTime = the_globals->getThermalTime(); + if( globals->haveThermalTime() ){ + info->thermalTime = globals->getThermalTime(); } // check for the temperature set flag - if( the_globals->haveTempSet() ) info->setTemp = the_globals->getTempSet(); + if( globals->haveTempSet() ) info->setTemp = globals->getTempSet(); // get some of the tricky things that may still be in the globals double boxVector[3]; - if( the_globals->haveBox() ){ - boxVector[0] = the_globals->getBox(); - boxVector[1] = the_globals->getBox(); - boxVector[2] = the_globals->getBox(); + if( globals->haveBox() ){ + boxVector[0] = globals->getBox(); + boxVector[1] = globals->getBox(); + boxVector[2] = globals->getBox(); info->setBox( boxVector ); } - else if( the_globals->haveDensity() ){ + else if( globals->haveDensity() ){ double vol; - vol = (double)tot_nmol / the_globals->getDensity(); + vol = (double)tot_nmol / globals->getDensity(); boxVector[0] = pow( vol, ( 1.0 / 3.0 ) ); boxVector[1] = boxVector[0]; boxVector[2] = boxVector[0]; @@ -925,29 +725,29 @@ void SimSetup::gatherInfo( void ){ info->setBox( boxVector ); } else{ - if( !the_globals->haveBoxX() ){ + if( !globals->haveBoxX() ){ sprintf( painCave.errMsg, "SimSetup error, no periodic BoxX size given.\n" ); painCave.isFatal = 1; simError(); } - boxVector[0] = the_globals->getBoxX(); + boxVector[0] = globals->getBoxX(); - if( !the_globals->haveBoxY() ){ + if( !globals->haveBoxY() ){ sprintf( painCave.errMsg, "SimSetup error, no periodic BoxY size given.\n" ); painCave.isFatal = 1; simError(); } - boxVector[1] = the_globals->getBoxY(); + boxVector[1] = globals->getBoxY(); - if( !the_globals->haveBoxZ() ){ + if( !globals->haveBoxZ() ){ sprintf( painCave.errMsg, "SimSetup error, no periodic BoxZ size given.\n" ); painCave.isFatal = 1; simError(); } - boxVector[2] = the_globals->getBoxZ(); + boxVector[2] = globals->getBoxZ(); info->setBox( boxVector ); } @@ -977,15 +777,16 @@ void SimSetup::finalInfoCheck( void ){ } #ifdef IS_MPI - int myUse = usesDipoles - MPI_Allreduce( &myUse, &UsesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD ); + int myUse = usesDipoles; + MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD ); #endif //is_mpi + double theEcr, theEst; - if (the_globals->getUseRF() ) { + if (globals->getUseRF() ) { info->useReactionField = 1; - if( !the_globals->haveECR() ){ + if( !globals->haveECR() ){ sprintf( painCave.errMsg, "SimSetup Warning: using default value of 1/2 the smallest " "box length for the electrostaticCutoffRadius.\n" @@ -993,27 +794,29 @@ 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 = the_globals->getECR(); + theEcr = globals->getECR(); } - if( !the_globals->haveEST() ){ + if( !globals->haveEST() ){ sprintf( painCave.errMsg, "SimSetup Warning: using default value of 0.05 * the " "electrostaticCutoffRadius for the electrostaticSkinThickness\n" ); painCave.isFatal = 0; simError(); - info->est = 0.05 * info->ecr; + theEst = 0.05 * theEcr; } else { - info->est = the_globals->getEST(); + theEst= globals->getEST(); } + + info->setEcr( theEcr, theEst ); - if(!the_globals->haveDielectric() ){ + if(!globals->haveDielectric() ){ sprintf( painCave.errMsg, "SimSetup Error: You are trying to use Reaction Field without" "setting a dielectric constant!\n" @@ -1021,39 +824,41 @@ void SimSetup::finalInfoCheck( void ){ painCave.isFatal = 1; simError(); } - info->dielectric = the_globals->getDielectric(); + info->dielectric = globals->getDielectric(); } else { if (usesDipoles) { - if( !the_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; + 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->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 = the_globals->getECR(); + theEcr = globals->getECR(); } - if( !the_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; + if( !globals->haveEST() ){ + 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 = the_globals->getEST(); + theEst= globals->getEST(); } + + info->setEcr( theEcr, theEst ); } } @@ -1066,17 +871,17 @@ void SimSetup::initSystemCoords( void ){ void SimSetup::initSystemCoords( void ){ - if( the_globals->haveInitialConfig() ){ + if( globals->haveInitialConfig() ){ InitializeFromFile* fileInit; #ifdef IS_MPI // is_mpi if( worldRank == 0 ){ #endif //is_mpi - fileInit = new InitializeFromFile( the_globals->getInitialConfig() ); + fileInit = new InitializeFromFile( globals->getInitialConfig() ); #ifdef IS_MPI }else fileInit = new InitializeFromFile( NULL ); #endif - fileInit->read_xyz( info ); // default velocities on + fileInit->readInit( info ); // default velocities on delete fileInit; } @@ -1113,8 +918,8 @@ void SimSetup::makeOutNames( void ){ if( worldRank == 0 ){ #endif // is_mpi - if( the_globals->haveFinalConfig() ){ - strcpy( info->finalName, the_globals->getFinalConfig() ); + if( globals->haveFinalConfig() ){ + strcpy( info->finalName, globals->getFinalConfig() ); } else{ strcpy( info->finalName, inFileName ); @@ -1197,6 +1002,8 @@ void SimSetup::sysObjectsCreation( void ){ void SimSetup::sysObjectsCreation( void ){ + int i; + // create the forceField createFF(); @@ -1219,6 +1026,13 @@ void SimSetup::sysObjectsCreation( void ){ makeSysArrays(); + // make and initialize the molecules (all but atomic coordinates) + + makeMolecules(); + info->identArray = new int[info->n_atoms]; + for(i=0; in_atoms; i++){ + info->identArray[i] = the_atoms[i]->getIdent(); + } @@ -1237,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"); @@ -1254,6 +1072,8 @@ void SimSetup::compList( void ){ void SimSetup::compList( void ){ + int i; + comp_stamps = new MoleculeStamp*[n_components]; // make an array of molecule stamps that match the components used. @@ -1279,7 +1099,7 @@ void SimSetup::compList( void ){ // extract the component from the list; - currentStamp = the_stamps->extractMolStamp( id ); + currentStamp = stamps->extractMolStamp( id ); if( currentStamp == NULL ){ sprintf( painCave.errMsg, "SimSetup error: Component \"%s\" was not found in the " @@ -1303,7 +1123,9 @@ void SimSetup::calcSysValues( void ){ } void SimSetup::calcSysValues( void ){ + int i, j, k; + tot_atoms = 0; tot_bonds = 0; tot_bends = 0; @@ -1333,6 +1155,7 @@ void SimSetup::mpiMolDivide( void ){ void SimSetup::mpiMolDivide( void ){ + int i, j, k; int localMol, allMol; int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; @@ -1402,7 +1225,9 @@ void SimSetup::makeSysArrays( void ){ void SimSetup::makeSysArrays( void ){ + int i, j, k; + // create the atom and short range interaction arrays Atom::createArrays(info->n_atoms); @@ -1476,3 +1301,270 @@ void SimSetup::makeSysArrays( void ){ the_ff->setSimInfo( info ); } + +void SimSetup::makeIntegrator( void ){ + + 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->setTargetTemp(globals->getTargetTemp()); + + if (globals->haveTauThermostat()) + myNVT->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 NPTi_ENS: + myNPTi = new NPTi( info, the_ff ); + myNPTi->setTargetTemp( globals->getTargetTemp() ); + + if (globals->haveTargetPressure()) + myNPTi->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() ) + myNPTi->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() ) + myNPTi->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 NPTf_ENS: + myNPTf = new NPTf( info, the_ff ); + myNPTf->setTargetTemp( globals->getTargetTemp()); + + if (globals->haveTargetPressure()) + myNPTf->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() ) + myNPTf->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() ) + myNPTf->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 NPTim_ENS: + myNPTim = new NPTim( info, the_ff ); + myNPTim->setTargetTemp( globals->getTargetTemp()); + + if (globals->haveTargetPressure()) + myNPTim->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() ) + myNPTim->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() ) + myNPTim->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 NPTfm_ENS: + myNPTfm = new NPTfm( info, the_ff ); + myNPTfm->setTargetTemp( globals->getTargetTemp()); + + if (globals->haveTargetPressure()) + myNPTfm->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() ) + myNPTfm->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() ) + myNPTfm->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 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"); + painCave.isFatal = 1; + simError(); + } + +} + +void SimSetup::initFortran( void ){ + + info->refreshSim(); + + if( !strcmp( info->mixingRule, "standard") ){ + the_ff->initForceField( LB_MIXING_RULE ); + } + else if( !strcmp( info->mixingRule, "explicit") ){ + the_ff->initForceField( EXPLICIT_MIXING_RULE ); + } + else{ + sprintf( painCave.errMsg, + "SimSetup Error: unknown mixing rule -> \"%s\"\n", + info->mixingRule ); + painCave.isFatal = 1; + simError(); + } + + +#ifdef IS_MPI + strcpy( checkPointMsg, + "Successfully intialized the mixingRule for Fortran." ); + MPIcheckPoint(); +#endif // is_mpi + +}