Merge pull request #109 from electronic-structure/develop

Develop

src/Band/diag_full_potential.hpp

@@ -257,7 +257,7 @@ inline void Band::get_singular_components(K_point* kp__) const
     int n = ncomp;
 
     if (ctx_.control().verbosity_ >= 3 && kp__->comm().rank() == 0) {
-        DUMP("iterative solver tolerance: %18.12f", ctx_.iterative_solver_tolerance());
+        printf("iterative solver tolerance: %18.12f\n", ctx_.iterative_solver_tolerance());
     }
 
     #ifdef __PRINT_MEMORY_USAGE
@@ -273,13 +273,43 @@ inline void Band::get_singular_components(K_point* kp__) const
         /* apply Hamiltonian and overlap operators to the new basis functions */
         apply_fv_o(kp__, true, true, N, n, phi, ophi);
 
+        if (ctx_.control().verification_ >= 1) {
+            dmatrix<double_complex> tmp;
+            set_subspace_mtrx(1, 0, N + n, phi, ophi, ovlp, tmp);
+
+            if (ctx_.control().verification_ >= 2) {
+                ovlp.serialize("overlap", N + n);
+            }
+
+            double max_diff = check_hermitian(ovlp, N + n);
+            if (max_diff > 1e-12) {
+                std::stringstream s;
+                s << "overlap matrix is not hermitian, max_err = " << max_diff;
+                TERMINATE(s);
+            }
+        }
+
         orthogonalize(N, n, phi, ophi, ovlp, res);
 
         /* setup eigen-value problem
          * N is the number of previous basis functions
          * n is the number of new basis functions */
         set_subspace_mtrx(1, N, n, phi, ophi, ovlp, ovlp_old);
 
+        if (ctx_.control().verification_ >= 1) {
+
+            if (ctx_.control().verification_ >= 2) {
+                ovlp.serialize("overlap", N + n);
+            }
+
+            double max_diff = check_hermitian(ovlp, N + n);
+            if (max_diff > 1e-12) {
+                std::stringstream s;
+                s << "overlap matrix is not hermitian, max_err = " << max_diff;
+                TERMINATE(s);
+            }
+        }
+
         /* increase size of the variation space */
         N += n;
 
@@ -293,19 +323,27 @@ inline void Band::get_singular_components(K_point* kp__) const
             s << "error in diagonalziation";
             TERMINATE(s);
         }
+
+        for (auto e: eval) {
+            if (e < 0) {
+                std::stringstream s;
+                s << "overlap matrix is not positively defined";
+                TERMINATE(s);
+            }
+        }
 
         if (ctx_.control().verbosity_ >= 3 && kp__->comm().rank() == 0) {
-            DUMP("step: %i, current subspace size: %i, maximum subspace size: %i", k, N, num_phi);
+            printf("step: %i, current subspace size: %i, maximum subspace size: %i\n", k, N, num_phi);
             if (ctx_.control().verbosity_ >= 4) {
                 for (int i = 0; i < ncomp; i++) {
-                    DUMP("eval[%i]=%20.16f, diff=%20.16f", i, eval[i], std::abs(eval[i] - eval_old[i]));
+                    printf("eval[%i]=%20.16f, diff=%20.16f\n", i, eval[i], std::abs(eval[i] - eval_old[i]));
                 }
             }
         }
 
         /* don't compute residuals on last iteration */
         if (k != itso.num_steps_ - 1) {
-            /* get new preconditionined residuals, and also hpsi and opsi as a by-product */
+            /* get new preconditionined residuals, and also opsi and psi as a by-product */
             n = residuals(kp__, 0, N, ncomp, eval, eval_old, evec, ophi, phi, opsi, psi, res, o_diag, diag1);
         }
 
@@ -322,7 +360,11 @@ inline void Band::get_singular_components(K_point* kp__) const
             }
             else { /* otherwise, set Psi as a new trial basis */
                 if (ctx_.control().verbosity_ >= 3 && kp__->comm().rank() == 0) {
-                    DUMP("subspace size limit reached");
+                    printf("subspace size limit reached\n");
+                }
+
+                if (itso.converge_by_energy_) {
+                    transform(ctx_.processing_unit(), ophi, 0, N, evec, 0, 0, opsi, 0, ncomp);
                 }
 
                 ovlp_old.zero();
@@ -331,6 +373,7 @@ inline void Band::get_singular_components(K_point* kp__) const
                 }
                 /* update basis functions */
                 phi.copy_from(psi, 0, ncomp, ctx_.processing_unit());
