TABLE OF CONTENTS


ABINIT/symspgr [ Functions ]

[ Top ] [ Functions ]

NAME

 symspgr

FUNCTION

 Using the type of each symmetry operation (found in symplanes.f and symaxes.f):
 proper symmetries 1,2,2_1,3,3_1,3_2,4,4_1,4_2,4_3,6,6_1,...6_5
 improper symmetries -1,m,a,b,c,d,n,g,-3,-4,-6 ,
 build an array with the number of such operations. then, call symlist.f to identify the space group.
 The identification is not unambiguous still ...

COPYRIGHT

 Copyright (C) 2000-2018 ABINIT group (RC, XG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

 bravais(11): bravais(1)=iholohedry
              bravais(2)=center
              bravais(3:11)=coordinates of rprimd in the axes
              of the conventional bravais lattice (*2 if center/=0)
 nsym=actual number of symmetries
 symrel(3,3,nsym)= nsym symmetry operations in real space in terms
   of primitive translations
 tnons(3,nsym)=nonsymmorphic translations for each symmetry (would
   be 0 0 0 each for a symmorphic space group)

OUTPUT

 spgroup=symmetry space group number

NOTES

 It is assumed that the symmetry operations will be entered in the
 symrel tnons arrays, for the PRIMITIVE cell. The matrix of transformation
 from the primitive cell to the conventional cell is described
 in the array "bravais" (see symlatt.F90).
 The present routine first make the transformation from the
 primitive coordinates to the conventional ones, then eventually
 generate additional symmetries, taking into account the
 centering translations.
 Then, the order and determinant of each symmetry operation
 is determined.

 For proper symmetries (rotations), the
 associated translation is also determined.
 However, left or right handed screw rotations are
 not (presently) distinguished, and will be attributed equally
 to left or right.

 For the detailed description of the labelling of the axes,
 see symaxes.f and symplanes.f

PARENTS

      symanal

CHILDREN

      spgdata,symcharac,symdet,symlist_bcc,symlist_fcc,symlist_others
      symlist_prim,symrelrot,wrtout,xred2xcart

SOURCE

 64 #if defined HAVE_CONFIG_H
 65 #include "config.h"
 66 #endif
 67 
 68 #include "abi_common.h"
 69 
 70 subroutine symspgr(bravais,nsym,spgroup,symrel,tnons,tolsym)
 71 
 72  use defs_basis
 73  use m_errors
 74  use m_profiling_abi
 75  use m_numeric_tools, only : OPERATOR(.x.)
 76 
 77  use m_geometry,   only : xred2xcart
 78 
 79 !This section has been created automatically by the script Abilint (TD).
 80 !Do not modify the following lines by hand.
 81 #undef ABI_FUNC
 82 #define ABI_FUNC 'symspgr'
 83  use interfaces_14_hidewrite
 84  use interfaces_41_geometry, except_this_one => symspgr
 85 !End of the abilint section
 86 
 87  implicit none
 88 
 89 !Arguments ------------------------------------
 90 !scalars
 91  integer,intent(in) :: nsym
 92  integer,intent(out) :: spgroup
 93  real(dp),intent(in) :: tolsym
 94 !arrays
 95  integer,intent(in) :: bravais(11),symrel(3,3,nsym)
 96  real(dp),intent(inout) :: tnons(3,nsym)
 97 
 98 !Local variables-------------------------------
 99 !scalars
100 ! logical,parameter :: verbose=.FALSE.
101  integer :: additional_info,brvltt,center,direction=0,found,iholohedry,ii
102  integer :: ishift,isym,jj,nshift,nsymconv,spgaxor,spgorig,sporder
103  character(len=1) :: brvsb
104  character(len=15) :: intsb,ptintsb,ptschsb,schsb
105  character(len=35) :: intsbl
106  character(len=500) :: message
107  character(len=128) :: label
108 !arrays
109  integer :: n_axes(31),n_axest(31),prime(5),test_direction(3),symrel_uni(3,3)
110  integer :: uniaxis(3),uniaxis_try(3)
111  integer,allocatable :: determinant(:),symrelconv(:,:,:),t_axes(:)
112  real(dp) :: axes(3,3),rprimdconv(3,3),trialt(3),vect(3,3)
113  real(dp),allocatable :: shift(:,:),tnonsconv(:,:)
114 
115 !**************************************************************************
116 
117  DBG_ENTER("COLL")
118 
119 !Initialize brvltt, from bravais(2) and bravais(1)
120  center=bravais(2)
121  iholohedry=bravais(1)
122  brvltt=1
123  if(center==-1)brvltt=2  ! Inner centering
124  if(center==-3)brvltt=3  ! Face centering
125  if(center==1)brvltt=5  ! A-Face centering
126  if(center==2)brvltt=6  ! B-Face centering
127  if(center==3)brvltt=4  ! C-Face centering
128  if(iholohedry==5)brvltt=7  ! Rhombohedral
129 
130 !Produce the symmetry operations, in the axis of the conventional cell
131  nsymconv=nsym
132  if(center/=0)nsymconv=2*nsymconv
133  if(center==-3)nsymconv=4*nsym
134  ABI_ALLOCATE(symrelconv,(3,3,nsymconv))
135  ABI_ALLOCATE(tnonsconv,(3,nsymconv))
136 
137 !Produce symrel and tnons in conventional axes,
138 !name them symrelconv and tnonsconv
139  rprimdconv(:,1)=bravais(3:5)
140  rprimdconv(:,2)=bravais(6:8)
141  rprimdconv(:,3)=bravais(9:11)
142 
143  if(center/=0)rprimdconv(:,:)=rprimdconv(:,:)*half
144 
145  axes(:,:)=zero
146  axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one
147  symrelconv(:,:,1:nsym)=symrel(:,:,1:nsym)
148 !Note that the number of symmetry operations is still nsym
149  call symrelrot(nsym,rprimdconv,axes,symrelconv,tolsym)
150  call xred2xcart(nsym,rprimdconv,tnonsconv,tnons)
151 !Gives the associated translation, with components in the
152 !interval ]-0.5,0.5] .
153  tnonsconv(:,1:nsym)=tnonsconv(:,1:nsym)-nint(tnonsconv(:,1:nsym)-tol6)
154 
155 !If the Bravais lattice is centered, duplicate or quadruplicate
156 !the number of symmetry operations, using the Bravais
157 !lattice shifts
158  nshift=1
159  if(center/=0)nshift=2
160  if(center==-3)nshift=4
161  ABI_ALLOCATE(shift,(3,nshift))
162  shift(:,1)=zero
163  if(center/=0 .and. center/=-3)then
164    shift(:,2)=half
165    if(center==1)shift(1,2)=zero
166    if(center==2)shift(2,2)=zero
167    if(center==3)shift(3,2)=zero
168  else if(center==-3)then
169    shift(:,2)=half ; shift(1,2)=zero
170    shift(:,3)=half ; shift(2,3)=zero
171    shift(:,4)=half ; shift(3,4)=zero
172  end if ! center/=0 or -3
173  if(nshift/=1)then
174    do ishift=2,nshift
175      symrelconv(:,:,(ishift-1)*nsym+1:ishift*nsym)=symrelconv(:,:,1:nsym)
176      do isym=1,nsym
177        tnonsconv(:,(ishift-1)*nsym+isym)=tnonsconv(:,isym)+shift(:,ishift)
178      end do
179    end do ! ishift
180  end if ! nshift/=1
181 
182 !At this stage, all the symmetry operations are available,
183 !expressed in the conventional axis, and also include
184 !the Bravais lattive translations, and associated operations...
185 
186  n_axes(:)=0
187 
188  ABI_ALLOCATE(determinant,(nsymconv))
189 
190 !Get the determinant
191  call symdet(determinant,nsymconv,symrelconv)
192 
193 !Get the order of each the symmetry operation, as well as the maximal order
194 !Also, examine whether each symmetry operation is the inversion, or a root
195 !of the inversion (like -3)
196 !Decide which kind of point symmetry operation it is
197 !Finally assign tnonsconv order and decide the space symmetry operation
198 
199  ABI_ALLOCATE(t_axes,(nsymconv))
200 
201  do isym=1,nsymconv
202 
203    call symcharac(center, determinant(isym), iholohedry, isym, label, &
204    symrelconv(:,:,isym), tnonsconv(:,isym), t_axes(isym))
205    if (t_axes(isym) == -1) then
206      write(message, '(a,a,i3,a,3(a,3i4,a),a,3es22.12,a,a,3es22.12)' )ch10,&
207 &     ' symspgr : problem with isym=',isym,ch10,&
208 &     '  symrelconv(:,1,isym)=',symrelconv(:,1,isym),ch10,&
209 &     '  symrelconv(:,2,isym)=',symrelconv(:,2,isym),ch10,&
210 &     '  symrelconv(:,3,isym)=',symrelconv(:,3,isym),ch10,&
211 &     '  tnonsconv(:,isym)=',tnonsconv(:,isym),ch10,&
212 &     '  trialt(:)=',trialt(:)
213      call wrtout(std_out,message,'COLL')
214      write(message, '(a,i4,2a)' )&
215 &     'The space symmetry operation number',isym,ch10,&
216 &     'is not a (translated) root of unity'
217      MSG_BUG(message)
218    else if (t_axes(isym) == -2) then
219      write(message, '(a,i0,a)' )'The symmetry operation number ',isym,' is not a root of unity'
220      MSG_BUG(message)
221    end if
222 
223    n_axes(t_axes(isym))=n_axes(t_axes(isym))+1
224 
225  end do ! isym=1,nsymconv
226 
227  if (sum(n_axes)-nsymconv/=0) then
228    write(message, '(7a)' )&
229 &   'Not all the symmetries have been recognized. ',ch10,&
230 &   'This might be due either to an error in the input file',ch10,&
231 &   'or to a BUG in ABINIT',ch10,&
232 &   'Please contact the ABINIT group.'
233    MSG_WARNING(message)
234  end if
235 
236 !DEBUG
237 !write(std_out,*)' symspgr : brvltt,nsymconv=',brvltt,nsymconv
238 !write(std_out,*)' n_axes(1:10)=',n_axes(1:10)
239 !write(std_out,*)' n_axes(11:20)=',n_axes(11:20)
240 !write(std_out,*)' n_axes(21:31)=',n_axes(21:31)
241 !ENDDEBUG
242 
243 !Treat cases in which the space group cannot be identified on the
244 !basis of n_axes one need additional informations
245  if(brvltt==1)then
246 !  If the bravais lattice is primitive
247    if(nsymconv==4)then
248      n_axest=(/0,0,0,0,0,0,0,1,1,0,  0,0,0,0,0,2,0,0,0,0,  0,0,0,0,0,0,0,0,0,0,0/)
249      if(sum((n_axes-n_axest)**2)==0)then    ! Spgroup 27 (Pcc2) or 32 (Pba2)
250        write(std_out,*)' symspgr: 27 or 32'
251        additional_info=2
252 !      Select binary axis
253        do isym=1,nsymconv
254          if(t_axes(isym)==8)then
255 !          Find direction of binary axis
256            if(symrelconv(1,1,isym)==1)direction=1
257            if(symrelconv(2,2,isym)==1)direction=2
258            if(symrelconv(3,3,isym)==1)direction=3
259          end if
260        end do
261 !      Examine the projection of the translation vector of the a, b or c mirror planes
262 !      onto the binary axis
263        do isym=1,nsymconv
264          if(t_axes(isym)==16)then
265            if(abs(tnonsconv(direction,isym))>tol8)additional_info=1
266          end if
267        end do
268      end if
269    else if(nsymconv==8)then
270      n_axest=(/0,0,0,0,1,0,0,1,1,0,  0,0,0,0,1,2,0,0,0,2,  0,0,0,0,0,0,0,0,0,0,0/)
271      if(sum((n_axes-n_axest)**2)==0)then    ! Spgroup 55 (Pbam) or 57 (Pbcm)
272        write(std_out,*)' symspgr: 55 or 57'
273        additional_info=1
274 !      Select mirror plane m
275        do isym=1,nsymconv
276          if(t_axes(isym)==15)then
277 !          Find direction of mirror plane
278            if(symrelconv(1,1,isym)==-1)direction=1
279            if(symrelconv(2,2,isym)==-1)direction=2
280            if(symrelconv(3,3,isym)==-1)direction=3
281          end if
282        end do
283 !      Examine the projection of the translation vector of the a, b, or c mirror planes
284 !      onto the binary axis
285        do isym=1,nsymconv
286          if(t_axes(isym)==16)then
287            if(abs(tnonsconv(direction,isym))>tol8)additional_info=2
288          end if
289        end do
290      end if
291      n_axest=(/0,0,0,0,1,0,0,1,1,0,  0,0,0,0,0,2,0,1,0,2,  0,0,0,0,0,0,0,0,0,0,0/)
292      if(sum((n_axes-n_axest)**2)==0)then    ! Spgroup 56 (Pccn) or 60 (Pbcn)
293        write(std_out,*)' symspgr: 56 or 60'
294        additional_info=1
295 !      Select mirror plane n
296        do isym=1,nsymconv
297          if(t_axes(isym)==18)then
298 !          Find direction of mirror plane
299            if(symrelconv(1,1,isym)==-1)direction=1
300            if(symrelconv(2,2,isym)==-1)direction=2
301            if(symrelconv(3,3,isym)==-1)direction=3
302          end if
303        end do
304 !      Examine the projection of the translation vector of the a, b, or c mirror planes
305 !      onto the binary axis
306        do isym=1,nsymconv
307          if(t_axes(isym)==16)then
308            if(abs(tnonsconv(direction,isym))<tol8)additional_info=2
309          end if
310        end do
311      end if
312    end if
313  else if(brvltt==2)then
314 !  In the few next lines, use additional_info as a flag
315    additional_info=0
316 !  If the bravais lattice is inner-centered
317    if(nsymconv==8)then
318 !    Test spgroup 23 (I222) or 24 (I2_{1}2_{1}2_{1})
319      n_axest=(/0,0,0,0,0,0,1,1,3,0,  0,0,0,0,0,0,0,0,0,3,  0,0,0,0,0,0,0,0,0,0,0/)
320      if(sum((n_axes-n_axest)**2)==0) additional_info=1
321    else if(nsymconv==24)then
322 !    Test spgroup 197 (I23) or 199 (I2_{1}3)
323      n_axest=(/0,0,0,0,0,0,1,1,3,16, 0,0,0,0,0,0,0,0,0,3,  0,0,0,0,0,0,0,0,0,0,0/)
324      if(sum((n_axes-n_axest)**2)==0) additional_info=1
325    end if
326    if(additional_info==1)then
327      write(std_out,*)' symspgr: (23 or 24) or (197 or 199)'
328 !    Select the three binary axes (they might be 2 or 2_1 !)
329      test_direction(:)=0
330      do isym=1,nsymconv
331        if(t_axes(isym)==20)then
332 !        Find direction of axis
333          do direction=1,3
334            if(symrelconv(direction,direction,isym)==1)then
335              test_direction(direction)=1
336              if(abs(tnonsconv(direction,isym))<tol8)then
337                vect(:,direction)=tnonsconv(:,isym)
338              else
339                vect(:,direction)=tnonsconv(:,isym)+half
340              end if
341              vect(:,direction)=vect(:,direction)-nint(vect(:,direction)-tol8)
342              vect(direction,direction)=zero
343            end if
344          end do ! direction=1,3
345        end if ! if binary axis
346      end do ! isym
347      if(test_direction(1)/=1 .or. test_direction(2)/=1 .and. test_direction(3)/=1)then
348        write(message, '(5a,3i4)' )&
349 &       'For space groups 23, 24, 197 or 197, the three binary axes',ch10,&
350 &       'are not equally partitioned along the x, y and z directions',ch10,&
351 &       'test_direction(1:3)=',test_direction(:)
352        MSG_BUG(message)
353      end if
354      additional_info=1
355      if(abs(vect(1,2)-vect(1,3))>tol8 .or. &
356 &     abs(vect(2,1)-vect(2,3))>tol8 .or. &
357 &     abs(vect(3,1)-vect(3,2))>tol8) additional_info=2
358    end if ! additional informations are needed
359  end if ! brvltt==1
360 
361  if (brvltt==0 .or. brvltt==1) then ! Primitive
362    call symlist_prim(additional_info,nsymconv,n_axes,spgroup)
363  else if(brvltt==2)then
364    call symlist_bcc(additional_info,nsymconv,n_axes,spgroup)
365  else if(brvltt==3)then
366    call symlist_fcc(nsymconv,n_axes,spgroup)
367  else
368    call symlist_others(brvltt,nsymconv,n_axes,spgroup)
369  end if
370 
371  if(spgroup==0) then
372    write(message, '(a,a,a,a,a)' )&
373 &   'Could not find the space group.',ch10,&
374 &   'This often happens when the user selects a restricted set of symmetries ',ch10,&
375 &   'in the input file, instead of letting the code automatically find symmetries.'
376    MSG_WARNING(message)
377  end if
378 
379  spgorig=1 ; spgaxor=1
380  call spgdata(brvsb,intsb,intsbl,ptintsb,ptschsb,schsb,spgaxor,spgroup,sporder,spgorig)
381 
382  if(spgroup/=0)then
383    write(message, '(a,i4,2x,a,a,a,a,a)' ) ' symspgr : spgroup=',spgroup,trim(brvsb),trim(intsb),'   (=',trim(schsb),')'
384    call wrtout(std_out,message,'COLL')
385  end if
386 
387  if(bravais(1)==7)then
388    write(message, '(a)' ) ' symspgr : optical characteristics = isotropic '
389    call wrtout(std_out,message,'COLL')
390  else if(bravais(1)==4 .or. bravais(1)==5 .or. bravais(1)==6)then
391    write(message, '(a)' ) ' symspgr : optical characteristics = uniaxial '
392    call wrtout(std_out,message,'COLL')
393 !  Identify the first symmetry operation that is order 3, 4 or 6
394    found=0
395    do isym=1,nsym
396 !    Proper rotations
397      if( minval( abs( t_axes(isym)-(/10,12,14,22,23,24,25,26,27,28,29,30,31/) ))==0) then
398        found=1 ; exit
399 !   Improper symmetry operations
400      else if( minval( abs( t_axes(isym)-(/1,2,3/) ))==0) then
401        found=-1 ; exit
402      end if
403    end do
404    if(found==-1 .or. found==1)then
405      symrel_uni=symrel(:,:,isym)
406      if(found==-1)symrel_uni=-symrel_uni
407 !    Now, symrel_uni is a rotation of order 3, 4, 6, for which the axis must be identified
408 !    It is actually the only eigenvector with eigenvalue 1. It can be found by cross products
409 !    Subtract the unit matrix.
410      do ii=1,3
411        symrel_uni(ii,ii)=symrel_uni(ii,ii)-1
412      end do
413      found=0
414      do ii=1,3
415        jj=ii+1 ; if(jj==4)jj=1
416 !      Cross product
417        uniaxis=symrel_uni(ii,:).x.symrel_uni(jj,:)
418        if(sum(uniaxis**2)/=0)then
419          found=1 ; exit
420        end if
421      end do
422      if(found==1)then
423 !      Try to reduce the length, by an integer factor (try only primes 2, 3, 5, 7, 11)
424        prime=(/2,3,5,7,11/)
425        ii=1
426        do while (ii<6)
427          uniaxis_try=uniaxis/prime(ii)
428          if(sum(abs(uniaxis_try*prime(ii)-uniaxis))==0)then
429            uniaxis=uniaxis_try
430          else
431            ii=ii+1
432          end if
433        end do
434        write(message, '(a,3i4)' ) ' Optical axis (in reduced coordinates, real space ) :',uniaxis
435      end if
436    end if
437    if(found==0)then
438      write(message, '(a)' ) ' However, the axis has not been found. Sorry for this.'
439    end if
440    call wrtout(std_out,message,'COLL')
441  end if
442 
443  ABI_DEALLOCATE(determinant)
444  ABI_DEALLOCATE(shift)
445  ABI_DEALLOCATE(symrelconv)
446  ABI_DEALLOCATE(tnonsconv)
447  ABI_DEALLOCATE(t_axes)
448 
449  DBG_EXIT("COLL")
450 
451 end subroutine symspgr