force.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (c) 2017-2025  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.
//
// ---------------------------------------------------------------------
 
 
#ifndef force_H_
#define force_H_
#include "vselfBinsManager.h"
#include "dftParameters.h"
 
#include "constants.h"
#include "headers.h"
#include "meshMovementGaussian.h"
#include <dftd.h>
#include <oncvClass.h>
#include <AtomicCenteredNonLocalOperator.h>
#include <feevaluationWrapper.h>
#include <groupSymmetry.h>
namespace dftfe
{
  // forward declaration
  template <dftfe::utils::MemorySpace memory>
  class dftClass;
 
  /**
   * @brief computes configurational forces in KSDFT
   *
   * This class computes and stores the configurational forces corresponding to
   * geometry optimization. It uses the formulation in the paper by Motamarri
   * et.al. (https://link.aps.org/doi/10.1103/PhysRevB.97.165132) which provides
   * an unified approach to atomic forces corresponding to internal atomic
   * relaxation and cell stress corresponding to cell relaxation.
   *
   * @author Sambit Das
   */
  template <dftfe::utils::MemorySpace memorySpace>
  class forceClass
  {
    friend class dftClass<memorySpace>;
 
  public:
    /** @brief Constructor.
     *
     *  @param _dftPtr pointer to dftClass
     *  @param mpi_comm_parent parent mpi_communicator
     *  @param mpi_comm_domain domain decomposition mpi_communicator
     */
    forceClass(dftClass<memorySpace> *_dftPtr,
               const MPI_Comm        &mpi_comm_parent,
               const MPI_Comm        &mpi_comm_domain,
               const dftParameters   &dftParams);
 
    /** @brief initializes data structures inside forceClass assuming unmoved triangulation.
     *
     *  initUnmoved is the first step of the initialization/reinitialization of
     * force class when starting from a new unmoved triangulation. It creates
     * the dofHandler with linear finite elements and three components
     * corresponding to the three force components. It also creates the
     * corresponding constraint matrices which is why an unmoved triangulation
     * is necessary. Finally this function also initializes the gaussianMovePar
     * data member.
     *
     *  @param triangulation reference to unmoved triangulation where the mesh nodes have not
     *  been manually moved.
     *  @param isElectrostaticsMesh boolean parameter specifying whether this triangulatio is to be used for
     *  for the electrostatics part of the configurational force.
     *  @return void.
     */
    void
    initUnmoved(const dealii::Triangulation<3, 3>      &triangulation,
                const dealii::Triangulation<3, 3>      &serialTriangulation,
                const std::vector<std::vector<double>> &domainBoundingVectors,
                const bool                              isElectrostaticsMesh);
 
    /** @brief initializes data structures inside forceClass which depend on the moved mesh.
     *
     *  initMoved is the second step (first step call initUnmoved) of the
     * initialization/reinitialization of force class when starting from a new
     * mesh, and the first step when recomputing forces on a perturbed mesh.
     * initMoved assumes that the triangulation whose reference was passed to
     * the forceClass object in the initUnmoved call now has its nodes moved
     * such that all atomic positions lie on nodes.
     *
     *  @return void.
     */
    void
    initMoved(
      std::vector<const dealii::DoFHandler<3> *> &dofHandlerVectorMatrixFree,
      std::vector<const dealii::AffineConstraints<double> *>
                &constraintsVectorMatrixFree,
      const bool isElectrostaticsMesh);
 
    /** @brief initializes and precomputes pseudopotential related data structuers required for configurational force
     *  and stress computation.
     *
     *  This function is only activated for pseudopotential calculations and is
     * currently called when initializing/reinitializing the dftClass object.
     * This function initializes and precomputes the pseudopotential
     * datastructures for local and non-local parts. Separate internal function
     * calls are made for KB and ONCV projectors.
     *
     *  @return void.
     */
    void
    initPseudoData();
 
