Feature: output mag force in spin constraint calculation (#4037)

* Feature: output mag force in spin constraint calculation * add ut for print_Mag_Force
deepmodeling · Apr 22, 2024 · fbb6f2c · fbb6f2c
1 parent 6840f52
commit fbb6f2c
Show file tree

Hide file tree

Showing 7 changed files with 69 additions and 9 deletions.
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
@@ -1131,6 +1131,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
     {
         SpinConstrain<TK, psi::DEVICE_CPU>& sc = SpinConstrain<TK, psi::DEVICE_CPU>::getScInstance();
         sc.cal_MW(istep, &(this->LM), true);
+        sc.print_Mag_Force();
     }
 
     if (!GlobalV::CAL_FORCE && !GlobalV::CAL_STRESS)

diff --git a/source/module_hamilt_lcao/module_deltaspin/lambda_loop_helper.cpp b/source/module_hamilt_lcao/module_deltaspin/lambda_loop_helper.cpp
@@ -4,8 +4,8 @@
 template <>
 void SpinConstrain<std::complex<double>, psi::DEVICE_CPU>::print_termination()
 {
-    print_2d("after-optimization spin: (print in the inner loop): ", this->Mi_, this->nspin_);
-    print_2d("after-optimization lambda: (print in the inner loop): ", this->lambda_, this->nspin_);
+    print_2d("after-optimization spin (uB): (print in the inner loop): ", this->Mi_, this->nspin_);
+    print_2d("after-optimization lambda (Ry/uB): (print in the inner loop): ", this->lambda_, this->nspin_);
     std::cout << "Inner optimization for lambda ends." << std::endl;
     std::cout << "===============================================================================" << std::endl;
 }

diff --git a/source/module_hamilt_lcao/module_deltaspin/spin_constrain.cpp b/source/module_hamilt_lcao/module_deltaspin/spin_constrain.cpp
@@ -1,4 +1,5 @@
 #include "spin_constrain.h"
+#include "module_base/formatter_fmt.h"
 
 #include <cmath>
 
@@ -541,20 +542,66 @@ void SpinConstrain<FPTYPE, Device>::print_Mi(bool print)
     int nat = this->get_nat();
     if (print)
     {
+        formatter::Fmt fmt;
+        fmt.set_width(20);
+        fmt.set_precision(10);
+        fmt.set_fillChar(' ');
+        fmt.set_fixed(false);
+        fmt.set_right(true);
+        fmt.set_error(false);
+        std::cout << "Total Magnetism (uB): " << std::endl;
         for (int iat = 0; iat < nat; ++iat)
         {
             if (this->nspin_ == 2)
             {
-                std::cout << "Total Magnetism on atom: " << iat << " " << std::setprecision(10) << " (" << Mi_[iat].z << ")" << std::endl;
+                std::cout << "ATOM " << std::left << std::setw(6) << iat << fmt.format(Mi_[iat].z) << std::endl;
             }
             else if (this->nspin_ ==4)
             {
-                std::cout << "Total Magnetism on atom: " << iat << " " << std::setprecision(10) << " (" << Mi_[iat].x
-                        << ", " << Mi_[iat].y << ", " << Mi_[iat].z << ")" << std::endl;
+                std::cout << "ATOM " << std::left << std::setw(6) << iat << fmt.format(Mi_[iat].x) << fmt.format(Mi_[iat].y) << fmt.format(Mi_[iat].z) << std::endl;
             }
         }
     }
 }
 
+/// print magnetic force (defined as \frac{\delta{L}}/{\delta{Mi}} = -lambda[iat])
+template <typename FPTYPE, typename Device>
+void SpinConstrain<FPTYPE, Device>::print_Mag_Force()
+{
+    this->check_atomCounts();
+    int nat = this->get_nat();
+    formatter::Fmt fmt;
+    fmt.set_width(20);
+    fmt.set_precision(10);
+    fmt.set_fillChar(' ');
+    fmt.set_fixed(false);
+    fmt.set_right(true);
+    fmt.set_error(false);
+    std::cout << "Final optimal lambda (Ry/uB): " << std::endl;
+    for (int iat = 0; iat < nat; ++iat)
+    {
+        if (this->nspin_ == 2)
+        {
+            std::cout << "ATOM " << std::left << std::setw(6) << iat << fmt.format(lambda_[iat].z) << std::endl;
+        }
+        else if (this->nspin_ ==4)
+        {
+            std::cout << "ATOM " << std::left << std::setw(6) << iat << fmt.format(lambda_[iat].x) << fmt.format(lambda_[iat].y) << fmt.format(lambda_[iat].z) << std::endl;
+        }
+    }
+    std::cout << "Magnetic force (Ry/uB): " << std::endl;
+    for (int iat = 0; iat < nat; ++iat)
+    {
+        if (this->nspin_ == 2)
+        {
+            std::cout << "ATOM " << std::left << std::setw(6) << iat << fmt.format(-lambda_[iat].z) << std::endl;
+        }
+        else if (this->nspin_ ==4)
+        {
+            std::cout << "ATOM " << std::left << std::setw(6) << iat << fmt.format(-lambda_[iat].x) << fmt.format(-lambda_[iat].y) << fmt.format(-lambda_[iat].z) << std::endl;
+        }
+    }
+}
+
 template class SpinConstrain<std::complex<double>, psi::DEVICE_CPU>;
 template class SpinConstrain<double, psi::DEVICE_CPU>;
