src/brains/SimCreator.cpp

 /*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
/**
 * @file SimCreator.cpp
 * @author tlin
 * @date 11/03/2004
 * @time 13:51am
 * @version 1.0
 */

#include <sprng.h>

#include "brains/MoleculeCreator.hpp"
#include "brains/SimCreator.hpp"
#include "brains/SimSnapshotManager.hpp"
#include "io/DumpReader.hpp"
#include "io/parse_me.h"
#include "UseTheForce/ForceFieldFactory.hpp"
#include "utils/simError.h"
#include "utils/StringUtils.hpp"
#ifdef IS_MPI
#include "io/mpiBASS.h"
#include "math/randomSPRNG.hpp"
#endif

namespace oopse {

void SimCreator::parseFile(const std::string mdFileName,  MakeStamps* stamps, Globals* simParams){

#ifdef IS_MPI

    if (worldRank == 0) {
#endif // is_mpi

        simParams->initalize();
        set_interface_stamps(stamps, simParams);

#ifdef IS_MPI

        mpiEventInit();

#endif

        yacc_BASS(mdFileName.c_str());

#ifdef IS_MPI

        throwMPIEvent(NULL);
    } else {
        set_interface_stamps(stamps, simParams);
        mpiEventInit();
        MPIcheckPoint();
        mpiEventLoop();
    }

#endif

}

SimInfo*  SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
    
    MakeStamps * stamps = new MakeStamps();

    Globals * simParams = new Globals();

    //parse meta-data file
    parseFile(mdFileName, stamps, simParams);

    //create the force field
    ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
                          simParams->getForceField());
    
    if (ff == NULL) {
        sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
                simParams->getForceField());
        painCave.isFatal = 1;
        simError();
    }

    if (simParams->haveForceFieldFileName()) {
        ff->setForceFieldFileName(simParams->getForceFieldFileName());
    }
    
    std::string forcefieldFileName;
    forcefieldFileName = ff->getForceFieldFileName();

    if (simParams->haveForceFieldVariant()) {
        //If the force field has variant, the variant force field name will be
        //Base.variant.frc. For exampel EAM.u6.frc
        
        std::string variant = simParams->getForceFieldVariant();

        std::string::size_type pos = forcefieldFileName.rfind(".frc");
        variant = "." + variant;
        if (pos != std::string::npos) {
            forcefieldFileName.insert(pos, variant);
        } else {
            //If the default force field file name does not containt .frc suffix, just append the .variant
            forcefieldFileName.append(variant);
        }
    } 
    
    ff->parse(forcefieldFileName);
    
    //extract the molecule stamps
    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
    compList(stamps, simParams, moleculeStampPairs);

    //create SimInfo
    SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);

    //gather parameters (SimCreator only retrieves part of the parameters)
    gatherParameters(info, mdFileName);

    //divide the molecules and determine the global index of molecules
#ifdef IS_MPI
    divideMolecules(info);
#endif 

    //create the molecules
    createMolecules(info);


    //allocate memory for DataStorage(circular reference, need to break it)
    info->setSnapshotManager(new SimSnapshotManager(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. 
    setGlobalIndex(info);

    //Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
    //atoms don't have the global index yet  (their global index are all initialized to -1).
    //Therefore we have to call addExcludePairs explicitly here. A way to work around is that
    //we can determine the beginning global indices of atoms before they get created.
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
        info->addExcludePairs(mol);
    }
    

    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
    //eta, chi for NPT integrator)
    if (loadInitCoords)
        loadCoordinates(info);    
    
    return info;
}

