--- trunk/OOPSE/libmdtools/SimSetup.cpp 2003/03/26 22:02:36 414 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2003/04/07 21:42:19 474 @@ -66,6 +66,7 @@ void SimSetup::createSim( void ){ MakeStamps *the_stamps; Globals* the_globals; + ExtendedSystem* the_extendedsystem; int i, j; // get the stamps and globals; @@ -80,18 +81,40 @@ void SimSetup::createSim( void ){ // get the ones we know are there, yet still may need some work. n_components = the_globals->getNComponents(); strcpy( force_field, the_globals->getForceField() ); + + // get the ensemble and set up an extended system if we need it: strcpy( ensemble, the_globals->getEnsemble() ); + if( !strcasecmp( ensemble, "NPT" ) ) { + the_extendedsystem = new ExtendedSystem( simnfo ); + the_extendedsystem->setTargetTemp(the_globals->getTargetTemp()); + the_extendedsystem->setTargetPressure(the_globals->getTargetPressure()); + the_extendedsystem->setQmass(the_globals->getQmass()); + the_extendedsystem->setTauRelax(the_globals->getTauRelax()); + } else if ( !strcasecmp( ensemble, "NVT") ) { + the_extendedsystem = new ExtendedSystem( simnfo ); + the_extendedsystem->setTargetTemp(the_globals->getTargetTemp()); + the_extendedsystem->setQmass(the_globals->getQmass()); + } else if ( !strcasecmp( ensemble, "NVE") ) { + } else { + sprintf( painCave.errMsg, + "SimSetup Warning. Unrecognized Ensemble -> %s, " + "reverting to NVE for this simulation.\n", + ensemble ); + painCave.isFatal = 0; + simError(); + strcpy( ensemble, "NVE" ); + } strcpy( simnfo->ensemble, ensemble ); strcpy( simnfo->mixingRule, the_globals->getMixingRule() ); simnfo->usePBC = the_globals->getPBC(); - - - 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(); + int usesDipoles = 0; + if( !strcmp( force_field, "TraPPE_Ex" ) ){ + the_ff = new TraPPE_ExFF(); + usesDipoles = 1; + } + else if( !strcasecmp( force_field, "LJ" ) ) the_ff = new LJ_FF(); else{ sprintf( painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n", @@ -242,12 +265,13 @@ void SimSetup::createSim( void ){ globalIndex = mpiSim->divideLabor(); - - // set up the local variables int localMol, allMol; int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; + + int* mol2proc = mpiSim->getMolToProcMap(); + int* molCompType = mpiSim->getMolComponentType(); allMol = 0; localMol = 0; @@ -259,8 +283,7 @@ void SimSetup::createSim( void ){ for( j=0; jgetMyMolStart() <= allMol && - allMol <= mpiSim->getMyMolEnd() ){ + if( mol2proc[j] == worldRank ){ local_atoms += comp_stamps[i]->getNAtoms(); local_bonds += comp_stamps[i]->getNBonds(); @@ -279,7 +302,7 @@ void SimSetup::createSim( void ){ if( local_atoms != simnfo->n_atoms ){ sprintf( painCave.errMsg, "SimSetup error: mpiSim's localAtom (%d) and SimSetup's" - " localAtom (%d) are note equal.\n", + " localAtom (%d) are not equal.\n", simnfo->n_atoms, local_atoms ); painCave.isFatal = 1; @@ -304,13 +327,45 @@ void SimSetup::createSim( void ){ Atom::createArrays(simnfo->n_atoms); the_atoms = new Atom*[simnfo->n_atoms]; the_molecules = new Molecule[simnfo->n_mol]; + int molIndex; + // initialize the molecule's stampID's +#ifdef IS_MPI + + + molIndex = 0; + for(i=0; igetTotNmol(); i++){ + + if(mol2proc[i] == worldRank ){ + the_molecules[molIndex].setStampID( molCompType[i] ); + the_molecules[molIndex].setMyIndex( molIndex ); + molIndex++; + } + } + +#else // is_mpi + + molIndex = 0; + for(i=0; in_SRI ){ + Exclude::createArray(simnfo->n_SRI); the_excludes = new Exclude*[simnfo->n_SRI]; + for( int ex=0; exn_SRI; ex++) the_excludes[ex] = new Exclude(ex); simnfo->globalExcludes = new int; - simnfo->n_exclude = tot_SRI; + simnfo->n_exclude = simnfo->n_SRI; } else{ @@ -326,7 +381,7 @@ void SimSetup::createSim( void ){ // set the arrays into the SimInfo object simnfo->atoms = the_atoms; - simnfo->sr_interactions = the_sris; + simnfo->molecules = the_molecules; simnfo->nGlobalExcludes = 0; simnfo->excludes = the_excludes; @@ -383,25 +438,12 @@ void SimSetup::createSim( void ){ the_ff->setSimInfo( simnfo ); - makeAtoms(); + makeMolecules(); simnfo->identArray = new int[simnfo->n_atoms]; for(i=0; in_atoms; i++){ simnfo->identArray[i] = the_atoms[i]->getIdent(); } - if( tot_bonds ){ - makeBonds(); - } - - if( tot_bends ){ - makeBends(); - } - - if( tot_torsions ){ - makeTorsions(); - } - - if (the_globals->getUseRF() ) { simnfo->useReactionField = 1; @@ -443,11 +485,11 @@ void SimSetup::createSim( void ){ } simnfo->dielectric = the_globals->getDielectric(); } else { - if (simnfo->n_dipoles) { + if (usesDipoles) { if( !the_globals->haveECR() ){ sprintf( painCave.errMsg, - "SimSetup Warning: using default value of 1/2 the smallest" + "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; @@ -463,7 +505,7 @@ void SimSetup::createSim( void ){ if( !the_globals->haveEST() ){ sprintf( painCave.errMsg, - "SimSetup Warning: using default value of 5% of the" + "SimSetup Warning: using default value of 5%% of the " "electrostaticCutoffRadius for the " "electrostaticSkinThickness\n" ); @@ -638,15 +680,28 @@ void SimSetup::createSim( void ){ // new AllLong( simnfo ); - if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff ); - if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo, the_ff ); - if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo, the_ff ); - if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff ); + if( !strcmp( force_field, "TraPPE_Ex" ) ){ + new Symplectic(simnfo, the_ff, the_extendedsystem); + std::cerr << "called new Symplecic\n"; + fprintf( stderr, "called new Symplectic. stderr\n" ); + } + else if( !strcmp( force_field, "LJ" ) ){ + new Verlet( *simnfo, the_ff, the_extendedsystem ); + std::cerr << "called new Verlet\n"; + fprintf( stderr, "called new Verlet. stderr\n" ); + } + else { + std::cerr << "I'm a bug.\n"; + fprintf( stderr, "Ima bug. stderr %s\n", force_field); + } +#ifdef IS_MPI + mpiSim->mpiRefresh(); +#endif - // initialize the Fortran - + + simnfo->refreshSim(); if( !strcmp( simnfo->mixingRule, "standard") ){ @@ -683,14 +738,21 @@ void SimSetup::makeMolecules( void ){ BondStamp* currentBond; BendStamp* currentBend; TorsionStamp* currentTorsion; + + bond_pair* theBonds; + bend_set* theBends; + torsion_set* theTorsions; + //init the forceField paramters the_ff->readParams(); - // init the molecules + // init the atoms + double ux, uy, uz, u, uSqr; + atomOffset = 0; excludeOffset = 0; for(i=0; in_mol; i++){ @@ -707,7 +769,7 @@ void SimSetup::makeMolecules( void ){ info.myExcludes = &the_excludes[excludeOffset]; info.myBonds = new Bond*[info.nBonds]; info.myBends = new Bend*[info.nBends]; - info.myTorsions = new Torsions*[info.nTorsions]; + info.myTorsions = new Torsion*[info.nTorsions]; theBonds = new bond_pair[info.nBonds]; theBends = new bend_set[info.nBends]; @@ -717,7 +779,7 @@ void SimSetup::makeMolecules( void ){ for(j=0; jgetAtom( j ); + currentAtom = comp_stamps[stampID]->getAtom( j ); if( currentAtom->haveOrientation() ){ dAtom = new DirectionalAtom(j + atomOffset); @@ -758,8 +820,8 @@ void SimSetup::makeMolecules( void ){ theBonds[j].a = currentBond->getA() + atomOffset; theBonds[j].b = currentBond->getB() + atomOffset; - exI = theBonds[i].a; - exJ = theBonds[i].b; + exI = theBonds[j].a; + exJ = theBonds[j].b; // exclude_I must always be the smaller of the pair if( exI > exJ ){ @@ -775,6 +837,7 @@ void SimSetup::makeMolecules( void ){ the_excludes[j+excludeOffset]->setPair( exI, exJ ); #else // isn't MPI + the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); #endif //is_mpi } @@ -790,7 +853,7 @@ void SimSetup::makeMolecules( void ){ if( currentBend->haveExtras() ){ - extras = current_bend->getExtras(); + extras = currentBend->getExtras(); current_extra = extras; while( current_extra != NULL ){ @@ -812,7 +875,7 @@ void SimSetup::makeMolecules( void ){ default: sprintf( painCave.errMsg, - "SimSetup Error: ghostVectorSource was neiter a " + "SimSetup Error: ghostVectorSource was neither a " "double nor an int.\n" "-->Bend[%d] in %s\n", j, comp_stamps[stampID]->getID() ); @@ -906,368 +969,23 @@ void SimSetup::makeMolecules( void ){ the_molecules[i].initialize( info ); - atomOffset += info.nAtoms; - } - - // clean up the forcefield - - the_ff->cleanMe(); -} - - - -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; igetMyMolStart() <= molIndex && - molIndex <= mpiSim->getMyMolEnd() ){ -#endif // is_mpi - - molStart = index; - nMemb = comp_stamps[i]->getNAtoms(); - for( k=0; kgetNAtoms(); k++ ){ - - 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 ); - - // 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++; - } - } - -#ifdef IS_MPI - for( i=0; igetMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] ); - - delete[] globalIndex; - - mpiSim->mpiRefresh(); -#endif //IS_MPI - - the_ff->initializeAtoms(); -} - -void SimSetup::makeBonds( void ){ - - int i, j, k, index, offset, molIndex, exI, exJ, tempEx; - bond_pair* the_bonds; - BondStamp* current_bond; - - the_bonds = new bond_pair[tot_bonds]; - index = 0; - offset = 0; - molIndex = 0; - - for( i=0; igetMyMolStart() <= 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; - exI = the_bonds[index].a; - exJ = the_bonds[index].b; - // exclude_I must always be the smaller of the pair - if( exI > exJ ){ - tempEx = exI; - exI = exJ; - exJ = tempEx; - } - - -#ifdef IS_MPI - - the_excludes[index*2] = - the_atoms[exI]->getGlobalIndex() + 1; - the_excludes[index*2 + 1] = - the_atoms[exJ]->getGlobalIndex() + 1; - -#else // isn't MPI - - the_excludes[index*2] = exI + 1; - the_excludes[index*2 + 1] = exJ + 1; - // fortran index from 1 (hence the +1 in the indexing) -#endif //is_mpi - - // 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, molIndex, exI, exJ, tempEx; - bend_set* the_bends; - BendStamp* current_bend; - LinkedAssign* extras; - LinkedAssign* current_extra; - - - the_bends = new bend_set[tot_bends]; - index = 0; - offset = 0; - molIndex = 0; - for( i=0; igetMyMolStart() <= molIndex && - molIndex <= mpiSim->getMyMolEnd() ){ -#endif // is_mpi - - 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; - - if( current_bend->haveExtras() ){ - - extras = current_bend->getExtras(); - current_extra = extras; - - while( current_extra != NULL ){ - if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){ - - switch( current_extra->getType() ){ - - case 0: - the_bends[index].ghost = - current_extra->getInt() + offset; - the_bends[index].isGhost = 1; - break; - - case 1: - the_bends[index].ghost = - (int)current_extra->getDouble() + offset; - the_bends[index].isGhost = 1; - break; - - default: - sprintf( painCave.errMsg, - "SimSetup Error: ghostVectorSource was neiter a " - "double nor an int.\n" - "-->Bend[%d] in %s\n", - k, comp_stamps[i]->getID() ); - painCave.isFatal = 1; - simError(); - } - } - - else{ - - sprintf( painCave.errMsg, - "SimSetup Error: unhandled bend assignment:\n" - " -->%s in Bend[%d] in %s\n", - current_extra->getlhs(), - k, comp_stamps[i]->getID() ); - painCave.isFatal = 1; - simError(); - } - - current_extra = current_extra->getNext(); - } - } - - if( !the_bends[index].isGhost ){ - - exI = the_bends[index].a; - exJ = the_bends[index].c; - } - else{ - - exI = the_bends[index].a; - exJ = the_bends[index].b; - } - - // exclude_I must always be the smaller of the pair - if( exI > exJ ){ - tempEx = exI; - exI = exJ; - exJ = tempEx; - } - - -#ifdef IS_MPI - - the_excludes[(index + tot_bonds)*2] = - the_atoms[exI]->getGlobalIndex() + 1; - the_excludes[(index + tot_bonds)*2 + 1] = - the_atoms[exJ]->getGlobalIndex() + 1; - -#else // isn't MPI - - the_excludes[(index + tot_bonds)*2] = exI + 1; - the_excludes[(index + tot_bonds)*2 + 1] = exJ + 1; - // fortran index from 1 (hence the +1 in the indexing) -#endif //is_mpi - - - // increment the index and repeat; - index++; - } - offset += comp_stamps[i]->getNAtoms(); - -#ifdef IS_MPI - } -#endif //is_mpi - - molIndex++; - } + atomOffset += info.nAtoms; + delete[] theBonds; + delete[] theBends; + delete[] theTorsions; } #ifdef IS_MPI - sprintf( checkPointMsg, - "Successfully created the bends list.\n" ); + sprintf( checkPointMsg, "all molecules initialized succesfully" ); MPIcheckPoint(); #endif // is_mpi - - the_ff->initializeBends( the_bends ); -} + // clean up the forcefield + the_ff->calcRcut(); + the_ff->cleanMe(); -void SimSetup::makeTorsions( void ){ - - int i, j, k, index, offset, molIndex, exI, exJ, tempEx; - 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; - the_torsions[index].c = current_torsion->getC() + offset; - the_torsions[index].d = current_torsion->getD() + offset; - - exI = the_torsions[index].a; - exJ = the_torsions[index].d; - - - // exclude_I must always be the smaller of the pair - if( exI > exJ ){ - tempEx = exI; - exI = exJ; - exJ = tempEx; - } - - -#ifdef IS_MPI - - the_excludes[(index + tot_bonds + tot_bends)*2] = - the_atoms[exI]->getGlobalIndex() + 1; - the_excludes[(index + tot_bonds + tot_bends)*2 + 1] = - the_atoms[exJ]->getGlobalIndex() + 1; - -#else // isn't MPI - - the_excludes[(index + tot_bonds + tot_bends)*2] = exI + 1; - the_excludes[(index + tot_bonds + tot_bends)*2 + 1] = exJ + 1; - // fortran indexes from 1 (hence the +1 in the indexing) -#endif //is_mpi - - - // increment the index and repeat; - index++; - } - offset += comp_stamps[i]->getNAtoms(); - -#ifdef IS_MPI - } -#endif //is_mpi - - molIndex++; - } - } - - the_ff->initializeTorsions( the_torsions ); } void SimSetup::initFromBass( void ){