Skip to content

Introduce new constructor for class Ions #330

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
Mar 27, 2025
Merged

Introduce new constructor for class Ions #330

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
6 changes: 3 additions & 3 deletions src/ABPG.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ void ABPG<T>::setup(T& orbitals)
//
// orthof=true: wants orthonormalized updated wave functions
template <class T>
int ABPG<T>::updateWF(T& orbitals, Ions& /*ions*/, const double precond_factor,
int ABPG<T>::updateWF(T& orbitals, Ions& ions, const double precond_factor,
const bool /*orthof*/, T& work_orbitals, const bool accelerate,
const bool print_res, const double atol)
{
Expand All @@ -51,8 +51,8 @@ int ABPG<T>::updateWF(T& orbitals, Ions& /*ions*/, const double precond_factor,
T res("Residual", orbitals, false);

const bool check_res = (atol > 0.);
double normRes = mgmol_strategy_->computeResidual(
orbitals, work_orbitals, res, (print_res || check_res), check_res);
double normRes = mgmol_strategy_->computeResidual(orbitals, work_orbitals,
ions, res, (print_res || check_res), check_res);
if (normRes < atol && check_res)
{
abpg_nl_update_tm_.stop();
Expand Down
4 changes: 2 additions & 2 deletions src/DFTsolver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,8 +196,8 @@ double DFTsolver<OrbitalsType>::evaluateEnergy(

// Get the new total energy
const double ts = 0.5 * proj_matrices_->computeEntropy(); // in [Ha]
eks_history_[0]
= energy_->evaluateTotal(ts, proj_matrices_, orbitals, print_flag, os_);
eks_history_[0] = energy_->evaluateTotal(
ts, proj_matrices_, ions_, orbitals, print_flag, os_);

sum_eig_[1] = sum_eig_[0];
sum_eig_[0] = 2. * proj_matrices_->getEigSum(); // 2.*sum in [Ry]
Expand Down
8 changes: 4 additions & 4 deletions src/DavidsonSolver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -500,7 +500,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(

ts0 = evalEntropy(projmatrices, (ct.verbose > 1), os_);
e0 = energy_->evaluateTotal(
ts0, projmatrices, orbitals, printE, os_);
ts0, projmatrices, ions_, orbitals, printE, os_);

retval = checkConvergence(e0, outer_it, ct.conv_tol);
if (retval == 0 || (outer_it == ct.max_electronic_steps))
Expand Down Expand Up @@ -549,7 +549,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(

ts0 = evalEntropy(proj_mat2N_.get(), (ct.verbose > 1), os_);
e0 = energy_->evaluateTotal(
ts0, proj_mat2N_.get(), orbitals, printE, os_);
ts0, proj_mat2N_.get(), ions_, orbitals, printE, os_);
}

// 2N x 2N target...
Expand Down Expand Up @@ -621,8 +621,8 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(

const double ts1
= evalEntropy(proj_mat2N_.get(), (ct.verbose > 2), os_);
const double e1 = energy_->evaluateTotal(
ts1, proj_mat2N_.get(), orbitals, ct.verbose - 1, os_);
const double e1 = energy_->evaluateTotal(ts1, proj_mat2N_.get(),
ions_, orbitals, ct.verbose - 1, os_);

// line minimization
beta = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_);
Expand Down
29 changes: 14 additions & 15 deletions src/Energy.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,11 +25,10 @@
#define RY2HA 0.5

template <class T>
Energy<T>::Energy(const pb::Grid& mygrid, const Ions& ions,
const Potentials& pot, const Electrostatic& es, const Rho<T>& rho,
const XConGrid& xc, SpreadPenaltyInterface<T>* spread_penalty)
Energy<T>::Energy(const pb::Grid& mygrid, const Potentials& pot,
const Electrostatic& es, const Rho<T>& rho, const XConGrid& xc,
SpreadPenaltyInterface<T>* spread_penalty)
: mygrid_(mygrid),
ions_(ions),
pot_(pot),
es_(es),
rho_(rho),
Expand Down Expand Up @@ -59,7 +58,7 @@ double Energy<T>::getEVrhoRho() const
}

template <class T>
double Energy<T>::evaluateEnergyIonsInVext()
double Energy<T>::evaluateEnergyIonsInVext(Ions& ions)
{
double energy = 0.;

Expand All @@ -69,12 +68,12 @@ double Energy<T>::evaluateEnergyIonsInVext()
//(*MPIdata::sout)<<"Energy<T>::evaluateEnergyIonsInVext()"<<std::endl;
double position[3];
std::vector<double> positions;
positions.reserve(3 * ions_.local_ions().size());
positions.reserve(3 * ions.local_ions().size());

// loop over ions
int nions = 0;
std::vector<Ion*>::const_iterator ion = ions_.local_ions().begin();
while (ion != ions_.local_ions().end())
std::vector<Ion*>::const_iterator ion = ions.local_ions().begin();
while (ion != ions.local_ions().end())
{
(*ion)->getPosition(position);
positions.push_back(position[0]);
Expand All @@ -88,9 +87,9 @@ double Energy<T>::evaluateEnergyIonsInVext()
pot_.getValVext(positions, val);

// loop over ions again
ion = ions_.local_ions().begin();
ion = ions.local_ions().begin();
int ion_index = 0;
while (ion != ions_.local_ions().end())
while (ion != ions.local_ions().end())
{
const double z = (*ion)->getZion();
// int ion_index=(*ion)->index();
Expand All @@ -112,16 +111,16 @@ double Energy<T>::evaluateEnergyIonsInVext()

template <class T>
double Energy<T>::evaluateTotal(const double ts, // in [Ha]
ProjectedMatricesInterface* projmatrices, const T& phi, const int verbosity,
std::ostream& os)
ProjectedMatricesInterface* projmatrices, Ions& ions, const T& phi,
const int verbosity, std::ostream& os)
{
eval_te_tm_.start();

Control& ct = *(Control::instance());

const double eself = ions_.energySelf();
const double ediff = ions_.energyDiff(ct.bcPoisson);
const double eipot = evaluateEnergyIonsInVext();
const double eself = ions.energySelf();
const double ediff = ions.energyDiff(ct.bcPoisson);
const double eipot = evaluateEnergyIonsInVext(ions);

const double eigsum = 0.5 * projmatrices->getExpectationH();

Expand Down
12 changes: 5 additions & 7 deletions src/Energy.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#define MGMOL_ENERGY_H

#include "Grid.h"
#include "Ions.h"
#include "Rho.h"
#include "SpreadPenaltyInterface.h"
#include "Timer.h"
Expand All @@ -20,7 +21,6 @@
#include <vector>

class Potentials;
class Ions;
class Electrostatic;
class ProjectedMatricesInterface;
class XConGrid;
Expand All @@ -29,7 +29,6 @@ template <class T>
class Energy
{
const pb::Grid& mygrid_;
const Ions& ions_;
const Potentials& pot_;
const Electrostatic& es_;
const Rho<T>& rho_;
Expand All @@ -45,16 +44,15 @@ class Energy
double getEVrhoRho() const;

public:
Energy(const pb::Grid&, const Ions&, const Potentials&,
const Electrostatic&, const Rho<T>&, const XConGrid&,
SpreadPenaltyInterface<T>*);
Energy(const pb::Grid&, const Potentials&, const Electrostatic&,
const Rho<T>&, const XConGrid&, SpreadPenaltyInterface<T>*);

static Timer eval_te_tm() { return eval_te_tm_; }

double evaluateTotal(const double ts, ProjectedMatricesInterface*,
double evaluateTotal(const double ts, ProjectedMatricesInterface*, Ions&,
const T& phi, const int, std::ostream&);

double evaluateEnergyIonsInVext();
double evaluateEnergyIonsInVext(Ions&);

void saveVofRho();
};
Expand Down
4 changes: 2 additions & 2 deletions src/GrassmanLineMinimization.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ Timer GrassmanLineMinimization<T>::update_states_tm_("Grassman_update_states");
//
// orthof=true: wants orthonormalized updated wave functions
template <class T>
int GrassmanLineMinimization<T>::updateWF(T& orbitals, Ions& /*ions*/,
int GrassmanLineMinimization<T>::updateWF(T& orbitals, Ions& ions,
const double precond_factor, const bool orthof, T& work_orbitals,
const bool accelerate, const bool print_res, const double atol)
{
Expand Down Expand Up @@ -61,7 +61,7 @@ int GrassmanLineMinimization<T>::updateWF(T& orbitals, Ions& /*ions*/,
// Update wavefunctions
const bool check_res = (atol > 0.);
double normRes = mgmol_strategy_->computeResidual(orbitals, work_orbitals,
*new_grad_, (print_res || check_res), check_res);
ions, *new_grad_, (print_res || check_res), check_res);
if (normRes < atol && check_res)
{
nl_update_tm_.stop();
Expand Down
14 changes: 7 additions & 7 deletions src/HamiltonianMVPSolver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,8 +160,8 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(

// compute energy at origin
const int printE = (ct.verbose > 1) ? 1 : 0;
double e0
= energy_->evaluateTotal(ts0, projmatrices, orbitals, printE, os_);
double e0 = energy_->evaluateTotal(
ts0, projmatrices, ions_, orbitals, printE, os_);

//
// compute energy at end for new H
Expand Down Expand Up @@ -189,8 +189,8 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
projmatrices->setHB2H();

// compute energy at end (beta=1.)
double e1
= energy_->evaluateTotal(ts1, projmatrices, orbitals, printE, os_);
double e1 = energy_->evaluateTotal(
ts1, projmatrices, ions_, orbitals, printE, os_);

//
// evaluate energy at mid-point
Expand Down Expand Up @@ -226,8 +226,8 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
projmatrices->setHB2H();

// compute energy at midpoint
double ei
= energy_->evaluateTotal(tsi, projmatrices, orbitals, printE, os_);
double ei = energy_->evaluateTotal(
tsi, projmatrices, ions_, orbitals, printE, os_);

// line minimization
double beta
Expand Down Expand Up @@ -285,7 +285,7 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(

// compute energy at end (beta=1.)
ei = energy_->evaluateTotal(
tsi, projmatrices, orbitals, printE, os_);
tsi, projmatrices, ions_, orbitals, printE, os_);

// line minimization
beta = minQuadPolynomialFrom3values(
Expand Down
Loading