[Pw_forum] Error in divide_class ( with D_2 symmetry)

[Dal Corso Andrea]

On Wed, 2007-08-15 at 14:35 -0400, Tao Sun wrote:
> 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
>  -----------------------------------------------------  
> 

You are right, divide_class is wrong for D_2. Thank you for reporting
this. Actually in your system two of the three C_2 axes are not parallel
to x or y, so even correcting the bug will not solve the problem.

One possibility is to take C_2 parallel to the axis of 
the first rotation, C_2' parallel to the axis of the 
second rotation, and C_2'' parallel to the axis of the third rotation.
This can be done substituting the loop over classes with:

   DO iclass=2,nclass
      which_irr(iclass)=iclass
   ENDDO
 
For divide_class_so we can do the same thing. The loop over the
classes becomes:


   first=.true.
   DO iclass=2,nclass
      ts=tipo_sym(smat(1,1,elem(1,iclass)))
      IF (ts==1) THEN
         which_irr(iclass)=2
         first=.false.
      ELSE
         if (first) then
            which_irr(iclass)=iclass+1
         else
            which_irr(iclass)=iclass
         endif
      END IF
   ENDDO


Hope this helps. 

Andrea





> 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
> 
> 
> _______________________________________________
> Pw_forum mailing list
> Pw_forum at pwscf.org
> http://www.democritos.it/mailman/listinfo/pw_forum
-- 
Andrea Dal Corso                    Tel. 0039-040-3787428
SISSA, Via Beirut 2/4               Fax. 0039-040-3787528
34014 Trieste (Italy)               e-mail: dalcorso at sissa.it