Thank you for pointing these out. You may check my input file (epw4.in) in the Dropbox folder, I assume I'm only using standard inputs. From the source code, here's what I have so far:

The line that sets ibin to -332 should be line 520, which is:

Combining these two expressions, I get dbin to be approximately equal to (lambda_max*eps2)/(lambda_max+eps2). From constants_epw.f90 it looks like eps2 = 1.0E-2. If my math is correct, for dbin to be negative, we'd have to have (approximately) -1.0E-2 < lambda_max < 0. Of course, since lambdas should be nonnegative, this also doesn't make sense.

It seems to me there's no obvious variable where the problem may creep in. Do you agree with that? In order to pin down where the negative numbers come from, I should print some of these variables during the run. But I'm not experienced with Fortran code editing. Could you give me an example of such a line to print a variable to the output file? For instance, if I wanted to print ibin after line 520, what line would I need to add? Thank you very much for all your help!

You need to find the origin in the following ways:

First, you need to print out the values of relevant variables on each MPI task.
To do this, please proceed as follows:
(1) Change the line of "MANUAL_DFLAGS =" in make.inc in the QE root directory with that of "MANUAL_DFLAGS = -DDEBUG".
(2) Issue "make clean".
(3) Build again the code; make ... (you should not issue "./configure ..." again).

Second, check the values of nbin and dbin.
(There might be some typos in the following modifications.)

The output file looks the same (epw4.out3) but 127 out.0_* files were also created. Only one of them printed a warning message (out.0_32). The other 126 were identical to out.0_1.

As you can see in out.0_32, there are 3 negative elements of g2. Does everything else make sense? If so, the problem seems to be these negative elements of g2. Where do you think they might originate, and how can I pin that down? I greatly appreciate your help.

I think that something wrong happens in your case.

We need to get more information and we need the following modifications in the SUBROUTINE write_ephmat in io_eliashberg.f90.
(There might be some typos below.)

(before)
IF (ixkf(lower_bnd + ik - 1) > 0) THEN
IF (ixkqf(ixkf(lower_bnd + ik - 1), iq) > 0) THEN
DO imode = 1, nmodes ! phonon modes
wq = wf(imode, iq)
inv_wq = 1.0 / (two * wq)
!
DO ibnd = 1, nbndfst
IF (ABS(ekfs(ibnd, ixkf(lower_bnd + ik - 1)) - ef0) < fsthick) THEN
DO jbnd = 1, nbndfst
IF (ABS(ekfs(jbnd, ixkqf(ixkf(lower_bnd + ik - 1), iq)) - ef0) < fsthick) THEN
!
! here we take into account the zero-point DSQRT(hbar/2M\omega)
! with hbar = 1 and M already contained in the eigenmodes
! g2 is Ry^2, wkf must already account for the spin factor
!
IF (shortrange .AND. (ABS(xqf(1, iq)) > eps8 .OR. ABS(xqf(2, iq)) > eps8 &
.OR. ABS(xqf(3, iq)) > eps8 )) THEN
! SP: The abs has to be removed. Indeed the epf17 can be a pure imaginary
! number, in which case its square will be a negative number.
g2 = REAL((epf17(jbnd, ibnd, imode, ik)**two) * inv_wq)
ELSE
g2 = ABS(epf17(jbnd, ibnd, imode, ik))**two * inv_wq
ENDIF
ind(my_pool_id + 1) = ind(my_pool_id + 1) + 1
tmp_g2(ind(my_pool_id + 1)) = g2
ENDIF
ENDDO ! jbnd
ENDIF
ENDDO ! ibnd
ENDDO ! imode
ENDIF ! ixkqf
ENDIF ! ixkf

(after)
IF (ixkf(lower_bnd + ik - 1) > 0) THEN
IF (ixkqf(ixkf(lower_bnd + ik - 1), iq) > 0) THEN
DO imode = 1, nmodes ! phonon modes
wq = wf(imode, iq)
inv_wq = 1.0 / (two * wq)
!
DO ibnd = 1, nbndfst
IF (ABS(ekfs(ibnd, ixkf(lower_bnd + ik - 1)) - ef0) < fsthick) THEN
DO jbnd = 1, nbndfst
IF (ABS(ekfs(jbnd, ixkqf(ixkf(lower_bnd + ik - 1), iq)) - ef0) < fsthick) THEN
!
! here we take into account the zero-point DSQRT(hbar/2M\omega)
! with hbar = 1 and M already contained in the eigenmodes
! g2 is Ry^2, wkf must already account for the spin factor
!
IF (shortrange .AND. (ABS(xqf(1, iq)) > eps8 .OR. ABS(xqf(2, iq)) > eps8 &
.OR. ABS(xqf(3, iq)) > eps8 )) THEN
! SP: The abs has to be removed. Indeed the epf17 can be a pure imaginary
! number, in which case its square will be a negative number.
g2 = REAL((epf17(jbnd, ibnd, imode, ik)**two) * inv_wq)
ELSE
g2 = ABS(epf17(jbnd, ibnd, imode, ik))**two * inv_wq
ENDIF
!!!!!
IF (g2 < 0.d0) THEN
WRITE(stdout, *) 'warning: g2 is negative. g2=', g2
WRITE(stdout, *) 'inv_wq=', inv_wq
WRITE(stdout, *) 'iq=', iq
WRITE(stdout, *) 'ik=', ik
WRITE(stdout, *) 'imode=', imode
WRITE(stdout, *) 'ibnd=', ibnd
WRITE(stdout, *) 'jbnd=', jbnd
IF (shortrange) WRITE(stdout, *) 'shortrange is .true.'
ENDIF
!!!!!
ind(my_pool_id + 1) = ind(my_pool_id + 1) + 1
tmp_g2(ind(my_pool_id + 1)) = g2
ENDIF
ENDDO ! jbnd
ENDIF
ENDDO ! ibnd
ENDDO ! imode
ENDIF ! ixkqf
ENDIF ! ixkf

You just need to focus on the region sandwiched by "!!!!!" and you need to rerun the evaluation of e-ph vertex on fine grids.