Skip to content

Commit 9e7a043

Browse files
AsTonyshmentAsTonyshment
authored
Fix: Support out_current=2 with multiple MPI processes (deepmodeling#7624)
1 parent 38eb1db commit 9e7a043

4 files changed

Lines changed: 42 additions & 85 deletions

File tree

‎docs/advanced/input_files/input-main.md‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -4372,9 +4372,10 @@
43724372
### out_current
43734373

43744374
- **Type**: Integer
4375-
- **Description**: - 0: Do not output current.
4376-
- 1: Output current using the two-center integral, faster.
4377-
- 2: Output current using the matrix commutation, more precise.
4375+
- **Description**: Controls the current-density output method for LCAO RT-TDDFT.
4376+
- 0: Do not output current.
4377+
- 1: Explicitly construct the velocity operator from the momentum, vector-potential, and KB nonlocal-pseudopotential terms using two-center integral / spherical grid integral: $$\hat{v}_{\alpha}=-\mathrm{i}\nabla_{\alpha}+A_{\alpha}(t)+\mathrm{i}\left[\widetilde{V}_{\mathrm{NL}}^{\mathrm{KB}},r_{\alpha}\right],$$ where $\widetilde{V}_{\mathrm{NL}}^{\mathrm{KB}}=\mathrm{e}^{-\mathrm{i}\boldsymbol{A}(t)\cdot\boldsymbol{r}}\hat{V}_{\mathrm{NL}}^{\mathrm{KB}}\mathrm{e}^{\mathrm{i}\boldsymbol{A}(t)\cdot\boldsymbol{r}}$. $\boldsymbol{A}(t)$ is nonzero only for the velocity gauge (td_stype=1); otherwise $\boldsymbol{A}(t)=0$. Other nonlocal Hamiltonian terms (e.g., EXX) are not included explicitly.
4378+
- 2: Use the full Hamiltonian to construct the generalized velocity matrix in a nonorthogonal NAO basis: $$\widetilde{v}_{\alpha}=\partial_{\alpha}H+\mathrm{i}HS^{-1}\mathcal{R}_{\alpha}-\mathrm{i}\mathcal{R}_{\alpha}S^{-1}H-HS^{-1}\partial_{\alpha}S.$$ This includes all contributions available in the real-space Hamiltonian matrix when enabled. This method is more general but more expensive.
43784379
- **Default**: 0
43794380

43804381
### out_current_k

‎docs/parameters.yaml‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3468,9 +3468,10 @@ parameters:
34683468
category: "RT-TDDFT: Real-Time Time-Dependent Density Functional Theory"
34693469
type: Integer
34703470
description: |
3471+
Controls the current-density output method for LCAO RT-TDDFT.
34713472
* 0: Do not output current.
3472-
* 1: Output current using the two-center integral, faster.
3473-
* 2: Output current using the matrix commutation, more precise.
3473+
* 1: Explicitly construct the velocity operator from the momentum, vector-potential, and KB nonlocal-pseudopotential terms using two-center integral / spherical grid integral: $$\hat{v}_{\alpha}=-\mathrm{i}\nabla_{\alpha}+A_{\alpha}(t)+\mathrm{i}\left[\widetilde{V}_{\mathrm{NL}}^{\mathrm{KB}},r_{\alpha}\right],$$ where $\widetilde{V}_{\mathrm{NL}}^{\mathrm{KB}}=\mathrm{e}^{-\mathrm{i}\boldsymbol{A}(t)\cdot\boldsymbol{r}}\hat{V}_{\mathrm{NL}}^{\mathrm{KB}}\mathrm{e}^{\mathrm{i}\boldsymbol{A}(t)\cdot\boldsymbol{r}}$. $\boldsymbol{A}(t)$ is nonzero only for the velocity gauge (td_stype=1); otherwise $\boldsymbol{A}(t)=0$. Other nonlocal Hamiltonian terms (e.g., EXX) are not included explicitly.
3474+
* 2: Use the full Hamiltonian to construct the generalized velocity matrix in a nonorthogonal NAO basis: $$\widetilde{v}_{\alpha}=\partial_{\alpha}H+\mathrm{i}HS^{-1}\mathcal{R}_{\alpha}-\mathrm{i}\mathcal{R}_{\alpha}S^{-1}H-HS^{-1}\partial_{\alpha}S.$$ This includes all contributions available in the real-space Hamiltonian matrix when enabled. This method is more general but more expensive.
34743475
default_value: "0"
34753476
unit: ""
34763477
availability: ""

‎source/source_io/module_current/td_current_io_comm.cpp‎

Lines changed: 31 additions & 77 deletions
Original file line numberDiff line numberDiff line change
@@ -15,8 +15,8 @@
1515
#include "td_current_io.h"
1616
#ifdef __EXX
1717
#include "source_lcao/module_operator_lcao/op_exx_lcao.h"
18-
#include "source_lcao/module_ri/Exx_LRI_interface.h"
1918
#include "source_lcao/module_ri/Exx_LRI.h"
19+
#include "source_lcao/module_ri/Exx_LRI_interface.h"
2020
#endif
2121
#ifdef __LCAO
2222
template <typename TR, typename TA>
@@ -163,23 +163,19 @@ void ModuleIO::set_rR_from_hR(const UnitCell& ucell,
163163
const int N2 = iw2n2[iw2];
164164
const int m2 = iw2m2[iw2];
165165

166-
// std::cout<<"L1: "<<L1<<" L2: "<<L2<<" N1: "<<N1<<" N2: "<<N2<<" m1: "<<m1<<" m2:
167-
// "<<m2<<std::endl;
168166
ModuleBase::Vector3<double> tmp_r
169-
= r_calculator
170-
.get_psi_r_psi(tau1 * ucell.lat0, T1, L1, m1, N1, tau2 * ucell.lat0, T2, L2, m2, N2);
167+
= r_calculator.get_psi_r_psi(tau1 * ucell.lat0, T1, L1, m1, N1, tau2 * ucell.lat0, T2, L2, m2, N2);
171168
for (size_t i_alpha = 0; i_alpha != 3; ++i_alpha)
172169
{
173170
hamilt::BaseMatrix<double>* HlocR = rR[i_alpha]->find_matrix(iat1, iat2, r_index);
174171
if (HlocR != nullptr)
175172
{
176-
HlocR->add_element(iw1, iw2, tmp_r[i_alpha]);
173+
// Taoni fix 2026-07-12: HlocR uses local block indices, while row_indexes and col_indexes identify orbitals.
174+
for (int ipol = 0; ipol < npol; ++ipol)
175+
{
176+
HlocR->add_element(iw1l + ipol, iw2l + ipol, tmp_r[i_alpha]);
177+
}
177178
}
178-
// if (i_alpha == 2)
179-
// {
180-
// std::cout << "iw1: " << iw1 << " iw2: " << iw2 << " i_alpha: " << i_alpha
181-
// << " tmp_r: " << tmp_r[i_alpha] << std::endl;
182-
// }
183179
}
184180
}
185181
}
@@ -220,9 +216,8 @@ void ModuleIO::sum_HR(const UnitCell& ucell,
220216
{
221217
atoms_pos[iat] = RI_Util::Vector3_to_array3(ucell.atoms[ucell.iat2it[iat]].tau[ucell.iat2ia[iat]]);
222218
}
223-
const std::array<std::array<double, 3>, 3> latvec = {RI_Util::Vector3_to_array3(ucell.a1),
224-
RI_Util::Vector3_to_array3(ucell.a2),
225-
RI_Util::Vector3_to_array3(ucell.a3)};
219+
const std::array<std::array<double, 3>, 3> latvec
220+
= {RI_Util::Vector3_to_array3(ucell.a1), RI_Util::Vector3_to_array3(ucell.a2), RI_Util::Vector3_to_array3(ucell.a3)};
226221
cell_nearest.init(atoms_pos, latvec, Rs_period);
227222
hamilt::reallocate_hcontainer(ucell.nat, full_hR, Rs_period, &cell_nearest);
228223
}
@@ -349,20 +344,6 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell,
349344
std::complex<double>* r_is_h = new std::complex<double>[pv->nloc];
350345
std::complex<double>* h_is_ps = new std::complex<double>[pv->nloc];
351346

