Fix bug (deltaspin): restore LCAO spin rotation invariance

source/source_hsolver/module_genelpa/elpa_new_complex.cpp

@@ -79,7 +79,10 @@ int ELPA_Solver::generalized_eigenvector(std::complex<double>* A, std::complex<d
             t=-1;
             timer(myid, "A*U^-1", "2.1a", t);
         }
-        ScalapackConnector::gemm('C', 'N', nFull, nFull, nFull, 1.0, A, B, 0.0, zwork.data(), desc);
+        // LCAO provides the authoritative upper triangle of the Hermitian
+        // Hamiltonian. PZHEMM observes that contract, while PZGEMM would read
+        // an uninitialized/stale lower triangle through A^H.
+        ScalapackConnector::hemm('L', 'U', nFull, 1.0, A, B, 0.0, zwork.data(), desc);
         if(loglevel>1)
         {
             timer(myid, "A*U^-1", "2.1a", t);
@@ -105,6 +108,9 @@ int ELPA_Solver::generalized_eigenvector(std::complex<double>* A, std::complex<d
             t=-1;
             timer(myid, "B*A^T", "2.1b", t);
         }
+        // A now stores the general product H * U^-1, so its conjugate
+        // transpose must be formed with GEMM rather than treating A as
+        // Hermitian a second time.
         ScalapackConnector::gemm('N', 'C', nFull, nFull, nFull, 1.0, B, A, 0.0, zwork.data(), desc);
         if(loglevel>1)
         {

source/source_hsolver/test/diago_lcao_test.cpp

@@ -371,6 +371,33 @@ INSTANTIATE_TEST_SUITE_P(
         DiagoPrepare<std::complex<double>>(0, 0, 1, 0, "scalapack_gvx", "H-KPoints-Si2.dat", "S-KPoints-Si2.dat"),
         DiagoPrepare<std::complex<double>>(0, 0, 32, 0, "scalapack_gvx", "H-KPoints-Si64.dat", "S-KPoints-Si64.dat")));
 
+#ifdef __ELPA
+TEST(DiagoElpaComplexTest, UsesAuthoritativeUpperTriangle)
+{
+    std::stringstream out_info;
+    DiagoPrepare<std::complex<double>> dp(0, 0, 1, 0, "genelpa", "H-KPoints-Si2.dat", "S-KPoints-Si2.dat");
+    ASSERT_TRUE(dp.produce_HS());
+
+    if (dp.myrank == 0)
+    {
+        dp.diago_lapack();
+        for (int row = 1; row < dp.nlocal; ++row)
+        {
+            for (int col = 0; col < row; ++col)
+            {
+                dp.h[row * dp.nlocal + col] = std::complex<double>(17.0 + row + col, -13.0);
+            }
+        }
+    }
+
+    dp.diago();
+    if (dp.myrank == 0)
+    {
+        EXPECT_TRUE(dp.compare_eigen(out_info)) << out_info.str();
+    }
+}
+#endif
+
 int main(int argc, char** argv)
 {
     MPI_Init(&argc, &argv);

source/source_lcao/module_deltaspin/mi_tools.h

@@ -1,11 +1,12 @@
 #ifndef MI_TOOLS_H
 #define MI_TOOLS_H
 
+#include "source_base/vector3.h"
+
+#include <array>
 #include <complex>
 #include <vector>
 
-#include "source_base/vector3.h"
-
 /**
  * @file mi_tools.h
  * @brief Free-function utilities for computing atomic magnetic moments (Mi)
@@ -23,6 +24,21 @@
 namespace spinconstrain
 {
 
+/**
+ * @brief Convert a Cartesian Pauli vector to a spinor-space operator.
+ *
+ * For sigma_y = [[0, -i], [i, 0]], lambda dot sigma is stored in row-major
+ * order as {lambda_z, lambda_x - i lambda_y,
+ *           lambda_x + i lambda_y, -lambda_z}.
+ */
+inline std::array<std::complex<double>, 4> pauli_vector_to_spinor(const ModuleBase::Vector3<double>& lambda)
+{
+    return {{std::complex<double>(lambda.z, 0.0),
+             std::complex<double>(lambda.x, -lambda.y),
+             std::complex<double>(lambda.x, lambda.y),
+             std::complex<double>(-lambda.z, 0.0)}};
+}
+
 /**
  * @brief Convert spinor occupation matrix to magnetic moment vector using Pauli matrices.
  *

source/source_lcao/module_deltaspin/test/deltaspin_core_test.cpp

@@ -1,8 +1,10 @@
+#include "source_lcao/module_deltaspin/mi_tools.h"
+
 #include "gtest/gtest.h"
+#include <algorithm>
 #include <cmath>
 #include <complex>
 #include <vector>
-#include <algorithm>
 
 /***********************************************************************
  * Unit tests for DeltaSpin core algorithms.
@@ -92,6 +94,28 @@ TEST_F(PauliToMomentTest, GeneralCase_AllComponents)
     EXPECT_NEAR(M.z, 0.2, 1e-15);
 }
 
+TEST(PauliConventionTest, LambdaExpectationMatchesDotMoment)
+{
+    const double amplitude = 1.0 / std::sqrt(2.0);
+    const std::complex<double> spinor[2] = {{amplitude, 0.0}, {0.0, amplitude}};
+    const std::complex<double> occ[4] = {std::conj(spinor[0]) * spinor[0],
+                                         std::conj(spinor[0]) * spinor[1],
+                                         std::conj(spinor[1]) * spinor[0],
+                                         std::conj(spinor[1]) * spinor[1]};
+    const ModuleBase::Vector3<double> lambda(0.0, 2.0, 0.0);
+    const auto matrix = spinconstrain::pauli_vector_to_spinor(lambda);
+    const auto moment = spinconstrain::pauli_to_moment(occ, 1.0);
+
+    const std::complex<double> h_up = matrix[0] * spinor[0] + matrix[1] * spinor[1];
+    const std::complex<double> h_down = matrix[2] * spinor[0] + matrix[3] * spinor[1];
+    const double expectation = (std::conj(spinor[0]) * h_up + std::conj(spinor[1]) * h_down).real();
+    const double dot_moment = lambda.x * moment.x + lambda.y * moment.y + lambda.z * moment.z;
+
+    EXPECT_NEAR(matrix[1].imag(), -2.0, 1e-15);
+    EXPECT_NEAR(matrix[2].imag(), 2.0, 1e-15);
+    EXPECT_NEAR(expectation, dot_moment, 1e-15);
+}
+
 // =====================================================================
 // 2. calculate_delta_hcc: Pauli matrix expansion
 //

source/source_lcao/module_operator_lcao/dspin_lcao.cpp

@@ -1,10 +1,12 @@
 #include "dspin_lcao.h"
-#include "source_lcao/module_deltaspin/spin_constrain.h"
-#include "source_base/timer.h"
+
 #include "source_base/memory_recorder.h"
-#include "source_base/tool_title.h"
 #include "source_base/parallel_reduce.h"
+#include "source_base/timer.h"
+#include "source_base/tool_title.h"
 #include "source_io/module_parameter/parameter.h"
+#include "source_lcao/module_deltaspin/mi_tools.h"
+#include "source_lcao/module_deltaspin/spin_constrain.h"
 
 template <typename TK, typename TR>
 hamilt::DeltaSpin<hamilt::OperatorLCAO<TK, TR>>::DeltaSpin(HS_Matrix_K<TK>* hsk_in,
@@ -56,16 +58,13 @@ inline void cal_coeff_lambda(const std::vector<double>& current_lambda, std::vec
     coefficients[1] = -current_lambda[0];
 }
 inline void cal_coeff_lambda(const std::vector<double>& current_lambda, std::vector<std::complex<double>>& coefficients)
-{// {\lambda^{I,3}, \lambda^{I,1}+i\lambda^{I,2}, \lambda^{I,1}-i\lambda^{I,2}, -\lambda^{I,3}}
-    // The occupation matrix is built conj-first (occ[1]=conj(c_up)*c_dn, spin_constrain.cpp),
-    // and pauli_to_moment measures the physical +m_y from it (bare +Im). The lambda operator must
-    // drive that measured moment with the matching feedback sign, i.e. the pre-#7664 convention
-    // lambda_{ud}=lambda_x+i*lambda_y. #7664 flipped this together with the measurement (mirror-y);
-    // #7748 reverted the measurement but not this, leaving the constraint loop driving the y-mirror.
-    coefficients[0] = std::complex<double>(current_lambda[2], 0.0);
-    coefficients[1] = std::complex<double>(current_lambda[0] , current_lambda[1]);
-    coefficients[2] = std::complex<double>(current_lambda[0] , -1 * current_lambda[1]);
-    coefficients[3] = std::complex<double>(-1 * current_lambda[2], 0.0);
+{
+    const ModuleBase::Vector3<double> lambda(current_lambda[0], current_lambda[1], current_lambda[2]);
+    const auto spinor = spinconstrain::pauli_vector_to_spinor(lambda);
+    for (int is = 0; is < 4; ++is)
+    {
+        coefficients[is] = spinor[is];
+    }
 }
 
 template <typename TK, typename TR>
@@ -635,4 +634,4 @@ void hamilt::DeltaSpin<hamilt::OperatorLCAO<TK, TR>>::cal_PI_sub(
 
 template class hamilt::DeltaSpin<hamilt::OperatorLCAO<double, double>>;
 template class hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, double>>;
-template class hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>>;
+template class hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>>;