Error in anisotropic calculation

[mdogan]

Dear H. Lee,

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:

Code: Select all

ibin = NINT(lambda_eph / dbin) + 1

This means lambda_eph / dbin = -333 (approximately). So, either lambda_eph or dbin is negative. lambda_eph seems to be set in lines 511-517:

Code: Select all

lambda_eph = zero
DO imode = 1, nmodes
IF (wf(imode, iq0) > eps_acustic) THEN
    lambda_eph = lambda_eph + g2(ik, iq, ibnd, jbnd, imode) / wf(imode,iq0)
  ENDIF
ENDDO
lambda_eph = 2.d0 * lambda_eph * dosef

As you said, if dose, wf and g2 are nonnegative, lambda_eph should be as well.

dbin should be defined in lines 493-494:

Code: Select all

nbin = NINT(1.1d0 * MAXVAL(lambda_max(:)) / eps2) + 1
dbin = 1.1d0 * MAXVAL(lambda_max(:)) / DBLE(nbin)

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!

Best,
Mehmet

[hlee]

Dear Mehmet:

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.)

(before)

Code: Select all

    IF (iverbosity == 2) THEN
      nbin = NINT(1.1d0 * MAXVAL(lambda_max(:)) / eps2) + 1
      dbin = 1.1d0 * MAXVAL(lambda_max(:)) / DBLE(nbin)
      ALLOCATE(lambda_pairs(nbin), STAT = ierr)
      IF (ierr /= 0) CALL errore('evaluate_a2f_lambda', 'Error allocating lambda_pairs', 1)
      lambda_pairs(:) = zero
    ENDIF

(after)

Code: Select all

    IF (iverbosity == 2) THEN
      nbin = NINT(1.1d0 * MAXVAL(lambda_max(:)) / eps2) + 1
      dbin = 1.1d0 * MAXVAL(lambda_max(:)) / DBLE(nbin)
      ALLOCATE(lambda_pairs(nbin), STAT = ierr)
      IF (ierr /= 0) CALL errore('evaluate_a2f_lambda', 'Error allocating lambda_pairs', 1)
      lambda_pairs(:) = zero
      IF (nbin < 0) WRITE(stdout, '(a)') 'WARNING: nbin is negative!'
      IF (dbin < 0.d0) WRITE(stdout, '(a)') 'WARNING: dbin is negative!'
    ENDIF

Third, do the same thing for g2, wf, and dosef:

(before)

Code: Select all

                IF (iverbosity == 2) THEN
                  ibin = NINT(lambda_eph / dbin) + 1
                  weight2 =  w0g(ibnd, ik) * w0g(jbnd,ixkqf(ik, iq0))
                  lambda_pairs(ibin) = lambda_pairs(ibin) + weight2
                ENDIF

(after)

Code: Select all

                 IF (iverbosity == 2) THEN
                  ibin = NINT(lambda_eph / dbin) + 1
                  IF (ibin < 0) THEN
		     WRITE(stdout, '(a)') 'WARNING: ibin is negative!'
	             WRITE(stdout, *) 'dosef=', dosef
	             WRITE(stdout, *) 'g2=', g2(ik, iq, ibnd, jbnd, 1:nmodes)
	             WRITE(stdout, *) 'wf=', wf(1:nmodes,iq0)
                  ENDIF
                  weight2 =  w0g(ibnd, ik) * w0g(jbnd,ixkqf(ik, iq0))
                  lambda_pairs(ibin) = lambda_pairs(ibin) + weight2
                ENDIF

Lastly, you need to check the main output of epw.out as well as out.* files.

Sincerely,

H. Lee

[mdogan]

Dear H. Lee,

Thank you very much for your detailed guidance in debugging. Here is what happened. (Refer to the Dropbox folder [https://www.dropbox.com/sh/a12q7198eiq7 ... b06Ha?dl=0] for the files.)

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.

Best,
Mehmet

[hlee]

Dear Mehmet:

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.

Sincerely,

H. Lee