Skip to content

Commit

Permalink
add to hK passed in rather than LM.Hloc
Browse files Browse the repository at this point in the history
  • Loading branch information
@maki49
maki49 committed Dec 9, 2023
1 parent 69b0b93 commit c4eebcb
Show file tree
Hide file tree
Showing 5 changed files with 314 additions and 333 deletions.
70 changes: 26 additions & 44 deletions source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,39 +139,41 @@ void LCAO_Matrix::allocate_HS_R(const int &nnR)
// set 'dtype' matrix element (iw1_all, iw2_all) with
// an input value 'v'
//------------------------------------------------------
void LCAO_Matrix::set_HSgamma(
const int &iw1_all, // index i for atomic orbital (row)
const int &iw2_all, // index j for atomic orbital (column)
const double& v, // value for matrix element (i,j)
double* HSloc) //input pointer for store the matrix
template<typename T>
void LCAO_Matrix::set_mat2d(
const int& global_ir, // index i for atomic orbital (row)
const int& global_ic, // index j for atomic orbital (column)
const T& v, // value for matrix element (i,j)
const Parallel_Orbitals& pv,
T* HSloc) //input pointer for store the matrix
{
// use iw1_all and iw2_all to set Hloc
// becareful! The ir and ic may be < 0 !!!
const int ir = this->ParaV->global2local_row(iw1_all);
const int ic = this->ParaV->global2local_col(iw2_all);
const int ir = pv.global2local_row(global_ir);
const int ic = pv.global2local_col(global_ic);

//const int index = ir * ParaO.ncol + ic;
long index=0;

// save the matrix as column major format
if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER())
{
index=ic*this->ParaV->nrow+ir;
index = ic * pv.nrow + ir;
}
else
{
index=ir*this->ParaV->ncol+ic;
index = ir * pv.ncol + ic;
}

if( index >= this->ParaV->nloc)
if (index >= pv.nloc)
{
std::cout << " iw1_all = " << iw1_all << std::endl;
std::cout << " iw2_all = " << iw2_all << std::endl;
std::cout << " iw1_all = " << global_ir << std::endl;
std::cout << " iw2_all = " << global_ic << std::endl;
std::cout << " ir = " << ir << std::endl;
std::cout << " ic = " << ic << std::endl;
std::cout << " index = " << index << std::endl;
std::cout << " this->ParaV->nloc = " << this->ParaV->nloc << std::endl;
ModuleBase::WARNING_QUIT("LCAO_Matrix","set_HSgamma");
std::cout << " ParaV->nloc = " << pv.nloc << std::endl;
ModuleBase::WARNING_QUIT("LCAO_Matrix", "set_mat2d");
}

//using input pointer HSloc
Expand All @@ -180,41 +182,21 @@ void LCAO_Matrix::set_HSgamma(
return;
}

void LCAO_Matrix::set_HSk(const int &iw1_all, const int &iw2_all, const std::complex<double> &v, const char &dtype, const int spin)
void LCAO_Matrix::set_HSgamma(const int& iw1_all, const int& iw2_all, const double& v, double* HSloc)
{
LCAO_Matrix::set_mat2d<double>(iw1_all, iw2_all, v, *this->ParaV, HSloc);
return;
}
void LCAO_Matrix::set_HSk(const int& iw1_all, const int& iw2_all, const std::complex<double>& v, const char& dtype, const int spin)
{
// use iw1_all and iw2_all to set Hloc
// becareful! The ir and ic may < 0!!!!!!!!!!!!!!!!
const int ir = this->ParaV->global2local_row(iw1_all);
const int ic = this->ParaV->global2local_col(iw2_all);
//const int index = ir * this->ParaV->ncol + ic;
long index;
if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER())
{
index=ic*this->ParaV->nrow+ir;
}
else
{
index=ir*this->ParaV->ncol+ic;
}
assert(index < this->ParaV->nloc);
if (dtype=='S')//overlap Hamiltonian.
{
this->Sloc2[index] += v;
}
LCAO_Matrix::set_mat2d<std::complex<double>>(iw1_all, iw2_all, v, *this->ParaV, this->Sloc2.data());
else if (dtype=='T' || dtype=='N')// kinetic and nonlocal Hamiltonian.
{
this->Hloc_fixed2[index] += v; // because kinetic and nonlocal Hamiltonian matrices are already block-cycle staraged after caculated in lcao_nnr.cpp
// this statement will not be used.
}
LCAO_Matrix::set_mat2d<std::complex<double>>(iw1_all, iw2_all, v, *this->ParaV, this->Hloc_fixed2.data());
else if (dtype=='L') // Local potential Hamiltonian.
{
this->Hloc2[index] += v;
}
LCAO_Matrix::set_mat2d<std::complex<double>>(iw1_all, iw2_all, v, *this->ParaV, this->Hloc2.data());
else
{
ModuleBase::WARNING_QUIT("LCAO_Matrix","set_HSk");
}