    /** @brief computes the configurational force on all atoms corresponding to a Gaussian generator,
     *  which represents perturbation of the underlying space.
     *
     *  The Gaussian generator is taken to be exp(-d_gaussianConstant*r^2), r
     * being the distance from the atom. Currently d_gaussianConstant is
     * hardcoded to be 4.0. To get the computed atomic forces use getAtomsForces
     *
     *  @return void.
     */
    void
    computeAtomsForces(
      const dealii::MatrixFree<3, double>       &matrixFreeData,
      const std::shared_ptr<groupSymmetryClass> &groupSymmetryPtr,
      const dispersionCorrection                &dispersionCorr,
      const dftfe::uInt                          eigenDofHandlerIndex,
      const dftfe::uInt                          smearedChargeQuadratureId,
      const dftfe::uInt                          lpspQuadratureIdElectro,
      const dealii::MatrixFree<3, double>       &matrixFreeDataElectro,
      const dftfe::uInt                          phiTotDofHandlerIndexElectro,
      const distributedCPUVec<double>           &phiTotRhoOutElectro,
      const std::vector<
        dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
        &rhoOutValues,
      const std::vector<
        dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
        &gradRhoOutValues,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &rhoTotalOutValuesLpsp,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &gradRhoTotalOutValuesLpsp,
      const std::map<dealii::CellId, std::vector<double>> &rhoCoreValues,
      const std::map<dealii::CellId, std::vector<double>> &gradRhoCoreValues,
      const std::map<dealii::CellId, std::vector<double>> &hessianRhoCoreValues,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &gradRhoCoreAtoms,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                                                          &hessianRhoCoreAtoms,
      const std::map<dealii::CellId, std::vector<double>> &pseudoVLocElectro,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                                              &pseudoVLocAtomsElectro,
      const dealii::AffineConstraints<double> &hangingPlusPBCConstraintsElectro,
      const vselfBinsManager                  &vselfBinsManagerElectro);
 
    /** @brief returns a copy of the configurational force on all global atoms.
     *
     *  computeAtomsForces must be called prior to this function call.
     *
     *  @return std::vector<double> flattened array of the configurational force on all atoms,
     *  the three force components on each atom being the leading dimension.
     * Units- Hartree/Bohr
     */
    std::vector<double> &
    getAtomsForces();
 
    /** @brief prints the currently stored configurational forces on atoms and the Gaussian generator constant
     *  used to compute them.
     *
     *  @return void.
     */
    void
    printAtomsForces();
 
    /** @brief Update force generator Gaussian constant.
     *
     *  @return void.
     */
    // void updateGaussianConstant(const double newGaussianConstant);
 
    /** @brief computes the configurational stress on the domain corresponding to
     *  affine deformation of the periodic cell.
     *
     *  This function cannot be called for fully non-periodic calculations.
     *
     *  @return void.
     */
    void
    computeStress(
      const dealii::MatrixFree<3, double>       &matrixFreeData,
      const std::shared_ptr<groupSymmetryClass> &groupSymmetryPtr,
      const dispersionCorrection                &dispersionCorr,
      const dftfe::uInt                          eigenDofHandlerIndex,
      const dftfe::uInt                          smearedChargeQuadratureId,
      const dftfe::uInt                          lpspQuadratureIdElectro,
      const dealii::MatrixFree<3, double>       &matrixFreeDataElectro,
      const dftfe::uInt                          phiTotDofHandlerIndexElectro,
      const distributedCPUVec<double>           &phiTotRhoOutElectro,
      const std::vector<
        dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
        &rhoOutValues,
      const std::vector<
        dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
        &gradRhoOutValues,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &rhoTotalOutValuesLpsp,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &gradRhoTotalOutValuesLpsp,
      const std::map<dealii::CellId, std::vector<double>> &pseudoVLocElectro,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &pseudoVLocAtomsElectro,
      const std::map<dealii::CellId, std::vector<double>> &rhoCoreValues,
      const std::map<dealii::CellId, std::vector<double>> &gradRhoCoreValues,
      const std::map<dealii::CellId, std::vector<double>> &hessianRhoCoreValues,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &gradRhoCoreAtoms,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                                              &hessianRhoCoreAtoms,
      const dealii::AffineConstraints<double> &hangingPlusPBCConstraintsElectro,
      const vselfBinsManager                  &vselfBinsManagerElectro);
 
    /** @brief prints the currently stored configurational stress tensor.
     *
     *  @return void.
     */
    void
    printStress();
 
    /** @brief returns a copy of the current stress tensor value.
     *
     *  computeStress must be call prior to this function call.
     *
     *  @return dealii::Tensor<2,3,double>  second order stress Tensor in Hartree/Bohr^3
     */
    dealii::Tensor<2, 3, double> &
    getStress();
 
    /** @brief get the value of Gaussian generator parameter (d_gaussianConstant).
     * Gaussian generator: Gamma(r)= exp(-d_gaussianConstant*r^2).
     *
     */
    // double getGaussianGeneratorParameter() const;
 
  private:
    /** @brief Locates and stores the global dof indices of d_dofHandlerForce whose cooridinates match
     *  with the atomic positions.
     *
     *  @return void.
     */
    void
    locateAtomCoreNodesForce(const dealii::DoFHandler<3> &dofHandlerForce,
                             const dealii::IndexSet &locally_owned_dofsForce,
                             std::map<std::pair<dftfe::uInt, dftfe::uInt>,
                                      dftfe::uInt>  &atomsForceDofs);
 
