[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 656 by mmeineke, Tue Jul 29 16:32:37 2003 UTC

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

Diff Legend

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

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

Diff Legend

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