ModuleBase::WARNING_QUIT("LCAO_Matrix", "set_HSk");
return;
}

Expand Down
3 changes: 2 additions & 1 deletion source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -189,7 +189,8 @@ class LCAO_Matrix
double* DHloc_fixed_23;
double* DHloc_fixed_33;


template <typename T>
static void set_mat2d(const int& global_ir, const int& global_ic, const T& v, const Parallel_Orbitals& pv, T* mat);
void set_HSgamma(const int& iw1_all, const int& iw2_all, const double& v, double* HSloc);
void set_HSk(const int &iw1_all, const int &iw2_all, const std::complex<double> &v, const char &dtype, const int spin = 0);

Expand Down
22 changes: 14 additions & 8 deletions source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,15 +33,17 @@ void OperatorEXX<OperatorLCAO<double, double>>::contributeHk(int ik)
kv,
ik,
GlobalC::exx_info.info_global.hybrid_alpha,
*this->LM->Hexxd,
*this->LM);
*this->LM->Hexxd,
*this->LM->ParaV,
*this->hK);
else
RI_2D_Comm::add_Hexx(
kv,
ik,
GlobalC::exx_info.info_global.hybrid_alpha,
*this->LM->Hexxc,
*this->LM);
*this->LM->Hexxc,
*this->LM->ParaV,
*this->hK);
}
}

Expand All @@ -57,14 +59,16 @@ void OperatorEXX<OperatorLCAO<std::complex<double>, double>>::contributeHk(int i
ik,
GlobalC::exx_info.info_global.hybrid_alpha,
*this->LM->Hexxd,
*this->LM);
*this->LM->ParaV,
*this->hK);
else
RI_2D_Comm::add_Hexx(
kv,
ik,
GlobalC::exx_info.info_global.hybrid_alpha,
*this->LM->Hexxc,
*this->LM);
*this->LM->ParaV,
*this->hK);
}
}

Expand All @@ -80,14 +84,16 @@ void OperatorEXX<OperatorLCAO<std::complex<double>, std::complex<double>>>::cont
ik,
GlobalC::exx_info.info_global.hybrid_alpha,
*this->LM->Hexxd,
*this->LM);
*this->LM->ParaV,
*this->hK);
else
RI_2D_Comm::add_Hexx(
kv,
ik,
GlobalC::exx_info.info_global.hybrid_alpha,
*this->LM->Hexxc,
*this->LM);
*this->LM->ParaV,
*this->hK);
}
}

Expand Down
141 changes: 71 additions & 70 deletions source/module_ri/RI_2D_Comm.h
View file
Original file line number Diff line number Diff line change
@@ -1,71 +1,72 @@
//=======================
// AUTHOR : Peize Lin
// DATE : 2022-08-17
//=======================

#ifndef RI_2D_COMM_H
#define RI_2D_COMM_H

#include "module_basis/module_ao/parallel_orbitals.h"
#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h"
#include "module_cell/klist.h"

#include <RI/global/Tensor.h>

#include <array>
#include <vector>
#include <tuple>
#include <map>
#include <set>
#include <deque>