+                ophi.copy_from(opsi, 0, ncomp, ctx_.processing_unit());
                 /* number of basis functions that we already have */
                 N = ncomp;
             }
@@ -346,6 +389,10 @@ inline void Band::get_singular_components(K_point* kp__) const
     }
     #endif
 
+    if (ctx_.control().verbosity_ >= 2 && kp__->comm().rank() == 0) {
+        printf("lowest and highest eigen-values of the singluar components: %20.16f %20.16f\n", eval.front(), eval.back());
+    }
+
     kp__->comm().barrier();
 }
 

src/Band/diag_pseudo_potential.hpp

@@ -283,7 +283,7 @@ inline int Band::diag_pseudo_potential_davidson(K_point*       kp__,
                     break;
                 } else { /* otherwise, set Psi as a new trial basis */
                     if (ctx_.control().verbosity_ >= 3 && kp__->comm().rank() == 0) {
-                        DUMP("subspace size limit reached");
+                        printf("subspace size limit reached\n");
                     }
                     hmlt_old.zero();
                     for (int i = 0; i < num_bands; i++) {

src/Beta_projectors/beta_projectors.h

@@ -70,7 +70,7 @@ class Beta_projectors: public Beta_projectors_base<1>
                         int idxrf = atom_type.indexb(xi).idxrf;
 
                         auto z = std::pow(double_complex(0, -1), l) * fourpi / std::sqrt(ctx_.unit_cell().omega());
-                        pw_coeffs_t_[0](igkloc, atom_type.offset_lo() + xi) = z * gkvec_rlm[lm] * beta_radial_integrals.value(idxrf, iat, vs[0]);
+                        pw_coeffs_t_[0](igkloc, atom_type.offset_lo() + xi) = z * gkvec_rlm[lm] * beta_radial_integrals.value<int, int>(idxrf, iat, vs[0]);
                     }
                 }
             }

src/Beta_projectors/beta_projectors_strain_deriv.h

@@ -83,7 +83,7 @@ class Beta_projectors_strain_deriv : public Beta_projectors_base<9>
                                 if (l == 0) {
                                     auto z = fourpi / std::sqrt(ctx_.unit_cell().omega());
 
-                                    auto d1 = beta_ri0.value(idxrf, iat, gvs[0]) * (-p * y00);
+                                    auto d1 = beta_ri0.value<int, int>(idxrf, iat, gvs[0]) * (-p * y00);
 
                                     pw_coeffs_t_[mu + nu * 3](igkloc, atom_type.offset_lo() + xi) = z * d1;
                                 } else {
@@ -170,7 +170,7 @@ class Beta_projectors_strain_deriv : public Beta_projectors_base<9>
             for (int iat = 0; iat < ctx_.unit_cell().num_atom_types(); iat++) {
                 for (int l = 0; l <= lmax_beta_ + 2; l++) {
                     for (int j = 0; j < ctx_.unit_cell().atom_type(iat).mt_radial_basis_size(); j++) {
-                        tmp(j, l, iat) = beta_ri1.value(j, l, iat, gvs[0]);
+                        tmp(j, l, iat) = beta_ri1.value<int, int, int>(j, l, iat, gvs[0]);
                     }
                 }
             }
@@ -201,7 +201,7 @@ class Beta_projectors_strain_deriv : public Beta_projectors_base<9>
                                 //    gkvec_rlm[c.lm3] * c.coef * tmp(idxrf, l3, iat) * r_f[nu];
                             pw_coeffs_t_[mu + nu * 3](igkloc, atom_type.offset_lo() + xi) += z3 * z2;
 
-                            auto d2 = beta_ri0.value(idxrf, iat, gvs[0]) * (-p * gkvec_rlm[lm]);
+                            auto d2 = beta_ri0.value<int, int>(idxrf, iat, gvs[0]) * (-p * gkvec_rlm[lm]);
 
                             pw_coeffs_t_[mu + nu * 3](igkloc, atom_type.offset_lo() + xi) += z1 * d2 * zil[l];
                         }

src/Density/initial_density.hpp

