--- trunk/OOPSE/libmdtools/SimSetup.cpp 2003/03/26 20:22:02 407 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2003/03/28 19:30:59 434 @@ -88,9 +88,7 @@ void SimSetup::createSim( void ){ - 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(); + if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF(); else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF(); else{ sprintf( painCave.errMsg, @@ -241,13 +239,14 @@ 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 +258,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 +277,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,29 +302,61 @@ 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] ); + molIndex++; + } + } +#else // is_mpi + + molIndex = 0; + for(i=0; in_SRI ){ - the_sris = new SRI*[simnfo->n_SRI]; - the_excludes = new int[2 * simnfo->n_SRI]; + + std::cerr << "n_SRI = " << simnfo->n_SRI << "\n"; + + 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; } else{ - the_excludes = new int[2]; - the_excludes[0] = 0; - the_excludes[1] = 0; + Exclude::createArray( 1 ); + the_excludes = new Exclude*; + the_excludes[0] = new Exclude(0); + the_excludes[0]->setPair( 0,0 ); simnfo->globalExcludes = new int; simnfo->globalExcludes[0] = 0; - - simnfo->n_exclude = 1; + simnfo->n_exclude = 0; } // 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 +413,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; @@ -638,15 +655,16 @@ 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 ); - +#ifdef IS_MPI + mpiSim->mpiRefresh(); +#endif // initialize the Fortran - + + simnfo->refreshSim(); if( !strcmp( simnfo->mixingRule, "standard") ){ @@ -674,22 +692,32 @@ void SimSetup::makeMolecules( void ){ void SimSetup::makeMolecules( void ){ - int i, j, exI, exJ, tempEx, stampID, atomOffset; + int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; molInit info; DirectionalAtom* dAtom; + LinkedAssign* extras; + LinkedAssign* current_extra; AtomStamp* currentAtom; 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++){ stampID = the_molecules[i].getStampID(); @@ -698,11 +726,13 @@ void SimSetup::makeMolecules( void ){ 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; + 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 Torsions*[info.nTorsions]; + info.myTorsions = new Torsion*[info.nTorsions]; theBonds = new bond_pair[info.nBonds]; theBends = new bend_set[info.nBends]; @@ -712,7 +742,7 @@ void SimSetup::makeMolecules( void ){ for(j=0; jgetAtom( j ); + currentAtom = comp_stamps[stampID]->getAtom( j ); if( currentAtom->haveOrientation() ){ dAtom = new DirectionalAtom(j + atomOffset); @@ -747,7 +777,7 @@ void SimSetup::makeMolecules( void ){ } // make the bonds - for(j=0; jgetBond( j ); theBonds[j].a = currentBond->getA() + atomOffset; @@ -763,388 +793,159 @@ void SimSetup::makeMolecules( void ){ exJ = tempEx; } #ifdef IS_MPI + tempEx = exI; + exI = the_atoms[tempEx]->getGlobalIndex() + 1; + tempEx = exJ; + exJ = the_atoms[tempEx]->getGlobalIndex() + 1; - the_excludes[index*2] = - the_atoms[exI]->getGlobalIndex() + 1; - the_excludes[index*2 + 1] = - the_atoms[exJ]->getGlobalIndex() + 1; - + the_excludes[j+excludeOffset]->setPair( exI, exJ ); #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) - + the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); #endif //is_mpi - - } + } + excludeOffset += info.nBonds; - - - - - - - - - - - - - - -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++ ){ + //make the bends + for(j=0; jgetBend( j ); + theBends[j].a = currentBend->getA() + atomOffset; + theBends[j].b = currentBend->getB() + atomOffset; + theBends[j].c = currentBend->getC() + atomOffset; - current_atom = comp_stamps[i]->getAtom( k ); - if( current_atom->haveOrientation() ){ + if( currentBend->haveExtras() ){ - dAtom = new DirectionalAtom(index); - simnfo->n_oriented++; - the_atoms[index] = dAtom; + extras = currentBend->getExtras(); + current_extra = extras; - 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 ); + while( current_extra != NULL ){ + if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){ + + switch( current_extra->getType() ){ + + case 0: + theBends[j].ghost = + current_extra->getInt() + atomOffset; + theBends[j].isGhost = 1; + break; + + case 1: + theBends[j].ghost = + (int)current_extra->getDouble() + atomOffset; + theBends[j].isGhost = 1; + break; + + default: + sprintf( painCave.errMsg, + "SimSetup Error: ghostVectorSource was neither a " + "double nor an int.\n" + "-->Bend[%d] in %s\n", + j, comp_stamps[stampID]->getID() ); + painCave.isFatal = 1; + simError(); + } } + else{ - the_atoms[index] = new GeneralAtom(index); + + sprintf( painCave.errMsg, + "SimSetup Error: unhandled bend assignment:\n" + " -->%s in Bend[%d] in %s\n", + current_extra->getlhs(), + j, comp_stamps[stampID]->getID() ); + painCave.isFatal = 1; + simError(); } - the_atoms[index]->setType( current_atom->getType() ); - the_atoms[index]->setIndex( index ); - // increment the index and repeat; - index++; + current_extra = current_extra->getNext(); } + } + + if( !theBends[j].isGhost ){ + + exI = theBends[j].a; + exJ = theBends[j].c; + } + else{ - 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 + exI = theBends[j].a; + exJ = theBends[j].b; } -#endif //is_mpi - molIndex++; + // exclude_I must always be the smaller of the pair + if( exI > exJ ){ + tempEx = exI; + exI = exJ; + exJ = tempEx; + } +#ifdef IS_MPI + tempEx = exI; + exI = the_atoms[tempEx]->getGlobalIndex() + 1; + tempEx = exJ; + exJ = the_atoms[tempEx]->getGlobalIndex() + 1; + + the_excludes[j+excludeOffset]->setPair( exI, exJ ); +#else // isn't MPI + the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); +#endif //is_mpi } - } + excludeOffset += info.nBends; -#ifdef IS_MPI - for( i=0; igetMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] ); - - delete[] globalIndex; + for(j=0; jgetTorsion( j ); + theTorsions[j].a = currentTorsion->getA() + atomOffset; + theTorsions[j].b = currentTorsion->getB() + atomOffset; + theTorsions[j].c = currentTorsion->getC() + atomOffset; + theTorsions[j].d = currentTorsion->getD() + atomOffset; + + exI = theTorsions[j].a; + exJ = theTorsions[j].d; - 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; i exJ ){ + tempEx = exI; + exI = exJ; + exJ = tempEx; + } #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; - - 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; - + tempEx = exI; + exI = the_atoms[tempEx]->getGlobalIndex() + 1; + tempEx = exJ; + exJ = the_atoms[tempEx]->getGlobalIndex() + 1; + + the_excludes[j+excludeOffset]->setPair( exI, exJ ); #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) + the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) ); #endif //is_mpi - - // increment the index and repeat; - index++; - } - offset += comp_stamps[i]->getNAtoms(); - -#ifdef IS_MPI - } -#endif //is_mpi - - molIndex++; - } - } + } + excludeOffset += info.nTorsions; - the_ff->initializeBonds( the_bonds ); -} + + // send the arrays off to the forceField for init. -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; - } + 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 ); -#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++; - } + the_molecules[i].initialize( info ); + 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 ); -} - -void SimSetup::makeTorsions( void ){ + // clean up the forcefield + the_ff->calcRcut(); + the_ff->cleanMe(); - 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 ){