352-
// for (size_t i_alpha = 0; i_alpha != 3; ++i_alpha)
353-
// {
354-
// for (int i = 0; i < hR.size_atom_pairs(); ++i)
355-
// {
356-
// hamilt::AtomPair<double>& tmp = rR[i_alpha]->get_atom_pair(i);
357-
// std::cout<<"cal_velocity_basis_k: "<<rR[i_alpha]->size_atom_pairs()<<" R_size:
358-
// "<<tmp.get_R_size()<<std::endl; for(int ir = 0;ir < tmp.get_R_size(); ++ir )
359-
// {
360-
// const ModuleBase::Vector3<int> r_index = tmp.get_R_index(ir);
361-
// std::cout<<"r_index: "<<r_index<<std::endl;
362-
// }
363-
// }
364-
// }
365-
366347
for (size_t ik = 0; ik != kv.get_nks(); ++ik)
367348
{
368349
// set H(k), S(k)
@@ -456,7 +437,6 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell,
456437
module_rt::folding_partial_HR(ucell, sR, partial_sk, kv.kvec_d[ik], i_alpha, nrow, 1);
457438
}
458439
// 3.3 set r(k)
459-
// std::cout << "set r(k): " << "i_alpha: " << i_alpha << std::endl;
460440
ModuleBase::GlobalFunc::ZEROS(rk, pv->nloc);
461441
// folding_rR(rR[i_alpha], partial_sk, rk, pv, kv.kvec_d[ik], nrow, 1);
462442
if (elecstate::H_TDDFT_pw::stype == 2)
@@ -575,50 +555,14 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell,
575555
1,
576556
pv->desc);
577557
// 4.4 h_is_r will be changed to partial_Hk + IMAG_UNIT * (Hk * Sk_inv * rk)
578-
ScalapackConnector::geadd('N',
579-
nlocal,
580-
nlocal,
581-
one_real,
582-
partial_hk,
583-
1,
584-
1,
585-
pv->desc,
586-
one_imag,
587-
h_is_r,
588-
1,
589-
1,
590-
pv->desc);
558+
ScalapackConnector::geadd('N', nlocal, nlocal, one_real, partial_hk, 1, 1, pv->desc, one_imag, h_is_r, 1, 1, pv->desc);
591559
// 4.5 r_is_h will be changed to h_is_r - IMAG_UNIT * (rk * Sk_inv * Hk)
592-
ScalapackConnector::geadd('N',
593-
nlocal,
594-
nlocal,
595-
one_real,
596-
h_is_r,
597-
1,
598-
1,
599-
pv->desc,
600-
neg_one_imag,
601-
r_is_h,
602-
1,
603-
1,
604-
pv->desc);
560+
ScalapackConnector::geadd('N', nlocal, nlocal, one_real, h_is_r, 1, 1, pv->desc, neg_one_imag, r_is_h, 1, 1, pv->desc);
605561
// 4.6 h_is_ps will be changed to r_is_h - Hk * Sk_inv * partial_Sk
606-
ScalapackConnector::geadd('N',
607-
nlocal,
608-
nlocal,
609-
one_real,
610-
r_is_h,
611-
1,
612-
1,
613-
pv->desc,
614-
neg_one_real,
615-
h_is_ps,
616-
1,
617-
1,
618-
pv->desc);
562+
ScalapackConnector::geadd('N', nlocal, nlocal, one_real, r_is_h, 1, 1, pv->desc, neg_one_real, h_is_ps, 1, 1, pv->desc);
619563
// 5. copy h_is_ps to velocity_basis_k[ik][i_alpha]
620564
BlasConnector::copy(pv->nloc, h_is_ps, 1, velocity_basis_k[ik][i_alpha], 1);
621-
}
565+
}
622566
}
623567

