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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 407 by mmeineke, Wed Mar 26 20:22:02 2003 UTC vs.
Revision 561 by mmeineke, Fri Jun 20 20:29:36 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 NPT_ENS 2
20	+
21	+	#define FF_DUFF 0
22	+	#define FF_LJ 1
23	+
24	+
25		SimSetup::SimSetup(){
26		stamps = new MakeStamps();
27		globals = new Globals();
#	Line 66 \| Line 76 \| void SimSetup::createSim( void ){
76
77		MakeStamps *the_stamps;
78		Globals* the_globals;
79	<	int i, j;
79	>	int i, j, k, globalAtomIndex;
80	>
81	>	int ensembleCase;
82	>	int ffCase;
83	>
84	>	ensembleCase = -1;
85	>	ffCase = -1;
86
87		// get the stamps and globals;
88		the_stamps = stamps;
#	Line 80 \| Line 96 \| void SimSetup::createSim( void ){
96		// get the ones we know are there, yet still may need some work.
97		n_components = the_globals->getNComponents();
98		strcpy( force_field, the_globals->getForceField() );
99	+
100	+	if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
101	+	else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
102	+	else{
103	+	sprintf( painCave.errMsg,
104	+	"SimSetup Error. Unrecognized force field -> %s\n",
105	+	force_field );
106	+	painCave.isFatal = 1;
107	+	simError();
108	+	}
109	+
110	+	// get the ensemble:
111		strcpy( ensemble, the_globals->getEnsemble() );
112	+
113	+	if( !strcasecmp( ensemble, "NVE" )) ensembleCase = NVE_ENS;
114	+	else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
115	+	else if( !strcasecmp( ensemble, "NPT" )) ensembleCase = NPT_ENS;
116	+	else{
117	+	sprintf( painCave.errMsg,
118	+	"SimSetup Warning. Unrecognized Ensemble -> %s, "
119	+	"reverting to NVE for this simulation.\n",
120	+	ensemble );
121	+	painCave.isFatal = 0;
122	+	simError();
123	+	strcpy( ensemble, "NVE" );
124	+	ensembleCase = NVE_ENS;
125	+	}
126		strcpy( simnfo->ensemble, ensemble );
127
128	+
129	+	// if( !strcasecmp( ensemble, "NPT" ) ) {
130	+	// the_extendedsystem = new ExtendedSystem( simnfo );
131	+	// the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
132	+	// if (the_globals->haveTargetPressure())
133	+	// the_extendedsystem->setTargetPressure(the_globals->getTargetPressure());
134	+	// else {
135	+	// sprintf( painCave.errMsg,
136	+	// "SimSetup error: If you use the constant pressure\n"
137	+	// " ensemble, you must set targetPressure.\n"
138	+	// " This was found in the BASS file.\n");
139	+	// painCave.isFatal = 1;
140	+	// simError();
141	+	// }
142	+
143	+	// if (the_globals->haveTauThermostat())
144	+	// the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
145	+	// else if (the_globals->haveQmass())
146	+	// the_extendedsystem->setQmass(the_globals->getQmass());
147	+	// else {
148	+	// sprintf( painCave.errMsg,
149	+	// "SimSetup error: If you use one of the constant temperature\n"
150	+	// " ensembles, you must set either tauThermostat or qMass.\n"
151	+	// " Neither of these was found in the BASS file.\n");
152	+	// painCave.isFatal = 1;
153	+	// simError();
154	+	// }
155	+
156	+	// if (the_globals->haveTauBarostat())
157	+	// the_extendedsystem->setTauBarostat(the_globals->getTauBarostat());
158	+	// else {
159	+	// sprintf( painCave.errMsg,
160	+	// "SimSetup error: If you use the constant pressure\n"
161	+	// " ensemble, you must set tauBarostat.\n"
162	+	// " This was found in the BASS file.\n");
163	+	// painCave.isFatal = 1;
164	+	// simError();
165	+	// }
166	+
167	+	// } else if ( !strcasecmp( ensemble, "NVT") ) {
168	+	// the_extendedsystem = new ExtendedSystem( simnfo );
169	+	// the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
170	+
171	+	// if (the_globals->haveTauThermostat())
172	+	// the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
173	+	// else if (the_globals->haveQmass())
174	+	// the_extendedsystem->setQmass(the_globals->getQmass());
175	+	// else {
176	+	// sprintf( painCave.errMsg,
177	+	// "SimSetup error: If you use one of the constant temperature\n"
178	+	// " ensembles, you must set either tauThermostat or qMass.\n"
179	+	// " Neither of these was found in the BASS file.\n");
180	+	// painCave.isFatal = 1;
181	+	// simError();
182	+	// }
183	+
184		strcpy( simnfo->mixingRule, the_globals->getMixingRule() );
185		simnfo->usePBC = the_globals->getPBC();
186
187	+	int usesDipoles = 0;
188	+	switch( ffCase ){
189
190	+	case FF_DUFF:
191	+	the_ff = new DUFF();
192	+	usesDipoles = 1;
193	+	break;
194
195	<	if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
196	<	else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
197	<	else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
198	<	else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
199	<	else{
195	>	case FF_LJ:
196	>	the_ff = new LJFF();
197	>	break;
198	>
199	>	default:
200		sprintf( painCave.errMsg,
201	<	"SimSetup Error. Unrecognized force field -> %s\n",
98	<	force_field );
201	>	"SimSetup Error. Unrecognized force field in case statement.\n");
202		painCave.isFatal = 1;
203		simError();
204		}
#	Line 105 \| Line 208 \| void SimSetup::createSim( void ){
208		MPIcheckPoint();
209		#endif // is_mpi
210
108	–
109	–
211		// get the components and calculate the tot_nMol and indvidual n_mol
212		the_components = the_globals->getComponents();
213		components_nmol = new int[n_components];
#	Line 230 \| Line 331 \| void SimSetup::createSim( void ){
331		simnfo->n_torsions = tot_torsions;
332		simnfo->n_SRI = tot_SRI;
333		simnfo->n_mol = tot_nmol;
233	–
334
335	+	simnfo->molMembershipArray = new int[tot_atoms];
336	+
337		#ifdef IS_MPI
338
339		// divide the molecules among processors here.
340
341		mpiSim = new mpiSimulation( simnfo );
342
241	–
242	–
343		globalIndex = mpiSim->divideLabor();
344
245	–
246	–
345		// set up the local variables
346
347		int localMol, allMol;
348		int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
349	+
350	+	int* mol2proc = mpiSim->getMolToProcMap();
351	+	int* molCompType = mpiSim->getMolComponentType();
352
353		allMol = 0;
354		localMol = 0;
#	Line 255 \| Line 356 \| void SimSetup::createSim( void ){
356		local_bonds = 0;
357		local_bends = 0;
358		local_torsions = 0;
359	+	globalAtomIndex = 0;
360	+
361	+
362		for( i=0; i<n_components; i++ ){
363
364		for( j=0; j<components_nmol[i]; j++ ){
365
366	<	if( mpiSim->getMyMolStart() <= allMol &&
263	<	allMol <= mpiSim->getMyMolEnd() ){
366	>	if( mol2proc[allMol] == worldRank ){
367
368		local_atoms += comp_stamps[i]->getNAtoms();
369		local_bonds += comp_stamps[i]->getNBonds();
#	Line 268 \| Line 371 \| void SimSetup::createSim( void ){
371		local_torsions += comp_stamps[i]->getNTorsions();
372		localMol++;
373		}
374	<	allMol++;
374	>	for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
375	>	simnfo->molMembershipArray[globalAtomIndex] = allMol;
376	>	globalAtomIndex++;
377	>	}
378	>
379	>	allMol++;
380		}
381		}
382		local_SRI = local_bonds + local_bends + local_torsions;
383
276	–
384		simnfo->n_atoms = mpiSim->getMyNlocal();
385
386		if( local_atoms != simnfo->n_atoms ){
387		sprintf( painCave.errMsg,
388		"SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
389	<	" localAtom (%d) are note equal.\n",
389	>	" localAtom (%d) are not equal.\n",
390		simnfo->n_atoms,
391		local_atoms );
392		painCave.isFatal = 1;
#	Line 304 \| Line 411 \| void SimSetup::createSim( void ){
411		Atom::createArrays(simnfo->n_atoms);
412		the_atoms = new Atom*[simnfo->n_atoms];
413		the_molecules = new Molecule[simnfo->n_mol];
414	+	int molIndex;
415
416	+	// initialize the molecule's stampID's
417
418	<	if( simnfo->n_SRI ){
419	<	the_sris = new SRI*[simnfo->n_SRI];
420	<	the_excludes = new int[2 * simnfo->n_SRI];
418	>	#ifdef IS_MPI
419	>
420	>
421	>	molIndex = 0;
422	>	for(i=0; i<mpiSim->getTotNmol(); i++){
423	>
424	>	if(mol2proc[i] == worldRank ){
425	>	the_molecules[molIndex].setStampID( molCompType[i] );
426	>	the_molecules[molIndex].setMyIndex( molIndex );
427	>	the_molecules[molIndex].setGlobalIndex( i );
428	>	molIndex++;
429	>	}
430	>	}
431	>
432	>	#else // is_mpi
433	>
434	>	molIndex = 0;
435	>	globalAtomIndex = 0;
436	>	for(i=0; i<n_components; i++){
437	>	for(j=0; j<components_nmol[i]; j++ ){
438	>	the_molecules[molIndex].setStampID( i );
439	>	the_molecules[molIndex].setMyIndex( molIndex );
440	>	the_molecules[molIndex].setGlobalIndex( molIndex );
441	>	for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
442	>	simnfo->molMembershipArray[globalAtomIndex] = molIndex;
443	>	globalAtomIndex++;
444	>	}
445	>	molIndex++;
446	>	}
447	>	}
448	>
449	>
450	>	#endif // is_mpi
451	>
452	>
453	>	if( simnfo->n_SRI ){
454	>
455	>	Exclude::createArray(simnfo->n_SRI);
456	>	the_excludes = new Exclude*[simnfo->n_SRI];
457	>	for( int ex=0; ex<simnfo->n_SRI; ex++) the_excludes[ex] = new Exclude(ex);
458		simnfo->globalExcludes = new int;
459	<	simnfo->n_exclude = tot_SRI;
459	>	simnfo->n_exclude = simnfo->n_SRI;
460		}
461		else{
462
463	<	the_excludes = new int[2];
464	<	the_excludes[0] = 0;
465	<	the_excludes[1] = 0;
463	>	Exclude::createArray( 1 );
464	>	the_excludes = new Exclude*;
465	>	the_excludes[0] = new Exclude(0);
466	>	the_excludes[0]->setPair( 0,0 );
467		simnfo->globalExcludes = new int;
468		simnfo->globalExcludes[0] = 0;
469	<
323	<	simnfo->n_exclude = 1;
469	>	simnfo->n_exclude = 0;
470		}
471
472		// set the arrays into the SimInfo object
473
474		simnfo->atoms = the_atoms;
475	<	simnfo->sr_interactions = the_sris;
475	>	simnfo->molecules = the_molecules;
476		simnfo->nGlobalExcludes = 0;
477		simnfo->excludes = the_excludes;
478
#	Line 383 \| Line 529 \| void SimSetup::createSim( void ){
529
530		the_ff->setSimInfo( simnfo );
531
532	<	makeAtoms();
532	>	makeMolecules();
533		simnfo->identArray = new int[simnfo->n_atoms];
534		for(i=0; i<simnfo->n_atoms; i++){
535		simnfo->identArray[i] = the_atoms[i]->getIdent();
536		}
537
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	–
538		if (the_globals->getUseRF() ) {
539		simnfo->useReactionField = 1;
540
#	Line 443 \| Line 576 \| void SimSetup::createSim( void ){
576		}
577		simnfo->dielectric = the_globals->getDielectric();
578		} else {
579	<	if (simnfo->n_dipoles) {
579	>	if (usesDipoles) {
580
581		if( !the_globals->haveECR() ){
582		sprintf( painCave.errMsg,
583	<	"SimSetup Warning: using default value of 1/2 the smallest"
583	>	"SimSetup Warning: using default value of 1/2 the smallest "
584		"box length for the electrostaticCutoffRadius.\n"
585		"I hope you have a very fast processor!\n");
586		painCave.isFatal = 0;
#	Line 463 \| Line 596 \| void SimSetup::createSim( void ){
596
597		if( !the_globals->haveEST() ){
598		sprintf( painCave.errMsg,
599	<	"SimSetup Warning: using default value of 5% of the"
599	>	"SimSetup Warning: using default value of 5%% of the "
600		"electrostaticCutoffRadius for the "
601		"electrostaticSkinThickness\n"
602		);
#	Line 634 \| Line 767 \| void SimSetup::createSim( void ){
767		if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
768
769
770	<	// // make the longe range forces and the integrator
770	>	// make the integrator
771	>
772	>
773	>	NVT* myNVT = NULL;
774	>	switch( ensembleCase ){
775
776	<	// new AllLong( simnfo );
776	>	case NVE_ENS:
777	>	new NVE( simnfo, the_ff );
778	>	break;
779
780	<	if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
781	<	if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo, the_ff );
782	<	if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo, the_ff );
783	<	if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
780	>	case NVT_ENS:
781	>	myNVT = new NVT( simnfo, the_ff );
782	>	myNVT->setTargetTemp(the_globals->getTargetTemp());
783	>
784	>	if (the_globals->haveTauThermostat())
785	>	myNVT->setTauThermostat(the_globals->getTauThermostat());
786	>	else if (the_globals->haveQmass())
787	>	myNVT->setQmass(the_globals->getQmass());
788	>	else {
789	>	sprintf( painCave.errMsg,
790	>	"SimSetup error: If you use the NVT\n"
791	>	" ensemble, you must set either tauThermostat or qMass.\n"
792	>	" Neither of these was found in the BASS file.\n");
793	>	painCave.isFatal = 1;
794	>	simError();
795	>	}
796	>	break;
797	>
798	>	default:
799	>	sprintf( painCave.errMsg,
800	>	"SimSetup Error. Unrecognized ensemble in case statement.\n");
801	>	painCave.isFatal = 1;
802	>	simError();
803	>	}
804
805
806	+	#ifdef IS_MPI
807	+	mpiSim->mpiRefresh();
808	+	#endif
809
810		// initialize the Fortran
811	<
811	>
812	>
813		simnfo->refreshSim();
814
815		if( !strcmp( simnfo->mixingRule, "standard") ){
#	Line 674 \| Line 837 \| void SimSetup::makeMolecules( void ){
837
838		void SimSetup::makeMolecules( void ){
839
840	<	int i, j, exI, exJ, tempEx, stampID, atomOffset;
840	>	int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
841		molInit info;
842		DirectionalAtom* dAtom;
843	+	LinkedAssign* extras;
844	+	LinkedAssign* current_extra;
845		AtomStamp* currentAtom;
846		BondStamp* currentBond;
847		BendStamp* currentBend;
848		TorsionStamp* currentTorsion;
849	+
850	+	bond_pair* theBonds;
851	+	bend_set* theBends;
852	+	torsion_set* theTorsions;
853	+
854
855		//init the forceField paramters
856
857		the_ff->readParams();
858
859
860	<	// init the molecules
860	>	// init the atoms
861
862	+	double ux, uy, uz, u, uSqr;
863	+
864		atomOffset = 0;
865	+	excludeOffset = 0;
866		for(i=0; i<simnfo->n_mol; i++){
867
868		stampID = the_molecules[i].getStampID();
#	Line 698 \| Line 871 \| void SimSetup::makeMolecules( void ){
871		info.nBonds = comp_stamps[stampID]->getNBonds();
872		info.nBends = comp_stamps[stampID]->getNBends();
873		info.nTorsions = comp_stamps[stampID]->getNTorsions();
874	<
874	>	info.nExcludes = info.nBonds + info.nBends + info.nTorsions;
875	>
876		info.myAtoms = &the_atoms[atomOffset];
877	+	info.myExcludes = &the_excludes[excludeOffset];
878		info.myBonds = new Bond*[info.nBonds];
879		info.myBends = new Bend*[info.nBends];
880	<	info.myTorsions = new Torsions*[info.nTorsions];
880	>	info.myTorsions = new Torsion*[info.nTorsions];
881
882		theBonds = new bond_pair[info.nBonds];
883		theBends = new bend_set[info.nBends];
#	Line 712 \| Line 887 \| void SimSetup::makeMolecules( void ){
887
888		for(j=0; j<info.nAtoms; j++){
889
890	<	currentAtom = theComponents[stampID]->getAtom( j );
890	>	currentAtom = comp_stamps[stampID]->getAtom( j );
891		if( currentAtom->haveOrientation() ){
892
893		dAtom = new DirectionalAtom(j + atomOffset);
#	Line 747 \| Line 922 \| void SimSetup::makeMolecules( void ){
922		}
923
924		// make the bonds
925	<	for(j=0; j<nBonds; j++){
925	>	for(j=0; j<info.nBonds; j++){
926
927		currentBond = comp_stamps[stampID]->getBond( j );
928		theBonds[j].a = currentBond->getA() + atomOffset;
929		theBonds[j].b = currentBond->getB() + atomOffset;
930
931	<	exI = theBonds[i].a;
932	<	exJ = theBonds[i].b;
931	>	exI = theBonds[j].a;
932	>	exJ = theBonds[j].b;
933
934		// exclude_I must always be the smaller of the pair
935		if( exI > exJ ){
#	Line 763 \| Line 938 \| void SimSetup::makeMolecules( void ){
938		exJ = tempEx;
939		}
940		#ifdef IS_MPI
941	<
942	<	the_excludes[index*2] =
943	<	the_atoms[exI]->getGlobalIndex() + 1;
944	<	the_excludes[index*2 + 1] =
770	<	the_atoms[exJ]->getGlobalIndex() + 1;
941	>	tempEx = exI;
942	>	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
943	>	tempEx = exJ;
944	>	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
945
946	+	the_excludes[j+excludeOffset]->setPair( exI, exJ );
947		#else // isn't MPI
773	–
774	–	the_excludes[index*2] = exI + 1;
775	–	the_excludes[index*2 + 1] = exJ + 1;
776	–	// fortran index from 1 (hence the +1 in the indexing)
948
949	+	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
950		#endif //is_mpi
951	<
952	<	}
951	>	}
952	>	excludeOffset += info.nBonds;
953
954	<
955	<
956	<
957	<
958	<
959	<
960	<
789	<
790	<
791	<
792	<
793	<
794	<
795	<
796	<	void SimSetup::makeAtoms( void ){
797	<
798	<	int i, j, k, index;
799	<	double ux, uy, uz, uSqr, u;
800	<	AtomStamp* current_atom;
801	<
802	<	DirectionalAtom* dAtom;
803	<	int molIndex, molStart, molEnd, nMemb, lMolIndex;
804	<
805	<	lMolIndex = 0;
806	<	molIndex = 0;
807	<	index = 0;
808	<	for( i=0; i<n_components; i++ ){
809	<
810	<	for( j=0; j<components_nmol[i]; j++ ){
811	<
812	<	#ifdef IS_MPI
813	<	if( mpiSim->getMyMolStart() <= molIndex &&
814	<	molIndex <= mpiSim->getMyMolEnd() ){
815	<	#endif // is_mpi
816	<
817	<	molStart = index;
818	<	nMemb = comp_stamps[i]->getNAtoms();
819	<	for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
954	>	//make the bends
955	>	for(j=0; j<info.nBends; j++){
956	>
957	>	currentBend = comp_stamps[stampID]->getBend( j );
958	>	theBends[j].a = currentBend->getA() + atomOffset;
959	>	theBends[j].b = currentBend->getB() + atomOffset;
960	>	theBends[j].c = currentBend->getC() + atomOffset;
961
962	<	current_atom = comp_stamps[i]->getAtom( k );
822	<	if( current_atom->haveOrientation() ){
962	>	if( currentBend->haveExtras() ){
963
964	<	dAtom = new DirectionalAtom(index);
965	<	simnfo->n_oriented++;
826	<	the_atoms[index] = dAtom;
964	>	extras = currentBend->getExtras();
965	>	current_extra = extras;
966
967	<	ux = current_atom->getOrntX();
968	<	uy = current_atom->getOrntY();
969	<	uz = current_atom->getOrntZ();
970	<
971	<	uSqr = (ux * ux) + (uy * uy) + (uz * uz);
972	<
973	<	u = sqrt( uSqr );
974	<	ux = ux / u;
975	<	uy = uy / u;
976	<	uz = uz / u;
977	<
978	<	dAtom->setSUx( ux );
979	<	dAtom->setSUy( uy );
980	<	dAtom->setSUz( uz );
981	<	}
982	<	else{
983	<	the_atoms[index] = new GeneralAtom(index);
984	<	}
985	<	the_atoms[index]->setType( current_atom->getType() );
986	<	the_atoms[index]->setIndex( index );
987	<
988	<	// increment the index and repeat;
989	<	index++;
990	<	}
991	<
992	<	molEnd = index -1;
993	<	the_molecules[lMolIndex].setNMembers( nMemb );
855	<	the_molecules[lMolIndex].setStartAtom( molStart );
856	<	the_molecules[lMolIndex].setEndAtom( molEnd );
857	<	the_molecules[lMolIndex].setStampID( i );
858	<	lMolIndex++;
859	<
860	<	#ifdef IS_MPI
861	<	}
862	<	#endif //is_mpi
863	<
864	<	molIndex++;
865	<	}
866	<	}
867	<
868	<	#ifdef IS_MPI
869	<	for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] );
870	<
871	<	delete[] globalIndex;
872	<
873	<	mpiSim->mpiRefresh();
874	<	#endif //IS_MPI
967	>	while( current_extra != NULL ){
968	>	if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
969	>
970	>	switch( current_extra->getType() ){
971	>
972	>	case 0:
973	>	theBends[j].ghost =
974	>	current_extra->getInt() + atomOffset;
975	>	theBends[j].isGhost = 1;
976	>	break;
977	>
978	>	case 1:
979	>	theBends[j].ghost =
980	>	(int)current_extra->getDouble() + atomOffset;
981	>	theBends[j].isGhost = 1;
982	>	break;
983	>
984	>	default:
985	>	sprintf( painCave.errMsg,
986	>	"SimSetup Error: ghostVectorSource was neither a "
987	>	"double nor an int.\n"
988	>	"-->Bend[%d] in %s\n",
989	>	j, comp_stamps[stampID]->getID() );
990	>	painCave.isFatal = 1;
991	>	simError();
992	>	}
993	>	}
994
995	<	the_ff->initializeAtoms();
996	<	}
997	<
998	<	void SimSetup::makeBonds( void ){
999	<
1000	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1001	<	bond_pair* the_bonds;
1002	<	BondStamp* current_bond;
1003	<
885	<	the_bonds = new bond_pair[tot_bonds];
886	<	index = 0;
887	<	offset = 0;
888	<	molIndex = 0;
889	<
890	<	for( i=0; i<n_components; i++ ){
891	<
892	<	for( j=0; j<components_nmol[i]; j++ ){
893	<
894	<	#ifdef IS_MPI
895	<	if( mpiSim->getMyMolStart() <= molIndex &&
896	<	molIndex <= mpiSim->getMyMolEnd() ){
897	<	#endif // is_mpi
898	<
899	<	for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
900	<
901	<	current_bond = comp_stamps[i]->getBond( k );
902	<	the_bonds[index].a = current_bond->getA() + offset;
903	<	the_bonds[index].b = current_bond->getB() + offset;
904	<
905	<	exI = the_bonds[index].a;
906	<	exJ = the_bonds[index].b;
907	<
908	<	// exclude_I must always be the smaller of the pair
909	<	if( exI > exJ ){
910	<	tempEx = exI;
911	<	exI = exJ;
912	<	exJ = tempEx;
995	>	else{
996	>
997	>	sprintf( painCave.errMsg,
998	>	"SimSetup Error: unhandled bend assignment:\n"
999	>	" -->%s in Bend[%d] in %s\n",
1000	>	current_extra->getlhs(),
1001	>	j, comp_stamps[stampID]->getID() );
1002	>	painCave.isFatal = 1;
1003	>	simError();
1004		}
914	–
1005
1006	<	#ifdef IS_MPI
917	<
918	<	the_excludes[index*2] =
919	<	the_atoms[exI]->getGlobalIndex() + 1;
920	<	the_excludes[index*2 + 1] =
921	<	the_atoms[exJ]->getGlobalIndex() + 1;
922	<
923	<	#else // isn't MPI
924	<
925	<	the_excludes[index*2] = exI + 1;
926	<	the_excludes[index*2 + 1] = exJ + 1;
927	<	// fortran index from 1 (hence the +1 in the indexing)
928	<	#endif //is_mpi
929	<
930	<	// increment the index and repeat;
931	<	index++;
1006	>	current_extra = current_extra->getNext();
1007		}
1008	<	offset += comp_stamps[i]->getNAtoms();
1008	>	}
1009	>
1010	>	if( !theBends[j].isGhost ){
1011	>
1012	>	exI = theBends[j].a;
1013	>	exJ = theBends[j].c;
1014	>	}
1015	>	else{
1016
1017	<	#ifdef IS_MPI
1017	>	exI = theBends[j].a;
1018	>	exJ = theBends[j].b;
1019		}
937	–	#endif //is_mpi
1020
1021	<	molIndex++;
1022	<	}
1023	<	}
1021	>	// exclude_I must always be the smaller of the pair
1022	>	if( exI > exJ ){
1023	>	tempEx = exI;
1024	>	exI = exJ;
1025	>	exJ = tempEx;
1026	>	}
1027	>	#ifdef IS_MPI
1028	>	tempEx = exI;
1029	>	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
1030	>	tempEx = exJ;
1031	>	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
1032	>
1033	>	the_excludes[j+excludeOffset]->setPair( exI, exJ );
1034	>	#else // isn't MPI
1035	>	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
1036	>	#endif //is_mpi
1037	>	}
1038	>	excludeOffset += info.nBends;
1039
1040	<	the_ff->initializeBonds( the_bonds );
1041	<	}
1040	>	for(j=0; j<info.nTorsions; j++){
1041	>
1042	>	currentTorsion = comp_stamps[stampID]->getTorsion( j );
1043	>	theTorsions[j].a = currentTorsion->getA() + atomOffset;
1044	>	theTorsions[j].b = currentTorsion->getB() + atomOffset;
1045	>	theTorsions[j].c = currentTorsion->getC() + atomOffset;
1046	>	theTorsions[j].d = currentTorsion->getD() + atomOffset;
1047	>
1048	>	exI = theTorsions[j].a;
1049	>	exJ = theTorsions[j].d;
1050
1051	<	void SimSetup::makeBends( void ){
1051	>	// exclude_I must always be the smaller of the pair
1052	>	if( exI > exJ ){
1053	>	tempEx = exI;
1054	>	exI = exJ;
1055	>	exJ = tempEx;
1056	>	}
1057	>	#ifdef IS_MPI
1058	>	tempEx = exI;
1059	>	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
1060	>	tempEx = exJ;
1061	>	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
1062	>
1063	>	the_excludes[j+excludeOffset]->setPair( exI, exJ );
1064	>	#else // isn't MPI
1065	>	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
1066	>	#endif //is_mpi
1067	>	}
1068	>	excludeOffset += info.nTorsions;
1069
1070	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1071	<	bend_set* the_bends;
950	<	BendStamp* current_bend;
951	<	LinkedAssign* extras;
952	<	LinkedAssign* current_extra;
953	<
1070	>
1071	>	// send the arrays off to the forceField for init.
1072
1073	<	the_bends = new bend_set[tot_bends];
1074	<	index = 0;
1075	<	offset = 0;
1076	<	molIndex = 0;
959	<	for( i=0; i<n_components; i++ ){
960	<
961	<	for( j=0; j<components_nmol[i]; j++ ){
1073	>	the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
1074	>	the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
1075	>	the_ff->initializeBends( info.nBends, info.myBends, theBends );
1076	>	the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
1077
963	–	#ifdef IS_MPI
964	–	if( mpiSim->getMyMolStart() <= molIndex &&
965	–	molIndex <= mpiSim->getMyMolEnd() ){
966	–	#endif // is_mpi
1078
1079	<	for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
969	<
970	<	current_bend = comp_stamps[i]->getBend( k );
971	<	the_bends[index].a = current_bend->getA() + offset;
972	<	the_bends[index].b = current_bend->getB() + offset;
973	<	the_bends[index].c = current_bend->getC() + offset;
974	<
975	<	if( current_bend->haveExtras() ){
976	<
977	<	extras = current_bend->getExtras();
978	<	current_extra = extras;
979	<
980	<	while( current_extra != NULL ){
981	<	if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
982	<
983	<	switch( current_extra->getType() ){
984	<
985	<	case 0:
986	<	the_bends[index].ghost =
987	<	current_extra->getInt() + offset;
988	<	the_bends[index].isGhost = 1;
989	<	break;
990	<
991	<	case 1:
992	<	the_bends[index].ghost =
993	<	(int)current_extra->getDouble() + offset;
994	<	the_bends[index].isGhost = 1;
995	<	break;
996	<
997	<	default:
998	<	sprintf( painCave.errMsg,
999	<	"SimSetup Error: ghostVectorSource was neiter a "
1000	<	"double nor an int.\n"
1001	<	"-->Bend[%d] in %s\n",
1002	<	k, comp_stamps[i]->getID() );
1003	<	painCave.isFatal = 1;
1004	<	simError();
1005	<	}
1006	<	}
1007	<
1008	<	else{
1009	<
1010	<	sprintf( painCave.errMsg,
1011	<	"SimSetup Error: unhandled bend assignment:\n"
1012	<	" -->%s in Bend[%d] in %s\n",
1013	<	current_extra->getlhs(),
1014	<	k, comp_stamps[i]->getID() );
1015	<	painCave.isFatal = 1;
1016	<	simError();
1017	<	}
1018	<
1019	<	current_extra = current_extra->getNext();
1020	<	}
1021	<	}
1022	<
1023	<	if( !the_bends[index].isGhost ){
1024	<
1025	<	exI = the_bends[index].a;
1026	<	exJ = the_bends[index].c;
1027	<	}
1028	<	else{
1029	<
1030	<	exI = the_bends[index].a;
1031	<	exJ = the_bends[index].b;
1032	<	}
1033	<
1034	<	// exclude_I must always be the smaller of the pair
1035	<	if( exI > exJ ){
1036	<	tempEx = exI;
1037	<	exI = exJ;
1038	<	exJ = tempEx;
1039	<	}
1079	>	the_molecules[i].initialize( info );
1080
1081
1082	<	#ifdef IS_MPI
1083	<
1084	<	the_excludes[(index + tot_bonds)*2] =
1085	<	the_atoms[exI]->getGlobalIndex() + 1;
1046	<	the_excludes[(index + tot_bonds)*2 + 1] =
1047	<	the_atoms[exJ]->getGlobalIndex() + 1;
1048	<
1049	<	#else // isn't MPI
1050	<
1051	<	the_excludes[(index + tot_bonds)*2] = exI + 1;
1052	<	the_excludes[(index + tot_bonds)*2 + 1] = exJ + 1;
1053	<	// fortran index from 1 (hence the +1 in the indexing)
1054	<	#endif //is_mpi
1055	<
1056	<
1057	<	// increment the index and repeat;
1058	<	index++;
1059	<	}
1060	<	offset += comp_stamps[i]->getNAtoms();
1061	<
1062	<	#ifdef IS_MPI
1063	<	}
1064	<	#endif //is_mpi
1065	<
1066	<	molIndex++;
1067	<	}
1082	>	atomOffset += info.nAtoms;
1083	>	delete[] theBonds;
1084	>	delete[] theBends;
1085	>	delete[] theTorsions;
1086		}
1087
1088		#ifdef IS_MPI
1089	<	sprintf( checkPointMsg,
1072	<	"Successfully created the bends list.\n" );
1089	>	sprintf( checkPointMsg, "all molecules initialized succesfully" );
1090		MPIcheckPoint();
1091		#endif // is_mpi
1075	–
1092
1093	<	the_ff->initializeBends( the_bends );
1094	<	}
1093	>	// clean up the forcefield
1094	>	the_ff->calcRcut();
1095	>	the_ff->cleanMe();
1096
1080	–	void SimSetup::makeTorsions( void ){
1081	–
1082	–	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1083	–	torsion_set* the_torsions;
1084	–	TorsionStamp* current_torsion;
1085	–
1086	–	the_torsions = new torsion_set[tot_torsions];
1087	–	index = 0;
1088	–	offset = 0;
1089	–	molIndex = 0;
1090	–	for( i=0; i<n_components; i++ ){
1091	–
1092	–	for( j=0; j<components_nmol[i]; j++ ){
1093	–
1094	–	#ifdef IS_MPI
1095	–	if( mpiSim->getMyMolStart() <= molIndex &&
1096	–	molIndex <= mpiSim->getMyMolEnd() ){
1097	–	#endif // is_mpi
1098	–
1099	–	for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
1100	–
1101	–	current_torsion = comp_stamps[i]->getTorsion( k );
1102	–	the_torsions[index].a = current_torsion->getA() + offset;
1103	–	the_torsions[index].b = current_torsion->getB() + offset;
1104	–	the_torsions[index].c = current_torsion->getC() + offset;
1105	–	the_torsions[index].d = current_torsion->getD() + offset;
1106	–
1107	–	exI = the_torsions[index].a;
1108	–	exJ = the_torsions[index].d;
1109	–
1110	–
1111	–	// exclude_I must always be the smaller of the pair
1112	–	if( exI > exJ ){
1113	–	tempEx = exI;
1114	–	exI = exJ;
1115	–	exJ = tempEx;
1116	–	}
1117	–
1118	–
1119	–	#ifdef IS_MPI
1120	–
1121	–	the_excludes[(index + tot_bonds + tot_bends)*2] =
1122	–	the_atoms[exI]->getGlobalIndex() + 1;
1123	–	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] =
1124	–	the_atoms[exJ]->getGlobalIndex() + 1;
1125	–
1126	–	#else // isn't MPI
1127	–
1128	–	the_excludes[(index + tot_bonds + tot_bends)*2] = exI + 1;
1129	–	the_excludes[(index + tot_bonds + tot_bends)*2 + 1] = exJ + 1;
1130	–	// fortran indexes from 1 (hence the +1 in the indexing)
1131	–	#endif //is_mpi
1132	–
1133	–
1134	–	// increment the index and repeat;
1135	–	index++;
1136	–	}
1137	–	offset += comp_stamps[i]->getNAtoms();
1138	–
1139	–	#ifdef IS_MPI
1140	–	}
1141	–	#endif //is_mpi
1142	–
1143	–	molIndex++;
1144	–	}
1145	–	}
1146	–
1147	–	the_ff->initializeTorsions( the_torsions );
1097		}
1098
1099		void SimSetup::initFromBass( void ){

Diff Legend

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

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents): Revision 407 by mmeineke, Wed Mar 26 20:22:02 2003 UTC vs. Revision 561 by mmeineke, Fri Jun 20 20:29:36 2003 UTC

Diff Legend

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 407 by mmeineke, Wed Mar 26 20:22:02 2003 UTC vs.
Revision 561 by mmeineke, Fri Jun 20 20:29:36 2003 UTC