ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/interface_implementation/SimSetup.cpp
(Generate patch)

Comparing trunk/mdtools/interface_implementation/SimSetup.cpp (file contents):
Revision 113 by mmeineke, Mon Sep 23 15:12:56 2002 UTC vs.
Revision 253 by chuckv, Thu Jan 30 15:20:21 2003 UTC

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

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines