[ViewVC] Diff of: group/trunk/OOPSE/libmdtools/SimSetup.cpp

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 614 by mmeineke, Tue Jul 15 17:57:04 2003 UTC vs.
Revision 653 by chuckv, Fri Jul 25 20:00:17 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		stamps = new MakeStamps();
31		globals = new Globals();
#	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;
589	+	else if( !strcasecmp( force_field, "EAM" )) ffCase = FF_EAM;
590		else{
591		sprintf( painCave.errMsg,
592		"SimSetup Error. Unrecognized force field -> %s\n",
#	Line 815 \| Line 597 \| void SimSetup::gatherInfo( void ){
597
598		// get the ensemble
599
600	<	strcpy( ensemble, the_globals->getEnsemble() );
600	>	strcpy( ensemble, globals->getEnsemble() );
601
602		if( !strcasecmp( ensemble, "NVE" )) ensembleCase = NVE_ENS;
603		else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
#	Line 838 \| Line 620 \| void SimSetup::gatherInfo( void ){
620
621		// get the mixing rule
622
623	<	strcpy( info->mixingRule, the_globals->getMixingRule() );
624	<	info->usePBC = the_globals->getPBC();
623	>	strcpy( info->mixingRule, globals->getMixingRule() );
624	>	info->usePBC = globals->getPBC();
625
626
627		// get the components and calculate the tot_nMol and indvidual n_mol
628
629	<	the_components = the_globals->getComponents();
629	>	the_components = globals->getComponents();
630		components_nmol = new int[n_components];
631
632
633	<	if( !the_globals->haveNMol() ){
633	>	if( !globals->haveNMol() ){
634		// we don't have the total number of molecules, so we assume it is
635		// given in each component
636
#	Line 881 \| Line 663 \| void SimSetup::gatherInfo( void ){
663
664		// set the status, sample, and thermal kick times
665
666	<	if( the_globals->haveSampleTime() ){
667	<	info->sampleTime = the_globals->getSampleTime();
666	>	if( globals->haveSampleTime() ){
667	>	info->sampleTime = globals->getSampleTime();
668		info->statusTime = info->sampleTime;
669		info->thermalTime = info->sampleTime;
670		}
671		else{
672	<	info->sampleTime = the_globals->getRunTime();
672	>	info->sampleTime = globals->getRunTime();
673		info->statusTime = info->sampleTime;
674		info->thermalTime = info->sampleTime;
675		}
676
677	<	if( the_globals->haveStatusTime() ){
678	<	info->statusTime = the_globals->getStatusTime();
677	>	if( globals->haveStatusTime() ){
678	>	info->statusTime = globals->getStatusTime();
679		}
680
681	<	if( the_globals->haveThermalTime() ){
682	<	info->thermalTime = the_globals->getThermalTime();
681	>	if( globals->haveThermalTime() ){
682	>	info->thermalTime = globals->getThermalTime();
683		}
684
685		// check for the temperature set flag
686
687	<	if( the_globals->haveTempSet() ) info->setTemp = the_globals->getTempSet();
687	>	if( globals->haveTempSet() ) info->setTemp = globals->getTempSet();
688
689		// get some of the tricky things that may still be in the globals
690
691		double boxVector[3];
692	<	if( the_globals->haveBox() ){
693	<	boxVector[0] = the_globals->getBox();
694	<	boxVector[1] = the_globals->getBox();
695	<	boxVector[2] = the_globals->getBox();
692	>	if( globals->haveBox() ){
693	>	boxVector[0] = globals->getBox();
694	>	boxVector[1] = globals->getBox();
695	>	boxVector[2] = globals->getBox();
696
697		info->setBox( boxVector );
698		}
699	<	else if( the_globals->haveDensity() ){
699	>	else if( globals->haveDensity() ){
700
701		double vol;
702	<	vol = (double)tot_nmol / the_globals->getDensity();
702	>	vol = (double)tot_nmol / globals->getDensity();
703		boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
704		boxVector[1] = boxVector[0];
705		boxVector[2] = boxVector[0];
#	Line 925 \| Line 707 \| void SimSetup::gatherInfo( void ){
707		info->setBox( boxVector );
708		}
709		else{
710	<	if( !the_globals->haveBoxX() ){
710	>	if( !globals->haveBoxX() ){
711		sprintf( painCave.errMsg,
712		"SimSetup error, no periodic BoxX size given.\n" );
713		painCave.isFatal = 1;
714		simError();
715		}
716	<	boxVector[0] = the_globals->getBoxX();
716	>	boxVector[0] = globals->getBoxX();
717
718	<	if( !the_globals->haveBoxY() ){
718	>	if( !globals->haveBoxY() ){
719		sprintf( painCave.errMsg,
720		"SimSetup error, no periodic BoxY size given.\n" );
721		painCave.isFatal = 1;
722		simError();
723		}
724	<	boxVector[1] = the_globals->getBoxY();
724	>	boxVector[1] = globals->getBoxY();
725
726	<	if( !the_globals->haveBoxZ() ){
726	>	if( !globals->haveBoxZ() ){
727		sprintf( painCave.errMsg,
728		"SimSetup error, no periodic BoxZ size given.\n" );
729		painCave.isFatal = 1;
730		simError();
731		}
732	<	boxVector[2] = the_globals->getBoxZ();
732	>	boxVector[2] = globals->getBoxZ();
733
734		info->setBox( boxVector );
735		}
#	Line 977 \| Line 759 \| void SimSetup::finalInfoCheck( void ){
759		}
760
761		#ifdef IS_MPI
762	<	int myUse = usesDipoles
763	<	MPI_Allreduce( &myUse, &UsesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
762	>	int myUse = usesDipoles;
763	>	MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
764		#endif //is_mpi
765
766	+	double theEcr, theEst;
767
768	<	if (the_globals->getUseRF() ) {
768	>	if (globals->getUseRF() ) {
769		info->useReactionField = 1;
770
771	<	if( !the_globals->haveECR() ){
771	>	if( !globals->haveECR() ){
772		sprintf( painCave.errMsg,
773		"SimSetup Warning: using default value of 1/2 the smallest "
774		"box length for the electrostaticCutoffRadius.\n"
#	Line 993 \| Line 776 \| void SimSetup::finalInfoCheck( void ){
776		painCave.isFatal = 0;
777		simError();
778		double smallest;
779	<	smallest = info->boxLx;
780	<	if (info->boxLy <= smallest) smallest = info->boxLy;
781	<	if (info->boxLz <= smallest) smallest = info->boxLz;
782	<	info->ecr = 0.5 * smallest;
779	>	smallest = info->boxL[0];
780	>	if (info->boxL[1] <= smallest) smallest = info->boxL[1];
781	>	if (info->boxL[2] <= smallest) smallest = info->boxL[2];
782	>	theEcr = 0.5 * smallest;
783		} else {
784	<	info->ecr = the_globals->getECR();
784	>	theEcr = globals->getECR();
785		}
786
787	<	if( !the_globals->haveEST() ){
787	>	if( !globals->haveEST() ){
788		sprintf( painCave.errMsg,
789		"SimSetup Warning: using default value of 0.05 * the "
790		"electrostaticCutoffRadius for the electrostaticSkinThickness\n"
791		);
792		painCave.isFatal = 0;
793		simError();
794	<	info->est = 0.05 * info->ecr;
794	>	theEst = 0.05 * theEcr;
795		} else {
796	<	info->est = the_globals->getEST();
796	>	theEst= globals->getEST();
797		}
798	+
799	+	info->setEcr( theEcr, theEst );
800
801	<	if(!the_globals->haveDielectric() ){
801	>	if(!globals->haveDielectric() ){
802		sprintf( painCave.errMsg,
803		"SimSetup Error: You are trying to use Reaction Field without"
804		"setting a dielectric constant!\n"
#	Line 1021 \| Line 806 \| void SimSetup::finalInfoCheck( void ){
806		painCave.isFatal = 1;
807		simError();
808		}
809	<	info->dielectric = the_globals->getDielectric();
809	>	info->dielectric = globals->getDielectric();
810		}
811		else {
812		if (usesDipoles) {
813
814	<	if( !the_globals->haveECR() ){
815	<	sprintf( painCave.errMsg,
816	<	"SimSetup Warning: using default value of 1/2 the smallest "
817	<	"box length for the electrostaticCutoffRadius.\n"
818	<	"I hope you have a very fast processor!\n");
819	<	painCave.isFatal = 0;
820	<	simError();
821	<	double smallest;
822	<	smallest = info->boxLx;
823	<	if (info->boxLy <= smallest) smallest = info->boxLy;
824	<	if (info->boxLz <= smallest) smallest = info->boxLz;
825	<	info->ecr = 0.5 * smallest;
814	>	if( !globals->haveECR() ){
815	>	sprintf( painCave.errMsg,
816	>	"SimSetup Warning: using default value of 1/2 the smallest "
817	>	"box length for the electrostaticCutoffRadius.\n"
818	>	"I hope you have a very fast processor!\n");
819	>	painCave.isFatal = 0;
820	>	simError();
821	>	double smallest;
822	>	smallest = info->boxL[0];
823	>	if (info->boxL[1] <= smallest) smallest = info->boxL[1];
824	>	if (info->boxL[2] <= smallest) smallest = info->boxL[2];
825	>	theEcr = 0.5 * smallest;
826		} else {
827	<	info->ecr = the_globals->getECR();
827	>	theEcr = globals->getECR();
828		}
829
830	<	if( !the_globals->haveEST() ){
831	<	sprintf( painCave.errMsg,
832	<	"SimSetup Warning: using default value of 5%% of the "
833	<	"electrostaticCutoffRadius for the "
834	<	"electrostaticSkinThickness\n"
835	<	);
836	<	painCave.isFatal = 0;
837	<	simError();
838	<	info->est = 0.05 * info->ecr;
839	<	} else {
840	<	info->est = the_globals->getEST();
830	>	if( !globals->haveEST() ){
831	>	sprintf( painCave.errMsg,
832	>	"SimSetup Warning: using default value of 0.05 * the "
833	>	"electrostaticCutoffRadius for the "
834	>	"electrostaticSkinThickness\n"
835	>	);
836	>	painCave.isFatal = 0;
837	>	simError();
838	>	theEst = 0.05 * theEcr;
839	>	} else {
840	>	theEst= globals->getEST();
841		}
842	+
843	+	info->setEcr( theEcr, theEst );
844		}
845		}
846
#	Line 1066 \| Line 853 \| void SimSetup::initSystemCoords( void ){
853
854		void SimSetup::initSystemCoords( void ){
855
856	<	if( the_globals->haveInitialConfig() ){
856	>	if( globals->haveInitialConfig() ){
857
858		InitializeFromFile* fileInit;
859		#ifdef IS_MPI // is_mpi
860		if( worldRank == 0 ){
861		#endif //is_mpi
862	<	fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
862	>	fileInit = new InitializeFromFile( globals->getInitialConfig() );
863		#ifdef IS_MPI
864		}else fileInit = new InitializeFromFile( NULL );
865		#endif
866	<	fileInit->read_xyz( info ); // default velocities on
866	>	fileInit->readInit( info ); // default velocities on
867
868		delete fileInit;
869		}
#	Line 1113 \| Line 900 \| void SimSetup::makeOutNames( void ){
900		if( worldRank == 0 ){
901		#endif // is_mpi
902
903	<	if( the_globals->haveFinalConfig() ){
904	<	strcpy( info->finalName, the_globals->getFinalConfig() );
903	>	if( globals->haveFinalConfig() ){
904	>	strcpy( info->finalName, globals->getFinalConfig() );
905		}
906		else{
907		strcpy( info->finalName, inFileName );
#	Line 1197 \| Line 984 \| void SimSetup::sysObjectsCreation( void ){
984
985		void SimSetup::sysObjectsCreation( void ){
986
987	+	int i;
988	+
989		// create the forceField
990
991		createFF();
#	Line 1219 \| Line 1008 \| void SimSetup::sysObjectsCreation( void ){
1008
1009		makeSysArrays();
1010
1011	+	// make and initialize the molecules (all but atomic coordinates)
1012
1013	+	makeMolecules();
1014	+	info->identArray = new int[info->n_atoms];
1015	+	for(i=0; i<info->n_atoms; i++){
1016	+	info->identArray[i] = the_atoms[i]->getIdent();
1017	+	}
1018	+
1019
1020
1021		}
#	Line 1237 \| Line 1033 \| void SimSetup::createFF( void ){
1033		the_ff = new LJFF();
1034		break;
1035
1036	+	case FF_EAM:
1037	+	the_ff = new EAM_FF();
1038	+	break;
1039	+
1040		default:
1041		sprintf( painCave.errMsg,
1042		"SimSetup Error. Unrecognized force field in case statement.\n");
#	Line 1254 \| Line 1054 \| void SimSetup::compList( void ){
1054
1055		void SimSetup::compList( void ){
1056
1057	+	int i;
1058	+
1059		comp_stamps = new MoleculeStamp*[n_components];
1060
1061		// make an array of molecule stamps that match the components used.
#	Line 1279 \| Line 1081 \| void SimSetup::compList( void ){
1081
1082		// extract the component from the list;
1083
1084	<	currentStamp = the_stamps->extractMolStamp( id );
1084	>	currentStamp = stamps->extractMolStamp( id );
1085		if( currentStamp == NULL ){
1086		sprintf( painCave.errMsg,
1087		"SimSetup error: Component \"%s\" was not found in the "
#	Line 1303 \| Line 1105 \| void SimSetup::calcSysValues( void ){
1105		}
1106
1107		void SimSetup::calcSysValues( void ){
1108	+	int i, j, k;
1109
1110	+
1111		tot_atoms = 0;
1112		tot_bonds = 0;
1113		tot_bends = 0;
#	Line 1333 \| Line 1137 \| void SimSetup::mpiMolDivide( void ){
1137
1138		void SimSetup::mpiMolDivide( void ){
1139
1140	+	int i, j, k;
1141		int localMol, allMol;
1142		int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1143
#	Line 1402 \| Line 1207 \| void SimSetup::makeSysArrays( void ){
1207
1208
1209		void SimSetup::makeSysArrays( void ){
1210	+	int i, j, k;
1211
1212	+
1213		// create the atom and short range interaction arrays
1214
1215		Atom::createArrays(info->n_atoms);
#	Line 1474 \| Line 1281 \| void SimSetup::makeSysArrays( void ){
1281		info->excludes = the_excludes;
1282
1283		the_ff->setSimInfo( info );
1284	+
1285	+	}
1286	+
1287	+	void SimSetup::makeIntegrator( void ){
1288	+
1289	+	NVT<RealIntegrator>* myNVT = NULL;
1290	+	NPTi<RealIntegrator>* myNPTi = NULL;
1291	+	NPTf<RealIntegrator>* myNPTf = NULL;
1292	+	NPTim<RealIntegrator>* myNPTim = NULL;
1293	+	NPTfm<RealIntegrator>* myNPTfm = NULL;
1294	+
1295	+	switch( ensembleCase ){
1296	+
1297	+	case NVE_ENS:
1298	+	new NVE<RealIntegrator>( info, the_ff );
1299	+	break;
1300	+
1301	+	case NVT_ENS:
1302	+	myNVT = new NVT<RealIntegrator>( info, the_ff );
1303	+	myNVT->setTargetTemp(globals->getTargetTemp());
1304	+
1305	+	if (globals->haveTauThermostat())
1306	+	myNVT->setTauThermostat(globals->getTauThermostat());
1307	+
1308	+	else {
1309	+	sprintf( painCave.errMsg,
1310	+	"SimSetup error: If you use the NVT\n"
1311	+	" ensemble, you must set tauThermostat.\n");
1312	+	painCave.isFatal = 1;
1313	+	simError();
1314	+	}
1315	+	break;
1316	+
1317	+	case NPTi_ENS:
1318	+	myNPTi = new NPTi<RealIntegrator>( info, the_ff );
1319	+	myNPTi->setTargetTemp( globals->getTargetTemp() );
1320	+
1321	+	if (globals->haveTargetPressure())
1322	+	myNPTi->setTargetPressure(globals->getTargetPressure());
1323	+	else {
1324	+	sprintf( painCave.errMsg,
1325	+	"SimSetup error: If you use a constant pressure\n"
1326	+	" ensemble, you must set targetPressure in the BASS file.\n");
1327	+	painCave.isFatal = 1;
1328	+	simError();
1329	+	}
1330	+
1331	+	if( globals->haveTauThermostat() )
1332	+	myNPTi->setTauThermostat( globals->getTauThermostat() );
1333	+	else{
1334	+	sprintf( painCave.errMsg,
1335	+	"SimSetup error: If you use an NPT\n"
1336	+	" ensemble, you must set tauThermostat.\n");
1337	+	painCave.isFatal = 1;
1338	+	simError();
1339	+	}
1340
1341	+	if( globals->haveTauBarostat() )
1342	+	myNPTi->setTauBarostat( globals->getTauBarostat() );
1343	+	else{
1344	+	sprintf( painCave.errMsg,
1345	+	"SimSetup error: If you use an NPT\n"
1346	+	" ensemble, you must set tauBarostat.\n");
1347	+	painCave.isFatal = 1;
1348	+	simError();
1349	+	}
1350	+	break;
1351	+
1352	+	case NPTf_ENS:
1353	+	myNPTf = new NPTf<RealIntegrator>( info, the_ff );
1354	+	myNPTf->setTargetTemp( globals->getTargetTemp());
1355	+
1356	+	if (globals->haveTargetPressure())
1357	+	myNPTf->setTargetPressure(globals->getTargetPressure());
1358	+	else {
1359	+	sprintf( painCave.errMsg,
1360	+	"SimSetup error: If you use a constant pressure\n"
1361	+	" ensemble, you must set targetPressure in the BASS file.\n");
1362	+	painCave.isFatal = 1;
1363	+	simError();
1364	+	}
1365	+
1366	+	if( globals->haveTauThermostat() )
1367	+	myNPTf->setTauThermostat( globals->getTauThermostat() );
1368	+	else{
1369	+	sprintf( painCave.errMsg,
1370	+	"SimSetup error: If you use an NPT\n"
1371	+	" ensemble, you must set tauThermostat.\n");
1372	+	painCave.isFatal = 1;
1373	+	simError();
1374	+	}
1375	+
1376	+	if( globals->haveTauBarostat() )
1377	+	myNPTf->setTauBarostat( globals->getTauBarostat() );
1378	+	else{
1379	+	sprintf( painCave.errMsg,
1380	+	"SimSetup error: If you use an NPT\n"
1381	+	" ensemble, you must set tauBarostat.\n");
1382	+	painCave.isFatal = 1;
1383	+	simError();
1384	+	}
1385	+	break;
1386	+
1387	+	case NPTim_ENS:
1388	+	myNPTim = new NPTim<RealIntegrator>( info, the_ff );
1389	+	myNPTim->setTargetTemp( globals->getTargetTemp());
1390	+
1391	+	if (globals->haveTargetPressure())
1392	+	myNPTim->setTargetPressure(globals->getTargetPressure());
1393	+	else {
1394	+	sprintf( painCave.errMsg,
1395	+	"SimSetup error: If you use a constant pressure\n"
1396	+	" ensemble, you must set targetPressure in the BASS file.\n");
1397	+	painCave.isFatal = 1;
1398	+	simError();
1399	+	}
1400	+
1401	+	if( globals->haveTauThermostat() )
1402	+	myNPTim->setTauThermostat( globals->getTauThermostat() );
1403	+	else{
1404	+	sprintf( painCave.errMsg,
1405	+	"SimSetup error: If you use an NPT\n"
1406	+	" ensemble, you must set tauThermostat.\n");
1407	+	painCave.isFatal = 1;
1408	+	simError();
1409	+	}
1410	+
1411	+	if( globals->haveTauBarostat() )
1412	+	myNPTim->setTauBarostat( globals->getTauBarostat() );
1413	+	else{
1414	+	sprintf( painCave.errMsg,
1415	+	"SimSetup error: If you use an NPT\n"
1416	+	" ensemble, you must set tauBarostat.\n");
1417	+	painCave.isFatal = 1;
1418	+	simError();
1419	+	}
1420	+	break;
1421	+
1422	+	case NPTfm_ENS:
1423	+	myNPTfm = new NPTfm<RealIntegrator>( info, the_ff );
1424	+	myNPTfm->setTargetTemp( globals->getTargetTemp());
1425	+
1426	+	if (globals->haveTargetPressure())
1427	+	myNPTfm->setTargetPressure(globals->getTargetPressure());
1428	+	else {
1429	+	sprintf( painCave.errMsg,
1430	+	"SimSetup error: If you use a constant pressure\n"
1431	+	" ensemble, you must set targetPressure in the BASS file.\n");
1432	+	painCave.isFatal = 1;
1433	+	simError();
1434	+	}
1435	+
1436	+	if( globals->haveTauThermostat() )
1437	+	myNPTfm->setTauThermostat( globals->getTauThermostat() );
1438	+	else{
1439	+	sprintf( painCave.errMsg,
1440	+	"SimSetup error: If you use an NPT\n"
1441	+	" ensemble, you must set tauThermostat.\n");
1442	+	painCave.isFatal = 1;
1443	+	simError();
1444	+	}
1445	+
1446	+	if( globals->haveTauBarostat() )
1447	+	myNPTfm->setTauBarostat( globals->getTauBarostat() );
1448	+	else{
1449	+	sprintf( painCave.errMsg,
1450	+	"SimSetup error: If you use an NPT\n"
1451	+	" ensemble, you must set tauBarostat.\n");
1452	+	painCave.isFatal = 1;
1453	+	simError();
1454	+	}
1455	+	break;
1456	+
1457	+	default:
1458	+	sprintf( painCave.errMsg,
1459	+	"SimSetup Error. Unrecognized ensemble in case statement.\n");
1460	+	painCave.isFatal = 1;
1461	+	simError();
1462	+	}
1463	+
1464		}
1465	+
1466	+	void SimSetup::initFortran( void ){
1467	+
1468	+	info->refreshSim();
1469	+
1470	+	if( !strcmp( info->mixingRule, "standard") ){
1471	+	the_ff->initForceField( LB_MIXING_RULE );
1472	+	}
1473	+	else if( !strcmp( info->mixingRule, "explicit") ){
1474	+	the_ff->initForceField( EXPLICIT_MIXING_RULE );
1475	+	}
1476	+	else{
1477	+	sprintf( painCave.errMsg,
1478	+	"SimSetup Error: unknown mixing rule -> \"%s\"\n",
1479	+	info->mixingRule );
1480	+	painCave.isFatal = 1;
1481	+	simError();
1482	+	}
1483	+
1484	+
1485	+	#ifdef IS_MPI
1486	+	strcpy( checkPointMsg,
1487	+	"Successfully intialized the mixingRule for Fortran." );
1488	+	MPIcheckPoint();
1489	+	#endif // is_mpi
1490	+
1491	+	}

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents): Revision 614 by mmeineke, Tue Jul 15 17:57:04 2003 UTC vs. Revision 653 by chuckv, Fri Jul 25 20:00:17 2003 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 614 by mmeineke, Tue Jul 15 17:57:04 2003 UTC vs.
Revision 653 by chuckv, Fri Jul 25 20:00:17 2003 UTC