Skip to content

Commit

Permalink
Gauss Subdivided quadrature rule working
Browse files Browse the repository at this point in the history
  • Loading branch information
Avirup Sircar committed Oct 19, 2024
1 parent 99e4d30 commit 1d2b818
Show file tree
Hide file tree
Showing 15 changed files with 5,283 additions and 84 deletions.
327 changes: 327 additions & 0 deletions analysis/PseudoConverter/PSPConverter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -222,6 +222,330 @@ xmltoDensityFile(std::string file_path_in, std::string file_path_out)
}
file.close();
}

void
xmltoProjectorFile(std::string file_path_in, std::string file_path_out)
{
// List of momentum values
std::vector<std::string> tag_name_parent;
tag_name_parent.push_back("PP_NONLOCAL");
std::vector<int> ang_mom_list;
for (int i = 1; i < xmlNodeChildCount(tag_name_parent, file_path_in); i++)
{
std::string pp_beta_str = "PP_BETA.";
pp_beta_str += std::to_string(i);
std::vector<std::string> tag_name;
tag_name.push_back("PP_NONLOCAL");
tag_name.push_back(pp_beta_str);
std::vector<std::string> attr_type;
std::vector<std::string> attr_value;
XmlTagReaderAttr(tag_name, file_path_in, &attr_type, &attr_value);
unsigned int index = 0;
std::string to_search = "angular_momentum";
auto it = std::find(attr_type.begin(), attr_type.end(), to_search);
if (it == attr_type.end())
{
throw std::invalid_argument(
"angular momentum attribute not found");
}
else
{
index = std::distance(attr_type.begin(), it);
ang_mom_list.push_back(std::stoi(attr_value[index]));
}
}

// Unique angular momentum values
std::vector<int> ang_mom_unique_list;
auto is_unique =
std::adjacent_find(ang_mom_list.begin(), ang_mom_list.end()) ==
ang_mom_list.end();
if (!is_unique)
{
ang_mom_unique_list = ang_mom_list;
std::sort(ang_mom_unique_list.begin(), ang_mom_unique_list.end());
auto it =
ang_mom_unique_list.erase(std::unique(ang_mom_unique_list.begin(),
ang_mom_unique_list.end()));
ang_mom_unique_list.resize(distance(ang_mom_unique_list.begin(), it));
}
else
{
ang_mom_unique_list = ang_mom_list;
}

// Multiplicity of unique angular momentum values
std::vector<int> ang_mom_multiplicity_list;
for (int i = 0; i < ang_mom_unique_list.size(); i++)
{
int count = 0;
for (int j = 0; j < ang_mom_list.size(); j++)
{
if (ang_mom_list[j] == ang_mom_unique_list[i])
{
count++;
}
}
}
// Beta index for same angular momentum
std::vector<std::vector<int>> beta_index;
for (int i = 0; i < ang_mom_unique_list.size(); i++)
{
beta_index.push_back((std::vector<int>()));
for (int j = 0; j < ang_mom_list.size(); j++)
{
if (ang_mom_list[j] == ang_mom_unique_list[i])
{
beta_index[i].push_back(j + 1);
}
}
}

// Extracting radial coordinates
std::vector<double> radial_coord;
std::vector<std::string> radial_tag;
radial_tag.push_back("PP_MESH");
radial_tag.push_back("PP_R");
radial_coord = XmlTagReaderMain(radial_tag, file_path_in);

// Extracting projector data according to angular momentum
for (int i = 0; i < ang_mom_unique_list.size(); i++)
{
std::vector<std::vector<double>> beta_values;
std::string proj_str = "/proj_l";
proj_str += std::to_string(ang_mom_unique_list[i]);
proj_str += ".dat";
for (int j = 0; j < beta_index[i].size(); j++)
{
std::string pp_beta_str = "PP_BETA.";
pp_beta_str += std::to_string(beta_index[i][j]);
std::vector<std::string> beta_tag;
beta_tag.push_back("PP_NONLOCAL");
beta_tag.push_back(pp_beta_str);
beta_values.push_back(std::vector<double>());
beta_values[j] = XmlTagReaderMain(beta_tag, file_path_in);
std::vector<double> trial =
XmlTagReaderMain(beta_tag, file_path_in);
}

std::fstream file;
file.open(file_path_out + proj_str, std::ios::out);
file << std::setprecision(15);
if (file.is_open())
{
for (int l = 0; l < radial_coord.size(); l++)
{
if (l == 0)
{
file << radial_coord[0] << " ";
for (int m = 0; m < beta_values.size(); m++)
{
if (m != (beta_values.size() - 1))
file << beta_values[m][1] / radial_coord[1] << " ";
else
file << beta_values[m][1] / radial_coord[1]
<< std::endl;
}
}
else
{
file << radial_coord[l] << " ";
for (int m = 0; m < beta_values.size(); m++)
{
if (m != (beta_values.size() - 1))
file << beta_values[m][l] / radial_coord[l] << " ";
else
file << beta_values[m][l] / radial_coord[l]
<< std::endl;
}
}
}
}
file.close();
}
}

