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

Comparing branches/new_design/OOPSE-2.0/src/brains/SimCreator.cpp (file contents):
Revision 1712 by tim, Thu Nov 4 20:55:01 2004 UTC vs.
Revision 1727 by tim, Thu Nov 11 16:41:58 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 SimSetup::parseFile(char* mdfile,  MakeStamps* stamps, Globals* globals){
-<
+SimModel* SimCreator::createSim(const std::string& mdfile) {
-<
+    MakeStamps* stamps;
->
+SimInfo* SimCreator::createSim(const std::string& mdfile) {
->
+    MakeStamps* stamps = new MakeStamps();
+    Globals* globals = new Globals();
+    parseFile(mdfile, stamps, globals);
+    //create the force field
-<
+    ForceFiled* ff = ForceFiledFactory::getInstance()->createForceField(globals->getForceField());
->
+    ForceFiled* ff = ForceFieldFactory::getInstance()->createObject(globals->getForceField());
-<
+    //create SimModel
-<
+    model = new SimInfo();
->
+    //create SimInfo
->
+    SimInfo* info = new SimInfo();
->
+    info->setGlobals(globals);
->
+    info->setForceField(ff);
->
->
+    //extract the molecule stamps and add them into SimInfo
->
+    compList(stamps, info);
-<
+    //gather parameters
-<
+    gatherParameters(model);
->
+    //gather parameters (SimCreator only retrieves the parameters which will be used to create
->
+    // the simulation)
->
+    gatherParameters(info);
-–
-–
+    //divide the molecules and determine the global index of molecules
-<
->
+    //create the molecules
-+
+    createMolecules(info);
-<
+    //create atoms, bonds, bend, torsions, rigidbodies
-<
->
+    //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)
-<
+    model->setSnapshotManager(new SimSnapshotManager(model);
->
+    info->setSnapshotManager(new SimSnapshotManager(info);
-<
+    //load initial coordinates
-<
+    DumpReader reader();
-<
+    reader.readFrame(model);
->
+    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
->
+    //eta, chi for NPT integrator)
->
+    DumpReader reader(info->getInitFilename());
->
+    int nframes = reader->getNframes();
->
->
+    if (nframes > 0) {
->
+        reader.readFrame(info, nframes - 1);
->
+    } else {
->
+        //invalid initial coordinate file
->
->
+    }
->
+    //initialize fortran
-<
+    return model;
->
->
+    return info;
-<
+void SimCreator::gatherParameters(SimModel* model) {
-<
+    model->addProperty(new StringGenericData("Ensemble", globals->getForceFiled()));
-<
+    model->addProperty(new DoubleGenericData("dt"), globals->getDt());
->
+void SimCreator::gatherParameters(SimInfo* info) {
->
+    //model->addProperty(new StringGenericData("Ensemble", globals->getForceFiled()));
->
+    //model->addProperty(new DoubleGenericData("dt"), globals->getDt());
->
->
+    //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);
->
->
+    //
->
+    std::string prefix;
->
+#ifdef IS_MPI
->
+    if (worldRank == 0){
->
+#endif // is_mpi
->
->
+        if(globals->haveFinalConfig()) {
->
+            prefix = getPrefix(globals->getFinalConfig());
->
+        } else {
->
+            prefix = getPrefix(mdfile_);
->
+        }
->
->
+        info->setFinalConfFileName(prefix + ".eor");
->
+        info->setDumpFileName(prefix + ".dump");
->
+        info->setStatFileName(prefix + ".stat");
->
->
+#ifdef IS_MPI
->
+    }
->
+#endif
->
->
->
+#ifdef IS_MPI
-<
+void SimCreator::mpiMolDivide() {
->
+void SimCreator::divideMolecules(SimInfo* info){
->
+    int nComponents;
->
+    MoleculeStamp ** compStamps;
->
+    randomSPRNG * myRandom;
->
+    int * componentsNmol;
->
+    int * AtomsPerProc;
-+
+    double numerator;
-+
+    double denominator;
-+
+    double precast;
-+
+    double x;
-+
+    double y;
-+
+    double a;
-+
+    int old_atoms;
-+
+    int add_atoms;
-+
+    int new_atoms;
-+
+    int nTarget;
-+
+    int molIndex;
-+
+    int atomIndex;
-+
+    int done;
-+
+    int i;
-+
+    int j;
-+
+    int loops;
-+
+    int which_proc;
-+
+    int nmol_global,
-+
+    nmol_local;
-+
+    int natoms_global,
-+
-+
+    int baseSeed = info->getSeed();
-+
+    CutoffGroupStamp * cg;
-+
-+
+    nComponents = info->nComponents;
-+
+    compStamps = info->compStamps;
-+
+    componentsNmol = info->componentsNmol;
-+
+    AtomsPerProc = new int[parallelData->nProcessors];
-+
-+
+    parallelData->nMolGlobal = info->n_mol;
-+
-+
+    if (parallelData->nProcessors > parallelData->nMolGlobal) {
-+
+        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",
-+
+                parallelData->nProcessors,
-+
+                parallelData->nMolGlobal);
-+
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+    }
-+
-+
+    myRandom = new randomSPRNG(baseSeed);
-+
-+
+    a = 3.0 * parallelData->nMolGlobal / parallelData->nAtomsGlobal;
-+
-+
+    // Initialize things that we'll send out later:
-+
+    for( i = 0; i < parallelData->nProcessors; i++ ) {
-+
+        AtomsPerProc[i] = 0;
-+
+    }
-+
-+
+    for( i = 0; i < parallelData->nMolGlobal; i++ ) {
-+
+        // default to an error condition:
-+
+        MolToProcMap[i] = -1;
-+
+    }
-+
-+
+    if (parallelData->myNode == 0) {
-+
+        numerator = info->n_atoms;
-+
+        denominator = parallelData->nProcessors;
-+
+        precast = numerator / denominator;
-+
+        nTarget = (int)(precast + 0.5);
-+
-+
+        // Build the array of molecule component types first
-+
+        molIndex = 0;
-+
-+
+        for( i = 0; i < nComponents; i++ ) {
-+
+            for( j = 0; j < componentsNmol[i]; j++ ) {
-+
+                molIndex++;
-+
+            }
-+
+        }
-+
-+
+        for( i = 0; i < molIndex; i++ ) {
-+
+            done = 0;
-+
+            loops = 0;
-+
-+
+            while (!done) {
-+
+                loops++;
-+
-+
+                // Pick a processor at random
-+
-+
+                which_proc
-+
-+
+                = (int) (myRandom->getRandom() * parallelData->nProcessors);
-+
-+
+                // How many atoms does this processor have?
-+
-+
+                old_atoms = AtomsPerProc[which_proc];
-+
+                add_atoms = compStamps[MolComponentType[i]]->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, parallelData->nMolGlobal, MPI_INT, 0,
-+
+                  MPI_COMM_WORLD);
-+
+    } else {
-+
-+
+    // Listen to your marching orders from processor 0:
-+
-+
+        MPI_Bcast(MolToProcMap, parallelData->nMolGlobal, MPI_INT, 0,
-+
+                  MPI_COMM_WORLD); }
-+
-+
+    // Let's all check for sanity:
-+
-+
+    nmol_local = 0;
-+
-+
+    for( i = 0; i < parallelData->nMolGlobal; i++ ) {
-+
+        if (MolToProcMap[i] == parallelData->myNode) {
-+
+            nmol_local++;
-+
+        }
-+
+    }
-+
-+
+    MPI_Allreduce(&nmol_local, &nmol_global, 1, MPI_INT, MPI_SUM,
-+
+                  MPI_COMM_WORLD);
-+
-+
+    if (nmol_global != info->n_mol) {
-+
+        sprintf(painCave.errMsg,
-+
+                "The sum of all nmol_local, %d, did not equal the "
-+
+                    "total number of molecules, %d.\n",
-+
+                nmol_global,
-+
+                info->n_mol);
-+
-+
+        painCave.isFatal = 1;
-+
+        simError();
-+
+    }
-+
-+
+    sprintf(checkPointMsg,
-+
+            "Successfully divided the molecules among the processors.\n");
-+
+    MPIcheckPoint();
-+
-+
+#endif
-+
+Molecule* SimCreator::createMolecules(SimInfo* info) {
-+
+    MoleculeCreator molCreator;
-+
+    int stampId;
-+
-+
-+
+    for (int i = 0; i < info->getNGlobalMolecules(); i++){
-+
+        if (info->getMolToProc(i) == worldRank){
-+
-+
+            stampId = info->getMoleculeStampId(i);
-+
+            Molecule* mol = molCreator.createMolecule(info->getForceField(),
-+
+                                                                                info->getMoleculeStamp(stampId), stampId, i);
-+
-+
+            info->addMolecule(mol);
-+
+        }
-+
+    }
-+
-+
+}
-+
-+
+void SimSetup::compList(MakeStamps* stamps,SimInfo* info) {
-+
+    int i;
-+
+    char* id;
-+
+    MoleculeStamp* currentStamp;
-+
+    Component* the_components = info->getGlobals()->getComponents();
-+
+    int n_components = info->getGlobals()->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();
-+
+            }
-+
-+
+            info->addMoleculeStamp(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::setGlobalIndices(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;
-+
-+
+#ifdef 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);
-+
-+
+    tmpAtomInProc[myNode] = info->getNAtoms();
-+
+    tmpRigidBodiesInProc[myNode] = info->getNRigidBodiess();
-+
+    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroupss();
-+
-+
+    //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++);
-+
+        }
-+
-+
+    }
-+
+}
-+
+} //end namespace oopse

Diff Legend

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

Comparing branches/new_design/OOPSE-2.0/src/brains/SimCreator.cpp (file contents): Revision 1712 by tim, Thu Nov 4 20:55:01 2004 UTC vs. Revision 1727 by tim, Thu Nov 11 16:41:58 2004 UTC

Diff Legend

Comparing branches/new_design/OOPSE-2.0/src/brains/SimCreator.cpp (file contents):
Revision 1712 by tim, Thu Nov 4 20:55:01 2004 UTC vs.
Revision 1727 by tim, Thu Nov 11 16:41:58 2004 UTC