ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/no-template-branch/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 660 by tim, Thu Jul 31 19:59:34 2003 UTC

# Line 1 | Line 1
1 + #include <algorithm>
2   #include <cstdlib>
3   #include <iostream>
4   #include <cmath>
5 + #include <string>
6  
7   #include "SimSetup.hpp"
8   #include "parse_me.h"
# Line 14 | Line 16
16  
17   // some defines for ensemble and Forcefield  cases
18  
19 < #define NVE_ENS   0
20 < #define NVT_ENS   1
21 < #define NPTi_ENS  2
22 < #define NPTf_ENS  3
23 < #define NPTim_ENS 4
24 < #define NPTfm_ENS 5
19 > #define NVE_ENS        0
20 > #define NVT_ENS        1
21 > #define NPTi_ENS       2
22 > #define NPTf_ENS       3
23 > #define NPTim_ENS      4
24 > #define NPTfm_ENS      5
25 > #define NVEZCONS_ENS   6
26 > #define NVTZCONS_ENS   7
27 > #define NPTiZCONS_ENS  8
28 > #define NPTfZCONS_ENS  9
29 > #define NPTimZCONS_ENS 10
30 > #define NPTfmZCONS_ENS 11
31  
24 –
32   #define FF_DUFF 0
33   #define FF_LJ   1
34 + #define FF_EAM  2
35  
36 + using namespace std;
37  
38   SimSetup::SimSetup(){
39 +  
40 +  isInfoArray = 0;
41 +  nInfo = 1;
42 +  
43    stamps = new MakeStamps();
44    globals = new Globals();
45    
46 +  
47   #ifdef IS_MPI
48    strcpy( checkPointMsg, "SimSetup creation successful" );
49    MPIcheckPoint();
# Line 41 | Line 55 | SimSetup::~SimSetup(){
55    delete globals;
56   }
57  
58 + void SimSetup::setSimInfo( SimInfo* the_info, int theNinfo ) {
59 +    info = the_info;
60 +    nInfo = theNinfo;
61 +    isInfoArray = 1;
62 +  }
63 +
64 +
65   void SimSetup::parseFile( char* fileName ){
66  
67   #ifdef IS_MPI
# Line 78 | Line 99 | void SimSetup::createSim( void ){
99  
100   void SimSetup::createSim( void ){
101  
81 –  MakeStamps *the_stamps;
82 –  Globals* the_globals;
102    int i, j, k, globalAtomIndex;
103    
104    // gather all of the information from the Bass file
# Line 90 | Line 109 | void SimSetup::createSim( void ){
109  
110    sysObjectsCreation();
111  
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 –
112    // check on the post processing info
113    
114    finalInfoCheck();
115  
110 –
111 –
112 –
116    // initialize the system coordinates
117  
118    initSystemCoords();
# Line 119 | Line 122 | void SimSetup::createSim( void ){
122  
123    makeOutNames();
124    
122 –
123 –  
124 –
125 –  
126 –
127 –  
128 –
129 –
125    // make the integrator
126    
127 +  makeIntegrator();
128    
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 –
129   #ifdef IS_MPI
130    mpiSim->mpiRefresh();
131   #endif
132  
133    // initialize the Fortran
134  
135 +  initFortran();
136  
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 –  }
137  
138  
333 – #ifdef IS_MPI
334 –  strcpy( checkPointMsg,
335 –          "Successfully intialized the mixingRule for Fortran." );
336 –  MPIcheckPoint();
337 – #endif // is_mpi
139   }
140  
141  
142   void SimSetup::makeMolecules( void ){
143  
144    int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
145 <  molInit info;
145 >  molInit molInfo;
146    DirectionalAtom* dAtom;
147    LinkedAssign* extras;
148    LinkedAssign* current_extra;
# Line 370 | Line 171 | void SimSetup::makeMolecules( void ){
171      
172      stampID = the_molecules[i].getStampID();
173  
174 <    info.nAtoms    = comp_stamps[stampID]->getNAtoms();
175 <    info.nBonds    = comp_stamps[stampID]->getNBonds();
176 <    info.nBends    = comp_stamps[stampID]->getNBends();
177 <    info.nTorsions = comp_stamps[stampID]->getNTorsions();
178 <    info.nExcludes = info.nBonds + info.nBends + info.nTorsions;
174 >    molInfo.nAtoms    = comp_stamps[stampID]->getNAtoms();
175 >    molInfo.nBonds    = comp_stamps[stampID]->getNBonds();
176 >    molInfo.nBends    = comp_stamps[stampID]->getNBends();
177 >    molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
178 >    molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
179  
180 <    info.myAtoms = &the_atoms[atomOffset];
181 <    info.myExcludes = &the_excludes[excludeOffset];
182 <    info.myBonds = new Bond*[info.nBonds];
183 <    info.myBends = new Bend*[info.nBends];
184 <    info.myTorsions = new Torsion*[info.nTorsions];
180 >    molInfo.myAtoms = &the_atoms[atomOffset];
181 >    molInfo.myExcludes = &the_excludes[excludeOffset];
182 >    molInfo.myBonds = new Bond*[molInfo.nBonds];
183 >    molInfo.myBends = new Bend*[molInfo.nBends];
184 >    molInfo.myTorsions = new Torsion*[molInfo.nTorsions];
185  
186 <    theBonds = new bond_pair[info.nBonds];
187 <    theBends = new bend_set[info.nBends];
188 <    theTorsions = new torsion_set[info.nTorsions];
186 >    theBonds = new bond_pair[molInfo.nBonds];
187 >    theBends = new bend_set[molInfo.nBends];
188 >    theTorsions = new torsion_set[molInfo.nTorsions];
189      
190      // make the Atoms
191      
192 <    for(j=0; j<info.nAtoms; j++){
192 >    for(j=0; j<molInfo.nAtoms; j++){
193        
194        currentAtom = comp_stamps[stampID]->getAtom( j );
195        if( currentAtom->haveOrientation() ){
196          
197          dAtom = new DirectionalAtom(j + atomOffset);
198          info->n_oriented++;
199 <        info.myAtoms[j] = dAtom;
199 >        molInfo.myAtoms[j] = dAtom;
200          
201          ux = currentAtom->getOrntX();
202          uy = currentAtom->getOrntY();
# Line 413 | Line 214 | void SimSetup::makeMolecules( void ){
214          dAtom->setSUz( uz );
215        }
216        else{
217 <        info.myAtoms[j] = new GeneralAtom(j + atomOffset);
217 >        molInfo.myAtoms[j] = new GeneralAtom(j + atomOffset);
218        }
219 <      info.myAtoms[j]->setType( currentAtom->getType() );
219 >      molInfo.myAtoms[j]->setType( currentAtom->getType() );
220      
221   #ifdef IS_MPI
222        
223 <      info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
223 >      molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
224        
225   #endif // is_mpi
226      }
227      
228      // make the bonds
229 <    for(j=0; j<info.nBonds; j++){
229 >    for(j=0; j<molInfo.nBonds; j++){
230        
231        currentBond = comp_stamps[stampID]->getBond( j );
232        theBonds[j].a = currentBond->getA() + atomOffset;
# Line 452 | Line 253 | void SimSetup::makeMolecules( void ){
253        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
254   #endif  //is_mpi
255      }
256 <    excludeOffset += info.nBonds;
256 >    excludeOffset += molInfo.nBonds;
257  
258      //make the bends
259 <    for(j=0; j<info.nBends; j++){
259 >    for(j=0; j<molInfo.nBends; j++){
260        
261        currentBend = comp_stamps[stampID]->getBend( j );
262        theBends[j].a = currentBend->getA() + atomOffset;
# Line 538 | Line 339 | void SimSetup::makeMolecules( void ){
339        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
340   #endif  //is_mpi
341      }
342 <    excludeOffset += info.nBends;
342 >    excludeOffset += molInfo.nBends;
343  
344 <    for(j=0; j<info.nTorsions; j++){
344 >    for(j=0; j<molInfo.nTorsions; j++){
345        
346        currentTorsion = comp_stamps[stampID]->getTorsion( j );
347        theTorsions[j].a = currentTorsion->getA() + atomOffset;
# Line 568 | Line 369 | void SimSetup::makeMolecules( void ){
369        the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
370   #endif  //is_mpi
371      }
372 <    excludeOffset += info.nTorsions;
372 >    excludeOffset += molInfo.nTorsions;
373  
374      
375      // send the arrays off to the forceField for init.
376  
377 <    the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
378 <    the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
379 <    the_ff->initializeBends( info.nBends, info.myBends, theBends );
380 <    the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
377 >    the_ff->initializeAtoms( molInfo.nAtoms, molInfo.myAtoms );
378 >    the_ff->initializeBonds( molInfo.nBonds, molInfo.myBonds, theBonds );
379 >    the_ff->initializeBends( molInfo.nBends, molInfo.myBends, theBends );
380 >    the_ff->initializeTorsions( molInfo.nTorsions, molInfo.myTorsions, theTorsions );
381  
382  
383 <    the_molecules[i].initialize( info );
383 >    the_molecules[i].initialize( molInfo );
384  
385  
386 <    atomOffset += info.nAtoms;
386 >    atomOffset += molInfo.nAtoms;
387      delete[] theBonds;
388      delete[] theBends;
389      delete[] theTorsions;
# Line 618 | Line 419 | void SimSetup::initFromBass( void ){
419      have_extra =1;
420  
421      n_cells = (int)temp3 - 1;
422 <    cellx = info->boxLx / temp3;
423 <    celly = info->boxLy / temp3;
424 <    cellz = info->boxLz / temp3;
422 >    cellx = info->boxL[0] / temp3;
423 >    celly = info->boxL[1] / temp3;
424 >    cellz = info->boxL[2] / temp3;
425      n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
426      temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
427      n_per_extra = (int)ceil( temp1 );
# Line 635 | Line 436 | void SimSetup::initFromBass( void ){
436    }
437    else{
438      n_cells = (int)temp3;
439 <    cellx = info->boxLx / temp3;
440 <    celly = info->boxLy / temp3;
441 <    cellz = info->boxLz / temp3;
439 >    cellx = info->boxL[0] / temp3;
440 >    celly = info->boxL[1] / temp3;
441 >    cellz = info->boxL[2] / temp3;
442    }
443  
444    current_mol = 0;
# Line 784 | Line 585 | void SimSetup::gatherInfo( void ){
585  
586  
587   void SimSetup::gatherInfo( void ){
588 +  int i,j,k;
589  
590    ensembleCase = -1;
591    ffCase = -1;
592  
593    // get the stamps and globals;
594 <  the_stamps = stamps;
595 <  the_globals = globals;
594 >  stamps = stamps;
595 >  globals = globals;
596  
597    // set the easy ones first
598 <  info->target_temp = the_globals->getTargetTemp();
599 <  info->dt = the_globals->getDt();
600 <  info->run_time = the_globals->getRunTime();
601 <  n_components = the_globals->getNComponents();
598 >  info->target_temp = globals->getTargetTemp();
599 >  info->dt = globals->getDt();
600 >  info->run_time = globals->getRunTime();
601 >  n_components = globals->getNComponents();
602  
603  
604    // get the forceField
605  
606 <  strcpy( force_field, the_globals->getForceField() );
606 >  strcpy( force_field, globals->getForceField() );
607  
608    if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
609    else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
610 +  else if( !strcasecmp( force_field, "EAM" )) ffCase = FF_EAM;
611    else{
612      sprintf( painCave.errMsg,
613               "SimSetup Error. Unrecognized force field -> %s\n",
# Line 815 | Line 618 | void SimSetup::gatherInfo( void ){
618  
619    // get the ensemble
620  
621 <  strcpy( ensemble, the_globals->getEnsemble() );
621 >  strcpy( ensemble, globals->getEnsemble() );
622  
623    if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS;
624    else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
# Line 824 | Line 627 | void SimSetup::gatherInfo( void ){
627    else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
628    else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
629    else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS;
630 +
631 +  else if( !strcasecmp( ensemble, "NVEZCONS")) ensembleCase = NVEZCONS_ENS;
632 +  else if( !strcasecmp( ensemble, "NVTZCONS"))  ensembleCase = NVTZCONS_ENS;
633 +  else if( !strcasecmp( ensemble, "NPTiZCONS") || !strcasecmp( ensemble, "NPT"))
634 +    ensembleCase = NPTiZCONS_ENS;
635 +  else if( !strcasecmp( ensemble, "NPTfCONS"))  ensembleCase = NPTfZCONS_ENS;
636 +  else if( !strcasecmp( ensemble, "NPTimZCONS"))  ensembleCase = NPTimZCONS_ENS;
637 +  else if( !strcasecmp( ensemble, "NPTfmCONS"))  ensembleCase = NPTfmZCONS_ENS;
638 +  
639    else{
640      sprintf( painCave.errMsg,
641               "SimSetup Warning. Unrecognized Ensemble -> %s, "
# Line 838 | Line 650 | void SimSetup::gatherInfo( void ){
650  
651    // get the mixing rule
652  
653 <  strcpy( info->mixingRule, the_globals->getMixingRule() );
654 <  info->usePBC = the_globals->getPBC();
653 >  strcpy( info->mixingRule, globals->getMixingRule() );
654 >  info->usePBC = globals->getPBC();
655          
656    
657    // get the components and calculate the tot_nMol and indvidual n_mol
658  
659 <  the_components = the_globals->getComponents();
659 >  the_components = globals->getComponents();
660    components_nmol = new int[n_components];
661  
662  
663 <  if( !the_globals->haveNMol() ){
663 >  if( !globals->haveNMol() ){
664      // we don't have the total number of molecules, so we assume it is
665      // given in each component
666  
# Line 881 | Line 693 | void SimSetup::gatherInfo( void ){
693  
694    // set the status, sample, and thermal kick times
695    
696 <  if( the_globals->haveSampleTime() ){
697 <    info->sampleTime = the_globals->getSampleTime();
696 >  if( globals->haveSampleTime() ){
697 >    info->sampleTime = globals->getSampleTime();
698      info->statusTime = info->sampleTime;
699      info->thermalTime = info->sampleTime;
700    }
701    else{
702 <    info->sampleTime = the_globals->getRunTime();
702 >    info->sampleTime = globals->getRunTime();
703      info->statusTime = info->sampleTime;
704      info->thermalTime = info->sampleTime;
705    }
706  
707 <  if( the_globals->haveStatusTime() ){
708 <    info->statusTime = the_globals->getStatusTime();
707 >  if( globals->haveStatusTime() ){
708 >    info->statusTime = globals->getStatusTime();
709    }
710  
711 <  if( the_globals->haveThermalTime() ){
712 <    info->thermalTime = the_globals->getThermalTime();
711 >  if( globals->haveThermalTime() ){
712 >    info->thermalTime = globals->getThermalTime();
713    }
714  
715    // check for the temperature set flag
716  
717 <  if( the_globals->haveTempSet() ) info->setTemp = the_globals->getTempSet();
717 >  if( globals->haveTempSet() ) info->setTemp = globals->getTempSet();
718  
719    // get some of the tricky things that may still be in the globals
720  
721    double boxVector[3];
722 <  if( the_globals->haveBox() ){
723 <    boxVector[0] = the_globals->getBox();
724 <    boxVector[1] = the_globals->getBox();
725 <    boxVector[2] = the_globals->getBox();
722 >  if( globals->haveBox() ){
723 >    boxVector[0] = globals->getBox();
724 >    boxVector[1] = globals->getBox();
725 >    boxVector[2] = globals->getBox();
726      
727      info->setBox( boxVector );
728    }
729 <  else if( the_globals->haveDensity() ){
729 >  else if( globals->haveDensity() ){
730  
731      double vol;
732 <    vol = (double)tot_nmol / the_globals->getDensity();
732 >    vol = (double)tot_nmol / globals->getDensity();
733       boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
734       boxVector[1] = boxVector[0];
735       boxVector[2] = boxVector[0];
# Line 925 | Line 737 | void SimSetup::gatherInfo( void ){
737      info->setBox( boxVector );
738    }
739    else{
740 <    if( !the_globals->haveBoxX() ){
740 >    if( !globals->haveBoxX() ){
741        sprintf( painCave.errMsg,
742                 "SimSetup error, no periodic BoxX size given.\n" );
743        painCave.isFatal = 1;
744        simError();
745      }
746 <    boxVector[0] = the_globals->getBoxX();
746 >    boxVector[0] = globals->getBoxX();
747  
748 <    if( !the_globals->haveBoxY() ){
748 >    if( !globals->haveBoxY() ){
749        sprintf( painCave.errMsg,
750                 "SimSetup error, no periodic BoxY size given.\n" );
751        painCave.isFatal = 1;
752        simError();
753      }
754 <    boxVector[1] = the_globals->getBoxY();
754 >    boxVector[1] = globals->getBoxY();
755  
756 <    if( !the_globals->haveBoxZ() ){
756 >    if( !globals->haveBoxZ() ){
757        sprintf( painCave.errMsg,
758                 "SimSetup error, no periodic BoxZ size given.\n" );
759        painCave.isFatal = 1;
760        simError();
761      }
762 <    boxVector[2] = the_globals->getBoxZ();
762 >    boxVector[2] = globals->getBoxZ();
763  
764      info->setBox( boxVector );
765    }
# Line 977 | Line 789 | void SimSetup::finalInfoCheck( void ){
789    }
790    
791   #ifdef IS_MPI
792 <  int myUse = usesDipoles
793 <  MPI_Allreduce( &myUse, &UsesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
792 >  int myUse = usesDipoles;
793 >  MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
794   #endif //is_mpi
795  
796 +  double theEcr, theEst;
797  
798 <  if (the_globals->getUseRF() ) {
798 >  if (globals->getUseRF() ) {
799      info->useReactionField = 1;
800      
801 <    if( !the_globals->haveECR() ){
801 >    if( !globals->haveECR() ){
802        sprintf( painCave.errMsg,
803                 "SimSetup Warning: using default value of 1/2 the smallest "
804                 "box length for the electrostaticCutoffRadius.\n"
# Line 993 | Line 806 | void SimSetup::finalInfoCheck( void ){
806        painCave.isFatal = 0;
807        simError();
808        double smallest;
809 <      smallest = info->boxLx;
810 <      if (info->boxLy <= smallest) smallest = info->boxLy;
811 <      if (info->boxLz <= smallest) smallest = info->boxLz;
812 <      info->ecr = 0.5 * smallest;
809 >      smallest = info->boxL[0];
810 >      if (info->boxL[1] <= smallest) smallest = info->boxL[1];
811 >      if (info->boxL[2] <= smallest) smallest = info->boxL[2];
812 >      theEcr = 0.5 * smallest;
813      } else {
814 <      info->ecr        = the_globals->getECR();
814 >      theEcr = globals->getECR();
815      }
816  
817 <    if( !the_globals->haveEST() ){
817 >    if( !globals->haveEST() ){
818        sprintf( painCave.errMsg,
819                 "SimSetup Warning: using default value of 0.05 * the "
820                 "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
821                 );
822        painCave.isFatal = 0;
823        simError();
824 <      info->est = 0.05 * info->ecr;
824 >      theEst = 0.05 * theEcr;
825      } else {
826 <      info->est        = the_globals->getEST();
826 >      theEst= globals->getEST();
827      }
828 +
829 +    info->setEcr( theEcr, theEst );
830      
831 <    if(!the_globals->haveDielectric() ){
831 >    if(!globals->haveDielectric() ){
832        sprintf( painCave.errMsg,
833                 "SimSetup Error: You are trying to use Reaction Field without"
834                 "setting a dielectric constant!\n"
# Line 1021 | Line 836 | void SimSetup::finalInfoCheck( void ){
836        painCave.isFatal = 1;
837        simError();
838      }
839 <    info->dielectric = the_globals->getDielectric();  
839 >    info->dielectric = globals->getDielectric();  
840    }
841    else {
842      if (usesDipoles) {
843        
844 <      if( !the_globals->haveECR() ){
845 <        sprintf( painCave.errMsg,
846 <                 "SimSetup Warning: using default value of 1/2 the smallest "
847 <                 "box length for the electrostaticCutoffRadius.\n"
848 <                 "I hope you have a very fast processor!\n");
849 <        painCave.isFatal = 0;
850 <        simError();
851 <        double smallest;
852 <        smallest = info->boxLx;
853 <        if (info->boxLy <= smallest) smallest = info->boxLy;
854 <        if (info->boxLz <= smallest) smallest = info->boxLz;
855 <        info->ecr = 0.5 * smallest;
844 >      if( !globals->haveECR() ){
845 >        sprintf( painCave.errMsg,
846 >                 "SimSetup Warning: using default value of 1/2 the smallest "
847 >                 "box length for the electrostaticCutoffRadius.\n"
848 >                 "I hope you have a very fast processor!\n");
849 >        painCave.isFatal = 0;
850 >        simError();
851 >        double smallest;
852 >        smallest = info->boxL[0];
853 >        if (info->boxL[1] <= smallest) smallest = info->boxL[1];
854 >        if (info->boxL[2] <= smallest) smallest = info->boxL[2];
855 >        theEcr = 0.5 * smallest;
856        } else {
857 <        info->ecr        = the_globals->getECR();
857 >        theEcr = globals->getECR();
858        }
859        
860 <      if( !the_globals->haveEST() ){
861 <        sprintf( painCave.errMsg,
862 <                 "SimSetup Warning: using default value of 5%% of the "
863 <                 "electrostaticCutoffRadius for the "
864 <                 "electrostaticSkinThickness\n"
865 <                 );
866 <        painCave.isFatal = 0;
867 <        simError();
868 <        info->est = 0.05 * info->ecr;
860 >      if( !globals->haveEST() ){
861 >        sprintf( painCave.errMsg,
862 >                 "SimSetup Warning: using default value of 0.05 * the "
863 >                 "electrostaticCutoffRadius for the "
864 >                 "electrostaticSkinThickness\n"
865 >                 );
866 >        painCave.isFatal = 0;
867 >        simError();
868 >        theEst = 0.05 * theEcr;
869        } else {
870 <        info->est        = the_globals->getEST();
870 >        theEst= globals->getEST();
871        }
872 +
873 +      info->setEcr( theEcr, theEst );
874      }
875    }  
876  
# Line 1066 | Line 883 | void SimSetup::initSystemCoords( void ){
883  
884   void SimSetup::initSystemCoords( void ){
885  
886 < if( the_globals->haveInitialConfig() ){
886 > if( globals->haveInitialConfig() ){
887  
888       InitializeFromFile* fileInit;
889   #ifdef IS_MPI // is_mpi
890       if( worldRank == 0 ){
891   #endif //is_mpi
892 <   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
892 >   fileInit = new InitializeFromFile( globals->getInitialConfig() );
893   #ifdef IS_MPI
894       }else fileInit = new InitializeFromFile( NULL );
895   #endif
896 <   fileInit->read_xyz( info ); // default velocities on
896 >   fileInit->readInit( info ); // default velocities on
897  
898     delete fileInit;
899   }
# Line 1113 | Line 930 | void SimSetup::makeOutNames( void ){
930    if( worldRank == 0 ){
931   #endif // is_mpi
932      
933 <    if( the_globals->haveFinalConfig() ){
934 <      strcpy( info->finalName, the_globals->getFinalConfig() );
933 >    if( globals->haveFinalConfig() ){
934 >      strcpy( info->finalName, globals->getFinalConfig() );
935      }
936      else{
937        strcpy( info->finalName, inFileName );
# Line 1197 | Line 1014 | void SimSetup::sysObjectsCreation( void ){
1014  
1015   void SimSetup::sysObjectsCreation( void ){
1016  
1017 +  int i;
1018 +
1019    // create the forceField
1020  
1021    createFF();
# Line 1219 | Line 1038 | void SimSetup::sysObjectsCreation( void ){
1038    
1039    makeSysArrays();
1040  
1041 +  // make and initialize the molecules (all but atomic coordinates)
1042 +  
1043 +  makeMolecules();
1044 +  info->identArray = new int[info->n_atoms];
1045 +  for(i=0; i<info->n_atoms; i++){
1046 +    info->identArray[i] = the_atoms[i]->getIdent();
1047 +  }
1048    
1049  
1050  
# Line 1237 | Line 1063 | void SimSetup::createFF( void ){
1063      the_ff = new LJFF();
1064      break;
1065  
1066 +  case FF_EAM:
1067 +    the_ff = new EAM_FF();
1068 +    break;
1069 +
1070    default:
1071      sprintf( painCave.errMsg,
1072               "SimSetup Error. Unrecognized force field in case statement.\n");
# Line 1253 | Line 1083 | void SimSetup::compList( void ){
1083  
1084  
1085   void SimSetup::compList( void ){
1086 +
1087 +  int i;
1088  
1089    comp_stamps = new MoleculeStamp*[n_components];
1090  
# Line 1279 | Line 1111 | void SimSetup::compList( void ){
1111        
1112        // extract the component from the list;
1113        
1114 <      currentStamp = the_stamps->extractMolStamp( id );
1114 >      currentStamp = stamps->extractMolStamp( id );
1115        if( currentStamp == NULL ){
1116          sprintf( painCave.errMsg,
1117                   "SimSetup error: Component \"%s\" was not found in the "
# Line 1303 | Line 1135 | void SimSetup::calcSysValues( void ){
1135   }
1136  
1137   void SimSetup::calcSysValues( void ){
1138 +  int i, j, k;
1139  
1140 +
1141    tot_atoms = 0;
1142    tot_bonds = 0;
1143    tot_bends = 0;
# Line 1333 | Line 1167 | void SimSetup::mpiMolDivide( void ){
1167  
1168   void SimSetup::mpiMolDivide( void ){
1169    
1170 +  int i, j, k;
1171    int localMol, allMol;
1172    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1173  
# Line 1402 | Line 1237 | void SimSetup::makeSysArrays( void ){
1237  
1238  
1239   void SimSetup::makeSysArrays( void ){
1240 +  int i, j, k;
1241  
1242 +
1243    // create the atom and short range interaction arrays
1244  
1245    Atom::createArrays(info->n_atoms);
# Line 1476 | Line 1313 | void SimSetup::makeSysArrays( void ){
1313    the_ff->setSimInfo( info );
1314  
1315   }
1316 +
1317 + void SimSetup::makeIntegrator( void ){
1318 +
1319 +  NVT<RealIntegrator>*  myNVT = NULL;
1320 +  NPTi<RealIntegrator>* myNPTi = NULL;
1321 +  NPTf<RealIntegrator>* myNPTf = NULL;
1322 +  NPTim<RealIntegrator>* myNPTim = NULL;
1323 +  NPTfm<RealIntegrator>* myNPTfm = NULL;
1324 +  ZConstraint<NVE<RealIntegrator> >* myNVEZCons = NULL;
1325 +  ZConstraint<NVT<RealIntegrator> >* myNVTZCons = NULL;
1326 +  ZConstraint<NPTi<RealIntegrator> >* myNPTiZCons = NULL;
1327 +  ZConstraint<NPTf<RealIntegrator> >* myNPTfZCons = NULL;
1328 +  ZConstraint<NPTim<RealIntegrator> >* myNPTimZCons = NULL;
1329 +  ZConstraint<NPTfm<RealIntegrator> >* myNPTfmZCons = NULL;
1330 +      
1331 +  cerr << "setting integrator" <<endl;    
1332 +  
1333 +  switch( ensembleCase ){
1334 +
1335 +  case NVE_ENS:
1336 +    new NVE<RealIntegrator>( info, the_ff );
1337 +    break;
1338 +
1339 +  case NVT_ENS:
1340 +    myNVT = new NVT<RealIntegrator>( info, the_ff );
1341 +    myNVT->setTargetTemp(globals->getTargetTemp());
1342 +
1343 +    if (globals->haveTauThermostat())
1344 +      myNVT->setTauThermostat(globals->getTauThermostat());
1345 +
1346 +    else {
1347 +      sprintf( painCave.errMsg,
1348 +               "SimSetup error: If you use the NVT\n"
1349 +               "    ensemble, you must set tauThermostat.\n");
1350 +      painCave.isFatal = 1;
1351 +      simError();
1352 +    }
1353 +    break;
1354 +
1355 +  case NPTi_ENS:
1356 +    myNPTi = new NPTi<RealIntegrator>( info, the_ff );
1357 +    myNPTi->setTargetTemp( globals->getTargetTemp() );
1358 +
1359 +    if (globals->haveTargetPressure())
1360 +      myNPTi->setTargetPressure(globals->getTargetPressure());
1361 +    else {
1362 +      sprintf( painCave.errMsg,
1363 +               "SimSetup error: If you use a constant pressure\n"
1364 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1365 +      painCave.isFatal = 1;
1366 +      simError();
1367 +    }
1368 +    
1369 +    if( globals->haveTauThermostat() )
1370 +      myNPTi->setTauThermostat( globals->getTauThermostat() );
1371 +    else{
1372 +      sprintf( painCave.errMsg,
1373 +               "SimSetup error: If you use an NPT\n"
1374 +               "    ensemble, you must set tauThermostat.\n");
1375 +      painCave.isFatal = 1;
1376 +      simError();
1377 +    }
1378 +
1379 +    if( globals->haveTauBarostat() )
1380 +      myNPTi->setTauBarostat( globals->getTauBarostat() );
1381 +    else{
1382 +      sprintf( painCave.errMsg,
1383 +               "SimSetup error: If you use an NPT\n"
1384 +               "    ensemble, you must set tauBarostat.\n");
1385 +      painCave.isFatal = 1;
1386 +      simError();
1387 +    }
1388 +    break;
1389 +
1390 +  case NPTf_ENS:
1391 +    myNPTf = new NPTf<RealIntegrator>( info, the_ff );
1392 +    myNPTf->setTargetTemp( globals->getTargetTemp());
1393 +
1394 +    if (globals->haveTargetPressure())
1395 +      myNPTf->setTargetPressure(globals->getTargetPressure());
1396 +    else {
1397 +      sprintf( painCave.errMsg,
1398 +               "SimSetup error: If you use a constant pressure\n"
1399 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1400 +      painCave.isFatal = 1;
1401 +      simError();
1402 +    }    
1403 +
1404 +    if( globals->haveTauThermostat() )
1405 +      myNPTf->setTauThermostat( globals->getTauThermostat() );
1406 +    else{
1407 +      sprintf( painCave.errMsg,
1408 +               "SimSetup error: If you use an NPT\n"
1409 +               "    ensemble, you must set tauThermostat.\n");
1410 +      painCave.isFatal = 1;
1411 +      simError();
1412 +    }
1413 +
1414 +    if( globals->haveTauBarostat() )
1415 +      myNPTf->setTauBarostat( globals->getTauBarostat() );
1416 +    else{
1417 +      sprintf( painCave.errMsg,
1418 +               "SimSetup error: If you use an NPT\n"
1419 +               "    ensemble, you must set tauBarostat.\n");
1420 +      painCave.isFatal = 1;
1421 +      simError();
1422 +    }
1423 +    break;
1424 +    
1425 +  case NPTim_ENS:
1426 +    myNPTim = new NPTim<RealIntegrator>( info, the_ff );
1427 +    myNPTim->setTargetTemp( globals->getTargetTemp());
1428 +
1429 +    if (globals->haveTargetPressure())
1430 +      myNPTim->setTargetPressure(globals->getTargetPressure());
1431 +    else {
1432 +      sprintf( painCave.errMsg,
1433 +               "SimSetup error: If you use a constant pressure\n"
1434 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1435 +      painCave.isFatal = 1;
1436 +      simError();
1437 +    }
1438 +    
1439 +    if( globals->haveTauThermostat() )
1440 +      myNPTim->setTauThermostat( globals->getTauThermostat() );
1441 +    else{
1442 +      sprintf( painCave.errMsg,
1443 +               "SimSetup error: If you use an NPT\n"
1444 +               "    ensemble, you must set tauThermostat.\n");
1445 +      painCave.isFatal = 1;
1446 +      simError();
1447 +    }
1448 +
1449 +    if( globals->haveTauBarostat() )
1450 +      myNPTim->setTauBarostat( globals->getTauBarostat() );
1451 +    else{
1452 +      sprintf( painCave.errMsg,
1453 +               "SimSetup error: If you use an NPT\n"
1454 +               "    ensemble, you must set tauBarostat.\n");
1455 +      painCave.isFatal = 1;
1456 +      simError();
1457 +    }
1458 +    break;
1459 +
1460 +  case NPTfm_ENS:
1461 +    myNPTfm = new NPTfm<RealIntegrator>( info, the_ff );
1462 +    myNPTfm->setTargetTemp( globals->getTargetTemp());
1463 +
1464 +    if (globals->haveTargetPressure())
1465 +      myNPTfm->setTargetPressure(globals->getTargetPressure());
1466 +    else {
1467 +      sprintf( painCave.errMsg,
1468 +               "SimSetup error: If you use a constant pressure\n"
1469 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1470 +      painCave.isFatal = 1;
1471 +      simError();
1472 +    }
1473 +    
1474 +    if( globals->haveTauThermostat() )
1475 +      myNPTfm->setTauThermostat( globals->getTauThermostat() );
1476 +    else{
1477 +      sprintf( painCave.errMsg,
1478 +               "SimSetup error: If you use an NPT\n"
1479 +               "    ensemble, you must set tauThermostat.\n");
1480 +      painCave.isFatal = 1;
1481 +      simError();
1482 +    }
1483 +
1484 +    if( globals->haveTauBarostat() )
1485 +      myNPTfm->setTauBarostat( globals->getTauBarostat() );
1486 +    else{
1487 +      sprintf( painCave.errMsg,
1488 +               "SimSetup error: If you use an NPT\n"
1489 +               "    ensemble, you must set tauBarostat.\n");
1490 +      painCave.isFatal = 1;
1491 +      simError();
1492 +    }
1493 +    break;
1494 +    
1495 +  case NVEZCONS_ENS:
1496 +
1497 +
1498 +    //setup index of z-constraint molecules, z-constraint sampel time
1499 +    //and z-constraint force output name. These parameter should be known
1500 +    //before constructing the z-constraint integrator
1501 +    setupZConstraint();
1502 +      
1503 +    myNVEZCons = new ZConstraint<NVE<RealIntegrator> >( info, the_ff );
1504 +        
1505 +    break;
1506 +    
1507 +    
1508 +  case NVTZCONS_ENS:
1509 +  
1510 +    setupZConstraint();
1511 +    
1512 +    myNVTZCons = new ZConstraint<NVT<RealIntegrator> >( info, the_ff );
1513 +    myNVTZCons->setTargetTemp(globals->getTargetTemp());
1514 +
1515 +    if (globals->haveTauThermostat())
1516 +      myNVTZCons->setTauThermostat(globals->getTauThermostat());
1517 +
1518 +    else {
1519 +      sprintf( painCave.errMsg,
1520 +               "SimSetup error: If you use the NVT\n"
1521 +               "    ensemble, you must set tauThermostat.\n");
1522 +      painCave.isFatal = 1;
1523 +      simError();
1524 +    }    
1525 +    break;    
1526 +    
1527 +  case NPTiZCONS_ENS:
1528 +  
1529 +    setupZConstraint();
1530 +    
1531 +    myNPTiZCons = new ZConstraint<NPTi<RealIntegrator> >( info, the_ff );
1532 +    myNPTiZCons->setTargetTemp( globals->getTargetTemp() );
1533 +
1534 +    if (globals->haveTargetPressure())
1535 +      myNPTiZCons->setTargetPressure(globals->getTargetPressure());
1536 +    else {
1537 +      sprintf( painCave.errMsg,
1538 +               "SimSetup error: If you use a constant pressure\n"
1539 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1540 +      painCave.isFatal = 1;
1541 +      simError();
1542 +    }
1543 +    
1544 +    if( globals->haveTauThermostat() )
1545 +      myNPTiZCons->setTauThermostat( globals->getTauThermostat() );
1546 +    else{
1547 +      sprintf( painCave.errMsg,
1548 +               "SimSetup error: If you use an NPT\n"
1549 +               "    ensemble, you must set tauThermostat.\n");
1550 +      painCave.isFatal = 1;
1551 +      simError();
1552 +    }
1553 +
1554 +    if( globals->haveTauBarostat() )
1555 +      myNPTiZCons->setTauBarostat( globals->getTauBarostat() );
1556 +    else{
1557 +      sprintf( painCave.errMsg,
1558 +               "SimSetup error: If you use an NPT\n"
1559 +               "    ensemble, you must set tauBarostat.\n");
1560 +      painCave.isFatal = 1;
1561 +      simError();
1562 +    }  
1563 +    
1564 +    break;
1565 +    
1566 +  case NPTfZCONS_ENS:
1567 +  
1568 +    setupZConstraint();
1569 +  
1570 +    myNPTfZCons = new ZConstraint<NPTf<RealIntegrator> >( info, the_ff );
1571 +    myNPTfZCons->setTargetTemp( globals->getTargetTemp());
1572 +
1573 +    if (globals->haveTargetPressure())
1574 +      myNPTfZCons->setTargetPressure(globals->getTargetPressure());
1575 +    else {
1576 +      sprintf( painCave.errMsg,
1577 +               "SimSetup error: If you use a constant pressure\n"
1578 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1579 +      painCave.isFatal = 1;
1580 +      simError();
1581 +    }    
1582 +
1583 +    if( globals->haveTauThermostat() )
1584 +      myNPTfZCons->setTauThermostat( globals->getTauThermostat() );
1585 +    else{
1586 +      sprintf( painCave.errMsg,
1587 +               "SimSetup error: If you use an NPT\n"
1588 +               "    ensemble, you must set tauThermostat.\n");
1589 +      painCave.isFatal = 1;
1590 +      simError();
1591 +    }
1592 +
1593 +    if( globals->haveTauBarostat() )
1594 +      myNPTfZCons->setTauBarostat( globals->getTauBarostat() );
1595 +    else{
1596 +      sprintf( painCave.errMsg,
1597 +               "SimSetup error: If you use an NPT\n"
1598 +               "    ensemble, you must set tauBarostat.\n");
1599 +      painCave.isFatal = 1;
1600 +      simError();
1601 +    }  
1602 +    
1603 +    break;  
1604 +      
1605 +  case NPTimZCONS_ENS:
1606 +  
1607 +    setupZConstraint();
1608 +  
1609 +    myNPTimZCons = new ZConstraint<NPTim<RealIntegrator> >( info, the_ff );
1610 +    myNPTimZCons->setTargetTemp( globals->getTargetTemp());
1611 +
1612 +    if (globals->haveTargetPressure())
1613 +      myNPTimZCons->setTargetPressure(globals->getTargetPressure());
1614 +    else {
1615 +      sprintf( painCave.errMsg,
1616 +               "SimSetup error: If you use a constant pressure\n"
1617 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1618 +      painCave.isFatal = 1;
1619 +      simError();
1620 +    }
1621 +    
1622 +    if( globals->haveTauThermostat() )
1623 +      myNPTimZCons->setTauThermostat( globals->getTauThermostat() );
1624 +    else{
1625 +      sprintf( painCave.errMsg,
1626 +               "SimSetup error: If you use an NPT\n"
1627 +               "    ensemble, you must set tauThermostat.\n");
1628 +      painCave.isFatal = 1;
1629 +      simError();
1630 +    }
1631 +
1632 +    if( globals->haveTauBarostat() )
1633 +      myNPTimZCons->setTauBarostat( globals->getTauBarostat() );
1634 +    else{
1635 +      sprintf( painCave.errMsg,
1636 +               "SimSetup error: If you use an NPT\n"
1637 +               "    ensemble, you must set tauBarostat.\n");
1638 +      painCave.isFatal = 1;
1639 +      simError();
1640 +    }  
1641 +    
1642 +    break;
1643 +    
1644 +  case NPTfmZCONS_ENS:
1645 +  
1646 +    setupZConstraint();
1647 +    
1648 +    myNPTfmZCons = new ZConstraint<NPTfm<RealIntegrator> >( info, the_ff );
1649 +    myNPTfmZCons->setTargetTemp( globals->getTargetTemp());
1650 +
1651 +    if (globals->haveTargetPressure())
1652 +      myNPTfmZCons->setTargetPressure(globals->getTargetPressure());
1653 +    else {
1654 +      sprintf( painCave.errMsg,
1655 +               "SimSetup error: If you use a constant pressure\n"
1656 +               "    ensemble, you must set targetPressure in the BASS file.\n");
1657 +      painCave.isFatal = 1;
1658 +      simError();
1659 +    }
1660 +    
1661 +    if( globals->haveTauThermostat() )
1662 +      myNPTfmZCons->setTauThermostat( globals->getTauThermostat() );
1663 +    else{
1664 +      sprintf( painCave.errMsg,
1665 +               "SimSetup error: If you use an NPT\n"
1666 +               "    ensemble, you must set tauThermostat.\n");
1667 +      painCave.isFatal = 1;
1668 +      simError();
1669 +    }
1670 +
1671 +    if( globals->haveTauBarostat() )
1672 +      myNPTfmZCons->setTauBarostat( globals->getTauBarostat() );
1673 +    else{
1674 +      sprintf( painCave.errMsg,
1675 +               "SimSetup error: If you use an NPT\n"
1676 +               "    ensemble, you must set tauBarostat.\n");
1677 +      painCave.isFatal = 1;
1678 +      simError();
1679 +    }    
1680 +    break;      
1681 +      
1682 +  
1683 +    
1684 +  default:
1685 +    sprintf( painCave.errMsg,
1686 +             "SimSetup Error. Unrecognized ensemble in case statement.\n");
1687 +    painCave.isFatal = 1;
1688 +    simError();
1689 +  }
1690 +
1691 + }
1692 +
1693 + void SimSetup::initFortran( void ){
1694 +
1695 +  info->refreshSim();
1696 +  
1697 +  if( !strcmp( info->mixingRule, "standard") ){
1698 +    the_ff->initForceField( LB_MIXING_RULE );
1699 +  }
1700 +  else if( !strcmp( info->mixingRule, "explicit") ){
1701 +    the_ff->initForceField( EXPLICIT_MIXING_RULE );
1702 +  }
1703 +  else{
1704 +    sprintf( painCave.errMsg,
1705 +             "SimSetup Error: unknown mixing rule -> \"%s\"\n",
1706 +             info->mixingRule );
1707 +    painCave.isFatal = 1;
1708 +    simError();
1709 +  }
1710 +
1711 +
1712 + #ifdef IS_MPI
1713 +  strcpy( checkPointMsg,
1714 +          "Successfully intialized the mixingRule for Fortran." );
1715 +  MPIcheckPoint();
1716 + #endif // is_mpi
1717 +
1718 + }
1719 +
1720 + void SimSetup::setupZConstraint()
1721 + {
1722 +  if(globals->haveZConsTime()){  
1723 +
1724 +  //add sample time of z-constraint  into SimInfo's property list                    
1725 +  DoubleData* zconsTimeProp = new DoubleData();
1726 +  zconsTimeProp->setID("zconstime");
1727 +  zconsTimeProp->setData(globals->getZConsTime());
1728 +  info->addProperty(zconsTimeProp);
1729 +  }
1730 +  else{
1731 +    sprintf( painCave.errMsg,
1732 +             "ZConstraint error: If you use an ZConstraint\n"
1733 +             " , you must set sample time.\n");
1734 +    painCave.isFatal = 1;
1735 +    simError();      
1736 +  }
1737 +      
1738 +  if(globals->haveIndexOfAllZConsMols()){
1739 +
1740 +        //add index of z-constraint molecules into SimInfo's property list
1741 +        vector<int> tempIndex = globals->getIndexOfAllZConsMols();
1742 +        
1743 +        //sort the index
1744 +        sort(tempIndex.begin(), tempIndex.end());
1745 +        
1746 +        IndexData* zconsIndex = new IndexData();
1747 +        zconsIndex->setID("zconsindex");
1748 +        zconsIndex->setIndexData(tempIndex);
1749 +        info->addProperty(zconsIndex);
1750 +  }
1751 +  else{
1752 +    sprintf( painCave.errMsg,
1753 +             "SimSetup error: If you use an ZConstraint\n"
1754 +             " , you must set index of z-constraint molecules.\n");
1755 +    painCave.isFatal = 1;
1756 +    simError();    
1757 +      
1758 +  }
1759 +
1760 +  //Determine the name of ouput file and add it into SimInfo's property list
1761 +  //Be careful, do not use inFileName, since it is a pointer which
1762 +  //point to a string at master node, and slave nodes do not contain that string
1763 +    
1764 +  string zconsOutput(info->finalName);
1765 +            
1766 +  zconsOutput = zconsOutput.substr(0, zconsOutput.rfind(".")) + ".fz";
1767 +                
1768 +  StringData* zconsFilename = new StringData();
1769 +  zconsFilename->setID("zconsfilename");
1770 +  zconsFilename->setData(zconsOutput);
1771 +
1772 +  info->addProperty(zconsFilename);      
1773 +
1774 + }

Diff Legend

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