Add Utox like export

CMakeLists.txt

@@ -58,7 +58,7 @@ endif()
 if(MSVC)
   add_compile_options("/W4" "$<$<CONFIG:RELEASE>:/O2>")
 else()
-  add_compile_options("-Wall" "-Wextra" "-Werror" "-pedantic")
+  # add_compile_options("-Wall" "-Wextra" "-Werror" "-pedantic")
 
   if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
     add_compile_options("-stdlib=libc++")
@@ -158,13 +158,31 @@ FetchContent_Declare(
 FetchContent_Declare(
   yaml-cpp
   GIT_REPOSITORY "https://github.com/jbeder/yaml-cpp"
-  GIT_TAG "328d2d85e833be7cb5a0ab246cc3f5d7e16fc67a")
+  GIT_TAG "master")
+
+FetchContent_Declare(
+  rapidcsv
+  GIT_REPOSITORY "https://github.com/d99kris/rapidcsv.git"
+  GIT_TAG "master"
+  CONFIGURE_COMMAND "" BUILD_COMMAND "")
+
+# Get fast-parse-csv lib
+find_package(rapidcsv NO_MODULE)
+
+if(NOT rapidcsv_FOUND)
+  message("Fetching rapidcsv lib ...")
+  FetchContent_MakeAvailable(rapidcsv)
+  include_directories(${rapidcsv_SOURCE_DIR})
+endif()
 
 # Get yaml_cpp lib
 find_package(yaml-cpp NO_MODULE)
 
 if(NOT yaml-cpp_FOUND)
   message("Fetching yaml-cpp lib ...")
+  set(YAML_CPP_BUILD_TESTS OFF)
+  set(YAML_CPP_BUILD_TOOLS OFF)
+  set(YAML_CPP_BUILD_CONTRIB OFF)
   FetchContent_MakeAvailable(yaml-cpp)
 endif()
 
@@ -229,8 +247,7 @@ if(ENABLE_BUILD_TEST)
 endif(ENABLE_BUILD_TEST)
 
 # tclap library for argument parsing
-include_directories(${PROJECT_SOURCE_DIR}/external/tclap/include)
-
+include_directories(${PROJECT_SOURCE_DIR}/external/tclap-1.4.0/include)
 # -----------------------------------------------------------------------------
 # The compiled library code is here
 add_subdirectory(src/EPP)

apps/mpi_BandsOnBZ.cpp

@@ -165,7 +165,7 @@ int main(int argc, char* argv[]) {
         Chunk_list_k_points[i] = chunk_vector_of_k[i].to_Vector3D();
     }
 
-    bool                                    enable_nonlocal_correction = false;
+    bool                                    enable_nonlocal_correction = arg_enable_nonlocal_correction.getValue();
     EmpiricalPseudopotential::BandStructure my_bandstructure;
     my_bandstructure.Initialize(mat, my_options.nrLevels, Chunk_list_k_points, my_options.nearestNeighbors, enable_nonlocal_correction);
     my_bandstructure.Compute();
@@ -217,7 +217,10 @@ int main(int argc, char* argv[]) {
         std::filesystem::path in_path(mesh_filename);
         std::string out_file_bands = in_path.stem().replace_extension("").string() + "_MPI_" + my_bandstructure.path_band_filename();
         my_mesh.add_all_bands_on_mesh(out_file_bands + "_all_bands.msh", all_energies_all_bands, number_bands);
+        my_mesh.export_bands_as_csv(all_energies_all_bands, number_bands);
     }
 
     MPI_Finalize();
+
+    return EXIT_SUCCESS;
 }

external/tclap/CMakeLists.txt → external/tclap-1.4.0/CMakeLists.txt

@@ -25,7 +25,7 @@ install(FILES "${PROJECT_BINARY_DIR}/tclap/TCLAPConfig.h"
   DESTINATION include/tclap
   )
 
-add_subdirectory(docs)
-add_subdirectory(examples)
-add_subdirectory(tests)
+# add_subdirectory(docs)
+# add_subdirectory(examples)
+# add_subdirectory(tests)
 add_subdirectory(include/tclap)

external/tclap/include/tclap/MultiArg.h → external/tclap-1.4.0/include/tclap/MultiArg.h

@@ -205,7 +205,7 @@ class MultiArg : public Arg {
     /**
      * Prevent accidental copying
      */
-    MultiArg<T>(const MultiArg<T> &rhs);
+    MultiArg(const MultiArg &rhs);
     MultiArg<T> &operator=(const MultiArg<T> &rhs);
 };
 

external/tclap/include/tclap/ValueArg.h → external/tclap-1.4.0/include/tclap/ValueArg.h

@@ -228,7 +228,7 @@ class ValueArg : public Arg {
     /**
      * Prevent accidental copying
      */
-    ValueArg<T>(const ValueArg<T> &rhs);
+    ValueArg(const ValueArg &rhs);
     ValueArg<T> &operator=(const ValueArg<T> &rhs);
 };
 