@@ -22,7 +22,7 @@ inline void Density::initial_density_pseudo()
     Radial_integrals_rho_pseudo ri(unit_cell_, ctx_.pw_cutoff(), 20);
     auto v = ctx_.make_periodic_function<index_domain_t::local>([&ri](int iat, double g)
                                                                  {
-                                                                     return ri.value(iat, g);
+                                                                     return ri.value<int>(iat, g);
                                                                  });
 
     if (ctx_.control().print_checksum_ && ctx_.comm().rank() == 0) {
@@ -169,22 +169,34 @@ inline void Density::initial_density_full_pot()
     /* compute contribution from free atoms to the interstitial density */
     auto v = ctx_.make_periodic_function<index_domain_t::local>([&ri](int iat, double g)
                                                                  {
-                                                                     return ri.value(iat, g);
+                                                                     return ri.value<int>(iat, g);
                                                                  });
 
-    #ifdef __PRINT_OBJECT_CHECKSUM
-    double_complex z = mdarray<double_complex, 1>(&v[0], ctx_.gvec().num_gvec()).checksum();
-    DUMP("checksum(rho_pw): %18.10f %18.10f", z.real(), z.imag());
-    #endif
+    double v0{0};
+    if (ctx_.comm().rank() == 0) {
+        v0 = v[0].real();
+    }
+    ctx_.comm().bcast(&v0, 1, 0);
+
+    if (ctx_.control().print_checksum_) {
+        auto z = mdarray<double_complex, 1>(&v[0], ctx_.gvec().count()).checksum();
+        ctx_.comm().allreduce(&z, 1);
+        if (ctx_.comm().rank() == 0) {
+            print_checksum("rho_pw", z);
+        }
+    }
 
     /* set plane-wave coefficients of the charge density */
     std::memcpy(&rho_->f_pw_local(0), &v[0], ctx_.gvec().count() * sizeof(double_complex));
     /* convert charge deisnty to real space mesh */
     rho_->fft_transform(1);
 
-    #ifdef __PRINT_OBJECT_CHECKSUM
-    DUMP("checksum(rho_rg): %18.10f", rho_->checksum_rg());
-    #endif
+    if (ctx_.control().print_checksum_) {
+        auto cs = rho_->checksum_rg();
+        if (ctx_.comm().rank() == 0) {
+            print_checksum("rho_rg", cs);
+        }
+    }
 
     /* remove possible negative noise */
     for (int ir = 0; ir < ctx_.fft().local_size(); ir++) {
@@ -201,14 +213,16 @@ inline void Density::initial_density_full_pot()
         int ig = ctx_.gvec().offset() + igloc;
         /* index of the G-vector shell */
         int igsh = ctx_.gvec().shell(ig);
-        if (gsh_map.count(igsh) == 0) gsh_map[igsh] = std::vector<int>();
+        if (gsh_map.count(igsh) == 0) {
+            gsh_map[igsh] = std::vector<int>();
+        }
         gsh_map[igsh].push_back(igloc);
     }
 
     /* list of G-shells for the curent MPI rank */
-    std::vector<std::pair<int, std::vector<int> > > gsh_list;
+    std::vector<std::pair<int, std::vector<int>>> gsh_list;
     for (auto& i: gsh_map) {
-        gsh_list.push_back(std::pair<int, std::vector<int> >(i.first, i.second));
+        gsh_list.push_back(std::pair<int, std::vector<int>>(i.first, i.second));
     }
 
     int lmax = 1; //ctx_.lmax_rho();
@@ -274,10 +288,12 @@ inline void Density::initial_density_full_pot()
     ctx_.comm().allreduce(znulm.at<CPU>(), (int)znulm.size());
     t3.stop();
 
-    #ifdef __PRINT_OBJECT_CHECKSUM
-    double_complex z3 = znulm.checksum();
-    DUMP("checksum(znulm): %18.10f %18.10f", std::real(z3), std::imag(z3));
-    #endif
+    if (ctx_.control().print_checksum_) {
+        double_complex z3 = znulm.checksum();
+        if (ctx_.comm().rank() == 0) {
+            print_checksum("znulm", z3);
+        }
+    }
 
     sddk::timer t4("sirius::Density::initial_density|rholm");
 
@@ -303,9 +319,16 @@ inline void Density::initial_density_full_pot()
         }
         for (int ir = 0; ir < unit_cell_.atom(ia).num_mt_points(); ir++) {
             double x = unit_cell_.atom(ia).radial_grid(ir);
-            rhoylm(0, ir) += (v[0] - unit_cell_.atom(ia).type().free_atom_density(x)) / y00;
+            rhoylm(0, ir) += (v0 - unit_cell_.atom(ia).type().free_atom_density(x)) / y00;
         }
         auto rhorlm = convert(rhoylm);
