Skip to content

Commit

Permalink
Fix: error with DFTU force&stress with NSPIN=4 (#4825)
Browse files Browse the repository at this point in the history 
* Fix: error with DFTU force&stress with NSPIN=4

* Fix: add test

---------

Co-authored-by: dyzheng <[email protected]>
  • Loading branch information
@dyzheng
dyzheng and dyzheng authored Jul 29, 2024
1 parent 3845b76 commit d8cca53
Show file tree
Hide file tree
Showing 3 changed files with 47 additions and 18 deletions.
57 changes: 41 additions & 16 deletions source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_force_stress.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,18 +113,35 @@ void DFTU<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
}
// first iteration to calculate occupation matrix
std::vector<double> occ(tlp1 * tlp1 * this->nspin, 0);
for (int i = 0; i < occ.size(); i++)
if(this->nspin ==2)
{
const int is = i / (tlp1 * tlp1);
const int ii = i % (tlp1 * tlp1);
occ[i] = this->dftu->locale[iat0][target_L][0][is].c[ii];
for (int i = 0; i < occ.size(); i++)
{
const int is = i / (tlp1 * tlp1);
const int ii = i % (tlp1 * tlp1);
occ[i] = this->dftu->locale[iat0][target_L][0][is].c[ii];
}
}
else
{
for (int i = 0; i < occ.size(); i++)
{
occ[i] = this->dftu->locale[iat0][target_L][0][0].c[i];
}
}

// calculate VU
const double u_value = this->dftu->U[T0];
std::vector<double> VU(occ.size());
double eu_tmp = 0;
this->cal_v_of_u(occ, tlp1, u_value, &VU[0], eu_tmp);
if(this->nspin == 4)
{
for (int i = 0; i < VU.size(); i++)
{
VU[i] /= 2.0;
}
}

// second iteration to calculate force and stress
// calculate Force for atom J
Expand Down Expand Up @@ -255,15 +272,19 @@ void DFTU<OperatorLCAO<TK, TR>>::cal_force_IJR(const int& iat1,
const int m_size = int(sqrt(vu_in.size() / nspin));
const int m_size2 = m_size * m_size;
// step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4
std::vector<int> step_trace(npol, 0);
std::vector<int> step_trace(npol * npol, 0);
if (npol == 2) {
step_trace[1] = col_indexes.size() + 1;
step_trace[1] = 1;
step_trace[2] = col_indexes.size();
step_trace[3] = col_indexes.size() + 1;
}
double tmp[3];
// calculate the local matrix
for (int is = 0; is < nspin; is++)
{
double* dm_pointer = dmR_pointer[is]->get_pointer();
const int is0 = nspin==2 ? is : 0;
const int step_is = nspin==4 ? is : 0;
const double* dm_pointer = dmR_pointer[is0]->get_pointer();
for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol)
{
const std::vector<double>& nlm1 = nlm1_all.find(row_indexes[iw1l])->second;
Expand All @@ -277,11 +298,11 @@ void DFTU<OperatorLCAO<TK, TR>>::cal_force_IJR(const int& iat1,
{
for (int m2 = 0; m2 < m_size; m2++)
{
tmp[0] = vu_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size] * nlm2[m2] * dm_pointer[0];
tmp[0] = vu_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size] * nlm2[m2] * dm_pointer[step_trace[step_is]];
tmp[1]
= vu_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size * 2] * nlm2[m2] * dm_pointer[0];
= vu_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size * 2] * nlm2[m2] * dm_pointer[step_trace[step_is]];
tmp[2]
= vu_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size * 3] * nlm2[m2] * dm_pointer[0];
= vu_in[m1 * m_size + m2 + is * m_size2] * nlm1[m1 + m_size * 3] * nlm2[m2] * dm_pointer[step_trace[step_is]];
// force1 = - VU * <d phi_{I,R1}/d R1|chi_m> * <chi_m'|phi_{J,R2}>
// force2 = - VU * <phi_{I,R1}|d chi_m/d R0> * <chi_m'|phi_{J,R2>}
force1[0] += tmp[0];
Expand All @@ -292,7 +313,7 @@ void DFTU<OperatorLCAO<TK, TR>>::cal_force_IJR(const int& iat1,
force2[2] -= tmp[2];
}
}
dm_pointer++;
dm_pointer += npol;
}
dm_pointer += (npol - 1) * col_indexes.size();
}
Expand Down Expand Up @@ -324,14 +345,18 @@ void DFTU<OperatorLCAO<TK, TR>>::cal_stress_IJR(const int& iat1,
const int m_size = int(sqrt(vu_in.size() / nspin));
const int m_size2 = m_size * m_size;
// step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4
std::vector<int> step_trace(npol, 0);
std::vector<int> step_trace(npol * npol, 0);
if (npol == 2) {
step_trace[1] = col_indexes.size() + 1;
step_trace[1] = 1;
step_trace[2] = col_indexes.size();
step_trace[3] = col_indexes.size() + 1;
}
// calculate the local matrix
for (int is = 0; is < nspin; is++)
{
double* dm_pointer = dmR_pointer[is]->get_pointer();
const int is0 = nspin==2 ? is : 0;
const int step_is = nspin==4 ? is : 0;
const double* dm_pointer = dmR_pointer[is0]->get_pointer();
for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol)
{
const std::vector<double>& nlm1 = nlm1_all.find(row_indexes[iw1l])->second;
Expand All @@ -345,7 +370,7 @@ void DFTU<OperatorLCAO<TK, TR>>::cal_stress_IJR(const int& iat1,
{
for (int m2 = 0; m2 < m_size; m2++)
{
double tmp = vu_in[m1 * m_size + m2 + is * m_size2] * dm_pointer[0];
double tmp = vu_in[m1 * m_size + m2 + is * m_size2] * dm_pointer[step_trace[step_is]];
// std::cout<<__FILE__<<__LINE__<<" "<<tmp<<" "<<m1<<" "<<m2<<" "<<nlm1[m1 + m_size * 2]<<"
// "<<nlm2[m2 + m_size * 2]<<" "<<dis1.y<<" "<<dis2.y<<std::endl;
stress[0]
Expand All @@ -365,7 +390,7 @@ void DFTU<OperatorLCAO<TK, TR>>::cal_stress_IJR(const int& iat1,
+ nlm1[m1] * nlm2[m2 + m_size * 3] * dis2.z);
}
}
dm_pointer++;
dm_pointer += npol;
}
dm_pointer += (npol - 1) * col_indexes.size();
}
Expand Down
2 changes: 2 additions & 0 deletions tests/integrate/260_NO_DJ_PK_PU_SO/INPUT
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,8 @@ basis_type lcao
gamma_only 0
noncolin 1
lspinorb 1
cal_force 1
cal_stress 1

#Parameter DFT+U
dft_plus_u 1
Expand Down
6 changes: 4 additions & 2 deletions tests/integrate/260_NO_DJ_PK_PU_SO/result.ref
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
etotref -6789.2817503491560274
etotref -6789.2817503491551179
etotperatomref -3394.6408751746
totaltimeref 4.39
totalforceref 11.319866
totalstressref 4893.060045
totaltimeref 5.40

0 comments on commit d8cca53

Please sign in to comment.