Skip to content

Update DensityMatrix class #342

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jeanlucf22 merged 1 commit into from
Apr 11, 2025
Merged

Update DensityMatrix class #342

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions src/DavidsonSolver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -206,8 +206,8 @@ double DavidsonSolver<OrbitalsType, MatrixType>::evaluateDerivative(
const double dbeta = 0.0001;
*work2N_ = dm2Ninit;
work2N_->axpy(dbeta, delta_dm);
// proj_mat2N_->setDM(*work2N_,orbitals.getIterativeIndex());
proj_mat2N_->setDM(*work2N_, -1);

proj_mat2N_->setDM(*work2N_);
proj_mat2N_->computeOccupationsFromDM();

const double tsd0e = evalEntropy(proj_mat2N_.get(), false, os_);
Expand Down Expand Up @@ -250,7 +250,7 @@ void DavidsonSolver<OrbitalsType, MatrixType>::buildTarget2N_MVP(
// if( onpe0 )os_<<"Compute X2N..."<<endl;
proj_mat2N_->setHB2H();

proj_mat2N_->updateDM(orbitals_index);
proj_mat2N_->updateDM();

target = proj_mat2N_->dm();

Expand Down Expand Up @@ -584,7 +584,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
//
if (mmpi.PE0() && ct.verbose > 2)
os_ << "Target energy..." << std::endl;
proj_mat2N_->setDM(target, orbitals.getIterativeIndex());
proj_mat2N_->setDM(target);
proj_mat2N_->computeOccupationsFromDM();
double nel = proj_mat2N_->getNel();
if (mmpi.PE0() && ct.verbose > 2)
Expand Down Expand Up @@ -643,7 +643,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
// update DM
*work2N_ = dm2Ninit;
work2N_->axpy(beta, delta_dm);
proj_mat2N_->setDM(*work2N_, orbitals.getIterativeIndex());
proj_mat2N_->setDM(*work2N_);

if (inner_it < ct.dm_inner_steps - 1)
{
Expand Down Expand Up @@ -809,7 +809,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
= dynamic_cast<ProjectedMatrices<MatrixType>*>(
orbitals.getProjMatrices());
assert(pmat);
pmat->buildDM(new_occ, orbitals.getIterativeIndex());
pmat->buildDM(new_occ);

if (retval == 0) break;

Expand Down
81 changes: 39 additions & 42 deletions src/DensityMatrix.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,6 @@

const double factor_kernel4dot = 10.;

#define PROCRUSTES 0

#define MGMOL_DENSITYMATRIX_FAIL(X) \
{ \
std::cerr << "DensityMatrix failure:" << std::endl; \
Expand All @@ -36,7 +34,7 @@ const double factor_kernel4dot = 10.;
template <class MatrixType>
DensityMatrix<MatrixType>::DensityMatrix(const int ndim)
: dim_(ndim),
orbitals_index_(-1),
update_index_(-1),
occ_uptodate_(false),
uniform_occ_(false),
stripped_(false)
Expand Down Expand Up @@ -66,8 +64,8 @@ DensityMatrix<MatrixType>::~DensityMatrix()
}

template <class MatrixType>
void DensityMatrix<MatrixType>::build(const MatrixType& zmat,
const std::vector<double>& occ, const int new_orbitals_index)
void DensityMatrix<MatrixType>::build(
const MatrixType& zmat, const std::vector<double>& occ)
{
#ifdef PRINT_OPERATIONS
MGmol_MPI& mmpi = *(MGmol_MPI::instance());
Expand Down Expand Up @@ -98,32 +96,30 @@ void DensityMatrix<MatrixType>::build(const MatrixType& zmat,
work_->symm('r', 'l', 1., gamma, zmat, 0.);
kernel4dot_->gemm('n', 't', 1., *work_, zmat, 0.);

stripped_ = false;
orbitals_index_ = new_orbitals_index;
stripped_ = false;
update_index_++;
}

template <class MatrixType>
void DensityMatrix<MatrixType>::build(
const MatrixType& zmat, const int new_orbitals_index)
void DensityMatrix<MatrixType>::build(const MatrixType& zmat)
{
build(zmat, occupation_, new_orbitals_index);
build(zmat, occupation_);
}

// build diagonal matrix
template <class MatrixType>
void DensityMatrix<MatrixType>::build(
const std::vector<double>& occ, const int new_orbitals_index)
void DensityMatrix<MatrixType>::build(const std::vector<double>& occ)
{
assert(dm_ != nullptr);
assert(!occ.empty());

setOccupations(occ);

build(new_orbitals_index);
build();
}

template <class MatrixType>
void DensityMatrix<MatrixType>::build(const int new_orbitals_index)
void DensityMatrix<MatrixType>::build()
{
MGmol_MPI& mmpi = *(MGmol_MPI::instance());
#ifdef PRINT_OPERATIONS
Expand All @@ -148,13 +144,12 @@ void DensityMatrix<MatrixType>::build(const int new_orbitals_index)
* std::min(1., factor_kernel4dot * occupation_[i]));
kernel4dot_->setDiagonal(w);

stripped_ = false;
orbitals_index_ = new_orbitals_index;
stripped_ = false;
update_index_++;
}

template <class MatrixType>
void DensityMatrix<MatrixType>::setUniform(
const double nel, const int new_orbitals_index)
void DensityMatrix<MatrixType>::setUniform(const double nel)
{
assert(!occupation_.empty());

Expand All @@ -171,22 +166,13 @@ void DensityMatrix<MatrixType>::setUniform(

uniform_occ_ = true;

build(occupation_, new_orbitals_index);
}

template <class MatrixType>
void DensityMatrix<MatrixType>::buildFromBlock(const MatrixType& block00)
{
dm_->clear();
dm_->assign(block00, 0, 0);
dm_->print(std::cout, 0, 0, 25, 25);
build(occupation_);
}

template <class MatrixType>
void DensityMatrix<MatrixType>::rotate(
const MatrixType& rotation_matrix, const bool flag_eigen)
{

if (!flag_eigen)
{
MatrixType invU(rotation_matrix);
Expand All @@ -204,6 +190,8 @@ void DensityMatrix<MatrixType>::rotate(
invU.getrs('n', tmp, ipiv);

*dm_ = tmp;

update_index_++;
}
}

Expand Down Expand Up @@ -370,8 +358,7 @@ double DensityMatrix<MatrixType>::computeEntropy() const
}

template <class MatrixType>
void DensityMatrix<MatrixType>::setto2InvS(
const MatrixType& invS, const int orbitals_index)
void DensityMatrix<MatrixType>::setto2InvS(const MatrixType& invS)
{
*dm_ = invS;
dm_->scal(orbital_occupation_);
Expand All @@ -382,8 +369,8 @@ void DensityMatrix<MatrixType>::setto2InvS(
occupation_[st] = 1.;
occ_uptodate_ = true;
}
uniform_occ_ = false;
orbitals_index_ = orbitals_index;
uniform_occ_ = false;
update_index_++;
}

template <class MatrixType>
Expand All @@ -400,8 +387,7 @@ void DensityMatrix<MatrixType>::stripS(const MatrixType& ls)
}

template <class MatrixType>
void DensityMatrix<MatrixType>::dressUpS(
const MatrixType& ls, const int new_orbitals_index)
void DensityMatrix<MatrixType>::dressUpS(const MatrixType& ls)
{
assert(stripped_);

Expand All @@ -410,10 +396,11 @@ void DensityMatrix<MatrixType>::dressUpS(
*dm_ = *work_;
ls.trtrs('l', 't', 'n', *dm_);

orbitals_index_ = new_orbitals_index;
occ_uptodate_ = false;
uniform_occ_ = false;
stripped_ = false;
update_index_++;

occ_uptodate_ = false;
uniform_occ_ = false;
stripped_ = false;
}

// dm_ -> u*dm_*u^T
Expand All @@ -424,6 +411,8 @@ void DensityMatrix<MatrixType>::transform(const MatrixType& u)
{
work_->gemm('n', 't', 1., *dm_, u, 0.);
dm_->gemm('n', 'n', 1., u, *work_, 0.);

update_index_++;
}

template <class MatrixType>
Expand All @@ -434,13 +423,21 @@ double DensityMatrix<MatrixType>::getExpectation(const MatrixType& A)
}

template <class MatrixType>
void DensityMatrix<MatrixType>::mix(
const double mix, const MatrixType& matA, const int new_orbitals_index)
void DensityMatrix<MatrixType>::mix(const double mix, const MatrixType& matA)
{
dm_->scal(1. - mix);

dm_->axpy(mix, matA);
orbitals_index_ = new_orbitals_index;

update_index_++;
}

template <class MatrixType>
void DensityMatrix<MatrixType>::linearExtrapolate(const MatrixType& previous_dm)
{
dm_->scal(2.);
dm_->axpy(-1., previous_dm);

update_index_++;
}

template <class MatrixType>
Expand Down
40 changes: 22 additions & 18 deletions src/DensityMatrix.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,10 @@ class DensityMatrix
MatrixType* kernel4dot_;
MatrixType* work_;

int orbitals_index_;
/*!
* Keep track of changes, incremented every time dm_ is updated
*/
int update_index_;

bool occ_uptodate_;
bool uniform_occ_;
Expand All @@ -46,16 +49,16 @@ class DensityMatrix
DensityMatrix& operator=(const DensityMatrix&);
DensityMatrix(const DensityMatrix&);

void build(const int new_orbitals_index);
void build();

public:
DensityMatrix(const int ndim);

~DensityMatrix();

void setUniform(const double nel, const int new_orbitals_index);
void setUniform(const double nel);

int getOrbitalsIndex() const { return orbitals_index_; }
int getIndex() const { return update_index_; }

bool occupationsUptodate() const { return occ_uptodate_; }
bool fromUniformOccupations() const { return uniform_occ_; }
Expand Down Expand Up @@ -84,10 +87,10 @@ class DensityMatrix

const MatrixType& kernel4dot() const { return *kernel4dot_; }

void setMatrix(const MatrixType& mat, const int orbitals_index)
void setMatrix(const MatrixType& mat)
{
*dm_ = mat;
orbitals_index_ = orbitals_index;
*dm_ = mat;
update_index_++;

setDummyOcc();

Expand All @@ -112,7 +115,7 @@ class DensityMatrix
uniform_occ_ = false;
stripped_ = false;

orbitals_index_ = 0;
update_index_ = 0;
}

void getOccupations(std::vector<double>& occ) const
Expand All @@ -126,31 +129,32 @@ class DensityMatrix

void setOccupations(const std::vector<double>& occ);

void setto2InvS(const MatrixType& invS, const int orbitals_index);
void setto2InvS(const MatrixType& invS);

void stripS(const MatrixType& ls);
void dressUpS(const MatrixType& ls, const int new_orbitals_index);
void dressUpS(const MatrixType& ls);

// dm_ -> u*dm_*u^T
void transform(const MatrixType& u);

void buildFromBlock(const MatrixType& block00);

double computeEntropy() const;
void computeOccupations(const MatrixType& ls);
void build(const std::vector<double>& occ, const int new_orbitals_index);
void build(const MatrixType& z, const int new_orbitals_index);
void build(const MatrixType& z, const std::vector<double>& occ,
const int new_orbitals_index);
void build(const std::vector<double>& occ);
void build(const MatrixType& z);
void build(const MatrixType& z, const std::vector<double>& occ);

void rotate(const MatrixType& rotation_matrix, const bool flag_eigen);
void printOccupations(std::ostream& os) const;
void diagonalize(const MatrixType& ls, std::vector<double>& occ);
void diagonalize(
const char eigv, std::vector<double>& occ, MatrixType& vect);
double getExpectation(const MatrixType& A);
void mix(
const double mix, const MatrixType& matA, const int new_orbitals_index);
void mix(const double mix, const MatrixType& matA);

/*!
* dm <- dm + (dm-previous_dm) = 2.*dm - previous_dm
*/
void linearExtrapolate(const MatrixType& previous_dm);

int write(HDFrestart& h5f_file, std::string& name);
int read(HDFrestart& h5f_file, std::string& name);
Expand Down
15 changes: 6 additions & 9 deletions src/DensityMatrixSparse.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ DensityMatrixSparse::DensityMatrixSparse(
MGmol_MPI& mmpi = *(MGmol_MPI::instance());
orbital_occupation_ = mmpi.nspin() > 1 ? 1. : 2.;

orbitals_index_ = -1;
update_index_ = -1;

if (dim_ > 0)
{
Expand All @@ -44,27 +44,24 @@ DensityMatrixSparse::~DensityMatrixSparse()
}
}

void DensityMatrixSparse::setUniform(const double nel, const int orbitals_index)
void DensityMatrixSparse::setUniform(const double nel)
{
const double occ = (double)((double)nel / (double)dim_);
assert(occ < 1.01);
orbitals_index_ = orbitals_index;
update_index_++;
(*dm_).reset();
const double uval = (double)occ * orbital_occupation_;
for (std::vector<int>::const_iterator st = locvars_.begin();
st != locvars_.end(); ++st)
(*dm_).insertMatrixElement(*st, *st, uval, INSERT, true);

return;
}

void DensityMatrixSparse::setto2InvS(
const VariableSizeMatrix<sparserow>& invS, const int orbitals_index)
void DensityMatrixSparse::setto2InvS(const VariableSizeMatrix<sparserow>& invS)
{
*dm_ = invS;
dm_->scale(orbital_occupation_);

orbitals_index_ = orbitals_index;
update_index_++;
}
// build density matrix, given computed locally centered data
void DensityMatrixSparse::assembleMatrixFromCenteredData(
Expand Down Expand Up @@ -92,7 +89,7 @@ void DensityMatrixSparse::assembleMatrixFromCenteredData(
dtor_DM.updateLocalRows((*dm_));
gather_DM_tm_.stop();

orbitals_index_ = orbitals_index;
update_index_ = orbitals_index;
}
// compute trace of dot product dm_ . vsmat
double DensityMatrixSparse::getTraceDotProductWithMat(
Expand Down
Loading