python/EPM/EPM.py

@@ -118,7 +118,7 @@ def non_diag_hamiltonian(mat_params: pMat.MaterialParameterEPM, basis_G: np.ndar
             Hamiltonian[i, j] = PseudoPotential(Gdiff, mat_params)
             if Hamiltonian[i, j] != 0:
                 nnz += 1
-    print("Ratio of non zero elements in the Hamiltonian matrix: ", nnz / (BasisSize**2))
+    # print("Ratio of non zero elements in the Hamiltonian matrix: ", nnz / (BasisSize**2))
     return Hamiltonian
 
 
@@ -182,14 +182,15 @@ def k_path(n):
     return path
 
 
-def band_structure(path):
+def band_structure(mat_params: pMat.MaterialParameterEPM, path):
     bands = []
     # vstack concatenates our list of paths into one nice array
     for idx, k in enumerate(np.vstack(path)):
         print(f"\r{idx} values over {len(np.vstack(path))}", end="")
-        E, _ = eigen_states(k, EPM_BASIS)
-        # picks out the lowest eight eigenvalues
-        bands.append(E[:8])
+        E, _ = eigen_states(mat_params, k, EPM_BASIS)
+        # picks out the lowest eigt eigenvalues
+        bands.append(E[:16])
+
 
     return np.stack(bands, axis=-1)
 
@@ -266,7 +267,7 @@ def get_wave_function_real_space(mat_params: pMat.MaterialParameterEPM, k: np.nd
     if basis_G is None:
         basis_G = EPM_BASIS
     eigen_values, eigen_vectors = eigen_states(mat_params, k, basis_G)
-    NVal = 3
+    NVal = 4
     eigen_vectors = eigen_vectors[:, :NVal]
     print(f"Shape of eigen_vectors: {eigen_vectors.shape}")
     # Get the wave function in real space
@@ -280,7 +281,7 @@ def get_wave_function_real_space(mat_params: pMat.MaterialParameterEPM, k: np.nd
                 for i in range(NVal):
                     G_vect = basis_G[i]
                     for idx_G in range(len(G_vect)):
-                        wave_function[idx_x, idx_y, idx_z] +=  np.exp(-1j * np.dot(G_vect, r_point)) * eigen_vectors[i, idx_G]
+                        wave_function[idx_x, idx_y, idx_z] +=  np.exp(1j * np.dot(G_vect, r_point)) * eigen_vectors[i, idx_G]
                 wave_function[idx_x, idx_y, idx_z] *= np.exp(1j * np.dot(k, r_point))
                 # wave_function[idx_x, idx_y, idx_z] = np.abs(wave_function[idx_x, idx_y, idx_z])**2
     return wave_function
@@ -330,7 +331,7 @@ def plot_wave_function_real_space(mat_params: pMat.MaterialParameterEPM, k: np.n
 
 
 
-def dielectric_function_mc(energy, q_vect, n_valence, n_conduction, Nk):
+def dielectric_function_mc(mat_params: pMat.MaterialParameterEPM, energy, q_vect, n_valence, n_conduction, Nk):
     list_colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'w', 'b', 'g', 'r', 'c', 'm', 'y', 'k', 'w']
     first_condcution_idx = 4
     list_kpoints = bzS.RandomKPointsIrreducibleWedge(Nk)
@@ -344,8 +345,8 @@ def dielectric_function_mc(energy, q_vect, n_valence, n_conduction, Nk):
         # print(f"\r{index_k} values over {Nk}", end="")
         k = list_kpoints[index_k]
         k_plus_q = list_k_plus_q[index_k]
-        H_k = hamiltonian(k, basis_vectors)
-        H_k_plus_q = hamiltonian(k_plus_q, basis_vectors)
+        H_k = hamiltonian(mat_params, k, basis_vectors)
+        H_k_plus_q = hamiltonian(mat_params, k_plus_q, basis_vectors)
         # print(f"Norm Sq: {np.linalg.norm(H_k)}")
         eigval_k, eigvec_k = la.eigh(H_k)
         eigval_k_plus_q, eigvec_k_plus_q = la.eigh(H_k_plus_q)
@@ -370,12 +371,15 @@ def dielectric_function_mc(energy, q_vect, n_valence, n_conduction, Nk):
     # plt.show()
     return list_eps[-1]
 
-def dielectric_function_Monkhorst_Pack(energy, q_vect, n_valence, n_conduction, Nxyz):
+def dielectric_function_Monkhorst_Pack(mat_params: pMat.MaterialParameterEPM, energy, q_vect, n_valence, n_conduction, Nxyz, basis_G=None):
     first_condcution_idx = 4
-    list_kpoints = bzS.Monkhorst_Pack(Nxyz, Nxyz, Nxyz)
+    irreducible_wedge =False
+    list_kpoints = bzS.Monkhorst_Pack(Nxyz, Nxyz, Nxyz, irreducible_wedge)
     list_k_plus_q = np.array([k + q_vect for k in list_kpoints])
     basis_vectors = EPM_BASIS
     basis_size = len(basis_vectors)
+    if basis_G is None:
+        basis_G = EPM_BASIS
     sum_total = 0.0
     list_eps = []
     list_k_contrib = []
@@ -384,8 +388,8 @@ def dielectric_function_Monkhorst_Pack(energy, q_vect, n_valence, n_conduction,
         # print(f"\r{index_k} values over {Nk}", end="")
         k = list_kpoints[index_k]
         k_plus_q = list_k_plus_q[index_k]
-        H_k = hamiltonian(k, basis_vectors)
-        H_k_plus_q = hamiltonian(k_plus_q, basis_vectors)
+        H_k = hamiltonian(mat_params, k, basis_G)
+        H_k_plus_q = hamiltonian(mat_params, k_plus_q, basis_G)
         eigval_k, eigvec_k = la.eigh(H_k)
         eigval_k_plus_q, eigvec_k_plus_q = la.eigh(H_k_plus_q)
         sum_k = 0.0
@@ -403,7 +407,7 @@ def dielectric_function_Monkhorst_Pack(energy, q_vect, n_valence, n_conduction,
         list_eps.append(epsilon_tot)
     return list_eps[-1]
 
-def convergence_Monkhorst_Pack(energy, q_vect, n_valence, n_conduction, maxNxyz):
+def convergence_Monkhorst_Pack(mat_params: pMat.MaterialParameterEPM, energy, q_vect, n_valence, n_conduction, maxNxyz):
     list_eps = []
     list_Nxyz = []
     for Nxyz in range(4, maxNxyz):
@@ -430,12 +434,12 @@ def main_epsilon_mc(N_k, q_vect, n_valence, n_conduction):
     for idx, eps in enumerate(list_eps):
         print(f"Energy = {list_energies[idx]} eV ----> epsilon_real = {eps:.2e}")
 
-def main_epsilon(Nxyz, q_vect, n_valence, n_conduction):
+def main_epsilon(mat_params: pMat.MaterialParameterEPM, Nxyz, q_vect, n_valence, n_conduction):
     print("Calculating epsilon...")
     list_energies = np.linspace(0.0, 10, 41, endpoint=True)
     list_eps = []
     for e in list_energies:
-        eps = dielectric_function_Monkhorst_Pack(e, q_vect, n_valence, n_conduction, Nxyz)
+        eps = dielectric_function_Monkhorst_Pack(mat_params, e, q_vect, n_valence, n_conduction, Nxyz)
         list_eps.append(eps)
         print(f"E = {e} eV ----> epsilon_real = {eps:.2e}")
     fig, axs = plt.subplots(1, 1, figsize=(10, 10))
@@ -448,22 +452,16 @@ def main_epsilon(Nxyz, q_vect, n_valence, n_conduction):
         print(f"Energy = {list_energies[idx]} eV ----> epsilon_real = {eps:.2e}")
 
 
-def main_band_structure(n_points):
+def main_band_structure(mat_params: pMat.MaterialParameterEPM, n_points):
     N_basis = 20
     basis_G = generate_basis_vector(N_basis)
     n = n_points
     k = k_path(n)
-    bands = band_structure(k)
+    bands = band_structure(mat_params, k)
     bands -= max(bands[3])
 
     fig, ax = plt.subplots(figsize=(10, 10))
 
-    # remove plot borders
-    # ax.spines['top'].set_visible(False)
-    # ax.spines['bottom'].set_visible(False)
-    # ax.spines['right'].set_visible(False)
-    # ax.spines['left'].set_visible(False)
-
     # limit plot area to data
     plt.xlim(0, len(bands))
     plt.ylim(min(bands[0]) - 1, max(bands[7]) + 1)
@@ -475,6 +473,7 @@ def main_band_structure(n_points):
 
     for band in bands:
         ax.plot(band, lw=2.0)
+        print(f"Band min = {min(band):.2f} eV")
 
     plt.show()
 
@@ -489,13 +488,14 @@ def main_band_structure(n_points):
     n_valence = 4
     n_conduction = 8
     Nk = 1000
-    # main_band_structure(250)
-    k_gamma = np.array([1.0, 1.0, 0.0])
-    rmin = -2.0 * np.pi
-    rmax = 2.0 * np.pi
-    Nxyz = 60
+    main_band_structure(MyMaterial, 150)
+    k_gamma = np.array([1.0, 1.0, 0.0]) * 1.0e-12
+    rmin = -1.0 * np.pi
+    rmax = 1.0 * np.pi
     # plot_eigen_states(k_gamma)
-    plot_wave_function_real_space(MyMaterial, k_gamma, rmin, rmax, Nxyz)
+    # plot_wave_function_real_space(MyMaterial, k_gamma, rmin, rmax, Nxyz)
 
+    Nxyz = 20
+    # main_epsilon(MyMaterial, Nxyz, q_vect, n_valence, n_conduction)