Skip to content

Commit

Permalink
CGDESCENT and LBFGS: supress trivial ion update call and don't comput…
Browse files Browse the repository at this point in the history 
…e forces if only energy is needed.
  • Loading branch information
@dsambit
dsambit committed Nov 13, 2019
1 parent 123b91f commit aefae00
Show file tree
Hide file tree
Showing 9 changed files with 147 additions and 117 deletions.
21 changes: 18 additions & 3 deletions include/cg_descent_wrapper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,12 @@ inline double cgd_value(double* x, INT n)
0,
MPI_COMM_WORLD);

functionToBeMinimized->update(solutionInc);
double sumAbsDisp=0.0;
for(unsigned int i = 0; i < solutionInc.size(); ++i)
sumAbsDisp+=std::fabs(solutionInc[i]);

if (sumAbsDisp>1e-12)
functionToBeMinimized->update(solutionInc,false);

std::vector<double> funcValue;
functionToBeMinimized->value(funcValue);
Expand Down Expand Up @@ -82,7 +87,12 @@ inline void cgd_gradient(double* g, double* x, INT n)
0,
MPI_COMM_WORLD);

functionToBeMinimized->update(solutionInc);
double sumAbsDisp=0.0;
for(unsigned int i = 0; i < solutionInc.size(); ++i)
sumAbsDisp+=std::fabs(solutionInc[i]);

if (sumAbsDisp>1e-12)
functionToBeMinimized->update(solutionInc);

functionToBeMinimized->save();

Expand Down Expand Up @@ -119,7 +129,12 @@ inline double cgd_value_gradient(double* g, double* x, INT n)
0,
MPI_COMM_WORLD);

functionToBeMinimized->update(solutionInc);
double sumAbsDisp=0.0;
for(unsigned int i = 0; i < solutionInc.size(); ++i)
sumAbsDisp+=std::fabs(solutionInc[i]);

if (sumAbsDisp>1e-12)
functionToBeMinimized->update(solutionInc);

functionToBeMinimized->save();