void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {

    //setup seed for random number generator
    int seedValue;
    Globals * simParams = info->getSimParams();

    if (simParams->haveSeed()) {
        seedValue = simParams->getSeed();

        if (seedValue < 100000000 ) {
            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 simParams->haveSeed()

    info->setSeed(seedValue);


    //figure out the ouput file names
    std::string prefix;

#ifdef IS_MPI

    if (worldRank == 0) {
#endif // is_mpi

        if (simParams->haveFinalConfig()) {
            prefix = getPrefix(simParams->getFinalConfig());
        } else {
            prefix = 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;
    std::vector<int> atomsPerProc;
    int nGlobalMols = info->getNGlobalMolecules();
    std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
    
    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);

    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();
    }

    MTRand myRandom(info->getSeed(), nProcessors, worldRank);


    a = 3.0 * nGlobalMols / info->getNGlobalAtoms();

    //initialize atomsPerProc
    atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);

    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.rand() * 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[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
    } else {

        // Listen to your marching orders from processor 0:

        MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
    }

    info->setMolToProcMap(molToProcMap);
    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->getLocalIndexManager());

            info->addMolecule(mol);

#ifdef IS_MPI

        }

#endif

    } //end for(int i=0)   
}

void SimCreator::compList(MakeStamps *stamps, Globals* simParams,
                        std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
    int i;
    char * id;
    MoleculeStamp * currentStamp;
    Component** the_components = simParams->getComponents();
    int n_components = simParams->getNComponents();

    if (!simParams->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,
                        "SimCreator 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))->getStamp();

            if (currentStamp == NULL) {
                sprintf(painCave.errMsg,
                        "SimCreator error: Component \"%s\" was not found in the "
                            "list of declared molecules\n", id);

                painCave.isFatal = 1;
                simError();
            }

            moleculeStampPairs.push_back(
                std::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) {
    SimInfo::MoleculeIterator mi;
    Molecule::AtomIterator ai;
    Molecule::RigidBodyIterator ri;
    Molecule::CutoffGroupIterator ci;
    Molecule * mol;
    Atom * atom;
    RigidBody * rb;
    CutoffGroup * cg;
    int beginAtomIndex;
    int beginRigidBodyIndex;
    int beginCutoffGroupIndex;
    int nGlobalAtoms = info->getNGlobalAtoms();
    
#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 < myNode; i++) {
        beginAtomIndex += NumAtomsInProc[i];
        beginRigidBodyIndex += NumRigidBodiesInProc[i];
        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
    }

#endif

    //rigidbody's index begins right after atom's
    beginRigidBodyIndex += info->getNGlobalAtoms();

    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++);
        }
    }

    //fill globalGroupMembership
    std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
    for(mol = info->beginMolecule(mi); mol != NULL; mol = info->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();
            }

        }       
    }

#ifdef IS_MPI    
    // 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.
    std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
    MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
     info->setGlobalGroupMembership(tmpGroupMembership);
#else
    info->setGlobalGroupMembership(globalGroupMembership);
#endif

    //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();
        }
    }

#ifdef IS_MPI
    std::vector<int> tmpMolMembership(nGlobalAtoms, 0);

    MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    
    info->setGlobalMolMembership(tmpMolMembership);
#else
    info->setGlobalMolMembership(globalMolMembership);
#endif

}

void SimCreator::loadCoordinates(SimInfo* info) {
    Globals* simParams;
    simParams = info->getSimParams();
    
    if (!simParams->haveInitialConfig()) {
        sprintf(painCave.errMsg,
                "Cannot intialize a simulation without an initial configuration file.\n");
        painCave.isFatal = 1;;
        simError();
    }
        
    DumpReader reader(info, simParams->getInitialConfig());
    int nframes = reader.getNFrames();

    if (nframes > 0) {
        reader.readFrame(nframes - 1);
    } else {
        //invalid initial coordinate file
        sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
                simParams->getInitialConfig());
        painCave.isFatal = 1;
        simError();
    }

    //copy the current snapshot to previous snapshot
    info->getSnapshotManager()->advance();
}

} //end namespace oopse