    void
    createBinObjectsForce(
      const dealii::DoFHandler<3>             &dofHandler,
      const dealii::DoFHandler<3>             &dofHandlerForce,
      const dealii::AffineConstraints<double> &hangingPlusPBCConstraints,
      const vselfBinsManager                  &vselfBinsManager,
      std::vector<std::vector<dealii::DoFHandler<3>::active_cell_iterator>>
        &cellsVselfBallsDofHandler,
      std::vector<std::vector<dealii::DoFHandler<3>::active_cell_iterator>>
        &cellsVselfBallsDofHandlerForce,
      std::vector<std::map<dealii::CellId, dftfe::uInt>>
        &cellsVselfBallsClosestAtomIdDofHandler,
      std::map<dftfe::uInt, dftfe::uInt> &AtomIdBinIdLocalDofHandler,
      std::vector<std::map<dealii::DoFHandler<3>::active_cell_iterator,
                           std::vector<dftfe::uInt>>>
        &cellFacesVselfBallSurfacesDofHandler,
      std::vector<std::map<dealii::DoFHandler<3>::active_cell_iterator,
                           std::vector<dftfe::uInt>>>
        &cellFacesVselfBallSurfacesDofHandlerForce);
 
    void
    configForceLinFEInit(
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dealii::MatrixFree<3, double> &matrixFreeDataElectro);
 
    void
    configForceLinFEFinalize();
 
    void
    computeConfigurationalForceEEshelbyTensorFPSPFnlLinFE(
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    eigenDofHandlerIndex,
      const dftfe::uInt                    smearedChargeQuadratureId,
      const dftfe::uInt                    lpspQuadratureIdElectro,
      const dealii::MatrixFree<3, double> &matrixFreeDataElectro,
      const dftfe::uInt                    phiTotDofHandlerIndexElectro,
      const distributedCPUVec<double>     &phiTotRhoOutElectro,
      const std::vector<
        dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
        &rhoOutValues,
      const std::vector<
        dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
        &gradRhoOutValues,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &rhoTotalOutValuesLpsp,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &gradRhoTotalOutValuesLpsp,
      const std::map<dealii::CellId, std::vector<double>> &rhoCoreValues,
      const std::map<dealii::CellId, std::vector<double>> &gradRhoCoreValues,
      const std::map<dealii::CellId, std::vector<double>> &hessianRhoCoreValues,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &gradRhoCoreAtoms,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                                                          &hessianRhoCoreAtoms,
      const std::map<dealii::CellId, std::vector<double>> &pseudoVLocElectro,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                             &pseudoVLocAtomsElectro,
      const vselfBinsManager &vselfBinsManagerElectro);
 
    void
    computeConfigurationalForceEEshelbyEElectroPhiTot(
      const dealii::MatrixFree<3, double> &matrixFreeDataElectro,
      const dftfe::uInt                    phiTotDofHandlerIndexElectro,
      const dftfe::uInt                    smearedChargeQuadratureId,
      const dftfe::uInt                    lpspQuadratureIdElectro,
      const distributedCPUVec<double>     &phiTotRhoOutElectro,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &rhoTotalOutValues,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &rhoTotalOutValuesLpsp,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &gradRhoTotalOutValues,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &gradRhoTotalOutValuesLpsp,
      const std::map<dealii::CellId, std::vector<double>> &pseudoVLocElectro,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                             &pseudoVLocAtomsElectro,
      const vselfBinsManager &vselfBinsManagerElectro);
 
    void
    computeConfigurationalForcePhiExtLinFE();
 
    void
    computeConfigurationalForceEselfLinFE(
      const dealii::DoFHandler<3>         &dofHandlerElectro,
      const vselfBinsManager              &vselfBinsManagerElectro,
      const dealii::MatrixFree<3, double> &matrixFreeDataElectro,
      const dftfe::uInt                    smearedChargeQuadratureId);
 
    void
    computeConfigurationalForceEselfNoSurfaceLinFE();
 
    void
    computeConfigurationalForceTotalLinFE(
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    eigenDofHandlerIndex,
      const dftfe::uInt                    smearedChargeQuadratureId,
      const dftfe::uInt                    lpspQuadratureIdElectro,
      const dealii::MatrixFree<3, double> &matrixFreeDataElectro,
      const dftfe::uInt                    phiTotDofHandlerIndexElectro,
      const distributedCPUVec<double>     &phiTotRhoOutElectro,
      const std::vector<
        dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
        &rhoOutValues,
      const std::vector<
        dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
        &gradRhoOutValues,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &rhoTotalOutValuesLpsp,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &gradRhoTotalOutValuesLpsp,
      const std::map<dealii::CellId, std::vector<double>> &rhoCoreValues,
      const std::map<dealii::CellId, std::vector<double>> &gradRhoCoreValues,
      const std::map<dealii::CellId, std::vector<double>> &hessianRhoCoreValues,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &gradRhoCoreAtoms,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                                                          &hessianRhoCoreAtoms,
      const std::map<dealii::CellId, std::vector<double>> &pseudoVLocElectro,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                             &pseudoVLocAtomsElectro,
      const vselfBinsManager &vselfBinsManagerElectro);
 
