Line data Source code
1 : !!****m* ABINIT/m_symfind
2 : !! NAME
3 : !! m_symfind
4 : !!
5 : !! FUNCTION
6 : !! Symmetry finder high-level API.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2000-2026 ABINIT group (XG, RC)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_symfind
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_symlist
28 :
29 :
30 : use m_matrix, only : mati3inv, matr3inv, mati3det
31 : use m_symtk, only : chkprimit, symrelrot, symdet, symcharac, holocell, symatm, &
32 : smallprim, print_symmetries, sg_multable, symmetrize_tnons, symmetrize_xred
33 : use m_geometry, only : acrossb, xred2xcart
34 : use m_spgdata, only : getptgroupma, symptgroup, spgdata
35 :
36 : implicit none
37 :
38 : private
39 : !!***
40 :
41 : public :: symfind ! From the symmetries of the Bravais lattice,
42 : ! select those that leave invariant the system, and generate tnons. Not always robust.
43 : public :: symfind_expert ! Wrap symfind to provide robust determination of the symmetries,
44 : ! for which resymmetrization of atomic positions and tnons is needed.
45 : public :: symanal ! Find the space group from the list of symmetries and lattice parameters
46 : public :: symbrav ! Determine the Bravais information from the list of symmetry operations, and the lattice vectors.
47 : public :: symlatt ! Find the Bravais lattice and its symmetry operations (ptsymrel).
48 : ! From the unit cell vectors (rprimd) and the corresponding metric tensor.
49 :
50 : contains
51 : !!***
52 :
53 : !!****f* m_symfind/symfind
54 : !! NAME
55 : !! symfind
56 : !!
57 : !! FUNCTION
58 : !! Symmetry finder.
59 : !! From the symmetries of the Bravais lattice (ptsymrel),
60 : !! select those that leave invariant the system, and generate
61 : !! the corresponding tnons vectors and symafm information.
62 : !! Unlike symfind_expert, does NOT resymmetrize atomic positions and tnons for more robust determination of the symmetries.
63 : !! The algorithm is explained in T.G. Worlton and J.L. Warren, Comp. Phys. Comm. 3, 88 (1972) [[cite:Worlton1972]]
64 : !!
65 : !! INPUTS
66 : !! chrgat(natom) (optional)=target charge for each atom. Not always used, it depends on the value of constraint_kind
67 : !! invardir_red (optional)=reduced coordinates of an invariant direction (only acting with symrel - not tnons)
68 : !! invaraxial_red (optional)=reduced coordinates of an axial vector, similar to invardir_red, but transforms with an additional
69 : !! deteminant factor under symrel operations
70 : !! invar_z (optional)= if 1, the z direction must stay invariant for all symrel applied ;
71 : !! if 2, z must stay invariant and also there cannot be any associated tnons along z.
72 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space
73 : !! msym=default maximal number of symmetries
74 : !! natom=number of atoms in cell.
75 : !! nptsym=number of point symmetries of the Bravais lattice
76 : !! nspden= number of spin-density components. When 4, the three components of spinat are taken into account, instead of only z-component.
77 : !! nucdipmom(3,natom) (optional) array of nuclear dipole moments
78 : !! ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
79 : !! of the Bravais lattice in real space in terms
80 : !! of primitive translations.
81 : !! spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
82 : !! tolsym=tolerance for the symmetries
83 : !! typat(natom)=integer identifying type of atom.
84 : !! use_inversion=1 if inversion and improper rotations can be included in set of symmetries
85 : !! xred(3,natom)=reduced coordinates of atoms in terms of real space
86 : !! primitive translations
87 : !!
88 : !! OUTPUT
89 : !! ierr (optional)=if non-zero, the symmetry operations do not form a group
90 : !! nsym=actual number of symmetries
91 : !! symafm(1:msym)=(anti)ferromagnetic part of nsym symmetry operations
92 : !! symrel(3,3,1:msym)= nsym symmetry operations in real space in terms
93 : !! of primitive translations
94 : !! tnons(3,1:msym)=nonsymmorphic translations for each symmetry (would
95 : !! be 0 0 0 each for a symmorphic space group)
96 : !!
97 : !! SOURCE
98 :
99 18005 : subroutine symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
100 18005 : prtvol, ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
101 35586 : chrgat,ierr,nucdipmom,invardir_red,invaraxial_red,invar_z) ! Optional
102 :
103 : !Arguments ------------------------------------
104 : !scalars
105 : integer,intent(in) :: msym,natom,nptsym,nspden,use_inversion
106 : integer,intent(in) :: prtvol
107 : integer,optional,intent(in) :: invar_z
108 : integer,optional,intent(out) :: ierr
109 : integer,intent(out) :: nsym
110 : real(dp),intent(in) :: tolsym
111 : !arrays
112 : integer,intent(in) :: ptsymrel(3,3,msym),typat(natom)
113 : integer,intent(inout) :: symafm(msym),symrel(3,3,msym) !vz_i
114 : real(dp),intent(in) :: gprimd(3,3),spinat(3,natom),xred(3,natom)
115 : real(dp),optional,intent(in) :: invardir_red(3),invaraxial_red(3),chrgat(natom)
116 : real(dp),optional, intent(in) :: nucdipmom(3,natom)
117 : real(dp),intent(inout) :: tnons(3,msym) !vz_i
118 :
119 : !Local variables-------------------------------
120 : !scalars
121 : integer :: found3,foundcl,iatom,iatom0,iatom1,iatom2,iatom3,iclass,iclass0,ierr_,ii
122 : integer :: isym,jj,kk,natom0,nclass,ntrial,printed,trialafm,trialok
123 : integer :: mm(3,3), detR
124 : real(dp) :: det,diff1,diff2,diff3,diffr1,diffr2,diffr3,ndnorm,nucdipmomcl2,nucdipmomcl20
125 : real(dp) :: spinat2,spinatcl2,spinatcl20,tolsym2
126 : ! TRUE if antiferro symmetries are used with non-collinear magnetism.
127 : integer :: afm_noncoll=1 !For nspden==4. If 1, all symops are permitted ; if 0 symafm must be 1.
128 : !For nspden=4. If noncoll_orthorhombic1, require the symmetry operations to be a subset of the orthorhombic symmetries, except if all spinat=0..
129 : integer :: noncoll_orthorhombic=0
130 : logical :: test_sameabsspin,test_samechrg
131 : logical :: test_samenucdipmom
132 : character(len=500) :: msg
133 : !arrays
134 18005 : integer,allocatable :: class(:,:),natomcl(:),typecl(:)
135 : real(dp) :: diff(3),invardir_red_rot(3),invaraxial_red_rot(3),hand2(3),hand3(3),ndtest(3),rprimd(3,3),spinat0(3),xred0(3)
136 : !real(dp) :: symnucdipmom2(3)
137 : real(dp) :: symnucdipmom2cart(3,3),symnucdipmom2red(3,3)
138 : real(dp) :: symspinat1(3),symspinat2(3),symxred2(3),trialnons(3)
139 18005 : real(dp),allocatable :: chrgat_(:)
140 18005 : real(dp),allocatable :: chrgatcl(:)
141 18005 : real(dp),allocatable :: local_nucdipmom(:,:,:),nucdipmomcl(:,:),nucdipmomred(:,:,:)
142 18005 : real(dp),allocatable :: spinatcl(:,:),spinatred(:,:)
143 : !**************************************************************************
144 :
145 : !write(std_out,'(a)')' m_symfind%symfind : enter '; call flush(std_out)
146 :
147 54015 : ABI_MALLOC(local_nucdipmom,(3,3,natom))
148 596765 : local_nucdipmom(:,:,:) = zero
149 18005 : if(present(nucdipmom)) then
150 193265 : local_nucdipmom(1:3,1,:) = nucdipmom(1:3,:)
151 : end if
152 : !DEBUG
153 : ! if (prtvol>1) msg="remove me later"
154 : ! write(std_out,*)' symfind : enter'
155 : ! call flush(6)
156 : ! write(std_out,*)' ptsymrel matrices are :'
157 : ! do isym=1,nptsym
158 : ! write(std_out,'(i4,4x,9i4)' )isym,ptsymrel(:,:,isym)
159 : ! end do
160 : ! write(std_out,*)' symfind : natom=',natom
161 : ! do iatom=1,natom
162 : ! write(std_out,*)' atom number',iatom
163 : ! write(std_out,*)' typat =',typat(iatom)
164 : ! write(std_out,*)' spinat =',spinat(:,iatom)
165 : ! write(std_out,*)' xred =',xred(:,iatom)
166 : ! if (ANY(ABS(local_nucdipmom)>tol8)) then
167 : ! write(std_out,*)' nucdipmom = ',local_nucdipmom(:,1,iatom)
168 : ! end if
169 : ! end do
170 : ! write(std_out,*)' '
171 : ! call flush(6)
172 : !ENDDEBUG
173 :
174 : ABI_UNUSED(prtvol)
175 :
176 : !Find the number of classes of atoms (type, chrg and spinat must be identical,
177 : !spinat might differ by a sign, if aligned with the z direction, or,
178 : ! type and nucdipmom must be identical)
179 : !natomcl(iclass) will contain the number of atoms in the class
180 : !typecl(iclass) will contain the type of the atoms in the class
181 : !chrgcl(iclass) will contain the charge of the atoms in the class
182 : !spinatcl(1:3,iclass) will contain the spinat of the atoms in the class
183 : !class(1:natomclass(iclass),iclass) will contain the index of the
184 : !atoms belonging to the class
185 72020 : ABI_MALLOC(class,(natom+3,natom))
186 54015 : ABI_MALLOC(natomcl,(natom))
187 36010 : ABI_MALLOC(typecl,(natom))
188 54015 : ABI_MALLOC(chrgat_,(natom))
189 36010 : ABI_MALLOC(chrgatcl,(natom))
190 54015 : ABI_MALLOC(spinatcl,(3,natom))
191 36010 : ABI_MALLOC(nucdipmomcl,(3,natom))
192 :
193 18005 : tolsym2=tolsym**2
194 :
195 62525 : chrgat_(:)=zero
196 18005 : if(present(chrgat))then
197 61661 : chrgat_(:)=chrgat(:)
198 : endif
199 :
200 : ! for each nuclear dipole we need a local right handed coord system, so we can
201 : ! test later for whether a symmetry operation preserves the circulation induced
202 : ! by the dipole
203 62525 : do iatom=1, natom
204 178080 : ndnorm=sqrt(DOT_PRODUCT(local_nucdipmom(1:3,1,iatom),local_nucdipmom(1:3,1,iatom)))
205 :
206 : ! if nuclear dipole has effectively zero size, move on to the next atom
207 44520 : if (ndnorm < tol8) cycle
208 :
209 : ! for testing purposes, we care only about direction so renormalize to unity
210 496 : local_nucdipmom(1:3,1,iatom) = local_nucdipmom(1:3,1,iatom)/ndnorm
211 :
212 : ! make a random vector, each component is (0,1]
213 124 : call random_number(ndtest)
214 :
215 : ! vector 2 is constructed to be orthogonal to original nuclear dipole moment vector
216 124 : call acrossb(local_nucdipmom(1:3,1,iatom),ndtest(1:3),local_nucdipmom(1:3,2,iatom))
217 :
218 : ! vector 3 is orthogonal to 1 and 2, and 1,2,3 form a right-handed set
219 62525 : call acrossb(local_nucdipmom(1:3,1,iatom),local_nucdipmom(1:3,2,iatom),local_nucdipmom(1:3,3,iatom))
220 : end do
221 :
222 : ! need rprimd later to transform back to cart coords
223 18005 : call matr3inv(gprimd,rprimd)
224 :
225 : !write(std_out,'(a)')' m_symfind%symfind : before initialise with the first atom '; call flush(std_out)
226 :
227 : !Initialise with the first atom
228 18005 : nclass=1
229 18005 : natomcl(1)=1
230 18005 : typecl(1)=typat(1)
231 18005 : chrgatcl(1)=chrgat_(1)
232 72020 : spinatcl(:,1)=spinat(:,1)
233 72020 : nucdipmomcl(:,1)=local_nucdipmom(:,1,1)
234 18005 : class(1,1)=1
235 18005 : if(natom>1)then
236 40135 : do iatom=2,natom
237 : ! DEBUG
238 : ! write(std_out,*)' '
239 : ! write(std_out,*)' symfind : examine iatom=',iatom
240 : ! ENDDEBUG
241 38751 : foundcl=0
242 38751 : do iclass=1,nclass
243 : ! Compare the typat, chrg and spinat of atom iatom with existing ones.
244 : ! At this stage, admit either identical spinat, or spin-flip spinat.
245 38751 : if( typat(iatom)==typecl(iclass)) then
246 20612 : test_samechrg= (abs(chrgat_(iatom)-chrgatcl(iclass))<tolsym)
247 20612 : if(nspden/=4)then
248 19968 : test_sameabsspin=(abs(abs(spinat(3,iatom))-abs(spinatcl(3,iclass)))<tolsym)
249 : else if(nspden==4)then
250 644 : spinat2 =spinat(1,iatom)**2+spinat(2,iatom)**2+spinat(3,iatom)**2
251 644 : spinatcl2=spinatcl(1,iclass)**2+spinatcl(2,iclass)**2+spinatcl(3,iclass)**2
252 644 : test_sameabsspin=abs(spinat2-spinatcl2)<tolsym
253 : endif
254 : test_samenucdipmom= &
255 : & abs(local_nucdipmom(1,1,iatom)-nucdipmomcl(1,iclass))<tolsym .and. &
256 : & abs(local_nucdipmom(2,1,iatom)-nucdipmomcl(2,iclass))<tolsym .and. &
257 20612 : & abs(local_nucdipmom(3,1,iatom)-nucdipmomcl(3,iclass))<tolsym
258 : ! note in the following test, m_chkinp/chkinp has already prevented nucdipmom to be
259 : ! nonzero when spinat is nonzero
260 20612 : if( test_samechrg .and. test_sameabsspin .and. test_samenucdipmom ) then
261 : ! write(std_out,*)' symfind : find it belongs to class iclass=',iclass
262 : ! write(std_out,*)' symfind : spinat(:,iatom)=',spinat(:,iatom)
263 : ! write(std_out,*)' symfind : spinatcl(:,iclass)=',spinatcl(:,iclass)
264 : ! write(std_out,*)' symfind : test_sameabsspin=',test_sameabsspin
265 : ! write(std_out,*)' '
266 20530 : natomcl(iclass)=natomcl(iclass)+1
267 20530 : class(natomcl(iclass),iclass)=iatom
268 : foundcl=1
269 : exit
270 : end if
271 : end if
272 : end do
273 : ! If no class with these characteristics exist, create one
274 13620 : if(foundcl==0)then
275 5985 : nclass=nclass+1
276 5985 : natomcl(nclass)=1
277 5985 : typecl(nclass)=typat(iatom)
278 5985 : chrgatcl(nclass)=chrgat_(iatom)
279 23940 : spinatcl(:,nclass)=spinat(:,iatom)
280 23940 : nucdipmomcl(:,nclass)=local_nucdipmom(:,1,iatom)
281 5985 : class(1,nclass)=iatom
282 : end if
283 : end do
284 : end if
285 :
286 : !write(std_out,*)' '
287 : !write(std_out,*)' symfind : found ',nclass,' nclass of atoms'
288 : !do iclass=1,nclass
289 : !write(std_out,*)' class number',iclass
290 : !write(std_out,*)' natomcl =',natomcl(iclass)
291 : !write(std_out,*)' typecl =',typecl(iclass)
292 : !write(std_out,*)' spinatcl=',spinatcl(:,iclass)
293 : !write(std_out,*)' class =',(class(iatom,iclass),iatom=1,natomcl(iclass))
294 : !end do
295 : !write(std_out,*)' '
296 :
297 : !write(std_out,'(a)')' m_symfind%symfind : before select the class '; call flush(std_out)
298 :
299 : !Select the class with the least number of atoms, and non-zero spinat if any
300 : !It is important to select a magnetic class of atom, if any, otherwise
301 : !the determination of the initial (inclusive) set of symmetries takes only
302 : !non-magnetic symmetries, and not both magnetic and non-magnetic ones, see later.
303 : !On the contrary, the chrgat_ data does not play any role, it is invariant upon atomic-centered symmetries
304 18005 : iclass0=1
305 18005 : natom0=natomcl(1)
306 18005 : spinatcl20=spinatcl(1,1)**2+spinatcl(2,1)**2+spinatcl(3,1)**2
307 18005 : nucdipmomcl20=nucdipmomcl(1,1)**2+nucdipmomcl(2,1)**2+nucdipmomcl(3,1)**2
308 18005 : if(nclass>1)then
309 11602 : do iclass=2,nclass
310 5985 : spinatcl2=spinatcl(1,iclass)**2+spinatcl(2,iclass)**2+spinatcl(3,iclass)**2
311 5985 : nucdipmomcl2=nucdipmomcl(1,iclass)**2+nucdipmomcl(2,iclass)**2+nucdipmomcl(3,iclass)**2
312 : if( (natomcl(iclass)<natom0 &
313 : & .and. .not. (spinatcl20>tolsym .and. spinatcl2<tolsym) &
314 : & .and. .not. (nucdipmomcl20>tolsym .and. nucdipmomcl2<tolsym) ) &
315 : & .or. (spinatcl20<tolsym .and. spinatcl2>tolsym) &
316 11602 : & .or. (nucdipmomcl20<tolsym .and. nucdipmomcl2>tolsym)) then
317 5985 : iclass0=iclass
318 5985 : natom0=natomcl(iclass)
319 5985 : spinatcl20=spinatcl2
320 5985 : nucdipmomcl20=nucdipmomcl2
321 : end if
322 : end do
323 : end if
324 :
325 18005 : printed=0
326 :
327 : !If non-collinear spinat have to be used, transfer them in reduced coordinates
328 18005 : if (nspden==4) then
329 1068 : ABI_MALLOC(spinatred,(3,natom))
330 1866 : do iatom=1,natom
331 5862 : do ii=1,3
332 5328 : spinatred(1:3,iatom)=MATMUL(TRANSPOSE(gprimd),spinat(1:3,iatom))
333 : end do
334 : end do
335 : end if
336 :
337 : !write(std_out,*)' '
338 : !write(std_out,*)' symfind : has selected iclass0=',iclass0
339 : !write(std_out,*)' # iatom xred spinat (spinatred if nspden=4) '
340 : !do iatom0=1,natomcl(iclass0)
341 : !iatom=class(iatom0,iclass0)
342 : !if(nspden/=4)then
343 : ! write(std_out,'(2i4,6f10.4)' )iatom0,iatom,xred(:,iatom),spinat(:,iatom)
344 : !else if(nspden==4)then
345 : ! write(std_out,'(2i4,9f10.4)' )iatom0,iatom,xred(:,iatom),spinat(:,iatom),spinatred(:,iatom)
346 : !endif
347 : !end do
348 : !write(std_out,*)' '
349 :
350 : !represent nuclear dipole moments in reduced coords
351 36010 : ABI_MALLOC(nucdipmomred,(3,3,natom))
352 62525 : do iatom=1,natom
353 196085 : do ii=1,3
354 178080 : nucdipmomred(1:3,ii,iatom)=MATMUL(TRANSPOSE(gprimd),local_nucdipmom(1:3,ii,iatom))
355 : end do
356 : end do
357 :
358 : !write(std_out,'(a)')' m_symfind%symfind : before big loop '; call flush(std_out)
359 :
360 : !Big loop over each symmetry operation of the Bravais lattice
361 18005 : nsym=0
362 754453 : do isym=1,nptsym
363 :
364 : !write(std_out,'(a,i4)')' m_symfind%symfind : enter loop isym=',isym; call flush(std_out)
365 :
366 736448 : if(present(invardir_red))then
367 : ! ji: Check whether symmetry operation leaves invardir_red invariant
368 : invardir_red_rot(:) = ptsymrel(:,1,isym)*invardir_red(1) + &
369 : & ptsymrel(:,2,isym)*invardir_red(2) + &
370 2910688 : & ptsymrel(:,3,isym)*invardir_red(3)
371 2910688 : diff(:)=invardir_red(:)-invardir_red_rot(:)
372 727672 : if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 ) cycle
373 : endif
374 :
375 : ! check whether symmetry operation leaves invaraxial_red invariant (axial vector)
376 733488 : if(present(invaraxial_red))then
377 : !! CALL MatrixHyb_getDet(symref_hyb(isym), detR)
378 9421256 : mm(:,:) = ptsymrel(:,:,isym)
379 724712 : call mati3det(mm,detR)
380 : invaraxial_red_rot(:) = detR * (&
381 : ptsymrel(:,1,isym)*invaraxial_red(1) + &
382 : & ptsymrel(:,2,isym)*invaraxial_red(2) + &
383 2898848 : & ptsymrel(:,3,isym)*invaraxial_red(3) )
384 2898848 : diff(:)=invaraxial_red(:)-invaraxial_red_rot(:)
385 :
386 724712 : if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 ) cycle
387 : endif
388 :
389 : !write(std_out,'(a,i4)')' m_symfind%symfind : 1'; call flush(std_out)
390 :
391 : det=ptsymrel(1,1,isym)*ptsymrel(2,2,isym)*ptsymrel(3,3,isym)+&
392 : & ptsymrel(2,1,isym)*ptsymrel(3,2,isym)*ptsymrel(1,3,isym)+&
393 : & ptsymrel(1,2,isym)*ptsymrel(2,3,isym)*ptsymrel(3,1,isym) - &
394 : & (ptsymrel(3,1,isym)*ptsymrel(2,2,isym)*ptsymrel(1,3,isym)+&
395 : & ptsymrel(2,1,isym)*ptsymrel(1,2,isym)*ptsymrel(3,3,isym)+&
396 732704 : & ptsymrel(3,2,isym)*ptsymrel(2,3,isym)*ptsymrel(1,1,isym))
397 732704 : if(use_inversion==0 .and. det==-1) cycle
398 :
399 : !write(std_out,'(a,i4)')' m_symfind%symfind : 2'; call flush(std_out)
400 :
401 : ! jellium slab and spatially varying chemical potential cases:
402 : ! (actually, an inversion symmetry/mirror plane perpendicular to z symmetry operation might still be allowed... TO BE DONE !)
403 722436 : if(present(invar_z))then
404 713668 : if (invar_z/=0) then
405 : ! check whether symmetry operation produce a rotation only in the xy plane
406 : if (ptsymrel(1,3,isym)/=0 .or. ptsymrel(2,3,isym)/=0 .or. &
407 1856 : ptsymrel(3,1,isym)/=0 .or. ptsymrel(3,2,isym)/=0 ) cycle
408 : ! check whether symmetry operation does not change the z
409 832 : if( ptsymrel(3,3,isym)/=1 ) cycle
410 : end if
411 : end if
412 :
413 : !write(std_out,'(a,i4)')' m_symfind%symfind : 3'; call flush(std_out)
414 :
415 : ! If noncoll_orthorhombic=1, require orthorhombic operations of symmetries, except if spinat=0.
416 720996 : if (nspden==4 .and. noncoll_orthorhombic==1)then
417 0 : if(sum(abs(spinat(:,:)))>tol14)then
418 : if( ptsymrel(1,3,isym)/=0 .or. ptsymrel(2,3,isym)/=0 .or. &
419 : & ptsymrel(1,2,isym)/=0 .or. ptsymrel(3,2,isym)/=0 .or. &
420 0 : & ptsymrel(2,1,isym)/=0 .or. ptsymrel(3,2,isym)/=0 ) cycle
421 : endif
422 : endif
423 :
424 : !write(std_out,'(a,i4)')' m_symfind%symfind : 4'; call flush(std_out)
425 :
426 : ! Select a tentative set of associated translations
427 : ! First compute the symmetric of the first atom in the smallest class,
428 : ! using the point symmetry, and also the symmetric of spinat(red).
429 720996 : iatom0=class(1,iclass0)
430 : xred0(:)=ptsymrel(:,1,isym)*xred(1,iatom0)+ &
431 : ptsymrel(:,2,isym)*xred(2,iatom0)+ &
432 2883984 : ptsymrel(:,3,isym)*xred(3,iatom0)
433 720996 : if (nspden/=4) then
434 2831344 : spinat0(:)=spinat(:,iatom0)
435 : else
436 : spinat0(:)=det*(&
437 : ptsymrel(:,1,isym)*spinatred(1,iatom0)+ &
438 : & ptsymrel(:,2,isym)*spinatred(2,iatom0)+ &
439 52640 : & ptsymrel(:,3,isym)*spinatred(3,iatom0))
440 : ! spinat should be treated as an axial vector
441 : ! i.e. the improper part of a symm. op. has no effect on spinat
442 : endif
443 :
444 : !write(std_out,'(a,i4)')' m_symfind%symfind : 5'; call flush(std_out)
445 :
446 : ! From the set of possible images, deduce tentative translations,
447 : ! and magnetic factor then test whether it send each atom on a symmetric one
448 : ntrial=0
449 2012161 : do ii=1,natom0
450 : !write(std_out,'(a,2i4)')' symfind : loop isym,ii=',isym,ii
451 1273160 : iatom1=class(ii,iclass0)
452 :
453 : ! The tentative translation is found
454 5092640 : trialnons(:)=xred(:,iatom1)-xred0(:)
455 : ! Compare the spinat vectors
456 1273160 : if (nspden/=4) then
457 5013808 : symspinat1(:)=spinat(:,iatom1)
458 : else
459 78832 : symspinat1(:)=spinatred(:,iatom1)
460 : end if
461 :
462 : !write(std_out,'(a,6f10.4)')' symspinat1,spinat0=',symspinat1(:),spinat0(:)
463 :
464 1273160 : trialafm=1
465 5092640 : if(sum(abs(symspinat1(:)-spinat0(:)))>tolsym)then
466 16318 : trialafm=-1
467 16318 : if(nspden==4 .and. afm_noncoll==0)cycle
468 65272 : if(sum(abs(symspinat1(:)+spinat0(:)))>tolsym)cycle
469 : endif
470 :
471 5055360 : if(sum(abs(local_nucdipmom(:,1,iatom1)-local_nucdipmom(:,1,iatom0)))>tolsym)then
472 : write(msg,'(3a,3i5)')&
473 0 : 'Problem with matching the nuclear dipole moment within a class.',ch10,&
474 0 : 'isym,iatom0,iatom1=',isym,iatom0,iatom1
475 0 : ABI_ERROR(msg)
476 : end if
477 :
478 : ! jellium slab case: check whether symmetry operation has no translational
479 : ! component along z
480 1263840 : if(present(invar_z))then
481 1250772 : if( invar_z==2 .and. abs(trialnons(3)) > tolsym ) cycle
482 : endif
483 : trialok=1
484 :
485 : ! write(std_out, '(a,i3,a,i3,a,i3,a,3f12.4,i3)') ' Try isym=',isym,' sending iatom0 ',iatom0,' to iatom1 ',iatom1,' with trialnons(:),trialafm =',trialnons(:),trialafm
486 :
487 : ! Loop over all classes, then all atoms in the class,
488 : ! to find whether they have a symmetric
489 2070394 : do iclass=1,nclass
490 4162874 : do jj=1,natomcl(iclass)
491 :
492 3356320 : iatom2=class(jj,iclass)
493 : ! Generate the tentative symmetric position of iatom2
494 : symxred2(:)=ptsymrel(:,1,isym)*xred(1,iatom2)+ &
495 : & ptsymrel(:,2,isym)*xred(2,iatom2)+ &
496 13425280 : & ptsymrel(:,3,isym)*xred(3,iatom2)+ trialnons(:)
497 : ! Generate the tentative symmetric spinat of iatom2
498 3356320 : if (nspden/=4) then
499 13327616 : symspinat2(:)=trialafm*spinat(:,iatom2)
500 : else
501 : symspinat2(:)=trialafm*det*(ptsymrel(:,1,isym)*spinatred(1,iatom2)+ &
502 : & ptsymrel(:,2,isym)*spinatred(2,iatom2)+ &
503 97664 : & ptsymrel(:,3,isym)*spinatred(3,iatom2))
504 : end if
505 : ! Generate the tentative symmetric nucdipmom of iatom2
506 13425280 : do kk = 1, 3
507 : symnucdipmom2red(:,kk)=ptsymrel(:,1,isym)*nucdipmomred(1,kk,iatom2)+ &
508 : & ptsymrel(:,2,isym)*nucdipmomred(2,kk,iatom2)+ &
509 40275840 : & ptsymrel(:,3,isym)*nucdipmomred(3,kk,iatom2)
510 : ! transform back to cart coords for final comparison to nucdipmom
511 164459680 : symnucdipmom2cart(:,kk)=MATMUL(rprimd,symnucdipmom2red(:,kk))
512 : end do
513 :
514 : ! write(std_out,'(a,i4,a,3f8.4,a,3f8.4,a,3f8.4)')&
515 : ! ' Test iatom2=',iatom2,' at xred=',xred(:,iatom2),'. Is sent to',symxred2(:),' with symspinat2=',symspinat2(:)
516 : ! write(std_out,'(a,3f8.4)')' and nucdipmom2=',symnucdipmom2cart(:,1)
517 :
518 : ! Check whether there exists an atom of the same class at the
519 : ! same location, with the correct spinat and nuclear dipole moment circulation
520 28436716 : do kk=1,natomcl(iclass)
521 :
522 27760693 : found3=1
523 27760693 : iatom3=class(kk,iclass)
524 : ! Check the location
525 27760693 : diffr1=xred(1,iatom3)-symxred2(1)
526 27760693 : diff1=diffr1-nint(diffr1)
527 27760693 : if(diff1**2>tolsym2)then
528 : found3=0 ; cycle
529 : else
530 6100390 : diffr2=xred(2,iatom3)-symxred2(2)
531 6100390 : diff2=diffr2-nint(diffr2)
532 6100390 : if(diff2**2>tolsym2)then
533 : found3=0 ; cycle
534 : else
535 3300059 : diffr3=xred(3,iatom3)-symxred2(3)
536 3300059 : diff3=diffr3-nint(diffr3)
537 3300059 : if( (diff1**2+diff2**2+diff3**2) > tolsym**2 )then
538 : found3=0 ; cycle
539 : endif
540 : endif
541 : endif
542 : ! Check the spinat
543 2682689 : if (nspden/=4) then
544 10650892 : diff(:)=spinat(:,iatom3)-symspinat2(:)
545 : else
546 79864 : diff(:)=spinatred(:,iatom3)-symspinat2(:)
547 : end if
548 2682689 : if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 )then
549 : found3=0
550 : cycle
551 : endif
552 : ! Check the nucdipmom
553 : ! hand3 gives original circulation sense of nuclear dipole
554 2682497 : call acrossb(local_nucdipmom(1:3,2,iatom3),local_nucdipmom(1:3,3,iatom3),hand3)
555 :
556 : ! hand2 gives circulation sense of tentative, symmetry equivalent nuclear dipole
557 2682497 : call acrossb(symnucdipmom2cart(1:3,2),symnucdipmom2cart(1:3,3),hand2)
558 :
559 10729988 : diff(:)=hand3(:)-hand2(:)
560 10723804 : if( any(abs(diff)>tolsym) )found3=0
561 :
562 : !diff(:)=symnucdipmom2cart(:,1) - local_nucdipmom(:,1,iatom3)
563 : !if ( any(abs(diff)>tolsym) ) found3=0
564 :
565 676023 : if(found3==1)exit
566 : end do ! End loop over iatom3
567 :
568 4162874 : if(found3==0)then
569 : trialok=0
570 : exit
571 : end if
572 : end do ! End loop over iatom2
573 :
574 2070394 : if(trialok==0)exit
575 : end do ! End loop over all classes
576 :
577 : !write(std_out,*)' For trial isym=',isym,', trialok = ',trialok
578 : !write(std_out,*)' '
579 :
580 2000288 : if(trialok==1)then
581 587817 : nsym=nsym+1
582 587817 : if(nsym>msym)then
583 : write(msg,'(a,i0,2a,i0,4a)')&
584 0 : 'The number of symmetries (including non-symmorphic translations) is: ', nsym, ch10,&
585 0 : 'is larger than maxnsym: ',msym,ch10,&
586 0 : 'Action: increase maxnsym in the input, or take a cell that is primitive, ',ch10,&
587 0 : 'or at least smaller than the present one.'
588 0 : ABI_ERROR(msg)
589 : end if
590 : ntrial=ntrial+1
591 7641621 : symrel(:,:,nsym)=ptsymrel(:,:,isym)
592 587817 : symafm(nsym)=trialafm
593 2351268 : tnons(:,nsym)=trialnons(:)-nint(trialnons(:)-tolsym)
594 : end if
595 :
596 : end do ! End the loop on tentative translations
597 : end do ! End big loop over each symmetry operation of the Bravais lattice
598 :
599 : !write(std_out,'(a)')' m_symfind%symfind : after big loop, will call ABI_FREE '; call flush(std_out)
600 :
601 18005 : ABI_FREE(class)
602 18005 : ABI_FREE(natomcl)
603 18005 : ABI_FREE(chrgat_)
604 18005 : ABI_FREE(chrgatcl)
605 18005 : ABI_FREE(spinatcl)
606 18005 : ABI_FREE(typecl)
607 18005 : ABI_FREE(local_nucdipmom)
608 18005 : ABI_FREE(nucdipmomcl)
609 18005 : ABI_FREE(nucdipmomred)
610 18005 : if (nspden==4) then
611 534 : ABI_FREE(spinatred)
612 : end if
613 :
614 : !write(std_out,'(a,i6)')' m_symfind%symfind : call sg_multable, nsym= ',nsym; call flush(std_out)
615 :
616 : ! The algorithm in sg_multable is still cubic in nsym, so avoid calling it uselessly when nsym is too large
617 18005 : if(present(ierr) .or. nsym<=384)then
618 18005 : call sg_multable(nsym, symafm, symrel, ierr_, tnons=tnons, tnons_tol=tolsym)
619 : else
620 0 : ierr_=0
621 : endif
622 :
623 : !write(std_out,'(a)')' m_symfind%symfind : call print_symmetries, ierr_= ',ierr_; call flush(std_out)
624 :
625 18006 : if (ierr_/=0) call print_symmetries([std_out], nsym, symrel, tnons, symafm)
626 :
627 18005 : if(.not.present(ierr))then
628 9080 : ABI_CHECK(ierr_==0,"Error in group closure")
629 : else
630 8925 : ierr=ierr_
631 : endif
632 :
633 : !write(msg,'(a,I0,es16.6,a)')' symfind : exit, nsym, tolsym=',nsym,tolsym,ch10
634 : !write(msg,'(2a)') trim(msg),' symrel matrices, symafm and tnons are :'
635 : !call wrtout(std_out,msg)
636 : !do isym=1,nsym
637 : ! write(msg,'(i4,4x,3i4,2x,3i4,2x,3i4,4x,i4,4x,3f8.4)' ) isym,symrel(:,:,isym),symafm(isym),tnons(:,isym)
638 : ! call wrtout(std_out,msg)
639 : !end do
640 : !stop
641 :
642 : !write(std_out,'(a)')' m_symfind%symfind : exit '; call flush(std_out)
643 :
644 53803 : end subroutine symfind
645 : !!***
646 :
647 : !!****f* m_symfind/symfind_expert
648 : !! NAME
649 : !! symfind_expert
650 : !!
651 : !! FUNCTION
652 : !! Symmetry finder, with an added layer of robustness compared to symfind,
653 : !! and for which resymmetrization of atomic positions and tnons is needed..
654 : !! From the symmetries of the Bravais lattice (ptsymrel),
655 : !! select those that leave invariant the system, and generate
656 : !! the corresponding tnons vectors and symafm information.
657 : !! Unlike symfind_expert, does NOT resymmetrize atomic positions and tnons for more robust determination of the symmetries.
658 : !! The algorithm is explained in T.G. Worlton and J.L. Warren, Comp. Phys. Comm. 3, 88 (1972) [[cite:Worton1972]]
659 : !!
660 : !! INPUTS
661 : !! chrgat(natom) (optional)=target charge for each atom. Not always used, it depends on the value of constraint_kind
662 : !! invardir_red (optional)=reduced coordinates of an invariant direction (only acting with symrel - not tnons)
663 : !! invaraxial_red (optional)=reduced coordinates of an axial vector, similar to invardir_red, but transforms with an additional
664 : !! deteminant factor under symrel operations
665 : !! invar_z (optional)= if 1, the z direction must stay invariant for all symrel applied ;
666 : !! if 2, z must stay invariant and also there cannot be any associated tnons along z.
667 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space
668 : !! msym=default maximal number of symmetries
669 : !! natom=number of atoms in cell.
670 : !! nptsym=number of point symmetries of the Bravais lattice
671 : !! nspden= number of spin-density components. When 4, the three components of spinat are taken into account, instead of only z-component.
672 : !! nucdipmom(3,natom) (optional) array of nuclear dipole moments
673 : ! pawspnorb=flag: 1 if spin-orbit coupling is activated
674 : !! ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
675 : !! of the Bravais lattice in real space in terms of primitive translations.
676 : !! spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
677 : !! tolsym=tolerance for the symmetries
678 : !! typat(natom)=integer identifying type of atom.
679 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
680 : !!
681 : !! OUTPUT
682 : !! nsym=actual number of symmetries
683 : !! symafm(1:msym)=(anti)ferromagnetic part of nsym symmetry operations
684 : !! symrel(3,3,1:msym)= nsym symmetry operations in real space in terms of primitive translations
685 : !! tnons(3,1:msym)=nonsymmorphic translations for each symmetry (would be 0 0 0 each for a symmorphic space group)
686 : !!
687 : !! SIDE EFFECTS
688 : !! xred(3,natom)=reduced coordinates of atoms in terms of real space
689 : !! primitive translations. Might be changed during the resymmetrization.
690 : !!
691 : !! SOURCE
692 :
693 8924 : subroutine symfind_expert(gprimd,msym,natom,nptsym,nspden,nsym,&
694 8924 : pawspnorb,prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,usepaw,xred,&
695 : chrgat,nucdipmom,invardir_red,invaraxial_red,invar_z) ! Optional - although for the time being all are required ...
696 :
697 : !Arguments ------------------------------------
698 : !scalars
699 : integer,intent(in) :: msym,natom,nptsym,nspden,pawspnorb,usepaw
700 : integer,intent(in) :: prtvol
701 : integer,optional,intent(in) :: invar_z
702 : integer,intent(out) :: nsym
703 : real(dp),intent(in) :: tolsym
704 : !arrays
705 : integer,intent(in) :: ptsymrel(3,3,msym),typat(natom)
706 : integer,intent(inout) :: symafm(msym),symrel(3,3,msym) !vz_i
707 : real(dp),intent(in) :: gprimd(3,3),spinat(3,natom)
708 : real(dp),intent(inout) :: xred(3,natom)
709 : real(dp),optional,intent(in) :: invardir_red(3),invaraxial_red(3),chrgat(natom)
710 : real(dp),optional, intent(in) :: nucdipmom(3,natom)
711 : real(dp),intent(inout) :: tnons(3,msym) !vz_i
712 :
713 : !Local variables-------------------------------
714 : !scalars
715 : integer, save :: print_comment_tolsym=1
716 : integer :: fixed_mismatch,mismatch_fft_tnons
717 : integer :: ierr,isym,use_inversion
718 : character(len=1000) :: msg
719 : !arrays
720 8924 : integer,allocatable :: indsym(:,:,:),symrec(:,:,:)
721 8924 : real(dp),allocatable :: tnons_new(:,:)
722 : !**************************************************************************
723 :
724 : ! write(std_out,*)' m_symfind%symfind_expert : enter '
725 :
726 8924 : use_inversion=1
727 8924 : if (usepaw == 1 .and. (nspden==4.or.pawspnorb>0)) then
728 237 : ABI_COMMENT("Removing inversion and improper rotations from initial space group because of PAW + SOC")
729 : ! MMignolet: PAW can be used with inversion, however it results in seg faults in the dmft code. To enable when this is fixed...
730 237 : use_inversion=0
731 : end if
732 :
733 : ! write(std_out,*)' m_symfind%symfind_expert : before call symfind (1) '
734 :
735 : call symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
736 : prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
737 8924 : chrgat=chrgat,nucdipmom=nucdipmom,ierr=ierr,invardir_red=invardir_red,invaraxial_red=invaraxial_red,invar_z=invar_z)
738 :
739 : ! write(std_out,*)' m_symfind%symfind_expert : after call symfind (1) '
740 :
741 : !If the group closure is not obtained, which should be exceptional, try with a larger tolsym (three times larger)
742 8924 : if(ierr/=0)then
743 1 : ABI_WARNING('Will try to obtain group closure by using a tripled tolsym.')
744 : call symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
745 : prtvol,ptsymrel,spinat,symafm,symrel,tnons,three*tolsym,typat,use_inversion,xred,&
746 1 : chrgat=chrgat,nucdipmom=nucdipmom,ierr=ierr,invardir_red=invardir_red,invaraxial_red=invaraxial_red,invar_z=invar_z)
747 1 : ABI_CHECK(ierr==0,"Error in group closure")
748 1 : ABI_WARNING('Succeeded to obtain group closure by using a tripled tolsym.')
749 : endif
750 :
751 : ! If the tolerance on symmetries is bigger than 1.e-8, symmetrize tnons for gliding or screw operations,
752 : ! symmetrize the atomic positions and recompute the symmetry operations
753 8924 : if(tolsym>1.00001e-8)then
754 :
755 8868 : call symmetrize_tnons(nsym,symrel,tnons,tolsym)
756 35472 : ABI_MALLOC(indsym,(4,natom,nsym))
757 26604 : ABI_MALLOC(symrec,(3,3,nsym))
758 298124 : do isym=1,nsym
759 298124 : call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
760 : end do
761 8868 : call symatm(indsym,natom,nsym,symrec,tnons,tolsym,typat,xred)
762 8868 : call symmetrize_xred(natom,nsym,symrel,tnons,xred,indsym=indsym)
763 8868 : ABI_FREE(indsym)
764 8868 : ABI_FREE(symrec)
765 :
766 8868 : if(print_comment_tolsym==1)then
767 : write(msg,'(a,es12.3,18a)')&
768 1360 : 'The tolerance on symmetries =',tolsym,' is bigger than 1.0e-8.',ch10,&
769 1360 : 'In order to avoid spurious effects, the atomic coordinates have been',ch10,&
770 1360 : 'symmetrized before storing them in the dataset internal variable.',ch10,&
771 1360 : 'So, do not be surprised by the fact that your input variables (xcart, xred, ...)',ch10,&
772 1360 : 'do not correspond exactly to the ones echoed by ABINIT, the latter being used to do the calculations.',ch10,&
773 1360 : 'This is not a problem per se.',ch10,&
774 1360 : 'Still, in order to avoid this symmetrization (e.g. for specific debugging/development),',&
775 1360 : ' decrease tolsym to 1.0e-8 or lower,',ch10,&
776 1360 : 'or (much preferred) use input primitive vectors that are accurate to better than 1.0e-8.',ch10,&
777 2720 : 'This message will only be printed once, even if there are other datasets where tolsym is bigger than 1.0e-8.'
778 1360 : ABI_COMMENT(msg)
779 1360 : print_comment_tolsym=0
780 : endif
781 :
782 : !write(std_out,*)' m_symfind%symfind_expert : before call symfind (3) '
783 :
784 : call symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
785 : prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
786 8868 : chrgat=chrgat,nucdipmom=nucdipmom,invardir_red=invardir_red,invaraxial_red=invaraxial_red,invar_z=invar_z)
787 :
788 : ! write(std_out,*)' m_symfind%symfind_expert : after call symfind (3) '
789 :
790 : !Needs one more resymmetrization, for the tnons
791 26604 : ABI_MALLOC(tnons_new,(3,nsym))
792 :
793 : call symmetrize_xred(natom,nsym,symrel,tnons,xred,&
794 8868 : fixed_mismatch=fixed_mismatch,mismatch_fft_tnons=mismatch_fft_tnons,tnons_new=tnons_new,tolsym=tolsym)
795 1165892 : tnons(:,1:nsym)=tnons_new(:,:)
796 26604 : ABI_FREE(tnons_new)
797 : end if ! tolsym >1.00001e-8
798 :
799 : ! write(std_out,*)' m_symfind%symfind_expert : exit '
800 :
801 8924 : end subroutine symfind_expert
802 : !!***
803 :
804 : !!****f* m_symfind/symanal
805 : !! NAME
806 : !! symanal
807 : !!
808 : !! FUNCTION
809 : !! Find the space group, Bravais lattice, including Shubnikov characteristics
810 : !! from the list of symmetries (including magnetic characteristics), and lattice parameters
811 : !! Warning: the recognition of the space group might not yet work for the
812 : !! Shubnikov group of type IV
813 : !!
814 : !! INPUTS
815 : !! chkprim= if 1 then stop if the cell is not primitive
816 : !! msym=default maximal number of symmetries
817 : !! nsym=actual number of symmetries
818 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
819 : !! symafm(1:msym)=(anti)ferromagnetic part of symmetry operations
820 : !! symrel(3,3,1:msym)=symmetry operations in real space in terms
821 : !! of primitive translations
822 : !! tnons(3,1:msym)=nonsymmorphic translations for symmetry operations
823 : !! tolsym=tolerance for the symmetry operations
824 : !! [verbose]= if true, will list the symmetry operation labels
825 : !!
826 : !! OUTPUT
827 : !! bravais(11)=characteristics of Bravais lattice (see symlatt.F90)
828 : !! genafm(3)=magnetic translation generator (in case of Shubnikov group type IV)
829 : !! ptgroupma = magnetic point group number
830 : !! spgroup=symmetry space group
831 : !!
832 : !! SOURCE
833 :
834 11934 : subroutine symanal(bravais,chkprim,genafm,msym,nsym,ptgroupma,rprimd,spgroup,symafm,symrel,tnons,tolsym, &
835 : verbose) ! optional
836 :
837 : !Arguments ------------------------------------
838 : !scalars
839 : integer,intent(in) :: chkprim,msym,nsym
840 : integer,intent(out) :: ptgroupma,spgroup
841 : real(dp),intent(in) :: tolsym
842 : logical,optional,intent(in) :: verbose
843 : !arrays
844 : integer,intent(out) :: bravais(11)
845 : integer,intent(in) :: symafm(msym),symrel(3,3,msym)
846 : real(dp),intent(in) :: rprimd(3,3), tnons(3,msym)
847 : real(dp),intent(out) :: genafm(3)
848 :
849 : !Local variables-------------------------------
850 : !scalars
851 : integer, parameter :: maxsym=192
852 : ! In this routine, maxsym is used either to determine the ptsymrel from the rprimd (routine symlatt),
853 : ! so, it might be up to 192 = 4*48 for FCC, and also to define the maximum number of symmetry operation labels,
854 : ! but only in case the cell is primitive, which gives the same upper bound. Thus in this routine, msym might
855 : ! be equal to nsym.
856 : integer :: iholohedry_nomagn,isym,isym_nomagn,multi, nptsym,nsym_nomagn,shubnikov
857 : logical :: verbose_
858 : character(len=5) :: ptgroup,ptgroupha
859 : character(len=500) :: msg
860 : !arrays
861 : integer :: identity(3,3)
862 11934 : integer,allocatable :: ptsymrel(:,:,:),symrel_nomagn(:,:,:)
863 11934 : real(dp),allocatable :: tnons_nomagn(:,:)
864 : character(len=128) :: labels(maxsym)
865 : ! *************************************************************************
866 :
867 : !write(std_out,*)' symanal : enter'
868 : !write(std_out,*)' symanal : chkprim =',chkprim
869 : !write(std_out,*)' symanal : nsym=',nsym
870 : !do isym=1,nsym
871 : ! write(std_out,*)' symanal : symrel=',symrel(1:3,1:3,isym)
872 : !enddo
873 :
874 11934 : verbose_=.false.
875 11934 : if (present(verbose)) verbose_=verbose
876 :
877 : !This routine finds the Bravais characteristics, without actually looking at the symmetry operations.
878 11934 : ABI_MALLOC(ptsymrel,(3,3,maxsym))
879 11934 : call symlatt(bravais,dev_null,maxsym,nptsym,ptsymrel,rprimd,tolsym)
880 11934 : ABI_FREE(ptsymrel)
881 :
882 : ! Check whether the cell is primitive or not.
883 11934 : call chkprimit(chkprim,multi,nsym,symafm,symrel)
884 :
885 11934 : spgroup=0 ; ptgroupma=0 ; genafm(:)=zero
886 :
887 11934 : if (multi>1) then ! Modify bravais if the cell is not primitive ; no determination of the space group
888 289 : bravais(1)=-bravais(1)
889 : else
890 :
891 : ! The cell is primitive, so that the space group can be
892 : ! determined. Need to distinguish Fedorov and Shubnikov groups.
893 : ! Do not distinguish Shubnikov types I and II.
894 : ! Also identify genafm, in case of Shubnikov type IV
895 11645 : identity(:,:)=reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/))
896 11645 : shubnikov=1
897 297715 : do isym=1,nsym
898 297715 : if(symafm(isym)==-1)then
899 5740 : shubnikov=3
900 74620 : if(sum(abs(symrel(:,:,isym)-identity(:,:)))==0)then
901 607 : shubnikov=4
902 2428 : genafm(:)=tnons(:,isym)
903 : !write(std_out,*)' isym=',isym
904 : !write(std_out,*)' symrel(:,:,isym)',symrel(:,:,isym)
905 : !write(std_out,*)' tnons(:,isym)',tnons(:,isym)
906 : !write(std_out,*)' symafm(isym)',symafm(isym)
907 : exit
908 : end if
909 : end if
910 : end do
911 :
912 11645 : if(shubnikov/=1)then
913 1505 : if(shubnikov==3)write(msg, '(a)' )' Shubnikov space group type III'
914 1505 : if(shubnikov==4)write(msg, '(a)' )' Shubnikov space group type IV'
915 1505 : call wrtout(std_out,msg)
916 : end if
917 :
918 11645 : if(shubnikov==1 .or. shubnikov==3)then
919 : ! Find the correct Bravais characteristics and point group
920 : ! Should also be used for Shubnikov groups of type IV ...
921 11038 : call symbrav(bravais,msym,nsym,ptgroup,rprimd,symrel,tolsym)
922 :
923 : ! Find the space group
924 11038 : call symspgr(bravais,labels,nsym,spgroup,symrel,tnons,tolsym)
925 :
926 11038 : if (verbose_) then
927 2275 : do isym=1,nsym
928 2203 : write(msg,'(a,i3,2a)')' symanal : the symmetry operation no. ',isym,' is ',trim(labels(isym))
929 2275 : call wrtout(std_out,msg)
930 : enddo
931 : endif
932 :
933 : end if
934 :
935 11038 : if (shubnikov/=1) then
936 :
937 : ! Determine nonmagnetic symmetry operations
938 1505 : nsym_nomagn=nsym/2
939 4515 : ABI_MALLOC(symrel_nomagn,(3,3,nsym_nomagn))
940 4515 : ABI_MALLOC(tnons_nomagn,(3,nsym_nomagn))
941 1505 : isym_nomagn=0
942 23661 : do isym=1,nsym
943 23661 : if(symafm(isym)==1)then
944 11078 : isym_nomagn=isym_nomagn+1
945 144014 : symrel_nomagn(:,:,isym_nomagn)=symrel(:,:,isym)
946 44312 : tnons_nomagn(:,isym_nomagn)=tnons(:,isym)
947 : end if
948 : end do
949 :
950 1505 : if (shubnikov==3) then
951 :
952 : ! write(std_out,*)' symanal : will enter symbrav with halved symmetry set'
953 : ! write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
954 : ! do isym=1,nsym_nomagn
955 : ! write(std_out,'(i3,2x,9i3,3es12.2,i3)')isym,symrel_nomagn(:,:,isym),tnons_nomagn(:,isym)
956 : ! end do
957 :
958 : ! Find the point group of the halved symmetry set
959 898 : call symptgroup(iholohedry_nomagn,nsym_nomagn,ptgroupha,symrel_nomagn)
960 :
961 : ! Deduce the magnetic point group (ptgroupma) from ptgroup and ptgroupha
962 898 : call getptgroupma(ptgroup,ptgroupha,ptgroupma)
963 :
964 607 : else if(shubnikov==4)then
965 :
966 : ! Find the Fedorov space group of the halved symmetry set
967 607 : call symspgr(bravais,labels,nsym_nomagn,spgroup,symrel_nomagn,tnons_nomagn,tolsym)
968 :
969 : ! The magnetic translation generator genafm has already been determined
970 : ! write(std_out,*)' genafm =',genafm, ' spgroup=',spgroup
971 607 : if (verbose_) then
972 0 : write(msg, '(a)' )' Select only the non-magnetic symmetry operations '
973 0 : call wrtout(std_out,msg)
974 :
975 0 : do isym=1,nsym
976 0 : if(symafm(isym)==1)then
977 0 : isym_nomagn=isym_nomagn+1
978 0 : write(msg,'(a,i3,2a)')' symspgr : the symmetry operation no. ',isym,' is ',trim(labels(isym_nomagn))
979 0 : call wrtout(std_out,msg)
980 : endif
981 : enddo
982 : endif
983 :
984 : end if
985 :
986 1505 : ABI_FREE(symrel_nomagn)
987 1505 : ABI_FREE(tnons_nomagn)
988 : end if ! Shubnikov groups
989 :
990 : end if
991 :
992 11934 : end subroutine symanal
993 : !!***
994 :
995 : !!****f* m_symfind/symbrav
996 : !! NAME
997 : !! symbrav
998 : !!
999 : !! FUNCTION
1000 : !! From the list of symmetry operations, and the lattice vectors,
1001 : !! determine the Bravais information (including the holohedry, the centering,
1002 : !! the coordinate of the primitive vectors in the conventional vectors),
1003 : !! as well as the point group.
1004 : !!
1005 : !! INPUTS
1006 : !! msym=dimension of symrel
1007 : !! nsym=actual number of symmetries
1008 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1009 : !! symrel(3,3,msym)=symmetry operations in real space in terms
1010 : !! of primitive translations
1011 : !! tolsym=tolerance for the symmetries
1012 : !!
1013 : !! OUTPUT
1014 : !! bravais(11): bravais(1)=iholohedry
1015 : !! bravais(2)=center
1016 : !! bravais(3:11)=coordinates of rprimd in the axes
1017 : !! of the conventional bravais lattice (*2 if center/=0)
1018 : !! ptgroup=symmetry point group
1019 : !! [axis(3)]=Invariant axis in the conventional vector coordinates
1020 : !! Set to (/0,0,0/) if the lattice belongs to the same holohedry as the lattice+atoms (+electric field + ...).
1021 : !!
1022 : !! SOURCE
1023 :
1024 11038 : subroutine symbrav(bravais,msym,nsym,ptgroup,rprimd,symrel,tolsym,axis)
1025 :
1026 : !Arguments ------------------------------------
1027 : !scalars
1028 : integer,intent(in) :: msym,nsym
1029 : real(dp),intent(in) :: tolsym
1030 : character(len=5),intent(out) :: ptgroup
1031 : !arrays
1032 : integer,intent(in) :: symrel(3,3,msym)
1033 : integer,optional,intent(out) :: axis(3)
1034 : integer,intent(out) :: bravais(11)
1035 : real(dp),intent(in) :: rprimd(3,3)
1036 :
1037 : !Local variables-------------------------------
1038 : !scalars
1039 : integer :: iaxis,ii,bravais1now,ideform,iholohedry,invariant,isym
1040 : integer :: jaxis,next_stage,nptsym,problem,maxsym
1041 : integer, parameter :: naxes_ortho=22, naxes_hexa=7
1042 : real(dp) :: norm,scprod
1043 : character(len=500) :: msg
1044 : !arrays
1045 : integer :: identity(3,3),axis_trial(3),hexa_axes(3,naxes_hexa),ortho_axes(3,naxes_ortho)
1046 11038 : integer,allocatable :: ptsymrel(:,:,:),symrelconv(:,:,:)
1047 : real(dp) :: axes(3,3),axis_cart(3),axis_red(3)
1048 : real(dp) :: rprimdconv(3,3),rprimdtry(3,3),rprimdnow(3,3)
1049 : real(dp) :: rprimdconv_invt(3,3)
1050 : !**************************************************************************
1051 :
1052 : !DEBUG
1053 : !write(std_out,*)' symbrav : enter '
1054 : !call flush(std_out)
1055 : !ENDDEBUG
1056 :
1057 : identity(:,:)=0
1058 11038 : identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
1059 :
1060 11038 : ortho_axes(:,:)=0
1061 11038 : ortho_axes(1,1)=1
1062 11038 : ortho_axes(2,2)=1
1063 11038 : ortho_axes(3,3)=1
1064 44152 : ortho_axes(:,4)=(/0,1,1/)
1065 44152 : ortho_axes(:,5)=(/1,0,1/)
1066 44152 : ortho_axes(:,6)=(/1,1,0/)
1067 44152 : ortho_axes(:,7)=(/0,1,-1/)
1068 44152 : ortho_axes(:,8)=(/-1,0,1/)
1069 44152 : ortho_axes(:,9)=(/1,-1,0/)
1070 44152 : ortho_axes(:,10)=(/0,1,2/)
1071 44152 : ortho_axes(:,11)=(/2,0,1/)
1072 44152 : ortho_axes(:,12)=(/1,2,0/)
1073 44152 : ortho_axes(:,13)=(/1,1,1/)
1074 44152 : ortho_axes(:,14)=(/-1,1,1/)
1075 44152 : ortho_axes(:,15)=(/1,-1,1/)
1076 44152 : ortho_axes(:,16)=(/1,1,-1/)
1077 44152 : ortho_axes(:,17)=(/2,1,1/)
1078 44152 : ortho_axes(:,18)=(/1,2,1/)
1079 44152 : ortho_axes(:,19)=(/1,1,2/)
1080 44152 : ortho_axes(:,20)=(/2,1,-1/)
1081 44152 : ortho_axes(:,21)=(/-1,2,1/)
1082 44152 : ortho_axes(:,22)=(/1,-1,2/)
1083 :
1084 11038 : hexa_axes(:,:)=0
1085 11038 : hexa_axes(1,1)=1
1086 11038 : hexa_axes(2,2)=1
1087 11038 : hexa_axes(3,3)=1
1088 44152 : hexa_axes(:,4)=(/1,-1,0/)
1089 44152 : hexa_axes(:,5)=(/2,1,0/)
1090 44152 : hexa_axes(:,6)=(/1,1,0/)
1091 44152 : hexa_axes(:,7)=(/1,2,0/)
1092 :
1093 : !Determine the point group from the list of symmetry operations.
1094 : !Also determine the holohedry, up to one undeterminacy : hR versus hP
1095 11038 : call symptgroup(iholohedry,nsym,ptgroup,symrel)
1096 :
1097 : !DEBUG
1098 : !write(std_out,*)' symbrav, after symptgroup: nsym=',nsym
1099 : !call flush(std_out)
1100 : !write(std_out,*)' symbrav: symrel='
1101 : !do isym=1,nsym
1102 : ! write(std_out,'(9i4)')symrel(:,:,isym)
1103 : !enddo
1104 : !write(std_out,*)' symbrav: iholohedry=',iholohedry
1105 : !call flush(std_out)
1106 : !ENDDEBUG
1107 :
1108 : !Loop over trial deformations
1109 : !This is needed in case the Bravais lattice determination from the lattice vectors
1110 : !has a higher holohedry than the real one, in which the symmetry
1111 : !operations for the atoms (or electric field, etc) are taken into account
1112 11038 : iaxis=0
1113 11038 : invariant=0
1114 11038 : next_stage=0
1115 11038 : rprimdnow(:,:)=rprimd(:,:)
1116 11038 : rprimdtry(:,:)=rprimd(:,:)
1117 33114 : ABI_MALLOC(symrelconv,(3,3,nsym))
1118 :
1119 : !At most will have to try naxes_ortho*5 deformations (naxes_ortho axes, five stages)
1120 : !First, test whether the current recognition of Bravais lattice is problematic (iholohedry differs from bravais(1)).
1121 : !Then, if there is a problem, test different deformations of rprimd, one after the other.
1122 : !For each try, there is a new rprimdtry from the different set of deformation, that is generated later in the loop
1123 : !Also bravais1now and rprimdnow might changed (and progressively lowered).
1124 : !The latter change induces at most 5 stages for the computation (cubic->tetragonal->orthorhombic->monoclinic->triclinic).
1125 : !After an upgrade of bravais1now and rprimdnow, one has to restart the full set of deformations.
1126 : !Not sure that this procedure resolves all cases, but seems to work on >40000 inaccurate POSCAR files.
1127 15432 : do ideform=1,naxes_ortho*5
1128 :
1129 15432 : maxsym=max(192,msym)
1130 46296 : ABI_MALLOC(ptsymrel,(3,3,maxsym))
1131 15432 : call symlatt(bravais,std_out,maxsym,nptsym,ptsymrel,rprimdtry,tolsym)
1132 15432 : ABI_FREE(ptsymrel)
1133 :
1134 : !DEBUG
1135 : !write(std_out,*)' symbrav: inside loop with ideform,iaxis=',ideform,iaxis
1136 : !write(std_out,'(a,9f12.6)')' rprimdtry=',rprimdtry(:,:)
1137 : !write(std_out,'(a,2i4)')' bravais(1:2)=',bravais(1:2)
1138 : !call flush(std_out)
1139 : !ENDDEBUG
1140 :
1141 :
1142 : ! Examine the agreement with bravais(1)
1143 : ! Warning : might change Bravais lattice hR to hP, if hexagonal axes
1144 15432 : problem=0
1145 23334 : select case (bravais(1))
1146 : case (7)
1147 7902 : if(iholohedry<6)problem=1
1148 7902 : if(iholohedry==6)problem=2
1149 : case (6)
1150 668 : if(iholohedry<4)problem=1
1151 668 : if(iholohedry==7 .or. iholohedry==4)problem=2
1152 : ! Here, change hR into hP
1153 668 : if(iholohedry==5)iholohedry=6
1154 : case (5)
1155 806 : if(iholohedry<4)problem=1
1156 806 : if(iholohedry==7 .or. iholohedry==6 .or. iholohedry==4)problem=2
1157 : case (4)
1158 2040 : if(iholohedry<4)problem=1
1159 2040 : if(iholohedry>4)problem=2
1160 : case (3)
1161 898 : if(iholohedry<3)problem=1
1162 898 : if(iholohedry>3)problem=2
1163 : case (2)
1164 1721 : if(iholohedry<2)problem=1
1165 1721 : if(iholohedry>2)problem=2
1166 : case (1)
1167 15432 : if(iholohedry>1)problem=2
1168 : end select
1169 :
1170 : ! This is the usual situation, in which the lattice belong to the same holohedry
1171 : ! as the lattice+atoms (+electric field + ...)
1172 15432 : if(problem==0)exit
1173 :
1174 4394 : if(problem==2)then
1175 0 : if(iaxis==0)then
1176 : write(msg, '(3a,i3,3a,i3,7a)' )&
1177 0 : 'The Bravais lattice determined only from the primitive',ch10,&
1178 0 : 'vectors (rprim or angdeg), bravais(1)=',bravais(1),', is not compatible',ch10,&
1179 0 : 'with the real one, iholohedry=',iholohedry,', obtained by taking into',ch10,&
1180 0 : 'account the symmetry operations. This might be due to an insufficient',ch10,&
1181 0 : 'number of digits in the specification of rprim (at least 10),',ch10,&
1182 0 : 'or to an erroneous rprim or angdeg. If this is not the case, then ...'
1183 0 : ABI_BUG(msg)
1184 : end if
1185 0 : if(iaxis==1)then
1186 : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1187 0 : 'Could not succeed to determine the bravais lattice',ch10,&
1188 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1189 0 : 'bravais(1)=',bravais(1),ch10,&
1190 0 : 'iholohedry=',iholohedry
1191 0 : ABI_BUG(msg)
1192 : end if
1193 : ! Try to increase tolsym to find the Bravais lattice.
1194 : maxsym=max(192,msym)
1195 0 : ABI_MALLOC(ptsymrel,(3,3,maxsym))
1196 : !DEBUG
1197 : ! write(6,*)' symbrav : will call symlatt, 3*tolsym=',3*tolsym
1198 : !ENDDEBUG
1199 0 : call symlatt(bravais,std_out,maxsym,nptsym,ptsymrel,rprimdtry,3*tolsym)
1200 0 : ABI_FREE(ptsymrel)
1201 0 : if(bravais(1)==iholohedry)then
1202 : ! Succeeded
1203 : exit
1204 : else
1205 : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1206 0 : 'Could not succeed to determine the bravais lattice, even after considering a larger tolsym',ch10,&
1207 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1208 0 : 'bravais(1)=',bravais(1),ch10,&
1209 0 : 'iholohedry=',iholohedry
1210 0 : ABI_BUG(msg)
1211 : end if
1212 : end if
1213 :
1214 4394 : if(problem==1)then ! One is left with the problem=1 case, basically iholohedry is lower than bravais(1)
1215 4394 : if(iaxis==0)then
1216 : write(msg, '(a,a,a,i3,a,a,a,i3,a,a,a)' )&
1217 2878 : 'The Bravais lattice determined only from the primitive',ch10,&
1218 2878 : 'vectors, bravais(1)=',bravais(1),', is more symmetric',ch10,&
1219 2878 : 'than the real one, iholohedry=',iholohedry,', obtained by taking into',ch10,&
1220 5756 : 'account the atomic positions. Start deforming the primitive vector set.'
1221 2878 : ABI_COMMENT(msg)
1222 : next_stage=1
1223 : else if(iaxis/=0)then
1224 1516 : if(bravais(1)<bravais1now)then
1225 : write(msg, '(3a,i3,3a,i3,2a)' )&
1226 1301 : 'The Bravais lattice determined from modified primitive',ch10,&
1227 1301 : 'vectors, bravais(1)=',bravais(1),', has a lower symmetry than before,',ch10,&
1228 1301 : 'but is still more symmetric than the real one, iholohedry=',iholohedry,ch10,&
1229 2602 : 'obtained by taking into account the atomic positions.'
1230 1301 : ABI_COMMENT(msg)
1231 : next_stage=1
1232 215 : else if(iaxis==1)then
1233 : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1234 0 : 'Could not succeed to determine the bravais lattice',ch10,&
1235 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1236 0 : 'bravais(1)=',bravais(1),ch10,&
1237 0 : 'iholohedry=',iholohedry
1238 0 : ABI_BUG(msg)
1239 : end if
1240 : end if
1241 : end if ! problem==1
1242 :
1243 : ! One is here when problem=1 (iholohedry < bravais(1)) and either
1244 : ! - iaxis=0 (no deformation has been tried yet),
1245 : ! - some deformation iaxis has been tried giving bravais(1), but iholohedry < bravais(1) < bravais1now
1246 : ! - some deformation iaxis has been tried giving bravais(1), but iholohedry < bravais(1) = bravais1now and iaxis/=1 .
1247 : ! Also, note that next_stage is still 0 when bravais(1)=bravais1now .
1248 : ! The loop has been ended (so, the search failed) when bravais(1)=bravais1now and iaxis==1.
1249 :
1250 : if(next_stage==1)then
1251 4179 : bravais1now=bravais(1)
1252 4179 : rprimdnow(:,:)=rprimdtry(:,:)
1253 : ! Generate the symmetry operations in the conventional vector coordinates
1254 16716 : rprimdconv(:,1)=bravais(3:5)
1255 16716 : rprimdconv(:,2)=bravais(6:8)
1256 16716 : rprimdconv(:,3)=bravais(9:11)
1257 4179 : axes(:,:)=zero
1258 4179 : axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one
1259 209397 : symrelconv(:,:,1:nsym)=symrel(:,:,1:nsym)
1260 4179 : call symrelrot(nsym,rprimdconv,axes,symrelconv,tolsym)
1261 4179 : if(bravais(1)/=6)then
1262 4162 : iaxis=naxes_ortho+1
1263 : else
1264 17 : iaxis=naxes_hexa+1
1265 : end if
1266 : next_stage=0
1267 : !DEBUG
1268 : ! write(std_out,*)' symbrav: next stage, bravais(1), bravais(2) and symrelconv'
1269 : ! write(std_out,'(a,2i4)')' bravais(1:2)=',bravais(1:2)
1270 : ! write(std_out,'(a,9f12.6)')' rprimdconv=',rprimdconv(:,:)
1271 : ! do isym=1,nsym
1272 : ! write(std_out,'(9i4)')symrelconv(:,:,isym)
1273 : ! enddo
1274 : ! call flush(std_out)
1275 : !ENDDEBUG
1276 : end if
1277 :
1278 : ! Go to the next iaxis that will be left invariant
1279 4394 : iaxis=iaxis-1
1280 24832 : do jaxis=iaxis,1,-1
1281 24832 : if(bravais(1)/=6)then
1282 99252 : axis_trial(:)=ortho_axes(:,jaxis)
1283 : else
1284 76 : axis_trial(:)=hexa_axes(:,jaxis)
1285 : end if
1286 : ! DEBUG
1287 : ! write(std_out,*)' symbrav : ixaxis, trial jaxis=',iaxis,jaxis
1288 : ! write(std_out,*)' axis_trial=',axis_trial
1289 : ! ENDDEBUG
1290 24832 : invariant=1
1291 : ! Examine whether all symmetry operations leave the axis invariant (might be reversed, though)
1292 224594 : do isym=1,nsym
1293 1797858 : if(sum(abs(matmul(symrelconv(:,:,isym),axis_trial)+(-axis_trial(:))))/=0 .and. &
1294 7016502 : sum(abs(matmul(symrelconv(:,:,isym),axis_trial)+axis_trial(:)))/=0 )invariant=0
1295 : end do
1296 24832 : if(invariant==1)then
1297 4394 : iaxis=jaxis
1298 : !DEBUG
1299 : ! write(msg, '(2a,i3)' )ch10,' symbrav : found invariant axis, jaxis=',jaxis
1300 : ! call wrtout(std_out,msg)
1301 : !ENDDEBUG
1302 4394 : exit
1303 : end if
1304 : end do
1305 :
1306 4394 : if(invariant==0)then
1307 : ! Not a single axis was invariant with respect to all operations ?!
1308 : ! do isym=1,nsym; write(std_out, '(a,10i4)' )' isym,symrelconv=',isym,symrelconv(:,:,isym); enddo
1309 : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1310 0 : 'Could not succeed to determine the bravais lattice (not a single invariant)',ch10,&
1311 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1312 0 : 'bravais(1)=',bravais(1),ch10,&
1313 0 : 'iholohedry=',iholohedry
1314 0 : ABI_BUG(msg)
1315 : end if
1316 :
1317 4394 : call matr3inv(rprimdconv,rprimdconv_invt)
1318 : axis_red(:)=axis_trial(1)*rprimdconv_invt(1,:)+ &
1319 : & axis_trial(2)*rprimdconv_invt(2,:)+ &
1320 17576 : & axis_trial(3)*rprimdconv_invt(3,:)
1321 : axis_cart(:)=axis_red(1)*rprimdnow(:,1)+ &
1322 : & axis_red(2)*rprimdnow(:,2)+ &
1323 17576 : & axis_red(3)*rprimdnow(:,3)
1324 17576 : norm=sum(axis_cart(:)**2)
1325 : !DEBUG
1326 : ! write(6,*)' axis_trial =',axis_trial
1327 : ! write(6,*)' axis_red =',axis_red
1328 : ! write(6,*)' axis_cart =',axis_cart
1329 : ! write(6,*)' rprimdnow=',rprimdnow
1330 : !ENDDEBUG
1331 : ! Expand by a uniform, quite arbitrary, dilatation, along the invariant axis
1332 : ! Note : make these dilatation different, according to ideform
1333 : ! XG 20151221 : Still, the interplay between the size of the deformation and the tolsym is not easy to address.
1334 : ! Indeed the deformation must be sufficiently large to be perceived by symlatt as a real breaking of the
1335 : ! symmetry of the lattice. In order to deal with all the small values od tolsym, it has been set at a minimum of tol3,
1336 : ! but it must also be larger than tolsym. Moreover, for some axis choice, the deformation is not aligned with the axis, decreasing
1337 : ! the effective deformation length. An additional factor of three is thus included, actually increased to six just to be sure...
1338 44046 : do ii=1,3
1339 13182 : scprod=axis_cart(1)*rprimdnow(1,ii)+axis_cart(2)*rprimdnow(2,ii)+axis_cart(3)*rprimdnow(3,ii)
1340 57122 : rprimdtry(:,ii)=rprimdnow(:,ii)+ideform*(max(tol3,six*tolsym)-tol6)*scprod/norm*axis_cart(:)
1341 : end do
1342 :
1343 : !DEBUG
1344 : ! write(6,*)' rprimdtry=',rprimdtry
1345 : !ENDDEBUG
1346 :
1347 : end do ! ideform
1348 :
1349 11038 : if(bravais(1)/=iholohedry)then
1350 : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1351 0 : 'Despite efforts, Could not succeed to determine the bravais lattice :',ch10,&
1352 0 : 'bravais(1)=',bravais(1),ch10,&
1353 0 : 'iholohedry=',iholohedry
1354 0 : ABI_BUG(msg)
1355 : end if
1356 :
1357 11038 : ABI_FREE(symrelconv)
1358 :
1359 11038 : if (PRESENT(axis)) then ! Return symmetry axis.
1360 0 : axis=(/0,0,0/)
1361 0 : if (iaxis/=0) then
1362 0 : if(bravais(1)/=6)then
1363 0 : axis=ortho_axes(:,iaxis)
1364 : else
1365 0 : axis=hexa_axes(:,iaxis)
1366 : end if
1367 : end if
1368 : end if
1369 :
1370 : !DEBUG
1371 : !write(std_out,'(a)')' symbrav : exit '
1372 : !ENDDEBUG
1373 :
1374 11038 : end subroutine symbrav
1375 : !!***
1376 :
1377 : !!****f* m_symfind/symspgr
1378 : !! NAME
1379 : !! symspgr
1380 : !!
1381 : !! FUNCTION
1382 : !! Find the type of each symmetry operation (calling symcharac):
1383 : !! proper symmetries 1,2,2_1,3,3_1,3_2,4,4_1,4_2,4_3,6,6_1,...6_5
1384 : !! improper symmetries -1,m,a,b,c,d,n,g,-3,-4,-6 ,
1385 : !! Then, build an array with the number of such operations.
1386 : !! Then, call symlist to identify the space group.
1387 : !!
1388 : !! INPUTS
1389 : !! bravais(11): bravais(1)=iholohedry
1390 : !! bravais(2)=center
1391 : !! bravais(3:11)=coordinates of rprimd in the axes
1392 : !! of the conventional bravais lattice (*2 if center/=0)
1393 : !! nsym=actual number of symmetries
1394 : !! symrel(3,3,nsym)= nsym symmetry operations in real space in terms
1395 : !! of primitive translations
1396 : !! tnons(3,nsym)=nonsymmorphic translations for each symmetry (would
1397 : !! be 0 0 0 each for a symmorphic space group)
1398 : !!
1399 : !! OUTPUT
1400 : !! labels(maxsym=192)= labels of the symmetry operations
1401 : !! spgroup=symmetry space group number
1402 : !!
1403 : !! NOTES
1404 : !! It is assumed that the symmetry operations will be entered in the
1405 : !! symrel tnons arrays, for the PRIMITIVE cell. The matrix of transformation
1406 : !! from the primitive cell to the conventional cell is described
1407 : !! in the array "bravais" (see symlatt.F90).
1408 : !! The present routine first make the transformation from the
1409 : !! primitive coordinates to the conventional ones, then eventually
1410 : !! generate additional symmetries, taking into account the
1411 : !! centering translations.
1412 : !! Then, the order and determinant of each symmetry operation
1413 : !! is determined.
1414 : !!
1415 : !! For proper symmetries (rotations), the
1416 : !! associated translation is also determined.
1417 : !! However, left or right handed screw rotations are
1418 : !! not (presently) distinguished, and will be attributed equally
1419 : !! to left or right.
1420 : !!
1421 : !! For the detailed description of the labelling of the axes, see symaxes.f and symplanes.f
1422 : !!
1423 : !! SOURCE
1424 :
1425 11645 : subroutine symspgr(bravais,labels,nsym,spgroup,symrel,tnons,tolsym)
1426 :
1427 : use m_numeric_tools, only : OPERATOR(.x.)
1428 :
1429 : !Arguments ------------------------------------
1430 : !scalars
1431 : integer,intent(in) :: nsym
1432 : integer,intent(out) :: spgroup
1433 : real(dp),intent(in) :: tolsym
1434 : !arrays
1435 : integer,intent(in) :: bravais(11),symrel(3,3,nsym)
1436 : real(dp),intent(in) :: tnons(3,nsym)
1437 : character(len=128),intent(out) :: labels(192) ! 192 = maxsym
1438 :
1439 : !Local variables-------------------------------
1440 : !scalars
1441 : ! logical,parameter :: verbose=.FALSE.
1442 : integer :: additional_info,brvltt,center,direction=0,found,iholohedry,ii
1443 : integer :: ishift,isym,jj,nshift,nsymconv,spgaxor,spgorig,sporder
1444 : character(len=1) :: brvsb
1445 : character(len=15) :: intsb,ptintsb,ptschsb,schsb
1446 : character(len=35) :: intsbl
1447 : character(len=500) :: msg
1448 : !arrays
1449 : integer :: ivec1(3), ivec2(3)
1450 : integer :: n_axes(31),n_axest(31),prime(5),test_direction(3),symrel_uni(3,3)
1451 : integer :: uniaxis(3),uniaxis_try(3)
1452 11645 : integer,allocatable :: determinant(:),symrelconv(:,:,:),t_axes(:)
1453 : real(dp) :: axes(3,3),rprimdconv(3,3),vect(3,3)
1454 11645 : real(dp),allocatable :: shift(:,:),tnonsconv(:,:)
1455 : !**************************************************************************
1456 :
1457 : DBG_ENTER("COLL")
1458 :
1459 : !Initialize brvltt, from bravais(2) and bravais(1)
1460 11645 : center=bravais(2)
1461 11645 : iholohedry=bravais(1)
1462 11645 : brvltt=1
1463 11645 : if(center==-1)brvltt=2 ! Inner centering
1464 11645 : if(center==-3)brvltt=3 ! Face centering
1465 11645 : if(center==1)brvltt=5 ! A-Face centering
1466 11645 : if(center==2)brvltt=6 ! B-Face centering
1467 11645 : if(center==3)brvltt=4 ! C-Face centering
1468 11645 : if(iholohedry==5)brvltt=7 ! Rhombohedral
1469 :
1470 : !Produce the symmetry operations, in the axis of the conventional cell
1471 11645 : nsymconv=nsym
1472 11645 : if(center/=0)nsymconv=2*nsymconv
1473 11645 : if(center==-3)nsymconv=4*nsym
1474 34935 : ABI_MALLOC(symrelconv,(3,3,nsymconv))
1475 34935 : ABI_MALLOC(tnonsconv,(3,nsymconv))
1476 :
1477 : !Produce symrel and tnons in conventional axes,
1478 : !name them symrelconv and tnonsconv
1479 46580 : rprimdconv(:,1)=bravais(3:5)
1480 46580 : rprimdconv(:,2)=bravais(6:8)
1481 46580 : rprimdconv(:,3)=bravais(9:11)
1482 :
1483 79241 : if(center/=0)rprimdconv(:,:)=rprimdconv(:,:)*half
1484 :
1485 11645 : axes(:,:)=zero
1486 11645 : axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one
1487 3748326 : symrelconv(:,:,1:nsym)=symrel(:,:,1:nsym)
1488 : !Note that the number of symmetry operations is still nsym
1489 11645 : call symrelrot(nsym,rprimdconv,axes,symrelconv,tolsym)
1490 :
1491 11645 : call xred2xcart(nsym,rprimdconv,tnonsconv,tnons)
1492 : !Gives the associated translation, with components in the interval ]-0.5,0.5] .
1493 1161393 : tnonsconv(:,1:nsym)=tnonsconv(:,1:nsym)-nint(tnonsconv(:,1:nsym)-tol6)
1494 :
1495 : !If the Bravais lattice is centered, duplicate or quadruplicate
1496 : !the number of symmetry operations, using the Bravais lattice shifts
1497 11645 : nshift=1
1498 11645 : if(center/=0)nshift=2
1499 11645 : if(center==-3)nshift=4
1500 23290 : ABI_MALLOC(shift,(3,nshift))
1501 46580 : shift(:,1)=zero
1502 11645 : if(center/=0 .and. center/=-3)then
1503 5164 : shift(:,2)=half
1504 1291 : if(center==1)shift(1,2)=zero
1505 1291 : if(center==2)shift(2,2)=zero
1506 1291 : if(center==3)shift(3,2)=zero
1507 10354 : else if(center==-3)then
1508 17368 : shift(:,2)=half ; shift(1,2)=zero
1509 17368 : shift(:,3)=half ; shift(2,3)=zero
1510 17368 : shift(:,4)=half ; shift(3,4)=zero
1511 : end if ! center/=0 or -3
1512 11645 : if(nshift/=1)then
1513 19950 : do ishift=2,nshift
1514 14757877 : symrelconv(:,:,(ishift-1)*nsym+1:ishift*nsym)=symrelconv(:,:,1:nsym)
1515 587010 : do isym=1,nsym
1516 2282557 : tnonsconv(:,(ishift-1)*nsym+isym)=tnonsconv(:,isym)+shift(:,ishift)
1517 : end do
1518 : end do ! ishift
1519 : end if ! nshift/=1
1520 :
1521 : !At this stage, all the symmetry operations are available,
1522 : !expressed in the conventional axis, and also include
1523 : !the Bravais lattive translations, and associated operations...
1524 :
1525 11645 : n_axes(:)=0
1526 :
1527 34935 : ABI_MALLOC(determinant,(nsymconv))
1528 :
1529 : !Get the determinant
1530 11645 : call symdet(determinant,nsymconv,symrelconv)
1531 :
1532 : !Get the order of each the symmetry operation, as well as the maximal order
1533 : !Also, examine whether each symmetry operation is the inversion, or a root of the inversion (like -3)
1534 : !Decide which kind of point symmetry operation it is
1535 : !Finally assign tnonsconv order and decide the space symmetry operation
1536 :
1537 23290 : ABI_MALLOC(t_axes,(nsymconv))
1538 :
1539 866142 : do isym=1,nsymconv
1540 :
1541 : ! Note : nsymconv might be bigger than 192, but only for non-primitive cells, in which case labels will not be echoed anywhere.
1542 : ! 192 is the fixed dimension of labels, so this avoids possible memory problems.
1543 : call symcharac(center, determinant(isym), iholohedry, isym, labels(mod(isym-1,192)+1), &
1544 854497 : symrelconv(:,:,isym), tnonsconv(:,isym), t_axes(isym))
1545 854497 : if (t_axes(isym) == -1) then
1546 0 : write(msg, '(a,a,i3,a,3(a,3i4,a),a,3es22.12,a,a,3es22.12)' )ch10,&
1547 0 : ' symspgr: problem with isym=',isym,ch10,&
1548 0 : ' symrelconv(:,1,isym)=',symrelconv(:,1,isym),ch10,&
1549 0 : ' symrelconv(:,2,isym)=',symrelconv(:,2,isym),ch10,&
1550 0 : ' symrelconv(:,3,isym)=',symrelconv(:,3,isym),ch10,&
1551 0 : ' tnonsconv(:,isym)=',tnonsconv(:,isym)
1552 0 : call wrtout(std_out,msg)
1553 0 : write(msg, '(a,i0,2a)' )'The space symmetry operation number',isym,ch10,'is not a (translated) root of unity'
1554 0 : ABI_BUG(msg)
1555 854497 : else if (t_axes(isym) == -2) then
1556 0 : write(msg, '(a,i0,a)' )'The symmetry operation number ',isym,' is not a root of unity'
1557 0 : ABI_BUG(msg)
1558 : end if
1559 :
1560 866142 : n_axes(t_axes(isym))=n_axes(t_axes(isym))+1
1561 : end do ! isym=1,nsymconv
1562 :
1563 372640 : if (sum(n_axes)-nsymconv/=0) then
1564 : write(msg, '(7a)' )&
1565 0 : 'Not all the symmetries have been recognized. ',ch10,&
1566 0 : 'This might be due either to an error in the input file',ch10,&
1567 0 : 'or to a BUG in ABINIT',ch10,&
1568 0 : 'Please contact the ABINIT group.'
1569 0 : ABI_WARNING(msg)
1570 : end if
1571 :
1572 : !write(std_out,*)' symspgr : brvltt,nsymconv=',brvltt,nsymconv
1573 : !write(std_out,*)' n_axes(1:10)=',n_axes(1:10)
1574 : !write(std_out,*)' n_axes(11:20)=',n_axes(11:20)
1575 : !write(std_out,*)' n_axes(21:31)=',n_axes(21:31)
1576 :
1577 : !Treat cases in which the space group cannot be identified on the
1578 : !basis of n_axes one need additional information
1579 11645 : if(brvltt==1)then
1580 : ! If the bravais lattice is primitive
1581 5204 : if(nsymconv==4)then
1582 309 : 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/)
1583 9888 : if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 27 (Pcc2) or 32 (Pba2)
1584 16 : write(std_out,*)' symspgr: 27 or 32'
1585 16 : additional_info=2
1586 : ! Select binary axis
1587 80 : do isym=1,nsymconv
1588 80 : if(t_axes(isym)==8)then
1589 : ! Find direction of binary axis
1590 16 : if(symrelconv(1,1,isym)==1)direction=1
1591 16 : if(symrelconv(2,2,isym)==1)direction=2
1592 16 : if(symrelconv(3,3,isym)==1)direction=3
1593 : end if
1594 : end do
1595 : ! Examine the projection of the translation vector of the a, b or c mirror planes onto the binary axis
1596 80 : do isym=1,nsymconv
1597 80 : if(t_axes(isym)==16)then
1598 32 : if(abs(tnonsconv(direction,isym))>tol8)additional_info=1
1599 : end if
1600 : end do
1601 : end if
1602 4895 : else if(nsymconv==8)then
1603 879 : 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/)
1604 28128 : if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 55 (Pbam) or 57 (Pbcm)
1605 26 : write(std_out,*)' symspgr: 55 or 57'
1606 26 : additional_info=1
1607 : ! Select mirror plane m
1608 234 : do isym=1,nsymconv
1609 234 : if(t_axes(isym)==15)then
1610 : ! Find direction of mirror plane
1611 26 : if(symrelconv(1,1,isym)==-1)direction=1
1612 26 : if(symrelconv(2,2,isym)==-1)direction=2
1613 26 : if(symrelconv(3,3,isym)==-1)direction=3
1614 : end if
1615 : end do
1616 : ! Examine the projection of the translation vector of the a, b, or c mirror planes
1617 : ! onto the binary axis
1618 234 : do isym=1,nsymconv
1619 234 : if(t_axes(isym)==16)then
1620 52 : if(abs(tnonsconv(direction,isym))>tol8)additional_info=2
1621 : end if
1622 : end do
1623 : end if
1624 879 : 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/)
1625 28128 : if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 56 (Pccn) or 60 (Pbcn)
1626 26 : write(std_out,*)' symspgr: 56 or 60'
1627 26 : additional_info=1
1628 : ! Select mirror plane n
1629 234 : do isym=1,nsymconv
1630 234 : if(t_axes(isym)==18)then
1631 : ! Find direction of mirror plane
1632 26 : if(symrelconv(1,1,isym)==-1)direction=1
1633 26 : if(symrelconv(2,2,isym)==-1)direction=2
1634 26 : if(symrelconv(3,3,isym)==-1)direction=3
1635 : end if
1636 : end do
1637 : ! Examine the projection of the translation vector of the a, b, or c mirror planes
1638 : ! onto the binary axis
1639 234 : do isym=1,nsymconv
1640 234 : if(t_axes(isym)==16)then
1641 52 : if(abs(tnonsconv(direction,isym))<tol8)additional_info=2
1642 : end if
1643 : end do
1644 : end if
1645 : end if
1646 6441 : else if(brvltt==2)then
1647 : ! In the few next lines, use additional_info as a flag
1648 698 : additional_info=0
1649 : ! If the bravais lattice is inner-centered
1650 698 : if(nsymconv==8)then
1651 : ! Test spgroup 23 (I222) or 24 (I2_{1}2_{1}2_{1})
1652 105 : 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/)
1653 3360 : if(sum((n_axes-n_axest)**2)==0) additional_info=1
1654 593 : else if(nsymconv==24)then
1655 : ! Test spgroup 197 (I23) or 199 (I2_{1}3)
1656 2 : 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/)
1657 64 : if(sum((n_axes-n_axest)**2)==0) additional_info=1
1658 : end if
1659 698 : if(additional_info==1)then
1660 8 : write(std_out,*)' symspgr: (23 or 24) or (197 or 199)'
1661 : ! Select the three binary axes (they might be 2 or 2_1 !)
1662 8 : test_direction(:)=0
1663 104 : do isym=1,nsymconv
1664 104 : if(t_axes(isym)==20)then
1665 : ! Find direction of axis
1666 96 : do direction=1,3
1667 96 : if(symrelconv(direction,direction,isym)==1)then
1668 24 : test_direction(direction)=1
1669 24 : if(abs(tnonsconv(direction,isym))<tol8)then
1670 0 : vect(:,direction)=tnonsconv(:,isym)
1671 : else
1672 96 : vect(:,direction)=tnonsconv(:,isym)+half
1673 : end if
1674 96 : vect(:,direction)=vect(:,direction)-nint(vect(:,direction)-tol8)
1675 24 : vect(direction,direction)=zero
1676 : end if
1677 : end do ! direction=1,3
1678 : end if ! if binary axis
1679 : end do ! isym
1680 8 : if(test_direction(1)/=1 .or. test_direction(2)/=1 .and. test_direction(3)/=1)then
1681 : write(msg, '(5a,3i4)' )&
1682 0 : 'For space groups 23, 24, 197 or 197, the three binary axes',ch10,&
1683 0 : 'are not equally partitioned along the x, y and z directions',ch10,&
1684 0 : 'test_direction(1:3)=',test_direction(:)
1685 0 : ABI_BUG(msg)
1686 : end if
1687 8 : additional_info=1
1688 8 : if(abs(vect(1,2)-vect(1,3))>tol8 .or. abs(vect(2,1)-vect(2,3))>tol8 .or. &
1689 4 : abs(vect(3,1)-vect(3,2))>tol8) additional_info=2
1690 : end if ! additional information are needed
1691 : end if ! brvltt==1
1692 :
1693 11645 : if (brvltt==0 .or. brvltt==1) then ! Primitive
1694 5204 : call symlist_prim(additional_info,nsymconv,n_axes,spgroup)
1695 6441 : else if(brvltt==2)then
1696 698 : call symlist_bcc(additional_info,nsymconv,n_axes,spgroup)
1697 5743 : else if(brvltt==3)then
1698 4342 : call symlist_fcc(nsymconv,n_axes,spgroup)
1699 : else
1700 1401 : call symlist_others(brvltt,nsymconv,n_axes,spgroup)
1701 : end if
1702 :
1703 11645 : if(spgroup==0) then
1704 : write(msg, '(5a)' )&
1705 17 : 'Could not find the space group.',ch10,&
1706 17 : 'This often happens when the user selects a restricted set of symmetries ',ch10,&
1707 34 : 'in the input file, instead of letting the code automatically find symmetries.'
1708 17 : ABI_WARNING(msg)
1709 : end if
1710 :
1711 11645 : spgorig=1 ; spgaxor=1
1712 11645 : call spgdata(brvsb,intsb,intsbl,ptintsb,ptschsb,schsb,spgaxor,spgroup,sporder,spgorig)
1713 :
1714 11645 : if(spgroup/=0)then
1715 11628 : write(msg, '(a,i4,2x,a,a,a,a,a)' ) ' symspgr: spgroup=',spgroup,trim(brvsb),trim(intsb),' (=',trim(schsb),')'
1716 11628 : call wrtout(std_out,msg)
1717 : end if
1718 :
1719 11645 : if(bravais(1)==7)then
1720 5589 : write(msg, '(a)' ) ' symspgr: optical characteristics = isotropic '
1721 5589 : call wrtout(std_out,msg)
1722 6056 : else if(bravais(1)==4 .or. bravais(1)==5 .or. bravais(1)==6)then
1723 3053 : write(msg, '(a)' ) ' symspgr: optical characteristics = uniaxial '
1724 3053 : call wrtout(std_out,msg)
1725 : ! Identify the first symmetry operation that is order 3, 4 or 6
1726 3053 : found=0
1727 16561 : do isym=1,nsym
1728 231529 : if( minval( abs( t_axes(isym)-(/10,12,14,22,23,24,25,26,27,28,29,30,31/) ))==0) then
1729 : ! Proper rotations
1730 : found=1 ; exit
1731 :
1732 58004 : else if( minval( abs( t_axes(isym)-(/1,2,3/) ))==0) then
1733 : ! Improper symmetry operations
1734 993 : found=-1 ; exit
1735 : end if
1736 : end do
1737 3053 : if(found==-1 .or. found==1)then
1738 39364 : symrel_uni=symrel(:,:,isym)
1739 14944 : if(found==-1)symrel_uni=-symrel_uni
1740 : ! Now, symrel_uni is a rotation of order 3, 4, 6, for which the axis must be identified
1741 : ! It is actually the only eigenvector with eigenvalue 1. It can be found by cross products
1742 : ! Subtract the unit matrix.
1743 12112 : do ii=1,3
1744 12112 : symrel_uni(ii,ii)=symrel_uni(ii,ii)-1
1745 : end do
1746 3641 : found=0
1747 3641 : do ii=1,3
1748 3641 : jj=ii+1 ; if(jj==4)jj=1
1749 : ! Cross product
1750 25487 : ivec1 = symrel_uni(ii,:); ivec2 = symrel_uni(jj,:)
1751 3641 : uniaxis = ivec1 .x. ivec2
1752 14564 : if(sum(uniaxis**2)/=0)then
1753 : found=1 ; exit
1754 : end if
1755 : end do
1756 3028 : if(found==1)then
1757 : ! Try to reduce the length, by an integer factor (try only primes 2, 3, 5, 7, 11)
1758 3028 : prime=(/2,3,5,7,11/)
1759 3028 : ii=1
1760 19736 : do while (ii<6)
1761 66832 : uniaxis_try=uniaxis/prime(ii)
1762 69860 : if(sum(abs(uniaxis_try*prime(ii)-uniaxis))==0)then
1763 1568 : uniaxis=uniaxis_try
1764 : else
1765 15140 : ii=ii+1
1766 : end if
1767 : end do
1768 3028 : write(msg, '(a,3i4)' ) ' Optical axis (in reduced coordinates, real space ) :',uniaxis
1769 : end if
1770 : end if
1771 : if(found==0)then
1772 25 : write(msg, '(a)' ) ' However, the axis has not been found. Sorry for this.'
1773 : end if
1774 3053 : call wrtout(std_out,msg)
1775 : end if
1776 :
1777 11645 : ABI_FREE(determinant)
1778 11645 : ABI_FREE(shift)
1779 11645 : ABI_FREE(symrelconv)
1780 11645 : ABI_FREE(tnonsconv)
1781 11645 : ABI_FREE(t_axes)
1782 :
1783 : DBG_EXIT("COLL")
1784 :
1785 11645 : end subroutine symspgr
1786 : !!***
1787 :
1788 : !!****f* m_symfind/symlatt
1789 : !! NAME
1790 : !! symlatt
1791 : !!
1792 : !! FUNCTION
1793 : !! From the unit cell vectors (rprimd) and the corresponding metric tensor,
1794 : !! find the Bravais lattice and its symmetry operations (ptsymrel).
1795 : !! 1) Find the shortest possible primitive vectors for the lattice
1796 : !! 2) Determines the holohedral group of the lattice, and the
1797 : !! axes to be used for the conventional cell
1798 : !! (this is a delicate part, in which the centering of the
1799 : !! reduced cell must be taken into account)
1800 : !! The idea is to determine the basis vectors of the conventional
1801 : !! cell from the reduced cell basis vectors.
1802 : !! 3) Generate the symmetry operations of the holohedral group
1803 : !!
1804 : !! INPUTS
1805 : !! iout=unit number of output file
1806 : !! msym=default maximal number of symmetries. WARNING : cannot be simply set to nsym, because
1807 : !! the number of symmetries found here will likely be bigger than sym !
1808 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1809 : !! tolsym=tolerance for the symmetries
1810 : !!
1811 : !! OUTPUT
1812 : !! bravais(11): bravais(1)=iholohedry
1813 : !! bravais(2)=center
1814 : !! bravais(3:11)=coordinates of rprim in the axes
1815 : !! of the conventional bravais lattice (*2 if center/=0)
1816 : !! nptsym=number of point symmetries of the Bravais lattice
1817 : !! ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
1818 : !! of the Bravais lattice in real space in terms of primitive translations.
1819 : !!
1820 : !! NOTES
1821 : !! WARNING: bravais(1) might be given a negative value in another
1822 : !! routine, if the cell is non-primitive.
1823 : !! The holohedral groups are numbered as follows
1824 : !! (see international tables for crystallography (1983), p. 13)
1825 : !! iholohedry=1 triclinic 1bar
1826 : !! iholohedry=2 monoclinic 2/m
1827 : !! iholohedry=3 orthorhombic mmm
1828 : !! iholohedry=4 tetragonal 4/mmm
1829 : !! iholohedry=5 trigonal 3bar m
1830 : !! iholohedry=6 hexagonal 6/mmm
1831 : !! iholohedry=7 cubic m3bar m
1832 : !! Centering
1833 : !! center=0 no centering
1834 : !! center=-1 body-centered
1835 : !! center=-3 face-centered
1836 : !! center=1 A-face centered
1837 : !! center=2 B-face centered
1838 : !! center=3 C-face centered
1839 : !!
1840 : !! SOURCE
1841 :
1842 39896 : subroutine symlatt(bravais,iout,msym,nptsym,ptsymrel,rprimd,tolsym)
1843 :
1844 : !Arguments ------------------------------------
1845 : !scalars
1846 : integer,intent(in) :: iout,msym
1847 : integer,intent(out) :: nptsym
1848 : real(dp),intent(in) :: tolsym
1849 : !arrays
1850 : integer,intent(out) :: bravais(11),ptsymrel(3,3,msym)
1851 : real(dp),intent(in) :: rprimd(3,3)
1852 :
1853 : !Local variables-------------------------------
1854 : !scalars
1855 : integer,parameter :: mgen=4
1856 : integer :: center,fact,found,foundc,ia,iaxis1,iaxis2
1857 : integer :: isign1,isign2,ib,icase,igen,iholohedry,ii,index,isym
1858 : integer :: itrial,jj,jsym,ngen=0,orthogonal,sign12,sign13,sign23,sumsign
1859 : real(dp) :: determinant,norm2a,norm2b,norm2c,norm2trial,reduceda,reducedb,sca
1860 : real(dp) :: scalarprod,scb,trace,trace_best,val
1861 : character(len=500) :: msg
1862 : !arrays
1863 : integer,parameter :: list_holo(7)=(/7,6,4,3,5,2,1/)
1864 : integer :: ang90(3),equal(3),gen(3,3,mgen),gen2xy(3,3),gen2y(3,3),gen2z(3,3)
1865 : integer :: gen3(3,3),gen6(3,3),icoord(3,3),identity(3,3),nvecta(3),nvectb(3)
1866 : integer :: order(mgen)
1867 : real(dp) :: axes(3,3),axesinvt(3,3),axes_best(3,3),axes_try(3,3)
1868 : real(dp) :: cell_base(3,3),coord(3,3),metmin(3,3)
1869 : real(dp) :: minim(3,3),scprods(3,3),vecta(3),vectb(3),vectc(3),vin1(3),vin2(3),vext(3)
1870 : !**************************************************************************
1871 :
1872 : !DEBUG
1873 : !write(std_out,'(a,es14.6)') ' m_symfind%symlatt : enter, tolsym= ',tolsym
1874 : !call flush(std_out)
1875 : !ENDDEBUG
1876 :
1877 39896 : identity(:,:)=0 ; identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
1878 39896 : nvecta(1)=2 ; nvectb(1)=3
1879 39896 : nvecta(2)=1 ; nvectb(2)=3
1880 39896 : nvecta(3)=1 ; nvectb(3)=2
1881 :
1882 : !--------------------------------------------------------------------------
1883 : !Reduce the input vectors to a set of minimal vectors
1884 39896 : call smallprim(metmin,minim,rprimd)
1885 :
1886 : !DEBUG
1887 : !write(std_out,*)' symlatt : minim(:,1)=',minim(:,1)
1888 : !write(std_out,*)' symlatt : minim(:,2)=',minim(:,2)
1889 : !write(std_out,*)' symlatt : minim(:,3)=',minim(:,3)
1890 : !call flush(std_out)
1891 : !ENDDEBUG
1892 :
1893 : !--------------------------------------------------------------------------
1894 : !Examine the angles and vector lengths
1895 39896 : ang90(:)=0
1896 39896 : if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1
1897 39896 : if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1
1898 39896 : if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1
1899 39896 : equal(:)=0
1900 39896 : if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1
1901 39896 : if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1
1902 39896 : if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1
1903 :
1904 : !DEBUG
1905 : !write(std_out,*)' ang90=',ang90(:)
1906 : !write(std_out,*)' equal=',equal(:)
1907 : !call flush(std_out)
1908 : !ENDDEBUG
1909 :
1910 : !-----------------------------------------------------------------------
1911 : !Identification of the centering
1912 :
1913 39896 : foundc=0
1914 : !Default values
1915 39896 : fact=1 ; center=0
1916 39896 : cell_base(:,:)=minim(:,:)
1917 :
1918 : !Examine each holohedral group
1919 : !This search is ordered : should not be happy with tetragonal,
1920 : !while there is FCC ...
1921 125500 : do index=1,6
1922 :
1923 : ! If the holohedry is already found, exit
1924 121802 : if(foundc==1)exit
1925 :
1926 : ! Initialize the target holohedry
1927 85604 : iholohedry=list_holo(index)
1928 :
1929 : ! DEBUG
1930 : ! write(std_out,*)
1931 : ! write(std_out,*)' symlatt : trial holohedry',iholohedry
1932 : ! ENDDEBUG
1933 :
1934 85604 : orthogonal=0
1935 85604 : if(iholohedry==7 .or. iholohedry==4 .or. iholohedry==3)orthogonal=1
1936 :
1937 : ! Now, will examine different working hypothesis.
1938 : ! The set of these hypothesis is thought to cover all possible cases ...
1939 :
1940 : ! Working hypothesis : the basis is orthogonal
1941 85604 : if(ang90(1)+ang90(2)+ang90(3)==3 .and. orthogonal==1)then
1942 21452 : fact=1 ; center=0
1943 21452 : cell_base(:,:)=minim(:,:)
1944 : ! Checks that the basis vectors are OK for the target holohedry
1945 21452 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1946 : end if
1947 :
1948 : ! Select one trial direction
1949 300247 : do itrial=1,3
1950 :
1951 : ! If the holohedry is already found, exit
1952 209399 : if(foundc==1)exit
1953 :
1954 174747 : ia=nvecta(itrial) ; ib=nvectb(itrial)
1955 :
1956 : ! This is in case of hexagonal holohedry
1957 174747 : if(foundc==0 .and. iholohedry==6 .and. ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then
1958 7401 : reduceda=metmin(ib,ia)/metmin(ia,ia)
1959 7401 : fact=1 ; center=0
1960 7401 : if(abs(reduceda+0.5d0)<tolsym)then
1961 8508 : cell_base(:,1)=minim(:,ia)
1962 8508 : cell_base(:,2)=minim(:,ib)
1963 8508 : cell_base(:,3)=minim(:,itrial)
1964 : !DEBUG
1965 : ! write(std_out,*)' cell_base(:,1)=',cell_base(:,1)
1966 : ! write(std_out,*)' cell_base(:,2)=',cell_base(:,2)
1967 : ! write(std_out,*)' cell_base(:,3)=',cell_base(:,3)
1968 : !ENDDEBUG
1969 : ! Checks that the basis vectors are OK for the target holohedry
1970 2127 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1971 5274 : else if(abs(reduceda-0.5d0)<tolsym)then
1972 276 : cell_base(:,1)=minim(:,ia)
1973 276 : cell_base(:,2)=minim(:,ib)-minim(:,ia)
1974 276 : cell_base(:,3)=minim(:,itrial)
1975 : ! Checks that the basis vectors are OK for the target holohedry
1976 69 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1977 : end if
1978 : end if
1979 :
1980 : ! Working hypothesis : the conventional cell is orthogonal,
1981 : ! and the two other vectors are axes of the conventional cell
1982 174747 : if(foundc==0 .and. orthogonal==1 .and. ang90(itrial)==1)then
1983 :
1984 : ! Compute the reduced coordinate of trial vector in the basis
1985 : ! of the two other vectors
1986 37462 : reduceda=metmin(itrial,ia)/metmin(ia,ia)
1987 37462 : reducedb=metmin(itrial,ib)/metmin(ib,ib)
1988 149848 : cell_base(:,ia)=minim(:,ia)
1989 149848 : cell_base(:,ib)=minim(:,ib)
1990 37462 : if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. &
1991 37462 : & ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym) )then
1992 4596 : if(abs(abs(reduceda)-0.5d0)<tolsym)center=ib
1993 4596 : if(abs(abs(reducedb)-0.5d0)<tolsym)center=ia
1994 4596 : fact=2
1995 : cell_base(:,itrial)= &
1996 18384 : & (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0
1997 4596 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1998 32866 : else if( abs(abs(reduceda)-0.5d0)<tolsym .and.&
1999 : & abs(abs(reducedb)-0.5d0)<tolsym ) then
2000 553 : fact=2 ; center=-1
2001 : cell_base(:,itrial)= &
2002 2212 : & (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0
2003 553 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2004 : end if
2005 : end if
2006 :
2007 : ! Working hypothesis : the conventional cell is orthogonal, and
2008 : ! the trial vector is one of the future axes, and the face perpendicular to it is centered
2009 : if(foundc==0 .and. iholohedry==3 .and. &
2010 174747 : & ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then
2011 460 : fact=2 ; center=itrial
2012 1840 : cell_base(:,ia)=minim(:,ia)+minim(:,ib)
2013 1840 : cell_base(:,ib)=minim(:,ia)-minim(:,ib)
2014 1840 : cell_base(:,itrial)=minim(:,itrial)
2015 : ! Checks that the basis vectors are OK for the target holohedry
2016 460 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2017 : end if
2018 :
2019 : ! DEBUG
2020 : ! write(std_out,*)' after test_b, foundc=',foundc
2021 : ! ENDDEBUG
2022 :
2023 : ! Working hypothesis : the conventional cell is orthogonal, and
2024 : ! the trial vector is one of the future axes
2025 174747 : if(foundc==0 .and. orthogonal==1)then
2026 : ! Compute the projection of the two other vectors on the trial vector
2027 105793 : reduceda=metmin(itrial,ia)/metmin(itrial,itrial)
2028 105793 : reducedb=metmin(itrial,ib)/metmin(itrial,itrial)
2029 : ! If both projections are half-integer, one might have found an axis
2030 105793 : if( abs(abs(reduceda)-0.5d0)<tolsym .and.&
2031 : & abs(abs(reducedb)-0.5d0)<tolsym ) then
2032 95468 : vecta(:)=minim(:,ia)-reduceda*minim(:,itrial)
2033 95468 : vectb(:)=minim(:,ib)-reducedb*minim(:,itrial)
2034 23867 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2035 23867 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2036 23867 : scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
2037 : ! Note the order of selection : body-centered is prefered
2038 : ! over face centered, which is correct for the tetragonal case
2039 23867 : if(abs(norm2a-norm2b)<tolsym*half*(norm2a+norm2b))then
2040 : ! The lattice is body centered
2041 17042 : fact=2 ; center=-1
2042 68168 : cell_base(:,ia)=vecta(:)+vectb(:)
2043 68168 : cell_base(:,ib)=vecta(:)-vectb(:)
2044 68168 : cell_base(:,itrial)=minim(:,itrial)
2045 17042 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2046 6825 : else if(abs(scalarprod)<tolsym*half*(norm2a+norm2b))then
2047 : ! The lattice is face centered
2048 0 : fact=2 ; center=-3
2049 0 : cell_base(:,ia)=2.0d0*vecta(:)
2050 0 : cell_base(:,ib)=2.0d0*vectb(:)
2051 0 : cell_base(:,itrial)=minim(:,itrial)
2052 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2053 : end if
2054 : end if
2055 : end if
2056 :
2057 : ! DEBUG
2058 : ! write(std_out,*)' after test_c, foundc=',foundc
2059 : ! ENDDEBUG
2060 :
2061 : ! Working hypothesis : the conventional cell is orthogonal,
2062 : ! and body centered with no basis vector being an axis,
2063 : ! in which case the basis vectors must be equal (even for orthorhombic)
2064 : if(foundc==0 .and. orthogonal==1 .and. &
2065 174747 : & equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then
2066 : ! Compute the combination of the two other vectors
2067 92488 : vecta(:)=minim(:,ia)+minim(:,ib)
2068 92488 : vectb(:)=minim(:,ia)-minim(:,ib)
2069 23122 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2070 23122 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2071 : ! Project the trial vector on the first of the two vectors
2072 : reduceda=( minim(1,itrial)*vecta(1)+ &
2073 : & minim(2,itrial)*vecta(2)+ &
2074 23122 : & minim(3,itrial)*vecta(3) )/norm2a
2075 : reducedb=( minim(1,itrial)*vectb(1)+ &
2076 : & minim(2,itrial)*vectb(2)+ &
2077 23122 : & minim(3,itrial)*vectb(3) )/norm2b
2078 23122 : if( abs(abs(reduceda)-0.5d0)<tolsym )then
2079 : ! The first vector is an axis
2080 3653 : fact=2 ; center=-1
2081 14612 : cell_base(:,ia)=vecta(:)
2082 14612 : vecta(:)=minim(:,itrial)-reduceda*vecta(:)
2083 14612 : vectb(:)=0.5d0*vectb(:)
2084 14612 : cell_base(:,ib)=vecta(:)+vectb(:)
2085 14612 : cell_base(:,itrial)=vecta(:)-vectb(:)
2086 3653 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2087 19469 : else if( abs(abs(reducedb)-0.5d0)<tolsym )then
2088 : ! The second vector is an axis
2089 84 : fact=2 ; center=-1
2090 336 : cell_base(:,ib)=vectb(:)
2091 336 : vectb(:)=minim(:,itrial)-reducedb*vectb(:)
2092 336 : vecta(:)=0.5d0*vecta(:)
2093 336 : cell_base(:,ia)=vectb(:)+vecta(:)
2094 336 : cell_base(:,itrial)=vectb(:)-vecta(:)
2095 84 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2096 : end if
2097 : end if
2098 :
2099 : ! Working hypothesis : the conventional cell is orthogonal,
2100 : ! and face centered, in the case where two minimal vectors are equal
2101 174747 : if(foundc==0 .and. orthogonal==1 .and. equal(itrial)==1 ) then
2102 : ! Compute the combination of these two vectors
2103 155408 : vecta(:)=minim(:,ia)+minim(:,ib)
2104 155408 : vectb(:)=minim(:,ia)-minim(:,ib)
2105 38852 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2106 38852 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2107 : ! Project the trial vector on the two vectors
2108 : reduceda=( minim(1,itrial)*vecta(1)+ &
2109 : & minim(2,itrial)*vecta(2)+ &
2110 38852 : & minim(3,itrial)*vecta(3) )/norm2a
2111 : reducedb=( minim(1,itrial)*vectb(1)+ &
2112 : & minim(2,itrial)*vectb(2)+ &
2113 38852 : & minim(3,itrial)*vectb(3) )/norm2b
2114 38852 : if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. &
2115 38852 : & ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym) )then
2116 978 : fact=2 ; center=-3
2117 : cell_base(:,itrial)= &
2118 3912 : & (minim(:,itrial)-reduceda*vecta(:)-reducedb*vectb(:) )*2.0d0
2119 3912 : cell_base(:,ia)=vecta(:)
2120 3912 : cell_base(:,ib)=vectb(:)
2121 978 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2122 : end if
2123 : end if
2124 :
2125 : ! Working hypothesis : the conventional cell is orthogonal,
2126 : ! face centered, but no two vectors are on the same "square"
2127 174747 : if(foundc==0 .and. orthogonal==1)then
2128 : ! Compute the combination of these two vectors
2129 414892 : vecta(:)=minim(:,ia)+minim(:,ib)
2130 414892 : vectb(:)=minim(:,ia)-minim(:,ib)
2131 103723 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2132 103723 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2133 : ! The trial vector length must be equal to one of these lengths
2134 103723 : if(abs(metmin(itrial,itrial)-norm2a)<tolsym*norm2a)then
2135 1211 : fact=2 ; center=-3
2136 4844 : cell_base(:,ia)=vecta(:)+minim(:,itrial)
2137 4844 : cell_base(:,ib)=vecta(:)-minim(:,itrial)
2138 : ! Project vectb perpendicular to cell_base(:,ia) and cell_base(:,ib)
2139 1211 : norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2
2140 1211 : norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2
2141 : reduceda=( cell_base(1,ia)*vectb(1)+ &
2142 : & cell_base(2,ia)*vectb(2)+ &
2143 1211 : & cell_base(3,ia)*vectb(3) )/norm2a
2144 : reducedb=( cell_base(1,ib)*vectb(1)+ &
2145 : & cell_base(2,ib)*vectb(2)+ &
2146 1211 : & cell_base(3,ib)*vectb(3) )/norm2b
2147 1211 : if( abs(abs(reduceda)-0.5d0)<tolsym .and. &
2148 : & abs(abs(reducedb)-0.5d0)<tolsym )then
2149 0 : cell_base(:,itrial)=vectb(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib)
2150 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2151 : end if
2152 102512 : else if(abs(metmin(itrial,itrial)-norm2b)<tolsym*norm2b)then
2153 17144 : fact=2 ; center=-3
2154 68576 : cell_base(:,ia)=vectb(:)+minim(:,itrial)
2155 68576 : cell_base(:,ib)=vectb(:)-minim(:,itrial)
2156 : ! Project vecta perpendicular to cell_base(:,ia) and cell_base(:,ib)
2157 17144 : norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2
2158 17144 : norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2
2159 : reduceda=( cell_base(1,ia)*vecta(1)+ &
2160 : & cell_base(2,ia)*vecta(2)+ &
2161 17144 : & cell_base(3,ia)*vecta(3) )/norm2a
2162 : reducedb=( cell_base(1,ib)*vecta(1)+ &
2163 : & cell_base(2,ib)*vecta(2)+ &
2164 17144 : & cell_base(3,ib)*vecta(3) )/norm2b
2165 17144 : if( abs(abs(reduceda)-0.5d0)<tolsym .and. &
2166 : & abs(abs(reducedb)-0.5d0)<tolsym )then
2167 65624 : cell_base(:,itrial)=vecta(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib)
2168 16406 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2169 : end if
2170 : end if
2171 : end if
2172 :
2173 : ! Working hypothesis : the cell is rhombohedral, and
2174 : ! the three minimal vectors have same length and same absolute scalar product
2175 : if(foundc==0 .and. iholohedry==5 .and. &
2176 174747 : & equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then
2177 347 : if( abs(abs(metmin(1,2))-abs(metmin(1,3)))<tolsym*metmin(1,1) .and. &
2178 347 : & abs(abs(metmin(1,2))-abs(metmin(2,3)))<tolsym*metmin(2,2) )then
2179 323 : fact=1 ; center=0
2180 323 : cell_base(:,:)=minim(:,:)
2181 : ! One might have to change the sign of one of the vectors
2182 323 : sign12=1 ; sign13=1 ; sign23=1
2183 323 : if(metmin(1,2)<0.0d0)sign12=-1
2184 323 : if(metmin(1,3)<0.0d0)sign13=-1
2185 323 : if(metmin(2,3)<0.0d0)sign23=-1
2186 323 : sumsign=sign12+sign13+sign23
2187 323 : if(sumsign==-1)then
2188 50 : if(sign12==1)cell_base(:,3)=-cell_base(:,3)
2189 50 : if(sign13==1)cell_base(:,2)=-cell_base(:,2)
2190 200 : if(sign23==1)cell_base(:,1)=-cell_base(:,1)
2191 273 : else if(sumsign==1)then
2192 4 : if(sign12==-1)cell_base(:,3)=-cell_base(:,3)
2193 4 : if(sign13==-1)cell_base(:,2)=-cell_base(:,2)
2194 16 : if(sign23==-1)cell_base(:,1)=-cell_base(:,1)
2195 : end if
2196 323 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2197 : end if
2198 : end if
2199 :
2200 : ! DEBUG
2201 : ! write(std_out,*)' after test_3a, foundc=',foundc
2202 : ! write(std_out,*)' after test_3a, itrial=',itrial
2203 : ! write(std_out,*)' after test_3a, equal(:)=',equal(:)
2204 : ! ENDDEBUG
2205 :
2206 : ! Working hypothesis : the cell is rhombohedral, one vector
2207 : ! is parallel to the trigonal axis
2208 174747 : if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 )then
2209 20496 : vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib)
2210 2928 : norm2trial=minim(1,itrial)**2+minim(2,itrial)**2+minim(3,itrial)**2
2211 : reduceda=( minim(1,itrial)*vecta(1)+ &
2212 : & minim(2,itrial)*vecta(2)+ &
2213 2928 : & minim(3,itrial)*vecta(3) )/norm2trial
2214 : reducedb=( minim(1,itrial)*vectb(1)+ &
2215 : & minim(2,itrial)*vectb(2)+ &
2216 2928 : & minim(3,itrial)*vectb(3) )/norm2trial
2217 : ! DEBUG
2218 : ! write(std_out,*)' reduceda,reducedb=',reduceda,reducedb
2219 : ! ENDDEBUG
2220 2928 : if(abs(abs(reduceda)-1.0d0/3.0d0)<tolsym .and. &
2221 : & abs(abs(reducedb)-1.0d0/3.0d0)<tolsym ) then
2222 : ! Possibly change of sign to make positive the scalar product with
2223 : ! the vector parallel to the trigonal axis
2224 82 : if(reduceda<zero)vecta(:)=-vecta(:)
2225 130 : if(reducedb<zero)vectb(:)=-vectb(:)
2226 : ! Projection on the orthogonal plane
2227 148 : vecta(:)=vecta(:)-abs(reduceda)*cell_base(:,itrial)
2228 148 : vectb(:)=vectb(:)-abs(reducedb)*cell_base(:,itrial)
2229 : ! These two vectors should have an angle of 120 degrees
2230 37 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2231 37 : scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
2232 : ! DEBUG
2233 : ! write(std_out,*)' norm2a,scalarprod=',norm2a,scalarprod
2234 : ! ENDDEBUG
2235 37 : if(abs(two*scalarprod+norm2a)<tolsym*norm2a)then
2236 12 : fact=1 ; center=0
2237 12 : if(scalarprod>0.0d0)vectb(:)=-vectb(:)
2238 : ! Now vecta and vectb have an angle of 120 degrees
2239 48 : cell_base(:,1)=cell_base(:,itrial)/3.0d0+vecta(:)
2240 48 : cell_base(:,2)=cell_base(:,itrial)/3.0d0+vectb(:)
2241 48 : cell_base(:,3)=cell_base(:,itrial)/3.0d0-vecta(:)-vectb(:)
2242 : ! DEBUG
2243 : ! write(std_out,*)' cell_base(:,1)=',cell_base(:,1)
2244 : ! write(std_out,*)' cell_base(:,2)=',cell_base(:,2)
2245 : ! write(std_out,*)' cell_base(:,3)=',cell_base(:,3)
2246 : ! ENDDEBUG
2247 12 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2248 : end if
2249 : end if
2250 : end if
2251 :
2252 : ! Working hypothesis : the cell is rhombohedral, one vector
2253 : ! is in the plane perpendicular to the trigonal axis
2254 174747 : if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then
2255 11664 : vecta(:)=minim(:,ia)+minim(:,ib)
2256 11664 : vectb(:)=minim(:,ia)-minim(:,ib)
2257 2916 : norm2trial=cell_base(1,itrial)**2+cell_base(2,itrial)**2+cell_base(3,itrial)**2
2258 2916 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2259 2916 : norm2b=vecta(1)**2+vecta(2)**2+vecta(3)**2
2260 : reduceda=( cell_base(1,itrial)*vecta(1)+ &
2261 : & cell_base(2,itrial)*vecta(2)+ &
2262 2916 : & cell_base(3,itrial)*vecta(3) )/norm2trial
2263 : reducedb=( cell_base(1,itrial)*vectb(1)+ &
2264 : & cell_base(2,itrial)*vectb(2)+ &
2265 2916 : & cell_base(3,itrial)*vectb(3) )/norm2trial
2266 2916 : if(abs(norm2trial-norm2a)<tolsym*norm2a .and. &
2267 : & abs(abs(2*reduceda)-norm2trial)<tolsym*norm2trial )then
2268 0 : fact=1 ; center=0
2269 0 : cell_base(:,1)=minim(:,ia)
2270 0 : cell_base(:,2)=-minim(:,ib)
2271 0 : cell_base(:,3)=-minim(:,ib)+2*reduceda*minim(:,itrial)
2272 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2273 2916 : else if (abs(norm2trial-norm2b)<tolsym*norm2b .and. &
2274 : & abs(abs(2*reducedb)-norm2trial)<tolsym*norm2trial )then
2275 0 : fact=1 ; center=0
2276 0 : cell_base(:,1)=minim(:,ia)
2277 0 : cell_base(:,2)=minim(:,ib)
2278 0 : cell_base(:,3)=minim(:,ib)+2*reducedb*minim(:,itrial)
2279 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2280 : end if
2281 : end if
2282 :
2283 : ! Working hypothesis : the cell is rhombohedral, two vectors
2284 : ! are in the plane perpendicular to the trigonal axis
2285 174747 : if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then
2286 20412 : vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib)
2287 2916 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2288 2916 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2289 2916 : scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
2290 2916 : if(abs(abs(2*scalarprod)-norm2a)<tolsym*norm2a)then
2291 : ! This is in order to have 120 angle between vecta and vectb
2292 1949 : if(scalarprod>0.0d0)vectb(:)=-vectb(:)
2293 : reduceda=( cell_base(1,itrial)*vecta(1)+ &
2294 : & cell_base(2,itrial)*vecta(2)+ &
2295 1238 : & cell_base(3,itrial)*vecta(3) )/norm2a
2296 : reducedb=( cell_base(1,itrial)*vectb(1)+ &
2297 : & cell_base(2,itrial)*vectb(2)+ &
2298 1238 : & cell_base(3,itrial)*vectb(3) )/norm2b
2299 1238 : fact=1 ; center=0
2300 4952 : cell_base(:,1)=minim(:,itrial)
2301 1238 : if(abs(reduceda-0.5d0)<tolsym .and. abs(reducedb)<tolsym )then
2302 900 : cell_base(:,2)=minim(:,itrial)-vecta(:)
2303 900 : cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:)
2304 225 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2305 1013 : else if(abs(reduceda-0.5d0)<tolsym.and. abs(reducedb+0.5d0)<tolsym )then
2306 24 : cell_base(:,2)=minim(:,itrial)-vecta(:)
2307 24 : cell_base(:,3)=minim(:,itrial)+vectb(:)
2308 6 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2309 1007 : else if(abs(reduceda)<tolsym .and. abs(reducedb+0.5d0)<tolsym )then
2310 0 : cell_base(:,2)=minim(:,itrial)+vectb(:)
2311 0 : cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:)
2312 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2313 1007 : else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb)<tolsym)then
2314 24 : cell_base(:,2)=minim(:,itrial)+vecta(:)
2315 24 : cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:)
2316 6 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2317 1001 : else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb-0.5d0)<tolsym)then
2318 4004 : cell_base(:,2)=minim(:,itrial)+vecta(:)
2319 4004 : cell_base(:,3)=minim(:,itrial)-vectb(:)
2320 1001 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2321 0 : else if(abs(reduceda)<tolsym .and. abs(reducedb-0.5d0)<tolsym )then
2322 0 : cell_base(:,2)=minim(:,itrial)-vectb(:)
2323 0 : cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:)
2324 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2325 : end if
2326 : end if
2327 : end if
2328 :
2329 : ! Working hypothesis : monoclinic holohedry, primitive. Then, two angles are 90 degrees
2330 : if(foundc==0 .and. iholohedry==2 .and. &
2331 174747 : & ang90(ia)==1 .and. ang90(ib)==1 ) then
2332 297 : fact=1 ; center=0
2333 1188 : cell_base(:,1)=minim(:,ia)
2334 1188 : cell_base(:,2)=minim(:,itrial)
2335 1188 : cell_base(:,3)=minim(:,ib)
2336 : ! Checks that the basis vectors are OK for the target holohedry
2337 297 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2338 : end if
2339 :
2340 : ! Monoclinic holohedry, one-face-centered cell
2341 : ! Working hypothesis, two vectors have equal length.
2342 1048482 : do icase=1,5
2343 1048482 : if(foundc==0 .and. iholohedry==2 .and. equal(itrial)==1 ) then
2344 14692 : vecta(:)=cell_base(:,ia)+cell_base(:,ib)
2345 14692 : vectb(:)=cell_base(:,ia)-cell_base(:,ib)
2346 3673 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2347 3673 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2348 : ! The minim(:,trial) vector belongs to the
2349 : ! plane parallel to the cell_base(:,ia),cell_base(:,ib) plane
2350 : ! In that plane, must try minim(:,itrial),
2351 : ! as well as the 4 different combinations of
2352 : ! minim(:,itrial) with the vectors in the plane
2353 6850 : if(icase==1)vectc(:)=minim(:,itrial)
2354 5644 : if(icase==2)vectc(:)=minim(:,itrial)+cell_base(:,ia)
2355 5638 : if(icase==3)vectc(:)=minim(:,itrial)+cell_base(:,ib)
2356 5626 : if(icase==4)vectc(:)=minim(:,itrial)-cell_base(:,ia)
2357 5626 : if(icase==5)vectc(:)=minim(:,itrial)-cell_base(:,ib)
2358 3673 : norm2c=vectc(1)**2+vectc(2)**2+vectc(3)**2
2359 : sca=vectc(1)*vecta(1)+&
2360 : & vectc(2)*vecta(2)+&
2361 3673 : & vectc(3)*vecta(3)
2362 : scb=vectc(1)*vectb(1)+&
2363 : & vectc(2)*vectb(2)+&
2364 3673 : & vectc(3)*vectb(3)
2365 : ! DEBUG
2366 : ! write(std_out,*)' symlatt : test iholohedry=2, sca,scb=',sca,scb
2367 : ! ENDDEBUG
2368 3673 : if(abs(sca)<tolsym*sqrt(norm2c*norm2a) .or. abs(scb)<tolsym*sqrt(norm2c*norm2b))then
2369 408 : fact=2 ; center=3
2370 : ! The itrial direction is centered
2371 1632 : cell_base(:,3)=vectc(:)
2372 408 : if(abs(sca)<tolsym*sqrt(norm2c*norm2a))then
2373 560 : cell_base(:,2)=vecta(:)
2374 560 : cell_base(:,1)=vectb(:)
2375 140 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2376 268 : else if(abs(scb)<tolsym*sqrt(norm2c*norm2b))then
2377 1072 : cell_base(:,2)=vectb(:)
2378 1072 : cell_base(:,1)=vecta(:)
2379 268 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2380 : end if
2381 : end if
2382 : end if
2383 : end do ! icase=1,5
2384 :
2385 : ! Monoclinic holohedry, one-face-centered cell, but non equivalent.
2386 : ! This case, one pair of vectors is orthogonal
2387 174747 : if(foundc==0 .and. iholohedry==2 .and. ang90(itrial)==1) then
2388 4520 : vecta(:)=minim(:,ia)
2389 4520 : vectb(:)=minim(:,ib)
2390 1130 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2391 1130 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2392 : ! Project the trial vector on the two vectors
2393 : reduceda=( minim(1,itrial)*vecta(1)+ &
2394 : & minim(2,itrial)*vecta(2)+ &
2395 1130 : & minim(3,itrial)*vecta(3) )/norm2a
2396 : reducedb=( minim(1,itrial)*vectb(1)+ &
2397 : & minim(2,itrial)*vectb(2)+ &
2398 1130 : & minim(3,itrial)*vectb(3) )/norm2b
2399 1130 : if(abs(abs(reduceda)-0.5d0)<tolsym .or. abs(abs(reducedb)-0.5d0)<tolsym) then
2400 713 : fact=2 ; center=3
2401 713 : if(abs(abs(reduceda)-0.5d0)<tolsym)then
2402 20 : cell_base(:,2)=vecta(:)
2403 20 : cell_base(:,3)=vectb(:)
2404 20 : cell_base(:,1)=2*(minim(:,itrial)-reduceda*vecta(:))
2405 5 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2406 708 : else if(abs(abs(reducedb)-0.5d0)<tolsym)then
2407 2832 : cell_base(:,2)=vectb(:)
2408 2832 : cell_base(:,3)=vecta(:)
2409 2832 : cell_base(:,1)=2*(minim(:,itrial)-reducedb*vectb(:))
2410 708 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2411 : end if
2412 : end if
2413 : end if
2414 :
2415 : ! Monoclinic holohedry, one-face-centered cell, but non equivalent.
2416 : ! This case, no pair of vectors is orthogonal, no pair of vector of equal lentgh
2417 260351 : if(foundc==0 .and. iholohedry==2)then
2418 : ! Try to find a vector that belongs to the mediator plane, or the binary vector.
2419 : ! There must be one such vector, if centered monoclinic and no pair of vectors of equal length,
2420 : ! either among the three vectors, or among one of their differences or sums.
2421 : ! And there must be, among the two other vectors, one vector whose projection
2422 : ! on this vector is half the length of this vector.
2423 28768 : vecta(:)=minim(:,ia)
2424 28768 : vectb(:)=minim(:,ib)
2425 : ! Try the different possibilities for the vector on which the projection will be half ...
2426 31919 : do ii=1,5
2427 48956 : if(ii==1)vectc(:)=minim(:,itrial)
2428 46817 : if(ii==2)vectc(:)=minim(:,itrial)+vecta(:)
2429 41126 : if(ii==3)vectc(:)=minim(:,itrial)-vecta(:)
2430 41126 : if(ii==4)vectc(:)=minim(:,itrial)+vectb(:)
2431 41015 : if(ii==5)vectc(:)=minim(:,itrial)-vectb(:)
2432 27380 : norm2trial=vectc(1)**2+vectc(2)**2+vectc(3)**2
2433 : ! Project the two vectors on the trial vector
2434 27380 : reduceda=( vectc(1)*vecta(1)+vectc(2)*vecta(2)+vectc(3)*vecta(3) )/norm2trial
2435 27380 : reducedb=( vectc(1)*vectb(1)+vectc(2)*vectb(2)+vectc(3)*vectb(3) )/norm2trial
2436 27380 : found=0
2437 27380 : if(abs(abs(reduceda)-0.5d0)<tolsym)then
2438 1967 : vin1(:)=vectc(:)
2439 7868 : vin2(:)=2.0d0*(vecta(:)-reduceda*vectc(:))
2440 1967 : vext(:)=vectb(:)
2441 : found=1
2442 25413 : else if(abs(abs(reducedb)-0.5d0)<tolsym)then
2443 686 : vin1(:)=vectc(:)
2444 2744 : vin2(:)=2.0d0*(vectb(:)-reduceda*vectc(:))
2445 686 : vext(:)=vecta(:)
2446 : found=1
2447 : end if
2448 4539 : if(found==1)exit
2449 : end do
2450 : ! Now, vin1 and vin2 are perpendicular to each other, and in the plane that contains the binary vector.
2451 : ! One of them must be the binary vector if any.
2452 : ! On the other hand, vext is out-of-plane. Might belong to the mediator plane or not.
2453 : ! If C monoclinc, then the projection of this vext on the binary vector will be either 0 or +1/2 or -1/2.
2454 : ! The binary axis must be stored in cell_base(:,2) for conventional C-cell
2455 7192 : if(found==1)then
2456 2653 : found=0
2457 :
2458 : ! Test vin1 being the binary axis
2459 2653 : norm2trial=vin1(1)**2+vin1(2)**2+vin1(3)**2
2460 2653 : reduceda=(vext(1)*vin1(1)+vext(2)*vin1(2)+vext(3)*vin1(3))/norm2trial
2461 2653 : if(abs(reduceda)<tolsym)then ! vin1 is the binary axis and vext is in the mediator plane
2462 4868 : found=1
2463 4868 : cell_base(:,1)=vin2(:)
2464 4868 : cell_base(:,2)=vin1(:)
2465 4868 : cell_base(:,3)=vext(:)
2466 1436 : else if(abs(abs(reduceda)-0.5d0)<tolsym)then ! vin1 is the binary axis and vext has +1/2 or -1/2 as projection
2467 144 : found=1
2468 144 : cell_base(:,1)=vin2(:)
2469 144 : cell_base(:,2)=vin1(:)
2470 144 : cell_base(:,3)=vext(:)-reduceda*vin1(:)+vin2(:)*half
2471 : else
2472 : ! Test vin2 being the binary axis
2473 1400 : norm2trial=vin2(1)**2+vin2(2)**2+vin2(3)**2
2474 1400 : reduceda=(vext(1)*vin2(1)+vext(2)*vin2(2)+vext(3)*vin2(3))/norm2trial
2475 1400 : if(abs(reduceda)<tolsym)then ! vin2 is the binary axis and vext is in the mediator plane
2476 328 : found=1
2477 328 : cell_base(:,1)=vin1(:)
2478 328 : cell_base(:,2)=vin2(:)
2479 328 : cell_base(:,3)=vext(:)
2480 1318 : else if(abs(abs(reduceda)-0.5d0)<tolsym)then ! vin2 is the binary axis and vext has +1/2 or -1/2 as projection
2481 0 : found=1
2482 0 : cell_base(:,1)=vin1(:)
2483 0 : cell_base(:,2)=vin2(:)
2484 0 : cell_base(:,3)=vext(:)-reduceda*vin2(:)+vin1(:)*half
2485 : end if
2486 : end if
2487 :
2488 : if(found==1)then
2489 1335 : fact=2 ; center=3
2490 1335 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2491 : end if
2492 : end if
2493 : end if
2494 :
2495 : end do ! Do-loop on three different directions
2496 : end do ! Do-loop on different target holohedries
2497 :
2498 39896 : if(foundc==0)then
2499 1601 : iholohedry=1 ; fact=1 ; center=0
2500 1601 : cell_base(:,:)=minim(:,:)
2501 : end if
2502 :
2503 : !DEBUG
2504 : !write(std_out,*)' symlatt : done with centering tests, foundc=',foundc
2505 : !write(std_out,*)' center=',center
2506 : !write(std_out,*)' iholohedry=',iholohedry
2507 : !call flush(std_out)
2508 : !ENDDEBUG
2509 :
2510 : !--------------------------------------------------------------------------
2511 : !Final check on the Bravais lattice, using the basis vectors
2512 :
2513 : !Recompute the metric tensor
2514 39896 : if(foundc==1)then
2515 153180 : do ii=1,3
2516 : metmin(:,ii)=cell_base(1,:)*cell_base(1,ii)+&
2517 : & cell_base(2,:)*cell_base(2,ii)+&
2518 497835 : & cell_base(3,:)*cell_base(3,ii)
2519 : end do
2520 : end if
2521 :
2522 : !Examine the angles and vector lengths
2523 39896 : ang90(:)=0
2524 39896 : if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1
2525 39896 : if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1
2526 39896 : if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1
2527 39896 : equal(:)=0
2528 39896 : if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1
2529 39896 : if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1
2530 39896 : if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1
2531 :
2532 : !DEBUG
2533 : !write(std_out,*)' symlatt : recompute the metric tensor '
2534 : !write(std_out,*)' ang90=',ang90
2535 : !write(std_out,*)' equal=',equal
2536 : !call flush(std_out)
2537 : !ENDDEBUG
2538 :
2539 : !The axes will be aligned with the previously determined
2540 : !basis vectors, EXCEPT for the tetragonal cell, see later
2541 39896 : axes(:,:)=cell_base(:,:)
2542 :
2543 39896 : found=0
2544 : !Check orthogonal conventional cells
2545 39896 : if(ang90(1)+ang90(2)+ang90(3)==3)then
2546 :
2547 : ! Cubic system
2548 32432 : if(equal(1)+equal(2)+equal(3)==3)then
2549 : ! However, one-face centered is not admitted
2550 24402 : if(center==0 .or. center==-1 .or. center==-3)then
2551 16502 : iholohedry=7 ; found=1
2552 7900 : if(center==0)then
2553 6709 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cP (primitive cubic)'
2554 17693 : else if(center==-1)then
2555 1191 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cI (body-centered cubic)'
2556 : else if(center==-3)then
2557 16502 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cF (face-centered cubic)'
2558 : end if
2559 : end if
2560 : end if
2561 :
2562 : ! Tetragonal system
2563 8030 : if(found==0 .and. &
2564 : & (equal(1)==1 .or. equal(2)==1 .or. equal(3)==1) )then
2565 : ! However, one-face centered or face-centered is not admitted
2566 5429 : if(center==0 .or. center==-1)then
2567 5351 : iholohedry=4 ; found=1
2568 5351 : if(equal(1)==1)then
2569 24540 : axes(:,3)=cell_base(:,1) ; axes(:,1)=cell_base(:,2) ; axes(:,2)=cell_base(:,3)
2570 2897 : else if(equal(2)==1)then
2571 1260 : axes(:,3)=cell_base(:,2) ; axes(:,2)=cell_base(:,1) ; axes(:,1)=cell_base(:,3)
2572 2771 : else if(equal(3)==1)then
2573 2771 : axes(:,:)=cell_base(:,:)
2574 : end if
2575 5351 : if(center==0)then
2576 4745 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is tP (primitive tetragonal)'
2577 : else if(center==-1)then
2578 606 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is tI (body-centered tetragonal)'
2579 : end if
2580 : end if
2581 : end if
2582 :
2583 : ! Orthorhombic system
2584 : if(found==0)then
2585 2679 : iholohedry=3 ; found=1
2586 2679 : axes(:,:)=cell_base(:,:)
2587 2679 : if(center==0)then
2588 1751 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oP (primitive orthorhombic)'
2589 928 : else if(center==-1)then
2590 302 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oI (body-centered orthorhombic)'
2591 626 : else if(center==1 .or. center==2 .or. center==3)then
2592 526 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oC (one-face-centered orthorhombic)'
2593 100 : else if(center==-3)then
2594 100 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oF (face-centered orthorhombic)'
2595 : end if
2596 : end if
2597 :
2598 : else
2599 :
2600 : ! Hexagonal system
2601 7464 : 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
2602 2193 : iholohedry=6 ; found=1
2603 2193 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is hP (primitive hexagonal)'
2604 : end if
2605 :
2606 : ! Rhombohedral system
2607 : if(found==0 .and. equal(1)+equal(2)+equal(3)==3 .and. &
2608 5271 : & abs(metmin(2,1)-metmin(3,2))<tolsym*metmin(2,2) .and. &
2609 : & abs(metmin(2,1)-metmin(3,1))<tolsym*metmin(1,1) )then
2610 1573 : iholohedry=5 ; found=1
2611 1573 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is hR (rhombohedral)'
2612 : end if
2613 :
2614 : ! Monoclinic system
2615 3698 : if(found==0 .and. ang90(1)+ang90(2)+ang90(3)==2 )then
2616 2097 : iholohedry=2 ; found=1
2617 2097 : if(center==0)then
2618 297 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is mP (primitive monoclinic)'
2619 1800 : else if(center==3)then
2620 1800 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is mC (one-face-centered monoclinic)'
2621 : end if
2622 : end if
2623 :
2624 : ! Triclinic system
2625 : if(found==0)then
2626 1601 : iholohedry=1 ; found=1
2627 1601 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is aP (primitive triclinic)'
2628 : end if
2629 :
2630 : end if
2631 :
2632 39896 : call wrtout(iout,msg)
2633 :
2634 : !DEBUG
2635 : !write(std_out,*)' symlatt : after checking conventional orthogonal cell '
2636 : !call flush(std_out)
2637 : !ENDDEBUG
2638 :
2639 : !--------------------------------------------------------------------------
2640 : !Make sure that axes form a right-handed coordinate system
2641 : !(Note : this should be done in the body of the routine,
2642 : !by making changes that leave the sign of the mixed product of the three
2643 : !vectors invariant)
2644 : determinant=axes(1,1)*axes(2,2)*axes(3,3) &
2645 : & +axes(1,2)*axes(2,3)*axes(3,1) &
2646 : & +axes(1,3)*axes(3,2)*axes(2,1) &
2647 : & -axes(1,1)*axes(3,2)*axes(2,3) &
2648 : & -axes(1,3)*axes(2,2)*axes(3,1) &
2649 39896 : & -axes(1,2)*axes(2,1)*axes(3,3)
2650 39896 : if(determinant<zero)then
2651 70356 : axes(:,:)=-axes(:,:)
2652 : end if
2653 :
2654 : !DEBUG
2655 : !write(std_out,'(a,i4)')' symlatt : before itrial do loop, iholohedry= ',iholohedry
2656 : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',&
2657 : !& rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
2658 : !call flush(std_out)
2659 : !ENDDEBUG
2660 :
2661 : !--------------------------------------------------------------------------
2662 : !Prefer symmetry axes on the same side as the primitive axes,
2663 : !when the changes are allowed
2664 63189 : do itrial=1,100
2665 :
2666 : ! DEBUG
2667 : ! write(std_out,'(a)')' '
2668 : ! write(std_out,'(a,i5)')' symlatt : itrial do loop, itrial= ',itrial
2669 : ! write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',&
2670 : ! & axes(:,1),ch10,axes(:,2),ch10,axes(:,3)
2671 : ! call flush(std_out)
2672 : ! ENDDEBUG
2673 :
2674 252756 : do ia=1,3
2675 : scprods(ia,:)=axes(1,ia)*rprimd(1,:)+&
2676 : & axes(2,ia)*rprimd(2,:)+&
2677 758268 : & axes(3,ia)*rprimd(3,:)
2678 758268 : norm2trial=sum(axes(:,ia)**2)
2679 821457 : scprods(ia,:)=scprods(ia,:)/sqrt(norm2trial)
2680 : end do
2681 252756 : do ia=1,3
2682 758268 : norm2trial=sum(rprimd(:,ia)**2)
2683 821457 : scprods(:,ia)=scprods(:,ia)/sqrt(norm2trial)
2684 : end do
2685 :
2686 : !DEBUG
2687 : ! write(std_out,'(a,3f12.6)')' diagonal scalar products ',scprods(1,1),scprods(2,2),scprods(3,3)
2688 : ! call flush(std_out)
2689 : !ENDDEBUG
2690 :
2691 : ! One should now try all the generators of the
2692 : ! proper rotations of each Bravais lattice, coupled with change of
2693 : ! signs of each vector. This is not done systematically in what follows ...
2694 : ! Here, the third axis is left unchanged
2695 63189 : if(iholohedry/=5)then
2696 58811 : if(scprods(1,1)<-tolsym .and. scprods(2,2)<-tolsym)then
2697 16408 : axes(:,1)=-axes(:,1) ; axes(:,2)=-axes(:,2)
2698 : cycle
2699 : end if
2700 : end if
2701 : ! The first (or second) axis is left unchanged
2702 60845 : if(iholohedry/=5 .and. iholohedry/=6)then
2703 54274 : if(scprods(2,2)<-tolsym .and. scprods(3,3)<-tolsym)then
2704 8757 : axes(:,2)=-axes(:,2) ; axes(:,3)=-axes(:,3)
2705 : cycle
2706 : end if
2707 53023 : if(scprods(1,1)<-tolsym .and. scprods(3,3)<-tolsym)then
2708 3318 : axes(:,1)=-axes(:,1) ; axes(:,3)=-axes(:,3)
2709 : cycle
2710 : end if
2711 : end if
2712 : ! Permutation of the three axis
2713 59120 : if(iholohedry==5 .or. iholohedry==7)then
2714 45199 : trace=scprods(1,1)+scprods(2,2)+scprods(3,3)
2715 45199 : if(trace+tolsym< scprods(1,2)+scprods(2,3)+scprods(3,1))then
2716 3591 : vecta(:)=axes(:,1) ; axes(:,1)=axes(:,3)
2717 3591 : axes(:,3)=axes(:,2); axes(:,2)=vecta(:)
2718 : cycle
2719 : end if
2720 44686 : if(trace+tolsym < scprods(1,3)+scprods(2,1)+scprods(3,2))then
2721 121933 : vecta(:)=axes(:,1) ; axes(:,1)=axes(:,2)
2722 121933 : axes(:,2)=axes(:,3); axes(:,3)=vecta(:)
2723 : cycle
2724 : end if
2725 : ! This case is observed when the three new vectors
2726 : ! are pointing opposite to the three original vectors
2727 : ! One takes their opposite, then switch two of them, then process
2728 : ! them again in the loop
2729 354471 : if(sum(scprods(:,:))<-tolsym)then
2730 5168 : axes(:,1)=-axes(:,1)
2731 5168 : vecta(:)=-axes(:,2)
2732 5168 : axes(:,2)=-axes(:,3)
2733 5168 : axes(:,3)=vecta(:)
2734 : cycle
2735 : end if
2736 : end if
2737 :
2738 : ! Actually, for iholohedry==7, can test specifically all possibilities
2739 : ! and take the best one.
2740 : ! Not activated, because changing the order of symmetries in many tests !
2741 : if(iholohedry==7 .and. .false.)then
2742 : ! if(iholohedry==7)then
2743 :
2744 : !DEBUG
2745 : !write(std_out,'(a,a)')ch10,' enter search of all possibilities for iholohedry==7 '
2746 : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',&
2747 : !& axes(:,1),ch10,axes(:,2),ch10,axes(:,3)
2748 : !call flush(std_out)
2749 : !ENDDEBUG
2750 :
2751 : do iaxis1=1,3
2752 : axes_try(:,1)=axes(:,iaxis1)
2753 : do iaxis2=1,2
2754 : if(iaxis1==1)axes_try(:,2)=axes(:,1+iaxis2)
2755 : if(iaxis1==2)axes_try(:,2)=axes(:,2*iaxis2-1)
2756 : if(iaxis1==3)axes_try(:,2)=axes(:,iaxis2)
2757 : if(iaxis2==1)axes_try(:,3)=axes(:,3)
2758 : if(iaxis2==2)axes_try(:,3)=axes(:,1)
2759 : if(iaxis1==1.and.iaxis2==2)axes_try(:,3)=axes(:,2)
2760 : if(iaxis1==3.and.iaxis2==1)axes_try(:,3)=axes(:,2)
2761 : do isign1=1,-1,-2
2762 : axes_try(:,1)=-axes_try(:,1)
2763 : do isign2=1,-1,-2
2764 : axes_try(:,2)=-axes_try(:,2)
2765 : determinant=axes_try(1,1)*axes_try(2,2)*axes_try(3,3) &
2766 : & +axes_try(1,2)*axes_try(2,3)*axes_try(3,1) &
2767 : & +axes_try(1,3)*axes_try(3,2)*axes_try(2,1) &
2768 : & -axes_try(1,1)*axes_try(3,2)*axes_try(2,3) &
2769 : & -axes_try(1,3)*axes_try(2,2)*axes_try(3,1) &
2770 : & -axes_try(1,2)*axes_try(2,1)*axes_try(3,3)
2771 : if(determinant<zero)axes_try(:,3)=-axes_try(:,3)
2772 : do ia=1,3
2773 : scprods(ia,:)=axes_try(1,ia)*rprimd(1,:)+&
2774 : & axes_try(2,ia)*rprimd(2,:)+&
2775 : & axes_try(3,ia)*rprimd(3,:)
2776 : norm2trial=sum(axes_try(:,ia)**2)
2777 : scprods(ia,:)=scprods(ia,:)/sqrt(norm2trial)
2778 : end do
2779 : do ia=1,3
2780 : norm2trial=sum(rprimd(:,ia)**2)
2781 : scprods(:,ia)=scprods(:,ia)/sqrt(norm2trial)
2782 : end do
2783 : trace=scprods(1,1)+scprods(2,2)+scprods(3,3)
2784 : if(iaxis1==1.and.iaxis2==1.and.isign1==1.and.isign2==1)then
2785 : trace_best=trace
2786 : axes_best=axes_try
2787 : else if (trace>trace_best+tolsym)then
2788 : trace_best=trace
2789 : axes_best=axes_try
2790 : endif
2791 : enddo ! isign2
2792 : enddo ! isign1
2793 : enddo ! iaxes2
2794 : enddo ! iaxes1
2795 : axes=axes_best
2796 : endif ! iholohedry=7
2797 0 : exit
2798 : end do
2799 :
2800 : !--------------------------------------------------------------------------
2801 :
2802 : !DEBUG
2803 : !write(std_out,'(a,a)')ch10,' after order/sign optimization do-loop '
2804 : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',&
2805 : !& rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
2806 : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',&
2807 : !& axes(:,1),ch10,axes(:,2),ch10,axes(:,3)
2808 : !call flush(std_out)
2809 : !ENDDEBUG
2810 :
2811 : !Compute the coordinates of rprimd in the system defined by axes(:,:)
2812 39896 : call matr3inv(axes,axesinvt)
2813 159584 : do ii=1,3
2814 : coord(:,ii)=rprimd(1,ii)*axesinvt(1,:)+ &
2815 : & rprimd(2,ii)*axesinvt(2,:)+ &
2816 518648 : & rprimd(3,ii)*axesinvt(3,:)
2817 : end do
2818 :
2819 : !Check that the coordinates are integers, or half-integer in
2820 : !the case there is a centering, and generate integer coordinates
2821 159584 : do ii=1,3
2822 518648 : do jj=1,3
2823 359064 : val=coord(ii,jj)*fact
2824 359064 : if(abs(val-nint(val))>fact*two*tolsym)then
2825 : write(msg,'(4a,a,3es18.10,a,a,3es18.10,a,a,3es18.10,a,a,i4)')&
2826 0 : 'One of the coordinates of rprimd in axes is non-integer,',ch10,&
2827 0 : 'or non-half-integer (if centering), within 2*tolsym.',ch10,&
2828 0 : 'coord=',coord(:,1),ch10,&
2829 0 : ' ',coord(:,2),ch10,&
2830 0 : ' ',coord(:,3),ch10,&
2831 0 : 'fact=',fact
2832 0 : ABI_BUG(msg)
2833 : end if
2834 478752 : icoord(ii,jj)=nint(val)
2835 : end do
2836 : end do
2837 :
2838 : !Store the bravais lattice characteristics
2839 39896 : bravais(1)=iholohedry
2840 39896 : bravais(2)=center
2841 159584 : bravais(3:5)=icoord(1:3,1)
2842 159584 : bravais(6:8)=icoord(1:3,2)
2843 159584 : bravais(9:11)=icoord(1:3,3)
2844 :
2845 : !--------------------------------------------------------------------------
2846 : !Initialize the set of symmetries
2847 : !Bravais lattices are always invariant under identity and inversion
2848 :
2849 : !Identity and inversion
2850 997400 : ptsymrel(:,:,1)=identity(:,:) ; ptsymrel(:,:,2)=-identity(:,:)
2851 39896 : nptsym=2
2852 :
2853 : !Keep this for IFCv70 compiler
2854 : if(nptsym/=2)then
2855 : write(msg,'(a,a,a,a)')ch10,&
2856 : ' symlatt : BUG -',ch10,&
2857 : ' Crazy error, compiler bug '
2858 : call wrtout(std_out,msg)
2859 : end if
2860 :
2861 : !--------------------------------------------------------------------------
2862 : !Initialize some generators
2863 : !gen6 is defined in a coordinated system with gamma=120 degrees
2864 39896 : gen6(:,:)=0 ; gen6(3,3)=1 ; gen6(1,1)=1 ; gen6(1,2)=-1 ; gen6(2,1)=1
2865 39896 : gen3(:,:)=0 ; gen3(1,2)=1 ; gen3(2,3)=1 ; gen3(3,1)=1
2866 39896 : gen2xy(:,:)=0 ; gen2xy(2,1)=1 ; gen2xy(1,2)=1; gen2xy(3,3)=1
2867 39896 : gen2y(:,:)=0 ; gen2y(1,1)=-1; gen2y(2,2)=1 ; gen2y(3,3)=-1
2868 39896 : gen2z(:,:)=0 ; gen2z(1,1)=-1; gen2z(2,2)=-1; gen2z(3,3)=1
2869 :
2870 : !--------------------------------------------------------------------------
2871 :
2872 : !Define the generators for each holohedry (inversion is already included)
2873 39896 : if(iholohedry==6)then
2874 2193 : ngen=2
2875 28509 : gen(:,:,1)=gen2xy(:,:) ; order(1)=2
2876 28509 : gen(:,:,2)=gen6(:,:) ; order(2)=6
2877 37703 : else if(iholohedry==5)then
2878 1573 : ngen=2
2879 20449 : gen(:,:,1)=gen2xy(:,:) ; order(1)=2
2880 20449 : gen(:,:,2)=gen3(:,:) ; order(2)=3
2881 : else
2882 469690 : gen(:,:,1)=gen2y(:,:) ; order(1)=2
2883 469690 : gen(:,:,2)=gen2z(:,:) ; order(2)=2
2884 469690 : gen(:,:,3)=gen2xy(:,:) ; order(3)=2
2885 469690 : gen(:,:,4)=gen3(:,:) ; order(4)=3
2886 36130 : if(iholohedry<=4)ngen=iholohedry-1
2887 36130 : if(iholohedry==7)ngen=4
2888 : end if
2889 :
2890 : !Build the point symmetry operations from generators, in the reduced system
2891 : !of coordinates defined by axes(:,:)
2892 39896 : if(ngen/=0)then
2893 166943 : do igen=1,ngen
2894 1410298 : do isym=1+nptsym,order(igen)*nptsym
2895 1281650 : jsym=isym-nptsym
2896 5255248 : do ii=1,3
2897 : ptsymrel(:,ii,isym)=gen(:,1,igen)*ptsymrel(1,ii,jsym)+ &
2898 : & gen(:,2,igen)*ptsymrel(2,ii,jsym)+ &
2899 16661450 : & gen(:,3,igen)*ptsymrel(3,ii,jsym)
2900 : end do
2901 : end do
2902 166943 : nptsym=order(igen)*nptsym
2903 :
2904 : end do
2905 : end if
2906 :
2907 : !--------------------------------------------------------------------------
2908 :
2909 : !Transform symmetry matrices in the system defined by rprimd
2910 39896 : call symrelrot(nptsym,axes,rprimd,ptsymrel,tolsym)
2911 :
2912 : !DEBUG
2913 : !write(std_out,'(a)') ' symlatt : exit '
2914 : !call flush(std_out)
2915 : !stop
2916 : !ENDDEBUG
2917 :
2918 39896 : end subroutine symlatt
2919 : !!***
2920 :
2921 199762 : end module m_symfind
2922 : !!***
|
