--- branches/mmeineke/OOPSE/libmdtools/SimSetup.cpp 2003/03/21 17:42:12 377 +++ 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 @@ -333,30 +333,6 @@ void SimSetup::createSim( void ){ // get some of the tricky things that may still be in the globals - if( simnfo->n_dipoles ){ - - if( !the_globals->haveRRF() ){ - sprintf( painCave.errMsg, - "SimSetup Error, system has dipoles, but no rRF was set.\n"); - painCave.isFatal = 1; - simError(); - } - if( !the_globals->haveDielectric() ){ - 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(); @@ -426,9 +402,84 @@ void SimSetup::createSim( void ){ } - + if (the_globals->getUseRF() ) { + simnfo->useReactionField = 1; + + if( !the_globals->haveECR() ){ + sprintf( painCave.errMsg, + "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; + simError(); + double smallest; + smallest = simnfo->box_x; + if (simnfo->box_y <= smallest) smallest = simnfo->box_y; + if (simnfo->box_z <= smallest) smallest = simnfo->box_z; + simnfo->ecr = 0.5 * smallest; + } else { + simnfo->ecr = the_globals->getECR(); + } + if( !the_globals->haveEST() ){ + sprintf( painCave.errMsg, + "SimSetup Warning: using default value of 0.05 * the " + "electrostaticCutoffRadius for the electrostaticSkinThickness\n" + ); + painCave.isFatal = 0; + simError(); + simnfo->est = 0.05 * simnfo->ecr; + } else { + simnfo->est = the_globals->getEST(); + } + + if(!the_globals->haveDielectric() ){ + sprintf( painCave.errMsg, + "SimSetup Error: You are trying to use Reaction Field without" + "setting a dielectric constant!\n" + ); + painCave.isFatal = 1; + simError(); + } + simnfo->dielectric = the_globals->getDielectric(); + } else { + if (simnfo->n_dipoles) { + + if( !the_globals->haveECR() ){ + sprintf( painCave.errMsg, + "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; + simError(); + double smallest; + smallest = simnfo->box_x; + if (simnfo->box_y <= smallest) smallest = simnfo->box_y; + if (simnfo->box_z <= smallest) smallest = simnfo->box_z; + simnfo->ecr = 0.5 * smallest; + } else { + simnfo->ecr = the_globals->getECR(); + } + + if( !the_globals->haveEST() ){ + sprintf( painCave.errMsg, + "SimSetup Warning: using default value of 5% of the" + "electrostaticCutoffRadius for the " + "electrostaticSkinThickness\n" + ); + painCave.isFatal = 0; + simError(); + simnfo->est = 0.05 * simnfo->ecr; + } else { + simnfo->est = the_globals->getEST(); + } + } + } +#ifdef IS_MPI + strcpy( checkPointMsg, "electrostatic parameters check out" ); + MPIcheckPoint(); +#endif // is_mpi if( the_globals->haveInitialConfig() ){ @@ -617,8 +668,253 @@ void SimSetup::createSim( void ){ strcpy( checkPointMsg, "Successfully intialized the mixingRule for Fortran." ); 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 ){