<commit_before>#include <sirius.h>

using namespace sirius;

void init_wf(K_point* kp__, Wave_functions& phi__, int num_bands__, int num_mag_dims__)
{
    std::vector<double> tmp(0xFFFF);
    for (int i = 0; i < 0xFFFF; i++) {
        tmp[i] = utils::random<double>();
    }

    phi__.pw_coeffs(0).prime().zero();

    //#pragma omp parallel for schedule(static)
    for (int i = 0; i < num_bands__; i++) {
        for (int igk_loc = 0; igk_loc < kp__->num_gkvec_loc(); igk_loc++) {
            /* global index of G+k vector */
            int igk = kp__->idxgk(igk_loc);
            if (igk == 0) {
                phi__.pw_coeffs(0).prime(igk_loc, i) = 1.0;
            }
            if (igk == i + 1) {
                phi__.pw_coeffs(0).prime(igk_loc, i) = 0.5;
            }
            if (igk == i + 2) {
                phi__.pw_coeffs(0).prime(igk_loc, i) = 0.25;
            }
            if (igk == i + 3) {
                phi__.pw_coeffs(0).prime(igk_loc, i) = 0.125;
            }
            phi__.pw_coeffs(0).prime(igk_loc, i) += tmp[(igk + i) & 0xFFFF] * 1e-5;
        }
    }

    if (num_mag_dims__ == 3) {
        /* make pure spinor up- and dn- wave functions */
        phi__.copy_from(device_t::CPU, num_bands__, phi__, 0, 0, 1, num_bands__);
    }
}

void test_davidson(cmd_args const& args__)
{
    auto pu        = get_device_t(args__.value<std::string>("device", "CPU"));
    auto pw_cutoff = args__.value<double>("pw_cutoff", 30);
    auto gk_cutoff = args__.value<double>("gk_cutoff", 10);
    auto N         = args__.value<int>("N", 1);
    auto mpi_grid  = args__.value<std::vector<int>>("mpi_grid", {1, 1});
    auto solver    = args__.value<std::string>("solver", "lapack");

    utils::timer t1("test_davidson|setup");

    /* create simulation context */
    Simulation_context ctx(
        "{"
        "   \"parameters\" : {"
        "        \"electronic_structure_method\" : \"pseudopotential\""
        "    },"
        "   \"control\" : {"
        "       \"verification\" : 0"
        "    }"
        "}");

    /* add a new atom type to the unit cell */
    auto& atype = ctx.unit_cell().add_atom_type("Cu");
    /* set pseudo charge */
    atype.zn(11);
    /* set radial grid */
    atype.set_radial_grid(radial_grid_t::lin_exp, 1000, 0.0, 100.0, 6);
    /* cutoff at ~1 a.u. */
    int icut = atype.radial_grid().index_of(1.0);
    double rcut = atype.radial_grid(icut);
    /* create beta radial function */
    std::vector<double> beta(icut + 1);
    std::vector<double> beta1(icut + 1);
    for (int l = 0; l <= 2; l++) {
        for (int i = 0; i <= icut; i++) {
            double x = atype.radial_grid(i);
            beta[i] = utils::confined_polynomial(x, rcut, l, l + 1, 0);
            beta1[i] = utils::confined_polynomial(x, rcut, l, l + 2, 0);
        }
        /* add radial function for l */
        atype.add_beta_radial_function(l, beta);
        atype.add_beta_radial_function(l, beta1);
    }
    /* set local part of potential */
    std::vector<double> vloc(atype.radial_grid().num_points());
    for (int i = 0; i < atype.radial_grid().num_points(); i++) {
        //double x =  atype.radial_grid(i);
        vloc[i] = 0.0; //-atype.zn() / (std::exp(-x * (x + 1)) + x);
    }
    atype.local_potential(vloc);
    /* set Dion matrix */
    int nbf = atype.num_beta_radial_functions();
    matrix<double> dion(nbf, nbf);
    dion.zero();
    //for (int i = 0; i < nbf; i++) {
    //    dion(i, i) = -10.0;
    //}
    atype.d_mtrx_ion(dion);
    /* set atomic density */
    std::vector<double> arho(atype.radial_grid().num_points());
    for (int i = 0; i < atype.radial_grid().num_points(); i++) {
        double x = atype.radial_grid(i);
        arho[i] = 2 * atype.zn() * std::exp(-x * x) * x;
    }
    atype.ps_total_charge_density(arho);

    /* lattice constant */
    double a{5};
    /* set lattice vectors */
    ctx.unit_cell().set_lattice_vectors({{a * N, 0, 0},
                                         {0, a * N, 0},
                                         {0, 0, a * N}});
    /* add atoms */
    double p = 1.0 / N;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            for (int k = 0; k < N; k++) {
                ctx.unit_cell().add_atom("Cu", {i * p, j * p, k * p});
            }
        }
    }

    /* initialize the context */
    ctx.set_verbosity(1);
    ctx.pw_cutoff(pw_cutoff);
    ctx.gk_cutoff(gk_cutoff);
    ctx.set_processing_unit(pu);
    ctx.set_mpi_grid_dims(mpi_grid);
    ctx.gen_evp_solver_name(solver);
    ctx.std_evp_solver_name(solver);

    t1.stop();

    ctx.set_verbosity(1);
    ctx.iterative_solver_tolerance(1e-12);
    //ctx.set_iterative_solver_type("exact");
    ctx.initialize();

    Density rho(ctx);
    rho.initial_density();
    rho.zero();<commit_msg>change naming convention for setter/getter functions<commit_after>#include <sirius.h>

using namespace sirius;