    void
    FPSPLocalGammaAtomsElementalContribution(
      std::map<dftfe::uInt, std::vector<double>>
                                          &forceContributionFPSPLocalGammaAtoms,
      dealii::FEValues<3>                 &feValues,
      dealii::FEFaceValues<3>             &feFaceValues,
      FEEvaluationWrapperClass<3>         &forceEval,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    phiTotDofHandlerIndexElectro,
      const dftfe::uInt                    cell,
      const dealii::AlignedVector<dealii::VectorizedArray<double>> &rhoQuads,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>> &gradRhoQuads,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                             &pseudoVLocAtoms,
      const vselfBinsManager &vselfBinsManager,
      const std::vector<std::map<dealii::CellId, dftfe::uInt>>
        &cellsVselfBallsClosestAtomIdDofHandler);
 
    void
    FPhiTotSmearedChargesGammaAtomsElementalContribution(
      std::map<dftfe::uInt, std::vector<double>>
                                  &forceContributionSmearedChargesGammaAtoms,
      FEEvaluationWrapperClass<3> &forceEval,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    cell,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>> &gradPhiTotQuads,
      const std::vector<dftfe::uInt> &nonTrivialAtomIdsMacroCell,
      const std::map<dealii::CellId, std::vector<dftfe::Int>>
        &bQuadAtomIdsAllAtoms,
      const dealii::AlignedVector<dealii::VectorizedArray<double>>
        &smearedbQuads);
 
    void
    FVselfSmearedChargesGammaAtomsElementalContribution(
      std::map<dftfe::uInt, std::vector<double>>
                                  &forceContributionSmearedChargesGammaAtoms,
      FEEvaluationWrapperClass<3> &forceEval,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    cell,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>>
                                     &gradVselfBinQuads,
      const std::vector<dftfe::uInt> &nonTrivialAtomIdsMacroCell,
      const std::map<dealii::CellId, std::vector<dftfe::Int>>
        &bQuadAtomIdsAllAtoms,
      const dealii::AlignedVector<dealii::VectorizedArray<double>>
        &smearedbQuads);
 
    void
    FNonlinearCoreCorrectionGammaAtomsElementalContribution(
      std::map<dftfe::uInt, std::vector<double>>
        &forceContributionFNonlinearCoreCorrectionGammaAtoms,
      FEEvaluationWrapperClass<3>         &forceEval,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    cell,
      const dealii::AlignedVector<dealii::VectorizedArray<double>> &vxcQuads,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &gradRhoCoreAtoms);
 
 
    void
    FNonlinearCoreCorrectionGammaAtomsElementalContribution(
      std::map<dftfe::uInt, std::vector<double>>
        &forceContributionFNonlinearCoreCorrectionGammaAtoms,
      FEEvaluationWrapperClass<3>         &forceEval,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    cell,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>> &derExcGradRho,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &hessianRhoCoreAtoms);
 
    void
    FNonlinearCoreCorrectionGammaAtomsElementalContributionSpinPolarized(
      std::map<dftfe::uInt, std::vector<double>>
        &forceContributionFNonlinearCoreCorrectionGammaAtoms,
      FEEvaluationWrapperClass<3>         &forceEval,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    cell,
      const dealii::AlignedVector<dealii::VectorizedArray<double>>
        &vxcQuadsSpin0,
      const dealii::AlignedVector<dealii::VectorizedArray<double>>
        &vxcQuadsSpin1,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>>
        &derExcGradRhoSpin0,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>>
        &derExcGradRhoSpin1,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &gradRhoCoreAtoms,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                &hessianRhoCoreAtoms,
      const bool isXCGGA = false);
 
    void
    distributeForceContributionFPSPLocalGammaAtoms(
      const std::map<dftfe::uInt, std::vector<double>>
        &forceContributionFPSPLocalGammaAtoms,
      const std::map<std::pair<dftfe::uInt, dftfe::uInt>, dftfe::uInt>
                                              &atomsForceDofs,
      const dealii::AffineConstraints<double> &constraintsNoneForce,
      distributedCPUVec<double>               &configForceVectorLinFE);
 
