spin_constrain.h

Go to the documentation of this file.

#ifndef SPIN_CONSTRAIN_H
#define SPIN_CONSTRAIN_H
 
#include <map>
#include <vector>
 
#include "source_base/constants.h"
#include "source_base/tool_quit.h"
#include "source_base/tool_title.h"
#include "source_base/vector3.h"
#include "source_basis/module_ao/parallel_orbitals.h"
#include "source_cell/klist.h"
#include "source_cell/unitcell.h"
#include "source_hamilt/operator.h"
#include "source_estate/elecstate.h"
 
#ifdef __LCAO
#include "source_estate/module_dm/density_matrix.h" // mohan add 2025-11-02
#endif
 
namespace spinconstrain
{
 
struct ScAtomData;
 
template <typename TK>
class SpinConstrain
{
public:
  void init_sc(double sc_thr_in,
               int nsc_in,
               int nsc_min_in,
               double alpha_trial_in,
               double sccut_in,
               double sc_drop_thr_in,
               const UnitCell& ucell,
               Parallel_Orbitals* ParaV_in,
               int nspin_in,
               const K_Vectors& kv_in,
               void* p_hamilt_in,
               void* psi_in,
#ifdef __LCAO
               elecstate::DensityMatrix<TK, double> *dm_in, // mohan add 2025-11-02
#endif
               elecstate::ElecState* pelec_in,
               ModulePW::PW_Basis_K* pw_wfc_in = nullptr);
 
  void cal_mi_lcao(const int& step, bool print = false);
 
  void cal_mi_pw();
 
  void cal_mw_from_lambda(int i_step, 
          const ModuleBase::Vector3<double>* delta_lambda = nullptr);
 
  double cal_escon();
 
  double get_escon();
 
  void run_lambda_loop(int outer_step, 
          bool rerun = true);
 
  void update_psi_charge(const ModuleBase::Vector3<double>* delta_lambda, bool pw_solve = true);
 
  void calculate_delta_hcc(std::complex<double>* h_tmp, 
          const std::complex<double>* becp_k, 
          const ModuleBase::Vector3<double>* delta_lambda, 
          const int nbands, const int nkb, const int* nh_iat);
 
  bool check_rms_stop(int outer_step, int i_step, double rms_error, double duration, double total_duration);
 
  void check_restriction(const std::vector<ModuleBase::Vector3<double>>& search, double& alpha_trial);
 
  bool check_gradient_decay(std::vector<ModuleBase::Vector3<double>> new_spin,
                            std::vector<ModuleBase::Vector3<double>> old_spin,
                            std::vector<ModuleBase::Vector3<double>> new_delta_lambda,
                            std::vector<ModuleBase::Vector3<double>> old_delta_lambda,
                            bool print = false);
  double cal_alpha_opt(std::vector<ModuleBase::Vector3<double>> spin,
                       std::vector<ModuleBase::Vector3<double>> spin_plus,
                       const double alpha_trial);
  void print_header();
  void print_termination();
 
  void print_Mi(std::ofstream& ofs_running);
 
  void print_Mag_Force(std::ofstream& ofs_running);
 
  bool higher_mag_prec = false;
 
public:
    Parallel_Orbitals *ParaV = nullptr;
    //--------------------------------------------------------------------------------
    // pointers for solve Hamiltonian to get new Magnetization from Lambda
    void* p_hamilt = nullptr;
    void* psi = nullptr;
    elecstate::ElecState* pelec = nullptr;
    ModulePW::PW_Basis_K* pw_wfc_ = nullptr;
#ifdef __LCAO
    elecstate::DensityMatrix<TK, double>* dm_;
#endif
    double tpiba = 0.0; 
    const double meV_to_Ry = 7.349864435130999e-05;
    K_Vectors kv_;
    //--------------------------------------------------------------------------------
 
