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 "brains/SimCreator.hpp"

namespace oopse {

void SimSetup::parseFile(const std::string& mdFileName, MakeStamps *stamps, Globals *globals) {

#ifdef IS_MPI

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

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

#ifdef IS_MPI

        mpiEventInit();

#endif

        yacc_BASS(mdFileName.c_str());

#ifdef IS_MPI

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

#endif

}

SimInfo*  SimCreator::createSim(const std::string & mdFileName) {

    mdFileName_ = mdFileName;
    
    MakeStamps * stamps = new MakeStamps();

    Globals * globals = new Globals();

    //parse meta-data file
    parseFile(mdFileName_, stamps, globals);

    //create the force field
    ForceFiled * ff = ForceFieldFactory::getInstance()->createObject(
                          globals->getForceField());

    if (ff == NULL) {
        sprintf(painCave.errMsg, "ForceFiled Factory can not create %s force field\n",
                globals->getForceField());
        painCave.isFatal = 1;
        simError();
    }

    //extract the molecule stamps
    std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
    compList(stamps, globals, moleculeStampPairs);

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

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

    //divide the molecules and determine the global index of molecules
    divideMolecules(info);

    //create the molecules
    createMolecules(info);

    //set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
    //global index will never change again). Local indices of atoms and rigidbodies are already set by 
    //MoleculeCreator class which actually delegates the responsibility to LocalIndexManager. 
    setGlobalIndices(info);

    //allocate memory for DataStorage(circular reference, need to break it)
    info->setSnapshotManager(new SimSnapshotManager(info));

    //initialize fortran -- setup the cutoff 
    initFortran(info);

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

void SimCreator::gatherParameters(SimInfo *info) {

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

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

        if (seedValue / 1000000000 == 0) {
            sprintf(painCave.errMsg,
                    "Seed for sprng library should contain at least 9 digits\n"
                        "OOPSE will generate a seed for user\n");

            painCave.isFatal = 0;
            simError();

            //using seed generated by system instead of invalid seed set by user 

#ifndef IS_MPI

            seedValue = make_sprng_seed();

#else

            if (worldRank == 0) {
                seedValue = make_sprng_seed();
            }

            MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);

#endif

        } //end if (seedValue /1000000000 == 0)
    } else {

#ifndef IS_MPI

        seedValue = make_sprng_seed();

#else

        if (worldRank == 0) {
            seedValue = make_sprng_seed();
        }

        MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);

#endif

    } //end of globals->haveSeed()

    info->setSeed(seedValue);


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

#ifdef IS_MPI

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

        if (globals->haveFinalConfig()) {
            prefix = StringUtils::getPrefix(globals->getFinalConfig());
        } else {
            prefix = StringUtils::getPrefix(mdfile_);
        }

        info->setFinalConfigFileName(prefix + ".eor");
        info->setDumpFileName(prefix + ".dump");
        info->setStatFileName(prefix + ".stat");

#ifdef IS_MPI

    }

#endif

}

#ifdef IS_MPI

