[Pw_forum] Bug in QE 5.0 with ibrav=0 and using a instead of celldm

[Giovanni Pizzi]

Dear all,
I believe there is a bug in QE 5.0 (pw.x) that occurs when using ibrav=0 
and defining 'a' instead of 'celldm(1)'.
In this case, celldm is not set and remains equal to zero. (This 
behavior was different in 4.3.2, where instead the celldm(1) value was 
set from the 'a' value.)

This doesn't seem to create problems in the pw.x calculation.
However, in the matdyn files generated by ph.x, the celldm(1) is set to 
zero.
When the matdyn files are then given to q2r.x, the code crashes, 
complaining of a wrong celldm value.

In particular this is what happens in q2r.x:
1. the first q point is read (matdyn1) using the routine 
read_dyn_from_file; celldm(1) is set to 0
2. the code calls the routine latgen: this routine sets celldm(1) to the 
a value different from zero (moreover this value may be wrong if 'a' in 
the input file was not set to the length of the first lattice vector)
3. the code reads the second q point (matdyn2) and compares the celldm 
value read from file (which is zero) with the old value (which was set 
to a value different from zero at step 2): being different, the code stops.

I think that to correct it, one should replace the following lines 
(cell_base.f90, lines 152 and following)
      IF( celldm( 1 ) /= 0.0_DP ) THEN
         alat = celldm( 1 )
      ELSE IF ( a /= 0.0_dp ) THEN
         alat = a / bohr_radius_angs
      ELSE
         alat = SQRT ( at(1,1)**2+at(2,1)**2+at(3,1)**2 )
      END IF
      at(:,:) = at(:,:) / alat

with
      IF( celldm( 1 ) /= 0.0_DP ) THEN
         alat = celldm( 1 )
      ELSE IF ( a /= 0.0_dp ) THEN
         alat = a / bohr_radius_angs
     celldm( 1 ) = alat
      ELSE
         alat = SQRT ( at(1,1)**2+at(2,1)**2+at(3,1)**2 )
     celldm( 1 ) = alat
      END IF
      at(:,:) = at(:,:) / alat

Let me know what you think.

Best,
Giovanni Pizzi



-- 
Giovanni Pizzi
Post-doctoral Research Scientist
EPFL STI IMX THEOS
MXC 319 (Bâtiment MXC)
Station 12
CH-1015 Lausanne (Switzerland)
Phone: +41 21 69 31159