Line data Source code
1 : !!****m* ABINIT/m_symtk
2 : !! NAME
3 : !! m_symtk
4 : !!
5 : !! FUNCTION
6 : !! Low-level tools related to symmetries
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (RC, XG, GMR, MG, JWZ)
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_symtk
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 : use m_linalg_interfaces
28 :
29 : use m_fstrings, only : sjoin, ltoa
30 : use m_matrix, only : mati3inv, mati3det, matr3inv
31 : use m_numeric_tools, only : isinteger, wrap2_pmhalf
32 : use m_hide_lapack, only : matrginv
33 :
34 : implicit none
35 :
36 : private
37 : !!***
38 :
39 : public :: symdet ! Compute determinant of each input symmetry matrix sym(3,3,i)
40 : public :: chkgrp ! Checks that a set of input symmetries constitutes a group.
41 : public :: sg_multable ! Checks that a set of input symmetries constitutes a group.
42 : ! TODO: This improved version should replace chkgrp.
43 : public :: chkorthsy ! Check the orthogonality of the symmetry operations
44 : public :: chkprimit ! Check whether the cell is primitive or not.
45 : public :: symrelrot ! Transform symmetry matrices to new coordinate system.
46 : public :: littlegroup_q ! Determines the symmetry operations by which reciprocal vector q is preserved.
47 : public :: matpointsym ! Symmetrizes a 3x3 input matrix using the point symmetry of the input atom
48 : public :: holocell ! Examine whether the trial conventional cell described by cell_base
49 : ! is coherent with the required holohedral group.
50 : public :: symmetrize_rprimd ! Generate new rprimd on the basis of the expected characteristics of the conventional cell
51 : public :: symmetrize_tnons ! Enforce accurate tnons for glide and screw symmetries
52 : public :: symmetrize_xred ! Symmetrize atomic coordinates using input symmetry matrices symrel
53 : public :: symchk ! Symmetry checker for atomic coordinates.
54 : public :: symatm ! Build indsym table describing the action of the symmetry operations on the atomic positions.
55 : public :: symcharac ! Get the type of axis for the symmetry.
56 : public :: smallprim ! Find the smallest possible primitive vectors for an input lattice
57 : public :: print_symmetries ! Helper function to print symmetries in a nice format.
58 : public :: rot2str ! Return string with info on rotation.
59 : public :: sym_order ! Return the order n of a (possibly non-symmorphic) operation (rot^n = identity, S^n(r) = r + T).
60 : !!***
61 :
62 : contains
63 : !!***
64 :
65 : !!****f* m_symtk/symdet
66 : !! NAME
67 : !! symdet
68 : !!
69 : !! FUNCTION
70 : !! Compute determinant of each input symmetry matrix sym(3,3,i)
71 : !! and check that the determinant is always +/- 1. Integer arithmetic.
72 : !!
73 : !! INPUTS
74 : !! nsym=number of symmetry operations
75 : !! sym(3,3,nsym)=integer symmetry array
76 : !!
77 : !! OUTPUT
78 : !! determinant(nsym)=determinant of each symmetry operation
79 : !!
80 : !! SOURCE
81 :
82 33011 : subroutine symdet(determinant, nsym, sym)
83 :
84 : !Arguments ------------------------------------
85 : !scalars
86 : integer,intent(in) :: nsym
87 : !arrays
88 : integer,intent(in) :: sym(3,3,nsym)
89 : integer,intent(out) :: determinant(nsym)
90 :
91 : !Local variables-------------------------------
92 : !scalars
93 : integer :: det,isym
94 : character(len=500) :: msg
95 : ! *************************************************************************
96 :
97 1308559 : do isym=1,nsym
98 1275548 : call mati3det(sym(:,:,isym),det)
99 1275548 : determinant(isym)=det
100 1308559 : if (abs(det)/=1) then
101 : write(msg,'(2(a,i0), 5a)')&
102 0 : 'Abs(determinant) for symmetry number ',isym,' is ',det,' .',ch10,&
103 0 : 'For a legitimate symmetry, abs(determinant) must be 1.',ch10,&
104 0 : 'Action: check your symmetry operations (symrel) in input file.'
105 0 : ABI_ERROR(msg)
106 : end if
107 : end do
108 :
109 33011 : end subroutine symdet
110 : !!***
111 :
112 : !!****f* m_symtk/chkgrp
113 : !! NAME
114 : !! chkgrp
115 : !!
116 : !! FUNCTION
117 : !! Checks that a set of input symmetries constitutes a group.
118 : !!
119 : !! INPUTS
120 : !! nsym = number of symmetry operations
121 : !! symafm = (anti)ferromagnetic part of symmetry operations
122 : !! symrel = 3D matrix containg symmetry operations
123 : !!
124 : !! OUTPUT
125 : !! ierr=Status error.
126 : !!
127 : !! TODO
128 : !! SHOULD ALSO CHECK THE tnons !
129 : !!
130 : !! SOURCE
131 :
132 0 : subroutine chkgrp(nsym, symafm, symrel, ierr)
133 :
134 : !Arguments ------------------------------------
135 : !scalars
136 : integer,intent(in) :: nsym
137 : integer,intent(out) :: ierr
138 : !arrays
139 : integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
140 :
141 : !Local variables-------------------------------
142 : !scalars
143 : integer :: isym, jsym, ksym, symafmchk, testeq, print_warning
144 : logical :: found_inv
145 : character(len=500) :: msg
146 : !arrays
147 : integer :: chk(3,3)
148 : ! *************************************************************************
149 :
150 : !write(std_out,*)' chkgrp : enter'
151 : !write(std_out,*)' isym symrel symafm '
152 : !do isym=1,nsym
153 : !write(std_out,'(i3,a,9i3,a,i3)' )isym,' ',symrel(:,:,isym),' ',symafm(isym)
154 : !end do
155 :
156 0 : ierr = 0
157 0 : print_warning = 1
158 :
159 : ! 1) Identity must be the first symmetry.
160 0 : if (any(symrel(:,:,1) /= identity_3d .or. symafm(1)/=1 )) then
161 0 : ABI_WARNING("First operation must be the identity operator")
162 0 : ierr = ierr + 1
163 : end if
164 :
165 : ! 2) The inverse of each element must belong to the group.
166 0 : do isym=1,nsym
167 0 : call mati3inv(symrel(:,:,isym), chk)
168 0 : chk = transpose(chk)
169 0 : found_inv = .FALSE.
170 0 : do jsym=1,nsym
171 0 : if (all(symrel(:,:,jsym) == chk) .and. (symafm(jsym) * symafm(isym) == 1)) then
172 : found_inv = .TRUE.; EXIT
173 : end if
174 : end do
175 :
176 0 : if (.not. found_inv) then
177 : write(msg,'(a,i0,2a)')&
178 0 : "Cannot find the inverse of symmetry operation ",isym,ch10,"Input symmetries do not form a group!"
179 0 : ABI_WARNING(msg)
180 0 : ierr = ierr + 1
181 : end if
182 : end do
183 :
184 : ! Check closure under composition.
185 0 : do isym=1,nsym
186 0 : do jsym=1,nsym
187 :
188 : ! Compute the product of the two symmetries
189 0 : chk = MATMUL(symrel(:,:,jsym), symrel(:,:,isym))
190 0 : symafmchk = symafm(jsym) * symafm(isym)
191 :
192 : ! Check that product array is one of the original symmetries.
193 0 : do ksym=1,nsym
194 0 : testeq = 1
195 0 : if ( ANY(chk/=symrel(:,:,ksym) )) testeq = 0
196 : #if 0
197 : ! FIXME this check make v4/t26 and v4/t27 fails.
198 : ! The rotational part is in the group but with different magnetic part!
199 : if (symafmchk /= symafm(ksym)) testeq=0
200 : #endif
201 0 : if (testeq==1) exit ! The test is positive
202 : end do
203 :
204 0 : if (testeq == 0 .and. print_warning == 1) then
205 : ! The test is negative
206 : write(msg, '(a,2i3,a,9a)' )&
207 0 : 'Product of symmetries',isym,jsym,' is not in group.',ch10,&
208 0 : 'This indicates that the input symmetry elements',ch10,&
209 0 : 'do not possess closure under group composition.',ch10,&
210 0 : 'ABINIT might stop with an ERROR after trying to correct and making a few more checks.',ch10,&
211 0 : 'Action: check symrel, symafm and possibly atomic positions, and fix them.'
212 0 : ABI_WARNING(msg)
213 0 : ierr = ierr + 1
214 0 : print_warning = 0
215 : end if
216 :
217 : end do ! jsym
218 : end do ! isym
219 :
220 0 : end subroutine chkgrp
221 : !!***
222 :
223 : !!****f* m_symtk/sg_multable
224 : !! NAME
225 : !! sg_multable
226 : !!
227 : !! FUNCTION
228 : !! Checks that a set of input symmetries constitutes a group.
229 : !! Treat reasonably well large set of symmetries, where pure translations are present.
230 : !! The translations are optional. This allows to test symrec.
231 : !!
232 : !! INPUTS
233 : !! nsym=number of symmetry operations
234 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
235 : !! symrel(3,3,nsym)=symmetry operations in real space.
236 : !! [tnons(3,nsym)]=Fractional translations.
237 : !! [tnons_tol]= tolerance on the match for tnons
238 : !!
239 : !! OUTPUT
240 : !! ierr=Status error. A non-zero value signals failure.
241 : !! [multable(4,nsym,nsym)]= Optional output.
242 : !! multable(1,sym1,sym2) gives the index of the symmetry product S1 * S2 in the symrel array. 0 if not found.
243 : !! multable(2:4,sym1,sym2)= the lattice vector that has to added to the fractional translation
244 : !! of the operation of index multable(1,sym1,sym2) to obtain the fractional translation of the product S1 * S2.
245 : !! [toinv(4,nsym)]= Optional output.
246 : !! toinv(1,sym1)=Gives the index of the inverse of the symmetry operation.
247 : !! S1 * S1^{-1} = {E, L} with E the identity and L a real-space lattice vector.
248 : !! toinv(2:4,sym1)=The lattice vector L
249 : !! Note that toinv can be easily obtained from multable but sometimes we do not need the full table.
250 : !!
251 : !! TODO
252 : !! This improved version should replace chkgrp.
253 : !!
254 : !! SOURCE
255 :
256 92940 : subroutine sg_multable(nsym, symafm, symrel, ierr, &
257 35995 : tnons, tnons_tol, multable, toinv) ! optional
258 :
259 : !Arguments ------------------------------------
260 : !scalars
261 : integer,intent(in) :: nsym
262 : integer,intent(out) :: ierr
263 : real(dp),optional,intent(in) :: tnons_tol
264 : !arrays
265 : integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
266 : integer,optional,intent(out) :: multable(4,nsym,nsym), toinv(4,nsym)
267 : real(dp),optional,intent(in) :: tnons(3,nsym)
268 :
269 : !Local variables-------------------------------
270 : !scalars
271 : integer :: echo,found,ilist_symrel,nptsymm,prd_symafm,prd_ptsymm,ptsymm1,ptsymm2,ptsymm3, sym1,sym2,sym3
272 : real(dp) :: tnons_tol_
273 : logical :: found_inv,iseq
274 : character(len=500) :: msg
275 : !arrays
276 : integer :: nlist_symrel(48),prd_symrel(3,3),ptmultable(48,48),ptsymrel(3,3,48)
277 46470 : integer,allocatable :: ptsymm(:),list_symrel(:,:)
278 : real(dp) :: prd_tnons(3)
279 46470 : real(dp),allocatable :: tnons_(:,:)
280 : ! *************************************************************************
281 :
282 : !write(std_out,*)' m_symtk%sg_multable : enter, nsym= ',nsym
283 46470 : ierr = 0
284 :
285 139410 : ABI_MALLOC(tnons_,(3,nsym))
286 :
287 4358082 : tnons_=zero
288 3749841 : if(present(tnons)) tnons_=tnons
289 :
290 46470 : tnons_tol_=tol5
291 46470 : if(present(tnons_tol)) tnons_tol_=tnons_tol
292 :
293 : ! 1) Identity must be the first symmetry. Do not check if tnons_ == 0 as cell might not be primitive.
294 604110 : if (any(symrel(:,:,1) /= identity_3d .or. symafm(1) /= 1)) then
295 0 : ABI_WARNING("First operation must be the identity operator")
296 0 : ierr = ierr + 1
297 : end if
298 :
299 : ! 2) The inverse of each element must belong to the group.
300 46470 : echo = 1
301 1124369 : do sym1=1,nsym
302 1077900 : found_inv = .FALSE.
303 30484571 : do sym2=1,nsym
304 1219382800 : prd_symrel = matmul(symrel(:,:,sym1), symrel(:,:,sym2))
305 853567960 : prd_tnons = tnons_(:,sym1) + matmul(symrel(:,:,sym1), tnons_(:,sym2))
306 30484570 : prd_symafm = symafm(sym1)*symafm(sym2)
307 66904161 : if (all(prd_symrel == identity_3d) .and. isinteger(prd_tnons, tnons_tol_) .and. prd_symafm == 1) then
308 1077899 : found_inv = .TRUE.
309 1077899 : if (present(toinv)) then
310 0 : toinv(1, sym1) = sym2; toinv(2:4, sym1) = nint(prd_tnons)
311 : end if
312 1077899 : exit
313 : end if
314 : end do
315 :
316 46469 : if (.not. found_inv) then
317 : if(echo == 1) then
318 1 : write(msg,'(a,i0,2a)')"Cannot find the inverse of symmetry operation ",sym1,ch10,"Input symmetries do not form a group "
319 1 : ABI_WARNING(msg)
320 1 : echo = 0
321 : endif
322 1 : ierr = ierr + 1
323 1 : exit
324 : end if
325 : end do
326 :
327 : ! 3)
328 : !In order to avoid potential cubic scaling with number of atoms, in exotic cases, with large prefactor,
329 : !set up lookup table for the point symmetry part of the symmetry operations.
330 : !Still cubic, but with a reduced prefactor. To fully eliminate cubic scaling, should
331 : !set up lookup table for the tnons_ as well.
332 :
333 139410 : ABI_MALLOC(list_symrel,(nsym,48))
334 139410 : ABI_MALLOC(ptsymm,(nsym))
335 :
336 46470 : nlist_symrel(:)=0
337 : ! Initialize with the first symmetry operation
338 604110 : ptsymrel(1:3,1:3,1)=symrel(:,:,1)
339 46470 : ptsymm(1)=1
340 46470 : nptsymm=1
341 46470 : list_symrel(1,1)=1
342 46470 : nlist_symrel(1)=1
343 : !If more than one symmetry operation, then loop on the other ones, find whether the ptsymm has already been found,
344 : !or create one new item in the list
345 46470 : if(nsym/=1)then
346 1070096 : do sym1=2,nsym
347 : found=0
348 22428303 : do ptsymm2=1,nptsymm
349 42076570 : if(all(symrel(:,:,list_symrel(1,ptsymm2)) == symrel(:,:,sym1)))then
350 115453 : ptsymm(sym1)=ptsymm2 ; found=1
351 115453 : nlist_symrel(ptsymm2)=nlist_symrel(ptsymm2)+1
352 115453 : list_symrel(nlist_symrel(ptsymm2),ptsymm2)=sym1
353 115453 : cycle
354 : endif
355 : enddo
356 1070096 : if(found==0)then
357 915980 : nptsymm=nptsymm+1
358 : !write(std_out,*)' current value of nptsymm, sym1=',nptsymm, sym1
359 11907740 : ptsymrel(1:3,1:3,nptsymm)=symrel(:,:,sym1)
360 915980 : ptsymm(sym1)=nptsymm
361 915980 : nlist_symrel(nptsymm)=1
362 915980 : list_symrel(1,nptsymm)=sym1
363 : endif
364 : enddo
365 : endif
366 :
367 : !Check that each point symmetry is associated to the same number of translations
368 38663 : if(nptsymm/=1)then
369 993284 : do ptsymm1=1,nptsymm
370 993284 : if(nlist_symrel(ptsymm1)/=nlist_symrel(1))then
371 : write(msg, '(9a)' )&
372 0 : 'The number of translations (and possibly symafm) associated to the same symrel',ch10,&
373 0 : 'is not the same for all point symmetries',ch10,&
374 0 : 'This indicates that the input symmetry elements',ch10,&
375 0 : 'do not possess closure under group composition.',ch10,&
376 0 : 'Action: check symrel, symafm and fix them.'
377 0 : ABI_WARNING(msg)
378 0 : echo = 0
379 0 : ierr = ierr + 1
380 0 : if (present(multable)) then
381 0 : multable(1,:,:) = 0; multable(2:4,:,:) = huge(0)
382 : end if
383 : exit
384 : endif
385 : enddo
386 : endif
387 :
388 : !write(std_out,*)' final value of nptsymm=',nptsymm
389 :
390 : ! 4) Check closure under composition and construct multiplication table of ptsymrel
391 46470 : echo = 1
392 1008916 : do ptsymm1=1,nptsymm
393 962447 : sym1=list_symrel(1,ptsymm1)
394 39126382 : do ptsymm2=1,nptsymm
395 38163936 : sym2=list_symrel(1,ptsymm2)
396 : ! Compute the product of the two symmetries.
397 1526557440 : prd_symrel = matmul(symrel(:,:,sym1), symrel(:,:,sym2))
398 : ! Check that product array is one of the original point symmetries.
399 : iseq= .false.
400 880541709 : do ptsymm3=1,nptsymm
401 2009955119 : iseq = all(prd_symrel == symrel(:,:,list_symrel(1,ptsymm3)))
402 880541709 : if (iseq) then
403 38163935 : ptmultable(ptsymm1,ptsymm2) = ptsymm3; exit
404 : endif
405 : end do
406 :
407 39126382 : if (.not. iseq .and. echo == 1) then
408 : if (echo == 1)then
409 : ! The test is negative
410 1 : prd_symafm = symafm(sym1) * symafm(sym2)
411 28 : prd_tnons = tnons_(:, sym1) + matmul(symrel(:,:,sym1), tnons_(:,sym2))
412 : write(msg, '(a,2(i0,1x),2a,3i3,f11.6,i3,a,2(3i3,f11.6,a),5a)' )&
413 1 : 'Product of symmetries:',sym1,sym2,' is not in group.',ch10,&
414 1 : prd_symrel(1,1:3),prd_tnons(1),prd_symafm,ch10,&
415 1 : prd_symrel(2,1:3),prd_tnons(2),ch10,&
416 1 : prd_symrel(3,1:3),prd_tnons(3),ch10,&
417 1 : 'This indicates that the input symmetry elements',ch10,&
418 1 : 'do not possess closure under group composition.',ch10,&
419 2 : 'Action: check symrel, symafm and fix them.'
420 1 : ABI_WARNING(msg)
421 1 : echo = 0
422 : endif
423 1 : ierr = ierr + 1
424 1 : if (present(multable)) then
425 0 : multable(1, sym1, sym2) = 0; multable(2:4, sym1, sym2) = huge(0)
426 : end if
427 : exit
428 : end if
429 :
430 : end do ! ptsymm2
431 :
432 : !write(std_out,*)' ptmultable for ptsymm1=',ptsymm1,' by batch of 16 values '
433 : !write(std_out,'(16i3)')ptmultable(ptsymm1,1:16)
434 : !write(std_out,'(16i3)')ptmultable(ptsymm1,17:32)
435 : !write(std_out,'(16i3)')ptmultable(ptsymm1,33:48)
436 :
437 1008916 : if (echo == 0) exit
438 : end do ! ptsymm1
439 :
440 : ! 5)
441 : ! Check closure under composition and construct multiplication table.
442 : ! However, does this only if the ptgroup has been successfull.
443 46470 : if(echo/=0 .and. ierr==0)then
444 1124367 : do sym1=1,nsym
445 1077898 : ptsymm1=ptsymm(sym1)
446 61867348 : do sym2=1,nsym
447 60789450 : ptsymm2=ptsymm(sym2)
448 :
449 : !The equal number of translations for each point symmetry has been checked earlier.
450 : !If the full table is not requested, it is now sufficient to check that
451 : !the product of all symmetry operations sym1 with a pure translation (ptsymm=1), or with one of the instances
452 : !for each point symmetries is indeed present in the table.
453 : !This is done to save CPU time when the number of symmetry operations is bigger than 384.
454 :
455 60789450 : if (nsym >384 .and. .not.(present(multable))) then
456 0 : if(ptsymm2/=1 .and. sym2/=list_symrel(1,ptsymm2)) cycle
457 : end if
458 :
459 : !if(ptsymm2<1 .or. ptsymm2>48)then
460 : !write(std_out,*)' sym1,sym2,ptsymm1,ptsymm2=',sym1,sym2,ptsymm1,ptsymm2
461 : !endif
462 :
463 : ! Compute the product of the two symmetries. Convention {A,a} {B,b} = {AB, a + Ab}
464 : ! prd_symrel = matmul(symrel(:,:,sym1), symrel(:,:,sym2))
465 60789450 : prd_ptsymm=ptmultable(ptsymm1,ptsymm2)
466 790262850 : prd_symrel=ptsymrel(:,:,prd_ptsymm)
467 60789450 : prd_symafm = symafm(sym1) * symafm(sym2)
468 1702104600 : prd_tnons = tnons_(:, sym1) + matmul(symrel(:,:,sym1), tnons_(:,sym2))
469 : !write(std_out,*)' prd_ptsymm,prdsymrel=',prd_ptsymm,prd_symrel
470 :
471 : ! Check that product array is one of the original symmetries.
472 : ! Only explore those symmetries that have a symrel that is the product of the two symrel of sym1 and sym2.
473 60789450 : iseq = .False.
474 90945936 : do ilist_symrel=1,nlist_symrel(prd_ptsymm)
475 90945936 : sym3=list_symrel(ilist_symrel,prd_ptsymm)
476 90945936 : iseq = isinteger(prd_tnons(1) - tnons_(1,sym3), tnons_tol_)
477 90945936 : if(iseq)then
478 71070192 : iseq = isinteger(prd_tnons(2) - tnons_(2,sym3), tnons_tol_)
479 71070192 : if(iseq)then
480 61487848 : iseq = isinteger(prd_tnons(3) - tnons_(3,sym3), tnons_tol_)
481 61487848 : if(iseq)then
482 60805986 : iseq = (prd_symafm == symafm(sym3))
483 60805986 : if(iseq)then
484 : ! The test is positive
485 60789450 : if (present(multable)) then
486 0 : multable(1,sym1,sym2) = sym3; multable(2:4,sym1,sym2) = nint(prd_tnons - tnons_(:,sym3))
487 : end if
488 : exit
489 : endif
490 : endif
491 : endif
492 : endif
493 : end do
494 61867348 : if (.not. iseq .and. echo == 1) then
495 : if (echo == 1)then
496 : ! The test is negative
497 : write(msg, '(a,2(i0,1x),2a,3i3,f11.6,i3,a,2(3i3,f11.6,a),5a)' )&
498 0 : 'Product of symmetries:',sym1,sym2,' is not in group.',ch10,&
499 0 : prd_symrel(1,1:3),prd_tnons(1),prd_symafm,ch10,&
500 0 : prd_symrel(2,1:3),prd_tnons(2),ch10,&
501 0 : prd_symrel(3,1:3),prd_tnons(3),ch10,&
502 0 : 'This indicates that the input symmetry elements',ch10,&
503 0 : 'do not possess closure under group composition.',ch10,&
504 0 : 'Action: check symrel, symafm and fix them.'
505 0 : ABI_WARNING(msg)
506 : echo = 0
507 : endif
508 0 : ierr = ierr + 1
509 0 : if (present(multable)) then
510 0 : multable(1, sym1, sym2) = 0; multable(2:4, sym1, sym2) = huge(0)
511 : end if
512 0 : exit
513 : end if
514 : end do ! sym2
515 1124367 : if (echo == 0) exit
516 : end do ! sym1
517 : else
518 1 : if (present(multable)) then
519 0 : do sym1=1,nsym
520 0 : do sym2=1,nsym
521 0 : multable(1, sym1, sym2) = 0; multable(2:4, sym1, sym2) = huge(0)
522 : enddo
523 : enddo
524 : endif
525 : endif
526 :
527 46470 : ABI_FREE(list_symrel)
528 46470 : ABI_FREE(ptsymm)
529 46470 : ABI_FREE(tnons_)
530 : !write(std_out,*)' m_symtk%sg_multable : exit '
531 :
532 82465 : end subroutine sg_multable
533 : !!***
534 :
535 : !!****f* m_symtk/chkorthsy
536 : !! NAME
537 : !! chkorthsy
538 : !!
539 : !! FUNCTION
540 : !! Check the orthogonality of the symmetry operations
541 : !! (lengths and absolute values of scalar products should be preserved)
542 : !!
543 : !! INPUTS
544 : !! gprimd(3,3)=dimensional primitive transl. for reciprocal space (bohr**-1)
545 : !! rmet=Real space metric.
546 : !! nsym=actual number of symmetries
547 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
548 : !! symrel(3,3,1:nsym)=symmetry operations in real space in terms of primitive translations
549 : !! tolsym=defines the tolerance on the orthogonality, after multiplication by 2.
550 : !!
551 : !! SIDE EFFECTS
552 : !! iexit= if 0 at input, will do the check, and stop if there is a problem, return 0 if no problem
553 : !! if 1 at input, will always output, return 0 if no problem, -1 if there is a problem,
554 : !! also, suppresses printing of problem
555 : !!
556 : !! SOURCE
557 :
558 17072 : subroutine chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel,tolsym)
559 :
560 : !Arguments ------------------------------------
561 : !scalars
562 : integer,intent(in) :: nsym
563 : integer,intent(inout) :: iexit
564 : real(dp),intent(in) :: tolsym
565 : !arrays
566 : integer,intent(in) :: symrel(3,3,nsym)
567 : real(dp),intent(in) :: gprimd(3,3),rmet(3,3),rprimd(3,3)
568 :
569 : !Local variables-------------------------------
570 : !scalars
571 : integer :: ii,isym,jj
572 : real(dp) :: residual,rmet2
573 : character(len=500) :: msg
574 : !arrays
575 : real(dp) :: prods(3,3),rmet_sym(3,3),rprimd_sym(3,3)
576 : ! *************************************************************************
577 :
578 : !write(std_out,'(a,i3)') ' chkorthsy : enter, iexit= ',iexit
579 : !write(std_out,'(a,i3)') ' nsym=',nsym
580 : !do isym=1,nsym
581 : ! write(std_out,'(9i4)')symrel(:,:,isym)
582 : !enddo
583 : !write(std_out, '(a)') ' Matrix rprimd :'
584 : !do ii=1,3
585 : ! write(std_out, '(3es16.8)')rprimd(:,ii)
586 : !enddo
587 : !write(std_out, '(a)') ' Matrix rmet :'
588 : !do ii=1,3
589 : ! write(std_out, '(3es16.8)')rmet(:,ii)
590 : !enddo
591 :
592 17072 : rmet2=zero
593 68288 : do ii=1,3
594 221936 : do jj=1,3
595 204864 : rmet2=rmet2+rmet(ii,jj)**2
596 : end do
597 : end do
598 :
599 : !Loop over all symmetry operations
600 388213 : do isym=1,nsym
601 :
602 : !write(std_out,'(a,a,i4)') ch10,' Check for isym=',isym
603 : ! Compute symmetric of primitive vectors under point symmetry operations
604 1484680 : do ii=1,3
605 : rprimd_sym(:,ii)=symrel(1,ii,isym)*rprimd(:,1)+&
606 : symrel(2,ii,isym)*rprimd(:,2)+&
607 4825210 : symrel(3,ii,isym)*rprimd(:,3)
608 : end do
609 :
610 : ! If the new lattice is the same as the original one, the lengths and angles are preserved.
611 1484680 : do ii=1,3
612 : rmet_sym(ii,:)=rprimd_sym(1,ii)*rprimd_sym(1,:)+&
613 : rprimd_sym(2,ii)*rprimd_sym(2,:)+&
614 4825210 : rprimd_sym(3,ii)*rprimd_sym(3,:)
615 : end do
616 :
617 371170 : residual=zero
618 1484680 : do ii=1,3
619 4825210 : do jj=1,3
620 4454040 : residual=residual+(rmet_sym(ii,jj)-rmet(ii,jj))**2
621 : end do
622 : end do
623 :
624 371170 : if(sqrt(residual) > four*tolsym*sqrt(rmet2))then
625 29 : if(iexit==0)then
626 0 : write(std_out, '(a)') ' Matrix rprimd :'
627 0 : do ii=1,3
628 0 : write(std_out, '(3es16.8)')rprimd(:,ii)
629 : enddo
630 0 : write(std_out, '(a)') ' Matrix rmet :'
631 0 : do ii=1,3
632 0 : write(std_out, '(3es16.8)')rmet(:,ii)
633 : enddo
634 0 : write(std_out, '(a)') ' Matrix rprimd_sym :'
635 0 : do ii=1,3
636 0 : write(std_out, '(3es16.8)')rprimd_sym(:,ii)
637 : enddo
638 0 : write(std_out, '(a)') ' Matrix rmet_sym :'
639 0 : do ii=1,3
640 0 : write(std_out, '(3es16.8)')rmet_sym(:,ii)
641 : enddo
642 0 : write(std_out, '(a)') ' Matrix rmet_sym-rmet :'
643 0 : do ii=1,3
644 0 : write(std_out, '(3es16.8)')(rmet_sym(:,ii)-rmet(:,ii))
645 : enddo
646 : write(msg, '(a,i0,5a,es12.4,a,es12.4,6a)' )&
647 0 : 'The symmetry operation number ',isym,' does not preserve',ch10,&
648 0 : 'vector lengths and angles.',ch10,&
649 0 : 'The value of the square root of residual is: ',sqrt(residual),&
650 0 : ' that is greater than threshold:', four*tolsym*sqrt(rmet2),ch10,&
651 0 : 'Action: modify rprim, acell and/or symrel so that',ch10,&
652 0 : 'vector lengths and angles are preserved.',ch10,&
653 0 : 'Beware, the tolerance on symmetry operations is very small.'
654 0 : ABI_ERROR(msg)
655 : else
656 29 : iexit=-1
657 : end if
658 : end if
659 :
660 : ! Also, the scalar product of rprimd_sym and gprimd must give integer numbers
661 1484680 : do ii=1,3
662 : prods(ii,:)=rprimd_sym(1,ii)*gprimd(1,:)+ &
663 : rprimd_sym(2,ii)*gprimd(2,:)+ &
664 4825210 : rprimd_sym(3,ii)*gprimd(3,:)
665 : end do
666 :
667 1484680 : do ii=1,3
668 4825210 : do jj=1,3
669 3340530 : residual=prods(ii,jj)-anint(prods(ii,jj))
670 4454040 : if(abs(residual)>two*tolsym)then
671 0 : if(iexit==0)then
672 : write(msg, '(a,i0,5a,es12.4,a,es12.4,4a)' )&
673 0 : 'The symmetry operation number ',isym,' generates',ch10,&
674 0 : 'a different lattice.',ch10,&
675 0 : 'The value of the residual is: ',residual, 'that is greater than the threshold:', two*tolsym, ch10,&
676 0 : 'Action: modify rprim, acell and/or symrel so that',ch10,&
677 0 : 'the lattice is preserved.'
678 0 : ABI_ERROR(msg)
679 : else
680 0 : iexit=-1
681 : end if
682 : end if
683 : end do
684 : end do
685 :
686 388213 : if(iexit==-1) exit
687 : end do ! isym
688 :
689 17072 : if(iexit==1)iexit=0
690 :
691 : !write(std_out,'(a)') ' chkorthsy : exit '
692 :
693 17072 : end subroutine chkorthsy
694 : !!***
695 :
696 : !!****f* m_symtk/chkprimit
697 : !! NAME
698 : !! chkprimit
699 : !!
700 : !! FUNCTION
701 : !! Check whether the cell is primitive or not. If chkprim/=0 and the cell is non-primitive, stops.
702 : !!
703 : !! INPUTS
704 : !! chkprim= if non-zero, check that the unit cell is primitive.
705 : !! nsym=actual number of symmetries.
706 : !! symafm(nsym)= (anti)ferromagnetic part of symmetry operations.
707 : !! symrel(3,3,nsym)= nsym symmetry operations in real space in terms of primitive translations.
708 : !!
709 : !! OUTPUT
710 : !! multi=multiplicity of the unit cell
711 : !! translation(nsym)= (optional) set to 1 if the symetry operation is a pure translation
712 : !!
713 : !! SOURCE
714 :
715 4335 : subroutine chkprimit(chkprim, multi, nsym, symafm, symrel, is_translation)
716 :
717 : !Arguments ------------------------------------
718 : !scalars
719 : integer,intent(in) :: chkprim,nsym
720 : integer,intent(out) :: multi
721 : !arrays
722 : integer,intent(in) :: symafm(nsym),symrel(3,3,nsym)
723 : integer,intent(out),optional :: is_translation(nsym)
724 :
725 : !Local variables-------------------------------
726 : !scalars
727 : integer :: isym
728 : character(len=500) :: msg
729 : !**************************************************************************
730 :
731 158445 : if(present(is_translation)) is_translation(:)=0
732 :
733 : !Loop over each symmetry operation of the Bravais lattice
734 : !Find whether it is the identity, or a pure translation, without change of sign of the spin
735 16264 : multi=0
736 480812 : do isym=1,nsym
737 464548 : if( abs(symrel(1,1,isym)-1)+&
738 : abs(symrel(2,2,isym)-1)+&
739 : abs(symrel(3,3,isym)-1)+&
740 : abs(symrel(1,2,isym))+abs(symrel(2,1,isym))+&
741 : abs(symrel(2,3,isym))+abs(symrel(3,2,isym))+&
742 : abs(symrel(3,1,isym))+abs(symrel(1,3,isym))+&
743 16264 : abs(symafm(isym)-1) == 0 )then
744 17232 : multi=multi+1
745 17232 : if(present(is_translation))then
746 4632 : is_translation(isym)=1
747 : endif
748 : end if
749 : end do
750 :
751 : !Check whether the cell is primitive
752 16264 : if(multi>1)then
753 406 : if(chkprim>0)then
754 : write(msg,'(a,a,a,i0,a,a,a,a,a,a,a,a,a)')&
755 0 : 'According to the symmetry finder, the unit cell is',ch10,&
756 0 : 'NOT primitive. The multiplicity is ',multi,' .',ch10,&
757 0 : 'The use of non-primitive unit cells is allowed',ch10,&
758 0 : 'only when the current chkprim is 0.',ch10,&
759 0 : 'Action: either change your unit cell (rprim or angdeg),',ch10,&
760 0 : 'or set chkprim to 0.'
761 0 : ABI_ERROR(msg)
762 406 : else if(chkprim==0)then
763 : write(msg,'(3a,i0,a,a,a)')&
764 289 : 'According to the symmetry finder, the unit cell is',ch10,&
765 289 : 'not primitive, with multiplicity= ',multi,'.',ch10,&
766 578 : 'This is allowed, as the current chkprim is 0.'
767 289 : ABI_COMMENT(msg)
768 : end if
769 : end if
770 :
771 16264 : end subroutine chkprimit
772 : !!***
773 :
774 : !!****f* m_symtk/symrelrot
775 : !! NAME
776 : !! symrelrot
777 : !!
778 : !! FUNCTION
779 : !! Transform the symmetry matrices symrel expressed in the coordinate system rprimd,
780 : !! to symmetry matrices symrel expressed in the new coordinate system rprimd_new
781 : !!
782 : !! INPUTS
783 : !! nsym=number of symmetries
784 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
785 : !! rprimd_new(3,3)=new dimensional primitive translations for real space (bohr)
786 : !!
787 : !! SIDE EFFECTS
788 : !! Input/Output
789 : !! ierr= (at input) if present, will deal with error code outside of the routine.
790 : !! (at output) return 0 if no problem, 1 otherwise
791 : !! symrel(3,3,nsym)=symmetry operations in real space in terms
792 : !! of primitive translations rprimd at input and rprimd_new at output
793 : !!
794 : !! SOURCE
795 :
796 56201 : subroutine symrelrot(nsym, rprimd, rprimd_new, symrel, tolsym, ierr)
797 :
798 : !Arguments ------------------------------------
799 : !scalars
800 : integer,intent(in) :: nsym
801 : integer,intent(inout),optional :: ierr
802 : real(dp),intent(in) :: tolsym
803 : !arrays
804 : integer,intent(inout) :: symrel(3,3,nsym)
805 : real(dp),intent(in) :: rprimd(3,3),rprimd_new(3,3)
806 :
807 : !Local variables-------------------------------
808 : !scalars
809 : integer :: ierr_,ii,isym,jj
810 : real(dp) :: val
811 : character(len=500) :: msg
812 : !arrays
813 112402 : integer :: symrel_tmp(3,3,nsym)
814 : real(dp) :: coord(3,3),coordinvt(3,3),matr1(3,3),matr2(3,3),rprimd_invt(3,3)
815 : !**************************************************************************
816 :
817 56201 : ierr_=0
818 :
819 : !Compute the coordinates of rprimd_new in the system defined by rprimd(:,:)
820 56201 : call matr3inv(rprimd,rprimd_invt)
821 224804 : do ii=1,3
822 : coord(:,ii)=rprimd_new(1,ii)*rprimd_invt(1,:)+ &
823 : rprimd_new(2,ii)*rprimd_invt(2,:)+ &
824 730613 : rprimd_new(3,ii)*rprimd_invt(3,:)
825 : end do
826 :
827 : !Transform symmetry matrices in the system defined by rprimd_new
828 56201 : call matr3inv(coord,coordinvt)
829 1751874 : do isym=1,nsym
830 6782692 : do ii=1,3
831 : matr1(:,ii)=symrel(:,1,isym)*coord(1,ii)+&
832 : symrel(:,2,isym)*coord(2,ii)+&
833 22043749 : symrel(:,3,isym)*coord(3,ii)
834 : end do
835 6782692 : do ii=1,3
836 : matr2(:,ii)=coordinvt(1,:)*matr1(1,ii)+&
837 : coordinvt(2,:)*matr1(2,ii)+&
838 22043749 : coordinvt(3,:)*matr1(3,ii)
839 : end do
840 :
841 : !write(std_out, '(a,10i4)')' symrelrot : isym, symrel=',isym,symrel(:,:,isym)
842 : !write(std_out, '(a,9es16.6)')' transformed to ', matr2(:,:)
843 :
844 : ! Check that the new symmetry matrices are made of integers, and store them
845 6838893 : do ii=1,3
846 22043749 : do jj=1,3
847 15261057 : val=matr2(ii,jj)
848 : ! Need to allow for ten times tolsym, in case of centered Bravais lattices (but do it for all lattices ...)
849 15261057 : if(abs(val-nint(val))>ten*tolsym)then
850 6163 : ierr_=1
851 6163 : if(.not.(present(ierr))) then
852 : write(msg,'(2a,a,i3,a,a,3es14.6,a,a,3es14.6,a,a,3es14.6)')&
853 0 : 'One of the components of symrel is non-integer within 10*tolsym,',ch10,&
854 0 : ' for isym=',isym,ch10,&
855 0 : ' symrel=',matr2(:,1),ch10,&
856 0 : ' ',matr2(:,2),ch10,&
857 0 : ' ',matr2(:,3)
858 0 : ABI_ERROR_CLASS(msg, "TolSymError")
859 : endif
860 : end if
861 20348076 : symrel_tmp(ii,jj,isym)=nint(val)
862 : end do
863 : end do
864 : end do ! isym
865 :
866 : ! Upgrade symrel only if there is no error
867 21946290 : if(ierr_==0) symrel(:,:,:)=symrel_tmp(:,:,:)
868 :
869 56201 : if(present(ierr)) ierr=ierr_
870 :
871 56201 : end subroutine symrelrot
872 : !!***
873 :
874 : !!****f* m_symtk/littlegroup_q
875 : !! NAME
876 : !! littlegroup_q
877 : !!
878 : !! FUNCTION
879 : !! Determines the symmetry operations by which the reciprocal vector qpt is preserved,
880 : !! modulo a primitive reciprocal lattice vector, and the time-reversal symmetry.
881 : !!
882 : !! INPUTS
883 : !! nsym=number of space group symmetries
884 : !! qpt(3)= vector in reciprocal space
885 : !! symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
886 : !! [prtvol]=integer flag defining the verbosity of output. =0 if no output is provided.
887 : !! use_sym= integer flag. If 1 provide output of electron-phonon "gkk" matrix elements, for further
888 : !! treatment by mrggkk utility or anaddb utility. If 0 no output is provided.
889 : !!
890 : !! OUTPUT
891 : !! symq(4,2,nsym)= three first numbers define the G vector;
892 : !! fourth number is zero if the q-vector is not preserved, 1 otherwise
893 : !! second index is one without time-reversal symmetry, two with time-reversal symmetry
894 : !! timrev=1 if the time-reversal symmetry preserves the wavevector, modulo a reciprocal lattice vector (in principle, see below).
895 : !!
896 : !! NOTES
897 : !! The condition is: $q = O S(q) - G$
898 : !! with O being either the identity or the time reversal symmetry (= inversion in reciprocal space)
899 : !! and G being a primitive vector of the reciprocal lattice.
900 : !! If the time-reversal (alone) also preserves q, modulo a lattice vector, then timrev is set to 1, otherwise 0.
901 : !!
902 : !! TODO
903 : !! timrev is put to 1 only for Gamma. Better handling should be provided in further version.
904 : !!
905 : !! SOURCE
906 :
907 314067 : subroutine littlegroup_q(nsym, qpt, symq, symrec, symafm, timrev, &
908 : prtvol, use_sym) ! optional
909 :
910 : !Arguments -------------------------------
911 : !scalars
912 : integer,intent(in) :: nsym
913 : integer,intent(in),optional :: prtvol,use_sym
914 : integer,intent(out) :: timrev
915 : !arrays
916 : integer,intent(in) :: symrec(3,3,nsym), symafm(nsym)
917 : integer,intent(out) :: symq(4,2,nsym)
918 : real(dp),intent(in) :: qpt(3)
919 :
920 : !Local variables -------------------------
921 : !scalars
922 : integer :: ii,isign,isym,itirev,my_prtvol
923 : real(dp),parameter :: tol=2.d-8
924 : !real(dp),parameter :: tol=tol4
925 : real(dp) :: reduce
926 : character(len=500) :: msg
927 : !arrays
928 : real(dp) :: difq(3),qsym(3),shift(3)
929 : ! *********************************************************************
930 :
931 314067 : my_prtvol=0; if (PRESENT(prtvol)) my_prtvol=prtvol
932 :
933 : ! Initialise the array symq
934 97805582 : symq = 0
935 :
936 : isym = symafm(1) ! just to fool abirules and use symafm for the moment
937 :
938 9176932 : do isym=1,nsym
939 : ! if (symafm(isym) /= 1) cycle ! skip afm symops
940 : ! TODO: check how much of the afm syms are coded in the rf part of the code. cf
941 : ! test v3 / 12
942 26902662 : do itirev=1,2
943 17725730 : isign=3-2*itirev ! isign is 1 without time-reversal, -1 with time-reversal
944 :
945 : ! Get the symmetric of the vector
946 70902920 : do ii=1,3
947 : qsym(ii)=qpt(1)*isign*symrec(ii,1,isym)&
948 : +qpt(2)*isign*symrec(ii,2,isym)&
949 70902920 : +qpt(3)*isign*symrec(ii,3,isym)
950 : end do
951 :
952 : ! Get the difference between the symmetric and the original vector
953 17725730 : symq(4,itirev,isym)=1
954 70902920 : do ii=1,3
955 53177190 : difq(ii)=qsym(ii)-qpt(ii)
956 : ! Project modulo 1 in the interval ]-1/2,1/2] such that difq = reduce + shift
957 53177190 : call wrap2_pmhalf(difq(ii),reduce,shift(ii))
958 70902920 : if(abs(reduce)>tol)symq(4,itirev,isym)=0
959 : end do
960 :
961 : ! SP: When prtgkk is asked (GKK matrix element will be output), one has to
962 : ! disable symmetries. There is otherwise a gauge problem with the unperturbed
963 : ! and the perturbed wavefunctions. This leads to a +- 5% increase in computational
964 : ! cost but provide the correct GKKs (i.e. the same as without the use of symmetries.)
965 :
966 17725730 : if (PRESENT(use_sym)) then
967 2376 : if (use_sym == 0) then
968 2376 : symq(4,itirev,isym)=0
969 2376 : symq(4,itirev,1)=1
970 : end if
971 : end if
972 :
973 : ! If the operation succeded, change shift from real(dp) to integer, then exit loop
974 26588595 : if(symq(4,itirev,isym)/=0)then
975 1302689 : if (my_prtvol>0) then
976 11599 : if(itirev==1)write(msg,'(a,i4,a)')' littlegroup_q: found symmetry',isym,' preserves q '
977 11599 : if(itirev==2)write(msg,'(a,i4,a)')' littlegroup_q: found symmetry ',isym,' + TimeReversal preserves q '
978 11599 : call wrtout(std_out,msg)
979 : end if
980 : ! Uses the mathematical function NINT = nearest integer
981 5210756 : do ii=1,3
982 5210756 : symq(ii,itirev,isym)=nint(shift(ii))
983 : end do
984 : end if
985 :
986 : end do !itirev
987 : end do !isym
988 :
989 : ! Test time-reversal symmetry
990 314067 : timrev=1
991 1256268 : do ii=1,3
992 : ! Unfortunately, this version does not work yet ...
993 : ! call wrap2_pmhalf(2*qpt(ii),reduce,shift(ii))
994 : ! if(abs(reduce)>tol)timrev=0
995 : ! So, this is left ...
996 1256268 : if(abs(qpt(ii))>tol)timrev=0
997 : end do
998 :
999 314067 : if(timrev==1.and.my_prtvol>0)then
1000 : write(msg, '(3a)' )&
1001 88 : ' littlegroup_q: able to use time-reversal symmetry. ',ch10,&
1002 176 : ' (except for gamma, not yet able to use time-reversal symmetry)'
1003 88 : call wrtout(std_out,msg)
1004 : end if
1005 :
1006 314067 : end subroutine littlegroup_q
1007 : !!***
1008 :
1009 : !!****f* m_symtk/matpointsym
1010 : !! NAME
1011 : !! matpointsym
1012 : !!
1013 : !! FUNCTION
1014 : !! For given order of point group, symmetrizes a 3x3 input matrix using the
1015 : !! point symmetry of the input atom
1016 : !!
1017 : !! INPUTS
1018 : !! iatom=index of atom to symmetrize around
1019 : !! natom=number of atoms in cell
1020 : !! nsym=order of group
1021 : !! rprimd(3,3)= real space primitive vectors
1022 : !! symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
1023 : !! tnons(3,nsym) = nonsymmorphic translations
1024 : !! xred(3,natom)=locations of atoms in reduced coordinates
1025 : !!
1026 : !! SIDE EFFECTS
1027 : !! mat3(3,3) = matrix to be symmetrized, in cartesian frame
1028 : !!
1029 : !! SOURCE
1030 :
1031 129 : subroutine matpointsym(iatom,mat3,natom,nsym,rprimd,symrel,tnons,xred)
1032 :
1033 : !Arguments ------------------------------------
1034 : !scalars
1035 : integer,intent(in) :: iatom,natom,nsym
1036 : !arrays
1037 : integer,intent(in) :: symrel(3,3,nsym)
1038 : real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym),xred(3,natom)
1039 : real(dp),intent(inout) :: mat3(3,3)
1040 :
1041 : !Local variables-------------------------------
1042 : !scalars
1043 : integer :: cell_index,cell_indexp,ii,isym,nsym_point
1044 : real(dp) :: xreddiff
1045 : !arrays
1046 : integer :: symrel_it(3,3)
1047 : real(dp) :: mat3_tri(3,3),mat3_tri_sym(3,3),rprimd_inv(3,3),tmp_mat(3,3)
1048 : real(dp) :: xredp(3)
1049 : !**************************************************************************
1050 :
1051 : !copy rprimd input and construct inverse
1052 129 : rprimd_inv = rprimd
1053 129 : call matrginv(rprimd_inv,3,3)
1054 :
1055 : !transform input mat3 to triclinic frame with rprimd^{-1} * mat3 * rprimd
1056 129 : call dgemm('N','N',3,3,3,one,rprimd_inv,3,mat3,3,zero,tmp_mat,3)
1057 129 : call dgemm('N','N',3,3,3,one,tmp_mat,3,rprimd,3,zero,mat3_tri,3)
1058 :
1059 : !loop over symmetry elements to obtain symmetrized input matrix
1060 129 : mat3_tri_sym = zero
1061 129 : nsym_point = 0
1062 2649 : do isym = 1, nsym
1063 :
1064 : ! skip any nonsymmorphic symmetry elements, want to consider point elements only
1065 10080 : if(dot_product(tnons(:,isym),tnons(:,isym))>tol8) cycle
1066 :
1067 : ! for current symmetry element, find transformed reduced coordinates of target atom
1068 : ! via xredp = symrel * xred
1069 17160 : call dgemv('N',3,3,one,dble(symrel(:,:,isym)),3,xred(:,iatom),1,zero,xredp,1)
1070 :
1071 :
1072 : ! shift xredp into the same unit cell as xred, for comparison
1073 : ! label cells as 0..1:0 1..2:1 2..3:2 and -1..0:-1 -2..-1:-2 and so forth
1074 5280 : do ii = 1, 3
1075 :
1076 3960 : cell_index = int(xred(ii,iatom))
1077 3960 : if(xred(ii,iatom) < zero) cell_index = cell_index - 1
1078 3960 : cell_indexp = int(xredp(ii))
1079 3960 : if(xredp(ii) < zero) cell_indexp = cell_indexp - 1
1080 :
1081 5036 : do while (cell_indexp < cell_index)
1082 1076 : xredp(ii) = xredp(ii)+one
1083 5036 : cell_indexp = cell_indexp + 1
1084 : end do
1085 5636 : do while (cell_indexp > cell_index)
1086 356 : xredp(ii) = xredp(ii)-one
1087 356 : cell_indexp = cell_indexp - 1
1088 : end do
1089 :
1090 : end do
1091 :
1092 : ! now compare xredp to xred
1093 5280 : xreddiff = dot_product(xredp-xred(:,iatom),xredp-xred(:,iatom))
1094 :
1095 1449 : if (xreddiff < tol8) then
1096 :
1097 : ! accumulate symrel^{-1}*mat3_tri*symrel into mat3_tri_sym iff xredp = xred + L,
1098 : ! where is a lattice vector, so symrel leaves the target atom invariant
1099 :
1100 : ! mati3inv gives the inverse transpose of symrel
1101 818 : call mati3inv(symrel(:,:,isym),symrel_it)
1102 10634 : call dgemm('N','N',3,3,3,one,mat3_tri,3,dble(symrel(:,:,isym)),3,zero,tmp_mat,3)
1103 10634 : call dgemm('T','N',3,3,3,one,dble(symrel_it),3,tmp_mat,3,one,mat3_tri_sym,3)
1104 818 : nsym_point = nsym_point + 1
1105 : end if
1106 :
1107 : end do
1108 :
1109 : !normalize by number of point symmetry operations
1110 1677 : mat3_tri_sym = mat3_tri_sym/dble(nsym_point)
1111 :
1112 : !transform mat3_tri_sym to cartesian frame with rprimd * mat3_tri_sym * rprimd^{-1}
1113 :
1114 129 : call dgemm('N','N',3,3,3,one,mat3_tri_sym,3,rprimd_inv,3,zero,tmp_mat,3)
1115 129 : call dgemm('N','N',3,3,3,one,rprimd,3,tmp_mat,3,zero,mat3,3)
1116 :
1117 129 : end subroutine matpointsym
1118 : !!***
1119 :
1120 : !!****f* m_symtk/holocell
1121 : !! NAME
1122 : !! holocell
1123 : !!
1124 : !! FUNCTION
1125 : !! Examine whether the trial conventional cell described by cell_base
1126 : !! is coherent with the required holohedral group.
1127 : !! Possibly enforce the holohedry and modify the basis vectors.
1128 : !! Note: for iholohedry=4, the tetragonal axis is not required to be along the C axis.
1129 : !!
1130 : !! INPUTS
1131 : !! enforce= if 0, only check; if =1, enforce exactly the holohedry
1132 : !! iholohedry=required holohegral group (uses its absolute value, since when the multiplicity of the cell is
1133 : !! more than one, the sign of iholohedry is changed).
1134 : !! iholohedry=1 triclinic 1bar
1135 : !! iholohedry=2 monoclinic 2/m
1136 : !! iholohedry=3 orthorhombic mmm
1137 : !! iholohedry=4 tetragonal 4/mmm
1138 : !! iholohedry=5 trigonal 3bar m
1139 : !! iholohedry=6 hexagonal 6/mmm
1140 : !! iholohedry=7 cubic m3bar m
1141 : !! tolsym=tolerance for the symmetry operations
1142 : !!
1143 : !! OUTPUT
1144 : !! foundc=1 if the basis vectors supports the required holohedry ; =0 otherwise
1145 : !!
1146 : !! SIDE EFFECTS
1147 : !! cell_base(3,3)=basis vectors of the conventional cell (changed if enforce==1, otherwise unchanged)
1148 : !!
1149 : !! SOURCE
1150 :
1151 71745 : subroutine holocell(cell_base,enforce,foundc,iholohedry,tolsym)
1152 :
1153 : !Arguments ------------------------------------
1154 : !scalars
1155 : integer,intent(in) :: enforce,iholohedry
1156 : integer,intent(out) :: foundc
1157 : real(dp),intent(in) :: tolsym
1158 : !arrays
1159 : real(dp),intent(inout) :: cell_base(3,3)
1160 :
1161 : !Local variables ------------------------------
1162 : !scalars
1163 : integer :: allequal,ii,orth
1164 : real(dp):: aa,scprod1
1165 : character(len=500) :: msg
1166 : !arrays
1167 : integer :: ang90(3),equal(3)
1168 : real(dp) :: length(3),metric(3,3),norm(3),rbasis(3,3),rconv(3,3),rconv_new(3,3)
1169 : real(dp) :: rnormalized(3,3),symmetrized_length(3)
1170 : !**************************************************************************
1171 :
1172 71745 : if(abs(iholohedry)<1 .or. abs(iholohedry)>7)then
1173 0 : write(msg, '(a,i0)' )'Abs(iholohedry) should be between 1 and 7, while iholohedry=',iholohedry
1174 0 : ABI_BUG(msg)
1175 : end if
1176 :
1177 286980 : do ii=1,3
1178 : metric(:,ii)=cell_base(1,:)*cell_base(1,ii)+&
1179 : & cell_base(2,:)*cell_base(2,ii)+&
1180 932685 : & cell_base(3,:)*cell_base(3,ii)
1181 : end do
1182 :
1183 : !Examine the angles and vector lengths
1184 71745 : ang90(:)=0
1185 71745 : if(metric(1,2)**2<tolsym**2*metric(1,1)*metric(2,2))ang90(3)=1
1186 71745 : if(metric(1,3)**2<tolsym**2*metric(1,1)*metric(3,3))ang90(2)=1
1187 71745 : if(metric(2,3)**2<tolsym**2*metric(2,2)*metric(3,3))ang90(1)=1
1188 71745 : orth=0
1189 71745 : if(ang90(1)==1 .and. ang90(2)==1 .and. ang90(3)==1) orth=1
1190 71745 : equal(:)=0
1191 71745 : if(abs(metric(1,1)-metric(2,2))<tolsym*half*(metric(1,1)+metric(2,2)))equal(3)=1
1192 71745 : if(abs(metric(1,1)-metric(3,3))<tolsym*half*(metric(1,1)+metric(3,3)))equal(2)=1
1193 71745 : if(abs(metric(2,2)-metric(3,3))<tolsym*half*(metric(2,2)+metric(3,3)))equal(1)=1
1194 71745 : allequal=0
1195 71745 : if(equal(1)==1 .and. equal(2)==1 .and. equal(3)==1) allequal=1
1196 :
1197 : !DEBUG
1198 : !write(std_out,*)' holocell : enforce, iholohedry=',enforce, iholohedry
1199 : !write(std_out,*)' holocell : ang90=',ang90
1200 : !write(std_out,*)' holocell : equal=',equal
1201 : !!write(std_out,*)' holocell : tolsym=',tolsym
1202 : !!write(std_out,*)' holocell : metric(1,1)=',metric(1,1)
1203 : !!write(std_out,*)' holocell : metric(1,2)=',metric(1,2)
1204 : !ENDDEBUG
1205 :
1206 71745 : foundc=0
1207 71745 : if(abs(iholohedry)==1) foundc=1
1208 71745 : if(abs(iholohedry)==2 .and. ang90(1)+ang90(3)==2 ) foundc=1
1209 71745 : if(abs(iholohedry)==3 .and. orth==1) foundc=1
1210 71745 : if(abs(iholohedry)==4 .and. orth==1 .and. &
1211 5357 : & (equal(3)==1 .or. equal(2)==1 .or. equal(1)==1) ) foundc=1
1212 : if(abs(iholohedry)==5 .and. allequal==1 .and. &
1213 71745 : & (abs(metric(1,2)-metric(2,3))<tolsym*metric(2,2)) .and. &
1214 1573 : & (abs(metric(1,2)-metric(1,3))<tolsym*metric(1,1)) ) foundc=1
1215 : if(abs(iholohedry)==6 .and. equal(3)==1 .and. &
1216 71745 : & ang90(1)==1 .and. ang90(2)==1 .and. &
1217 2193 : & abs(2*metric(1,2)+metric(1,1))<tolsym*metric(1,1) ) foundc=1
1218 71745 : if(abs(iholohedry)==7 .and. orth==1 .and. allequal==1) foundc=1
1219 :
1220 : !write(std_out, '(a,2i4)' )' holocell : foundc, enforce=',foundc,enforce
1221 :
1222 : !-------------------------------------------------------------------------------------
1223 : !Possibly enforce the holohedry (if it is to be enforced !)
1224 :
1225 71745 : if(foundc==0.and.enforce==1.and.abs(iholohedry)/=1)then
1226 :
1227 : ! Copy the cell_base vectors, and possibly fix the tetragonal axis to be the c-axis
1228 : ! XG20201016 WARNING : in principle, one should NOT use the 'equal' information, since precisely this enforcement
1229 : ! has the aim to reinstall the symmetries while they are broken !!
1230 29 : if(abs(iholohedry)==4.and.equal(1)==1)then
1231 0 : rconv(:,3)=cell_base(:,1) ; rconv(:,1)=cell_base(:,2) ; rconv(:,2)=cell_base(:,3)
1232 29 : else if (abs(iholohedry)==4.and.equal(2)==1)then
1233 0 : rconv(:,3)=cell_base(:,2) ; rconv(:,2)=cell_base(:,1) ; rconv(:,1)=cell_base(:,3)
1234 : else
1235 29 : rconv(:,:)=cell_base(:,:)
1236 : end if
1237 :
1238 : ! Compute the length of the three conventional vectors
1239 116 : length(1)=sqrt(sum(rconv(:,1)**2))
1240 116 : length(2)=sqrt(sum(rconv(:,2)**2))
1241 116 : length(3)=sqrt(sum(rconv(:,3)**2))
1242 :
1243 : ! Take care of the first conventional vector aligned with rbasis(:,3) (or aligned with the trigonal axis if rhombohedral)
1244 : ! and choice of the first normalized direction
1245 29 : if(abs(iholohedry)==5)then
1246 4 : rbasis(:,3)=third*(rconv(:,1)+rconv(:,2)+rconv(:,3))
1247 : else
1248 112 : rbasis(:,3)=rconv(:,3)
1249 : end if
1250 116 : norm(3)=sqrt(sum(rbasis(:,3)**2))
1251 116 : rnormalized(:,3)=rbasis(:,3)/norm(3)
1252 :
1253 : ! Projection of the first conventional vector perpendicular to rbasis(:,3)
1254 : ! and choice of the first normalized direction
1255 116 : scprod1=sum(rnormalized(:,3)*rconv(:,1))
1256 116 : rbasis(:,1)=rconv(:,1)-rnormalized(:,3)*scprod1
1257 116 : norm(1)=sqrt(sum(rbasis(:,1)**2))
1258 116 : rnormalized(:,1)=rbasis(:,1)/norm(1)
1259 :
1260 : ! Generation of the second vector, perpendicular to the third and first
1261 29 : rnormalized(1,2)=rnormalized(2,3)*rnormalized(3,1)-rnormalized(3,3)*rnormalized(2,1)
1262 29 : rnormalized(2,2)=rnormalized(3,3)*rnormalized(1,1)-rnormalized(1,3)*rnormalized(3,1)
1263 29 : rnormalized(3,2)=rnormalized(1,3)*rnormalized(2,1)-rnormalized(2,3)*rnormalized(1,1)
1264 :
1265 : ! Compute the vectors of the conventional cell, on the basis of iholohedry
1266 : if(abs(iholohedry)==2)then
1267 28 : rconv_new(:,3)=rconv(:,3)
1268 28 : rconv_new(:,1)=rconv(:,1)
1269 28 : rconv_new(:,2)=rnormalized(:,2)*length(2) ! Now, the y axis is perpendicular to the two others, that have not been changed
1270 : else if(abs(iholohedry)==3.or.abs(iholohedry)==4.or.abs(iholohedry)==7)then
1271 14 : if(abs(iholohedry)==7)then
1272 35 : symmetrized_length(1:3)=sum(length(:))*third
1273 9 : else if(abs(iholohedry)==4)then
1274 5 : symmetrized_length(3)=length(3)
1275 15 : symmetrized_length(1:2)=half*(length(1)+length(2))
1276 4 : else if(abs(iholohedry)==3)then
1277 4 : symmetrized_length(:)=length(:)
1278 : end if
1279 56 : do ii=1,3
1280 182 : rconv_new(:,ii)=rnormalized(:,ii)*symmetrized_length(ii)
1281 : end do
1282 : else if(abs(iholohedry)==5)then
1283 : ! In the normalized basis, they have coordinates (a,0,c), and (-a/2,+-sqrt(3)/2*a,c)
1284 : ! c is known, but a is computed from the knowledge of the average length of the initial vectors
1285 4 : aa=sqrt(sum(length(:)**2)*third-norm(3)**2)
1286 4 : rconv_new(:,1)=aa*rnormalized(:,1)+rbasis(:,3)
1287 4 : rconv_new(:,2)=aa*half*(-rnormalized(:,1)+sqrt(three)*rnormalized(:,2))+rbasis(:,3)
1288 4 : rconv_new(:,3)=aa*half*(-rnormalized(:,1)-sqrt(three)*rnormalized(:,2))+rbasis(:,3)
1289 : else if(abs(iholohedry)==6)then
1290 :
1291 : ! In the normalized basis, they have coordinates (a,0,0), (-a/2,+-sqrt(3)/2*a,0), and (0,0,c)
1292 : ! c is known, but a is computed from the knowledge of the average length of the initial vectors
1293 7 : aa=half*(length(1)+length(2))
1294 28 : rconv_new(:,1)=aa*rnormalized(:,1)
1295 28 : rconv_new(:,2)=aa*half*(-rnormalized(:,1)+sqrt(three)*rnormalized(:,2))
1296 28 : rconv_new(:,3)=rconv(:,3)
1297 : end if
1298 :
1299 : ! Copy back the cell_base vectors
1300 29 : if(abs(iholohedry)==4.and.equal(1)==1)then
1301 0 : cell_base(:,3)=rconv_new(:,2) ; cell_base(:,2)=rconv_new(:,1) ; cell_base(:,1)=rconv_new(:,3)
1302 29 : else if (abs(iholohedry)==4.and.equal(2)==1)then
1303 0 : cell_base(:,3)=rconv_new(:,1) ; cell_base(:,1)=rconv_new(:,2) ; cell_base(:,2)=rconv_new(:,3)
1304 : else
1305 29 : cell_base(:,:)=rconv_new(:,:)
1306 : end if
1307 :
1308 : end if
1309 :
1310 71745 : end subroutine holocell
1311 : !!***
1312 :
1313 : !!****f* m_symtk/symmetrize_rprimd
1314 : !! NAME
1315 : !! symmetrize_rprimd
1316 : !!
1317 : !! FUNCTION
1318 : !! Supposing the input rprimd does not preserve the length and angles
1319 : !! following the symmetries, will generates a new set rprimd,
1320 : !! on the basis of the expected characteristics of the conventional cell, as specified in bravais(:)
1321 : !!
1322 : !! INPUTS
1323 : !! bravais(11): bravais(1)=iholohedry
1324 : !! bravais(2)=center
1325 : !! bravais(3:11)=coordinates of rprimd in the axes
1326 : !! of the conventional bravais lattice (*2 if center/=0)
1327 : !! nsym=actual number of symmetries
1328 : !! symrel(3,3,1:nsym)=symmetry operations in real space in terms of primitive translations
1329 : !! tolsym=tolerance for the symmetry operations (only for checking purposes, the new set rprimd will
1330 : !! be coherent with the symmetry operations at a much accurate level).
1331 : !!
1332 : !! SIDE EFFECTS
1333 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1334 : !!
1335 : !! SOURCE
1336 :
1337 29 : subroutine symmetrize_rprimd(bravais,nsym,rprimd,symrel,tolsym)
1338 :
1339 : !Arguments ------------------------------------
1340 : !scalars
1341 : integer,intent(in) :: nsym
1342 : real(dp),intent(in) :: tolsym
1343 : !arrays
1344 : integer,intent(in) :: bravais(11),symrel(3,3,nsym)
1345 : real(dp),intent(inout) :: rprimd(3,3)
1346 :
1347 : !Local variables-------------------------------
1348 : !scalars
1349 : integer :: foundc,iexit,ii,jj
1350 : real(dp):: rprimd_maxabs
1351 : !character(len=500) :: msg
1352 : !arrays
1353 : real(dp):: aa(3,3),ait(3,3),cell_base(3,3),gprimd(3,3),rmet(3,3),rprimd_new(3,3)
1354 : ! *************************************************************************
1355 :
1356 : !write(std_out,'(a)') ' symmetrize_rprimd : enter '
1357 :
1358 : !Build the conventional cell basis vectors in cartesian coordinates
1359 116 : aa(:,1)=bravais(3:5)
1360 116 : aa(:,2)=bravais(6:8)
1361 116 : aa(:,3)=bravais(9:11)
1362 : !Inverse transpose
1363 29 : call matr3inv(aa,ait)
1364 116 : do ii=1,3
1365 377 : cell_base(:,ii)=ait(ii,1)*rprimd(:,1)+ait(ii,2)*rprimd(:,2)+ait(ii,3)*rprimd(:,3)
1366 : end do
1367 :
1368 : !write(std_out,'(a)') ' before holocell, cell_base ='
1369 : !do ii=1,3
1370 : ! write(std_out,'(3es16.8)') cell_base(:,ii)
1371 : !enddo
1372 :
1373 : !Enforce the proper holohedry on the conventional cell vectors.
1374 29 : call holocell(cell_base,1,foundc,bravais(1),tolsym)
1375 :
1376 : !write(std_out,'(a)') ' after holocell, cell_base ='
1377 : !do ii=1,3
1378 : ! write(std_out,'(3es16.8)') cell_base(:,ii)
1379 : !enddo
1380 :
1381 : !Reconstruct the dimensional primitive vectors
1382 116 : do ii=1,3
1383 377 : rprimd_new(:,ii)=aa(1,ii)*cell_base(:,1)+aa(2,ii)*cell_base(:,2)+aa(3,ii)*cell_base(:,3)
1384 : end do
1385 :
1386 : !Suppress meaningless values
1387 377 : rprimd_maxabs=maxval(abs(rprimd_new))
1388 116 : do ii=1,3
1389 377 : do jj=1,3
1390 348 : if(abs(rprimd(ii,jj))<tol12*rprimd_maxabs)rprimd(ii,jj)=zero
1391 : enddo
1392 : enddo
1393 :
1394 29 : rprimd(:,:)=rprimd_new(:,:)
1395 :
1396 : !Check whether the symmetry operations are consistent with the lattice vectors
1397 1160 : rmet = MATMUL(TRANSPOSE(rprimd), rprimd)
1398 29 : call matr3inv(rprimd, gprimd)
1399 : !call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1400 29 : iexit=0
1401 :
1402 29 : call chkorthsy(gprimd,iexit,nsym,rmet,rprimd,symrel,tolsym)
1403 :
1404 : !write(std_out,'(a)') ' symmetrize_rprimd : exit '
1405 :
1406 29 : end subroutine symmetrize_rprimd
1407 : !!***
1408 :
1409 : !!****f* m_symtk/symmetrize_tnons
1410 : !! NAME
1411 : !! symmetrize_tnons
1412 : !!
1413 : !! FUNCTION
1414 : !! Given the order of a symmetry operation, make sure that tnons is
1415 : !! such that applying "order" times the symmetry operation
1416 : !! generate the unity operation, accurately.
1417 : !!
1418 : !! INPUTS
1419 : !! nsym=actual number of symmetries
1420 : !! symrel(3,3,1:nsym)=symmetry operations in real space in terms of primitive translations
1421 : !! tolsym=tolerance that was used to determine the symmetry operations
1422 : !!
1423 : !! SIDE EFFECTS
1424 : !! tnons(3,1:nsym)= non-symmorphic translation vectors
1425 : !!
1426 : !! SOURCE
1427 :
1428 8863 : subroutine symmetrize_tnons(nsym,symrel,tnons,tolsym)
1429 :
1430 : !Arguments ------------------------------------
1431 : !scalars
1432 : integer,intent(in) :: nsym
1433 : real(dp),intent(in) :: tolsym
1434 : !arrays
1435 : integer,intent(in) :: symrel(3,3,nsym)
1436 : real(dp),intent(inout) :: tnons(3,nsym)
1437 :
1438 : !Local variables-------------------------------
1439 : !scalars
1440 : integer :: iorder,isym,order
1441 : !character(len=500) :: msg
1442 : !arrays
1443 : integer :: symrel_mult(3,3)
1444 : integer :: unitmat(3,3)
1445 : real(dp):: tnons_mult(3)
1446 : ! *************************************************************************
1447 :
1448 : !write(std_out,'(a)') ' symmetrize_tnons : enter '
1449 :
1450 8863 : unitmat=0
1451 8863 : unitmat(1,1)=1 ; unitmat(2,2)=1 ; unitmat(3,3)=1
1452 :
1453 297879 : do isym=1,nsym
1454 : !write(std_out,'(a,i4,9i3,3es16.6)') ' isym,symrel,tnons=',isym,symrel(:,:,isym),tnons(:,isym)
1455 :
1456 4624256 : symrel_mult(:,:)=symrel(:,:,isym) ; tnons_mult(:)=tnons(:,isym)
1457 648234 : order=0
1458 : !Determine the order of the operation
1459 648234 : do iorder=1,48
1460 25929360 : symrel_mult(:,:)=matmul(symrel(:,:,isym),symrel_mult(:,:))
1461 18150552 : tnons_mult(:)=matmul(symrel(:,:,isym),tnons_mult(:))+tnons(:,isym)
1462 8427042 : if(sum(abs(symrel_mult-unitmat))==0)then
1463 : if(abs(tnons_mult(1)-nint(tnons_mult(1)))<tolsym*iorder .and. &
1464 293515 : abs(tnons_mult(2)-nint(tnons_mult(2)))<tolsym*iorder .and. &
1465 : abs(tnons_mult(3)-nint(tnons_mult(3)))<tolsym*iorder)then
1466 : !The order has been found
1467 289016 : order=iorder+1
1468 : !Now, adjust the tnons vector, in order to obtain the exact identity operation at order "order"
1469 1156064 : tnons_mult(:)=(tnons_mult(:)-nint(tnons_mult(:)))/(dble(order))
1470 289016 : if(abs(tnons_mult(1))>1.00001e-8) tnons(1,isym)=tnons(1,isym)-tnons_mult(1)
1471 289016 : if(abs(tnons_mult(2))>1.00001e-8) tnons(2,isym)=tnons(2,isym)-tnons_mult(2)
1472 289016 : if(abs(tnons_mult(3))>1.00001e-8) tnons(3,isym)=tnons(3,isym)-tnons_mult(3)
1473 289016 : exit
1474 : endif
1475 : endif
1476 :
1477 : enddo ! iorder
1478 :
1479 297879 : if(order==0)then
1480 0 : ABI_BUG("Was unable to find order of operation")
1481 : endif
1482 : enddo
1483 :
1484 : !write(std_out,'(a)') ' symmetrize_tnons : exit '
1485 :
1486 8863 : end subroutine symmetrize_tnons
1487 : !!***
1488 :
1489 : !!****f* m_symtk/symmetrize_xred
1490 : !! NAME
1491 : !! symmetrize_xred
1492 : !!
1493 : !! FUNCTION
1494 : !! Symmetrize atomic coordinates.
1495 : !! Two tasks can be executed :
1496 : !! A. If optional argument indsym is present.
1497 : !! Using input symmetry matrices symrel
1498 : !! which are expressed in terms of the basis of real space primitive
1499 : !! translations (array elements are integers), use indsym to make all corresponding atoms coordinate
1500 : !! fullfil exactly symmetry operations.
1501 : !! Input array indsym(4,isym,iatom) gives label of atom into which iatom
1502 : !! is rotated by INVERSE of symmetry element isym and also gives primitive
1503 : !! translation to get back to unit cell.
1504 : !! This version uses improvement in algorithm suggested by Andrew
1505 : !! Horsfield (see symatm.f).
1506 : !! B. If optional argument tolsym AND tnons_new are defined.
1507 : !! Might also adjust xred in order for tnons to be aligned with the FFT grids.
1508 : !! This will deliver new tnons_new.
1509 : !! NOTE : Actually, should make two separate routines !
1510 : !!
1511 : !! INPUTS
1512 : !! indsym(4,nsym,natom)=(optional) indirect indexing array giving label of atom
1513 : !! into which iatom is rotated by symmetry element isym
1514 : !! natom=number of atoms
1515 : !! nsym=number of symmetries in group
1516 : !! symrel(3,3,nsym)=symmetry matrices in terms of real space
1517 : !! primitive translations
1518 : !! tnons(3,nsym)=nonsymmorphic translations for symmetries
1519 : !! tolsym=(optional) tolerance on symmetries. When defined, one will try to align the symmetry operations with the FFT grid,
1520 : !! if the modification is less than tolsym. Take tolsym equal to 1 to deliver possibly large changes of xred,
1521 : !! giving suggestions of xred modifications, to be proposed to users.
1522 : !!
1523 : !! OUTPUT
1524 : !! tnons_new(3,nsym)=(optional)nonsymmorphic translations for symmetries
1525 : !!
1526 : !! SIDE EFFECTS
1527 : !! Input/Output
1528 : !! xred(3,natom)=
1529 : !! (input) atomic coordinates in terms of real space translations
1530 : !! (output) symmetrized atomic coordinates in terms
1531 : !! of real space translations
1532 : !! fixed_mismatch=(optional) At input, needs to be present for tnons_new to be computed
1533 : !! At output : 1 if there is a mismatch and this mismatch has been fixed, 0 otherwise
1534 : !! mismatch_fft_tnons=(optional) At input, needs to be present for tnons_new to be computed
1535 : !! Atd output : non-zero if there is a mismatch between the fft grid and the tnons, gives the number
1536 : !! of the first symmetry operation for which there is such a mismatch. Zero otherwise.
1537 : !!
1538 : !! SOURCE
1539 :
1540 28190 : subroutine symmetrize_xred(natom,nsym,symrel,tnons,xred,fixed_mismatch,indsym,mismatch_fft_tnons,tnons_new,tolsym)
1541 :
1542 : !Arguments ------------------------------------
1543 : !scalars
1544 : integer,intent(in) :: natom,nsym
1545 : integer,intent(out),optional :: fixed_mismatch,mismatch_fft_tnons
1546 : !arrays
1547 : integer,intent(in),optional :: indsym(4,nsym,natom)
1548 : integer,intent(in) :: symrel(3,3,nsym)
1549 : real(dp),intent(in) :: tnons(3,nsym)
1550 : real(dp),intent(in),optional :: tolsym
1551 : real(dp),intent(out),optional :: tnons_new(3,nsym)
1552 : real(dp),intent(inout) :: xred(3,natom)
1553 :
1554 : !Local variables-------------------------------
1555 : !scalars
1556 : integer :: iatom,ib,ii,info,irank,isym,isym2
1557 : integer :: jj,mismatch_fft_tnons_current
1558 : real(dp) :: diff
1559 : logical :: dissimilar
1560 : !arrays
1561 : real(dp) :: delta(4),fc(3),mat(3,3),mult(4)=(/eight,nine,ten,three*four/)
1562 : real(dp) :: sgval(3),tsum(3),tt(3),work(15),xredshift(3,1)
1563 43512 : real(dp),allocatable :: xredsym(:,:)
1564 : real(dp) :: transl(3) ! translation vector
1565 : ! *************************************************************************
1566 :
1567 : !Check whether group contains more than identity, if not then simply return after possible copying.
1568 43512 : if(present(tnons_new))then
1569 1791642 : tnons_new(:,1:nsym)=tnons(:,1:nsym)
1570 : endif
1571 43512 : if(present(fixed_mismatch))fixed_mismatch=0
1572 43512 : if(present(mismatch_fft_tnons))mismatch_fft_tnons=0
1573 :
1574 43512 : if (nsym>1) then
1575 : ! write(std_out,*)
1576 : ! write(std_out,'(a,i4)') 'symmetrize_xred: enter, nsym=',nsym
1577 : ! do iatom=1,natom
1578 : ! write(std_out,'(a,i4,3es16.6)') 'iatom,xred=',iatom,xred(:,iatom)
1579 : ! enddo
1580 : ! do isym=1,nsym
1581 : ! write(std_out,'(a,i4,9i3,3es16.6)') 'isym,symrel,tnons',isym,symrel(:,:,isym),tnons(:,isym)
1582 : ! enddo
1583 : ! write(std_out,*)' present(tnons_new),present(tolsym)=',present(tnons_new),present(tolsym)
1584 : ! write(std_out,*)
1585 :
1586 87639 : ABI_MALLOC(xredsym,(3,natom))
1587 390037 : xredsym(:,:)=xred(:,1:natom)
1588 :
1589 29213 : if(present(indsym))then
1590 :
1591 : ! Loop over atoms to determine new, symmetrized positions.
1592 51438 : do iatom=1,natom
1593 35998 : tsum(:)=0.0d0
1594 : !
1595 : ! Loop over symmetries
1596 1020600 : do isym=1,nsym
1597 : ! atom ib is atom into which iatom is rotated by inverse of
1598 : ! symmetry isym (inverse of symrel(mu,nu,isym))
1599 984602 : ib=indsym(4,isym,iatom)
1600 : ! Find the reduced coordinates after translation=t(indsym)+transl
1601 3938408 : fc(:)=xred(:,ib)+dble(indsym(1:3,isym,iatom))
1602 : ! Compute [S * (x(indsym)+transl) ] + tnonsymmorphic
1603 : tt(:)=dble(symrel(:,1,isym))*fc(1)+&
1604 : dble(symrel(:,2,isym))*fc(2)+&
1605 3938408 : dble(symrel(:,3,isym))*fc(3)+ tnons(:,isym)
1606 :
1607 : ! Average over nominally equivalent atomic positions
1608 3974406 : tsum(:)=tsum(:)+tt(:)
1609 : end do ! isym
1610 : !
1611 : ! Set symmetrized result to sum over number of terms
1612 159432 : xredsym(:,iatom)=tsum(:)/dble(nsym)
1613 :
1614 : end do ! iatom
1615 : endif ! present(indsym)
1616 :
1617 : !do iatom=1,natom
1618 : ! write(std_out,'(a,i4,3es20.10)') 'iatom,xredsym=',iatom,xredsym(:,iatom)
1619 : !enddo
1620 :
1621 : ! Loop over symmetry operations to determine possibly new tnons, as well as symmetrized positions.
1622 29213 : if(present(tolsym) .and. present(tnons_new) .and. present(fixed_mismatch) .and. present(mismatch_fft_tnons))then
1623 13773 : fixed_mismatch=0
1624 13773 : mismatch_fft_tnons=0
1625 : !The use of tolsym here is only to favor 0.5 over -0.5
1626 1783897 : tnons_new(:,:)=tnons(:,:)-nint(tnons(:,:)-tolsym)
1627 455467 : do isym=1,nsym
1628 : mismatch_fft_tnons_current=0
1629 1766884 : do ii=1,3
1630 6625815 : delta(:)=tnons(ii,isym)*mult(:)
1631 6625815 : delta(:)=delta(:)-nint(delta(:))
1632 : ! Is there is a mismatch between FFT and isym for all multipliers ?
1633 1982134 : if( all(abs(delta(:))>tol8*mult(:)) ) mismatch_fft_tnons_current=1
1634 : enddo
1635 : !Declare the first symmetry operation that induces a problem
1636 441721 : if(mismatch_fft_tnons_current>0 .and. (mismatch_fft_tnons==0)) mismatch_fft_tnons=isym
1637 :
1638 : ! However, also try to propose a solution.
1639 13773 : if(mismatch_fft_tnons_current==1)then
1640 : ! Compute the pseudo-inverse of symrel-1, then multiply tnons
1641 23329 : mat(:,:)=zero; mat(1,1)=one; mat(2,2)=one; mat(3,3)=one
1642 : ! This is symrel-1
1643 303277 : mat(:,:)=symrel(:,:,isym)-mat(:,:)
1644 93316 : do ii=1,3
1645 : !Select the smallest modification tnons
1646 349935 : delta(:)=tnons(ii,isym)*mult(:)
1647 349935 : delta(:)=(delta(:)-nint(delta(:)))/mult(:)
1648 69987 : xredshift(ii,1)=delta(1)
1649 303277 : do jj=2,4
1650 279948 : if(abs(delta(jj))<abs(xredshift(ii,1)))xredshift(ii,1)=delta(jj)
1651 : enddo
1652 : enddo
1653 23329 : call dgelss(3,3,1,mat,3,xredshift(:,1),3,sgval,tol5,irank,work,15,info)
1654 :
1655 : ! xredshift(:,1) is now the tentative shift, to be tested for all symmetries
1656 26238 : if( all(abs(xredshift(:,1))<tolsym) )then
1657 29 : fixed_mismatch=1
1658 913 : do isym2=1, nsym
1659 27466 : tnons_new(:,isym2)=tnons(:,isym2)+xredshift(:,1)-matmul(symrel(:,:,isym2),xredshift(:,1))
1660 3544 : do ii=1,3
1661 : !tnons might now be slighly non-zero. Set to zero such values
1662 2658 : if(abs(tnons_new(ii,isym2))<tol6**2)tnons_new(ii,isym2)=zero
1663 13290 : delta(:)=tnons_new(ii,isym2)*mult(:)
1664 13290 : delta(:)=delta(:)-nint(delta(:))
1665 : ! Is the mismatch between FFT and symmetries still present for all the multipliers ??
1666 3552 : if( all(abs(delta(:))>tol8*mult(:)) ) fixed_mismatch=0
1667 : enddo
1668 913 : if(fixed_mismatch==0)exit
1669 : enddo
1670 : endif
1671 23329 : if(fixed_mismatch==1)exit
1672 : endif ! mismatch_fft_tnons_current==1
1673 : end do ! isym
1674 13773 : if(mismatch_fft_tnons/=0)then
1675 1011 : if(fixed_mismatch==1)then
1676 196 : do iatom=1,natom
1677 703 : xredsym(:,iatom)=xredsym(:,iatom)+xredshift(:,1)
1678 : enddo
1679 : endif
1680 : endif
1681 :
1682 : ! write(std_out,*) ' '
1683 : ! write(std_out,*) ' mismatch_fft_tnons, fixed_mismatch=',mismatch_fft_tnons, fixed_mismatch
1684 : ! do iatom=1,natom
1685 : ! write(std_out,'(a,i4,3es20.10)') ' iatom,xredsym=',iatom,xredsym(:,iatom)
1686 : ! enddo
1687 :
1688 : endif ! present(tolsym) .and. present(tnons_new)
1689 :
1690 : ! --------------------------------------------------------------
1691 : ! Will update the atomic positions only if it is worth to do so.
1692 :
1693 116852 : transl(:)=xredsym(:,1)-nint(xredsym(:,1))
1694 :
1695 : ! Compute the smallest translation to an integer
1696 90206 : do jj=2,natom
1697 273185 : do ii=1,3
1698 182979 : diff=xredsym(ii,jj)-nint(xredsym(ii,jj))
1699 243972 : if (diff<transl(ii)) transl(ii)=diff
1700 : end do
1701 : end do
1702 :
1703 : ! Test if the translation on each direction is small
1704 : ! Tolerance 1E-13
1705 116852 : do ii=1,3
1706 116852 : if (abs(transl(ii))>1e-13) transl(ii)=0.0
1707 : end do
1708 :
1709 : ! Execute translation
1710 119419 : do jj=1,natom
1711 390037 : do ii=1,3
1712 360824 : xredsym(ii,jj)=xredsym(ii,jj)-transl(ii)
1713 : end do
1714 : end do
1715 :
1716 : ! Test if xredsym is too similar to xred
1717 : ! Tolerance 1E-15
1718 : dissimilar=.FALSE.
1719 119419 : do jj=1,natom
1720 390037 : do ii=1,3
1721 360824 : if (abs(xredsym(ii,jj)-xred(ii,jj))>1E-15) dissimilar=.TRUE.
1722 : end do
1723 : end do
1724 :
1725 35209 : if (dissimilar) xred(:,:)=xredsym(:,:)
1726 29213 : ABI_FREE(xredsym)
1727 :
1728 : ! End condition of nsym/=1
1729 : end if
1730 :
1731 : !write(std_out,*)
1732 : !write(std_out,'(a)') 'symmetrize_xred : exit'
1733 : !do iatom=1,natom
1734 : ! write(std_out,'(a,i4,3es20.10)') 'iatom,xred=',iatom,xred(:,iatom)
1735 : !enddo
1736 : !if(present(tnons_new))then
1737 : ! do isym=1,nsym
1738 : ! write(std_out,'(a,i4,9i3,3es20.10)') 'isym,symrel,tnons_new',isym,symrel(:,:,isym),tnons_new(:,isym)
1739 : ! enddo
1740 : !endif
1741 : !write(std_out,*)
1742 :
1743 43512 : end subroutine symmetrize_xred
1744 : !!***
1745 :
1746 : !!****f* m_symtk/symchk
1747 : !! NAME
1748 : !! symchk
1749 : !!
1750 : !! FUNCTION
1751 : !! Symmetry checker for atomic coordinates.
1752 : !! Checks for translated atomic coordinate tratom(3) to agree
1753 : !! with some coordinate xred(3,iatom) where atomic types agree too.
1754 : !! All coordinates are "reduced", i.e. given in terms of primitive
1755 : !! reciprocal translations.
1756 : !!
1757 : !! INPUTS
1758 : !! natom=number of atoms in unit cell
1759 : !! tratom(3)=reduced coordinates for a single atom which presumably
1760 : !! result from the application of a symmetry operation to an atomic
1761 : !! coordinate
1762 : !! trtypat=type of atom (integer) translated to tratom
1763 : !! typat(natom)=types of all atoms in unit cell (integer)
1764 : !! xred(3,natom)=reduced coordinates for all atoms in unit cell
1765 : !!
1766 : !! OUTPUT
1767 : !! difmin(3)=minimum difference between apparently equivalent atoms
1768 : !! (give value separately for each coordinate)--note that value
1769 : !! may be NEGATIVE so take abs later if needed
1770 : !! eatom=atom label of atom which is SAME as tratom to within a primitive
1771 : !! cell translation ("equivalent atom")
1772 : !! transl(3)=primitive cell translation to make iatom same as tratom (integers)
1773 : !!
1774 : !! SOURCE
1775 :
1776 3116497 : subroutine symchk(difmin,eatom,natom,tratom,transl,trtypat,typat,xred)
1777 :
1778 : !Arguments ------------------------------------
1779 : !scalars
1780 : integer,intent(in) :: natom,trtypat
1781 : integer,intent(out) :: eatom
1782 : !arrays
1783 : integer,intent(in) :: typat(natom)
1784 : integer,intent(out) :: transl(3)
1785 : real(dp),intent(in) :: tratom(3),xred(3,natom)
1786 : real(dp),intent(out) :: difmin(3)
1787 :
1788 : !Local variables-------------------------------
1789 : !scalars
1790 : integer :: iatom,jatom,trans1,trans2,trans3
1791 : real(dp) :: test,test1,test2,test3,testmn
1792 : ! *************************************************************************
1793 :
1794 : ! write(std_out,'(a,a,i4,3f18.12)') ch10,' symchk : enter, trtypat,tratom=',trtypat,tratom
1795 :
1796 : !Start testmn out at large value
1797 3116497 : testmn=1000000.d0
1798 :
1799 : !Loop through atoms. when types agree, check for agreement after primitive translation
1800 3116497 : jatom=1
1801 8802923 : do iatom=1,natom
1802 8798161 : if (trtypat/=typat(iatom)) cycle
1803 :
1804 : ! Check all three components
1805 7254844 : test1=tratom(1)-xred(1,iatom)
1806 7254844 : test2=tratom(2)-xred(2,iatom)
1807 7254844 : test3=tratom(3)-xred(3,iatom)
1808 : ! Find nearest integer part of difference
1809 7254844 : trans1=nint(test1)
1810 7254844 : trans2=nint(test2)
1811 7254844 : trans3=nint(test3)
1812 : ! Check whether, after translation, they agree
1813 7254844 : test1=test1-dble(trans1)
1814 7254844 : test2=test2-dble(trans2)
1815 7254844 : test3=test3-dble(trans3)
1816 7254844 : test=abs(test1)+abs(test2)+abs(test3)
1817 7259606 : if (test<tol10) then
1818 : ! Note that abs() is not taken here
1819 3111735 : difmin(1)=test1
1820 3111735 : difmin(2)=test2
1821 3111735 : difmin(3)=test3
1822 3111735 : jatom=iatom
1823 3111735 : transl(1)=trans1
1824 3111735 : transl(2)=trans2
1825 3111735 : transl(3)=trans3
1826 : ! Break out of loop when agreement is within tolerance
1827 3111735 : exit
1828 : else
1829 : ! Keep track of smallest difference if greater than tol10
1830 4143109 : if (test<testmn) then
1831 1610680 : testmn=test
1832 : ! Note that abs() is not taken here
1833 1610680 : difmin(1)=test1
1834 1610680 : difmin(2)=test2
1835 1610680 : difmin(3)=test3
1836 1610680 : jatom=iatom
1837 1610680 : transl(1)=trans1
1838 1610680 : transl(2)=trans2
1839 1610680 : transl(3)=trans3
1840 : end if
1841 : end if
1842 :
1843 : end do ! End loop over iatom. Note a "cycle" and an "exit" inside the loop
1844 :
1845 3116497 : eatom=jatom
1846 :
1847 3116497 : end subroutine symchk
1848 : !!***
1849 :
1850 : !!****f* m_symtk/symatm
1851 : !! NAME
1852 : !! symatm
1853 : !!
1854 : !! FUNCTION
1855 : !! For each symmetry operation, find the number of the position to
1856 : !! which each atom is sent in the unit cell by the INVERSE of the
1857 : !! symmetry operation inv(symrel); i.e. this is the atom which, when acted
1858 : !! upon by the given symmetry element isym, gets transformed into atom iatom.
1859 : !!
1860 : !! This routine uses the fact that inv(symrel)=trans(symrec),
1861 : !! the inverse of the symmetry operation expressed in the basis of real
1862 : !! space primitive translations equals the transpose of the same symmetry
1863 : !! operation expressed in the basis of reciprocal space primitive transl:
1864 : !!
1865 : !! $ xred(nu,indsym(4,isym,ia)) = symrec(mu,nu,isym)*(xred(mu,ia)-tnons(mu,isym)) - transl(mu)$
1866 : !!
1867 : !! where $transl$ is also a set of integers and
1868 : !! where translation transl places coordinates within unit cell (note sign).
1869 : !! Note that symrec is the set of arrays which are actually input here.
1870 : !! These arrays have integer elements.
1871 : !! tnons is the nonsymmorphic translation or else is zero.
1872 : !! If nsym=1 (i.e. only the identity symmetry is present) then
1873 : !! indsym merely takes each atom into itself.
1874 : !! The array of integer translations "transl" gets included within array "indsym" as seen below.
1875 : !! This routine has been improved using ideas of p. 649 of notes,
1876 : !! implementing suggestion of Andrew Horsfield: replace search for
1877 : !! equivalent atoms using direct primitive cell translations by
1878 : !! use of dot product relation which must produce an integer.
1879 : !! Relation:
1880 : !!
1881 : !! $[inv(S(i)) * (x(a)-tnons(i)) - x(inv(S)(i,a))] = integer$
1882 : !!
1883 : !! where S(i) is the symmetry matrix in real space, tnons=nonsymmorphic translation
1884 : !! (may be 0 0 0), and $x(inv(S)(i,a))$ is sought atom into which $x(a)$ gets
1885 : !! rotated by $inv(S)$. Integer gives primitive translation coordinates to get
1886 : !! back to original unit cell.
1887 : !! Equivalent to $S*t(b)+tnons-x(a)=another$ $integer$ for $x(b)=x(inv(S))$.
1888 : !!
1889 : !! INPUTS
1890 : !! natom=number of atoms in cell.
1891 : !! nsym=number of space group symmetries.
1892 : !! symrec(3,3,nsym)=symmetries expressed in terms of their action on
1893 : !! reciprocal space primitive translations (integer).
1894 : !! tnons(3,nsym)=nonsymmorphic translations for each symmetry (would
1895 : !! be 0 0 0 each for a symmorphic space group)
1896 : !! typat(natom)=integer identifying type of atom.
1897 : !! xred(3,natom)=reduced coordinates of atoms in terms of real space
1898 : !! primitive translations
1899 : !! tolsym=tolerance for the symmetries
1900 : !! [print_indsym]: Print indsym table to std_out if the number of atoms is smaller that print_indsym
1901 : !! Default: -1 i.e. no output is provided.
1902 : !!
1903 : !! OUTPUT
1904 : !! indsym(4,nsym,natom)=indirect indexing array described above: for each
1905 : !! isym,iatom, fourth element is label of atom into
1906 : !! which iatom is sent by INVERSE of symmetry operation
1907 : !! isym; first three elements are the primitive
1908 : !! translations which must be subtracted after the
1909 : !! transformation to get back to the original unit cell.
1910 : !!
1911 : !! SOURCE
1912 :
1913 50679 : subroutine symatm(indsym, natom, nsym, symrec, tnons, tolsym, typat, xred, print_indsym)
1914 :
1915 : !Arguments ------------------------------------
1916 : !scalars
1917 : integer,intent(in) :: natom,nsym
1918 : integer,optional,intent(in) :: print_indsym
1919 : real(dp), intent(in) :: tolsym
1920 : !arrays
1921 : integer,intent(in) :: symrec(3,3,nsym),typat(natom)
1922 : integer,intent(out) :: indsym(4,nsym,natom)
1923 : real(dp),intent(in) :: tnons(3,nsym),xred(3,natom)
1924 :
1925 : !Local variables-------------------------------
1926 : !scalars
1927 : integer :: eatom,errout,iatom,ii,isym,mu,print_indsym_
1928 : real(dp) :: difmax,err
1929 : character(len=500) :: msg
1930 : !arrays
1931 : integer :: transl(3)
1932 : real(dp) :: difmin(3),tratom(3)
1933 : ! *************************************************************************
1934 :
1935 : !write(std_out,'(a,i4,es12.4)')' symatm : enter, nsym,tolsym=',nsym,tolsym
1936 : !write(std_out,'(a,es12.4)')' symatm : xred='
1937 : !do ii=1,natom
1938 : ! write(std_out,'(i4,3es18.10)')ii,xred(1:3,ii)
1939 : !enddo
1940 : !write(std_out,'(a,es12.4)')' symatm : isym,symrec,tnons='
1941 : !do isym=1,nsym
1942 : ! write(std_out,'(i6,9i4,3es18.10)')isym,symrec(:,:,isym),tnons(1:3,isym)
1943 : !enddo
1944 :
1945 50679 : err=zero
1946 50679 : errout=0
1947 :
1948 1224837 : do isym=1,nsym
1949 3828118 : do iatom=1,natom
1950 :
1951 : ! Apply inverse transformation to original coordinates.
1952 : ! Note TRANSPOSE of symrec that corresponds to symrel^{-1}
1953 :
1954 10413124 : do mu=1,3
1955 : tratom(mu) = dble(symrec(1,mu,isym))*(xred(1,iatom)-tnons(1,isym))&
1956 : & +dble(symrec(2,mu,isym))*(xred(2,iatom)-tnons(2,isym))&
1957 10413124 : & +dble(symrec(3,mu,isym))*(xred(3,iatom)-tnons(3,isym))
1958 : end do
1959 :
1960 : ! Find symmetrically equivalent atom
1961 2603281 : call symchk(difmin,eatom,natom,tratom,transl,typat(iatom),typat,xred)
1962 :
1963 : ! Put information into array indsym: translations and label
1964 2603281 : indsym(1,isym,iatom)=transl(1)
1965 2603281 : indsym(2,isym,iatom)=transl(2)
1966 2603281 : indsym(3,isym,iatom)=transl(3)
1967 2603281 : indsym(4,isym,iatom)=eatom
1968 :
1969 : ! Keep track of maximum difference between transformed coordinates and
1970 : ! nearest "target" coordinate
1971 2603281 : difmax=max(abs(difmin(1)),abs(difmin(2)),abs(difmin(3)))
1972 2603281 : err=max(err,difmax)
1973 :
1974 2603281 : if(errout==3)then
1975 : write(msg, '(a)' )&
1976 1 : ' Suppress warning about finding symmetrically equivalent atoms, as mentioned already three times.'
1977 1 : ABI_WARNING(msg)
1978 1 : errout=errout+1
1979 : endif
1980 :
1981 6380720 : if (difmax>tolsym .and. errout<3) then ! Print warnings if differences exceed tolerance
1982 : write(msg, '(3a,i3,a,i6,a,i3,a,a,3f18.12,3a,es12.4)' )&
1983 3 : ' Trouble finding symmetrically equivalent atoms',ch10,&
1984 3 : ' Applying inv of symm number',isym,' to atom number',iatom,' of typat',typat(iatom),ch10,&
1985 3 : ' gives tratom=',tratom(1:3),'.',ch10,&
1986 6 : ' This is further away from every atom in crystal than the allowed tolerance, tolsym=',tolsym
1987 3 : ABI_WARNING(msg)
1988 :
1989 : write(msg, '(a,3i3,a,a,3i3,a,a,3i3)' ) &
1990 3 : ' The inverse symmetry matrix is',symrec(1,1:3,isym),ch10,&
1991 3 : ' ',symrec(2,1:3,isym),ch10,&
1992 6 : ' ',symrec(3,1:3,isym)
1993 3 : call wrtout(std_out,msg)
1994 3 : write(msg, '(a,3f18.12)' )' and the nonsymmorphic transl. tnons =',(tnons(mu,isym),mu=1,3)
1995 :
1996 3 : call wrtout(std_out,msg)
1997 : write(msg, '(a,1p,3es12.4,a,a,i5)' ) &
1998 3 : ' The nearest coordinate differs by',difmin(1:3),ch10,&
1999 6 : ' for indsym(nearest atom)=',indsym(4,isym,iatom)
2000 3 : call wrtout(std_out,msg)
2001 :
2002 : ! Use errout to reduce volume of error diagnostic output
2003 3 : if (errout==0) then
2004 1 : write(msg,'(6a)') ch10,&
2005 1 : ' This indicates that when symatm attempts to find atoms symmetrically',ch10, &
2006 1 : ' related to a given atom, the nearest candidate is further away than some',ch10,&
2007 2 : ' tolerance. Should check atomic coordinates and symmetry group input data.'
2008 1 : call wrtout(std_out,msg)
2009 : end if
2010 3 : errout=errout+1
2011 :
2012 : end if !difmax>tol
2013 : end do !iatom
2014 : end do !isym
2015 :
2016 : ! MG: Do not change this behaviour. symatm is called many times in the EPH code in which we have tons of q-points
2017 : ! and it's really annoying to see this output repeated over and over again.
2018 : ! If you need to print the indsym table at the beginning of the calculation, find the call to symatm
2019 : ! and pass the optional argument print_indsym_ or use `abitk crystal_print FILE --prtvol 1`
2020 50679 : print_indsym_ = -1; if (present(print_indsym)) print_indsym_ = print_indsym
2021 50679 : if (natom <= print_indsym_) then
2022 29525 : do iatom=1,natom
2023 20084 : write(msg, '(a,i0,a)' )' symatm: atom number ',iatom,' is reached starting at atom'
2024 20084 : call wrtout(std_out,msg)
2025 54562 : do ii=1,(nsym-1)/24+1
2026 25037 : if(natom<100)then
2027 25037 : write(msg, '(1x,24i3)' ) (indsym(4,isym,iatom),isym=1+(ii-1)*24,min(nsym,ii*24))
2028 : else
2029 0 : write(msg, '(1x,24i6)' ) (indsym(4,isym,iatom),isym=1+(ii-1)*24,min(nsym,ii*24))
2030 : end if
2031 45121 : call wrtout(std_out,msg)
2032 : end do
2033 : end do
2034 : end if
2035 :
2036 50679 : if (err>tolsym) then
2037 1 : write(msg, '(1x,a,1p,e14.5,a,e12.4)' )'symatm: maximum (delta t)=',err,' is larger than tol=',tolsym
2038 1 : ABI_WARNING(msg)
2039 : end if
2040 :
2041 : ! Stop execution if error is really big
2042 50679 : if (err>0.01d0) then
2043 : write(msg,'(5a)')&
2044 0 : 'Largest error (above) is so large (0.01) that either input atomic coordinates (xred)',ch10,&
2045 0 : 'are wrong or space group symmetry data is wrong.',ch10,&
2046 0 : 'Action: correct your input file.'
2047 0 : ABI_ERROR(msg)
2048 : end if
2049 :
2050 50679 : end subroutine symatm
2051 : !!***
2052 :
2053 : !!****f* m_symtk/symcharac
2054 : !! NAME
2055 : !! symcharac
2056 : !!
2057 : !! FUNCTION
2058 : !! Get the type of axis for the symmetry.
2059 : !!
2060 : !! INPUTS
2061 : !! center=bravais(2)
2062 : !! determinant=the value of the determinant of sym
2063 : !! iholohedry=bravais(1)
2064 : !! isym=number of the symmetry operation that is currently analyzed
2065 : !! order=the order of the symmetry
2066 : !! symrel(3,3)= the symmetry matrix
2067 : !! tnons(3)=nonsymmorphic translations
2068 : !!
2069 : !! OUTPUT
2070 : !! label=a human readable text for the characteristic of the symmetry
2071 : !! type_axis=an identifier for the type of symmetry
2072 : !!
2073 : !! SOURCE
2074 :
2075 853537 : subroutine symcharac(center, determinant, iholohedry, isym, label, symrel, tnons, type_axis)
2076 :
2077 : !Arguments ------------------------------------
2078 : !scalars
2079 : integer, intent(in) :: determinant, center, iholohedry, isym
2080 : integer, intent(out) :: type_axis
2081 : character(len=128), intent(out) :: label
2082 : !arrays
2083 : integer,intent(in) :: symrel(3,3)
2084 : real(dp),intent(in) :: tnons(3)
2085 :
2086 : !Local variables-------------------------------
2087 : !scalars
2088 : logical,parameter :: verbose=.FALSE.
2089 : integer :: tnons_order, identified, ii, order, iorder
2090 : character(len=500) :: msg
2091 : !arrays
2092 : integer :: identity(3,3),matrix(3,3),trial(3,3)
2093 : real(dp) :: reduced(3),trialt(3)
2094 : !**************************************************************************
2095 :
2096 853537 : identity(:,:)=0
2097 853537 : identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
2098 853537 : trial(:,:)=identity(:,:)
2099 853537 : matrix(:,:)=symrel(:,:)
2100 :
2101 853537 : order=0
2102 2350958 : do iorder=1,6
2103 94038320 : trial=matmul(matrix,trial)
2104 30562454 : if(sum((trial-identity)**2)==0)then
2105 717065 : order=iorder
2106 853537 : exit
2107 : end if
2108 21240609 : if(sum((trial+identity)**2)==0)then
2109 136472 : order=iorder
2110 136472 : exit
2111 : end if
2112 : end do
2113 :
2114 853537 : if(order==0)then
2115 0 : type_axis = -2
2116 0 : return
2117 : end if
2118 :
2119 : ! Determination of the characteristics of proper symmetries (rotations)
2120 853537 : if (determinant==1) then
2121 :
2122 : ! Determine the translation vector associated to the rotations
2123 : ! and its order: apply the symmetry operation
2124 : ! then analyse the resulting vector.
2125 431036 : identified=0
2126 431036 : trialt(:)=zero
2127 1606196 : do ii=1,order
2128 33335516 : trialt(:)=matmul(symrel(:,:),trialt(:))+tnons(:)
2129 : end do
2130 : ! Gives the associated translation, with components in the interval [-0.5,0.5] .
2131 1724144 : reduced(:)=trialt(:)-nint(trialt(:)-tol6)
2132 :
2133 1724144 : if(sum(abs(reduced(:)))<tol6)identified=1
2134 2137532 : if( (center==1 .or. center==-3) .and. sum(abs(reduced(:)-(/zero,half,half/)))<tol6 )identified=2
2135 2137597 : if( (center==2 .or. center==-3) .and. sum(abs(reduced(:)-(/half,zero,half/)))<tol6 )identified=3
2136 2137075 : if( (center==3 .or. center==-3) .and. sum(abs(reduced(:)-(/half,half,zero/)))<tol6 )identified=4
2137 1724144 : if(center==-1.and. sum(abs(reduced(:)-(/half,half,half/)))<tol6 )identified=5
2138 :
2139 : ! If the symmetry operation has not been identified, there is a problem ...
2140 428066 : if(identified==0) then
2141 0 : type_axis = -1
2142 0 : return
2143 : end if
2144 :
2145 : ! Compute the translation vector associated with one rotation
2146 1724144 : trialt(:)=trialt(:)/order
2147 1724144 : trialt(:)=trialt(:)-nint(trialt(:)-tol6)
2148 :
2149 : ! Analyse the resulting vector.
2150 : identified=0
2151 786690 : do ii=1,order
2152 3146760 : reduced(:)=ii*trialt(:)-nint(ii*trialt(:)-tol6)
2153 3146760 : if(sum(abs(reduced(:)))<tol6)identified=1
2154 3910265 : if( (center==1 .or. center==-3) .and. sum(abs(reduced(:)-(/zero,half,half/)))<tol6 )identified=2
2155 3910330 : if( (center==2 .or. center==-3) .and. sum(abs(reduced(:)-(/half,zero,half/)))<tol6 )identified=3
2156 3909810 : if( (center==3 .or. center==-3) .and. sum(abs(reduced(:)-(/half,half,zero/)))<tol6 )identified=4
2157 3146760 : if(center==-1.and. sum(abs(reduced(:)-(/half,half,half/)))<tol6 )identified=5
2158 :
2159 783720 : if(identified/=0)then
2160 431036 : tnons_order=ii
2161 431036 : exit
2162 : end if
2163 : end do ! ii
2164 :
2165 : ! Determinant (here=+1, as we are dealing with proper symmetry operations),
2166 : ! order, tnons_order and identified are enough to determine the kind of symmetry operation
2167 :
2168 25942 : select case(order)
2169 : case (1) ! point symmetry 1
2170 25942 : if(identified==1) then
2171 11640 : type_axis=8 ! 1
2172 11640 : write(label,'(a)') 'the identity'
2173 : else
2174 14302 : type_axis=7 ! t
2175 14302 : write(label,'(a)') 'a pure translation '
2176 : end if
2177 :
2178 405094 : if (verbose) then
2179 : write(msg,'(a,i3,2a)')' symspgr : the symmetry operation no. ',isym,' is ',trim(label)
2180 : call wrtout(std_out,msg)
2181 : end if
2182 :
2183 : case (2,3,4,6) ! point symmetry 2,3,4,6 - rotations
2184 431036 : call symaxes(center,iholohedry,isym,symrel,label,order,tnons_order,trialt,type_axis)
2185 : end select
2186 :
2187 422501 : else if (determinant==-1)then
2188 :
2189 : ! Now, take care of the improper symmetry operations.
2190 : ! Their treatment is relatively easy, except for the mirror planes
2191 17874 : select case(order)
2192 : case (1) ! point symmetry 1
2193 17874 : type_axis=5 ! -1
2194 17874 : write(label,'(a)') 'an inversion'
2195 : case (2) ! point symmetry 2 - planes
2196 290029 : call symplanes(center,iholohedry,isym,symrel,tnons,label,type_axis)
2197 : case (3) ! point symmetry 3
2198 118598 : type_axis=3 ! -3
2199 118598 : write(label,'(a)') 'a -3 axis '
2200 : case (4) ! point symmetry 1
2201 114160 : type_axis=2 ! -4
2202 114160 : write(label,'(a)') 'a -4 axis '
2203 : case (6) ! point symmetry 1
2204 438 : type_axis=1 ! -6
2205 422501 : write(label,'(a)') 'a -6 axis '
2206 : end select
2207 :
2208 : if (order /= 2 .and. verbose) then
2209 : write(msg,'(a,i3,2a)')' symspgr : the symmetry operation no. ',isym,' is ',trim(label)
2210 : call wrtout(std_out,msg)
2211 : end if
2212 :
2213 : end if ! determinant==1 or -1
2214 :
2215 : end subroutine symcharac
2216 : !!***
2217 :
2218 : !!****f* m_symtk/symaxes
2219 : !! NAME
2220 : !! symaxes
2221 : !!
2222 : !! FUNCTION
2223 : !! Determines the type of symmetry operation, for
2224 : !! the proper symmetries 2,2_1,3,3_1,3_2,4,4_1,4_2,4_3,6,6_1,...6_5
2225 : !!
2226 : !! INPUTS
2227 : !! center=type of bravais lattice centering
2228 : !! center=0 no centering
2229 : !! center=-1 body-centered
2230 : !! center=-3 face-centered
2231 : !! center=1 A-face centered
2232 : !! center=2 B-face centered
2233 : !! center=3 C-face centered
2234 : !! iholohedry=type of holohedry
2235 : !! iholohedry=1 triclinic 1bar
2236 : !! iholohedry=2 monoclinic 2/m
2237 : !! iholohedry=3 orthorhombic mmm
2238 : !! iholohedry=4 tetragonal 4/mmm
2239 : !! iholohedry=5 trigonal 3bar m (rhombohedral Bravais latt)
2240 : !! iholohedry=6 hexagonal 6/mmm
2241 : !! iholohedry=7 cubic m3bar m
2242 : !! isym=number of the symmetry operation that is currently analyzed
2243 : !! isymrelconv=symrel matrix for the particular operation, in conv. axes
2244 : !! ordersym=order of the symmetry operation
2245 : !! tnons_order=order of the screw translation
2246 : !! trialt(3)=screw translation associated with the symmetry operation
2247 : !! in conventional axes (all components in the range ]-1/2,1/2] )
2248 : !!
2249 : !! OUTPUT
2250 : !! label=a user friendly label for the rotation
2251 : !! type_axis=type of the symmetry operation
2252 : !!
2253 : !! NOTES
2254 : !! It is assumed that the symmetry operations will be entered in the
2255 : !! symrel tnonsconv arrays, for the CONVENTIONAL cell.
2256 : !! For proper symmetries (rotations), the
2257 : !! associated translation is determined.
2258 : !!
2259 : !! There is a subtlety with translations associated with rotations:
2260 : !! all the rotations with axis parallel to the one analysed do not all have the
2261 : !! same translation characteristics. This is clearly seen
2262 : !! in the extended Hermann-Mauguin symbols, see the internationa table for crystallography, chapter 4.
2263 : !! In the treatment that we adopt, one will distinguish
2264 : !! the cases of primitive Bravais lattices, and centered bravais lattices. In the latter case, in the present routine,
2265 : !! at the exception of the trigonal axis for the cubic system, we explicitely generate the correct ratio of different
2266 : !! translations, so that their type can be explicitely assigned,
2267 : !! without confusion. By contrast, for primitive lattices,
2268 : !! the "tnons" that has been transmitted to the present routine
2269 : !! might be one of the few possible translations vectors,
2270 : !! nearly at random. We deal with this case by the explicit
2271 : !! examination of the system classes, and the identification
2272 : !! of such a possibility. In particular:
2273 : !! (1) for the trigonal axis in the rhombohedral Bravais lattice,
2274 : !! or in the cubic system, there is an equal number of 3, 3_1,
2275 : !! and 3_2 axes parallel to each other, in a cell that
2276 : !! is primitive (as well as conventional). In this particular case,
2277 : !! in the present
2278 : !! routine, all 3, 3_1 and 3_2 axes are assigned to be 3 axes,
2279 : !! independently of the centering.
2280 : !! (2) for the 4- or 6- axes, no confusion is possible :
2281 : !! in the primitive cell, there is only one possible translation,
2282 : !! while in the centered cells, the correct ratio of translation
2283 : !! vectors will be generated
2284 : !! (3) for the binary axes, there is no problem when the cell
2285 : !! is centered, but there are problems
2286 : !! (3a) for the tP Bravais lattice, for an axis in a tertiary direction,
2287 : !! (see the description of the lattice symmetry directions
2288 : !! table 2.4.1 of the international tables for crystallography),
2289 : !! where the family of axes is made equally of 2 and 2_1 axis.
2290 : !! In this case, we attribute the binary axis to the specific class
2291 : !! of "tertiary 2-axis". We keep track of the 2 or 2_1
2292 : !! characteristics of all other binary axes
2293 : !! (3b) for the tI Bravais lattice, in all the directions,
2294 : !! there is an equal number of 2 and 2_1 axes. We distinguish
2295 : !! the primary and secondary family from the tertiary family.
2296 : !! (3c) for the hP Bravais lattice, each binary axis can present
2297 : !! no translation or be a screw axis (in the same direction).
2298 : !! For primary axes, one need the "2" and "2_1" classification,
2299 : !! while for secondary and tertiary axes, the associated
2300 : !! translation vector will have not importance.
2301 : !! However, one will need to distinguish secondary from
2302 : !! tertiary, and these from primary axes.
2303 : !! So, this is the most complicated case, for binary axes,
2304 : !! with the following sets of binary axes : "2", "2_1",
2305 : !! "secondary 2" and "tertiary 2".
2306 : !! (3d) for the hR Bravais lattice, each binary axis can present
2307 : !! no translation or be a screw axis (in the same direction).
2308 : !! There is no distinction between tertiary axes and other, so that
2309 : !! we simply assign a binary axis to "2-axis"
2310 : !! (3e) for the cP lattice, the binary axes along tertiary directions
2311 : !! can also have different translation vectors, while for the primary
2312 : !! direction, there is no such ambiguity. So, we will attribute
2313 : !! tertiary 2 axis to the "tertiary 2-axis" set (there are always 6),
2314 : !! and attribute 2 and 2_1 primary axes to the corresponding sets.
2315 : !!
2316 : !! SOURCE
2317 :
2318 405094 : subroutine symaxes(center,iholohedry,isym,isymrelconv,label,ordersym,tnons_order,trialt,type_axis)
2319 :
2320 : !Arguments ------------------------------------
2321 : !scalars
2322 : integer,intent(in) :: center,iholohedry,isym,ordersym,tnons_order
2323 : integer,intent(out) :: type_axis
2324 : character(len=128),intent(out) :: label
2325 : !arrays
2326 : integer,intent(in) :: isymrelconv(3,3)
2327 : real(dp),intent(in) :: trialt(3)
2328 :
2329 : !Local variables-------------------------------
2330 : !scalars
2331 : logical,parameter :: verbose=.FALSE.
2332 : character(len=500) :: msg
2333 : integer :: direction,directiontype
2334 : real(dp),parameter :: nzero=1.0d-6
2335 : !**************************************************************************
2336 :
2337 : !write(std_out,*)' symaxes : enter, isym=',isym
2338 : !write(std_out,*)' symaxes : iholohedry, ',iholohedry
2339 : !write(std_out,*)' symaxes : center, ',center
2340 :
2341 564666 : select case(ordersym)
2342 : case (2) ! point symmetry 2
2343 : ! Must characterize directiontype for cP, tP, tI, and hP Bravais lattices
2344 159572 : directiontype=1
2345 159572 : if( iholohedry==4 .or. iholohedry==7) then ! tP or cP Bravais lattices
2346 151876 : if(abs(isymrelconv(1,1))+ &
2347 : abs(isymrelconv(2,2))+ &
2348 : abs(isymrelconv(3,3)) ==1) directiontype=3
2349 7696 : else if(iholohedry==6)then ! hP Bravais lattice
2350 29406 : if(sum(isymrelconv(:,:))/=-1 )directiontype=2
2351 58812 : if(sum(isymrelconv(:,:))==0 .or. sum(isymrelconv(:,:))==-3 ) directiontype=3
2352 : ! directiontype=1 corresponds to a primary axis
2353 : ! directiontype=2 corresponds to a tertiary axis
2354 : ! directiontype=3 corresponds to a secondary axis
2355 : end if
2356 :
2357 : ! DEBUG
2358 : ! write(std_out,*)' directiontype=',directiontype
2359 : ! write(std_out,'(a,3i6)' )' isymrelconv(1:3)=',isymrelconv(:,1)
2360 : ! write(std_out,'(a,3i6)' )' isymrelconv(4:6)=',isymrelconv(:,2)
2361 : ! write(std_out,'(a,3i6)' )' isymrelconv(7:9)=',isymrelconv(:,3)
2362 : ! write(std_out,'(a,i)' )' tnons_order=',tnons_order
2363 : ! ENDDEBUG
2364 :
2365 : ! Now, classify the 2 axes
2366 1377 : if(directiontype==2)then
2367 1065 : type_axis=4 ! secondary 2 (only in the hP Bravais latt case)
2368 1065 : write(label,'(a)') 'a secondary 2-axis '
2369 :
2370 158507 : else if(directiontype==3 .and. iholohedry==4)then
2371 2672 : type_axis=21 ! tertiary 2
2372 2672 : write(label,'(a)') 'a tertiary 2-axis '
2373 89199 : else if(directiontype==3 .and. center==0 .and. (iholohedry==6.or.iholohedry==7) )then
2374 6447 : type_axis=21 ! tertiary 2
2375 6447 : write(label,'(a)') 'a tertiary 2-axis '
2376 149388 : else if(tnons_order==1 .or. (iholohedry==4 .and. center==-1) .or. iholohedry==5)then
2377 80116 : type_axis=9 ! 2
2378 80116 : write(label,'(a)') 'a 2-axis '
2379 : else
2380 69272 : type_axis=20 ! 2_1
2381 69272 : write(label,'(a)') 'a 2_1-axis '
2382 : end if
2383 :
2384 : case (3) ! point symmetry 3
2385 153262 : if(tnons_order==1)then
2386 99406 : type_axis=10 ! 3
2387 99406 : write(label,'(a)') 'a 3-axis '
2388 53856 : else if(iholohedry==5 .or. iholohedry==7)then
2389 : ! This is a special situation : in the same family of parallel 3-axis,
2390 : ! one will have an equal number of 3, 3_1 and 3_2 axes, so that
2391 : ! it is non-sense to try to classify one of them.
2392 53498 : type_axis=10 ! 3, 3_1 or 3_2, undistinguishable
2393 53498 : write(label,'(a)') 'a 3, 3_1 or 3_2 axis '
2394 : else
2395 : ! write(std_out,*)'isymrelconv=',isymrelconv(:,:)
2396 : ! write(std_out,*)'trialt=',trialt(:)
2397 : ! Must recognize 3_1 or 3_2
2398 358 : if(isymrelconv(1,1)==0)then ! 3+
2399 179 : if(abs(trialt(3)-third)<nzero)type_axis=22 ! 3_1
2400 179 : if(abs(trialt(3)+third)<nzero)type_axis=23 ! 3_2
2401 179 : else if(isymrelconv(1,1)==-1)then ! 3-
2402 179 : if(abs(trialt(3)-third)<nzero)type_axis=23 ! 3_2
2403 179 : if(abs(trialt(3)+third)<nzero)type_axis=22 ! 3_1
2404 : end if
2405 358 : write(label,'(a)') 'a 3_1 or 3_2-axis '
2406 : end if
2407 :
2408 : case (4) ! point symmetry 4
2409 91636 : if(tnons_order==1)then
2410 22366 : type_axis=12 ! 4
2411 22366 : write(label,'(a)') 'a 4-axis '
2412 69270 : else if(tnons_order==2)then
2413 15464 : type_axis=25 ! 4_2
2414 15464 : write(label,'(a)') 'a 4_2-axis '
2415 53806 : else if(center/=0)then
2416 53692 : type_axis=24 ! 4_1 or 4_3
2417 53692 : write(label,'(a)') 'a 4_1 or 4_3-axis '
2418 : else
2419 : ! write(std_out,*)'isymrelconv=',isymrelconv(:,:)
2420 : ! write(std_out,*)'trialt=',trialt(:)
2421 : ! Must recognize 4_1 or 4_3, along the three primary directions
2422 456 : do direction=1,3
2423 456 : if(isymrelconv(direction,direction)==1)then !
2424 : if( (direction==1 .and. isymrelconv(2,3)==-1) .or. &
2425 114 : (direction==2 .and. isymrelconv(3,1)==-1) .or. &
2426 : (direction==3 .and. isymrelconv(1,2)==-1) )then ! 4+
2427 57 : if(abs(trialt(direction)-quarter)<nzero)type_axis=24 ! 4_1
2428 57 : if(abs(trialt(direction)+quarter)<nzero)type_axis=26 ! 4_3
2429 : else if( (direction==1 .and. isymrelconv(2,3)==1) .or. &
2430 57 : (direction==2 .and. isymrelconv(3,1)==1) .or. &
2431 : (direction==3 .and. isymrelconv(1,2)==1) )then ! 4-
2432 57 : if(abs(trialt(direction)-quarter)<nzero)type_axis=26 ! 4_3
2433 57 : if(abs(trialt(direction)+quarter)<nzero)type_axis=24 ! 4_1
2434 : end if
2435 : end if
2436 : end do
2437 114 : write(label,'(a)') 'a 4_1 or 4_3-axis '
2438 : end if
2439 :
2440 : case (6) ! point symmetry 6
2441 405094 : if(tnons_order==1)then
2442 162 : type_axis=14 ! 6
2443 162 : write(label,'(a)') 'a 6-axis '
2444 462 : else if(tnons_order==2)then
2445 374 : type_axis=29 ! 6_3
2446 374 : write(label,'(a)') 'a 6_3-axis '
2447 88 : else if(tnons_order==3)then
2448 : !write(std_out,*)'isymrelconv=',isymrelconv(:,:)
2449 : !write(std_out,*)'trialt=',trialt(:)
2450 : !Must recognize 6_2 or 6_4
2451 46 : if(isymrelconv(1,1)==1)then ! 6+
2452 23 : if(abs(trialt(3)-third)<nzero)type_axis=28 ! 6_2
2453 23 : if(abs(trialt(3)+third)<nzero)type_axis=30 ! 6_4
2454 23 : else if(isymrelconv(1,1)==0)then ! 6-
2455 23 : if(abs(trialt(3)-third)<nzero)type_axis=30 ! 6_4
2456 23 : if(abs(trialt(3)+third)<nzero)type_axis=28 ! 6_2
2457 : end if
2458 46 : write(label,'(a)') 'a 6_2 or 6_4-axis '
2459 : else
2460 : !write(std_out,*)'isymrelconv=',isymrelconv(:,:)
2461 : !write(std_out,*)'trialt=',trialt(:)
2462 : !Must recognize 6_1 or 6_5
2463 42 : if(isymrelconv(1,1)==1)then ! 6+
2464 21 : if(abs(trialt(3)-sixth)<nzero)type_axis=27 ! 6_1
2465 21 : if(abs(trialt(3)+sixth)<nzero)type_axis=31 ! 6_5
2466 21 : else if(isymrelconv(1,1)==0)then ! 6-
2467 21 : if(abs(trialt(3)-sixth)<nzero)type_axis=31 ! 6_5
2468 21 : if(abs(trialt(3)+sixth)<nzero)type_axis=27 ! 6_1
2469 : end if
2470 42 : write(label,'(a)') 'a 6_1 or 6_5-axis '
2471 : end if
2472 :
2473 : end select
2474 :
2475 : if (verbose) then
2476 : write(msg,'(a,i3,a,a)')' symaxes : the symmetry operation no. ',isym,' is ', trim(label)
2477 : call wrtout(std_out,msg)
2478 : end if
2479 :
2480 405094 : end subroutine symaxes
2481 : !!***
2482 :
2483 : !!****f* m_symtk/symplanes
2484 : !! NAME
2485 : !! symplanes
2486 : !!
2487 : !! FUNCTION
2488 : !! Determines the type of symmetry mirror planes: m,a,b,c,d,n,g.
2489 : !! This is used (see symlist.f) to identify the space group.
2490 : !!
2491 : !! INPUTS
2492 : !! center=type of bravais lattice centering
2493 : !! center=0 no centering
2494 : !! center=-1 body-centered
2495 : !! center=-3 face-centered
2496 : !! center=1 A-face centered
2497 : !! center=2 B-face centered
2498 : !! center=3 C-face centered
2499 : !! iholohedry=type of holohedry
2500 : !! iholohedry=1 triclinic 1bar
2501 : !! iholohedry=2 monoclinic 2/m
2502 : !! iholohedry=3 orthorhombic mmm
2503 : !! iholohedry=4 tetragonal 4/mmm
2504 : !! iholohedry=5 trigonal 3bar m
2505 : !! iholohedry=6 hexagonal 6/mmm
2506 : !! iholohedry=7 cubic m3bar m
2507 : !! isym=number of the symmetry operation that is currently analyzed
2508 : !! isymrelconv=symrel matrix for the particular operation, in conv. coord.
2509 : !! itnonsconv=tnons vector for the particular operation, in conv. coord
2510 : !!
2511 : !! OUTPUT
2512 : !! label=user friendly label of the plane
2513 : !! type_axis=type of the symmetry operation
2514 : !!
2515 : !! NOTES
2516 : !! One follows the
2517 : !! conventions explained in table 1.3 of the international tables for
2518 : !! crystallography. In the case of the rhombohedral system,
2519 : !! one takes into account the first footnote of this table 1.3 .
2520 : !! In general, we will assign the different symmetries to
2521 : !! the following numbers : m -> 15 , (a, b or c) -> 16,
2522 : !! d -> 17, n -> 18 , g -> 19
2523 : !! However, there is the same problem as for binary axes,
2524 : !! namely, for parallel mirror planes, one can find different
2525 : !! translation vectors, and these might be found at random,
2526 : !! depending on the input tnons.
2527 : !! (1) In the tP case, one will distinguish tertiary
2528 : !! mirror plane, for which it is important to know whether they are
2529 : !! m or c (for tertiary planes in tP, g is equivalent to m and n is equivalent to c).
2530 : !! On the other hand, it is important to distinguish among
2531 : !! primary and secondary mirror planes, those that are m,(a or b),c, or n.
2532 : !! To summarize, the number of the symmetry will be :
2533 : !! m (primary, secondary or tertiary) -> 15 ,
2534 : !! secondary (a or b) -> 16, secondary c -> 17,
2535 : !! primary or secondary n -> 18 , tertiary c -> 19
2536 : !! (2) In the tI case, one will distinguish tertiary
2537 : !! mirror plane, for which it is important to know whether they are
2538 : !! m or d (for tertiary planes in tI, c is equivalent to m.
2539 : !! On the other hand, it is important to distinguish among
2540 : !! primary and secondary mirror planes, those that are m (equivalent to n),
2541 : !! or a,b or c.
2542 : !! To summarize, the number of the symmetry will be :
2543 : !! m (primary, secondary, tertiary) -> 15 ,
2544 : !! a,b or c (primary or secondary) -> 16, tertiary d -> 17
2545 : !! (3) For hP and hR, a m plane is always coupled to a a or b plane,
2546 : !! while a c plane is always coupled to an n plane. On the other
2547 : !! hand, it is important to distinguish between primary or secondary
2548 : !! mirror planes, and tertiary mirror planes. So we will keep the
2549 : !! following sets : m non-tertiary (that includes a or b non-tertiary) -> 15,
2550 : !! c non-tertiary (that includes n non-tertiary) -> 16,
2551 : !! m tertiary (that includes a or b non-tertiary) -> 17,
2552 : !! c tertiary (that includes n non-tertiary) -> 18.
2553 : !! For hR, all mirror planes are secondary.
2554 : !! (4) For the cP lattice, in the same spirit, one can see that
2555 : !! the tertiary m and g mirror planes are to be classified as "m" -> 15,
2556 : !! while n, a and c are to be classified as "n" -> 18. There is no need
2557 : !! to distinguish between primary, secondary or tertiary axes.
2558 : !!
2559 : !! SOURCE
2560 :
2561 171431 : subroutine symplanes(center,iholohedry,isym,isymrelconv,itnonsconv,label,type_axis)
2562 :
2563 : !Arguments ------------------------------------
2564 : !scalars
2565 : integer,intent(in) :: center,iholohedry,isym
2566 : integer,intent(out) :: type_axis
2567 : character(len = 128), intent(out) :: label
2568 : !arrays
2569 : integer,intent(in) :: isymrelconv(3,3)
2570 : real(dp),intent(in) :: itnonsconv(3)
2571 :
2572 : !Local variables-------------------------------
2573 : !scalars
2574 : logical,parameter :: verbose=.FALSE.
2575 : character(len=500) :: msg
2576 : integer :: directiontype,sum_elements
2577 : real(dp),parameter :: nzero=1.0d-6
2578 : !arrays
2579 : integer :: identity(3,3),mirrormxy(3,3),mirrormyz(3,3),mirrormzx(3,3)
2580 : integer :: mirrorx(3,3),mirrorxy(3,3),mirrory(3,3),mirroryz(3,3),mirrorz(3,3)
2581 : integer :: mirrorzx(3,3)
2582 : real(dp) :: trialt(3)
2583 : ! real(dp) :: itnonsconv2(3),trialt2(3)
2584 : !**************************************************************************
2585 :
2586 : !write(std_out,*)' symplanes : enter'
2587 : !write(std_out,*)' center,iholohedry,isym,isymrelconv,itnonsconv=',center,iholohedry,isym,isymrelconv,itnonsconv
2588 :
2589 171431 : identity(:,:)=0
2590 171431 : identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
2591 :
2592 : !Will be a mirror plane, but one must characterize
2593 : !(1) the type of plane (primary, secondary or tertiary)
2594 : !(2) the gliding vector. One now defines a few matrices.
2595 171431 : mirrorx(:,:)=identity(:,:) ; mirrorx(1,1)=-1
2596 171431 : mirrory(:,:)=identity(:,:) ; mirrory(2,2)=-1
2597 171431 : mirrorz(:,:)=identity(:,:) ; mirrorz(3,3)=-1
2598 171431 : mirrorxy(:,:)=0 ; mirrorxy(1,2)=1 ; mirrorxy(2,1)=1 ; mirrorxy(3,3)=1
2599 171431 : mirrorzx(:,:)=0 ; mirrorzx(1,3)=1 ; mirrorzx(3,1)=1 ; mirrorzx(2,2)=1
2600 171431 : mirroryz(:,:)=0 ; mirroryz(2,3)=1 ; mirroryz(3,2)=1 ; mirroryz(1,1)=1
2601 171431 : mirrormxy(:,:)=0 ; mirrormxy(1,2)=-1 ; mirrormxy(2,1)=-1 ; mirrormxy(3,3)=1
2602 171431 : mirrormzx(:,:)=0 ; mirrormzx(1,3)=-1 ; mirrormzx(3,1)=-1 ; mirrormzx(2,2)=1
2603 171431 : mirrormyz(:,:)=0 ; mirrormyz(2,3)=-1 ; mirrormyz(3,2)=-1 ; mirrormyz(1,1)=1
2604 :
2605 : !Determine the type of plane. At the end,
2606 : !directiontype=1 will correspond to a primary axis (or equivalent
2607 : !axes for orthorhombic)
2608 : !directiontype=2 will correspond to a secondary axis
2609 : !directiontype=3 will correspond to a tertiary axis
2610 : !See table 2.4.1, 11.2 and 11.3 of the international tables for crystallography
2611 171431 : directiontype=0
2612 : !The sum of elements of the matrices allow to characterize them
2613 2228603 : sum_elements=sum(isymrelconv(:,:))
2614 :
2615 171431 : if(sum_elements==1)then
2616 : ! The mirror plane perpendicular to the c axis is always primary
2617 688194 : if( sum(abs(isymrelconv(:,:)-mirrorz(:,:)))==0 )then
2618 : directiontype=1
2619 : ! All the other planes with a symrel matrix whose sum of elements is 1
2620 : ! are a or b planes. They are primary or
2621 : ! secondary planes, depending the holohedry.
2622 462839 : else if(sum(isymrelconv(:,:))==1)then
2623 35603 : if( iholohedry==2 .or. iholohedry==3 .or. iholohedry==7 )then
2624 : directiontype=1
2625 : else if(iholohedry==4 .or. iholohedry==6)then
2626 171431 : directiontype=2
2627 : end if
2628 : end if
2629 : end if
2630 :
2631 : !All the planes with a symrel matrix whose sum of elements
2632 : !is 2 are secondary planes (table 11.3).
2633 171431 : if( sum_elements==2 ) directiontype=2
2634 :
2635 : !The planes with a symrel matrix whose sum of elements
2636 : !is 3 or 0 are tertiary planes
2637 171431 : if( sum_elements==3 .or. sum_elements==0 )directiontype=3
2638 :
2639 : !One is left with sum_elements=-1, tertiary for tetragonal
2640 : !or cubic, secondary for hexagonal
2641 171431 : if( sum_elements==-1)then
2642 57549 : if(iholohedry==4 .or. iholohedry==7)directiontype=3
2643 57549 : if(iholohedry==6)directiontype=2
2644 : end if
2645 :
2646 :
2647 : !Now, determine the gliding vector
2648 : !First, apply the symmetry operation
2649 : !to itnonsconv, in order to get the translation vector
2650 : !under the application of twice the symmetry operation
2651 4800068 : trialt(:)=matmul(isymrelconv(:,:),itnonsconv(:)) +itnonsconv(:)
2652 : !Get the translation associated with one application,
2653 : !and force its components to be in the interval ]-0.5,0.5] .
2654 685724 : trialt(:)=trialt(:)*half
2655 685724 : trialt(:)=trialt(:)-nint(trialt(:)-nzero)
2656 :
2657 : !If there is a glide vector for the initial choice of itnonsconv,
2658 : !it might be that it disappears if itnonsconv is translated by a
2659 : !lattice vector of the conventional cell
2660 : !if(trialt(1)**2+trialt(2)**2+trialt(3)**2>tol5)then
2661 : !do ii=1,3
2662 : !itnonsconv2(:)=itnonsconv(:)
2663 : !itnonsconv2(ii)=itnonsconv(ii)+one
2664 : !trialt2(:)=matmul(isymrelconv(:,:),itnonsconv2(:)) +itnonsconv2(:)
2665 : !trialt2(:)=trialt2(:)*half
2666 : !trialt2(:)=trialt2(:)-nint(trialt2(:)-nzero)
2667 : !if(trialt2(1)**2+trialt2(2)**2+trialt2(3)**2<tol5)then
2668 : !trialt(:)=trialt2(:)
2669 : !endif
2670 : !enddo
2671 : !endif
2672 :
2673 171431 : write(msg,'(a)') ' symplanes...'
2674 :
2675 : !Must use the convention of table 1.3 of the international
2676 : !tables for crystallography, see also pp 788 and 789.
2677 : !Often, one needs to specialize the selection according
2678 : !to the Bravais lattice or the system.
2679 :
2680 685724 : if(sum(abs(trialt(:)))<nzero .and. iholohedry/=6)then
2681 61705 : type_axis=15 ! m
2682 61705 : write(label,'(a)') 'a mirror plane'
2683 109726 : else if(iholohedry==4 .and. center==0)then ! primitive tetragonal
2684 :
2685 931 : if(directiontype==1)then
2686 121 : type_axis=18 ! primary n
2687 121 : write(label,'(a)') 'a primary n plane'
2688 810 : else if(directiontype==2)then
2689 3332 : if(sum(abs(trialt(:)-(/half,zero,zero/)))<nzero .or. sum(abs(trialt(:)-(/zero,half,zero/)))<nzero)then
2690 130 : type_axis=16 ! secondary a or b
2691 130 : write(label,'(a)') 'a secondary a or b plane'
2692 1384 : else if(sum(abs(trialt(:)-(/zero,zero,half/)))<nzero)then
2693 130 : type_axis=17 ! secondary c
2694 130 : write(label,'(a)') 'a secondary c plane'
2695 : else
2696 216 : type_axis=18 ! secondary n
2697 216 : write(label,'(a)') 'a secondary n plane'
2698 : end if ! directiontype==2
2699 334 : else if(directiontype==3)then
2700 334 : if( abs(trialt(3))<nzero )then
2701 74 : type_axis=15 ! tertiary m
2702 74 : write(label,'(a)') 'a tertiary m plane'
2703 260 : else if( abs(trialt(3)-half)<nzero )then
2704 260 : type_axis=19 ! tertiary c
2705 260 : write(label,'(a)') 'a tertiary c plane'
2706 : end if
2707 : end if
2708 :
2709 108795 : else if(iholohedry==4 .and. center==-1)then ! inner tetragonal
2710 :
2711 1189 : if(directiontype==1 .or. directiontype==2)then
2712 : if(sum(abs(trialt(:)-(/half,zero,zero/)))<nzero .or. &
2713 6730 : sum(abs(trialt(:)-(/zero,half,zero/)))<nzero .or. &
2714 : sum(abs(trialt(:)-(/zero,zero,half/)))<nzero )then
2715 218 : type_axis=16 ! a, b, or c
2716 218 : write(label,'(a)') 'an a, b or c plane'
2717 : else if(sum(abs(trialt(:)-(/half,half,zero/)))<nzero .or. &
2718 4550 : sum(abs(trialt(:)-(/zero,half,half/)))<nzero .or. &
2719 : sum(abs(trialt(:)-(/half,zero,half/)))<nzero )then
2720 455 : type_axis=15 ! n plane, equivalent to m
2721 455 : write(label,'(a)') 'a m plane'
2722 : end if ! directiontype==1 or 2
2723 516 : else if(directiontype==3)then
2724 516 : if( abs(trialt(3))<nzero .or. abs(trialt(3)-half)<nzero )then
2725 304 : type_axis=15 ! tertiary c, equivalent to m
2726 304 : write(label,'(a)') 'a tertiary m plane'
2727 : else
2728 212 : type_axis=17 ! tertiary d
2729 212 : write(label,'(a)') 'a tertiary d plane'
2730 : end if
2731 : end if
2732 :
2733 107606 : else if(iholohedry==5)then ! hR
2734 :
2735 484 : if( abs(sum(abs(trialt(:)))-one) < nzero) then
2736 0 : type_axis=15 ! secondary m
2737 0 : write(label,'(a)') 'a secondary m plane'
2738 847 : else if( abs(sum(abs(trialt(:)))-half) < nzero .or. abs(sum(abs(trialt(:)))-three*half) < nzero )then
2739 121 : type_axis=16 ! secondary c
2740 121 : write(label,'(a)') 'a secondary c plane'
2741 : end if
2742 :
2743 107485 : else if(iholohedry==6)then ! hP
2744 :
2745 2235 : if(directiontype==1)then
2746 219 : if( abs(trialt(3)) <nzero )then
2747 219 : type_axis=15 ! primary m
2748 219 : write(label,'(a)') 'a primary m plane'
2749 : end if
2750 2016 : else if(directiontype==2)then
2751 1128 : if( abs(trialt(3)) <nzero )then
2752 951 : type_axis=15 ! secondary m
2753 951 : write(label,'(a)') 'a secondary m plane'
2754 177 : else if( abs(trialt(3)-half) < nzero ) then
2755 177 : type_axis=16 ! secondary c
2756 177 : write(label,'(a)') 'a secondary c plane'
2757 : end if
2758 888 : else if(directiontype==3)then
2759 888 : if( abs(trialt(3)) <nzero )then
2760 291 : type_axis=17 ! tertiary m
2761 291 : write(label,'(a)') 'a tertiary m plane'
2762 597 : else if( abs(trialt(3)-half) < nzero ) then
2763 597 : type_axis=18 ! tertiary c
2764 597 : write(label,'(a)') 'a tertiary c plane'
2765 : end if
2766 : end if ! directiontype
2767 :
2768 : ! else if(iholohedry==7 .and. center==0)then ! cP
2769 105250 : else if(iholohedry==7)then ! cP
2770 :
2771 102965 : if(directiontype==1)then
2772 : if((sum(abs(isymrelconv(:,:)-mirrorx(:,:)))==0 .and. &
2773 : sum(abs(two*abs(trialt(:))-(/zero,half,half/)))<nzero ).or. &
2774 : (sum(abs(isymrelconv(:,:)-mirrory(:,:)))==0 .and. &
2775 1709268 : sum(abs(two*abs(trialt(:))-(/half,zero,half/)))<nzero ).or. &
2776 : (sum(abs(isymrelconv(:,:)-mirrorz(:,:)))==0 .and. &
2777 : sum(abs(two*abs(trialt(:))-(/half,half,zero/)))<nzero ) ) then
2778 26544 : type_axis=17 ! d
2779 26544 : write(label,'(a)') 'a d plane'
2780 : else
2781 10614 : type_axis=18 ! primary n
2782 10614 : write(label,'(a)') 'a primary n plane'
2783 : end if
2784 65807 : else if(directiontype==3)then
2785 263228 : if(sum(abs(two*abs(trialt(:))-(/half,half,half/)))<nzero )then
2786 72 : type_axis=17 ! d
2787 72 : write(label,'(a)') 'a d plane'
2788 460145 : else if( abs(sum(abs(trialt(:)))-half) < nzero .or. abs(sum(abs(trialt(:)))-three*half) < nzero ) then
2789 1998 : type_axis=18 ! tertiary n
2790 1998 : write(label,'(a)') 'a tertiary n plane'
2791 254948 : else if( abs(sum(abs(trialt(:)))-one) < nzero )then
2792 63737 : type_axis=15 ! tertiary m
2793 63737 : write(label,'(a)') 'a tertiary m plane'
2794 : end if
2795 : end if
2796 :
2797 : ! Now, treat all other cases (including other centered Bravais lattices)
2798 : else if(sum(abs(trialt(:)-(/half,zero,zero/)))<nzero .or. &
2799 22850 : sum(abs(trialt(:)-(/zero,half,zero/)))<nzero .or. &
2800 : sum(abs(trialt(:)-(/zero,zero,half/)))<nzero )then
2801 1417 : type_axis=16 ! a, b or c
2802 1417 : write(label,'(a)') 'an a,b, or c plane'
2803 8680 : else if( (directiontype==1 .or. directiontype==2) .and. &
2804 : (sum(abs(trialt(:)-(/half,half,zero/)))<nzero .or. &
2805 : sum(abs(trialt(:)-(/zero,half,half/)))<nzero .or. &
2806 : sum(abs(trialt(:)-(/half,zero,half/)))<nzero ) )then
2807 776 : type_axis=18 ! n
2808 776 : write(label,'(a)') 'an n plane'
2809 368 : else if( directiontype==3 .and. sum(abs(trialt(:)-(/half,half,half/)))<nzero )then
2810 0 : type_axis=18 ! n
2811 0 : write(label,'(a)') 'an n plane'
2812 : else if((sum(abs(isymrelconv(:,:)-mirrorx(:,:)))==0 .and. &
2813 : sum(abs(two*abs(trialt(:))-(/zero,half,half/)))<nzero ).or. &
2814 : (sum(abs(isymrelconv(:,:)-mirrory(:,:)))==0 .and. &
2815 4232 : sum(abs(two*abs(trialt(:))-(/half,zero,half/)))<nzero ).or. &
2816 : (sum(abs(isymrelconv(:,:)-mirrorz(:,:)))==0 .and. &
2817 : sum(abs(two*abs(trialt(:))-(/half,half,zero/)))<nzero ) ) then
2818 92 : type_axis=17 ! d
2819 92 : write(label,'(a)') 'a d plane'
2820 0 : else if( directiontype==3 .and. sum(abs(two*abs(trialt(:))-(/half,half,half/)))<nzero)then
2821 0 : type_axis=17 ! d
2822 0 : write(label,'(a)') 'a d plane'
2823 : else
2824 0 : type_axis=19 ! g (all other planes with
2825 : ! unconventional glide vector)
2826 0 : write(label,'(a)') 'a g plane'
2827 : end if
2828 :
2829 : if (verbose) then
2830 : write(msg,'(a,i3,a,a)')' symplanes : the symmetry operation no. ',isym,' is ', trim(label)
2831 : call wrtout(std_out,msg)
2832 : end if
2833 :
2834 171431 : end subroutine symplanes
2835 : !!***
2836 :
2837 : !!****f* m_symtk/smallprim
2838 : !!
2839 : !! NAME
2840 : !! smallprim
2841 : !!
2842 : !! FUNCTION
2843 : !! Find the smallest possible primitive vectors for an input lattice
2844 : !! This algorithm is not as restrictive as the conditions mentioned at p.740
2845 : !! of the international tables for crystallography (1983).
2846 : !! The final vectors form a right-handed basis, while their
2847 : !! sign and ordering is chosen such as to maximize the overlap
2848 : !! with the original vectors in order.
2849 : !!
2850 : !! INPUTS
2851 : !! rprimd(3,3)=primitive vectors
2852 : !!
2853 : !! OUTPUT
2854 : !! metmin(3,3)=metric for the new (minimal) primitive vectors
2855 : !! minim(3,3)=minimal primitive translations
2856 : !!
2857 : !! NOTES
2858 : !! The routine might as well be defined without
2859 : !! metmin as argument, but it is more convenient to have it
2860 : !!
2861 : !! SOURCE
2862 :
2863 53184 : subroutine smallprim(metmin,minim,rprimd)
2864 :
2865 : !Arguments ------------------------------------
2866 : !arrays
2867 : real(dp),intent(in) :: rprimd(3,3)
2868 : real(dp),intent(out) :: metmin(3,3),minim(3,3)
2869 :
2870 : !Local variables-------------------------------
2871 : !scalars
2872 : integer :: ia,ib,ii,ilong,itrial,minimal
2873 : integer :: iiter, maxiter = 100000
2874 : real(dp) :: determinant,length2,metsum
2875 : character(len=500) :: msg
2876 : !arrays
2877 : integer :: nvecta(3),nvectb(3)
2878 : real(dp) :: rmet(3,3),scprod(3),tmpvect(3)
2879 : !**************************************************************************
2880 :
2881 : !call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2882 2127360 : rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
2883 :
2884 53184 : nvecta(1)=2 ; nvectb(1)=3
2885 53184 : nvecta(2)=1 ; nvectb(2)=3
2886 53184 : nvecta(3)=1 ; nvectb(3)=2
2887 :
2888 53184 : minim(:,:)=rprimd(:,:)
2889 53184 : metmin(:,:)=rmet(:,:)
2890 :
2891 : !DEBUG
2892 : !write(std_out,*)' smallprim : starting values, rprim '
2893 : !write(std_out,'(3f16.8)' )rprimd(:,1)
2894 : !write(std_out,'(3f16.8)' )rprimd(:,2)
2895 : !write(std_out,'(3f16.8)' )rprimd(:,3)
2896 : !write(std_out,*)' smallprim : starting values, rmet '
2897 : !write(std_out,'(3f16.8)' )rmet(:,1)
2898 : !write(std_out,'(3f16.8)' )rmet(:,2)
2899 : !write(std_out,'(3f16.8)' )rmet(:,3)
2900 : !call flush(std_out)
2901 : !ENDDEBUG
2902 :
2903 : !Note this loop without index
2904 65720 : do iiter = 1, maxiter
2905 :
2906 : ! Will exit if minimal=1 is still valid after a trial
2907 : ! to reduce the vectors of each of the three pairs
2908 : minimal=1
2909 :
2910 262880 : do itrial=1,3
2911 :
2912 197160 : ia=nvecta(itrial) ; ib=nvectb(itrial)
2913 : ! Make sure the scalar product is negative
2914 197160 : if(metmin(ia,ib)>tol8)then
2915 345836 : minim(:,ia)=-minim(:,ia)
2916 86459 : metmin(ia,ib)=-metmin(ia,ib) ; metmin(ib,ia)=-metmin(ib,ia)
2917 86459 : metmin(ia,itrial)=-metmin(ia,itrial)
2918 86459 : metmin(itrial,ia)=-metmin(itrial,ia)
2919 : end if
2920 : ! Compute the length of the sum vector
2921 197160 : length2=metmin(ia,ia)+2*metmin(ia,ib)+metmin(ib,ib)
2922 : ! Replace the first vector by the sum vector if the latter is smaller
2923 262880 : if(length2/metmin(ia,ia) < one-tol8)then
2924 73668 : minim(:,ia)=minim(:,ia)+minim(:,ib)
2925 18417 : metmin(ia,ia)=length2
2926 18417 : metmin(ia,ib)=metmin(ia,ib)+metmin(ib,ib)
2927 18417 : metmin(ia,itrial)=metmin(ia,itrial)+metmin(ib,itrial)
2928 18417 : metmin(ib,ia)=metmin(ia,ib)
2929 18417 : metmin(itrial,ia)=metmin(ia,itrial)
2930 18417 : minimal=0
2931 : ! Replace the second vector by the sum vector if the latter is smaller
2932 178743 : else if(length2/metmin(ib,ib) < one-tol8)then
2933 13088 : minim(:,ib)=minim(:,ia)+minim(:,ib)
2934 3272 : metmin(ib,ib)=length2
2935 3272 : metmin(ia,ib)=metmin(ia,ib)+metmin(ia,ia)
2936 3272 : metmin(itrial,ib)=metmin(itrial,ib)+metmin(itrial,ia)
2937 3272 : metmin(ib,ia)=metmin(ia,ib)
2938 3272 : metmin(ib,itrial)=metmin(itrial,ib)
2939 3272 : minimal=0
2940 : end if
2941 :
2942 : end do
2943 :
2944 65720 : if(minimal==1)exit
2945 : end do
2946 :
2947 53184 : if (iiter >= maxiter) then
2948 0 : write(msg,'(a,i0,a)') 'the loop has failed to find a set of minimal vectors in ',maxiter,' iterations.'
2949 0 : ABI_BUG(msg)
2950 : end if
2951 :
2952 : !DEBUG
2953 : !write(std_out,*)' smallprim : after pair optimization '
2954 : !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' minim =',ch10,minim(:,1),ch10,minim(:,2),ch10,minim(:,3)
2955 : !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' metmin =',ch10,metmin(:,1),ch10,metmin(:,2),ch10,metmin(:,3)
2956 : !write(std_out,*)' smallprim : will start triplet optimization ',ch10
2957 : !call flush(std_out)
2958 : !ENDDEBUG
2959 :
2960 : !At this stage, the three vectors have angles between each other that are
2961 : !comprised between 90 and 120 degrees. It might still be that minus the vector
2962 : !that is the sum of the three vectors is smaller than the longest of these vectors
2963 53523 : do iiter = 1, maxiter
2964 :
2965 : ! Will exit if minimal=1 is still valid after a trial
2966 : ! to replace the longest of the three vectors by one of the triplet sum of the three vectors, with plus or minus sign
2967 :
2968 : ! Find longest of the three vectors
2969 53523 : ilong=1
2970 53523 : if( metmin(2,2)/metmin(1,1) > one + tol8 )ilong=2
2971 53523 : if( metmin(3,3)/metmin(ilong,ilong) > one + tol8)ilong=3
2972 :
2973 : ! Try combination with all same signs
2974 53523 : minimal=1
2975 695799 : metsum=sum(metmin(:,:))
2976 53523 : itrial=0
2977 53523 : if( metsum/metmin(ilong,ilong) < one - tol8)then
2978 : ! Better combination indeed ...
2979 1356 : minim(:,ilong)=minim(:,1)+minim(:,2)+minim(:,3)
2980 13560 : metmin=MATMUL(TRANSPOSE(minim),minim)
2981 : minimal=0
2982 : else
2983 : ! Try combinations with sign of itrial different from others
2984 53184 : metsum=two*(metmin(1,1)+metmin(2,2)+metmin(3,3))-metsum
2985 212736 : do itrial=1,3
2986 159552 : ia=nvecta(itrial) ; ib=nvectb(itrial)
2987 212736 : if( (metsum+four*metmin(ia,ib))/metmin(ilong,ilong) < one - tol8)then
2988 : ! Better combination indeed ...
2989 0 : metsum=metsum+four*metmin(ia,ib)
2990 0 : minim(:,ilong)=-minim(:,itrial)+minim(:,ia)+minim(:,ib)
2991 0 : metmin=MATMUL(TRANSPOSE(minim),minim)
2992 : minimal=0
2993 : exit
2994 : endif
2995 : enddo
2996 : endif
2997 :
2998 : !DEBUG
2999 : !write(std_out,*)' smallprim : triplet optimization, iiter,ilong,itrial= ',iiter,ilong,itrial
3000 : !write(std_out,*)' smallprim : predict met for the new vector=',metsum
3001 : !call flush(std_out)
3002 : !ENDDEBUG
3003 :
3004 : ! do itrial=1,3
3005 : ! ia=nvecta(itrial) ; ib=nvectb(itrial)
3006 : ! if(metmin(ia,ia)/metsum > one + tol8)then
3007 : ! minim(:,ia)=-minim(:,1)-minim(:,2)-minim(:,3)
3008 : ! metmin(ia,ib)=-sum(metmin(:,ib))
3009 : ! metmin(ia,itrial)=-sum(metmin(:,itrial))
3010 : ! metmin(ia,ia)=metsum
3011 : ! metmin(ib,ia)=metmin(ia,ib)
3012 : ! metmin(itrial,ia)=metmin(ia,itrial)
3013 : ! minimal=0
3014 : ! end if
3015 : ! end do
3016 :
3017 : !DEBUG
3018 : !write(std_out,*)' smallprim : found better primitive vector using triplets, itrial= ',itrial
3019 : !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' minim =',ch10,minim(:,1),ch10,minim(:,2),ch10,minim(:,3)
3020 : !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' metmin =',ch10,metmin(:,1),ch10,metmin(:,2),ch10,metmin(:,3)
3021 : !write(std_out,*)' smallprim : will continue triplet optimization ',ch10
3022 : !call flush(std_out)
3023 : !ENDDEBUG
3024 :
3025 53523 : if(minimal==1)exit
3026 :
3027 : end do
3028 :
3029 53184 : if (iiter >= maxiter) then
3030 0 : write(msg, '(a,i0,a)') 'the second loop has failed to find a set of minimal vectors in ',maxiter, 'iterations.'
3031 0 : ABI_BUG(msg)
3032 : end if
3033 :
3034 : !DEBUG
3035 : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
3036 : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' minim =',minim(:,1),ch10,minim(:,2),ch10,minim(:,3)
3037 : !ENDDEBUG
3038 :
3039 : !DEBUG
3040 : !Change sign of the third vector if not right-handed basis
3041 : !determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+&
3042 : !& minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+&
3043 : !& minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3))
3044 : !write(std_out,*)' smallprim: determinant=',determinant
3045 : !ENDDEBUG
3046 :
3047 : !Choose the first vector
3048 : !Compute the scalar product of the three minimal vectors with the first original vector
3049 53184 : scprod(:)=zero
3050 212736 : do ii=1,3
3051 691392 : scprod(:)=scprod(:)+minim(ii,:)*rprimd(ii,1)
3052 : end do
3053 : !Determine the vector with the maximal absolute overlap
3054 53184 : itrial=1
3055 53184 : if(abs(scprod(2))>abs(scprod(1))+tol8)itrial=2
3056 53184 : if(abs(scprod(3))>abs(scprod(itrial))+tol8)itrial=3
3057 : !Switch the vectors if needed
3058 53184 : if(itrial/=1)then
3059 11492 : tmpvect(:)=minim(:,1)
3060 11492 : minim(:,1)=minim(:,itrial)
3061 11492 : minim(:,itrial)=tmpvect(:)
3062 : end if
3063 : !Choose the sign
3064 64120 : if(scprod(itrial)<tol8)minim(:,1)=-minim(:,1)
3065 :
3066 : !DEBUG
3067 : !Change sign of the third vector if not right-handed basis
3068 : !determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+&
3069 : !& minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+&
3070 : !& minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3))
3071 : !write(std_out,*)' smallprim: determinant=',determinant
3072 : !ENDDEBUG
3073 :
3074 : !Choose the second vector
3075 : !Compute the scalar product of the second and third minimal vectors with the second original vector
3076 159552 : scprod(2:3)=zero
3077 212736 : do ii=1,3
3078 531840 : scprod(2:3)=scprod(2:3)+minim(ii,2:3)*rprimd(ii,2)
3079 : end do
3080 : !Determine the vector with the maximal absolute overlap
3081 53184 : itrial=2
3082 53184 : if(abs(scprod(3))>abs(scprod(2))+tol8)itrial=3
3083 : !Switch the vectors if needed
3084 : if(itrial/=2)then
3085 4944 : tmpvect(:)=minim(:,2)
3086 4944 : minim(:,2)=minim(:,itrial)
3087 4944 : minim(:,itrial)=tmpvect(:)
3088 : end if
3089 : !Choose the sign
3090 140844 : if(scprod(itrial)<tol8)minim(:,2)=-minim(:,2)
3091 :
3092 : !Change sign of the third vector if not right-handed basis
3093 : determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+&
3094 : & minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+&
3095 53184 : & minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3))
3096 55796 : if(determinant<-tol8)minim(:,3)=-minim(:,3)
3097 53184 : if(abs(determinant)<tol8)then
3098 0 : ABI_BUG('minim gives vanishing unit cell volume.')
3099 : end if
3100 :
3101 : !Final computation of metmin
3102 212736 : do ii=1,3
3103 691392 : metmin(ii,:)=minim(1,ii)*minim(1,:)+ minim(2,ii)*minim(2,:)+ minim(3,ii)*minim(3,:)
3104 : end do
3105 :
3106 : !DEBUG
3107 : !write(std_out,'(2a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',ch10,rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
3108 : !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' minim =',ch10,minim(:,1),ch10,minim(:,2),ch10,minim(:,3)
3109 : !write(std_out,'(2a,3es16.8,a,3es16.8,a,3es16.8)')' metmin =',ch10,metmin(:,1),ch10,metmin(:,2),ch10,metmin(:,3)
3110 : !write(std_out,'(a)')' smallprim : exit '
3111 : !call flush(std_out)
3112 : !ENDDEBUG
3113 :
3114 53184 : end subroutine smallprim
3115 : !!***
3116 :
3117 : !!****f* m_symtk/print_symmetries
3118 : !! NAME
3119 : !! print_symmetries
3120 : !!
3121 : !! FUNCTION
3122 : !! Helper function to print the set of symmetries.
3123 : !!
3124 : !! SOURCE
3125 :
3126 1 : subroutine print_symmetries(units, nsym, symrel, tnons, symafm)
3127 :
3128 : !Arguments ------------------------------------
3129 : !scalars
3130 : integer,intent(in) :: nsym, units(:)
3131 : !arrays
3132 : integer,intent(in) :: symrel(3,3,nsym),symafm(nsym)
3133 : real(dp),intent(in) :: tnons(3,nsym)
3134 :
3135 : !Local variables-------------------------------
3136 : integer :: isym, isymin, isymend, ii, jj
3137 : character(len=500) :: msg
3138 : ! *********************************************************************
3139 :
3140 1 : write(msg,'(2a)')ch10,' Symmetry operations in real space (Rotation tnons AFM)'
3141 1 : call wrtout(units, msg)
3142 :
3143 3 : do isymin=1,nsym,4
3144 2 : isymend=isymin+3
3145 2 : if (isymend>nsym) isymend=nsym
3146 8 : do ii=1,3
3147 66 : write(msg,'(4(3i3,f11.6,i3,3x))')((symrel(ii,jj,isym),jj=1,3),tnons(ii,isym),symafm(isym),isym=isymin,isymend)
3148 8 : call wrtout(units,msg)
3149 : end do
3150 2 : write(msg,'(a)')ch10
3151 3 : call wrtout(units, msg)
3152 : end do
3153 :
3154 1 : end subroutine print_symmetries
3155 : !!***
3156 :
3157 : !---------------------------------------------------------------
3158 : ! Main driver: analyze a rotation given in reduced coordinates.
3159 : !
3160 : ! rred : 3x3 rotation matrix in reduced coordinates (should have
3161 : ! integer entries for a crystallographic symmetry operation
3162 : ! label : output descriptive string
3163 : ! rprim : OPTIONAL 3x3 matrix whose COLUMNS are the lattice vectors
3164 : ! a1, a2, a3 in Cartesian coordinates. If present, a
3165 : ! Cartesian axis and unit vector are also reported.
3166 : !---------------------------------------------------------------
3167 : ! Given a rotation (point-symmetry) operation expressed as a 3x3 matrix
3168 : ! in REDUCED (fractional/lattice) coordinates, this module determines:
3169 : !
3170 : ! - whether the operation is proper (det = +1) or improper (det = -1)
3171 : ! - the rotation order n and angle (in degrees)
3172 : ! - the rotation axis, expressed in reduced coordinates (small integer
3173 : ! triplet) and, optionally, in Cartesian coordinates if the lattice
3174 : ! vectors are supplied
3175 : ! - a crystallographic-style label: 1, 2, 3, 4, 6 (proper)
3176 : ! -1, m, -3, -4, -6 (improper)
3177 : !
3178 : ! Key fact used: trace and determinant are invariant under a similarity
3179 : ! transform (R_cart = A * R_red * A^-1), so det/trace/angle can be
3180 : ! obtained directly from the reduced-coordinate matrix without ever
3181 : ! forming the Cartesian matrix. Only the axis direction needs the
3182 : ! lattice vectors to be expressed in real space.
3183 :
3184 0 : subroutine rot2str(rred, label, rprim)
3185 : integer, intent(in) :: rred(3,3)
3186 : character(len=*), intent(out) :: label
3187 : real(dp), intent(in), optional :: rprim(3,3)
3188 :
3189 : real(dp), parameter :: TOL = 1.0e-4_dp
3190 : real(dp) :: trR, trP, costh, theta_deg
3191 : real(dp) :: P(3,3), axis_red(3)
3192 : integer :: detR, order
3193 : logical :: proper, is_mirror, is_inversion, is_identity
3194 : character(len=64) :: axisstr, anglestr, cartstr
3195 : character(len=8) :: ordlab
3196 :
3197 0 : call mati3det(rred, detR)
3198 0 : trR = rred(1,1) + rred(2,2) + rred(3,3)
3199 :
3200 : ! Proper part P of the operation: P = R if proper, P = -R if improper.
3201 : ! P always has det(P) = +1 and represents a pure rotation.
3202 0 : proper = (detR > zero)
3203 0 : if (proper) then
3204 0 : P = rred
3205 : else
3206 0 : P = -rred
3207 : end if
3208 0 : trP = P(1,1) + P(2,2) + P(3,3)
3209 :
3210 0 : costh = (trP - 1.0_dp) * 0.5_dp
3211 0 : costh = max(-1.0_dp, min(1.0_dp, costh))
3212 0 : theta_deg = acos(costh) * 180.0_dp / pi
3213 :
3214 0 : if (theta_deg < TOL) then
3215 0 : order = 1
3216 : else
3217 0 : order = nint(360.0_dp / theta_deg)
3218 : end if
3219 :
3220 0 : is_identity = (proper .and. order == 1)
3221 0 : is_inversion = ((.not. proper) .and. order == 1)
3222 0 : is_mirror = ((.not. proper) .and. order == 2)
3223 :
3224 : ! --- rotation axis (reduced coordinates) ---
3225 0 : if (is_identity .or. is_inversion) then
3226 0 : axis_red = zero
3227 0 : axisstr = "(none - no unique axis)"
3228 : else
3229 0 : call rotation_axis(P, axis_red)
3230 0 : write(axisstr, '(A,3(F7.3,1X),A)') "[ ", axis_red, "] (reduced coords)"
3231 : end if
3232 :
3233 : ! --- optional Cartesian axis ---
3234 0 : cartstr = ""
3235 0 : if (present(rprim) .and. .not. (is_identity .or. is_inversion)) then
3236 : block
3237 : real(dp) :: axc(3), nrm
3238 0 : axc = matmul(rprim, axis_red)
3239 0 : nrm = sqrt(sum(axc**2))
3240 0 : if (nrm > TOL) axc = axc / nrm
3241 0 : write(cartstr, '(A,3(F7.4,1X),A)') ", cart axis [ ", axc, "]"
3242 : end block
3243 : end if
3244 :
3245 : ! --- crystallographic label ---
3246 0 : write(ordlab,'(I0)') order
3247 0 : if (is_identity) then
3248 0 : label = "1 (identity, proper, angle=0.0 deg)"
3249 0 : else if (is_inversion) then
3250 0 : label = "-1 (inversion center, improper, angle=0.0 deg)"
3251 0 : else if (is_mirror) then
3252 0 : write(anglestr,'(F6.2)') theta_deg
3253 0 : label = "m (mirror plane, improper, normal "//trim(axisstr)//trim(cartstr)//")"
3254 0 : else if (proper) then
3255 0 : write(anglestr,'(F6.2)') theta_deg
3256 : label = trim(ordlab)//"-fold proper rotation, axis "//trim(axisstr)// &
3257 0 : trim(cartstr)//", angle = "//trim(adjustl(anglestr))//" deg"
3258 : else
3259 0 : write(anglestr,'(F6.2)') theta_deg
3260 : label = "-"//trim(ordlab)//" (roto-inversion, improper), axis "// &
3261 0 : trim(axisstr)//trim(cartstr)//", angle = "// trim(adjustl(anglestr))//" deg"
3262 : end if
3263 :
3264 0 : end subroutine rot2str
3265 :
3266 : !---------------------------------------------------------------
3267 : ! Extract the rotation axis of a proper rotation matrix P (det=+1)
3268 : ! by computing the adjugate of M = P - I. Since M is singular
3269 : ! (rank <= 2 for any rotation other than identity), every column
3270 : ! of adj(M) is proportional to the null vector of M, i.e. to the
3271 : ! rotation axis. We pick the column of largest norm for numerical
3272 : ! robustness, then reduce it to small integers via the GCD.
3273 : !---------------------------------------------------------------
3274 0 : subroutine rotation_axis(P, axis)
3275 : real(dp), intent(in) :: P(3,3)
3276 : real(dp), intent(out) :: axis(3)
3277 : real(dp) :: M(3,3), adj(3,3), nrm(3)
3278 : integer :: i, jbest
3279 : real(dp), parameter :: TOL = 1.0e-4_dp
3280 : real(dp) :: best
3281 :
3282 0 : M = P
3283 0 : M(1,1) = M(1,1) - 1.0_dp
3284 0 : M(2,2) = M(2,2) - 1.0_dp
3285 0 : M(3,3) = M(3,3) - 1.0_dp
3286 :
3287 0 : adj(1,1) = M(2,2)*M(3,3) - M(2,3)*M(3,2)
3288 0 : adj(1,2) = M(1,3)*M(3,2) - M(1,2)*M(3,3)
3289 0 : adj(1,3) = M(1,2)*M(2,3) - M(1,3)*M(2,2)
3290 0 : adj(2,1) = M(2,3)*M(3,1) - M(2,1)*M(3,3)
3291 0 : adj(2,2) = M(1,1)*M(3,3) - M(1,3)*M(3,1)
3292 0 : adj(2,3) = M(1,3)*M(2,1) - M(1,1)*M(2,3)
3293 0 : adj(3,1) = M(2,1)*M(3,2) - M(2,2)*M(3,1)
3294 0 : adj(3,2) = M(1,2)*M(3,1) - M(1,1)*M(3,2)
3295 0 : adj(3,3) = M(1,1)*M(2,2) - M(1,2)*M(2,1)
3296 :
3297 0 : do i = 1, 3
3298 0 : nrm(i) = sqrt(adj(1,i)**2 + adj(2,i)**2 + adj(3,i)**2)
3299 : end do
3300 :
3301 0 : jbest = maxloc(nrm, dim=1)
3302 0 : best = nrm(jbest)
3303 :
3304 0 : if (best < TOL) then
3305 : ! Degenerate fallback (M ~ 0, e.g. numerical issues): just
3306 : ! return a zero vector; caller already filters identity/inversion.
3307 0 : axis = zero
3308 0 : return
3309 : end if
3310 :
3311 0 : axis = adj(:, jbest)
3312 0 : call reduce_to_small_integers(axis)
3313 :
3314 : end subroutine rotation_axis
3315 :
3316 : !---------------------------------------------------------------
3317 : ! Rescale a (near-)integer vector by the GCD of its rounded
3318 : ! components, and fix an overall sign convention (first nonzero
3319 : ! component positive) so the axis is reported in a canonical form.
3320 : !---------------------------------------------------------------
3321 0 : subroutine reduce_to_small_integers(v)
3322 : real(dp), intent(inout) :: v(3)
3323 : integer :: iv(3), g, i
3324 :
3325 0 : iv = nint(v)
3326 0 : if (all(iv == 0)) return
3327 :
3328 0 : g = 0
3329 0 : do i = 1, 3
3330 0 : g = igcd(g, abs(iv(i)))
3331 : end do
3332 0 : if (g > 0) iv = iv / g
3333 :
3334 0 : do i = 1, 3
3335 0 : if (iv(i) /= 0) then
3336 0 : if (iv(i) < 0) iv = -iv
3337 : exit
3338 : end if
3339 : end do
3340 :
3341 0 : v = real(iv, dp)
3342 : end subroutine reduce_to_small_integers
3343 : !!***
3344 :
3345 0 : recursive function igcd(a, b) result(g)
3346 : integer, intent(in) :: a, b
3347 : integer :: g
3348 0 : if (b == 0) then
3349 0 : g = a
3350 : else
3351 0 : g = igcd(b, mod(a, b))
3352 : end if
3353 0 : end function igcd
3354 :
3355 : !!****f* m_symtk/sym_order
3356 : !! NAME
3357 : !! sym_order
3358 : !!
3359 : !! FUNCTION
3360 : !! Return the order n of a (possibly non-symmorphic) space-group operation {rot|tnons},
3361 : !! i.e. the smallest integer such that rot^n = identity. By the crystallographic
3362 : !! restriction theorem, n must be one of {1, 2, 3, 4, 6}: the routine aborts if none
3363 : !! of these values gives the identity, as this signals that rot is not a valid
3364 : !! crystallographic point-group operation.
3365 : !!
3366 : !! If tnons is given, also compute the cumulative translation T such that
3367 : !! applying the operation n times gives:
3368 : !!
3369 : !! S^n(r) = r + T, with T = [I + rot + rot^2 + ... + rot^(n-1)] . tnons
3370 : !!
3371 : !! For symmorphic operations (or when tnons is a lattice vector times a screw/glide
3372 : !! fraction that closes exactly), T reduces to a lattice vector.
3373 : !!
3374 : !! INPUTS
3375 : !! rot(3,3)=Rotation matrix in reduced coordinates (e.g. symrel(:,:,isym)).
3376 : !! tnons(3)=Optional non-symmorphic translation in reduced coordinates (e.g. tnons(:,isym)).
3377 : !!
3378 : !! OUTPUT
3379 : !! n=Order of the operation, one of {1, 2, 3, 4, 6}.
3380 : !! isproper=.True. if rot is a proper rotation (det=+1), .False. if
3381 : !! improper (det=-1, e.g. mirror, inversion, rotoinversion).
3382 : !! trans(3)=Cumulative translation T (see above).
3383 : !! msg and ierr= Error messate and exit status.
3384 : !!
3385 : !! SOURCE
3386 :
3387 0 : subroutine sym_order(rot, tnons, n, isproper, trans, msg, ierr)
3388 :
3389 : !Arguments ------------------------------------
3390 : integer,intent(in) :: rot(3,3)
3391 : real(dp),intent(in) :: tnons(3)
3392 : integer,intent(out) :: n
3393 : logical,intent(out) :: isproper
3394 : integer,intent(out) :: trans(3), ierr
3395 : character(len=*),intent(out) :: msg
3396 :
3397 : !Local variables-------------------------------
3398 : integer,parameter :: norders = 5
3399 : integer,parameter :: allowed_orders(norders) = (/1, 2, 3, 4, 6/)
3400 : integer,parameter :: identity(3,3) = reshape((/1,0,0, 0,1,0, 0,0,1/), (/3,3/))
3401 : integer :: io, k, det
3402 : integer :: rot_k(3,3)
3403 : real(dp) :: trans_dp(3)
3404 : logical :: found
3405 : ! *********************************************************************
3406 :
3407 0 : ierr = 0; msg = ""
3408 :
3409 0 : call mati3det(rot, det)
3410 0 : isproper = (det == 1)
3411 :
3412 0 : found = .False.
3413 0 : do io=1,norders
3414 0 : n = allowed_orders(io)
3415 0 : rot_k = identity
3416 0 : do k=1,n
3417 0 : rot_k = matmul(rot, rot_k)
3418 : end do
3419 0 : if (all(rot_k == identity)) then
3420 : found = .True.
3421 : exit
3422 : end if
3423 : end do
3424 :
3425 0 : if (.not. found) then
3426 0 : msg = "Rotation order does not belong to {1, 2, 3, 4, 6}: this is not a valid crystallographic point-group operation!"
3427 0 : ierr = 1; return
3428 : end if
3429 :
3430 0 : trans_dp = zero
3431 0 : rot_k = identity
3432 0 : do k=1,n
3433 0 : trans_dp = trans_dp + matmul(rot_k, tnons)
3434 0 : rot_k = matmul(rot, rot_k)
3435 : end do
3436 :
3437 : ! Make sure trans_dp is integer.
3438 0 : trans = nint(trans_dp)
3439 0 : if (.not. isinteger(trans_dp, tol=tol5)) then
3440 0 : ierr = 2
3441 0 : msg = sjoin("T = [I + rot + rot^2 + ... + rot^(n-1)] . tnons is not integer:", ltoa(trans_dp))
3442 : end if
3443 :
3444 0 : end subroutine sym_order
3445 : !!***
3446 :
3447 : end module m_symtk
3448 : !!***
|