Revision:	381
Committed:	Tue Mar 1 14:45:45 2005 UTC (20 years, 8 months ago) by tim
File size:	21567 byte(s)
Log Message:	adding MersenneTwister random number generator
#	Content
1	/*
2	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	*
4	* The University of Notre Dame grants you ("Licensee") a
5	* non-exclusive, royalty free, license to use, modify and
6	* redistribute this software in source and binary code form, provided
7	* that the following conditions are met:
8	*
9	* 1. Acknowledgement of the program authors must be made in any
10	* publication of scientific results based in part on use of the
11	* program. An acceptable form of acknowledgement is citation of
12	* the article in which the program was described (Matthew
13	* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	* Parallel Simulation Engine for Molecular Dynamics,"
16	* J. Comput. Chem. 26, pp. 252-271 (2005))
17	*
18	* 2. Redistributions of source code must retain the above copyright
19	* notice, this list of conditions and the following disclaimer.
20	*
21	* 3. Redistributions in binary form must reproduce the above copyright
22	* notice, this list of conditions and the following disclaimer in the
23	* documentation and/or other materials provided with the
24	* distribution.
25	*
26	* This software is provided "AS IS," without a warranty of any
27	* kind. All express or implied conditions, representations and
28	* warranties, including any implied warranty of merchantability,
29	* fitness for a particular purpose or non-infringement, are hereby
30	* excluded. The University of Notre Dame and its licensors shall not
31	* be liable for any damages suffered by licensee as a result of
32	* using, modifying or distributing the software or its
33	* derivatives. In no event will the University of Notre Dame or its
34	* licensors be liable for any lost revenue, profit or data, or for
35	* direct, indirect, special, consequential, incidental or punitive
36	* damages, however caused and regardless of the theory of liability,
37	* arising out of the use of or inability to use software, even if the
38	* University of Notre Dame has been advised of the possibility of
39	* such damages.
40	*/
41
42	/**
43	* @file SimCreator.cpp
44	* @author tlin
45	* @date 11/03/2004
46	* @time 13:51am
47	* @version 1.0
48	*/
49
50	#include <sprng.h>
51
52	#include "brains/MoleculeCreator.hpp"
53	#include "brains/SimCreator.hpp"
54	#include "brains/SimSnapshotManager.hpp"
55	#include "io/DumpReader.hpp"
56	#include "io/parse_me.h"
57	#include "UseTheForce/ForceFieldFactory.hpp"
58	#include "utils/simError.h"
59	#include "utils/StringUtils.hpp"
60	#ifdef IS_MPI
61	#include "io/mpiBASS.h"
62	#include "math/randomSPRNG.hpp"
63	#endif
64
65	namespace oopse {
66
67	void SimCreator::parseFile(const std::string mdFileName, MakeStamps* stamps, Globals* simParams){
68
69	#ifdef IS_MPI
70
71	if (worldRank == 0) {
72	#endif // is_mpi
73
74	simParams->initalize();
75	set_interface_stamps(stamps, simParams);
76
77	#ifdef IS_MPI
78
79	mpiEventInit();
80
81	#endif
82
83	yacc_BASS(mdFileName.c_str());
84
85	#ifdef IS_MPI
86
87	throwMPIEvent(NULL);
88	} else {
89	set_interface_stamps(stamps, simParams);
90	mpiEventInit();
91	MPIcheckPoint();
92	mpiEventLoop();
93	}
94
95	#endif
96
97	}
98
99	SimInfo* SimCreator::createSim(const std::string & mdFileName, bool loadInitCoords) {
100
101	MakeStamps * stamps = new MakeStamps();
102
103	Globals * simParams = new Globals();
104
105	//parse meta-data file
106	parseFile(mdFileName, stamps, simParams);
107
108	//create the force field
109	ForceField * ff = ForceFieldFactory::getInstance()->createForceField(
110	simParams->getForceField());
111
112	if (ff == NULL) {
113	sprintf(painCave.errMsg, "ForceField Factory can not create %s force field\n",
114	simParams->getForceField());
115	painCave.isFatal = 1;
116	simError();
117	}
118
119	if (simParams->haveForceFieldFileName()) {
120	ff->setForceFieldFileName(simParams->getForceFieldFileName());
121	}
122
123	std::string forcefieldFileName;
124	forcefieldFileName = ff->getForceFieldFileName();
125
126	if (simParams->haveForceFieldVariant()) {
127	//If the force field has variant, the variant force field name will be
128	//Base.variant.frc. For exampel EAM.u6.frc
129
130	std::string variant = simParams->getForceFieldVariant();
131
132	std::string::size_type pos = forcefieldFileName.rfind(".frc");
133	variant = "." + variant;
134	if (pos != std::string::npos) {
135	forcefieldFileName.insert(pos, variant);
136	} else {
137	//If the default force field file name does not containt .frc suffix, just append the .variant
138	forcefieldFileName.append(variant);
139	}
140	}
141
142	ff->parse(forcefieldFileName);
143
144	//extract the molecule stamps
145	std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
146	compList(stamps, simParams, moleculeStampPairs);
147
148	//create SimInfo
149	SimInfo * info = new SimInfo(moleculeStampPairs, ff, simParams);
150
151	//gather parameters (SimCreator only retrieves part of the parameters)
152	gatherParameters(info, mdFileName);
153
154	//divide the molecules and determine the global index of molecules
155	#ifdef IS_MPI
156	divideMolecules(info);
157	#endif
158
159	//create the molecules
160	createMolecules(info);
161
162
163	//allocate memory for DataStorage(circular reference, need to break it)
164	info->setSnapshotManager(new SimSnapshotManager(info));
165
166	//set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
167	//global index will never change again). Local indices of atoms and rigidbodies are already set by
168	//MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
169	setGlobalIndex(info);
170
171	//Alought addExculdePairs is called inside SimInfo's addMolecule method, at that point
172	//atoms don't have the global index yet (their global index are all initialized to -1).
173	//Therefore we have to call addExcludePairs explicitly here. A way to work around is that
174	//we can determine the beginning global indices of atoms before they get created.
175	SimInfo::MoleculeIterator mi;
176	Molecule* mol;
177	for (mol= info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
178	info->addExcludePairs(mol);
179	}
180
181
182	//load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
183	//eta, chi for NPT integrator)
184	if (loadInitCoords)
185	loadCoordinates(info);
186
187	return info;
188	}
189
190	void SimCreator::gatherParameters(SimInfo *info, const std::string& mdfile) {
191
192	//setup seed for random number generator
193	int seedValue;
194	Globals * simParams = info->getSimParams();
195
196	if (simParams->haveSeed()) {
197	seedValue = simParams->getSeed();
198
199	if (seedValue < 100000000 ) {
200	sprintf(painCave.errMsg,
201	"Seed for sprng library should contain at least 9 digits\n"
202	"OOPSE will generate a seed for user\n");
203
204	painCave.isFatal = 0;
205	simError();
206
207	//using seed generated by system instead of invalid seed set by user
208
209	#ifndef IS_MPI
210
211	seedValue = make_sprng_seed();
212
213	#else
214
215	if (worldRank == 0) {
216	seedValue = make_sprng_seed();
217	}
218
219	MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
220
221	#endif
222
223	} //end if (seedValue /1000000000 == 0)
224	} else {
225
226	#ifndef IS_MPI
227
228	seedValue = make_sprng_seed();
229
230	#else
231
232	if (worldRank == 0) {
233	seedValue = make_sprng_seed();
234	}
235
236	MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
237
238	#endif
239
240	} //end of simParams->haveSeed()
241
242	info->setSeed(seedValue);
243
244
245	//figure out the ouput file names
246	std::string prefix;
247
248	#ifdef IS_MPI
249
250	if (worldRank == 0) {
251	#endif // is_mpi
252
253	if (simParams->haveFinalConfig()) {
254	prefix = getPrefix(simParams->getFinalConfig());
255	} else {
256	prefix = getPrefix(mdfile);
257	}
258
259	info->setFinalConfigFileName(prefix + ".eor");
260	info->setDumpFileName(prefix + ".dump");
261	info->setStatFileName(prefix + ".stat");
262
263	#ifdef IS_MPI
264
265	}
266
267	#endif
268
269	}
270
271	#ifdef IS_MPI
272	void SimCreator::divideMolecules(SimInfo *info) {
273	double numerator;
274	double denominator;
275	double precast;
276	double x;
277	double y;
278	double a;
279	int old_atoms;
280	int add_atoms;
281	int new_atoms;
282	int nTarget;
283	int done;
284	int i;
285	int j;
286	int loops;
287	int which_proc;
288	int nProcessors;
289	std::vector<int> atomsPerProc;
290	int nGlobalMols = info->getNGlobalMolecules();
291	std::vector<int> molToProcMap(nGlobalMols, -1); // default to an error condition:
292
293	MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
294
295	if (nProcessors > nGlobalMols) {
296	sprintf(painCave.errMsg,
297	"nProcessors (%d) > nMol (%d)\n"
298	"\tThe number of processors is larger than\n"
299	"\tthe number of molecules. This will not result in a \n"
300	"\tusable division of atoms for force decomposition.\n"
301	"\tEither try a smaller number of processors, or run the\n"
302	"\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
303
304	painCave.isFatal = 1;
305	simError();
306	}
307
308	MTRand myRandom(info->getSeed(), nProcessors, worldRank);
309
310
311	a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
312
313	//initialize atomsPerProc
314	atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
315
316	if (worldRank == 0) {
317	numerator = info->getNGlobalAtoms();
318	denominator = nProcessors;
319	precast = numerator / denominator;
320	nTarget = (int)(precast + 0.5);
321
322	for(i = 0; i < nGlobalMols; i++) {
323	done = 0;
324	loops = 0;
325
326	while (!done) {
327	loops++;
328
329	// Pick a processor at random
330
331	which_proc = (int) (myRandom.rand() * nProcessors);
332
333	//get the molecule stamp first
334	int stampId = info->getMoleculeStampId(i);
335	MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
336
337	// How many atoms does this processor have so far?
338	old_atoms = atomsPerProc[which_proc];
339	add_atoms = moleculeStamp->getNAtoms();
340	new_atoms = old_atoms + add_atoms;
341
342	// If we've been through this loop too many times, we need
343	// to just give up and assign the molecule to this processor
344	// and be done with it.
345
346	if (loops > 100) {
347	sprintf(painCave.errMsg,
348	"I've tried 100 times to assign molecule %d to a "
349	" processor, but can't find a good spot.\n"
350	"I'm assigning it at random to processor %d.\n",
351	i, which_proc);
352
353	painCave.isFatal = 0;
354	simError();
355
356	molToProcMap[i] = which_proc;
357	atomsPerProc[which_proc] += add_atoms;
358
359	done = 1;
360	continue;
361	}
362
363	// If we can add this molecule to this processor without sending
364	// it above nTarget, then go ahead and do it:
365
366	if (new_atoms <= nTarget) {
367	molToProcMap[i] = which_proc;
368	atomsPerProc[which_proc] += add_atoms;
369
370	done = 1;
371	continue;
372	}
373
374	// The only situation left is when new_atoms > nTarget. We
375	// want to accept this with some probability that dies off the
376	// farther we are from nTarget
377
378	// roughly: x = new_atoms - nTarget
379	// Pacc(x) = exp(- a * x)
380	// where a = penalty / (average atoms per molecule)
381
382	x = (double)(new_atoms - nTarget);
383	y = myRandom.getRandom();
384
385	if (y < exp(- a * x)) {
386	molToProcMap[i] = which_proc;
387	atomsPerProc[which_proc] += add_atoms;
388
389	done = 1;
390	continue;
391	} else {
392	continue;
393	}
394	}
395	}
396
397	// Spray out this nonsense to all other processors:
398
399	MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
400	} else {
401
402	// Listen to your marching orders from processor 0:
403
404	MPI_Bcast(&molToProcMap[0], nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
405	}
406
407	info->setMolToProcMap(molToProcMap);
408	sprintf(checkPointMsg,
409	"Successfully divided the molecules among the processors.\n");
410	MPIcheckPoint();
411	}
412
413	#endif
414
415	void SimCreator::createMolecules(SimInfo *info) {
416	MoleculeCreator molCreator;
417	int stampId;
418
419	for(int i = 0; i < info->getNGlobalMolecules(); i++) {
420
421	#ifdef IS_MPI
422
423	if (info->getMolToProc(i) == worldRank) {
424	#endif
425
426	stampId = info->getMoleculeStampId(i);
427	Molecule * mol = molCreator.createMolecule(info->getForceField(), info->getMoleculeStamp(stampId),
428	stampId, i, info->getLocalIndexManager());
429
430	info->addMolecule(mol);
431
432	#ifdef IS_MPI
433
434	}
435
436	#endif
437
438	} //end for(int i=0)
439	}
440
441	void SimCreator::compList(MakeStamps stamps, Globals simParams,
442	std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
443	int i;
444	char * id;
445	MoleculeStamp * currentStamp;
446	Component** the_components = simParams->getComponents();
447	int n_components = simParams->getNComponents();
448
449	if (!simParams->haveNMol()) {
450	// we don't have the total number of molecules, so we assume it is
451	// given in each component
452
453	for(i = 0; i < n_components; i++) {
454	if (!the_components[i]->haveNMol()) {
455	// we have a problem
456	sprintf(painCave.errMsg,
457	"SimCreator Error. No global NMol or component NMol given.\n"
458	"\tCannot calculate the number of atoms.\n");
459
460	painCave.isFatal = 1;
461	simError();
462	}
463
464	id = the_components[i]->getType();
465	currentStamp = (stamps->extractMolStamp(id))->getStamp();
466
467	if (currentStamp == NULL) {
468	sprintf(painCave.errMsg,
469	"SimCreator error: Component \"%s\" was not found in the "
470	"list of declared molecules\n", id);
471
472	painCave.isFatal = 1;
473	simError();
474	}
475
476	moleculeStampPairs.push_back(
477	std::make_pair(currentStamp, the_components[i]->getNMol()));
478	} //end for (i = 0; i < n_components; i++)
479	} else {
480	sprintf(painCave.errMsg, "SimSetup error.\n"
481	"\tSorry, the ability to specify total"
482	" nMols and then give molfractions in the components\n"
483	"\tis not currently supported."
484	" Please give nMol in the components.\n");
485
486	painCave.isFatal = 1;
487	simError();
488	}
489
490	#ifdef IS_MPI
491
492	strcpy(checkPointMsg, "Component stamps successfully extracted\n");
493	MPIcheckPoint();
494
495	#endif // is_mpi
496
497	}
498
499	void SimCreator::setGlobalIndex(SimInfo *info) {
500	SimInfo::MoleculeIterator mi;
501	Molecule::AtomIterator ai;
502	Molecule::RigidBodyIterator ri;
503	Molecule::CutoffGroupIterator ci;
504	Molecule * mol;
505	Atom * atom;
506	RigidBody * rb;
507	CutoffGroup * cg;
508	int beginAtomIndex;
509	int beginRigidBodyIndex;
510	int beginCutoffGroupIndex;
511	int nGlobalAtoms = info->getNGlobalAtoms();
512
513	#ifndef IS_MPI
514
515	beginAtomIndex = 0;
516	beginRigidBodyIndex = 0;
517	beginCutoffGroupIndex = 0;
518
519	#else
520
521	int nproc;
522	int myNode;
523
524	myNode = worldRank;
525	MPI_Comm_size(MPI_COMM_WORLD, &nproc);
526
527	std::vector < int > tmpAtomsInProc(nproc, 0);
528	std::vector < int > tmpRigidBodiesInProc(nproc, 0);
529	std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
530	std::vector < int > NumAtomsInProc(nproc, 0);
531	std::vector < int > NumRigidBodiesInProc(nproc, 0);
532	std::vector < int > NumCutoffGroupsInProc(nproc, 0);
533
534	tmpAtomsInProc[myNode] = info->getNAtoms();
535	tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
536	tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
537
538	//do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
539	MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
540	MPI_SUM, MPI_COMM_WORLD);
541	MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
542	MPI_INT, MPI_SUM, MPI_COMM_WORLD);
543	MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
544	MPI_INT, MPI_SUM, MPI_COMM_WORLD);
545
546	beginAtomIndex = 0;
547	beginRigidBodyIndex = 0;
548	beginCutoffGroupIndex = 0;
549
550	for(int i = 0; i < myNode; i++) {
551	beginAtomIndex += NumAtomsInProc[i];
552	beginRigidBodyIndex += NumRigidBodiesInProc[i];
553	beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
554	}
555
556	#endif
557
558	//rigidbody's index begins right after atom's
559	beginRigidBodyIndex += info->getNGlobalAtoms();
560
561	for(mol = info->beginMolecule(mi); mol != NULL;
562	mol = info->nextMolecule(mi)) {
563
564	//local index(index in DataStorge) of atom is important
565	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
566	atom->setGlobalIndex(beginAtomIndex++);
567	}
568
569	for(rb = mol->beginRigidBody(ri); rb != NULL;
570	rb = mol->nextRigidBody(ri)) {
571	rb->setGlobalIndex(beginRigidBodyIndex++);
572	}
573
574	//local index of cutoff group is trivial, it only depends on the order of travesing
575	for(cg = mol->beginCutoffGroup(ci); cg != NULL;
576	cg = mol->nextCutoffGroup(ci)) {
577	cg->setGlobalIndex(beginCutoffGroupIndex++);
578	}
579	}
580
581	//fill globalGroupMembership
582	std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
583	for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
584	for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
585
586	for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
587	globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
588	}
589
590	}
591	}
592
593	#ifdef IS_MPI
594	// Since the globalGroupMembership has been zero filled and we've only
595	// poked values into the atoms we know, we can do an Allreduce
596	// to get the full globalGroupMembership array (We think).
597	// This would be prettier if we could use MPI_IN_PLACE like the MPI-2
598	// docs said we could.
599	std::vector<int> tmpGroupMembership(nGlobalAtoms, 0);
600	MPI_Allreduce(&globalGroupMembership[0], &tmpGroupMembership[0], nGlobalAtoms,
601	MPI_INT, MPI_SUM, MPI_COMM_WORLD);
602	info->setGlobalGroupMembership(tmpGroupMembership);
603	#else
604	info->setGlobalGroupMembership(globalGroupMembership);
605	#endif
606
607	//fill molMembership
608	std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
609
610	for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
611
612	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
613	globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
614	}
615	}
616
617	#ifdef IS_MPI
618	std::vector<int> tmpMolMembership(nGlobalAtoms, 0);
619
620	MPI_Allreduce(&globalMolMembership[0], &tmpMolMembership[0], nGlobalAtoms,
621	MPI_INT, MPI_SUM, MPI_COMM_WORLD);
622
623	info->setGlobalMolMembership(tmpMolMembership);
624	#else
625	info->setGlobalMolMembership(globalMolMembership);
626	#endif
627
628	}
629
630	void SimCreator::loadCoordinates(SimInfo* info) {
631	Globals* simParams;
632	simParams = info->getSimParams();
633
634	if (!simParams->haveInitialConfig()) {
635	sprintf(painCave.errMsg,
636	"Cannot intialize a simulation without an initial configuration file.\n");
637	painCave.isFatal = 1;;
638	simError();
639	}
640
641	DumpReader reader(info, simParams->getInitialConfig());
642	int nframes = reader.getNFrames();
643
644	if (nframes > 0) {
645	reader.readFrame(nframes - 1);
646	} else {
647	//invalid initial coordinate file
648	sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
649	simParams->getInitialConfig());
650	painCave.isFatal = 1;
651	simError();
652	}
653
654	//copy the current snapshot to previous snapshot
655	info->getSnapshotManager()->advance();
656	}
657
658	} //end namespace oopse
659
660