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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 414 by mmeineke, Wed Mar 26 22:02:36 2003 UTC vs.
Revision 605 by gezelter, Tue Jul 15 03:27:24 2003 UTC

#	Line 12 \| Line 12 \| SimSetup::SimSetup(){
12		#include "mpiSimulation.hpp"
13		#endif
14
15	+	// some defines for ensemble and Forcefield cases
16	+
17	+	#define NVE_ENS 0
18	+	#define NVT_ENS 1
19	+	#define NPTi_ENS 2
20	+	#define NPTf_ENS 3
21	+	#define NPTim_ENS 4
22	+	#define NPTfm_ENS 5
23	+
24	+
25	+	#define FF_DUFF 0
26	+	#define FF_LJ 1
27	+
28	+
29		SimSetup::SimSetup(){
30		stamps = new MakeStamps();
31		globals = new Globals();
#	Line 66 \| Line 80 \| void SimSetup::createSim( void ){
80
81		MakeStamps *the_stamps;
82		Globals* the_globals;
83	<	int i, j;
83	>	int i, j, k, globalAtomIndex;
84	>
85	>	int ensembleCase;
86	>	int ffCase;
87	>
88	>	ensembleCase = -1;
89	>	ffCase = -1;
90
91		// get the stamps and globals;
92		the_stamps = stamps;
#	Line 80 \| Line 100 \| void SimSetup::createSim( void ){
100		// get the ones we know are there, yet still may need some work.
101		n_components = the_globals->getNComponents();
102		strcpy( force_field, the_globals->getForceField() );
103	+
104	+	if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
105	+	else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
106	+	else{
107	+	sprintf( painCave.errMsg,
108	+	"SimSetup Error. Unrecognized force field -> %s\n",
109	+	force_field );
110	+	painCave.isFatal = 1;
111	+	simError();
112	+	}
113	+
114	+	// get the ensemble:
115		strcpy( ensemble, the_globals->getEnsemble() );
116	+
117	+	if( !strcasecmp( ensemble, "NVE" )) ensembleCase = NVE_ENS;
118	+	else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
119	+	else if( !strcasecmp( ensemble, "NPTi" ) \|\| !strcasecmp( ensemble, "NPT") )
120	+	ensembleCase = NPTi_ENS;
121	+	else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
122	+	else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
123	+	else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS;
124	+	else{
125	+	sprintf( painCave.errMsg,
126	+	"SimSetup Warning. Unrecognized Ensemble -> %s, "
127	+	"reverting to NVE for this simulation.\n",
128	+	ensemble );
129	+	painCave.isFatal = 0;
130	+	simError();
131	+	strcpy( ensemble, "NVE" );
132	+	ensembleCase = NVE_ENS;
133	+	}
134		strcpy( simnfo->ensemble, ensemble );
135
136	+
137	+	// if( !strcasecmp( ensemble, "NPT" ) ) {
138	+	// the_extendedsystem = new ExtendedSystem( simnfo );
139	+	// the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
140	+	// if (the_globals->haveTargetPressure())
141	+	// the_extendedsystem->setTargetPressure(the_globals->getTargetPressure());
142	+	// else {
143	+	// sprintf( painCave.errMsg,
144	+	// "SimSetup error: If you use the constant pressure\n"
145	+	// " ensemble, you must set targetPressure.\n"
146	+	// " This was found in the BASS file.\n");
147	+	// painCave.isFatal = 1;
148	+	// simError();
149	+	// }
150	+
151	+	// if (the_globals->haveTauThermostat())
152	+	// the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
153	+	// else if (the_globals->haveQmass())
154	+	// the_extendedsystem->setQmass(the_globals->getQmass());
155	+	// else {
156	+	// sprintf( painCave.errMsg,
157	+	// "SimSetup error: If you use one of the constant temperature\n"
158	+	// " ensembles, you must set either tauThermostat or qMass.\n"
159	+	// " Neither of these was found in the BASS file.\n");
160	+	// painCave.isFatal = 1;
161	+	// simError();
162	+	// }
163	+
164	+	// if (the_globals->haveTauBarostat())
165	+	// the_extendedsystem->setTauBarostat(the_globals->getTauBarostat());
166	+	// else {
167	+	// sprintf( painCave.errMsg,
168	+	// "SimSetup error: If you use the constant pressure\n"
169	+	// " ensemble, you must set tauBarostat.\n"
170	+	// " This was found in the BASS file.\n");
171	+	// painCave.isFatal = 1;
172	+	// simError();
173	+	// }
174	+
175	+	// } else if ( !strcasecmp( ensemble, "NVT") ) {
176	+	// the_extendedsystem = new ExtendedSystem( simnfo );
177	+	// the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
178	+
179	+	// if (the_globals->haveTauThermostat())
180	+	// the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
181	+	// else if (the_globals->haveQmass())
182	+	// the_extendedsystem->setQmass(the_globals->getQmass());
183	+	// else {
184	+	// sprintf( painCave.errMsg,
185	+	// "SimSetup error: If you use one of the constant temperature\n"
186	+	// " ensembles, you must set either tauThermostat or qMass.\n"
187	+	// " Neither of these was found in the BASS file.\n");
188	+	// painCave.isFatal = 1;
189	+	// simError();
190	+	// }
191	+
192		strcpy( simnfo->mixingRule, the_globals->getMixingRule() );
193		simnfo->usePBC = the_globals->getPBC();
194
195	+	int usesDipoles = 0;
196	+	switch( ffCase ){
197	+
198	+	case FF_DUFF:
199	+	the_ff = new DUFF();
200	+	usesDipoles = 1;
201	+	break;
202
203	+	case FF_LJ:
204	+	the_ff = new LJFF();
205	+	break;
206
207	<	if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
92	<	else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
93	<	else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
94	<	else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
95	<	else{
207	>	default:
208		sprintf( painCave.errMsg,
209	<	"SimSetup Error. Unrecognized force field -> %s\n",
98	<	force_field );
209	>	"SimSetup Error. Unrecognized force field in case statement.\n");
210		painCave.isFatal = 1;
211		simError();
212		}
#	Line 105 \| Line 216 \| void SimSetup::createSim( void ){
216		MPIcheckPoint();
217		#endif // is_mpi
218
108	–
109	–
219		// get the components and calculate the tot_nMol and indvidual n_mol
220		the_components = the_globals->getComponents();
221		components_nmol = new int[n_components];
#	Line 230 \| Line 339 \| void SimSetup::createSim( void ){
339		simnfo->n_torsions = tot_torsions;
340		simnfo->n_SRI = tot_SRI;
341		simnfo->n_mol = tot_nmol;
233	–
342
343	+	simnfo->molMembershipArray = new int[tot_atoms];
344	+
345		#ifdef IS_MPI
346
347		// divide the molecules among processors here.
348
349		mpiSim = new mpiSimulation( simnfo );
350
241	–
242	–
351		globalIndex = mpiSim->divideLabor();
352
245	–
246	–
353		// set up the local variables
354
355		int localMol, allMol;
356		int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
357	+
358	+	int* mol2proc = mpiSim->getMolToProcMap();
359	+	int* molCompType = mpiSim->getMolComponentType();
360
361		allMol = 0;
362		localMol = 0;
#	Line 255 \| Line 364 \| void SimSetup::createSim( void ){
364		local_bonds = 0;
365		local_bends = 0;
366		local_torsions = 0;
367	+	globalAtomIndex = 0;
368	+
369	+
370		for( i=0; i<n_components; i++ ){
371
372		for( j=0; j<components_nmol[i]; j++ ){
373
374	<	if( mpiSim->getMyMolStart() <= allMol &&
263	<	allMol <= mpiSim->getMyMolEnd() ){
374	>	if( mol2proc[allMol] == worldRank ){
375
376		local_atoms += comp_stamps[i]->getNAtoms();
377		local_bonds += comp_stamps[i]->getNBonds();
#	Line 268 \| Line 379 \| void SimSetup::createSim( void ){
379		local_torsions += comp_stamps[i]->getNTorsions();
380		localMol++;
381		}
382	<	allMol++;
382	>	for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
383	>	simnfo->molMembershipArray[globalAtomIndex] = allMol;
384	>	globalAtomIndex++;
385	>	}
386	>
387	>	allMol++;
388		}
389		}
390		local_SRI = local_bonds + local_bends + local_torsions;
391
276	–
392		simnfo->n_atoms = mpiSim->getMyNlocal();
393
394		if( local_atoms != simnfo->n_atoms ){
395		sprintf( painCave.errMsg,
396		"SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
397	<	" localAtom (%d) are note equal.\n",
397	>	" localAtom (%d) are not equal.\n",
398		simnfo->n_atoms,
399		local_atoms );
400		painCave.isFatal = 1;
#	Line 304 \| Line 419 \| void SimSetup::createSim( void ){
419		Atom::createArrays(simnfo->n_atoms);
420		the_atoms = new Atom*[simnfo->n_atoms];
421		the_molecules = new Molecule[simnfo->n_mol];
422	+	int molIndex;
423
424	+	// initialize the molecule's stampID's
425
426	+	#ifdef IS_MPI
427	+
428	+
429	+	molIndex = 0;
430	+	for(i=0; i<mpiSim->getTotNmol(); i++){
431	+
432	+	if(mol2proc[i] == worldRank ){
433	+	the_molecules[molIndex].setStampID( molCompType[i] );
434	+	the_molecules[molIndex].setMyIndex( molIndex );
435	+	the_molecules[molIndex].setGlobalIndex( i );
436	+	molIndex++;
437	+	}
438	+	}
439	+
440	+	#else // is_mpi
441	+
442	+	molIndex = 0;
443	+	globalAtomIndex = 0;
444	+	for(i=0; i<n_components; i++){
445	+	for(j=0; j<components_nmol[i]; j++ ){
446	+	the_molecules[molIndex].setStampID( i );
447	+	the_molecules[molIndex].setMyIndex( molIndex );
448	+	the_molecules[molIndex].setGlobalIndex( molIndex );
449	+	for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
450	+	simnfo->molMembershipArray[globalAtomIndex] = molIndex;
451	+	globalAtomIndex++;
452	+	}
453	+	molIndex++;
454	+	}
455	+	}
456	+
457	+
458	+	#endif // is_mpi
459	+
460	+
461		if( simnfo->n_SRI ){
462	+
463		Exclude::createArray(simnfo->n_SRI);
464		the_excludes = new Exclude*[simnfo->n_SRI];
465	+	for( int ex=0; ex<simnfo->n_SRI; ex++) the_excludes[ex] = new Exclude(ex);
466		simnfo->globalExcludes = new int;
467	<	simnfo->n_exclude = tot_SRI;
467	>	simnfo->n_exclude = simnfo->n_SRI;
468		}
469		else{
470
#	Line 326 \| Line 480 \| void SimSetup::createSim( void ){
480		// set the arrays into the SimInfo object
481
482		simnfo->atoms = the_atoms;
483	<	simnfo->sr_interactions = the_sris;
483	>	simnfo->molecules = the_molecules;
484		simnfo->nGlobalExcludes = 0;
485		simnfo->excludes = the_excludes;
486
487
488		// get some of the tricky things that may still be in the globals
489
490	<
490	>	double boxVector[3];
491		if( the_globals->haveBox() ){
492	<	simnfo->box_x = the_globals->getBox();
493	<	simnfo->box_y = the_globals->getBox();
494	<	simnfo->box_z = the_globals->getBox();
492	>	boxVector[0] = the_globals->getBox();
493	>	boxVector[1] = the_globals->getBox();
494	>	boxVector[2] = the_globals->getBox();
495	>
496	>	simnfo->setBox( boxVector );
497		}
498		else if( the_globals->haveDensity() ){
499
500		double vol;
501		vol = (double)tot_nmol / the_globals->getDensity();
502	<	simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
503	<	simnfo->box_y = simnfo->box_x;
504	<	simnfo->box_z = simnfo->box_x;
502	>	boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
503	>	boxVector[1] = boxVector[0];
504	>	boxVector[2] = boxVector[0];
505	>
506	>	simnfo->setBox( boxVector );
507		}
508		else{
509		if( !the_globals->haveBoxX() ){
#	Line 354 \| Line 512 \| void SimSetup::createSim( void ){
512		painCave.isFatal = 1;
513		simError();
514		}
515	<	simnfo->box_x = the_globals->getBoxX();
515	>	boxVector[0] = the_globals->getBoxX();
516
517		if( !the_globals->haveBoxY() ){
518		sprintf( painCave.errMsg,
#	Line 362 \| Line 520 \| void SimSetup::createSim( void ){
520		painCave.isFatal = 1;
521		simError();
522		}
523	<	simnfo->box_y = the_globals->getBoxY();
523	>	boxVector[1] = the_globals->getBoxY();
524
525		if( !the_globals->haveBoxZ() ){
526		sprintf( painCave.errMsg,
#	Line 370 \| Line 528 \| void SimSetup::createSim( void ){
528		painCave.isFatal = 1;
529		simError();
530		}
531	<	simnfo->box_z = the_globals->getBoxZ();
531	>	boxVector[2] = the_globals->getBoxZ();
532	>
533	>	simnfo->setBox( boxVector );
534		}
535
536		#ifdef IS_MPI
#	Line 383 \| Line 543 \| void SimSetup::createSim( void ){
543
544		the_ff->setSimInfo( simnfo );
545
546	<	makeAtoms();
546	>	makeMolecules();
547		simnfo->identArray = new int[simnfo->n_atoms];
548		for(i=0; i<simnfo->n_atoms; i++){
549		simnfo->identArray[i] = the_atoms[i]->getIdent();
550		}
551
392	–	if( tot_bonds ){
393	–	makeBonds();
394	–	}
395	–
396	–	if( tot_bends ){
397	–	makeBends();
398	–	}
399	–
400	–	if( tot_torsions ){
401	–	makeTorsions();
402	–	}
403	–
404	–
552		if (the_globals->getUseRF() ) {
553		simnfo->useReactionField = 1;
554
#	Line 413 \| Line 560 \| void SimSetup::createSim( void ){
560		painCave.isFatal = 0;
561		simError();
562		double smallest;
563	<	smallest = simnfo->box_x;
564	<	if (simnfo->box_y <= smallest) smallest = simnfo->box_y;
565	<	if (simnfo->box_z <= smallest) smallest = simnfo->box_z;
563	>	smallest = simnfo->boxLx;
564	>	if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy;
565	>	if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz;
566		simnfo->ecr = 0.5 * smallest;
567		} else {
568		simnfo->ecr = the_globals->getECR();
#	Line 443 \| Line 590 \| void SimSetup::createSim( void ){
590		}
591		simnfo->dielectric = the_globals->getDielectric();
592		} else {
593	<	if (simnfo->n_dipoles) {
593	>	if (usesDipoles) {
594
595		if( !the_globals->haveECR() ){
596		sprintf( painCave.errMsg,
597	<	"SimSetup Warning: using default value of 1/2 the smallest"
597	>	"SimSetup Warning: using default value of 1/2 the smallest "
598		"box length for the electrostaticCutoffRadius.\n"
599		"I hope you have a very fast processor!\n");
600		painCave.isFatal = 0;
601		simError();
602		double smallest;
603	<	smallest = simnfo->box_x;
604	<	if (simnfo->box_y <= smallest) smallest = simnfo->box_y;
605	<	if (simnfo->box_z <= smallest) smallest = simnfo->box_z;
603	>	smallest = simnfo->boxLx;
604	>	if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy;
605	>	if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz;
606		simnfo->ecr = 0.5 * smallest;
607		} else {
608		simnfo->ecr = the_globals->getECR();
#	Line 463 \| Line 610 \| void SimSetup::createSim( void ){
610
611		if( !the_globals->haveEST() ){
612		sprintf( painCave.errMsg,
613	<	"SimSetup Warning: using default value of 5% of the"
613	>	"SimSetup Warning: using default value of 5%% of the "
614		"electrostaticCutoffRadius for the "
615		"electrostaticSkinThickness\n"
616		);
#	Line 634 \| Line 781 \| void SimSetup::createSim( void ){
781		if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
782
783
784	<	// // make the longe range forces and the integrator
785	<
786	<	// new AllLong( simnfo );
787	<
788	<	if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
789	<	if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo, the_ff );
790	<	if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo, the_ff );
791	<	if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
784	>	// make the integrator
785	>
786	>
787	>	NVT* myNVT = NULL;
788	>	NPTi* myNPTi = NULL;
789	>	NPTf* myNPTf = NULL;
790	>	NPTim* myNPTim = NULL;
791	>	NPTfm* myNPTfm = NULL;
792
793	+	switch( ensembleCase ){
794
795	+	case NVE_ENS:
796	+	new NVE( simnfo, the_ff );
797	+	break;
798
799	+	case NVT_ENS:
800	+	myNVT = new NVT( simnfo, the_ff );
801	+	myNVT->setTargetTemp(the_globals->getTargetTemp());
802	+
803	+	if (the_globals->haveTauThermostat())
804	+	myNVT->setTauThermostat(the_globals->getTauThermostat());
805	+
806	+	else {
807	+	sprintf( painCave.errMsg,
808	+	"SimSetup error: If you use the NVT\n"
809	+	" ensemble, you must set tauThermostat.\n");
810	+	painCave.isFatal = 1;
811	+	simError();
812	+	}
813	+	break;
814	+
815	+	case NPTi_ENS:
816	+	myNPTi = new NPTi( simnfo, the_ff );
817	+	myNPTi->setTargetTemp( the_globals->getTargetTemp());
818	+
819	+	if (the_globals->haveTargetPressure())
820	+	myNPTi->setTargetPressure(the_globals->getTargetPressure());
821	+	else {
822	+	sprintf( painCave.errMsg,
823	+	"SimSetup error: If you use a constant pressure\n"
824	+	" ensemble, you must set targetPressure in the BASS file.\n");
825	+	painCave.isFatal = 1;
826	+	simError();
827	+	}
828	+
829	+	if( the_globals->haveTauThermostat() )
830	+	myNPTi->setTauThermostat( the_globals->getTauThermostat() );
831	+	else{
832	+	sprintf( painCave.errMsg,
833	+	"SimSetup error: If you use an NPT\n"
834	+	" ensemble, you must set tauThermostat.\n");
835	+	painCave.isFatal = 1;
836	+	simError();
837	+	}
838	+
839	+	if( the_globals->haveTauBarostat() )
840	+	myNPTi->setTauBarostat( the_globals->getTauBarostat() );
841	+	else{
842	+	sprintf( painCave.errMsg,
843	+	"SimSetup error: If you use an NPT\n"
844	+	" ensemble, you must set tauBarostat.\n");
845	+	painCave.isFatal = 1;
846	+	simError();
847	+	}
848	+	break;
849	+
850	+	case NPTf_ENS:
851	+	myNPTf = new NPTf( simnfo, the_ff );
852	+	myNPTf->setTargetTemp( the_globals->getTargetTemp());
853	+
854	+	if (the_globals->haveTargetPressure())
855	+	myNPTf->setTargetPressure(the_globals->getTargetPressure());
856	+	else {
857	+	sprintf( painCave.errMsg,
858	+	"SimSetup error: If you use a constant pressure\n"
859	+	" ensemble, you must set targetPressure in the BASS file.\n");
860	+	painCave.isFatal = 1;
861	+	simError();
862	+	}
863	+
864	+	if( the_globals->haveTauThermostat() )
865	+	myNPTf->setTauThermostat( the_globals->getTauThermostat() );
866	+	else{
867	+	sprintf( painCave.errMsg,
868	+	"SimSetup error: If you use an NPT\n"
869	+	" ensemble, you must set tauThermostat.\n");
870	+	painCave.isFatal = 1;
871	+	simError();
872	+	}
873	+
874	+	if( the_globals->haveTauBarostat() )
875	+	myNPTf->setTauBarostat( the_globals->getTauBarostat() );
876	+	else{
877	+	sprintf( painCave.errMsg,
878	+	"SimSetup error: If you use an NPT\n"
879	+	" ensemble, you must set tauBarostat.\n");
880	+	painCave.isFatal = 1;
881	+	simError();
882	+	}
883	+	break;
884	+
885	+	case NPTim_ENS:
886	+	myNPTim = new NPTim( simnfo, the_ff );
887	+	myNPTim->setTargetTemp( the_globals->getTargetTemp());
888	+
889	+	if (the_globals->haveTargetPressure())
890	+	myNPTim->setTargetPressure(the_globals->getTargetPressure());
891	+	else {
892	+	sprintf( painCave.errMsg,
893	+	"SimSetup error: If you use a constant pressure\n"
894	+	" ensemble, you must set targetPressure in the BASS file.\n");
895	+	painCave.isFatal = 1;
896	+	simError();
897	+	}
898	+
899	+	if( the_globals->haveTauThermostat() )
900	+	myNPTim->setTauThermostat( the_globals->getTauThermostat() );
901	+	else{
902	+	sprintf( painCave.errMsg,
903	+	"SimSetup error: If you use an NPT\n"
904	+	" ensemble, you must set tauThermostat.\n");
905	+	painCave.isFatal = 1;
906	+	simError();
907	+	}
908	+
909	+	if( the_globals->haveTauBarostat() )
910	+	myNPTim->setTauBarostat( the_globals->getTauBarostat() );
911	+	else{
912	+	sprintf( painCave.errMsg,
913	+	"SimSetup error: If you use an NPT\n"
914	+	" ensemble, you must set tauBarostat.\n");
915	+	painCave.isFatal = 1;
916	+	simError();
917	+	}
918	+	break;
919	+
920	+	case NPTfm_ENS:
921	+	myNPTfm = new NPTfm( simnfo, the_ff );
922	+	myNPTfm->setTargetTemp( the_globals->getTargetTemp());
923	+
924	+	if (the_globals->haveTargetPressure())
925	+	myNPTfm->setTargetPressure(the_globals->getTargetPressure());
926	+	else {
927	+	sprintf( painCave.errMsg,
928	+	"SimSetup error: If you use a constant pressure\n"
929	+	" ensemble, you must set targetPressure in the BASS file.\n");
930	+	painCave.isFatal = 1;
931	+	simError();
932	+	}
933	+
934	+	if( the_globals->haveTauThermostat() )
935	+	myNPTfm->setTauThermostat( the_globals->getTauThermostat() );
936	+	else{
937	+	sprintf( painCave.errMsg,
938	+	"SimSetup error: If you use an NPT\n"
939	+	" ensemble, you must set tauThermostat.\n");
940	+	painCave.isFatal = 1;
941	+	simError();
942	+	}
943	+
944	+	if( the_globals->haveTauBarostat() )
945	+	myNPTfm->setTauBarostat( the_globals->getTauBarostat() );
946	+	else{
947	+	sprintf( painCave.errMsg,
948	+	"SimSetup error: If you use an NPT\n"
949	+	" ensemble, you must set tauBarostat.\n");
950	+	painCave.isFatal = 1;
951	+	simError();
952	+	}
953	+	break;
954	+
955	+	default:
956	+	sprintf( painCave.errMsg,
957	+	"SimSetup Error. Unrecognized ensemble in case statement.\n");
958	+	painCave.isFatal = 1;
959	+	simError();
960	+	}
961	+
962	+
963	+	#ifdef IS_MPI
964	+	mpiSim->mpiRefresh();
965	+	#endif
966	+
967		// initialize the Fortran
968	<
968	>
969	>
970		simnfo->refreshSim();
971
972		if( !strcmp( simnfo->mixingRule, "standard") ){
#	Line 683 \| Line 1003 \| void SimSetup::makeMolecules( void ){
1003		BondStamp* currentBond;
1004		BendStamp* currentBend;
1005		TorsionStamp* currentTorsion;
1006	+
1007	+	bond_pair* theBonds;
1008	+	bend_set* theBends;
1009	+	torsion_set* theTorsions;
1010	+
1011
1012		//init the forceField paramters
1013
1014		the_ff->readParams();
1015
1016
1017	<	// init the molecules
1017	>	// init the atoms
1018
1019	+	double ux, uy, uz, u, uSqr;
1020	+
1021		atomOffset = 0;
1022		excludeOffset = 0;
1023		for(i=0; i<simnfo->n_mol; i++){
#	Line 707 \| Line 1034 \| void SimSetup::makeMolecules( void ){
1034		info.myExcludes = &the_excludes[excludeOffset];
1035		info.myBonds = new Bond*[info.nBonds];
1036		info.myBends = new Bend*[info.nBends];
1037	<	info.myTorsions = new Torsions*[info.nTorsions];
1037	>	info.myTorsions = new Torsion*[info.nTorsions];
1038
1039		theBonds = new bond_pair[info.nBonds];
1040		theBends = new bend_set[info.nBends];
#	Line 717 \| Line 1044 \| void SimSetup::makeMolecules( void ){
1044
1045		for(j=0; j<info.nAtoms; j++){
1046
1047	<	currentAtom = theComponents[stampID]->getAtom( j );
1047	>	currentAtom = comp_stamps[stampID]->getAtom( j );
1048		if( currentAtom->haveOrientation() ){
1049
1050		dAtom = new DirectionalAtom(j + atomOffset);
#	Line 758 \| Line 1085 \| void SimSetup::makeMolecules( void ){
1085		theBonds[j].a = currentBond->getA() + atomOffset;
1086		theBonds[j].b = currentBond->getB() + atomOffset;
1087
1088	<	exI = theBonds[i].a;
1089	<	exJ = theBonds[i].b;
1088	>	exI = theBonds[j].a;
1089	>	exJ = theBonds[j].b;
1090
1091		// exclude_I must always be the smaller of the pair
1092		if( exI > exJ ){
#	Line 775 \| Line 1102 \| void SimSetup::makeMolecules( void ){
1102
1103		the_excludes[j+excludeOffset]->setPair( exI, exJ );
1104		#else // isn't MPI
1105	+
1106		the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
1107		#endif //is_mpi
1108		}
#	Line 790 \| Line 1118 \| void SimSetup::makeMolecules( void ){
1118
1119		if( currentBend->haveExtras() ){
1120
1121	<	extras = current_bend->getExtras();
1121	>	extras = currentBend->getExtras();
1122		current_extra = extras;
1123
1124		while( current_extra != NULL ){
#	Line 812 \| Line 1140 \| void SimSetup::makeMolecules( void ){
1140
1141		default:
1142		sprintf( painCave.errMsg,
1143	<	"SimSetup Error: ghostVectorSource was neiter a "
1143	>	"SimSetup Error: ghostVectorSource was neither a "
1144		"double nor an int.\n"
1145		"-->Bend[%d] in %s\n",
1146		j, comp_stamps[stampID]->getID() );
#	Line 906 \| Line 1234 \| void SimSetup::makeMolecules( void ){
1234
1235
1236		the_molecules[i].initialize( info );
909	–	atomOffset += info.nAtoms;
910	–	}
911	–
912	–	// clean up the forcefield
913	–
914	–	the_ff->cleanMe();
915	–	}
916	–
917	–
1237
919	–	void SimSetup::makeAtoms( void ){
1238
1239	<	int i, j, k, index;
1240	<	double ux, uy, uz, uSqr, u;
1241	<	AtomStamp* current_atom;
1242	<
925	<	DirectionalAtom* dAtom;
926	<	int molIndex, molStart, molEnd, nMemb, lMolIndex;
927	<
928	<	lMolIndex = 0;
929	<	molIndex = 0;
930	<	index = 0;
931	<	for( i=0; i<n_components; i++ ){
932	<
933	<	for( j=0; j<components_nmol[i]; j++ ){
934	<
935	<	#ifdef IS_MPI
936	<	if( mpiSim->getMyMolStart() <= molIndex &&
937	<	molIndex <= mpiSim->getMyMolEnd() ){
938	<	#endif // is_mpi
939	<
940	<	molStart = index;
941	<	nMemb = comp_stamps[i]->getNAtoms();
942	<	for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
943	<
944	<	current_atom = comp_stamps[i]->getAtom( k );
945	<	if( current_atom->haveOrientation() ){
946	<
947	<	dAtom = new DirectionalAtom(index);
948	<	simnfo->n_oriented++;
949	<	the_atoms[index] = dAtom;
950	<
951	<	ux = current_atom->getOrntX();
952	<	uy = current_atom->getOrntY();
953	<	uz = current_atom->getOrntZ();
954	<
955	<	uSqr = (ux * ux) + (uy * uy) + (uz * uz);
956	<
957	<	u = sqrt( uSqr );
958	<	ux = ux / u;
959	<	uy = uy / u;
960	<	uz = uz / u;
961	<
962	<	dAtom->setSUx( ux );
963	<	dAtom->setSUy( uy );
964	<	dAtom->setSUz( uz );
965	<	}
966	<	else{
967	<	the_atoms[index] = new GeneralAtom(index);
968	<	}
969	<	the_atoms[index]->setType( current_atom->getType() );
970	<	the_atoms[index]->setIndex( index );
971	<
972	<	// increment the index and repeat;
973	<	index++;
974	<	}
975	<
976	<	molEnd = index -1;
977	<	the_molecules[lMolIndex].setNMembers( nMemb );
978	<	the_molecules[lMolIndex].setStartAtom( molStart );
979	<	the_molecules[lMolIndex].setEndAtom( molEnd );
980	<	the_molecules[lMolIndex].setStampID( i );
981	<	lMolIndex++;
982	<
983	<	#ifdef IS_MPI
984	<	}
985	<	#endif //is_mpi
986	<
987	<	molIndex++;
988	<	}
989	<	}
990	<
991	<	#ifdef IS_MPI
992	<	for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] );
993	<
994	<	delete[] globalIndex;
995	<
996	<	mpiSim->mpiRefresh();
997	<	#endif //IS_MPI
998	<
999	<	the_ff->initializeAtoms();
1000	<	}
1001	<
1002	<	void SimSetup::makeBonds( void ){
1003	<
1004	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1005	<	bond_pair* the_bonds;
1006	<	BondStamp* current_bond;
1007	<
1008	<	the_bonds = new bond_pair[tot_bonds];
1009	<	index = 0;
1010	<	offset = 0;
1011	<	molIndex = 0;
1012	<
1013	<	for( i=0; i<n_components; i++ ){
1014	<
1015	<	for( j=0; j<components_nmol[i]; j++ ){
1016	<
1017	<	#ifdef IS_MPI
1018	<	if( mpiSim->getMyMolStart() <= molIndex &&
1019	<	molIndex <= mpiSim->getMyMolEnd() ){
1020	<	#endif // is_mpi
1021	<
1022	<	for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
1023	<
1024	<	current_bond = comp_stamps[i]->getBond( k );
1025	<	the_bonds[index].a = current_bond->getA() + offset;
1026	<	the_bonds[index].b = current_bond->getB() + offset;
1027	<
1028	<	exI = the_bonds[index].a;
1029	<	exJ = the_bonds[index].b;
1030	<
1031	<	// exclude_I must always be the smaller of the pair
1032	<	if( exI > exJ ){
1033	<	tempEx = exI;
1034	<	exI = exJ;
1035	<	exJ = tempEx;
1036	<	}
1037	<
1038	<
1039	<	#ifdef IS_MPI
1040	<
1041	<	the_excludes[index*2] =
1042	<	the_atoms[exI]->getGlobalIndex() + 1;
1043	<	the_excludes[index*2 + 1] =
1044	<	the_atoms[exJ]->getGlobalIndex() + 1;
1045	<
1046	<	#else // isn't MPI
1047	<
1048	<	the_excludes[index*2] = exI + 1;
1049	<	the_excludes[index*2 + 1] = exJ + 1;
1050	<	// fortran index from 1 (hence the +1 in the indexing)
1051	<	#endif //is_mpi
1052	<
1053	<	// increment the index and repeat;
1054	<	index++;
1055	<	}
1056	<	offset += comp_stamps[i]->getNAtoms();
1057	<
1058	<	#ifdef IS_MPI
1059	<	}
1060	<	#endif //is_mpi
1061	<
1062	<	molIndex++;
1063	<	}
1064	<	}
1065	<
1066	<	the_ff->initializeBonds( the_bonds );
1067	<	}
1068	<
1069	<	void SimSetup::makeBends( void ){
1070	<
1071	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1072	<	bend_set* the_bends;
1073	<	BendStamp* current_bend;
1074	<	LinkedAssign* extras;
1075	<	LinkedAssign* current_extra;
1076	<
1077	<
1078	<	the_bends = new bend_set[tot_bends];
1079	<	index = 0;
1080	<	offset = 0;
1081	<	molIndex = 0;
1082	<	for( i=0; i<n_components; i++ ){
1083	<
1084	<	for( j=0; j<components_nmol[i]; j++ ){
1085	<
1086	<	#ifdef IS_MPI
1087	<	if( mpiSim->getMyMolStart() <= molIndex &&
1088	<	molIndex <= mpiSim->getMyMolEnd() ){
1089	<	#endif // is_mpi
1090	<
1091	<	for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
1092	<
1093	<	current_bend = comp_stamps[i]->getBend( k );
1094	<	the_bends[index].a = current_bend->getA() + offset;
1095	<	the_bends[index].b = current_bend->getB() + offset;
1096	<	the_bends[index].c = current_bend->getC() + offset;
1097	<
1098	<	if( current_bend->haveExtras() ){
1099	<
1100	<	extras = current_bend->getExtras();
1101	<	current_extra = extras;
1102	<
1103	<	while( current_extra != NULL ){
1104	<	if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
1105	<
1106	<	switch( current_extra->getType() ){
1107	<
1108	<	case 0:
1109	<	the_bends[index].ghost =
1110	<	current_extra->getInt() + offset;
1111	<	the_bends[index].isGhost = 1;
1112	<	break;
1113	<
1114	<	case 1:
1115	<	the_bends[index].ghost =
1116	<	(int)current_extra->getDouble() + offset;
1117	<	the_bends[index].isGhost = 1;
1118	<	break;
1119	<
1120	<	default:
1121	<	sprintf( painCave.errMsg,
1122	<	"SimSetup Error: ghostVectorSource was neiter a "
1123	<	"double nor an int.\n"
1124	<	"-->Bend[%d] in %s\n",
1125	<	k, comp_stamps[i]->getID() );
1126	<	painCave.isFatal = 1;
1127	<	simError();
1128	<	}
1129	<	}
1130	<
1131	<	else{
1132	<
1133	<	sprintf( painCave.errMsg,
1134	<	"SimSetup Error: unhandled bend assignment:\n"
1135	<	" -->%s in Bend[%d] in %s\n",
1136	<	current_extra->getlhs(),
1137	<	k, comp_stamps[i]->getID() );
1138	<	painCave.isFatal = 1;
1139	<	simError();
1140	<	}
1141	<
1142	<	current_extra = current_extra->getNext();
1143	<	}
1144	<	}
1145	<
1146	<	if( !the_bends[index].isGhost ){
1147	<
1148	<	exI = the_bends[index].a;
1149	<	exJ = the_bends[index].c;
1150	<	}
1151	<	else{
1152	<
1153	<	exI = the_bends[index].a;
1154	<	exJ = the_bends[index].b;
1155	<	}
1156	<
1157	<	// exclude_I must always be the smaller of the pair
1158	<	if( exI > exJ ){
1159	<	tempEx = exI;
1160	<	exI = exJ;
1161	<	exJ = tempEx;
1162	<	}
1163	<
1164	<
1165	<	#ifdef IS_MPI
1166	<
1167	<	the_excludes[(index + tot_bonds)*2] =
1168	<	the_atoms[exI]->getGlobalIndex() + 1;
1169	<	the_excludes[(index + tot_bonds)*2 + 1] =
1170	<	the_atoms[exJ]->getGlobalIndex() + 1;
1171	<
1172	<	#else // isn't MPI
1173	<
1174	<	the_excludes[(index + tot_bonds)*2] = exI + 1;
1175	<	the_excludes[(index + tot_bonds)*2 + 1] = exJ + 1;
1176	<	// fortran index from 1 (hence the +1 in the indexing)
1177	<	#endif //is_mpi
1178	<
1179	<
1180	<	// increment the index and repeat;
1181	<	index++;
1182	<	}
1183	<	offset += comp_stamps[i]->getNAtoms();
1184	<
1185	<	#ifdef IS_MPI
1186	<	}
1187	<	#endif //is_mpi
1188	<
1189	<	molIndex++;
1190	<	}
1239	>	atomOffset += info.nAtoms;
1240	>	delete[] theBonds;
1241	>	delete[] theBends;
1242	>	delete[] theTorsions;
1243		}
1244
1245		#ifdef IS_MPI
1246	<	sprintf( checkPointMsg,
1195	<	"Successfully created the bends list.\n" );
1246	>	sprintf( checkPointMsg, "all molecules initialized succesfully" );
1247		MPIcheckPoint();
1248		#endif // is_mpi
1198	–
1249
1250	<	the_ff->initializeBends( the_bends );
1251	<	}
1250	>	// clean up the forcefield
1251	>	the_ff->calcRcut();
1252	>	the_ff->cleanMe();
1253
1203	–	void SimSetup::makeTorsions( void ){
1204	–
1205	–	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1206	–	torsion_set* the_torsions;
1207	–	TorsionStamp* current_torsion;
1208	–
1209	–	the_torsions = new torsion_set[tot_torsions];
1210	–	index = 0;
1211	–	offset = 0;
1212	–	molIndex = 0;
1213	–	for( i=0; i<n_components; i++ ){
1214	–
1215	–	for( j=0; j<components_nmol[i]; j++ ){
1216	–
1217	–	#ifdef IS_MPI
1218	–	if( mpiSim->getMyMolStart() <= molIndex &&
1219	–	molIndex <= mpiSim->getMyMolEnd() ){
1220	–	#endif // is_mpi
1221	–
1222	–	for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
1223	–
1224	–	current_torsion = comp_stamps[i]->getTorsion( k );
1225	–	the_torsions[index].a = current_torsion->getA() + offset;
1226	–	the_torsions[index].b = current_torsion->getB() + offset;
1227	–	the_torsions[index].c = current_torsion->getC() + offset;
1228	–	the_torsions[index].d = current_torsion->getD() + offset;
1229	–
1230	–	exI = the_torsions[index].a;
1231	–	exJ = the_torsions[index].d;
1232	–
1233	–
1234	–	// exclude_I must always be the smaller of the pair
1235	–	if( exI > exJ ){
1236	–	tempEx = exI;
1237	–	exI = exJ;
1238	–	exJ = tempEx;
1239	–	}
1240	–
1241	–
1242	–	#ifdef IS_MPI
1243	–
1244	–	the_excludes[(index + tot_bonds + tot_bends)*2] =
1245	–	the_atoms[exI]->getGlobalIndex() + 1;
1246	–	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] =
1247	–	the_atoms[exJ]->getGlobalIndex() + 1;
1248	–
1249	–	#else // isn't MPI
1250	–
1251	–	the_excludes[(index + tot_bonds + tot_bends)*2] = exI + 1;
1252	–	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] = exJ + 1;
1253	–	// fortran indexes from 1 (hence the +1 in the indexing)
1254	–	#endif //is_mpi
1255	–
1256	–
1257	–	// increment the index and repeat;
1258	–	index++;
1259	–	}
1260	–	offset += comp_stamps[i]->getNAtoms();
1261	–
1262	–	#ifdef IS_MPI
1263	–	}
1264	–	#endif //is_mpi
1265	–
1266	–	molIndex++;
1267	–	}
1268	–	}
1269	–
1270	–	the_ff->initializeTorsions( the_torsions );
1254		}
1255
1256		void SimSetup::initFromBass( void ){
#	Line 1289 \| Line 1272 \| void SimSetup::initFromBass( void ){
1272		have_extra =1;
1273
1274		n_cells = (int)temp3 - 1;
1275	<	cellx = simnfo->box_x / temp3;
1276	<	celly = simnfo->box_y / temp3;
1277	<	cellz = simnfo->box_z / temp3;
1275	>	cellx = simnfo->boxLx / temp3;
1276	>	celly = simnfo->boxLy / temp3;
1277	>	cellz = simnfo->boxLz / temp3;
1278		n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
1279		temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
1280		n_per_extra = (int)ceil( temp1 );
#	Line 1306 \| Line 1289 \| void SimSetup::initFromBass( void ){
1289		}
1290		else{
1291		n_cells = (int)temp3;
1292	<	cellx = simnfo->box_x / temp3;
1293	<	celly = simnfo->box_y / temp3;
1294	<	cellz = simnfo->box_z / temp3;
1292	>	cellx = simnfo->boxLx / temp3;
1293	>	celly = simnfo->boxLy / temp3;
1294	>	cellz = simnfo->boxLz / temp3;
1295		}
1296
1297		current_mol = 0;

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents): Revision 414 by mmeineke, Wed Mar 26 22:02:36 2003 UTC vs. Revision 605 by gezelter, Tue Jul 15 03:27:24 2003 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 414 by mmeineke, Wed Mar 26 22:02:36 2003 UTC vs.
Revision 605 by gezelter, Tue Jul 15 03:27:24 2003 UTC