void
xmltoDenomFile(std::string file_path_in, std::string file_path_out)
{
// Extracting Diagonal Matrix
std::vector<double> diagonal_mat;
std::vector<std::string> diagonal_tag;
diagonal_tag.push_back("PP_NONLOCAL");
diagonal_tag.push_back("PP_DIJ");
diagonal_mat = XmlTagReaderMain(diagonal_tag, file_path_in);

std::vector<std::string> tag_name_parent;
tag_name_parent.push_back("PP_NONLOCAL");
int n = xmlNodeChildCount(tag_name_parent, file_path_in) - 1;

// Writing the denom.dat
std::fstream file;
file.open(file_path_out + "/denom.dat", std::ios::out);
file << std::setprecision(12);
if (file.is_open())
{
for (int l = 0; l < diagonal_mat.size(); l++)
{
if (l != 0 & (l % n == 0))
file << std::endl;
file << diagonal_mat[l] / 2 << " ";
}
}
file.close();
}

void
xmltoSummaryFile(std::string file_path_in, std::string file_path_out)
{
// List of momentum values
std::vector<std::string> tag_name_parent;
tag_name_parent.push_back("PP_NONLOCAL");
std::vector<int> ang_mom_list;
for (int i = 1; i < xmlNodeChildCount(tag_name_parent, file_path_in); i++)
{
std::string pp_beta_str = "PP_BETA.";
pp_beta_str += std::to_string(i);
std::vector<std::string> tag_name;
tag_name.push_back("PP_NONLOCAL");
tag_name.push_back(pp_beta_str);
std::vector<std::string> attr_type;
std::vector<std::string> attr_value;
XmlTagReaderAttr(tag_name, file_path_in, &attr_type, &attr_value);
unsigned int index = 0;
std::string to_search = "angular_momentum";
auto it = std::find(attr_type.begin(), attr_type.end(), to_search);
if (it == attr_type.end())
{
throw std::invalid_argument(
"angular momentum attribute not found");
}
else
{
index = std::distance(attr_type.begin(), it);
ang_mom_list.push_back(std::stoi(attr_value[index]));
}
}
// Unique angular momentum values
std::vector<int> ang_mom_unique_list;
auto is_unique =
std::adjacent_find(ang_mom_list.begin(), ang_mom_list.end()) ==
ang_mom_list.end();
if (!is_unique)
{
ang_mom_unique_list = ang_mom_list;
std::sort(ang_mom_unique_list.begin(), ang_mom_unique_list.end());
auto it =
ang_mom_unique_list.erase(std::unique(ang_mom_unique_list.begin(),
ang_mom_unique_list.end()));
ang_mom_unique_list.resize(distance(ang_mom_unique_list.begin(), it));
}
else
{
ang_mom_unique_list = ang_mom_list;
}

// Multiplicity of unique angular momentum values
std::vector<int> ang_mom_multiplicity_list;
for (int i = 0; i < ang_mom_unique_list.size(); i++)
{
int count = 0;
for (int j = 0; j < ang_mom_list.size(); j++)
{
if (ang_mom_list[j] == ang_mom_unique_list[i])
{
count++;
}
}
ang_mom_multiplicity_list.push_back(count);
}
int row_index = 0;
int index = 0;
std::vector<std::vector<int>> out_proj_arr;
for (int i = 0; i < ang_mom_unique_list.size(); i++)
{
int l = ang_mom_unique_list[i];
for (int j = 0; j < ang_mom_multiplicity_list[i]; j++)
{
int m = -l;
for (int k = 0; k < 2 * l + 1; k++)
{
out_proj_arr.push_back((std::vector<int>()));
out_proj_arr[row_index].push_back(index);
out_proj_arr[row_index].push_back(l);
out_proj_arr[row_index].push_back(m);
m++;
row_index++;
}
index++;
}
}
// Writing the supplementary

std::fstream file;
file_path_out.append("/PseudoAtomDat");
file.open(file_path_out, std::ios::out);
if (file.is_open())
{
// Total projector
file << out_proj_arr.size() << std::endl;
// Projector data
int m = out_proj_arr.size();
int n = m==0 ? 0 : out_proj_arr[0].size();
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
file << out_proj_arr[i][j] << " ";
file << std::endl;
}

for (int i = 0; i < ang_mom_unique_list.size(); i++)
{
std::string proj_str = "proj_l";
proj_str += std::to_string(ang_mom_unique_list[i]);
file << proj_str << ".dat" << std::endl;
file << ang_mom_multiplicity_list[i] << std::endl;
}
// Name for D_ij file
file << "denom.dat" << std::endl;
}
// Orbitals
std::vector<std::string> pswfc_tag;
pswfc_tag.push_back("PP_PSWFC");
for (int i = 1; i <= xmlNodeChildCount(pswfc_tag, file_path_in); i++)
{
// Reading chi data
std::string pp_chi_str = "PP_CHI.";
pp_chi_str += std::to_string(i);
std::vector<std::string> chi_tag;
chi_tag.push_back("PP_PSWFC");
chi_tag.push_back(pp_chi_str);
std::vector<std::string> attr_type;
std::vector<std::string> attr_value;
XmlTagReaderAttr(chi_tag, file_path_in, &attr_type, &attr_value);
unsigned int index = 0;
std::string to_search = "label";
auto it =
std::find(attr_type.begin(), attr_type.end(), to_search);
if (it == attr_type.end())
{
throw std::invalid_argument(
"orbital label attribute not found");
}
else
{
index = std::distance(attr_type.begin(), it);
}
std::string orbital_string = attr_value[index];
for (auto &w : orbital_string)
{
w = tolower(w);
}
file << orbital_string + ".dat" << std::endl;
}
file.close();
}