Expand Down
4 changes: 2 additions & 2 deletions include/dft.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -153,7 +153,7 @@ namespace dftfe {
/**
* @brief Kohn-Sham ground-state solve using SCF iteration
*/
void solve();
void solve(const bool computeForces=true);

/**
* @brief Number of Kohn-Sham eigen values to be computed
Expand Down Expand Up @@ -501,7 +501,7 @@ namespace dftfe {
* However, it works for time reversal symmetry.
*
*/
void computeElectrostaticEnergyHRefined();
void computeElectrostaticEnergyHRefined(const bool computeForces=true);
void computeElectrostaticEnergyPRefined();
/**
*@brief Computes Fermi-energy obtained by imposing constraint on the number of electrons
Expand Down
3 changes: 2 additions & 1 deletion include/geoOptCell.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,8 @@ namespace dftfe {
*
* @param solution delta(epsilon).
*/
void update(const std::vector<double> & solution);
void update(const std::vector<double> & solution,
const bool computeForces=true);

/**
* @brief create checkpoint file for current domain bounding vectors and atomic coordinates.
Expand Down
3 changes: 2 additions & 1 deletion include/geoOptIon.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,8 @@ namespace dftfe {
* @param solution displacement of the atoms with respect to their current position.
* The size of the solution vector is equal to the number of unknowns.
*/
void update(const std::vector<double> & solution);
void update(const std::vector<double> & solution,
const bool computeForces=true);

/**
* @brief create checkpoint file for current domain bounding vectors and atomic coordinates.
Expand Down
3 changes: 2 additions & 1 deletion include/nonlinearSolverProblem.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,8 @@ namespace dftfe {
*
* @param solution Updated solution.
*/
virtual void update(const std::vector<double> & solution) = 0;
virtual void update(const std::vector<double> & solution,
const bool computeForces=true) = 0;

/**
* @brief Obtain current solution.
Expand Down
110 changes: 58 additions & 52 deletions src/dft/dft.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -913,7 +913,7 @@ namespace dftfe {
//dft solve
//
template<unsigned int FEOrder>
void dftClass<FEOrder>::solve()
void dftClass<FEOrder>::solve(const bool computeForces)
{

//
Expand Down Expand Up @@ -1852,31 +1852,34 @@ namespace dftfe {

computing_timer.enter_section("Ion force computation");
computingTimerStandard.enter_section("Ion force computation");
forcePtr->computeAtomsForces(matrix_free_data,
eigenDofHandlerIndex,
phiExtDofHandlerIndex,
phiTotDofHandlerIndex,
d_phiTotRhoIn,
d_phiTotRhoOut,
d_phiExt,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager,
matrix_free_data,
phiTotDofHandlerIndex,
phiExtDofHandlerIndex,
d_phiTotRhoOut,
d_phiExt,
*rhoOutValues,
*gradRhoOutValues,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager);
forcePtr->printAtomsForces();
if (computeForces)
{
forcePtr->computeAtomsForces(matrix_free_data,
eigenDofHandlerIndex,
phiExtDofHandlerIndex,
phiTotDofHandlerIndex,
d_phiTotRhoIn,
d_phiTotRhoOut,
d_phiExt,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager,
matrix_free_data,
phiTotDofHandlerIndex,
phiExtDofHandlerIndex,
d_phiTotRhoOut,
d_phiExt,
*rhoOutValues,
*gradRhoOutValues,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager);
forcePtr->printAtomsForces();
}
computingTimerStandard.exit_section("Ion force computation");
computing_timer.exit_section("Ion force computation");
}
Expand All @@ -1888,38 +1891,41 @@ namespace dftfe {

computing_timer.enter_section("Cell stress computation");
computingTimerStandard.enter_section("Cell stress computation");
forcePtr->computeStress(matrix_free_data,
eigenDofHandlerIndex,
phiExtDofHandlerIndex,
phiTotDofHandlerIndex,
d_phiTotRhoIn,
d_phiTotRhoOut,
d_phiExt,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager,
matrix_free_data,
phiTotDofHandlerIndex,
phiExtDofHandlerIndex,
d_phiTotRhoOut,
d_phiExt,
*rhoOutValues,
*gradRhoOutValues,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager);
forcePtr->printStress();
if (computeForces)
{
forcePtr->computeStress(matrix_free_data,
eigenDofHandlerIndex,
phiExtDofHandlerIndex,
phiTotDofHandlerIndex,
d_phiTotRhoIn,
d_phiTotRhoOut,
d_phiExt,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager,
matrix_free_data,
phiTotDofHandlerIndex,
phiExtDofHandlerIndex,
d_phiTotRhoOut,
d_phiExt,
*rhoOutValues,
*gradRhoOutValues,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager);
forcePtr->printStress();
}
computingTimerStandard.exit_section("Cell stress computation");
computing_timer.exit_section("Cell stress computation");
}
#endif

if(dftParameters::electrostaticsHRefinement)
computeElectrostaticEnergyHRefined();
computeElectrostaticEnergyHRefined(computeForces);

if(dftParameters::electrostaticsPRefinement)
computeElectrostaticEnergyPRefined();
Expand Down
112 changes: 59 additions & 53 deletions src/dft/electrostaticHRefinedEnergy.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@


template<unsigned int FEOrder>
void dftClass<FEOrder>::computeElectrostaticEnergyHRefined()
void dftClass<FEOrder>::computeElectrostaticEnergyHRefined(const bool computeForces)
{
computing_timer.enter_section("h refinement electrostatics");
computingTimerStandard.enter_section("h refinement electrostatics");
Expand Down Expand Up @@ -333,7 +333,7 @@ void dftClass<FEOrder>::computeElectrostaticEnergyHRefined()
delzRhoNodalFieldRefined.update_ghost_values();
}


//
//move the refined mesh so that it forms exact subdivison of coarse moved mesh
//
Expand Down Expand Up @@ -363,7 +363,7 @@ void dftClass<FEOrder>::computeElectrostaticEnergyHRefined()
}




//
//call init for the force computation subsequently
Expand Down Expand Up @@ -694,31 +694,34 @@ void dftClass<FEOrder>::computeElectrostaticEnergyHRefined()

computing_timer.enter_section("Ion force computation");
computingTimerStandard.enter_section("Ion force computation");
forcePtr->computeAtomsForces(matrix_free_data,
eigenDofHandlerIndex,
phiExtDofHandlerIndex,
phiTotDofHandlerIndex,
d_phiTotRhoIn,
d_phiTotRhoOut,
d_phiExt,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager,
matrixFreeDataHRefined,
phiTotDofHandlerIndexHRefined,
phiExtDofHandlerIndexHRefined,
phiTotRhoOutHRefined,
phiExtHRefined,
rhoOutHRefinedQuadValues,
gradRhoOutHRefinedQuadValues,
pseudoVLocHRefined,
gradPseudoVLocHRefined,
gradPseudoVLocAtomsHRefined,
onlyHangingNodeConstraints,
vselfBinsManagerHRefined);
forcePtr->printAtomsForces();
if (computeForces)
{
forcePtr->computeAtomsForces(matrix_free_data,
eigenDofHandlerIndex,
phiExtDofHandlerIndex,
phiTotDofHandlerIndex,
d_phiTotRhoIn,
d_phiTotRhoOut,
d_phiExt,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager,
matrixFreeDataHRefined,
phiTotDofHandlerIndexHRefined,
phiExtDofHandlerIndexHRefined,
phiTotRhoOutHRefined,
phiExtHRefined,
rhoOutHRefinedQuadValues,
gradRhoOutHRefinedQuadValues,
pseudoVLocHRefined,
gradPseudoVLocHRefined,
gradPseudoVLocAtomsHRefined,
onlyHangingNodeConstraints,
vselfBinsManagerHRefined);
forcePtr->printAtomsForces();
}
computingTimerStandard.exit_section("Ion force computation");
computing_timer.exit_section("Ion force computation");
}
Expand All @@ -728,31 +731,34 @@ void dftClass<FEOrder>::computeElectrostaticEnergyHRefined()

computing_timer.enter_section("Cell stress computation");
computingTimerStandard.enter_section("Cell stress computation");
forcePtr->computeStress(matrix_free_data,
eigenDofHandlerIndex,
phiExtDofHandlerIndex,
phiTotDofHandlerIndex,
d_phiTotRhoIn,
d_phiTotRhoOut,
d_phiExt,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager,
matrixFreeDataHRefined,
phiTotDofHandlerIndexHRefined,
phiExtDofHandlerIndexHRefined,
phiTotRhoOutHRefined,
phiExtHRefined,
rhoOutHRefinedQuadValues,
gradRhoOutHRefinedQuadValues,
pseudoVLocHRefined,
gradPseudoVLocHRefined,
gradPseudoVLocAtomsHRefined,
onlyHangingNodeConstraints,
vselfBinsManagerHRefined);
forcePtr->printStress();
if (computeForces)
{
forcePtr->computeStress(matrix_free_data,
eigenDofHandlerIndex,
phiExtDofHandlerIndex,
phiTotDofHandlerIndex,
d_phiTotRhoIn,
d_phiTotRhoOut,
d_phiExt,
d_pseudoVLoc,
d_gradPseudoVLoc,
d_gradPseudoVLocAtoms,
d_noConstraints,
d_vselfBinsManager,
matrixFreeDataHRefined,
phiTotDofHandlerIndexHRefined,
phiExtDofHandlerIndexHRefined,
phiTotRhoOutHRefined,
phiExtHRefined,
rhoOutHRefinedQuadValues,
gradRhoOutHRefinedQuadValues,
pseudoVLocHRefined,
gradPseudoVLocHRefined,
gradPseudoVLocAtomsHRefined,
onlyHangingNodeConstraints,
vselfBinsManagerHRefined);
forcePtr->printStress();
}
computingTimerStandard.exit_section("Cell stress computation");
computing_timer.exit_section("Cell stress computation");
}
Expand Down
4 changes: 2 additions & 2 deletions src/geoOpt/geoOptCell.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -309,7 +309,7 @@ void geoOptCell<FEOrder>::precondition(std::vector<double> & s,
}

template<unsigned int FEOrder>
void geoOptCell<FEOrder>::update(const std::vector<double> & solution)
void geoOptCell<FEOrder>::update(const std::vector<double> & solution, const bool computeForces)
{


Expand Down Expand Up @@ -379,7 +379,7 @@ void geoOptCell<FEOrder>::update(const std::vector<double> & solution)
d_totalUpdateCalls+=1;
dftPtr->deformDomain(deformationGradient);

dftPtr->solve();
dftPtr->solve(computeForces);
// if ion optimization is on, then for every cell relaxation also relax the atomic forces
if (dftParameters::isIonOpt)
{
Expand Down
Loading

0 comments on commit aefae00

Please sign in to comment.