--- trunk/OOPSE/libmdtools/SimSetup.cpp 2003/03/26 22:02:36 414 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2003/03/27 20:36:16 424 @@ -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, @@ -248,6 +246,9 @@ void SimSetup::createSim( void ){ 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 +260,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 +279,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,8 +304,36 @@ 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 ){ Exclude::createArray(simnfo->n_SRI); the_excludes = new Exclude*[simnfo->n_SRI]; @@ -326,7 +354,6 @@ void SimSetup::createSim( void ){ // set the arrays into the SimInfo object simnfo->atoms = the_atoms; - simnfo->sr_interactions = the_sris; simnfo->nGlobalExcludes = 0; simnfo->excludes = the_excludes; @@ -383,25 +410,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,8 +652,6 @@ 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 ); @@ -910,366 +922,10 @@ void SimSetup::makeMolecules( void ){ } // clean up the forcefield - + the_ff->calcRcut(); 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++; - } - } - -#ifdef IS_MPI - sprintf( checkPointMsg, - "Successfully created the bends list.\n" ); - MPIcheckPoint(); -#endif // is_mpi - - - the_ff->initializeBends( the_bends ); -} - -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 ){ int i, j, k;