@@ -95,15 +95,15 @@ namespace dftfe
9595 d_domainBoundingVectorsInverse.copyFrom (inv3 (d_domainBoundingVectors));
9696 if (d_isGroupSymmetry)
9797 {
98- const dftfe::Int max_size = 500 ;
98+ const dftfe::Int max_size = 2000 ;
9999 int rotation[max_size][3 ][3 ];
100100 double translation[max_size][3 ];
101101 double lattice[3 ][3 ];
102102 double position[d_numAtoms][3 ];
103103 int types[d_numAtoms];
104104 for (dftfe::uInt iVec = 0 ; iVec < 3 ; ++iVec)
105105 for (dftfe::uInt jDim = 0 ; jDim < 3 ; ++jDim)
106- lattice[iVec][jDim ] = d_domainBoundingVectors[3 * iVec + jDim];
106+ lattice[jDim][iVec ] = d_domainBoundingVectors[3 * iVec + jDim];
107107 for (dftfe::uInt iAtom = 0 ; iAtom < d_numAtoms; ++iAtom)
108108 {
109109 types[iAtom] = atomLocationsFractional[iAtom][0 ];
@@ -118,7 +118,7 @@ namespace dftfe
118118 position,
119119 types,
120120 d_numAtoms,
121- 1e-5 );
121+ 1e-8 );
122122 else
123123 {
124124 int equivalent_atoms[d_numAtoms];
@@ -136,25 +136,25 @@ namespace dftfe
136136 types,
137137 spins,
138138 d_numAtoms,
139- 1e-5 );
139+ 1e-8 );
140140 }
141141 d_symmMat.reserve (d_numSymm);
142142 d_symmMatInverse.reserve (d_numSymm);
143143 d_translation.reserve (d_numSymm);
144144 dftfe::uInt numSymm = 0 ;
145145 for (dftfe::uInt iSymm = 0 ; iSymm < d_numSymm; ++iSymm)
146- if (std::abs (translation[iSymm][0 ]) < 1e-8 &&
147- std::abs (translation[iSymm][1 ]) < 1e-8 &&
148- std::abs (translation[iSymm][2 ]) < 1e-8 )
149- {
150- d_symmMat.push_back (std::vector<double >(9 , 0.0 ));
151- d_translation.push_back (std::vector<double >(3 , 0.0 ));
152- for (dftfe::uInt jDim = 0 ; jDim < 3 ; ++jDim)
146+ {
147+ d_symmMat.push_back (std::vector<double >(9 , 0.0 ));
148+ d_translation.push_back (std::vector<double >(3 , 0.0 ));
149+ for (dftfe::uInt jDim = 0 ; jDim < 3 ; ++jDim)
150+ {
151+ d_translation.back ()[jDim] = translation[iSymm][jDim];
153152 for (dftfe::uInt kDim = 0 ; kDim < 3 ; ++kDim )
154153 d_symmMat.back ()[kDim * 3 + jDim] =
155154 static_cast <double >(rotation[iSymm][jDim][kDim ]);
156- d_symmMatInverse.push_back (inv3 (d_symmMat.back ()));
157- }
155+ }
156+ d_symmMatInverse.push_back (inv3 (d_symmMat.back ()));
157+ }
158158 d_symmMat.shrink_to_fit ();
159159 d_symmMatInverse.shrink_to_fit ();
160160 d_translation.shrink_to_fit ();
@@ -262,7 +262,7 @@ namespace dftfe
262262 for (int kDim = 0 ; kDim < 3 ; ++kDim )
263263 dot += m[iDim * 3 + kDim ] * m[jDim * 3 + kDim ];
264264 double orthoVal = iDim == jDim ? 1.0 : 0.0 ;
265- if (std::fabs (dot - orthoVal) > 1e-6 )
265+ if (std::fabs (dot - orthoVal) > 1e-8 )
266266 return false ;
267267 }
268268 }
@@ -295,7 +295,7 @@ namespace dftfe
295295 std::vector<dftfe::Int> kPointSymmetryMap (numKPoints, -1 );
296296 auto wrap = [](double x) {
297297 double r = std::remainder (x, 1.0 );
298- return (r <= - 0.5 ? 0.5 : r);
298+ return (r >= 0.5 ? r - 1.0 : r);
299299 };
300300 auto periodicDist = [](double a, double b) noexcept {
301301 double d = std::fabs (a - b);
@@ -384,7 +384,9 @@ namespace dftfe
384384 requiredPointCoordinates.clear ();
385385 requiredPointCoordinates.reserve (d_numSymm * nodalCoordinates.size ());
386386 localDoFIndexToPointIndexMap.clear ();
387- localDoFIndexToPointIndexMap.resize (d_numSymm);
387+ localDoFIndexToPointIndexMap.resize (
388+ d_numSymm, std::vector<dftfe::uInt>(dofHandler.n_locally_owned_dofs ()));
389+ std::map<std::array<std::int64_t , 3 >, dftfe::uInt> pointToPointIndexMap;
388390 for (dftfe::uInt iSymm = 0 ; iSymm < d_numSymm; ++iSymm)
389391 for (dealii::IndexSet::ElementIterator it = locallyOwnedNodes.begin ();
390392 it != locallyOwnedNodes.end ();
@@ -409,7 +411,7 @@ namespace dftfe
409411 for (dftfe::uInt jDim = 0 ; jDim < 3 ; ++jDim)
410412 transformedNodeCoordinatesFrac[iDim] +=
411413 d_symmMatInverse[iSymm][jDim * 3 + iDim] *
412- currentNodeCoordinatesFrac[jDim];
414+ ( currentNodeCoordinatesFrac[jDim] - d_translation[iSymm][jDim]) ;
413415 for (dftfe::uInt iDim = 0 ; iDim < 3 ; ++iDim)
414416 transformedNodeCoordinatesFrac[iDim] =
415417 transformedNodeCoordinatesFrac[iDim] -
@@ -421,10 +423,30 @@ namespace dftfe
421423 transformedNodeCoordinatesCart[iDim] +=
422424 d_domainBoundingVectors[3 * jDim + iDim] *
423425 transformedNodeCoordinatesFrac[jDim];
424- requiredPointCoordinates.push_back (transformedNodeCoordinatesCart);
425- localDoFIndexToPointIndexMap[iSymm][locallyOwnedNodes
426- .index_within_set (*it)] =
427- requiredPointCoordinates.size () - 1 ;
426+ std::array<std::int64_t , 3 > roundedCoords;
427+ roundedCoords[0 ] =
428+ std::llround (transformedNodeCoordinatesCart[0 ] * 1e8 );
429+ roundedCoords[1 ] =
430+ std::llround (transformedNodeCoordinatesCart[1 ] * 1e8 );
431+ roundedCoords[2 ] =
432+ std::llround (transformedNodeCoordinatesCart[2 ] * 1e8 );
433+ auto pointIterator = pointToPointIndexMap.find (roundedCoords);
434+ if (pointIterator != pointToPointIndexMap.end ())
435+ {
436+ localDoFIndexToPointIndexMap[iSymm][locallyOwnedNodes
437+ .index_within_set (*it)] =
438+ pointIterator->second ;
439+ }
440+ else
441+ {
442+ requiredPointCoordinates.push_back (
443+ transformedNodeCoordinatesCart);
444+ localDoFIndexToPointIndexMap[iSymm][locallyOwnedNodes
445+ .index_within_set (*it)] =
446+ requiredPointCoordinates.size () - 1 ;
447+ pointToPointIndexMap[roundedCoords] =
448+ requiredPointCoordinates.size () - 1 ;
449+ }
428450 }
429451 requiredPointCoordinates.shrink_to_fit ();
430452 remotePointCache.reinit (requiredPointCoordinates,
@@ -454,14 +476,15 @@ namespace dftfe
454476 for (dftfe::uInt iSymm = 0 ; iSymm < d_numSymm; ++iSymm)
455477 for (dftfe::uInt iPoint = 0 ; iPoint < numPoints; ++iPoint)
456478 {
457- std::vector<double > transformedPoint = d_translation[iSymm] ;
479+ std::vector<double > transformedPoint = { 0.0 , 0.0 , 0.0 } ;
458480 for (dftfe::uInt jDim = 0 ; jDim < 3 ; ++jDim)
459- transformedPoint[jDim] += d_symmMatInverse[iSymm][0 * 3 + jDim] *
460- globalPointCoords[3 * iPoint + 0 ] +
461- d_symmMatInverse[iSymm][1 * 3 + jDim] *
462- globalPointCoords[3 * iPoint + 1 ] +
463- d_symmMatInverse[iSymm][2 * 3 + jDim] *
464- globalPointCoords[3 * iPoint + 2 ];
481+ transformedPoint[jDim] =
482+ d_symmMatInverse[iSymm][0 * 3 + jDim] *
483+ (globalPointCoords[3 * iPoint + 0 ] - d_translation[iSymm][0 ]) +
484+ d_symmMatInverse[iSymm][1 * 3 + jDim] *
485+ (globalPointCoords[3 * iPoint + 1 ] - d_translation[iSymm][1 ]) +
486+ d_symmMatInverse[iSymm][2 * 3 + jDim] *
487+ (globalPointCoords[3 * iPoint + 2 ] - d_translation[iSymm][2 ]);
465488 for (dftfe::uInt jDim = 0 ; jDim < 3 ; ++jDim)
466489 transformedPoint[jDim] =
467490 transformedPoint[jDim] - std::floor (transformedPoint[jDim]);
@@ -500,8 +523,7 @@ namespace dftfe
500523 for (dftfe::uInt iDoF = 0 ; iDoF < scalarField.locally_owned_size ();
501524 ++iDoF)
502525 scalarField.local_element (iDoF) +=
503- pointValues[localDoFIndexToPointIndexMap[iSymm].find (iDoF)->second ] /
504- d_numSymm;
526+ pointValues[localDoFIndexToPointIndexMap[iSymm][iDoF]] / d_numSymm;
505527 }
506528
507529 void
0 commit comments