refactoring

src/inih/cpp/INIReader.cpp

@@ -81,13 +81,13 @@ void INIReader::dump_compilation_info() const {
 void INIReader::dump_env_info() const {
   auto log = spdlog::get("loggers");
   log->debug("-----------enviroment variables------------");
-  const std::vector<std::string> env_variables = {
+  const std::vector<std::string> env_variables{
       "OMP_DYNAMIC",  "OMP_SCHEDULE",  "OMP_NUM_THREADS", "MKL_NUM_THREADS",
       "KMP_AFFINITY", "OMP_PROC_BIND", "OMP_PLACES",      "GOMP_CPU_AFFINITY"};
-  const std::vector<std::string> slurm_vars = {
+  const std::vector<std::string> slurm_vars{
       "SLURM_JOB_ID", "SLURM_SUBMIT_DIR", "SLURM_NTASKS", "SLURM_JOB_NODELIST"};
-  const std::vector<std::string> pbs_vars = {"PBS_JOBID", "PBS_O_WORKDIR",
-                                             "PBS_NP", "PBS_NODEFILE"};
+  const std::vector<std::string> pbs_vars{"PBS_JOBID", "PBS_O_WORKDIR",
+                                          "PBS_NP", "PBS_NODEFILE"};
 
   if (getenv(slurm_vars.at(0).c_str())) {
     for (const auto &var : slurm_vars) {
@@ -173,7 +173,8 @@ void INIReader::dump_parsed() const {
                   "supported parameters.",
                   param_in_file);
       } else {
-        log->warn("Unrecognized parameter in the input file: " + param_in_file);
+        log->warn("Unrecognized parameter in the input file: {}",
+                  param_in_file);
       }
     }
   }

src/isolated.cpp

@@ -19,7 +19,7 @@ std::vector<double> nonlinear_fit(const double &opt_tol,
   try {
     opt.optimize(fit_parameters, fit_MSE);
   } catch (const std::exception &e) {
-    log->error("Nonlinear fitting failed: " + std::string(e.what()));
+    log->error("Nonlinear fitting failed: {}", std::string(e.what()));
   }
 
   return fit_parameters;

src/slabcc.cpp

@@ -3,15 +3,12 @@
 // See the accompanying LICENSE.txt file for terms.
 
 #include "isolated.hpp"
+#include "slabcc_consts.hpp"
 #include "slabcc_model.hpp"
 #include "stdafx.h"
 #include "vasp.hpp"
 
 int verbosity_level = 0;
-const double ang_to_bohr = 1e-10 / arma::datum::a_0; // 1.88972612546
-const double Hartree_to_eV = arma::datum::R_inf * arma::datum::h *
-                             arma::datum::c_0 / arma::datum::eV *
-                             2; // 27.2113860193
 
 int main(int argc, char *argv[]) {
   slabcc_model model;
@@ -144,24 +141,24 @@ int main(int argc, char *argv[]) {
   check_slabcc_compatiblity(Neutral_supercell, Charged_supercell);
 
   // cell vectors of the CHGCAR and LOCPOT files (bohr)
-  const arma::mat33 input_cell_vectors = abs(Neutral_supercell.cell_vectors) *
-                                         Neutral_supercell.scaling *
-                                         ang_to_bohr;
+  const arma::mat33 input_cell_vectors =
+      arma::abs(Neutral_supercell.cell_vectors) * Neutral_supercell.scaling *
+      ang_to_bohr;
   const arma::urowvec3 input_grid_size = SizeVec(Neutral_supercell.charge);
   model.init_supercell(input_cell_vectors, input_grid_size);
 
   const arma::rowvec3 relative_shift = 0.5 - slabcenter;
   model.rounded_relative_shift =
-      round(model.cell_grid % relative_shift) / model.cell_grid;
+      arma::round(model.cell_grid % relative_shift) / model.cell_grid;
 
   model.interfaces = fmod(
       model.interfaces + model.rounded_relative_shift(normal_direction), 1);
 
   if (!output_diffs_only) {
     Neutral_supercell.shift(model.rounded_relative_shift);
     Charged_supercell.shift(model.rounded_relative_shift);
-    model.charge_position +=
-        repmat(model.rounded_relative_shift, model.charge_position.n_rows, 1);
+    model.charge_position += arma::repmat(model.rounded_relative_shift,
+                                          model.charge_position.n_rows, 1);
     model.charge_position = fmod_p(model.charge_position, 1);
   }
   log->debug("Slab normal direction index (0-2): {}", model.normal_direction);
@@ -236,9 +233,9 @@ int main(int argc, char *argv[]) {
     const arma::rowvec3 optimization_grid_size =
         opt_grid_x * arma::conv_to<arma::rowvec>::from(model.cell_grid);
     const arma::urowvec3 optimization_grid = {
-        (arma::uword)optimization_grid_size(0),
-        (arma::uword)optimization_grid_size(1),
-        (arma::uword)optimization_grid_size(2)};
+        static_cast<arma::uword>(optimization_grid_size(0)),
+        static_cast<arma::uword>(optimization_grid_size(1)),
+        static_cast<arma::uword>(optimization_grid_size(2))};
     model.change_grid(optimization_grid);
     model.update_V_target();
     model.optimize(opt_algo, opt_tol, max_eval, max_time,
@@ -405,9 +402,9 @@ int main(int argc, char *argv[]) {
     const arma::rowvec3 extrapolation_grid_size =
         extrapol_grid_x * arma::conv_to<arma::rowvec>::from(model.cell_grid);
     const arma::urowvec3 extrapolation_grid = {
-        (arma::uword)extrapolation_grid_size(0),
-        (arma::uword)extrapolation_grid_size(1),
-        (arma::uword)extrapolation_grid_size(2)};
+        static_cast<arma::uword>(extrapolation_grid_size(0)),
+        static_cast<arma::uword>(extrapolation_grid_size(1)),
+        static_cast<arma::uword>(extrapolation_grid_size(2))};
     model.change_grid(extrapolation_grid);
     model.adjust_extrapolation_grid(extrapol_steps_num, extrapol_steps_size);
     if (as_size(model.cell_grid) !=
@@ -475,10 +472,11 @@ int main(int argc, char *argv[]) {
       E_isolated = cs.at(0) + (cs.at(1) - madelung_term) / cs.at(3);
       E_correction = E_isolated - EperModel0 - model.total_charge * dV;
     } else { // bulk and slab models
-      const arma::colvec pols = polyfit(sizes, Es, 1);
-      const arma::colvec evals = polyval(pols, sizes.t());
-      const auto linearfit_MSE = accu(square(evals.t() - Es)) / Es.n_elem * 100;
-      const arma::rowvec slopes = diff(Es) / diff(sizes);
+      const arma::colvec pols = arma::polyfit(sizes, Es, 1);
+      const arma::colvec evals = arma::polyval(pols, sizes.t());
+      const auto linearfit_MSE =
+          arma::accu(arma::square(evals.t() - Es)) / Es.n_elem * 100;
+      const arma::rowvec slopes = arma::diff(Es) / arma::diff(sizes);
       const auto extrapol_error_periodic =
           abs(slopes(0) - slopes(slopes.n_elem - 1));
       log->debug("--------------------------------------------------------");

src/slabcc_consts.hpp

@@ -0,0 +1,9 @@
+#pragma once
+
+#include <armadillo>
+
+const double ang_to_bohr = 1e-10 / arma::datum::a_0; // 1.88972612546
+const double Hartree_to_eV = arma::datum::R_inf * arma::datum::h *
+                             arma::datum::c_0 / arma::datum::eV *
+                             2; // 27.2113860193
+const double PI = arma::datum::pi;

src/slabcc_input.cpp

@@ -6,7 +6,7 @@
 
 void input_data::verify() const {
   auto log = spdlog::get("loggers");
-  charge_sigma = abs(charge_sigma);
+  charge_sigma = arma::abs(charge_sigma);
   max_eval = abs(max_eval);
   max_time = abs(max_time);
   interfaces = fmod_p(interfaces, 1);
@@ -101,16 +101,17 @@ void input_data::verify() const {
 
   if (charge_sigma.n_elem == 1) {
     charge_sigma = charge_sigma(0) * arma::ones(charge_number, 3);
-    log->debug("Only one charge_sigma is defined!");
-
+    if (trivariate) {
+      log->debug("Only one charge_sigma is defined!");
+    }
   } else if (sigma_rows == charge_number) {
     if (sigma_cols == 3) {
       if (trivariate) {
         log->debug("All the charge_sigma values are properly defined!");
       } else {
         const arma::mat isotropic_sig = arma::repmat(charge_sigma.col(0), 1, 3);
         const arma::mat sig_diff = arma::abs(isotropic_sig - charge_sigma);
-        if (any(vectorise(sig_diff) > 0.01)) {
+        if (arma::any(arma::vectorise(sig_diff) > 0.01)) {
           charge_sigma = isotropic_sig;
           log->warn("charge_sigma is not defined properly! charge_sigma={} "
                     "will be used.",

src/slabcc_math.cpp

@@ -112,7 +112,6 @@ arma::vec planar_average(const arma::uword &direction,
 }
 
 arma::cx_vec fft(arma::vec X) {
-  // TODO: should come up with a better solution than reinterpret_cast
   arma::cx_vec out(X.n_elem);
   fftw_plan plan = fftw_plan_dft_r2c_1d(
       X.n_elem, X.memptr(), reinterpret_cast<fftw_complex *>(out.memptr()),

src/slabcc_math.hpp

@@ -5,12 +5,11 @@
 #pragma once
 #include "arma_io.hpp"
 #include "general_io.hpp"
+#include "slabcc_consts.hpp"
 #include "spline.hpp"
 #include <armadillo>
 #include <fftw3.h>
 
-#define PI arma::datum::pi
-
 // These functions extend the functionality of the included Armadillo library by
 // providing:
 // (arma,iostream) << and >> operators for reading from iostreams to armadillo

src/slabcc_model.cpp

@@ -123,8 +123,6 @@ void slabcc_model::gaussian_charges_gen() {
         0, cell_vectors_lengths(2) - cell_vectors_lengths(2) / cell_grid(2),
         cell_grid(2));
 
-    CHG = arma::zeros<arma::cx_cube>(as_size(cell_grid));
-
     for (arma::uword i = 0; i < charge_fraction.n_elem; ++i) {
       // shift the axis reference to position of the Gaussian charge center
       arma::rowvec x =
@@ -191,6 +189,8 @@ void slabcc_model::gaussian_charges_gen() {
       // checked!
 
       const double Q = charge_fraction(i) * defect_charge;
+      CHG = arma::zeros<arma::cx_cube>(as_size(cell_grid));
+
       if (trivariate_charge) {
         CHG += arma::cx_cube(
             Q / (pow(2 * PI, 1.5) * arma::prod(charge_sigma.row(i))) *
@@ -222,11 +222,11 @@ slabcc_model::data_packer(opt_switches optimize) const {
   const arma::rowvec3 relative_move_step =
       move_step * ang_to_bohr / cell_vectors_lengths;
   //------interfaces----
-  std::vector<double> step_size = {relative_move_step(normal_direction),
-                                   relative_move_step(normal_direction)};
-  std::vector<double> opt_param = {interfaces(0), interfaces(1)};
-  std::vector<double> low_b = {0, 0}; // lower bounds
-  std::vector<double> upp_b = {1, 1}; // upper bounds
+  std::vector<double> step_size{relative_move_step(normal_direction),
+                                relative_move_step(normal_direction)};
+  std::vector<double> opt_param{interfaces(0), interfaces(1)};
+  std::vector<double> low_b{0, 0}; // lower bounds
+  std::vector<double> upp_b{1, 1}; // upper bounds
   if (!optimize.interfaces) {
     low_b = opt_param;
     upp_b = opt_param;
@@ -458,11 +458,7 @@ bool slabcc_model::had_discretization_error() {
       log->debug("Model charge error on the former grid size: {}",
                  new_charge_error);
       last_charge_error = new_charge_error;
-      const arma::rowvec3 new_grid_size =
-          1.5 * arma::conv_to<arma::rowvec>::from(cell_grid);
-      const arma::urowvec3 new_grid = {(arma::uword)new_grid_size(0),
-                                       (arma::uword)new_grid_size(1),
-                                       (arma::uword)new_grid_size(2)};
+      const arma::urowvec3 new_grid = cell_grid + cell_grid / 2;
       change_grid(new_grid);
       log->debug("New model charge grid size: {}", to_string(cell_grid));
       return true;
@@ -843,8 +839,9 @@ void slabcc_model::verify_CHG(const arma::cube &defect_charge) {
 
     // find the grid index of the interfaces
     arma::rowvec2 interfaces_index = cell_grid(normal_direction) * interfaces;
-    arma::urowvec2 interfaces_grid_i = {(arma::uword)interfaces_index(0),
-                                        (arma::uword)interfaces_index(1)};
+    arma::urowvec2 interfaces_grid_i = {
+        static_cast<arma::uword>(interfaces_index(0)),
+        static_cast<arma::uword>(interfaces_index(1))};
     interfaces_grid_i = arma::sort(interfaces_grid_i);
     std::vector<arma::span> spans = {
         arma::span(), arma::span(),
@@ -861,8 +858,8 @@ void slabcc_model::verify_CHG(const arma::cube &defect_charge) {
     arma::rowvec2 defect_interfaces_index =
         defect_grid(normal_direction) * interfaces;
     arma::urowvec2 defect_interfaces_grid_i = {
-        (arma::uword)defect_interfaces_index(0),
-        (arma::uword)defect_interfaces_index(1)};
+        static_cast<arma::uword>(defect_interfaces_index(0)),
+        static_cast<arma::uword>(defect_interfaces_index(1))};
     defect_interfaces_grid_i = sort(defect_interfaces_grid_i);
     std::vector<arma::span> defect_spans = {
         arma::span(), arma::span(),

src/slabcc_model.hpp

@@ -1,11 +1,9 @@
 #pragma once
+#include "slabcc_consts.hpp"
 #include "slabcc_input.hpp"
 #include "slabcc_math.hpp"
 #include "vasp.hpp"
 
-extern const double Hartree_to_eV;
-extern const double ang_to_bohr;
-
 struct opt_variable {
   arma::rowvec2 &interfaces;
   arma::mat &charge_sigma, &charge_rotations;

src/vasp.hpp

@@ -3,9 +3,8 @@
 // See the accompanying LICENSE.txt file for terms.
 
 #pragma once
+#include "slabcc_consts.hpp"
 #include "slabcc_math.hpp"
-// TODO: CONSTEXPER
-extern const double ang_to_bohr;
 
 // These functions are processing the VASP files:
 // read/write POSCAR, CHGCAR, LOCPOT files