624568
delete[] hk;
@@ -648,7 +592,6 @@ void ModuleIO::cal_velocity_matrix(const psi::Psi<std::complex<double>>* psi,
648592
const char C_char = 'C';
649593
const std::complex<double> one_real = ModuleBase::ONE;
650594
const std::complex<double> zero_complex = ModuleBase::ZERO;
651-
const double zero_double = 0.0;
652595
const int nlocal = PARAM.globalv.nlocal;
653596
const int nbands = PARAM.inp.nbands;
654597
std::complex<double>* vk_c = new std::complex<double>[pv->ncol_bands * pv->nrow_bands]; // local one
@@ -708,9 +651,12 @@ void ModuleIO::cal_velocity_matrix(const psi::Psi<std::complex<double>>* psi,
708651
{
709652
for (int ic = 0; ic < PARAM.inp.nbands; ++ic)
710653
{
711-
const int irc = ic * pv->nrow + ir;
712654
if (pv->in_this_processor(ir, ic))
713655
{
656+
// Taoni fix 2026-07-12: vk_c follows the local block-cyclic layout described by desc_Eij.
657+
const int local_row = pv->global2local_row(ir);
658+
const int local_col = pv->global2local_col(ic);
659+
const int irc = local_col * pv->nrow + local_row;
714660
velocity_k[ik][i_alpha](ir, ic) = vk_c[irc];
715661
}
716662
}
@@ -766,13 +712,22 @@ void ModuleIO::cal_current_comm_k(const UnitCell& ucell,
766712

767713
// sum n and m for current_k
768714
for (size_t ik = 0; ik != kv.get_nks(); ++ik)
715+
{
769716
for (size_t i_alpha = 0; i_alpha != 3; ++i_alpha)
770717
{
771718
for (size_t ib = 0; ib != PARAM.inp.nbands; ++ib)
772719
{
773720
current_k[ik][i_alpha] -= pelec->wg(ik, ib) * velocity_k[ik][i_alpha](ib, ib).real() / 2.0; // for unit
774721
}
775722
}
723+
}
724+
// Taoni fix 2026-07-12: Reduce the current_k values across all MPI processes to get the total current for each k-point.
725+
for (size_t ik = 0; ik != kv.get_nks(); ++ik)
726+
{
727+
Parallel_Reduce::reduce_all(current_k[ik].x);
728+
Parallel_Reduce::reduce_all(current_k[ik].y);
729+
Parallel_Reduce::reduce_all(current_k[ik].z);
730+
}
776731
for (size_t i_alpha = 0; i_alpha < 3; ++i_alpha)
777732
{
778733
delete rR[i_alpha];
@@ -824,14 +779,14 @@ void ModuleIO::write_current(const UnitCell& ucell,
824779
{
825780
if (GlobalV::MY_RANK == 0 && TD_info::out_current_k)
826781
{
827-
std::string filename = PARAM.globalv.global_out_dir + "currents" + std::to_string(is) + "k"
828-
+ std::to_string(ik) + "comm.txt";
782+
std::string filename
783+
= PARAM.globalv.global_out_dir + "currents" + std::to_string(is) + "k" + std::to_string(ik) + "comm.txt";
829784
std::ofstream fout;
830785
fout.open(filename, std::ios::app);
831786
fout << std::setprecision(16);
832787
fout << std::scientific;
833-
fout << istep << " " << current_k[ik][0] / omega << " " << current_k[ik][1] / omega << " "
834-
<< current_k[ik][2] / omega << std::endl;
788+
fout << istep << " " << current_k[ik][0] / omega << " " << current_k[ik][1] / omega << " " << current_k[ik][2] / omega
789+
<< std::endl;
835790
fout.close();
836791
}
837792
}
@@ -853,8 +808,7 @@ void ModuleIO::write_current(const UnitCell& ucell,
853808
fout.open(filename, std::ios::app);
854809
fout << std::setprecision(16);
855810
fout << std::scientific;
856-
fout << istep << " " << current_total[0] / omega << " " << current_total[1] / omega << " "
857-
<< current_total[2] / omega << std::endl;
811+
fout << istep << " " << current_total[0] / omega << " " << current_total[1] / omega << " " << current_total[2] / omega << std::endl;
858812
fout.close();
859813
}
860814

‎source/source_io/module_parameter/read_input_item_output.cpp‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1483,9 +1483,10 @@ In molecular dynamics calculations, the output frequency is controlled by out_fr
14831483
item.annotation = "output current or not";
14841484
item.category = "RT-TDDFT: Real-Time Time-Dependent Density Functional Theory";
14851485
item.type = "Integer";
1486-
item.description = R"(* 0: Do not output current.
1487-
* 1: Output current using the two-center integral, faster.
1488-
* 2: Output current using the matrix commutation, more precise.)";
1486+
item.description = R"(Controls the current-density output method for LCAO RT-TDDFT.
1487+
* 0: Do not output current.
1488+
* 1: Explicitly construct the velocity operator from the momentum, vector-potential, and KB nonlocal-pseudopotential terms using two-center integral / spherical grid integral: $$\hat{v}_{\alpha}=-\mathrm{i}\nabla_{\alpha}+A_{\alpha}(t)+\mathrm{i}\left[\widetilde{V}_{\mathrm{NL}}^{\mathrm{KB}},r_{\alpha}\right],$$ where $\widetilde{V}_{\mathrm{NL}}^{\mathrm{KB}}=\mathrm{e}^{-\mathrm{i}\boldsymbol{A}(t)\cdot\boldsymbol{r}}\hat{V}_{\mathrm{NL}}^{\mathrm{KB}}\mathrm{e}^{\mathrm{i}\boldsymbol{A}(t)\cdot\boldsymbol{r}}$. $\boldsymbol{A}(t)$ is nonzero only for the velocity gauge (td_stype=1); otherwise $\boldsymbol{A}(t)=0$. Other nonlocal Hamiltonian terms (e.g., EXX) are not included explicitly.
1489+
* 2: Use the full Hamiltonian to construct the generalized velocity matrix in a nonorthogonal NAO basis: $$\widetilde{v}_{\alpha}=\partial_{\alpha}H+\mathrm{i}HS^{-1}\mathcal{R}_{\alpha}-\mathrm{i}\mathcal{R}_{\alpha}S^{-1}H-HS^{-1}\partial_{\alpha}S.$$ This includes all contributions available in the real-space Hamiltonian matrix when enabled. This method is more general but more expensive.)";
14891490
item.default_value = "0";
14901491
item.unit = "";
14911492
item.availability = "";

0 commit comments

Comments
 (0)