void SimCreator::divideMolecules(SimInfo *info) {
    double numerator;
    double denominator;
    double precast;
    double x;
    double y;
    double a;
    int old_atoms;
    int add_atoms;
    int new_atoms;
    int nTarget;
    int done;
    int i;
    int j;
    int loops;
    int which_proc;
    int nProcessors;
    int nGlobalMols;
    std::vector < int > atomsPerProc;
    randomSPRNG myRandom(info->getSeed());
    int * molToProcMap = info->getMolToProcMapPointer();

    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
    nGlobalMols = info->getNGlobalMolecules();

    if (nProcessors > nGlobalMols) {
        sprintf(painCave.errMsg,
                "nProcessors (%d) > nMol (%d)\n"
                    "\tThe number of processors is larger than\n"
                    "\tthe number of molecules.  This will not result in a \n"
                    "\tusable division of atoms for force decomposition.\n"
                    "\tEither try a smaller number of processors, or run the\n"
                    "\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);

        painCave.isFatal = 1;
        simError();
    }

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

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

    for(i = 0; i < nGlobalMols; i++) {
        // default to an error condition:
        molToProcMap[i] = -1;
    }

    if (worldRank == 0) {
        numerator = info->getNGlobalAtoms();
        denominator = nProcessors;
        precast = numerator / denominator;
        nTarget = (int)(precast + 0.5);

        for(i = 0; i < nGlobalMols; i++) {
            done = 0;
            loops = 0;

            while (!done) {
                loops++;

                // Pick a processor at random

                which_proc = (int) (myRandom.getRandom() * nProcessors);

                //get the molecule stamp first
                int stampId = info->getMoleculeStampId(i);
                MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);

                // How many atoms does this processor have so far?
                old_atoms = atomsPerProc[which_proc];
                add_atoms = moleculeStamp->getNAtoms();
                new_atoms = old_atoms + add_atoms;

                // If we've been through this loop too many times, we need
                // to just give up and assign the molecule to this processor
                // and be done with it. 

                if (loops > 100) {
                    sprintf(painCave.errMsg,
                            "I've tried 100 times to assign molecule %d to a "
                                " processor, but can't find a good spot.\n"
                                "I'm assigning it at random to processor %d.\n",
                            i, which_proc);

                    painCave.isFatal = 0;
                    simError();

                    molToProcMap[i] = which_proc;
                    atomsPerProc[which_proc] += add_atoms;

                    done = 1;
                    continue;
                }

                // If we can add this molecule to this processor without sending
                // it above nTarget, then go ahead and do it:

                if (new_atoms <= nTarget) {
                    molToProcMap[i] = which_proc;
                    atomsPerProc[which_proc] += add_atoms;

                    done = 1;
                    continue;
                }

                // The only situation left is when new_atoms > nTarget.  We
                // want to accept this with some probability that dies off the
                // farther we are from nTarget

                // roughly:  x = new_atoms - nTarget
                //           Pacc(x) = exp(- a * x)
                // where a = penalty / (average atoms per molecule)

                x = (double)(new_atoms - nTarget);
                y = myRandom.getRandom();

                if (y < exp(- a * x)) {
                    molToProcMap[i] = which_proc;
                    atomsPerProc[which_proc] += add_atoms;

                    done = 1;
                    continue;
                } else {
                    continue;
                }
            }
        }

        // Spray out this nonsense to all other processors:

        MPI_Bcast(molToProcMap, nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
    } else {

        // Listen to your marching orders from processor 0:

        MPI_Bcast(molToProcMap, nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
    }

    sprintf(checkPointMsg,
            "Successfully divided the molecules among the processors.\n");
    MPIcheckPoint();
}

#endif

void SimCreator::createMolecules(SimInfo *info) {
    MoleculeCreator molCreator;
    int stampId;

    for(int i = 0; i < info->getNGlobalMolecules(); i++) {

#ifdef IS_MPI

        if (info->getMolToProc(i) == worldRank) {
#endif

            stampId = info->getMoleculeStampId(i);
            Molecule * mol = molCreator.createMolecule(info->getForceField(),
                                                info->getMoleculeStamp(stampId), stampId, i);

            info->addMolecule(mol);

#ifdef IS_MPI

        }

#endif

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

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

    if (!globals->haveNMol()) {
        // we don't have the total number of molecules, so we assume it is
        // given in each component

        for(i = 0; i < n_components; i++) {
            if (!the_components[i]->haveNMol()) {
                // we have a problem
                sprintf(painCave.errMsg,
                        "SimSetup Error. No global NMol or component NMol given.\n"
                            "\tCannot calculate the number of atoms.\n");

                painCave.isFatal = 1;
                simError();
            }

            id = the_components[i]->getType();
            currentStamp = stamps->extractMolStamp(id);

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

                painCave.isFatal = 1;
                simError();
            }

            moleculeStampPairs.push_back(
                make_pair(currentStamp, the_components[i]->getNMol));
        } //end for (i = 0; i < n_components; i++)
    } else {
        sprintf(painCave.errMsg, "SimSetup error.\n"
                                     "\tSorry, the ability to specify total"
                                     " nMols and then give molfractions in the components\n"
                                     "\tis not currently supported."
                                     " Please give nMol in the components.\n");

        painCave.isFatal = 1;
        simError();
    }

#ifdef IS_MPI

    strcpy(checkPointMsg, "Component stamps successfully extracted\n");
    MPIcheckPoint();

#endif // is_mpi

}

void SimCreator::setGlobalIndex(SimInfo *info) {
    typename SimInfo::MoleculeIterator mi;
    typename Molecule::AtomIterator ai;
    typename Molecule::RigidBodyIterator ri;
    typename Molecule::CutoffGroupIterator ci;
    Molecule * mol;
    Atom * atom;
    RigidBody * rb;
    CutoffGroup * cg;
    int beginAtomIndex;
    int beginRigidBodyIndex;
    int beginCutoffGroupIndex;

#ifndef IS_MPI

    beginAtomIndex = 0;
    beginRigidBodyIndex = 0;
    beginCutoffGroupIndex = 0;

#else

    int nproc;
    int myNode;

    myNode = worldRank;
    MPI_Comm_size(MPI_COMM_WORLD, &nproc);

    std::vector < int > tmpAtomsInProc(nproc, 0);
    std::vector < int > tmpRigidBodiesInProc(nproc, 0);
    std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
    std::vector < int > NumAtomsInProc(nproc, 0);
    std::vector < int > NumRigidBodiesInProc(nproc, 0);
    std::vector < int > NumCutoffGroupsInProc(nproc, 0);

    tmpAtomsInProc[myNode] = info->getNAtoms();
    tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();

    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
                  MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);

    beginAtomIndex = 0;
    beginRigidBodyIndex = 0;
    beginCutoffGroupIndex = 0;

    for(int i = 0; i < nproc; i++) {
        beginAtomIndex += NumAtomsInProc[i];
        beginRigidBodyIndex += NumRigidBodiesInProc[i];
        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
    }

#endif

    for(mol = info->beginMolecule(mi); mol != NULL;
        mol = info->nextMolecule(mi)) {

        //local index(index in DataStorge) of atom is important
        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
            atom->setGlobalIndex(beginAtomIndex++);
        }

        for(rb = mol->beginRigidBody(ri); rb != NULL;
            rb = mol->nextRigidBody(ri)) {
            rb->setGlobalIndex(beginRigidBodyIndex++);
        }

        //local index of cutoff group is trivial, it only depends on the order of travesing
        for(cg = mol->beginCutoffGroup(ci); cg != NULL;
            cg = mol->nextCutoffGroup(ci)) {
            cg->setGlobalIndex(beginCutoffGroupIndex++);
        }
    }

    //fill globalGroupMembership
    std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {

            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
                globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
            }

        }       
    }

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

    MPI_Allreduce(&globalGroupMembership[0],
                  info->getGlobalGroupMembershipPointer(),
                  info->getNGlobalAtoms(),
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
#else
    std::copy(globalGroupMembership.begin(), globalGroupMembership.end(), 
        info->getGlobalGroupMembershipPointer());
#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
    MPI_Allreduce(&globalMolMembership[0],
                  info->getGlobalMolMembershipPointer(),
                  info->getNGlobalAtoms(),
                  MPI_INT, MPI_SUM, MPI_COMM_WORLD);
#else
    std::copy(globalMolMembership.begin(), globalMolMembership.end(), 
        info->getGlobalMolMembershipPointer());
#endif

}


