[ViewVC] Diff of: group/branches/new_design/OOPSE-3.0/src/brains/SimCreator.cpp

Comparing branches/new_design/OOPSE-3.0/src/brains/SimCreator.cpp (file contents):
Revision 1703 by tim, Wed Nov 3 19:56:02 2004 UTC vs.
Revision 1733 by tim, Fri Nov 12 06:19:04 2004 UTC

+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ */
-–
-–
+ /**
-–
+  * @file SimCreator.cpp
-–
+  * @author tlin
-–
+  * @date 11/03/2004
-–
+  * @time 13:51am
-–
+  * @version 1.0
-–
+  */
-<
+#include "brains/SimCreator.hpp"
->
+/**
->
+ * @file SimCreator.cpp
->
+ * @author tlin
->
+ * @date 11/03/2004
->
+ * @time 13:51am
->
+ * @version 1.0
->
+ */
->
->
+#include "brains/SimCreator.hpp"
->
+namespace oopse {
-–
+void SimCreator::parseFile(const string & mdfile) {
-+
+void SimSetup::parseFile(const std::string& mdFileName, MakeStamps *stamps, Globals *globals) {
-+
-+
+#ifdef IS_MPI
-+
-+
+    if (worldRank == 0) {
-+
+#endif // is_mpi
-+
-+
+        globals->initalize();
-+
+        set_interface_stamps(stamps, globals);
-+
-+
+#ifdef IS_MPI
-+
-+
+        mpiEventInit();
-+
-+
+#endif
-+
-+
+        yacc_BASS(mdFileName.c_str());
-+
-+
+#ifdef IS_MPI
-+
-+
+        throwMPIEvent(NULL);
-+
+    } else {
-+
+        set_interface_stamps(stamps, globals);
-+
+        mpiEventInit();
-+
+        MPIcheckPoint();
-+
+        mpiEventLoop();
-+
+    }
-+
-+
+#endif
-+
-+
+SimInfo*  SimCreator::createSim(const std::string & mdFileName) {
-+
-+
+    mdFileName_ = mdFileName;
-+
-+
+    MakeStamps * stamps = new MakeStamps();
-+
-+
+    Globals * globals = new Globals();
-+
-+
+    //parse meta-data file
-+
+    parseFile(mdFileName_, stamps, globals);
-+
-+
+    //create the force field
-+
+    ForceFiled * ff = ForceFieldFactory::getInstance()->createObject(
-+
+                          globals->getForceField());
-+
-+
+    if (ff == NULL) {
-+
+        sprintf(painCave.errMsg, "ForceFiled Factory can not create %s force field\n",
-+
+                globals->getForceField());
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+    }
-+
-+
+    //extract the molecule stamps
-+
+    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
-+
+    compList(stamps, globals, moleculeStampPairs);
-+
-+
+    //create SimInfo
-+
+    SimInfo * info = new SimInfo(moleculeStampPairs, ff, globals);
-+
-+
+    //gather parameters (SimCreator only retrieves part of the parameters)
-+
+    gatherParameters(info);
-+
-+
+    //divide the molecules and determine the global index of molecules
-+
+    divideMolecules(info);
-+
-+
+    //create the molecules
-+
+    createMolecules(info);
-+
-+
+    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
-+
+    //global index will never change again). Local indices of atoms and rigidbodies are already set by
-+
+    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
-+
+    setGlobalIndices(info);
-+
-+
+    //allocate memory for DataStorage(circular reference, need to break it)
-+
+    info->setSnapshotManager(new SimSnapshotManager(info));
-+
-+
+    //initialize fortran -- setup the cutoff
-+
+    initFortran(info);
-+
-+
+    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
-+
+    //eta, chi for NPT integrator)
-+
+    loadCoordinates(info);
-+
+    return info;
-+
+}
-+
-+
+void SimCreator::gatherParameters(SimInfo *info) {
-+
-+
+    //setup seed for random number generator
-+
+    int seedValue;
-+
+    Globals * globals = info->getGlobals();
-+
-+
+    if (globals->haveSeed()) {
-+
+        seedValue = globals->getSeed();
-+
-+
+        if (seedValue / 1000000000 == 0) {
-+
+            sprintf(painCave.errMsg,
-+
+                    "Seed for sprng library should contain at least 9 digits\n"
-+
+                        "OOPSE will generate a seed for user\n");
-+
-+
+            painCave.isFatal = 0;
-+
+            simError();
-+
-+
+            //using seed generated by system instead of invalid seed set by user
-+
-+
+#ifndef IS_MPI
-+
-+
+            seedValue = make_sprng_seed();
-+
-+
+#else
-+
-+
+            if (worldRank == 0) {
-+
+                seedValue = make_sprng_seed();
-+
+            }
-+
-+
+            MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
-+
-+
+#endif
-+
-+
+        } //end if (seedValue /1000000000 == 0)
-+
+    } else {
-+
-+
+#ifndef IS_MPI
-+
-+
+        seedValue = make_sprng_seed();
-+
-+
+#else
-+
-+
+        if (worldRank == 0) {
-+
+            seedValue = make_sprng_seed();
-+
+        }
-+
-+
+        MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
-+
-+
+#endif
-+
-+
+    } //end of globals->haveSeed()
-+
-+
+    info->setSeed(seedValue);
-+
-+
-+
+    //figure out the ouput file names
-+
+    std::string prefix;
-+
-+
+#ifdef IS_MPI
-+
-+
+    if (worldRank == 0) {
-+
+#endif // is_mpi
-+
-+
+        if (globals->haveFinalConfig()) {
-+
+            prefix = StringUtils::getPrefix(globals->getFinalConfig());
-+
+        } else {
-+
+            prefix = StringUtils::getPrefix(mdfile_);
-+
+        }
-+
-+
+        info->setFinalConfigFileName(prefix + ".eor");
-+
+        info->setDumpFileName(prefix + ".dump");
-+
+        info->setStatFileName(prefix + ".stat");
-+
-+
+#ifdef IS_MPI
-+
-+
+    }
-+
-+
+#endif
-+
-+
+}
-+
-+
+#ifdef IS_MPI
-+
-+
+void SimCreator::divideMolecules(SimInfo *info) {
-+
+    double numerator;
-+
+    double denominator;
-+
+    double precast;
-+
+    double x;
-+
+    double y;
-+
+    double a;
-+
+    int old_atoms;
-+
+    int add_atoms;
-+
+    int new_atoms;
-+
+    int nTarget;
-+
+    int done;
-+
+    int i;
-+
+    int j;
-+
+    int loops;
-+
+    int which_proc;
-+
+    int nProcessors;
-+
+    int nGlobalMols;
-+
+    std::vector < int > atomsPerProc;
-+
+    randomSPRNG myRandom(info->getSeed());
-+
+    int * molToProcMap = info->getMolToProcMapPointer();
-+
-+
+    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
-+
+    nGlobalMols = info->getNGlobalMolecules();
-+
-+
+    if (nProcessors > nGlobalMols) {
-+
+        sprintf(painCave.errMsg,
-+
+                "nProcessors (%d) > nMol (%d)\n"
-+
+                    "\tThe number of processors is larger than\n"
-+
+                    "\tthe number of molecules.  This will not result in a \n"
-+
+                    "\tusable division of atoms for force decomposition.\n"
-+
+                    "\tEither try a smaller number of processors, or run the\n"
-+
+                    "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
-+
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+    }
-+
-+
+    a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
-+
-+
+    //initialize atomsPerProc
-+
+    atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
-+
-+
+    for(i = 0; i < nGlobalMols; i++) {
-+
+        // default to an error condition:
-+
+        molToProcMap[i] = -1;
-+
+    }
-+
-+
+    if (worldRank == 0) {
-+
+        numerator = info->getNGlobalAtoms();
-+
+        denominator = nProcessors;
-+
+        precast = numerator / denominator;
-+
+        nTarget = (int)(precast + 0.5);
-+
-+
+        for(i = 0; i < nGlobalMols; i++) {
-+
+            done = 0;
-+
+            loops = 0;
-+
-+
+            while (!done) {
-+
+                loops++;
-+
-+
+                // Pick a processor at random
-+
-+
+                which_proc = (int) (myRandom.getRandom() * nProcessors);
-+
-+
+                //get the molecule stamp first
-+
+                int stampId = info->getMoleculeStampId(i);
-+
+                MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
-+
-+
+                // How many atoms does this processor have so far?
-+
+                old_atoms = atomsPerProc[which_proc];
-+
+                add_atoms = moleculeStamp->getNAtoms();
-+
+                new_atoms = old_atoms + add_atoms;
-+
-+
+                // If we've been through this loop too many times, we need
-+
+                // to just give up and assign the molecule to this processor
-+
+                // and be done with it.
-+
-+
+                if (loops > 100) {
-+
+                    sprintf(painCave.errMsg,
-+
+                            "I've tried 100 times to assign molecule %d to a "
-+
+                                " processor, but can't find a good spot.\n"
-+
+                                "I'm assigning it at random to processor %d.\n",
-+
+                            i, which_proc);
-+
-+
+                    painCave.isFatal = 0;
-+
+                    simError();
-+
-+
+                    molToProcMap[i] = which_proc;
-+
+                    atomsPerProc[which_proc] += add_atoms;
-+
-+
+                    done = 1;
-+
+                    continue;
-+
+                }
-+
-+
+                // If we can add this molecule to this processor without sending
-+
+                // it above nTarget, then go ahead and do it:
-+
-+
+                if (new_atoms <= nTarget) {
-+
+                    molToProcMap[i] = which_proc;
-+
+                    atomsPerProc[which_proc] += add_atoms;
-+
-+
+                    done = 1;
-+
+                    continue;
-+
+                }
-+
-+
+                // The only situation left is when new_atoms > nTarget.  We
-+
+                // want to accept this with some probability that dies off the
-+
+                // farther we are from nTarget
-+
-+
+                // roughly:  x = new_atoms - nTarget
-+
+                //           Pacc(x) = exp(- a * x)
-+
+                // where a = penalty / (average atoms per molecule)
-+
-+
+                x = (double)(new_atoms - nTarget);
-+
+                y = myRandom.getRandom();
-+
-+
+                if (y < exp(- a * x)) {
-+
+                    molToProcMap[i] = which_proc;
-+
+                    atomsPerProc[which_proc] += add_atoms;
-+
-+
+                    done = 1;
-+
+                    continue;
-+
+                } else {
-+
+                    continue;
-+
+                }
-+
+            }
-+
+        }
-+
-+
+        // Spray out this nonsense to all other processors:
-+
-+
+        MPI_Bcast(molToProcMap, nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
-+
+    } else {
-+
-+
+        // Listen to your marching orders from processor 0:
-+
-+
+        MPI_Bcast(molToProcMap, nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
-+
+    }
-+
-+
+    sprintf(checkPointMsg,
-+
+            "Successfully divided the molecules among the processors.\n");
-+
+    MPIcheckPoint();
-+
+}
-+
-+
+#endif
-+
-+
+void SimCreator::createMolecules(SimInfo *info) {
-+
+    MoleculeCreator molCreator;
-+
+    int stampId;
-+
-+
+    for(int i = 0; i < info->getNGlobalMolecules(); i++) {
-+
-+
+#ifdef IS_MPI
-+
-+
+        if (info->getMolToProc(i) == worldRank) {
-+
+#endif
-+
-+
+            stampId = info->getMoleculeStampId(i);
-+
+            Molecule * mol = molCreator.createMolecule(info->getForceField(),
-+
+                                                info->getMoleculeStamp(stampId), stampId, i);
-+
-+
+            info->addMolecule(mol);
-+
-+
+#ifdef IS_MPI
-+
-+
+        }
-+
-+
+#endif
-+
-+
+    } //end for(int i=0)
-+
+}
-+
-+
+void SimSetup::compList(MakeStamps *stamps, Globals* globals,
-+
+                        std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
-+
+    int i;
-+
+    char * id;
-+
+    MoleculeStamp * currentStamp;
-+
+    Component * the_components = globals->getComponents();
-+
+    int n_components = globals->getNComponents();
-+
-+
+    if (!globals->haveNMol()) {
-+
+        // we don't have the total number of molecules, so we assume it is
-+
+        // given in each component
-+
-+
+        for(i = 0; i < n_components; i++) {
-+
+            if (!the_components[i]->haveNMol()) {
-+
+                // we have a problem
-+
+                sprintf(painCave.errMsg,
-+
+                        "SimSetup Error. No global NMol or component NMol given.\n"
-+
+                            "\tCannot calculate the number of atoms.\n");
-+
-+
+                painCave.isFatal = 1;
-+
+                simError();
-+
+            }
-+
-+
+            id = the_components[i]->getType();
-+
+            currentStamp = stamps->extractMolStamp(id);
-+
-+
+            if (currentStamp == NULL) {
-+
+                sprintf(painCave.errMsg,
-+
+                        "SimSetup error: Component \"%s\" was not found in the "
-+
+                            "list of declared molecules\n", id);
-+
-+
+                painCave.isFatal = 1;
-+
+                simError();
-+
+            }
-+
-+
+            moleculeStampPairs.push_back(
-+
+                make_pair(currentStamp, the_components[i]->getNMol));
-+
+        } //end for (i = 0; i < n_components; i++)
-+
+    } else {
-+
+        sprintf(painCave.errMsg, "SimSetup error.\n"
-+
+                                     "\tSorry, the ability to specify total"
-+
+                                     " nMols and then give molfractions in the components\n"
-+
+                                     "\tis not currently supported."
-+
+                                     " Please give nMol in the components.\n");
-+
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+    }
-+
-+
+#ifdef IS_MPI
-+
-+
+    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
-+
+    MPIcheckPoint();
-+
-+
+#endif // is_mpi
-+
-+
+}
-+
-+
+void SimCreator::setGlobalIndex(SimInfo *info) {
-+
+    typename SimInfo::MoleculeIterator mi;
-+
+    typename Molecule::AtomIterator ai;
-+
+    typename Molecule::RigidBodyIterator ri;
-+
+    typename Molecule::CutoffGroupIterator ci;
-+
+    Molecule * mol;
-+
+    Atom * atom;
-+
+    RigidBody * rb;
-+
+    CutoffGroup * cg;
-+
+    int beginAtomIndex;
-+
+    int beginRigidBodyIndex;
-+
+    int beginCutoffGroupIndex;
-+
-+
+#ifndef IS_MPI
-+
-+
+    beginAtomIndex = 0;
-+
+    beginRigidBodyIndex = 0;
-+
+    beginCutoffGroupIndex = 0;
-+
-+
+#else
-+
-+
+    int nproc;
-+
+    int myNode;
-+
-+
+    myNode = worldRank;
-+
+    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
-+
-+
+    std::vector < int > tmpAtomsInProc(nproc, 0);
-+
+    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
-+
+    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
-+
+    std::vector < int > NumAtomsInProc(nproc, 0);
-+
+    std::vector < int > NumRigidBodiesInProc(nproc, 0);
-+
+    std::vector < int > NumCutoffGroupsInProc(nproc, 0);
-+
-+
+    tmpAtomsInProc[myNode] = info->getNAtoms();
-+
+    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
-+
+    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
-+
-+
+    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
-+
+    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
-+
+                  MPI_SUM, MPI_COMM_WORLD);
-+
+    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
-+
+                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-+
+    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
-+
+                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-+
-+
+    beginAtomIndex = 0;
-+
+    beginRigidBodyIndex = 0;
-+
+    beginCutoffGroupIndex = 0;
-+
-+
+    for(int i = 0; i < nproc; i++) {
-+
+        beginAtomIndex += NumAtomsInProc[i];
-+
+        beginRigidBodyIndex += NumRigidBodiesInProc[i];
-+
+        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
-+
+    }
-+
-+
+#endif
-+
-+
+    for(mol = info->beginMolecule(mi); mol != NULL;
-+
+        mol = info->nextMolecule(mi)) {
-+
-+
+        //local index(index in DataStorge) of atom is important
-+
+        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
-+
+            atom->setGlobalIndex(beginAtomIndex++);
-+
+        }
-+
-+
+        for(rb = mol->beginRigidBody(ri); rb != NULL;
-+
+            rb = mol->nextRigidBody(ri)) {
-+
+            rb->setGlobalIndex(beginRigidBodyIndex++);
-+
+        }
-+
-+
+        //local index of cutoff group is trivial, it only depends on the order of travesing
-+
+        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
-+
+            cg = mol->nextCutoffGroup(ci)) {
-+
+            cg->setGlobalIndex(beginCutoffGroupIndex++);
-+
+        }
-+
+    }
-+
-+
+    std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
-+
+    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
-+
+        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
-+
-+
+            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
-+
+                globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
-+
+            }
-+
-+
+        }
-+
+    }
-+
-+
+    // Since the globalGroupMembership has been zero filled and we've only
-+
+    // poked values into the atoms we know, we can do an Allreduce
-+
+    // to get the full globalGroupMembership array (We think).
-+
+    // This would be prettier if we could use MPI_IN_PLACE like the MPI-2
-+
+    // docs said we could.
-+
-+
+    MPI_Allreduce(&globalGroupMembership[0],
-+
+                  info->getGlobalGroupMembershipPointer(),
-+
+                  info->getNGlobalAtoms(),
-+
+                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-+
-+
+    //fill molMembership
-+
+    std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
-+
-+
+    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
-+
-+
+        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
-+
+            globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
-+
+        }
-+
-+
+    MPI_Allreduce(&globalMolMembership[0],
-+
+                  info->getGlobalMolMembershipPointer(),
-+
+                  info->getNGlobalAtoms(),
-+
+                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-+
-+
+}
-+
-+
-+
+void SimCreator::initFortran(SimInfo* info) {
-+
+    info->update();
-+
-+
+    std::vector<AtomType*> atomTypes;
-+
+    std::vector<AtomType*>::iterator i;
-+
+    std::vector<double> cutoffRadius;
-+
+    ForceField* ff;
-+
-+
+    //get the unique atom types
-+
+    atomTypes = info->getUniqueAtomTypes();
-+
-+
+    ff = info->getForceField();
-+
-+
+    //query the max cutoff radius among these atom types
-+
+    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
-+
+        cutoffRadius.push_back(ff->(*i));
-+
+    }
-+
-+
+    double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end());
-+
-+
+#ifdef IS_MPI
-+
+    //pick the max cutoff radius among the processors
-+
+#endif
-+
-+
-+
-+
+}
-+
-+
+void SimCreator::loadCoordinates(SimInfo* info) {
-+
+    Globals* globals;
-+
+    globals = info->getGlobals();
-+
-+
+    if (!globals->haveInitialConfig()) {
-+
+        sprintf(painCave.errMsg,
-+
+                "Cannot intialize a simulation without an initial configuration file.\n");
-+
+        painCave.isFatal = 1;;
-+
+        simError();
-+
+    }
-+
-+
+    DumpReader reader(globals->getInitialConfig());
-+
+    int nframes = reader->getNframes();
-+
-+
+    if (nframes > 0) {
-+
+        reader.readFrame(info, nframes - 1);
-+
+    } else {
-+
+    //invalid initial coordinate file
-+
+        sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
-+
+                globals->getInitialConfig());
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+    }
-+
-+
+}
-+
+} //end namespace oopse

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing branches/new_design/OOPSE-3.0/src/brains/SimCreator.cpp (file contents): Revision 1703 by tim, Wed Nov 3 19:56:02 2004 UTC vs. Revision 1733 by tim, Fri Nov 12 06:19:04 2004 UTC

Diff Legend

Comparing branches/new_design/OOPSE-3.0/src/brains/SimCreator.cpp (file contents):
Revision 1703 by tim, Wed Nov 3 19:56:02 2004 UTC vs.
Revision 1733 by tim, Fri Nov 12 06:19:04 2004 UTC