diff --git a/source/module_hamilt_lcao/module_deltaspin/spin_constrain.h b/source/module_hamilt_lcao/module_deltaspin/spin_constrain.h
@@ -87,6 +87,9 @@ class SpinConstrain
   /// print mi
   void print_Mi(bool print = false);
 
+  /// print magnetic force, defined as \frac{\delta{L}}/{\delta{Mi}} = -lambda[iat])
+  void print_Mag_Force();
+
   /// collect_mw from matrix multiplication result
   void collect_MW(ModuleBase::matrix& MecMulP, const ModuleBase::ComplexMatrix& mud, int nw, int isk);
 

diff --git a/source/module_hamilt_lcao/module_deltaspin/test/cal_mw_helper_test.cpp b/source/module_hamilt_lcao/module_deltaspin/test/cal_mw_helper_test.cpp
@@ -66,7 +66,8 @@ TEST_F(SpinConstrainTest, CalculateMW)
     testing::internal::CaptureStdout();
     sc.print_Mi(true);
     std::string output = testing::internal::GetCapturedStdout();
-    EXPECT_THAT(output, testing::HasSubstr("Total Magnetism on atom: 0  (2, 3, 4)"));
+    EXPECT_THAT(output, testing::HasSubstr("Total Magnetism (uB):"));
+    EXPECT_THAT(output, testing::HasSubstr("ATOM 0         2.0000000000e+00    3.0000000000e+00    4.0000000000e+00"));
 }
 
 TEST_F(SpinConstrainTest, CollectMW)
@@ -141,7 +142,8 @@ TEST_F(SpinConstrainTest, CalculateMWS2)
     testing::internal::CaptureStdout();
     sc.print_Mi(true);
     std::string output = testing::internal::GetCapturedStdout();
-    EXPECT_THAT(output, testing::HasSubstr("Total Magnetism on atom: 0  (-1)"));
+    EXPECT_THAT(output, testing::HasSubstr("Total Magnetism (uB):"));
+    EXPECT_THAT(output, testing::HasSubstr("ATOM 0        -1.0000000000e+00"));
 }
 
 TEST_F(SpinConstrainTest, CollectMWS2)

diff --git a/source/module_hamilt_lcao/module_deltaspin/test/lambda_loop_helper_test.cpp b/source/module_hamilt_lcao/module_deltaspin/test/lambda_loop_helper_test.cpp
@@ -40,10 +40,15 @@ TEST_F(SpinConstrainTest, PrintTermination)
     sc.set_sc_lambda(sc_lambda.data(), 1);
     testing::internal::CaptureStdout();
     sc.print_termination();
+    sc.print_Mag_Force();
     std::string output = testing::internal::GetCapturedStdout();
     EXPECT_THAT(output, testing::HasSubstr("Inner optimization for lambda ends."));
     EXPECT_THAT(output, testing::HasSubstr("ATOM 1         0.0000000000e+00    0.0000000000e+00    0.0000000000e+00"));
     EXPECT_THAT(output, testing::HasSubstr("ATOM 1         1.0000000000e+00    2.0000000000e+00    3.0000000000e+00"));
+    EXPECT_THAT(output, testing::HasSubstr("Final optimal lambda (Ry/uB):"));
+    EXPECT_THAT(output, testing::HasSubstr("ATOM 1         1.0000000000e+00    2.0000000000e+00    3.0000000000e+00"));
+    EXPECT_THAT(output, testing::HasSubstr("Magnetic force (Ry/uB):"));
+    EXPECT_THAT(output, testing::HasSubstr("ATOM 0        -1.0000000000e+00   -2.0000000000e+00   -3.0000000000e+00"));
 }
 
 TEST_F(SpinConstrainTest, CheckRmsStop)

diff --git a/source/module_hamilt_lcao/module_deltaspin/test/spin_constrain_test.cpp b/source/module_hamilt_lcao/module_deltaspin/test/spin_constrain_test.cpp
@@ -398,10 +398,12 @@ TYPED_TEST(SpinConstrainTest, PrintMi)
     testing::internal::CaptureStdout();
     this->sc.print_Mi(true);
     std::string output = testing::internal::GetCapturedStdout();
-    EXPECT_THAT(output, testing::HasSubstr("Total Magnetism on atom: 0  (0, 0, 0)"));
+    EXPECT_THAT(output, testing::HasSubstr("Total Magnetism (uB):"));
+    EXPECT_THAT(output, testing::HasSubstr("ATOM 0         0.0000000000e+00    0.0000000000e+00    0.0000000000e+00"));
     this->sc.set_nspin(2);
      testing::internal::CaptureStdout();
     this->sc.print_Mi(true);
     output = testing::internal::GetCapturedStdout();
-    EXPECT_THAT(output, testing::HasSubstr("Total Magnetism on atom: 0  (0)"));
+    EXPECT_THAT(output, testing::HasSubstr("Total Magnetism (uB):"));
+    EXPECT_THAT(output, testing::HasSubstr("ATOM 0         0.0000000000e+00"));
 }