[ViewVC] Diff of: OpenMD/trunk/src/parallel/ForceMatrixDecomposition.cpp

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1571 by gezelter, Fri May 27 16:45:44 2011 UTC vs.
Revision 1584 by gezelter, Fri Jun 17 20:16:35 2011 UTC

#	Line 57 \| Line 57 \| namespace OpenMD {
57		storageLayout_ = sman_->getStorageLayout();
58		ff_ = info_->getForceField();
59		nLocal_ = snap_->getNumberOfAtoms();
60	–	nGroups_ = snap_->getNumberOfCutoffGroups();
60
61	+	nGroups_ = info_->getNLocalCutoffGroups();
62	+	cerr << "in dId, nGroups = " << nGroups_ << "\n";
63		// gather the information for atomtype IDs (atids):
64	<	identsLocal = info_->getIdentArray();
64	>	idents = info_->getIdentArray();
65		AtomLocalToGlobal = info_->getGlobalAtomIndices();
66		cgLocalToGlobal = info_->getGlobalGroupIndices();
67		vector<int> globalGroupMembership = info_->getGlobalGroupMembership();
68	<	vector<RealType> massFactorsLocal = info_->getMassFactors();
68	>	massFactors = info_->getMassFactors();
69	>
70		PairList excludes = info_->getExcludedInteractions();
71		PairList oneTwo = info_->getOneTwoInteractions();
72		PairList oneThree = info_->getOneThreeInteractions();
73		PairList oneFour = info_->getOneFourInteractions();
72	–	vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0);
74
75		#ifdef IS_MPI
76
#	Line 77 \| Line 78 \| namespace OpenMD {
78		AtomCommRealRow = new Communicator<Row,RealType>(nLocal_);
79		AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal_);
80		AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal_);
81	+	AtomCommPotRow = new Communicator<Row,potVec>(nLocal_);
82
83		AtomCommIntColumn = new Communicator<Column,int>(nLocal_);
84		AtomCommRealColumn = new Communicator<Column,RealType>(nLocal_);
85		AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal_);
86		AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal_);
87	+	AtomCommPotColumn = new Communicator<Column,potVec>(nLocal_);
88
89		cgCommIntRow = new Communicator<Row,int>(nGroups_);
90		cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups_);
#	Line 102 \| Line 105 \| namespace OpenMD {
105		cgRowData.setStorageLayout(DataStorage::dslPosition);
106		cgColData.resize(nGroupsInCol_);
107		cgColData.setStorageLayout(DataStorage::dslPosition);
108	+
109	+	identsRow.resize(nAtomsInRow_);
110	+	identsCol.resize(nAtomsInCol_);
111
112	<	vector<vector<RealType> > pot_row(N_INTERACTION_FAMILIES,
113	<	vector<RealType> (nAtomsInRow_, 0.0));
108	<	vector<vector<RealType> > pot_col(N_INTERACTION_FAMILIES,
109	<	vector<RealType> (nAtomsInCol_, 0.0));
112	>	AtomCommIntRow->gather(idents, identsRow);
113	>	AtomCommIntColumn->gather(idents, identsCol);
114
111	–	identsRow.reserve(nAtomsInRow_);
112	–	identsCol.reserve(nAtomsInCol_);
113	–
114	–	AtomCommIntRow->gather(identsLocal, identsRow);
115	–	AtomCommIntColumn->gather(identsLocal, identsCol);
116	–
115		AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal);
116		AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal);
117
118		cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal);
119		cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal);
120
121	<	AtomCommRealRow->gather(massFactorsLocal, massFactorsRow);
122	<	AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol);
121	>	AtomCommRealRow->gather(massFactors, massFactorsRow);
122	>	AtomCommRealColumn->gather(massFactors, massFactorsCol);
123
124		groupListRow_.clear();
125	<	groupListRow_.reserve(nGroupsInRow_);
125	>	groupListRow_.resize(nGroupsInRow_);
126		for (int i = 0; i < nGroupsInRow_; i++) {
127		int gid = cgRowToGlobal[i];
128		for (int j = 0; j < nAtomsInRow_; j++) {
#	Line 135 \| Line 133 \| namespace OpenMD {
133		}
134
135		groupListCol_.clear();
136	<	groupListCol_.reserve(nGroupsInCol_);
136	>	groupListCol_.resize(nGroupsInCol_);
137		for (int i = 0; i < nGroupsInCol_; i++) {
138		int gid = cgColToGlobal[i];
139		for (int j = 0; j < nAtomsInCol_; j++) {
#	Line 145 \| Line 143 \| namespace OpenMD {
143		}
144		}
145
146	<	skipsForRowAtom.clear();
147	<	skipsForRowAtom.reserve(nAtomsInRow_);
146	>	skipsForAtom.clear();
147	>	skipsForAtom.resize(nAtomsInRow_);
148	>	toposForAtom.clear();
149	>	toposForAtom.resize(nAtomsInRow_);
150	>	topoDist.clear();
151	>	topoDist.resize(nAtomsInRow_);
152		for (int i = 0; i < nAtomsInRow_; i++) {
153		int iglob = AtomRowToGlobal[i];
152	–	for (int j = 0; j < nAtomsInCol_; j++) {
153	–	int jglob = AtomColToGlobal[j];
154	–	if (excludes.hasPair(iglob, jglob))
155	–	skipsForRowAtom[i].push_back(j);
156	–	}
157	–	}
154
159	–	toposForRowAtom.clear();
160	–	toposForRowAtom.reserve(nAtomsInRow_);
161	–	for (int i = 0; i < nAtomsInRow_; i++) {
162	–	int iglob = AtomRowToGlobal[i];
163	–	int nTopos = 0;
155		for (int j = 0; j < nAtomsInCol_; j++) {
156	<	int jglob = AtomColToGlobal[j];
156	>	int jglob = AtomColToGlobal[j];
157	>
158	>	if (excludes.hasPair(iglob, jglob))
159	>	skipsForAtom[i].push_back(j);
160	>
161		if (oneTwo.hasPair(iglob, jglob)) {
162	<	toposForRowAtom[i].push_back(j);
163	<	topoDistRow[i][nTopos] = 1;
164	<	nTopos++;
162	>	toposForAtom[i].push_back(j);
163	>	topoDist[i].push_back(1);
164	>	} else {
165	>	if (oneThree.hasPair(iglob, jglob)) {
166	>	toposForAtom[i].push_back(j);
167	>	topoDist[i].push_back(2);
168	>	} else {
169	>	if (oneFour.hasPair(iglob, jglob)) {
170	>	toposForAtom[i].push_back(j);
171	>	topoDist[i].push_back(3);
172	>	}
173	>	}
174		}
171	–	if (oneThree.hasPair(iglob, jglob)) {
172	–	toposForRowAtom[i].push_back(j);
173	–	topoDistRow[i][nTopos] = 2;
174	–	nTopos++;
175	–	}
176	–	if (oneFour.hasPair(iglob, jglob)) {
177	–	toposForRowAtom[i].push_back(j);
178	–	topoDistRow[i][nTopos] = 3;
179	–	nTopos++;
180	–	}
175		}
176		}
177
178		#endif
179
180		groupList_.clear();
181	<	groupList_.reserve(nGroups_);
181	>	groupList_.resize(nGroups_);
182		for (int i = 0; i < nGroups_; i++) {
183		int gid = cgLocalToGlobal[i];
184		for (int j = 0; j < nLocal_; j++) {
185		int aid = AtomLocalToGlobal[j];
186	<	if (globalGroupMembership[aid] == gid)
186	>	if (globalGroupMembership[aid] == gid) {
187		groupList_[i].push_back(j);
188	+	}
189		}
190		}
191
192	<	skipsForLocalAtom.clear();
193	<	skipsForLocalAtom.reserve(nLocal_);
192	>	skipsForAtom.clear();
193	>	skipsForAtom.resize(nLocal_);
194	>	toposForAtom.clear();
195	>	toposForAtom.resize(nLocal_);
196	>	topoDist.clear();
197	>	topoDist.resize(nLocal_);
198
199		for (int i = 0; i < nLocal_; i++) {
200		int iglob = AtomLocalToGlobal[i];
201	+
202		for (int j = 0; j < nLocal_; j++) {
203	<	int jglob = AtomLocalToGlobal[j];
203	>	int jglob = AtomLocalToGlobal[j];
204	>
205		if (excludes.hasPair(iglob, jglob))
206	<	skipsForLocalAtom[i].push_back(j);
206	>	skipsForAtom[i].push_back(j);
207	>
208	>	if (oneTwo.hasPair(iglob, jglob)) {
209	>	toposForAtom[i].push_back(j);
210	>	topoDist[i].push_back(1);
211	>	} else {
212	>	if (oneThree.hasPair(iglob, jglob)) {
213	>	toposForAtom[i].push_back(j);
214	>	topoDist[i].push_back(2);
215	>	} else {
216	>	if (oneFour.hasPair(iglob, jglob)) {
217	>	toposForAtom[i].push_back(j);
218	>	topoDist[i].push_back(3);
219	>	}
220	>	}
221	>	}
222		}
223		}
224	+
225	+	createGtypeCutoffMap();
226	+	}
227	+
228	+	void ForceMatrixDecomposition::createGtypeCutoffMap() {
229
230	<	toposForLocalAtom.clear();
231	<	toposForLocalAtom.reserve(nLocal_);
232	<	for (int i = 0; i < nLocal_; i++) {
233	<	int iglob = AtomLocalToGlobal[i];
234	<	int nTopos = 0;
235	<	for (int j = 0; j < nLocal_; j++) {
236	<	int jglob = AtomLocalToGlobal[j];
237	<	if (oneTwo.hasPair(iglob, jglob)) {
238	<	toposForLocalAtom[i].push_back(j);
239	<	topoDistLocal[i][nTopos] = 1;
240	<	nTopos++;
241	<	}
242	<	if (oneThree.hasPair(iglob, jglob)) {
243	<	toposForLocalAtom[i].push_back(j);
244	<	topoDistLocal[i][nTopos] = 2;
245	<	nTopos++;
246	<	}
247	<	if (oneFour.hasPair(iglob, jglob)) {
248	<	toposForLocalAtom[i].push_back(j);
249	<	topoDistLocal[i][nTopos] = 3;
250	<	nTopos++;
230	>	RealType tol = 1e-6;
231	>	RealType rc;
232	>	int atid;
233	>	set<AtomType*> atypes = info_->getSimulatedAtomTypes();
234	>	vector<RealType> atypeCutoff;
235	>	atypeCutoff.resize( atypes.size() );
236	>
237	>	for (set<AtomType*>::iterator at = atypes.begin();
238	>	at != atypes.end(); ++at){
239	>	atid = (*at)->getIdent();
240	>
241	>	if (userChoseCutoff_)
242	>	atypeCutoff[atid] = userCutoff_;
243	>	else
244	>	atypeCutoff[atid] = interactionMan_->getSuggestedCutoffRadius(*at);
245	>	}
246	>
247	>	vector<RealType> gTypeCutoffs;
248	>
249	>	// first we do a single loop over the cutoff groups to find the
250	>	// largest cutoff for any atypes present in this group.
251	>	#ifdef IS_MPI
252	>	vector<RealType> groupCutoffRow(nGroupsInRow_, 0.0);
253	>	groupRowToGtype.resize(nGroupsInRow_);
254	>	for (int cg1 = 0; cg1 < nGroupsInRow_; cg1++) {
255	>	vector<int> atomListRow = getAtomsInGroupRow(cg1);
256	>	for (vector<int>::iterator ia = atomListRow.begin();
257	>	ia != atomListRow.end(); ++ia) {
258	>	int atom1 = (*ia);
259	>	atid = identsRow[atom1];
260	>	if (atypeCutoff[atid] > groupCutoffRow[cg1]) {
261	>	groupCutoffRow[cg1] = atypeCutoff[atid];
262		}
263	+	}
264	+
265	+	bool gTypeFound = false;
266	+	for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
267	+	if (abs(groupCutoffRow[cg1] - gTypeCutoffs[gt]) < tol) {
268	+	groupRowToGtype[cg1] = gt;
269	+	gTypeFound = true;
270	+	}
271	+	}
272	+	if (!gTypeFound) {
273	+	gTypeCutoffs.push_back( groupCutoffRow[cg1] );
274	+	groupRowToGtype[cg1] = gTypeCutoffs.size() - 1;
275	+	}
276	+
277	+	}
278	+	vector<RealType> groupCutoffCol(nGroupsInCol_, 0.0);
279	+	groupColToGtype.resize(nGroupsInCol_);
280	+	for (int cg2 = 0; cg2 < nGroupsInCol_; cg2++) {
281	+	vector<int> atomListCol = getAtomsInGroupColumn(cg2);
282	+	for (vector<int>::iterator jb = atomListCol.begin();
283	+	jb != atomListCol.end(); ++jb) {
284	+	int atom2 = (*jb);
285	+	atid = identsCol[atom2];
286	+	if (atypeCutoff[atid] > groupCutoffCol[cg2]) {
287	+	groupCutoffCol[cg2] = atypeCutoff[atid];
288	+	}
289	+	}
290	+	bool gTypeFound = false;
291	+	for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
292	+	if (abs(groupCutoffCol[cg2] - gTypeCutoffs[gt]) < tol) {
293	+	groupColToGtype[cg2] = gt;
294	+	gTypeFound = true;
295	+	}
296	+	}
297	+	if (!gTypeFound) {
298	+	gTypeCutoffs.push_back( groupCutoffCol[cg2] );
299	+	groupColToGtype[cg2] = gTypeCutoffs.size() - 1;
300	+	}
301	+	}
302	+	#else
303	+
304	+	vector<RealType> groupCutoff(nGroups_, 0.0);
305	+	groupToGtype.resize(nGroups_);
306	+
307	+	cerr << "nGroups = " << nGroups_ << "\n";
308	+	for (int cg1 = 0; cg1 < nGroups_; cg1++) {
309	+
310	+	groupCutoff[cg1] = 0.0;
311	+	vector<int> atomList = getAtomsInGroupRow(cg1);
312	+
313	+	for (vector<int>::iterator ia = atomList.begin();
314	+	ia != atomList.end(); ++ia) {
315	+	int atom1 = (*ia);
316	+	atid = idents[atom1];
317	+	if (atypeCutoff[atid] > groupCutoff[cg1]) {
318	+	groupCutoff[cg1] = atypeCutoff[atid];
319	+	}
320	+	}
321	+
322	+	bool gTypeFound = false;
323	+	for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
324	+	if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) {
325	+	groupToGtype[cg1] = gt;
326	+	gTypeFound = true;
327	+	}
328	+	}
329	+	if (!gTypeFound) {
330	+	gTypeCutoffs.push_back( groupCutoff[cg1] );
331	+	groupToGtype[cg1] = gTypeCutoffs.size() - 1;
332		}
333		}
334	+	#endif
335	+
336	+	cerr << "gTypeCutoffs.size() = " << gTypeCutoffs.size() << "\n";
337	+	// Now we find the maximum group cutoff value present in the simulation
338	+
339	+	RealType groupMax = *max_element(gTypeCutoffs.begin(), gTypeCutoffs.end());
340	+
341	+	#ifdef IS_MPI
342	+	MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, MPI::MAX);
343	+	#endif
344	+
345	+	RealType tradRcut = groupMax;
346	+
347	+	for (int i = 0; i < gTypeCutoffs.size(); i++) {
348	+	for (int j = 0; j < gTypeCutoffs.size(); j++) {
349	+	RealType thisRcut;
350	+	switch(cutoffPolicy_) {
351	+	case TRADITIONAL:
352	+	thisRcut = tradRcut;
353	+	break;
354	+	case MIX:
355	+	thisRcut = 0.5 * (gTypeCutoffs[i] + gTypeCutoffs[j]);
356	+	break;
357	+	case MAX:
358	+	thisRcut = max(gTypeCutoffs[i], gTypeCutoffs[j]);
359	+	break;
360	+	default:
361	+	sprintf(painCave.errMsg,
362	+	"ForceMatrixDecomposition::createGtypeCutoffMap "
363	+	"hit an unknown cutoff policy!\n");
364	+	painCave.severity = OPENMD_ERROR;
365	+	painCave.isFatal = 1;
366	+	simError();
367	+	break;
368	+	}
369	+
370	+	pair<int,int> key = make_pair(i,j);
371	+	gTypeCutoffMap[key].first = thisRcut;
372	+
373	+	if (thisRcut > largestRcut_) largestRcut_ = thisRcut;
374	+
375	+	gTypeCutoffMap[key].second = thisRcut*thisRcut;
376	+
377	+	gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2);
378	+
379	+	// sanity check
380	+
381	+	if (userChoseCutoff_) {
382	+	if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) {
383	+	sprintf(painCave.errMsg,
384	+	"ForceMatrixDecomposition::createGtypeCutoffMap "
385	+	"user-specified rCut (%lf) does not match computed group Cutoff\n", userCutoff_);
386	+	painCave.severity = OPENMD_ERROR;
387	+	painCave.isFatal = 1;
388	+	simError();
389	+	}
390	+	}
391	+	}
392	+	}
393		}
394	<
394	>
395	>
396	>	groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) {
397	>	int i, j;
398	>	#ifdef IS_MPI
399	>	i = groupRowToGtype[cg1];
400	>	j = groupColToGtype[cg2];
401	>	#else
402	>	i = groupToGtype[cg1];
403	>	j = groupToGtype[cg2];
404	>	#endif
405	>	return gTypeCutoffMap[make_pair(i,j)];
406	>	}
407	>
408	>	int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) {
409	>	for (int j = 0; j < toposForAtom[atom1].size(); j++) {
410	>	if (toposForAtom[atom1][j] == atom2)
411	>	return topoDist[atom1][j];
412	>	}
413	>	return 0;
414	>	}
415	>
416	>	void ForceMatrixDecomposition::zeroWorkArrays() {
417	>	pairwisePot = 0.0;
418	>	embeddingPot = 0.0;
419	>
420	>	#ifdef IS_MPI
421	>	if (storageLayout_ & DataStorage::dslForce) {
422	>	fill(atomRowData.force.begin(), atomRowData.force.end(), V3Zero);
423	>	fill(atomColData.force.begin(), atomColData.force.end(), V3Zero);
424	>	}
425	>
426	>	if (storageLayout_ & DataStorage::dslTorque) {
427	>	fill(atomRowData.torque.begin(), atomRowData.torque.end(), V3Zero);
428	>	fill(atomColData.torque.begin(), atomColData.torque.end(), V3Zero);
429	>	}
430	>
431	>	fill(pot_row.begin(), pot_row.end(),
432	>	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
433	>
434	>	fill(pot_col.begin(), pot_col.end(),
435	>	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
436	>
437	>	if (storageLayout_ & DataStorage::dslParticlePot) {
438	>	fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0);
439	>	fill(atomColData.particlePot.begin(), atomColData.particlePot.end(), 0.0);
440	>	}
441	>
442	>	if (storageLayout_ & DataStorage::dslDensity) {
443	>	fill(atomRowData.density.begin(), atomRowData.density.end(), 0.0);
444	>	fill(atomColData.density.begin(), atomColData.density.end(), 0.0);
445	>	}
446	>
447	>	if (storageLayout_ & DataStorage::dslFunctional) {
448	>	fill(atomRowData.functional.begin(), atomRowData.functional.end(), 0.0);
449	>	fill(atomColData.functional.begin(), atomColData.functional.end(), 0.0);
450	>	}
451	>
452	>	if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
453	>	fill(atomRowData.functionalDerivative.begin(),
454	>	atomRowData.functionalDerivative.end(), 0.0);
455	>	fill(atomColData.functionalDerivative.begin(),
456	>	atomColData.functionalDerivative.end(), 0.0);
457	>	}
458	>
459	>	#else
460	>
461	>	if (storageLayout_ & DataStorage::dslParticlePot) {
462	>	fill(snap_->atomData.particlePot.begin(),
463	>	snap_->atomData.particlePot.end(), 0.0);
464	>	}
465	>
466	>	if (storageLayout_ & DataStorage::dslDensity) {
467	>	fill(snap_->atomData.density.begin(),
468	>	snap_->atomData.density.end(), 0.0);
469	>	}
470	>	if (storageLayout_ & DataStorage::dslFunctional) {
471	>	fill(snap_->atomData.functional.begin(),
472	>	snap_->atomData.functional.end(), 0.0);
473	>	}
474	>	if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
475	>	fill(snap_->atomData.functionalDerivative.begin(),
476	>	snap_->atomData.functionalDerivative.end(), 0.0);
477	>	}
478	>	#endif
479	>
480	>	}
481	>
482	>
483		void ForceMatrixDecomposition::distributeData() {
484		snap_ = sman_->getCurrentSnapshot();
485		storageLayout_ = sman_->getStorageLayout();
#	Line 267 \| Line 515 \| namespace OpenMD {
515		#endif
516		}
517
518	+	/* collects information obtained during the pre-pair loop onto local
519	+	* data structures.
520	+	*/
521		void ForceMatrixDecomposition::collectIntermediateData() {
522		snap_ = sman_->getCurrentSnapshot();
523		storageLayout_ = sman_->getStorageLayout();
#	Line 278 \| Line 529 \| namespace OpenMD {
529		snap_->atomData.density);
530
531		int n = snap_->atomData.density.size();
532	<	std::vector<RealType> rho_tmp(n, 0.0);
532	>	vector<RealType> rho_tmp(n, 0.0);
533		AtomCommRealColumn->scatter(atomColData.density, rho_tmp);
534		for (int i = 0; i < n; i++)
535		snap_->atomData.density[i] += rho_tmp[i];
536		}
537		#endif
538		}
539	<
539	>
540	>	/*
541	>	* redistributes information obtained during the pre-pair loop out to
542	>	* row and column-indexed data structures
543	>	*/
544		void ForceMatrixDecomposition::distributeIntermediateData() {
545		snap_ = sman_->getCurrentSnapshot();
546		storageLayout_ = sman_->getStorageLayout();
#	Line 343 \| Line 598 \| namespace OpenMD {
598
599		nLocal_ = snap_->getNumberOfAtoms();
600
601	<	vector<vector<RealType> > pot_temp(N_INTERACTION_FAMILIES,
602	<	vector<RealType> (nLocal_, 0.0));
601	>	vector<potVec> pot_temp(nLocal_,
602	>	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
603	>
604	>	// scatter/gather pot_row into the members of my column
605	>
606	>	AtomCommPotRow->scatter(pot_row, pot_temp);
607	>
608	>	for (int ii = 0; ii < pot_temp.size(); ii++ )
609	>	pairwisePot += pot_temp[ii];
610
611	<	for (int i = 0; i < N_INTERACTION_FAMILIES; i++) {
612	<	AtomCommRealRow->scatter(pot_row[i], pot_temp[i]);
613	<	for (int ii = 0; ii < pot_temp[i].size(); ii++ ) {
614	<	pot_local[i] += pot_temp[i][ii];
615	<	}
616	<	}
611	>	fill(pot_temp.begin(), pot_temp.end(),
612	>	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
613	>
614	>	AtomCommPotColumn->scatter(pot_col, pot_temp);
615	>
616	>	for (int ii = 0; ii < pot_temp.size(); ii++ )
617	>	pairwisePot += pot_temp[ii];
618		#endif
619	+
620		}
621
622		int ForceMatrixDecomposition::getNAtomsInRow() {
#	Line 427 \| Line 691 \| namespace OpenMD {
691		#ifdef IS_MPI
692		return massFactorsRow[atom1];
693		#else
694	<	return massFactorsLocal[atom1];
694	>	cerr << "mfs = " << massFactors.size() << " atom1 = " << atom1 << "\n";
695	>	return massFactors[atom1];
696		#endif
697		}
698
#	Line 435 \| Line 700 \| namespace OpenMD {
700		#ifdef IS_MPI
701		return massFactorsCol[atom2];
702		#else
703	<	return massFactorsLocal[atom2];
703	>	return massFactors[atom2];
704		#endif
705
706		}
#	Line 453 \| Line 718 \| namespace OpenMD {
718		return d;
719		}
720
721	<	vector<int> ForceMatrixDecomposition::getSkipsForRowAtom(int atom1) {
722	<	#ifdef IS_MPI
458	<	return skipsForRowAtom[atom1];
459	<	#else
460	<	return skipsForLocalAtom[atom1];
461	<	#endif
721	>	vector<int> ForceMatrixDecomposition::getSkipsForAtom(int atom1) {
722	>	return skipsForAtom[atom1];
723		}
724
725		/**
726	<	* there are a number of reasons to skip a pair or a particle mostly
727	<	* we do this to exclude atoms who are involved in short range
728	<	* interactions (bonds, bends, torsions), but we also need to
729	<	* exclude some overcounted interactions that result from the
730	<	* parallel decomposition.
726	>	* There are a number of reasons to skip a pair or a
727	>	* particle. Mostly we do this to exclude atoms who are involved in
728	>	* short range interactions (bonds, bends, torsions), but we also
729	>	* need to exclude some overcounted interactions that result from
730	>	* the parallel decomposition.
731		*/
732		bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) {
733		int unique_id_1, unique_id_2;
#	Line 491 \| Line 752 \| namespace OpenMD {
752		unique_id_2 = atom2;
753		#endif
754
755	<	#ifdef IS_MPI
756	<	for (vector<int>::iterator i = skipsForRowAtom[atom1].begin();
496	<	i != skipsForRowAtom[atom1].end(); ++i) {
755	>	for (vector<int>::iterator i = skipsForAtom[atom1].begin();
756	>	i != skipsForAtom[atom1].end(); ++i) {
757		if ( (*i) == unique_id_2 ) return true;
498	–	}
499	–	#else
500	–	for (vector<int>::iterator i = skipsForLocalAtom[atom1].begin();
501	–	i != skipsForLocalAtom[atom1].end(); ++i) {
502	–	if ( (*i) == unique_id_2 ) return true;
503	–	}
504	–	#endif
505	–	}
506	–
507	–	int ForceMatrixDecomposition::getTopoDistance(int atom1, int atom2) {
508	–
509	–	#ifdef IS_MPI
510	–	for (int i = 0; i < toposForRowAtom[atom1].size(); i++) {
511	–	if ( toposForRowAtom[atom1][i] == atom2 ) return topoDistRow[atom1][i];
512	–	}
513	–	#else
514	–	for (int i = 0; i < toposForLocalAtom[atom1].size(); i++) {
515	–	if ( toposForLocalAtom[atom1][i] == atom2 ) return topoDistLocal[atom1][i];
758		}
517	–	#endif
759
760	<	// zero is default for unconnected (i.e. normal) pair interactions
520	<	return 0;
760	>	return false;
761		}
762
763	+
764		void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){
765		#ifdef IS_MPI
766		atomRowData.force[atom1] += fg;
#	Line 537 \| Line 778 \| namespace OpenMD {
778		}
779
780		// filling interaction blocks with pointers
781	<	InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) {
782	<	InteractionData idat;
542	<
781	>	void ForceMatrixDecomposition::fillInteractionData(InteractionData &idat,
782	>	int atom1, int atom2) {
783		#ifdef IS_MPI
784
785		idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
786		ff_->getAtomType(identsCol[atom2]) );
787	<
787	>
788		if (storageLayout_ & DataStorage::dslAmat) {
789		idat.A1 = &(atomRowData.aMat[atom1]);
790		idat.A2 = &(atomColData.aMat[atom2]);
#	Line 565 \| Line 805 \| namespace OpenMD {
805		idat.rho2 = &(atomColData.density[atom2]);
806		}
807
808	+	if (storageLayout_ & DataStorage::dslFunctional) {
809	+	idat.frho1 = &(atomRowData.functional[atom1]);
810	+	idat.frho2 = &(atomColData.functional[atom2]);
811	+	}
812	+
813		if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
814		idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]);
815		idat.dfrho2 = &(atomColData.functionalDerivative[atom2]);
816		}
817
818	+	if (storageLayout_ & DataStorage::dslParticlePot) {
819	+	idat.particlePot1 = &(atomRowData.particlePot[atom1]);
820	+	idat.particlePot2 = &(atomColData.particlePot[atom2]);
821	+	}
822	+
823		#else
824
825	<	idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]),
826	<	ff_->getAtomType(identsLocal[atom2]) );
825	>	idat.atypes = make_pair( ff_->getAtomType(idents[atom1]),
826	>	ff_->getAtomType(idents[atom2]) );
827
828		if (storageLayout_ & DataStorage::dslAmat) {
829		idat.A1 = &(snap_->atomData.aMat[atom1]);
#	Line 590 \| Line 840 \| namespace OpenMD {
840		idat.t2 = &(snap_->atomData.torque[atom2]);
841		}
842
843	<	if (storageLayout_ & DataStorage::dslDensity) {
843	>	if (storageLayout_ & DataStorage::dslDensity) {
844		idat.rho1 = &(snap_->atomData.density[atom1]);
845		idat.rho2 = &(snap_->atomData.density[atom2]);
846		}
847
848	+	if (storageLayout_ & DataStorage::dslFunctional) {
849	+	idat.frho1 = &(snap_->atomData.functional[atom1]);
850	+	idat.frho2 = &(snap_->atomData.functional[atom2]);
851	+	}
852	+
853		if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
854		idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]);
855		idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
856		}
857	+
858	+	if (storageLayout_ & DataStorage::dslParticlePot) {
859	+	idat.particlePot1 = &(snap_->atomData.particlePot[atom1]);
860	+	idat.particlePot2 = &(snap_->atomData.particlePot[atom2]);
861	+	}
862	+
863		#endif
603	–	return idat;
864		}
865
866	<	InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){
866	>
867	>	void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat, int atom1, int atom2) {
868	>	#ifdef IS_MPI
869	>	pot_row[atom1] += 0.5 * *(idat.pot);
870	>	pot_col[atom2] += 0.5 * *(idat.pot);
871
872	<	InteractionData idat;
872	>	atomRowData.force[atom1] += *(idat.f1);
873	>	atomColData.force[atom2] -= *(idat.f1);
874	>	#else
875	>	pairwisePot += *(idat.pot);
876	>
877	>	snap_->atomData.force[atom1] += *(idat.f1);
878	>	snap_->atomData.force[atom2] -= *(idat.f1);
879	>	#endif
880	>
881	>	}
882	>
883	>
884	>	void ForceMatrixDecomposition::fillSkipData(InteractionData &idat,
885	>	int atom1, int atom2) {
886		#ifdef IS_MPI
887		idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
888		ff_->getAtomType(identsCol[atom2]) );
#	Line 614 \| Line 891 \| namespace OpenMD {
891		idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
892		idat.eFrame2 = &(atomColData.electroFrame[atom2]);
893		}
894	+
895		if (storageLayout_ & DataStorage::dslTorque) {
896		idat.t1 = &(atomRowData.torque[atom1]);
897		idat.t2 = &(atomColData.torque[atom2]);
898		}
899	<	if (storageLayout_ & DataStorage::dslForce) {
900	<	idat.t1 = &(atomRowData.force[atom1]);
901	<	idat.t2 = &(atomColData.force[atom2]);
899	>
900	>	if (storageLayout_ & DataStorage::dslSkippedCharge) {
901	>	idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]);
902	>	idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]);
903		}
904		#else
905	<	idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]),
906	<	ff_->getAtomType(identsLocal[atom2]) );
905	>	idat.atypes = make_pair( ff_->getAtomType(idents[atom1]),
906	>	ff_->getAtomType(idents[atom2]) );
907
908		if (storageLayout_ & DataStorage::dslElectroFrame) {
909		idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
910		idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
911		}
912	+
913		if (storageLayout_ & DataStorage::dslTorque) {
914		idat.t1 = &(snap_->atomData.torque[atom1]);
915		idat.t2 = &(snap_->atomData.torque[atom2]);
916		}
917	<	if (storageLayout_ & DataStorage::dslForce) {
918	<	idat.t1 = &(snap_->atomData.force[atom1]);
919	<	idat.t2 = &(snap_->atomData.force[atom2]);
917	>
918	>	if (storageLayout_ & DataStorage::dslSkippedCharge) {
919	>	idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]);
920	>	idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]);
921		}
922		#endif
923		}
924
925	+
926	+	void ForceMatrixDecomposition::unpackSkipData(InteractionData &idat, int atom1, int atom2) {
927	+	#ifdef IS_MPI
928	+	pot_row[atom1] += 0.5 * *(idat.pot);
929	+	pot_col[atom2] += 0.5 * *(idat.pot);
930	+	#else
931	+	pairwisePot += *(idat.pot);
932	+	#endif
933	+
934	+	}
935	+
936	+
937		/*
938		* buildNeighborList
939		*
#	Line 650 \| Line 943 \| namespace OpenMD {
943		vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() {
944
945		vector<pair<int, int> > neighborList;
946	+	groupCutoffs cuts;
947		#ifdef IS_MPI
948		cellListRow_.clear();
949		cellListCol_.clear();
#	Line 657 \| Line 951 \| namespace OpenMD {
951		cellList_.clear();
952		#endif
953
954	<	// dangerous to not do error checking.
661	<	RealType rCut_;
662	<
663	<	RealType rList_ = (rCut_ + skinThickness_);
954	>	RealType rList_ = (largestRcut_ + skinThickness_);
955		RealType rl2 = rList_ * rList_;
956		Snapshot* snap_ = sman_->getCurrentSnapshot();
957		Mat3x3d Hmat = snap_->getHmat();
#	Line 676 \| Line 967 \| namespace OpenMD {
967		Vector3d rs, scaled, dr;
968		Vector3i whichCell;
969		int cellIndex;
970	+	int nCtot = nCells_.x() * nCells_.y() * nCells_.z();
971
972		#ifdef IS_MPI
973	+	cellListRow_.resize(nCtot);
974	+	cellListCol_.resize(nCtot);
975	+	#else
976	+	cellList_.resize(nCtot);
977	+	#endif
978	+
979	+	#ifdef IS_MPI
980		for (int i = 0; i < nGroupsInRow_; i++) {
981		rs = cgRowData.position[i];
982	+
983		// scaled positions relative to the box vectors
984		scaled = invHmat * rs;
985	+
986		// wrap the vector back into the unit box by subtracting integer box
987		// numbers
988	<	for (int j = 0; j < 3; j++)
988	>	for (int j = 0; j < 3; j++) {
989		scaled[j] -= roundMe(scaled[j]);
990	+	scaled[j] += 0.5;
991	+	}
992
993		// find xyz-indices of cell that cutoffGroup is in.
994		whichCell.x() = nCells_.x() * scaled.x();
#	Line 694 \| Line 997 \| namespace OpenMD {
997
998		// find single index of this cell:
999		cellIndex = Vlinear(whichCell, nCells_);
1000	+
1001		// add this cutoff group to the list of groups in this cell;
1002		cellListRow_[cellIndex].push_back(i);
1003		}
1004
1005		for (int i = 0; i < nGroupsInCol_; i++) {
1006		rs = cgColData.position[i];
1007	+
1008		// scaled positions relative to the box vectors
1009		scaled = invHmat * rs;
1010	+
1011		// wrap the vector back into the unit box by subtracting integer box
1012		// numbers
1013	<	for (int j = 0; j < 3; j++)
1013	>	for (int j = 0; j < 3; j++) {
1014		scaled[j] -= roundMe(scaled[j]);
1015	+	scaled[j] += 0.5;
1016	+	}
1017
1018		// find xyz-indices of cell that cutoffGroup is in.
1019		whichCell.x() = nCells_.x() * scaled.x();
#	Line 714 \| Line 1022 \| namespace OpenMD {
1022
1023		// find single index of this cell:
1024		cellIndex = Vlinear(whichCell, nCells_);
1025	+
1026		// add this cutoff group to the list of groups in this cell;
1027		cellListCol_[cellIndex].push_back(i);
1028		}
1029		#else
1030		for (int i = 0; i < nGroups_; i++) {
1031		rs = snap_->cgData.position[i];
1032	+
1033		// scaled positions relative to the box vectors
1034		scaled = invHmat * rs;
1035	+
1036		// wrap the vector back into the unit box by subtracting integer box
1037		// numbers
1038	<	for (int j = 0; j < 3; j++)
1038	>	for (int j = 0; j < 3; j++) {
1039		scaled[j] -= roundMe(scaled[j]);
1040	+	scaled[j] += 0.5;
1041	+	}
1042
1043		// find xyz-indices of cell that cutoffGroup is in.
1044		whichCell.x() = nCells_.x() * scaled.x();
#	Line 733 \| Line 1046 \| namespace OpenMD {
1046		whichCell.z() = nCells_.z() * scaled.z();
1047
1048		// find single index of this cell:
1049	<	cellIndex = Vlinear(whichCell, nCells_);
1049	>	cellIndex = Vlinear(whichCell, nCells_);
1050	>
1051		// add this cutoff group to the list of groups in this cell;
1052		cellList_[cellIndex].push_back(i);
1053		}
1054		#endif
1055
742	–
743	–
1056		for (int m1z = 0; m1z < nCells_.z(); m1z++) {
1057		for (int m1y = 0; m1y < nCells_.y(); m1y++) {
1058		for (int m1x = 0; m1x < nCells_.x(); m1x++) {
#	Line 785 \| Line 1097 \| namespace OpenMD {
1097		if (m2 != m1 \|\| cgColToGlobal[(j2)] < cgRowToGlobal[(j1)]) {
1098		dr = cgColData.position[(j2)] - cgRowData.position[(j1)];
1099		snap_->wrapVector(dr);
1100	<	if (dr.lengthSquare() < rl2) {
1100	>	cuts = getGroupCutoffs( (j1), (j2) );
1101	>	if (dr.lengthSquare() < cuts.third) {
1102		neighborList.push_back(make_pair((j1), (j2)));
1103		}
1104		}
1105		}
1106		}
1107		#else
1108	+
1109		for (vector<int>::iterator j1 = cellList_[m1].begin();
1110		j1 != cellList_[m1].end(); ++j1) {
1111		for (vector<int>::iterator j2 = cellList_[m2].begin();
1112		j2 != cellList_[m2].end(); ++j2) {
1113	<
1113	>
1114		// Always do this if we're in different cells or if
1115		// we're in the same cell and the global index of the
1116		// j2 cutoff group is less than the j1 cutoff group
#	Line 804 \| Line 1118 \| namespace OpenMD {
1118		if (m2 != m1 \|\| (j2) < (j1)) {
1119		dr = snap_->cgData.position[(j2)] - snap_->cgData.position[(j1)];
1120		snap_->wrapVector(dr);
1121	<	if (dr.lengthSquare() < rl2) {
1121	>	cuts = getGroupCutoffs( (j1), (j2) );
1122	>	if (dr.lengthSquare() < cuts.third) {
1123		neighborList.push_back(make_pair((j1), (j2)));
1124		}
1125		}
#	Line 815 \| Line 1130 \| namespace OpenMD {
1130		}
1131		}
1132		}
1133	<
1133	>
1134		// save the local cutoff group positions for the check that is
1135		// done on each loop:
1136		saved_CG_positions_.clear();
1137		for (int i = 0; i < nGroups_; i++)
1138		saved_CG_positions_.push_back(snap_->cgData.position[i]);
1139	<
1139	>
1140		return neighborList;
1141		}
1142		} //end namespace OpenMD

Diff Legend

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

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents): Revision 1571 by gezelter, Fri May 27 16:45:44 2011 UTC vs. Revision 1584 by gezelter, Fri Jun 17 20:16:35 2011 UTC

Diff Legend

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1571 by gezelter, Fri May 27 16:45:44 2011 UTC vs.
Revision 1584 by gezelter, Fri Jun 17 20:16:35 2011 UTC