    void
    accumulateForceContributionGammaAtomsFloating(
      const std::map<dftfe::uInt, std::vector<double>>
                          &forceContributionLocalGammaAtoms,
      std::vector<double> &accumForcesVector);
 
 
    void
    FnlGammaAtomsElementalContribution(
      std::map<dftfe::uInt, std::vector<double>>
                                          &forceContributionFnlGammaAtoms,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      FEEvaluationWrapperClass<3>         &forceEvalNLP,
      const std::shared_ptr<
        AtomicCenteredNonLocalOperator<dataTypes::number, memorySpace>>
                                      nonLocalOp,
      dftfe::uInt                     numNonLocalAtomsCurrentProcess,
      const std::vector<dftfe::Int>  &globalChargeIdNonLocalAtoms,
      const std::vector<dftfe::uInt> &numberPseudoWaveFunctionsPerAtom,
      const dftfe::uInt               cell,
      const std::map<dealii::CellId, dftfe::uInt> &cellIdToCellNumberMap,
#ifdef USE_COMPLEX
      const std::vector<dataTypes::number>
        &projectorKetTimesPsiTimesVTimesPartOccContractionPsiQuadsFlattened,
#endif
      const std::vector<dataTypes::number> &zetaDeltaVQuadsFlattened,
      const std::vector<dataTypes::number> &
        projectorKetTimesPsiTimesVTimesPartOccContractionGradPsiQuadsFlattened);
 
 
    void
    FnlGammaxElementalContribution(
      dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>> &FVectQuads,
      const dealii::MatrixFree<3, double>                      &matrixFreeData,
      const dftfe::uInt                                         numQuadPoints,
      const std::shared_ptr<
        AtomicCenteredNonLocalOperator<dataTypes::number, memorySpace>>
                                      nonLocalOp,
      const dftfe::uInt               numNonLocalAtomsCurrentProcess,
      const std::vector<dftfe::Int>  &globalChargeIdNonLocalAtoms,
      const std::vector<dftfe::uInt> &numberPseudoWaveFunctionsPerAtom,
      const dftfe::uInt               cell,
      const std::map<dealii::CellId, dftfe::uInt> &cellIdToCellNumberMap,
      const std::vector<dataTypes::number>        &zetaDeltaVQuadsFlattened,
      const std::vector<dataTypes::number> &
        projectorKetTimesPsiTimesVTimesPartOccContractionGradPsiQuadsFlattened);
 
    void
    distributeForceContributionFnlGammaAtoms(
      const std::map<dftfe::uInt, std::vector<double>>
        &forceContributionFnlGammaAtoms);
 
    void
    stressEnlElementalContribution(
      dealii::Tensor<2, 3, double>        &stressContribution,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    numQuadPoints,
      const std::vector<double>           &jxwQuadsSubCells,
      const dftfe::uInt                    cell,
      const dftfe::uInt                    numNonLocalAtomsCurrentProcess,
      const std::shared_ptr<
        AtomicCenteredNonLocalOperator<dataTypes::number, memorySpace>>
                                      nonLocalOp,
      const std::vector<dftfe::uInt> &numberPseudoWaveFunctionsPerAtom,
      const std::map<dealii::CellId, dftfe::uInt> &cellIdToCellNumberMap,
      const std::vector<dataTypes::number> &zetalmDeltaVlProductDistImageAtoms,
#ifdef USE_COMPLEX
      const std::vector<dataTypes::number>
        &projectorKetTimesPsiTimesVTimesPartOccContractionPsiQuadsFlattened,
#endif
      const std::vector<dataTypes::number>
        &projectorKetTimesPsiTimesVTimesPartOccContractionGradPsiQuadsFlattened,
      const bool isSpinPolarized);
 
    void
    computeAtomsForcesGaussianGenerator(
      bool allowGaussianOverlapOnAtoms = false);
 
    void
    computeFloatingAtomsForces();
 
    void
    computeStressEself(
      const dealii::DoFHandler<3>         &dofHandlerElectro,
      const vselfBinsManager              &vselfBinsManagerElectro,
      const dealii::MatrixFree<3, double> &matrixFreeDataElectro,
      const dftfe::uInt                    smearedChargeQuadratureId);
 
    void
    computeStressEEshelbyEPSPEnlEk(
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    eigenDofHandlerIndex,
      const dftfe::uInt                    smearedChargeQuadratureId,
      const dftfe::uInt                    lpspQuadratureIdElectro,
      const dealii::MatrixFree<3, double> &matrixFreeDataElectro,
      const dftfe::uInt                    phiTotDofHandlerIndexElectro,
      const distributedCPUVec<double>     &phiTotRhoOutElectro,
      const std::vector<
        dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
        &rhoOutValues,
      const std::vector<
        dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>>
        &gradRhoOutValues,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &rhoTotalOutValuesLpsp,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &gradRhoTotalOutValuesLpsp,
      const std::map<dealii::CellId, std::vector<double>> &pseudoVLocElectro,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &pseudoVLocAtomsElectro,
      const std::map<dealii::CellId, std::vector<double>> &rhoCoreValues,
      const std::map<dealii::CellId, std::vector<double>> &gradRhoCoreValues,
      const std::map<dealii::CellId, std::vector<double>> &hessianRhoCoreValues,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &gradRhoCoreAtoms,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                             &hessianRhoCoreAtoms,
      const vselfBinsManager &vselfBinsManagerElectro);
 
