Length gague rt-tddft terminated during current calculation

[ESROAMER]

Describe the bug

Running length gague rt-tddft with out_current 1 will terminated the program at the first time step during current calculation. Here is an example.
ex.zip
The program runs normally when using 1 thread and 20 processes, but it crashes when using 1 thread and 28 processes. The same problem does not occur in the velocity gauge.

Expected behavior

No response

To Reproduce

No response

Environment

No response

Additional Context

No response

Task list for Issue attackers (only for developers)

• Verify the issue is not a duplicate.
• Describe the bug.
• Steps to reproduce.
• Expected behavior.
• Error message.
• Environment details.
• Additional context.
• Assign a priority level (low, medium, high, urgent).
• Assign the issue to a team member.
• Label the issue with relevant tags.
• Identify possible related issues.
• Create a unit test or automated test to reproduce the bug (if applicable).
• Fix the bug.
• Test the fix.
• Update documentation (if necessary).
• Close the issue and inform the reporter (if applicable).

[AsTonyshment]

After some painful debugging, I found that the issue might be in line 288 of the td_current.cpp file. There is no issue under the velocity gauge, but when using the length gauge with 28 MPI processes in parallel (actually for any np > 20), the call this->current_term[i]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]) returns nullptr, which does not happen under the velocity gauge. I added debugging information before this line, and the output is:

Loop i = 0
Error: find_matrix() returned nullptr for iat1 = 1, iat2 = 3, R_vector=(0, 0, 0), thread_id = 0

But currently I have no idea why this is happening.

[AsTonyshment]

Adding a non-null pointer check can bypass this issue, i.e., make the following modifications at line 288 and line 420 in the td_current.cpp file:

Change

abacus-develop/source/module_hamilt_lcao/module_tddft/td_current.cpp

Line 288 in 70b7f58

tmp_c[i] = this->current_term[i]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2])->get_pointer();

to

auto matrix_ptr = this->current_term[i]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]);
if (matrix_ptr != nullptr)
    tmp_c[i] = matrix_ptr->get_pointer();

Change

abacus-develop/source/module_hamilt_lcao/module_tddft/td_current.cpp

Line 420 in 70b7f58

tmp_c[i] = this->current_term[i]->find_matrix(iat1, iat2, R_index2)->get_pointer();

to

auto matrix_ptr = this->current_term[i]->find_matrix(iat1, iat2, R_index2);
if (matrix_ptr != nullptr)
    tmp_c[i] = matrix_ptr->get_pointer();

However, I'm not sure if this approach is valid. What do you think @ESROAMER?

[ESROAMER]

Adding a non-null pointer check can bypass this issue, i.e., make the following modifications at line 288 and line 420 in the td_current.cpp file:

Change

abacus-develop/source/module_hamilt_lcao/module_tddft/td_current.cpp

Line 288 in 70b7f58

tmp_c[i] = this->current_term[i]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2])->get_pointer();

to
auto matrix_ptr = this->current_term[i]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]);
if (matrix_ptr != nullptr)
tmp_c[i] = matrix_ptr->get_pointer();
Change

abacus-develop/source/module_hamilt_lcao/module_tddft/td_current.cpp

Line 420 in 70b7f58

tmp_c[i] = this->current_term[i]->find_matrix(iat1, iat2, R_index2)->get_pointer();

to
auto matrix_ptr = this->current_term[i]->find_matrix(iat1, iat2, R_index2);
if (matrix_ptr != nullptr)
tmp_c[i] = matrix_ptr->get_pointer();
However, I'm not sure if this approach is valid. What do you think @ESROAMER?

I suggest testing an example of a molecule with an applied electric field for a few steps. If the multi-process results are consistent with the single-core results under this scheme, then I think it should be valid.