--- trunk/OOPSE/libmdtools/SimInfo.cpp 2004/04/28 22:06:29 1139 +++ trunk/OOPSE/libmdtools/SimInfo.cpp 2004/06/02 14:21:54 1218 @@ -42,11 +42,10 @@ SimInfo::SimInfo(){ thermalTime = 0.0; currentTime = 0.0; rCut = 0.0; - ecr = 0.0; - est = 0.0; + rSw = 0.0; haveRcut = 0; - haveEcr = 0; + haveRsw = 0; boxIsInit = 0; resetTime = 1e99; @@ -63,8 +62,11 @@ SimInfo::SimInfo(){ useReactionField = 0; useGB = 0; useEAM = 0; - useMolecularCutoffs = 0; + useSolidThermInt = 0; + useLiquidThermInt = 0; + haveCutoffGroups = false; + excludes = Exclude::Instance(); myConfiguration = new SimState(); @@ -72,6 +74,8 @@ SimInfo::SimInfo(){ has_minimizer = false; the_minimizer =NULL; + ngroup = 0; + wrapMeSimInfo( this ); } @@ -84,7 +88,7 @@ SimInfo::~SimInfo(){ for(i = properties.begin(); i != properties.end(); i++) delete (*i).second; - + } void SimInfo::setBox(double newBox[3]) { @@ -187,25 +191,27 @@ void SimInfo::calcHmatInv( void ) { if( oldOrtho != orthoRhombic ){ - if( orthoRhombic ){ + if( orthoRhombic ) { sprintf( painCave.errMsg, - "OOPSE is switching from the default Non-Orthorhombic\n" + "\n\tOOPSE is switching from the default Non-Orthorhombic\n" "\tto the faster Orthorhombic periodic boundary computations.\n" "\tThis is usually a good thing, but if you wan't the\n" "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n" "\tvariable ( currently set to %G ) smaller.\n", orthoTolerance); + painCave.severity = OOPSE_INFO; simError(); } else { sprintf( painCave.errMsg, - "OOPSE is switching from the faster Orthorhombic to the more\n" + "\n\tOOPSE is switching from the faster Orthorhombic to the more\n" "\tflexible Non-Orthorhombic periodic boundary computations.\n" "\tThis is usually because the box has deformed under\n" "\tNPTf integration. If you wan't to live on the edge with\n" "\tthe Orthorhombic computations, make the orthoBoxTolerance\n" "\tvariable ( currently set to %G ) larger.\n", orthoTolerance); + painCave.severity = OOPSE_WARNING; simError(); } } @@ -437,66 +443,58 @@ void SimInfo::refreshSim(){ //fInfo.SIM_uses_RF = 0; fInfo.SIM_uses_GB = useGB; fInfo.SIM_uses_EAM = useEAM; - fInfo.SIM_uses_molecular_cutoffs = useMolecularCutoffs; n_exclude = excludes->getSize(); excl = excludes->getFortranArray(); - + #ifdef IS_MPI - n_global = mpiSim->getTotAtoms(); + n_global = mpiSim->getNAtomsGlobal(); #else n_global = n_atoms; #endif - + isError = 0; - + + getFortranGroupArrays(this, FglobalGroupMembership, mfact); + //it may not be a good idea to pass the address of first element in vector + //since c++ standard does not require vector to be stored continuously in meomory + //Most of the compilers will organize the memory of vector continuously setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, - &nGlobalExcludes, globalExcludes, molMembershipArray, - &isError ); + &nGlobalExcludes, globalExcludes, molMembershipArray, + &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError); if( isError ){ - + sprintf( painCave.errMsg, - "There was an error setting the simulation information in fortran.\n" ); + "There was an error setting the simulation information in fortran.\n" ); painCave.isFatal = 1; simError(); } - + #ifdef IS_MPI sprintf( checkPointMsg, "succesfully sent the simulation information to fortran.\n"); MPIcheckPoint(); #endif // is_mpi - + this->ndf = this->getNDF(); this->ndfRaw = this->getNDFraw(); this->ndfTrans = this->getNDFtranslational(); } void SimInfo::setDefaultRcut( double theRcut ){ - + haveRcut = 1; rCut = theRcut; - - ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; - - notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); -} - -void SimInfo::setDefaultEcr( double theEcr ){ - - haveEcr = 1; - ecr = theEcr; + rList = rCut + 1.0; - ( rCut > ecr )? rList = rCut + 1.0: rList = ecr + 1.0; - - notifyFortranCutOffs( &rCut, &rList, &ecr, &est ); + notifyFortranCutOffs( &rCut, &rSw, &rList ); } -void SimInfo::setDefaultEcr( double theEcr, double theEst ){ +void SimInfo::setDefaultRcut( double theRcut, double theRsw ){ - est = theEst; - setDefaultEcr( theEcr ); + rSw = theRsw; + setDefaultRcut( theRcut ); } @@ -508,8 +506,8 @@ void SimInfo::checkCutOffs( void ){ if( rCut > maxCutoff ){ sprintf( painCave.errMsg, - "LJrcut is too large for the current periodic box.\n" - "\tCurrent Value of LJrcut = %G at time %G\n " + "\n\tcutoffRadius is too large for the current periodic box.\n" + "\tCurrent Value of cutoffRadius = %G at time %G\n " "\tThis is larger than half of at least one of the\n" "\tperiodic box vectors. Right now, the Box matrix is:\n" "\n" @@ -520,35 +518,16 @@ void SimInfo::checkCutOffs( void ){ Hmat[0][0], Hmat[0][1], Hmat[0][2], Hmat[1][0], Hmat[1][1], Hmat[1][2], Hmat[2][0], Hmat[2][1], Hmat[2][2]); + painCave.severity = OOPSE_ERROR; painCave.isFatal = 1; simError(); - } - - if( haveEcr ){ - if( ecr > maxCutoff ){ - sprintf( painCave.errMsg, - "electrostaticCutoffRadius is too large for the current\n" - "\tperiodic box.\n\n" - "\tCurrent Value of ECR = %G at time %G\n " - "\tThis is larger than half of at least one of the\n" - "\tperiodic box vectors. Right now, the Box matrix is:\n" - "\n" - "\t[ %G %G %G ]\n" - "\t[ %G %G %G ]\n" - "\t[ %G %G %G ]\n", - ecr, currentTime, - Hmat[0][0], Hmat[0][1], Hmat[0][2], - Hmat[1][0], Hmat[1][1], Hmat[1][2], - Hmat[2][0], Hmat[2][1], Hmat[2][2]); - painCave.isFatal = 1; - simError(); - } - } + } } else { // initialize this stuff before using it, OK? sprintf( painCave.errMsg, - "Trying to check cutoffs without a box.\n" + "\n\tTrying to check cutoffs without a box.\n" "\tOOPSE should have better programmers than that.\n" ); + painCave.severity = OOPSE_ERROR; painCave.isFatal = 1; simError(); } @@ -591,3 +570,54 @@ GenericData* SimInfo::getProperty(const string& propNa return NULL; } + +void SimInfo::getFortranGroupArrays(SimInfo* info, + vector& FglobalGroupMembership, + vector& mfact){ + + Molecule* myMols; + Atom** myAtoms; + int numAtom; + double mtot; + int numMol; + int numCutoffGroups; + CutoffGroup* myCutoffGroup; + vector::iterator iterCutoff; + Atom* cutoffAtom; + vector::iterator iterAtom; + int atomIndex; + double totalMass; + + mfact.clear(); + FglobalGroupMembership.clear(); + + + // Fix the silly fortran indexing problem +#ifdef IS_MPI + numAtom = mpiSim->getNAtomsGlobal(); +#else + numAtom = n_atoms; +#endif + for (int i = 0; i < numAtom; i++) + FglobalGroupMembership.push_back(globalGroupMembership[i] + 1); + + + myMols = info->molecules; + numMol = info->n_mol; + for(int i = 0; i < numMol; i++){ + numCutoffGroups = myMols[i].getNCutoffGroups(); + for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); + myCutoffGroup != NULL; + myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){ + + totalMass = myCutoffGroup->getMass(); + + for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); + cutoffAtom != NULL; + cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){ + mfact.push_back(cutoffAtom->getMass()/totalMass); + } + } + } + +}