[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 490 by gezelter, Fri Apr 11 15:16:59 2003 UTC

#	Line 66 \| Line 66 \| void SimSetup::createSim( void ){
66
67		MakeStamps *the_stamps;
68		Globals* the_globals;
69	<	int i, j;
69	>	ExtendedSystem* the_extendedsystem;
70	>	int i, j, k, globalAtomIndex;
71
72		// get the stamps and globals;
73		the_stamps = stamps;
#	Line 80 \| Line 81 \| void SimSetup::createSim( void ){
81		// get the ones we know are there, yet still may need some work.
82		n_components = the_globals->getNComponents();
83		strcpy( force_field, the_globals->getForceField() );
84	+
85	+	// get the ensemble and set up an extended system if we need it:
86		strcpy( ensemble, the_globals->getEnsemble() );
87	+	if( !strcasecmp( ensemble, "NPT" ) ) {
88	+	the_extendedsystem = new ExtendedSystem( simnfo );
89	+	the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
90	+	if (the_globals->haveTargetPressure())
91	+	the_extendedsystem->setTargetPressure(the_globals->getTargetPressure());
92	+	else {
93	+	sprintf( painCave.errMsg,
94	+	"SimSetup error: If you use the constant pressure\n"
95	+	" ensemble, you must set targetPressure.\n"
96	+	" This was found in the BASS file.\n");
97	+	painCave.isFatal = 1;
98	+	simError();
99	+	}
100	+
101	+	if (the_globals->haveTauThermostat())
102	+	the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
103	+	else if (the_globals->haveQmass())
104	+	the_extendedsystem->setQmass(the_globals->getQmass());
105	+	else {
106	+	sprintf( painCave.errMsg,
107	+	"SimSetup error: If you use one of the constant temperature\n"
108	+	" ensembles, you must set either tauThermostat or qMass.\n"
109	+	" Neither of these was found in the BASS file.\n");
110	+	painCave.isFatal = 1;
111	+	simError();
112	+	}
113	+
114	+	if (the_globals->haveTauBarostat())
115	+	the_extendedsystem->setTauBarostat(the_globals->getTauBarostat());
116	+	else {
117	+	sprintf( painCave.errMsg,
118	+	"SimSetup error: If you use the constant pressure\n"
119	+	" ensemble, you must set tauBarostat.\n"
120	+	" This was found in the BASS file.\n");
121	+	painCave.isFatal = 1;
122	+	simError();
123	+	}
124	+
125	+	} else if ( !strcasecmp( ensemble, "NVT") ) {
126	+	the_extendedsystem = new ExtendedSystem( simnfo );
127	+	the_extendedsystem->setTargetTemp(the_globals->getTargetTemp());
128	+
129	+	if (the_globals->haveTauThermostat())
130	+	the_extendedsystem->setTauThermostat(the_globals->getTauThermostat());
131	+	else if (the_globals->haveQmass())
132	+	the_extendedsystem->setQmass(the_globals->getQmass());
133	+	else {
134	+	sprintf( painCave.errMsg,
135	+	"SimSetup error: If you use one of the constant temperature\n"
136	+	" ensembles, you must set either tauThermostat or qMass.\n"
137	+	" Neither of these was found in the BASS file.\n");
138	+	painCave.isFatal = 1;
139	+	simError();
140	+	}
141	+
142	+	} else if ( !strcasecmp( ensemble, "NVE") ) {
143	+	} else {
144	+	sprintf( painCave.errMsg,
145	+	"SimSetup Warning. Unrecognized Ensemble -> %s, "
146	+	"reverting to NVE for this simulation.\n",
147	+	ensemble );
148	+	painCave.isFatal = 0;
149	+	simError();
150	+	strcpy( ensemble, "NVE" );
151	+	}
152		strcpy( simnfo->ensemble, ensemble );
153
154		strcpy( simnfo->mixingRule, the_globals->getMixingRule() );
155		simnfo->usePBC = the_globals->getPBC();
156
157	<
158	<
159	<	if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
160	<	else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
161	<	else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
162	<	else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
157	>	int usesDipoles = 0;
158	>	if( !strcmp( force_field, "TraPPE_Ex" ) ){
159	>	the_ff = new TraPPE_ExFF();
160	>	usesDipoles = 1;
161	>	}
162	>	else if( !strcasecmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
163		else{
164		sprintf( painCave.errMsg,
165		"SimSetup Error. Unrecognized force field -> %s\n",
#	Line 230 \| Line 298 \| void SimSetup::createSim( void ){
298		simnfo->n_torsions = tot_torsions;
299		simnfo->n_SRI = tot_SRI;
300		simnfo->n_mol = tot_nmol;
301	+	simnfo->molMembershipArray = new int[tot_atoms];
302
234	–
303		#ifdef IS_MPI
304
305		// divide the molecules among processors here.
306
307		mpiSim = new mpiSimulation( simnfo );
308
241	–
242	–
309		globalIndex = mpiSim->divideLabor();
310
245	–
246	–
311		// set up the local variables
312
313		int localMol, allMol;
314		int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
315	+
316	+	int* mol2proc = mpiSim->getMolToProcMap();
317	+	int* molCompType = mpiSim->getMolComponentType();
318
319		allMol = 0;
320		localMol = 0;
#	Line 255 \| Line 322 \| void SimSetup::createSim( void ){
322		local_bonds = 0;
323		local_bends = 0;
324		local_torsions = 0;
325	+	globalAtomIndex = 0;
326	+
327	+
328		for( i=0; i<n_components; i++ ){
329
330		for( j=0; j<components_nmol[i]; j++ ){
331
332	<	if( mpiSim->getMyMolStart() <= allMol &&
263	<	allMol <= mpiSim->getMyMolEnd() ){
332	>	if( mol2proc[allMol] == worldRank ){
333
334		local_atoms += comp_stamps[i]->getNAtoms();
335		local_bonds += comp_stamps[i]->getNBonds();
#	Line 268 \| Line 337 \| void SimSetup::createSim( void ){
337		local_torsions += comp_stamps[i]->getNTorsions();
338		localMol++;
339		}
340	<	allMol++;
340	>	for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
341	>	simnfo->molMembershipArray[globalAtomIndex] = allMol;
342	>	globalAtomIndex++;
343	>	}
344	>
345	>	allMol++;
346		}
347		}
348		local_SRI = local_bonds + local_bends + local_torsions;
349
350	+	if (worldRank != 0) {
351	+	for (i =0; i < tot_atoms; i++){
352	+	std::cerr << "i = " << i << " molMembershipArray[i] = " << simnfo->molMembershipArray[i] << "\n";
353	+	}
354	+	}
355
356		simnfo->n_atoms = mpiSim->getMyNlocal();
357
358		if( local_atoms != simnfo->n_atoms ){
359		sprintf( painCave.errMsg,
360		"SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
361	<	" localAtom (%d) are note equal.\n",
361	>	" localAtom (%d) are not equal.\n",
362		simnfo->n_atoms,
363		local_atoms );
364		painCave.isFatal = 1;
#	Line 304 \| Line 383 \| void SimSetup::createSim( void ){
383		Atom::createArrays(simnfo->n_atoms);
384		the_atoms = new Atom*[simnfo->n_atoms];
385		the_molecules = new Molecule[simnfo->n_mol];
386	+	int molIndex;
387
388	+	// initialize the molecule's stampID's
389
390	<	if( simnfo->n_SRI ){
391	<	the_sris = new SRI*[simnfo->n_SRI];
392	<	the_excludes = new int[2 * simnfo->n_SRI];
390	>	#ifdef IS_MPI
391	>
392	>
393	>	molIndex = 0;
394	>	for(i=0; i<mpiSim->getTotNmol(); i++){
395	>
396	>	if(mol2proc[i] == worldRank ){
397	>	the_molecules[molIndex].setStampID( molCompType[i] );
398	>	the_molecules[molIndex].setMyIndex( molIndex );
399	>	the_molecules[molIndex].setGlobalIndex( i );
400	>	molIndex++;
401	>	}
402	>	}
403	>
404	>	#else // is_mpi
405	>
406	>	molIndex = 0;
407	>	globalAtomIndex = 0;
408	>	for(i=0; i<n_components; i++){
409	>	for(j=0; j<components_nmol[i]; j++ ){
410	>	the_molecules[molIndex].setStampID( i );
411	>	the_molecules[molIndex].setMyIndex( molIndex );
412	>	the_molecules[molIndex].setGlobalIndex( molIndex );
413	>	for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
414	>	simnfo->molMembershipArray[globalAtomIndex] = molIndex;
415	>	globalAtomIndex++;
416	>	}
417	>	molIndex++;
418	>	}
419	>	}
420	>
421	>
422	>	#endif // is_mpi
423	>
424	>
425	>	if( simnfo->n_SRI ){
426	>
427	>	Exclude::createArray(simnfo->n_SRI);
428	>	the_excludes = new Exclude*[simnfo->n_SRI];
429	>	for( int ex=0; ex<simnfo->n_SRI; ex++) the_excludes[ex] = new Exclude(ex);
430		simnfo->globalExcludes = new int;
431	<	simnfo->n_exclude = tot_SRI;
431	>	simnfo->n_exclude = simnfo->n_SRI;
432		}
433		else{
434
435	<	the_excludes = new int[2];
436	<	the_excludes[0] = 0;
437	<	the_excludes[1] = 0;
435	>	Exclude::createArray( 1 );
436	>	the_excludes = new Exclude*;
437	>	the_excludes[0] = new Exclude(0);
438	>	the_excludes[0]->setPair( 0,0 );
439		simnfo->globalExcludes = new int;
440		simnfo->globalExcludes[0] = 0;
441	<
323	<	simnfo->n_exclude = 1;
441	>	simnfo->n_exclude = 0;
442		}
443
444		// set the arrays into the SimInfo object
445
446		simnfo->atoms = the_atoms;
447	<	simnfo->sr_interactions = the_sris;
447	>	simnfo->molecules = the_molecules;
448		simnfo->nGlobalExcludes = 0;
449		simnfo->excludes = the_excludes;
450
#	Line 383 \| Line 501 \| void SimSetup::createSim( void ){
501
502		the_ff->setSimInfo( simnfo );
503
504	<	makeAtoms();
504	>	makeMolecules();
505		simnfo->identArray = new int[simnfo->n_atoms];
506		for(i=0; i<simnfo->n_atoms; i++){
507		simnfo->identArray[i] = the_atoms[i]->getIdent();
508		}
509
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	–
510		if (the_globals->getUseRF() ) {
511		simnfo->useReactionField = 1;
512
#	Line 443 \| Line 548 \| void SimSetup::createSim( void ){
548		}
549		simnfo->dielectric = the_globals->getDielectric();
550		} else {
551	<	if (simnfo->n_dipoles) {
551	>	if (usesDipoles) {
552
553		if( !the_globals->haveECR() ){
554		sprintf( painCave.errMsg,
555	<	"SimSetup Warning: using default value of 1/2 the smallest"
555	>	"SimSetup Warning: using default value of 1/2 the smallest "
556		"box length for the electrostaticCutoffRadius.\n"
557		"I hope you have a very fast processor!\n");
558		painCave.isFatal = 0;
#	Line 463 \| Line 568 \| void SimSetup::createSim( void ){
568
569		if( !the_globals->haveEST() ){
570		sprintf( painCave.errMsg,
571	<	"SimSetup Warning: using default value of 5% of the"
571	>	"SimSetup Warning: using default value of 5%% of the "
572		"electrostaticCutoffRadius for the "
573		"electrostaticSkinThickness\n"
574		);
#	Line 638 \| Line 743 \| void SimSetup::createSim( void ){
743
744		// new AllLong( simnfo );
745
746	<	if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
747	<	if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo, the_ff );
748	<	if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo, the_ff );
749	<	if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
746	>
747	>	if( !strcmp( force_field, "TraPPE_Ex" ) ){
748	>	new Symplectic(simnfo, the_ff, the_extendedsystem);
749	>	}
750	>	else if( !strcmp( force_field, "LJ" ) ){
751	>	new Verlet( *simnfo, the_ff, the_extendedsystem );
752	>	}
753	>	else {
754	>	std::cerr << "I'm a bug.\n";
755	>	fprintf( stderr, "Ima bug. stderr %s\n", force_field);
756	>	}
757	>	#ifdef IS_MPI
758	>	mpiSim->mpiRefresh();
759	>	#endif
760
761	+	// initialize the Fortran
762
763
648	–	// initialize the Fortran
649	–
764		simnfo->refreshSim();
765
766		if( !strcmp( simnfo->mixingRule, "standard") ){
#	Line 674 \| Line 788 \| void SimSetup::makeMolecules( void ){
788
789		void SimSetup::makeMolecules( void ){
790
791	<	int i, j, exI, exJ, tempEx, stampID, atomOffset;
791	>	int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
792		molInit info;
793		DirectionalAtom* dAtom;
794	+	LinkedAssign* extras;
795	+	LinkedAssign* current_extra;
796		AtomStamp* currentAtom;
797		BondStamp* currentBond;
798		BendStamp* currentBend;
799		TorsionStamp* currentTorsion;
800	+
801	+	bond_pair* theBonds;
802	+	bend_set* theBends;
803	+	torsion_set* theTorsions;
804	+
805
806		//init the forceField paramters
807
808		the_ff->readParams();
809
810
811	<	// init the molecules
811	>	// init the atoms
812
813	+	double ux, uy, uz, u, uSqr;
814	+
815		atomOffset = 0;
816	+	excludeOffset = 0;
817		for(i=0; i<simnfo->n_mol; i++){
818
819		stampID = the_molecules[i].getStampID();
#	Line 698 \| Line 822 \| void SimSetup::makeMolecules( void ){
822		info.nBonds = comp_stamps[stampID]->getNBonds();
823		info.nBends = comp_stamps[stampID]->getNBends();
824		info.nTorsions = comp_stamps[stampID]->getNTorsions();
825	<
825	>	info.nExcludes = info.nBonds + info.nBends + info.nTorsions;
826	>
827		info.myAtoms = &the_atoms[atomOffset];
828	+	info.myExcludes = &the_excludes[excludeOffset];
829		info.myBonds = new Bond*[info.nBonds];
830		info.myBends = new Bend*[info.nBends];
831	<	info.myTorsions = new Torsions*[info.nTorsions];
831	>	info.myTorsions = new Torsion*[info.nTorsions];
832
833		theBonds = new bond_pair[info.nBonds];
834		theBends = new bend_set[info.nBends];
#	Line 712 \| Line 838 \| void SimSetup::makeMolecules( void ){
838
839		for(j=0; j<info.nAtoms; j++){
840
841	<	currentAtom = theComponents[stampID]->getAtom( j );
841	>	currentAtom = comp_stamps[stampID]->getAtom( j );
842		if( currentAtom->haveOrientation() ){
843
844		dAtom = new DirectionalAtom(j + atomOffset);
#	Line 747 \| Line 873 \| void SimSetup::makeMolecules( void ){
873		}
874
875		// make the bonds
876	<	for(j=0; j<nBonds; j++){
876	>	for(j=0; j<info.nBonds; j++){
877
878		currentBond = comp_stamps[stampID]->getBond( j );
879		theBonds[j].a = currentBond->getA() + atomOffset;
880		theBonds[j].b = currentBond->getB() + atomOffset;
881
882	<	exI = theBonds[i].a;
883	<	exJ = theBonds[i].b;
882	>	exI = theBonds[j].a;
883	>	exJ = theBonds[j].b;
884
885		// exclude_I must always be the smaller of the pair
886		if( exI > exJ ){
#	Line 763 \| Line 889 \| void SimSetup::makeMolecules( void ){
889		exJ = tempEx;
890		}
891		#ifdef IS_MPI
892	+	tempEx = exI;
893	+	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
894	+	tempEx = exJ;
895	+	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
896
897	<	the_excludes[index*2] =
768	<	the_atoms[exI]->getGlobalIndex() + 1;
769	<	the_excludes[index*2 + 1] =
770	<	the_atoms[exJ]->getGlobalIndex() + 1;
771	<
897	>	the_excludes[j+excludeOffset]->setPair( exI, exJ );
898		#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)
899
900	+	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
901		#endif //is_mpi
902	<
903	<	}
902	>	}
903	>	excludeOffset += info.nBonds;
904
905	<
906	<
907	<
908	<
909	<
910	<
911	<
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++ ){
905	>	//make the bends
906	>	for(j=0; j<info.nBends; j++){
907	>
908	>	currentBend = comp_stamps[stampID]->getBend( j );
909	>	theBends[j].a = currentBend->getA() + atomOffset;
910	>	theBends[j].b = currentBend->getB() + atomOffset;
911	>	theBends[j].c = currentBend->getC() + atomOffset;
912
913	<	current_atom = comp_stamps[i]->getAtom( k );
822	<	if( current_atom->haveOrientation() ){
913	>	if( currentBend->haveExtras() ){
914
915	<	dAtom = new DirectionalAtom(index);
916	<	simnfo->n_oriented++;
826	<	the_atoms[index] = dAtom;
915	>	extras = currentBend->getExtras();
916	>	current_extra = extras;
917
918	<	ux = current_atom->getOrntX();
919	<	uy = current_atom->getOrntY();
920	<	uz = current_atom->getOrntZ();
921	<
922	<	uSqr = (ux * ux) + (uy * uy) + (uz * uz);
923	<
924	<	u = sqrt( uSqr );
925	<	ux = ux / u;
926	<	uy = uy / u;
927	<	uz = uz / u;
928	<
929	<	dAtom->setSUx( ux );
930	<	dAtom->setSUy( uy );
931	<	dAtom->setSUz( uz );
932	<	}
933	<	else{
934	<	the_atoms[index] = new GeneralAtom(index);
935	<	}
936	<	the_atoms[index]->setType( current_atom->getType() );
937	<	the_atoms[index]->setIndex( index );
938	<
939	<	// increment the index and repeat;
940	<	index++;
941	<	}
942	<
943	<	molEnd = index -1;
944	<	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
875	<
876	<	the_ff->initializeAtoms();
877	<	}
878	<
879	<	void SimSetup::makeBonds( void ){
880	<
881	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
882	<	bond_pair* the_bonds;
883	<	BondStamp* current_bond;
884	<
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++ ){
918	>	while( current_extra != NULL ){
919	>	if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
920	>
921	>	switch( current_extra->getType() ){
922	>
923	>	case 0:
924	>	theBends[j].ghost =
925	>	current_extra->getInt() + atomOffset;
926	>	theBends[j].isGhost = 1;
927	>	break;
928	>
929	>	case 1:
930	>	theBends[j].ghost =
931	>	(int)current_extra->getDouble() + atomOffset;
932	>	theBends[j].isGhost = 1;
933	>	break;
934	>
935	>	default:
936	>	sprintf( painCave.errMsg,
937	>	"SimSetup Error: ghostVectorSource was neither a "
938	>	"double nor an int.\n"
939	>	"-->Bend[%d] in %s\n",
940	>	j, comp_stamps[stampID]->getID() );
941	>	painCave.isFatal = 1;
942	>	simError();
943	>	}
944	>	}
945
946	<	current_bond = comp_stamps[i]->getBond( k );
947	<	the_bonds[index].a = current_bond->getA() + offset;
948	<	the_bonds[index].b = current_bond->getB() + offset;
949	<
950	<	exI = the_bonds[index].a;
951	<	exJ = the_bonds[index].b;
952	<
953	<	// exclude_I must always be the smaller of the pair
954	<	if( exI > exJ ){
910	<	tempEx = exI;
911	<	exI = exJ;
912	<	exJ = tempEx;
946	>	else{
947	>
948	>	sprintf( painCave.errMsg,
949	>	"SimSetup Error: unhandled bend assignment:\n"
950	>	" -->%s in Bend[%d] in %s\n",
951	>	current_extra->getlhs(),
952	>	j, comp_stamps[stampID]->getID() );
953	>	painCave.isFatal = 1;
954	>	simError();
955		}
914	–
956
957	<	#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++;
957	>	current_extra = current_extra->getNext();
958		}
959	<	offset += comp_stamps[i]->getNAtoms();
959	>	}
960	>
961	>	if( !theBends[j].isGhost ){
962	>
963	>	exI = theBends[j].a;
964	>	exJ = theBends[j].c;
965	>	}
966	>	else{
967
968	<	#ifdef IS_MPI
968	>	exI = theBends[j].a;
969	>	exJ = theBends[j].b;
970		}
937	–	#endif //is_mpi
971
972	<	molIndex++;
973	<	}
974	<	}
972	>	// exclude_I must always be the smaller of the pair
973	>	if( exI > exJ ){
974	>	tempEx = exI;
975	>	exI = exJ;
976	>	exJ = tempEx;
977	>	}
978	>	#ifdef IS_MPI
979	>	tempEx = exI;
980	>	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
981	>	tempEx = exJ;
982	>	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
983	>
984	>	the_excludes[j+excludeOffset]->setPair( exI, exJ );
985	>	#else // isn't MPI
986	>	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
987	>	#endif //is_mpi
988	>	}
989	>	excludeOffset += info.nBends;
990
991	<	the_ff->initializeBonds( the_bonds );
992	<	}
991	>	for(j=0; j<info.nTorsions; j++){
992	>
993	>	currentTorsion = comp_stamps[stampID]->getTorsion( j );
994	>	theTorsions[j].a = currentTorsion->getA() + atomOffset;
995	>	theTorsions[j].b = currentTorsion->getB() + atomOffset;
996	>	theTorsions[j].c = currentTorsion->getC() + atomOffset;
997	>	theTorsions[j].d = currentTorsion->getD() + atomOffset;
998	>
999	>	exI = theTorsions[j].a;
1000	>	exJ = theTorsions[j].d;
1001
1002	<	void SimSetup::makeBends( void ){
1002	>	// exclude_I must always be the smaller of the pair
1003	>	if( exI > exJ ){
1004	>	tempEx = exI;
1005	>	exI = exJ;
1006	>	exJ = tempEx;
1007	>	}
1008	>	#ifdef IS_MPI
1009	>	tempEx = exI;
1010	>	exI = the_atoms[tempEx]->getGlobalIndex() + 1;
1011	>	tempEx = exJ;
1012	>	exJ = the_atoms[tempEx]->getGlobalIndex() + 1;
1013	>
1014	>	the_excludes[j+excludeOffset]->setPair( exI, exJ );
1015	>	#else // isn't MPI
1016	>	the_excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
1017	>	#endif //is_mpi
1018	>	}
1019	>	excludeOffset += info.nTorsions;
1020
1021	<	int i, j, k, index, offset, molIndex, exI, exJ, tempEx;
1022	<	bend_set* the_bends;
950	<	BendStamp* current_bend;
951	<	LinkedAssign* extras;
952	<	LinkedAssign* current_extra;
953	<
1021	>
1022	>	// send the arrays off to the forceField for init.
1023
1024	<	the_bends = new bend_set[tot_bends];
1025	<	index = 0;
1026	<	offset = 0;
1027	<	molIndex = 0;
959	<	for( i=0; i<n_components; i++ ){
1024	>	the_ff->initializeAtoms( info.nAtoms, info.myAtoms );
1025	>	the_ff->initializeBonds( info.nBonds, info.myBonds, theBonds );
1026	>	the_ff->initializeBends( info.nBends, info.myBends, theBends );
1027	>	the_ff->initializeTorsions( info.nTorsions, info.myTorsions, theTorsions );
1028
961	–	for( j=0; j<components_nmol[i]; j++ ){
1029
1030	<	#ifdef IS_MPI
964	<	if( mpiSim->getMyMolStart() <= molIndex &&
965	<	molIndex <= mpiSim->getMyMolEnd() ){
966	<	#endif // is_mpi
967	<
968	<	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	<	}
1030	>	the_molecules[i].initialize( info );
1031
1032
1033	<	#ifdef IS_MPI
1034	<
1035	<	the_excludes[(index + tot_bonds)*2] =
1036	<	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	<	}
1033	>	atomOffset += info.nAtoms;
1034	>	delete[] theBonds;
1035	>	delete[] theBends;
1036	>	delete[] theTorsions;
1037		}
1038
1039		#ifdef IS_MPI
1040	<	sprintf( checkPointMsg,
1072	<	"Successfully created the bends list.\n" );
1040	>	sprintf( checkPointMsg, "all molecules initialized succesfully" );
1041		MPIcheckPoint();
1042		#endif // is_mpi
1075	–
1043
1044	<	the_ff->initializeBends( the_bends );
1045	<	}
1044	>	// clean up the forcefield
1045	>	the_ff->calcRcut();
1046	>	the_ff->cleanMe();
1047
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 );
1048		}
1049
1050		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 490 by gezelter, Fri Apr 11 15:16:59 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 490 by gezelter, Fri Apr 11 15:16:59 2003 UTC