Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Crystal-symmetry routines originating from the K290/ACMI code.
10 : !>
11 : !> GROUP1, PGL1, ATFTM1, and ROT1 are based on the Worlton-Warren
12 : !> Computer Physics Communications package (1971, 1974).
13 : ! **************************************************************************************************
14 : MODULE kpsym
15 :
16 : USE kinds, ONLY: dp
17 : USE mathlib, ONLY: invmat
18 : USE string_utilities, ONLY: xstring
19 : #include "./base/base_uses.f90"
20 :
21 : IMPLICIT NONE
22 : PRIVATE
23 :
24 : PUBLIC :: GROUP1s
25 :
26 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpsym'
27 :
28 : ! **************************************************************************************************
29 :
30 : CONTAINS
31 :
32 : ! **************************************************************************************************
33 : !> \brief ...
34 : !> \param iout ...
35 : !> \param a1 ...
36 : !> \param a2 ...
37 : !> \param a3 ...
38 : !> \param nat ...
39 : !> \param ty ...
40 : !> \param x ...
41 : !> \param b1 ...
42 : !> \param b2 ...
43 : !> \param b3 ...
44 : !> \param ihg ...
45 : !> \param ihc ...
46 : !> \param isy ...
47 : !> \param li ...
48 : !> \param nc ...
49 : !> \param indpg ...
50 : !> \param ib ...
51 : !> \param ntvec ...
52 : !> \param v ...
53 : !> \param f0 ...
54 : !> \param r ...
55 : !> \param tvec ...
56 : !> \param origin ...
57 : !> \param rx ...
58 : !> \param isc ...
59 : !> \param delta ...
60 : ! **************************************************************************************************
61 958 : SUBROUTINE group1s(iout, a1, a2, a3, nat, ty, x, b1, b2, b3, &
62 : ihg, ihc, isy, li, nc, indpg, ib, ntvec, &
63 958 : v, f0, r, tvec, origin, rx, isc, delta)
64 : ! ==--------------------------------------------------------------==
65 : ! == WRITTEN ON SEPTEMBER 10TH - FROM THE ACMI COMPLEX ==
66 : ! == (WORLTON AND WARREN, COMPUT.PHYS.COMMUN. 8,71-84 (1974)) ==
67 : ! == (AND 3,88-117 (1972)) ==
68 : ! == BASIC CRYSTALLOGRAPHIC INFORMATION ==
69 : ! == ABOUT A GIVEN CRYSTAL STRUCTURE. ==
70 : ! == SUBROUTINES NEEDED: PGL1,ATFTM1,ROT1,RLV3 ==
71 : ! ==--------------------------------------------------------------==
72 : ! == INPUT DATA: ==
73 : ! == IOUT ... NUMBER OF THE OUTPUT UNIT FOR ON-LINE PRINTING ==
74 : ! == OF VARIOUS MESSAGES ==
75 : ! == IF IOUT<=0 NO MESSAGE ==
76 : ! == A1,A2,A3 .. ELEMENTARY TRANSLATIONS OF THE LATTICE, IN SOME ==
77 : ! == UNIT OF LENGTH ==
78 : ! == NAT .... NUMBER OF ATOMS IN THE UNIT CELL ==
79 : ! == ALL THE DIMENSIONS ARE SET FOR NAT <= 20 ==
80 : ! == TY ..... INTEGERS DISTINGUISHING BETWEEN THE ATOMS OF ==
81 : ! == DIFFERENT TYPE. TY(I) IS THE TYPE OF THE I-TH ATOM ==
82 : ! == OF THE BASIS ==
83 : ! == X ...... CARTESIAN COORDINATES OF THE NAT ATOMS OF THE BASIS ==
84 : ! == DELTA... REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
85 : ! ==--------------------------------------------------------------==
86 : ! == OUTPUT DATA: ==
87 : ! == B1,B2,B3 .. RECIPROCAL LATTICE VECTORS, NOT MULTIPLIED BY ==
88 : ! == ANY 2PI, IN UNITS RECIPROCAL TO THOSE OF A1,A2,A3 ==
89 : ! == IHG .... POINT GROUP OF THE PRIMITIVE LATTICE, HOLOHEDRAL ==
90 : ! == GROUP NUMBER: ==
91 : ! == IHG=1 STANDS FOR TRICLINIC SYSTEM ==
92 : ! == IHG=2 STANDS FOR MONOCLINIC SYSTEM ==
93 : ! == IHG=3 STANDS FOR ORTHORHOMBIC SYSTEM ==
94 : ! == IHG=4 STANDS FOR TETRAGONAL SYSTEM ==
95 : ! == IHG=5 STANDS FOR CUBIC SYSTEM ==
96 : ! == IHG=6 STANDS FOR TRIGONAL SYSTEM ==
97 : ! == IHG=7 STANDS FOR HEXAGONAL SYSTEM ==
98 : ! == IHC .... CODE DISTINGUISHING BETWEEN HEXAGONAL AND CUBIC ==
99 : ! == GROUPS ==
100 : ! == IHC=0 STANDS FOR HEXAGONAL GROUPS ==
101 : ! == IHC=1 STANDS FOR CUBIC GROUPS ==
102 : ! == ISY .... CODE INDICATING WHETHER THE SPACE GROUP IS ==
103 : ! == SYMMORPHIC OR NONSYMMORPHIC ==
104 : ! == ISY= 0 NONSYMMORPHIC GROUP ==
105 : ! == ISY= 1 SYMMORPHIC GROUP ==
106 : ! == ISY=-1 SYMMORPHIC GROUP WITH NON-STANDARD ORIGIN ==
107 : ! == ISY=-2 UNDETERMINED (NORMALLY NEVER) ==
108 : ! == THE GROUP IS CONSIDERED SYMMORPHIC IF FOR EACH ==
109 : ! == OPERATION OF THE POINT GROUP THE SUM OF THE 3 ==
110 : ! == COMPONENTS OF ABS(V(N)) (NONPRIMITIVE TRANSLATION, ==
111 : ! == SEE BELOW) IS LT. 0.0001 ==
112 : ! == ORIGIN STANDARD ORIGIN IF SYMMORPHIC (CRYSTAL COORDINATES) ==
113 : ! == LI ..... CODE INDICATING WHETHER THE POINT GROUP ==
114 : ! == OF THE CRYSTAL CONTAINS INVERSION OR NOT ==
115 : ! == (OPERATIONS 13 OR 25 IN RESPECTIVELY HEXAGONAL ==
116 : ! == OR CUBIC GROUPS). ==
117 : ! == LI=0 MEANS: DOES NOT CONTAIN INVERSION ==
118 : ! == LI>0 MEANS: THERE IS INVERSION IN THE POINT ==
119 : ! == GROUP OF THE CRYSTAL ==
120 : ! == NC ..... TOTAL NUMBER OF ELEMENTS IN THE POINT GROUP OF THE ==
121 : ! == CRYSTAL ==
122 : ! == INDPG .. POINT GROUP INDEX (DETERMINED IF SYMMORPHIC GROUP) ==
123 : ! == IB ..... LIST OF THE ROTATIONS CONSTITUTING THE POINT GROUP ==
124 : ! == OF THE CRYSTAL. THE NUMBERING IS THAT DEFINED IN ==
125 : ! == WORLTON AND WARREN, I.E. THE ONE MATERIALIZED IN THE==
126 : ! == ARRAY R (SEE BELOW) ==
127 : ! == ONLY THE FIRST NC ELEMENTS OF THE ARRAY IB ARE ==
128 : ! == MEANINGFUL ==
129 : ! == NTVEC .. NUMBER OF TRANSLATIONAL VECTORS ==
130 : ! == ASSOCIATED WITH IDENTITY OPERATOR I.E. ==
131 : ! == GIVES THE NUMBER OF IDENTICAL PRIMITIVE CELLS ==
132 : ! == V ...... NONPRIMITIVE TRANSLATIONS (IN THE CASE OF NONSYMMOR-==
133 : ! == PHIC GROUPS). V(I,N) IS THE I-TH COMPONENT ==
134 : ! == OF THE TRANSLATION CONNECTED WITH THE N-TH ELEMENT ==
135 : ! == OF THE POINT GROUP (I.E. WITH THE ROTATION ==
136 : ! == NUMBER IB(N) ). ==
137 : ! == ATTENTION: V(I) ARE NOT CARTESIAN COMPONENTS, ==
138 : ! == THEY REFER TO THE SYSTEM A1,A2,A3. ==
139 : ! == F0 ..... THE FUNCTION DEFINED IN MARADUDIN, IPATOVA BY ==
140 : ! == EQ. (3.2.12): ATOM TRANSFORMATION TABLE. ==
141 : ! == THE ELEMENT F0(N,KAPA) MEANS THAT THE N-TH ==
142 : ! == OPERATION OF THE SPACE GROUP (I.E. OPERATION NUMBER ==
143 : ! == IB(N), TOGETHER WITH AN EVENTUAL NONPRIMITIVE ==
144 : ! == TRANSLATION V(N)) TRANSFERS THE ATOM KAPA INTO THE ==
145 : ! == ATOM F0(N,KAPA). ==
146 : ! == THE 49TH LINE GIVES EQUIVALENT ATOMS FOR ==
147 : ! == FRACTIONAl TRANSLATIONS ASSOCIATED WITH IDENTITY ==
148 : ! == R ...... LIST OF THE 3 X 3 ROTATION MATRICES ==
149 : ! == (XYZ REPRESENTATION OF THE O(H) OR D(6)H GROUPS) ==
150 : ! == ALL 48 OR 24 MATRICES ARE LISTED. ==
151 : ! == FOLLOW NOTATION OF WORLTON-WARREN(1972) ==
152 : ! == TVEC .. LIST OF NTVEC TRANSLATIONAL VECTORS ==
153 : ! == ASSOCIATED WITH IDENTITY OPERATOR ==
154 : ! == TVEC(1:3,1) = \(0,0,0\) ==
155 : ! == (CRYSTAL COORDINATES) ==
156 : ! == RX ..... SCRATCH ARRAY ==
157 : ! == ISC .... SCRATCH ARRAY ==
158 : ! ==--------------------------------------------------------------==
159 : ! == PRINTED OUTPUT: ==
160 : ! == PROGRAM PRINTS THE TYPE OF THE LATTICE (IHG, IN WORDS), ==
161 : ! == LISTS THE OPERATIONS OF THE POINT GROUP OF THE ==
162 : ! == CRYSTAL, INDICATES WHETHER THE SPACE GROUP IS SYMMORPHIC OR ==
163 : ! == NONSYMMORPHIC AND WHETHER THE POINT GROUP OF THE CRYSTAL ==
164 : ! == CONTAINS INVERSION. ==
165 : ! ==--------------------------------------------------------------==
166 : INTEGER :: iout
167 : REAL(dp) :: a1(3), a2(3), a3(3)
168 : INTEGER :: nat, ty(nat)
169 : REAL(dp) :: x(3, nat), b1(3), b2(3), b3(3)
170 : INTEGER :: ihg, ihc, isy, li, nc, indpg, ib(48), &
171 : ntvec
172 : REAL(dp) :: v(3, 48)
173 : INTEGER :: f0(49, nat)
174 : REAL(dp) :: r(3, 3, 48), tvec(3, nat), origin(3), &
175 : rx(3, nat)
176 : INTEGER :: isc(nat)
177 : REAL(dp) :: delta
178 :
179 : INTEGER :: i, ncprim
180 : REAL(dp) :: a(3, 3), ai(3, 3), ap(3, 3), api(3, 3)
181 :
182 3832 : DO i = 1, 3
183 2874 : a(i, 1) = a1(i)
184 2874 : a(i, 2) = a2(i)
185 3832 : a(i, 3) = a3(i)
186 : END DO
187 : ! ==--------------------------------------------------------------==
188 : ! == A(I,J) IS THE I-TH CARTESIAN COMPONENT OF THE J-TH PRIMITIVE ==
189 : ! == TRANSLATION VECTOR OF THE DIRECT LATTICE ==
190 : ! == TY(I) IS AN INTEGER DISTINGUISHING ATOMS OF DIFFERENT TYPE, ==
191 : ! == I.E., DIFFERENT ATOMIC SPECIES ==
192 : ! == X(J,I) IS THE J-TH CARTESIAN COMPONENT OF THE POSITION ==
193 : ! == VECTOR FOR THE I-TH ATOM IN THE UNIT CELL. ==
194 : ! ==--------------------------------------------------------------==
195 : ! ==DETERMINE PRIMITIVE LATTICE VECTORS FOR THE RECIPROCAL LATTICE==
196 : ! ==--------------------------------------------------------------==
197 958 : CALL calbrec(a, ai)
198 3832 : DO i = 1, 3
199 2874 : b1(i) = ai(1, i)
200 2874 : b2(i) = ai(2, i)
201 3832 : b3(i) = ai(3, i)
202 : END DO
203 : ! ==--------------------------------------------------------------==
204 : ! Determination of the translation vectors associated with
205 : ! the Identity matrix i.e. if the cell is duplicated
206 : ! Give also the ``primitive lattice''
207 958 : CALL primlatt(a, ai, ap, api, nat, ty, x, ntvec, tvec, f0, isc, delta)
208 : ! ==--------------------------------------------------------------==
209 : ! Determination of the holohedral group (and crystal system)
210 958 : CALL pgl1(ap, api, ihc, nc, ib, ihg, r, delta)
211 958 : IF (ntvec > 1) THEN
212 : ! All rotations found by PGL1 have axes in x, y or z cart. axis
213 : ! So we have too check if we do not loose symmetry
214 186 : ncprim = nc
215 : ! The hexagonal system is found if the z axis is the sixfold axis
216 186 : CALL pgl1(a, ai, ihc, nc, ib, ihg, r, delta)
217 186 : IF (ncprim > nc) THEN
218 : ! More symmetry with
219 0 : CALL pgl1(ap, api, ihc, nc, ib, ihg, r, delta)
220 : END IF
221 : END IF
222 :
223 : ! Determination of the space group
224 : CALL atftm1(iout, r, v, x, f0, origin, ib, ty, nat, ihg, ihc, rx, &
225 958 : nc, indpg, ntvec, a, ai, li, isy, isc, delta)
226 :
227 958 : IF (iout > 0) THEN
228 335 : IF (li > 0) THEN
229 : IF (iout > 0) THEN
230 : WRITE (iout, '(1X,A)') &
231 146 : 'KPSYM| THE POINT GROUP OF THE CRYSTAL CONTAINS THE INVERSION'
232 : END IF
233 : END IF
234 335 : IF (iout > 0) THEN
235 335 : WRITE (iout, *)
236 : END IF
237 : END IF
238 :
239 958 : END SUBROUTINE group1s
240 : ! **************************************************************************************************
241 : !> \brief ...
242 : !> \param a ...
243 : !> \param ai ...
244 : ! **************************************************************************************************
245 1322 : SUBROUTINE calbrec(a, ai)
246 : ! ==--------------------------------------------------------------==
247 : ! == CALCULATE RECIPROCAL VECTOR BASIS (AI(1:3,1:3)) ==
248 : ! == INPUT: ==
249 : ! == A(3,3) A(I,J) IS THE I-TH CARTESIAN COMPONENT ==
250 : ! == OF THE J-TH PRIMITIVE TRANSLATION VECTOR OF ==
251 : ! == THE DIRECT LATTICE ==
252 : ! == OUTPUT: ==
253 : ! == AI(3,3) RECIPROCAL VECTOR BASIS ==
254 : ! ==--------------------------------------------------------------==
255 : REAL(dp) :: a(3, 3), ai(3, 3)
256 :
257 : INTEGER :: i, il, iu, j, jl, ju
258 : REAL(dp) :: det
259 :
260 : det = a(1, 1)*a(2, 2)*a(3, 3) + a(2, 1)*a(1, 3)*a(3, 2) + &
261 : a(3, 1)*a(1, 2)*a(2, 3) - a(1, 1)*a(2, 3)*a(3, 2) - &
262 1322 : A(2, 1)*A(1, 2)*A(3, 3) - A(3, 1)*A(1, 3)*A(2, 2)
263 1322 : det = 1._dp/det
264 5288 : DO i = 1, 3
265 3966 : il = 1
266 3966 : iu = 3
267 3966 : IF (i == 1) il = 2
268 2644 : IF (i == 3) iu = 2
269 17186 : DO j = 1, 3
270 11898 : jl = 1
271 11898 : ju = 3
272 11898 : IF (j == 1) jl = 2
273 7932 : IF (j == 3) ju = 2
274 : ai(j, i) = (-1._dp)**(i + j)*det* &
275 15864 : (A(IL, JL)*A(IU, JU) - A(IL, JU)*A(IU, JL))
276 : END DO
277 : END DO
278 : ! ==--------------------------------------------------------------==
279 1322 : RETURN
280 : END SUBROUTINE calbrec
281 : ! ==================================================================
282 : ! **************************************************************************************************
283 : !> \brief ...
284 : !> \param a ...
285 : !> \param ai ...
286 : !> \param ap ...
287 : !> \param api ...
288 : !> \param nat ...
289 : !> \param ty ...
290 : !> \param x ...
291 : !> \param ntvec ...
292 : !> \param tvec ...
293 : !> \param f0 ...
294 : !> \param isc ...
295 : !> \param delta ...
296 : ! **************************************************************************************************
297 958 : SUBROUTINE primlatt(a, ai, ap, api, nat, ty, x, ntvec, tvec, f0, isc, delta)
298 : ! ==--------------------------------------------------------------==
299 : ! == DETERMINATION OF THE TRANSLATION VECTORS ASSOCIATED WITH ==
300 : ! == THE IDENTITY SYMMETRY I.E. IF THE CELL IS DUPLICATED ==
301 : ! == GIVE ALSO THE PRIMITIVE DIRECT AND RECIPROCAL LATTICE VECTOR ==
302 : ! ==--------------------------------------------------------------==
303 : ! == INPUT: ==
304 : ! == A(3,3) A(I,J) IS THE I-TH CARTESIAN COMPONENT ==
305 : ! == OF THE J-TH TRANSLATION VECTOR OF ==
306 : ! == THE DIRECT LATTICE ==
307 : ! == AI(3,3) RECIPROCAL VECTOR BASIS (CARTESIAN) ==
308 : ! == NAT NUMBER OF ATOMS ==
309 : ! == TY(NAT) TYPE OF ATOMS ==
310 : ! == X(3,NAT) ATOMIC COORDINATES IN CARTESIAN COORDINATES ==
311 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
312 : ! == OUTPUT: ==
313 : ! == AP(3,3) COMPONENTS OF THE PRIMITIVE TRANSLATION VECTORS ==
314 : ! == API(3,3) PRIMITIVE RECIPROCAL BASIS VECTORS ==
315 : ! == BOTH BAISI ARE IN CARTESIAN COORDINATES ==
316 : ! == NTVEC NUMBER OF TRANSLATION VECTORS (FRACTIONNAL) ==
317 : ! == TVEC(3,NTVEC) COMPONENTS OF TRANSLATIONAL VECTORS ==
318 : ! == (CRYSTAL COORDINATES) ==
319 : ! == F0(49,NAT) GIVES INEQUIVALENT ATOM FOR EACH ATOM ==
320 : ! == THE 49-TH LINE ==
321 : ! == ISC(NAT) SCRATCH ARRAY ==
322 : ! ==--------------------------------------------------------------==
323 : REAL(dp) :: a(3, 3), ai(3, 3), ap(3, 3), api(3, 3)
324 : INTEGER :: nat, ty(nat)
325 : REAL(dp) :: x(3, nat)
326 : INTEGER :: ntvec
327 : REAL(dp) :: tvec(3, nat)
328 : INTEGER :: f0(49, nat), isc(nat)
329 : REAL(dp) :: delta
330 :
331 : INTEGER :: i, il, iv, j, k2
332 : LOGICAL :: oksym
333 : REAL(dp) :: vr(3), xb(3)
334 :
335 : ! Variables
336 : ! ==--------------------------------------------------------------==
337 : ! First we check if there exist fractional translational vectors
338 : ! associated with Identity operation i.e.
339 : ! if the cell is duplicated or not.
340 :
341 958 : ntvec = 1
342 958 : tvec(1, 1) = 0._dp
343 958 : tvec(2, 1) = 0._dp
344 958 : tvec(3, 1) = 0._dp
345 6050 : DO i = 1, nat
346 6050 : f0(49, i) = i
347 : END DO
348 5092 : DO k2 = 2, nat
349 4134 : IF (ty(1) /= ty(k2)) CYCLE
350 14344 : DO i = 1, 3
351 14344 : xb(i) = x(i, k2) - x(i, 1)
352 : END DO
353 : ! A fractional translation vector VR is defined.
354 3586 : CALL rlv3(ai, xb, vr, il, delta)
355 3586 : CALL checkrlv3(1, nat, ty, x, x, vr, f0, ai, isc, .TRUE., oksym, delta)
356 4544 : IF (oksym) THEN
357 : ! A fractional translational vector is found
358 542 : ntvec = ntvec + 1
359 : ! F0(49,1:NAT) gives number of equivalent atoms
360 : ! and has atom indexes of inequivalent atoms (for translation)
361 4830 : DO i = 1, nat
362 4830 : IF (f0(49, i) > f0(1, i)) f0(49, i) = f0(1, i)
363 : END DO
364 2168 : DO i = 1, 3
365 2168 : tvec(i, ntvec) = vr(i)
366 : END DO
367 : END IF
368 : END DO
369 : ! ==-------------------------------------------------------------==
370 3832 : DO i = 1, 3
371 2874 : ap(1, i) = a(1, i)
372 2874 : ap(2, i) = a(2, i)
373 2874 : ap(3, i) = a(3, i)
374 2874 : api(1, i) = ai(1, i)
375 2874 : api(2, i) = ai(2, i)
376 3832 : api(3, i) = ai(3, i)
377 : END DO
378 958 : IF (ntvec == 1) THEN
379 : ! The current cell is definitely a primitive one
380 : ! Copy A and AI to AP and API
381 : ELSE
382 : ! We are looking for the primitive lattice vector basis set
383 : ! AP is our current lattice vector basis
384 728 : DO iv = 2, ntvec
385 : ! TVEC in cartesian coordinates
386 2168 : DO i = 1, 3
387 : xb(i) = tvec(1, iv)*a(i, 1) &
388 : + TVEC(2, IV)*A(I, 2) &
389 2168 : + TVEC(3, IV)*A(I, 3)
390 : END DO
391 : ! We calculare TVEC in AP basis
392 542 : CALL rlv3(api, xb, vr, il, delta)
393 1440 : DO i = 1, 3
394 1254 : IF (ABS(vr(i)) > delta) THEN
395 364 : il = NINT(1._dp/ABS(vr(i)))
396 364 : IF (il > 1) THEN
397 : ! We replace AP(1:3,I) by TVEC(1:3,IV)
398 1456 : DO j = 1, 3
399 1456 : ap(j, i) = xb(j)
400 : END DO
401 : ! Calculate new API
402 364 : CALL calbrec(ap, api)
403 364 : EXIT
404 : END IF
405 : END IF
406 : END DO
407 : END DO
408 : END IF
409 : ! ==--------------------------------------------------------------==
410 958 : RETURN
411 : END SUBROUTINE primlatt
412 : ! ==================================================================
413 : ! **************************************************************************************************
414 : !> \brief ...
415 : !> \param a ...
416 : !> \param ai ...
417 : !> \param ihc ...
418 : !> \param nc ...
419 : !> \param ib ...
420 : !> \param ihg ...
421 : !> \param r ...
422 : !> \param delta ...
423 : ! **************************************************************************************************
424 1144 : SUBROUTINE pgl1(a, ai, ihc, nc, ib, ihg, r, delta)
425 : ! ==--------------------------------------------------------------==
426 : ! == WRITTEN ON SEPTEMBER 11TH, 1979 - FROM ACMI COMPLEX ==
427 : ! == AUXILIARY SUBROUTINE TO GROUP1 ==
428 : ! == SUBROUTINE PGL DETERMINES THE POINT GROUP OF THE LATTICE ==
429 : ! == AND THE CRYSTAL SYSTEM. ==
430 : ! == SUBROUTINES NEEDED: ROT1, RLV3 ==
431 : ! ==--------------------------------------------------------------==
432 : ! == WARNING: FOR THE HEXAGONAL SYSTEM, THE 3RD AXIS SUPPOSE ==
433 : ! == TO BE THE SIX-FOLD AXIS ==
434 : ! ==--------------------------------------------------------------==
435 : ! == INPUT: ==
436 : ! == A ..... DIRECT LATTICE VECTORS ==
437 : ! == AI .... RECIPROCAL LATTICE VECTORS ==
438 : ! == DELTA.. REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
439 : ! ==--------------------------------------------------------------==
440 : ! == OUTPUT: ==
441 : ! == IHC .... CODE DISTINGUISHING BETWEEN HEXAGONAL AND CUBIC ==
442 : ! == GROUPS ==
443 : ! == IHC=0 STANDS FOR HEXAGONAL GROUPS ==
444 : ! == IHC=1 STANDS FOR CUBIC GROUPS ==
445 : ! == NC .... NUMBER OF ROTATIONS IN THE POINT GROUP ==
446 : ! == IB .... SET OF ROTATION ==
447 : ! == IHG .... POINT GROUP OF THE PRIMITIVE LATTICE, HOLOHEDRAL ==
448 : ! == GROUP NUMBER: ==
449 : ! == IHG=1 STANDS FOR TRICLINIC SYSTEM ==
450 : ! == IHG=2 STANDS FOR MONOCLINIC SYSTEM ==
451 : ! == IHG=3 STANDS FOR ORTHORHOMBIC SYSTEM ==
452 : ! == IHG=4 STANDS FOR TETRAGONAL SYSTEM ==
453 : ! == IHG=5 STANDS FOR CUBIC SYSTEM ==
454 : ! == IHG=6 STANDS FOR TRIGONAL SYSTEM ==
455 : ! == IHG=7 STANDS FOR HEXAGONAL SYSTEM ==
456 : ! == R ...... LIST OF THE 3 X 3 ROTATION MATRICES ==
457 : ! == (XYZ REPRESENTATION OF THE O(H) OR D(6)H GROUPS) ==
458 : ! == ALL 48 OR 24 MATRICES ARE LISTED. ==
459 : ! == FOLLOW NOTATION OF WORLTON-WARREN(1972) ==
460 : ! ==--------------------------------------------------------------==
461 : REAL(dp) :: a(3, 3), ai(3, 3)
462 : INTEGER :: ihc, nc, ib(48), ihg
463 : REAL(dp) :: r(3, 3, 48), delta
464 :
465 : INTEGER :: i, j, k, lx, n, nr
466 : REAL(dp) :: tr, vr(3), xa(3)
467 :
468 2216 : DO ihc = 0, 1
469 : ! IHC is 0 for hexagonal groups and 1 for cubic groups.
470 2216 : IF (ihc == 0) THEN
471 : nr = 24
472 : ELSE
473 1072 : nr = 48
474 : END IF
475 2216 : nc = 0
476 : ! Constructs rotation operations.
477 2216 : CALL rot1(ihc, r)
478 81128 : loop_rotation: DO n = 1, nr
479 78912 : ib(n) = 0
480 : ! Rotate the A1,2,3 vectors by rotation No. N
481 215968 : DO k = 1, 3
482 690048 : DO i = 1, 3
483 517536 : xa(i) = 0._dp
484 2242656 : DO j = 1, 3
485 2070144 : xa(i) = xa(i) + r(i, j, n)*a(j, k)
486 : END DO
487 : END DO
488 172512 : CALL rlv3(ai, xa, vr, lx, delta)
489 172512 : tr = 0._dp
490 690048 : DO i = 1, 3
491 690048 : tr = tr + ABS(vr(i))
492 : END DO
493 : ! If VR.ne.0, then XA cannot be a multiple of a lattice vector
494 215968 : IF (tr > delta) CYCLE loop_rotation
495 : END DO
496 43456 : nc = nc + 1
497 81128 : ib(nc) = n
498 : END DO loop_rotation
499 : ! ==------------------------------------------------------------==
500 : ! IHG stands for holohedral group number.
501 2216 : IF (ihc == 0) THEN
502 : ! Hexagonal group:
503 1144 : IF (nc == 12) ihg = 6
504 1144 : IF (nc > 12) ihg = 7
505 1144 : IF (nc >= 12) RETURN
506 : ! Too few operations, try cubic group: (IHC=1,NR=48)
507 : ELSE
508 : ! Cubic group:
509 1072 : IF (nc < 4) ihg = 1
510 1072 : IF (nc == 4) ihg = 2
511 1072 : IF (nc > 4) ihg = 3
512 1072 : IF (nc == 16) ihg = 4
513 1072 : IF (nc > 16) ihg = 5
514 1072 : RETURN
515 : END IF
516 : END DO
517 : ! ==--------------------------------------------------------------==
518 : RETURN
519 : END SUBROUTINE pgl1
520 : ! ==================================================================
521 : ! **************************************************************************************************
522 : !> \brief ...
523 : !> \param ai ...
524 : !> \param xb ...
525 : !> \param vr ...
526 : !> \param il ...
527 : !> \param delta ...
528 : ! **************************************************************************************************
529 1466376 : SUBROUTINE rlv3(ai, xb, vr, il, delta)
530 : ! ==--------------------------------------------------------------==
531 : ! == WRITTEN ON SEPTEMBER 11TH, 1979 - FROM ACMI COMPLEX ==
532 : ! == AUXILIARY SUBROUTINE TO GROUP1 ==
533 : ! == SUBROUTINE RLV REMOVES A DIRECT LATTICE VECTOR ==
534 : ! == FROM XB LEAVING THE REMAINDER IN VR. ==
535 : ! == IF A NONZERO LATTICE VECTOR WAS REMOVED, IL IS MADE NONZERO. ==
536 : ! == VR STANDS FOR V-REFERENCE. ==
537 : ! ==--------------------------------------------------------------==
538 : ! == INPUT: ==
539 : ! == AI(I,J) ARE THE RECIPROCAL LATTICE VECTORS, ==
540 : ! == B(I) = AI(I,J),J=1,2,3 ==
541 : ! == XB(1:3) VECTOR IN CARTESIAN COORDINATES ==
542 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
543 : ! == OUTPUT: ==
544 : ! == VR IS NOT GIVEN IN CARTESIAN COORDINATES BUT ==
545 : ! == IN THE SYSTEM A1,A2,A3 (CRYSTAL COORDINATES) ==
546 : ! == AND BETWEEN -1/2 AND 1/2 ==
547 : ! == IL ABS OF VR ==
548 : ! == K.K., 23.10.1979 ==
549 : ! ==--------------------------------------------------------------==
550 : REAL(dp) :: ai(3, 3), xb(3), vr(3)
551 : INTEGER :: il
552 : REAL(dp) :: delta
553 :
554 : INTEGER :: i
555 : REAL(dp) :: ts
556 :
557 1466376 : il = 0
558 5865504 : DO i = 1, 3
559 5865504 : vr(i) = 0._dp
560 : END DO
561 1466376 : ts = ABS(xb(1)) + ABS(xb(2)) + ABS(xb(3))
562 1466376 : IF (ts <= delta) RETURN
563 5491552 : DO i = 1, 3
564 4118664 : vr(i) = vr(i) + ai(i, 1)*xb(1) + ai(i, 2)*xb(2) + ai(i, 3)*xb(3)
565 4118664 : il = il + NINT(ABS(vr(i)))
566 : ! Change in order to have correct determination of origin and
567 : ! symmorphic group (T.D 30/03/98)
568 : ! VR(I) = - MOD(real(VR(I),kind=dp),1._dp)
569 5491552 : vr(i) = NINT(vr(i)) - vr(i)
570 : END DO
571 : ! ==--------------------------------------------------------------==
572 : RETURN
573 : END SUBROUTINE rlv3
574 : ! ==================================================================
575 : ! **************************************************************************************************
576 : !> \brief ...
577 : !> \param iout ...
578 : !> \param r ...
579 : !> \param v ...
580 : !> \param x ...
581 : !> \param f0 ...
582 : !> \param origin ...
583 : !> \param ib ...
584 : !> \param ty ...
585 : !> \param nat ...
586 : !> \param ihg ...
587 : !> \param ihc ...
588 : !> \param rx ...
589 : !> \param nc ...
590 : !> \param indpg ...
591 : !> \param ntvec ...
592 : !> \param a ...
593 : !> \param ai ...
594 : !> \param li ...
595 : !> \param isy ...
596 : !> \param isc ...
597 : !> \param delta ...
598 : ! **************************************************************************************************
599 958 : SUBROUTINE atftm1(iout, r, v, x, f0, origin, ib, ty, nat, ihg, ihc, &
600 958 : rx, nc, indpg, ntvec, a, ai, li, isy, isc, delta)
601 : ! ==--------------------------------------------------------------==
602 : ! == WRITTEN ON SEPTEMBER 11TH, 1979 - FROM ACMI COMPLEX ==
603 : ! == AUXILIARY SUBROUTINE TO GROUP1 ==
604 : ! == SUBROUTINE ATFTMT DETERMINES ==
605 : ! == THE POINT GROUP OF THE CRYSTAL, ==
606 : ! == THE ATOM TRANSFORMATION TABLE,F0, ==
607 : ! == THE FRACTIONAL TRANSLATIONS,V, ==
608 : ! == ASSOCIATED WITH EACH ROTATION. ==
609 : ! == SUBROUTINES NEEDED: RLV3 CHECKRLV3 SYMMORPHIC XSTRING ==
610 : ! == MAY 14TH,1998: A LOT OF CHANGES (ARGUMENTS) ==
611 : ! == BETTER DETERMINATION OF V ==
612 : ! == SEP 15TH,1998: DETERMINATION OF FRACTIONAL TRANSLATIONAL VEC.==
613 : ! ==--------------------------------------------------------------==
614 : ! == INPUT: ==
615 : ! == IOUT Logical file number (output) ==
616 : ! == If IOUT<=0 no message ==
617 : ! == IHG Holohedral group number (determined by PGL1) ==
618 : ! == IHC Code distinguishing between hexagonal and cubic groups==
619 : ! == IHC=0 stands for hexagonal groups ==
620 : ! == IHC=1 stands for cubic groups ==
621 : ! == NC Number of rotation operations ==
622 : ! == NAT Number of atoms (used in the routine) ==
623 : ! == X Coordinates of atoms (cartesian) ==
624 : ! == TY Type of atoms ==
625 : ! == R Sets of transformation operations (cartesian) ==
626 : ! == IB Index giving NC operations in R ==
627 : ! == AI Reciprocal lattice vectors ==
628 : ! == NTVEC Number of translational vectors ==
629 : ! == associated with Identity ==
630 : ! == if primitive cell NTVEC=1, TVEC=(0,0,0) ==
631 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
632 : ! == OUTPUT: ==
633 : ! == RX(3,NAT) Scratch array ==
634 : ! == ISC(NAT) Scratch array ==
635 : ! == NC is modified (number of symmetry operations) ==
636 : ! == INDPG Point group index ==
637 : ! == V(3,48) The fractional translations associated ==
638 : ! == with each rotation (crystal coordinates) ==
639 : ! == F0(1:48,NAT) ==
640 : ! == The atom transformation table for rotation (48,NAT) ==
641 : ! == ORIGIN Standard origin if symmorphic (crystal coordinates) ==
642 : ! == ISY = 1 Isommorphic group ==
643 : ! == =-1 Isommorphic group with non-standard origin ==
644 : ! == = 0 Non-Isommorphic group ==
645 : ! == =-2 Undetermined (normally never) ==
646 : ! == LI ..... Code indicating whether the point group ==
647 : ! == of the crystal contains inversion or not ==
648 : ! == (operations 13 or 25 in respectively hexagonal ==
649 : ! == or cubic groups). ==
650 : ! == LI=0 : does not contain inversion ==
651 : ! == LI>0 : there is inversion in the point ==
652 : ! == group of the crystal ==
653 : ! ==--------------------------------------------------------------==
654 : ! INDPG group indpg group indpg group indpg group ==
655 : ! == 1 1 (c1) 9 3m (c3v) 17 4/mmm(d4h) 25 222(d2) ==
656 : ! == 2 <1>(ci) 10 <3>m(d3d) 18 6 (c6) 26 mm2(c2v) ==
657 : ! == 3 2 (c2) 11 4 (c4) 19 <6>(c3h) 27 mmm(d2h) ==
658 : ! == 4 m (c1h) 12 <4>(s4) 20 6/m(c6h) 28 23 (t) ==
659 : ! == 5 2/m(c2h) 13 4/m(c4h) 21 622(d6) 29 m3 (th) ==
660 : ! == 6 3 (c3) 14 422(d4) 22 6mm(c6v) 30 432(o) ==
661 : ! == 7 <3>(c3i) 15 4mm(c4v) 23 <6>m2(d3h) 31 <4>3m(td) ==
662 : ! == 8 32 (d3) 16 <4>2m(d2d) 24 6/mmm(d6h) 32 m3m(oh) ==
663 : ! ==--------------------------------------------------------------==
664 : ! rname_cubic: Name of 48 rotations (convention Warren-Worlton)
665 : INTEGER :: iout
666 : REAL(dp) :: r(3, 3, 48), v(3, 48), origin(3)
667 : INTEGER :: ib(48), nat, ty(nat), f0(49, nat)
668 : REAL(dp) :: x(3, nat)
669 : INTEGER :: ihg, ihc
670 : REAL(dp) :: rx(3, nat)
671 : INTEGER :: nc, indpg, ntvec
672 : REAL(dp) :: a(3, 3), ai(3, 3)
673 : INTEGER :: li, isy, isc(nat)
674 : REAL(dp) :: delta
675 :
676 : CHARACTER(len=10), DIMENSION(48), PARAMETER :: rname_cubic = [' 1 ', ' 2[ 10 0] ', &
677 : ' 2[ 01 0] ', ' 2[ 00 1] ', ' 3[-1-1-1]', ' 3[ 11-1] ', ' 3[-11 1] ', ' 3[ 1-11] ', &
678 : ' 3[ 11 1] ', ' 3[-11-1] ', ' 3[-1-11] ', ' 3[ 1-1-1]', ' 2[-11 0] ', ' 4[ 00 1] ', &
679 : ' 4[ 00-1] ', ' 2[ 11 0] ', ' 2[ 0-11] ', ' 2[ 01 1] ', ' 4[ 10 0] ', ' 4[-10 0] ', &
680 : ' 2[-10 1] ', ' 4[ 0-10] ', ' 2[ 10 1] ', ' 4[ 01 0] ', '-1 ', '-2[ 10 0] ', &
681 : '-2[ 01 0] ', '-2[ 00 1] ', '-3[-1-1-1]', '-3[ 11-1] ', '-3[-11 1] ', '-3[ 1-11] ', &
682 : '-3[ 11 1] ', '-3[-11-1] ', '-3[-1-11] ', '-3[ 1-1-1]', '-2[-11 0] ', '-4[ 00 1] ', &
683 : '-4[ 00-1] ', '-2[ 11 0] ', '-2[ 0-11] ', '-2[ 01 1] ', '-4[ 10 0] ', '-4[-10 0] ', &
684 : '-2[-10 1] ', '-4[ 0-10] ', '-2[ 10 1] ', '-4[ 01 0] ']
685 : CHARACTER(len=11), DIMENSION(24), PARAMETER :: rname_hexai = [' 1 ', ' 6[ 00 1] ', &
686 : ' 3[ 00 1] ', ' 2[ 00 1] ', ' 3[ 00 -1] ', ' 6[ 00 -1] ', ' 2[ 01 0] ', ' 2[-11 0] ', &
687 : ' 2[ 10 0] ', ' 2[ 21 0] ', ' 2[ 11 0] ', ' 2[ 12 0] ', '-1 ', '-6[ 00 1] ', &
688 : '-3[ 00 1] ', '-2[ 00 1] ', '-3[ 00 -1] ', '-6[ 00 -1] ', '-2[ 01 0] ', '-2[-11 0] ', &
689 : '-2[ 10 0] ', '-2[ 21 0] ', '-2[ 11 0] ', '-2[ 12 0] ']
690 : CHARACTER(len=12), DIMENSION(7), PARAMETER :: icst = ['TRICLINIC ', 'MONOCLINIC ', &
691 : 'ORTHORHOMBIC', 'TETRAGONAL ', 'CUBIC ', 'TRIGONAL ', 'HEXAGONAL ']
692 : CHARACTER(len=3), DIMENSION(32), PARAMETER :: pgrd = ['c1 ', 'ci ', 'c2 ', 'c1h', 'c2h', &
693 : 'c3 ', 'c3i', 'd3 ', 'c3v', 'd3 ', 'c4 ', 's4 ', 'c4h', 'd4 ', 'c4v', 'd2d', 'd4h', 'c6 ',&
694 : 'c3h', 'c6h', 'd6 ', 'c6v', 'd3h', 'd6h', 'd2 ', 'c2v', 'd2h', 't ', 'th ', 'o ', 'td ',&
695 : 'oh ']
696 : CHARACTER(len=5), DIMENSION(32), PARAMETER :: pgrp = [' 1', ' <1>', ' 2', ' m', &
697 : ' 2/m', ' 3', ' <3>', ' 32', ' 3m', ' <3>m', ' 4', ' <4>', ' 4/m', ' 422', &
698 : ' 4mm', '<4>2m', '4/mmm', ' 6', ' <6>', ' 6/m', ' 622', ' 6mm', '<6>m2', '6/mmm', &
699 : ' 222', ' mm2', ' mmm', ' 23', ' m3', ' 432', '<4>3m', ' m3m']
700 :
701 : INTEGER :: i, iis(48), il, info, j, k, k2, l, n, &
702 : nca, ni
703 : LOGICAL :: nodupli, oksym
704 : REAL(dp) :: vc(3, 48), vr(3), vs, xb(3)
705 :
706 958 : nodupli = ntvec == 1
707 958 : nca = 0
708 46942 : DO n = 1, 48
709 46942 : iis(n) = 0
710 : END DO
711 : ! Calculate translational vector for each operation
712 : ! and atom transformation table.
713 28790 : DO n = 1, nc
714 27832 : l = ib(n)
715 27832 : iis(l) = 1
716 169016 : DO k = 1, nat
717 592568 : DO i = 1, 3
718 564736 : rx(i, k) = r(i, 1, l)*x(1, k) + r(i, 2, l)*x(2, k) + r(i, 3, l)*x(3, k)
719 : END DO
720 : END DO
721 111328 : DO k = 1, 3
722 111328 : vr(k) = 0._dp
723 : END DO
724 : ! First we determine for VR=(/0,0,0/)
725 : ! IMPORTANT IF NOT UNIQUE ATOMS FOR DETERMINATION OF SYMMORPHIC
726 27832 : CALL checkrlv3(n, nat, ty, rx, x, vr, f0, ai, isc, nodupli, oksym, delta)
727 27832 : IF (.NOT. oksym) THEN
728 : ! Now we try other possible VR
729 : ! F0(49,1:NAT) has only inequivalent atom indexes for translation
730 96884 : DO k2 = 1, nat
731 85296 : IF (f0(49, k2) < k2) CYCLE
732 74928 : IF (ty(1) /= ty(k2)) CYCLE
733 260912 : DO i = 1, 3
734 260912 : xb(i) = rx(i, 1) - x(i, k2)
735 : END DO
736 : ! A translation vector VR is defined.
737 65228 : CALL rlv3(ai, xb, vr, il, delta)
738 : ! ==----------------------------------------------------------==
739 : ! == SUBROUTINE RLV3 REMOVES A DIRECT LATTICE VECTOR FROM XB ==
740 : ! == LEAVING THE REMAINDER IN VR. IF A NONZERO LATTICE ==
741 : ! == VECTOR WAS REMOVED, IL IS MADE NONZERO. ==
742 : ! == VR STANDS FOR V-REFERENCE. ==
743 : ! == VR IS NOT GIVEN IN CARTESIAN COORDINATES BUT ==
744 : ! == IN THE SYSTEM A1,A2,A3. K.K., 23.10.1979 ==
745 : ! ==----------------------------------------------------------==
746 65228 : CALL checkrlv3(n, nat, ty, rx, x, vr, f0, ai, isc, nodupli, oksym, delta)
747 76816 : IF (oksym) EXIT
748 : END DO
749 15508 : IF (.NOT. oksym) THEN
750 11588 : iis(l) = 0
751 11588 : CYCLE
752 : END IF
753 : END IF
754 16244 : nca = nca + 1
755 65934 : DO i = 1, 3
756 64976 : v(i, nca) = vr(i)
757 : END DO
758 : ! ==------------------------------------------------------------==
759 : ! == V(I,N) IS THE I-TH COMPONENT OF THE FRACTIONAL ==
760 : ! == TRANSLATION ASSOCIATED WITH THE ROTATION N. ==
761 : ! == ATTENTION: V(I) ARE NOT CARTESIAN COMPONENTS, THEY ARE ==
762 : ! == GIVEN IN THE SYSTEM A1,A2,A3. ==
763 : ! == K.K., 23.10. 1979 ==
764 : ! ==------------------------------------------------------------==
765 : END DO
766 : ! Remove unused operations
767 958 : i = 0
768 958 : ni = 13
769 958 : IF (ihg < 6) ni = 25
770 958 : li = 0
771 28790 : DO n = 1, nc
772 27832 : l = ib(n)
773 27832 : IF (iis(l) == 0) CYCLE
774 16244 : i = i + 1
775 16244 : ib(i) = ib(n)
776 16244 : IF (ib(i) == ni) li = i
777 92034 : DO k = 1, nat
778 102664 : f0(i, k) = f0(n, k)
779 : END DO
780 : END DO
781 : ! ==--------------------------------------------------------------==
782 958 : nc = i
783 958 : vs = 0._dp
784 17202 : DO n = 1, nc
785 17202 : vs = vs + ABS(v(1, n)) + ABS(v(2, n)) + ABS(v(3, n))
786 : END DO
787 : ! THE ORIGINAL VALUE DELTA=0.0001 WAS MODIFIED
788 : ! BY K.K. , SEPTEMBER 1979 TO 0.0005
789 : ! AND RETURNED TO 0.0001 BY RJN OCT 1987
790 958 : IF (vs > delta) THEN
791 302 : isy = 0
792 : ELSE
793 656 : isy = 1
794 : END IF
795 : ! ==--------------------------------------------------------------==
796 : ! Determination of the point group
797 : ! (Thierry Deutsch - 1998 [Maybe not complete!!])
798 958 : IF (ihg < 6) THEN
799 886 : IF (nc == 0) THEN
800 0 : IF (iout > 0) THEN
801 0 : WRITE (iout, '(" ATFTM1! IHG=",A," NC=",I2)') icst(ihg), nC
802 : END IF
803 0 : CPABORT('ATFTM1: NUMBER OF ROTATION NULL')
804 : ! Triclinic system
805 886 : ELSE IF (nc == 1) THEN
806 : ! IB=1
807 184 : indpg = 1 ! 1 (c1)
808 702 : ELSE IF (nc == 2 .AND. ib(2) == 25) THEN
809 : ! IB=125
810 16 : indpg = 2 ! <1>(ci)
811 686 : ELSE IF (nc == 2 .AND. ( &
812 : ib(2) == 4 .OR. & ! 2[001]
813 : ib(2) == 2 .OR. & ! 2[100]
814 : ib(2) == 3)) THEN ! 2[010]
815 : ! Monoclinic system
816 : ! IB=14 (z-axis) OR
817 : ! IB=12 (x-axis) OR
818 : ! IB=13 (y-axis)
819 8 : indpg = 3 ! 2 (c2)
820 678 : ELSE IF (nc == 2 .AND. ( &
821 : ib(2) == 28 .OR. &
822 : ib(2) == 26 .OR. &
823 : ib(2) == 27)) THEN
824 : ! IB=128 (z-axis) OR
825 : ! IB=126 (x-axis) OR
826 : ! IB=127 (y-axis)
827 60 : indpg = 4 ! m (c1h)
828 618 : ELSE IF (nc == 4 .AND. ( &
829 : ib(4) == 28 .OR. & ! 2[001]
830 : ib(4) == 27 .OR. & ! 2[010]
831 : ib(4) == 26 .OR. & ! 2[100]
832 : ib(4) == 37 .OR. & ! -2[-110]
833 : ib(4) == 40)) THEN ! 2[110]
834 : ! IB=1 425 28 (z-axis) OR
835 : ! IB=1 225 26 (x-axis) OR
836 : ! IB=1 325 27 (y-axis) OR
837 : ! IB=113 2537 (-xy-axis)OR
838 : ! IB=116 2540 (xy-axis)
839 138 : indpg = 5 ! 2/m(c2h)
840 480 : ELSE IF (nc == 4 .AND. ( &
841 : ib(4) == 15 .OR. &
842 : ib(4) == 20 .OR. &
843 : ib(4) == 24)) THEN
844 : ! Tetragonal system
845 : ! IB=14 1415 (z-axis) OR
846 : ! IB=12 1920 (x-axis) OR
847 : ! IB=13 2224 (y-axis)
848 2 : indpg = 11 ! 4 (c4)
849 478 : ELSE IF (nc == 4 .AND. ( &
850 : ib(4) == 39 .OR. &
851 : ib(4) == 44 .OR. &
852 : ib(4) == 48)) THEN
853 : ! IB=14 3839 (z-axis) OR
854 : ! IB=12 4344 (x-axis) OR
855 : ! IB=13 4648 (y-axis)
856 0 : indpg = 12 ! <4>(s4)
857 478 : ELSE IF (nc == 8 .AND. ( &
858 : (ib(3) == 14 .AND. ib(8) == 39) .OR. &
859 : (ib(3) == 19 .AND. ib(8) == 44) .OR. &
860 : (ib(3) == 22 .AND. ib(8) == 48))) THEN
861 : ! IB=14 1415 2825 3839 (z-axis) OR
862 : ! IB=12 1920 2625 4344 (x-axis) OR
863 : ! IB=13 2224 2725 4648 (y-axis)
864 0 : indpg = 13 ! 422(d4)
865 478 : ELSE IF (nc == 8 .AND. ib(4) == 4 .AND. ( &
866 : ib(8) == 16 .OR. &
867 : ib(8) == 20 .OR. &
868 : ib(8) == 24)) THEN
869 : ! IB=12 3 413 1415 16 (z-axis) OR
870 : ! IB=12 3 417 1920 18 (x-axis) OR
871 : ! IB=12 3 421 2224 23 (y-axis)
872 12 : indpg = 14 ! 4/m(c4h)
873 466 : ELSE IF (nc == 8 .AND. ( &
874 : ib(8) == 40 .OR. &
875 : ib(8) == 42 .OR. &
876 : ib(8) == 47)) THEN
877 : ! IB=14 1415 2627 3740 (z-axis) OR
878 : ! IB=12 1920 2827 4142 (x-axis) OR
879 : ! IB=13 2224 2628 4547 (y-axis)
880 2 : indpg = 15 ! 4mm(c4v)
881 464 : ELSE IF (nc == 8 .AND. ( &
882 : (ib(3) == 13 .AND. ib(8) == 39) .OR. &
883 : (ib(3) == 17 .AND. ib(8) == 44) .OR. &
884 : (ib(3) == 21 .AND. ib(8) == 48))) THEN
885 : ! IB=14 1316 2627 3839 (z-axis) OR
886 : ! IB=12 1718 2827 4344 (x-axis) OR
887 : ! IB=13 2123 2628 4648 (y-axis)
888 0 : indpg = 16 ! <4>2m(d2d)
889 464 : ELSE IF (nc == 16 .AND. ( &
890 : ib(16) == 40 .OR. &
891 : ib(16) == 44 .OR. &
892 : ib(16) == 48)) THEN
893 : ! IB=12 3 413 1415 1625 2627 2837 3839 40 (z-axis) OR
894 : ! IB=12 3 417 1920 1825 2627 2841 4344 42 (x-axis) OR
895 : ! IB=12 3 421 2224 2325 2627 2845 4648 47 (y-axis)
896 48 : indpg = 17 ! 4/mmm(d4h)
897 416 : ELSE IF (nc == 4 .AND. (ib(4) == 4)) THEN
898 : ! Orthorhombic system
899 : ! IB=12 3 4
900 0 : indpg = 25 ! 222(d2)
901 416 : ELSE IF (nc == 4 .AND. ( &
902 : ib(4) == 27 .OR. &
903 : ib(4) == 28)) THEN
904 : ! IB=13 2627 (z-axis) OR
905 : ! IB=12 2728 (x-axis) OR
906 : ! IB=14 2628 (y-axis) OR
907 0 : indpg = 26 ! mm2(c2v)
908 416 : ELSE IF (nc == 8) THEN
909 : ! IB=12 3 425 2627 28
910 0 : indpg = 27 ! mmm(d2h)
911 416 : ELSE IF (nc == 12 .AND. ( &
912 : ib(12) == 12 .OR. &
913 : ib(12) == 47 .OR. &
914 : ib(12) == 45)) THEN
915 : ! Cubic system
916 : ! IB=12 3 4 5 6 7 8 910 1112 OR
917 : ! IB=15 1113 1823 2530 3537 4247 OR
918 : ! IB=18 1016 1821 2532 3440 4245
919 0 : indpg = 28 ! 23 (t)
920 416 : ELSE IF (nc == 24 .AND. ib(24) == 36) THEN
921 : ! IB= 1 2 3 4 5 6 7 8 910 1112
922 : ! 2526 2728 2930 3132 3334 3536
923 0 : indpg = 29 ! m3 (th)
924 416 : ELSE IF (nc == 24 .AND. ib(24) == 24) THEN
925 : ! IB=12 3 45 6 78 9 1011 12
926 : ! 1314 1516 1718 1920 2122 2324
927 0 : indpg = 30 ! 432 (o)
928 416 : ELSE IF (nc == 24 .AND. ib(24) == 48) THEN
929 : ! IB=12 3 45 6 78 9 1011 12
930 : ! 3738 3940 4142 4345 4647 48
931 0 : indpg = 31 ! <4>3m(td)
932 416 : ELSE IF (nc == 48) THEN
933 : ! IB=1..48
934 272 : indpg = 32 ! m3m(oh)
935 : ELSE
936 : ! WRITE(6,'(" ATFTM1! IHG=",A," NC=",I2)') ICST(IHG),NC
937 : ! WRITE(6,'(" ATFTM1!",19I3)') (IB(I),I=1,NC)
938 : ! WRITE(6,'(" ATFTM1! THIS CASE IS UNKNOWN IN THE DATABASE")')
939 : ! Probably a sub-group of 32
940 144 : indpg = -32
941 : END IF
942 : ELSE IF (ihg >= 6) THEN
943 72 : IF (nc == 0) THEN
944 0 : IF (iout > 0) THEN
945 0 : WRITE (iout, '(" ATFTM1! IHG=",A," NC=",I2)') icst(ihg), nC
946 : END IF
947 0 : CPABORT('ATFTM1: NUMBER OF ROTATION NULL')
948 : ! Triclinic system
949 72 : ELSE IF (nc == 1) THEN
950 : ! IB=1
951 0 : indpg = 1 ! 1 (c1)
952 72 : ELSE IF (nc == 2 .AND. ib(2) == 13) THEN
953 : ! IB=113
954 0 : indpg = 2 ! <1>(ci)
955 72 : ELSE IF (nc == 2 .AND. ( &
956 : ib(2) == 4)) THEN ! 2[001]
957 : ! Monoclinic system
958 : ! IB=1 4
959 0 : indpg = 3 ! 2 (c2)
960 72 : ELSE IF (nc == 2 .AND. ( &
961 : ib(2) == 16)) THEN
962 : ! IB=116
963 0 : indpg = 4 ! m (c1h)
964 72 : ELSE IF (nc == 4 .AND. ( &
965 : ib(4) == 24 .OR. &
966 : ib(4) == 20)) THEN
967 : ! IB=112 1324 OR
968 : ! IB=1 813 20
969 0 : indpg = 5 ! 2/m(c2h)
970 72 : ELSE IF (nc == 3 .AND. ib(3) == 5) THEN
971 : ! Trigonal system
972 : ! IB=13 5
973 4 : indpg = 6 ! 3 (c3)
974 68 : ELSE IF (nc == 6 .AND. ib(6) == 17) THEN
975 : ! IB=113 1517 35
976 0 : indpg = 7 ! <3>(c3i)
977 68 : ELSE IF (nc == 6 .AND. ib(6) == 11) THEN
978 : ! IB=17 9 1135
979 0 : indpg = 8 ! 32 (d3)
980 68 : ELSE IF (nc == 6 .AND. ib(6) == 23) THEN
981 : ! IB=13 5 1921 23
982 0 : indpg = 9 ! 3m (c3v)
983 68 : ELSE IF (nc == 12 .AND. ib(12) == 23) THEN
984 : ! IB=13 5 79 1113 1517 1921 23
985 64 : indpg = 10 ! <3>m(d3d)
986 4 : ELSE IF (nc == 6 .AND. ib(6) == 6) THEN
987 : ! Hexagonal system
988 : ! IB=12 3 45 6
989 0 : indpg = 18 ! 6 (c6)
990 4 : ELSE IF (nc == 6 .AND. ib(6) == 18) THEN
991 : ! IB=13 5 1416 18
992 0 : indpg = 19 ! <6>(c3h)
993 4 : ELSE IF (nc == 12 .AND. ib(12) == 18) THEN
994 : ! IB=12 3 45 6 1314 1516 1718
995 0 : indpg = 20 ! 6/m(c6h)
996 4 : ELSE IF (nc == 12 .AND. ib(12) == 12) THEN
997 : ! IB=12 3 45 6 78 9 1011 12
998 0 : indpg = 21 ! 622(d6)
999 4 : ELSE IF (nc == 12 .AND. ib(2) == 2 .AND. ib(12) == 24) THEN
1000 : ! IB=12 3 45 6 1920 2122 2324
1001 0 : indpg = 22 ! 6mm(c6v)
1002 4 : ELSE IF (nc == 12 .AND. ib(2) == 3 .AND. ib(12) == 24) THEN
1003 : ! IB=13 5 79 1114 1618 2022 24
1004 0 : indpg = 23 ! <6>m2(d3h)
1005 4 : ELSE IF (nc == 24) THEN
1006 : ! IB=1..24
1007 4 : indpg = 24 ! 6/mmm(d6h)
1008 : ELSE
1009 : ! Probably a sub-group of 24
1010 : ! WRITE(6,'(" ATFTM1! IHG=",A," NC=",I2)') ICST(IHG),NC
1011 : ! WRITE(6,'(" ATFTM1!",48I3)') (IB(I),I=1,NC)
1012 : ! WRITE(6,'(" ATFTM1! THIS CASE IS UNKNOWN IN THE DATABASE")')
1013 0 : indpg = -24
1014 : END IF
1015 : END IF
1016 : ! ==--------------------------------------------------------------==
1017 : ! == Determination if the space group is symmorphic or not ==
1018 : ! ==--------------------------------------------------------------==
1019 958 : IF (isy /= 1) THEN
1020 : ! Transform V in cartesian coordinates
1021 7930 : DO n = 1, nc
1022 7628 : vc(1, n) = a(1, 1)*v(1, n) + a(1, 2)*v(2, n) + a(1, 3)*v(3, n)
1023 7628 : vc(2, n) = a(2, 1)*v(1, n) + a(2, 2)*v(2, n) + a(2, 3)*v(3, n)
1024 7930 : vc(3, n) = a(3, 1)*v(1, n) + a(3, 2)*v(2, n) + a(3, 3)*v(3, n)
1025 : END DO
1026 302 : CALL symmorphic(nc, ib, r, vc, ai, info, origin, delta)
1027 302 : IF (info == 1) THEN
1028 20 : CALL rlv3(ai, origin, xb, il, delta)
1029 : ! !!!RLV3 determines -XB in crystal coordinates
1030 : ! !!We want between 0.0 and 1.0
1031 80 : DO i = 1, 3
1032 80 : IF (-xb(i) >= 0._dp) THEN
1033 36 : origin(i) = -xb(i)
1034 : ELSE
1035 24 : origin(i) = 1._dp - xb(i)
1036 : END IF
1037 : END DO
1038 80 : DO i = 1, 3
1039 80 : xb(i) = a(i, 1)*origin(1) + a(i, 2)*origin(2) + a(i, 3)*origin(3)
1040 : END DO
1041 20 : isy = -1
1042 282 : ELSE IF (info == 0) THEN
1043 282 : isy = 0
1044 : ELSE
1045 0 : isy = -2
1046 : END IF
1047 : ELSE
1048 2624 : DO i = 1, 3
1049 2624 : origin(i) = 0._dp
1050 : END DO
1051 : END IF
1052 : ! ==--------------------------------------------------------------==
1053 : ! == Output ==
1054 : ! ==--------------------------------------------------------------==
1055 958 : IF (iout > 0) THEN
1056 : IF (iout > 0) THEN
1057 335 : WRITE (iout, *)
1058 : END IF
1059 335 : CALL xstring(icst(ihg), i, j)
1060 335 : IF ((ihg == 7 .AND. nc == 24) .OR. &
1061 : (ihg == 5 .AND. nc == 48)) THEN
1062 70 : IF (iout > 0) THEN
1063 : WRITE (iout, '(A,A,A)') &
1064 70 : ' KPSYM| THE POINT GROUP OF THE CRYSTAL IS THE FULL ', &
1065 70 : icst(ihg) (i:j), &
1066 140 : ' GROUP'
1067 : END IF
1068 : ELSE
1069 265 : IF (iout > 0) THEN
1070 : WRITE (iout, '(A,A,A,I2,A)') &
1071 265 : ' KPSYM| THE CRYSTAL SYSTEM IS ', &
1072 265 : icst(ihg) (i:j), &
1073 530 : ' WITH ', nc, ' OPERATIONS:'
1074 : END IF
1075 265 : IF (ihc == 0) THEN
1076 10 : IF (iout > 0) THEN
1077 121 : WRITE (iout, '( 5(5(A13),/))') (rname_hexai(ib(i)), i=1, nc)
1078 : END IF
1079 : ELSE
1080 255 : IF (iout > 0) THEN
1081 1138 : WRITE (iout, '(10(5(A13),/))') (rname_cubic(ib(i)), i=1, nc)
1082 : END IF
1083 : END IF
1084 : END IF
1085 : ! ==------------------------------------------------------------==
1086 335 : IF (isy == 1) THEN
1087 231 : IF (iout > 0) THEN
1088 : WRITE (iout, '(A)') &
1089 231 : ' KPSYM| THE SPACE GROUP OF THE CRYSTAL IS SYMMORPHIC'
1090 : END IF
1091 104 : ELSE IF (isy == -1) THEN
1092 8 : IF (iout > 0) THEN
1093 : WRITE (iout, '(A)') &
1094 8 : ' KPSYM| THE SPACE GROUP OF THE CRYSTAL IS SYMMORPHIC'
1095 : END IF
1096 8 : IF (iout > 0) THEN
1097 : WRITE (iout, '(A,A,/,T3,3F10.6,3X,3F10.6)') &
1098 8 : ' KPSYM| THE STANDARD ORIGIN OF COORDINATES IS: ', &
1099 16 : '[CARTESIAN] [CRYSTAL]', xb, origin
1100 : END IF
1101 96 : ELSE IF (isy == 0) THEN
1102 96 : IF (iout > 0) THEN
1103 : WRITE (iout, '(A,/,3X,A,F15.6,A)') &
1104 96 : ' KPSYM| THE SPACE GROUP IS NON-SYMMORPHIC,', &
1105 192 : ' (SUM OF TRANSLATION VECTORS=', vs, ')'
1106 : END IF
1107 0 : ELSE IF (isy == -2) THEN
1108 0 : IF (iout > 0) THEN
1109 : WRITE (iout, '(A,A)') &
1110 0 : ' KPSYM| CANNOT DETERMINE IF THE SPACE GROUP IS', &
1111 0 : ' SYMMORPHIC OR NOT'
1112 : END IF
1113 0 : IF (iout > 0) THEN
1114 : WRITE (iout, '(A,/,A,/,3X,A,F15.6,A)') &
1115 0 : ' KPSYM| THE SPACE GROUP IS NON-SYMMORPHIC,', &
1116 0 : ' KPSYM| OR ELSE A NON STANDARD ORIGIN OF COORDINATES WAS USED.', &
1117 0 : ' KPSYM| (SUM OF TRANSLATION VECTORS=', vs, ')'
1118 : END IF
1119 : END IF
1120 335 : IF (indpg > 0) THEN
1121 288 : CALL xstring(pgrp(indpg), i, j)
1122 288 : CALL xstring(pgrd(indpg), k, l)
1123 288 : IF (iout > 0) THEN
1124 : WRITE (iout, '(A,A,"(",A,")",T56,"[INDEX=",I2,"]")') &
1125 288 : ' KPSYM| THE POINT GROUP OF THE CRYSTAL IS ', pgrp(indpg) (i:j), &
1126 576 : pgrd(indpg) (k:l), indpg
1127 : END IF
1128 : ELSE
1129 47 : CALL xstring(pgrp(-indpg), i, j)
1130 47 : CALL xstring(pgrd(-indpg), k, l)
1131 47 : IF (iout > 0) THEN
1132 : WRITE (iout, '(A,I2,A,A,"(",A,")",T56,"[INDEX=",I2,"]")') &
1133 47 : ' KPSYM| POINT GROUP: GROUP ORDER=', nc, &
1134 47 : ' SUBGROUP OF ', pgrp(-indpg) (i:j), &
1135 94 : pgrd(-indpg) (k:l), -indpg
1136 : END IF
1137 : END IF
1138 335 : IF (ntvec == 1) THEN
1139 275 : IF (iout > 0) THEN
1140 : WRITE (iout, '(A,T60,I6)') &
1141 275 : ' KPSYM| NUMBER OF PRIMITIVE CELL:', ntvec
1142 : END IF
1143 : ELSE
1144 60 : IF (iout > 0) THEN
1145 : WRITE (iout, '(A,T60,I6)') &
1146 60 : ' KPSYM| NUMBER OF PRIMITIVE CELLS:', ntvec
1147 : END IF
1148 : END IF
1149 : END IF
1150 :
1151 958 : END SUBROUTINE atftm1
1152 :
1153 : ! **************************************************************************************************
1154 : !> \brief ...
1155 : !> \param n ...
1156 : !> \param nat ...
1157 : !> \param ty ...
1158 : !> \param rx ...
1159 : !> \param x ...
1160 : !> \param vr ...
1161 : !> \param f0 ...
1162 : !> \param ai ...
1163 : !> \param isc ...
1164 : !> \param nodupli ...
1165 : !> \param oksym ...
1166 : !> \param delta ...
1167 : ! **************************************************************************************************
1168 96646 : SUBROUTINE checkrlv3(n, nat, ty, rx, x, vr, f0, ai, isc, &
1169 : nodupli, oksym, delta)
1170 : ! ==--------------------------------------------------------------==
1171 : ! == WRITTEN IN MAY 14TH, 1998 (T.D.) ==
1172 : ! == CHECK IF RX+VR GIVES THE SAME LATTICE AS X ==
1173 : ! == BUILD THE ATOM TRANSFORMATION TABLE ==
1174 : ! ==--------------------------------------------------------------==
1175 : ! == INPUT: ==
1176 : ! == N ROTATION NUMBER (INDEX USED IN F0 BETWEEN 1 AND 48) ==
1177 : ! == NAT NUMBER OF ATOMS ==
1178 : ! == TY(1:NAT) TYPE OF ATOMS ==
1179 : ! == RX(1:3,1:NAT) ATOMIC COORDINATES FROM Nth ROTATION (CART.) ==
1180 : ! == X(1:3,1:NAT) ATOMIC COORDINATES (CARTESIAN) ==
1181 : ! == VR(1:3) TRANSLATION VECTOR (CRYSTAL COOR.) ==
1182 : ! == AI(1:3,1:3) LATTICE RECIPROCAL VECTORS ==
1183 : ! == NODUPLI .TRUE., THE CELL IS A PRIMITIVE ONE ==
1184 : ! == WE CAN SPEED UP ==
1185 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
1186 : ! == OUTPUT: ==
1187 : ! == F0(1:49,1:NAT) ATOM TRANSFORMATION TABLE ==
1188 : ! == F0 IS THE FUNCTION DEFINED IN MARADUDIN AND VOSK0 ==
1189 : ! == BY EQ.(2.35). ==
1190 : ! == IT DEFINES THE ATOM TRANSFORMATION TABLE ==
1191 : ! == OKSYM TRUE IF RX+VR = X ==
1192 : ! == ISC(1:NAT) SCRATCH ARRAY ==
1193 : ! == USED TO SPEED UP THE ROUTINE ==
1194 : ! == EACH ATOM IS ONLY ONCE AN IMAGE ==
1195 : ! == IF NO DUPLICATION OF THE CELL ==
1196 : ! ==--------------------------------------------------------------==
1197 : INTEGER :: n, nat, ty(nat)
1198 : REAL(dp) :: rx(3, nat), x(3, nat), vr(3)
1199 : INTEGER :: f0(49, nat)
1200 : REAL(dp) :: ai(3, 3)
1201 : INTEGER :: isc(nat)
1202 : LOGICAL :: nodupli, oksym
1203 : REAL(dp) :: delta
1204 :
1205 : INTEGER :: ia, ib, il
1206 : REAL(dp) :: tol, vt(3), xb(3)
1207 :
1208 : ! Fractional residuals at the tolerance boundary accumulate Cartesian
1209 : ! rotation and lattice-conversion roundoff. Account for that roundoff so
1210 : ! equivalent operations are not accepted or rejected by a few ulps.
1211 : tol = delta + 32.0_dp*EPSILON(1.0_dp)*MAX(1.0_dp, &
1212 7572020 : MAXVAL(SUM(ABS(ai), DIM=2))*MAX(MAXVAL(ABS(rx)), MAXVAL(ABS(x))))
1213 :
1214 874018 : DO ia = 1, nat
1215 874018 : isc(ia) = 0
1216 : END DO
1217 : ! Now we check if ROT(N)+VR gives a correct symmetry.
1218 305104 : atom: DO ia = 1, nat
1219 1601278 : DO ib = 1, nat
1220 1601278 : IF (ty(ia) == ty(ib) .AND. isc(ib) == 0) THEN
1221 1224424 : xb(1) = rx(1, ia) - x(1, ib)
1222 1224424 : xb(2) = rx(2, ia) - x(2, ib)
1223 1224424 : xb(3) = rx(3, ia) - x(3, ib)
1224 1224424 : CALL rlv3(ai, xb, vt, il, delta)
1225 : ! VT STANDS FOR V-TEST
1226 1990218 : oksym = ALL(ABS((vr - vt) - ANINT(vr - vt)) <= tol)
1227 1224424 : IF (oksym) THEN
1228 208458 : IF (nodupli) isc(ib) = 1
1229 208458 : f0(n, ia) = ib
1230 : ! IR+VR is the good one: another symmetry operation
1231 : ! Next atom
1232 : CYCLE atom
1233 : END IF
1234 : END IF
1235 : END DO
1236 : ! VR is not the correct translation vector
1237 96646 : RETURN
1238 : END DO atom
1239 : END SUBROUTINE checkrlv3
1240 : ! ==================================================================
1241 : ! **************************************************************************************************
1242 : !> \brief ...
1243 : !> \param nc ...
1244 : !> \param ib ...
1245 : !> \param r ...
1246 : !> \param v ...
1247 : !> \param ai ...
1248 : !> \param info ...
1249 : !> \param origin ...
1250 : !> \param delta ...
1251 : ! **************************************************************************************************
1252 302 : SUBROUTINE symmorphic(nc, ib, r, v, ai, info, origin, delta)
1253 : ! ==--------------------------------------------------------------==
1254 : ! == Check if the group is symmorphic with a non-standard origin ==
1255 : ! == WARNING: If there are equivalent atoms, this routine could ==
1256 : ! == not determine if the space group is symmorphic ==
1257 : ! == So you have to check if the solution V=0 works (see ATFTM1) ==
1258 : ! ==--------------------------------------------------------------==
1259 : ! == INPUT: ==
1260 : ! == NC Number of operations ==
1261 : ! == IB(NC) Index of operation in R ==
1262 : ! == R(3,3,48) Rotations ==
1263 : ! == V(3,NC) Fractional translations related to R(3,3,IB(NC)) ==
1264 : ! == R AND V ARE IN CARTESIAN COORDINATES ==
1265 : ! == AI(I,J) ARE THE RECIPROCAL LATTICE VECTORS, ==
1266 : ! == B(I) = AI(I,J),J=1,2,3 ==
1267 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
1268 : ! == ==
1269 : ! == OUTPUT: ==
1270 : ! == ORIGIN(1:3) Give standard origin (cartesian coordinates) ==
1271 : ! == Give the standard origin with smallest coordinates==
1272 : ! == if NTVEC /= 1 ==
1273 : ! == INFO = 1 The group is symmorphic ==
1274 : ! == INFO = 0 The group is not symmorphic ==
1275 : ! == INFO =-1 The routine cannot determine ==
1276 : ! ==--------------------------------------------------------------==
1277 : INTEGER :: nc, ib(nc)
1278 : REAL(dp) :: r(3, 3, 48), v(3, nc), ai(3, 3)
1279 : INTEGER :: info
1280 : REAL(dp) :: origin(3), delta
1281 :
1282 : INTEGER :: i, i1, ierror, igood(3), il, imissing2, &
1283 : imissing3, iok(3), ionly, ir, j, j1
1284 : REAL(dp) :: diag, dif, r2(2, 2), r3(3, 3), vr(3), &
1285 : xb(3)
1286 :
1287 : ! Variables
1288 : ! ==--------------------------------------------------------------==
1289 : ! Find a point A / V_R = (1-R).OA
1290 :
1291 1208 : DO i = 1, 3
1292 1208 : iok(i) = 0
1293 : END DO
1294 1208 : DO i = 1, 3
1295 1208 : origin(i) = 0._dp
1296 : END DO
1297 2622 : DO ir = 1, nc
1298 2594 : dif = v(1, ir)*v(1, ir) + v(2, ir)*v(2, ir) + v(3, ir)*v(3, ir)
1299 2622 : IF (dif > delta*delta) THEN
1300 2752 : DO i = 1, 3
1301 2752 : igood(i) = 1
1302 : END DO
1303 : ! V is non-zero. Construct matrix 1-R
1304 2752 : DO i = 1, 3
1305 8256 : DO j = 1, 3
1306 8256 : r3(i, j) = -r(i, j, ib(ir))
1307 : END DO
1308 2752 : r3(i, i) = 1 + r3(i, i)
1309 : END DO
1310 688 : CALL invmat(r3, ierror)
1311 688 : IF (ierror == 0) THEN
1312 : ! The matrix 3x3 has an inverse.
1313 464 : DO i = 1, 3
1314 : vr(i) = r3(i, 1)*v(1, ir) &
1315 : + r3(i, 2)*v(2, ir) &
1316 464 : + r3(i, 3)*v(3, ir)
1317 : END DO
1318 : ELSE
1319 : ! IERROR gives the column which causes some trouble
1320 : ! Construct matrix 1-R with 2x2
1321 572 : igood(ierror) = 0
1322 572 : imissing3 = ierror
1323 572 : i1 = 0
1324 2288 : DO i = 1, 3
1325 2288 : IF (i /= ierror) THEN
1326 1144 : i1 = i1 + 1
1327 1144 : j1 = 0
1328 4576 : DO j = 1, 3
1329 4576 : IF (j /= ierror) THEN
1330 2288 : j1 = j1 + 1
1331 2288 : r2(i1, j1) = -r(i, j, ib(ir))
1332 : END IF
1333 : END DO
1334 1144 : r2(i1, i1) = 1 + r2(i1, i1)
1335 : END IF
1336 : END DO
1337 572 : CALL invmat(r2, ierror)
1338 572 : IF (ierror == 0) THEN
1339 : ! The matrix 2X2 has an inverse.
1340 : ! Solve Vxy = (1-R).OAxy + OAz R3z (z is IMISSING3)
1341 : i1 = 0
1342 1936 : DO i = 1, 3
1343 1936 : IF (igood(i) == 1) THEN
1344 968 : i1 = i1 + 1
1345 968 : vr(i) = 0._dp
1346 968 : j1 = 0
1347 3872 : DO j = 1, 3
1348 3872 : IF (igood(j) == 1) THEN
1349 1936 : j1 = j1 + 1
1350 : vr(i) = vr(i) + r2(i1, j1)*(v(j, ir) + &
1351 1936 : origin(imissing3)*r(j, imissing3, ib(ir)))
1352 : END IF
1353 : END DO
1354 : ELSE
1355 484 : vr(i) = origin(i)
1356 : END IF
1357 : END DO
1358 : ELSE
1359 : ! Construct matrix 1-R with 1x1
1360 : i1 = 0
1361 352 : DO i = 1, 3
1362 352 : IF (i /= imissing3) THEN
1363 176 : i1 = i1 + 1
1364 176 : IF (i1 == ierror) THEN
1365 88 : igood(i) = 0
1366 88 : imissing2 = i
1367 : ELSE
1368 : ionly = i
1369 : END IF
1370 : END IF
1371 : END DO
1372 88 : diag = (1 - r(ionly, ionly, ib(ir)))
1373 88 : IF (ABS(diag) > delta) THEN
1374 : vr(ionly) = 1._dp/diag*(v(ionly, ir) + &
1375 : origin(imissing3)*r(ionly, imissing3, ib(ir)) + &
1376 88 : origin(imissing2)*r(ionly, imissing2, ib(ir)))
1377 : ELSE
1378 0 : vr(ionly) = origin(ionly)
1379 0 : igood(ionly) = 0
1380 : END IF
1381 88 : vr(imissing3) = origin(imissing3)
1382 88 : vr(imissing2) = origin(imissing2)
1383 : END IF
1384 : END IF
1385 : ! ==----------------------------------------------------------==
1386 : ! Compare VR with ORIGIN
1387 688 : dif = 0._dp
1388 : ! If NTVEC /=1 there are NTVEC possible standard origins
1389 2752 : DO i = 1, 3
1390 2752 : IF (iok(i) == 1) THEN
1391 860 : dif = dif + ABS(origin(i) - vr(i))
1392 : END IF
1393 : END DO
1394 688 : IF (dif > delta) THEN
1395 : ! Non-symmorphic
1396 274 : info = 0
1397 274 : RETURN
1398 : ELSE
1399 1656 : DO i = 1, 3
1400 1656 : IF (iok(i) /= 1 .AND. igood(i) == 1) THEN
1401 708 : iok(i) = 1
1402 708 : origin(i) = vr(i)
1403 : END IF
1404 : END DO
1405 : END IF
1406 : END IF
1407 : END DO
1408 : ! ==--------------------------------------------------------------==
1409 28 : IF (iok(1) == 0 .AND. iok(2) == 0 .AND. iok(3) == 0) THEN
1410 : ! Cannot not determine
1411 0 : info = -1
1412 0 : RETURN
1413 : END IF
1414 : ! The group is symmorphic
1415 28 : info = 1
1416 : ! Check
1417 84 : DO ir = 1, nc
1418 256 : DO i = 1, 3
1419 : vr(i) = r(i, 1, ib(ir))*origin(1) &
1420 : + r(i, 2, ib(ir))*origin(2) &
1421 192 : + r(i, 3, ib(ir))*origin(3)
1422 256 : vr(i) = (origin(i) - vr(i)) - v(i, ir)
1423 : END DO
1424 64 : CALL rlv3(ai, vr, xb, il, delta)
1425 64 : dif = ABS(xb(1)) + ABS(xb(2)) + ABS(xb(3))
1426 84 : IF (dif > delta) THEN
1427 : ! Non-symmorphic
1428 8 : info = 0
1429 8 : RETURN
1430 : END IF
1431 : END DO
1432 : ! ==--------------------------------------------------------------==
1433 : RETURN
1434 : END SUBROUTINE symmorphic
1435 : ! ==================================================================
1436 : ! **************************************************************************************************
1437 : !> \brief ...
1438 : !> \param ihc ...
1439 : !> \param r ...
1440 : ! **************************************************************************************************
1441 2216 : SUBROUTINE rot1(ihc, r)
1442 : ! ==--------------------------------------------------------------==
1443 : ! == WRITTEN ON FEBRUARY 17TH, 1976 ==
1444 : ! == GENERATION OF THE X,Y,Z-TRANSFORMATION MATRICES 3X3 ==
1445 : ! == FOR HEXAGONAL AND CUBIC GROUPS ==
1446 : ! == SUBROUTINES NEEDED -- NONE ==
1447 : ! ==--------------------------------------------------------------==
1448 : ! == THIS IS IDENTICAL WITH THE SUBROUTINE ROT OF WORLTON-WARREN ==
1449 : ! == (IN THE AC-COMPLEX), ONLY THE WAY OF TRANSFERRING THE DATA ==
1450 : ! == WAS CHANGED ==
1451 : ! ==--------------------------------------------------------------==
1452 : ! == INPUT DATA: ==
1453 : ! == IHC SWITCH DETERMINING IF WE DESIRE ==
1454 : ! == THE HEXAGONAL GROUP(IHC=0) OR THE CUBIC GROUP (IHC=1) ==
1455 : ! == OUTPUT DATA: ==
1456 : ! == R...THE 3X3 MATRICES OF THE DESIRED COORDINATE REPRESENTATION==
1457 : ! == THEIR NUMBERING CORRESPONDS TO THE SYMMETRY ELEMENTS AS ==
1458 : ! == LISTE IN WORLTON-WARREN ==
1459 : ! == (COMPUT. PHYS. COMM. 3(1972) 88--117) ==
1460 : ! == FOR IHC=0 THE FIRST 24 MATRICES OF THE ARRAY R REPRESENT ==
1461 : ! == THE FULL HEXAGONAL GROUP D(6H) ==
1462 : ! == FOR IHC=1 THE FIRST 48 MATRICES OF THE ARRAY R REPRESENT ==
1463 : ! == THE FULL CUBIC GROUP O(H) ==
1464 : ! ==--------------------------------------------------------------==
1465 : INTEGER :: ihc
1466 : REAL(dp) :: r(3, 3, 48)
1467 :
1468 : INTEGER :: i, j, k, n, nv
1469 : REAL(dp) :: c, s
1470 :
1471 8864 : DO j = 1, 3
1472 28808 : DO i = 1, 3
1473 983904 : DO n = 1, 48
1474 977256 : r(i, j, n) = 0._dp
1475 : END DO
1476 : END DO
1477 : END DO
1478 2216 : IF (ihc == 0) THEN
1479 : ! ==------------------------------------------------------------==
1480 : ! DEFINE THE GENERATORS FOR THE ROTATION MATRICES--HEXAGONAL GROUP
1481 : ! ==------------------------------------------------------------==
1482 1144 : c = 0.5_dp
1483 1144 : s = 0.5_dp*SQRT(3.0_dp)
1484 1144 : r(1, 1, 2) = c
1485 1144 : r(1, 2, 2) = -s
1486 1144 : r(2, 1, 2) = s
1487 1144 : r(2, 2, 2) = c
1488 1144 : r(1, 1, 7) = -c
1489 1144 : r(1, 2, 7) = -s
1490 1144 : r(2, 1, 7) = -s
1491 1144 : r(2, 2, 7) = c
1492 8008 : DO n = 1, 6
1493 6864 : r(3, 3, n) = 1._dp
1494 6864 : r(3, 3, n + 18) = 1._dp
1495 6864 : r(3, 3, n + 6) = -1._dp
1496 8008 : r(3, 3, n + 12) = -1._dp
1497 : END DO
1498 : ! ==------------------------------------------------------------==
1499 : ! == GENERATE THE REST OF THE ROTATION MATRICES ==
1500 : ! ==------------------------------------------------------------==
1501 3432 : DO i = 1, 2
1502 2288 : r(i, i, 1) = 1._dp
1503 8008 : DO j = 1, 2
1504 4576 : r(i, j, 6) = r(j, i, 2)
1505 16016 : DO k = 1, 2
1506 9152 : r(i, j, 3) = r(i, j, 3) + r(i, k, 2)*r(k, j, 2)
1507 9152 : r(i, j, 8) = r(i, j, 8) + r(i, k, 2)*r(k, j, 7)
1508 13728 : r(i, j, 12) = r(i, j, 12) + r(i, k, 7)*r(k, j, 2)
1509 : END DO
1510 : END DO
1511 : END DO
1512 3432 : DO i = 1, 2
1513 8008 : DO j = 1, 2
1514 4576 : r(i, j, 5) = r(j, i, 3)
1515 16016 : DO k = 1, 2
1516 9152 : r(i, j, 4) = r(i, j, 4) + r(i, k, 2)*r(k, j, 3)
1517 9152 : r(i, j, 9) = r(i, j, 9) + r(i, k, 2)*r(k, j, 8)
1518 9152 : r(i, j, 10) = r(i, j, 10) + r(i, k, 12)*r(k, j, 3)
1519 13728 : r(i, j, 11) = r(i, j, 11) + r(i, k, 12)*r(k, j, 2)
1520 : END DO
1521 : END DO
1522 : END DO
1523 14872 : DO n = 1, 12
1524 13728 : nv = n + 12
1525 42328 : DO i = 1, 2
1526 96096 : DO j = 1, 2
1527 82368 : r(i, j, nv) = -r(i, j, n)
1528 : END DO
1529 : END DO
1530 : END DO
1531 : ELSE
1532 : ! ==------------------------------------------------------------==
1533 : ! == DEFINE THE GENERATORS FOR THE ROTATION MATRICES-CUBIC GROUP==
1534 : ! ==------------------------------------------------------------==
1535 1072 : r(1, 3, 9) = 1._dp
1536 1072 : r(2, 1, 9) = 1._dp
1537 1072 : r(3, 2, 9) = 1._dp
1538 1072 : r(1, 1, 19) = 1._dp
1539 1072 : r(2, 3, 19) = -1._dp
1540 1072 : r(3, 2, 19) = 1._dp
1541 4288 : DO i = 1, 3
1542 3216 : r(i, i, 1) = 1._dp
1543 13936 : DO j = 1, 3
1544 9648 : r(i, j, 20) = r(j, i, 19)
1545 9648 : r(i, j, 5) = r(j, i, 9)
1546 41808 : DO k = 1, 3
1547 28944 : r(i, j, 2) = r(i, j, 2) + r(i, k, 19)*r(k, j, 19)
1548 28944 : r(i, j, 16) = r(i, j, 16) + r(i, k, 9)*r(k, j, 19)
1549 38592 : r(i, j, 23) = r(i, j, 23) + r(i, k, 19)*r(k, j, 9)
1550 : END DO
1551 : END DO
1552 : END DO
1553 4288 : DO i = 1, 3
1554 13936 : DO j = 1, 3
1555 41808 : DO k = 1, 3
1556 28944 : r(i, j, 6) = r(i, j, 6) + r(i, k, 2)*r(k, j, 5)
1557 28944 : r(i, j, 7) = r(i, j, 7) + r(i, k, 16)*r(k, j, 23)
1558 28944 : r(i, j, 8) = r(i, j, 8) + r(i, k, 5)*r(k, j, 2)
1559 28944 : r(i, j, 10) = r(i, j, 10) + r(i, k, 2)*r(k, j, 9)
1560 28944 : r(i, j, 11) = r(i, j, 11) + r(i, k, 9)*r(k, j, 2)
1561 28944 : r(i, j, 12) = r(i, j, 12) + r(i, k, 23)*r(k, j, 16)
1562 28944 : r(i, j, 14) = r(i, j, 14) + r(i, k, 16)*r(k, j, 2)
1563 28944 : r(i, j, 15) = r(i, j, 15) + r(i, k, 2)*r(k, j, 16)
1564 28944 : r(i, j, 22) = r(i, j, 22) + r(i, k, 23)*r(k, j, 2)
1565 38592 : r(i, j, 24) = r(i, j, 24) + r(i, k, 2)*r(k, j, 23)
1566 : END DO
1567 : END DO
1568 : END DO
1569 4288 : DO i = 1, 3
1570 13936 : DO j = 1, 3
1571 41808 : DO k = 1, 3
1572 28944 : r(i, j, 3) = r(i, j, 3) + r(i, k, 5)*r(k, j, 12)
1573 28944 : r(i, j, 4) = r(i, j, 4) + r(i, k, 5)*r(k, j, 10)
1574 28944 : r(i, j, 13) = r(i, j, 13) + r(i, k, 23)*r(k, j, 11)
1575 28944 : r(i, j, 17) = r(i, j, 17) + r(i, k, 16)*r(k, j, 12)
1576 28944 : r(i, j, 18) = r(i, j, 18) + r(i, k, 16)*r(k, j, 10)
1577 38592 : r(i, j, 21) = r(i, j, 21) + r(i, k, 12)*r(k, j, 15)
1578 : END DO
1579 : END DO
1580 : END DO
1581 26800 : DO n = 1, 24
1582 25728 : nv = n + 24
1583 25728 : r(1, 1, nv) = -r(1, 1, n)
1584 25728 : r(1, 2, nv) = -r(1, 2, n)
1585 25728 : r(1, 3, nv) = -r(1, 3, n)
1586 25728 : r(2, 1, nv) = -r(2, 1, n)
1587 25728 : r(2, 2, nv) = -r(2, 2, n)
1588 25728 : r(2, 3, nv) = -r(2, 3, n)
1589 25728 : r(3, 1, nv) = -r(3, 1, n)
1590 25728 : r(3, 2, nv) = -r(3, 2, n)
1591 26800 : r(3, 3, nv) = -r(3, 3, n)
1592 : END DO
1593 : END IF
1594 : ! ==--------------------------------------------------------------==
1595 2216 : RETURN
1596 : END SUBROUTINE rot1
1597 : ! ==================================================================
1598 : END MODULE kpsym
|