+        if (ctx_.control().print_checksum_) {
+            std::stringstream s;
+            s << "rhorlm(" << ia << ")";
+            auto cs = rhorlm.checksum();
+            print_checksum(s.str(), cs);
+        }
+
         for (int ir = 0; ir < unit_cell_.atom(ia).num_mt_points(); ir++) {
             for (int lm = 0; lm < lmmax; lm++) {
                 rho_->f_mt<index_domain_t::local>(lm, ir, ialoc) = rhorlm(lm, ir);

src/Geometry/forces.h

@@ -326,7 +326,7 @@ class Forces_PS
                 vector3d<double> gvec_cart = gvecs.gvec_cart(ig);
 
                 /* scalar part of a force without multipying by G-vector */
-                double_complex z = fact * fourpi * ri.value(iat, gvecs.gvec_len(ig)) *
+                double_complex z = fact * fourpi * ri.value<int>(iat, gvecs.gvec_len(ig)) *
                         std::conj(xc_pot->f_pw_local(igloc) * ctx_.gvec_phase_factor( ig, ia));
 
                 /* get force components multiplying by cartesian G-vector */
@@ -450,7 +450,7 @@ class Forces_PS
                 vector3d<double> gvec_cart = gvec.gvec_cart(ig);
 
                 /* scalar part of a force without multipying by G-vector */
-                double_complex z = fact * fourpi * ri.value(iat, gvec.gvec_len(ig)) *
+                double_complex z = fact * fourpi * ri.value<int>(iat, gvec.gvec_len(ig)) *
                         std::conj(dveff.f_pw_local(igloc) * ctx_.gvec_phase_factor(ig, ia));
 
                 /* get force components multiplying by cartesian G-vector */

src/Geometry/stress.h

@@ -808,7 +808,7 @@ class Stress {
 
         auto drhoc = ctx_.make_periodic_function<index_domain_t::local>([&ri_dg](int iat, double g)
                                                                         {
-                                                                            return ri_dg.value(iat, g);
+                                                                            return ri_dg.value<int>(iat, g);
                                                                         });
         double sdiag{0};
         int ig0 = (ctx_.comm().rank() == 0) ? 1 : 0;

src/K_point/initialize.hpp

@@ -142,15 +142,15 @@ inline void K_point::initialize()
             } else {
                 int ncomp = ctx_.iterative_solver_input().num_singular_;
                 if (ncomp < 0) {
-                    ncomp = ctx_.num_fv_states();
+                    ncomp = ctx_.num_fv_states() / 2;
                 }
 
                 singular_components_ = std::unique_ptr<wave_functions>(new wave_functions(ctx_.processing_unit(), gkvec(), ncomp));
                 singular_components_->pw_coeffs().prime().zero();
                 /* starting guess for wave-functions */
                 for (int i = 0; i < ncomp; i++) {
-                    for (int igloc = 0; igloc < gkvec().gvec_count(comm().rank()); igloc++) {
-                        int ig = igloc + gkvec().gvec_offset(comm().rank());
+                    for (int igloc = 0; igloc < gkvec().count(); igloc++) {
+                        int ig = igloc + gkvec().offset();
                         if (ig == i) {
                             singular_components_->pw_coeffs().prime(igloc, i) = 1.0;
                         }
@@ -160,7 +160,13 @@ inline void K_point::initialize()
                         if (ig == i + 2) {
                             singular_components_->pw_coeffs().prime(igloc, i) = 0.125;
                         }
-                        singular_components_->pw_coeffs().prime(igloc, i) += 0.01 * type_wrapper<double_complex>::random();
+                        //singular_components_->pw_coeffs().prime(igloc, i) += 0.01 * type_wrapper<double_complex>::random();
+                    }
+                }
+                if (ctx_.control().print_checksum_) {
+                    auto cs = singular_components_->checksum_pw(0, ncomp, ctx_.processing_unit());
+                    if (comm().rank() == 0) {
+                        print_checksum("singular_components", cs);
                     }
                 }
             }

src/SDDK/Makefile

@@ -4,7 +4,7 @@ include ../../make.inc
 
 SRC = sddk_api.cpp
 
-OBJ = $(SRC:.cpp=.o)
+OBJ = $(SRC:.cpp=.o) sddk.o
 
 ifneq (,$(findstring D__GPU,$(CXX_OPT)))
 	OBJ := $(OBJ) \

src/SDDK/matrix_storage.hpp

@@ -402,7 +402,7 @@ class matrix_storage<T, matrix_storage_t::slab>
         return extra_;
     }
 
-  template <sddk::memory_t mem_type>
+    template <sddk::memory_t mem_type>
     inline void zero(int i0__, int n__)
     {
         prime_.template zero<mem_type>(i0__ * num_rows_loc(), n__ * num_rows_loc());