    void
    computeStressEEshelbyEElectroPhiTot(
      const dealii::MatrixFree<3, double> &matrixFreeDataElectro,
      const dftfe::uInt                    phiTotDofHandlerIndexElectro,
      const dftfe::uInt                    smearedChargeQuadratureId,
      const dftfe::uInt                    lpspQuadratureIdElectro,
      const distributedCPUVec<double>     &phiTotRhoOutElectro,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &rhoTotalOutValues,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &rhoTotalOutValuesLpsp,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &gradRhoTotalOutValuesElectro,
      const dftfe::utils::MemoryStorage<double, dftfe::utils::MemorySpace::HOST>
        &gradRhoTotalOutValuesElectroLpsp,
      const std::map<dealii::CellId, std::vector<double>> &pseudoVLocElectro,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                             &pseudoVLocAtomsElectro,
      const vselfBinsManager &vselfBinsManagerElectro);
 
    void
    addEPSPStressContribution(
      dealii::FEValues<3>                 &feValues,
      dealii::FEFaceValues<3>             &feFaceValues,
      FEEvaluationWrapperClass<3>         &forceEval,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    phiTotDofHandlerIndexElectro,
      const dftfe::uInt                    cell,
      const dealii::AlignedVector<dealii::VectorizedArray<double>> &rhoQuads,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>> &gradRhoQuads,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                             &pseudoVLocAtoms,
      const vselfBinsManager &vselfBinsManager,
      const std::vector<std::map<dealii::CellId, dftfe::uInt>>
        &cellsVselfBallsClosestAtomIdDofHandler);
 
    void
    addENonlinearCoreCorrectionStressContribution(
      FEEvaluationWrapperClass<3>         &forceEval,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    cell,
      const dealii::AlignedVector<dealii::VectorizedArray<double>> &vxcQuads,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>> &derExcGradRho,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &gradRhoCoreAtoms,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &hessianRhoCoreAtoms);
 
    void
    addENonlinearCoreCorrectionStressContributionSpinPolarized(
      FEEvaluationWrapperClass<3>         &forceEval,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    cell,
      const dealii::AlignedVector<dealii::VectorizedArray<double>>
        &vxcQuadsSpin0,
      const dealii::AlignedVector<dealii::VectorizedArray<double>>
        &vxcQuadsSpin1,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>>
        &derExcGradRhoSpin0,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>>
        &derExcGradRhoSpin1,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
        &gradRhoCoreAtoms,
      const std::map<dftfe::uInt, std::map<dealii::CellId, std::vector<double>>>
                &hessianRhoCoreAtoms,
      const bool isXCGGA = false);
 
    void
    addEPhiTotSmearedStressContribution(
      FEEvaluationWrapperClass<3>         &forceEval,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    cell,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>> &gradPhiTotQuads,
      const std::vector<dftfe::uInt> &nonTrivialAtomImageIdsMacroCell,
      const std::map<dealii::CellId, std::vector<dftfe::Int>>
        &bQuadAtomIdsAllAtomsImages,
      const dealii::AlignedVector<dealii::VectorizedArray<double>>
        &smearedbQuads);
 
    void
    addEVselfSmearedStressContribution(
      FEEvaluationWrapperClass<3>         &forceEval,
      const dealii::MatrixFree<3, double> &matrixFreeData,
      const dftfe::uInt                    cell,
      const dealii::AlignedVector<
        dealii::Tensor<1, 3, dealii::VectorizedArray<double>>> &gradVselfQuads,
      const std::vector<dftfe::uInt> &nonTrivialAtomImageIdsMacroCell,
      const std::map<dealii::CellId, std::vector<dftfe::Int>>
        &bQuadAtomIdsAllAtomsImages,
      const dealii::AlignedVector<dealii::VectorizedArray<double>>
        &smearedbQuads);
 
    void
    computeElementalNonLocalPseudoOVDataForce();
 
 
    /// Parallel distributed vector field which stores the configurational force
    /// for each fem node corresponding to linear shape function generator (see
    /// equations 52-53 in
    /// (https://link.aps.org/doi/10.1103/PhysRevB.97.165132)). This vector
    /// doesn't contain contribution from terms which have sums over k points.
    distributedCPUVec<double> d_configForceVectorLinFE;
 
