abacus-develop/source/source_esolver/esolver_of.h at 1a6a3bc57909c5e9c43ec2a4507646a2f7f8127f · deepmodeling/abacus-develop · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
#ifndef ESOLVER_OF_H
#define ESOLVER_OF_H

#include "esolver_fp.h"
#include "source_base/opt_DCsrch.h"
#include "source_base/opt_TN.hpp"
#include "source_pw/module_ofdft/kedf_manager.h"
#include "source_psi/psi.h"

namespace ModuleESolver
{
class ESolver_OF : public ESolver_FP
{
  public:
    ESolver_OF();
    ~ESolver_OF();

    virtual void before_all_runners(UnitCell& ucell, const Input_para& inp) override;

    virtual void runner(UnitCell& ucell, const int istep) override;

    virtual void after_all_runners(UnitCell& ucell) override;

    virtual double cal_energy() override;

    virtual void cal_force(UnitCell& ucell, ModuleBase::matrix& force) override;

    virtual void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) override;

  protected:
    // ======================= variables ==========================
    // ---------- the kinetic energy density functionals ----------
    KEDF_Manager* kedf_manager_ = nullptr; // KEDF manager, which will be initialized in before_all_runners

    // ----------------- the optimization methods ------------------
    ModuleBase::Opt_CG* opt_cg_ = nullptr;
    ModuleBase::Opt_TN* opt_tn_ = nullptr;
    ModuleBase::Opt_DCsrch* opt_dcsrch_ = nullptr;
    ModuleBase::Opt_CG* opt_cg_mag_ = nullptr; // for spin2 case, under testing

    // ----------------- necessary parameters from INPUT ------------
    std::string of_kinetic_ = "wt";  // Kinetic energy functional, such as TF, VW, WT
    std::string of_method_ = "tn";   // optimization method, include cg1, cg2, tn (default), bfgs
    std::string of_conv_ = "energy"; // select the convergence criterion, potential, energy (default), or both
    double of_tole_ = 2e-6; // tolerance of the energy change (in Ry) for determining the convergence, default=2e-6 Ry
    double of_tolp_ = 1e-5; // tolerance of potential for determining the convergence, default=1e-5 in a.u.
    int max_iter_ = 50;     // scf_nmax

    // ------------------ parameters from other module --------------
    double dV_ = 0;           // volume of one grid point in real space
    double* nelec_ = nullptr; // number of electrons with each spin

    // ----- parameters and arrays used in density optimization -----
    int iter_ = 0;                                // iteration number
    double** pdirect_ = nullptr;                  // optimization direction of phi, which is sqrt(rho)
    std::complex<double>** precip_dir_ = nullptr; // direction in reciprocal space, used when of_full_pw=false.
    double* theta_ = nullptr;                     // step length
    double** pdEdphi_ = nullptr;                  // dE/dphi
    double** pdLdphi_ = nullptr;                  // dL/dphi
    double** pphi_ = nullptr;                     // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi().
    char* task_ = nullptr;                        // used in line search
    int tn_spin_flag_ = -1;                       // spin flag used in cal_potential, which will be called by opt_tn
    int max_dcsrch_ = 200;                        // max no. of line search
    int flag_ = -1;                               // flag of TN
    Charge<double>* ptemp_rho_ = nullptr;                 // used in line search
    psi::Psi<double>* psi_ = nullptr;             // sqrt(rho)

    // ----------------- used for convergence check -------------------
    double energy_llast_ = 0;
    double energy_last_ = 0;
    double energy_current_ = 0;
    double normdLdphi_llast_ = 100;
    double normdLdphi_last_ = 100;
    double normdLdphi_ = 100.;

    // ==================== main process of OFDFT ======================
    void before_opt(const int istep, UnitCell& ucell);
    void update_potential(UnitCell& ucell);
    void optimize(UnitCell& ucell);
    void update_rho();
    bool check_exit(bool& conv_esolver);
    void after_opt(const int istep, UnitCell& ucell, const bool conv_esolver);

    // ============================ tools ===============================
    // --------------------- initialize ---------------------------------
    void init_elecstate(UnitCell& ucell);
    void allocate_array();

    // --------------------- calculate physical qualities ---------------
    std::function<void(double*, double*)> bound_cal_potential_;
    void cal_potential_wrapper(double* ptemp_phi, double* rdLdphi);
    void cal_potential(double* ptemp_phi, double* rdLdphi, UnitCell& ucell);
    void cal_dEdtheta(double** ptemp_phi, Charge<double>* temp_rho, UnitCell& ucell, double* ptheta, double* rdEdtheta);
    double cal_mu(double* pphi, double* pdEdphi, double nelec);

    // --------------------- determine the optimization direction -------
    void adjust_direction();
    void check_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell);
    void test_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell);

    // --------------------- interface to blas --------------------------
    double inner_product(double* pa, double* pb, int length, double dV = 1)
    {
        double innerproduct = BlasConnector::dot(length, pa, 1, pb, 1);
        innerproduct *= dV;
        return innerproduct;
    }

    // ---------------------- interfaces to optimization methods --------
    void init_opt();
    void get_direction(UnitCell& ucell);
    void get_step_length(double* dEdtheta, double** ptemp_phi, UnitCell& ucell);
};
} // namespace ModuleESolver

#endif