[ViewVC] Diff of: group/branches/no-template-branch/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 658 by tim, Thu Jul 31 15:35:07 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 20 \| Line 22
22		#define NPTf_ENS 3
23		#define NPTim_ENS 4
24		#define NPTfm_ENS 5
25	+	#define NVEZCONS_ENS 6
26
27
28		#define FF_DUFF 0
29		#define FF_LJ 1
30	+	#define FF_EAM 2
31
32	+	using namespace std;
33
34		SimSetup::SimSetup(){
35	+
36	+	isInfoArray = 0;
37	+	nInfo = 1;
38	+
39		stamps = new MakeStamps();
40		globals = new Globals();
41
42	+
43		#ifdef IS_MPI
44		strcpy( checkPointMsg, "SimSetup creation successful" );
45		MPIcheckPoint();
#	Line 41 \| Line 51 \| SimSetup::~SimSetup(){
51		delete globals;
52		}
53
54	+	void SimSetup::setSimInfo( SimInfo* the_info, int theNinfo ) {
55	+	info = the_info;
56	+	nInfo = theNinfo;
57	+	isInfoArray = 1;
58	+	}
59	+
60	+
61		void SimSetup::parseFile( char* fileName ){
62
63		#ifdef IS_MPI
#	Line 78 \| Line 95 \| void SimSetup::createSim( void ){
95
96		void SimSetup::createSim( void ){
97
81	–	MakeStamps *the_stamps;
82	–	Globals* the_globals;
98		int i, j, k, globalAtomIndex;
99
100		// gather all of the information from the Bass file
#	Line 90 \| Line 105 \| void SimSetup::createSim( void ){
105
106		sysObjectsCreation();
107
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	–
108		// check on the post processing info
109
110		finalInfoCheck();
111
110	–
111	–
112	–
112		// initialize the system coordinates
113
114		initSystemCoords();
#	Line 119 \| Line 118 \| void SimSetup::createSim( void ){
118
119		makeOutNames();
120
122	–
123	–
124	–
125	–
126	–
127	–
128	–
129	–
121		// make the integrator
122
123	+	makeIntegrator();
124
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	–
125		#ifdef IS_MPI
126		mpiSim->mpiRefresh();
127		#endif
128
129		// initialize the Fortran
130
131	+	initFortran();
132
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	–	}
133
134
333	–	#ifdef IS_MPI
334	–	strcpy( checkPointMsg,
335	–	"Successfully intialized the mixingRule for Fortran." );
336	–	MPIcheckPoint();
337	–	#endif // is_mpi
135		}
136
137
138		void SimSetup::makeMolecules( void ){
139
140		int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
141	<	molInit info;
141	>	molInit molInfo;
142		DirectionalAtom* dAtom;
143		LinkedAssign* extras;
144		LinkedAssign* current_extra;
#	Line 370 \| Line 167 \| void SimSetup::makeMolecules( void ){
167
168		stampID = the_molecules[i].getStampID();
169
170	<	info.nAtoms = comp_stamps[stampID]->getNAtoms();
171	<	info.nBonds = comp_stamps[stampID]->getNBonds();
172	<	info.nBends = comp_stamps[stampID]->getNBends();
173	<	info.nTorsions = comp_stamps[stampID]->getNTorsions();
174	<	info.nExcludes = info.nBonds + info.nBends + info.nTorsions;
170	>	molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
171	>	molInfo.nBonds = comp_stamps[stampID]->getNBonds();
172	>	molInfo.nBends = comp_stamps[stampID]->getNBends();
173	>	molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
174	>	molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
175
176	<	info.myAtoms = &the_atoms[atomOffset];
177	<	info.myExcludes = &the_excludes[excludeOffset];
178	<	info.myBonds = new Bond*[info.nBonds];
179	<	info.myBends = new Bend*[info.nBends];
180	<	info.myTorsions = new Torsion*[info.nTorsions];
176	>	molInfo.myAtoms = &the_atoms[atomOffset];
177	>	molInfo.myExcludes = &the_excludes[excludeOffset];
178	>	molInfo.myBonds = new Bond*[molInfo.nBonds];
179	>	molInfo.myBends = new Bend*[molInfo.nBends];
180	>	molInfo.myTorsions = new Torsion*[molInfo.nTorsions];
181
182	<	theBonds = new bond_pair[info.nBonds];
183	<	theBends = new bend_set[info.nBends];
184	<	theTorsions = new torsion_set[info.nTorsions];
182	>	theBonds = new bond_pair[molInfo.nBonds];
183	>	theBends = new bend_set[molInfo.nBends];
184	>	theTorsions = new torsion_set[molInfo.nTorsions];
185
186		// make the Atoms
187
188	<	for(j=0; j<info.nAtoms; j++){
188	>	for(j=0; j<molInfo.nAtoms; j++){
189
190		currentAtom = comp_stamps[stampID]->getAtom( j );
191		if( currentAtom->haveOrientation() ){
192
193		dAtom = new DirectionalAtom(j + atomOffset);
194		info->n_oriented++;
195	<	info.myAtoms[j] = dAtom;
195	>	molInfo.myAtoms[j] = dAtom;
196
197		ux = currentAtom->getOrntX();
198		uy = currentAtom->getOrntY();
#	Line 413 \| Line 210 \| void SimSetup::makeMolecules( void ){
210		dAtom->setSUz( uz );
211		}
212		else{
213	<	info.myAtoms[j] = new GeneralAtom(j + atomOffset);
213	>	molInfo.myAtoms[j] = new GeneralAtom(j + atomOffset);
214		}
215	<	info.myAtoms[j]->setType( currentAtom->getType() );
215	>	molInfo.myAtoms[j]->setType( currentAtom->getType() );
216
217		#ifdef IS_MPI
218
219	<	info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
219	>	molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
220
221		#endif // is_mpi
222		}
223
224		// make the bonds
225	<	for(j=0; j<info.nBonds; j++){
225	>	for(j=0; j<molInfo.nBonds; j++){
226
227		currentBond = comp_stamps[stampID]->getBond( j );
228		theBonds[j].a = currentBond->getA() + atomOffset;
#	Line 452 \| Line 249 \| void SimSetup::makeMolecules( void ){
249		the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
250		#endif //is_mpi
251		}
252	<	excludeOffset += info.nBonds;
252	>	excludeOffset += molInfo.nBonds;
253
254		//make the bends
255	<	for(j=0; j<info.nBends; j++){
255	>	for(j=0; j<molInfo.nBends; j++){
256
257		currentBend = comp_stamps[stampID]->getBend( j );
258		theBends[j].a = currentBend->getA() + atomOffset;
#	Line 538 \| Line 335 \| void SimSetup::makeMolecules( void ){
335		the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
336		#endif //is_mpi
337		}
338	<	excludeOffset += info.nBends;
338	>	excludeOffset += molInfo.nBends;
339
340	<	for(j=0; j<info.nTorsions; j++){
340	>	for(j=0; j<molInfo.nTorsions; j++){
341
342		currentTorsion = comp_stamps[stampID]->getTorsion( j );
343		theTorsions[j].a = currentTorsion->getA() + atomOffset;
#	Line 568 \| Line 365 \| void SimSetup::makeMolecules( void ){
365		the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
366		#endif //is_mpi
367		}
368	<	excludeOffset += info.nTorsions;
368	>	excludeOffset += molInfo.nTorsions;
369
370
371		// send the arrays off to the forceField for init.
372
373	<	the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
374	<	the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
375	<	the_ff->initializeBends( info.nBends, info.myBends, theBends );
376	<	the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
373	>	the_ff->initializeAtoms( molInfo.nAtoms, molInfo.myAtoms );
374	>	the_ff->initializeBonds( molInfo.nBonds, molInfo.myBonds, theBonds );
375	>	the_ff->initializeBends( molInfo.nBends, molInfo.myBends, theBends );
376	>	the_ff->initializeTorsions( molInfo.nTorsions, molInfo.myTorsions, theTorsions );
377
378
379	<	the_molecules[i].initialize( info );
379	>	the_molecules[i].initialize( molInfo );
380
381
382	<	atomOffset += info.nAtoms;
382	>	atomOffset += molInfo.nAtoms;
383		delete[] theBonds;
384		delete[] theBends;
385		delete[] theTorsions;
#	Line 618 \| Line 415 \| void SimSetup::initFromBass( void ){
415		have_extra =1;
416
417		n_cells = (int)temp3 - 1;
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		n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
422		temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
423		n_per_extra = (int)ceil( temp1 );
#	Line 635 \| Line 432 \| void SimSetup::initFromBass( void ){
432		}
433		else{
434		n_cells = (int)temp3;
435	<	cellx = info->boxLx / temp3;
436	<	celly = info->boxLy / temp3;
437	<	cellz = info->boxLz / temp3;
435	>	cellx = info->boxL[0] / temp3;
436	>	celly = info->boxL[1] / temp3;
437	>	cellz = info->boxL[2] / temp3;
438		}
439
440		current_mol = 0;
#	Line 784 \| Line 581 \| void SimSetup::gatherInfo( void ){
581
582
583		void SimSetup::gatherInfo( void ){
584	+	int i,j,k;
585
586		ensembleCase = -1;
587		ffCase = -1;
588
589		// get the stamps and globals;
590	<	the_stamps = stamps;
591	<	the_globals = globals;
590	>	stamps = stamps;
591	>	globals = globals;
592
593		// set the easy ones first
594	<	info->target_temp = the_globals->getTargetTemp();
595	<	info->dt = the_globals->getDt();
596	<	info->run_time = the_globals->getRunTime();
597	<	n_components = the_globals->getNComponents();
594	>	info->target_temp = globals->getTargetTemp();
595	>	info->dt = globals->getDt();
596	>	info->run_time = globals->getRunTime();
597	>	n_components = globals->getNComponents();
598
599
600		// get the forceField
601
602	<	strcpy( force_field, the_globals->getForceField() );
602	>	strcpy( force_field, globals->getForceField() );
603
604		if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
605		else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
606	+	else if( !strcasecmp( force_field, "EAM" )) ffCase = FF_EAM;
607		else{
608		sprintf( painCave.errMsg,
609		"SimSetup Error. Unrecognized force field -> %s\n",
#	Line 815 \| Line 614 \| void SimSetup::gatherInfo( void ){
614
615		// get the ensemble
616
617	<	strcpy( ensemble, the_globals->getEnsemble() );
617	>	strcpy( ensemble, globals->getEnsemble() );
618
619		if( !strcasecmp( ensemble, "NVE" )) ensembleCase = NVE_ENS;
620		else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
#	Line 824 \| Line 623 \| void SimSetup::gatherInfo( void ){
623		else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
624		else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
625		else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS;
626	+	else if( !strcasecmp( ensemble, "NVEZCONS")) ensembleCase = NVEZCONS_ENS;
627		else{
628		sprintf( painCave.errMsg,
629		"SimSetup Warning. Unrecognized Ensemble -> %s, "
#	Line 838 \| Line 638 \| void SimSetup::gatherInfo( void ){
638
639		// get the mixing rule
640
641	<	strcpy( info->mixingRule, the_globals->getMixingRule() );
642	<	info->usePBC = the_globals->getPBC();
641	>	strcpy( info->mixingRule, globals->getMixingRule() );
642	>	info->usePBC = globals->getPBC();
643
644
645		// get the components and calculate the tot_nMol and indvidual n_mol
646
647	<	the_components = the_globals->getComponents();
647	>	the_components = globals->getComponents();
648		components_nmol = new int[n_components];
649
650
651	<	if( !the_globals->haveNMol() ){
651	>	if( !globals->haveNMol() ){
652		// we don't have the total number of molecules, so we assume it is
653		// given in each component
654
#	Line 881 \| Line 681 \| void SimSetup::gatherInfo( void ){
681
682		// set the status, sample, and thermal kick times
683
684	<	if( the_globals->haveSampleTime() ){
685	<	info->sampleTime = the_globals->getSampleTime();
684	>	if( globals->haveSampleTime() ){
685	>	info->sampleTime = globals->getSampleTime();
686		info->statusTime = info->sampleTime;
687		info->thermalTime = info->sampleTime;
688		}
689		else{
690	<	info->sampleTime = the_globals->getRunTime();
690	>	info->sampleTime = globals->getRunTime();
691		info->statusTime = info->sampleTime;
692		info->thermalTime = info->sampleTime;
693		}
694
695	<	if( the_globals->haveStatusTime() ){
696	<	info->statusTime = the_globals->getStatusTime();
695	>	if( globals->haveStatusTime() ){
696	>	info->statusTime = globals->getStatusTime();
697		}
698
699	<	if( the_globals->haveThermalTime() ){
700	<	info->thermalTime = the_globals->getThermalTime();
699	>	if( globals->haveThermalTime() ){
700	>	info->thermalTime = globals->getThermalTime();
701		}
702
703		// check for the temperature set flag
704
705	<	if( the_globals->haveTempSet() ) info->setTemp = the_globals->getTempSet();
705	>	if( globals->haveTempSet() ) info->setTemp = globals->getTempSet();
706
707		// get some of the tricky things that may still be in the globals
708
709		double boxVector[3];
710	<	if( the_globals->haveBox() ){
711	<	boxVector[0] = the_globals->getBox();
712	<	boxVector[1] = the_globals->getBox();
713	<	boxVector[2] = the_globals->getBox();
710	>	if( globals->haveBox() ){
711	>	boxVector[0] = globals->getBox();
712	>	boxVector[1] = globals->getBox();
713	>	boxVector[2] = globals->getBox();
714
715		info->setBox( boxVector );
716		}
717	<	else if( the_globals->haveDensity() ){
717	>	else if( globals->haveDensity() ){
718
719		double vol;
720	<	vol = (double)tot_nmol / the_globals->getDensity();
720	>	vol = (double)tot_nmol / globals->getDensity();
721		boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
722		boxVector[1] = boxVector[0];
723		boxVector[2] = boxVector[0];
#	Line 925 \| Line 725 \| void SimSetup::gatherInfo( void ){
725		info->setBox( boxVector );
726		}
727		else{
728	<	if( !the_globals->haveBoxX() ){
728	>	if( !globals->haveBoxX() ){
729		sprintf( painCave.errMsg,
730		"SimSetup error, no periodic BoxX size given.\n" );
731		painCave.isFatal = 1;
732		simError();
733		}
734	<	boxVector[0] = the_globals->getBoxX();
734	>	boxVector[0] = globals->getBoxX();
735
736	<	if( !the_globals->haveBoxY() ){
736	>	if( !globals->haveBoxY() ){
737		sprintf( painCave.errMsg,
738		"SimSetup error, no periodic BoxY size given.\n" );
739		painCave.isFatal = 1;
740		simError();
741		}
742	<	boxVector[1] = the_globals->getBoxY();
742	>	boxVector[1] = globals->getBoxY();
743
744	<	if( !the_globals->haveBoxZ() ){
744	>	if( !globals->haveBoxZ() ){
745		sprintf( painCave.errMsg,
746		"SimSetup error, no periodic BoxZ size given.\n" );
747		painCave.isFatal = 1;
748		simError();
749		}
750	<	boxVector[2] = the_globals->getBoxZ();
750	>	boxVector[2] = globals->getBoxZ();
751
752		info->setBox( boxVector );
753		}
#	Line 977 \| Line 777 \| void SimSetup::finalInfoCheck( void ){
777		}
778
779		#ifdef IS_MPI
780	<	int myUse = usesDipoles
781	<	MPI_Allreduce( &myUse, &UsesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
780	>	int myUse = usesDipoles;
781	>	MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
782		#endif //is_mpi
783
784	+	double theEcr, theEst;
785
786	<	if (the_globals->getUseRF() ) {
786	>	if (globals->getUseRF() ) {
787		info->useReactionField = 1;
788
789	<	if( !the_globals->haveECR() ){
789	>	if( !globals->haveECR() ){
790		sprintf( painCave.errMsg,
791		"SimSetup Warning: using default value of 1/2 the smallest "
792		"box length for the electrostaticCutoffRadius.\n"
#	Line 993 \| Line 794 \| void SimSetup::finalInfoCheck( void ){
794		painCave.isFatal = 0;
795		simError();
796		double smallest;
797	<	smallest = info->boxLx;
798	<	if (info->boxLy <= smallest) smallest = info->boxLy;
799	<	if (info->boxLz <= smallest) smallest = info->boxLz;
800	<	info->ecr = 0.5 * smallest;
797	>	smallest = info->boxL[0];
798	>	if (info->boxL[1] <= smallest) smallest = info->boxL[1];
799	>	if (info->boxL[2] <= smallest) smallest = info->boxL[2];
800	>	theEcr = 0.5 * smallest;
801		} else {
802	<	info->ecr = the_globals->getECR();
802	>	theEcr = globals->getECR();
803		}
804
805	<	if( !the_globals->haveEST() ){
805	>	if( !globals->haveEST() ){
806		sprintf( painCave.errMsg,
807		"SimSetup Warning: using default value of 0.05 * the "
808		"electrostaticCutoffRadius for the electrostaticSkinThickness\n"
809		);
810		painCave.isFatal = 0;
811		simError();
812	<	info->est = 0.05 * info->ecr;
812	>	theEst = 0.05 * theEcr;
813		} else {
814	<	info->est = the_globals->getEST();
814	>	theEst= globals->getEST();
815		}
816	+
817	+	info->setEcr( theEcr, theEst );
818
819	<	if(!the_globals->haveDielectric() ){
819	>	if(!globals->haveDielectric() ){
820		sprintf( painCave.errMsg,
821		"SimSetup Error: You are trying to use Reaction Field without"
822		"setting a dielectric constant!\n"
#	Line 1021 \| Line 824 \| void SimSetup::finalInfoCheck( void ){
824		painCave.isFatal = 1;
825		simError();
826		}
827	<	info->dielectric = the_globals->getDielectric();
827	>	info->dielectric = globals->getDielectric();
828		}
829		else {
830		if (usesDipoles) {
831
832	<	if( !the_globals->haveECR() ){
833	<	sprintf( painCave.errMsg,
834	<	"SimSetup Warning: using default value of 1/2 the smallest "
835	<	"box length for the electrostaticCutoffRadius.\n"
836	<	"I hope you have a very fast processor!\n");
837	<	painCave.isFatal = 0;
838	<	simError();
839	<	double smallest;
840	<	smallest = info->boxLx;
841	<	if (info->boxLy <= smallest) smallest = info->boxLy;
842	<	if (info->boxLz <= smallest) smallest = info->boxLz;
843	<	info->ecr = 0.5 * smallest;
832	>	if( !globals->haveECR() ){
833	>	sprintf( painCave.errMsg,
834	>	"SimSetup Warning: using default value of 1/2 the smallest "
835	>	"box length for the electrostaticCutoffRadius.\n"
836	>	"I hope you have a very fast processor!\n");
837	>	painCave.isFatal = 0;
838	>	simError();
839	>	double smallest;
840	>	smallest = info->boxL[0];
841	>	if (info->boxL[1] <= smallest) smallest = info->boxL[1];
842	>	if (info->boxL[2] <= smallest) smallest = info->boxL[2];
843	>	theEcr = 0.5 * smallest;
844		} else {
845	<	info->ecr = the_globals->getECR();
845	>	theEcr = globals->getECR();
846		}
847
848	<	if( !the_globals->haveEST() ){
849	<	sprintf( painCave.errMsg,
850	<	"SimSetup Warning: using default value of 5%% of the "
851	<	"electrostaticCutoffRadius for the "
852	<	"electrostaticSkinThickness\n"
853	<	);
854	<	painCave.isFatal = 0;
855	<	simError();
856	<	info->est = 0.05 * info->ecr;
848	>	if( !globals->haveEST() ){
849	>	sprintf( painCave.errMsg,
850	>	"SimSetup Warning: using default value of 0.05 * the "
851	>	"electrostaticCutoffRadius for the "
852	>	"electrostaticSkinThickness\n"
853	>	);
854	>	painCave.isFatal = 0;
855	>	simError();
856	>	theEst = 0.05 * theEcr;
857		} else {
858	<	info->est = the_globals->getEST();
858	>	theEst= globals->getEST();
859		}
860	+
861	+	info->setEcr( theEcr, theEst );
862		}
863		}
864
#	Line 1066 \| Line 871 \| void SimSetup::initSystemCoords( void ){
871
872		void SimSetup::initSystemCoords( void ){
873
874	<	if( the_globals->haveInitialConfig() ){
874	>	if( globals->haveInitialConfig() ){
875
876		InitializeFromFile* fileInit;
877		#ifdef IS_MPI // is_mpi
878		if( worldRank == 0 ){
879		#endif //is_mpi
880	<	fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
880	>	fileInit = new InitializeFromFile( globals->getInitialConfig() );
881		#ifdef IS_MPI
882		}else fileInit = new InitializeFromFile( NULL );
883		#endif
884	<	fileInit->read_xyz( info ); // default velocities on
884	>	fileInit->readInit( info ); // default velocities on
885
886		delete fileInit;
887		}
#	Line 1113 \| Line 918 \| void SimSetup::makeOutNames( void ){
918		if( worldRank == 0 ){
919		#endif // is_mpi
920
921	<	if( the_globals->haveFinalConfig() ){
922	<	strcpy( info->finalName, the_globals->getFinalConfig() );
921	>	if( globals->haveFinalConfig() ){
922	>	strcpy( info->finalName, globals->getFinalConfig() );
923		}
924		else{
925		strcpy( info->finalName, inFileName );
#	Line 1197 \| Line 1002 \| void SimSetup::sysObjectsCreation( void ){
1002
1003		void SimSetup::sysObjectsCreation( void ){
1004
1005	+	int i;
1006	+
1007		// create the forceField
1008
1009		createFF();
#	Line 1219 \| Line 1026 \| void SimSetup::sysObjectsCreation( void ){
1026
1027		makeSysArrays();
1028
1029	+	// make and initialize the molecules (all but atomic coordinates)
1030	+
1031	+	makeMolecules();
1032	+	info->identArray = new int[info->n_atoms];
1033	+	for(i=0; i<info->n_atoms; i++){
1034	+	info->identArray[i] = the_atoms[i]->getIdent();
1035	+	}
1036
1037
1038
#	Line 1237 \| Line 1051 \| void SimSetup::createFF( void ){
1051		the_ff = new LJFF();
1052		break;
1053
1054	+	case FF_EAM:
1055	+	the_ff = new EAM_FF();
1056	+	break;
1057	+
1058		default:
1059		sprintf( painCave.errMsg,
1060		"SimSetup Error. Unrecognized force field in case statement.\n");
#	Line 1254 \| Line 1072 \| void SimSetup::compList( void ){
1072
1073		void SimSetup::compList( void ){
1074
1075	+	int i;
1076	+
1077		comp_stamps = new MoleculeStamp*[n_components];
1078
1079		// make an array of molecule stamps that match the components used.
#	Line 1279 \| Line 1099 \| void SimSetup::compList( void ){
1099
1100		// extract the component from the list;
1101
1102	<	currentStamp = the_stamps->extractMolStamp( id );
1102	>	currentStamp = stamps->extractMolStamp( id );
1103		if( currentStamp == NULL ){
1104		sprintf( painCave.errMsg,
1105		"SimSetup error: Component \"%s\" was not found in the "
#	Line 1303 \| Line 1123 \| void SimSetup::calcSysValues( void ){
1123		}
1124
1125		void SimSetup::calcSysValues( void ){
1126	+	int i, j, k;
1127
1128	+
1129		tot_atoms = 0;
1130		tot_bonds = 0;
1131		tot_bends = 0;
#	Line 1333 \| Line 1155 \| void SimSetup::mpiMolDivide( void ){
1155
1156		void SimSetup::mpiMolDivide( void ){
1157
1158	+	int i, j, k;
1159		int localMol, allMol;
1160		int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1161
#	Line 1402 \| Line 1225 \| void SimSetup::makeSysArrays( void ){
1225
1226
1227		void SimSetup::makeSysArrays( void ){
1228	+	int i, j, k;
1229
1230	+
1231		// create the atom and short range interaction arrays
1232
1233		Atom::createArrays(info->n_atoms);
#	Line 1476 \| Line 1301 \| void SimSetup::makeSysArrays( void ){
1301		the_ff->setSimInfo( info );
1302
1303		}
1304	+
1305	+	void SimSetup::makeIntegrator( void ){
1306	+
1307	+	NVT<RealIntegrator>* myNVT = NULL;
1308	+	NPTi<RealIntegrator>* myNPTi = NULL;
1309	+	NPTf<RealIntegrator>* myNPTf = NULL;
1310	+	NPTim<RealIntegrator>* myNPTim = NULL;
1311	+	NPTfm<RealIntegrator>* myNPTfm = NULL;
1312	+	ZConstraint<NVE<RealIntegrator> >* myNVEZCons = NULL;
1313	+
1314	+	cerr << "setting integrator" <<endl;
1315	+
1316	+	switch( ensembleCase ){
1317	+
1318	+	case NVE_ENS:
1319	+	new NVE<RealIntegrator>( info, the_ff );
1320	+	break;
1321	+
1322	+	case NVT_ENS:
1323	+	myNVT = new NVT<RealIntegrator>( info, the_ff );
1324	+	myNVT->setTargetTemp(globals->getTargetTemp());
1325	+
1326	+	if (globals->haveTauThermostat())
1327	+	myNVT->setTauThermostat(globals->getTauThermostat());
1328	+
1329	+	else {
1330	+	sprintf( painCave.errMsg,
1331	+	"SimSetup error: If you use the NVT\n"
1332	+	" ensemble, you must set tauThermostat.\n");
1333	+	painCave.isFatal = 1;
1334	+	simError();
1335	+	}
1336	+	break;
1337	+
1338	+	case NPTi_ENS:
1339	+	myNPTi = new NPTi<RealIntegrator>( info, the_ff );
1340	+	myNPTi->setTargetTemp( globals->getTargetTemp() );
1341	+
1342	+	if (globals->haveTargetPressure())
1343	+	myNPTi->setTargetPressure(globals->getTargetPressure());
1344	+	else {
1345	+	sprintf( painCave.errMsg,
1346	+	"SimSetup error: If you use a constant pressure\n"
1347	+	" ensemble, you must set targetPressure in the BASS file.\n");
1348	+	painCave.isFatal = 1;
1349	+	simError();
1350	+	}
1351	+
1352	+	if( globals->haveTauThermostat() )
1353	+	myNPTi->setTauThermostat( globals->getTauThermostat() );
1354	+	else{
1355	+	sprintf( painCave.errMsg,
1356	+	"SimSetup error: If you use an NPT\n"
1357	+	" ensemble, you must set tauThermostat.\n");
1358	+	painCave.isFatal = 1;
1359	+	simError();
1360	+	}
1361	+
1362	+	if( globals->haveTauBarostat() )
1363	+	myNPTi->setTauBarostat( globals->getTauBarostat() );
1364	+	else{
1365	+	sprintf( painCave.errMsg,
1366	+	"SimSetup error: If you use an NPT\n"
1367	+	" ensemble, you must set tauBarostat.\n");
1368	+	painCave.isFatal = 1;
1369	+	simError();
1370	+	}
1371	+	break;
1372	+
1373	+	case NPTf_ENS:
1374	+	myNPTf = new NPTf<RealIntegrator>( info, the_ff );
1375	+	myNPTf->setTargetTemp( globals->getTargetTemp());
1376	+
1377	+	if (globals->haveTargetPressure())
1378	+	myNPTf->setTargetPressure(globals->getTargetPressure());
1379	+	else {
1380	+	sprintf( painCave.errMsg,
1381	+	"SimSetup error: If you use a constant pressure\n"
1382	+	" ensemble, you must set targetPressure in the BASS file.\n");
1383	+	painCave.isFatal = 1;
1384	+	simError();
1385	+	}
1386	+
1387	+	if( globals->haveTauThermostat() )
1388	+	myNPTf->setTauThermostat( globals->getTauThermostat() );
1389	+	else{
1390	+	sprintf( painCave.errMsg,
1391	+	"SimSetup error: If you use an NPT\n"
1392	+	" ensemble, you must set tauThermostat.\n");
1393	+	painCave.isFatal = 1;
1394	+	simError();
1395	+	}
1396	+
1397	+	if( globals->haveTauBarostat() )
1398	+	myNPTf->setTauBarostat( globals->getTauBarostat() );
1399	+	else{
1400	+	sprintf( painCave.errMsg,
1401	+	"SimSetup error: If you use an NPT\n"
1402	+	" ensemble, you must set tauBarostat.\n");
1403	+	painCave.isFatal = 1;
1404	+	simError();
1405	+	}
1406	+	break;
1407	+
1408	+	case NPTim_ENS:
1409	+	myNPTim = new NPTim<RealIntegrator>( info, the_ff );
1410	+	myNPTim->setTargetTemp( globals->getTargetTemp());
1411	+
1412	+	if (globals->haveTargetPressure())
1413	+	myNPTim->setTargetPressure(globals->getTargetPressure());
1414	+	else {
1415	+	sprintf( painCave.errMsg,
1416	+	"SimSetup error: If you use a constant pressure\n"
1417	+	" ensemble, you must set targetPressure in the BASS file.\n");
1418	+	painCave.isFatal = 1;
1419	+	simError();
1420	+	}
1421	+
1422	+	if( globals->haveTauThermostat() )
1423	+	myNPTim->setTauThermostat( globals->getTauThermostat() );
1424	+	else{
1425	+	sprintf( painCave.errMsg,
1426	+	"SimSetup error: If you use an NPT\n"
1427	+	" ensemble, you must set tauThermostat.\n");
1428	+	painCave.isFatal = 1;
1429	+	simError();
1430	+	}
1431	+
1432	+	if( globals->haveTauBarostat() )
1433	+	myNPTim->setTauBarostat( globals->getTauBarostat() );
1434	+	else{
1435	+	sprintf( painCave.errMsg,
1436	+	"SimSetup error: If you use an NPT\n"
1437	+	" ensemble, you must set tauBarostat.\n");
1438	+	painCave.isFatal = 1;
1439	+	simError();
1440	+	}
1441	+	break;
1442	+
1443	+	case NPTfm_ENS:
1444	+	myNPTfm = new NPTfm<RealIntegrator>( info, the_ff );
1445	+	myNPTfm->setTargetTemp( globals->getTargetTemp());
1446	+
1447	+	if (globals->haveTargetPressure())
1448	+	myNPTfm->setTargetPressure(globals->getTargetPressure());
1449	+	else {
1450	+	sprintf( painCave.errMsg,
1451	+	"SimSetup error: If you use a constant pressure\n"
1452	+	" ensemble, you must set targetPressure in the BASS file.\n");
1453	+	painCave.isFatal = 1;
1454	+	simError();
1455	+	}
1456	+
1457	+	if( globals->haveTauThermostat() )
1458	+	myNPTfm->setTauThermostat( globals->getTauThermostat() );
1459	+	else{
1460	+	sprintf( painCave.errMsg,
1461	+	"SimSetup error: If you use an NPT\n"
1462	+	" ensemble, you must set tauThermostat.\n");
1463	+	painCave.isFatal = 1;
1464	+	simError();
1465	+	}
1466	+
1467	+	if( globals->haveTauBarostat() )
1468	+	myNPTfm->setTauBarostat( globals->getTauBarostat() );
1469	+	else{
1470	+	sprintf( painCave.errMsg,
1471	+	"SimSetup error: If you use an NPT\n"
1472	+	" ensemble, you must set tauBarostat.\n");
1473	+	painCave.isFatal = 1;
1474	+	simError();
1475	+	}
1476	+	break;
1477	+
1478	+	case NVEZCONS_ENS:
1479	+	{
1480	+
1481	+	if(globals->haveZConsTime()){
1482	+
1483	+	//add sample time of z-constraint into SimInfo's property list
1484	+	DoubleData* zconsTimeProp = new DoubleData();
1485	+	zconsTimeProp->setID("zconstime");
1486	+	zconsTimeProp->setData(globals->getZConsTime());
1487	+	info->addProperty(zconsTimeProp);
1488	+	}
1489	+	else{
1490	+	sprintf( painCave.errMsg,
1491	+	"ZConstraint error: If you use an ZConstraint\n"
1492	+	" , you must set sample time.\n");
1493	+	painCave.isFatal = 1;
1494	+	simError();
1495	+	}
1496	+
1497	+	if(globals->haveIndexOfAllZConsMols()){
1498	+
1499	+	//add index of z-constraint molecules into SimInfo's property list
1500	+	vector<int> tempIndex = globals->getIndexOfAllZConsMols();
1501	+	sort(tempIndex.begin(), tempIndex.end());
1502	+
1503	+	IndexData* zconsIndex = new IndexData();
1504	+	zconsIndex->setID("zconsindex");
1505	+	zconsIndex->setIndexData(tempIndex);
1506	+	info->addProperty(zconsIndex);
1507	+	}
1508	+	else{
1509	+	sprintf( painCave.errMsg,
1510	+	"SimSetup error: If you use an ZConstraint\n"
1511	+	" , you must set index of z-constraint molecules.\n");
1512	+	painCave.isFatal = 1;
1513	+	simError();
1514	+
1515	+	}
1516	+
1517	+	//Determine the name of ouput file and add it into SimInfo's property list
1518	+	//Be careful, do not use inFileName, since it is a pointer which
1519	+	//point to a string at master node, and slave nodes do not contain that string
1520	+
1521	+	string zconsOutput(info->finalName);
1522	+
1523	+	zconsOutput = zconsOutput.substr(0, zconsOutput.rfind(".")) + ".fz";
1524	+
1525	+	StringData* zconsFilename = new StringData();
1526	+	zconsFilename->setID("zconsfilename");
1527	+	zconsFilename->setData(zconsOutput);
1528	+
1529	+	info->addProperty(zconsFilename);
1530	+
1531	+	myNVEZCons = new ZConstraint<NVE<RealIntegrator> >( info, the_ff );
1532	+
1533	+	break;
1534	+	}
1535	+
1536	+	default:
1537	+	sprintf( painCave.errMsg,
1538	+	"SimSetup Error. Unrecognized ensemble in case statement.\n");
1539	+	painCave.isFatal = 1;
1540	+	simError();
1541	+	}
1542	+
1543	+	}
1544	+
1545	+	void SimSetup::initFortran( void ){
1546	+
1547	+	info->refreshSim();
1548	+
1549	+	if( !strcmp( info->mixingRule, "standard") ){
1550	+	the_ff->initForceField( LB_MIXING_RULE );
1551	+	}
1552	+	else if( !strcmp( info->mixingRule, "explicit") ){
1553	+	the_ff->initForceField( EXPLICIT_MIXING_RULE );
1554	+	}
1555	+	else{
1556	+	sprintf( painCave.errMsg,
1557	+	"SimSetup Error: unknown mixing rule -> \"%s\"\n",
1558	+	info->mixingRule );
1559	+	painCave.isFatal = 1;
1560	+	simError();
1561	+	}
1562	+
1563	+
1564	+	#ifdef IS_MPI
1565	+	strcpy( checkPointMsg,
1566	+	"Successfully intialized the mixingRule for Fortran." );
1567	+	MPIcheckPoint();
1568	+	#endif // is_mpi
1569	+
1570	+	}

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 658 by tim, Thu Jul 31 15:35:07 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 658 by tim, Thu Jul 31 15:35:07 2003 UTC