Skip to content

Commit

Permalink
Merged in kerkerFinal (pull request #326)
Browse files Browse the repository at this point in the history
Kerker mixing implemented using Helmholtz solve on 2p DoFHandler and nodal-mixing on 2p DofHandler

Approved-by: Sambit Das <[email protected]>
  • Loading branch information
phanimotamarri committed Oct 16, 2019
2 parents 8749215 + 0508add commit 89fcbdb
Show file tree
Hide file tree
Showing 26 changed files with 2,738 additions and 1,139 deletions.
1 change: 1 addition & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ SET(TARGET_SRC
./src/dft/vselfBinsManager.cc
./src/dft/energyCalculator.cc
./src/poisson/poissonSolverProblem.cc
./src/helmholtz/kerkerSolverProblem.cc
./src/dftOperator/kohnShamDFTOperator.cc
./src/dftOperator/operator.cc
./src/force/force.cc
Expand Down
11 changes: 7 additions & 4 deletions include/constants.h
Original file line number Diff line number Diff line change
Expand Up @@ -13,26 +13,29 @@
//
// ---------------------------------------------------------------------
//
// @author Sambit Das
// @author Sambit Das and Phani Motamarri
//

#ifndef constants_H_
#define constants_H_

namespace dftfe {
//
//Add prefic C_ to all constants
//Add prefix C_ to all constants
//

/// Boltzmann constant
const double C_kb =3.166811429e-06;
const double C_kb = 3.166811429e-06;

/// problem space dimensions
const int C_DIM=3;
const int C_DIM = 3;

/// 1d quadrature rule order
template <unsigned int FEOrder> constexpr unsigned int C_num1DQuad(){return FEOrder+1;}

//kerker Helmholtz solve polynomial Order
template <unsigned int FEOrder> constexpr unsigned int C_num1DKerkerPoly(){return FEOrder;}


/// 1d quadrature rule order for non-local part of pseudopotential
template <unsigned int FEOrder> constexpr unsigned int C_num1DQuadPSP()
Expand Down
3 changes: 2 additions & 1 deletion include/dealiiLinearSolver.h
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,8 @@ namespace dftfe {
void solve(dealiiLinearSolverProblem & problem,
const double absTolerance,
const unsigned int maxNumberIterations,
const unsigned int debugLevel = 0);
const unsigned int debugLevel = 0,
bool distributeFlag = true);

private:

Expand Down
53 changes: 49 additions & 4 deletions include/dft.h
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@
#include <triangulationManager.h>
#include <poissonSolverProblem.h>
#include <dealiiLinearSolver.h>
#include <kerkerSolverProblem.h>

#include <interpolation.h>
#include <xc.h>
Expand Down Expand Up @@ -314,6 +315,29 @@ namespace dftfe {
void initElectronicFields(const unsigned int usePreviousGroundStateFields=0);
void initPseudoPotentialAll();

/**
* create a dofHandler containing finite-element interpolating polynomial twice of the original polynomial
* required for Kerker mixing
* and initialize various objects related to this refined dofHandler
*/
void createpRefinedDofHandler(parallel::distributed::Triangulation<3> & triangulation);
void initpRefinedObjects();

/**
*@brief interpolate nodal data to quadrature data using FEEvaluation
*
*@param[in] matrixFreeData matrix free data object
*@param[in] nodalField nodal data to be interpolated
*@param[out] quadratureValueData to be computed at quadrature points
*@param[out] quadratureGradValueData to be computed at quadrature points
*@param[in] isEvaluateGradData denotes a flag to evaluate gradients or not
*/
void interpolateNodalDataToQuadratureData(dealii::MatrixFree<3,double> & matrixFreeData,
vectorType & nodalField,
std::map<dealii::CellId, std::vector<double> > & quadratureValueData,
std::map<dealii::CellId, std::vector<double> > & quadratureGradValueData,
const bool isEvaluateGradData);

/**
*@brief Finds the global dof ids of the nodes containing atoms.
*
Expand Down Expand Up @@ -355,7 +379,13 @@ namespace dftfe {
std::map<dealii::CellId, std::vector<double> > * _rhoValuesSpinPolarized,
std::map<dealii::CellId, std::vector<double> > * _gradRhoValuesSpinPolarized,
const bool isEvaluateGradRho,
const bool isConsiderSpectrumSplitting);
const bool isConsiderSpectrumSplitting,
const bool lobattoNodesFlag = false);

/**
*@brief computes density nodal data from wavefunctions
*/
void computeRhoNodalFromPSI(bool isConsiderSpectrumSplitting);


/**
Expand Down Expand Up @@ -422,6 +452,14 @@ namespace dftfe {
const std::map<dealii::CellId, std::vector<double> > *rhoQuadValues);


double totalCharge(const dealii::MatrixFree<3,double> & matrixFreeDataObject,
const vectorType & rhoNodalField);


double fieldl2Norm(const dealii::MatrixFree<3,double> & matrixFreeDataObject,
const vectorType & rhoNodalField);



/**
*@brief Computes net magnetization from the difference of local spin densities
Expand All @@ -447,6 +485,10 @@ namespace dftfe {
double mixing_anderson_spinPolarized();
double mixing_broyden();
double mixing_broyden_spinPolarized();
double nodalDensity_mixing_simple(kerkerSolverProblem<C_num1DKerkerPoly<FEOrder>()> & solverProblem,
dealiiLinearSolver & dealiiLinearSolver);
double nodalDensity_mixing_anderson(kerkerSolverProblem<C_num1DKerkerPoly<FEOrder>()> & solverProblem,
dealiiLinearSolver & dealiiLinearSolver);


/**
Expand Down Expand Up @@ -613,10 +655,10 @@ namespace dftfe {
* dealii based FE data structres
*/
FESystem<3> FE, FEEigen;
DoFHandler<3> dofHandler, dofHandlerEigen;
DoFHandler<3> dofHandler, dofHandlerEigen, d_dofHandlerPRefined;
unsigned int eigenDofHandlerIndex,phiExtDofHandlerIndex,phiTotDofHandlerIndex,forceDofHandlerIndex;
unsigned int densityDofHandlerIndex;
MatrixFree<3,double> matrix_free_data;
MatrixFree<3,double> matrix_free_data, d_matrixFreeDataPRefined;
std::map<types::global_dof_index, Point<3> > d_supportPoints, d_supportPointsEigen;
std::vector<const ConstraintMatrix * > d_constraintsVector;

Expand Down Expand Up @@ -676,7 +718,7 @@ namespace dftfe {
dftUtils::constraintMatrixInfo constraintsNoneDataInfo2;


ConstraintMatrix constraintsNone, constraintsNoneEigen, d_constraintsForTotalPotential, d_noConstraints;
ConstraintMatrix constraintsNone, constraintsNoneEigen, d_constraintsForTotalPotential, d_noConstraints, d_constraintsPRefined;


/**
Expand Down Expand Up @@ -704,6 +746,9 @@ namespace dftfe {
std::map<dealii::CellId, std::vector<double> > *rhoInValues, *rhoOutValues, *rhoInValuesSpinPolarized, *rhoOutValuesSpinPolarized;
std::deque<std::map<dealii::CellId,std::vector<double> >> rhoInVals, rhoOutVals, rhoInValsSpinPolarized, rhoOutValsSpinPolarized;

vectorType d_rhoInNodalValues, d_rhoOutNodalValues, d_preCondResidualVector;
std::deque<vectorType> d_rhoInNodalVals, d_rhoOutNodalVals;


std::map<dealii::CellId, std::vector<double> > * gradRhoInValues, *gradRhoInValuesSpinPolarized;
std::map<dealii::CellId, std::vector<double> > * gradRhoOutValues, *gradRhoOutValuesSpinPolarized;
Expand Down
8 changes: 4 additions & 4 deletions include/dftParameters.h
Original file line number Diff line number Diff line change
Expand Up @@ -35,18 +35,18 @@ namespace dftfe {
{

extern unsigned int finiteElementPolynomialOrder,n_refinement_steps,numberEigenValues,xc_id, spinPolarized, nkx,nky,nkz , offsetFlagX,offsetFlagY,offsetFlagZ;
extern unsigned int chebyshevOrder,numPass,numSCFIterations,maxLinearSolverIterations, mixingHistory, npool, numberWaveFunctionsForEstimate, numLevels;
extern unsigned int chebyshevOrder,numPass,numSCFIterations,maxLinearSolverIterations, mixingHistory, npool, numberWaveFunctionsForEstimate, numLevels, maxLinearSolverIterationsHelmholtz;

extern double radiusAtomBall, mixingParameter;
extern double lowerEndWantedSpectrum,relLinearSolverTolerance,selfConsistentSolverTolerance,TVal, start_magnetization;
extern double lowerEndWantedSpectrum,absLinearSolverTolerance,selfConsistentSolverTolerance,TVal, start_magnetization,absLinearSolverToleranceHelmholtz;

extern bool isPseudopotential, periodicX, periodicY, periodicZ, useSymm, timeReversal,pseudoTestsFlag, constraintMagnetization, writeDosFile, writeLdosFile,writeLocalizationLengths;
extern bool isPseudopotential, periodicX, periodicY, periodicZ, useSymm, timeReversal,pseudoTestsFlag, constraintMagnetization, writeDosFile, writeLdosFile,writeLocalizationLengths, pinnedNodeForPBC;
extern std::string meshFileName,coordinatesFile,domainBoundingVectorsFile,kPointDataFile, ionRelaxFlagsFile, orthogType, algoType, pseudoPotentialFile;

extern double outerAtomBallRadius, innerAtomBallRadius, meshSizeOuterDomain;
extern bool autoUserMeshParams;
extern double meshSizeInnerBall, meshSizeOuterBall;
extern double chebyshevTolerance, topfrac;
extern double chebyshevTolerance, topfrac, kerkerParameter;
extern std::string mixingMethod,ionOptSolver;


Expand Down
154 changes: 154 additions & 0 deletions include/kerkerSolverProblem.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
// ---------------------------------------------------------------------
//
// Copyright (c) 2019-2020 The Regents of the University of Michigan and DFT-FE authors.
//
// This file is part of the DFT-FE code.
//
// The DFT-FE code is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the DFT-FE distribution.
//
// ---------------------------------------------------------------------
//


#include <dealiiLinearSolverProblem.h>
#include <triangulationManager.h>

#ifndef kerkerSolverProblem_H_
#define kerkerSolverProblem_H_

namespace dftfe {

/**
* @brief poisson solver problem class template. template parameter FEOrder
* is the finite element polynomial order
*
* @author Phani Motamarri
*/
template<unsigned int FEOrder>
class kerkerSolverProblem: public dealiiLinearSolverProblem {

public:

/// Constructor
kerkerSolverProblem(const MPI_Comm &mpi_comm);




/**
* @brief initialize the matrix-free data structures
*
* @param matrixFreeData structure to hold quadrature rule, constraints vector and appropriate dofHandler
* @param constraintMatrix to hold constraints in the given problem
* @param x vector to be initialized using matrix-free object
*
*/
void init(dealii::MatrixFree<3,double> & matrixFreeData,
dealii::ConstraintMatrix & constraintMatrix,
vectorType & x);


/**
* @brief reinitialize data structures .
*
* @param x vector to store initial guess and solution
* @param gradResidualValues stores the gradient of difference of input electron-density and output electron-density
* @param kerkerMixingParameter used in Kerker mixing scheme which usually represents Thomas Fermi wavevector (k_{TF}**2).
*
*/
void reinit(vectorType & x,
const std::map<dealii::CellId,std::vector<double> > & gradResidualValues,
double kerkerMixingParameter);

/**
* @brief get the reference to x field
*
* @return reference to x field. Assumes x field data structure is already initialized
*/
vectorType & getX();

/**
* @brief Compute A matrix multipled by x.
*
*/
void vmult(vectorType &Ax,
const vectorType &x) const;

/**
* @brief Compute right hand side vector for the problem Ax = rhs.
*
* @param rhs vector for the right hand side values
*/
void computeRhs(vectorType & rhs);

/**
* @brief Jacobi preconditioning.
*
*/
void precondition_Jacobi(vectorType& dst,
const vectorType& src,
const double omega) const;

/**
* @brief distribute x to the constrained nodes.
*
*/
void distributeX();

/// function needed by dealii to mimic SparseMatrix for Jacobi preconditioning
void subscribe (const char *identifier=0) const{};

/// function needed by dealii to mimic SparseMatrix for Jacobi preconditioning
void unsubscribe (const char *identifier=0) const{};

/// function needed by dealii to mimic SparseMatrix
bool operator!= (double val) const {return true;};

private:

/**
* @brief required for the cell_loop operation in dealii's MatrixFree class
*
*/
void AX (const dealii::MatrixFree<3,double> & matrixFreeData,
vectorType &dst,
const vectorType &src,
const std::pair<unsigned int,unsigned int> &cell_range) const;


/**
* @brief Compute the diagonal of A.
*
*/
void computeDiagonalA();


/// storage for diagonal of the A matrix
vectorType d_diagonalA;


/// pointer to the x vector being solved for
vectorType * d_xPtr;

//kerker mixing constant
double d_gamma;

/// pointer to electron density cell and grad residual data
const std::map<dealii::CellId,std::vector<double> >* d_quadGradResidualValuesPtr;
const dealii::DoFHandler<3> * d_dofHandlerPRefinedPtr;
const dealii::ConstraintMatrix * d_constraintMatrixPRefinedPtr;
const dealii::MatrixFree<3,double> * d_matrixFreeDataPRefinedPtr;

const MPI_Comm mpi_communicator;
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
dealii::ConditionalOStream pcout;
};

}
#endif // kerkerSolverProblem_H_
3 changes: 2 additions & 1 deletion include/linearSolver.h
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,8 @@ namespace dftfe {
virtual void solve(dealiiLinearSolverProblem & problem,
const double relTolerance,
const unsigned int maxNumberIterations,
const unsigned int debugLevel = 0)=0;
const unsigned int debugLevel = 0,
bool distributeFlag = true)=0;

private:
};
Expand Down
2 changes: 1 addition & 1 deletion setup.sh
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ cxx_flagsRelease="-O3"
withELPA=ON

#Optmization flag: 1 for optimized mode and 0 for debug mode compilation
optimizedMode=1
optimizedMode=0

###########################################################################
#Usually, no changes are needed below this line
Expand Down
Loading

0 comments on commit 89fcbdb

Please sign in to comment.