I traced the bug to this line:
call MPI_BCAST( input%preconditioning_param, 1, MPI_DOUBLE_PRECISION, 0, mpi%mpi_comm, ierr )
Rank 0 executes this line fine, all others fail. I don't understand why. It seems ierr is an array and is used as a scalar but this seems not to be the issue.
I did some changes to the line after this:
CALL ylmnorm_init(max(atoms%lmaxd, 2*hybrid%lexp))
But I don't really see how this could the line before.