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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 394 by gezelter, Mon Mar 24 21:55:34 2003 UTC vs.
Revision 582 by mmeineke, Wed Jul 9 15:33:46 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	+
22	+
23	+	#define FF_DUFF 0
24	+	#define FF_LJ 1
25	+
26	+
27		SimSetup::SimSetup(){
28		stamps = new MakeStamps();
29		globals = new Globals();
#	Line 66 \| Line 78 \| void SimSetup::createSim( void ){
78
79		MakeStamps *the_stamps;
80		Globals* the_globals;
81	<	int i, j;
81	>	int i, j, k, globalAtomIndex;
82	>
83	>	int ensembleCase;
84	>	int ffCase;
85	>
86	>	ensembleCase = -1;
87	>	ffCase = -1;
88
89		// get the stamps and globals;
90		the_stamps = stamps;
#	Line 80 \| Line 98 \| void SimSetup::createSim( void ){
98		// get the ones we know are there, yet still may need some work.
99		n_components = the_globals->getNComponents();
100		strcpy( force_field, the_globals->getForceField() );
101	+
102	+	if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
103	+	else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
104	+	else{
105	+	sprintf( painCave.errMsg,
106	+	"SimSetup Error. Unrecognized force field -> %s\n",
107	+	force_field );
108	+	painCave.isFatal = 1;
109	+	simError();
110	+	}
111	+
112	+	// get the ensemble:
113		strcpy( ensemble, the_globals->getEnsemble() );
114	+
115	+	if( !strcasecmp( ensemble, "NVE" )) ensembleCase = NVE_ENS;
116	+	else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
117	+	else if( !strcasecmp( ensemble, "NPTi" ) \|\| !strcasecmp( ensemble, "NPT") )
118	+	ensembleCase = NPTi_ENS;
119	+	else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
120	+	else{
121	+	sprintf( painCave.errMsg,
122	+	"SimSetup Warning. Unrecognized Ensemble -> %s, "
123	+	"reverting to NVE for this simulation.\n",
124	+	ensemble );
125	+	painCave.isFatal = 0;
126	+	simError();
127	+	strcpy( ensemble, "NVE" );
128	+	ensembleCase = NVE_ENS;
129	+	}
130		strcpy( simnfo->ensemble, ensemble );
131
132	+
133	+	// if( !strcasecmp( ensemble, "NPT" ) ) {
134	+	// the_extendedsystem = new ExtendedSystem( simnfo );
135	+	// the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
136	+	// if (the_globals->haveTargetPressure())
137	+	// the_extendedsystem->setTargetPressure(the_globals->getTargetPressure());
138	+	// else {
139	+	// sprintf( painCave.errMsg,
140	+	// "SimSetup error: If you use the constant pressure\n"
141	+	// " ensemble, you must set targetPressure.\n"
142	+	// " This was found in the BASS file.\n");
143	+	// painCave.isFatal = 1;
144	+	// simError();
145	+	// }
146	+
147	+	// if (the_globals->haveTauThermostat())
148	+	// the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
149	+	// else if (the_globals->haveQmass())
150	+	// the_extendedsystem->setQmass(the_globals->getQmass());
151	+	// else {
152	+	// sprintf( painCave.errMsg,
153	+	// "SimSetup error: If you use one of the constant temperature\n"
154	+	// " ensembles, you must set either tauThermostat or qMass.\n"
155	+	// " Neither of these was found in the BASS file.\n");
156	+	// painCave.isFatal = 1;
157	+	// simError();
158	+	// }
159	+
160	+	// if (the_globals->haveTauBarostat())
161	+	// the_extendedsystem->setTauBarostat(the_globals->getTauBarostat());
162	+	// else {
163	+	// sprintf( painCave.errMsg,
164	+	// "SimSetup error: If you use the constant pressure\n"
165	+	// " ensemble, you must set tauBarostat.\n"
166	+	// " This was found in the BASS file.\n");
167	+	// painCave.isFatal = 1;
168	+	// simError();
169	+	// }
170	+
171	+	// } else if ( !strcasecmp( ensemble, "NVT") ) {
172	+	// the_extendedsystem = new ExtendedSystem( simnfo );
173	+	// the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
174	+
175	+	// if (the_globals->haveTauThermostat())
176	+	// the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
177	+	// else if (the_globals->haveQmass())
178	+	// the_extendedsystem->setQmass(the_globals->getQmass());
179	+	// else {
180	+	// sprintf( painCave.errMsg,
181	+	// "SimSetup error: If you use one of the constant temperature\n"
182	+	// " ensembles, you must set either tauThermostat or qMass.\n"
183	+	// " Neither of these was found in the BASS file.\n");
184	+	// painCave.isFatal = 1;
185	+	// simError();
186	+	// }
187	+
188		strcpy( simnfo->mixingRule, the_globals->getMixingRule() );
189		simnfo->usePBC = the_globals->getPBC();
190
191	+	int usesDipoles = 0;
192	+	switch( ffCase ){
193
194	+	case FF_DUFF:
195	+	the_ff = new DUFF();
196	+	usesDipoles = 1;
197	+	break;
198
199	<	if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
200	<	else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
201	<	else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
202	<	else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
203	<	else{
199	>	case FF_LJ:
200	>	the_ff = new LJFF();
201	>	break;
202	>
203	>	default:
204		sprintf( painCave.errMsg,
205	<	"SimSetup Error. Unrecognized force field -> %s\n",
98	<	force_field );
205	>	"SimSetup Error. Unrecognized force field in case statement.\n");
206		painCave.isFatal = 1;
207		simError();
208		}
#	Line 105 \| Line 212 \| void SimSetup::createSim( void ){
212		MPIcheckPoint();
213		#endif // is_mpi
214
108	–
109	–
215		// get the components and calculate the tot_nMol and indvidual n_mol
216		the_components = the_globals->getComponents();
217		components_nmol = new int[n_components];
#	Line 230 \| Line 335 \| void SimSetup::createSim( void ){
335		simnfo->n_torsions = tot_torsions;
336		simnfo->n_SRI = tot_SRI;
337		simnfo->n_mol = tot_nmol;
233	–
338
339	+	simnfo->molMembershipArray = new int[tot_atoms];
340	+
341		#ifdef IS_MPI
342
343		// divide the molecules among processors here.
344
345		mpiSim = new mpiSimulation( simnfo );
346
241	–
242	–
347		globalIndex = mpiSim->divideLabor();
348
245	–
246	–
349		// set up the local variables
350
351		int localMol, allMol;
352		int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
353	+
354	+	int* mol2proc = mpiSim->getMolToProcMap();
355	+	int* molCompType = mpiSim->getMolComponentType();
356
357		allMol = 0;
358		localMol = 0;
#	Line 255 \| Line 360 \| void SimSetup::createSim( void ){
360		local_bonds = 0;
361		local_bends = 0;
362		local_torsions = 0;
363	+	globalAtomIndex = 0;
364	+
365	+
366		for( i=0; i<n_components; i++ ){
367
368		for( j=0; j<components_nmol[i]; j++ ){
369
370	<	if( mpiSim->getMyMolStart() <= allMol &&
263	<	allMol <= mpiSim->getMyMolEnd() ){
370	>	if( mol2proc[allMol] == worldRank ){
371
372		local_atoms += comp_stamps[i]->getNAtoms();
373		local_bonds += comp_stamps[i]->getNBonds();
#	Line 268 \| Line 375 \| void SimSetup::createSim( void ){
375		local_torsions += comp_stamps[i]->getNTorsions();
376		localMol++;
377		}
378	<	allMol++;
378	>	for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
379	>	simnfo->molMembershipArray[globalAtomIndex] = allMol;
380	>	globalAtomIndex++;
381	>	}
382	>
383	>	allMol++;
384		}
385		}
386		local_SRI = local_bonds + local_bends + local_torsions;
387
276	–
388		simnfo->n_atoms = mpiSim->getMyNlocal();
389
390		if( local_atoms != simnfo->n_atoms ){
391		sprintf( painCave.errMsg,
392		"SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
393	<	" localAtom (%d) are note equal.\n",
393	>	" localAtom (%d) are not equal.\n",
394		simnfo->n_atoms,
395		local_atoms );
396		painCave.isFatal = 1;
#	Line 304 \| Line 415 \| void SimSetup::createSim( void ){
415		Atom::createArrays(simnfo->n_atoms);
416		the_atoms = new Atom*[simnfo->n_atoms];
417		the_molecules = new Molecule[simnfo->n_mol];
418	+	int molIndex;
419
420	+	// initialize the molecule's stampID's
421
422	<	if( simnfo->n_SRI ){
423	<	the_sris = new SRI*[simnfo->n_SRI];
424	<	the_excludes = new int[2 * simnfo->n_SRI];
422	>	#ifdef IS_MPI
423	>
424	>
425	>	molIndex = 0;
426	>	for(i=0; i<mpiSim->getTotNmol(); i++){
427	>
428	>	if(mol2proc[i] == worldRank ){
429	>	the_molecules[molIndex].setStampID( molCompType[i] );
430	>	the_molecules[molIndex].setMyIndex( molIndex );
431	>	the_molecules[molIndex].setGlobalIndex( i );
432	>	molIndex++;
433	>	}
434	>	}
435	>
436	>	#else // is_mpi
437	>
438	>	molIndex = 0;
439	>	globalAtomIndex = 0;
440	>	for(i=0; i<n_components; i++){
441	>	for(j=0; j<components_nmol[i]; j++ ){
442	>	the_molecules[molIndex].setStampID( i );
443	>	the_molecules[molIndex].setMyIndex( molIndex );
444	>	the_molecules[molIndex].setGlobalIndex( molIndex );
445	>	for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
446	>	simnfo->molMembershipArray[globalAtomIndex] = molIndex;
447	>	globalAtomIndex++;
448	>	}
449	>	molIndex++;
450	>	}
451	>	}
452	>
453	>
454	>	#endif // is_mpi
455	>
456	>
457	>	if( simnfo->n_SRI ){
458	>
459	>	Exclude::createArray(simnfo->n_SRI);
460	>	the_excludes = new Exclude*[simnfo->n_SRI];
461	>	for( int ex=0; ex<simnfo->n_SRI; ex++) the_excludes[ex] = new Exclude(ex);
462		simnfo->globalExcludes = new int;
463	<	simnfo->n_exclude = tot_SRI;
463	>	simnfo->n_exclude = simnfo->n_SRI;
464		}
465		else{
466
467	<	the_excludes = new int[2];
468	<	the_excludes[0] = 0;
469	<	the_excludes[1] = 0;
467	>	Exclude::createArray( 1 );
468	>	the_excludes = new Exclude*;
469	>	the_excludes[0] = new Exclude(0);
470	>	the_excludes[0]->setPair( 0,0 );
471		simnfo->globalExcludes = new int;
472		simnfo->globalExcludes[0] = 0;
473	<
323	<	simnfo->n_exclude = 1;
473	>	simnfo->n_exclude = 0;
474		}
475
476		// set the arrays into the SimInfo object
477
478		simnfo->atoms = the_atoms;
479	<	simnfo->sr_interactions = the_sris;
479	>	simnfo->molecules = the_molecules;
480		simnfo->nGlobalExcludes = 0;
481		simnfo->excludes = the_excludes;
482
483
484		// get some of the tricky things that may still be in the globals
485
486	<
486	>	double boxVector[3];
487		if( the_globals->haveBox() ){
488	<	simnfo->box_x = the_globals->getBox();
489	<	simnfo->box_y = the_globals->getBox();
490	<	simnfo->box_z = the_globals->getBox();
488	>	boxVector[0] = the_globals->getBox();
489	>	boxVector[1] = the_globals->getBox();
490	>	boxVector[2] = the_globals->getBox();
491	>
492	>	simnfo->setBox( boxVector );
493		}
494		else if( the_globals->haveDensity() ){
495
496		double vol;
497		vol = (double)tot_nmol / the_globals->getDensity();
498	<	simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
499	<	simnfo->box_y = simnfo->box_x;
500	<	simnfo->box_z = simnfo->box_x;
498	>	boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
499	>	boxVector[1] = boxVector[0];
500	>	boxVector[2] = boxVector[0];
501	>
502	>	simnfo->setBox( boxVector );
503		}
504		else{
505		if( !the_globals->haveBoxX() ){
#	Line 354 \| Line 508 \| void SimSetup::createSim( void ){
508		painCave.isFatal = 1;
509		simError();
510		}
511	<	simnfo->box_x = the_globals->getBoxX();
511	>	boxVector[0] = the_globals->getBoxX();
512
513		if( !the_globals->haveBoxY() ){
514		sprintf( painCave.errMsg,
#	Line 362 \| Line 516 \| void SimSetup::createSim( void ){
516		painCave.isFatal = 1;
517		simError();
518		}
519	<	simnfo->box_y = the_globals->getBoxY();
519	>	boxVector[1] = the_globals->getBoxY();
520
521		if( !the_globals->haveBoxZ() ){
522		sprintf( painCave.errMsg,
#	Line 370 \| Line 524 \| void SimSetup::createSim( void ){
524		painCave.isFatal = 1;
525		simError();
526		}
527	<	simnfo->box_z = the_globals->getBoxZ();
527	>	boxVector[2] = the_globals->getBoxZ();
528	>
529	>	simnfo->setBox( boxVector );
530		}
531
532		#ifdef IS_MPI
#	Line 383 \| Line 539 \| void SimSetup::createSim( void ){
539
540		the_ff->setSimInfo( simnfo );
541
542	<	makeAtoms();
542	>	makeMolecules();
543		simnfo->identArray = new int[simnfo->n_atoms];
544		for(i=0; i<simnfo->n_atoms; i++){
545		simnfo->identArray[i] = the_atoms[i]->getIdent();
546		}
547
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	–
548		if (the_globals->getUseRF() ) {
549		simnfo->useReactionField = 1;
550
#	Line 413 \| Line 556 \| void SimSetup::createSim( void ){
556		painCave.isFatal = 0;
557		simError();
558		double smallest;
559	<	smallest = simnfo->box_x;
560	<	if (simnfo->box_y <= smallest) smallest = simnfo->box_y;
561	<	if (simnfo->box_z <= smallest) smallest = simnfo->box_z;
559	>	smallest = simnfo->boxLx;
560	>	if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy;
561	>	if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz;
562		simnfo->ecr = 0.5 * smallest;
563		} else {
564		simnfo->ecr = the_globals->getECR();
#	Line 443 \| Line 586 \| void SimSetup::createSim( void ){
586		}
587		simnfo->dielectric = the_globals->getDielectric();
588		} else {
589	<	if (simnfo->n_dipoles) {
589	>	if (usesDipoles) {
590
591		if( !the_globals->haveECR() ){
592		sprintf( painCave.errMsg,
593	<	"SimSetup Warning: using default value of 1/2 the smallest"
593	>	"SimSetup Warning: using default value of 1/2 the smallest "
594		"box length for the electrostaticCutoffRadius.\n"
595		"I hope you have a very fast processor!\n");
596		painCave.isFatal = 0;
597		simError();
598		double smallest;
599	<	smallest = simnfo->box_x;
600	<	if (simnfo->box_y <= smallest) smallest = simnfo->box_y;
601	<	if (simnfo->box_z <= smallest) smallest = simnfo->box_z;
599	>	smallest = simnfo->boxLx;
600	>	if (simnfo->boxLy <= smallest) smallest = simnfo->boxLy;
601	>	if (simnfo->boxLz <= smallest) smallest = simnfo->boxLz;
602		simnfo->ecr = 0.5 * smallest;
603		} else {
604		simnfo->ecr = the_globals->getECR();
#	Line 463 \| Line 606 \| void SimSetup::createSim( void ){
606
607		if( !the_globals->haveEST() ){
608		sprintf( painCave.errMsg,
609	<	"SimSetup Warning: using default value of 5% of the"
609	>	"SimSetup Warning: using default value of 5%% of the "
610		"electrostaticCutoffRadius for the "
611		"electrostaticSkinThickness\n"
612		);
#	Line 634 \| Line 777 \| void SimSetup::createSim( void ){
777		if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
778
779
780	<	// // make the longe range forces and the integrator
780	>	// make the integrator
781	>
782	>
783	>	NVT* myNVT = NULL;
784	>	switch( ensembleCase ){
785
786	<	// new AllLong( simnfo );
786	>	case NVE_ENS:
787	>	new NVE( simnfo, the_ff );
788	>	break;
789
790	<	if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
791	<	if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo, the_ff );
792	<	if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo, the_ff );
793	<	if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
790	>	case NVT_ENS:
791	>	myNVT = new NVT( simnfo, the_ff );
792	>	myNVT->setTargetTemp(the_globals->getTargetTemp());
793	>
794	>	if (the_globals->haveTauThermostat())
795	>	myNVT->setTauThermostat(the_globals->getTauThermostat());
796	>	// else if (the_globals->haveQmass())
797	>	// myNVT->setQmass(the_globals->getQmass());
798	>	else {
799	>	sprintf( painCave.errMsg,
800	>	"SimSetup error: If you use the NVT\n"
801	>	" ensemble, you must set either tauThermostat or qMass.\n"
802	>	" Neither of these was found in the BASS file.\n");
803	>	painCave.isFatal = 1;
804	>	simError();
805	>	}
806	>	break;
807
808	+	default:
809	+	sprintf( painCave.errMsg,
810	+	"SimSetup Error. Unrecognized ensemble in case statement.\n");
811	+	painCave.isFatal = 1;
812	+	simError();
813	+	}
814
815
816	+	#ifdef IS_MPI
817	+	mpiSim->mpiRefresh();
818	+	#endif
819	+
820		// initialize the Fortran
821	<
821	>
822	>
823		simnfo->refreshSim();
824
825		if( !strcmp( simnfo->mixingRule, "standard") ){
#	Line 671 \| Line 844 \| void SimSetup::makeAtoms( void ){
844		#endif // is_mpi
845		}
846
674	–	void SimSetup::makeAtoms( void ){
847
848	<	int i, j, k, index;
677	<	double ux, uy, uz, uSqr, u;
678	<	AtomStamp* current_atom;
848	>	void SimSetup::makeMolecules( void ){
849
850	+	int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
851	+	molInit info;
852		DirectionalAtom* dAtom;
853	<	int molIndex, molStart, molEnd, nMemb, lMolIndex;
853	>	LinkedAssign* extras;
854	>	LinkedAssign* current_extra;
855	>	AtomStamp* currentAtom;
856	>	BondStamp* currentBond;
857	>	BendStamp* currentBend;
858	>	TorsionStamp* currentTorsion;
859
860	<	lMolIndex = 0;
861	<	molIndex = 0;
862	<	index = 0;
686	<	for( i=0; i<n_components; i++ ){
860	>	bond_pair* theBonds;
861	>	bend_set* theBends;
862	>	torsion_set* theTorsions;
863
864	<	for( j=0; j<components_nmol[i]; j++ ){
864	>
865	>	//init the forceField paramters
866
867	<	#ifdef IS_MPI
691	<	if( mpiSim->getMyMolStart() <= molIndex &&
692	<	molIndex <= mpiSim->getMyMolEnd() ){
693	<	#endif // is_mpi
867	>	the_ff->readParams();
868
869	<	molStart = index;
870	<	nMemb = comp_stamps[i]->getNAtoms();
697	<	for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
698	<
699	<	current_atom = comp_stamps[i]->getAtom( k );
700	<	if( current_atom->haveOrientation() ){
701	<
702	<	dAtom = new DirectionalAtom(index);
703	<	simnfo->n_oriented++;
704	<	the_atoms[index] = dAtom;
705	<
706	<	ux = current_atom->getOrntX();
707	<	uy = current_atom->getOrntY();
708	<	uz = current_atom->getOrntZ();
709	<
710	<	uSqr = (ux * ux) + (uy * uy) + (uz * uz);
711	<
712	<	u = sqrt( uSqr );
713	<	ux = ux / u;
714	<	uy = uy / u;
715	<	uz = uz / u;
716	<
717	<	dAtom->setSUx( ux );
718	<	dAtom->setSUy( uy );
719	<	dAtom->setSUz( uz );
720	<	}
721	<	else{
722	<	the_atoms[index] = new GeneralAtom(index);
723	<	}
724	<	the_atoms[index]->setType( current_atom->getType() );
725	<	the_atoms[index]->setIndex( index );
726	<
727	<	// increment the index and repeat;
728	<	index++;
729	<	}
730	<
731	<	molEnd = index -1;
732	<	the_molecules[lMolIndex].setNMembers( nMemb );
733	<	the_molecules[lMolIndex].setStartAtom( molStart );
734	<	the_molecules[lMolIndex].setEndAtom( molEnd );
735	<	the_molecules[lMolIndex].setStampID( i );
736	<	lMolIndex++;
869	>
870	>	// init the atoms
871
872	<	#ifdef IS_MPI
873	<	}
874	<	#endif //is_mpi
875	<
876	<	molIndex++;
743	<	}
744	<	}
745	<
746	<	#ifdef IS_MPI
747	<	for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] );
872	>	double ux, uy, uz, u, uSqr;
873	>
874	>	atomOffset = 0;
875	>	excludeOffset = 0;
876	>	for(i=0; i<simnfo->n_mol; i++){
877
878	<	delete[] globalIndex;
878	>	stampID = the_molecules[i].getStampID();
879
880	<	mpiSim->mpiRefresh();
881	<	#endif //IS_MPI
882	<
883	<	the_ff->initializeAtoms();
884	<	}
880	>	info.nAtoms = comp_stamps[stampID]->getNAtoms();
881	>	info.nBonds = comp_stamps[stampID]->getNBonds();
882	>	info.nBends = comp_stamps[stampID]->getNBends();
883	>	info.nTorsions = comp_stamps[stampID]->getNTorsions();
884	>	info.nExcludes = info.nBonds + info.nBends + info.nTorsions;
885
886	<	void SimSetup::makeBonds( void ){
887	<
888	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
889	<	bond_pair* the_bonds;
890	<	BondStamp* current_bond;
762	<
763	<	the_bonds = new bond_pair[tot_bonds];
764	<	index = 0;
765	<	offset = 0;
766	<	molIndex = 0;
767	<
768	<	for( i=0; i<n_components; i++ ){
769	<
770	<	for( j=0; j<components_nmol[i]; j++ ){
771	<
772	<	#ifdef IS_MPI
773	<	if( mpiSim->getMyMolStart() <= molIndex &&
774	<	molIndex <= mpiSim->getMyMolEnd() ){
775	<	#endif // is_mpi
776	<
777	<	for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
778	<
779	<	current_bond = comp_stamps[i]->getBond( k );
780	<	the_bonds[index].a = current_bond->getA() + offset;
781	<	the_bonds[index].b = current_bond->getB() + offset;
886	>	info.myAtoms = &the_atoms[atomOffset];
887	>	info.myExcludes = &the_excludes[excludeOffset];
888	>	info.myBonds = new Bond*[info.nBonds];
889	>	info.myBends = new Bend*[info.nBends];
890	>	info.myTorsions = new Torsion*[info.nTorsions];
891
892	<	exI = the_bonds[index].a;
893	<	exJ = the_bonds[index].b;
894	<
895	<	// exclude_I must always be the smaller of the pair
896	<	if( exI > exJ ){
897	<	tempEx = exI;
898	<	exI = exJ;
899	<	exJ = tempEx;
900	<	}
901	<
793	<
794	<	#ifdef IS_MPI
795	<
796	<	the_excludes[index*2] =
797	<	the_atoms[exI]->getGlobalIndex() + 1;
798	<	the_excludes[index*2 + 1] =
799	<	the_atoms[exJ]->getGlobalIndex() + 1;
800	<
801	<	#else // isn't MPI
802	<
803	<	the_excludes[index*2] = exI + 1;
804	<	the_excludes[index*2 + 1] = exJ + 1;
805	<	// fortran index from 1 (hence the +1 in the indexing)
806	<	#endif //is_mpi
807	<
808	<	// increment the index and repeat;
809	<	index++;
810	<	}
811	<	offset += comp_stamps[i]->getNAtoms();
892	>	theBonds = new bond_pair[info.nBonds];
893	>	theBends = new bend_set[info.nBends];
894	>	theTorsions = new torsion_set[info.nTorsions];
895	>
896	>	// make the Atoms
897	>
898	>	for(j=0; j<info.nAtoms; j++){
899	>
900	>	currentAtom = comp_stamps[stampID]->getAtom( j );
901	>	if( currentAtom->haveOrientation() ){
902
903	<	#ifdef IS_MPI
903	>	dAtom = new DirectionalAtom(j + atomOffset);
904	>	simnfo->n_oriented++;
905	>	info.myAtoms[j] = dAtom;
906	>
907	>	ux = currentAtom->getOrntX();
908	>	uy = currentAtom->getOrntY();
909	>	uz = currentAtom->getOrntZ();
910	>
911	>	uSqr = (ux * ux) + (uy * uy) + (uz * uz);
912	>
913	>	u = sqrt( uSqr );
914	>	ux = ux / u;
915	>	uy = uy / u;
916	>	uz = uz / u;
917	>
918	>	dAtom->setSUx( ux );
919	>	dAtom->setSUy( uy );
920	>	dAtom->setSUz( uz );
921		}
922	<	#endif //is_mpi
922	>	else{
923	>	info.myAtoms[j] = new GeneralAtom(j + atomOffset);
924	>	}
925	>	info.myAtoms[j]->setType( currentAtom->getType() );
926	>
927	>	#ifdef IS_MPI
928
929	<	molIndex++;
930	<	}
931	<	}
929	>	info.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
930	>
931	>	#endif // is_mpi
932	>	}
933	>
934	>	// make the bonds
935	>	for(j=0; j<info.nBonds; j++){
936	>
937	>	currentBond = comp_stamps[stampID]->getBond( j );
938	>	theBonds[j].a = currentBond->getA() + atomOffset;
939	>	theBonds[j].b = currentBond->getB() + atomOffset;
940
941	<	the_ff->initializeBonds( the_bonds );
942	<	}
941	>	exI = theBonds[j].a;
942	>	exJ = theBonds[j].b;
943
944	<	void SimSetup::makeBends( void ){
945	<
946	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
947	<	bend_set* the_bends;
948	<	BendStamp* current_bend;
949	<	LinkedAssign* extras;
830	<	LinkedAssign* current_extra;
831	<
832	<
833	<	the_bends = new bend_set[tot_bends];
834	<	index = 0;
835	<	offset = 0;
836	<	molIndex = 0;
837	<	for( i=0; i<n_components; i++ ){
838	<
839	<	for( j=0; j<components_nmol[i]; j++ ){
840	<
944	>	// exclude_I must always be the smaller of the pair
945	>	if( exI > exJ ){
946	>	tempEx = exI;
947	>	exI = exJ;
948	>	exJ = tempEx;
949	>	}
950		#ifdef IS_MPI
951	<	if( mpiSim->getMyMolStart() <= molIndex &&
952	<	molIndex <= mpiSim->getMyMolEnd() ){
953	<	#endif // is_mpi
951	>	tempEx = exI;
952	>	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
953	>	tempEx = exJ;
954	>	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
955	>
956	>	the_excludes[j+excludeOffset]->setPair( exI, exJ );
957	>	#else // isn't MPI
958
959	<	for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
959	>	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
960	>	#endif //is_mpi
961	>	}
962	>	excludeOffset += info.nBonds;
963	>
964	>	//make the bends
965	>	for(j=0; j<info.nBends; j++){
966	>
967	>	currentBend = comp_stamps[stampID]->getBend( j );
968	>	theBends[j].a = currentBend->getA() + atomOffset;
969	>	theBends[j].b = currentBend->getB() + atomOffset;
970	>	theBends[j].c = currentBend->getC() + atomOffset;
971
972	<	current_bend = comp_stamps[i]->getBend( k );
849	<	the_bends[index].a = current_bend->getA() + offset;
850	<	the_bends[index].b = current_bend->getB() + offset;
851	<	the_bends[index].c = current_bend->getC() + offset;
852	<
853	<	if( current_bend->haveExtras() ){
972	>	if( currentBend->haveExtras() ){
973
974	<	extras = current_bend->getExtras();
975	<	current_extra = extras;
974	>	extras = currentBend->getExtras();
975	>	current_extra = extras;
976
977	<	while( current_extra != NULL ){
978	<	if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
977	>	while( current_extra != NULL ){
978	>	if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
979
980	<	switch( current_extra->getType() ){
980	>	switch( current_extra->getType() ){
981	>
982	>	case 0:
983	>	theBends[j].ghost =
984	>	current_extra->getInt() + atomOffset;
985	>	theBends[j].isGhost = 1;
986	>	break;
987
988	<	case 0:
989	<	the_bends[index].ghost =
990	<	current_extra->getInt() + offset;
991	<	the_bends[index].isGhost = 1;
992	<	break;
868	<
869	<	case 1:
870	<	the_bends[index].ghost =
871	<	(int)current_extra->getDouble() + offset;
872	<	the_bends[index].isGhost = 1;
873	<	break;
874	<
875	<	default:
876	<	sprintf( painCave.errMsg,
877	<	"SimSetup Error: ghostVectorSource was neiter a "
878	<	"double nor an int.\n"
879	<	"-->Bend[%d] in %s\n",
880	<	k, comp_stamps[i]->getID() );
881	<	painCave.isFatal = 1;
882	<	simError();
883	<	}
884	<	}
988	>	case 1:
989	>	theBends[j].ghost =
990	>	(int)current_extra->getDouble() + atomOffset;
991	>	theBends[j].isGhost = 1;
992	>	break;
993
994	<	else{
995	<
996	<	sprintf( painCave.errMsg,
997	<	"SimSetup Error: unhandled bend assignment:\n"
998	<	" -->%s in Bend[%d] in %s\n",
999	<	current_extra->getlhs(),
1000	<	k, comp_stamps[i]->getID() );
1001	<	painCave.isFatal = 1;
894	<	simError();
895	<	}
896	<
897	<	current_extra = current_extra->getNext();
994	>	default:
995	>	sprintf( painCave.errMsg,
996	>	"SimSetup Error: ghostVectorSource was neither a "
997	>	"double nor an int.\n"
998	>	"-->Bend[%d] in %s\n",
999	>	j, comp_stamps[stampID]->getID() );
1000	>	painCave.isFatal = 1;
1001	>	simError();
1002		}
1003		}
1004
901	–	if( !the_bends[index].isGhost ){
902	–
903	–	exI = the_bends[index].a;
904	–	exJ = the_bends[index].c;
905	–	}
1005		else{
1006
1007	<	exI = the_bends[index].a;
1008	<	exJ = the_bends[index].b;
1007	>	sprintf( painCave.errMsg,
1008	>	"SimSetup Error: unhandled bend assignment:\n"
1009	>	" -->%s in Bend[%d] in %s\n",
1010	>	current_extra->getlhs(),
1011	>	j, comp_stamps[stampID]->getID() );
1012	>	painCave.isFatal = 1;
1013	>	simError();
1014		}
1015
1016	<	// exclude_I must always be the smaller of the pair
1017	<	if( exI > exJ ){
1018	<	tempEx = exI;
1019	<	exI = exJ;
1020	<	exJ = tempEx;
1021	<	}
1016	>	current_extra = current_extra->getNext();
1017	>	}
1018	>	}
1019	>
1020	>	if( !theBends[j].isGhost ){
1021	>
1022	>	exI = theBends[j].a;
1023	>	exJ = theBends[j].c;
1024	>	}
1025	>	else{
1026	>
1027	>	exI = theBends[j].a;
1028	>	exJ = theBends[j].b;
1029	>	}
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	>	#ifdef IS_MPI
1038	>	tempEx = exI;
1039	>	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
1040	>	tempEx = exJ;
1041	>	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
1042	>
1043	>	the_excludes[j+excludeOffset]->setPair( exI, exJ );
1044	>	#else // isn't MPI
1045	>	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
1046	>	#endif //is_mpi
1047	>	}
1048	>	excludeOffset += info.nBends;
1049
1050	+	for(j=0; j<info.nTorsions; j++){
1051	+
1052	+	currentTorsion = comp_stamps[stampID]->getTorsion( j );
1053	+	theTorsions[j].a = currentTorsion->getA() + atomOffset;
1054	+	theTorsions[j].b = currentTorsion->getB() + atomOffset;
1055	+	theTorsions[j].c = currentTorsion->getC() + atomOffset;
1056	+	theTorsions[j].d = currentTorsion->getD() + atomOffset;
1057	+
1058	+	exI = theTorsions[j].a;
1059	+	exJ = theTorsions[j].d;
1060
1061	+	// exclude_I must always be the smaller of the pair
1062	+	if( exI > exJ ){
1063	+	tempEx = exI;
1064	+	exI = exJ;
1065	+	exJ = tempEx;
1066	+	}
1067		#ifdef IS_MPI
1068	<
1069	<	the_excludes[(index + tot_bonds)*2] =
1070	<	the_atoms[exI]->getGlobalIndex() + 1;
1071	<	the_excludes[(index + tot_bonds)*2 + 1] =
1072	<	the_atoms[exJ]->getGlobalIndex() + 1;
1073	<
1068	>	tempEx = exI;
1069	>	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
1070	>	tempEx = exJ;
1071	>	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
1072	>
1073	>	the_excludes[j+excludeOffset]->setPair( exI, exJ );
1074		#else // isn't MPI
1075	<
929	<	the_excludes[(index + tot_bonds)*2] = exI + 1;
930	<	the_excludes[(index + tot_bonds)*2 + 1] = exJ + 1;
931	<	// fortran index from 1 (hence the +1 in the indexing)
1075	>	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
1076		#endif //is_mpi
933	–
934	–
935	–	// increment the index and repeat;
936	–	index++;
937	–	}
938	–	offset += comp_stamps[i]->getNAtoms();
939	–
940	–	#ifdef IS_MPI
941	–	}
942	–	#endif //is_mpi
943	–
944	–	molIndex++;
1077		}
1078	<	}
1078	>	excludeOffset += info.nTorsions;
1079	>
1080	>
1081	>	// send the arrays off to the forceField for init.
1082	>
1083	>	the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
1084	>	the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
1085	>	the_ff->initializeBends( info.nBends, info.myBends, theBends );
1086	>	the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
1087
948	–	#ifdef IS_MPI
949	–	sprintf( checkPointMsg,
950	–	"Successfully created the bends list.\n" );
951	–	MPIcheckPoint();
952	–	#endif // is_mpi
953	–
1088
1089	<	the_ff->initializeBends( the_bends );
956	<	}
1089	>	the_molecules[i].initialize( info );
1090
958	–	void SimSetup::makeTorsions( void ){
1091
1092	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1093	<	torsion_set* the_torsions;
1094	<	TorsionStamp* current_torsion;
1092	>	atomOffset += info.nAtoms;
1093	>	delete[] theBonds;
1094	>	delete[] theBends;
1095	>	delete[] theTorsions;
1096	>	}
1097
964	–	the_torsions = new torsion_set[tot_torsions];
965	–	index = 0;
966	–	offset = 0;
967	–	molIndex = 0;
968	–	for( i=0; i<n_components; i++ ){
969	–
970	–	for( j=0; j<components_nmol[i]; j++ ){
971	–
1098		#ifdef IS_MPI
1099	<	if( mpiSim->getMyMolStart() <= molIndex &&
1100	<	molIndex <= mpiSim->getMyMolEnd() ){
1101	<	#endif // is_mpi
1099	>	sprintf( checkPointMsg, "all molecules initialized succesfully" );
1100	>	MPIcheckPoint();
1101	>	#endif // is_mpi
1102
1103	<	for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
1103	>	// clean up the forcefield
1104	>	the_ff->calcRcut();
1105	>	the_ff->cleanMe();
1106
979	–	current_torsion = comp_stamps[i]->getTorsion( k );
980	–	the_torsions[index].a = current_torsion->getA() + offset;
981	–	the_torsions[index].b = current_torsion->getB() + offset;
982	–	the_torsions[index].c = current_torsion->getC() + offset;
983	–	the_torsions[index].d = current_torsion->getD() + offset;
984	–
985	–	exI = the_torsions[index].a;
986	–	exJ = the_torsions[index].d;
987	–
988	–
989	–	// exclude_I must always be the smaller of the pair
990	–	if( exI > exJ ){
991	–	tempEx = exI;
992	–	exI = exJ;
993	–	exJ = tempEx;
994	–	}
995	–
996	–
997	–	#ifdef IS_MPI
998	–
999	–	the_excludes[(index + tot_bonds + tot_bends)*2] =
1000	–	the_atoms[exI]->getGlobalIndex() + 1;
1001	–	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] =
1002	–	the_atoms[exJ]->getGlobalIndex() + 1;
1003	–
1004	–	#else // isn't MPI
1005	–
1006	–	the_excludes[(index + tot_bonds + tot_bends)*2] = exI + 1;
1007	–	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] = exJ + 1;
1008	–	// fortran indexes from 1 (hence the +1 in the indexing)
1009	–	#endif //is_mpi
1010	–
1011	–
1012	–	// increment the index and repeat;
1013	–	index++;
1014	–	}
1015	–	offset += comp_stamps[i]->getNAtoms();
1016	–
1017	–	#ifdef IS_MPI
1018	–	}
1019	–	#endif //is_mpi
1020	–
1021	–	molIndex++;
1022	–	}
1023	–	}
1024	–
1025	–	the_ff->initializeTorsions( the_torsions );
1107		}
1108
1109		void SimSetup::initFromBass( void ){
#	Line 1044 \| Line 1125 \| void SimSetup::initFromBass( void ){
1125		have_extra =1;
1126
1127		n_cells = (int)temp3 - 1;
1128	<	cellx = simnfo->box_x / temp3;
1129	<	celly = simnfo->box_y / temp3;
1130	<	cellz = simnfo->box_z / temp3;
1128	>	cellx = simnfo->boxLx / temp3;
1129	>	celly = simnfo->boxLy / temp3;
1130	>	cellz = simnfo->boxLz / temp3;
1131		n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
1132		temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
1133		n_per_extra = (int)ceil( temp1 );
#	Line 1061 \| Line 1142 \| void SimSetup::initFromBass( void ){
1142		}
1143		else{
1144		n_cells = (int)temp3;
1145	<	cellx = simnfo->box_x / temp3;
1146	<	celly = simnfo->box_y / temp3;
1147	<	cellz = simnfo->box_z / temp3;
1145	>	cellx = simnfo->boxLx / temp3;
1146	>	celly = simnfo->boxLy / temp3;
1147	>	cellz = simnfo->boxLz / temp3;
1148		}
1149
1150		current_mol = 0;

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents): Revision 394 by gezelter, Mon Mar 24 21:55:34 2003 UTC vs. Revision 582 by mmeineke, Wed Jul 9 15:33:46 2003 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 394 by gezelter, Mon Mar 24 21:55:34 2003 UTC vs.
Revision 582 by mmeineke, Wed Jul 9 15:33:46 2003 UTC