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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 412 by mmeineke, Wed Mar 26 21:50:33 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	<	if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
204	<	else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
205	<	else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
206	<	else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
207	<	else{
203	>	case FF_LJ:
204	>	the_ff = new LJFF();
205	>	break;
206	>
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
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	<	// new AllLong( simnfo );
793	>	switch( ensembleCase ){
794
795	<	if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
796	<	if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo, the_ff );
797	<	if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo, the_ff );
798	<	if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
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 897 \| Line 1225 \| void SimSetup::makeMolecules( void ){
1225		excludeOffset += info.nTorsions;
1226
1227
1228	<
901	<
902	<
903	<
904	<
905	<
906	<	void SimSetup::makeAtoms( void ){
907	<
908	<	int i, j, k, index;
909	<	double ux, uy, uz, uSqr, u;
910	<	AtomStamp* current_atom;
1228	>	// send the arrays off to the forceField for init.
1229
1230	<	DirectionalAtom* dAtom;
1231	<	int molIndex, molStart, molEnd, nMemb, lMolIndex;
1232	<
1233	<	lMolIndex = 0;
916	<	molIndex = 0;
917	<	index = 0;
918	<	for( i=0; i<n_components; i++ ){
919	<
920	<	for( j=0; j<components_nmol[i]; j++ ){
921	<
922	<	#ifdef IS_MPI
923	<	if( mpiSim->getMyMolStart() <= molIndex &&
924	<	molIndex <= mpiSim->getMyMolEnd() ){
925	<	#endif // is_mpi
926	<
927	<	molStart = index;
928	<	nMemb = comp_stamps[i]->getNAtoms();
929	<	for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
930	<
931	<	current_atom = comp_stamps[i]->getAtom( k );
932	<	if( current_atom->haveOrientation() ){
933	<
934	<	dAtom = new DirectionalAtom(index);
935	<	simnfo->n_oriented++;
936	<	the_atoms[index] = dAtom;
937	<
938	<	ux = current_atom->getOrntX();
939	<	uy = current_atom->getOrntY();
940	<	uz = current_atom->getOrntZ();
941	<
942	<	uSqr = (ux * ux) + (uy * uy) + (uz * uz);
943	<
944	<	u = sqrt( uSqr );
945	<	ux = ux / u;
946	<	uy = uy / u;
947	<	uz = uz / u;
948	<
949	<	dAtom->setSUx( ux );
950	<	dAtom->setSUy( uy );
951	<	dAtom->setSUz( uz );
952	<	}
953	<	else{
954	<	the_atoms[index] = new GeneralAtom(index);
955	<	}
956	<	the_atoms[index]->setType( current_atom->getType() );
957	<	the_atoms[index]->setIndex( index );
958	<
959	<	// increment the index and repeat;
960	<	index++;
961	<	}
962	<
963	<	molEnd = index -1;
964	<	the_molecules[lMolIndex].setNMembers( nMemb );
965	<	the_molecules[lMolIndex].setStartAtom( molStart );
966	<	the_molecules[lMolIndex].setEndAtom( molEnd );
967	<	the_molecules[lMolIndex].setStampID( i );
968	<	lMolIndex++;
1230	>	the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
1231	>	the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
1232	>	the_ff->initializeBends( info.nBends, info.myBends, theBends );
1233	>	the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
1234
970	–	#ifdef IS_MPI
971	–	}
972	–	#endif //is_mpi
973	–
974	–	molIndex++;
975	–	}
976	–	}
1235
1236	<	#ifdef IS_MPI
979	<	for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] );
980	<
981	<	delete[] globalIndex;
982	<
983	<	mpiSim->mpiRefresh();
984	<	#endif //IS_MPI
985	<
986	<	the_ff->initializeAtoms();
987	<	}
988	<
989	<	void SimSetup::makeBonds( void ){
990	<
991	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
992	<	bond_pair* the_bonds;
993	<	BondStamp* current_bond;
994	<
995	<	the_bonds = new bond_pair[tot_bonds];
996	<	index = 0;
997	<	offset = 0;
998	<	molIndex = 0;
999	<
1000	<	for( i=0; i<n_components; i++ ){
1001	<
1002	<	for( j=0; j<components_nmol[i]; j++ ){
1003	<
1004	<	#ifdef IS_MPI
1005	<	if( mpiSim->getMyMolStart() <= molIndex &&
1006	<	molIndex <= mpiSim->getMyMolEnd() ){
1007	<	#endif // is_mpi
1008	<
1009	<	for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
1010	<
1011	<	current_bond = comp_stamps[i]->getBond( k );
1012	<	the_bonds[index].a = current_bond->getA() + offset;
1013	<	the_bonds[index].b = current_bond->getB() + offset;
1014	<
1015	<	exI = the_bonds[index].a;
1016	<	exJ = the_bonds[index].b;
1017	<
1018	<	// exclude_I must always be the smaller of the pair
1019	<	if( exI > exJ ){
1020	<	tempEx = exI;
1021	<	exI = exJ;
1022	<	exJ = tempEx;
1023	<	}
1024	<
1025	<
1026	<	#ifdef IS_MPI
1027	<
1028	<	the_excludes[index*2] =
1029	<	the_atoms[exI]->getGlobalIndex() + 1;
1030	<	the_excludes[index*2 + 1] =
1031	<	the_atoms[exJ]->getGlobalIndex() + 1;
1032	<
1033	<	#else // isn't MPI
1034	<
1035	<	the_excludes[index*2] = exI + 1;
1036	<	the_excludes[index*2 + 1] = exJ + 1;
1037	<	// fortran index from 1 (hence the +1 in the indexing)
1038	<	#endif //is_mpi
1039	<
1040	<	// increment the index and repeat;
1041	<	index++;
1042	<	}
1043	<	offset += comp_stamps[i]->getNAtoms();
1044	<
1045	<	#ifdef IS_MPI
1046	<	}
1047	<	#endif //is_mpi
1048	<
1049	<	molIndex++;
1050	<	}
1051	<	}
1052	<
1053	<	the_ff->initializeBonds( the_bonds );
1054	<	}
1055	<
1056	<	void SimSetup::makeBends( void ){
1057	<
1058	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1059	<	bend_set* the_bends;
1060	<	BendStamp* current_bend;
1061	<	LinkedAssign* extras;
1062	<	LinkedAssign* current_extra;
1063	<
1236	>	the_molecules[i].initialize( info );
1237
1065	–	the_bends = new bend_set[tot_bends];
1066	–	index = 0;
1067	–	offset = 0;
1068	–	molIndex = 0;
1069	–	for( i=0; i<n_components; i++ ){
1238
1239	<	for( j=0; j<components_nmol[i]; j++ ){
1240	<
1241	<	#ifdef IS_MPI
1242	<	if( mpiSim->getMyMolStart() <= molIndex &&
1075	<	molIndex <= mpiSim->getMyMolEnd() ){
1076	<	#endif // is_mpi
1077	<
1078	<	for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
1079	<
1080	<	current_bend = comp_stamps[i]->getBend( k );
1081	<	the_bends[index].a = current_bend->getA() + offset;
1082	<	the_bends[index].b = current_bend->getB() + offset;
1083	<	the_bends[index].c = current_bend->getC() + offset;
1084	<
1085	<	if( current_bend->haveExtras() ){
1086	<
1087	<	extras = current_bend->getExtras();
1088	<	current_extra = extras;
1089	<
1090	<	while( current_extra != NULL ){
1091	<	if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
1092	<
1093	<	switch( current_extra->getType() ){
1094	<
1095	<	case 0:
1096	<	the_bends[index].ghost =
1097	<	current_extra->getInt() + offset;
1098	<	the_bends[index].isGhost = 1;
1099	<	break;
1100	<
1101	<	case 1:
1102	<	the_bends[index].ghost =
1103	<	(int)current_extra->getDouble() + offset;
1104	<	the_bends[index].isGhost = 1;
1105	<	break;
1106	<
1107	<	default:
1108	<	sprintf( painCave.errMsg,
1109	<	"SimSetup Error: ghostVectorSource was neiter a "
1110	<	"double nor an int.\n"
1111	<	"-->Bend[%d] in %s\n",
1112	<	k, comp_stamps[i]->getID() );
1113	<	painCave.isFatal = 1;
1114	<	simError();
1115	<	}
1116	<	}
1117	<
1118	<	else{
1119	<
1120	<	sprintf( painCave.errMsg,
1121	<	"SimSetup Error: unhandled bend assignment:\n"
1122	<	" -->%s in Bend[%d] in %s\n",
1123	<	current_extra->getlhs(),
1124	<	k, comp_stamps[i]->getID() );
1125	<	painCave.isFatal = 1;
1126	<	simError();
1127	<	}
1128	<
1129	<	current_extra = current_extra->getNext();
1130	<	}
1131	<	}
1132	<
1133	<	if( !the_bends[index].isGhost ){
1134	<
1135	<	exI = the_bends[index].a;
1136	<	exJ = the_bends[index].c;
1137	<	}
1138	<	else{
1139	<
1140	<	exI = the_bends[index].a;
1141	<	exJ = the_bends[index].b;
1142	<	}
1143	<
1144	<	// exclude_I must always be the smaller of the pair
1145	<	if( exI > exJ ){
1146	<	tempEx = exI;
1147	<	exI = exJ;
1148	<	exJ = tempEx;
1149	<	}
1150	<
1151	<
1152	<	#ifdef IS_MPI
1153	<
1154	<	the_excludes[(index + tot_bonds)*2] =
1155	<	the_atoms[exI]->getGlobalIndex() + 1;
1156	<	the_excludes[(index + tot_bonds)*2 + 1] =
1157	<	the_atoms[exJ]->getGlobalIndex() + 1;
1158	<
1159	<	#else // isn't MPI
1160	<
1161	<	the_excludes[(index + tot_bonds)*2] = exI + 1;
1162	<	the_excludes[(index + tot_bonds)*2 + 1] = exJ + 1;
1163	<	// fortran index from 1 (hence the +1 in the indexing)
1164	<	#endif //is_mpi
1165	<
1166	<
1167	<	// increment the index and repeat;
1168	<	index++;
1169	<	}
1170	<	offset += comp_stamps[i]->getNAtoms();
1171	<
1172	<	#ifdef IS_MPI
1173	<	}
1174	<	#endif //is_mpi
1175	<
1176	<	molIndex++;
1177	<	}
1239	>	atomOffset += info.nAtoms;
1240	>	delete[] theBonds;
1241	>	delete[] theBends;
1242	>	delete[] theTorsions;
1243		}
1244
1245		#ifdef IS_MPI
1246	<	sprintf( checkPointMsg,
1182	<	"Successfully created the bends list.\n" );
1246	>	sprintf( checkPointMsg, "all molecules initialized succesfully" );
1247		MPIcheckPoint();
1248		#endif // is_mpi
1185	–
1249
1250	<	the_ff->initializeBends( the_bends );
1251	<	}
1250	>	// clean up the forcefield
1251	>	the_ff->calcRcut();
1252	>	the_ff->cleanMe();
1253
1190	–	void SimSetup::makeTorsions( void ){
1191	–
1192	–	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1193	–	torsion_set* the_torsions;
1194	–	TorsionStamp* current_torsion;
1195	–
1196	–	the_torsions = new torsion_set[tot_torsions];
1197	–	index = 0;
1198	–	offset = 0;
1199	–	molIndex = 0;
1200	–	for( i=0; i<n_components; i++ ){
1201	–
1202	–	for( j=0; j<components_nmol[i]; j++ ){
1203	–
1204	–	#ifdef IS_MPI
1205	–	if( mpiSim->getMyMolStart() <= molIndex &&
1206	–	molIndex <= mpiSim->getMyMolEnd() ){
1207	–	#endif // is_mpi
1208	–
1209	–	for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
1210	–
1211	–	current_torsion = comp_stamps[i]->getTorsion( k );
1212	–	the_torsions[index].a = current_torsion->getA() + offset;
1213	–	the_torsions[index].b = current_torsion->getB() + offset;
1214	–	the_torsions[index].c = current_torsion->getC() + offset;
1215	–	the_torsions[index].d = current_torsion->getD() + offset;
1216	–
1217	–	exI = the_torsions[index].a;
1218	–	exJ = the_torsions[index].d;
1219	–
1220	–
1221	–	// exclude_I must always be the smaller of the pair
1222	–	if( exI > exJ ){
1223	–	tempEx = exI;
1224	–	exI = exJ;
1225	–	exJ = tempEx;
1226	–	}
1227	–
1228	–
1229	–	#ifdef IS_MPI
1230	–
1231	–	the_excludes[(index + tot_bonds + tot_bends)*2] =
1232	–	the_atoms[exI]->getGlobalIndex() + 1;
1233	–	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] =
1234	–	the_atoms[exJ]->getGlobalIndex() + 1;
1235	–
1236	–	#else // isn't MPI
1237	–
1238	–	the_excludes[(index + tot_bonds + tot_bends)*2] = exI + 1;
1239	–	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] = exJ + 1;
1240	–	// fortran indexes from 1 (hence the +1 in the indexing)
1241	–	#endif //is_mpi
1242	–
1243	–
1244	–	// increment the index and repeat;
1245	–	index++;
1246	–	}
1247	–	offset += comp_stamps[i]->getNAtoms();
1248	–
1249	–	#ifdef IS_MPI
1250	–	}
1251	–	#endif //is_mpi
1252	–
1253	–	molIndex++;
1254	–	}
1255	–	}
1256	–
1257	–	the_ff->initializeTorsions( the_torsions );
1254		}
1255
1256		void SimSetup::initFromBass( void ){
#	Line 1276 \| 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 1293 \| 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 412 by mmeineke, Wed Mar 26 21:50:33 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 412 by mmeineke, Wed Mar 26 21:50:33 2003 UTC vs.
Revision 605 by gezelter, Tue Jul 15 03:27:24 2003 UTC