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

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1613 by gezelter, Thu Aug 18 20:18:19 2011 UTC vs.
Revision 1825 by gezelter, Wed Jan 9 19:27:52 2013 UTC

#	Line 36 \| Line 36
36		* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37		* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38		* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	<	* [4] Vardeman & Gezelter, in progress (2009).
39	>	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40	>	* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41		*/
42		#include "parallel/ForceMatrixDecomposition.hpp"
43		#include "math/SquareMatrix3.hpp"
#	Line 94 \| Line 95 \| namespace OpenMD {
95		storageLayout_ = sman_->getStorageLayout();
96		ff_ = info_->getForceField();
97		nLocal_ = snap_->getNumberOfAtoms();
98	<
98	>
99		nGroups_ = info_->getNLocalCutoffGroups();
100		// gather the information for atomtype IDs (atids):
101		idents = info_->getIdentArray();
#	Line 108 \| Line 109 \| namespace OpenMD {
109		PairList* oneTwo = info_->getOneTwoInteractions();
110		PairList* oneThree = info_->getOneThreeInteractions();
111		PairList* oneFour = info_->getOneFourInteractions();
112	<
112	>
113	>	if (needVelocities_)
114	>	snap_->cgData.setStorageLayout(DataStorage::dslPosition \|
115	>	DataStorage::dslVelocity);
116	>	else
117	>	snap_->cgData.setStorageLayout(DataStorage::dslPosition);
118	>
119		#ifdef IS_MPI
120
121		MPI::Intracomm row = rowComm.getComm();
#	Line 144 \| Line 151 \| namespace OpenMD {
151		cgRowData.resize(nGroupsInRow_);
152		cgRowData.setStorageLayout(DataStorage::dslPosition);
153		cgColData.resize(nGroupsInCol_);
154	<	cgColData.setStorageLayout(DataStorage::dslPosition);
155	<
154	>	if (needVelocities_)
155	>	// we only need column velocities if we need them.
156	>	cgColData.setStorageLayout(DataStorage::dslPosition \|
157	>	DataStorage::dslVelocity);
158	>	else
159	>	cgColData.setStorageLayout(DataStorage::dslPosition);
160	>
161		identsRow.resize(nAtomsInRow_);
162		identsCol.resize(nAtomsInCol_);
163
#	Line 164 \| Line 176 \| namespace OpenMD {
176		pot_row.resize(nAtomsInRow_);
177		pot_col.resize(nAtomsInCol_);
178
179	+	expot_row.resize(nAtomsInRow_);
180	+	expot_col.resize(nAtomsInCol_);
181	+
182		AtomRowToGlobal.resize(nAtomsInRow_);
183		AtomColToGlobal.resize(nAtomsInCol_);
184		AtomPlanIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal);
#	Line 247 \| Line 262 \| namespace OpenMD {
262		for (int j = 0; j < nLocal_; j++) {
263		int jglob = AtomLocalToGlobal[j];
264
265	<	if (excludes->hasPair(iglob, jglob))
265	>	if (excludes->hasPair(iglob, jglob))
266		excludesForAtom[i].push_back(j);
267
253	–
268		if (oneTwo->hasPair(iglob, jglob)) {
269		toposForAtom[i].push_back(j);
270		topoDist[i].push_back(1);
#	Line 296 \| Line 310 \| namespace OpenMD {
310
311		RealType tol = 1e-6;
312		largestRcut_ = 0.0;
299	–	RealType rc;
313		int atid;
314		set<AtomType*> atypes = info_->getSimulatedAtomTypes();
315
#	Line 381 \| Line 394 \| namespace OpenMD {
394		}
395
396		bool gTypeFound = false;
397	<	for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
397	>	for (unsigned int gt = 0; gt < gTypeCutoffs.size(); gt++) {
398		if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) {
399		groupToGtype[cg1] = gt;
400		gTypeFound = true;
#	Line 406 \| Line 419 \| namespace OpenMD {
419
420		RealType tradRcut = groupMax;
421
422	<	for (int i = 0; i < gTypeCutoffs.size(); i++) {
423	<	for (int j = 0; j < gTypeCutoffs.size(); j++) {
422	>	for (unsigned int i = 0; i < gTypeCutoffs.size(); i++) {
423	>	for (unsigned int j = 0; j < gTypeCutoffs.size(); j++) {
424		RealType thisRcut;
425		switch(cutoffPolicy_) {
426		case TRADITIONAL:
#	Line 450 \| Line 463 \| namespace OpenMD {
463		}
464		}
465
453	–
466		groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) {
467		int i, j;
468		#ifdef IS_MPI
#	Line 464 \| Line 476 \| namespace OpenMD {
476		}
477
478		int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) {
479	<	for (int j = 0; j < toposForAtom[atom1].size(); j++) {
479	>	for (unsigned int j = 0; j < toposForAtom[atom1].size(); j++) {
480		if (toposForAtom[atom1][j] == atom2)
481		return topoDist[atom1][j];
482		}
#	Line 474 \| Line 486 \| namespace OpenMD {
486		void ForceMatrixDecomposition::zeroWorkArrays() {
487		pairwisePot = 0.0;
488		embeddingPot = 0.0;
489	+	excludedPot = 0.0;
490	+	excludedSelfPot = 0.0;
491
492		#ifdef IS_MPI
493		if (storageLayout_ & DataStorage::dslForce) {
#	Line 492 \| Line 506 \| namespace OpenMD {
506		fill(pot_col.begin(), pot_col.end(),
507		Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
508
509	+	fill(expot_row.begin(), expot_row.end(),
510	+	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
511	+
512	+	fill(expot_col.begin(), expot_col.end(),
513	+	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
514	+
515		if (storageLayout_ & DataStorage::dslParticlePot) {
516		fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(),
517		0.0);
#	Line 525 \| Line 545 \| namespace OpenMD {
545		atomColData.skippedCharge.end(), 0.0);
546		}
547
548	+	if (storageLayout_ & DataStorage::dslFlucQForce) {
549	+	fill(atomRowData.flucQFrc.begin(),
550	+	atomRowData.flucQFrc.end(), 0.0);
551	+	fill(atomColData.flucQFrc.begin(),
552	+	atomColData.flucQFrc.end(), 0.0);
553	+	}
554	+
555	+	if (storageLayout_ & DataStorage::dslElectricField) {
556	+	fill(atomRowData.electricField.begin(),
557	+	atomRowData.electricField.end(), V3Zero);
558	+	fill(atomColData.electricField.begin(),
559	+	atomColData.electricField.end(), V3Zero);
560	+	}
561	+
562		#endif
563		// even in parallel, we need to zero out the local arrays:
564
#	Line 537 \| Line 571 \| namespace OpenMD {
571		fill(snap_->atomData.density.begin(),
572		snap_->atomData.density.end(), 0.0);
573		}
574	+
575		if (storageLayout_ & DataStorage::dslFunctional) {
576		fill(snap_->atomData.functional.begin(),
577		snap_->atomData.functional.end(), 0.0);
578		}
579	+
580		if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
581		fill(snap_->atomData.functionalDerivative.begin(),
582		snap_->atomData.functionalDerivative.end(), 0.0);
583		}
584	+
585		if (storageLayout_ & DataStorage::dslSkippedCharge) {
586		fill(snap_->atomData.skippedCharge.begin(),
587		snap_->atomData.skippedCharge.end(), 0.0);
588		}
589	<
589	>
590	>	if (storageLayout_ & DataStorage::dslElectricField) {
591	>	fill(snap_->atomData.electricField.begin(),
592	>	snap_->atomData.electricField.end(), V3Zero);
593	>	}
594		}
595
596
#	Line 572 \| Line 613 \| namespace OpenMD {
613		cgPlanVectorColumn->gather(snap_->cgData.position,
614		cgColData.position);
615
616	+
617	+
618	+	if (needVelocities_) {
619	+	// gather up the atomic velocities
620	+	AtomPlanVectorColumn->gather(snap_->atomData.velocity,
621	+	atomColData.velocity);
622	+
623	+	cgPlanVectorColumn->gather(snap_->cgData.velocity,
624	+	cgColData.velocity);
625	+	}
626	+
627
628		// if needed, gather the atomic rotation matrices
629		if (storageLayout_ & DataStorage::dslAmat) {
#	Line 580 \| Line 632 \| namespace OpenMD {
632		AtomPlanMatrixColumn->gather(snap_->atomData.aMat,
633		atomColData.aMat);
634		}
635	<
636	<	// if needed, gather the atomic eletrostatic frames
637	<	if (storageLayout_ & DataStorage::dslElectroFrame) {
638	<	AtomPlanMatrixRow->gather(snap_->atomData.electroFrame,
639	<	atomRowData.electroFrame);
640	<	AtomPlanMatrixColumn->gather(snap_->atomData.electroFrame,
641	<	atomColData.electroFrame);
635	>
636	>	// if needed, gather the atomic eletrostatic information
637	>	if (storageLayout_ & DataStorage::dslDipole) {
638	>	AtomPlanVectorRow->gather(snap_->atomData.dipole,
639	>	atomRowData.dipole);
640	>	AtomPlanVectorColumn->gather(snap_->atomData.dipole,
641	>	atomColData.dipole);
642	>	}
643	>
644	>	if (storageLayout_ & DataStorage::dslQuadrupole) {
645	>	AtomPlanMatrixRow->gather(snap_->atomData.quadrupole,
646	>	atomRowData.quadrupole);
647	>	AtomPlanMatrixColumn->gather(snap_->atomData.quadrupole,
648	>	atomColData.quadrupole);
649	>	}
650	>
651	>	// if needed, gather the atomic fluctuating charge values
652	>	if (storageLayout_ & DataStorage::dslFlucQPosition) {
653	>	AtomPlanRealRow->gather(snap_->atomData.flucQPos,
654	>	atomRowData.flucQPos);
655	>	AtomPlanRealColumn->gather(snap_->atomData.flucQPos,
656	>	atomColData.flucQPos);
657		}
658
659		#endif
#	Line 611 \| Line 678 \| namespace OpenMD {
678		for (int i = 0; i < n; i++)
679		snap_->atomData.density[i] += rho_tmp[i];
680		}
681	+
682	+	if (storageLayout_ & DataStorage::dslElectricField) {
683	+
684	+	AtomPlanVectorRow->scatter(atomRowData.electricField,
685	+	snap_->atomData.electricField);
686	+
687	+	int n = snap_->atomData.electricField.size();
688	+	vector<Vector3d> field_tmp(n, V3Zero);
689	+	AtomPlanVectorColumn->scatter(atomColData.electricField,
690	+	field_tmp);
691	+	for (int i = 0; i < n; i++)
692	+	snap_->atomData.electricField[i] += field_tmp[i];
693	+	}
694		#endif
695		}
696
#	Line 690 \| Line 770 \| namespace OpenMD {
770
771		}
772
773	+	if (storageLayout_ & DataStorage::dslFlucQForce) {
774	+
775	+	int nq = snap_->atomData.flucQFrc.size();
776	+	vector<RealType> fqfrc_tmp(nq, 0.0);
777	+
778	+	AtomPlanRealRow->scatter(atomRowData.flucQFrc, fqfrc_tmp);
779	+	for (int i = 0; i < nq; i++) {
780	+	snap_->atomData.flucQFrc[i] += fqfrc_tmp[i];
781	+	fqfrc_tmp[i] = 0.0;
782	+	}
783	+
784	+	AtomPlanRealColumn->scatter(atomColData.flucQFrc, fqfrc_tmp);
785	+	for (int i = 0; i < nq; i++)
786	+	snap_->atomData.flucQFrc[i] += fqfrc_tmp[i];
787	+
788	+	}
789	+
790		nLocal_ = snap_->getNumberOfAtoms();
791
792		vector<potVec> pot_temp(nLocal_,
793		Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
794	+	vector<potVec> expot_temp(nLocal_,
795	+	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
796
797		// scatter/gather pot_row into the members of my column
798
799		AtomPlanPotRow->scatter(pot_row, pot_temp);
800	+	AtomPlanPotRow->scatter(expot_row, expot_temp);
801
802	<	for (int ii = 0; ii < pot_temp.size(); ii++ )
802	>	for (int ii = 0; ii < pot_temp.size(); ii++ )
803		pairwisePot += pot_temp[ii];
804	<
804	>
805	>	for (int ii = 0; ii < expot_temp.size(); ii++ )
806	>	excludedPot += expot_temp[ii];
807	>
808	>	if (storageLayout_ & DataStorage::dslParticlePot) {
809	>	// This is the pairwise contribution to the particle pot. The
810	>	// embedding contribution is added in each of the low level
811	>	// non-bonded routines. In single processor, this is done in
812	>	// unpackInteractionData, not in collectData.
813	>	for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
814	>	for (int i = 0; i < nLocal_; i++) {
815	>	// factor of two is because the total potential terms are divided
816	>	// by 2 in parallel due to row/ column scatter
817	>	snap_->atomData.particlePot[i] += 2.0 * pot_temp[i](ii);
818	>	}
819	>	}
820	>	}
821	>
822		fill(pot_temp.begin(), pot_temp.end(),
823		Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
824	+	fill(expot_temp.begin(), expot_temp.end(),
825	+	Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
826
827		AtomPlanPotColumn->scatter(pot_col, pot_temp);
828	+	AtomPlanPotColumn->scatter(expot_col, expot_temp);
829
830		for (int ii = 0; ii < pot_temp.size(); ii++ )
831		pairwisePot += pot_temp[ii];
832	+
833	+	for (int ii = 0; ii < expot_temp.size(); ii++ )
834	+	excludedPot += expot_temp[ii];
835	+
836	+	if (storageLayout_ & DataStorage::dslParticlePot) {
837	+	// This is the pairwise contribution to the particle pot. The
838	+	// embedding contribution is added in each of the low level
839	+	// non-bonded routines. In single processor, this is done in
840	+	// unpackInteractionData, not in collectData.
841	+	for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
842	+	for (int i = 0; i < nLocal_; i++) {
843	+	// factor of two is because the total potential terms are divided
844	+	// by 2 in parallel due to row/ column scatter
845	+	snap_->atomData.particlePot[i] += 2.0 * pot_temp[i](ii);
846	+	}
847	+	}
848	+	}
849
850	+	if (storageLayout_ & DataStorage::dslParticlePot) {
851	+	int npp = snap_->atomData.particlePot.size();
852	+	vector<RealType> ppot_temp(npp, 0.0);
853	+
854	+	// This is the direct or embedding contribution to the particle
855	+	// pot.
856	+
857	+	AtomPlanRealRow->scatter(atomRowData.particlePot, ppot_temp);
858	+	for (int i = 0; i < npp; i++) {
859	+	snap_->atomData.particlePot[i] += ppot_temp[i];
860	+	}
861	+
862	+	fill(ppot_temp.begin(), ppot_temp.end(), 0.0);
863	+
864	+	AtomPlanRealColumn->scatter(atomColData.particlePot, ppot_temp);
865	+	for (int i = 0; i < npp; i++) {
866	+	snap_->atomData.particlePot[i] += ppot_temp[i];
867	+	}
868	+	}
869	+
870		for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
871		RealType ploc1 = pairwisePot[ii];
872		RealType ploc2 = 0.0;
#	Line 718 \| Line 875 \| namespace OpenMD {
875		}
876
877		for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
878	<	RealType ploc1 = embeddingPot[ii];
878	>	RealType ploc1 = excludedPot[ii];
879		RealType ploc2 = 0.0;
880		MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
881	<	embeddingPot[ii] = ploc2;
881	>	excludedPot[ii] = ploc2;
882		}
883	+
884	+	// Here be dragons.
885	+	MPI::Intracomm col = colComm.getComm();
886	+
887	+	col.Allreduce(MPI::IN_PLACE,
888	+	&snap_->frameData.conductiveHeatFlux[0], 3,
889	+	MPI::REALTYPE, MPI::SUM);
890
891	+
892		#endif
893
894		}
895
896	+	/**
897	+	* Collects information obtained during the post-pair (and embedding
898	+	* functional) loops onto local data structures.
899	+	*/
900	+	void ForceMatrixDecomposition::collectSelfData() {
901	+	snap_ = sman_->getCurrentSnapshot();
902	+	storageLayout_ = sman_->getStorageLayout();
903	+
904	+	#ifdef IS_MPI
905	+	for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
906	+	RealType ploc1 = embeddingPot[ii];
907	+	RealType ploc2 = 0.0;
908	+	MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
909	+	embeddingPot[ii] = ploc2;
910	+	}
911	+	for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
912	+	RealType ploc1 = excludedSelfPot[ii];
913	+	RealType ploc2 = 0.0;
914	+	MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
915	+	excludedSelfPot[ii] = ploc2;
916	+	}
917	+	#endif
918	+
919	+	}
920	+
921	+
922	+
923		int ForceMatrixDecomposition::getNAtomsInRow() {
924		#ifdef IS_MPI
925		return nAtomsInRow_;
#	Line 768 \| Line 960 \| namespace OpenMD {
960		return d;
961		}
962
963	+	Vector3d ForceMatrixDecomposition::getGroupVelocityColumn(int cg2){
964	+	#ifdef IS_MPI
965	+	return cgColData.velocity[cg2];
966	+	#else
967	+	return snap_->cgData.velocity[cg2];
968	+	#endif
969	+	}
970
971	+	Vector3d ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){
972	+	#ifdef IS_MPI
973	+	return atomColData.velocity[atom2];
974	+	#else
975	+	return snap_->atomData.velocity[atom2];
976	+	#endif
977	+	}
978	+
979	+
980		Vector3d ForceMatrixDecomposition::getAtomToGroupVectorRow(int atom1, int cg1){
981
982		Vector3d d;
#	Line 834 \| Line 1042 \| namespace OpenMD {
1042		* We need to exclude some overcounted interactions that result from
1043		* the parallel decomposition.
1044		*/
1045	<	bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) {
1045	>	bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2, int cg1, int cg2) {
1046		int unique_id_1, unique_id_2;
1047	<
1047	>
1048		#ifdef IS_MPI
1049		// in MPI, we have to look up the unique IDs for each atom
1050		unique_id_1 = AtomRowToGlobal[atom1];
1051		unique_id_2 = AtomColToGlobal[atom2];
1052	<
1053	<	// this situation should only arise in MPI simulations
1052	>	// group1 = cgRowToGlobal[cg1];
1053	>	// group2 = cgColToGlobal[cg2];
1054	>	#else
1055	>	unique_id_1 = AtomLocalToGlobal[atom1];
1056	>	unique_id_2 = AtomLocalToGlobal[atom2];
1057	>	int group1 = cgLocalToGlobal[cg1];
1058	>	int group2 = cgLocalToGlobal[cg2];
1059	>	#endif
1060	>
1061		if (unique_id_1 == unique_id_2) return true;
1062	<
1062	>
1063	>	#ifdef IS_MPI
1064		// this prevents us from doing the pair on multiple processors
1065		if (unique_id_1 < unique_id_2) {
1066		if ((unique_id_1 + unique_id_2) % 2 == 0) return true;
1067		} else {
1068	<	if ((unique_id_1 + unique_id_2) % 2 == 1) return true;
1068	>	if ((unique_id_1 + unique_id_2) % 2 == 1) return true;
1069		}
1070	+	#endif
1071	+
1072	+	#ifndef IS_MPI
1073	+	if (group1 == group2) {
1074	+	if (unique_id_1 < unique_id_2) return true;
1075	+	}
1076		#endif
1077	+
1078		return false;
1079		}
1080
#	Line 871 \| Line 1094 \| namespace OpenMD {
1094
1095		for (vector<int>::iterator i = excludesForAtom[atom1].begin();
1096		i != excludesForAtom[atom1].end(); ++i) {
1097	<	if ( (*i) == atom2 ) return true;
1097	>	if ( (*i) == atom2 ) return true;
1098		}
1099
1100		return false;
#	Line 910 \| Line 1133 \| namespace OpenMD {
1133		idat.A2 = &(atomColData.aMat[atom2]);
1134		}
1135
913	–	if (storageLayout_ & DataStorage::dslElectroFrame) {
914	–	idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
915	–	idat.eFrame2 = &(atomColData.electroFrame[atom2]);
916	–	}
917	–
1136		if (storageLayout_ & DataStorage::dslTorque) {
1137		idat.t1 = &(atomRowData.torque[atom1]);
1138		idat.t2 = &(atomColData.torque[atom2]);
1139		}
1140
1141	+	if (storageLayout_ & DataStorage::dslDipole) {
1142	+	idat.dipole1 = &(atomRowData.dipole[atom1]);
1143	+	idat.dipole2 = &(atomColData.dipole[atom2]);
1144	+	}
1145	+
1146	+	if (storageLayout_ & DataStorage::dslQuadrupole) {
1147	+	idat.quadrupole1 = &(atomRowData.quadrupole[atom1]);
1148	+	idat.quadrupole2 = &(atomColData.quadrupole[atom2]);
1149	+	}
1150	+
1151		if (storageLayout_ & DataStorage::dslDensity) {
1152		idat.rho1 = &(atomRowData.density[atom1]);
1153		idat.rho2 = &(atomColData.density[atom2]);
#	Line 945 \| Line 1173 \| namespace OpenMD {
1173		idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]);
1174		}
1175
1176	<	#else
1176	>	if (storageLayout_ & DataStorage::dslFlucQPosition) {
1177	>	idat.flucQ1 = &(atomRowData.flucQPos[atom1]);
1178	>	idat.flucQ2 = &(atomColData.flucQPos[atom2]);
1179	>	}
1180
1181	+	#else
1182	+
1183		idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]);
951	–	//idat.atypes = make_pair( ff_->getAtomType(idents[atom1]),
952	–	// ff_->getAtomType(idents[atom2]) );
1184
1185		if (storageLayout_ & DataStorage::dslAmat) {
1186		idat.A1 = &(snap_->atomData.aMat[atom1]);
1187		idat.A2 = &(snap_->atomData.aMat[atom2]);
1188		}
1189
1190	<	if (storageLayout_ & DataStorage::dslElectroFrame) {
960	<	idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
961	<	idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
962	<	}
1190	>	RealType ct = dot(idat.A1->getColumn(2), idat.A2->getColumn(2));
1191
1192		if (storageLayout_ & DataStorage::dslTorque) {
1193		idat.t1 = &(snap_->atomData.torque[atom1]);
1194		idat.t2 = &(snap_->atomData.torque[atom2]);
1195		}
1196
1197	+	if (storageLayout_ & DataStorage::dslDipole) {
1198	+	idat.dipole1 = &(snap_->atomData.dipole[atom1]);
1199	+	idat.dipole2 = &(snap_->atomData.dipole[atom2]);
1200	+	}
1201	+
1202	+	if (storageLayout_ & DataStorage::dslQuadrupole) {
1203	+	idat.quadrupole1 = &(snap_->atomData.quadrupole[atom1]);
1204	+	idat.quadrupole2 = &(snap_->atomData.quadrupole[atom2]);
1205	+	}
1206	+
1207		if (storageLayout_ & DataStorage::dslDensity) {
1208		idat.rho1 = &(snap_->atomData.density[atom1]);
1209		idat.rho2 = &(snap_->atomData.density[atom2]);
#	Line 990 \| Line 1228 \| namespace OpenMD {
1228		idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]);
1229		idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]);
1230		}
1231	+
1232	+	if (storageLayout_ & DataStorage::dslFlucQPosition) {
1233	+	idat.flucQ1 = &(snap_->atomData.flucQPos[atom1]);
1234	+	idat.flucQ2 = &(snap_->atomData.flucQPos[atom2]);
1235	+	}
1236	+
1237		#endif
1238		}
1239
1240
1241		void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat, int atom1, int atom2) {
1242		#ifdef IS_MPI
1243	<	pot_row[atom1] += 0.5 * *(idat.pot);
1244	<	pot_col[atom2] += 0.5 * *(idat.pot);
1243	>	pot_row[atom1] += RealType(0.5) * *(idat.pot);
1244	>	pot_col[atom2] += RealType(0.5) * *(idat.pot);
1245	>	expot_row[atom1] += RealType(0.5) * *(idat.excludedPot);
1246	>	expot_col[atom2] += RealType(0.5) * *(idat.excludedPot);
1247
1248		atomRowData.force[atom1] += *(idat.f1);
1249		atomColData.force[atom2] -= *(idat.f1);
1250	+
1251	+	if (storageLayout_ & DataStorage::dslFlucQForce) {
1252	+	atomRowData.flucQFrc[atom1] -= *(idat.dVdFQ1);
1253	+	atomColData.flucQFrc[atom2] -= *(idat.dVdFQ2);
1254	+	}
1255	+
1256	+	if (storageLayout_ & DataStorage::dslElectricField) {
1257	+	atomRowData.electricField[atom1] += *(idat.eField1);
1258	+	atomColData.electricField[atom2] += *(idat.eField2);
1259	+	}
1260	+
1261		#else
1262		pairwisePot += *(idat.pot);
1263	+	excludedPot += *(idat.excludedPot);
1264
1265		snap_->atomData.force[atom1] += *(idat.f1);
1266		snap_->atomData.force[atom2] -= *(idat.f1);
1267	+
1268	+	if (idat.doParticlePot) {
1269	+	// This is the pairwise contribution to the particle pot. The
1270	+	// embedding contribution is added in each of the low level
1271	+	// non-bonded routines. In parallel, this calculation is done
1272	+	// in collectData, not in unpackInteractionData.
1273	+	snap_->atomData.particlePot[atom1] += (idat.vpair) *(idat.sw);
1274	+	snap_->atomData.particlePot[atom2] += (idat.vpair) *(idat.sw);
1275	+	}
1276	+
1277	+	if (storageLayout_ & DataStorage::dslFlucQForce) {
1278	+	snap_->atomData.flucQFrc[atom1] -= *(idat.dVdFQ1);
1279	+	snap_->atomData.flucQFrc[atom2] -= *(idat.dVdFQ2);
1280	+	}
1281	+
1282	+	if (storageLayout_ & DataStorage::dslElectricField) {
1283	+	snap_->atomData.electricField[atom1] += *(idat.eField1);
1284	+	snap_->atomData.electricField[atom2] += *(idat.eField2);
1285	+	}
1286	+
1287		#endif
1288
1289		}
#	Line 1030 \| Line 1308 \| namespace OpenMD {
1308		#endif
1309
1310		RealType rList_ = (largestRcut_ + skinThickness_);
1033	–	RealType rl2 = rList_ * rList_;
1311		Snapshot* snap_ = sman_->getCurrentSnapshot();
1312		Mat3x3d Hmat = snap_->getHmat();
1313		Vector3d Hx = Hmat.getColumn(0);
#	Line 1074 \| Line 1351 \| namespace OpenMD {
1351		for (int j = 0; j < 3; j++) {
1352		scaled[j] -= roundMe(scaled[j]);
1353		scaled[j] += 0.5;
1354	+	// Handle the special case when an object is exactly on the
1355	+	// boundary (a scaled coordinate of 1.0 is the same as
1356	+	// scaled coordinate of 0.0)
1357	+	if (scaled[j] >= 1.0) scaled[j] -= 1.0;
1358		}
1359
1360		// find xyz-indices of cell that cutoffGroup is in.
#	Line 1098 \| Line 1379 \| namespace OpenMD {
1379		for (int j = 0; j < 3; j++) {
1380		scaled[j] -= roundMe(scaled[j]);
1381		scaled[j] += 0.5;
1382	+	// Handle the special case when an object is exactly on the
1383	+	// boundary (a scaled coordinate of 1.0 is the same as
1384	+	// scaled coordinate of 0.0)
1385	+	if (scaled[j] >= 1.0) scaled[j] -= 1.0;
1386		}
1387
1388		// find xyz-indices of cell that cutoffGroup is in.
#	Line 1124 \| Line 1409 \| namespace OpenMD {
1409		for (int j = 0; j < 3; j++) {
1410		scaled[j] -= roundMe(scaled[j]);
1411		scaled[j] += 0.5;
1412	+	// Handle the special case when an object is exactly on the
1413	+	// boundary (a scaled coordinate of 1.0 is the same as
1414	+	// scaled coordinate of 0.0)
1415	+	if (scaled[j] >= 1.0) scaled[j] -= 1.0;
1416		}
1417
1418		// find xyz-indices of cell that cutoffGroup is in.
#	Line 1190 \| Line 1479 \| namespace OpenMD {
1479		}
1480		}
1481		#else
1193	–
1482		for (vector<int>::iterator j1 = cellList_[m1].begin();
1483		j1 != cellList_[m1].end(); ++j1) {
1484		for (vector<int>::iterator j2 = cellList_[m2].begin();
1485		j2 != cellList_[m2].end(); ++j2) {
1486	<
1486	>
1487		// Always do this if we're in different cells or if
1488	<	// we're in the same cell and the global index of the
1489	<	// j2 cutoff group is less than the j1 cutoff group
1490	<
1491	<	if (m2 != m1 \|\| (j2) < (j1)) {
1488	>	// we're in the same cell and the global index of
1489	>	// the j2 cutoff group is greater than or equal to
1490	>	// the j1 cutoff group. Note that Rappaport's code
1491	>	// has a "less than" conditional here, but that
1492	>	// deals with atom-by-atom computation. OpenMD
1493	>	// allows atoms within a single cutoff group to
1494	>	// interact with each other.
1495	>
1496	>
1497	>
1498	>	if (m2 != m1 \|\| (j2) >= (j1) ) {
1499	>
1500		dr = snap_->cgData.position[(j2)] - snap_->cgData.position[(j1)];
1501		snap_->wrapVector(dr);
1502		cuts = getGroupCutoffs( (j1), (j2) );
#	Line 1219 \| Line 1515 \| namespace OpenMD {
1515		// branch to do all cutoff group pairs
1516		#ifdef IS_MPI
1517		for (int j1 = 0; j1 < nGroupsInRow_; j1++) {
1518	<	for (int j2 = 0; j2 < nGroupsInCol_; j2++) {
1518	>	for (int j2 = 0; j2 < nGroupsInCol_; j2++) {
1519		dr = cgColData.position[j2] - cgRowData.position[j1];
1520		snap_->wrapVector(dr);
1521		cuts = getGroupCutoffs( j1, j2 );
#	Line 1227 \| Line 1523 \| namespace OpenMD {
1523		neighborList.push_back(make_pair(j1, j2));
1524		}
1525		}
1526	<	}
1526	>	}
1527		#else
1528	<	for (int j1 = 0; j1 < nGroups_ - 1; j1++) {
1529	<	for (int j2 = j1 + 1; j2 < nGroups_; j2++) {
1528	>	// include all groups here.
1529	>	for (int j1 = 0; j1 < nGroups_; j1++) {
1530	>	// include self group interactions j2 == j1
1531	>	for (int j2 = j1; j2 < nGroups_; j2++) {
1532		dr = snap_->cgData.position[j2] - snap_->cgData.position[j1];
1533		snap_->wrapVector(dr);
1534		cuts = getGroupCutoffs( j1, j2 );
1535		if (dr.lengthSquare() < cuts.third) {
1536		neighborList.push_back(make_pair(j1, j2));
1537		}
1538	<	}
1539	<	}
1538	>	}
1539	>	}
1540		#endif
1541		}
1542

Diff Legend

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

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents): Revision 1613 by gezelter, Thu Aug 18 20:18:19 2011 UTC vs. Revision 1825 by gezelter, Wed Jan 9 19:27:52 2013 UTC

Diff Legend

Comparing branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1613 by gezelter, Thu Aug 18 20:18:19 2011 UTC vs.
Revision 1825 by gezelter, Wed Jan 9 19:27:52 2013 UTC