Merge release

src/ABPG.cc

@@ -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)
 {
@@ -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();

src/DFTsolver.cc

@@ -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]

src/DavidsonSolver.cc

@@ -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))
@@ -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...
@@ -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_);

src/Energy.cc

@@ -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),
@@ -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.;
 
@@ -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]);
@@ -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();
@@ -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();
 

src/Energy.h

@@ -11,6 +11,7 @@
 #define MGMOL_ENERGY_H
 
 #include "Grid.h"
+#include "Ions.h"
 #include "Rho.h"
 #include "SpreadPenaltyInterface.h"
 #include "Timer.h"
@@ -20,7 +21,6 @@
 #include <vector>
 
 class Potentials;
-class Ions;
 class Electrostatic;
 class ProjectedMatricesInterface;
 class XConGrid;
@@ -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_;
@@ -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();
 };

src/GrassmanLineMinimization.cc

@@ -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)
 {
@@ -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();

src/HamiltonianMVPSolver.cc

@@ -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
@@ -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
@@ -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
@@ -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(