--- trunk/mdtools/interface_implementation/SimSetup.cpp 2002/07/09 18:40:59 11 +++ trunk/mdtools/interface_implementation/SimSetup.cpp 2002/12/19 21:59:51 215 @@ -6,10 +6,22 @@ #include "parse_me.h" #include "LRI.hpp" #include "Integrator.hpp" +#include "simError.h" +#ifdef IS_MPI +#include "mpiBASS.h" +#include "mpiSimulation.hpp" +#include "bassDiag.hpp" +#endif + SimSetup::SimSetup(){ stamps = new MakeStamps(); globals = new Globals(); + +#ifdef IS_MPI + strcpy( checkPointMsg, "SimSetup creation successful" ); + MPIcheckPoint(); +#endif // IS_MPI } SimSetup::~SimSetup(){ @@ -19,16 +31,50 @@ void SimSetup::parseFile( char* fileName ){ void SimSetup::parseFile( char* fileName ){ - inFileName = fileName; - set_interface_stamps( stamps, globals ); - yacc_BASS( fileName ); +#ifdef IS_MPI + if( worldRank == 0 ){ +#endif // is_mpi + + inFileName = fileName; + set_interface_stamps( stamps, globals ); + +#ifdef IS_MPI + mpiEventInit(); +#endif + + yacc_BASS( fileName ); + +#ifdef IS_MPI + throwMPIEvent(NULL); + } + else receiveParse(); +#endif + } +#ifdef IS_MPI +void SimSetup::receiveParse(void){ + + set_interface_stamps( stamps, globals ); + mpiEventInit(); + MPIcheckPoint(); + mpiEventLoop(); + +} + + +void SimSetup::testMe(void){ + bassDiag* dumpMe = new bassDiag(globals,stamps); + dumpMe->dumpStamps(); + delete dumpMe; +} +#endif + void SimSetup::createSim( void ){ MakeStamps *the_stamps; Globals* the_globals; - int i; + int i, j; // get the stamps and globals; the_stamps = stamps; @@ -43,33 +89,43 @@ void SimSetup::createSim( void ){ n_components = the_globals->getNComponents(); strcpy( force_field, the_globals->getForceField() ); strcpy( ensemble, the_globals->getEnsemble() ); - + if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF(); else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF(); else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF(); + else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF(); else{ - std::cerr<< "SimSetup Error. Unrecognized force field -> " - << force_field << "\n"; - exit(8); + sprintf( painCave.errMsg, + "SimSetup Error. Unrecognized force field -> %s\n", + force_field ); + painCave.isFatal = 1; + simError(); } +#ifdef IS_MPI + strcpy( checkPointMsg, "ForceField creation successful" ); + MPIcheckPoint(); +#endif // is_mpi + // get the components and calculate the tot_nMol and indvidual n_mol the_components = the_globals->getComponents(); components_nmol = new int[n_components]; comp_stamps = new MoleculeStamp*[n_components]; - + if( !the_globals->haveNMol() ){ - // we don't have the total number of molecules, so we assume it is + // we don't have the total number of molecules, so we assume it is // given in each component tot_nmol = 0; for( i=0; ihaveNMol() ){ // we have a problem - std::cerr << "SimSetup Error. No global NMol or component NMol" - << " given. Cannot calculate the number of atoms.\n"; - exit( 8 ); + sprintf( painCave.errMsg, + "SimSetup Error. No global NMol or component NMol" + " given. Cannot calculate the number of atoms.\n" ); + painCave.isFatal = 1; + simError(); } tot_nmol += the_components[i]->getNMol(); @@ -77,32 +133,81 @@ void SimSetup::createSim( void ){ } } else{ - std::cerr << "NOT A SUPPORTED FEATURE\n"; + sprintf( painCave.errMsg, + "SimSetup error.\n" + "\tSorry, the ability to specify total" + " nMols and then give molfractions in the components\n" + "\tis not currently supported." + " Please give nMol in the components.\n" ); + painCave.isFatal = 1; + simError(); -// tot_nmol = the_globals->getNMol(); -// //we have the total number of molecules, now we check for molfractions -// for( i=0; ihaveMolFraction() ){ - -// if( !the_components[i]->haveNMol() ){ -// //we have a problem -// std::cerr << "SimSetup error. Neither molFraction nor " -// << " nMol was given in component - + // tot_nmol = the_globals->getNMol(); + + // //we have the total number of molecules, now we check for molfractions + // for( i=0; ihaveMolFraction() ){ + + // if( !the_components[i]->haveNMol() ){ + // //we have a problem + // std::cerr << "SimSetup error. Neither molFraction nor " + // << " nMol was given in component + } +#ifdef IS_MPI + strcpy( checkPointMsg, "Have the number of components" ); + MPIcheckPoint(); +#endif // is_mpi + // make an array of molecule stamps that match the components used. + // also extract the used stamps out into a separate linked list + simnfo->nComponents = n_components; + simnfo->componentsNmol = components_nmol; + simnfo->compStamps = comp_stamps; + simnfo->headStamp = new LinkedMolStamp(); + + char* id; + LinkedMolStamp* headStamp = simnfo->headStamp; + LinkedMolStamp* currentStamp = NULL; for( i=0; igetMolecule( the_components[i]->getType() ); + id = the_components[i]->getType(); + comp_stamps[i] = NULL; + + // check to make sure the component isn't already in the list + + comp_stamps[i] = headStamp->match( id ); + if( comp_stamps[i] == NULL ){ + + // extract the component from the list; + + currentStamp = the_stamps->extractMolStamp( id ); + if( currentStamp == NULL ){ + sprintf( painCave.errMsg, + "SimSetup error: Component \"%s\" was not found in the " + "list of declared molecules\n" + id ); + painCave.isFatal = 1; + simError(); + } + + headStamp->add( currentStamp ); + comp_stamps[i] = headStamp->match( id ); + } } +#ifdef IS_MPI + strcpy( checkPointMsg, "Component stamps successfully extracted\n" ); + MPIcheckPoint(); +#endif // is_mpi + + // caclulate the number of atoms, bonds, bends and torsions tot_atoms = 0; @@ -111,31 +216,97 @@ void SimSetup::createSim( void ){ tot_torsions = 0; for( i=0; igetNAtoms(); - tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); - tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); + tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); + tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); + tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); } - + tot_SRI = tot_bonds + tot_bends + tot_torsions; - + simnfo->n_atoms = tot_atoms; simnfo->n_bonds = tot_bonds; simnfo->n_bends = tot_bends; simnfo->n_torsions = tot_torsions; simnfo->n_SRI = tot_SRI; + simnfo->n_mol = tot_nmol; - // create the atom and short range interaction arrays - the_atoms = new Atom*[tot_atoms]; - // the_molecules = new Molecule[tot_nmol]; +#ifdef IS_MPI + + // divide the molecules among processors here. + mpiSim = new mpiSimulation( simnfo ); - if( tot_SRI ){ - the_sris = new SRI*[tot_SRI]; - the_excludes = new ex_pair[tot_SRI]; + mpiSim->divideLabor(); + + // set up the local variables + + int localMol, allMol; + int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; + + allMol = 0; + localMol = 0; + local_atoms = 0; + local_bonds = 0; + local_bends = 0; + local_torsions = 0; + for( i=0; igetMyMolStart() <= allMol && + allMol <= mpiSim->getMyMolEnd() ){ + + local_atoms += comp_stamps[i]->getNAtoms(); + local_bonds += comp_stamps[i]->getNBonds(); + local_bends += comp_stamps[i]->getNBends(); + local_torsions += comp_stamps[i]->getNTorsions(); + localMol++; + } + allMol++; + } } + local_SRI = local_bonds + local_bends + local_torsions; + + simnfo->n_atoms = mpiSim->getMyNlocal(); + + if( local_atoms != simnfo->n_atoms ){ + sprintf( painCave.errMsg, + "SimSetup error: mpiSim's localAtom (%d) and SimSetup's" + " localAtom (%d) are note equal.\n", + simnfo->n_atoms, + local_atoms ); + painCave.isFatal = 1; + simError(); + } + + simnfo->n_bonds = local_bonds; + simnfo->n_bends = local_bends; + simnfo->n_torsions = local_torsions; + simnfo->n_SRI = local_SRI; + simnfo->n_mol = localMol; + + strcpy( checkPointMsg, "Passed nlocal consistency check." ); + MPIcheckPoint(); + + +#endif // is_mpi + + + // create the atom and short range interaction arrays + + Atom::createArrays(simnfo->n_atoms); + the_atoms = new Atom*[simnfo->n_atoms]; + the_molecules = new Molecule[simnfo->n_mol]; + + + if( simnfo->n_SRI ){ + the_sris = new SRI*[simnfo->n_SRI]; + the_excludes = new ex_pair[simnfo->n_SRI]; + } + // set the arrays into the SimInfo object simnfo->atoms = the_atoms; @@ -143,10 +314,11 @@ void SimSetup::createSim( void ){ simnfo->n_exclude = tot_SRI; simnfo->excludes = the_excludes; + // initialize the arrays - + the_ff->setSimInfo( simnfo ); - + makeAtoms(); if( tot_bonds ){ @@ -161,33 +333,41 @@ void SimSetup::createSim( void ){ makeTorsions(); } - // makeMolecules(); // get some of the tricky things that may still be in the globals if( simnfo->n_dipoles ){ if( !the_globals->haveRRF() ){ - std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n"; - exit(8); + sprintf( painCave.errMsg, + "SimSetup Error, system has dipoles, but no rRF was set.\n"); + painCave.isFatal = 1; + simError(); } if( !the_globals->haveDielectric() ){ - std::cerr << "SimSetup Error, system has dipoles, but no" - << " dielectric was set.\n"; - exit(8); + sprintf( painCave.errMsg, + "SimSetup Error, system has dipoles, but no" + " dielectric was set.\n" ); + painCave.isFatal = 1; + simError(); } simnfo->rRF = the_globals->getRRF(); simnfo->dielectric = the_globals->getDielectric(); } +#ifdef IS_MPI + strcpy( checkPointMsg, "rRf and dielectric check out" ); + MPIcheckPoint(); +#endif // is_mpi + if( the_globals->haveBox() ){ simnfo->box_x = the_globals->getBox(); simnfo->box_y = the_globals->getBox(); simnfo->box_z = the_globals->getBox(); } else if( the_globals->haveDensity() ){ - + double vol; vol = (double)tot_nmol / the_globals->getDensity(); simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) ); @@ -196,120 +376,173 @@ void SimSetup::createSim( void ){ } else{ if( !the_globals->haveBoxX() ){ - std::cerr << "SimSetup error, no periodic BoxX size given.\n"; - exit(8); + sprintf( painCave.errMsg, + "SimSetup error, no periodic BoxX size given.\n" ); + painCave.isFatal = 1; + simError(); } simnfo->box_x = the_globals->getBoxX(); if( !the_globals->haveBoxY() ){ - std::cerr << "SimSetup error, no periodic BoxY size given.\n"; - exit(8); + sprintf( painCave.errMsg, + "SimSetup error, no periodic BoxY size given.\n" ); + painCave.isFatal = 1; + simError(); } simnfo->box_y = the_globals->getBoxY(); if( !the_globals->haveBoxZ() ){ - std::cerr << "SimSetup error, no periodic BoxZ size given.\n"; - exit(8); + sprintf( painCave.errMsg, + "SimSetup error, no periodic BoxZ size given.\n" ); + painCave.isFatal = 1; + simError(); } simnfo->box_z = the_globals->getBoxZ(); } - - if( the_globals->haveInitialConfig() ){ - InitializeFromFile* fileInit; - fileInit = new InitializeFromFile( the_globals->getInitialConfig() ); - - fileInit->read_xyz( simnfo ); // default velocities on - delete fileInit; - } - else{ - initFromBass(); - } +#ifdef IS_MPI + strcpy( checkPointMsg, "Box size set up" ); + MPIcheckPoint(); +#endif // is_mpi - if( the_globals->haveFinalConfig() ){ - strcpy( simnfo->finalName, the_globals->getFinalConfig() ); - } - else{ - strcpy( simnfo->finalName, inFileName ); - char* endTest; - int nameLength = strlen( simnfo->finalName ); - endTest = &(simnfo->finalName[nameLength - 5]); - if( !strcmp( endTest, ".bass" ) ){ - strcpy( endTest, ".eor" ); + + + if( the_globals->haveInitialConfig() ){ + + InitializeFromFile* fileInit; +#ifdef IS_MPI // is_mpi + if( worldRank == 0 ){ +#endif //is_mpi + fileInit = new InitializeFromFile( the_globals->getInitialConfig() ); +#ifdef IS_MPI + }else fileInit = new InitializeFromFile( NULL ); +#endif + fileInit->read_xyz( simnfo ); // default velocities on + + delete fileInit; + } + else{ + +#ifdef IS_MPI + + // no init from bass + + sprintf( painCave.errMsg, + "Cannot intialize a parallel simulation without an initial configuration file.\n" ); + painCave.isFatal; + simError(); + +#else + + initFromBass(); + + +#endif + } + +#ifdef IS_MPI + strcpy( checkPointMsg, "Successfully read in the initial configuration" ); + MPIcheckPoint(); +#endif // is_mpi + + + + + + + +#ifdef IS_MPI + if( worldRank == 0 ){ +#endif // is_mpi + + if( the_globals->haveFinalConfig() ){ + strcpy( simnfo->finalName, the_globals->getFinalConfig() ); } - else if( !strcmp( endTest, ".BASS" ) ){ - strcpy( endTest, ".eor" ); - } else{ - endTest = &(simnfo->finalName[nameLength - 4]); - if( !strcmp( endTest, ".bss" ) ){ + strcpy( simnfo->finalName, inFileName ); + char* endTest; + int nameLength = strlen( simnfo->finalName ); + endTest = &(simnfo->finalName[nameLength - 5]); + if( !strcmp( endTest, ".bass" ) ){ strcpy( endTest, ".eor" ); } - else if( !strcmp( endTest, ".mdl" ) ){ + else if( !strcmp( endTest, ".BASS" ) ){ strcpy( endTest, ".eor" ); } else{ - strcat( simnfo->finalName, ".eor" ); + endTest = &(simnfo->finalName[nameLength - 4]); + if( !strcmp( endTest, ".bss" ) ){ + strcpy( endTest, ".eor" ); + } + else if( !strcmp( endTest, ".mdl" ) ){ + strcpy( endTest, ".eor" ); + } + else{ + strcat( simnfo->finalName, ".eor" ); + } } } - } - // make the sample and status out names - - strcpy( simnfo->sampleName, inFileName ); - char* endTest; - int nameLength = strlen( simnfo->sampleName ); - endTest = &(simnfo->sampleName[nameLength - 5]); - if( !strcmp( endTest, ".bass" ) ){ - strcpy( endTest, ".dump" ); - } - else if( !strcmp( endTest, ".BASS" ) ){ - strcpy( endTest, ".dump" ); - } - else{ - endTest = &(simnfo->sampleName[nameLength - 4]); - if( !strcmp( endTest, ".bss" ) ){ + // make the sample and status out names + + strcpy( simnfo->sampleName, inFileName ); + char* endTest; + int nameLength = strlen( simnfo->sampleName ); + endTest = &(simnfo->sampleName[nameLength - 5]); + if( !strcmp( endTest, ".bass" ) ){ strcpy( endTest, ".dump" ); } - else if( !strcmp( endTest, ".mdl" ) ){ + else if( !strcmp( endTest, ".BASS" ) ){ strcpy( endTest, ".dump" ); } else{ - strcat( simnfo->sampleName, ".dump" ); + endTest = &(simnfo->sampleName[nameLength - 4]); + if( !strcmp( endTest, ".bss" ) ){ + strcpy( endTest, ".dump" ); + } + else if( !strcmp( endTest, ".mdl" ) ){ + strcpy( endTest, ".dump" ); + } + else{ + strcat( simnfo->sampleName, ".dump" ); + } } - } - - strcpy( simnfo->statusName, inFileName ); - nameLength = strlen( simnfo->statusName ); - endTest = &(simnfo->statusName[nameLength - 5]); - if( !strcmp( endTest, ".bass" ) ){ - strcpy( endTest, ".stat" ); - } - else if( !strcmp( endTest, ".BASS" ) ){ - strcpy( endTest, ".stat" ); - } - else{ - endTest = &(simnfo->statusName[nameLength - 4]); - if( !strcmp( endTest, ".bss" ) ){ + + strcpy( simnfo->statusName, inFileName ); + nameLength = strlen( simnfo->statusName ); + endTest = &(simnfo->statusName[nameLength - 5]); + if( !strcmp( endTest, ".bass" ) ){ strcpy( endTest, ".stat" ); } - else if( !strcmp( endTest, ".mdl" ) ){ + else if( !strcmp( endTest, ".BASS" ) ){ strcpy( endTest, ".stat" ); } else{ - strcat( simnfo->statusName, ".stat" ); + endTest = &(simnfo->statusName[nameLength - 4]); + if( !strcmp( endTest, ".bss" ) ){ + strcpy( endTest, ".stat" ); + } + else if( !strcmp( endTest, ".mdl" ) ){ + strcpy( endTest, ".stat" ); + } + else{ + strcat( simnfo->statusName, ".stat" ); + } } + +#ifdef IS_MPI } +#endif // is_mpi // set the status, sample, and themal kick times - + if( the_globals->haveSampleTime() ){ - simnfo->sampleTime = the_globals->getSampleTime(); + simnfo->sampleTime = the_globals->getSampleTime(); simnfo->statusTime = simnfo->sampleTime; simnfo->thermalTime = simnfo->sampleTime; } else{ - simnfo->sampleTime = the_globals->getRunTime(); + simnfo->sampleTime = the_globals->getRunTime(); simnfo->statusTime = simnfo->sampleTime; simnfo->thermalTime = simnfo->sampleTime; } @@ -325,148 +558,202 @@ void SimSetup::createSim( void ){ // check for the temperature set flag if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet(); - - + + // make the longe range forces and the integrator - + new AllLong( simnfo ); - + if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo ); if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo ); if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo ); } void SimSetup::makeAtoms( void ){ - + int i, j, k, index; double ux, uy, uz, uSqr, u; AtomStamp* current_atom; DirectionalAtom* dAtom; + int molIndex, molStart, molEnd, nMemb, lMolIndex; + lMolIndex = 0; + molIndex = 0; index = 0; for( i=0; igetNAtoms(); k++ ){ - - current_atom = comp_stamps[i]->getAtom( k ); - if( current_atom->haveOrientation() ){ - dAtom = new DirectionalAtom; - simnfo->n_oriented++; - the_atoms[index] = dAtom; - - ux = current_atom->getOrntX(); - uy = current_atom->getOrntY(); - uz = current_atom->getOrntZ(); +#ifdef IS_MPI + if( mpiSim->getMyMolStart() <= molIndex && + molIndex <= mpiSim->getMyMolEnd() ){ +#endif // is_mpi + + molStart = index; + nMemb = comp_stamps[i]->getNAtoms(); + for( k=0; kgetNAtoms(); k++ ){ - uSqr = (ux * ux) + (uy * uy) + (uz * uz); + current_atom = comp_stamps[i]->getAtom( k ); + if( current_atom->haveOrientation() ){ + + dAtom = new DirectionalAtom(index); + simnfo->n_oriented++; + the_atoms[index] = dAtom; + + ux = current_atom->getOrntX(); + uy = current_atom->getOrntY(); + uz = current_atom->getOrntZ(); + + uSqr = (ux * ux) + (uy * uy) + (uz * uz); + + u = sqrt( uSqr ); + ux = ux / u; + uy = uy / u; + uz = uz / u; + + dAtom->setSUx( ux ); + dAtom->setSUy( uy ); + dAtom->setSUz( uz ); + } + else{ + the_atoms[index] = new GeneralAtom(index); + } + the_atoms[index]->setType( current_atom->getType() ); + the_atoms[index]->setIndex( index ); - u = sqrt( uSqr ); - ux = ux / u; - uy = uy / u; - uz = uz / u; - - dAtom->setSUx( ux ); - dAtom->setSUy( uy ); - dAtom->setSUz( uz ); + // increment the index and repeat; + index++; } - else{ - the_atoms[index] = new GeneralAtom; - } - the_atoms[index]->setType( current_atom->getType() ); - the_atoms[index]->setIndex( index ); - // increment the index and repeat; - index++; + molEnd = index -1; + the_molecules[lMolIndex].setNMembers( nMemb ); + the_molecules[lMolIndex].setStartAtom( molStart ); + the_molecules[lMolIndex].setEndAtom( molEnd ); + the_molecules[lMolIndex].setStampID( i ); + lMolIndex++; + +#ifdef IS_MPI } +#endif //is_mpi + + molIndex++; } } - + the_ff->initializeAtoms(); } void SimSetup::makeBonds( void ){ - int i, j, k, index, offset; + int i, j, k, index, offset, molIndex; bond_pair* the_bonds; BondStamp* current_bond; the_bonds = new bond_pair[tot_bonds]; index = 0; offset = 0; + molIndex = 0;g1 + for( i=0; igetNBonds(); k++ ){ - - current_bond = comp_stamps[i]->getBond( k ); - the_bonds[index].a = current_bond->getA() + offset; - the_bonds[index].b = current_bond->getB() + offset; - the_excludes[index].i = the_bonds[index].a; - the_excludes[index].j = the_bonds[index].b; - - // increment the index and repeat; - index++; - } - offset += comp_stamps[i]->getNAtoms(); - } +#ifdef IS_MPI + if( mpiSim->getMyMolStart() <= molIndex && + molIndex <= mpiSim->getMyMolEnd() ){ +#endif // is_mpi + + for( k=0; kgetNBonds(); k++ ){ + + current_bond = comp_stamps[i]->getBond( k ); + the_bonds[index].a = current_bond->getA() + offset; + the_bonds[index].b = current_bond->getB() + offset; + + the_excludes[index].i = the_bonds[index].a; + the_excludes[index].j = the_bonds[index].b; + + // increment the index and repeat; + index++; + } + offset += comp_stamps[i]->getNAtoms(); + +#ifdef IS_MPI + } +#endif is_mpi + + molIndex++; + } } - + the_ff->initializeBonds( the_bonds ); } void SimSetup::makeBends( void ){ - int i, j, k, index, offset; + int i, j, k, index, offset, molIndex; bend_set* the_bends; BendStamp* current_bend; the_bends = new bend_set[tot_bends]; index = 0; offset = 0; + molIndex = 0; for( i=0; igetNBends(); k++ ){ - - current_bend = comp_stamps[i]->getBend( k ); - the_bends[index].a = current_bend->getA() + offset; - the_bends[index].b = current_bend->getB() + offset; - the_bends[index].c = current_bend->getC() + offset; - the_excludes[index + tot_bonds].i = the_bends[index].a; - the_excludes[index + tot_bonds].j = the_bends[index].c; +#ifdef IS_MPI + if( mpiSim->getMyMolStart() <= molIndex && + molIndex <= mpiSim->getMyMolEnd() ){ +#endif // is_mpi - // increment the index and repeat; - index++; + for( k=0; kgetNBends(); k++ ){ + + current_bend = comp_stamps[i]->getBend( k ); + the_bends[index].a = current_bend->getA() + offset; + the_bends[index].b = current_bend->getB() + offset; + the_bends[index].c = current_bend->getC() + offset; + + the_excludes[index + tot_bonds].i = the_bends[index].a; + the_excludes[index + tot_bonds].j = the_bends[index].c; + + // increment the index and repeat; + index++; + } + offset += comp_stamps[i]->getNAtoms(); + +#ifdef IS_MPI } - offset += comp_stamps[i]->getNAtoms(); +#endif //is_mpi + + molIndex++; } } - + the_ff->initializeBends( the_bends ); } void SimSetup::makeTorsions( void ){ - int i, j, k, index, offset; + int i, j, k, index, offset, molIndex; torsion_set* the_torsions; TorsionStamp* current_torsion; the_torsions = new torsion_set[tot_torsions]; index = 0; offset = 0; + molIndex = 0; for( i=0; igetMyMolStart() <= molIndex && + molIndex <= mpiSim->getMyMolEnd() ){ +#endif // is_mpi + for( k=0; kgetNTorsions(); k++ ){ - + current_torsion = comp_stamps[i]->getTorsion( k ); the_torsions[index].a = current_torsion->getA() + offset; the_torsions[index].b = current_torsion->getB() + offset; @@ -480,17 +767,18 @@ void SimSetup::makeTorsions( void ){ index++; } offset += comp_stamps[i]->getNAtoms(); + +#ifdef IS_MPI + } +#endif //is_mpi + + molIndex++; } } - + the_ff->initializeTorsions( the_torsions ); } -void SimSetup::makeMolecules( void ){ - - //empy for now -} - void SimSetup::initFromBass( void ){ int i, j, k; @@ -518,8 +806,11 @@ void SimSetup::initFromBass( void ){ n_per_extra = (int)ceil( temp1 ); if( n_per_extra > 4){ - std::cerr << "THere has been an error in constructing the non-complete lattice.\n"; - exit(8); + sprintf( painCave.errMsg, + "SimSetup error. There has been an error in constructing" + " the non-complete lattice.\n" ); + painCave.isFatal = 1; + simError(); } } else{ @@ -528,28 +819,28 @@ void SimSetup::initFromBass( void ){ celly = simnfo->box_y / temp3; cellz = simnfo->box_z / temp3; } - + current_mol = 0; current_comp_mol = 0; current_comp = 0; current_atom_ndx = 0; - + for( i=0; i < n_cells ; i++ ){ for( j=0; j < n_cells; j++ ){ for( k=0; k < n_cells; k++ ){ - + makeElement( i * cellx, j * celly, k * cellz ); - + makeElement( i * cellx + 0.5 * cellx, j * celly + 0.5 * celly, k * cellz ); - + makeElement( i * cellx, j * celly + 0.5 * celly, k * cellz + 0.5 * cellz ); - + makeElement( i * cellx + 0.5 * cellx, j * celly, k * cellz + 0.5 * cellz ); @@ -559,41 +850,41 @@ void SimSetup::initFromBass( void ){ if( have_extra ){ done = 0; - + int start_ndx; for( i=0; i < (n_cells+1) && !done; i++ ){ for( j=0; j < (n_cells+1) && !done; j++ ){ - + if( i < n_cells ){ - + if( j < n_cells ){ start_ndx = n_cells; } else start_ndx = 0; } else start_ndx = 0; - + for( k=start_ndx; k < (n_cells+1) && !done; k++ ){ - + makeElement( i * cellx, j * celly, k * cellz ); done = ( current_mol >= tot_nmol ); - + if( !done && n_per_extra > 1 ){ makeElement( i * cellx + 0.5 * cellx, j * celly + 0.5 * celly, k * cellz ); done = ( current_mol >= tot_nmol ); } - + if( !done && n_per_extra > 2){ makeElement( i * cellx, j * celly + 0.5 * celly, k * cellz + 0.5 * cellz ); done = ( current_mol >= tot_nmol ); } - + if( !done && n_per_extra > 3){ makeElement( i * cellx + 0.5 * cellx, j * celly, @@ -604,8 +895,8 @@ void SimSetup::initFromBass( void ){ } } } - - + + for( i=0; in_atoms; i++ ){ simnfo->atoms[i]->set_vx( 0.0 ); simnfo->atoms[i]->set_vy( 0.0 ); @@ -621,25 +912,28 @@ void SimSetup::makeElement( double x, double y, double double rotMat[3][3]; for( k=0; kgetNAtoms(); k++ ){ - + current_atom = comp_stamps[current_comp]->getAtom( k ); if( !current_atom->havePosition() ){ - std::cerr << "Component " << comp_stamps[current_comp]->getID() - << ", atom " << current_atom->getType() - << " does not have a position specified.\n" - << "The initialization routine is unable to give a start" - << " position.\n"; - exit(8); + sprintf( painCave.errMsg, + "SimSetup:initFromBass error.\n" + "\tComponent %s, atom %s does not have a position specified.\n" + "\tThe initialization routine is unable to give a start" + " position.\n", + comp_stamps[current_comp]->getID(), + current_atom->getType() ); + painCave.isFatal = 1; + simError(); } - + the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() ); the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() ); the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() ); - + if( the_atoms[current_atom_ndx]->isDirectional() ){ - + dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx]; - + rotMat[0][0] = 1.0; rotMat[0][1] = 0.0; rotMat[0][2] = 0.0; @@ -657,12 +951,12 @@ void SimSetup::makeElement( double x, double y, double current_atom_ndx++; } - + current_mol++; current_comp_mol++; if( current_comp_mol >= components_nmol[current_comp] ){ - + current_comp_mol = 0; current_comp++; }