update nclibxc and add input parameters

Author: YuleZhang936

source/module_hamilt_general/module_xc/CMakeLists.txt

@@ -24,6 +24,13 @@ add_library(
     NCLibxc/LebedevGrid.cpp
     NCLibxc/FibonacciGrid.cpp
     NCLibxc/FibonacciGrid.h
+    NCLibxc/benchmark_tests.cpp
+    NCLibxc/gga.cpp
+    NCLibxc/lda.cpp
+    NCLibxc/mgga.cpp
+    NCLibxc/math.cpp
+    NCLibxc/print.cpp
+    NCLibxc/torque.cpp
 )
 
 if(ENABLE_COVERAGE)

source/module_hamilt_general/module_xc/NCLibxc/NCLibxc.cpp

@@ -0,0 +1,44 @@
+#include "NCLibxc.h"
+#include "interface_to_libxc.h"
+#include <vector>
+#include <xc.h>
+
+
+///////////////////////////////////////////////////////////////////////////////////lda
+
+// post-processing of libxc. get partial derivatives from libxc and integrate them to get the 0th,1st,2nd derivatives of the functional
+//Note that e is the exchange and correlation energy per electron per volume. You need to multiply by \rho before the integration.
+/*
+v1=0, v2=1, f1=0,0, f2=0,1, f3=1,1
+*/
+
+void NCLibxc::postlibxc_lda(int xc_id, const std::vector<double>& rho_up, const std::vector<double>& rho_down, 
+                   std::vector<double>& e, std::vector<double>& v1, std::vector<double>& v2, 
+                   std::vector<double>& f1, std::vector<double>& f2, std::vector<double>& f3)
+{
+    LibxcInterface libxc(xc_id, true); // xc_id now passed from the caller
+
+    std::vector<double> exc = libxc.lda_exc(rho_up, rho_down);
+    std::vector<double> vrho_1(rho_up.size()), vrho_2(rho_down.size());
+    libxc.lda_vxc(rho_up, rho_down, vrho_1, vrho_2);
+
+    std::vector<double> v2rho2_1(rho_up.size()), v2rho2_2(rho_down.size()), v2rho2_3(rho_up.size());
+    libxc.lda_fxc(rho_up, rho_down, v2rho2_1, v2rho2_2, v2rho2_3);
+
+    
+    e.resize(rho_up.size());
+    v1.resize(rho_up.size());
+    v2.resize(rho_down.size());
+    f1.resize(rho_up.size());
+    f2.resize(rho_down.size());
+    f3.resize(rho_up.size());
+
+    for (size_t i = 0; i < rho_up.size(); ++i) {
+        e[i] = exc[i] ; 
+        v1[i] = vrho_1[i] ; 
+        v2[i] = vrho_2[i] ; 
+        f1[i] = v2rho2_1[i] ; 
+        f2[i] = v2rho2_2[i] ; 
+        f3[i] = v2rho2_3[i] ; 
+    }
+}

source/module_hamilt_general/module_xc/NCLibxc/NCLibxc.h

@@ -0,0 +1,72 @@
+#include "NCLibxc.h"
+#include <iostream>
+#include <vector>
+#include <array>
+#include <cmath>
+#include <complex>
+#include <iomanip>
+#include <stdexcept>
+#include <fstream>
+
+
+///////////////////////////////////////////////////////////////////////////////////
+//Matrix
+
+
+// Pauli matrix initialization
+const Matrix2x2 NCLibxc::sigma_x = {{{std::complex<double>(0.0, 0.0), std::complex<double>(1.0, 0.0)},
+                            {std::complex<double>(1.0, 0.0), std::complex<double>(0.0, 0.0)}}};
+
+const Matrix2x2 NCLibxc::sigma_y = {{{std::complex<double>(0.0, 0.0), std::complex<double>(0.0, -1.0)},
+                            {std::complex<double>(0.0, 1.0), std::complex<double>(0.0, 0.0)}}};
+
+const Matrix2x2 NCLibxc::sigma_z = {{{std::complex<double>(1.0, 0.0), std::complex<double>(0.0, 0.0)},
+                            {std::complex<double>(0.0, 0.0), std::complex<double>(-1.0, 0.0)}}};
+
+// matrix addition
+Matrix2x2 NCLibxc::add(const Matrix2x2 &a, const Matrix2x2 &b)
+{
+    Matrix2x2 result;
+    for (int i = 0; i < 2; ++i)
+    {
+        for (int j = 0; j < 2; ++j)
+        {
+            result[i][j] = a[i][j] + b[i][j];
+        }
+    }
+    return result;
+}
+
+// matrix scalar multiplication
+Matrix2x2 NCLibxc::scalar_multiply(const std::complex<double> &scalar, const Matrix2x2 &matrix)
+{
+    Matrix2x2 result;
+    for (int i = 0; i < 2; ++i)
+    {
+        for (int j = 0; j < 2; ++j)
+        {
+            result[i][j] = scalar * matrix[i][j];
+        }
+    }
+    return result;
+}
+
+
+// construct a Pauli matrix based on x, y, z
+Matrix2x2 NCLibxc::construct_pauli_matrix(double x, double y, double z)
+{
+    Matrix2x2 pauli_matrix = add(add(scalar_multiply(x, sigma_x), scalar_multiply(y, sigma_y)), scalar_multiply(z, sigma_z));
+    return pauli_matrix;
+}
+
+// identity matrix
+Matrix2x2 NCLibxc::identity_matrix()
+{
+    Matrix2x2 result;
+    result[0][0] = 1.0;
+    result[0][1] = 0.0;
+    result[1][0] = 0.0;
+    result[1][1] = 1.0;
+    return result;
+}
+///////////////////////////////////////////////////////////////////////////////////