src/brains/SimCreator.cpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * 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 <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();
    }

    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;
    randomSPRNG myRandom(info->getSeed());
    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();
    }

    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.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[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

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