  public:
    static SpinConstrain& getScInstance();
    SpinConstrain(SpinConstrain const&) = delete;
    SpinConstrain(SpinConstrain&&) = delete;
    void set_atomCounts(const std::map<int, int>& atomCounts_in);
    const std::map<int, int>& get_atomCounts() const;
    void set_orbitalCounts(const std::map<int, int>& orbitalCounts_in);
    const std::map<int, int>& get_orbitalCounts() const;
    void set_lnchiCounts(const std::map<int, std::map<int, int>>& lnchiCounts_in);
    const std::map<int, std::map<int, int>>& get_lnchiCounts() const;
    void set_sc_lambda();
    void set_sc_lambda(const ModuleBase::Vector3<double>* lambda_in, int nat_in);
    void set_target_mag();
    void set_target_mag(const ModuleBase::Vector3<double>* target_mag_in, int nat_in);
    void set_target_mag(const std::vector<ModuleBase::Vector3<double>>& target_mag_in);
    void set_constrain();
    void set_constrain(const ModuleBase::Vector3<int>* constrain_in, int nat_in);
    const std::vector<ModuleBase::Vector3<double>>& get_sc_lambda() const;
    const std::vector<ModuleBase::Vector3<double>>& get_target_mag() const;
    const std::vector<ModuleBase::Vector3<int>>& get_constrain() const;
    int get_nat();
    int get_ntype();
    void check_atomCounts();
    int get_iat(int itype, int atom_index);
    void set_nspin(int nspin);
    int get_nspin();
    void zero_Mi();
    double get_decay_grad(int itype);
    void set_decay_grad();
    const std::vector<double>& get_decay_grad();
    void set_decay_grad(const double* decay_grad_in, int ntype_in);
    void set_sc_drop_thr(double sc_drop_thr_in);
    void set_input_parameters(double sc_thr_in,
                              int nsc_in,
                              int nsc_min_in,
                              double alpha_trial_in,
                              double sccut_in,
                              double sc_drop_thr_in);
    double get_sc_thr();
    int get_nsc();
    int get_nsc_min();
    double get_alpha_trial();
    double get_sccut();
    double get_sc_drop_thr();
    void set_ParaV(Parallel_Orbitals* ParaV_in);
    void set_solver_parameters(const K_Vectors& kv_in,
                               void* p_hamilt_in,
                               void* psi_in,
                               elecstate::ElecState* pelec_in);
 
  private:
    SpinConstrain(){};                               // Private constructor
    ~SpinConstrain(){};                              // Destructor
    SpinConstrain& operator=(SpinConstrain const&) = delete;  // Copy assign
    SpinConstrain& operator=(SpinConstrain &&) = delete;      // Move assign
    std::map<int, std::vector<ScAtomData>> ScData;
    std::map<int, double> ScDecayGrad; // in unit of uB^2/eV
    std::vector<double> decay_grad_;   // in unit of uB^2/Ry
    std::map<int, int> atomCounts;
    std::map<int, int> orbitalCounts;
    std::map<int, std::map<int, int>> lnchiCounts;
    std::vector<ModuleBase::Vector3<double>> lambda_; // in unit of Ry/uB in code, but in unit of meV/uB in input file
    std::vector<ModuleBase::Vector3<double>> target_mag_; // in unit of uB
    std::vector<ModuleBase::Vector3<double>> Mi_; // in unit of uB
    std::vector<std::string> atomLabels_;
    double escon_ = 0.0;
    int nspin_ = 0;
    int npol_ = 1;
    int nsc_;
    int nsc_min_;
    double sc_drop_thr_ = 1e-3;
    double sc_thr_; // in unit of uB
    double current_sc_thr_;
    std::vector<ModuleBase::Vector3<int>> constrain_;
    bool debug = false;
    double alpha_trial_; // in unit of Ry/uB^2 = 0.01 eV/uB^2
    double restrict_current_; // in unit of Ry/uB = 3 eV/uB
 
  public:
    void set_operator(hamilt::Operator<TK>* op_in);
    void set_mag_converged(bool is_Mi_converged_in){this->is_Mi_converged = is_Mi_converged_in;}
    bool mag_converged() const {return this->is_Mi_converged;}
  private:
    hamilt::Operator<TK>* p_operator = nullptr;
    bool is_Mi_converged = false;
 
    TK* sub_h_save;
    TK* sub_s_save;
    TK* becp_save;
};
 
 
struct ScAtomData {
    int index;
    std::vector<double> lambda;
    std::vector<double> target_mag;
    std::vector<int> constrain;
    int mag_type;
    double target_mag_val;
    double target_mag_angle1;
    double target_mag_angle2;
};
 
} // namespace spinconstrain
 
#endif // SPIN_CONSTRAIN_H