void SimCreator::initFortran(SimInfo* info) {

    //setup fortran simulation (Actually, in parallel version it will initialze the parallel part of fortran first)
    info->update();

    //notify fortran whether reaction field is used or not
    //deprecated
    int isError;
    isError = 0;
    initFortranFF( &useReactionField, &isError );

    if(isError){
        sprintf( painCave.errMsg,
        "SimCreator::initFortran() error: There was an error initializing the forceField in fortran.\n" );
        painCave.isFatal = 1;
        simError();
    }
    
    //figure out the cutoff radius and pass it to fortran
    setCutoff(info);
}

void SimCreator::setCutoff(SimInfo* info) {
    Globals* globals = info->getGlobals();
    double rcut;  //cutoff radius
    double rsw; //switching radius

    simtype* st = info->getSimType();
    
    if (st->SIM_uses_Charges | st->SIM_uses_Dipoles | st->SIM_uses_RF) {
        
        if (!globals->haveRcut()){
            sprintf(painCave.errMsg,
                "SimCreator Warning: No value was set for the cutoffRadius.\n"
                "\tOOPSE will use a default value of 15.0 angstroms"
                "\tfor the cutoffRadius.\n");
            painCave.isFatal = 0;
            simError();
            rcut = 15.0;
        } else{
            rcut = globals->getRcut();
        }

        if (!globals->haveRsw()){
            sprintf(painCave.errMsg,
                "SimCreator Warning: No value was set for switchingRadius.\n"
                "\tOOPSE will use a default value of\n"
                "\t0.95 * cutoffRadius for the switchingRadius\n");
            painCave.isFatal = 0;
            simError();
            rsw = 0.95 * rcut;
        } else{
            rsw = globals->getRsw();
        }

    } else {
        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
        
        if (globals->haveRcut()) {
            rcut = globals->getRcut();
        } else {
            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
            rcut = info->calcMaxCutoffRadius();
        }

        if (globals->haveRsw()) {
            rsw  = globals->getRsw()
        } else {
            rsw = rcut;
        }
    
    }
    
    //store them into Siminfo
    info->setRcut(rcut);
    info->setRsw(rsw);
    
    double rnblist = rcut + 1; // skin of neighbor list

    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
    notifyFortranCutoffs(&rcut, &rsw, &rnblist);
}

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

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

}

} //end namespace oopse
Revision:	1735
Committed:	Fri Nov 12 17:40:03 2004 UTC (19 years, 9 months ago) by tim
File size:	20551 byte(s)
Log Message:	SimCreator and SimInfo are ready for unit test
#	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			#include "brains/SimCreator.hpp"
35	tim	1733
36	tim	1703	namespace oopse {
37
38	tim	1733	void SimSetup::parseFile(const std::string& mdFileName, MakeStamps stamps, Globals globals) {
39
40	tim	1712	#ifdef IS_MPI
41	tim	1733
42			if (worldRank == 0) {
43	tim	1712	#endif // is_mpi
44
45	tim	1733	globals->initalize();
46			set_interface_stamps(stamps, globals);
47	tim	1712
48	tim	1733	#ifdef IS_MPI
49	tim	1712
50	tim	1733	mpiEventInit();
51
52	tim	1712	#endif
53
54	tim	1733	yacc_BASS(mdFileName.c_str());
55	tim	1712
56			#ifdef IS_MPI
57	tim	1733
58			throwMPIEvent(NULL);
59			} else {
60			set_interface_stamps(stamps, globals);
61			mpiEventInit();
62			MPIcheckPoint();
63			mpiEventLoop();
64			}
65
66	tim	1712	#endif
67
68	tim	1703	}
69
70	tim	1733	SimInfo* SimCreator::createSim(const std::string & mdFileName) {
71	tim	1712
72	tim	1733	mdFileName_ = mdFileName;
73	tim	1712
74	tim	1733	MakeStamps * stamps = new MakeStamps();
75	tim	1712
76	tim	1733	Globals * globals = new Globals();
77
78	tim	1712	//parse meta-data file
79	tim	1733	parseFile(mdFileName_, stamps, globals);
80	tim	1712
81			//create the force field
82	tim	1733	ForceFiled * ff = ForceFieldFactory::getInstance()->createObject(
83			globals->getForceField());
84
85			if (ff == NULL) {
86			sprintf(painCave.errMsg, "ForceFiled Factory can not create %s force field\n",
87			globals->getForceField());
88			painCave.isFatal = 1;
89			simError();
90			}
91
92			//extract the molecule stamps
93			std::vector < std::pair<MoleculeStamp *, int> > moleculeStampPairs;
94			compList(stamps, globals, moleculeStampPairs);
95
96	tim	1719	//create SimInfo
97	tim	1733	SimInfo * info = new SimInfo(moleculeStampPairs, ff, globals);
98	tim	1719
99	tim	1733	//gather parameters (SimCreator only retrieves part of the parameters)
100	tim	1719	gatherParameters(info);
101	tim	1733
102	tim	1712	//divide the molecules and determine the global index of molecules
103	tim	1733	divideMolecules(info);
104
105	tim	1712	//create the molecules
106	tim	1719	createMolecules(info);
107	tim	1712
108	tim	1725	//set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
109			//global index will never change again). Local indices of atoms and rigidbodies are already set by
110	tim	1733	//MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
111	tim	1725	setGlobalIndices(info);
112	tim	1733
113	tim	1712	//allocate memory for DataStorage(circular reference, need to break it)
114	tim	1733	info->setSnapshotManager(new SimSnapshotManager(info));
115
116			//initialize fortran -- setup the cutoff
117			initFortran(info);
118
119	tim	1725	//load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
120			//eta, chi for NPT integrator)
121	tim	1733	loadCoordinates(info);
122	tim	1719	return info;
123	tim	1733	}
124	tim	1712
125	tim	1733	void SimCreator::gatherParameters(SimInfo *info) {
126	tim	1725
127			//setup seed for random number generator
128			int seedValue;
129	tim	1733	Globals * globals = info->getGlobals();
130
131			if (globals->haveSeed()) {
132	tim	1725	seedValue = globals->getSeed();
133
134	tim	1733	if (seedValue / 1000000000 == 0) {
135	tim	1725	sprintf(painCave.errMsg,
136	tim	1733	"Seed for sprng library should contain at least 9 digits\n"
137			"OOPSE will generate a seed for user\n");
138
139	tim	1725	painCave.isFatal = 0;
140			simError();
141
142	tim	1733	//using seed generated by system instead of invalid seed set by user
143
144	tim	1725	#ifndef IS_MPI
145	tim	1733
146	tim	1725	seedValue = make_sprng_seed();
147	tim	1733
148			#else
149
150			if (worldRank == 0) {
151	tim	1725	seedValue = make_sprng_seed();
152			}
153	tim	1733
154			MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
155
156			#endif
157
158	tim	1725	} //end if (seedValue /1000000000 == 0)
159	tim	1733	} else {
160
161	tim	1725	#ifndef IS_MPI
162	tim	1733
163	tim	1725	seedValue = make_sprng_seed();
164	tim	1733
165			#else
166
167			if (worldRank == 0) {
168	tim	1725	seedValue = make_sprng_seed();
169			}
170	tim	1733
171			MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
172
173	tim	1725	#endif
174	tim	1733
175			} //end of globals->haveSeed()
176
177	tim	1725	info->setSeed(seedValue);
178
179	tim	1733
180			//figure out the ouput file names
181	tim	1725	std::string prefix;
182	tim	1733
183	tim	1725	#ifdef IS_MPI
184	tim	1733
185			if (worldRank == 0) {
186	tim	1725	#endif // is_mpi
187	tim	1733
188			if (globals->haveFinalConfig()) {
189			prefix = StringUtils::getPrefix(globals->getFinalConfig());
190	tim	1725	} else {
191	tim	1733	prefix = StringUtils::getPrefix(mdfile_);
192	tim	1725	}
193	tim	1733
194			info->setFinalConfigFileName(prefix + ".eor");
195	tim	1725	info->setDumpFileName(prefix + ".dump");
196			info->setStatFileName(prefix + ".stat");
197
198			#ifdef IS_MPI
199	tim	1733
200	tim	1725	}
201
202	tim	1733	#endif
203	tim	1725
204	tim	1712	}
205
206			#ifdef IS_MPI
207
208	tim	1733	void SimCreator::divideMolecules(SimInfo *info) {
209	tim	1727	double numerator;
210			double denominator;
211			double precast;
212			double x;
213			double y;
214			double a;
215			int old_atoms;
216			int add_atoms;
217			int new_atoms;
218			int nTarget;
219			int done;
220			int i;
221			int j;
222			int loops;
223			int which_proc;
224	tim	1733	int nProcessors;
225			int nGlobalMols;
226			std::vector < int > atomsPerProc;
227			randomSPRNG myRandom(info->getSeed());
228			int * molToProcMap = info->getMolToProcMapPointer();
229	tim	1727
230	tim	1733	MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
231			nGlobalMols = info->getNGlobalMolecules();
232	tim	1727
233	tim	1733	if (nProcessors > nGlobalMols) {
234	tim	1727	sprintf(painCave.errMsg,
235			"nProcessors (%d) > nMol (%d)\n"
236			"\tThe number of processors is larger than\n"
237			"\tthe number of molecules. This will not result in a \n"
238			"\tusable division of atoms for force decomposition.\n"
239			"\tEither try a smaller number of processors, or run the\n"
240	tim	1733	"\tsingle-processor version of OOPSE.\n", nProcessors, nGlobalMols);
241	tim	1727
242			painCave.isFatal = 1;
243			simError();
244			}
245
246	tim	1733	a = 3.0 * nGlobalMols / info->getNGlobalAtoms();
247	tim	1727
248	tim	1733	//initialize atomsPerProc
249			atomsPerProc.insert(atomsPerProc.end(), nProcessors, 0);
250	tim	1727
251	tim	1733	for(i = 0; i < nGlobalMols; i++) {
252	tim	1727	// default to an error condition:
253	tim	1733	molToProcMap[i] = -1;
254	tim	1727	}
255
256	tim	1733	if (worldRank == 0) {
257			numerator = info->getNGlobalAtoms();
258			denominator = nProcessors;
259	tim	1727	precast = numerator / denominator;
260			nTarget = (int)(precast + 0.5);
261
262	tim	1733	for(i = 0; i < nGlobalMols; i++) {
263	tim	1727	done = 0;
264			loops = 0;
265
266			while (!done) {
267			loops++;
268
269			// Pick a processor at random
270
271	tim	1733	which_proc = (int) (myRandom.getRandom() * nProcessors);
272	tim	1727
273	tim	1733	//get the molecule stamp first
274			int stampId = info->getMoleculeStampId(i);
275			MoleculeStamp * moleculeStamp = info->getMoleculeStamp(stampId);
276	tim	1727
277	tim	1733	// How many atoms does this processor have so far?
278			old_atoms = atomsPerProc[which_proc];
279			add_atoms = moleculeStamp->getNAtoms();
280	tim	1727	new_atoms = old_atoms + add_atoms;
281
282			// If we've been through this loop too many times, we need
283			// to just give up and assign the molecule to this processor
284			// and be done with it.
285
286			if (loops > 100) {
287	tim	1733	sprintf(painCave.errMsg,
288			"I've tried 100 times to assign molecule %d to a "
289			" processor, but can't find a good spot.\n"
290			"I'm assigning it at random to processor %d.\n",
291			i, which_proc);
292	tim	1727
293			painCave.isFatal = 0;
294			simError();
295
296	tim	1733	molToProcMap[i] = which_proc;
297			atomsPerProc[which_proc] += add_atoms;
298	tim	1727
299			done = 1;
300			continue;
301			}
302
303			// If we can add this molecule to this processor without sending
304			// it above nTarget, then go ahead and do it:
305
306			if (new_atoms <= nTarget) {
307	tim	1733	molToProcMap[i] = which_proc;
308			atomsPerProc[which_proc] += add_atoms;
309	tim	1727
310			done = 1;
311			continue;
312			}
313
314			// The only situation left is when new_atoms > nTarget. We
315			// want to accept this with some probability that dies off the
316			// farther we are from nTarget
317
318			// roughly: x = new_atoms - nTarget
319			// Pacc(x) = exp(- a * x)
320			// where a = penalty / (average atoms per molecule)
321
322			x = (double)(new_atoms - nTarget);
323	tim	1733	y = myRandom.getRandom();
324	tim	1727
325			if (y < exp(- a * x)) {
326	tim	1733	molToProcMap[i] = which_proc;
327			atomsPerProc[which_proc] += add_atoms;
328	tim	1727
329			done = 1;
330			continue;
331	tim	1733	} else {
332			continue;
333			}
334	tim	1727	}
335			}
336
337			// Spray out this nonsense to all other processors:
338
339	tim	1733	MPI_Bcast(molToProcMap, nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
340	tim	1727	} else {
341
342	tim	1733	// Listen to your marching orders from processor 0:
343	tim	1727
344	tim	1733	MPI_Bcast(molToProcMap, nGlobalMols, MPI_INT, 0, MPI_COMM_WORLD);
345	tim	1727	}
346
347			sprintf(checkPointMsg,
348			"Successfully divided the molecules among the processors.\n");
349			MPIcheckPoint();
350	tim	1712	}
351	tim	1733
352	tim	1713	#endif
353	tim	1712
354	tim	1733	void SimCreator::createMolecules(SimInfo *info) {
355	tim	1725	MoleculeCreator molCreator;
356			int stampId;
357	tim	1713
358	tim	1733	for(int i = 0; i < info->getNGlobalMolecules(); i++) {
359	tim	1713
360	tim	1733	#ifdef IS_MPI
361
362			if (info->getMolToProc(i) == worldRank) {
363			#endif
364
365	tim	1725	stampId = info->getMoleculeStampId(i);
366	tim	1733	Molecule * mol = molCreator.createMolecule(info->getForceField(),
367			info->getMoleculeStamp(stampId), stampId, i);
368
369	tim	1719	info->addMolecule(mol);
370	tim	1733
371			#ifdef IS_MPI
372
373	tim	1719	}
374	tim	1733
375			#endif
376
377			} //end for(int i=0)
378	tim	1713	}
379
380	tim	1733	void SimSetup::compList(MakeStamps stamps, Globals globals,
381			std::vector < std::pair<MoleculeStamp *, int> > &moleculeStampPairs) {
382	tim	1719	int i;
383	tim	1733	char * id;
384			MoleculeStamp * currentStamp;
385			Component * the_components = globals->getComponents();
386			int n_components = globals->getNComponents();
387	tim	1719
388	tim	1733	if (!globals->haveNMol()) {
389	tim	1719	// we don't have the total number of molecules, so we assume it is
390			// given in each component
391
392	tim	1733	for(i = 0; i < n_components; i++) {
393			if (!the_components[i]->haveNMol()) {
394	tim	1719	// we have a problem
395			sprintf(painCave.errMsg,
396	tim	1733	"SimSetup Error. No global NMol or component NMol given.\n"
397			"\tCannot calculate the number of atoms.\n");
398
399	tim	1719	painCave.isFatal = 1;
400			simError();
401			}
402
403			id = the_components[i]->getType();
404			currentStamp = stamps->extractMolStamp(id);
405	tim	1733
406			if (currentStamp == NULL) {
407	tim	1719	sprintf(painCave.errMsg,
408	tim	1733	"SimSetup error: Component \"%s\" was not found in the "
409			"list of declared molecules\n", id);
410
411	tim	1719	painCave.isFatal = 1;
412			simError();
413	tim	1733	}
414	tim	1719
415	tim	1733	moleculeStampPairs.push_back(
416			make_pair(currentStamp, the_components[i]->getNMol));
417	tim	1719	} //end for (i = 0; i < n_components; i++)
418	tim	1733	} else {
419			sprintf(painCave.errMsg, "SimSetup error.\n"
420			"\tSorry, the ability to specify total"
421			" nMols and then give molfractions in the components\n"
422			"\tis not currently supported."
423			" Please give nMol in the components.\n");
424
425	tim	1719	painCave.isFatal = 1;
426			simError();
427			}
428	tim	1733
429	tim	1719	#ifdef IS_MPI
430	tim	1733
431			strcpy(checkPointMsg, "Component stamps successfully extracted\n");
432			MPIcheckPoint();
433
434	tim	1719	#endif // is_mpi
435	tim	1733
436	tim	1719	}
437
438	tim	1733	void SimCreator::setGlobalIndex(SimInfo *info) {
439	tim	1725	typename SimInfo::MoleculeIterator mi;
440			typename Molecule::AtomIterator ai;
441			typename Molecule::RigidBodyIterator ri;
442			typename Molecule::CutoffGroupIterator ci;
443	tim	1733	Molecule * mol;
444			Atom * atom;
445			RigidBody * rb;
446			CutoffGroup * cg;
447	tim	1725	int beginAtomIndex;
448			int beginRigidBodyIndex;
449			int beginCutoffGroupIndex;
450	tim	1719
451	tim	1733	#ifndef IS_MPI
452
453	tim	1725	beginAtomIndex = 0;
454			beginRigidBodyIndex = 0;
455			beginCutoffGroupIndex = 0;
456	tim	1733
457	tim	1725	#else
458	tim	1733
459	tim	1725	int nproc;
460			int myNode;
461
462			myNode = worldRank;
463			MPI_Comm_size(MPI_COMM_WORLD, &nproc);
464
465	tim	1733	std::vector < int > tmpAtomsInProc(nproc, 0);
466			std::vector < int > tmpRigidBodiesInProc(nproc, 0);
467			std::vector < int > tmpCutoffGroupsInProc(nproc, 0);
468			std::vector < int > NumAtomsInProc(nproc, 0);
469			std::vector < int > NumRigidBodiesInProc(nproc, 0);
470			std::vector < int > NumCutoffGroupsInProc(nproc, 0);
471
472			tmpAtomsInProc[myNode] = info->getNAtoms();
473			tmpRigidBodiesInProc[myNode] = info->getNRigidBodies();
474			tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroups();
475
476	tim	1725	//do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
477	tim	1733	MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,
478			MPI_SUM, MPI_COMM_WORLD);
479			MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc,
480			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
481			MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc,
482			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
483	tim	1725
484			beginAtomIndex = 0;
485			beginRigidBodyIndex = 0;
486			beginCutoffGroupIndex = 0;
487
488	tim	1733	for(int i = 0; i < nproc; i++) {
489	tim	1725	beginAtomIndex += NumAtomsInProc[i];
490			beginRigidBodyIndex += NumRigidBodiesInProc[i];
491	tim	1733	beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
492	tim	1725	}
493	tim	1733
494	tim	1725	#endif
495
496	tim	1733	for(mol = info->beginMolecule(mi); mol != NULL;
497			mol = info->nextMolecule(mi)) {
498
499	tim	1725	//local index(index in DataStorge) of atom is important
500	tim	1733	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
501	tim	1725	atom->setGlobalIndex(beginAtomIndex++);
502			}
503
504	tim	1733	for(rb = mol->beginRigidBody(ri); rb != NULL;
505			rb = mol->nextRigidBody(ri)) {
506	tim	1725	rb->setGlobalIndex(beginRigidBodyIndex++);
507			}
508	tim	1733
509	tim	1725	//local index of cutoff group is trivial, it only depends on the order of travesing
510	tim	1733	for(cg = mol->beginCutoffGroup(ci); cg != NULL;
511			cg = mol->nextCutoffGroup(ci)) {
512			cg->setGlobalIndex(beginCutoffGroupIndex++);
513			}
514			}
515
516	tim	1735	//fill globalGroupMembership
517	tim	1733	std::vector<int> globalGroupMembership(info->getNGlobalAtoms(), 0);
518			for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
519	tim	1725	for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
520	tim	1733
521			for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
522			globalGroupMembership[atom->getGlobalIndex()] = cg->getGlobalIndex();
523			}
524
525			}
526			}
527	tim	1735
528			#ifdef IS_MPI
529	tim	1733	// Since the globalGroupMembership has been zero filled and we've only
530			// poked values into the atoms we know, we can do an Allreduce
531			// to get the full globalGroupMembership array (We think).
532			// This would be prettier if we could use MPI_IN_PLACE like the MPI-2
533			// docs said we could.
534
535			MPI_Allreduce(&globalGroupMembership[0],
536			info->getGlobalGroupMembershipPointer(),
537			info->getNGlobalAtoms(),
538			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
539	tim	1735	#else
540			std::copy(globalGroupMembership.begin(), globalGroupMembership.end(),
541			info->getGlobalGroupMembershipPointer());
542			#endif
543	tim	1733
544			//fill molMembership
545			std::vector<int> globalMolMembership(info->getNGlobalAtoms(), 0);
546
547			for(mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
548
549			for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
550			globalMolMembership[atom->getGlobalIndex()] = mol->getGlobalIndex();
551			}
552
553	tim	1735	#ifdef IS_MPI
554	tim	1733	MPI_Allreduce(&globalMolMembership[0],
555			info->getGlobalMolMembershipPointer(),
556			info->getNGlobalAtoms(),
557			MPI_INT, MPI_SUM, MPI_COMM_WORLD);
558	tim	1735	#else
559			std::copy(globalMolMembership.begin(), globalMolMembership.end(),
560			info->getGlobalMolMembershipPointer());
561			#endif
562
563	tim	1733	}
564
565
566			void SimCreator::initFortran(SimInfo* info) {
567	tim	1735
568			//setup fortran simulation (Actually, in parallel version it will initialze the parallel part of fortran first)
569	tim	1733	info->update();
570
571	tim	1735	//notify fortran whether reaction field is used or not
572			//deprecated
573			int isError;
574			isError = 0;
575			initFortranFF( &useReactionField, &isError );
576	tim	1733
577	tim	1735	if(isError){
578			sprintf( painCave.errMsg,
579			"SimCreator::initFortran() error: There was an error initializing the forceField in fortran.\n" );
580			painCave.isFatal = 1;
581			simError();
582			}
583
584			//figure out the cutoff radius and pass it to fortran
585			setCutoff(info);
586			}
587	tim	1733
588	tim	1735	void SimCreator::setCutoff(SimInfo* info) {
589			Globals* globals = info->getGlobals();
590			double rcut; //cutoff radius
591			double rsw; //switching radius
592	tim	1733
593	tim	1735	simtype* st = info->getSimType();
594
595			if (st->SIM_uses_Charges \| st->SIM_uses_Dipoles \| st->SIM_uses_RF) {
596
597			if (!globals->haveRcut()){
598			sprintf(painCave.errMsg,
599			"SimCreator Warning: No value was set for the cutoffRadius.\n"
600			"\tOOPSE will use a default value of 15.0 angstroms"
601			"\tfor the cutoffRadius.\n");
602			painCave.isFatal = 0;
603			simError();
604			rcut = 15.0;
605			} else{
606			rcut = globals->getRcut();
607			}
608	tim	1733
609	tim	1735	if (!globals->haveRsw()){
610			sprintf(painCave.errMsg,
611			"SimCreator Warning: No value was set for switchingRadius.\n"
612			"\tOOPSE will use a default value of\n"
613			"\t0.95 * cutoffRadius for the switchingRadius\n");
614			painCave.isFatal = 0;
615			simError();
616			rsw = 0.95 * rcut;
617			} else{
618			rsw = globals->getRsw();
619			}
620	tim	1733
621	tim	1735	} else {
622			// if charge, dipole or reaction field is not used and the cutofff radius is not specified in
623			//meta-data file, the maximum cutoff radius calculated from forcefiled will be used
624
625			if (globals->haveRcut()) {
626			rcut = globals->getRcut();
627			} else {
628			//set cutoff radius to the maximum cutoff radius based on atom types in the whole system
629			rcut = info->calcMaxCutoffRadius();
630			}
631	tim	1733
632	tim	1735	if (globals->haveRsw()) {
633			rsw = globals->getRsw()
634			} else {
635			rsw = rcut;
636			}
637	tim	1733
638	tim	1735	}
639	tim	1733
640	tim	1735	//store them into Siminfo
641			info->setRcut(rcut);
642			info->setRsw(rsw);
643
644			double rnblist = rcut + 1; // skin of neighbor list
645
646			//Pass these cutoff radius etc. to fortran. This function should be called once and only once
647			notifyFortranCutoffs(&rcut, &rsw, &rnblist);
648	tim	1733	}
649
650			void SimCreator::loadCoordinates(SimInfo* info) {
651			Globals* globals;
652			globals = info->getGlobals();
653
654			if (!globals->haveInitialConfig()) {
655			sprintf(painCave.errMsg,
656			"Cannot intialize a simulation without an initial configuration file.\n");
657			painCave.isFatal = 1;;
658			simError();
659			}
660	tim	1725
661	tim	1733	DumpReader reader(globals->getInitialConfig());
662			int nframes = reader->getNframes();
663
664			if (nframes > 0) {
665			reader.readFrame(info, nframes - 1);
666			} else {
667	tim	1735	//invalid initial coordinate file
668	tim	1733	sprintf(painCave.errMsg, "Initial configuration file %s should at least contain one frame\n",
669			globals->getInitialConfig());
670			painCave.isFatal = 1;
671			simError();
672	tim	1725	}
673	tim	1733
674	tim	1725	}
675
676	tim	1703	} //end namespace oopse