    /// Parallel distributed vector field which stores the configurational force
    /// for each fem node corresponding to linear shape function generator (see
    /// equations 52-53 in
    /// (https://link.aps.org/doi/10.1103/PhysRevB.97.165132)). This vector only
    /// containts contribution from the electrostatic part.
    distributedCPUVec<double> d_configForceVectorLinFEElectro;
 
#ifdef USE_COMPLEX
    /// Parallel distributed vector field which stores the configurational force
    /// for each fem node corresponding to linear shape function generator (see
    /// equations 52-53 in
    /// (https://link.aps.org/doi/10.1103/PhysRevB.97.165132)). This vector only
    /// contains contribution from terms which have sums over k points.
    distributedCPUVec<double> d_configForceVectorLinFEKPoints;
#endif
 
 
    std::vector<double> d_forceAtomsFloating;
 
#ifdef USE_COMPLEX
    std::vector<double> d_forceAtomsFloatingKPoints;
#endif
 
 
 
    /// Gaussian generator constant. Gaussian generator: Gamma(r)=
    /// exp(-d_gaussianConstant*r^2)
    /// FIXME: Until the hanging nodes surface integral issue is fixed use a
    /// value >=4.0
    // double d_gaussianConstant;
 
    /// Storage for configurational force on all global atoms.
    std::vector<double> d_globalAtomsForces;
 
    /// Storage for configurational stress tensor
    dealii::Tensor<2, 3, double> d_stress;
 
 
 
    /* Part of the stress tensor which is summed over k points.
     * It is a temporary data structure required for stress evaluation
     * (d_stress) when parallization over k points is on.
     */
    dealii::Tensor<2, 3, double> d_stressKPoints;
 
    /* Dont use true except for debugging forces only without mesh movement, as
     * gaussian ovelap on atoms for move mesh is by default set to false
     */
    const bool d_allowGaussianOverlapOnAtoms = false;
 
    /// pointer to dft class
    dftClass<memorySpace> *dftPtr;
 
    /// Finite element object for configurational force computation. Linear
    /// finite elements with three force field components are used.
    dealii::FESystem<3> FEForce;
 
    /* DofHandler on which we define the configurational force field. Each
     * geometric fem node has three dofs corresponding the the three force
     * components. The generator for the configurational force on the fem node
     * is the linear shape function attached to it. This DofHandler is based on
     * the same triangulation on which we solve the dft problem.
     */
    dealii::DoFHandler<3> d_dofHandlerForce;
 
    /* DofHandler on which we define the configurational force field from
     * electrostatic part (without psp). Each geometric fem node has three dofs
     * corresponding the the three force components. The generator for the
     * configurational force on the fem node is the linear shape function
     * attached to it. This DofHandler is based on the same triangulation on
     * which we solve the dft problem.
     */
    dealii::DoFHandler<3> d_dofHandlerForceElectro;
 
    /// Index of the d_dofHandlerForce in the MatrixFree object stored in
    /// dftClass. This is required to correctly use FEEvaluation class.
    dftfe::uInt d_forceDofHandlerIndex;
 
    /// Index of the d_dofHandlerForceElectro in the MatrixFree object stored in
    /// dftClass. This is required to correctly use FEEvaluation class.
    dftfe::uInt d_forceDofHandlerIndexElectro;
 
    /// dealii::IndexSet of locally owned dofs of in d_dofHandlerForce the
    /// current processor
    dealii::IndexSet d_locally_owned_dofsForce;
 
    /// dealii::IndexSet of locally owned dofs of in d_dofHandlerForceElectro
    /// the current processor
    dealii::IndexSet d_locally_owned_dofsForceElectro;
 
    /// dealii::IndexSet of locally relevant dofs of in d_dofHandlerForce the
    /// current processor
    dealii::IndexSet d_locally_relevant_dofsForce;
 
    /// dealii::IndexSet of locally relevant dofs of in d_dofHandlerForceElectro
    /// the current processor
    dealii::IndexSet d_locally_relevant_dofsForceElectro;
 
    /// Constraint matrix for hanging node and periodic constaints on
    /// d_dofHandlerForce.
    dealii::AffineConstraints<double> d_constraintsNoneForce;
 
    /// Constraint matrix for hanging node and periodic constaints on
    /// d_dofHandlerForceElectro.
    dealii::AffineConstraints<double> d_constraintsNoneForceElectro;
 
    /// Internal data: map < <atomId,force component>, globaldof in
    /// d_dofHandlerForce>
    std::map<std::pair<dftfe::uInt, dftfe::uInt>, dftfe::uInt> d_atomsForceDofs;
 
    /// Internal data: map < <atomId,force component>, globaldof in
    /// d_dofHandlerForceElectro>
    std::map<std::pair<dftfe::uInt, dftfe::uInt>, dftfe::uInt>
      d_atomsForceDofsElectro;
 
    /// Internal data: stores cell iterators of all cells in
    /// dftPtr->d_dofHandler which are part of the vself ball. Outer vector is
    /// over vself bins.
    std::vector<std::vector<dealii::DoFHandler<3>::active_cell_iterator>>
      d_cellsVselfBallsDofHandler;
 
