--- trunk/OOPSE/libmdtools/SimSetup.cpp 2003/03/24 21:55:34 394 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2003/03/26 22:02:36 414 @@ -307,20 +307,20 @@ void SimSetup::createSim( void ){ if( simnfo->n_SRI ){ - the_sris = new SRI*[simnfo->n_SRI]; - the_excludes = new int[2 * simnfo->n_SRI]; + Exclude::createArray(simnfo->n_SRI); + the_excludes = new Exclude*[simnfo->n_SRI]; 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 @@ -670,7 +670,252 @@ void SimSetup::createSim( void ){ MPIcheckPoint(); #endif // is_mpi } + + +void SimSetup::makeMolecules( void ){ + + 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; + + //init the forceField paramters + + the_ff->readParams(); + + + // init the molecules + atomOffset = 0; + excludeOffset = 0; + for(i=0; in_mol; i++){ + + 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; + + 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]; + + theBonds = new bond_pair[info.nBonds]; + theBends = new bend_set[info.nBends]; + theTorsions = new torsion_set[info.nTorsions]; + + // make the Atoms + + for(j=0; jgetAtom( j ); + if( currentAtom->haveOrientation() ){ + + dAtom = new DirectionalAtom(j + atomOffset); + simnfo->n_oriented++; + info.myAtoms[j] = dAtom; + + ux = currentAtom->getOrntX(); + uy = currentAtom->getOrntY(); + uz = currentAtom->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{ + info.myAtoms[j] = new GeneralAtom(j + atomOffset); + } + info.myAtoms[j]->setType( currentAtom->getType() ); + +#ifdef IS_MPI + + info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] ); + +#endif // is_mpi + } + + // make the bonds + for(j=0; jgetBond( j ); + theBonds[j].a = currentBond->getA() + atomOffset; + theBonds[j].b = currentBond->getB() + atomOffset; + + exI = theBonds[i].a; + exJ = theBonds[i].b; + + // 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.nBonds; + + //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; + + if( currentBend->haveExtras() ){ + + extras = current_bend->getExtras(); + current_extra = extras; + + 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 neiter a " + "double nor an int.\n" + "-->Bend[%d] in %s\n", + j, comp_stamps[stampID]->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(), + j, comp_stamps[stampID]->getID() ); + painCave.isFatal = 1; + simError(); + } + + current_extra = current_extra->getNext(); + } + } + + if( !theBends[j].isGhost ){ + + exI = theBends[j].a; + exJ = theBends[j].c; + } + else{ + + exI = theBends[j].a; + exJ = theBends[j].b; + } + + // 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; + + 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; + + // 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.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_molecules[i].initialize( info ); + atomOffset += info.nAtoms; + } + + // clean up the forcefield + + the_ff->cleanMe(); +} + + + void SimSetup::makeAtoms( void ){ int i, j, k, index;