Branch data 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 : 33061 : 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 [ + + ]: 1310065 : do isym=1,nsym
98 : 1277004 : call mati3det(sym(:,:,isym),det)
99 : 1277004 : determinant(isym)=det
100 [ - + ]: 1310065 : 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 : 33061 : 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 : 93072 : subroutine sg_multable(nsym, symafm, symrel, ierr, &
257 : 36041 : 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 : 46536 : integer,allocatable :: ptsymm(:),list_symrel(:,:)
278 : : real(dp) :: prd_tnons(3)
279 : 46536 : real(dp),allocatable :: tnons_(:,:)
280 : : ! *************************************************************************
281 : :
282 : : !write(std_out,*)' m_symtk%sg_multable : enter, nsym= ',nsym
283 : 46536 : ierr = 0
284 : :
285 [ - + + - : 139608 : ABI_MALLOC(tnons_,(3,nsym))
+ - + - -
+ - + ]
286 : :
287 [ + + + + ]: 4362660 : tnons_=zero
288 [ + + + - : 3753633 : if(present(tnons)) tnons_=tnons
- + - - +
+ + + ]
289 : :
290 : 46536 : tnons_tol_=tol5
291 [ + + ]: 46536 : 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 [ + + + + : 604968 : 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 : 46536 : echo = 1
301 [ + + ]: 1125563 : do sym1=1,nsym
302 : 1079028 : found_inv = .FALSE.
303 [ + + ]: 30510627 : do sym2=1,nsym
304 [ + + + + : 1220425040 : prd_symrel = matmul(symrel(:,:,sym1), symrel(:,:,sym2))
+ + ]
305 [ + + + + : 854297528 : prd_tnons = tnons_(:,sym1) + matmul(symrel(:,:,sym1), tnons_(:,sym2))
+ + + + +
+ ]
306 : 30510626 : prd_symafm = symafm(sym1)*symafm(sym2)
307 [ + + + + : 66957865 : if (all(prd_symrel == identity_3d) .and. isinteger(prd_tnons, tnons_tol_) .and. prd_symafm == 1) then
+ + + + +
+ + + ]
308 : 1079027 : found_inv = .TRUE.
309 [ - + ]: 1079027 : if (present(toinv)) then
310 [ # # ]: 0 : toinv(1, sym1) = sym2; toinv(2:4, sym1) = nint(prd_tnons)
311 : : end if
312 : 1079027 : exit
313 : : end if
314 : : end do
315 : :
316 : 46535 : 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 [ - + + - : 139608 : ABI_MALLOC(list_symrel,(nsym,48))
+ - + - -
+ - + ]
334 [ - + + - : 139608 : ABI_MALLOC(ptsymm,(nsym))
+ - + - -
+ - + ]
335 : :
336 : 46536 : nlist_symrel(:)=0
337 : : ! Initialize with the first symmetry operation
338 [ + + + + ]: 604968 : ptsymrel(1:3,1:3,1)=symrel(:,:,1)
339 : 46536 : ptsymm(1)=1
340 : 46536 : nptsymm=1
341 : 46536 : list_symrel(1,1)=1
342 : 46536 : 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 [ + + ]: 46536 : if(nsym/=1)then
346 [ + + ]: 1071184 : do sym1=2,nsym
347 : : found=0
348 [ + + ]: 22454293 : do ptsymm2=1,nptsymm
349 [ + + + + : 42119556 : 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 [ + + ]: 1071184 : if(found==0)then
357 : 917042 : nptsymm=nptsymm+1
358 : : !write(std_out,*)' current value of nptsymm, sym1=',nptsymm, sym1
359 [ + + + + ]: 11921546 : ptsymrel(1:3,1:3,nptsymm)=symrel(:,:,sym1)
360 : 917042 : ptsymm(sym1)=nptsymm
361 : 917042 : nlist_symrel(nptsymm)=1
362 : 917042 : 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 [ + + ]: 38689 : if(nptsymm/=1)then
369 [ + + ]: 994398 : do ptsymm1=1,nptsymm
370 [ - + ]: 994398 : 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 : 46536 : echo = 1
392 [ + + ]: 1010110 : do ptsymm1=1,nptsymm
393 : 963575 : sym1=list_symrel(1,ptsymm1)
394 [ + + ]: 39178494 : do ptsymm2=1,nptsymm
395 : 38214920 : sym2=list_symrel(1,ptsymm2)
396 : : ! Compute the product of the two symmetries.
397 [ + + + + : 1528596800 : prd_symrel = matmul(symrel(:,:,sym1), symrel(:,:,sym2))
+ + ]
398 : : ! Check that product array is one of the original point symmetries.
399 : : iseq= .false.
400 [ + + ]: 881784757 : do ptsymm3=1,nptsymm
401 [ + + + + : 2012615863 : iseq = all(prd_symrel == symrel(:,:,list_symrel(1,ptsymm3)))
+ + ]
402 [ + + ]: 881784757 : if (iseq) then
403 : 38214919 : ptmultable(ptsymm1,ptsymm2) = ptsymm3; exit
404 : : endif
405 : : end do
406 : :
407 [ + + ]: 39178494 : 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 [ + + ]: 1010110 : 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 [ + + + - ]: 46536 : if(echo/=0 .and. ierr==0)then
444 [ + + ]: 1125561 : do sym1=1,nsym
445 : 1079026 : ptsymm1=ptsymm(sym1)
446 [ + + ]: 61919460 : do sym2=1,nsym
447 : 60840434 : 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 [ - + - - ]: 60840434 : 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 : 60840434 : prd_ptsymm=ptmultable(ptsymm1,ptsymm2)
466 [ + + + + ]: 790925642 : prd_symrel=ptsymrel(:,:,prd_ptsymm)
467 : 60840434 : prd_symafm = symafm(sym1) * symafm(sym2)
468 [ + + + + : 1703532152 : 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 : 60840434 : iseq = .False.
474 [ - + ]: 90996920 : do ilist_symrel=1,nlist_symrel(prd_ptsymm)
475 : 90996920 : sym3=list_symrel(ilist_symrel,prd_ptsymm)
476 : 90996920 : iseq = isinteger(prd_tnons(1) - tnons_(1,sym3), tnons_tol_)
477 [ + + ]: 90996920 : if(iseq)then
478 : 71121176 : iseq = isinteger(prd_tnons(2) - tnons_(2,sym3), tnons_tol_)
479 [ + + ]: 71121176 : if(iseq)then
480 : 61538832 : iseq = isinteger(prd_tnons(3) - tnons_(3,sym3), tnons_tol_)
481 [ + + ]: 61538832 : if(iseq)then
482 : 60856970 : iseq = (prd_symafm == symafm(sym3))
483 [ + + ]: 60856970 : if(iseq)then
484 : : ! The test is positive
485 [ - + ]: 60840434 : 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 [ - + ]: 61919460 : 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 [ + - ]: 1125561 : 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 : 46536 : ABI_FREE(list_symrel)
528 : 46536 : ABI_FREE(ptsymm)
529 [ - + ]: 46536 : ABI_FREE(tnons_)
530 : : !write(std_out,*)' m_symtk%sg_multable : exit '
531 : :
532 [ + + - + ]: 82577 : 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 : 17105 : 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 : 17105 : rmet2=zero
593 [ + + ]: 68420 : do ii=1,3
594 [ + + ]: 222365 : do jj=1,3
595 : 205260 : rmet2=rmet2+rmet(ii,jj)**2
596 : : end do
597 : : end do
598 : :
599 : : !Loop over all symmetry operations
600 [ + + ]: 388702 : 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 [ + + ]: 1486504 : do ii=1,3
605 : : rprimd_sym(:,ii)=symrel(1,ii,isym)*rprimd(:,1)+&
606 : : symrel(2,ii,isym)*rprimd(:,2)+&
607 [ + + ]: 4831138 : 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 [ + + ]: 1486504 : do ii=1,3
612 : : rmet_sym(ii,:)=rprimd_sym(1,ii)*rprimd_sym(1,:)+&
613 : : rprimd_sym(2,ii)*rprimd_sym(2,:)+&
614 [ + + ]: 4831138 : rprimd_sym(3,ii)*rprimd_sym(3,:)
615 : : end do
616 : :
617 : 371626 : residual=zero
618 [ + + ]: 1486504 : do ii=1,3
619 [ + + ]: 4831138 : do jj=1,3
620 : 4459512 : residual=residual+(rmet_sym(ii,jj)-rmet(ii,jj))**2
621 : : end do
622 : : end do
623 : :
624 [ + + ]: 371626 : 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 [ + + ]: 1486504 : do ii=1,3
662 : : prods(ii,:)=rprimd_sym(1,ii)*gprimd(1,:)+ &
663 : : rprimd_sym(2,ii)*gprimd(2,:)+ &
664 [ + + ]: 4831138 : rprimd_sym(3,ii)*gprimd(3,:)
665 : : end do
666 : :
667 [ + + ]: 1486504 : do ii=1,3
668 [ + + ]: 4831138 : do jj=1,3
669 : 3344634 : residual=prods(ii,jj)-anint(prods(ii,jj))
670 [ - + ]: 4459512 : 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 [ + + ]: 388702 : if(iexit==-1) exit
687 : : end do ! isym
688 : :
689 [ + + ]: 17105 : if(iexit==1)iexit=0
690 : :
691 : : !write(std_out,'(a)') ' chkorthsy : exit '
692 : :
693 : 17105 : 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 : 4340 : 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 [ + + + + ]: 158710 : 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 : 16284 : multi=0
736 [ + + ]: 481336 : do isym=1,nsym
737 [ + + ]: 465052 : 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 : 16284 : abs(symafm(isym)-1) == 0 )then
744 : 17252 : multi=multi+1
745 [ + + ]: 17252 : if(present(is_translation))then
746 : 4637 : is_translation(isym)=1
747 : : endif
748 : : end if
749 : : end do
750 : :
751 : : !Check whether the cell is primitive
752 [ + + ]: 16284 : 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 [ + + ]: 16284 : 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 : 56297 : 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 : 112594 : 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 : 56297 : ierr_=0
818 : :
819 : : !Compute the coordinates of rprimd_new in the system defined by rprimd(:,:)
820 : 56297 : call matr3inv(rprimd,rprimd_invt)
821 [ + + ]: 225188 : do ii=1,3
822 : : coord(:,ii)=rprimd_new(1,ii)*rprimd_invt(1,:)+ &
823 : : rprimd_new(2,ii)*rprimd_invt(2,:)+ &
824 [ + + ]: 731861 : rprimd_new(3,ii)*rprimd_invt(3,:)
825 : : end do
826 : :
827 : : !Transform symmetry matrices in the system defined by rprimd_new
828 : 56297 : call matr3inv(coord,coordinvt)
829 [ + + ]: 1754506 : do isym=1,nsym
830 [ + + ]: 6792836 : do ii=1,3
831 : : matr1(:,ii)=symrel(:,1,isym)*coord(1,ii)+&
832 : : symrel(:,2,isym)*coord(2,ii)+&
833 [ + + ]: 22076717 : symrel(:,3,isym)*coord(3,ii)
834 : : end do
835 [ + + ]: 6792836 : do ii=1,3
836 : : matr2(:,ii)=coordinvt(1,:)*matr1(1,ii)+&
837 : : coordinvt(2,:)*matr1(2,ii)+&
838 [ + + ]: 22076717 : 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 [ + + ]: 6849133 : do ii=1,3
846 [ + + ]: 22076717 : do jj=1,3
847 : 15283881 : 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 [ + + ]: 15283881 : if(abs(val-nint(val))>ten*tolsym)then
850 : 6241 : ierr_=1
851 [ - + ]: 6241 : 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 : 20378508 : 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 [ + + + + : 21979354 : if(ierr_==0) symrel(:,:,:)=symrel_tmp(:,:,:)
+ + + + ]
868 : :
869 [ + + ]: 56297 : if(present(ierr)) ierr=ierr_
870 : :
871 : 56297 : 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 : 314075 : 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 [ + + ]: 314075 : my_prtvol=0; if (PRESENT(prtvol)) my_prtvol=prtvol
932 : :
933 : : ! Initialise the array symq
934 [ + + + + : 97805678 : symq = 0
+ + ]
935 : :
936 : : isym = symafm(1) ! just to fool abirules and use symafm for the moment
937 : :
938 [ + + ]: 9176948 : 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 [ + + ]: 26902694 : do itirev=1,2
943 : 17725746 : isign=3-2*itirev ! isign is 1 without time-reversal, -1 with time-reversal
944 : :
945 : : ! Get the symmetric of the vector
946 [ + + ]: 70902984 : do ii=1,3
947 : : qsym(ii)=qpt(1)*isign*symrec(ii,1,isym)&
948 : : +qpt(2)*isign*symrec(ii,2,isym)&
949 : 70902984 : +qpt(3)*isign*symrec(ii,3,isym)
950 : : end do
951 : :
952 : : ! Get the difference between the symmetric and the original vector
953 : 17725746 : symq(4,itirev,isym)=1
954 [ + + ]: 70902984 : do ii=1,3
955 : 53177238 : difq(ii)=qsym(ii)-qpt(ii)
956 : : ! Project modulo 1 in the interval ]-1/2,1/2] such that difq = reduce + shift
957 : 53177238 : call wrap2_pmhalf(difq(ii),reduce,shift(ii))
958 [ + + ]: 70902984 : 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 [ + + ]: 17725746 : 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 [ + + ]: 26588619 : if(symq(4,itirev,isym)/=0)then
975 [ + + ]: 1302705 : 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 [ + + ]: 5210820 : do ii=1,3
982 : 5210820 : 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 : 314075 : timrev=1
991 [ + + ]: 1256300 : 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 [ + + ]: 1256300 : if(abs(qpt(ii))>tol)timrev=0
997 : : end do
998 : :
999 [ + + + + ]: 314075 : 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 : 314075 : 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 : 71817 : 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 [ - + ]: 71817 : 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 [ + + ]: 287268 : do ii=1,3
1178 : : metric(:,ii)=cell_base(1,:)*cell_base(1,ii)+&
1179 : : & cell_base(2,:)*cell_base(2,ii)+&
1180 [ + + ]: 933621 : & cell_base(3,:)*cell_base(3,ii)
1181 : : end do
1182 : :
1183 : : !Examine the angles and vector lengths
1184 : 71817 : ang90(:)=0
1185 [ + + ]: 71817 : if(metric(1,2)**2<tolsym**2*metric(1,1)*metric(2,2))ang90(3)=1
1186 [ + + ]: 71817 : if(metric(1,3)**2<tolsym**2*metric(1,1)*metric(3,3))ang90(2)=1
1187 [ + + ]: 71817 : if(metric(2,3)**2<tolsym**2*metric(2,2)*metric(3,3))ang90(1)=1
1188 : 71817 : orth=0
1189 [ + + + + : 71817 : if(ang90(1)==1 .and. ang90(2)==1 .and. ang90(3)==1) orth=1
+ + ]
1190 : 71817 : equal(:)=0
1191 [ + + ]: 71817 : if(abs(metric(1,1)-metric(2,2))<tolsym*half*(metric(1,1)+metric(2,2)))equal(3)=1
1192 [ + + ]: 71817 : if(abs(metric(1,1)-metric(3,3))<tolsym*half*(metric(1,1)+metric(3,3)))equal(2)=1
1193 [ + + ]: 71817 : if(abs(metric(2,2)-metric(3,3))<tolsym*half*(metric(2,2)+metric(3,3)))equal(1)=1
1194 : 71817 : allequal=0
1195 [ + + + + : 71817 : 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 : 71817 : foundc=0
1207 [ - + ]: 71817 : if(abs(iholohedry)==1) foundc=1
1208 [ + + + + ]: 71817 : if(abs(iholohedry)==2 .and. ang90(1)+ang90(3)==2 ) foundc=1
1209 [ + + + + ]: 71817 : if(abs(iholohedry)==3 .and. orth==1) foundc=1
1210 [ + + + + : 71817 : if(abs(iholohedry)==4 .and. orth==1 .and. &
+ + + + +
+ ]
1211 : 5359 : & (equal(3)==1 .or. equal(2)==1 .or. equal(1)==1) ) foundc=1
1212 : : if(abs(iholohedry)==5 .and. allequal==1 .and. &
1213 [ + + + + : 71817 : & (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 [ + + + + : 71817 : & 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 [ + + + + ]: 71817 : 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 [ + + + + : 71817 : 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 : 71817 : 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 : 8870 : 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 : 8870 : unitmat=0
1451 : 8870 : unitmat(1,1)=1 ; unitmat(2,2)=1 ; unitmat(3,3)=1
1452 : :
1453 [ + + ]: 298142 : do isym=1,nsym
1454 : : !write(std_out,'(a,i4,9i3,3es16.6)') ' isym,symrel,tnons=',isym,symrel(:,:,isym),tnons(:,isym)
1455 : :
1456 [ + + + + : 4628352 : symrel_mult(:,:)=symrel(:,:,isym) ; tnons_mult(:)=tnons(:,isym)
+ + ]
1457 : 648826 : order=0
1458 : : !Determine the order of the operation
1459 [ - + ]: 648826 : do iorder=1,48
1460 [ + + + + : 25953040 : symrel_mult(:,:)=matmul(symrel(:,:,isym),symrel_mult(:,:))
+ + ]
1461 [ + + + + : 18167128 : tnons_mult(:)=matmul(symrel(:,:,isym),tnons_mult(:))+tnons(:,isym)
+ + + + +
+ ]
1462 [ + + + + : 8434738 : if(sum(abs(symrel_mult-unitmat))==0)then
+ + ]
1463 : : if(abs(tnons_mult(1)-nint(tnons_mult(1)))<tolsym*iorder .and. &
1464 [ + + + + : 293771 : 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 : 289272 : order=iorder+1
1468 : : !Now, adjust the tnons vector, in order to obtain the exact identity operation at order "order"
1469 [ + + ]: 1157088 : tnons_mult(:)=(tnons_mult(:)-nint(tnons_mult(:)))/(dble(order))
1470 [ + + ]: 289272 : if(abs(tnons_mult(1))>1.00001e-8) tnons(1,isym)=tnons(1,isym)-tnons_mult(1)
1471 [ + + ]: 289272 : if(abs(tnons_mult(2))>1.00001e-8) tnons(2,isym)=tnons(2,isym)-tnons_mult(2)
1472 [ + + ]: 289272 : if(abs(tnons_mult(3))>1.00001e-8) tnons(3,isym)=tnons(3,isym)-tnons_mult(3)
1473 : 289272 : exit
1474 : : endif
1475 : : endif
1476 : :
1477 : : enddo ! iorder
1478 : :
1479 [ - + ]: 298142 : 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 : 8870 : 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 : 28209 : 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 : 43550 : 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 [ + + ]: 43550 : if(present(tnons_new))then
1569 [ + + + + ]: 1793485 : tnons_new(:,1:nsym)=tnons(:,1:nsym)
1570 : : endif
1571 [ + + ]: 43550 : if(present(fixed_mismatch))fixed_mismatch=0
1572 [ + + ]: 43550 : if(present(mismatch_fft_tnons))mismatch_fft_tnons=0
1573 : :
1574 [ + + ]: 43550 : 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 [ - + + - : 87705 : ABI_MALLOC(xredsym,(3,natom))
+ - + - -
+ - + ]
1587 [ + + + + ]: 390219 : xredsym(:,:)=xred(:,1:natom)
1588 : :
1589 [ + + ]: 29235 : if(present(indsym))then
1590 : :
1591 : : ! Loop over atoms to determine new, symmetrized positions.
1592 [ + + ]: 51469 : do iatom=1,natom
1593 : 36018 : tsum(:)=0.0d0
1594 : : !
1595 : : ! Loop over symmetries
1596 [ + + ]: 1021500 : 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 : 985482 : ib=indsym(4,isym,iatom)
1600 : : ! Find the reduced coordinates after translation=t(indsym)+transl
1601 [ + + ]: 3941928 : 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 [ + + ]: 3941928 : dble(symrel(:,3,isym))*fc(3)+ tnons(:,isym)
1606 : :
1607 : : ! Average over nominally equivalent atomic positions
1608 [ + + ]: 3977946 : tsum(:)=tsum(:)+tt(:)
1609 : : end do ! isym
1610 : : !
1611 : : ! Set symmetrized result to sum over number of terms
1612 [ + + ]: 159523 : 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 [ + + - + ]: 29235 : if(present(tolsym) .and. present(tnons_new) .and. present(fixed_mismatch) .and. present(mismatch_fft_tnons))then
1623 : 13784 : fixed_mismatch=0
1624 : 13784 : mismatch_fft_tnons=0
1625 : : !The use of tolsym here is only to favor 0.5 over -0.5
1626 [ + + + + ]: 1785700 : tnons_new(:,:)=tnons(:,:)-nint(tnons(:,:)-tolsym)
1627 [ + + ]: 455926 : do isym=1,nsym
1628 : : mismatch_fft_tnons_current=0
1629 [ + + ]: 1768676 : do ii=1,3
1630 [ + + ]: 6632535 : delta(:)=tnons(ii,isym)*mult(:)
1631 [ + + ]: 6632535 : delta(:)=delta(:)-nint(delta(:))
1632 : : ! Is there is a mismatch between FFT and isym for all multipliers ?
1633 [ + + + + : 1983926 : if( all(abs(delta(:))>tol8*mult(:)) ) mismatch_fft_tnons_current=1
+ + ]
1634 : : enddo
1635 : : !Declare the first symmetry operation that induces a problem
1636 [ + + + + ]: 442169 : 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 : 13784 : 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 [ + + ]: 13784 : 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 [ + + ]: 116940 : transl(:)=xredsym(:,1)-nint(xredsym(:,1))
1694 : :
1695 : : ! Compute the smallest translation to an integer
1696 [ + + ]: 90246 : do jj=2,natom
1697 [ + + ]: 273279 : do ii=1,3
1698 : 183033 : diff=xredsym(ii,jj)-nint(xredsym(ii,jj))
1699 [ + + ]: 244044 : 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 [ + + ]: 116940 : do ii=1,3
1706 [ + + ]: 116940 : if (abs(transl(ii))>1e-13) transl(ii)=0.0
1707 : : end do
1708 : :
1709 : : ! Execute translation
1710 [ + + ]: 119481 : do jj=1,natom
1711 [ + + ]: 390219 : do ii=1,3
1712 : 360984 : 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 [ + + ]: 119481 : do jj=1,natom
1720 [ + + ]: 390219 : do ii=1,3
1721 [ + + ]: 360984 : if (abs(xredsym(ii,jj)-xred(ii,jj))>1E-15) dissimilar=.TRUE.
1722 : : end do
1723 : : end do
1724 : :
1725 [ + + + + : 35231 : if (dissimilar) xred(:,:)=xredsym(:,:)
+ + ]
1726 [ - + ]: 29235 : 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 [ + + ]: 43550 : 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 : 3154584 : 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 : 3154584 : testmn=1000000.d0
1798 : :
1799 : : !Loop through atoms. when types agree, check for agreement after primitive translation
1800 : 3154584 : jatom=1
1801 [ + + ]: 8860018 : do iatom=1,natom
1802 [ + + ]: 8855256 : if (trtypat/=typat(iatom)) cycle
1803 : :
1804 : : ! Check all three components
1805 : 7311939 : test1=tratom(1)-xred(1,iatom)
1806 : 7311939 : test2=tratom(2)-xred(2,iatom)
1807 : 7311939 : test3=tratom(3)-xred(3,iatom)
1808 : : ! Find nearest integer part of difference
1809 : 7311939 : trans1=nint(test1)
1810 : 7311939 : trans2=nint(test2)
1811 : 7311939 : trans3=nint(test3)
1812 : : ! Check whether, after translation, they agree
1813 : 7311939 : test1=test1-dble(trans1)
1814 : 7311939 : test2=test2-dble(trans2)
1815 : 7311939 : test3=test3-dble(trans3)
1816 : 7311939 : test=abs(test1)+abs(test2)+abs(test3)
1817 [ + + ]: 7316701 : if (test<tol10) then
1818 : : ! Note that abs() is not taken here
1819 : 3149822 : difmin(1)=test1
1820 : 3149822 : difmin(2)=test2
1821 : 3149822 : difmin(3)=test3
1822 : 3149822 : jatom=iatom
1823 : 3149822 : transl(1)=trans1
1824 : 3149822 : transl(2)=trans2
1825 : 3149822 : transl(3)=trans3
1826 : : ! Break out of loop when agreement is within tolerance
1827 : 3149822 : exit
1828 : : else
1829 : : ! Keep track of smallest difference if greater than tol10
1830 [ + + ]: 4162117 : if (test<testmn) then
1831 : 1629688 : testmn=test
1832 : : ! Note that abs() is not taken here
1833 : 1629688 : difmin(1)=test1
1834 : 1629688 : difmin(2)=test2
1835 : 1629688 : difmin(3)=test3
1836 : 1629688 : jatom=iatom
1837 : 1629688 : transl(1)=trans1
1838 : 1629688 : transl(2)=trans2
1839 : 1629688 : 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 : 3154584 : eatom=jatom
1846 : :
1847 : 3154584 : 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 : 51132 : 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 : 51132 : err=zero
1946 : 51132 : errout=0
1947 : :
1948 [ + + ]: 1244369 : do isym=1,nsym
1949 [ + + ]: 3885737 : do iatom=1,natom
1950 : :
1951 : : ! Apply inverse transformation to original coordinates.
1952 : : ! Note TRANSPOSE of symrec that corresponds to symrel^{-1}
1953 : :
1954 [ + + ]: 10565472 : 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 : 10565472 : & +dble(symrec(3,mu,isym))*(xred(3,iatom)-tnons(3,isym))
1958 : : end do
1959 : :
1960 : : ! Find symmetrically equivalent atom
1961 : 2641368 : call symchk(difmin,eatom,natom,tratom,transl,typat(iatom),typat,xred)
1962 : :
1963 : : ! Put information into array indsym: translations and label
1964 : 2641368 : indsym(1,isym,iatom)=transl(1)
1965 : 2641368 : indsym(2,isym,iatom)=transl(2)
1966 : 2641368 : indsym(3,isym,iatom)=transl(3)
1967 : 2641368 : indsym(4,isym,iatom)=eatom
1968 : :
1969 : : ! Keep track of maximum difference between transformed coordinates and
1970 : : ! nearest "target" coordinate
1971 : 2641368 : difmax=max(abs(difmin(1)),abs(difmin(2)),abs(difmin(3)))
1972 : 2641368 : err=max(err,difmax)
1973 : :
1974 [ + + ]: 2641368 : 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 [ + + + + ]: 6475973 : 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 [ + + ]: 51132 : print_indsym_ = -1; if (present(print_indsym)) print_indsym_ = print_indsym
2021 [ + + ]: 51132 : if (natom <= print_indsym_) then
2022 [ + + ]: 29569 : do iatom=1,natom
2023 : 20108 : write(msg, '(a,i0,a)' )' symatm: atom number ',iatom,' is reached starting at atom'
2024 : 20108 : call wrtout(std_out,msg)
2025 [ + + ]: 54638 : do ii=1,(nsym-1)/24+1
2026 [ + - ]: 25069 : if(natom<100)then
2027 : 25069 : 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 : 45177 : call wrtout(std_out,msg)
2032 : : end do
2033 : : end do
2034 : : end if
2035 : :
2036 [ + + ]: 51132 : 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 [ - + ]: 51132 : 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 : 51132 : 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 : 854521 : 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 : 854521 : identity(:,:)=0
2097 : 854521 : identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
2098 : 854521 : trial(:,:)=identity(:,:)
2099 : 854521 : matrix(:,:)=symrel(:,:)
2100 : :
2101 : 854521 : order=0
2102 [ - + ]: 2353690 : do iorder=1,6
2103 [ + + + + : 94147600 : trial=matmul(matrix,trial)
+ + ]
2104 [ + + + + : 30597970 : if(sum((trial-identity)**2)==0)then
+ + ]
2105 : 717867 : order=iorder
2106 : 854521 : exit
2107 : : end if
2108 [ + + + + : 21265699 : if(sum((trial+identity)**2)==0)then
+ + ]
2109 : 136654 : order=iorder
2110 : 136654 : exit
2111 : : end if
2112 : : end do
2113 : :
2114 [ - + ]: 854521 : 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 [ + + ]: 854521 : 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 : 431532 : identified=0
2126 : 431532 : trialt(:)=zero
2127 [ + + ]: 1608062 : do ii=1,order
2128 [ + + + + : 33374372 : trialt(:)=matmul(symrel(:,:),trialt(:))+tnons(:)
+ + + + +
+ ]
2129 : : end do
2130 : : ! Gives the associated translation, with components in the interval [-0.5,0.5] .
2131 [ + + ]: 1726128 : reduced(:)=trialt(:)-nint(trialt(:)-tol6)
2132 : :
2133 [ + + + + ]: 1726128 : if(sum(abs(reduced(:)))<tol6)identified=1
2134 [ + + + + : 2139987 : if( (center==1 .or. center==-3) .and. sum(abs(reduced(:)-(/zero,half,half/)))<tol6 )identified=2
+ + ]
2135 [ + + + + : 2140052 : if( (center==2 .or. center==-3) .and. sum(abs(reduced(:)-(/half,zero,half/)))<tol6 )identified=3
+ + ]
2136 [ + + + + : 2139530 : if( (center==3 .or. center==-3) .and. sum(abs(reduced(:)-(/half,half,zero/)))<tol6 )identified=4
+ + ]
2137 [ + + + + : 1726128 : 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 [ - + ]: 428562 : 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 [ + + ]: 1726128 : trialt(:)=trialt(:)/order
2147 [ + + ]: 1726128 : trialt(:)=trialt(:)-nint(trialt(:)-tol6)
2148 : :
2149 : : ! Analyse the resulting vector.
2150 : : identified=0
2151 [ + - ]: 787756 : do ii=1,order
2152 [ + + ]: 3151024 : reduced(:)=ii*trialt(:)-nint(ii*trialt(:)-tol6)
2153 [ + + + + ]: 3151024 : if(sum(abs(reduced(:)))<tol6)identified=1
2154 [ + + + + : 3915560 : if( (center==1 .or. center==-3) .and. sum(abs(reduced(:)-(/zero,half,half/)))<tol6 )identified=2
+ + ]
2155 [ + + + + : 3915625 : if( (center==2 .or. center==-3) .and. sum(abs(reduced(:)-(/half,zero,half/)))<tol6 )identified=3
+ + ]
2156 [ + + + + : 3915105 : if( (center==3 .or. center==-3) .and. sum(abs(reduced(:)-(/half,half,zero/)))<tol6 )identified=4
+ + ]
2157 [ + + + + : 3151024 : if(center==-1.and. sum(abs(reduced(:)-(/half,half,half/)))<tol6 )identified=5
+ + ]
2158 : :
2159 [ + + ]: 784786 : if(identified/=0)then
2160 : 431532 : tnons_order=ii
2161 : 431532 : 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 : 25972 : select case(order)
2169 : : case (1) ! point symmetry 1
2170 [ + + ]: 25972 : if(identified==1) then
2171 : 11655 : type_axis=8 ! 1
2172 : 11655 : write(label,'(a)') 'the identity'
2173 : : else
2174 : 14317 : type_axis=7 ! t
2175 : 14317 : write(label,'(a)') 'a pure translation '
2176 : : end if
2177 : :
2178 : 405560 : 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 [ + + - ]: 431532 : call symaxes(center,iholohedry,isym,symrel,label,order,tnons_order,trialt,type_axis)
2185 : : end select
2186 : :
2187 [ + - ]: 422989 : 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 : 17896 : select case(order)
2192 : : case (1) ! point symmetry 1
2193 : 17896 : type_axis=5 ! -1
2194 : 17896 : write(label,'(a)') 'an inversion'
2195 : : case (2) ! point symmetry 2 - planes
2196 : 290371 : call symplanes(center,iholohedry,isym,symrel,tnons,label,type_axis)
2197 : : case (3) ! point symmetry 3
2198 : 118758 : type_axis=3 ! -3
2199 : 118758 : write(label,'(a)') 'a -3 axis '
2200 : : case (4) ! point symmetry 1
2201 : 114284 : type_axis=2 ! -4
2202 : 114284 : write(label,'(a)') 'a -4 axis '
2203 : : case (6) ! point symmetry 1
2204 : 438 : type_axis=1 ! -6
2205 [ + + + + : 422989 : 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 : 405560 : 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 : 565314 : select case(ordersym)
2342 : : case (2) ! point symmetry 2
2343 : : ! Must characterize directiontype for cP, tP, tI, and hP Bravais lattices
2344 : 159754 : directiontype=1
2345 [ + + ]: 159754 : if( iholohedry==4 .or. iholohedry==7) then ! tP or cP Bravais lattices
2346 [ + + ]: 152058 : 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 [ + + + + ]: 158689 : 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 [ + + + - ]: 89319 : 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 [ + + + + : 149570 : else if(tnons_order==1 .or. (iholohedry==4 .and. center==-1) .or. iholohedry==5)then
+ + - + ]
2377 : 80208 : type_axis=9 ! 2
2378 : 80208 : write(label,'(a)') 'a 2-axis '
2379 : : else
2380 : 69362 : type_axis=20 ! 2_1
2381 : 69362 : write(label,'(a)') 'a 2_1-axis '
2382 : : end if
2383 : :
2384 : : case (3) ! point symmetry 3
2385 [ + + ]: 153422 : if(tnons_order==1)then
2386 : 99506 : type_axis=10 ! 3
2387 : 99506 : write(label,'(a)') 'a 3-axis '
2388 [ + + ]: 53916 : 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 : 53558 : type_axis=10 ! 3, 3_1 or 3_2, undistinguishable
2393 : 53558 : 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 [ + + ]: 91760 : if(tnons_order==1)then
2410 : 22370 : type_axis=12 ! 4
2411 : 22370 : write(label,'(a)') 'a 4-axis '
2412 [ + + ]: 69390 : else if(tnons_order==2)then
2413 : 15464 : type_axis=25 ! 4_2
2414 : 15464 : write(label,'(a)') 'a 4_2-axis '
2415 [ + + ]: 53926 : else if(center/=0)then
2416 : 53812 : type_axis=24 ! 4_1 or 4_3
2417 : 53812 : 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 [ + + + + : 405560 : 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 : 405560 : 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 : 171613 : 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 : 171613 : identity(:,:)=0
2590 : 171613 : 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 : 171613 : mirrorx(:,:)=identity(:,:) ; mirrorx(1,1)=-1
2596 : 171613 : mirrory(:,:)=identity(:,:) ; mirrory(2,2)=-1
2597 : 171613 : mirrorz(:,:)=identity(:,:) ; mirrorz(3,3)=-1
2598 : 171613 : mirrorxy(:,:)=0 ; mirrorxy(1,2)=1 ; mirrorxy(2,1)=1 ; mirrorxy(3,3)=1
2599 : 171613 : mirrorzx(:,:)=0 ; mirrorzx(1,3)=1 ; mirrorzx(3,1)=1 ; mirrorzx(2,2)=1
2600 : 171613 : mirroryz(:,:)=0 ; mirroryz(2,3)=1 ; mirroryz(3,2)=1 ; mirroryz(1,1)=1
2601 : 171613 : mirrormxy(:,:)=0 ; mirrormxy(1,2)=-1 ; mirrormxy(2,1)=-1 ; mirrormxy(3,3)=1
2602 : 171613 : mirrormzx(:,:)=0 ; mirrormzx(1,3)=-1 ; mirrormzx(3,1)=-1 ; mirrormzx(2,2)=1
2603 : 171613 : 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 : 171613 : directiontype=0
2612 : : !The sum of elements of the matrices allow to characterize them
2613 [ + + + + ]: 2230969 : sum_elements=sum(isymrelconv(:,:))
2614 : :
2615 [ + + ]: 171613 : if(sum_elements==1)then
2616 : : ! The mirror plane perpendicular to the c axis is always primary
2617 [ + + + + : 689000 : 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 [ + + + + : 463359 : else if(sum(isymrelconv(:,:))==1)then
+ - ]
2623 [ + + - ]: 35643 : if( iholohedry==2 .or. iholohedry==3 .or. iholohedry==7 )then
2624 : : directiontype=1
2625 : : else if(iholohedry==4 .or. iholohedry==6)then
2626 : 171613 : 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 [ + + ]: 171613 : 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 [ + + ]: 171613 : 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 [ + + ]: 171613 : if( sum_elements==-1)then
2642 [ + + ]: 57609 : if(iholohedry==4 .or. iholohedry==7)directiontype=3
2643 [ + + ]: 57609 : 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 [ + + + + : 4805164 : 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 [ + + ]: 686452 : trialt(:)=trialt(:)*half
2655 [ + + ]: 686452 : 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 : 171613 : 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 [ + + + + : 686452 : if(sum(abs(trialt(:)))<nzero .and. iholohedry/=6)then
+ + ]
2681 : 61752 : type_axis=15 ! m
2682 : 61752 : write(label,'(a)') 'a mirror plane'
2683 [ + + + + ]: 109861 : 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 [ + + + - ]: 108930 : 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 [ + + ]: 107741 : 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 [ + + ]: 107620 : 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 [ + + ]: 105385 : else if(iholohedry==7)then ! cP
2770 : :
2771 [ + + ]: 103100 : 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 [ + + + + : 1712028 : 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 : 26604 : type_axis=17 ! d
2779 : 26604 : 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 [ + - ]: 65882 : else if(directiontype==3)then
2785 [ + + + + ]: 263528 : 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 [ + + + + : 460670 : 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 [ + + + - ]: 255248 : else if( abs(sum(abs(trialt(:)))-one) < nzero )then
2792 : 63812 : type_axis=15 ! tertiary m
2793 : 63812 : 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 : 171613 : 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 : 53265 : 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 [ + + + + : 2130600 : rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
+ + ]
2883 : :
2884 : 53265 : nvecta(1)=2 ; nvectb(1)=3
2885 : 53265 : nvecta(2)=1 ; nvectb(2)=3
2886 : 53265 : nvecta(3)=1 ; nvectb(3)=2
2887 : :
2888 : 53265 : minim(:,:)=rprimd(:,:)
2889 : 53265 : 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 [ + - ]: 65819 : 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 [ + + ]: 263276 : do itrial=1,3
2911 : :
2912 : 197457 : ia=nvecta(itrial) ; ib=nvectb(itrial)
2913 : : ! Make sure the scalar product is negative
2914 [ + + ]: 197457 : if(metmin(ia,ib)>tol8)then
2915 [ + + ]: 346144 : minim(:,ia)=-minim(:,ia)
2916 : 86536 : metmin(ia,ib)=-metmin(ia,ib) ; metmin(ib,ia)=-metmin(ib,ia)
2917 : 86536 : metmin(ia,itrial)=-metmin(ia,itrial)
2918 : 86536 : metmin(itrial,ia)=-metmin(itrial,ia)
2919 : : end if
2920 : : ! Compute the length of the sum vector
2921 : 197457 : 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 [ + + ]: 263276 : if(length2/metmin(ia,ia) < one-tol8)then
2924 [ + + ]: 73812 : minim(:,ia)=minim(:,ia)+minim(:,ib)
2925 : 18453 : metmin(ia,ia)=length2
2926 : 18453 : metmin(ia,ib)=metmin(ia,ib)+metmin(ib,ib)
2927 : 18453 : metmin(ia,itrial)=metmin(ia,itrial)+metmin(ib,itrial)
2928 : 18453 : metmin(ib,ia)=metmin(ia,ib)
2929 : 18453 : metmin(itrial,ia)=metmin(ia,itrial)
2930 : 18453 : minimal=0
2931 : : ! Replace the second vector by the sum vector if the latter is smaller
2932 [ + + ]: 179004 : 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 [ + + ]: 65819 : if(minimal==1)exit
2945 : : end do
2946 : :
2947 [ - + ]: 53265 : 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 [ + - ]: 53604 : 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 : 53604 : ilong=1
2970 [ + + ]: 53604 : if( metmin(2,2)/metmin(1,1) > one + tol8 )ilong=2
2971 [ + + ]: 53604 : if( metmin(3,3)/metmin(ilong,ilong) > one + tol8)ilong=3
2972 : :
2973 : : ! Try combination with all same signs
2974 : 53604 : minimal=1
2975 [ + + + + ]: 696852 : metsum=sum(metmin(:,:))
2976 : 53604 : itrial=0
2977 [ + + ]: 53604 : 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 : 53265 : metsum=two*(metmin(1,1)+metmin(2,2)+metmin(3,3))-metsum
2985 [ + + ]: 213060 : do itrial=1,3
2986 : 159795 : ia=nvecta(itrial) ; ib=nvectb(itrial)
2987 [ - + ]: 213060 : 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 [ + + ]: 53604 : if(minimal==1)exit
3026 : :
3027 : : end do
3028 : :
3029 [ - + ]: 53265 : 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 : 53265 : scprod(:)=zero
3050 [ + + ]: 213060 : do ii=1,3
3051 [ + + ]: 692445 : scprod(:)=scprod(:)+minim(ii,:)*rprimd(ii,1)
3052 : : end do
3053 : : !Determine the vector with the maximal absolute overlap
3054 : 53265 : itrial=1
3055 [ + + ]: 53265 : if(abs(scprod(2))>abs(scprod(1))+tol8)itrial=2
3056 [ + + ]: 53265 : if(abs(scprod(3))>abs(scprod(itrial))+tol8)itrial=3
3057 : : !Switch the vectors if needed
3058 [ + + ]: 53265 : 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 [ + + + + ]: 64201 : 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 [ + + ]: 159795 : scprod(2:3)=zero
3077 [ + + ]: 213060 : do ii=1,3
3078 [ + + ]: 532650 : 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 : 53265 : itrial=2
3082 [ + + ]: 53265 : 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 [ + + + + ]: 140985 : 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 : 53265 : & minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3))
3096 [ + + + + ]: 55877 : if(determinant<-tol8)minim(:,3)=-minim(:,3)
3097 [ - + ]: 53265 : 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 [ + + ]: 213060 : do ii=1,3
3103 [ + + ]: 692445 : 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 : 53265 : 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 : : !!***
|
