[Pw_forum] Error in divide_class ( with D_2 symmetry)

[Tao Sun]

Hi,

When I tried to compute phonon at q=(-0.01000100010001, -1.000100010001,
0.5) in a shear-strained fcc primitive cell

CELL_PARAMETERS (alat)
 -0.5 0.005  0.5
 -0.005 0.5  0.5
 -0.505 0.505 0.0

using espresso-3.2.2

I met the following error message:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     from divide_class : error #         1
     wrong axis
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     stopping ...

This error message comes from PW/divide_class.f90, in the D_2
symmetry part. between line 150-169:
 
-----------------------------------------------------
ELSEIF (code_group==8) THEN
!
!  D_2
!
!   versor find the rotation axis.
!   is_axis(ax,ind) is true if ax the axis given by ind
!   (1 asse x, 2 asse y, 3 asse z).
!
   DO iclass=2,nclass
      CALL versor(smat(1,1,elem(1,iclass)),ax)
      IF (is_axis(ax,3)) THEN
         which_irr(iclass)=2
      ELSEIF (is_axis(ax,2)) THEN
         which_irr(iclass)=3
      ELSEIF (is_axis(ax,2)) THEN
         which_irr(iclass)=4
      ELSE
         CALL errore('divide_class','wrong axis',1)
      ENDIF
   ENDDO
 -----------------------------------------------------  

The duplicated condition (is_axis(ax,2)) looks suspicious,
but when I changed one of them into is_axis(ax,1), the
error is still there.

If I remove the CALL errore(... sentence, the code can 
compute the eigen-values, except gives weird symmetry
label to each mode, as:

---------------------------------------------------------
     Mode symmetry, D_2  (222)  point group:

     omega(  1 -  1) =     81.09814  [cm-1]   -->   ?
     omega(  1 -  1) =     81.09814  [cm-1]   --> A
     omega(  1 -  1) =     81.09814  [cm-1]   -->   ?
     omega(  1 -  1) =     81.09814  [cm-1]   --> B_1
     omega(  2 -  2) =    111.55732  [cm-1]   -->   ?
     omega(  2 -  2) =    111.55732  [cm-1]   --> B_2
     omega(  2 -  2) =    111.55732  [cm-1]   -->   ?
     omega(  2 -  2) =    111.55732  [cm-1]   --> B_3
     omega(  3 -  3) =    124.59709  [cm-1]   -->   ?
     omega(  3 -  3) =    124.59709  [cm-1]   --> B_2
     omega(  3 -  3) =    124.59709  [cm-1]   -->   ?
     omega(  3 -  3) =    124.59709  [cm-1]   --> B_3
------------------------------------------------------

This error is independent of the system, we observe the
same behavior in two different fcc crystals(MgO and Pt).

BTW, similar problem exists in PW/divide_class_so.f90
while q points of other symmetry type(e.g. D_2h) seems OK.

It will be great if somebody can help us fix the problem.

Cheers.

Tao Sun

Graduate Student
Dept. of Physics & Astronomy
SUNY at Stony Brook