    /// Internal data: stores cell iterators of all cells in d_dofHandlerForce
    /// which are part of the vself ball. Outer vector is over vself bins.
    std::vector<std::vector<dealii::DoFHandler<3>::active_cell_iterator>>
      d_cellsVselfBallsDofHandlerForce;
 
    /// Internal data: stores map of vself ball cell Id  to the closest atom Id
    /// of that cell. Outer vector over vself bins.
    std::vector<std::map<dealii::CellId, dftfe::uInt>>
      d_cellsVselfBallsClosestAtomIdDofHandler;
 
    /// Internal data: stores the map of atom Id (only in the local processor)
    /// to the vself bin Id.
    std::map<dftfe::uInt, dftfe::uInt> d_AtomIdBinIdLocalDofHandler;
 
    /* Internal data: stores the face ids of dftPtr->d_dofHandler (single
     * component field) on which to evaluate the vself ball surface integral in
     * the configurational force expression. Outer vector is over the vself
     * bins. Inner map is between the cell iterator and the vector of face ids
     * to integrate on for that cell iterator.
     */
    std::vector<std::map<dealii::DoFHandler<3>::active_cell_iterator,
                         std::vector<dftfe::uInt>>>
      d_cellFacesVselfBallSurfacesDofHandler;
 
    /* Internal data: stores the face ids of d_dofHandlerForce (three component
     * field) on which to evaluate the vself ball surface integral in the
     * configurational force expression. Outer vector is over the vself bins.
     * Inner map is between the cell iterator and the vector of face ids to
     * integrate on for that cell iterator.
     */
    std::vector<std::map<dealii::DoFHandler<3>::active_cell_iterator,
                         std::vector<dftfe::uInt>>>
      d_cellFacesVselfBallSurfacesDofHandlerForce;
 
    /// Internal data: stores cell iterators of all cells in
    /// dftPtr->d_dofHandler which are part of the vself ball. Outer vector is
    /// over vself bins.
    std::vector<std::vector<dealii::DoFHandler<3>::active_cell_iterator>>
      d_cellsVselfBallsDofHandlerElectro;
 
    /// Internal data: stores cell iterators of all cells in d_dofHandlerForce
    /// which are part of the vself ball. Outer vector is over vself bins.
    std::vector<std::vector<dealii::DoFHandler<3>::active_cell_iterator>>
      d_cellsVselfBallsDofHandlerForceElectro;
 
    /// Internal data: stores map of vself ball cell Id  to the closest atom Id
    /// of that cell. Outer vector over vself bins.
    std::vector<std::map<dealii::CellId, dftfe::uInt>>
      d_cellsVselfBallsClosestAtomIdDofHandlerElectro;
 
    /// Internal data: stores the map of atom Id (only in the local processor)
    /// to the vself bin Id.
    std::map<dftfe::uInt, dftfe::uInt> d_AtomIdBinIdLocalDofHandlerElectro;
 
    /* Internal data: stores the face ids of dftPtr->d_dofHandler (single
     * component field) on which to evaluate the vself ball surface integral in
     * the configurational force expression. Outer vector is over the vself
     * bins. Inner map is between the cell iterator and the vector of face ids
     * to integrate on for that cell iterator.
     */
    std::vector<std::map<dealii::DoFHandler<3>::active_cell_iterator,
                         std::vector<dftfe::uInt>>>
      d_cellFacesVselfBallSurfacesDofHandlerElectro;
 
    /* Internal data: stores the face ids of d_dofHandlerForce (three component
     * field) on which to evaluate the vself ball surface integral in the
     * configurational force expression. Outer vector is over the vself bins.
     * Inner map is between the cell iterator and the vector of face ids to
     * integrate on for that cell iterator.
     */
    std::vector<std::map<dealii::DoFHandler<3>::active_cell_iterator,
                         std::vector<dftfe::uInt>>>
      d_cellFacesVselfBallSurfacesDofHandlerForceElectro;
 
    std::map<dealii::CellId, dealii::DoFHandler<3>::active_cell_iterator>
      d_cellIdToActiveCellIteratorMapDofHandlerRhoNodalElectro;
 
    std::vector<distributedCPUVec<double>> d_gaussianWeightsVecAtoms;
 
 
 
    /// domain decomposition mpi_communicator
    const MPI_Comm d_mpiCommParent;
 
    /// domain decomposition mpi_communicator
    const MPI_Comm mpi_communicator;
 
    const dftParameters &d_dftParams;
 
    /// number of mpi processes in the current pool
    const dftfe::uInt n_mpi_processes;
 
    /// current mpi process id in the current pool
    const dftfe::uInt this_mpi_process;
 
    /// conditional stream object to enable printing only on root processor
    /// across all pools
    dealii::ConditionalOStream pcout;
  };
 
} // namespace dftfe
#endif