1!{\src2tex{textfont=tt}}
2!!****f* ABINIT/symlatt
3!! NAME
4!! symlatt
5!!
6!! FUNCTION
7!! From the unit cell vectors (rprimd) and the corresponding metric tensor,
8!! find the Bravais lattice and its symmetry operations (ptsymrel).
9!! 1) Find the shortest possible primitive vectors for the lattice
10!! 2) Determines the holohedral group of the lattice, and the
11!!    axes to be used for the conventional cell
12!!    (this is a delicate part, in which the centering of the
13!!    reduced cell must be taken into account)
14!!    The idea is to determine the basis vectors of the conventional
15!!    cell from the reduced cell basis vectors.
16!! 3) Generate the symmetry operations of the holohedral group
17!!
18!! COPYRIGHT
19!! Copyright (C) 2000-2016 ABINIT group (XG)
20!! This file is distributed under the terms of the
21!! GNU General Public License, see ~abinit/COPYING
22!! or http://www.gnu.org/copyleft/gpl.txt .
23!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
24!!
25!! INPUTS
26!! msym=default maximal number of symmetries
27!! rprimd(3,3)=dimensional primitive translations for real space (bohr)
28!! tolsym=tolerance for the symmetries
29!!
30!! OUTPUT
31!!  bravais(11): bravais(1)=iholohedry
32!!               bravais(2)=center
33!!               bravais(3:11)=coordinates of rprim in the axes
34!!               of the conventional bravais lattice (*2 if center/=0)
35!! nptsym=number of point symmetries of the Bravais lattice
36!! ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
37!! of the Bravais lattice in real space in terms
38!! of primitive translations.
39!!
40!! NOTES
41!! WARNING: bravais(1) might be given a negative value in another
42!! routine, if the cell is non-primitive.
43!! The holohedral groups are numbered as follows
44!! (see international tables for crystallography (1983), p. 13)
45!! iholohedry=1   triclinic      1bar
46!! iholohedry=2   monoclinic     2/m
47!! iholohedry=3   orthorhombic   mmm
48!! iholohedry=4   tetragonal     4/mmm
49!! iholohedry=5   trigonal       3bar m
50!! iholohedry=6   hexagonal      6/mmm
51!! iholohedry=7   cubic          m3bar m
52!! Centering
53!! center=0        no centering
54!! center=-1       body-centered
55!! center=-3       face-centered
56!! center=1        A-face centered
57!! center=2        B-face centered
58!! center=3        C-face centered
59!!
60!! PARENTS
61!!      ingeo,inqpt,m_ab7_symmetry,m_use_ga,symanal,symbrav,thmeig
62!!
63!! CHILDREN
64!!      holocell,matr3inv,smallprim,symrelrot,wrtout
65!!
66!! SOURCE
67
68#if defined HAVE_CONFIG_H
69#include "config.h"
70#endif
71
72#include "abi_common.h"
73
74
75subroutine symlatt(bravais,msym,nptsym,ptsymrel,rprimd,tolsym)
76
77 use defs_basis
78 use m_errors
79 use m_profiling_abi
80
81!This section has been created automatically by the script Abilint (TD).
82!Do not modify the following lines by hand.
83#undef ABI_FUNC
84#define ABI_FUNC 'symlatt'
85 use interfaces_14_hidewrite
86 use interfaces_32_util
87 use interfaces_41_geometry, except_this_one => symlatt
88!End of the abilint section
89
90 implicit none
91
92!Arguments ------------------------------------
93!scalars
94 integer,intent(in) :: msym
95 integer,intent(out) :: nptsym
96 real(dp),intent(in) :: tolsym
97!arrays
98 integer,intent(out) :: bravais(11),ptsymrel(3,3,msym)
99 real(dp),intent(in) :: rprimd(3,3)
100
101!Local variables-------------------------------
102!scalars
103 integer,parameter :: mgen=4
104 integer :: center,fact,found,foundc,ia,ib,icase,igen,iholohedry,ii,index,isym
105 integer :: itrial,jj,jsym,ngen=0,orthogonal,sign12,sign13,sign23,sumsign
106 real(dp) :: determinant,norm2a,norm2b,norm2c,norm2trial,reduceda,reducedb,sca
107 real(dp) :: scalarprod,scb,trace,val
108 character(len=500) :: message
109!arrays
110 integer,parameter :: list_holo(7)=(/7,6,4,3,5,2,1/)
111 integer :: ang90(3),equal(3),gen(3,3,mgen),gen2xy(3,3),gen2y(3,3),gen2z(3,3)
112 integer :: gen3(3,3),gen6(3,3),icoord(3,3),identity(3,3),nvecta(3),nvectb(3)
113 integer :: order(mgen)
114 real(dp) :: axes(3,3),axesinvt(3,3),cell_base(3,3),coord(3,3),metmin(3,3)
115 real(dp) :: minim(3,3),scprods(3,3),vecta(3),vectb(3),vectc(3),vin1(3),vin2(3),vext(3)
116
117!**************************************************************************
118
119 identity(:,:)=0 ; identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
120 nvecta(1)=2 ; nvectb(1)=3
121 nvecta(2)=1 ; nvectb(2)=3
122 nvecta(3)=1 ; nvectb(3)=2
123
124!--------------------------------------------------------------------------
125!Reduce the input vectors to a set of minimal vectors
126 call smallprim(metmin,minim,rprimd)
127
128!DEBUG
129!write(std_out,*)' symlatt : minim(:,1)=',minim(:,1)
130!write(std_out,*)' symlatt : minim(:,2)=',minim(:,2)
131!write(std_out,*)' symlatt : minim(:,3)=',minim(:,3)
132!ENDDEBUG
133
134!--------------------------------------------------------------------------
135!Examine the angles and vector lengths
136 ang90(:)=0
137 if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1
138 if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1
139 if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1
140 equal(:)=0
141 if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1
142 if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1
143 if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1
144
145!DEBUG
146!write(std_out,*)' ang90=',ang90(:)
147!write(std_out,*)' equal=',equal(:)
148!ENDDEBUG
149
150!-----------------------------------------------------------------------
151!Identification of the centering
152
153 foundc=0
154!Default values
155 fact=1 ; center=0
156 cell_base(:,:)=minim(:,:)
157
158!Examine each holohedral group
159!This search is ordered : should not be happy with tetragonal,
160!while there is FCC ...
161 do index=1,6
162
163!  If the holohedry is already found, exit
164   if(foundc==1)exit
165
166!  Initialize the target holohedry
167   iholohedry=list_holo(index)
168
169!  DEBUG
170!  write(std_out,*)' symlatt : trial holohedry',iholohedry
171!  ENDDEBUG
172
173   orthogonal=0
174   if(iholohedry==7 .or. iholohedry==4 .or. iholohedry==3)orthogonal=1
175
176!  Now, will examine different working hypothesis.
177!  The set of these hypothesis is thought to cover all possible cases ...
178
179!  Working hypothesis : the basis is orthogonal
180   if(ang90(1)+ang90(2)+ang90(3)==3 .and. orthogonal==1)then
181     fact=1 ; center=0
182     cell_base(:,:)=minim(:,:)
183!    Checks that the basis vectors are OK for the target holohedry
184     call holocell(cell_base,0,foundc,iholohedry,tolsym)
185   end if
186
187!  Select one trial direction
188   do itrial=1,3
189
190!    If the holohedry is already found, exit
191     if(foundc==1)exit
192
193     ia=nvecta(itrial) ; ib=nvectb(itrial)
194
195!    This is in case of hexagonal holohedry
196     if(foundc==0 .and. iholohedry==6 .and. ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then
197       reduceda=metmin(ib,ia)/metmin(ia,ia)
198       fact=1 ; center=0
199       if(abs(reduceda+0.5d0)<tolsym)then
200         cell_base(:,1)=minim(:,ia)
201         cell_base(:,2)=minim(:,ib)
202         cell_base(:,3)=minim(:,itrial)
203!        Checks that the basis vectors are OK for the target holohedry
204         call holocell(cell_base,0,foundc,iholohedry,tolsym)
205       else if(abs(reduceda-0.5d0)<tolsym)then
206         cell_base(:,1)=minim(:,ia)
207         cell_base(:,2)=minim(:,ib)-minim(:,ia)
208         cell_base(:,3)=minim(:,itrial)
209!        Checks that the basis vectors are OK for the target holohedry
210         call holocell(cell_base,0,foundc,iholohedry,tolsym)
211       end if
212     end if
213
214!    Working hypothesis : the conventional cell is orthogonal,
215!    and the two other vectors are axes of the conventional cell
216     if(foundc==0 .and. orthogonal==1 .and. ang90(itrial)==1)then
217
218!      Compute the reduced coordinate of trial vector in the basis
219!      of the two other vectors
220       reduceda=metmin(itrial,ia)/metmin(ia,ia)
221       reducedb=metmin(itrial,ib)/metmin(ib,ib)
222       cell_base(:,ia)=minim(:,ia)
223       cell_base(:,ib)=minim(:,ib)
224       if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. &
225&       ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym)       )then
226         if(abs(abs(reduceda)-0.5d0)<tolsym)center=ib
227         if(abs(abs(reducedb)-0.5d0)<tolsym)center=ia
228         fact=2
229         cell_base(:,itrial)= &
230&         (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0
231         call holocell(cell_base,0,foundc,iholohedry,tolsym)
232       else if( abs(abs(reduceda)-0.5d0)<tolsym .and.&
233&         abs(abs(reducedb)-0.5d0)<tolsym       ) then
234         fact=2 ; center=-1
235         cell_base(:,itrial)= &
236&         (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0
237         call holocell(cell_base,0,foundc,iholohedry,tolsym)
238       end if
239     end if
240
241!    Working hypothesis : the conventional cell is orthogonal, and
242!    the trial vector is one of the future axes, and the face perpendicular to it is centered
243     if(foundc==0 .and. iholohedry==3 .and. &
244&     ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then
245       fact=2 ; center=itrial
246       cell_base(:,ia)=minim(:,ia)+minim(:,ib)
247       cell_base(:,ib)=minim(:,ia)-minim(:,ib)
248       cell_base(:,itrial)=minim(:,itrial)
249!      Checks that the basis vectors are OK for the target holohedry
250       call holocell(cell_base,0,foundc,iholohedry,tolsym)
251     end if
252
253!    DEBUG
254!    write(std_out,*)' after test_b, foundc=',foundc
255!    ENDDEBUG
256
257!    Working hypothesis : the conventional cell is orthogonal, and
258!    the trial vector is one of the future axes
259     if(foundc==0 .and. orthogonal==1)then
260!      Compute the projection of the two other vectors on the trial vector
261       reduceda=metmin(itrial,ia)/metmin(itrial,itrial)
262       reducedb=metmin(itrial,ib)/metmin(itrial,itrial)
263!      If both projections are half-integer, one might have found an axis
264       if( abs(abs(reduceda)-0.5d0)<tolsym .and.&
265&       abs(abs(reducedb)-0.5d0)<tolsym       ) then
266         vecta(:)=minim(:,ia)-reduceda*minim(:,itrial)
267         vectb(:)=minim(:,ib)-reducedb*minim(:,itrial)
268         norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
269         norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
270         scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
271!        Note the order of selection : body-centered is prefered
272!        over face centered, which is correct for the tetragonal case
273         if(abs(norm2a-norm2b)<tolsym*half*(norm2a+norm2b))then
274!          The lattice is body centered
275           fact=2 ; center=-1
276           cell_base(:,ia)=vecta(:)+vectb(:)
277           cell_base(:,ib)=vecta(:)-vectb(:)
278           cell_base(:,itrial)=minim(:,itrial)
279           call holocell(cell_base,0,foundc,iholohedry,tolsym)
280         else if(abs(scalarprod)<tolsym*half*(norm2a+norm2b))then
281!          The lattice is face centered
282           fact=2 ; center=-3
283           cell_base(:,ia)=2.0d0*vecta(:)
284           cell_base(:,ib)=2.0d0*vectb(:)
285           cell_base(:,itrial)=minim(:,itrial)
286           call holocell(cell_base,0,foundc,iholohedry,tolsym)
287         end if
288       end if
289     end if
290
291!    DEBUG
292!    write(std_out,*)' after test_c, foundc=',foundc
293!    ENDDEBUG
294
295!    Working hypothesis : the conventional cell is orthogonal,
296!    and body centered with no basis vector being an axis,
297!    in which case the basis vectors must be equal (even for orthorhombic)
298     if(foundc==0 .and. orthogonal==1 .and. &
299&     equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then
300!      Compute the combination of the two other vectors
301       vecta(:)=minim(:,ia)+minim(:,ib)
302       vectb(:)=minim(:,ia)-minim(:,ib)
303       norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
304       norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
305!      Project the trial vector on the first of the two vectors
306       reduceda=( minim(1,itrial)*vecta(1)+       &
307&       minim(2,itrial)*vecta(2)+       &
308&       minim(3,itrial)*vecta(3) )/norm2a
309       reducedb=( minim(1,itrial)*vectb(1)+       &
310&       minim(2,itrial)*vectb(2)+       &
311&       minim(3,itrial)*vectb(3) )/norm2b
312       if( abs(abs(reduceda)-0.5d0)<tolsym )then
313!        The first vector is an axis
314         fact=2 ; center=-1
315         cell_base(:,ia)=vecta(:)
316         vecta(:)=minim(:,itrial)-reduceda*vecta(:)
317         vectb(:)=0.5d0*vectb(:)
318         cell_base(:,ib)=vecta(:)+vectb(:)
319         cell_base(:,itrial)=vecta(:)-vectb(:)
320         call holocell(cell_base,0,foundc,iholohedry,tolsym)
321       else if( abs(abs(reducedb)-0.5d0)<tolsym )then
322!        The second vector is an axis
323         fact=2 ; center=-1
324         cell_base(:,ib)=vectb(:)
325         vectb(:)=minim(:,itrial)-reducedb*vectb(:)
326         vecta(:)=0.5d0*vecta(:)
327         cell_base(:,ia)=vectb(:)+vecta(:)
328         cell_base(:,itrial)=vectb(:)-vecta(:)
329         call holocell(cell_base,0,foundc,iholohedry,tolsym)
330       end if
331     end if
332
333!    Working hypothesis : the conventional cell is orthogonal,
334!    and face centered, in the case where two minimal vectors are equal
335     if(foundc==0 .and. orthogonal==1 .and. equal(itrial)==1 ) then
336!      Compute the combination of these two vectors
337       vecta(:)=minim(:,ia)+minim(:,ib)
338       vectb(:)=minim(:,ia)-minim(:,ib)
339       norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
340       norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
341!      Project the trial vector on the two vectors
342       reduceda=( minim(1,itrial)*vecta(1)+       &
343&       minim(2,itrial)*vecta(2)+       &
344&       minim(3,itrial)*vecta(3) )/norm2a
345       reducedb=( minim(1,itrial)*vectb(1)+       &
346&       minim(2,itrial)*vectb(2)+       &
347&       minim(3,itrial)*vectb(3) )/norm2b
348       if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. &
349&       ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym)       )then
350         fact=2 ; center=-3
351         cell_base(:,itrial)= &
352&         (minim(:,itrial)-reduceda*vecta(:)-reducedb*vectb(:) )*2.0d0
353         cell_base(:,ia)=vecta(:)
354         cell_base(:,ib)=vectb(:)
355         call holocell(cell_base,0,foundc,iholohedry,tolsym)
356       end if
357     end if
358
359!    Working hypothesis : the conventional cell is orthogonal,
360!    face centered, but no two vectors are on the same "square"
361     if(foundc==0 .and. orthogonal==1)then
362!      Compute the combination of these two vectors
363       vecta(:)=minim(:,ia)+minim(:,ib)
364       vectb(:)=minim(:,ia)-minim(:,ib)
365       norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
366       norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
367!      The trial vector length must be equal to one of these lengths
368       if(abs(metmin(itrial,itrial)-norm2a)<tolsym*norm2a)then
369         fact=2 ; center=-3
370         cell_base(:,ia)=vecta(:)+minim(:,itrial)
371         cell_base(:,ib)=vecta(:)-minim(:,itrial)
372!        Project vectb perpendicular to cell_base(:,ia) and cell_base(:,ib)
373         norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2
374         norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2
375         reduceda=( cell_base(1,ia)*vectb(1)+       &
376&         cell_base(2,ia)*vectb(2)+       &
377&         cell_base(3,ia)*vectb(3) )/norm2a
378         reducedb=( cell_base(1,ib)*vectb(1)+       &
379&         cell_base(2,ib)*vectb(2)+       &
380&         cell_base(3,ib)*vectb(3) )/norm2b
381         if( abs(abs(reduceda)-0.5d0)<tolsym .and.         &
382&         abs(abs(reducedb)-0.5d0)<tolsym      )then
383           cell_base(:,itrial)=vectb(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib)
384           call holocell(cell_base,0,foundc,iholohedry,tolsym)
385         end if
386       else if(abs(metmin(itrial,itrial)-norm2b)<tolsym*norm2b)then
387         fact=2 ; center=-3
388         cell_base(:,ia)=vectb(:)+minim(:,itrial)
389         cell_base(:,ib)=vectb(:)-minim(:,itrial)
390!        Project vecta perpendicular to cell_base(:,ia) and cell_base(:,ib)
391         norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2
392         norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2
393         reduceda=( cell_base(1,ia)*vecta(1)+       &
394&         cell_base(2,ia)*vecta(2)+       &
395&         cell_base(3,ia)*vecta(3) )/norm2a
396         reducedb=( cell_base(1,ib)*vecta(1)+       &
397&         cell_base(2,ib)*vecta(2)+       &
398&         cell_base(3,ib)*vecta(3) )/norm2b
399         if( abs(abs(reduceda)-0.5d0)<tolsym .and.         &
400&         abs(abs(reducedb)-0.5d0)<tolsym      )then
401           cell_base(:,itrial)=vecta(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib)
402           call holocell(cell_base,0,foundc,iholohedry,tolsym)
403         end if
404       end if
405     end if
406
407!    Working hypothesis : the cell is rhombohedral, and
408!    the three minimal vectors have same length and same absolute scalar product
409     if(foundc==0 .and. iholohedry==5 .and. &
410&     equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then
411       if( abs(abs(metmin(1,2))-abs(metmin(1,3)))<tolsym*metmin(1,1) .and.     &
412&       abs(abs(metmin(1,2))-abs(metmin(2,3)))<tolsym*metmin(2,2)      )then
413         fact=1 ; center=0
414         cell_base(:,:)=minim(:,:)
415!        One might have to change the sign of one of the vectors
416         sign12=1 ; sign13=1 ; sign23=1
417         if(metmin(1,2)<0.0d0)sign12=-1
418         if(metmin(1,3)<0.0d0)sign13=-1
419         if(metmin(2,3)<0.0d0)sign23=-1
420         sumsign=sign12+sign13+sign23
421         if(sumsign==-1)then
422           if(sign12==1)cell_base(:,3)=-cell_base(:,3)
423           if(sign13==1)cell_base(:,2)=-cell_base(:,2)
424           if(sign23==1)cell_base(:,1)=-cell_base(:,1)
425         else if(sumsign==1)then
426           if(sign12==-1)cell_base(:,3)=-cell_base(:,3)
427           if(sign13==-1)cell_base(:,2)=-cell_base(:,2)
428           if(sign23==-1)cell_base(:,1)=-cell_base(:,1)
429         end if
430         call holocell(cell_base,0,foundc,iholohedry,tolsym)
431       end if
432     end if
433
434!    DEBUG
435!    write(std_out,*)' after test_3a, foundc=',foundc
436!    write(std_out,*)' after test_3a, itrial=',itrial
437!    write(std_out,*)' after test_3a, equal(:)=',equal(:)
438!    ENDDEBUG
439
440!    Working hypothesis : the cell is rhombohedral, one vector
441!    is parallel to the trigonal axis
442     if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 )then
443       vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib)
444       norm2trial=minim(1,itrial)**2+minim(2,itrial)**2+minim(3,itrial)**2
445       reduceda=( minim(1,itrial)*vecta(1)+       &
446&       minim(2,itrial)*vecta(2)+       &
447&       minim(3,itrial)*vecta(3) )/norm2trial
448       reducedb=( minim(1,itrial)*vectb(1)+       &
449&       minim(2,itrial)*vectb(2)+       &
450&       minim(3,itrial)*vectb(3) )/norm2trial
451!      DEBUG
452!      write(std_out,*)' reduceda,reducedb=',reduceda,reducedb
453!      ENDDEBUG
454       if(abs(abs(reduceda)-1.0d0/3.0d0)<tolsym .and.      &
455&       abs(abs(reducedb)-1.0d0/3.0d0)<tolsym      ) then
456!        Possibly change of sign to make positive the scalar product with
457!        the vector parallel to the trigonal axis
458         if(reduceda<zero)vecta(:)=-vecta(:)
459         if(reducedb<zero)vectb(:)=-vectb(:)
460!        Projection on the orthogonal plane
461         vecta(:)=vecta(:)-abs(reduceda)*cell_base(:,itrial)
462         vectb(:)=vectb(:)-abs(reducedb)*cell_base(:,itrial)
463!        These two vectors should have an angle of 120 degrees
464         norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
465         scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
466!        DEBUG
467!        write(std_out,*)' norm2a,scalarprod=',norm2a,scalarprod
468!        ENDDEBUG
469         if(abs(two*scalarprod+norm2a)<tolsym*norm2a)then
470           fact=1 ; center=0
471           if(scalarprod>0.0d0)vectb(:)=-vectb(:)
472!          Now vecta and vectb have an angle of 120 degrees
473           cell_base(:,1)=cell_base(:,itrial)/3.0d0+vecta(:)
474           cell_base(:,2)=cell_base(:,itrial)/3.0d0+vectb(:)
475           cell_base(:,3)=cell_base(:,itrial)/3.0d0-vecta(:)-vectb(:)
476!          DEBUG
477!          write(std_out,*)' cell_base(:,1)=',cell_base(:,1)
478!          write(std_out,*)' cell_base(:,2)=',cell_base(:,2)
479!          write(std_out,*)' cell_base(:,3)=',cell_base(:,3)
480!          ENDDEBUG
481           call holocell(cell_base,0,foundc,iholohedry,tolsym)
482         end if
483       end if
484     end if
485
486!    Working hypothesis : the cell is rhombohedral, one vector
487!    is in the plane perpendicular to the trigonal axis
488     if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then
489       vecta(:)=minim(:,ia)+minim(:,ib)
490       vectb(:)=minim(:,ia)-minim(:,ib)
491       norm2trial=cell_base(1,itrial)**2+cell_base(2,itrial)**2+cell_base(3,itrial)**2
492       norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
493       norm2b=vecta(1)**2+vecta(2)**2+vecta(3)**2
494       reduceda=( cell_base(1,itrial)*vecta(1)+       &
495&       cell_base(2,itrial)*vecta(2)+       &
496&       cell_base(3,itrial)*vecta(3) )/norm2trial
497       reducedb=( cell_base(1,itrial)*vectb(1)+       &
498&       cell_base(2,itrial)*vectb(2)+       &
499&       cell_base(3,itrial)*vectb(3) )/norm2trial
500       if(abs(norm2trial-norm2a)<tolsym*norm2a .and. &
501&       abs(abs(2*reduceda)-norm2trial)<tolsym*norm2trial    )then
502         fact=1 ; center=0
503         cell_base(:,1)=minim(:,ia)
504         cell_base(:,2)=-minim(:,ib)
505         cell_base(:,3)=-minim(:,ib)+2*reduceda*minim(:,itrial)
506         call holocell(cell_base,0,foundc,iholohedry,tolsym)
507       else if (abs(norm2trial-norm2b)<tolsym*norm2b .and. &
508&         abs(abs(2*reducedb)-norm2trial)<tolsym*norm2trial    )then
509         fact=1 ; center=0
510         cell_base(:,1)=minim(:,ia)
511         cell_base(:,2)=minim(:,ib)
512         cell_base(:,3)=minim(:,ib)+2*reducedb*minim(:,itrial)
513         call holocell(cell_base,0,foundc,iholohedry,tolsym)
514       end if
515     end if
516
517!    Working hypothesis : the cell is rhombohedral, two vectors
518!    are in the plane perpendicular to the trigonal axis
519     if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then
520       vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib)
521       norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
522       norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
523       scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
524       if(abs(abs(2*scalarprod)-norm2a)<tolsym*norm2a)then
525!        This is in order to have 120 angle between vecta and vectb
526         if(scalarprod>0.0d0)vectb(:)=-vectb(:)
527         reduceda=( cell_base(1,itrial)*vecta(1)+        &
528&         cell_base(2,itrial)*vecta(2)+        &
529&         cell_base(3,itrial)*vecta(3) )/norm2a
530         reducedb=( cell_base(1,itrial)*vectb(1)+        &
531&         cell_base(2,itrial)*vectb(2)+        &
532&         cell_base(3,itrial)*vectb(3) )/norm2b
533         fact=1 ; center=0
534         cell_base(:,1)=minim(:,itrial)
535         if(abs(reduceda-0.5d0)<tolsym .and. abs(reducedb)<tolsym )then
536           cell_base(:,2)=minim(:,itrial)-vecta(:)
537           cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:)
538           call holocell(cell_base,0,foundc,iholohedry,tolsym)
539         else if(abs(reduceda-0.5d0)<tolsym.and. abs(reducedb+0.5d0)<tolsym )then
540           cell_base(:,2)=minim(:,itrial)-vecta(:)
541           cell_base(:,3)=minim(:,itrial)+vectb(:)
542           call holocell(cell_base,0,foundc,iholohedry,tolsym)
543         else if(abs(reduceda)<tolsym .and. abs(reducedb+0.5d0)<tolsym )then
544           cell_base(:,2)=minim(:,itrial)+vectb(:)
545           cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:)
546           call holocell(cell_base,0,foundc,iholohedry,tolsym)
547         else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb)<tolsym)then
548           cell_base(:,2)=minim(:,itrial)+vecta(:)
549           cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:)
550           call holocell(cell_base,0,foundc,iholohedry,tolsym)
551         else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb-0.5d0)<tolsym)then
552           cell_base(:,2)=minim(:,itrial)+vecta(:)
553           cell_base(:,3)=minim(:,itrial)-vectb(:)
554           call holocell(cell_base,0,foundc,iholohedry,tolsym)
555         else if(abs(reduceda)<tolsym .and. abs(reducedb-0.5d0)<tolsym )then
556           cell_base(:,2)=minim(:,itrial)-vectb(:)
557           cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:)
558           call holocell(cell_base,0,foundc,iholohedry,tolsym)
559         end if
560       end if
561     end if
562
563!    Working hypothesis : monoclinic holohedry, primitive. Then, two angles are 90 degrees
564     if(foundc==0 .and. iholohedry==2 .and. &
565&     ang90(ia)==1 .and. ang90(ib)==1 ) then
566       fact=1 ; center=0
567       cell_base(:,1)=minim(:,ia)
568       cell_base(:,2)=minim(:,itrial)
569       cell_base(:,3)=minim(:,ib)
570!      Checks that the basis vectors are OK for the target holohedry
571       call holocell(cell_base,0,foundc,iholohedry,tolsym)
572     end if
573
574!    Monoclinic holohedry, one-face-centered cell
575!    Working hypothesis, two vectors have equal length.
576     do icase=1,5
577       if(foundc==0 .and. iholohedry==2 .and. equal(itrial)==1 ) then
578         vecta(:)=cell_base(:,ia)+cell_base(:,ib)
579         vectb(:)=cell_base(:,ia)-cell_base(:,ib)
580         norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
581         norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
582!        The minim(:,trial) vector belongs to the
583!        plane parallel to the cell_base(:,ia),cell_base(:,ib) plane
584!        In that plane, must try minim(:,itrial),
585!        as well as the 4 different combinations of
586!        minim(:,itrial) with the vectors in the plane
587         if(icase==1)vectc(:)=minim(:,itrial)
588         if(icase==2)vectc(:)=minim(:,itrial)+cell_base(:,ia)
589         if(icase==3)vectc(:)=minim(:,itrial)+cell_base(:,ib)
590         if(icase==4)vectc(:)=minim(:,itrial)-cell_base(:,ia)
591         if(icase==5)vectc(:)=minim(:,itrial)-cell_base(:,ib)
592         norm2c=vectc(1)**2+vectc(2)**2+vectc(3)**2
593         sca=vectc(1)*vecta(1)+&
594&         vectc(2)*vecta(2)+&
595&         vectc(3)*vecta(3)
596         scb=vectc(1)*vectb(1)+&
597&         vectc(2)*vectb(2)+&
598&         vectc(3)*vectb(3)
599!        DEBUG
600!        write(std_out,*)' symlatt : test iholohedry=2, sca,scb=',sca,scb
601!        ENDDEBUG
602         if(abs(sca)<tolsym*sqrt(norm2c*norm2a) .or. abs(scb)<tolsym*sqrt(norm2c*norm2b))then
603           fact=2 ; center=3
604!          The itrial direction is centered
605           cell_base(:,3)=vectc(:)
606           if(abs(sca)<tolsym*sqrt(norm2c*norm2a))then
607             cell_base(:,2)=vecta(:)
608             cell_base(:,1)=vectb(:)
609             call holocell(cell_base,0,foundc,iholohedry,tolsym)
610           else if(abs(scb)<tolsym*sqrt(norm2c*norm2b))then
611             cell_base(:,2)=vectb(:)
612             cell_base(:,1)=vecta(:)
613             call holocell(cell_base,0,foundc,iholohedry,tolsym)
614           end if
615         end if
616       end if
617     end do ! icase=1,5
618
619!    Monoclinic holohedry, one-face-centered cell, but non equivalent.
620!    This case, one pair of vectors is orthogonal
621     if(foundc==0 .and. iholohedry==2 .and. ang90(itrial)==1) then
622       vecta(:)=minim(:,ia)
623       vectb(:)=minim(:,ib)
624       norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
625       norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
626!      Project the trial vector on the two vectors
627       reduceda=( minim(1,itrial)*vecta(1)+       &
628&       minim(2,itrial)*vecta(2)+       &
629&       minim(3,itrial)*vecta(3) )/norm2a
630       reducedb=( minim(1,itrial)*vectb(1)+       &
631&       minim(2,itrial)*vectb(2)+       &
632&       minim(3,itrial)*vectb(3) )/norm2b
633       if(abs(abs(reduceda)-0.5d0)<tolsym .or. abs(abs(reducedb)-0.5d0)<tolsym) then
634         fact=2 ; center=3
635         if(abs(abs(reduceda)-0.5d0)<tolsym)then
636           cell_base(:,2)=vecta(:)
637           cell_base(:,3)=vectb(:)
638           cell_base(:,1)=2*(minim(:,itrial)-reduceda*vecta(:))
639           call holocell(cell_base,0,foundc,iholohedry,tolsym)
640         else if(abs(abs(reducedb)-0.5d0)<tolsym)then
641           cell_base(:,2)=vectb(:)
642           cell_base(:,3)=vecta(:)
643           cell_base(:,1)=2*(minim(:,itrial)-reducedb*vectb(:))
644           call holocell(cell_base,0,foundc,iholohedry,tolsym)
645         end if
646       end if
647     end if
648
649!    Monoclinic holohedry, one-face-centered cell, but non equivalent.
650!    This case, no pair of vectors is orthogonal, no pair of vector of equal lentgh
651     if(foundc==0 .and. iholohedry==2)then
652!      Try to find a vector that belongs to the mediator plane, or the binary vector.
653!      There must be one such vector, if centered monoclinic and no pair of vectors of equal length,
654!      either among the three vectors, or among one of their differences or sums.
655!      And there must be, among the two other vectors, one vector whose projection
656!      on this vector is half the length of this vector.
657       vecta(:)=minim(:,ia)
658       vectb(:)=minim(:,ib)
659!      Try the different possibilities for the vector on which the projection will be half ...
660       do ii=1,5
661         if(ii==1)vectc(:)=minim(:,itrial)
662         if(ii==2)vectc(:)=minim(:,itrial)+vecta(:)
663         if(ii==3)vectc(:)=minim(:,itrial)-vecta(:)
664         if(ii==4)vectc(:)=minim(:,itrial)+vectb(:)
665         if(ii==5)vectc(:)=minim(:,itrial)-vectb(:)
666         norm2trial=vectc(1)**2+vectc(2)**2+vectc(3)**2
667!        Project the two vectors on the trial vector
668         reduceda=( vectc(1)*vecta(1)+vectc(2)*vecta(2)+vectc(3)*vecta(3) )/norm2trial
669         reducedb=( vectc(1)*vectb(1)+vectc(2)*vectb(2)+vectc(3)*vectb(3) )/norm2trial
670         found=0
671         if(abs(abs(reduceda)-0.5d0)<tolsym)then
672           vin1(:)=vectc(:)
673           vin2(:)=2.0d0*(vecta(:)-reduceda*vectc(:))
674           vext(:)=vectb(:)
675           found=1
676         else if(abs(abs(reducedb)-0.5d0)<tolsym)then
677           vin1(:)=vectc(:)
678           vin2(:)=2.0d0*(vectb(:)-reduceda*vectc(:))
679           vext(:)=vecta(:)
680           found=1
681         end if
682         if(found==1)exit
683       end do
684!      Now, vin1 and vin2 are perpendicular to each other, and in the plane that contains the binary vector.
685!      One of them must be the binary vector if any.
686!      On the other hand, vext is out-of-plane. Might belong to the mediator plane or not.
687!      If C monoclinc, then the projection of this vext on the binary vector will be either 0 or +1/2 or -1/2.
688!      The binary axis must be stored in cell_base(:,2) for conventional C-cell
689       if(found==1)then
690         found=0
691
692!        Test vin1 being the binary axis
693         norm2trial=vin1(1)**2+vin1(2)**2+vin1(3)**2
694         reduceda=(vext(1)*vin1(1)+vext(2)*vin1(2)+vext(3)*vin1(3))/norm2trial
695         if(abs(reduceda)<tolsym)then  ! vin1 is the binary axis and vext is in the mediator plane
696           found=1
697           cell_base(:,1)=vin2(:)
698           cell_base(:,2)=vin1(:)
699           cell_base(:,3)=vext(:)
700         else if(abs(abs(reduceda)-0.5d0)<tolsym)then  ! vin1 is the binary axis and vext has +1/2 or -1/2 as projection
701           found=1
702           cell_base(:,1)=vin2(:)
703           cell_base(:,2)=vin1(:)
704           cell_base(:,3)=vext(:)-reduceda*vin1(:)+vin2(:)*half
705         else
706!          Test vin2 being the binary axis
707           norm2trial=vin2(1)**2+vin2(2)**2+vin2(3)**2
708           reduceda=(vext(1)*vin2(1)+vext(2)*vin2(2)+vext(3)*vin2(3))/norm2trial
709           if(abs(reduceda)<tolsym)then  ! vin2 is the binary axis and vext is in the mediator plane
710             found=1
711             cell_base(:,1)=vin1(:)
712             cell_base(:,2)=vin2(:)
713             cell_base(:,3)=vext(:)
714           else if(abs(abs(reduceda)-0.5d0)<tolsym)then  ! vin2 is the binary axis and vext has +1/2 or -1/2 as projection
715             found=1
716             cell_base(:,1)=vin1(:)
717             cell_base(:,2)=vin2(:)
718             cell_base(:,3)=vext(:)-reduceda*vin2(:)+vin1(:)*half
719           end if
720         end if
721
722         if(found==1)then
723           fact=2 ; center=3
724           call holocell(cell_base,0,foundc,iholohedry,tolsym)
725         end if
726       end if
727     end if
728
729   end do ! Do-loop on three different directions
730 end do !  Do-loop on different target holohedries
731
732 if(foundc==0)then
733   iholohedry=1 ; fact=1 ; center=0
734   cell_base(:,:)=minim(:,:)
735 end if
736
737!DEBUG
738!write(std_out,*)' symlatt : done with centering tests, foundc=',foundc
739!write(std_out,*)'  center=',center
740!write(std_out,*)'  iholohedry=',iholohedry
741!ENDDEBUG
742
743!--------------------------------------------------------------------------
744!Final check on the Bravais lattice, using the basis vectors
745
746!Recompute the metric tensor
747 if(foundc==1)then
748   do ii=1,3
749     metmin(:,ii)=cell_base(1,:)*cell_base(1,ii)+&
750&     cell_base(2,:)*cell_base(2,ii)+&
751&     cell_base(3,:)*cell_base(3,ii)
752   end do
753 end if
754
755!Examine the angles and vector lengths
756 ang90(:)=0
757 if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1
758 if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1
759 if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1
760 equal(:)=0
761 if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1
762 if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1
763 if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1
764
765!DEBUG
766!write(std_out,*)' symlatt : recompute the  metric tensor '
767!write(std_out,*)'  ang90=',ang90
768!write(std_out,*)'  equal=',equal
769!ENDDEBUG
770
771!The axes will be aligned with the previously determined
772!basis vectors, EXCEPT for the tetragonal cell, see later
773 axes(:,:)=cell_base(:,:)
774
775 found=0
776!Check orthogonal conventional cells
777 if(ang90(1)+ang90(2)+ang90(3)==3)then
778
779!  Cubic system
780   if(equal(1)+equal(2)+equal(3)==3)then
781!    However, one-face centered is not admitted
782     if(center==0 .or. center==-1 .or. center==-3)then
783       iholohedry=7 ; found=1
784       if(center==0)then
785         write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is cP (primitive cubic)'
786       else if(center==-1)then
787         write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is cI (body-centered cubic)'
788       else if(center==-3)then
789         write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is cF (face-centered cubic)'
790       end if
791     end if
792   end if
793
794!  Tetragonal system
795   if(found==0 .and. &
796&   (equal(1)==1 .or. equal(2)==1 .or. equal(3)==1) )then
797!    However, one-face centered or face-centered is not admitted
798     if(center==0 .or. center==-1)then
799       iholohedry=4 ; found=1
800       if(equal(1)==1)then
801         axes(:,3)=cell_base(:,1) ; axes(:,1)=cell_base(:,2) ; axes(:,2)=cell_base(:,3)
802       else if(equal(2)==1)then
803         axes(:,3)=cell_base(:,2) ; axes(:,2)=cell_base(:,1) ; axes(:,1)=cell_base(:,3)
804       else if(equal(3)==1)then
805         axes(:,:)=cell_base(:,:)
806       end if
807       if(center==0)then
808         write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is tP (primitive tetragonal)'
809       else if(center==-1)then
810         write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is tI (body-centered tetragonal)'
811       end if
812     end if
813   end if
814
815!  Orthorhombic system
816   if(found==0)then
817     iholohedry=3 ; found=1
818     axes(:,:)=cell_base(:,:)
819     if(center==0)then
820       write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oP (primitive orthorhombic)'
821     else if(center==-1)then
822       write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oI (body-centered orthorhombic)'
823     else if(center==1 .or. center==2 .or. center==3)then
824       write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oC (one-face-centered orthorhombic)'
825     else if(center==-3)then
826       write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is oF (face-centered orthorhombic)'
827     end if
828   end if
829
830 else
831
832!  Hexagonal system
833   if(found==0 .and. ang90(1)==1 .and. ang90(2)==1 .and. equal(3)==1 .and. (2*metmin(2,1)+metmin(1,1))<tolsym*metmin(1,1))then
834     iholohedry=6 ; found=1
835     write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is hP (primitive hexagonal)'
836   end if
837
838!  Rhombohedral system
839   if(found==0 .and. equal(1)+equal(2)+equal(3)==3 .and.       &
840&   abs(metmin(2,1)-metmin(3,2))<tolsym*metmin(2,2)             .and.       &
841&   abs(metmin(2,1)-metmin(3,1))<tolsym*metmin(1,1) )then
842     iholohedry=5 ; found=1
843     write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is hR (rhombohedral)'
844   end if
845
846!  Monoclinic system
847   if(found==0 .and. ang90(1)+ang90(2)+ang90(3)==2 )then
848     iholohedry=2 ; found=1
849     if(center==0)then
850       write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is mP (primitive monoclinic)'
851     else if(center==3)then
852       write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is mC (one-face-centered monoclinic)'
853     end if
854   end if
855
856!  Triclinic system
857   if(found==0)then
858     iholohedry=1 ; found=1
859     write(message,'(a,a)')ch10,' symlatt : the Bravais lattice is aP (primitive triclinic)'
860   end if
861
862 end if
863
864 call wrtout(std_out,message,'COLL')
865
866!--------------------------------------------------------------------------
867!Make sure that axes form a right-handed coordinate system
868!(Note : this should be done in the body of the routine,
869!by making changes that leave the sign of the mixed product of the three
870!vectors invariant)
871 determinant=axes(1,1)*axes(2,2)*axes(3,3) &
872& +axes(1,2)*axes(2,3)*axes(3,1) &
873& +axes(1,3)*axes(3,2)*axes(2,1) &
874& -axes(1,1)*axes(3,2)*axes(2,3) &
875& -axes(1,3)*axes(2,2)*axes(3,1) &
876& -axes(1,2)*axes(2,1)*axes(3,3)
877 if(determinant<0.0d0)then
878   axes(:,:)=-axes(:,:)
879 end if
880
881!--------------------------------------------------------------------------
882!Prefer symmetry axes on the same side as the primitive axes,
883!when the changes are allowed
884 do
885   do ia=1,3
886     scprods(ia,:)=axes(1,ia)*rprimd(1,:)+&
887&     axes(2,ia)*rprimd(2,:)+&
888&     axes(3,ia)*rprimd(3,:)
889     norm2trial=sum(axes(:,ia)**2)
890     scprods(ia,:)=scprods(ia,:)/sqrt(norm2trial)
891   end do
892   do ia=1,3
893     norm2trial=sum(rprimd(:,ia)**2)
894     scprods(:,ia)=scprods(:,ia)/sqrt(norm2trial)
895   end do
896
897!  One should now try all the generators of the
898!  proper rotations of each Bravais lattice, coupled with change of
899!  signs of each vector. This is not done systematically in what follows ...
900!  Here, the third axis is left unchanged
901   if(iholohedry/=5)then
902     if(scprods(1,1)<-tolsym .and. scprods(2,2)<-tolsym)then
903       axes(:,1)=-axes(:,1) ; axes(:,2)=-axes(:,2)
904       cycle
905     end if
906   end if
907!  The first (or second) axis is left unchanged
908   if(iholohedry/=5 .and. iholohedry/=6)then
909     if(scprods(2,2)<-tolsym .and. scprods(3,3)<-tolsym)then
910       axes(:,2)=-axes(:,2) ; axes(:,3)=-axes(:,3)
911       cycle
912     end if
913     if(scprods(1,1)<-tolsym .and. scprods(3,3)<-tolsym)then
914       axes(:,1)=-axes(:,1) ; axes(:,3)=-axes(:,3)
915       cycle
916     end if
917   end if
918!  Permutation of the three axis
919   if(iholohedry==5 .or. iholohedry==7)then
920     trace=scprods(1,1)+scprods(2,2)+scprods(3,3)
921     if(trace+tolsym< scprods(1,2)+scprods(2,3)+scprods(3,1))then
922       vecta(:)=axes(:,1) ; axes(:,1)=axes(:,3)
923       axes(:,3)=axes(:,2); axes(:,2)=vecta(:)
924       cycle
925     end if
926     if(trace+tolsym < scprods(1,3)+scprods(2,1)+scprods(3,2))then
927       vecta(:)=axes(:,1) ; axes(:,1)=axes(:,2)
928       axes(:,2)=axes(:,3); axes(:,3)=vecta(:)
929       cycle
930     end if
931!    This case is observed when the three new vectors
932!    are pointing opposite to the three original vectors
933!    One takes their opposite, then switch to of them, then process
934!    them again in the loop
935     if(sum(scprods(:,:))<tolsym)then
936       axes(:,1)=-axes(:,1)
937       vecta(:)=-axes(:,2)
938       axes(:,2)=-axes(:,3)
939       axes(:,3)=vecta(:)
940       cycle
941     end if
942   end if
943   exit
944!  Other cases might be coded ...
945 end do
946
947!--------------------------------------------------------------------------
948
949!DEBUG
950!write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',&
951!&  rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
952!write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes  =',&
953!&  axes(:,1),ch10,axes(:,2),ch10,axes(:,3)
954!ENDDEBUG
955
956!Compute the coordinates of rprimd in the system defined by axes(:,:)
957 call matr3inv(axes,axesinvt)
958 do ii=1,3
959   coord(:,ii)=rprimd(1,ii)*axesinvt(1,:)+ &
960&   rprimd(2,ii)*axesinvt(2,:)+ &
961&   rprimd(3,ii)*axesinvt(3,:)
962 end do
963
964!Check that the coordinates are integers, or half-integer in
965!the case there is a centering, and generate integer coordinates
966 do ii=1,3
967   do jj=1,3
968     val=coord(ii,jj)*fact
969     if(abs(val-nint(val))>fact*tolsym)then
970       write(message,'(4a,a,3es18.10,a,a,3es18.10,a,a,3es18.10,a,a,i4)')&
971&       'One of the coordinates of rprimd in axes is non-integer,',ch10,&
972&       'or non-half-integer (if centering).',ch10,&
973&       'coord=',coord(:,1),ch10,&
974&       '      ',coord(:,2),ch10,&
975&       '      ',coord(:,3),ch10,&
976&       'fact=',fact
977       MSG_BUG(message)
978     end if
979     icoord(ii,jj)=nint(val)
980   end do
981 end do
982
983!Store the bravais lattice characteristics
984 bravais(1)=iholohedry
985 bravais(2)=center
986 bravais(3:5)=icoord(1:3,1)
987 bravais(6:8)=icoord(1:3,2)
988 bravais(9:11)=icoord(1:3,3)
989
990!--------------------------------------------------------------------------
991!Initialize the set of symmetries
992!Bravais lattices are always invariant under identity and inversion
993
994!Identity and inversion
995 ptsymrel(:,:,1)=identity(:,:) ; ptsymrel(:,:,2)=-identity(:,:)
996 nptsym=2
997
998!Keep this for IFCv70 compiler
999 if(nptsym/=2)then
1000   write(message,'(a,a,a,a)')ch10,&
1001&   ' symlatt : BUG -',ch10,&
1002&   '  Crazy error, compiler bug '
1003   call wrtout(std_out,message,'COLL')
1004 end if
1005
1006!--------------------------------------------------------------------------
1007!Initialize some generators
1008!gen6 is defined in a coordinated system with gamma=120 degrees
1009 gen6(:,:)=0  ; gen6(3,3)=1  ; gen6(1,1)=1  ; gen6(1,2)=-1 ; gen6(2,1)=1
1010 gen3(:,:)=0  ; gen3(1,2)=1  ; gen3(2,3)=1  ; gen3(3,1)=1
1011 gen2xy(:,:)=0 ; gen2xy(2,1)=1 ; gen2xy(1,2)=1; gen2xy(3,3)=1
1012 gen2y(:,:)=0 ; gen2y(1,1)=-1; gen2y(2,2)=1 ; gen2y(3,3)=-1
1013 gen2z(:,:)=0 ; gen2z(1,1)=-1; gen2z(2,2)=-1; gen2z(3,3)=1
1014
1015!--------------------------------------------------------------------------
1016
1017!Define the generators for each holohedry (inversion is already included)
1018 if(iholohedry==6)then
1019   ngen=2
1020   gen(:,:,1)=gen2xy(:,:) ; order(1)=2
1021   gen(:,:,2)=gen6(:,:)   ; order(2)=6
1022 else if(iholohedry==5)then
1023   ngen=2
1024   gen(:,:,1)=gen2xy(:,:) ; order(1)=2
1025   gen(:,:,2)=gen3(:,:)   ; order(2)=3
1026 else
1027   gen(:,:,1)=gen2y(:,:)  ; order(1)=2
1028   gen(:,:,2)=gen2z(:,:)  ; order(2)=2
1029   gen(:,:,3)=gen2xy(:,:) ; order(3)=2
1030   gen(:,:,4)=gen3(:,:)   ; order(4)=3
1031   if(iholohedry<=4)ngen=iholohedry-1
1032   if(iholohedry==7)ngen=4
1033 end if
1034
1035!Build the point symmetry operations from generators, in the reduced system
1036!of coordinates defined by axes(:,:)
1037 if(ngen/=0)then
1038   do igen=1,ngen
1039     do isym=1+nptsym,order(igen)*nptsym
1040       jsym=isym-nptsym
1041       do ii=1,3
1042         ptsymrel(:,ii,isym)=gen(:,1,igen)*ptsymrel(1,ii,jsym)+ &
1043&         gen(:,2,igen)*ptsymrel(2,ii,jsym)+ &
1044&         gen(:,3,igen)*ptsymrel(3,ii,jsym)
1045       end do
1046     end do
1047     nptsym=order(igen)*nptsym
1048
1049   end do
1050 end if
1051
1052!--------------------------------------------------------------------------
1053
1054!Transform symmetry matrices in the system defined by rprimd
1055 call symrelrot(nptsym,axes,rprimd,ptsymrel,tolsym)
1056
1057end subroutine symlatt
1058!!***
1059