[Pw_forum] Error in divide_class ( with D_2 symmetry)

[Tao Sun]

Dear Andrea,

Thanks. After your modification, the code ran smoothly
and yielded the following:

   Mode symmetry, D_2  (222)  point group:

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

I guess I need to look at the eigenvectors with respect to rotation
axis to judge whether these symmetry assignments are correct. Anyway
they are not very relevant to our current project and I hope you 
won't mind I didn't check this.

Cheers.

Tao  

On Thu, 16 Aug 2007, Dal Corso Andrea wrote:

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