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

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1568 by gezelter, Wed May 25 16:20:37 2011 UTC vs.
Revision 1582 by gezelter, Tue Jun 14 20:41:44 2011 UTC

#	Line 42 \| Line 42
42		#include "math/SquareMatrix3.hpp"
43		#include "nonbonded/NonBondedInteraction.hpp"
44		#include "brains/SnapshotManager.hpp"
45	+	#include "brains/PairList.hpp"
46
47		using namespace std;
48		namespace OpenMD {
#	Line 54 \| Line 55 \| namespace OpenMD {
55		void ForceMatrixDecomposition::distributeInitialData() {
56		snap_ = sman_->getCurrentSnapshot();
57		storageLayout_ = sman_->getStorageLayout();
58	+	ff_ = info_->getForceField();
59		nLocal_ = snap_->getNumberOfAtoms();
58	–	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();
65	+	AtomLocalToGlobal = info_->getGlobalAtomIndices();
66	+	cgLocalToGlobal = info_->getGlobalGroupIndices();
67	+	vector<int> globalGroupMembership = info_->getGlobalGroupMembership();
68	+	massFactors = info_->getMassFactors();
69	+	PairList excludes = info_->getExcludedInteractions();
70	+	PairList oneTwo = info_->getOneTwoInteractions();
71	+	PairList oneThree = info_->getOneThreeInteractions();
72	+	PairList oneFour = info_->getOneFourInteractions();
73	+
74		#ifdef IS_MPI
75
76		AtomCommIntRow = new Communicator<Row,int>(nLocal_);
77		AtomCommRealRow = new Communicator<Row,RealType>(nLocal_);
78		AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal_);
79		AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal_);
80	+	AtomCommPotRow = new Communicator<Row,potVec>(nLocal_);
81
82		AtomCommIntColumn = new Communicator<Column,int>(nLocal_);
83		AtomCommRealColumn = new Communicator<Column,RealType>(nLocal_);
84		AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal_);
85		AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal_);
86	+	AtomCommPotColumn = new Communicator<Column,potVec>(nLocal_);
87
88		cgCommIntRow = new Communicator<Row,int>(nGroups_);
89		cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups_);
#	Line 88 \| Line 104 \| namespace OpenMD {
104		cgRowData.setStorageLayout(DataStorage::dslPosition);
105		cgColData.resize(nGroupsInCol_);
106		cgColData.setStorageLayout(DataStorage::dslPosition);
107	+
108	+	identsRow.resize(nAtomsInRow_);
109	+	identsCol.resize(nAtomsInCol_);
110
92	–	vector<vector<RealType> > pot_row(N_INTERACTION_FAMILIES,
93	–	vector<RealType> (nAtomsInRow_, 0.0));
94	–	vector<vector<RealType> > pot_col(N_INTERACTION_FAMILIES,
95	–	vector<RealType> (nAtomsInCol_, 0.0));
96	–
97	–
98	–	vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0);
99	–
100	–	// gather the information for atomtype IDs (atids):
101	–	vector<int> identsLocal = info_->getIdentArray();
102	–	identsRow.reserve(nAtomsInRow_);
103	–	identsCol.reserve(nAtomsInCol_);
104	–
111		AtomCommIntRow->gather(identsLocal, identsRow);
112		AtomCommIntColumn->gather(identsLocal, identsCol);
113
108	–	AtomLocalToGlobal = info_->getGlobalAtomIndices();
114		AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal);
115		AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal);
116
112	–	cgLocalToGlobal = info_->getGlobalGroupIndices();
117		cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal);
118		cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal);
119
120	<	// still need:
121	<	// topoDist
122	<	// exclude
120	>	AtomCommRealRow->gather(massFactors, massFactorsRow);
121	>	AtomCommRealColumn->gather(massFactors, massFactorsCol);
122	>
123	>	groupListRow_.clear();
124	>	groupListRow_.resize(nGroupsInRow_);
125	>	for (int i = 0; i < nGroupsInRow_; i++) {
126	>	int gid = cgRowToGlobal[i];
127	>	for (int j = 0; j < nAtomsInRow_; j++) {
128	>	int aid = AtomRowToGlobal[j];
129	>	if (globalGroupMembership[aid] == gid)
130	>	groupListRow_[i].push_back(j);
131	>	}
132	>	}
133	>
134	>	groupListCol_.clear();
135	>	groupListCol_.resize(nGroupsInCol_);
136	>	for (int i = 0; i < nGroupsInCol_; i++) {
137	>	int gid = cgColToGlobal[i];
138	>	for (int j = 0; j < nAtomsInCol_; j++) {
139	>	int aid = AtomColToGlobal[j];
140	>	if (globalGroupMembership[aid] == gid)
141	>	groupListCol_[i].push_back(j);
142	>	}
143	>	}
144	>
145	>	skipsForAtom.clear();
146	>	skipsForAtom.resize(nAtomsInRow_);
147	>	toposForAtom.clear();
148	>	toposForAtom.resize(nAtomsInRow_);
149	>	topoDist.clear();
150	>	topoDist.resize(nAtomsInRow_);
151	>	for (int i = 0; i < nAtomsInRow_; i++) {
152	>	int iglob = AtomRowToGlobal[i];
153	>
154	>	for (int j = 0; j < nAtomsInCol_; j++) {
155	>	int jglob = AtomColToGlobal[j];
156	>
157	>	if (excludes.hasPair(iglob, jglob))
158	>	skipsForAtom[i].push_back(j);
159	>
160	>	if (oneTwo.hasPair(iglob, jglob)) {
161	>	toposForAtom[i].push_back(j);
162	>	topoDist[i].push_back(1);
163	>	} else {
164	>	if (oneThree.hasPair(iglob, jglob)) {
165	>	toposForAtom[i].push_back(j);
166	>	topoDist[i].push_back(2);
167	>	} else {
168	>	if (oneFour.hasPair(iglob, jglob)) {
169	>	toposForAtom[i].push_back(j);
170	>	topoDist[i].push_back(3);
171	>	}
172	>	}
173	>	}
174	>	}
175	>	}
176	>
177		#endif
178	+
179	+	groupList_.clear();
180	+	groupList_.resize(nGroups_);
181	+	for (int i = 0; i < nGroups_; i++) {
182	+	int gid = cgLocalToGlobal[i];
183	+	for (int j = 0; j < nLocal_; j++) {
184	+	int aid = AtomLocalToGlobal[j];
185	+	if (globalGroupMembership[aid] == gid) {
186	+	groupList_[i].push_back(j);
187	+	}
188	+	}
189	+	}
190	+
191	+	skipsForAtom.clear();
192	+	skipsForAtom.resize(nLocal_);
193	+	toposForAtom.clear();
194	+	toposForAtom.resize(nLocal_);
195	+	topoDist.clear();
196	+	topoDist.resize(nLocal_);
197	+
198	+	for (int i = 0; i < nLocal_; i++) {
199	+	int iglob = AtomLocalToGlobal[i];
200	+
201	+	for (int j = 0; j < nLocal_; j++) {
202	+	int jglob = AtomLocalToGlobal[j];
203	+
204	+	if (excludes.hasPair(iglob, jglob))
205	+	skipsForAtom[i].push_back(j);
206	+
207	+	if (oneTwo.hasPair(iglob, jglob)) {
208	+	toposForAtom[i].push_back(j);
209	+	topoDist[i].push_back(1);
210	+	} else {
211	+	if (oneThree.hasPair(iglob, jglob)) {
212	+	toposForAtom[i].push_back(j);
213	+	topoDist[i].push_back(2);
214	+	} else {
215	+	if (oneFour.hasPair(iglob, jglob)) {
216	+	toposForAtom[i].push_back(j);
217	+	topoDist[i].push_back(3);
218	+	}
219	+	}
220	+	}
221	+	}
222	+	}
223	+
224	+	createGtypeCutoffMap();
225		}
226	+
227	+	void ForceMatrixDecomposition::createGtypeCutoffMap() {
228	+
229	+	RealType tol = 1e-6;
230	+	RealType rc;
231	+	int atid;
232	+	set<AtomType*> atypes = info_->getSimulatedAtomTypes();
233	+	vector<RealType> atypeCutoff;
234	+	atypeCutoff.resize( atypes.size() );
235	+
236	+	for (set<AtomType*>::iterator at = atypes.begin();
237	+	at != atypes.end(); ++at){
238	+	rc = interactionMan_->getSuggestedCutoffRadius(*at);
239	+	atid = (*at)->getIdent();
240	+	atypeCutoff[atid] = rc;
241	+	}
242	+
243	+	vector<RealType> gTypeCutoffs;
244	+
245	+	// first we do a single loop over the cutoff groups to find the
246	+	// largest cutoff for any atypes present in this group.
247	+	#ifdef IS_MPI
248	+	vector<RealType> groupCutoffRow(nGroupsInRow_, 0.0);
249	+	groupRowToGtype.resize(nGroupsInRow_);
250	+	for (int cg1 = 0; cg1 < nGroupsInRow_; cg1++) {
251	+	vector<int> atomListRow = getAtomsInGroupRow(cg1);
252	+	for (vector<int>::iterator ia = atomListRow.begin();
253	+	ia != atomListRow.end(); ++ia) {
254	+	int atom1 = (*ia);
255	+	atid = identsRow[atom1];
256	+	if (atypeCutoff[atid] > groupCutoffRow[cg1]) {
257	+	groupCutoffRow[cg1] = atypeCutoff[atid];
258	+	}
259	+	}
260	+
261	+	bool gTypeFound = false;
262	+	for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
263	+	if (abs(groupCutoffRow[cg1] - gTypeCutoffs[gt]) < tol) {
264	+	groupRowToGtype[cg1] = gt;
265	+	gTypeFound = true;
266	+	}
267	+	}
268	+	if (!gTypeFound) {
269	+	gTypeCutoffs.push_back( groupCutoffRow[cg1] );
270	+	groupRowToGtype[cg1] = gTypeCutoffs.size() - 1;
271	+	}
272	+
273	+	}
274	+	vector<RealType> groupCutoffCol(nGroupsInCol_, 0.0);
275	+	groupColToGtype.resize(nGroupsInCol_);
276	+	for (int cg2 = 0; cg2 < nGroupsInCol_; cg2++) {
277	+	vector<int> atomListCol = getAtomsInGroupColumn(cg2);
278	+	for (vector<int>::iterator jb = atomListCol.begin();
279	+	jb != atomListCol.end(); ++jb) {
280	+	int atom2 = (*jb);
281	+	atid = identsCol[atom2];
282	+	if (atypeCutoff[atid] > groupCutoffCol[cg2]) {
283	+	groupCutoffCol[cg2] = atypeCutoff[atid];
284	+	}
285	+	}
286	+	bool gTypeFound = false;
287	+	for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
288	+	if (abs(groupCutoffCol[cg2] - gTypeCutoffs[gt]) < tol) {
289	+	groupColToGtype[cg2] = gt;
290	+	gTypeFound = true;
291	+	}
292	+	}
293	+	if (!gTypeFound) {
294	+	gTypeCutoffs.push_back( groupCutoffCol[cg2] );
295	+	groupColToGtype[cg2] = gTypeCutoffs.size() - 1;
296	+	}
297	+	}
298	+	#else
299	+
300	+	vector<RealType> groupCutoff(nGroups_, 0.0);
301	+	groupToGtype.resize(nGroups_);
302	+
303	+	cerr << "nGroups = " << nGroups_ << "\n";
304	+	for (int cg1 = 0; cg1 < nGroups_; cg1++) {
305	+
306	+	groupCutoff[cg1] = 0.0;
307	+	vector<int> atomList = getAtomsInGroupRow(cg1);
308	+
309	+	for (vector<int>::iterator ia = atomList.begin();
310	+	ia != atomList.end(); ++ia) {
311	+	int atom1 = (*ia);
312	+	atid = identsLocal[atom1];
313	+	if (atypeCutoff[atid] > groupCutoff[cg1]) {
314	+	groupCutoff[cg1] = atypeCutoff[atid];
315	+	}
316	+	}
317	+
318	+	bool gTypeFound = false;
319	+	for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
320	+	if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) {
321	+	groupToGtype[cg1] = gt;
322	+	gTypeFound = true;
323	+	}
324	+	}
325	+	if (!gTypeFound) {
326	+	gTypeCutoffs.push_back( groupCutoff[cg1] );
327	+	groupToGtype[cg1] = gTypeCutoffs.size() - 1;
328	+	}
329	+	}
330	+	#endif
331	+
332	+	cerr << "gTypeCutoffs.size() = " << gTypeCutoffs.size() << "\n";
333	+	// Now we find the maximum group cutoff value present in the simulation
334	+
335	+	RealType groupMax = *max_element(gTypeCutoffs.begin(), gTypeCutoffs.end());
336	+
337	+	#ifdef IS_MPI
338	+	MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, MPI::MAX);
339	+	#endif
340
341	+	RealType tradRcut = groupMax;
342
343	+	for (int i = 0; i < gTypeCutoffs.size(); i++) {
344	+	for (int j = 0; j < gTypeCutoffs.size(); j++) {
345	+	RealType thisRcut;
346	+	switch(cutoffPolicy_) {
347	+	case TRADITIONAL:
348	+	thisRcut = tradRcut;
349	+	break;
350	+	case MIX:
351	+	thisRcut = 0.5 * (gTypeCutoffs[i] + gTypeCutoffs[j]);
352	+	break;
353	+	case MAX:
354	+	thisRcut = max(gTypeCutoffs[i], gTypeCutoffs[j]);
355	+	break;
356	+	default:
357	+	sprintf(painCave.errMsg,
358	+	"ForceMatrixDecomposition::createGtypeCutoffMap "
359	+	"hit an unknown cutoff policy!\n");
360	+	painCave.severity = OPENMD_ERROR;
361	+	painCave.isFatal = 1;
362	+	simError();
363	+	break;
364	+	}
365
366	+	pair<int,int> key = make_pair(i,j);
367	+	gTypeCutoffMap[key].first = thisRcut;
368	+
369	+	if (thisRcut > largestRcut_) largestRcut_ = thisRcut;
370	+
371	+	gTypeCutoffMap[key].second = thisRcut*thisRcut;
372	+
373	+	gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2);
374	+
375	+	// sanity check
376	+
377	+	if (userChoseCutoff_) {
378	+	if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) {
379	+	sprintf(painCave.errMsg,
380	+	"ForceMatrixDecomposition::createGtypeCutoffMap "
381	+	"user-specified rCut does not match computed group Cutoff\n");
382	+	painCave.severity = OPENMD_ERROR;
383	+	painCave.isFatal = 1;
384	+	simError();
385	+	}
386	+	}
387	+	}
388	+	}
389	+	}
390	+
391	+
392	+	groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) {
393	+	int i, j;
394	+	#ifdef IS_MPI
395	+	i = groupRowToGtype[cg1];
396	+	j = groupColToGtype[cg2];
397	+	#else
398	+	i = groupToGtype[cg1];
399	+	j = groupToGtype[cg2];
400	+	#endif
401	+	return gTypeCutoffMap[make_pair(i,j)];
402	+	}
403	+
404	+	int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) {
405	+	for (int j = 0; j < toposForAtom[atom1].size(); j++) {
406	+	if (toposForAtom[atom1][j] == atom2)
407	+	return topoDist[atom1][j];
408	+	}
409	+	return 0;
410	+	}
411	+
412	+	void ForceMatrixDecomposition::zeroWorkArrays() {
413	+
414	+	for (int j = 0; j < N_INTERACTION_FAMILIES; j++) {
415	+	longRangePot_[j] = 0.0;
416	+	}
417	+
418	+	#ifdef IS_MPI
419	+	if (storageLayout_ & DataStorage::dslForce) {
420	+	fill(atomRowData.force.begin(), atomRowData.force.end(), V3Zero);
421	+	fill(atomColData.force.begin(), atomColData.force.end(), V3Zero);
422	+	}
423	+
424	+	if (storageLayout_ & DataStorage::dslTorque) {
425	+	fill(atomRowData.torque.begin(), atomRowData.torque.end(), V3Zero);
426	+	fill(atomColData.torque.begin(), atomColData.torque.end(), V3Zero);
427	+	}
428	+
429	+	fill(pot_row.begin(), pot_row.end(),
430	+	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
431	+
432	+	fill(pot_col.begin(), pot_col.end(),
433	+	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
434	+
435	+	pot_local = 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 156 \| 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 167 \| 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 232 \| 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	>	pot_local += 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	>	pot_local += pot_temp[ii];
618	>
619		#endif
620		}
621
622	+	int ForceMatrixDecomposition::getNAtomsInRow() {
623	+	#ifdef IS_MPI
624	+	return nAtomsInRow_;
625	+	#else
626	+	return nLocal_;
627	+	#endif
628	+	}
629	+
630	+	/**
631	+	* returns the list of atoms belonging to this group.
632	+	*/
633	+	vector<int> ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){
634	+	#ifdef IS_MPI
635	+	return groupListRow_[cg1];
636	+	#else
637	+	return groupList_[cg1];
638	+	#endif
639	+	}
640	+
641	+	vector<int> ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){
642	+	#ifdef IS_MPI
643	+	return groupListCol_[cg2];
644	+	#else
645	+	return groupList_[cg2];
646	+	#endif
647	+	}
648
649		Vector3d ForceMatrixDecomposition::getIntergroupVector(int cg1, int cg2){
650		Vector3d d;
#	Line 285 \| Line 686 \| namespace OpenMD {
686		snap_->wrapVector(d);
687		return d;
688		}
689	+
690	+	RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) {
691	+	#ifdef IS_MPI
692	+	return massFactorsRow[atom1];
693	+	#else
694	+	return massFactors[atom1];
695	+	#endif
696	+	}
697	+
698	+	RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) {
699	+	#ifdef IS_MPI
700	+	return massFactorsCol[atom2];
701	+	#else
702	+	return massFactors[atom2];
703	+	#endif
704	+
705	+	}
706
707		Vector3d ForceMatrixDecomposition::getInteratomicVector(int atom1, int atom2){
708		Vector3d d;
#	Line 297 \| Line 715 \| namespace OpenMD {
715
716		snap_->wrapVector(d);
717		return d;
718	+	}
719	+
720	+	vector<int> ForceMatrixDecomposition::getSkipsForAtom(int atom1) {
721	+	return skipsForAtom[atom1];
722	+	}
723	+
724	+	/**
725	+	* There are a number of reasons to skip a pair or a
726	+	* particle. Mostly we do this to exclude atoms who are involved in
727	+	* short range interactions (bonds, bends, torsions), but we also
728	+	* need to exclude some overcounted interactions that result from
729	+	* the parallel decomposition.
730	+	*/
731	+	bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) {
732	+	int unique_id_1, unique_id_2;
733	+
734	+	#ifdef IS_MPI
735	+	// in MPI, we have to look up the unique IDs for each atom
736	+	unique_id_1 = AtomRowToGlobal[atom1];
737	+	unique_id_2 = AtomColToGlobal[atom2];
738	+
739	+	// this situation should only arise in MPI simulations
740	+	if (unique_id_1 == unique_id_2) return true;
741	+
742	+	// this prevents us from doing the pair on multiple processors
743	+	if (unique_id_1 < unique_id_2) {
744	+	if ((unique_id_1 + unique_id_2) % 2 == 0) return true;
745	+	} else {
746	+	if ((unique_id_1 + unique_id_2) % 2 == 1) return true;
747	+	}
748	+	#else
749	+	// in the normal loop, the atom numbers are unique
750	+	unique_id_1 = atom1;
751	+	unique_id_2 = atom2;
752	+	#endif
753	+
754	+	for (vector<int>::iterator i = skipsForAtom[atom1].begin();
755	+	i != skipsForAtom[atom1].end(); ++i) {
756	+	if ( (*i) == unique_id_2 ) return true;
757	+	}
758	+
759		}
760
761	+
762		void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){
763		#ifdef IS_MPI
764		atomRowData.force[atom1] += fg;
#	Line 316 \| Line 776 \| namespace OpenMD {
776		}
777
778		// filling interaction blocks with pointers
779	<	InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) {
780	<	InteractionData idat;
321	<
779	>	void ForceMatrixDecomposition::fillInteractionData(InteractionData &idat,
780	>	int atom1, int atom2) {
781		#ifdef IS_MPI
782	+
783	+	idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
784	+	ff_->getAtomType(identsCol[atom2]) );
785	+
786		if (storageLayout_ & DataStorage::dslAmat) {
787		idat.A1 = &(atomRowData.aMat[atom1]);
788		idat.A2 = &(atomColData.aMat[atom2]);
#	Line 340 \| Line 803 \| namespace OpenMD {
803		idat.rho2 = &(atomColData.density[atom2]);
804		}
805
806	+	if (storageLayout_ & DataStorage::dslFunctional) {
807	+	idat.frho1 = &(atomRowData.functional[atom1]);
808	+	idat.frho2 = &(atomColData.functional[atom2]);
809	+	}
810	+
811		if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
812		idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]);
813		idat.dfrho2 = &(atomColData.functionalDerivative[atom2]);
814		}
815	+
816	+	if (storageLayout_ & DataStorage::dslParticlePot) {
817	+	idat.particlePot1 = &(atomRowData.particlePot[atom1]);
818	+	idat.particlePot2 = &(atomColData.particlePot[atom2]);
819	+	}
820	+
821		#else
822	+
823	+	idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]),
824	+	ff_->getAtomType(identsLocal[atom2]) );
825	+
826		if (storageLayout_ & DataStorage::dslAmat) {
827		idat.A1 = &(snap_->atomData.aMat[atom1]);
828		idat.A2 = &(snap_->atomData.aMat[atom2]);
#	Line 365 \| Line 843 \| namespace OpenMD {
843		idat.rho2 = &(snap_->atomData.density[atom2]);
844		}
845
846	+	if (storageLayout_ & DataStorage::dslFunctional) {
847	+	idat.frho1 = &(snap_->atomData.functional[atom1]);
848	+	idat.frho2 = &(snap_->atomData.functional[atom2]);
849	+	}
850	+
851		if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
852		idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]);
853		idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
854		}
855	+
856	+	if (storageLayout_ & DataStorage::dslParticlePot) {
857	+	idat.particlePot1 = &(snap_->atomData.particlePot[atom1]);
858	+	idat.particlePot2 = &(snap_->atomData.particlePot[atom2]);
859	+	}
860	+
861		#endif
373	–	return idat;
862		}
863
864	<	InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){
864	>
865	>	void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat, int atom1, int atom2) {
866	>	#ifdef IS_MPI
867	>	pot_row[atom1] += 0.5 * *(idat.pot);
868	>	pot_col[atom2] += 0.5 * *(idat.pot);
869
870	<	InteractionData idat;
870	>	atomRowData.force[atom1] += *(idat.f1);
871	>	atomColData.force[atom2] -= *(idat.f1);
872	>	#else
873	>	longRangePot_ += *(idat.pot);
874	>
875	>	snap_->atomData.force[atom1] += *(idat.f1);
876	>	snap_->atomData.force[atom2] -= *(idat.f1);
877	>	#endif
878	>
879	>	}
880	>
881	>
882	>	void ForceMatrixDecomposition::fillSkipData(InteractionData &idat,
883	>	int atom1, int atom2) {
884		#ifdef IS_MPI
885	+	idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
886	+	ff_->getAtomType(identsCol[atom2]) );
887	+
888		if (storageLayout_ & DataStorage::dslElectroFrame) {
889		idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
890		idat.eFrame2 = &(atomColData.electroFrame[atom2]);
#	Line 385 \| Line 893 \| namespace OpenMD {
893		idat.t1 = &(atomRowData.torque[atom1]);
894		idat.t2 = &(atomColData.torque[atom2]);
895		}
388	–	if (storageLayout_ & DataStorage::dslForce) {
389	–	idat.t1 = &(atomRowData.force[atom1]);
390	–	idat.t2 = &(atomColData.force[atom2]);
391	–	}
896		#else
897	+	idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]),
898	+	ff_->getAtomType(identsLocal[atom2]) );
899	+
900		if (storageLayout_ & DataStorage::dslElectroFrame) {
901		idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
902		idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
#	Line 398 \| Line 905 \| namespace OpenMD {
905		idat.t1 = &(snap_->atomData.torque[atom1]);
906		idat.t2 = &(snap_->atomData.torque[atom2]);
907		}
908	<	if (storageLayout_ & DataStorage::dslForce) {
402	<	idat.t1 = &(snap_->atomData.force[atom1]);
403	<	idat.t2 = &(snap_->atomData.force[atom2]);
404	<	}
405	<	#endif
406	<
908	>	#endif
909		}
910
409	–
410	–
411	–
911		/*
912		* buildNeighborList
913		*
#	Line 418 \| Line 917 \| namespace OpenMD {
917		vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() {
918
919		vector<pair<int, int> > neighborList;
920	+	groupCutoffs cuts;
921		#ifdef IS_MPI
922		cellListRow_.clear();
923		cellListCol_.clear();
#	Line 425 \| Line 925 \| namespace OpenMD {
925		cellList_.clear();
926		#endif
927
928	<	// dangerous to not do error checking.
429	<	RealType rCut_;
430	<
431	<	RealType rList_ = (rCut_ + skinThickness_);
928	>	RealType rList_ = (largestRcut_ + skinThickness_);
929		RealType rl2 = rList_ * rList_;
930		Snapshot* snap_ = sman_->getCurrentSnapshot();
931		Mat3x3d Hmat = snap_->getHmat();
#	Line 444 \| Line 941 \| namespace OpenMD {
941		Vector3d rs, scaled, dr;
942		Vector3i whichCell;
943		int cellIndex;
944	+	int nCtot = nCells_.x() * nCells_.y() * nCells_.z();
945
946		#ifdef IS_MPI
947	+	cellListRow_.resize(nCtot);
948	+	cellListCol_.resize(nCtot);
949	+	#else
950	+	cellList_.resize(nCtot);
951	+	#endif
952	+
953	+	#ifdef IS_MPI
954		for (int i = 0; i < nGroupsInRow_; i++) {
955		rs = cgRowData.position[i];
956	+
957		// scaled positions relative to the box vectors
958		scaled = invHmat * rs;
959	+
960		// wrap the vector back into the unit box by subtracting integer box
961		// numbers
962	<	for (int j = 0; j < 3; j++)
962	>	for (int j = 0; j < 3; j++) {
963		scaled[j] -= roundMe(scaled[j]);
964	+	scaled[j] += 0.5;
965	+	}
966
967		// find xyz-indices of cell that cutoffGroup is in.
968		whichCell.x() = nCells_.x() * scaled.x();
#	Line 462 \| Line 971 \| namespace OpenMD {
971
972		// find single index of this cell:
973		cellIndex = Vlinear(whichCell, nCells_);
974	+
975		// add this cutoff group to the list of groups in this cell;
976		cellListRow_[cellIndex].push_back(i);
977		}
978
979		for (int i = 0; i < nGroupsInCol_; i++) {
980		rs = cgColData.position[i];
981	+
982		// scaled positions relative to the box vectors
983		scaled = invHmat * rs;
984	+
985		// wrap the vector back into the unit box by subtracting integer box
986		// numbers
987	<	for (int j = 0; j < 3; j++)
987	>	for (int j = 0; j < 3; j++) {
988		scaled[j] -= roundMe(scaled[j]);
989	+	scaled[j] += 0.5;
990	+	}
991
992		// find xyz-indices of cell that cutoffGroup is in.
993		whichCell.x() = nCells_.x() * scaled.x();
#	Line 482 \| Line 996 \| namespace OpenMD {
996
997		// find single index of this cell:
998		cellIndex = Vlinear(whichCell, nCells_);
999	+
1000		// add this cutoff group to the list of groups in this cell;
1001		cellListCol_[cellIndex].push_back(i);
1002		}
1003		#else
1004		for (int i = 0; i < nGroups_; i++) {
1005		rs = snap_->cgData.position[i];
1006	+
1007		// scaled positions relative to the box vectors
1008		scaled = invHmat * rs;
1009	+
1010		// wrap the vector back into the unit box by subtracting integer box
1011		// numbers
1012	<	for (int j = 0; j < 3; j++)
1012	>	for (int j = 0; j < 3; j++) {
1013		scaled[j] -= roundMe(scaled[j]);
1014	+	scaled[j] += 0.5;
1015	+	}
1016
1017		// find xyz-indices of cell that cutoffGroup is in.
1018		whichCell.x() = nCells_.x() * scaled.x();
#	Line 501 \| Line 1020 \| namespace OpenMD {
1020		whichCell.z() = nCells_.z() * scaled.z();
1021
1022		// find single index of this cell:
1023	<	cellIndex = Vlinear(whichCell, nCells_);
1023	>	cellIndex = Vlinear(whichCell, nCells_);
1024	>
1025		// add this cutoff group to the list of groups in this cell;
1026		cellList_[cellIndex].push_back(i);
1027		}
1028		#endif
1029
510	–
511	–
1030		for (int m1z = 0; m1z < nCells_.z(); m1z++) {
1031		for (int m1y = 0; m1y < nCells_.y(); m1y++) {
1032		for (int m1x = 0; m1x < nCells_.x(); m1x++) {
#	Line 553 \| Line 1071 \| namespace OpenMD {
1071		if (m2 != m1 \|\| cgColToGlobal[(j2)] < cgRowToGlobal[(j1)]) {
1072		dr = cgColData.position[(j2)] - cgRowData.position[(j1)];
1073		snap_->wrapVector(dr);
1074	<	if (dr.lengthSquare() < rl2) {
1074	>	cuts = getGroupCutoffs( (j1), (j2) );
1075	>	if (dr.lengthSquare() < cuts.third) {
1076		neighborList.push_back(make_pair((j1), (j2)));
1077		}
1078		}
1079		}
1080		}
1081		#else
1082	+
1083		for (vector<int>::iterator j1 = cellList_[m1].begin();
1084		j1 != cellList_[m1].end(); ++j1) {
1085		for (vector<int>::iterator j2 = cellList_[m2].begin();
1086		j2 != cellList_[m2].end(); ++j2) {
1087	<
1087	>
1088		// Always do this if we're in different cells or if
1089		// we're in the same cell and the global index of the
1090		// j2 cutoff group is less than the j1 cutoff group
#	Line 572 \| Line 1092 \| namespace OpenMD {
1092		if (m2 != m1 \|\| (j2) < (j1)) {
1093		dr = snap_->cgData.position[(j2)] - snap_->cgData.position[(j1)];
1094		snap_->wrapVector(dr);
1095	<	if (dr.lengthSquare() < rl2) {
1095	>	cuts = getGroupCutoffs( (j1), (j2) );
1096	>	if (dr.lengthSquare() < cuts.third) {
1097		neighborList.push_back(make_pair((j1), (j2)));
1098		}
1099		}
#	Line 583 \| Line 1104 \| namespace OpenMD {
1104		}
1105		}
1106		}
1107	<
1107	>
1108		// save the local cutoff group positions for the check that is
1109		// done on each loop:
1110		saved_CG_positions_.clear();
1111		for (int i = 0; i < nGroups_; i++)
1112		saved_CG_positions_.push_back(snap_->cgData.position[i]);
1113	<
1113	>
1114		return neighborList;
1115		}
1116		} //end namespace OpenMD

Diff Legend

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

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents): Revision 1568 by gezelter, Wed May 25 16:20:37 2011 UTC vs. Revision 1582 by gezelter, Tue Jun 14 20:41:44 2011 UTC

Diff Legend

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1568 by gezelter, Wed May 25 16:20:37 2011 UTC vs.
Revision 1582 by gezelter, Tue Jun 14 20:41:44 2011 UTC