void init_wf(K_point* kp__, Wave_functions& phi__, int num_bands__, int num_mag_dims__)
{
    std::vector<double> tmp(0xFFFF);
    for (int i = 0; i < 0xFFFF; i++) {
        tmp[i] = utils::random<double>();
    }

    phi__.pw_coeffs(0).prime().zero();

    //#pragma omp parallel for schedule(static)
    for (int i = 0; i < num_bands__; i++) {
        for (int igk_loc = 0; igk_loc < kp__->num_gkvec_loc(); igk_loc++) {
            /* global index of G+k vector */
            int igk = kp__->idxgk(igk_loc);
            if (igk == 0) {
                phi__.pw_coeffs(0).prime(igk_loc, i) = 1.0;
            }
            if (igk == i + 1) {
                phi__.pw_coeffs(0).prime(igk_loc, i) = 0.5;
            }
            if (igk == i + 2) {
                phi__.pw_coeffs(0).prime(igk_loc, i) = 0.25;
            }
            if (igk == i + 3) {
                phi__.pw_coeffs(0).prime(igk_loc, i) = 0.125;
            }
            phi__.pw_coeffs(0).prime(igk_loc, i) += tmp[(igk + i) & 0xFFFF] * 1e-5;
        }
    }

    if (num_mag_dims__ == 3) {
        /* make pure spinor up- and dn- wave functions */
        phi__.copy_from(device_t::CPU, num_bands__, phi__, 0, 0, 1, num_bands__);
    }
}

void test_davidson(cmd_args const& args__)
{
    auto pu        = get_device_t(args__.value<std::string>("device", "CPU"));
    auto pw_cutoff = args__.value<double>("pw_cutoff", 30);
    auto gk_cutoff = args__.value<double>("gk_cutoff", 10);
    auto N         = args__.value<int>("N", 1);
    auto mpi_grid  = args__.value<std::vector<int>>("mpi_grid", {1, 1});
    auto solver    = args__.value<std::string>("solver", "lapack");

    utils::timer t1("test_davidson|setup");

    /* create simulation context */
    Simulation_context ctx(
        "{"
        "   \"parameters\" : {"
        "        \"electronic_structure_method\" : \"pseudopotential\""
        "    },"
        "   \"control\" : {"
        "       \"verification\" : 0"
        "    }"
        "}");

    /* add a new atom type to the unit cell */
    auto& atype = ctx.unit_cell().add_atom_type("Cu");
    /* set pseudo charge */
    atype.zn(11);
    /* set radial grid */
    atype.set_radial_grid(radial_grid_t::lin_exp, 1000, 0.0, 100.0, 6);
    /* cutoff at ~1 a.u. */
    int icut = atype.radial_grid().index_of(1.0);
    double rcut = atype.radial_grid(icut);
    /* create beta radial function */
    std::vector<double> beta(icut + 1);
    std::vector<double> beta1(icut + 1);
    for (int l = 0; l <= 2; l++) {
        for (int i = 0; i <= icut; i++) {
            double x = atype.radial_grid(i);
            beta[i] = utils::confined_polynomial(x, rcut, l, l + 1, 0);
            beta1[i] = utils::confined_polynomial(x, rcut, l, l + 2, 0);
        }
        /* add radial function for l */
        atype.add_beta_radial_function(l, beta);
        atype.add_beta_radial_function(l, beta1);
    }
    /* set local part of potential */
    std::vector<double> vloc(atype.radial_grid().num_points());
    for (int i = 0; i < atype.radial_grid().num_points(); i++) {
        //double x =  atype.radial_grid(i);
        vloc[i] = 0.0; //-atype.zn() / (std::exp(-x * (x + 1)) + x);
    }
    atype.local_potential(vloc);
    /* set Dion matrix */
    int nbf = atype.num_beta_radial_functions();
    matrix<double> dion(nbf, nbf);
    dion.zero();
    //for (int i = 0; i < nbf; i++) {
    //    dion(i, i) = -10.0;
    //}
    atype.d_mtrx_ion(dion);
    /* set atomic density */
    std::vector<double> arho(atype.radial_grid().num_points());
    for (int i = 0; i < atype.radial_grid().num_points(); i++) {
        double x = atype.radial_grid(i);
        arho[i] = 2 * atype.zn() * std::exp(-x * x) * x;
    }
    atype.ps_total_charge_density(arho);

    /* lattice constant */
    double a{5};
    /* set lattice vectors */
    ctx.unit_cell().set_lattice_vectors({{a * N, 0, 0},
                                         {0, a * N, 0},
                                         {0, 0, a * N}});
    /* add atoms */
    double p = 1.0 / N;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            for (int k = 0; k < N; k++) {
                ctx.unit_cell().add_atom("Cu", {i * p, j * p, k * p});
            }
        }
    }

    /* initialize the context */
    ctx.verbosity(1);
    ctx.pw_cutoff(pw_cutoff);
    ctx.gk_cutoff(gk_cutoff);
    ctx.set_processing_unit(pu);
    ctx.set_mpi_grid_dims(mpi_grid);
    ctx.gen_evp_solver_name(solver);
    ctx.std_evp_solver_name(solver);

    t1.stop();

    ctx.verbosity(1);
    ctx.iterative_solver_tolerance(1e-12);
    //ctx.set_iterative_solver_type("exact");
    ctx.initialize();

    Density rho(ctx);
    rho.initial_density();
    rho.zero();
{'fraction_non_alphanumeric': 0.1712421009012742, 'fraction_numerical': 0.02030456852791878, 'mean_word_length': 2.0502369668246447, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7529615', 'n_tokens_mistral': 3791, 'n_tokens_neox': 3583, 'n_words': 744}