operator_force_stress_utils.h

Go to the documentation of this file.

#ifndef OPERATOR_FORCE_STRESS_UTILS_H
#define OPERATOR_FORCE_STRESS_UTILS_H
 
#include "source_base/matrix.h"
#include "source_cell/unitcell.h"
#include <vector>
 
namespace OperatorForceStress {
 
inline void setup_step_trace(int npol, int col_size, std::vector<int>& step_trace)
{
    step_trace.resize(npol * npol, 0);
    if (npol == 2)
    {
        step_trace[1] = 1;
        step_trace[2] = col_size;
        step_trace[3] = col_size + 1;
    }
}
 
struct OrbitalQuantumNumbers
{
    int L;  
    int N;  
    int m;  
    int M;  
};
 
inline OrbitalQuantumNumbers get_orbital_qn(const Atom& atom, int iw)
{
    OrbitalQuantumNumbers qn;
    qn.L = atom.iw2l[iw];
    qn.N = atom.iw2n[iw];
    qn.m = atom.iw2m[iw];
    qn.M = (qn.m % 2 == 0) ? -qn.m / 2 : (qn.m + 1) / 2;
    return qn;
}
 
template <typename T>
inline double get_real_part(const T& val)
{
    return val.real();
}
 
template <>
inline double get_real_part<double>(const double& val)
{
    return val;
}
 
inline void rearrange_stress_matrix(ModuleBase::matrix& stress)
{
    stress.c[8] = stress.c[5]; // stress(2,2)
    stress.c[7] = stress.c[4]; // stress(2,1)
    stress.c[6] = stress.c[2]; // stress(2,0)
    stress.c[5] = stress.c[4]; // stress(1,2)
    stress.c[4] = stress.c[3]; // stress(1,1)
    stress.c[3] = stress.c[1]; // stress(1,0)
}
 
void finalize_force_stress(
    bool cal_force,
    bool cal_stress,
    const UnitCell* ucell,
    const std::vector<double>& stress_tmp,
    ModuleBase::matrix& force,
    ModuleBase::matrix& stress,
    double force_factor = 2.0,
    double stress_factor = 2.0);
 
} // namespace OperatorForceStress
 
#endif // OPERATOR_FORCE_STRESS_UTILS_H