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 626 by mmeineke, Wed Jul 16 21:30:56 2003 UTC

# Line 78 | Line 78 | void SimSetup::createSim( void ){
78  
79   void SimSetup::createSim( void ){
80  
81 –  MakeStamps *the_stamps;
82 –  Globals* the_globals;
81    int i, j, k, globalAtomIndex;
82    
83    // gather all of the information from the Bass file
# Line 90 | Line 88 | void SimSetup::createSim( void ){
88  
89    sysObjectsCreation();
90  
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 –
91    // check on the post processing info
92    
93    finalInfoCheck();
94  
110 –
111 –
112 –
95    // initialize the system coordinates
96  
97    initSystemCoords();
# Line 119 | Line 101 | void SimSetup::createSim( void ){
101  
102    makeOutNames();
103    
122 –
123 –  
124 –
125 –  
126 –
127 –  
128 –
129 –
104    // make the integrator
105    
106 +  makeIntegrator();
107    
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 –
108   #ifdef IS_MPI
109    mpiSim->mpiRefresh();
110   #endif
111  
112    // initialize the Fortran
113  
114 +  initFortran();
115  
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 –  }
116  
117  
333 – #ifdef IS_MPI
334 –  strcpy( checkPointMsg,
335 –          "Successfully intialized the mixingRule for Fortran." );
336 –  MPIcheckPoint();
337 – #endif // is_mpi
118   }
119  
120  
121   void SimSetup::makeMolecules( void ){
122  
123    int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
124 <  molInit info;
124 >  molInit molInfo;
125    DirectionalAtom* dAtom;
126    LinkedAssign* extras;
127    LinkedAssign* current_extra;
# Line 370 | Line 150 | void SimSetup::makeMolecules( void ){
150      
151      stampID = the_molecules[i].getStampID();
152  
153 <    info.nAtoms    = comp_stamps[stampID]->getNAtoms();
154 <    info.nBonds    = comp_stamps[stampID]->getNBonds();
155 <    info.nBends    = comp_stamps[stampID]->getNBends();
156 <    info.nTorsions = comp_stamps[stampID]->getNTorsions();
157 <    info.nExcludes = info.nBonds + info.nBends + info.nTorsions;
153 >    molInfo.nAtoms    = comp_stamps[stampID]->getNAtoms();
154 >    molInfo.nBonds    = comp_stamps[stampID]->getNBonds();
155 >    molInfo.nBends    = comp_stamps[stampID]->getNBends();
156 >    molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
157 >    molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
158  
159 <    info.myAtoms = &the_atoms[atomOffset];
160 <    info.myExcludes = &the_excludes[excludeOffset];
161 <    info.myBonds = new Bond*[info.nBonds];
162 <    info.myBends = new Bend*[info.nBends];
163 <    info.myTorsions = new Torsion*[info.nTorsions];
159 >    molInfo.myAtoms = &the_atoms[atomOffset];
160 >    molInfo.myExcludes = &the_excludes[excludeOffset];
161 >    molInfo.myBonds = new Bond*[molInfo.nBonds];
162 >    molInfo.myBends = new Bend*[molInfo.nBends];
163 >    molInfo.myTorsions = new Torsion*[molInfo.nTorsions];
164  
165 <    theBonds = new bond_pair[info.nBonds];
166 <    theBends = new bend_set[info.nBends];
167 <    theTorsions = new torsion_set[info.nTorsions];
165 >    theBonds = new bond_pair[molInfo.nBonds];
166 >    theBends = new bend_set[molInfo.nBends];
167 >    theTorsions = new torsion_set[molInfo.nTorsions];
168      
169      // make the Atoms
170      
171 <    for(j=0; j<info.nAtoms; j++){
171 >    for(j=0; j<molInfo.nAtoms; j++){
172        
173        currentAtom = comp_stamps[stampID]->getAtom( j );
174        if( currentAtom->haveOrientation() ){
175          
176          dAtom = new DirectionalAtom(j + atomOffset);
177          info->n_oriented++;
178 <        info.myAtoms[j] = dAtom;
178 >        molInfo.myAtoms[j] = dAtom;
179          
180          ux = currentAtom->getOrntX();
181          uy = currentAtom->getOrntY();
# Line 413 | Line 193 | void SimSetup::makeMolecules( void ){
193          dAtom->setSUz( uz );
194        }
195        else{
196 <        info.myAtoms[j] = new GeneralAtom(j + atomOffset);
196 >        molInfo.myAtoms[j] = new GeneralAtom(j + atomOffset);
197        }
198 <      info.myAtoms[j]->setType( currentAtom->getType() );
198 >      molInfo.myAtoms[j]->setType( currentAtom->getType() );
199      
200   #ifdef IS_MPI
201        
202 <      info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
202 >      molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
203        
204   #endif // is_mpi
205      }
206      
207      // make the bonds
208 <    for(j=0; j<info.nBonds; j++){
208 >    for(j=0; j<molInfo.nBonds; j++){
209        
210        currentBond = comp_stamps[stampID]->getBond( j );
211        theBonds[j].a = currentBond->getA() + atomOffset;
# Line 452 | Line 232 | void SimSetup::makeMolecules( void ){
232        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
233   #endif  //is_mpi
234      }
235 <    excludeOffset += info.nBonds;
235 >    excludeOffset += molInfo.nBonds;
236  
237      //make the bends
238 <    for(j=0; j<info.nBends; j++){
238 >    for(j=0; j<molInfo.nBends; j++){
239        
240        currentBend = comp_stamps[stampID]->getBend( j );
241        theBends[j].a = currentBend->getA() + atomOffset;
# Line 538 | Line 318 | void SimSetup::makeMolecules( void ){
318        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
319   #endif  //is_mpi
320      }
321 <    excludeOffset += info.nBends;
321 >    excludeOffset += molInfo.nBends;
322  
323 <    for(j=0; j<info.nTorsions; j++){
323 >    for(j=0; j<molInfo.nTorsions; j++){
324        
325        currentTorsion = comp_stamps[stampID]->getTorsion( j );
326        theTorsions[j].a = currentTorsion->getA() + atomOffset;
# Line 568 | Line 348 | void SimSetup::makeMolecules( void ){
348        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
349   #endif  //is_mpi
350      }
351 <    excludeOffset += info.nTorsions;
351 >    excludeOffset += molInfo.nTorsions;
352  
353      
354      // send the arrays off to the forceField for init.
355  
356 <    the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
357 <    the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
358 <    the_ff->initializeBends( info.nBends, info.myBends, theBends );
359 <    the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
356 >    the_ff->initializeAtoms( molInfo.nAtoms, molInfo.myAtoms );
357 >    the_ff->initializeBonds( molInfo.nBonds, molInfo.myBonds, theBonds );
358 >    the_ff->initializeBends( molInfo.nBends, molInfo.myBends, theBends );
359 >    the_ff->initializeTorsions( molInfo.nTorsions, molInfo.myTorsions, theTorsions );
360  
361  
362 <    the_molecules[i].initialize( info );
362 >    the_molecules[i].initialize( molInfo );
363  
364  
365 <    atomOffset += info.nAtoms;
365 >    atomOffset += molInfo.nAtoms;
366      delete[] theBonds;
367      delete[] theBends;
368      delete[] theTorsions;
# Line 618 | Line 398 | void SimSetup::initFromBass( void ){
398      have_extra =1;
399  
400      n_cells = (int)temp3 - 1;
401 <    cellx = info->boxLx / temp3;
402 <    celly = info->boxLy / temp3;
403 <    cellz = info->boxLz / temp3;
401 >    cellx = info->boxL[0] / temp3;
402 >    celly = info->boxL[1] / temp3;
403 >    cellz = info->boxL[2] / temp3;
404      n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
405      temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
406      n_per_extra = (int)ceil( temp1 );
# Line 635 | Line 415 | void SimSetup::initFromBass( void ){
415    }
416    else{
417      n_cells = (int)temp3;
418 <    cellx = info->boxLx / temp3;
419 <    celly = info->boxLy / temp3;
420 <    cellz = info->boxLz / temp3;
418 >    cellx = info->boxL[0] / temp3;
419 >    celly = info->boxL[1] / temp3;
420 >    cellz = info->boxL[2] / temp3;
421    }
422  
423    current_mol = 0;
# Line 784 | Line 564 | void SimSetup::gatherInfo( void ){
564  
565  
566   void SimSetup::gatherInfo( void ){
567 +  int i,j,k;
568  
569    ensembleCase = -1;
570    ffCase = -1;
571  
572    // get the stamps and globals;
573 <  the_stamps = stamps;
574 <  the_globals = globals;
573 >  stamps = stamps;
574 >  globals = globals;
575  
576    // set the easy ones first
577 <  info->target_temp = the_globals->getTargetTemp();
578 <  info->dt = the_globals->getDt();
579 <  info->run_time = the_globals->getRunTime();
580 <  n_components = the_globals->getNComponents();
577 >  info->target_temp = globals->getTargetTemp();
578 >  info->dt = globals->getDt();
579 >  info->run_time = globals->getRunTime();
580 >  n_components = globals->getNComponents();
581  
582  
583    // get the forceField
584  
585 <  strcpy( force_field, the_globals->getForceField() );
585 >  strcpy( force_field, globals->getForceField() );
586  
587    if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
588    else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
# Line 815 | Line 596 | void SimSetup::gatherInfo( void ){
596  
597    // get the ensemble
598  
599 <  strcpy( ensemble, the_globals->getEnsemble() );
599 >  strcpy( ensemble, globals->getEnsemble() );
600  
601    if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS;
602    else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
# Line 838 | Line 619 | void SimSetup::gatherInfo( void ){
619  
620    // get the mixing rule
621  
622 <  strcpy( info->mixingRule, the_globals->getMixingRule() );
623 <  info->usePBC = the_globals->getPBC();
622 >  strcpy( info->mixingRule, globals->getMixingRule() );
623 >  info->usePBC = globals->getPBC();
624          
625    
626    // get the components and calculate the tot_nMol and indvidual n_mol
627  
628 <  the_components = the_globals->getComponents();
628 >  the_components = globals->getComponents();
629    components_nmol = new int[n_components];
630  
631  
632 <  if( !the_globals->haveNMol() ){
632 >  if( !globals->haveNMol() ){
633      // we don't have the total number of molecules, so we assume it is
634      // given in each component
635  
# Line 881 | Line 662 | void SimSetup::gatherInfo( void ){
662  
663    // set the status, sample, and thermal kick times
664    
665 <  if( the_globals->haveSampleTime() ){
666 <    info->sampleTime = the_globals->getSampleTime();
665 >  if( globals->haveSampleTime() ){
666 >    info->sampleTime = globals->getSampleTime();
667      info->statusTime = info->sampleTime;
668      info->thermalTime = info->sampleTime;
669    }
670    else{
671 <    info->sampleTime = the_globals->getRunTime();
671 >    info->sampleTime = globals->getRunTime();
672      info->statusTime = info->sampleTime;
673      info->thermalTime = info->sampleTime;
674    }
675  
676 <  if( the_globals->haveStatusTime() ){
677 <    info->statusTime = the_globals->getStatusTime();
676 >  if( globals->haveStatusTime() ){
677 >    info->statusTime = globals->getStatusTime();
678    }
679  
680 <  if( the_globals->haveThermalTime() ){
681 <    info->thermalTime = the_globals->getThermalTime();
680 >  if( globals->haveThermalTime() ){
681 >    info->thermalTime = globals->getThermalTime();
682    }
683  
684    // check for the temperature set flag
685  
686 <  if( the_globals->haveTempSet() ) info->setTemp = the_globals->getTempSet();
686 >  if( globals->haveTempSet() ) info->setTemp = globals->getTempSet();
687  
688    // get some of the tricky things that may still be in the globals
689  
690    double boxVector[3];
691 <  if( the_globals->haveBox() ){
692 <    boxVector[0] = the_globals->getBox();
693 <    boxVector[1] = the_globals->getBox();
694 <    boxVector[2] = the_globals->getBox();
691 >  if( globals->haveBox() ){
692 >    boxVector[0] = globals->getBox();
693 >    boxVector[1] = globals->getBox();
694 >    boxVector[2] = globals->getBox();
695      
696      info->setBox( boxVector );
697    }
698 <  else if( the_globals->haveDensity() ){
698 >  else if( globals->haveDensity() ){
699  
700      double vol;
701 <    vol = (double)tot_nmol / the_globals->getDensity();
701 >    vol = (double)tot_nmol / globals->getDensity();
702       boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
703       boxVector[1] = boxVector[0];
704       boxVector[2] = boxVector[0];
# Line 925 | Line 706 | void SimSetup::gatherInfo( void ){
706      info->setBox( boxVector );
707    }
708    else{
709 <    if( !the_globals->haveBoxX() ){
709 >    if( !globals->haveBoxX() ){
710        sprintf( painCave.errMsg,
711                 "SimSetup error, no periodic BoxX size given.\n" );
712        painCave.isFatal = 1;
713        simError();
714      }
715 <    boxVector[0] = the_globals->getBoxX();
715 >    boxVector[0] = globals->getBoxX();
716  
717 <    if( !the_globals->haveBoxY() ){
717 >    if( !globals->haveBoxY() ){
718        sprintf( painCave.errMsg,
719                 "SimSetup error, no periodic BoxY size given.\n" );
720        painCave.isFatal = 1;
721        simError();
722      }
723 <    boxVector[1] = the_globals->getBoxY();
723 >    boxVector[1] = globals->getBoxY();
724  
725 <    if( !the_globals->haveBoxZ() ){
725 >    if( !globals->haveBoxZ() ){
726        sprintf( painCave.errMsg,
727                 "SimSetup error, no periodic BoxZ size given.\n" );
728        painCave.isFatal = 1;
729        simError();
730      }
731 <    boxVector[2] = the_globals->getBoxZ();
731 >    boxVector[2] = globals->getBoxZ();
732  
733      info->setBox( boxVector );
734    }
# Line 977 | Line 758 | void SimSetup::finalInfoCheck( void ){
758    }
759    
760   #ifdef IS_MPI
761 <  int myUse = usesDipoles
762 <  MPI_Allreduce( &myUse, &UsesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
761 >  int myUse = usesDipoles;
762 >  MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
763   #endif //is_mpi
764  
765 +  double theEcr, theEst;
766  
767 <  if (the_globals->getUseRF() ) {
767 >  if (globals->getUseRF() ) {
768      info->useReactionField = 1;
769      
770 <    if( !the_globals->haveECR() ){
770 >    if( !globals->haveECR() ){
771        sprintf( painCave.errMsg,
772                 "SimSetup Warning: using default value of 1/2 the smallest "
773                 "box length for the electrostaticCutoffRadius.\n"
# Line 993 | Line 775 | void SimSetup::finalInfoCheck( void ){
775        painCave.isFatal = 0;
776        simError();
777        double smallest;
778 <      smallest = info->boxLx;
779 <      if (info->boxLy <= smallest) smallest = info->boxLy;
780 <      if (info->boxLz <= smallest) smallest = info->boxLz;
781 <      info->ecr = 0.5 * smallest;
778 >      smallest = info->boxL[0];
779 >      if (info->boxL[1] <= smallest) smallest = info->boxL[1];
780 >      if (info->boxL[2] <= smallest) smallest = info->boxL[2];
781 >      theEcr = 0.5 * smallest;
782      } else {
783 <      info->ecr        = the_globals->getECR();
783 >      theEcr = globals->getECR();
784      }
785  
786 <    if( !the_globals->haveEST() ){
786 >    if( !globals->haveEST() ){
787        sprintf( painCave.errMsg,
788                 "SimSetup Warning: using default value of 0.05 * the "
789                 "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
790                 );
791        painCave.isFatal = 0;
792        simError();
793 <      info->est = 0.05 * info->ecr;
793 >      theEst = 0.05 * theEcr;
794      } else {
795 <      info->est        = the_globals->getEST();
795 >      theEst= globals->getEST();
796      }
797 +
798 +    info->setEcr( theEcr, theEst );
799      
800 <    if(!the_globals->haveDielectric() ){
800 >    if(!globals->haveDielectric() ){
801        sprintf( painCave.errMsg,
802                 "SimSetup Error: You are trying to use Reaction Field without"
803                 "setting a dielectric constant!\n"
# Line 1021 | Line 805 | void SimSetup::finalInfoCheck( void ){
805        painCave.isFatal = 1;
806        simError();
807      }
808 <    info->dielectric = the_globals->getDielectric();  
808 >    info->dielectric = globals->getDielectric();  
809    }
810    else {
811      if (usesDipoles) {
812        
813 <      if( !the_globals->haveECR() ){
814 <        sprintf( painCave.errMsg,
815 <                 "SimSetup Warning: using default value of 1/2 the smallest "
816 <                 "box length for the electrostaticCutoffRadius.\n"
817 <                 "I hope you have a very fast processor!\n");
818 <        painCave.isFatal = 0;
819 <        simError();
820 <        double smallest;
821 <        smallest = info->boxLx;
822 <        if (info->boxLy <= smallest) smallest = info->boxLy;
823 <        if (info->boxLz <= smallest) smallest = info->boxLz;
824 <        info->ecr = 0.5 * smallest;
813 >      if( !globals->haveECR() ){
814 >        sprintf( painCave.errMsg,
815 >                 "SimSetup Warning: using default value of 1/2 the smallest "
816 >                 "box length for the electrostaticCutoffRadius.\n"
817 >                 "I hope you have a very fast processor!\n");
818 >        painCave.isFatal = 0;
819 >        simError();
820 >        double smallest;
821 >        smallest = info->boxL[0];
822 >        if (info->boxL[1] <= smallest) smallest = info->boxL[1];
823 >        if (info->boxL[2] <= smallest) smallest = info->boxL[2];
824 >        theEcr = 0.5 * smallest;
825        } else {
826 <        info->ecr        = the_globals->getECR();
826 >        theEcr = globals->getECR();
827        }
828        
829 <      if( !the_globals->haveEST() ){
830 <        sprintf( painCave.errMsg,
831 <                 "SimSetup Warning: using default value of 5%% of the "
832 <                 "electrostaticCutoffRadius for the "
833 <                 "electrostaticSkinThickness\n"
834 <                 );
835 <        painCave.isFatal = 0;
836 <        simError();
837 <        info->est = 0.05 * info->ecr;
829 >      if( !globals->haveEST() ){
830 >        sprintf( painCave.errMsg,
831 >                 "SimSetup Warning: using default value of 0.05 * the "
832 >                 "electrostaticCutoffRadius for the "
833 >                 "electrostaticSkinThickness\n"
834 >                 );
835 >        painCave.isFatal = 0;
836 >        simError();
837 >        theEst = 0.05 * theEcr;
838        } else {
839 <        info->est        = the_globals->getEST();
839 >        theEst= globals->getEST();
840        }
841 +
842 +      info->setEcr( theEcr, theEst );
843      }
844    }  
845  
# Line 1066 | Line 852 | void SimSetup::initSystemCoords( void ){
852  
853   void SimSetup::initSystemCoords( void ){
854  
855 < if( the_globals->haveInitialConfig() ){
855 > if( globals->haveInitialConfig() ){
856  
857       InitializeFromFile* fileInit;
858   #ifdef IS_MPI // is_mpi
859       if( worldRank == 0 ){
860   #endif //is_mpi
861 <   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
861 >   fileInit = new InitializeFromFile( globals->getInitialConfig() );
862   #ifdef IS_MPI
863       }else fileInit = new InitializeFromFile( NULL );
864   #endif
# Line 1113 | Line 899 | void SimSetup::makeOutNames( void ){
899    if( worldRank == 0 ){
900   #endif // is_mpi
901      
902 <    if( the_globals->haveFinalConfig() ){
903 <      strcpy( info->finalName, the_globals->getFinalConfig() );
902 >    if( globals->haveFinalConfig() ){
903 >      strcpy( info->finalName, globals->getFinalConfig() );
904      }
905      else{
906        strcpy( info->finalName, inFileName );
# Line 1197 | Line 983 | void SimSetup::sysObjectsCreation( void ){
983  
984   void SimSetup::sysObjectsCreation( void ){
985  
986 +  int i;
987 +
988    // create the forceField
989  
990    createFF();
# Line 1219 | Line 1007 | void SimSetup::sysObjectsCreation( void ){
1007    
1008    makeSysArrays();
1009  
1010 +  // make and initialize the molecules (all but atomic coordinates)
1011    
1012 +  makeMolecules();
1013 +  info->identArray = new int[info->n_atoms];
1014 +  for(i=0; i<info->n_atoms; i++){
1015 +    info->identArray[i] = the_atoms[i]->getIdent();
1016 +  }
1017 +  
1018  
1019  
1020   }
# Line 1253 | Line 1048 | void SimSetup::compList( void ){
1048  
1049  
1050   void SimSetup::compList( void ){
1051 +
1052 +  int i;
1053  
1054    comp_stamps = new MoleculeStamp*[n_components];
1055  
# Line 1279 | Line 1076 | void SimSetup::compList( void ){
1076        
1077        // extract the component from the list;
1078        
1079 <      currentStamp = the_stamps->extractMolStamp( id );
1079 >      currentStamp = stamps->extractMolStamp( id );
1080        if( currentStamp == NULL ){
1081          sprintf( painCave.errMsg,
1082                   "SimSetup error: Component \"%s\" was not found in the "
# Line 1303 | Line 1100 | void SimSetup::calcSysValues( void ){
1100   }
1101  
1102   void SimSetup::calcSysValues( void ){
1103 +  int i, j, k;
1104  
1105 +
1106    tot_atoms = 0;
1107    tot_bonds = 0;
1108    tot_bends = 0;
# Line 1333 | Line 1132 | void SimSetup::mpiMolDivide( void ){
1132  
1133   void SimSetup::mpiMolDivide( void ){
1134    
1135 +  int i, j, k;
1136    int localMol, allMol;
1137    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1138  
# Line 1402 | Line 1202 | void SimSetup::makeSysArrays( void ){
1202  
1203  
1204   void SimSetup::makeSysArrays( void ){
1205 +  int i, j, k;
1206  
1207 +
1208    // create the atom and short range interaction arrays
1209  
1210    Atom::createArrays(info->n_atoms);
# Line 1476 | Line 1278 | void SimSetup::makeSysArrays( void ){
1278    the_ff->setSimInfo( info );
1279  
1280   }
1281 +
1282 + void SimSetup::makeIntegrator( void ){
1283 +
1284 +  NVT*  myNVT = NULL;
1285 +  NPTi* myNPTi = NULL;
1286 +  NPTf* myNPTf = NULL;
1287 +  NPTim* myNPTim = NULL;
1288 +  NPTfm* myNPTfm = NULL;
1289 +
1290 +  switch( ensembleCase ){
1291 +
1292 +  case NVE_ENS:
1293 +    new NVE( info, the_ff );
1294 +    break;
1295 +
1296 +  case NVT_ENS:
1297 +    myNVT = new NVT( info, the_ff );
1298 +    myNVT->setTargetTemp(globals->getTargetTemp());
1299 +
1300 +    if (globals->haveTauThermostat())
1301 +      myNVT->setTauThermostat(globals->getTauThermostat());
1302 +
1303 +    else {
1304 +      sprintf( painCave.errMsg,
1305 +               "SimSetup error: If you use the NVT\n"
1306 +               "    ensemble, you must set tauThermostat.\n");
1307 +      painCave.isFatal = 1;
1308 +      simError();
1309 +    }
1310 +    break;
1311 +
1312 +  case NPTi_ENS:
1313 +    myNPTi = new NPTi( info, the_ff );
1314 +    myNPTi->setTargetTemp( globals->getTargetTemp() );
1315 +
1316 +    if (globals->haveTargetPressure())
1317 +      myNPTi->setTargetPressure(globals->getTargetPressure());
1318 +    else {
1319 +      sprintf( painCave.errMsg,
1320 +               "SimSetup error: If you use a constant pressure\n"
1321 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1322 +      painCave.isFatal = 1;
1323 +      simError();
1324 +    }
1325 +    
1326 +    if( globals->haveTauThermostat() )
1327 +      myNPTi->setTauThermostat( globals->getTauThermostat() );
1328 +    else{
1329 +      sprintf( painCave.errMsg,
1330 +               "SimSetup error: If you use an NPT\n"
1331 +               "    ensemble, you must set tauThermostat.\n");
1332 +      painCave.isFatal = 1;
1333 +      simError();
1334 +    }
1335 +
1336 +    if( globals->haveTauBarostat() )
1337 +      myNPTi->setTauBarostat( globals->getTauBarostat() );
1338 +    else{
1339 +      sprintf( painCave.errMsg,
1340 +               "SimSetup error: If you use an NPT\n"
1341 +               "    ensemble, you must set tauBarostat.\n");
1342 +      painCave.isFatal = 1;
1343 +      simError();
1344 +    }
1345 +    break;
1346 +
1347 +  case NPTf_ENS:
1348 +    myNPTf = new NPTf( info, the_ff );
1349 +    myNPTf->setTargetTemp( globals->getTargetTemp());
1350 +
1351 +    if (globals->haveTargetPressure())
1352 +      myNPTf->setTargetPressure(globals->getTargetPressure());
1353 +    else {
1354 +      sprintf( painCave.errMsg,
1355 +               "SimSetup error: If you use a constant pressure\n"
1356 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1357 +      painCave.isFatal = 1;
1358 +      simError();
1359 +    }    
1360 +
1361 +    if( globals->haveTauThermostat() )
1362 +      myNPTf->setTauThermostat( globals->getTauThermostat() );
1363 +    else{
1364 +      sprintf( painCave.errMsg,
1365 +               "SimSetup error: If you use an NPT\n"
1366 +               "    ensemble, you must set tauThermostat.\n");
1367 +      painCave.isFatal = 1;
1368 +      simError();
1369 +    }
1370 +
1371 +    if( globals->haveTauBarostat() )
1372 +      myNPTf->setTauBarostat( globals->getTauBarostat() );
1373 +    else{
1374 +      sprintf( painCave.errMsg,
1375 +               "SimSetup error: If you use an NPT\n"
1376 +               "    ensemble, you must set tauBarostat.\n");
1377 +      painCave.isFatal = 1;
1378 +      simError();
1379 +    }
1380 +    break;
1381 +    
1382 +  case NPTim_ENS:
1383 +    myNPTim = new NPTim( info, the_ff );
1384 +    myNPTim->setTargetTemp( globals->getTargetTemp());
1385 +
1386 +    if (globals->haveTargetPressure())
1387 +      myNPTim->setTargetPressure(globals->getTargetPressure());
1388 +    else {
1389 +      sprintf( painCave.errMsg,
1390 +               "SimSetup error: If you use a constant pressure\n"
1391 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1392 +      painCave.isFatal = 1;
1393 +      simError();
1394 +    }
1395 +    
1396 +    if( globals->haveTauThermostat() )
1397 +      myNPTim->setTauThermostat( globals->getTauThermostat() );
1398 +    else{
1399 +      sprintf( painCave.errMsg,
1400 +               "SimSetup error: If you use an NPT\n"
1401 +               "    ensemble, you must set tauThermostat.\n");
1402 +      painCave.isFatal = 1;
1403 +      simError();
1404 +    }
1405 +
1406 +    if( globals->haveTauBarostat() )
1407 +      myNPTim->setTauBarostat( globals->getTauBarostat() );
1408 +    else{
1409 +      sprintf( painCave.errMsg,
1410 +               "SimSetup error: If you use an NPT\n"
1411 +               "    ensemble, you must set tauBarostat.\n");
1412 +      painCave.isFatal = 1;
1413 +      simError();
1414 +    }
1415 +    break;
1416 +
1417 +  case NPTfm_ENS:
1418 +    myNPTfm = new NPTfm( info, the_ff );
1419 +    myNPTfm->setTargetTemp( globals->getTargetTemp());
1420 +
1421 +    if (globals->haveTargetPressure())
1422 +      myNPTfm->setTargetPressure(globals->getTargetPressure());
1423 +    else {
1424 +      sprintf( painCave.errMsg,
1425 +               "SimSetup error: If you use a constant pressure\n"
1426 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1427 +      painCave.isFatal = 1;
1428 +      simError();
1429 +    }
1430 +    
1431 +    if( globals->haveTauThermostat() )
1432 +      myNPTfm->setTauThermostat( globals->getTauThermostat() );
1433 +    else{
1434 +      sprintf( painCave.errMsg,
1435 +               "SimSetup error: If you use an NPT\n"
1436 +               "    ensemble, you must set tauThermostat.\n");
1437 +      painCave.isFatal = 1;
1438 +      simError();
1439 +    }
1440 +
1441 +    if( globals->haveTauBarostat() )
1442 +      myNPTfm->setTauBarostat( globals->getTauBarostat() );
1443 +    else{
1444 +      sprintf( painCave.errMsg,
1445 +               "SimSetup error: If you use an NPT\n"
1446 +               "    ensemble, you must set tauBarostat.\n");
1447 +      painCave.isFatal = 1;
1448 +      simError();
1449 +    }
1450 +    break;
1451 +
1452 +  default:
1453 +    sprintf( painCave.errMsg,
1454 +             "SimSetup Error. Unrecognized ensemble in case statement.\n");
1455 +    painCave.isFatal = 1;
1456 +    simError();
1457 +  }
1458 +
1459 + }
1460 +
1461 + void SimSetup::initFortran( void ){
1462 +
1463 +  info->refreshSim();
1464 +  
1465 +  if( !strcmp( info->mixingRule, "standard") ){
1466 +    the_ff->initForceField( LB_MIXING_RULE );
1467 +  }
1468 +  else if( !strcmp( info->mixingRule, "explicit") ){
1469 +    the_ff->initForceField( EXPLICIT_MIXING_RULE );
1470 +  }
1471 +  else{
1472 +    sprintf( painCave.errMsg,
1473 +             "SimSetup Error: unknown mixing rule -> \"%s\"\n",
1474 +             info->mixingRule );
1475 +    painCave.isFatal = 1;
1476 +    simError();
1477 +  }
1478 +
1479 +
1480 + #ifdef IS_MPI
1481 +  strcpy( checkPointMsg,
1482 +          "Successfully intialized the mixingRule for Fortran." );
1483 +  MPIcheckPoint();
1484 + #endif // is_mpi
1485 +
1486 + }

Diff Legend

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