void
xmltoOrbitalFile(std::string file_path_in, std::string file_path_out)
{ // Extracting radial coordinates
Expand Down Expand Up @@ -301,7 +625,10 @@ pseudoPotentialToDftefeParser(const std::string file_path_in,
{
xmltoLocalPotential(file_path_in, file_path_out);
xmltoDensityFile(file_path_in, file_path_out);
xmltoDenomFile(file_path_in, file_path_out);
xmltoProjectorFile(file_path_in, file_path_out);
xmltoOrbitalFile(file_path_in, file_path_out);
xmltoSummaryFile(file_path_in, file_path_out);
}

int main(int argc, char** argv)
Expand Down
22 changes: 22 additions & 0 deletions analysis/classicalEnrichmentComparison/PSP/KSDFTClassical/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@
cmake_minimum_required(VERSION 3.20)
project(poisson_problem_comprison)
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

set(MAIN_PROJECT_DIR $ENV{DFTEFE_PATH})
message(${MAIN_PROJECT_DIR})

include_directories(${MAIN_PROJECT_DIR}/src)
add_subdirectory(${MAIN_PROJECT_DIR}/src/utils ${MAIN_PROJECT_DIR}/test/basis/lib/buildUtils)
add_subdirectory(${MAIN_PROJECT_DIR}/src/basis ${MAIN_PROJECT_DIR}/test/basis/lib/buildBasis)
add_subdirectory(${MAIN_PROJECT_DIR}/src/atoms ${MAIN_PROJECT_DIR}/test/basis/lib/buildAtoms)
add_subdirectory(${MAIN_PROJECT_DIR}/src/quadrature ${MAIN_PROJECT_DIR}/test/basis/lib/buildQuadrature)
add_subdirectory(${MAIN_PROJECT_DIR}/src/linearAlgebra ${MAIN_PROJECT_DIR}/test/linearAlgebra/lib/buildLinearAlgebra)
add_subdirectory(${MAIN_PROJECT_DIR}/src/electrostatics ${MAIN_PROJECT_DIR}/test/electrostatics/lib/electrostatics)
add_subdirectory(${MAIN_PROJECT_DIR}/src/ksdft ${MAIN_PROJECT_DIR}/test/ksdft/lib/ksdft)

if(ENABLE_MPI)
add_compile_definitions(DFTEFE_WITH_MPI)
add_executable(TestKohnShamDft ${MAIN_PROJECT_DIR}/analysis/classicalEnrichmentComparison/PSP/KSDFTClassical/TestKohnShamDft.cpp)
target_link_libraries(TestKohnShamDft PUBLIC dft-efe-basis dft-efe-electrostatics dft-efe-ksdft dft-efe-utils dft-efe-quadrature dft-efe-linalg dft-efe-atoms)
endif()
Loading

0 comments on commit 1d2b818

Please sign in to comment.