namespace RI_2D_Comm
{
using TA = int;
using Tcell = int;
static const size_t Ndim = 3;
using TC = std::array<Tcell,Ndim>;
using TAC = std::pair<TA,TC>;

//public:
template<typename Tdata, typename Tmatrix>
extern std::vector<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>>
split_m2D_ktoR(const K_Vectors &kv, const std::vector<const Tmatrix*> &mks_2D, const Parallel_Orbitals &pv);

// judge[is] = {s0, s1}
extern std::vector<std::tuple<std::set<TA>, std::set<TA>>>
get_2D_judge(const Parallel_Orbitals &pv);

template<typename Tdata>
extern void add_Hexx(
const K_Vectors &kv,
const int ik,
const double alpha,
const std::vector<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>> &Hs,
LCAO_Matrix &lm);

template<typename Tdata>
extern std::vector<std::vector<Tdata>> Hexxs_to_Hk(
const K_Vectors &kv,
const Parallel_Orbitals &pv,
const std::vector< std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>> &Hexxs,
const int ik);
template<typename Tdata>
std::vector<std::vector<Tdata>> pulay_mixing(
const Parallel_Orbitals &pv,
std::deque<std::vector<std::vector<Tdata>>> &Hk_seq,
const std::vector<std::vector<Tdata>> &Hk_new,
const double mixing_beta,
const std::string mixing_mode);

//private:
extern std::vector<int> get_ik_list(const K_Vectors &kv, const int is_k);
extern inline std::tuple<int,int,int> get_iat_iw_is_block(const int iwt);
extern inline int get_is_block(const int is_k, const int is_row_b, const int is_col_b);
extern inline std::tuple<int,int> split_is_block(const int is_b);
extern inline int get_iwt(const int iat, const int iw_b, const int is_b);
}

#include "RI_2D_Comm.hpp"

//=======================
// AUTHOR : Peize Lin
// DATE : 2022-08-17
//=======================

#ifndef RI_2D_COMM_H
#define RI_2D_COMM_H

#include "module_basis/module_ao/parallel_orbitals.h"
#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h"
#include "module_cell/klist.h"

#include <RI/global/Tensor.h>

#include <array>
#include <vector>
#include <tuple>
#include <map>
#include <set>
#include <deque>

namespace RI_2D_Comm
{
using TA = int;
using Tcell = int;
static const size_t Ndim = 3;
using TC = std::array<Tcell,Ndim>;
using TAC = std::pair<TA,TC>;

//public:
template<typename Tdata, typename Tmatrix>
extern std::vector<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>>
split_m2D_ktoR(const K_Vectors &kv, const std::vector<const Tmatrix*> &mks_2D, const Parallel_Orbitals &pv);

// judge[is] = {s0, s1}
extern std::vector<std::tuple<std::set<TA>, std::set<TA>>>
get_2D_judge(const Parallel_Orbitals &pv);

template<typename Tdata, typename TK>
extern void add_Hexx(
const K_Vectors& kv,
const int ik,
const double alpha,
const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>>& Hs,
const Parallel_Orbitals& pv,
std::vector<TK>& Hloc);

template<typename Tdata>
extern std::vector<std::vector<Tdata>> Hexxs_to_Hk(
const K_Vectors &kv,
const Parallel_Orbitals &pv,
const std::vector< std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>> &Hexxs,
const int ik);
template<typename Tdata>
std::vector<std::vector<Tdata>> pulay_mixing(
const Parallel_Orbitals &pv,
std::deque<std::vector<std::vector<Tdata>>> &Hk_seq,
const std::vector<std::vector<Tdata>> &Hk_new,
const double mixing_beta,
const std::string mixing_mode);

//private:
extern std::vector<int> get_ik_list(const K_Vectors &kv, const int is_k);
extern inline std::tuple<int,int,int> get_iat_iw_is_block(const int iwt);
extern inline int get_is_block(const int is_k, const int is_row_b, const int is_col_b);
extern inline std::tuple<int,int> split_is_block(const int is_b);
extern inline int get_iwt(const int iat, const int iw_b, const int is_b);
}

#include "RI_2D_Comm.hpp"

#endif
Loading

0 comments on commit c4eebcb

Please sign in to comment.