[Octopus-users] A bug in paramagnetic_current?

[mouseyan]

Hej, Developers,

    I found that when using time evolution enabling parallel in
states, the output current density is different when changing the
number of nodes used for states.

The original code in subroutine states_paramagnetic_current is:

#if defined(HAVE_MPI)
    if(st%parallel_in_states) then
      ALLOCATE(red(NP_PART, NDIM, st%d%nspin), NP_PART*NDIM*st%d%nspin)
      call MPI_Allreduce(jp(1, 1, 1), red(1, 1, 1), NP*NDIM*st%d%nspin,       &
        MPI_FLOAT, MPI_SUM, st%mpi_grp%comm, mpi_err)
      jp = red
      deallocate(red)
    end if
#endif

I have modified it, shown in the following.

#if defined(HAVE_MPI)
    if(st%parallel_in_states) then
      ALLOCATE(red(1:NP),NP)
      do i = 1, st%d%nspin
         do k = 1, NDIM
           call MPI_Allreduce(jp(1, k, i), red(1), NP,       &
                MPI_FLOAT, MPI_SUM, st%mpi_grp%comm, mpi_err)
           call lalg_copy(NP, red, jp(:, k, i))
         enddo
      enddo
      deallocate(red)
    end if
#endif

The problem was solved and the continuity equation is satisfied.
Please check that. I hope it will help.

All the best