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(char* mdfile,  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(mdfile);

#ifdef IS_MPI
    throwMPIEvent(NULL);
  }
  else{
  set_interface_stamps(stamps, globals);
  mpiEventInit();
  MPIcheckPoint();
  mpiEventLoop();
  }
#endif

}


SimInfo* SimCreator::createSim(const std::string& mdfile) {
    MakeStamps* stamps = new MakeStamps();
    
    Globals* globals = new Globals();

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

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

    //extract the molecule stamps and add them into SimInfo
    compList(stamps, info);
    
    //gather parameters (SimCreator only retrieves the parameters which will be used to create
    // the simulation)
    gatherParameters(info);
    
    //divide the molecules and determine the global index of molecules
    
    //create the molecules
    createMolecules(info);

    //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);
    
    //load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
    //eta, chi for NPT integrator)
    DumpReader reader(info->getInitFilename());
    int nframes = reader->getNframes();

    if (nframes > 0) {
        reader.readFrame(info, nframes - 1);
    } else {
        //invalid initial coordinate file
        
    }

    
    //initialize fortran    
    
    return info;
}    

void SimCreator::gatherParameters(SimInfo* info) {
    //model->addProperty(new StringGenericData("Ensemble", globals->getForceFiled()));
    //model->addProperty(new DoubleGenericData("dt"), globals->getDt());

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

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

    //using seed generated by system instead of invalid seed set by user 
#ifndef IS_MPI
            seedValue = make_sprng_seed();
#else 
            if (worldRank == 0){
                seedValue = make_sprng_seed();
            }
            MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);   
#endif      
        } //end if (seedValue /1000000000 == 0)
    } else{
    
#ifndef IS_MPI
        seedValue = make_sprng_seed();
#else 
        if (worldRank == 0){
            seedValue = make_sprng_seed();
        }
        MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);   
#endif
    }//end of globals->haveSeed()
    info->setSeed(seedValue);

    //
    std::string prefix;
#ifdef IS_MPI
    if (worldRank == 0){
#endif // is_mpi
      
        if(globals->haveFinalConfig()) {
            prefix = getPrefix(globals->getFinalConfig());   
        } else {
            prefix = getPrefix(mdfile_);
        }
        
        info->setFinalConfFileName(prefix + ".eor");
        info->setDumpFileName(prefix + ".dump");
        info->setStatFileName(prefix + ".stat");

#ifdef IS_MPI
    }
#endif    


}

#ifdef IS_MPI
void SimCreator::mpiMolDivide(){
 
    mpiSim = new mpiSimulation(info);

    mpiSim->divideLabor();

}
#endif

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

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

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


void SimSetup::compList(MakeStamps* stamps,SimInfo* info) {
    int i;
    char* id;
    MoleculeStamp* currentStamp;
    Component* the_components = info->getGlobals()->getComponents();
    int n_components = info->getGlobals()->getNComponents();

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

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

            id = the_components[i]->getType();
            currentStamp = stamps->extractMolStamp(id);
            if (currentStamp == NULL){
                sprintf(painCave.errMsg,
                    "SimSetup error: Component \"%s\" was not found in the "
                    "list of declared molecules\n",
                    id);
                painCave.isFatal = 1;
                simError();
            }      

            info->addMoleculeStamp(currentStamp, the_components[i]->getNMol);
          
        } //end for (i = 0; i < n_components; i++)
        
    } else{
        sprintf(painCave.errMsg,
                "SimSetup error.\n"
                "\tSorry, the ability to specify total"
                " nMols and then give molfractions in the components\n"
                "\tis not currently supported."
                " Please give nMol in the components.\n");
        painCave.isFatal = 1;
        simError();
    }
  
#ifdef IS_MPI
  strcpy(checkPointMsg, "Component stamps successfully extracted\n");
  MPIcheckPoint();
#endif // is_mpi
}

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

#ifdef IS_MPI
    beginAtomIndex = 0;
    beginRigidBodyIndex = 0;
    beginCutoffGroupIndex = 0;
#else
    int nproc;
    int myNode;

    myNode = worldRank;
    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
    
    std::vector<int> tmpAtomsInProc(nproc, 0);
    std::vector<int> tmpRigidBodiesInProc(nproc, 0);
    std::vector<int> tmpCutoffGroupsInProc(nproc, 0);    
    std::vector<int> NumAtomsInProc(nproc, 0);
    std::vector<int> NumRigidBodiesInProc(nproc, 0);
    std::vector<int> NumCutoffGroupsInProc(nproc, 0);

    tmpAtomInProc[myNode] = info->getNAtoms();
    tmpRigidBodiesInProc[myNode] = info->getNRigidBodiess();
    tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroupss();
    
    //do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
    MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
    MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc, MPI_INT,MPI_SUM, MPI_COMM_WORLD);

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

    for (int i = 0; i < nproc; i++) {
        beginAtomIndex += NumAtomsInProc[i];
        beginRigidBodyIndex += NumRigidBodiesInProc[i];
        beginCutoffGroupIndex += NumCutoffGroupsInProc[i];        
    }
    
#endif
    
    for (mol = info->beginMolecule(mi); mol != NULL; mol  = info->nextMolecule(mi)) {

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

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

} //end namespace oopse
Revision:	1725
Committed:	Wed Nov 10 22:01:06 2004 UTC (19 years, 8 months ago) by tim
File size:	10375 byte(s)
Log Message:	another painful day (1) SimCreator, SimInfo, mpiSimulation (2) DumpReader, DumpWriter (InitializeFrom File will be removed) (3) ForceField (at least LJ) and BondType, BendType, TorsionType (4)Integrator (5)oopse.cpp (6)visitors & Dump2XYZ (7)SimpleBuilder (8)Constraint & ZConstraint
#	Content
1	/*
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	/**
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	namespace oopse {
36
37	void SimSetup::parseFile(char* mdfile, MakeStamps* stamps, Globals* globals){
38	#ifdef IS_MPI
39	if (worldRank == 0){
40	#endif // is_mpi
41
42
43	globals->initalize();
44	set_interface_stamps(stamps, globals);
45
46	#ifdef IS_MPI
47	mpiEventInit();
48	#endif
49
50	yacc_BASS(mdfile);
51
52	#ifdef IS_MPI
53	throwMPIEvent(NULL);
54	}
55	else{
56	set_interface_stamps(stamps, globals);
57	mpiEventInit();
58	MPIcheckPoint();
59	mpiEventLoop();
60	}
61	#endif
62
63	}
64
65
66	SimInfo* SimCreator::createSim(const std::string& mdfile) {
67	MakeStamps* stamps = new MakeStamps();
68
69	Globals* globals = new Globals();
70
71	//parse meta-data file
72	parseFile(mdfile, stamps, globals);
73
74	//create the force field
75	ForceFiled* ff = ForceFieldFactory::getInstance()->createObject(globals->getForceField());
76
77	//create SimInfo
78	SimInfo* info = new SimInfo();
79	info->setGlobals(globals);
80	info->setForceField(ff);
81
82	//extract the molecule stamps and add them into SimInfo
83	compList(stamps, info);
84
85	//gather parameters (SimCreator only retrieves the parameters which will be used to create
86	// the simulation)
87	gatherParameters(info);
88
89	//divide the molecules and determine the global index of molecules
90
91	//create the molecules
92	createMolecules(info);
93
94	//set the global index of atoms, rigidbodies and cutoffgroups (only need to be set once, the
95	//global index will never change again). Local indices of atoms and rigidbodies are already set by
96	//MoleculeCreator class which actually delegates the responsibility to LocalIndexManager.
97	setGlobalIndices(info);
98
99	//allocate memory for DataStorage(circular reference, need to break it)
100	info->setSnapshotManager(new SimSnapshotManager(info);
101
102	//load initial coordinates, some extra information are pushed into SimInfo's property map ( such as
103	//eta, chi for NPT integrator)
104	DumpReader reader(info->getInitFilename());
105	int nframes = reader->getNframes();
106
107	if (nframes > 0) {
108	reader.readFrame(info, nframes - 1);
109	} else {
110	//invalid initial coordinate file
111
112	}
113
114
115	//initialize fortran
116
117	return info;
118	}
119
120	void SimCreator::gatherParameters(SimInfo* info) {
121	//model->addProperty(new StringGenericData("Ensemble", globals->getForceFiled()));
122	//model->addProperty(new DoubleGenericData("dt"), globals->getDt());
123
124	//setup seed for random number generator
125	int seedValue;
126	Globals* globals = info->getGlobals();
127
128	if (globals->haveSeed()){
129	seedValue = globals->getSeed();
130
131	if(seedValue / 1000000000 == 0){
132	sprintf(painCave.errMsg,
133	"Seed for sprng library should contain at least 9 digits\n"
134	"OOPSE will generate a seed for user\n");
135	painCave.isFatal = 0;
136	simError();
137
138	//using seed generated by system instead of invalid seed set by user
139	#ifndef IS_MPI
140	seedValue = make_sprng_seed();
141	#else
142	if (worldRank == 0){
143	seedValue = make_sprng_seed();
144	}
145	MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
146	#endif
147	} //end if (seedValue /1000000000 == 0)
148	} else{
149
150	#ifndef IS_MPI
151	seedValue = make_sprng_seed();
152	#else
153	if (worldRank == 0){
154	seedValue = make_sprng_seed();
155	}
156	MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);
157	#endif
158	}//end of globals->haveSeed()
159	info->setSeed(seedValue);
160
161	//
162	std::string prefix;
163	#ifdef IS_MPI
164	if (worldRank == 0){
165	#endif // is_mpi
166
167	if(globals->haveFinalConfig()) {
168	prefix = getPrefix(globals->getFinalConfig());
169	} else {
170	prefix = getPrefix(mdfile_);
171	}
172
173	info->setFinalConfFileName(prefix + ".eor");
174	info->setDumpFileName(prefix + ".dump");
175	info->setStatFileName(prefix + ".stat");
176
177	#ifdef IS_MPI
178	}
179	#endif
180
181
182
183	}
184
185	#ifdef IS_MPI
186	void SimCreator::mpiMolDivide(){
187
188	mpiSim = new mpiSimulation(info);
189
190	mpiSim->divideLabor();
191
192	}
193	#endif
194
195	Molecule* SimCreator::createMolecules(SimInfo* info) {
196	MoleculeCreator molCreator;
197	int stampId;
198
199
200	for (int i = 0; i < info->getNGlobalMolecules(); i++){
201	if (info->getMolToProc(i) == worldRank){
202
203	stampId = info->getMoleculeStampId(i);
204	Molecule* mol = molCreator.createMolecule(info->getForceField(),
205	info->getMoleculeStamp(stampId), stampId, i);
206
207	info->addMolecule(mol);
208	}
209	}
210
211	}
212
213
214	void SimSetup::compList(MakeStamps* stamps,SimInfo* info) {
215	int i;
216	char* id;
217	MoleculeStamp* currentStamp;
218	Component* the_components = info->getGlobals()->getComponents();
219	int n_components = info->getGlobals()->getNComponents();
220
221	if (!globals->haveNMol()){
222	// we don't have the total number of molecules, so we assume it is
223	// given in each component
224
225	for (i = 0; i < n_components; i++){
226	if (!the_components[i]->haveNMol()){
227	// we have a problem
228	sprintf(painCave.errMsg,
229	"SimSetup Error. No global NMol or component NMol given.\n"
230	"\tCannot calculate the number of atoms.\n");
231	painCave.isFatal = 1;
232	simError();
233	}
234
235	id = the_components[i]->getType();
236	currentStamp = stamps->extractMolStamp(id);
237	if (currentStamp == NULL){
238	sprintf(painCave.errMsg,
239	"SimSetup error: Component \"%s\" was not found in the "
240	"list of declared molecules\n",
241	id);
242	painCave.isFatal = 1;
243	simError();
244	}
245
246	info->addMoleculeStamp(currentStamp, the_components[i]->getNMol);
247
248	} //end for (i = 0; i < n_components; i++)
249
250	} else{
251	sprintf(painCave.errMsg,
252	"SimSetup error.\n"
253	"\tSorry, the ability to specify total"
254	" nMols and then give molfractions in the components\n"
255	"\tis not currently supported."
256	" Please give nMol in the components.\n");
257	painCave.isFatal = 1;
258	simError();
259	}
260
261	#ifdef IS_MPI
262	strcpy(checkPointMsg, "Component stamps successfully extracted\n");
263	MPIcheckPoint();
264	#endif // is_mpi
265	}
266
267	void SimCreator::setGlobalIndices(SimInfo* info) {
268	typename SimInfo::MoleculeIterator mi;
269	typename Molecule::AtomIterator ai;
270	typename Molecule::RigidBodyIterator ri;
271	typename Molecule::CutoffGroupIterator ci;
272	Molecule* mol;
273	Atom* atom;
274	RigidBody* rb;
275	CutoffGroup* cg;
276	int beginAtomIndex;
277	int beginRigidBodyIndex;
278	int beginCutoffGroupIndex;
279
280	#ifdef IS_MPI
281	beginAtomIndex = 0;
282	beginRigidBodyIndex = 0;
283	beginCutoffGroupIndex = 0;
284	#else
285	int nproc;
286	int myNode;
287
288	myNode = worldRank;
289	MPI_Comm_size(MPI_COMM_WORLD, &nproc);
290
291	std::vector<int> tmpAtomsInProc(nproc, 0);
292	std::vector<int> tmpRigidBodiesInProc(nproc, 0);
293	std::vector<int> tmpCutoffGroupsInProc(nproc, 0);
294	std::vector<int> NumAtomsInProc(nproc, 0);
295	std::vector<int> NumRigidBodiesInProc(nproc, 0);
296	std::vector<int> NumCutoffGroupsInProc(nproc, 0);
297
298	tmpAtomInProc[myNode] = info->getNAtoms();
299	tmpRigidBodiesInProc[myNode] = info->getNRigidBodiess();
300	tmpCutoffGroupsInProc[myNode] = info->getNCutoffGroupss();
301
302	//do MPI_ALLREDUCE to exchange the total number of atoms, rigidbodies and cutoff groups
303	MPI_Allreduce(&tmpAtomsInProc[0], &NumAtomsInProc[0], nproc, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
304	MPI_Allreduce(&tmpRigidBodiesInProc[0], &NumRigidBodiesInProc[0], nproc, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
305	MPI_Allreduce(&tmpCutoffGroupsInProc[0], &NumCutoffGroupsInProc[0], nproc, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
306
307	beginAtomIndex = 0;
308	beginRigidBodyIndex = 0;
309	beginCutoffGroupIndex = 0;
310
311	for (int i = 0; i < nproc; i++) {
312	beginAtomIndex += NumAtomsInProc[i];
313	beginRigidBodyIndex += NumRigidBodiesInProc[i];
314	beginCutoffGroupIndex += NumCutoffGroupsInProc[i];
315	}
316
317	#endif
318
319	for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
320
321	//local index(index in DataStorge) of atom is important
322	for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
323	atom->setGlobalIndex(beginAtomIndex++);
324	}
325
326	for (rb = mol->beginRigidBody(ri); rb != NULL; rb = mol->nextRigidBody(ri)) {
327	rb->setGlobalIndex(beginRigidBodyIndex++);
328	}
329
330	//local index of cutoff group is trivial, it only depends on the order of travesing
331	for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
332	cg->setGlobalIndex(beginCutoffGroupIndex++);
333	}
334
335	}
336	}
337
338	} //end namespace oopse