# | Line 6 | Line 6 | |
---|---|---|
6 | #include "parse_me.h" | |
7 | #include "LRI.hpp" | |
8 | #include "Integrator.hpp" | |
9 | + | #include "simError.h" |
10 | ||
11 | + | #ifdef IS_MPI |
12 | + | #include "mpiBASS.h" |
13 | + | #include "bassDiag.hpp" |
14 | + | #endif |
15 | + | |
16 | SimSetup::SimSetup(){ | |
17 | stamps = new MakeStamps(); | |
18 | globals = new Globals(); | |
19 | + | |
20 | + | #ifdef IS_MPI |
21 | + | strcpy( checkPointMsg, "SimSetup creation successful" ); |
22 | + | MPIcheckPoint(); |
23 | + | #endif // IS_MPI |
24 | } | |
25 | ||
26 | SimSetup::~SimSetup(){ | |
# | Line 19 | Line 30 | void SimSetup::parseFile( char* fileName ){ | |
30 | ||
31 | void SimSetup::parseFile( char* fileName ){ | |
32 | ||
33 | < | inFileName = fileName; |
34 | < | set_interface_stamps( stamps, globals ); |
35 | < | yacc_BASS( fileName ); |
33 | > | #ifdef IS_MPI |
34 | > | if( worldRank == 0 ){ |
35 | > | #endif // is_mpi |
36 | > | |
37 | > | inFileName = fileName; |
38 | > | set_interface_stamps( stamps, globals ); |
39 | > | |
40 | > | #ifdef IS_MPI |
41 | > | mpiEventInit(); |
42 | > | #endif |
43 | > | |
44 | > | yacc_BASS( fileName ); |
45 | > | |
46 | > | #ifdef IS_MPI |
47 | > | throwMPIEvent(NULL); |
48 | > | } |
49 | > | else receiveParse(); |
50 | > | #endif |
51 | > | |
52 | } | |
53 | ||
54 | + | #ifdef IS_MPI |
55 | + | void SimSetup::receiveParse(void){ |
56 | + | |
57 | + | set_interface_stamps( stamps, globals ); |
58 | + | mpiEventInit(); |
59 | + | MPIcheckPoint(); |
60 | + | mpiEventLoop(); |
61 | + | |
62 | + | } |
63 | + | |
64 | + | |
65 | + | void SimSetup::testMe(void){ |
66 | + | bassDiag* dumpMe = new bassDiag(globals,stamps); |
67 | + | dumpMe->dumpStamps(); |
68 | + | delete dumpMe; |
69 | + | } |
70 | + | #endif |
71 | + | |
72 | void SimSetup::createSim( void ){ | |
73 | ||
74 | MakeStamps *the_stamps; | |
75 | Globals* the_globals; | |
76 | < | int i; |
76 | > | int i, j; |
77 | ||
78 | // get the stamps and globals; | |
79 | the_stamps = stamps; | |
# | Line 43 | Line 88 | void SimSetup::createSim( void ){ | |
88 | n_components = the_globals->getNComponents(); | |
89 | strcpy( force_field, the_globals->getForceField() ); | |
90 | strcpy( ensemble, the_globals->getEnsemble() ); | |
91 | < | |
91 | > | |
92 | if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF(); | |
93 | else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF(); | |
94 | else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF(); | |
95 | else{ | |
96 | < | std::cerr<< "SimSetup Error. Unrecognized force field -> " |
97 | < | << force_field << "\n"; |
98 | < | exit(8); |
96 | > | sprintf( painCave.errMsg, |
97 | > | "SimSetup Error. Unrecognized force field -> %s\n", |
98 | > | force_field ); |
99 | > | painCave.isFatal = 1; |
100 | > | simError(); |
101 | } | |
102 | ||
103 | + | #ifdef IS_MPI |
104 | + | strcpy( checkPointMsg, "ForceField creation successful" ); |
105 | + | MPIcheckPoint(); |
106 | + | #endif // is_mpi |
107 | + | |
108 | // get the components and calculate the tot_nMol and indvidual n_mol | |
109 | the_components = the_globals->getComponents(); | |
110 | components_nmol = new int[n_components]; | |
111 | comp_stamps = new MoleculeStamp*[n_components]; | |
112 | < | |
112 | > | |
113 | if( !the_globals->haveNMol() ){ | |
114 | < | // we don't have the total number of molecules, so we assume it is |
114 | > | // we don't have the total number of molecules, so we assume it is |
115 | // given in each component | |
116 | ||
117 | tot_nmol = 0; | |
118 | for( i=0; i<n_components; i++ ){ | |
119 | < | |
119 | > | |
120 | if( !the_components[i]->haveNMol() ){ | |
121 | // we have a problem | |
122 | < | std::cerr << "SimSetup Error. No global NMol or component NMol" |
123 | < | << " given. Cannot calculate the number of atoms.\n"; |
124 | < | exit( 8 ); |
122 | > | sprintf( painCave.errMsg, |
123 | > | "SimSetup Error. No global NMol or component NMol" |
124 | > | " given. Cannot calculate the number of atoms.\n" ); |
125 | > | painCave.isFatal = 1; |
126 | > | simError(); |
127 | } | |
128 | ||
129 | tot_nmol += the_components[i]->getNMol(); | |
# | Line 77 | Line 131 | void SimSetup::createSim( void ){ | |
131 | } | |
132 | } | |
133 | else{ | |
134 | < | std::cerr << "NOT A SUPPORTED FEATURE\n"; |
134 | > | sprintf( painCave.errMsg, |
135 | > | "SimSetup error.\n" |
136 | > | "\tSorry, the ability to specify total" |
137 | > | " nMols and then give molfractions in the components\n" |
138 | > | "\tis not currently supported." |
139 | > | " Please give nMol in the components.\n" ); |
140 | > | painCave.isFatal = 1; |
141 | > | simError(); |
142 | ||
82 | – | // tot_nmol = the_globals->getNMol(); |
143 | ||
144 | < | // //we have the total number of molecules, now we check for molfractions |
145 | < | // for( i=0; i<n_components; i++ ){ |
146 | < | |
147 | < | // if( !the_components[i]->haveMolFraction() ){ |
148 | < | |
149 | < | // if( !the_components[i]->haveNMol() ){ |
150 | < | // //we have a problem |
151 | < | // std::cerr << "SimSetup error. Neither molFraction nor " |
152 | < | // << " nMol was given in component |
153 | < | |
144 | > | // tot_nmol = the_globals->getNMol(); |
145 | > | |
146 | > | // //we have the total number of molecules, now we check for molfractions |
147 | > | // for( i=0; i<n_components; i++ ){ |
148 | > | |
149 | > | // if( !the_components[i]->haveMolFraction() ){ |
150 | > | |
151 | > | // if( !the_components[i]->haveNMol() ){ |
152 | > | // //we have a problem |
153 | > | // std::cerr << "SimSetup error. Neither molFraction nor " |
154 | > | // << " nMol was given in component |
155 | > | |
156 | } | |
157 | ||
158 | + | #ifdef IS_MPI |
159 | + | strcpy( checkPointMsg, "Have the number of components" ); |
160 | + | MPIcheckPoint(); |
161 | + | #endif // is_mpi |
162 | + | |
163 | // make an array of molecule stamps that match the components used. | |
164 | + | // also extract the used stamps out into a separate linked list |
165 | ||
166 | + | simnfo->nComponents = n_components; |
167 | + | simnfo->componentsNmol = components_nmol; |
168 | + | simnfo->compStamps = comp_stamps; |
169 | + | simnfo->headStamp = new LinkedMolStamp(); |
170 | + | |
171 | + | char* id; |
172 | + | LinkedMolStamp* headStamp = simnfo->headStamp; |
173 | + | LinkedMolStamp* currentStamp = NULL; |
174 | for( i=0; i<n_components; i++ ){ | |
175 | ||
176 | < | comp_stamps[i] = |
177 | < | the_stamps->getMolecule( the_components[i]->getType() ); |
176 | > | id = the_components[i]->getType(); |
177 | > | comp_stamps[i] = NULL; |
178 | > | |
179 | > | // check to make sure the component isn't already in the list |
180 | > | |
181 | > | comp_stamps[i] = headStamp->match( id ); |
182 | > | if( comp_stamps[i] == NULL ){ |
183 | > | |
184 | > | // extract the component from the list; |
185 | > | |
186 | > | currentStamp = the_stamps->extractMolStamp( id ); |
187 | > | if( currentStamp == NULL ){ |
188 | > | sprintf( painCave.errMsg, |
189 | > | "SimSetup error: Component \"%s\" was not found in the " |
190 | > | "list of declared molecules\n" |
191 | > | id ); |
192 | > | painCave.isFatal = 1; |
193 | > | simError(); |
194 | > | } |
195 | > | |
196 | > | headStamp->add( currentStamp ); |
197 | > | comp_stamps[i] = headStamp->match( id ); |
198 | > | } |
199 | } | |
200 | ||
201 | + | #ifdef IS_MPI |
202 | + | strcpy( checkPointMsg, "Component stamps successfully extracted\n" ); |
203 | + | MPIcheckPoint(); |
204 | + | #endif // is_mpi |
205 | ||
206 | ||
207 | + | |
208 | + | |
209 | // caclulate the number of atoms, bonds, bends and torsions | |
210 | ||
211 | tot_atoms = 0; | |
# | Line 111 | Line 214 | void SimSetup::createSim( void ){ | |
214 | tot_torsions = 0; | |
215 | for( i=0; i<n_components; i++ ){ | |
216 | ||
217 | < | tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); |
218 | < | tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); |
219 | < | tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); |
217 | > | tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms(); |
218 | > | tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds(); |
219 | > | tot_bends += components_nmol[i] * comp_stamps[i]->getNBends(); |
220 | tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions(); | |
221 | } | |
222 | < | |
222 | > | |
223 | tot_SRI = tot_bonds + tot_bends + tot_torsions; | |
224 | < | |
224 | > | |
225 | simnfo->n_atoms = tot_atoms; | |
226 | simnfo->n_bonds = tot_bonds; | |
227 | simnfo->n_bends = tot_bends; | |
228 | simnfo->n_torsions = tot_torsions; | |
229 | simnfo->n_SRI = tot_SRI; | |
230 | + | simnfo->n_mol = tot_nmol; |
231 | ||
128 | – | // create the atom and short range interaction arrays |
232 | ||
233 | < | the_atoms = new Atom*[tot_atoms]; |
234 | < | // the_molecules = new Molecule[tot_nmol]; |
233 | > | #ifdef IS_MPI |
234 | > | |
235 | > | // divide the molecules among processors here. |
236 | ||
237 | + | mpiSimulation* mpiSim = new mpiSimulation( simnfo ); |
238 | ||
239 | < | if( tot_SRI ){ |
240 | < | the_sris = new SRI*[tot_SRI]; |
241 | < | the_excludes = new ex_pair[tot_SRI]; |
239 | > | mpiSim->divideLabor(); |
240 | > | |
241 | > | // set up the local variables |
242 | > | |
243 | > | int localMol, allMol; |
244 | > | int local_atoms, local_bonds, local_bends, local_torsions, local_SRI; |
245 | > | |
246 | > | allMol = 0; |
247 | > | localMol = 0; |
248 | > | local_atoms = 0; |
249 | > | local_bonds = 0; |
250 | > | local_bends = 0; |
251 | > | local_torsions = 0; |
252 | > | for( i=0; i<n_components; i++ ){ |
253 | > | |
254 | > | for( j=0; j<components_nmol[i]; j++ ){ |
255 | > | |
256 | > | if( mpiSim->getMyMolStart() <= allMol && |
257 | > | allMol <= mpiSim->getMyMolEnd() ){ |
258 | > | |
259 | > | local_atoms += comp_stamps[i]->getNAtoms(); |
260 | > | local_bonds += comp_stamps[i]->getNBonds(); |
261 | > | local_bends += comp_stamps[i]->getNBends(); |
262 | > | local_torsions += comp_stamps[i]->getNTorsions(); |
263 | > | localMol++; |
264 | > | } |
265 | > | allMol++; |
266 | > | } |
267 | } | |
268 | + | local_SRI = local_bonds + local_bends + local_torsions; |
269 | + | |
270 | ||
271 | + | simnfo->n_atoms = mpiSim->getMyNlocal(); |
272 | + | |
273 | + | if( local_atoms != simnfo->n_atoms ){ |
274 | + | sprintf( painCave.errMsg, |
275 | + | "SimSetup error: mpiSim's localAtom (%d) and SimSetup's" |
276 | + | " localAtom (%d) are note equal.\n", |
277 | + | simnfo->n_atoms, |
278 | + | local_atoms ); |
279 | + | painCave.isFatal = 1; |
280 | + | simError(); |
281 | + | } |
282 | + | |
283 | + | simnfo->n_bonds = local_bonds; |
284 | + | simnfo->n_bends = local_bends; |
285 | + | simnfo->n_torsions = local_torsions; |
286 | + | simnfo->n_SRI = local_SRI; |
287 | + | simnfo->n_mol = localMol; |
288 | + | |
289 | + | strcpy( checkPointMsg, "Passed nlocal consistency check." ); |
290 | + | MPIcheckPoint(); |
291 | + | |
292 | + | |
293 | + | #endif // is_mpi |
294 | + | |
295 | + | |
296 | + | // create the atom and short range interaction arrays |
297 | + | |
298 | + | Atom::createArrays(simnfo->n_atoms); |
299 | + | the_atoms = new Atom*[simnfo->n_atoms]; |
300 | + | the_molecules = new Molecule[simnfo->n_mol]; |
301 | + | |
302 | + | |
303 | + | if( simnfo->n_SRI ){ |
304 | + | the_sris = new SRI*[simnfo->n_SRI]; |
305 | + | the_excludes = new ex_pair[simnfo->n_SRI]; |
306 | + | } |
307 | + | |
308 | // set the arrays into the SimInfo object | |
309 | ||
310 | simnfo->atoms = the_atoms; | |
# | Line 143 | Line 312 | void SimSetup::createSim( void ){ | |
312 | simnfo->n_exclude = tot_SRI; | |
313 | simnfo->excludes = the_excludes; | |
314 | ||
315 | + | |
316 | // initialize the arrays | |
317 | < | |
317 | > | |
318 | the_ff->setSimInfo( simnfo ); | |
319 | < | |
319 | > | |
320 | makeAtoms(); | |
321 | ||
322 | if( tot_bonds ){ | |
# | Line 161 | Line 331 | void SimSetup::createSim( void ){ | |
331 | makeTorsions(); | |
332 | } | |
333 | ||
164 | – | // makeMolecules(); |
334 | ||
335 | // get some of the tricky things that may still be in the globals | |
336 | ||
337 | if( simnfo->n_dipoles ){ | |
338 | ||
339 | if( !the_globals->haveRRF() ){ | |
340 | < | std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n"; |
341 | < | exit(8); |
340 | > | sprintf( painCave.errMsg, |
341 | > | "SimSetup Error, system has dipoles, but no rRF was set.\n"); |
342 | > | painCave.isFatal = 1; |
343 | > | simError(); |
344 | } | |
345 | if( !the_globals->haveDielectric() ){ | |
346 | < | std::cerr << "SimSetup Error, system has dipoles, but no" |
347 | < | << " dielectric was set.\n"; |
348 | < | exit(8); |
346 | > | sprintf( painCave.errMsg, |
347 | > | "SimSetup Error, system has dipoles, but no" |
348 | > | " dielectric was set.\n" ); |
349 | > | painCave.isFatal = 1; |
350 | > | simError(); |
351 | } | |
352 | ||
353 | simnfo->rRF = the_globals->getRRF(); | |
354 | simnfo->dielectric = the_globals->getDielectric(); | |
355 | } | |
356 | ||
357 | + | #ifdef IS_MPI |
358 | + | strcpy( checkPointMsg, "rRf and dielectric check out" ); |
359 | + | MPIcheckPoint(); |
360 | + | #endif // is_mpi |
361 | + | |
362 | if( the_globals->haveBox() ){ | |
363 | simnfo->box_x = the_globals->getBox(); | |
364 | simnfo->box_y = the_globals->getBox(); | |
365 | simnfo->box_z = the_globals->getBox(); | |
366 | } | |
367 | else if( the_globals->haveDensity() ){ | |
368 | < | |
368 | > | |
369 | double vol; | |
370 | vol = (double)tot_nmol / the_globals->getDensity(); | |
371 | simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) ); | |
# | Line 196 | Line 374 | void SimSetup::createSim( void ){ | |
374 | } | |
375 | else{ | |
376 | if( !the_globals->haveBoxX() ){ | |
377 | < | std::cerr << "SimSetup error, no periodic BoxX size given.\n"; |
378 | < | exit(8); |
377 | > | sprintf( painCave.errMsg, |
378 | > | "SimSetup error, no periodic BoxX size given.\n" ); |
379 | > | painCave.isFatal = 1; |
380 | > | simError(); |
381 | } | |
382 | simnfo->box_x = the_globals->getBoxX(); | |
383 | ||
384 | if( !the_globals->haveBoxY() ){ | |
385 | < | std::cerr << "SimSetup error, no periodic BoxY size given.\n"; |
386 | < | exit(8); |
385 | > | sprintf( painCave.errMsg, |
386 | > | "SimSetup error, no periodic BoxY size given.\n" ); |
387 | > | painCave.isFatal = 1; |
388 | > | simError(); |
389 | } | |
390 | simnfo->box_y = the_globals->getBoxY(); | |
391 | ||
392 | if( !the_globals->haveBoxZ() ){ | |
393 | < | std::cerr << "SimSetup error, no periodic BoxZ size given.\n"; |
394 | < | exit(8); |
393 | > | sprintf( painCave.errMsg, |
394 | > | "SimSetup error, no periodic BoxZ size given.\n" ); |
395 | > | painCave.isFatal = 1; |
396 | > | simError(); |
397 | } | |
398 | simnfo->box_z = the_globals->getBoxZ(); | |
399 | } | |
216 | – | |
217 | – | if( the_globals->haveInitialConfig() ){ |
218 | – | InitializeFromFile* fileInit; |
219 | – | fileInit = new InitializeFromFile( the_globals->getInitialConfig() ); |
220 | – | |
221 | – | fileInit->read_xyz( simnfo ); // default velocities on |
400 | ||
401 | < | delete fileInit; |
402 | < | } |
403 | < | else{ |
404 | < | initFromBass(); |
227 | < | } |
401 | > | #ifdef IS_MPI |
402 | > | strcpy( checkPointMsg, "Box size set up" ); |
403 | > | MPIcheckPoint(); |
404 | > | #endif // is_mpi |
405 | ||
406 | < | if( the_globals->haveFinalConfig() ){ |
407 | < | strcpy( simnfo->finalName, the_globals->getFinalConfig() ); |
408 | < | } |
409 | < | else{ |
410 | < | strcpy( simnfo->finalName, inFileName ); |
411 | < | char* endTest; |
412 | < | int nameLength = strlen( simnfo->finalName ); |
413 | < | endTest = &(simnfo->finalName[nameLength - 5]); |
414 | < | if( !strcmp( endTest, ".bass" ) ){ |
415 | < | strcpy( endTest, ".eor" ); |
406 | > | |
407 | > | |
408 | > | if( the_globals->haveInitialConfig() ){ |
409 | > | |
410 | > | InitializeFromFile* fileInit; |
411 | > | #ifdef IS_MPI // is_mpi |
412 | > | if( worldRank == 0 ){ |
413 | > | #endif //is_mpi |
414 | > | fileInit = new InitializeFromFile( the_globals->getInitialConfig() ); |
415 | > | #ifdef IS_MPI |
416 | > | }else fileInit = new InitializeFromFile( NULL ); |
417 | > | #endif |
418 | > | fileInit->read_xyz( simnfo ); // default velocities on |
419 | > | |
420 | > | delete fileInit; |
421 | > | } |
422 | > | else{ |
423 | > | |
424 | > | #ifdef IS_MPI |
425 | > | |
426 | > | // no init from bass |
427 | > | |
428 | > | sprintf( painCave.errMsg, |
429 | > | "Cannot intialize a parallel simulation without an initial configuration file.\n" ); |
430 | > | painCave.isFatal; |
431 | > | simError(); |
432 | > | |
433 | > | #else |
434 | > | |
435 | > | initFromBass(); |
436 | > | |
437 | > | |
438 | > | #endif |
439 | > | } |
440 | > | |
441 | > | #ifdef IS_MPI |
442 | > | strcpy( checkPointMsg, "Successfully read in the initial configuration" ); |
443 | > | MPIcheckPoint(); |
444 | > | #endif // is_mpi |
445 | > | |
446 | > | |
447 | > | |
448 | > | |
449 | > | |
450 | > | |
451 | > | |
452 | > | #ifdef IS_MPI |
453 | > | if( worldRank == 0 ){ |
454 | > | #endif // is_mpi |
455 | > | |
456 | > | if( the_globals->haveFinalConfig() ){ |
457 | > | strcpy( simnfo->finalName, the_globals->getFinalConfig() ); |
458 | } | |
240 | – | else if( !strcmp( endTest, ".BASS" ) ){ |
241 | – | strcpy( endTest, ".eor" ); |
242 | – | } |
459 | else{ | |
460 | < | endTest = &(simnfo->finalName[nameLength - 4]); |
461 | < | if( !strcmp( endTest, ".bss" ) ){ |
460 | > | strcpy( simnfo->finalName, inFileName ); |
461 | > | char* endTest; |
462 | > | int nameLength = strlen( simnfo->finalName ); |
463 | > | endTest = &(simnfo->finalName[nameLength - 5]); |
464 | > | if( !strcmp( endTest, ".bass" ) ){ |
465 | strcpy( endTest, ".eor" ); | |
466 | } | |
467 | < | else if( !strcmp( endTest, ".mdl" ) ){ |
467 | > | else if( !strcmp( endTest, ".BASS" ) ){ |
468 | strcpy( endTest, ".eor" ); | |
469 | } | |
470 | else{ | |
471 | < | strcat( simnfo->finalName, ".eor" ); |
471 | > | endTest = &(simnfo->finalName[nameLength - 4]); |
472 | > | if( !strcmp( endTest, ".bss" ) ){ |
473 | > | strcpy( endTest, ".eor" ); |
474 | > | } |
475 | > | else if( !strcmp( endTest, ".mdl" ) ){ |
476 | > | strcpy( endTest, ".eor" ); |
477 | > | } |
478 | > | else{ |
479 | > | strcat( simnfo->finalName, ".eor" ); |
480 | > | } |
481 | } | |
482 | } | |
255 | – | } |
483 | ||
484 | < | // make the sample and status out names |
485 | < | |
486 | < | strcpy( simnfo->sampleName, inFileName ); |
487 | < | char* endTest; |
488 | < | int nameLength = strlen( simnfo->sampleName ); |
489 | < | endTest = &(simnfo->sampleName[nameLength - 5]); |
490 | < | if( !strcmp( endTest, ".bass" ) ){ |
264 | < | strcpy( endTest, ".dump" ); |
265 | < | } |
266 | < | else if( !strcmp( endTest, ".BASS" ) ){ |
267 | < | strcpy( endTest, ".dump" ); |
268 | < | } |
269 | < | else{ |
270 | < | endTest = &(simnfo->sampleName[nameLength - 4]); |
271 | < | if( !strcmp( endTest, ".bss" ) ){ |
484 | > | // make the sample and status out names |
485 | > | |
486 | > | strcpy( simnfo->sampleName, inFileName ); |
487 | > | char* endTest; |
488 | > | int nameLength = strlen( simnfo->sampleName ); |
489 | > | endTest = &(simnfo->sampleName[nameLength - 5]); |
490 | > | if( !strcmp( endTest, ".bass" ) ){ |
491 | strcpy( endTest, ".dump" ); | |
492 | } | |
493 | < | else if( !strcmp( endTest, ".mdl" ) ){ |
493 | > | else if( !strcmp( endTest, ".BASS" ) ){ |
494 | strcpy( endTest, ".dump" ); | |
495 | } | |
496 | else{ | |
497 | < | strcat( simnfo->sampleName, ".dump" ); |
497 | > | endTest = &(simnfo->sampleName[nameLength - 4]); |
498 | > | if( !strcmp( endTest, ".bss" ) ){ |
499 | > | strcpy( endTest, ".dump" ); |
500 | > | } |
501 | > | else if( !strcmp( endTest, ".mdl" ) ){ |
502 | > | strcpy( endTest, ".dump" ); |
503 | > | } |
504 | > | else{ |
505 | > | strcat( simnfo->sampleName, ".dump" ); |
506 | > | } |
507 | } | |
508 | < | } |
509 | < | |
510 | < | strcpy( simnfo->statusName, inFileName ); |
511 | < | nameLength = strlen( simnfo->statusName ); |
512 | < | endTest = &(simnfo->statusName[nameLength - 5]); |
285 | < | if( !strcmp( endTest, ".bass" ) ){ |
286 | < | strcpy( endTest, ".stat" ); |
287 | < | } |
288 | < | else if( !strcmp( endTest, ".BASS" ) ){ |
289 | < | strcpy( endTest, ".stat" ); |
290 | < | } |
291 | < | else{ |
292 | < | endTest = &(simnfo->statusName[nameLength - 4]); |
293 | < | if( !strcmp( endTest, ".bss" ) ){ |
508 | > | |
509 | > | strcpy( simnfo->statusName, inFileName ); |
510 | > | nameLength = strlen( simnfo->statusName ); |
511 | > | endTest = &(simnfo->statusName[nameLength - 5]); |
512 | > | if( !strcmp( endTest, ".bass" ) ){ |
513 | strcpy( endTest, ".stat" ); | |
514 | } | |
515 | < | else if( !strcmp( endTest, ".mdl" ) ){ |
515 | > | else if( !strcmp( endTest, ".BASS" ) ){ |
516 | strcpy( endTest, ".stat" ); | |
517 | } | |
518 | else{ | |
519 | < | strcat( simnfo->statusName, ".stat" ); |
519 | > | endTest = &(simnfo->statusName[nameLength - 4]); |
520 | > | if( !strcmp( endTest, ".bss" ) ){ |
521 | > | strcpy( endTest, ".stat" ); |
522 | > | } |
523 | > | else if( !strcmp( endTest, ".mdl" ) ){ |
524 | > | strcpy( endTest, ".stat" ); |
525 | > | } |
526 | > | else{ |
527 | > | strcat( simnfo->statusName, ".stat" ); |
528 | > | } |
529 | } | |
530 | + | |
531 | + | #ifdef IS_MPI |
532 | } | |
533 | + | #endif // is_mpi |
534 | ||
535 | // set the status, sample, and themal kick times | |
536 | < | |
536 | > | |
537 | if( the_globals->haveSampleTime() ){ | |
538 | < | simnfo->sampleTime = the_globals->getSampleTime(); |
538 | > | simnfo->sampleTime = the_globals->getSampleTime(); |
539 | simnfo->statusTime = simnfo->sampleTime; | |
540 | simnfo->thermalTime = simnfo->sampleTime; | |
541 | } | |
542 | else{ | |
543 | < | simnfo->sampleTime = the_globals->getRunTime(); |
543 | > | simnfo->sampleTime = the_globals->getRunTime(); |
544 | simnfo->statusTime = simnfo->sampleTime; | |
545 | simnfo->thermalTime = simnfo->sampleTime; | |
546 | } | |
# | Line 325 | Line 556 | void SimSetup::createSim( void ){ | |
556 | // check for the temperature set flag | |
557 | ||
558 | if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet(); | |
559 | < | |
560 | < | |
559 | > | |
560 | > | |
561 | // make the longe range forces and the integrator | |
562 | < | |
562 | > | |
563 | new AllLong( simnfo ); | |
564 | < | |
564 | > | |
565 | if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo ); | |
566 | if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo ); | |
567 | if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo ); | |
568 | } | |
569 | ||
570 | void SimSetup::makeAtoms( void ){ | |
571 | < | |
571 | > | |
572 | int i, j, k, index; | |
573 | double ux, uy, uz, uSqr, u; | |
574 | AtomStamp* current_atom; | |
575 | DirectionalAtom* dAtom; | |
576 | + | int molIndex, molStart, molEnd, nMemb, lMolIndex; |
577 | ||
578 | + | lMolIndex = 0; |
579 | + | molIndex = 0; |
580 | index = 0; | |
581 | for( i=0; i<n_components; i++ ){ | |
582 | < | |
582 | > | |
583 | for( j=0; j<components_nmol[i]; j++ ){ | |
350 | – | |
351 | – | for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){ |
352 | – | |
353 | – | current_atom = comp_stamps[i]->getAtom( k ); |
354 | – | if( current_atom->haveOrientation() ){ |
584 | ||
585 | < | dAtom = new DirectionalAtom; |
586 | < | simnfo->n_oriented++; |
587 | < | the_atoms[index] = dAtom; |
588 | < | |
589 | < | ux = current_atom->getOrntX(); |
590 | < | uy = current_atom->getOrntY(); |
591 | < | uz = current_atom->getOrntZ(); |
585 | > | #ifdef IS_MPI |
586 | > | if( simnfo->mpiSim->getMyMolStart() <= molIndex && |
587 | > | molIndex <= simnfo->mpiSim->getMyMolEnd() ){ |
588 | > | #endif // is_mpi |
589 | > | |
590 | > | molStart = index; |
591 | > | nMemb = comp_stamps[i]->getNAtoms(); |
592 | > | for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){ |
593 | ||
594 | < | uSqr = (ux * ux) + (uy * uy) + (uz * uz); |
594 | > | current_atom = comp_stamps[i]->getAtom( k ); |
595 | > | if( current_atom->haveOrientation() ){ |
596 | > | |
597 | > | dAtom = new DirectionalAtom(index); |
598 | > | simnfo->n_oriented++; |
599 | > | the_atoms[index] = dAtom; |
600 | > | |
601 | > | ux = current_atom->getOrntX(); |
602 | > | uy = current_atom->getOrntY(); |
603 | > | uz = current_atom->getOrntZ(); |
604 | > | |
605 | > | uSqr = (ux * ux) + (uy * uy) + (uz * uz); |
606 | > | |
607 | > | u = sqrt( uSqr ); |
608 | > | ux = ux / u; |
609 | > | uy = uy / u; |
610 | > | uz = uz / u; |
611 | > | |
612 | > | dAtom->setSUx( ux ); |
613 | > | dAtom->setSUy( uy ); |
614 | > | dAtom->setSUz( uz ); |
615 | > | } |
616 | > | else{ |
617 | > | the_atoms[index] = new GeneralAtom(index); |
618 | > | } |
619 | > | the_atoms[index]->setType( current_atom->getType() ); |
620 | > | the_atoms[index]->setIndex( index ); |
621 | ||
622 | < | u = sqrt( uSqr ); |
623 | < | ux = ux / u; |
368 | < | uy = uy / u; |
369 | < | uz = uz / u; |
370 | < | |
371 | < | dAtom->setSUx( ux ); |
372 | < | dAtom->setSUy( uy ); |
373 | < | dAtom->setSUz( uz ); |
622 | > | // increment the index and repeat; |
623 | > | index++; |
624 | } | |
375 | – | else{ |
376 | – | the_atoms[index] = new GeneralAtom; |
377 | – | } |
378 | – | the_atoms[index]->setType( current_atom->getType() ); |
379 | – | the_atoms[index]->setIndex( index ); |
625 | ||
626 | < | // increment the index and repeat; |
627 | < | index++; |
626 | > | molEnd = index -1; |
627 | > | the_molecules[lMolIndex].setNMembers( nMemb ); |
628 | > | the_molecules[lMolIndex].setStartAtom( molStart ); |
629 | > | the_molecules[lMolIndex].setEndAtom( molEnd ); |
630 | > | the_molecules[lMolIndex].setStampID( i ); |
631 | > | lMolIndex++; |
632 | > | |
633 | > | #ifdef IS_MPI |
634 | } | |
635 | + | #endif //is_mpi |
636 | + | |
637 | + | molIndex++; |
638 | } | |
639 | } | |
640 | < | |
640 | > | |
641 | the_ff->initializeAtoms(); | |
642 | } | |
643 | ||
644 | void SimSetup::makeBonds( void ){ | |
645 | ||
646 | < | int i, j, k, index, offset; |
646 | > | int i, j, k, index, offset, molIndex; |
647 | bond_pair* the_bonds; | |
648 | BondStamp* current_bond; | |
649 | ||
650 | the_bonds = new bond_pair[tot_bonds]; | |
651 | index = 0; | |
652 | offset = 0; | |
653 | + | molIndex = 0; |
654 | for( i=0; i<n_components; i++ ){ | |
655 | < | |
655 | > | |
656 | for( j=0; j<components_nmol[i]; j++ ){ | |
402 | – | |
403 | – | for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){ |
404 | – | |
405 | – | current_bond = comp_stamps[i]->getBond( k ); |
406 | – | the_bonds[index].a = current_bond->getA() + offset; |
407 | – | the_bonds[index].b = current_bond->getB() + offset; |
657 | ||
658 | < | the_excludes[index].i = the_bonds[index].a; |
659 | < | the_excludes[index].j = the_bonds[index].b; |
660 | < | |
661 | < | // increment the index and repeat; |
662 | < | index++; |
663 | < | } |
664 | < | offset += comp_stamps[i]->getNAtoms(); |
665 | < | } |
658 | > | #ifdef IS_MPI |
659 | > | if( simnfo->mpiSim->getMyMolStart() <= molIndex && |
660 | > | molIndex <= simnfo->mpiSim->getMyMolEnd() ){ |
661 | > | #endif // is_mpi |
662 | > | |
663 | > | for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){ |
664 | > | |
665 | > | current_bond = comp_stamps[i]->getBond( k ); |
666 | > | the_bonds[index].a = current_bond->getA() + offset; |
667 | > | the_bonds[index].b = current_bond->getB() + offset; |
668 | > | |
669 | > | the_excludes[index].i = the_bonds[index].a; |
670 | > | the_excludes[index].j = the_bonds[index].b; |
671 | > | |
672 | > | // increment the index and repeat; |
673 | > | index++; |
674 | > | } |
675 | > | offset += comp_stamps[i]->getNAtoms(); |
676 | > | |
677 | > | #ifdef IS_MPI |
678 | > | } |
679 | > | #endif is_mpi |
680 | > | |
681 | > | molIndex++; |
682 | > | } |
683 | } | |
684 | < | |
684 | > | |
685 | the_ff->initializeBonds( the_bonds ); | |
686 | } | |
687 | ||
688 | void SimSetup::makeBends( void ){ | |
689 | ||
690 | < | int i, j, k, index, offset; |
690 | > | int i, j, k, index, offset, molIndex; |
691 | bend_set* the_bends; | |
692 | BendStamp* current_bend; | |
693 | ||
694 | the_bends = new bend_set[tot_bends]; | |
695 | index = 0; | |
696 | offset = 0; | |
697 | + | molIndex = 0; |
698 | for( i=0; i<n_components; i++ ){ | |
699 | < | |
699 | > | |
700 | for( j=0; j<components_nmol[i]; j++ ){ | |
434 | – | |
435 | – | for( k=0; k<comp_stamps[i]->getNBends(); k++ ){ |
436 | – | |
437 | – | current_bend = comp_stamps[i]->getBend( k ); |
438 | – | the_bends[index].a = current_bend->getA() + offset; |
439 | – | the_bends[index].b = current_bend->getB() + offset; |
440 | – | the_bends[index].c = current_bend->getC() + offset; |
701 | ||
702 | < | the_excludes[index + tot_bonds].i = the_bends[index].a; |
703 | < | the_excludes[index + tot_bonds].j = the_bends[index].c; |
702 | > | #ifdef IS_MPI |
703 | > | if( simnfo->mpiSim->getMyMolStart() <= molIndex && |
704 | > | molIndex <= simnfo->mpiSim->getMyMolEnd() ){ |
705 | > | #endif // is_mpi |
706 | ||
707 | < | // increment the index and repeat; |
708 | < | index++; |
707 | > | for( k=0; k<comp_stamps[i]->getNBends(); k++ ){ |
708 | > | |
709 | > | current_bend = comp_stamps[i]->getBend( k ); |
710 | > | the_bends[index].a = current_bend->getA() + offset; |
711 | > | the_bends[index].b = current_bend->getB() + offset; |
712 | > | the_bends[index].c = current_bend->getC() + offset; |
713 | > | |
714 | > | the_excludes[index + tot_bonds].i = the_bends[index].a; |
715 | > | the_excludes[index + tot_bonds].j = the_bends[index].c; |
716 | > | |
717 | > | // increment the index and repeat; |
718 | > | index++; |
719 | > | } |
720 | > | offset += comp_stamps[i]->getNAtoms(); |
721 | > | |
722 | > | #ifdef IS_MPI |
723 | } | |
724 | < | offset += comp_stamps[i]->getNAtoms(); |
724 | > | #endif //is_mpi |
725 | > | |
726 | > | molIndex++; |
727 | } | |
728 | } | |
729 | < | |
729 | > | |
730 | the_ff->initializeBends( the_bends ); | |
731 | } | |
732 | ||
733 | void SimSetup::makeTorsions( void ){ | |
734 | ||
735 | < | int i, j, k, index, offset; |
735 | > | int i, j, k, index, offset, molIndex; |
736 | torsion_set* the_torsions; | |
737 | TorsionStamp* current_torsion; | |
738 | ||
739 | the_torsions = new torsion_set[tot_torsions]; | |
740 | index = 0; | |
741 | offset = 0; | |
742 | + | molIndex = 0; |
743 | for( i=0; i<n_components; i++ ){ | |
744 | < | |
744 | > | |
745 | for( j=0; j<components_nmol[i]; j++ ){ | |
746 | < | |
746 | > | |
747 | > | #ifdef IS_MPI |
748 | > | if( simnfo->mpiSim->getMyMolStart() <= molIndex && |
749 | > | molIndex <= simnfo->mpiSim->getMyMolEnd() ){ |
750 | > | #endif // is_mpi |
751 | > | |
752 | for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){ | |
753 | < | |
753 | > | |
754 | current_torsion = comp_stamps[i]->getTorsion( k ); | |
755 | the_torsions[index].a = current_torsion->getA() + offset; | |
756 | the_torsions[index].b = current_torsion->getB() + offset; | |
# | Line 480 | Line 764 | void SimSetup::makeTorsions( void ){ | |
764 | index++; | |
765 | } | |
766 | offset += comp_stamps[i]->getNAtoms(); | |
767 | + | |
768 | + | #ifdef IS_MPI |
769 | + | } |
770 | + | #endif //is_mpi |
771 | + | |
772 | + | molIndex++; |
773 | } | |
774 | } | |
775 | < | |
775 | > | |
776 | the_ff->initializeTorsions( the_torsions ); | |
777 | } | |
778 | ||
489 | – | void SimSetup::makeMolecules( void ){ |
490 | – | |
491 | – | //empy for now |
492 | – | } |
493 | – | |
779 | void SimSetup::initFromBass( void ){ | |
780 | ||
781 | int i, j, k; | |
# | Line 518 | Line 803 | void SimSetup::initFromBass( void ){ | |
803 | n_per_extra = (int)ceil( temp1 ); | |
804 | ||
805 | if( n_per_extra > 4){ | |
806 | < | std::cerr << "THere has been an error in constructing the non-complete lattice.\n"; |
807 | < | exit(8); |
806 | > | sprintf( painCave.errMsg, |
807 | > | "SimSetup error. There has been an error in constructing" |
808 | > | " the non-complete lattice.\n" ); |
809 | > | painCave.isFatal = 1; |
810 | > | simError(); |
811 | } | |
812 | } | |
813 | else{ | |
# | Line 528 | Line 816 | void SimSetup::initFromBass( void ){ | |
816 | celly = simnfo->box_y / temp3; | |
817 | cellz = simnfo->box_z / temp3; | |
818 | } | |
819 | < | |
819 | > | |
820 | current_mol = 0; | |
821 | current_comp_mol = 0; | |
822 | current_comp = 0; | |
823 | current_atom_ndx = 0; | |
824 | < | |
824 | > | |
825 | for( i=0; i < n_cells ; i++ ){ | |
826 | for( j=0; j < n_cells; j++ ){ | |
827 | for( k=0; k < n_cells; k++ ){ | |
828 | < | |
828 | > | |
829 | makeElement( i * cellx, | |
830 | j * celly, | |
831 | k * cellz ); | |
832 | < | |
832 | > | |
833 | makeElement( i * cellx + 0.5 * cellx, | |
834 | j * celly + 0.5 * celly, | |
835 | k * cellz ); | |
836 | < | |
836 | > | |
837 | makeElement( i * cellx, | |
838 | j * celly + 0.5 * celly, | |
839 | k * cellz + 0.5 * cellz ); | |
840 | < | |
840 | > | |
841 | makeElement( i * cellx + 0.5 * cellx, | |
842 | j * celly, | |
843 | k * cellz + 0.5 * cellz ); | |
# | Line 559 | Line 847 | void SimSetup::initFromBass( void ){ | |
847 | ||
848 | if( have_extra ){ | |
849 | done = 0; | |
850 | < | |
850 | > | |
851 | int start_ndx; | |
852 | for( i=0; i < (n_cells+1) && !done; i++ ){ | |
853 | for( j=0; j < (n_cells+1) && !done; j++ ){ | |
854 | < | |
854 | > | |
855 | if( i < n_cells ){ | |
856 | < | |
856 | > | |
857 | if( j < n_cells ){ | |
858 | start_ndx = n_cells; | |
859 | } | |
860 | else start_ndx = 0; | |
861 | } | |
862 | else start_ndx = 0; | |
863 | < | |
863 | > | |
864 | for( k=start_ndx; k < (n_cells+1) && !done; k++ ){ | |
865 | < | |
865 | > | |
866 | makeElement( i * cellx, | |
867 | j * celly, | |
868 | k * cellz ); | |
869 | done = ( current_mol >= tot_nmol ); | |
870 | < | |
870 | > | |
871 | if( !done && n_per_extra > 1 ){ | |
872 | makeElement( i * cellx + 0.5 * cellx, | |
873 | j * celly + 0.5 * celly, | |
874 | k * cellz ); | |
875 | done = ( current_mol >= tot_nmol ); | |
876 | } | |
877 | < | |
877 | > | |
878 | if( !done && n_per_extra > 2){ | |
879 | makeElement( i * cellx, | |
880 | j * celly + 0.5 * celly, | |
881 | k * cellz + 0.5 * cellz ); | |
882 | done = ( current_mol >= tot_nmol ); | |
883 | } | |
884 | < | |
884 | > | |
885 | if( !done && n_per_extra > 3){ | |
886 | makeElement( i * cellx + 0.5 * cellx, | |
887 | j * celly, | |
# | Line 604 | Line 892 | void SimSetup::initFromBass( void ){ | |
892 | } | |
893 | } | |
894 | } | |
895 | < | |
896 | < | |
895 | > | |
896 | > | |
897 | for( i=0; i<simnfo->n_atoms; i++ ){ | |
898 | simnfo->atoms[i]->set_vx( 0.0 ); | |
899 | simnfo->atoms[i]->set_vy( 0.0 ); | |
# | Line 621 | Line 909 | void SimSetup::makeElement( double x, double y, double | |
909 | double rotMat[3][3]; | |
910 | ||
911 | for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){ | |
912 | < | |
912 | > | |
913 | current_atom = comp_stamps[current_comp]->getAtom( k ); | |
914 | if( !current_atom->havePosition() ){ | |
915 | < | std::cerr << "Component " << comp_stamps[current_comp]->getID() |
916 | < | << ", atom " << current_atom->getType() |
917 | < | << " does not have a position specified.\n" |
918 | < | << "The initialization routine is unable to give a start" |
919 | < | << " position.\n"; |
920 | < | exit(8); |
915 | > | sprintf( painCave.errMsg, |
916 | > | "SimSetup:initFromBass error.\n" |
917 | > | "\tComponent %s, atom %s does not have a position specified.\n" |
918 | > | "\tThe initialization routine is unable to give a start" |
919 | > | " position.\n", |
920 | > | comp_stamps[current_comp]->getID(), |
921 | > | current_atom->getType() ); |
922 | > | painCave.isFatal = 1; |
923 | > | simError(); |
924 | } | |
925 | < | |
925 | > | |
926 | the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() ); | |
927 | the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() ); | |
928 | the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() ); | |
929 | < | |
929 | > | |
930 | if( the_atoms[current_atom_ndx]->isDirectional() ){ | |
931 | < | |
931 | > | |
932 | dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx]; | |
933 | < | |
933 | > | |
934 | rotMat[0][0] = 1.0; | |
935 | rotMat[0][1] = 0.0; | |
936 | rotMat[0][2] = 0.0; | |
# | Line 657 | Line 948 | void SimSetup::makeElement( double x, double y, double | |
948 | ||
949 | current_atom_ndx++; | |
950 | } | |
951 | < | |
951 | > | |
952 | current_mol++; | |
953 | current_comp_mol++; | |
954 | ||
955 | if( current_comp_mol >= components_nmol[current_comp] ){ | |
956 | < | |
956 | > | |
957 | current_comp_mol = 0; | |
958 | current_comp++; | |
959 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |