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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 614 by mmeineke, Tue Jul 15 17:57:04 2003 UTC vs.
Revision 656 by mmeineke, Tue Jul 29 16:32:37 2003 UTC

# Line 24 | Line 24
24  
25   #define FF_DUFF 0
26   #define FF_LJ   1
27 + #define FF_EAM  2
28  
28 –
29   SimSetup::SimSetup(){
30 +  
31 +  isInfoArray = 0;
32 +  nInfo = 1;
33 +  
34    stamps = new MakeStamps();
35    globals = new Globals();
36    
37 +  
38   #ifdef IS_MPI
39    strcpy( checkPointMsg, "SimSetup creation successful" );
40    MPIcheckPoint();
# Line 41 | Line 46 | SimSetup::~SimSetup(){
46    delete globals;
47   }
48  
49 + void SimSetup::setSimInfo( SimInfo* the_info, int theNinfo ) {
50 +    info = the_info;
51 +    nInfo = theNinfo;
52 +    isInfoArray = 1;
53 +  }
54 +
55 +
56   void SimSetup::parseFile( char* fileName ){
57  
58   #ifdef IS_MPI
# Line 78 | Line 90 | void SimSetup::createSim( void ){
90  
91   void SimSetup::createSim( void ){
92  
81 –  MakeStamps *the_stamps;
82 –  Globals* the_globals;
93    int i, j, k, globalAtomIndex;
94    
95    // gather all of the information from the Bass file
# Line 90 | Line 100 | void SimSetup::createSim( void ){
100  
101    sysObjectsCreation();
102  
93 –
94 –
95 –  // initialize the arrays
96 –
97 –
98 –
99 –  makeMolecules();
100 –  info->identArray = new int[info->n_atoms];
101 –  for(i=0; i<info->n_atoms; i++){
102 –    info->identArray[i] = the_atoms[i]->getIdent();
103 –  }
104 –  
105 –
103    // check on the post processing info
104    
105    finalInfoCheck();
106  
110 –
111 –
112 –
107    // initialize the system coordinates
108  
109    initSystemCoords();
# Line 119 | Line 113 | void SimSetup::createSim( void ){
113  
114    makeOutNames();
115    
122 –
123 –  
124 –
125 –  
126 –
127 –  
128 –
129 –
116    // make the integrator
117    
118 +  makeIntegrator();
119    
133 –  NVT*  myNVT = NULL;
134 –  NPTi* myNPTi = NULL;
135 –  NPTf* myNPTf = NULL;
136 –  NPTim* myNPTim = NULL;
137 –  NPTfm* myNPTfm = NULL;
138 –
139 –  switch( ensembleCase ){
140 –
141 –  case NVE_ENS:
142 –    new NVE( info, the_ff );
143 –    break;
144 –
145 –  case NVT_ENS:
146 –    myNVT = new NVT( info, the_ff );
147 –    myNVT->setTargetTemp(the_globals->getTargetTemp());
148 –
149 –    if (the_globals->haveTauThermostat())
150 –      myNVT->setTauThermostat(the_globals->getTauThermostat());
151 –
152 –    else {
153 –      sprintf( painCave.errMsg,
154 –               "SimSetup error: If you use the NVT\n"
155 –               "    ensemble, you must set tauThermostat.\n");
156 –      painCave.isFatal = 1;
157 –      simError();
158 –    }
159 –    break;
160 –
161 –  case NPTi_ENS:
162 –    myNPTi = new NPTi( info, the_ff );
163 –    myNPTi->setTargetTemp( the_globals->getTargetTemp() );
164 –
165 –    if (the_globals->haveTargetPressure())
166 –      myNPTi->setTargetPressure(the_globals->getTargetPressure());
167 –    else {
168 –      sprintf( painCave.errMsg,
169 –               "SimSetup error: If you use a constant pressure\n"
170 –               "    ensemble, you must set targetPressure in the BASS file.\n");
171 –      painCave.isFatal = 1;
172 –      simError();
173 –    }
174 –    
175 –    if( the_globals->haveTauThermostat() )
176 –      myNPTi->setTauThermostat( the_globals->getTauThermostat() );
177 –    else{
178 –      sprintf( painCave.errMsg,
179 –               "SimSetup error: If you use an NPT\n"
180 –               "    ensemble, you must set tauThermostat.\n");
181 –      painCave.isFatal = 1;
182 –      simError();
183 –    }
184 –
185 –    if( the_globals->haveTauBarostat() )
186 –      myNPTi->setTauBarostat( the_globals->getTauBarostat() );
187 –    else{
188 –      sprintf( painCave.errMsg,
189 –               "SimSetup error: If you use an NPT\n"
190 –               "    ensemble, you must set tauBarostat.\n");
191 –      painCave.isFatal = 1;
192 –      simError();
193 –    }
194 –    break;
195 –
196 –  case NPTf_ENS:
197 –    myNPTf = new NPTf( info, the_ff );
198 –    myNPTf->setTargetTemp( the_globals->getTargetTemp());
199 –
200 –    if (the_globals->haveTargetPressure())
201 –      myNPTf->setTargetPressure(the_globals->getTargetPressure());
202 –    else {
203 –      sprintf( painCave.errMsg,
204 –               "SimSetup error: If you use a constant pressure\n"
205 –               "    ensemble, you must set targetPressure in the BASS file.\n");
206 –      painCave.isFatal = 1;
207 –      simError();
208 –    }    
209 –
210 –    if( the_globals->haveTauThermostat() )
211 –      myNPTf->setTauThermostat( the_globals->getTauThermostat() );
212 –    else{
213 –      sprintf( painCave.errMsg,
214 –               "SimSetup error: If you use an NPT\n"
215 –               "    ensemble, you must set tauThermostat.\n");
216 –      painCave.isFatal = 1;
217 –      simError();
218 –    }
219 –
220 –    if( the_globals->haveTauBarostat() )
221 –      myNPTf->setTauBarostat( the_globals->getTauBarostat() );
222 –    else{
223 –      sprintf( painCave.errMsg,
224 –               "SimSetup error: If you use an NPT\n"
225 –               "    ensemble, you must set tauBarostat.\n");
226 –      painCave.isFatal = 1;
227 –      simError();
228 –    }
229 –    break;
230 –    
231 –  case NPTim_ENS:
232 –    myNPTim = new NPTim( info, the_ff );
233 –    myNPTim->setTargetTemp( the_globals->getTargetTemp());
234 –
235 –    if (the_globals->haveTargetPressure())
236 –      myNPTim->setTargetPressure(the_globals->getTargetPressure());
237 –    else {
238 –      sprintf( painCave.errMsg,
239 –               "SimSetup error: If you use a constant pressure\n"
240 –               "    ensemble, you must set targetPressure in the BASS file.\n");
241 –      painCave.isFatal = 1;
242 –      simError();
243 –    }
244 –    
245 –    if( the_globals->haveTauThermostat() )
246 –      myNPTim->setTauThermostat( the_globals->getTauThermostat() );
247 –    else{
248 –      sprintf( painCave.errMsg,
249 –               "SimSetup error: If you use an NPT\n"
250 –               "    ensemble, you must set tauThermostat.\n");
251 –      painCave.isFatal = 1;
252 –      simError();
253 –    }
254 –
255 –    if( the_globals->haveTauBarostat() )
256 –      myNPTim->setTauBarostat( the_globals->getTauBarostat() );
257 –    else{
258 –      sprintf( painCave.errMsg,
259 –               "SimSetup error: If you use an NPT\n"
260 –               "    ensemble, you must set tauBarostat.\n");
261 –      painCave.isFatal = 1;
262 –      simError();
263 –    }
264 –    break;
265 –
266 –  case NPTfm_ENS:
267 –    myNPTfm = new NPTfm( info, the_ff );
268 –    myNPTfm->setTargetTemp( the_globals->getTargetTemp());
269 –
270 –    if (the_globals->haveTargetPressure())
271 –      myNPTfm->setTargetPressure(the_globals->getTargetPressure());
272 –    else {
273 –      sprintf( painCave.errMsg,
274 –               "SimSetup error: If you use a constant pressure\n"
275 –               "    ensemble, you must set targetPressure in the BASS file.\n");
276 –      painCave.isFatal = 1;
277 –      simError();
278 –    }
279 –    
280 –    if( the_globals->haveTauThermostat() )
281 –      myNPTfm->setTauThermostat( the_globals->getTauThermostat() );
282 –    else{
283 –      sprintf( painCave.errMsg,
284 –               "SimSetup error: If you use an NPT\n"
285 –               "    ensemble, you must set tauThermostat.\n");
286 –      painCave.isFatal = 1;
287 –      simError();
288 –    }
289 –
290 –    if( the_globals->haveTauBarostat() )
291 –      myNPTfm->setTauBarostat( the_globals->getTauBarostat() );
292 –    else{
293 –      sprintf( painCave.errMsg,
294 –               "SimSetup error: If you use an NPT\n"
295 –               "    ensemble, you must set tauBarostat.\n");
296 –      painCave.isFatal = 1;
297 –      simError();
298 –    }
299 –    break;
300 –
301 –  default:
302 –    sprintf( painCave.errMsg,
303 –             "SimSetup Error. Unrecognized ensemble in case statement.\n");
304 –    painCave.isFatal = 1;
305 –    simError();
306 –  }
307 –
308 –
120   #ifdef IS_MPI
121    mpiSim->mpiRefresh();
122   #endif
123  
124    // initialize the Fortran
125  
126 +  initFortran();
127  
316 –  info->refreshSim();
317 –  
318 –  if( !strcmp( info->mixingRule, "standard") ){
319 –    the_ff->initForceField( LB_MIXING_RULE );
320 –  }
321 –  else if( !strcmp( info->mixingRule, "explicit") ){
322 –    the_ff->initForceField( EXPLICIT_MIXING_RULE );
323 –  }
324 –  else{
325 –    sprintf( painCave.errMsg,
326 –             "SimSetup Error: unknown mixing rule -> \"%s\"\n",
327 –             info->mixingRule );
328 –    painCave.isFatal = 1;
329 –    simError();
330 –  }
128  
129  
333 – #ifdef IS_MPI
334 –  strcpy( checkPointMsg,
335 –          "Successfully intialized the mixingRule for Fortran." );
336 –  MPIcheckPoint();
337 – #endif // is_mpi
130   }
131  
132  
133   void SimSetup::makeMolecules( void ){
134  
135    int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
136 <  molInit info;
136 >  molInit molInfo;
137    DirectionalAtom* dAtom;
138    LinkedAssign* extras;
139    LinkedAssign* current_extra;
# Line 370 | Line 162 | void SimSetup::makeMolecules( void ){
162      
163      stampID = the_molecules[i].getStampID();
164  
165 <    info.nAtoms    = comp_stamps[stampID]->getNAtoms();
166 <    info.nBonds    = comp_stamps[stampID]->getNBonds();
167 <    info.nBends    = comp_stamps[stampID]->getNBends();
168 <    info.nTorsions = comp_stamps[stampID]->getNTorsions();
169 <    info.nExcludes = info.nBonds + info.nBends + info.nTorsions;
165 >    molInfo.nAtoms    = comp_stamps[stampID]->getNAtoms();
166 >    molInfo.nBonds    = comp_stamps[stampID]->getNBonds();
167 >    molInfo.nBends    = comp_stamps[stampID]->getNBends();
168 >    molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
169 >    molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
170  
171 <    info.myAtoms = &the_atoms[atomOffset];
172 <    info.myExcludes = &the_excludes[excludeOffset];
173 <    info.myBonds = new Bond*[info.nBonds];
174 <    info.myBends = new Bend*[info.nBends];
175 <    info.myTorsions = new Torsion*[info.nTorsions];
171 >    molInfo.myAtoms = &the_atoms[atomOffset];
172 >    molInfo.myExcludes = &the_excludes[excludeOffset];
173 >    molInfo.myBonds = new Bond*[molInfo.nBonds];
174 >    molInfo.myBends = new Bend*[molInfo.nBends];
175 >    molInfo.myTorsions = new Torsion*[molInfo.nTorsions];
176  
177 <    theBonds = new bond_pair[info.nBonds];
178 <    theBends = new bend_set[info.nBends];
179 <    theTorsions = new torsion_set[info.nTorsions];
177 >    theBonds = new bond_pair[molInfo.nBonds];
178 >    theBends = new bend_set[molInfo.nBends];
179 >    theTorsions = new torsion_set[molInfo.nTorsions];
180      
181      // make the Atoms
182      
183 <    for(j=0; j<info.nAtoms; j++){
183 >    for(j=0; j<molInfo.nAtoms; j++){
184        
185        currentAtom = comp_stamps[stampID]->getAtom( j );
186        if( currentAtom->haveOrientation() ){
187          
188          dAtom = new DirectionalAtom(j + atomOffset);
189          info->n_oriented++;
190 <        info.myAtoms[j] = dAtom;
190 >        molInfo.myAtoms[j] = dAtom;
191          
192          ux = currentAtom->getOrntX();
193          uy = currentAtom->getOrntY();
# Line 413 | Line 205 | void SimSetup::makeMolecules( void ){
205          dAtom->setSUz( uz );
206        }
207        else{
208 <        info.myAtoms[j] = new GeneralAtom(j + atomOffset);
208 >        molInfo.myAtoms[j] = new GeneralAtom(j + atomOffset);
209        }
210 <      info.myAtoms[j]->setType( currentAtom->getType() );
210 >      molInfo.myAtoms[j]->setType( currentAtom->getType() );
211      
212   #ifdef IS_MPI
213        
214 <      info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
214 >      molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
215        
216   #endif // is_mpi
217      }
218      
219      // make the bonds
220 <    for(j=0; j<info.nBonds; j++){
220 >    for(j=0; j<molInfo.nBonds; j++){
221        
222        currentBond = comp_stamps[stampID]->getBond( j );
223        theBonds[j].a = currentBond->getA() + atomOffset;
# Line 452 | Line 244 | void SimSetup::makeMolecules( void ){
244        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
245   #endif  //is_mpi
246      }
247 <    excludeOffset += info.nBonds;
247 >    excludeOffset += molInfo.nBonds;
248  
249      //make the bends
250 <    for(j=0; j<info.nBends; j++){
250 >    for(j=0; j<molInfo.nBends; j++){
251        
252        currentBend = comp_stamps[stampID]->getBend( j );
253        theBends[j].a = currentBend->getA() + atomOffset;
# Line 538 | Line 330 | void SimSetup::makeMolecules( void ){
330        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
331   #endif  //is_mpi
332      }
333 <    excludeOffset += info.nBends;
333 >    excludeOffset += molInfo.nBends;
334  
335 <    for(j=0; j<info.nTorsions; j++){
335 >    for(j=0; j<molInfo.nTorsions; j++){
336        
337        currentTorsion = comp_stamps[stampID]->getTorsion( j );
338        theTorsions[j].a = currentTorsion->getA() + atomOffset;
# Line 568 | Line 360 | void SimSetup::makeMolecules( void ){
360        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
361   #endif  //is_mpi
362      }
363 <    excludeOffset += info.nTorsions;
363 >    excludeOffset += molInfo.nTorsions;
364  
365      
366      // send the arrays off to the forceField for init.
367  
368 <    the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
369 <    the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
370 <    the_ff->initializeBends( info.nBends, info.myBends, theBends );
371 <    the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
368 >    the_ff->initializeAtoms( molInfo.nAtoms, molInfo.myAtoms );
369 >    the_ff->initializeBonds( molInfo.nBonds, molInfo.myBonds, theBonds );
370 >    the_ff->initializeBends( molInfo.nBends, molInfo.myBends, theBends );
371 >    the_ff->initializeTorsions( molInfo.nTorsions, molInfo.myTorsions, theTorsions );
372  
373  
374 <    the_molecules[i].initialize( info );
374 >    the_molecules[i].initialize( molInfo );
375  
376  
377 <    atomOffset += info.nAtoms;
377 >    atomOffset += molInfo.nAtoms;
378      delete[] theBonds;
379      delete[] theBends;
380      delete[] theTorsions;
# Line 618 | Line 410 | void SimSetup::initFromBass( void ){
410      have_extra =1;
411  
412      n_cells = (int)temp3 - 1;
413 <    cellx = info->boxLx / temp3;
414 <    celly = info->boxLy / temp3;
415 <    cellz = info->boxLz / temp3;
413 >    cellx = info->boxL[0] / temp3;
414 >    celly = info->boxL[1] / temp3;
415 >    cellz = info->boxL[2] / temp3;
416      n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
417      temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
418      n_per_extra = (int)ceil( temp1 );
# Line 635 | Line 427 | void SimSetup::initFromBass( void ){
427    }
428    else{
429      n_cells = (int)temp3;
430 <    cellx = info->boxLx / temp3;
431 <    celly = info->boxLy / temp3;
432 <    cellz = info->boxLz / temp3;
430 >    cellx = info->boxL[0] / temp3;
431 >    celly = info->boxL[1] / temp3;
432 >    cellz = info->boxL[2] / temp3;
433    }
434  
435    current_mol = 0;
# Line 784 | Line 576 | void SimSetup::gatherInfo( void ){
576  
577  
578   void SimSetup::gatherInfo( void ){
579 +  int i,j,k;
580  
581    ensembleCase = -1;
582    ffCase = -1;
583  
584    // get the stamps and globals;
585 <  the_stamps = stamps;
586 <  the_globals = globals;
585 >  stamps = stamps;
586 >  globals = globals;
587  
588    // set the easy ones first
589 <  info->target_temp = the_globals->getTargetTemp();
590 <  info->dt = the_globals->getDt();
591 <  info->run_time = the_globals->getRunTime();
592 <  n_components = the_globals->getNComponents();
589 >  info->target_temp = globals->getTargetTemp();
590 >  info->dt = globals->getDt();
591 >  info->run_time = globals->getRunTime();
592 >  n_components = globals->getNComponents();
593  
594  
595    // get the forceField
596  
597 <  strcpy( force_field, the_globals->getForceField() );
597 >  strcpy( force_field, globals->getForceField() );
598  
599    if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
600    else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
601 +  else if( !strcasecmp( force_field, "EAM" )) ffCase = FF_EAM;
602    else{
603      sprintf( painCave.errMsg,
604               "SimSetup Error. Unrecognized force field -> %s\n",
# Line 815 | Line 609 | void SimSetup::gatherInfo( void ){
609  
610    // get the ensemble
611  
612 <  strcpy( ensemble, the_globals->getEnsemble() );
612 >  strcpy( ensemble, globals->getEnsemble() );
613  
614    if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS;
615    else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
# Line 838 | Line 632 | void SimSetup::gatherInfo( void ){
632  
633    // get the mixing rule
634  
635 <  strcpy( info->mixingRule, the_globals->getMixingRule() );
636 <  info->usePBC = the_globals->getPBC();
635 >  strcpy( info->mixingRule, globals->getMixingRule() );
636 >  info->usePBC = globals->getPBC();
637          
638    
639    // get the components and calculate the tot_nMol and indvidual n_mol
640  
641 <  the_components = the_globals->getComponents();
641 >  the_components = globals->getComponents();
642    components_nmol = new int[n_components];
643  
644  
645 <  if( !the_globals->haveNMol() ){
645 >  if( !globals->haveNMol() ){
646      // we don't have the total number of molecules, so we assume it is
647      // given in each component
648  
# Line 881 | Line 675 | void SimSetup::gatherInfo( void ){
675  
676    // set the status, sample, and thermal kick times
677    
678 <  if( the_globals->haveSampleTime() ){
679 <    info->sampleTime = the_globals->getSampleTime();
678 >  if( globals->haveSampleTime() ){
679 >    info->sampleTime = globals->getSampleTime();
680      info->statusTime = info->sampleTime;
681      info->thermalTime = info->sampleTime;
682    }
683    else{
684 <    info->sampleTime = the_globals->getRunTime();
684 >    info->sampleTime = globals->getRunTime();
685      info->statusTime = info->sampleTime;
686      info->thermalTime = info->sampleTime;
687    }
688  
689 <  if( the_globals->haveStatusTime() ){
690 <    info->statusTime = the_globals->getStatusTime();
689 >  if( globals->haveStatusTime() ){
690 >    info->statusTime = globals->getStatusTime();
691    }
692  
693 <  if( the_globals->haveThermalTime() ){
694 <    info->thermalTime = the_globals->getThermalTime();
693 >  if( globals->haveThermalTime() ){
694 >    info->thermalTime = globals->getThermalTime();
695    }
696  
697    // check for the temperature set flag
698  
699 <  if( the_globals->haveTempSet() ) info->setTemp = the_globals->getTempSet();
699 >  if( globals->haveTempSet() ) info->setTemp = globals->getTempSet();
700  
701    // get some of the tricky things that may still be in the globals
702  
703    double boxVector[3];
704 <  if( the_globals->haveBox() ){
705 <    boxVector[0] = the_globals->getBox();
706 <    boxVector[1] = the_globals->getBox();
707 <    boxVector[2] = the_globals->getBox();
704 >  if( globals->haveBox() ){
705 >    boxVector[0] = globals->getBox();
706 >    boxVector[1] = globals->getBox();
707 >    boxVector[2] = globals->getBox();
708      
709      info->setBox( boxVector );
710    }
711 <  else if( the_globals->haveDensity() ){
711 >  else if( globals->haveDensity() ){
712  
713      double vol;
714 <    vol = (double)tot_nmol / the_globals->getDensity();
714 >    vol = (double)tot_nmol / globals->getDensity();
715       boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
716       boxVector[1] = boxVector[0];
717       boxVector[2] = boxVector[0];
# Line 925 | Line 719 | void SimSetup::gatherInfo( void ){
719      info->setBox( boxVector );
720    }
721    else{
722 <    if( !the_globals->haveBoxX() ){
722 >    if( !globals->haveBoxX() ){
723        sprintf( painCave.errMsg,
724                 "SimSetup error, no periodic BoxX size given.\n" );
725        painCave.isFatal = 1;
726        simError();
727      }
728 <    boxVector[0] = the_globals->getBoxX();
728 >    boxVector[0] = globals->getBoxX();
729  
730 <    if( !the_globals->haveBoxY() ){
730 >    if( !globals->haveBoxY() ){
731        sprintf( painCave.errMsg,
732                 "SimSetup error, no periodic BoxY size given.\n" );
733        painCave.isFatal = 1;
734        simError();
735      }
736 <    boxVector[1] = the_globals->getBoxY();
736 >    boxVector[1] = globals->getBoxY();
737  
738 <    if( !the_globals->haveBoxZ() ){
738 >    if( !globals->haveBoxZ() ){
739        sprintf( painCave.errMsg,
740                 "SimSetup error, no periodic BoxZ size given.\n" );
741        painCave.isFatal = 1;
742        simError();
743      }
744 <    boxVector[2] = the_globals->getBoxZ();
744 >    boxVector[2] = globals->getBoxZ();
745  
746      info->setBox( boxVector );
747    }
# Line 977 | Line 771 | void SimSetup::finalInfoCheck( void ){
771    }
772    
773   #ifdef IS_MPI
774 <  int myUse = usesDipoles
775 <  MPI_Allreduce( &myUse, &UsesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
774 >  int myUse = usesDipoles;
775 >  MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
776   #endif //is_mpi
777  
778 +  double theEcr, theEst;
779  
780 <  if (the_globals->getUseRF() ) {
780 >  if (globals->getUseRF() ) {
781      info->useReactionField = 1;
782      
783 <    if( !the_globals->haveECR() ){
783 >    if( !globals->haveECR() ){
784        sprintf( painCave.errMsg,
785                 "SimSetup Warning: using default value of 1/2 the smallest "
786                 "box length for the electrostaticCutoffRadius.\n"
# Line 993 | Line 788 | void SimSetup::finalInfoCheck( void ){
788        painCave.isFatal = 0;
789        simError();
790        double smallest;
791 <      smallest = info->boxLx;
792 <      if (info->boxLy <= smallest) smallest = info->boxLy;
793 <      if (info->boxLz <= smallest) smallest = info->boxLz;
794 <      info->ecr = 0.5 * smallest;
791 >      smallest = info->boxL[0];
792 >      if (info->boxL[1] <= smallest) smallest = info->boxL[1];
793 >      if (info->boxL[2] <= smallest) smallest = info->boxL[2];
794 >      theEcr = 0.5 * smallest;
795      } else {
796 <      info->ecr        = the_globals->getECR();
796 >      theEcr = globals->getECR();
797      }
798  
799 <    if( !the_globals->haveEST() ){
799 >    if( !globals->haveEST() ){
800        sprintf( painCave.errMsg,
801                 "SimSetup Warning: using default value of 0.05 * the "
802                 "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
803                 );
804        painCave.isFatal = 0;
805        simError();
806 <      info->est = 0.05 * info->ecr;
806 >      theEst = 0.05 * theEcr;
807      } else {
808 <      info->est        = the_globals->getEST();
808 >      theEst= globals->getEST();
809      }
810 +
811 +    info->setEcr( theEcr, theEst );
812      
813 <    if(!the_globals->haveDielectric() ){
813 >    if(!globals->haveDielectric() ){
814        sprintf( painCave.errMsg,
815                 "SimSetup Error: You are trying to use Reaction Field without"
816                 "setting a dielectric constant!\n"
# Line 1021 | Line 818 | void SimSetup::finalInfoCheck( void ){
818        painCave.isFatal = 1;
819        simError();
820      }
821 <    info->dielectric = the_globals->getDielectric();  
821 >    info->dielectric = globals->getDielectric();  
822    }
823    else {
824      if (usesDipoles) {
825        
826 <      if( !the_globals->haveECR() ){
827 <        sprintf( painCave.errMsg,
828 <                 "SimSetup Warning: using default value of 1/2 the smallest "
829 <                 "box length for the electrostaticCutoffRadius.\n"
830 <                 "I hope you have a very fast processor!\n");
831 <        painCave.isFatal = 0;
832 <        simError();
833 <        double smallest;
834 <        smallest = info->boxLx;
835 <        if (info->boxLy <= smallest) smallest = info->boxLy;
836 <        if (info->boxLz <= smallest) smallest = info->boxLz;
837 <        info->ecr = 0.5 * smallest;
826 >      if( !globals->haveECR() ){
827 >        sprintf( painCave.errMsg,
828 >                 "SimSetup Warning: using default value of 1/2 the smallest "
829 >                 "box length for the electrostaticCutoffRadius.\n"
830 >                 "I hope you have a very fast processor!\n");
831 >        painCave.isFatal = 0;
832 >        simError();
833 >        double smallest;
834 >        smallest = info->boxL[0];
835 >        if (info->boxL[1] <= smallest) smallest = info->boxL[1];
836 >        if (info->boxL[2] <= smallest) smallest = info->boxL[2];
837 >        theEcr = 0.5 * smallest;
838        } else {
839 <        info->ecr        = the_globals->getECR();
839 >        theEcr = globals->getECR();
840        }
841        
842 <      if( !the_globals->haveEST() ){
843 <        sprintf( painCave.errMsg,
844 <                 "SimSetup Warning: using default value of 5%% of the "
845 <                 "electrostaticCutoffRadius for the "
846 <                 "electrostaticSkinThickness\n"
847 <                 );
848 <        painCave.isFatal = 0;
849 <        simError();
850 <        info->est = 0.05 * info->ecr;
842 >      if( !globals->haveEST() ){
843 >        sprintf( painCave.errMsg,
844 >                 "SimSetup Warning: using default value of 0.05 * the "
845 >                 "electrostaticCutoffRadius for the "
846 >                 "electrostaticSkinThickness\n"
847 >                 );
848 >        painCave.isFatal = 0;
849 >        simError();
850 >        theEst = 0.05 * theEcr;
851        } else {
852 <        info->est        = the_globals->getEST();
852 >        theEst= globals->getEST();
853        }
854 +
855 +      info->setEcr( theEcr, theEst );
856      }
857    }  
858  
# Line 1066 | Line 865 | void SimSetup::initSystemCoords( void ){
865  
866   void SimSetup::initSystemCoords( void ){
867  
868 < if( the_globals->haveInitialConfig() ){
868 > if( globals->haveInitialConfig() ){
869  
870       InitializeFromFile* fileInit;
871   #ifdef IS_MPI // is_mpi
872       if( worldRank == 0 ){
873   #endif //is_mpi
874 <   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
874 >   fileInit = new InitializeFromFile( globals->getInitialConfig() );
875   #ifdef IS_MPI
876       }else fileInit = new InitializeFromFile( NULL );
877   #endif
878 <   fileInit->read_xyz( info ); // default velocities on
878 >   fileInit->readInit( info ); // default velocities on
879  
880     delete fileInit;
881   }
# Line 1113 | Line 912 | void SimSetup::makeOutNames( void ){
912    if( worldRank == 0 ){
913   #endif // is_mpi
914      
915 <    if( the_globals->haveFinalConfig() ){
916 <      strcpy( info->finalName, the_globals->getFinalConfig() );
915 >    if( globals->haveFinalConfig() ){
916 >      strcpy( info->finalName, globals->getFinalConfig() );
917      }
918      else{
919        strcpy( info->finalName, inFileName );
# Line 1197 | Line 996 | void SimSetup::sysObjectsCreation( void ){
996  
997   void SimSetup::sysObjectsCreation( void ){
998  
999 +  int i;
1000 +
1001    // create the forceField
1002  
1003    createFF();
# Line 1219 | Line 1020 | void SimSetup::sysObjectsCreation( void ){
1020    
1021    makeSysArrays();
1022  
1023 +  // make and initialize the molecules (all but atomic coordinates)
1024    
1025 +  makeMolecules();
1026 +  info->identArray = new int[info->n_atoms];
1027 +  for(i=0; i<info->n_atoms; i++){
1028 +    info->identArray[i] = the_atoms[i]->getIdent();
1029 +  }
1030 +  
1031  
1032  
1033   }
# Line 1237 | Line 1045 | void SimSetup::createFF( void ){
1045      the_ff = new LJFF();
1046      break;
1047  
1048 +  case FF_EAM:
1049 +    the_ff = new EAM_FF();
1050 +    break;
1051 +
1052    default:
1053      sprintf( painCave.errMsg,
1054               "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1254 | Line 1066 | void SimSetup::compList( void ){
1066  
1067   void SimSetup::compList( void ){
1068  
1069 +  int i;
1070 +
1071    comp_stamps = new MoleculeStamp*[n_components];
1072  
1073    // make an array of molecule stamps that match the components used.
# Line 1279 | Line 1093 | void SimSetup::compList( void ){
1093        
1094        // extract the component from the list;
1095        
1096 <      currentStamp = the_stamps->extractMolStamp( id );
1096 >      currentStamp = stamps->extractMolStamp( id );
1097        if( currentStamp == NULL ){
1098          sprintf( painCave.errMsg,
1099                   "SimSetup error: Component \"%s\" was not found in the "
# Line 1303 | Line 1117 | void SimSetup::calcSysValues( void ){
1117   }
1118  
1119   void SimSetup::calcSysValues( void ){
1120 +  int i, j, k;
1121 +
1122  
1123    tot_atoms = 0;
1124    tot_bonds = 0;
# Line 1333 | Line 1149 | void SimSetup::mpiMolDivide( void ){
1149  
1150   void SimSetup::mpiMolDivide( void ){
1151    
1152 +  int i, j, k;
1153    int localMol, allMol;
1154    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1155  
# Line 1402 | Line 1219 | void SimSetup::makeSysArrays( void ){
1219  
1220  
1221   void SimSetup::makeSysArrays( void ){
1222 +  int i, j, k;
1223  
1224 +
1225    // create the atom and short range interaction arrays
1226  
1227    Atom::createArrays(info->n_atoms);
# Line 1474 | Line 1293 | void SimSetup::makeSysArrays( void ){
1293    info->excludes = the_excludes;
1294  
1295    the_ff->setSimInfo( info );
1296 +
1297 + }
1298 +
1299 + void SimSetup::makeIntegrator( void ){
1300 +
1301 +  NVT<RealIntegrator>*  myNVT = NULL;
1302 +  NPTi<RealIntegrator>* myNPTi = NULL;
1303 +  NPTf<RealIntegrator>* myNPTf = NULL;
1304 +  NPTim<RealIntegrator>* myNPTim = NULL;
1305 +  NPTfm<RealIntegrator>* myNPTfm = NULL;
1306 +
1307 +  switch( ensembleCase ){
1308 +
1309 +  case NVE_ENS:
1310 +    new NVE<RealIntegrator>( info, the_ff );
1311 +    break;
1312 +
1313 +  case NVT_ENS:
1314 +    myNVT = new NVT<RealIntegrator>( info, the_ff );
1315 +    myNVT->setTargetTemp(globals->getTargetTemp());
1316 +
1317 +    if (globals->haveTauThermostat())
1318 +      myNVT->setTauThermostat(globals->getTauThermostat());
1319 +
1320 +    else {
1321 +      sprintf( painCave.errMsg,
1322 +               "SimSetup error: If you use the NVT\n"
1323 +               "    ensemble, you must set tauThermostat.\n");
1324 +      painCave.isFatal = 1;
1325 +      simError();
1326 +    }
1327 +    break;
1328 +
1329 +  case NPTi_ENS:
1330 +    myNPTi = new NPTi<RealIntegrator>( info, the_ff );
1331 +    myNPTi->setTargetTemp( globals->getTargetTemp() );
1332 +
1333 +    if (globals->haveTargetPressure())
1334 +      myNPTi->setTargetPressure(globals->getTargetPressure());
1335 +    else {
1336 +      sprintf( painCave.errMsg,
1337 +               "SimSetup error: If you use a constant pressure\n"
1338 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1339 +      painCave.isFatal = 1;
1340 +      simError();
1341 +    }
1342 +    
1343 +    if( globals->haveTauThermostat() )
1344 +      myNPTi->setTauThermostat( globals->getTauThermostat() );
1345 +    else{
1346 +      sprintf( painCave.errMsg,
1347 +               "SimSetup error: If you use an NPT\n"
1348 +               "    ensemble, you must set tauThermostat.\n");
1349 +      painCave.isFatal = 1;
1350 +      simError();
1351 +    }
1352  
1353 +    if( globals->haveTauBarostat() )
1354 +      myNPTi->setTauBarostat( globals->getTauBarostat() );
1355 +    else{
1356 +      sprintf( painCave.errMsg,
1357 +               "SimSetup error: If you use an NPT\n"
1358 +               "    ensemble, you must set tauBarostat.\n");
1359 +      painCave.isFatal = 1;
1360 +      simError();
1361 +    }
1362 +    break;
1363 +
1364 +  case NPTf_ENS:
1365 +    myNPTf = new NPTf<RealIntegrator>( info, the_ff );
1366 +    myNPTf->setTargetTemp( globals->getTargetTemp());
1367 +
1368 +    if (globals->haveTargetPressure())
1369 +      myNPTf->setTargetPressure(globals->getTargetPressure());
1370 +    else {
1371 +      sprintf( painCave.errMsg,
1372 +               "SimSetup error: If you use a constant pressure\n"
1373 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1374 +      painCave.isFatal = 1;
1375 +      simError();
1376 +    }    
1377 +
1378 +    if( globals->haveTauThermostat() )
1379 +      myNPTf->setTauThermostat( globals->getTauThermostat() );
1380 +    else{
1381 +      sprintf( painCave.errMsg,
1382 +               "SimSetup error: If you use an NPT\n"
1383 +               "    ensemble, you must set tauThermostat.\n");
1384 +      painCave.isFatal = 1;
1385 +      simError();
1386 +    }
1387 +
1388 +    if( globals->haveTauBarostat() )
1389 +      myNPTf->setTauBarostat( globals->getTauBarostat() );
1390 +    else{
1391 +      sprintf( painCave.errMsg,
1392 +               "SimSetup error: If you use an NPT\n"
1393 +               "    ensemble, you must set tauBarostat.\n");
1394 +      painCave.isFatal = 1;
1395 +      simError();
1396 +    }
1397 +    break;
1398 +    
1399 +  case NPTim_ENS:
1400 +    myNPTim = new NPTim<RealIntegrator>( info, the_ff );
1401 +    myNPTim->setTargetTemp( globals->getTargetTemp());
1402 +
1403 +    if (globals->haveTargetPressure())
1404 +      myNPTim->setTargetPressure(globals->getTargetPressure());
1405 +    else {
1406 +      sprintf( painCave.errMsg,
1407 +               "SimSetup error: If you use a constant pressure\n"
1408 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1409 +      painCave.isFatal = 1;
1410 +      simError();
1411 +    }
1412 +    
1413 +    if( globals->haveTauThermostat() )
1414 +      myNPTim->setTauThermostat( globals->getTauThermostat() );
1415 +    else{
1416 +      sprintf( painCave.errMsg,
1417 +               "SimSetup error: If you use an NPT\n"
1418 +               "    ensemble, you must set tauThermostat.\n");
1419 +      painCave.isFatal = 1;
1420 +      simError();
1421 +    }
1422 +
1423 +    if( globals->haveTauBarostat() )
1424 +      myNPTim->setTauBarostat( globals->getTauBarostat() );
1425 +    else{
1426 +      sprintf( painCave.errMsg,
1427 +               "SimSetup error: If you use an NPT\n"
1428 +               "    ensemble, you must set tauBarostat.\n");
1429 +      painCave.isFatal = 1;
1430 +      simError();
1431 +    }
1432 +    break;
1433 +
1434 +  case NPTfm_ENS:
1435 +    myNPTfm = new NPTfm<RealIntegrator>( info, the_ff );
1436 +    myNPTfm->setTargetTemp( globals->getTargetTemp());
1437 +
1438 +    if (globals->haveTargetPressure())
1439 +      myNPTfm->setTargetPressure(globals->getTargetPressure());
1440 +    else {
1441 +      sprintf( painCave.errMsg,
1442 +               "SimSetup error: If you use a constant pressure\n"
1443 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1444 +      painCave.isFatal = 1;
1445 +      simError();
1446 +    }
1447 +    
1448 +    if( globals->haveTauThermostat() )
1449 +      myNPTfm->setTauThermostat( globals->getTauThermostat() );
1450 +    else{
1451 +      sprintf( painCave.errMsg,
1452 +               "SimSetup error: If you use an NPT\n"
1453 +               "    ensemble, you must set tauThermostat.\n");
1454 +      painCave.isFatal = 1;
1455 +      simError();
1456 +    }
1457 +
1458 +    if( globals->haveTauBarostat() )
1459 +      myNPTfm->setTauBarostat( globals->getTauBarostat() );
1460 +    else{
1461 +      sprintf( painCave.errMsg,
1462 +               "SimSetup error: If you use an NPT\n"
1463 +               "    ensemble, you must set tauBarostat.\n");
1464 +      painCave.isFatal = 1;
1465 +      simError();
1466 +    }
1467 +    break;
1468 +
1469 +  default:
1470 +    sprintf( painCave.errMsg,
1471 +             "SimSetup Error. Unrecognized ensemble in case statement.\n");
1472 +    painCave.isFatal = 1;
1473 +    simError();
1474 +  }
1475 +
1476   }
1477 +
1478 + void SimSetup::initFortran( void ){
1479 +
1480 +  info->refreshSim();
1481 +  
1482 +  if( !strcmp( info->mixingRule, "standard") ){
1483 +    the_ff->initForceField( LB_MIXING_RULE );
1484 +  }
1485 +  else if( !strcmp( info->mixingRule, "explicit") ){
1486 +    the_ff->initForceField( EXPLICIT_MIXING_RULE );
1487 +  }
1488 +  else{
1489 +    sprintf( painCave.errMsg,
1490 +             "SimSetup Error: unknown mixing rule -> \"%s\"\n",
1491 +             info->mixingRule );
1492 +    painCave.isFatal = 1;
1493 +    simError();
1494 +  }
1495 +
1496 +
1497 + #ifdef IS_MPI
1498 +  strcpy( checkPointMsg,
1499 +          "Successfully intialized the mixingRule for Fortran." );
1500 +  MPIcheckPoint();
1501 + #endif // is_mpi
1502 +
1503 + }

Diff Legend

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