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 K-points and crystal symmetry routines based on
10 : ! K290 code:
11 : ! Written on September 12th, 1979.
12 : ! IBM-retouched on October 27th, 1980.
13 : ! Generation of special points modified on 26-May-82 by ohn.
14 : ! Retouched on January 8th, 1997
15 : ! Integration in CPMD-FEMD Program by Thierry Deutsch
16 : ! ==--------------------------------------------------------------==
17 : ! Playing with special points and creation of 'CRYSTALLOGRAPHIC'
18 : ! File for band structure calculations.
19 : ! Generation of special points in k-space for an arbitrary lattice,
20 : ! Following the method Monkhorst,Pack, Phys. Rev. B13 (1976) 5188
21 : ! Modified by Macdonald, Phys. Rev. B18 (1978) 5897
22 : ! Modified also by Ole Holm Nielsen ("SYMMETRIZATION")
23 : ! ==--------------------------------------------------------------==
24 : ! (GROUP1, PGL1, ATFTM1, ROT1 FROM THE
25 : ! "COMPUTER PHYSICS COMMUNICATIONS" PACKAGE "ACMI" - (1971,1974)
26 : ! Worlton-Warren).
27 : ! **************************************************************************************************
28 : MODULE kpsym
29 :
30 : USE kinds, ONLY: dp
31 : USE mathlib, ONLY: invmat
32 : USE string_utilities, ONLY: xstring
33 : #include "./base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 : PRIVATE
37 :
38 : PUBLIC :: K290s, GROUP1s
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpsym'
41 :
42 : ! **************************************************************************************************
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief ...
48 : !> \param iout ...
49 : !> \param nat ...
50 : !> \param nkpoint ...
51 : !> \param nsp ...
52 : !> \param iq1 ...
53 : !> \param iq2 ...
54 : !> \param iq3 ...
55 : !> \param istriz ...
56 : !> \param a1 ...
57 : !> \param a2 ...
58 : !> \param a3 ...
59 : !> \param alat ...
60 : !> \param strain ...
61 : !> \param xkapa ...
62 : !> \param rx ...
63 : !> \param tvec ...
64 : !> \param ty ...
65 : !> \param isc ...
66 : !> \param f0 ...
67 : !> \param ntvec ...
68 : !> \param wvk0 ...
69 : !> \param wvkl ...
70 : !> \param lwght ...
71 : !> \param lrot ...
72 : !> \param nhash ...
73 : !> \param includ ...
74 : !> \param list ...
75 : !> \param rlist ...
76 : !> \param delta ...
77 : ! **************************************************************************************************
78 960 : SUBROUTINE k290s(iout, nat, nkpoint, nsp, iq1, iq2, iq3, istriz, &
79 960 : a1, a2, a3, alat, strain, xkapa, rx, tvec, &
80 960 : ty, isc, f0, ntvec, wvk0, wvkl, lwght, lrot, &
81 960 : nhash, includ, list, rlist, delta)
82 : ! ==================================================================
83 : ! WRITTEN ON SEPTEMBER 12TH, 1979.
84 : ! IBM-RETOUCHED ON OCTOBER 27TH, 1980.
85 : ! Tsukuba-retouched on March 19th, 2008.
86 : ! GENERATION OF SPECIAL POINTS MODIFIED ON 26-MAY-82 BY OHN.
87 : ! RETOUCHED ON JANUARY 8TH, 1997
88 : ! INTEGRATION IN CPMD-FEMD PROGRAM BY THIERRY DEUTSCH
89 : ! ==--------------------------------------------------------------==
90 : ! PLAYING WITH SPECIAL POINTS AND CREATION OF 'CRYSTALLOGRAPHIC'
91 : ! FILE FOR BAND STRUCTURE CALCULATIONS.
92 : ! GENERATION OF SPECIAL POINTS IN K-SPACE FOR AN ARBITRARY LATTICE,
93 : ! FOLLOWING THE METHOD MONKHORST,PACK, PHYS. REV. B13 (1976) 5188
94 : ! MODIFIED BY MACDONALD, PHYS. REV. B18 (1978) 5897
95 : ! MODIFIED ALSO BY OLE HOLM NIELSEN ("SYMMETRIZATION")
96 : ! ==--------------------------------------------------------------==
97 : ! TESTING THEIR EFFICIENCY AND PREPARATION OF THE
98 : ! "STRUCTURAL" FILE FOR RUNNING THE
99 : ! SELF-CONSISTENT BAND STRUCTURE PROGRAMS.
100 : ! IN THE CASES WHERE THE POINT GROUP OF THE CRYSTAL DOES NOT
101 : ! CONTAIN INVERSION, THE LATTER IS ARTIFICIALLY ADDED, IN ORDER
102 : ! TO MAKE USE OF THE HERMITICITY OF THE HAMILTONIAN
103 : ! ==--------------------------------------------------------------==
104 : ! == INPUT: ==
105 : ! == IOUT LOGIC FILE NUMBER ==
106 : ! == NAT NUMBER OF ATOMS ==
107 : ! == NKPOINT MAXIMAL NUMBER OF K POINTS ==
108 : ! == NSP NUMBER OF SPECIES ==
109 : ! == IQ1,IQ2,IQ3 THE MONKHORST-PACK MESH PARAMETERS ==
110 : ! == ISTRIZ SWITCH FOR SYMMETRIZATION ==
111 : ! == A1(3),A2(3),A3(3) LATTICE VECTORS ==
112 : ! == ALAT LATTICE CONSTANT ==
113 : ! == STRAIN(3,3) STRAIN APPLIED TO LATTICE IN ORDER ==
114 : ! == TO HAVE K POINTS WITH SYMMETRY OF STRAINED LATTICE ==
115 : ! == XKAPA(3,NAT) ATOMS COORDINATES ==
116 : ! == TY(NAT) TYPES OF ATOMS ==
117 : ! == WVK0(3) SHIFT FOR K POINTS MESh (MACDONALD ARTICLE) ==
118 : ! == NHASH SIZE OF THE HASH TABLES (LIST) ==
119 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
120 : ! == K-VECTOR < DELTA IS CONSIDERED ZERO ==
121 : ! == OUTPUT: ==
122 : ! == RX(3,NAT) SCRATCH ARRAY USED BY GROUP1 ROUTINE ==
123 : ! == TVEC(1:3,1:NTVEC) TRANSLATION VECTORS (SEE NTVEC) ==
124 : ! == ISC(NAT) SCRATCH ARRAY USED BY GROUP1 ROUTINE ==
125 : ! == F0(49,NAT) ATOM TRANSFORMATION TABLE ==
126 : ! == IF NTVEC/=1 THE 49TH GIVES INEQUIVALENT ATOMS ==
127 : ! == NTVEC NUMBER OF TRANSLATION VECTORS (IF NOT PRIMITIVE CELL)==
128 : ! == WVKL(3,NKPOINT) SPECIAL KPOINTS GENERATED ==
129 : ! == LWGHT(NKPOINT) WEIGHT FOR EACH K POINT ==
130 : ! == LROT(48,NKPOINT) SYMMETRY OPERATION FOR EACH K POINTS ==
131 : ! == INCLUD(NKPOINT) SCRATCH ARRAY USED BY SPPT2 ==
132 : ! == LIST(NKPOINT+NHASH) HASH TABLE USED BY SPPT2 ==
133 : ! == RLIST(3,NKPOINT) SCRATCH ARRAY USED BY SPPT2 ==
134 : ! ==--------------------------------------------------------------==
135 : ! SUBROUTINES NEEDED:
136 : ! SPPT2, GROUP1, PGL1, ATFTM1, ROT1, STRUCT,
137 : ! BZRDUC, INBZ, MESH, BZDEFI
138 : ! (GROUP1, PGL1, ATFTM1, ROT1 FROM THE
139 : ! "COMPUTER PHYSICS COMMUNICATIONS" PACKAGE "ACMI" - (1971,1974)
140 : ! WORLTON-WARREN).
141 : ! ==================================================================
142 : INTEGER :: iout, nat, nkpoint, nsp, iq1, iq2, iq3, &
143 : istriz
144 : REAL(KIND=dp) :: a1(3), a2(3), a3(3), alat, strain(6), &
145 : xkapa(3, nat), rx(3, nat), tvec(3, nat)
146 : INTEGER :: ty(nat), isc(nat), f0(49, nat), ntvec
147 : REAL(KIND=dp) :: wvk0(3), wvkl(3, nkpoint)
148 : INTEGER :: lwght(nkpoint), lrot(48, nkpoint), &
149 : nhash, includ(nkpoint), &
150 : list(nkpoint + nhash)
151 : REAL(KIND=dp) :: rlist(3, nkpoint), delta
152 :
153 : CHARACTER(len=10), DIMENSION(48), PARAMETER :: rname_cubic = [' 1 ', ' 2[ 10 0] ', &
154 : ' 2[ 01 0] ', ' 2[ 00 1] ', ' 3[-1-1-1]', ' 3[ 11-1] ', ' 3[-11 1] ', ' 3[ 1-11] ', &
155 : ' 3[ 11 1] ', ' 3[-11-1] ', ' 3[-1-11] ', ' 3[ 1-1-1]', ' 2[-11 0] ', ' 4[ 00 1] ', &
156 : ' 4[ 00-1] ', ' 2[ 11 0] ', ' 2[ 0-11] ', ' 2[ 01 1] ', ' 4[ 10 0] ', ' 4[-10 0] ', &
157 : ' 2[-10 1] ', ' 4[ 0-10] ', ' 2[ 10 1] ', ' 4[ 01 0] ', '-1 ', '-2[ 10 0] ', &
158 : '-2[ 01 0] ', '-2[ 00 1] ', '-3[-1-1-1]', '-3[ 11-1] ', '-3[-11 1] ', '-3[ 1-11] ', &
159 : '-3[ 11 1] ', '-3[-11-1] ', '-3[-1-11] ', '-3[ 1-1-1]', '-2[-11 0] ', '-4[ 00 1] ', &
160 : '-4[ 00-1] ', '-2[ 11 0] ', '-2[ 0-11] ', '-2[ 01 1] ', '-4[ 10 0] ', '-4[-10 0] ', &
161 : '-2[-10 1] ', '-4[ 0-10] ', '-2[ 10 1] ', '-4[ 01 0] ']
162 : CHARACTER(len=11), DIMENSION(24), PARAMETER :: rname_hexai = [' 1 ', ' 6[ 00 1] ', &
163 : ' 3[ 00 1] ', ' 2[ 00 1] ', ' 3[ 00 -1] ', ' 6[ 00 -1] ', ' 2[ 01 0] ', ' 2[-11 0] ', &
164 : ' 2[ 10 0] ', ' 2[ 21 0] ', ' 2[ 11 0] ', ' 2[ 12 0] ', '-1 ', '-6[ 00 1] ', &
165 : '-3[ 00 1] ', '-2[ 00 1] ', '-3[ 00 -1] ', '-6[ 00 -1] ', '-2[ 01 0] ', '-2[-11 0] ', &
166 : '-2[ 10 0] ', '-2[ 21 0] ', '-2[ 11 0] ', '-2[ 12 0] ']
167 : CHARACTER(len=12), DIMENSION(7), PARAMETER :: icst = ['TRICLINIC ', 'MONOCLINIC ', &
168 : 'ORTHORHOMBIC', 'TETRAGONAL ', 'CUBIC ', 'TRIGONAL ', 'HEXAGONAL ']
169 :
170 : INTEGER :: i, ib(48), ib0(48), ihc, ihc0, ihg, ihg0, indpg, indpg0, invadd, istrin, iswght, &
171 : isy, isy0, itype, j, k, l, li, li0, lmax, n, nc, nc0, ntot, ntvec0
172 : INTEGER, DIMENSION(49, 1) :: f00
173 : LOGICAL :: located_type
174 : REAL(KIND=dp) :: a01(3), a02(3), a03(3), b01(3), b02(3), b03(3), b1(3), b2(3), b3(3), &
175 : dtotstr, origin(3), origin0(3), proj1, proj2, proj3, r(3, 3, 48), r0(3, 3, 48), totstr, &
176 : tvec0(3, 1), volum, vv0(3)
177 : REAL(KIND=dp), DIMENSION(3, 1) :: x0
178 : REAL(KIND=dp), DIMENSION(3, 48) :: v, v0
179 :
180 960 : f00 = 0
181 960 : x0 = 0._dp
182 960 : v = 0._dp
183 960 : v0 = 0._dp
184 : ! ==--------------------------------------------------------------==
185 : ! READ IN LATTICE STRUCTURE
186 : ! ==--------------------------------------------------------------==
187 3840 : DO i = 1, 3
188 2880 : a01(i) = a1(i)/alat
189 2880 : a02(i) = a2(i)/alat
190 3840 : a03(i) = a3(i)/alat
191 : END DO
192 960 : IF (iout > 0) THEN
193 335 : WRITE (iout, '(" KPSYM| NUMBER OF ATOMS (STRUCT):",I6)') nat
194 : END IF
195 960 : IF (iout > 0) THEN
196 335 : WRITE (iout, '(" KPSYM|",10X,"K TYPE",14X,"X(K)")')
197 : END IF
198 960 : itype = 0
199 6062 : DO i = 1, nat
200 : ! Assign an atomic type (for internal purposes)
201 5102 : located_type = .FALSE.
202 5102 : IF (i /= 1) THEN
203 9294 : DO j = 1, (i - 1)
204 9294 : IF (ty(j) == ty(i)) THEN
205 : ! Type located
206 : located_type = .TRUE.
207 : EXIT
208 : END IF
209 : END DO
210 : ! New type
211 : END IF
212 4142 : IF (.NOT. located_type) THEN
213 1190 : itype = itype + 1
214 1190 : IF (itype > nsp) THEN
215 0 : IF (iout > 0) THEN
216 : WRITE (iout, '(A,I4,")")') &
217 0 : ' KPSYM| NUMBER OF ATOMIC TYPES EXCEEDS DIMENSION (NSP=)', &
218 0 : nsp
219 : END IF
220 0 : IF (iout > 0) THEN
221 : WRITE (iout, '(" KPSYM| THE ARRAY TY IS:",/,9(1X,10I7,/))') &
222 0 : (ty(j), j=1, nat)
223 : END IF
224 0 : CPABORT('K290: FATAL ERROR')
225 : END IF
226 : END IF
227 6062 : IF (iout > 0) THEN
228 : WRITE (iout, '(" KPSYM|",6X,I5,I6,3F10.5)') &
229 1733 : i, ty(i), (xkapa(j, i), j=1, 3)
230 : END IF
231 : END DO
232 : ! ==--------------------------------------------------------------==
233 : ! IS THE STRAIN SIGNIFICANT ?
234 : ! ==--------------------------------------------------------------==
235 960 : dtotstr = delta*delta
236 960 : totstr = 0._dp
237 960 : istrin = 0
238 6720 : DO i = 1, 6
239 6720 : totstr = totstr + ABS(strain(i))
240 : END DO
241 960 : IF (totstr > dtotstr) istrin = 1
242 : ! ==--------------------------------------------------------------==
243 : ! Volume of the cell.
244 : volum = a1(1)*a2(2)*a3(3) + a2(1)*a3(2)*a1(3) + &
245 : a3(1)*a1(2)*a2(3) - a1(3)*a2(2)*a3(1) - &
246 960 : A2(3)*A3(2)*A1(1) - A3(3)*A1(2)*A2(1)
247 960 : volum = ABS(volum)
248 960 : b1(1) = (a2(2)*a3(3) - a2(3)*a3(2))/volum
249 960 : b1(2) = (a2(3)*a3(1) - a2(1)*a3(3))/volum
250 960 : b1(3) = (a2(1)*a3(2) - a2(2)*a3(1))/volum
251 960 : b2(1) = (a3(2)*a1(3) - a3(3)*a1(2))/volum
252 960 : b2(2) = (a3(3)*a1(1) - a3(1)*a1(3))/volum
253 960 : b2(3) = (a3(1)*a1(2) - a3(2)*a1(1))/volum
254 960 : b3(1) = (a1(2)*a2(3) - a1(3)*a2(2))/volum
255 960 : b3(2) = (a1(3)*a2(1) - a1(1)*a2(3))/volum
256 960 : b3(3) = (a1(1)*a2(2) - a1(2)*a2(1))/volum
257 : ! ==--------------------------------------------------------------==
258 3840 : DO i = 1, 3
259 2880 : b01(i) = b1(i)*alat
260 2880 : b02(i) = b2(i)*alat
261 3840 : b03(i) = b3(i)*alat
262 : END DO
263 : ! ==--------------------------------------------------------------==
264 : ! == GROUP-THEORY ANALYSIS OF LATTICE ==
265 : ! ==--------------------------------------------------------------==
266 : CALL group1s(iout, a1, a2, a3, nat, ty, xkapa, b1, b2, b3, &
267 : ihg, ihc, isy, li, nc, indpg, ib, ntvec, &
268 960 : v, f0, r, tvec, origin, rx, isc, delta)
269 : ! ==--------------------------------------------------------------==
270 31134 : DO n = nc + 1, 48
271 31134 : ib(n) = 0
272 : END DO
273 : ! ==--------------------------------------------------------------==
274 960 : invadd = 0
275 960 : IF (li == 0) THEN
276 466 : IF (iout > 0) THEN
277 : WRITE (iout, '(A,/,A,/,A)') &
278 189 : ' KPSYM| ALTHOUGH THE POINT GROUP OF THE CRYSTAL DOES NOT', &
279 189 : ' KPSYM| CONTAIN INVERSION, THE SPECIAL POINT GENERATION ALGORITHM', &
280 378 : ' KPSYM| WILL CONSIDER IT AS A SYMMETRY OPERATION'
281 : END IF
282 466 : invadd = 1
283 : END IF
284 : ! ==--------------------------------------------------------------==
285 : ! == CRYSTALLOGRAPHIC DATA ==
286 : ! ==--------------------------------------------------------------==
287 960 : IF (iout > 0) THEN
288 335 : WRITE (iout, '(/," KPSYM| CRYSTALLOGRAPHIC DATA:")')
289 335 : WRITE (iout, '(4X,"A1",3F10.5,10X,"B1",3F10.5)') a1, b1
290 335 : WRITE (iout, '(4X,"A2",3F10.5,10X,"B2",3F10.5)') a2, b2
291 335 : WRITE (iout, '(4X,"A3",3F10.5,10X,"B3",3F10.5)') a3, b3
292 : END IF
293 : ! ==--------------------------------------------------------------==
294 : ! == GROUP-THEORETICAL INFORMATION ==
295 : ! ==--------------------------------------------------------------==
296 960 : IF (iout > 0) THEN
297 335 : WRITE (iout, '(/," KPSYM| GROUP-THEORETICAL INFORMATION:")')
298 : END IF
299 : ! IHG .... Point group of the primitive lattice, holohedral
300 960 : IF (iout > 0) THEN
301 : WRITE (iout, &
302 : '(" KPSYM| POINT GROUP OF THE PRIMITIVE LATTICE: ",A," SYSTEM")') &
303 335 : icst(ihg)
304 : END IF
305 : ! IHC .... Code distinguishing between hexagonal and cubic groups
306 : ! ISY .... Code indicating whether the space group is symmorphic
307 960 : IF (isy == 0) THEN
308 288 : IF (iout > 0) THEN
309 96 : WRITE (iout, '(" KPSYM|",4X,"NONSYMMORPHIC GROUP")')
310 : END IF
311 672 : ELSE IF (isy == 1) THEN
312 652 : IF (iout > 0) THEN
313 231 : WRITE (iout, '(" KPSYM|",4X,"SYMMORPHIC GROUP")')
314 : END IF
315 20 : ELSE IF (isy == -1) THEN
316 20 : IF (iout > 0) THEN
317 8 : WRITE (iout, '(" KPSYM|",4X,"SYMMORPHIC GROUP WITH NON-STANDARD ORIGIN")')
318 : END IF
319 0 : ELSE IF (isy == -2) THEN
320 0 : IF (iout > 0) THEN
321 0 : WRITE (iout, '(" KPSYM|",4X,"NONSYMMORPHIC GROUP???")')
322 : END IF
323 : END IF
324 : ! LI ..... Inversions symmetry
325 960 : IF (li == 0) THEN
326 466 : IF (iout > 0) THEN
327 189 : WRITE (iout, '(" KPSYM|",4X,"NO INVERSION SYMMETRY")')
328 : END IF
329 494 : ELSE IF (li > 0) THEN
330 494 : IF (iout > 0) THEN
331 146 : WRITE (iout, '(" KPSYM|",4X,"INVERSION SYMMETRY")')
332 : END IF
333 : END IF
334 : ! NC ..... Total number of elements in the point group
335 960 : IF (iout > 0) THEN
336 : WRITE (iout, &
337 335 : '(" KPSYM|",4X,"TOTAL NUMBER OF ELEMENTS IN THE POINT GROUP:",I3)') nc
338 : END IF
339 960 : IF (iout > 0) THEN
340 : WRITE (iout, '(" KPSYM|",4X,"TO SUM UP: (",I1,5I3,")")') &
341 335 : ihg, ihc, isy, li, nc, indpg
342 : END IF
343 : ! IB ..... List of the rotations constituting the point group
344 960 : IF (iout > 0) THEN
345 335 : WRITE (iout, '(/," KPSYM|",4X,"LIST OF THE ROTATIONS:")')
346 : END IF
347 960 : IF (iout > 0) THEN
348 335 : WRITE (iout, '(7X,12I4)') (ib(i), i=1, nc)
349 : END IF
350 : ! V ...... Nonprimitive translations (for nonsymmorphic groups)
351 960 : IF (isy <= 0) THEN
352 308 : IF (iout > 0) THEN
353 104 : WRITE (iout, '(/," KPSYM|",4X,"NONPRIMITIVE TRANSLATIONS:")')
354 : END IF
355 308 : IF (iout > 0) THEN
356 : WRITE (iout, '(A,A)') &
357 104 : ' ROT V IN THE BASIS A1, A2, A3 ', &
358 208 : 'V IN CARTESIAN COORDINATES'
359 : END IF
360 : ! Cartesian components of nonprimitive translation.
361 7926 : DO i = 1, nc
362 30472 : DO j = 1, 3
363 30472 : vv0(j) = v(1, i)*a1(j) + v(2, i)*a2(j) + v(3, i)*a3(j)
364 : END DO
365 7926 : IF (iout > 0) THEN
366 : WRITE (iout, '(1X,I3,3F10.5,3X,3F10.5)') &
367 2148 : ib(i), (v(j, i), j=1, 3), vv0
368 : END IF
369 : END DO
370 : END IF
371 : ! F0 ..... The function defined in Maradudin, Ipatova by
372 : ! eq. (3.2.12): atom transformation table.
373 960 : IF (iout > 0) THEN
374 : WRITE (iout, &
375 335 : '(/," KPSYM|",4X,"ATOM TRANSFORMATION TABLE (MARADUDIN,VOSKO):")')
376 : END IF
377 960 : IF (iout > 0) THEN
378 335 : WRITE (iout, '(5(4X,"R AT->AT"))')
379 : END IF
380 960 : IF (iout > 0) THEN
381 335 : WRITE (iout, '(I5," [Identity]")') 1
382 : END IF
383 15906 : DO k = 2, nc
384 84360 : DO j = 1, nat
385 69414 : IF (iout > 0) THEN
386 18819 : WRITE (iout, '(I5,2I4)', advance="no") ib(k), j, f0(k, j)
387 : END IF
388 84360 : IF ((MOD(j, 5) == 0) .AND. iout > 0) THEN
389 2060 : WRITE (iout, *)
390 : END IF
391 : END DO
392 15906 : IF ((MOD(j - 1, 5) /= 0) .AND. iout > 0) THEN
393 4019 : WRITE (iout, *)
394 : END IF
395 : END DO
396 : ! R ...... List of the 3 x 3 rotation matrices
397 960 : IF (iout > 0) THEN
398 335 : WRITE (iout, '(/," KPSYM|",4X,"LIST OF THE 3 X 3 ROTATION MATRICES:")')
399 : END IF
400 960 : IF (ihc == 0) THEN
401 874 : DO k = 1, nc
402 874 : IF (iout > 0) THEN
403 : WRITE (iout, &
404 : '(4X,I3," (",I2,": ",A11,")",2(3F14.6,/,25X),3F14.6)') &
405 1443 : k, ib(k), rname_hexai(ib(k)), ((r(i, j, ib(k)), j=1, 3), i=1, 3)
406 : END IF
407 : END DO
408 : ELSE
409 15992 : DO k = 1, nc
410 15992 : IF (iout > 0) THEN
411 : WRITE (iout, &
412 : '(4X,I3," (",I2,": ",A10,") ",2(3F14.6,/,25X),3F14.6)') &
413 55159 : k, ib(k), rname_cubic(ib(k)), ((r(i, j, ib(k)), j=1, 3), i=1, 3)
414 : END IF
415 : END DO
416 : END IF
417 : ! ==--------------------------------------------------------------==
418 : ! GENERATE THE BRAVAIS LATTICE
419 : ! ==--------------------------------------------------------------==
420 : CALL group1s(iout, a01, a02, a03, 1, ty, x0, b01, b02, b03, &
421 : ihg0, ihc0, isy0, li0, nc0, indpg0, ib0, ntvec0, &
422 960 : v0, f00, r0, tvec0, origin0, rx, isc, delta)
423 : ! ==--------------------------------------------------------------==
424 : ! It is assumed that the same 'type' of symmetry operations
425 : ! (cubic/hexagonal) will apply to the crystal as well as the Bravais
426 : ! lattice.
427 : ! ==--------------------------------------------------------------==
428 960 : IF (iout > 0) THEN
429 : WRITE (iout, '(/,1X,19("*"),A,25("*"))') &
430 335 : ' GENERATION OF SPECIAL POINTS '
431 : END IF
432 : ! Parameter Q of Monkhorst and Pack, generalized for 3 axes B1,2,3
433 960 : IF (iout > 0) THEN
434 : WRITE (iout, '(A,/,1X,3I5)') &
435 335 : ' KPSYM| MONKHORST-PACK PARAMETERS (GENERALIZED) IQ1,IQ2,IQ3:', &
436 670 : iq1, iq2, iq3
437 : END IF
438 : ! WVK0 is the shift of the whole mesh (see Macdonald)
439 960 : IF (iout > 0) THEN
440 : WRITE (iout, '(A,/,1X,3F10.5)') &
441 335 : ' KPSYM| CONSTANT VECTOR SHIFT (MACDONALD) OF THIS MESH:', wvk0
442 : END IF
443 960 : IF (ABS(iq1) + ABS(iq2) + ABS(iq3) == 0) RETURN
444 960 : IF (ABS(istriz) /= 1) THEN
445 0 : IF (iout > 0) THEN
446 0 : WRITE (iout, '(" KPSYM| INVALID SWITCH FOR SYMMETRIZATION",I10)') istriz
447 : END IF
448 0 : IF (iout > 0) THEN
449 0 : WRITE (iout, '(" KPSYM| INVALID SWITCH FOR SYMMETRIZATION",I10)') istriz
450 : END IF
451 0 : CPABORT('K290. ISTRIZ WRONG ARGUMENT')
452 : END IF
453 960 : IF (iout > 0) THEN
454 335 : WRITE (iout, '(" KPSYM| SYMMETRIZATION SWITCH: ",I3)', advance="no") istriz
455 : END IF
456 960 : IF (istriz == 1) THEN
457 960 : IF (iout > 0) THEN
458 335 : WRITE (iout, '(" (SYMMETRIZATION OF MONKHORST-PACK MESH)")')
459 : END IF
460 : ELSE
461 0 : IF (iout > 0) THEN
462 0 : WRITE (iout, '(" (NO SYMMETRIZATION OF MONKHORST-PACK MESH)")')
463 : END IF
464 : END IF
465 : ! Set to 0.
466 135820 : DO i = 1, nkpoint
467 135820 : lwght(i) = 0
468 : END DO
469 : ! ==--------------------------------------------------------------==
470 : ! == Generation of the points (they are not multiplied ==
471 : ! == by 2*Pi because B1,2,3 were not,either) ==
472 : ! ==--------------------------------------------------------------==
473 960 : IF (nc > nc0) THEN
474 : ! Due to non-use of primitive cell, the crystal has more
475 : ! rotations than Bravais lattice.
476 : ! We use only the rotations for Bravais lattices
477 0 : IF (ntvec == 1) THEN
478 0 : IF (iout > 0) THEN
479 0 : WRITE (iout, *) ' KPSYM| NUMBER OF ROTATIONS FOR BRAVAIS LATTICE', nc0
480 : END IF
481 0 : IF (iout > 0) THEN
482 0 : WRITE (iout, *) ' KPSYM| NUMBER OF ROTATIONS FOR CRYSTAL LATTICE', nc
483 : END IF
484 0 : IF (iout > 0) THEN
485 0 : WRITE (iout, *) ' KPSYM| NO DUPLICATION FOUND'
486 : END IF
487 0 : CPABORT('SOMETHING IS WRONG IN GROUP DETERMINATION')
488 : END IF
489 0 : nc = nc0
490 0 : DO i = 1, nc0
491 0 : ib(i) = ib0(i)
492 : END DO
493 0 : IF (iout > 0) THEN
494 0 : WRITE (iout, '(/,1X,20("! "),"WARNING",20("!"))')
495 : END IF
496 0 : IF (iout > 0) THEN
497 : WRITE (iout, '(A)') &
498 0 : ' KPSYM| THE CRYSTAL HAS MORE SYMMETRY THAN THE BRAVAIS LATTICE'
499 : END IF
500 0 : IF (iout > 0) THEN
501 : WRITE (iout, '(A)') &
502 0 : ' KPSYM| BECAUSE THIS IS NOT A PRIMITIVE CELL'
503 : END IF
504 0 : IF (iout > 0) THEN
505 : WRITE (iout, '(A)') &
506 0 : ' KPSYM| USE ONLY SYMMETRY FROM BRAVAIS LATTICE'
507 : END IF
508 0 : IF (iout > 0) THEN
509 0 : WRITE (iout, '(1X,20("! "),"WARNING",20("!"),/)')
510 : END IF
511 : END IF
512 : CALL sppt2(iout, iq1, iq2, iq3, wvk0, nkpoint, &
513 : a01, a02, a03, b01, b02, b03, &
514 : invadd, nc, ib, r, ntot, wvkl, lwght, lrot, nc0, ib0, istriz, &
515 960 : nhash, includ, list, rlist, delta)
516 : ! ==--------------------------------------------------------------==
517 : ! == Check on error signals ==
518 : ! ==--------------------------------------------------------------==
519 960 : IF (iout > 0) THEN
520 335 : WRITE (iout, '(/," KPSYM|",1X,I5," SPECIAL POINTS GENERATED")') ntot
521 : END IF
522 960 : IF (ntot == 0) THEN
523 : RETURN
524 960 : ELSE IF (ntot < 0) THEN
525 0 : IF (iout > 0) THEN
526 0 : WRITE (iout, '(A,I5,/,A,/,A)') ' KPSYM| DIMENSION NKPOINT =', nkpoint, &
527 0 : ' KPSYM| INSUFFICIENT FOR ACCOMMODATING ALL THE SPECIAL POINTS', &
528 0 : ' KPSYM| WHAT FOLLOWS IS AN INCOMPLETE LIST'
529 : END IF
530 0 : ntot = ABS(ntot)
531 : END IF
532 : ! Before using the list WVKL as wave vectors, they have to be
533 : ! multiplied by 2*Pi
534 : ! The list of weights LWGHT is not normalized
535 960 : iswght = 0
536 3470 : DO i = 1, ntot
537 3470 : iswght = iswght + lwght(i)
538 : END DO
539 960 : IF (iout > 0) THEN
540 : WRITE (iout, '(8X,A,T33,A,4X,A)') &
541 335 : 'WAVEVECTOR K', 'WEIGHT', 'UNFOLDING ROTATIONS'
542 : END IF
543 : ! Set near-zeroes equal to zero:
544 3470 : DO l = 1, ntot
545 10040 : DO i = 1, 3
546 10040 : IF (ABS(wvkl(i, l)) < delta) wvkl(i, l) = 0._dp
547 : END DO
548 2510 : IF (istrin /= 0) THEN
549 : ! Express special points in (unstrained) basis.
550 0 : proj1 = 0._dp
551 0 : proj2 = 0._dp
552 0 : proj3 = 0._dp
553 0 : DO i = 1, 3
554 0 : proj1 = proj1 + wvkl(i, l)*a01(i)
555 0 : proj2 = proj2 + wvkl(i, l)*a02(i)
556 0 : proj3 = proj3 + wvkl(i, l)*a03(i)
557 : END DO
558 0 : DO i = 1, 3
559 0 : wvkl(i, l) = proj1*b1(i) + proj2*b2(i) + proj3*b3(i)
560 : END DO
561 : END IF
562 2510 : lmax = lwght(l)
563 2510 : IF (iout > 0) THEN
564 : WRITE (iout, fmt='(1X,I5,3F8.4,I8,T42,12I3)') &
565 846 : l, (wvkl(i, l), i=1, 3), lwght(l), (lrot(i, l), i=1, MIN(lmax, 12))
566 : END IF
567 3650 : DO j = 13, lmax, 12
568 2690 : IF (iout > 0) THEN
569 : WRITE (iout, fmt='(T42,12I3)') &
570 61 : (lrot(i, l), i=j, MIN(lmax, j - 1 + 12))
571 : END IF
572 : END DO
573 : END DO
574 960 : IF (iout > 0) THEN
575 335 : WRITE (iout, '(24X,"TOTAL:",I8)') iswght
576 : END IF
577 : END SUBROUTINE k290s
578 : ! **************************************************************************************************
579 :
580 : ! **************************************************************************************************
581 : !> \brief ...
582 : !> \param iout ...
583 : !> \param a1 ...
584 : !> \param a2 ...
585 : !> \param a3 ...
586 : !> \param nat ...
587 : !> \param ty ...
588 : !> \param x ...
589 : !> \param b1 ...
590 : !> \param b2 ...
591 : !> \param b3 ...
592 : !> \param ihg ...
593 : !> \param ihc ...
594 : !> \param isy ...
595 : !> \param li ...
596 : !> \param nc ...
597 : !> \param indpg ...
598 : !> \param ib ...
599 : !> \param ntvec ...
600 : !> \param v ...
601 : !> \param f0 ...
602 : !> \param r ...
603 : !> \param tvec ...
604 : !> \param origin ...
605 : !> \param rx ...
606 : !> \param isc ...
607 : !> \param delta ...
608 : ! **************************************************************************************************
609 2880 : SUBROUTINE group1s(iout, a1, a2, a3, nat, ty, x, b1, b2, b3, &
610 : ihg, ihc, isy, li, nc, indpg, ib, ntvec, &
611 2880 : v, f0, r, tvec, origin, rx, isc, delta)
612 : ! ==--------------------------------------------------------------==
613 : ! == WRITTEN ON SEPTEMBER 10TH - FROM THE ACMI COMPLEX ==
614 : ! == (WORLTON AND WARREN, COMPUT.PHYS.COMMUN. 8,71-84 (1974)) ==
615 : ! == (AND 3,88-117 (1972)) ==
616 : ! == BASIC CRYSTALLOGRAPHIC INFORMATION ==
617 : ! == ABOUT A GIVEN CRYSTAL STRUCTURE. ==
618 : ! == SUBROUTINES NEEDED: PGL1,ATFTM1,ROT1,RLV3 ==
619 : ! ==--------------------------------------------------------------==
620 : ! == INPUT DATA: ==
621 : ! == IOUT ... NUMBER OF THE OUTPUT UNIT FOR ON-LINE PRINTING ==
622 : ! == OF VARIOUS MESSAGES ==
623 : ! == IF IOUT<=0 NO MESSAGE ==
624 : ! == A1,A2,A3 .. ELEMENTARY TRANSLATIONS OF THE LATTICE, IN SOME ==
625 : ! == UNIT OF LENGTH ==
626 : ! == NAT .... NUMBER OF ATOMS IN THE UNIT CELL ==
627 : ! == ALL THE DIMENSIONS ARE SET FOR NAT <= 20 ==
628 : ! == TY ..... INTEGERS DISTINGUISHING BETWEEN THE ATOMS OF ==
629 : ! == DIFFERENT TYPE. TY(I) IS THE TYPE OF THE I-TH ATOM ==
630 : ! == OF THE BASIS ==
631 : ! == X ...... CARTESIAN COORDINATES OF THE NAT ATOMS OF THE BASIS ==
632 : ! == DELTA... REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
633 : ! ==--------------------------------------------------------------==
634 : ! == OUTPUT DATA: ==
635 : ! == B1,B2,B3 .. RECIPROCAL LATTICE VECTORS, NOT MULTIPLIED BY ==
636 : ! == ANY 2PI, IN UNITS RECIPROCAL TO THOSE OF A1,A2,A3 ==
637 : ! == IHG .... POINT GROUP OF THE PRIMITIVE LATTICE, HOLOHEDRAL ==
638 : ! == GROUP NUMBER: ==
639 : ! == IHG=1 STANDS FOR TRICLINIC SYSTEM ==
640 : ! == IHG=2 STANDS FOR MONOCLINIC SYSTEM ==
641 : ! == IHG=3 STANDS FOR ORTHORHOMBIC SYSTEM ==
642 : ! == IHG=4 STANDS FOR TETRAGONAL SYSTEM ==
643 : ! == IHG=5 STANDS FOR CUBIC SYSTEM ==
644 : ! == IHG=6 STANDS FOR TRIGONAL SYSTEM ==
645 : ! == IHG=7 STANDS FOR HEXAGONAL SYSTEM ==
646 : ! == IHC .... CODE DISTINGUISHING BETWEEN HEXAGONAL AND CUBIC ==
647 : ! == GROUPS ==
648 : ! == IHC=0 STANDS FOR HEXAGONAL GROUPS ==
649 : ! == IHC=1 STANDS FOR CUBIC GROUPS ==
650 : ! == ISY .... CODE INDICATING WHETHER THE SPACE GROUP IS ==
651 : ! == SYMMORPHIC OR NONSYMMORPHIC ==
652 : ! == ISY= 0 NONSYMMORPHIC GROUP ==
653 : ! == ISY= 1 SYMMORPHIC GROUP ==
654 : ! == ISY=-1 SYMMORPHIC GROUP WITH NON-STANDARD ORIGIN ==
655 : ! == ISY=-2 UNDETERMINED (NORMALLY NEVER) ==
656 : ! == THE GROUP IS CONSIDERED SYMMORPHIC IF FOR EACH ==
657 : ! == OPERATION OF THE POINT GROUP THE SUM OF THE 3 ==
658 : ! == COMPONENTS OF ABS(V(N)) (NONPRIMITIVE TRANSLATION, ==
659 : ! == SEE BELOW) IS LT. 0.0001 ==
660 : ! == ORIGIN STANDARD ORIGIN IF SYMMORPHIC (CRYSTAL COORDINATES) ==
661 : ! == LI ..... CODE INDICATING WHETHER THE POINT GROUP ==
662 : ! == OF THE CRYSTAL CONTAINS INVERSION OR NOT ==
663 : ! == (OPERATIONS 13 OR 25 IN RESPECTIVELY HEXAGONAL ==
664 : ! == OR CUBIC GROUPS). ==
665 : ! == LI=0 MEANS: DOES NOT CONTAIN INVERSION ==
666 : ! == LI>0 MEANS: THERE IS INVERSION IN THE POINT ==
667 : ! == GROUP OF THE CRYSTAL ==
668 : ! == NC ..... TOTAL NUMBER OF ELEMENTS IN THE POINT GROUP OF THE ==
669 : ! == CRYSTAL ==
670 : ! == INDPG .. POINT GROUP INDEX (DETERMINED IF SYMMORPHIC GROUP) ==
671 : ! == IB ..... LIST OF THE ROTATIONS CONSTITUTING THE POINT GROUP ==
672 : ! == OF THE CRYSTAL. THE NUMBERING IS THAT DEFINED IN ==
673 : ! == WORLTON AND WARREN, I.E. THE ONE MATERIALIZED IN THE==
674 : ! == ARRAY R (SEE BELOW) ==
675 : ! == ONLY THE FIRST NC ELEMENTS OF THE ARRAY IB ARE ==
676 : ! == MEANINGFUL ==
677 : ! == NTVEC .. NUMBER OF TRANSLATIONAL VECTORS ==
678 : ! == ASSOCIATED WITH IDENTITY OPERATOR I.E. ==
679 : ! == GIVES THE NUMBER OF IDENTICAL PRIMITIVE CELLS ==
680 : ! == V ...... NONPRIMITIVE TRANSLATIONS (IN THE CASE OF NONSYMMOR-==
681 : ! == PHIC GROUPS). V(I,N) IS THE I-TH COMPONENT ==
682 : ! == OF THE TRANSLATION CONNECTED WITH THE N-TH ELEMENT ==
683 : ! == OF THE POINT GROUP (I.E. WITH THE ROTATION ==
684 : ! == NUMBER IB(N) ). ==
685 : ! == ATTENTION: V(I) ARE NOT CARTESIAN COMPONENTS, ==
686 : ! == THEY REFER TO THE SYSTEM A1,A2,A3. ==
687 : ! == F0 ..... THE FUNCTION DEFINED IN MARADUDIN, IPATOVA BY ==
688 : ! == EQ. (3.2.12): ATOM TRANSFORMATION TABLE. ==
689 : ! == THE ELEMENT F0(N,KAPA) MEANS THAT THE N-TH ==
690 : ! == OPERATION OF THE SPACE GROUP (I.E. OPERATION NUMBER ==
691 : ! == IB(N), TOGETHER WITH AN EVENTUAL NONPRIMITIVE ==
692 : ! == TRANSLATION V(N)) TRANSFERS THE ATOM KAPA INTO THE ==
693 : ! == ATOM F0(N,KAPA). ==
694 : ! == THE 49TH LINE GIVES EQUIVALENT ATOMS FOR ==
695 : ! == FRACTIONAl TRANSLATIONS ASSOCIATED WITH IDENTITY ==
696 : ! == R ...... LIST OF THE 3 X 3 ROTATION MATRICES ==
697 : ! == (XYZ REPRESENTATION OF THE O(H) OR D(6)H GROUPS) ==
698 : ! == ALL 48 OR 24 MATRICES ARE LISTED. ==
699 : ! == FOLLOW NOTATION OF WORLTON-WARREN(1972) ==
700 : ! == TVEC .. LIST OF NTVEC TRANSLATIONAL VECTORS ==
701 : ! == ASSOCIATED WITH IDENTITY OPERATOR ==
702 : ! == TVEC(1:3,1) = \(0,0,0\) ==
703 : ! == (CRYSTAL COORDINATES) ==
704 : ! == RX ..... SCRATCH ARRAY ==
705 : ! == ISC .... SCRATCH ARRAY ==
706 : ! ==--------------------------------------------------------------==
707 : ! == PRINTED OUTPUT: ==
708 : ! == PROGRAM PRINTS THE TYPE OF THE LATTICE (IHG, IN WORDS), ==
709 : ! == LISTS THE OPERATIONS OF THE POINT GROUP OF THE ==
710 : ! == CRYSTAL, INDICATES WHETHER THE SPACE GROUP IS SYMMORPHIC OR ==
711 : ! == NONSYMMORPHIC AND WHETHER THE POINT GROUP OF THE CRYSTAL ==
712 : ! == CONTAINS INVERSION. ==
713 : ! ==--------------------------------------------------------------==
714 : INTEGER :: iout
715 : REAL(dp) :: a1(3), a2(3), a3(3)
716 : INTEGER :: nat, ty(nat)
717 : REAL(dp) :: x(3, nat), b1(3), b2(3), b3(3)
718 : INTEGER :: ihg, ihc, isy, li, nc, indpg, ib(48), &
719 : ntvec
720 : REAL(dp) :: v(3, 48)
721 : INTEGER :: f0(49, nat)
722 : REAL(dp) :: r(3, 3, 48), tvec(3, nat), origin(3), &
723 : rx(3, nat)
724 : INTEGER :: isc(nat)
725 : REAL(dp) :: delta
726 :
727 : INTEGER :: i, ncprim
728 : REAL(dp) :: a(3, 3), ai(3, 3), ap(3, 3), api(3, 3)
729 :
730 11520 : DO i = 1, 3
731 8640 : a(i, 1) = a1(i)
732 8640 : a(i, 2) = a2(i)
733 11520 : a(i, 3) = a3(i)
734 : END DO
735 : ! ==--------------------------------------------------------------==
736 : ! == A(I,J) IS THE I-TH CARTESIAN COMPONENT OF THE J-TH PRIMITIVE ==
737 : ! == TRANSLATION VECTOR OF THE DIRECT LATTICE ==
738 : ! == TY(I) IS AN INTEGER DISTINGUISHING ATOMS OF DIFFERENT TYPE, ==
739 : ! == I.E., DIFFERENT ATOMIC SPECIES ==
740 : ! == X(J,I) IS THE J-TH CARTESIAN COMPONENT OF THE POSITION ==
741 : ! == VECTOR FOR THE I-TH ATOM IN THE UNIT CELL. ==
742 : ! ==--------------------------------------------------------------==
743 : ! ==DETERMINE PRIMITIVE LATTICE VECTORS FOR THE RECIPROCAL LATTICE==
744 : ! ==--------------------------------------------------------------==
745 2880 : CALL calbrec(a, ai)
746 11520 : DO i = 1, 3
747 8640 : b1(i) = ai(1, i)
748 8640 : b2(i) = ai(2, i)
749 11520 : b3(i) = ai(3, i)
750 : END DO
751 : ! ==--------------------------------------------------------------==
752 : ! Determination of the translation vectors associated with
753 : ! the Identity matrix i.e. if the cell is duplicated
754 : ! Give also the ``primitive lattice''
755 2880 : CALL primlatt(a, ai, ap, api, nat, ty, x, ntvec, tvec, f0, isc, delta)
756 : ! ==--------------------------------------------------------------==
757 : ! Determination of the holohedral group (and crystal system)
758 2880 : CALL pgl1(ap, api, ihc, nc, ib, ihg, r, delta)
759 2880 : IF (ntvec > 1) THEN
760 : ! All rotations found by PGL1 have axes in x, y or z cart. axis
761 : ! So we have too check if we do not loose symmetry
762 372 : ncprim = nc
763 : ! The hexagonal system is found if the z axis is the sixfold axis
764 372 : CALL pgl1(a, ai, ihc, nc, ib, ihg, r, delta)
765 372 : IF (ncprim > nc) THEN
766 : ! More symmetry with
767 0 : CALL pgl1(ap, api, ihc, nc, ib, ihg, r, delta)
768 : END IF
769 : END IF
770 :
771 : ! Determination of the space group
772 : CALL atftm1(iout, r, v, x, f0, origin, ib, ty, nat, ihg, ihc, rx, &
773 2880 : nc, indpg, ntvec, a, ai, li, isy, isc, delta)
774 :
775 2880 : IF (iout > 0) THEN
776 670 : IF (li > 0) THEN
777 : IF (iout > 0) THEN
778 : WRITE (iout, '(1X,A)') &
779 481 : 'KPSYM| THE POINT GROUP OF THE CRYSTAL CONTAINS THE INVERSION'
780 : END IF
781 : END IF
782 670 : IF (iout > 0) THEN
783 670 : WRITE (iout, *)
784 : END IF
785 : END IF
786 :
787 2880 : END SUBROUTINE group1s
788 : ! **************************************************************************************************
789 : !> \brief ...
790 : !> \param a ...
791 : !> \param ai ...
792 : ! **************************************************************************************************
793 3608 : SUBROUTINE calbrec(a, ai)
794 : ! ==--------------------------------------------------------------==
795 : ! == CALCULATE RECIPROCAL VECTOR BASIS (AI(1:3,1:3)) ==
796 : ! == INPUT: ==
797 : ! == A(3,3) A(I,J) IS THE I-TH CARTESIAN COMPONENT ==
798 : ! == OF THE J-TH PRIMITIVE TRANSLATION VECTOR OF ==
799 : ! == THE DIRECT LATTICE ==
800 : ! == OUTPUT: ==
801 : ! == AI(3,3) RECIPROCAL VECTOR BASIS ==
802 : ! ==--------------------------------------------------------------==
803 : REAL(dp) :: a(3, 3), ai(3, 3)
804 :
805 : INTEGER :: i, il, iu, j, jl, ju
806 : REAL(dp) :: det
807 :
808 : det = a(1, 1)*a(2, 2)*a(3, 3) + a(2, 1)*a(1, 3)*a(3, 2) + &
809 : a(3, 1)*a(1, 2)*a(2, 3) - a(1, 1)*a(2, 3)*a(3, 2) - &
810 3608 : A(2, 1)*A(1, 2)*A(3, 3) - A(3, 1)*A(1, 3)*A(2, 2)
811 3608 : det = 1._dp/det
812 14432 : DO i = 1, 3
813 10824 : il = 1
814 10824 : iu = 3
815 10824 : IF (i == 1) il = 2
816 7216 : IF (i == 3) iu = 2
817 46904 : DO j = 1, 3
818 32472 : jl = 1
819 32472 : ju = 3
820 32472 : IF (j == 1) jl = 2
821 21648 : IF (j == 3) ju = 2
822 : ai(j, i) = (-1._dp)**(i + j)*det* &
823 43296 : (A(IL, JL)*A(IU, JU) - A(IL, JU)*A(IU, JL))
824 : END DO
825 : END DO
826 : ! ==--------------------------------------------------------------==
827 3608 : RETURN
828 : END SUBROUTINE calbrec
829 : ! ==================================================================
830 : ! **************************************************************************************************
831 : !> \brief ...
832 : !> \param a ...
833 : !> \param ai ...
834 : !> \param ap ...
835 : !> \param api ...
836 : !> \param nat ...
837 : !> \param ty ...
838 : !> \param x ...
839 : !> \param ntvec ...
840 : !> \param tvec ...
841 : !> \param f0 ...
842 : !> \param isc ...
843 : !> \param delta ...
844 : ! **************************************************************************************************
845 2880 : SUBROUTINE primlatt(a, ai, ap, api, nat, ty, x, ntvec, tvec, f0, isc, delta)
846 : ! ==--------------------------------------------------------------==
847 : ! == DETERMINATION OF THE TRANSLATION VECTORS ASSOCIATED WITH ==
848 : ! == THE IDENTITY SYMMETRY I.E. IF THE CELL IS DUPLICATED ==
849 : ! == GIVE ALSO THE PRIMITIVE DIRECT AND RECIPROCAL LATTICE VECTOR ==
850 : ! ==--------------------------------------------------------------==
851 : ! == INPUT: ==
852 : ! == A(3,3) A(I,J) IS THE I-TH CARTESIAN COMPONENT ==
853 : ! == OF THE J-TH TRANSLATION VECTOR OF ==
854 : ! == THE DIRECT LATTICE ==
855 : ! == AI(3,3) RECIPROCAL VECTOR BASIS (CARTESIAN) ==
856 : ! == NAT NUMBER OF ATOMS ==
857 : ! == TY(NAT) TYPE OF ATOMS ==
858 : ! == X(3,NAT) ATOMIC COORDINATES IN CARTESIAN COORDINATES ==
859 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
860 : ! == OUTPUT: ==
861 : ! == AP(3,3) COMPONENTS OF THE PRIMITIVE TRANSLATION VECTORS ==
862 : ! == API(3,3) PRIMITIVE RECIPROCAL BASIS VECTORS ==
863 : ! == BOTH BAISI ARE IN CARTESIAN COORDINATES ==
864 : ! == NTVEC NUMBER OF TRANSLATION VECTORS (FRACTIONNAL) ==
865 : ! == TVEC(3,NTVEC) COMPONENTS OF TRANSLATIONAL VECTORS ==
866 : ! == (CRYSTAL COORDINATES) ==
867 : ! == F0(49,NAT) GIVES INEQUIVALENT ATOM FOR EACH ATOM ==
868 : ! == THE 49-TH LINE ==
869 : ! == ISC(NAT) SCRATCH ARRAY ==
870 : ! ==--------------------------------------------------------------==
871 : REAL(dp) :: a(3, 3), ai(3, 3), ap(3, 3), api(3, 3)
872 : INTEGER :: nat, ty(nat)
873 : REAL(dp) :: x(3, nat)
874 : INTEGER :: ntvec
875 : REAL(dp) :: tvec(3, nat)
876 : INTEGER :: f0(49, nat), isc(nat)
877 : REAL(dp) :: delta
878 :
879 : INTEGER :: i, il, iv, j, k2
880 : LOGICAL :: oksym
881 : REAL(dp) :: vr(3), xb(3)
882 :
883 : ! Variables
884 : ! ==--------------------------------------------------------------==
885 : ! First we check if there exist fractional translational vectors
886 : ! associated with Identity operation i.e.
887 : ! if the cell is duplicated or not.
888 :
889 2880 : ntvec = 1
890 2880 : tvec(1, 1) = 0._dp
891 2880 : tvec(2, 1) = 0._dp
892 2880 : tvec(3, 1) = 0._dp
893 14044 : DO i = 1, nat
894 14044 : f0(49, i) = i
895 : END DO
896 11164 : DO k2 = 2, nat
897 8284 : IF (ty(1) /= ty(k2)) CYCLE
898 28752 : DO i = 1, 3
899 28752 : xb(i) = x(i, k2) - x(i, 1)
900 : END DO
901 : ! A fractional translation vector VR is defined.
902 7188 : CALL rlv3(ai, xb, vr, il, delta)
903 7188 : CALL checkrlv3(1, nat, ty, x, x, vr, f0, ai, isc, .TRUE., oksym, delta)
904 10068 : IF (oksym) THEN
905 : ! A fractional translational vector is found
906 1084 : ntvec = ntvec + 1
907 : ! F0(49,1:NAT) gives number of equivalent atoms
908 : ! and has atom indexes of inequivalent atoms (for translation)
909 9660 : DO i = 1, nat
910 9660 : IF (f0(49, i) > f0(1, i)) f0(49, i) = f0(1, i)
911 : END DO
912 4336 : DO i = 1, 3
913 4336 : tvec(i, ntvec) = vr(i)
914 : END DO
915 : END IF
916 : END DO
917 : ! ==-------------------------------------------------------------==
918 11520 : DO i = 1, 3
919 8640 : ap(1, i) = a(1, i)
920 8640 : ap(2, i) = a(2, i)
921 8640 : ap(3, i) = a(3, i)
922 8640 : api(1, i) = ai(1, i)
923 8640 : api(2, i) = ai(2, i)
924 11520 : api(3, i) = ai(3, i)
925 : END DO
926 2880 : IF (ntvec == 1) THEN
927 : ! The current cell is definitely a primitive one
928 : ! Copy A and AI to AP and API
929 : ELSE
930 : ! We are looking for the primitive lattice vector basis set
931 : ! AP is our current lattice vector basis
932 1456 : DO iv = 2, ntvec
933 : ! TVEC in cartesian coordinates
934 4336 : DO i = 1, 3
935 : xb(i) = tvec(1, iv)*a(i, 1) &
936 : + TVEC(2, IV)*A(I, 2) &
937 4336 : + TVEC(3, IV)*A(I, 3)
938 : END DO
939 : ! We calculare TVEC in AP basis
940 1084 : CALL rlv3(api, xb, vr, il, delta)
941 2880 : DO i = 1, 3
942 2508 : IF (ABS(vr(i)) > delta) THEN
943 728 : il = NINT(1._dp/ABS(vr(i)))
944 728 : IF (il > 1) THEN
945 : ! We replace AP(1:3,I) by TVEC(1:3,IV)
946 2912 : DO j = 1, 3
947 2912 : ap(j, i) = xb(j)
948 : END DO
949 : ! Calculate new API
950 728 : CALL calbrec(ap, api)
951 728 : EXIT
952 : END IF
953 : END IF
954 : END DO
955 : END DO
956 : END IF
957 : ! ==--------------------------------------------------------------==
958 2880 : RETURN
959 : END SUBROUTINE primlatt
960 : ! ==================================================================
961 : ! **************************************************************************************************
962 : !> \brief ...
963 : !> \param a ...
964 : !> \param ai ...
965 : !> \param ihc ...
966 : !> \param nc ...
967 : !> \param ib ...
968 : !> \param ihg ...
969 : !> \param r ...
970 : !> \param delta ...
971 : ! **************************************************************************************************
972 3252 : SUBROUTINE pgl1(a, ai, ihc, nc, ib, ihg, r, delta)
973 : ! ==--------------------------------------------------------------==
974 : ! == WRITTEN ON SEPTEMBER 11TH, 1979 - FROM ACMI COMPLEX ==
975 : ! == AUXILIARY SUBROUTINE TO GROUP1 ==
976 : ! == SUBROUTINE PGL DETERMINES THE POINT GROUP OF THE LATTICE ==
977 : ! == AND THE CRYSTAL SYSTEM. ==
978 : ! == SUBROUTINES NEEDED: ROT1, RLV3 ==
979 : ! ==--------------------------------------------------------------==
980 : ! == WARNING: FOR THE HEXAGONAL SYSTEM, THE 3RD AXIS SUPPOSE ==
981 : ! == TO BE THE SIX-FOLD AXIS ==
982 : ! ==--------------------------------------------------------------==
983 : ! == INPUT: ==
984 : ! == A ..... DIRECT LATTICE VECTORS ==
985 : ! == AI .... RECIPROCAL LATTICE VECTORS ==
986 : ! == DELTA.. REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
987 : ! ==--------------------------------------------------------------==
988 : ! == OUTPUT: ==
989 : ! == IHC .... CODE DISTINGUISHING BETWEEN HEXAGONAL AND CUBIC ==
990 : ! == GROUPS ==
991 : ! == IHC=0 STANDS FOR HEXAGONAL GROUPS ==
992 : ! == IHC=1 STANDS FOR CUBIC GROUPS ==
993 : ! == NC .... NUMBER OF ROTATIONS IN THE POINT GROUP ==
994 : ! == IB .... SET OF ROTATION ==
995 : ! == IHG .... POINT GROUP OF THE PRIMITIVE LATTICE, HOLOHEDRAL ==
996 : ! == GROUP NUMBER: ==
997 : ! == IHG=1 STANDS FOR TRICLINIC SYSTEM ==
998 : ! == IHG=2 STANDS FOR MONOCLINIC SYSTEM ==
999 : ! == IHG=3 STANDS FOR ORTHORHOMBIC SYSTEM ==
1000 : ! == IHG=4 STANDS FOR TETRAGONAL SYSTEM ==
1001 : ! == IHG=5 STANDS FOR CUBIC SYSTEM ==
1002 : ! == IHG=6 STANDS FOR TRIGONAL SYSTEM ==
1003 : ! == IHG=7 STANDS FOR HEXAGONAL SYSTEM ==
1004 : ! == R ...... LIST OF THE 3 X 3 ROTATION MATRICES ==
1005 : ! == (XYZ REPRESENTATION OF THE O(H) OR D(6)H GROUPS) ==
1006 : ! == ALL 48 OR 24 MATRICES ARE LISTED. ==
1007 : ! == FOLLOW NOTATION OF WORLTON-WARREN(1972) ==
1008 : ! ==--------------------------------------------------------------==
1009 : REAL(dp) :: a(3, 3), ai(3, 3)
1010 : INTEGER :: ihc, nc, ib(48), ihg
1011 : REAL(dp) :: r(3, 3, 48), delta
1012 :
1013 : INTEGER :: i, j, k, lx, n, nr
1014 : REAL(dp) :: tr, vr(3), xa(3)
1015 :
1016 6294 : DO ihc = 0, 1
1017 : ! IHC is 0 for hexagonal groups and 1 for cubic groups.
1018 6294 : IF (ihc == 0) THEN
1019 : nr = 24
1020 : ELSE
1021 3042 : nr = 48
1022 : END IF
1023 6294 : nc = 0
1024 : ! Constructs rotation operations.
1025 6294 : CALL rot1(ihc, r)
1026 230358 : loop_rotation: DO n = 1, nr
1027 224064 : ib(n) = 0
1028 : ! Rotate the A1,2,3 vectors by rotation No. N
1029 610424 : DO k = 1, 3
1030 1953408 : DO i = 1, 3
1031 1465056 : xa(i) = 0._dp
1032 6348576 : DO j = 1, 3
1033 5860224 : xa(i) = xa(i) + r(i, j, n)*a(j, k)
1034 : END DO
1035 : END DO
1036 488352 : CALL rlv3(ai, xa, vr, lx, delta)
1037 488352 : tr = 0._dp
1038 1953408 : DO i = 1, 3
1039 1953408 : tr = tr + ABS(vr(i))
1040 : END DO
1041 : ! If VR.ne.0, then XA cannot be a multiple of a lattice vector
1042 610424 : IF (tr > delta) CYCLE loop_rotation
1043 : END DO
1044 122072 : nc = nc + 1
1045 230358 : ib(nc) = n
1046 : END DO loop_rotation
1047 : ! ==------------------------------------------------------------==
1048 : ! IHG stands for holohedral group number.
1049 6294 : IF (ihc == 0) THEN
1050 : ! Hexagonal group:
1051 3252 : IF (nc == 12) ihg = 6
1052 3252 : IF (nc > 12) ihg = 7
1053 3252 : IF (nc >= 12) RETURN
1054 : ! Too few operations, try cubic group: (IHC=1,NR=48)
1055 : ELSE
1056 : ! Cubic group:
1057 3042 : IF (nc < 4) ihg = 1
1058 3042 : IF (nc == 4) ihg = 2
1059 3042 : IF (nc > 4) ihg = 3
1060 3042 : IF (nc == 16) ihg = 4
1061 3042 : IF (nc > 16) ihg = 5
1062 3042 : RETURN
1063 : END IF
1064 : END DO
1065 : ! ==--------------------------------------------------------------==
1066 : RETURN
1067 : END SUBROUTINE pgl1
1068 : ! ==================================================================
1069 : ! **************************************************************************************************
1070 : !> \brief ...
1071 : !> \param ai ...
1072 : !> \param xb ...
1073 : !> \param vr ...
1074 : !> \param il ...
1075 : !> \param delta ...
1076 : ! **************************************************************************************************
1077 3111728 : SUBROUTINE rlv3(ai, xb, vr, il, delta)
1078 : ! ==--------------------------------------------------------------==
1079 : ! == WRITTEN ON SEPTEMBER 11TH, 1979 - FROM ACMI COMPLEX ==
1080 : ! == AUXILIARY SUBROUTINE TO GROUP1 ==
1081 : ! == SUBROUTINE RLV REMOVES A DIRECT LATTICE VECTOR ==
1082 : ! == FROM XB LEAVING THE REMAINDER IN VR. ==
1083 : ! == IF A NONZERO LATTICE VECTOR WAS REMOVED, IL IS MADE NONZERO. ==
1084 : ! == VR STANDS FOR V-REFERENCE. ==
1085 : ! ==--------------------------------------------------------------==
1086 : ! == INPUT: ==
1087 : ! == AI(I,J) ARE THE RECIPROCAL LATTICE VECTORS, ==
1088 : ! == B(I) = AI(I,J),J=1,2,3 ==
1089 : ! == XB(1:3) VECTOR IN CARTESIAN COORDINATES ==
1090 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
1091 : ! == OUTPUT: ==
1092 : ! == VR IS NOT GIVEN IN CARTESIAN COORDINATES BUT ==
1093 : ! == IN THE SYSTEM A1,A2,A3 (CRYSTAL COORDINATES) ==
1094 : ! == AND BETWEEN -1/2 AND 1/2 ==
1095 : ! == IL ABS OF VR ==
1096 : ! == K.K., 23.10.1979 ==
1097 : ! ==--------------------------------------------------------------==
1098 : REAL(dp) :: ai(3, 3), xb(3), vr(3)
1099 : INTEGER :: il
1100 : REAL(dp) :: delta
1101 :
1102 : INTEGER :: i
1103 : REAL(dp) :: ts
1104 :
1105 3111728 : il = 0
1106 12446912 : DO i = 1, 3
1107 12446912 : vr(i) = 0._dp
1108 : END DO
1109 3111728 : ts = ABS(xb(1)) + ABS(xb(2)) + ABS(xb(3))
1110 3111728 : IF (ts <= delta) RETURN
1111 11578336 : DO i = 1, 3
1112 8683752 : vr(i) = vr(i) + ai(i, 1)*xb(1) + ai(i, 2)*xb(2) + ai(i, 3)*xb(3)
1113 8683752 : il = il + NINT(ABS(vr(i)))
1114 : ! Change in order to have correct determination of origin and
1115 : ! symmorphic group (T.D 30/03/98)
1116 : ! VR(I) = - MOD(real(VR(I),kind=dp),1._dp)
1117 11578336 : vr(i) = NINT(vr(i)) - vr(i)
1118 : END DO
1119 : ! ==--------------------------------------------------------------==
1120 : RETURN
1121 : END SUBROUTINE rlv3
1122 : ! ==================================================================
1123 : ! **************************************************************************************************
1124 : !> \brief ...
1125 : !> \param iout ...
1126 : !> \param r ...
1127 : !> \param v ...
1128 : !> \param x ...
1129 : !> \param f0 ...
1130 : !> \param origin ...
1131 : !> \param ib ...
1132 : !> \param ty ...
1133 : !> \param nat ...
1134 : !> \param ihg ...
1135 : !> \param ihc ...
1136 : !> \param rx ...
1137 : !> \param nc ...
1138 : !> \param indpg ...
1139 : !> \param ntvec ...
1140 : !> \param a ...
1141 : !> \param ai ...
1142 : !> \param li ...
1143 : !> \param isy ...
1144 : !> \param isc ...
1145 : !> \param delta ...
1146 : ! **************************************************************************************************
1147 2880 : SUBROUTINE atftm1(iout, r, v, x, f0, origin, ib, ty, nat, ihg, ihc, &
1148 2880 : rx, nc, indpg, ntvec, a, ai, li, isy, isc, delta)
1149 : ! ==--------------------------------------------------------------==
1150 : ! == WRITTEN ON SEPTEMBER 11TH, 1979 - FROM ACMI COMPLEX ==
1151 : ! == AUXILIARY SUBROUTINE TO GROUP1 ==
1152 : ! == SUBROUTINE ATFTMT DETERMINES ==
1153 : ! == THE POINT GROUP OF THE CRYSTAL, ==
1154 : ! == THE ATOM TRANSFORMATION TABLE,F0, ==
1155 : ! == THE FRACTIONAL TRANSLATIONS,V, ==
1156 : ! == ASSOCIATED WITH EACH ROTATION. ==
1157 : ! == SUBROUTINES NEEDED: RLV3 CHECKRLV3 SYMMORPHIC XSTRING ==
1158 : ! == MAY 14TH,1998: A LOT OF CHANGES (ARGUMENTS) ==
1159 : ! == BETTER DETERMINATION OF V ==
1160 : ! == SEP 15TH,1998: DETERMINATION OF FRACTIONAL TRANSLATIONAL VEC.==
1161 : ! ==--------------------------------------------------------------==
1162 : ! == INPUT: ==
1163 : ! == IOUT Logical file number (output) ==
1164 : ! == If IOUT<=0 no message ==
1165 : ! == IHG Holohedral group number (determined by PGL1) ==
1166 : ! == IHC Code distinguishing between hexagonal and cubic groups==
1167 : ! == IHC=0 stands for hexagonal groups ==
1168 : ! == IHC=1 stands for cubic groups ==
1169 : ! == NC Number of rotation operations ==
1170 : ! == NAT Number of atoms (used in the routine) ==
1171 : ! == X Coordinates of atoms (cartesian) ==
1172 : ! == TY Type of atoms ==
1173 : ! == R Sets of transformation operations (cartesian) ==
1174 : ! == IB Index giving NC operations in R ==
1175 : ! == AI Reciprocal lattice vectors ==
1176 : ! == NTVEC Number of translational vectors ==
1177 : ! == associated with Identity ==
1178 : ! == if primitive cell NTVEC=1, TVEC=(0,0,0) ==
1179 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
1180 : ! == OUTPUT: ==
1181 : ! == RX(3,NAT) Scratch array ==
1182 : ! == ISC(NAT) Scratch array ==
1183 : ! == NC is modified (number of symmetry operations) ==
1184 : ! == INDPG Point group index ==
1185 : ! == V(3,48) The fractional translations associated ==
1186 : ! == with each rotation (crystal coordinates) ==
1187 : ! == F0(1:48,NAT) ==
1188 : ! == The atom transformation table for rotation (48,NAT) ==
1189 : ! == ORIGIN Standard origin if symmorphic (crystal coordinates) ==
1190 : ! == ISY = 1 Isommorphic group ==
1191 : ! == =-1 Isommorphic group with non-standard origin ==
1192 : ! == = 0 Non-Isommorphic group ==
1193 : ! == =-2 Undetermined (normally never) ==
1194 : ! == LI ..... Code indicating whether the point group ==
1195 : ! == of the crystal contains inversion or not ==
1196 : ! == (operations 13 or 25 in respectively hexagonal ==
1197 : ! == or cubic groups). ==
1198 : ! == LI=0 : does not contain inversion ==
1199 : ! == LI>0 : there is inversion in the point ==
1200 : ! == group of the crystal ==
1201 : ! ==--------------------------------------------------------------==
1202 : ! INDPG group indpg group indpg group indpg group ==
1203 : ! == 1 1 (c1) 9 3m (c3v) 17 4/mmm(d4h) 25 222(d2) ==
1204 : ! == 2 <1>(ci) 10 <3>m(d3d) 18 6 (c6) 26 mm2(c2v) ==
1205 : ! == 3 2 (c2) 11 4 (c4) 19 <6>(c3h) 27 mmm(d2h) ==
1206 : ! == 4 m (c1h) 12 <4>(s4) 20 6/m(c6h) 28 23 (t) ==
1207 : ! == 5 2/m(c2h) 13 4/m(c4h) 21 622(d6) 29 m3 (th) ==
1208 : ! == 6 3 (c3) 14 422(d4) 22 6mm(c6v) 30 432(o) ==
1209 : ! == 7 <3>(c3i) 15 4mm(c4v) 23 <6>m2(d3h) 31 <4>3m(td) ==
1210 : ! == 8 32 (d3) 16 <4>2m(d2d) 24 6/mmm(d6h) 32 m3m(oh) ==
1211 : ! ==--------------------------------------------------------------==
1212 : ! rname_cubic: Name of 48 rotations (convention Warren-Worlton)
1213 : INTEGER :: iout
1214 : REAL(dp) :: r(3, 3, 48), v(3, 48), origin(3)
1215 : INTEGER :: ib(48), nat, ty(nat), f0(49, nat)
1216 : REAL(dp) :: x(3, nat)
1217 : INTEGER :: ihg, ihc
1218 : REAL(dp) :: rx(3, nat)
1219 : INTEGER :: nc, indpg, ntvec
1220 : REAL(dp) :: a(3, 3), ai(3, 3)
1221 : INTEGER :: li, isy, isc(nat)
1222 : REAL(dp) :: delta
1223 :
1224 : CHARACTER(len=10), DIMENSION(48), PARAMETER :: rname_cubic = [' 1 ', ' 2[ 10 0] ', &
1225 : ' 2[ 01 0] ', ' 2[ 00 1] ', ' 3[-1-1-1]', ' 3[ 11-1] ', ' 3[-11 1] ', ' 3[ 1-11] ', &
1226 : ' 3[ 11 1] ', ' 3[-11-1] ', ' 3[-1-11] ', ' 3[ 1-1-1]', ' 2[-11 0] ', ' 4[ 00 1] ', &
1227 : ' 4[ 00-1] ', ' 2[ 11 0] ', ' 2[ 0-11] ', ' 2[ 01 1] ', ' 4[ 10 0] ', ' 4[-10 0] ', &
1228 : ' 2[-10 1] ', ' 4[ 0-10] ', ' 2[ 10 1] ', ' 4[ 01 0] ', '-1 ', '-2[ 10 0] ', &
1229 : '-2[ 01 0] ', '-2[ 00 1] ', '-3[-1-1-1]', '-3[ 11-1] ', '-3[-11 1] ', '-3[ 1-11] ', &
1230 : '-3[ 11 1] ', '-3[-11-1] ', '-3[-1-11] ', '-3[ 1-1-1]', '-2[-11 0] ', '-4[ 00 1] ', &
1231 : '-4[ 00-1] ', '-2[ 11 0] ', '-2[ 0-11] ', '-2[ 01 1] ', '-4[ 10 0] ', '-4[-10 0] ', &
1232 : '-2[-10 1] ', '-4[ 0-10] ', '-2[ 10 1] ', '-4[ 01 0] ']
1233 : CHARACTER(len=11), DIMENSION(24), PARAMETER :: rname_hexai = [' 1 ', ' 6[ 00 1] ', &
1234 : ' 3[ 00 1] ', ' 2[ 00 1] ', ' 3[ 00 -1] ', ' 6[ 00 -1] ', ' 2[ 01 0] ', ' 2[-11 0] ', &
1235 : ' 2[ 10 0] ', ' 2[ 21 0] ', ' 2[ 11 0] ', ' 2[ 12 0] ', '-1 ', '-6[ 00 1] ', &
1236 : '-3[ 00 1] ', '-2[ 00 1] ', '-3[ 00 -1] ', '-6[ 00 -1] ', '-2[ 01 0] ', '-2[-11 0] ', &
1237 : '-2[ 10 0] ', '-2[ 21 0] ', '-2[ 11 0] ', '-2[ 12 0] ']
1238 : CHARACTER(len=12), DIMENSION(7), PARAMETER :: icst = ['TRICLINIC ', 'MONOCLINIC ', &
1239 : 'ORTHORHOMBIC', 'TETRAGONAL ', 'CUBIC ', 'TRIGONAL ', 'HEXAGONAL ']
1240 : CHARACTER(len=3), DIMENSION(32), PARAMETER :: pgrd = ['c1 ', 'ci ', 'c2 ', 'c1h', 'c2h', &
1241 : 'c3 ', 'c3i', 'd3 ', 'c3v', 'd3 ', 'c4 ', 's4 ', 'c4h', 'd4 ', 'c4v', 'd2d', 'd4h', 'c6 ',&
1242 : 'c3h', 'c6h', 'd6 ', 'c6v', 'd3h', 'd6h', 'd2 ', 'c2v', 'd2h', 't ', 'th ', 'o ', 'td ',&
1243 : 'oh ']
1244 : CHARACTER(len=5), DIMENSION(32), PARAMETER :: pgrp = [' 1', ' <1>', ' 2', ' m', &
1245 : ' 2/m', ' 3', ' <3>', ' 32', ' 3m', ' <3>m', ' 4', ' <4>', ' 4/m', ' 422', &
1246 : ' 4mm', '<4>2m', '4/mmm', ' 6', ' <6>', ' 6/m', ' 622', ' 6mm', '<6>m2', '6/mmm', &
1247 : ' 222', ' mm2', ' mmm', ' 23', ' m3', ' 432', '<4>3m', ' m3m']
1248 :
1249 : INTEGER :: i, iis(48), il, info, j, k, k2, l, n, &
1250 : nca, ni
1251 : LOGICAL :: nodupli, oksym
1252 : REAL(dp) :: vc(3, 48), vr(3), vs, xb(3)
1253 :
1254 2880 : nodupli = ntvec == 1
1255 2880 : nca = 0
1256 141120 : DO n = 1, 48
1257 141120 : iis(n) = 0
1258 : END DO
1259 : ! Calculate translational vector for each operation
1260 : ! and atom transformation table.
1261 86736 : DO n = 1, nc
1262 83856 : l = ib(n)
1263 83856 : iis(l) = 1
1264 395376 : DO k = 1, nat
1265 1329936 : DO i = 1, 3
1266 1246080 : rx(i, k) = r(i, 1, l)*x(1, k) + r(i, 2, l)*x(2, k) + r(i, 3, l)*x(3, k)
1267 : END DO
1268 : END DO
1269 335424 : DO k = 1, 3
1270 335424 : vr(k) = 0._dp
1271 : END DO
1272 : ! First we determine for VR=(/0,0,0/)
1273 : ! IMPORTANT IF NOT UNIQUE ATOMS FOR DETERMINATION OF SYMMORPHIC
1274 83856 : CALL checkrlv3(n, nat, ty, rx, x, vr, f0, ai, isc, nodupli, oksym, delta)
1275 83856 : IF (.NOT. oksym) THEN
1276 : ! Now we try other possible VR
1277 : ! F0(49,1:NAT) has only inequivalent atom indexes for translation
1278 196524 : DO k2 = 1, nat
1279 172432 : IF (f0(49, k2) < k2) CYCLE
1280 151696 : IF (ty(1) /= ty(k2)) CYCLE
1281 529184 : DO i = 1, 3
1282 529184 : xb(i) = rx(i, 1) - x(i, k2)
1283 : END DO
1284 : ! A translation vector VR is defined.
1285 132296 : CALL rlv3(ai, xb, vr, il, delta)
1286 : ! ==----------------------------------------------------------==
1287 : ! == SUBROUTINE RLV3 REMOVES A DIRECT LATTICE VECTOR FROM XB ==
1288 : ! == LEAVING THE REMAINDER IN VR. IF A NONZERO LATTICE ==
1289 : ! == VECTOR WAS REMOVED, IL IS MADE NONZERO. ==
1290 : ! == VR STANDS FOR V-REFERENCE. ==
1291 : ! == VR IS NOT GIVEN IN CARTESIAN COORDINATES BUT ==
1292 : ! == IN THE SYSTEM A1,A2,A3. K.K., 23.10.1979 ==
1293 : ! ==----------------------------------------------------------==
1294 132296 : CALL checkrlv3(n, nat, ty, rx, x, vr, f0, ai, isc, nodupli, oksym, delta)
1295 156388 : IF (oksym) EXIT
1296 : END DO
1297 31912 : IF (.NOT. oksym) THEN
1298 24092 : iis(l) = 0
1299 24092 : CYCLE
1300 : END IF
1301 : END IF
1302 59764 : nca = nca + 1
1303 241936 : DO i = 1, 3
1304 239056 : v(i, nca) = vr(i)
1305 : END DO
1306 : ! ==------------------------------------------------------------==
1307 : ! == V(I,N) IS THE I-TH COMPONENT OF THE FRACTIONAL ==
1308 : ! == TRANSLATION ASSOCIATED WITH THE ROTATION N. ==
1309 : ! == ATTENTION: V(I) ARE NOT CARTESIAN COMPONENTS, THEY ARE ==
1310 : ! == GIVEN IN THE SYSTEM A1,A2,A3. ==
1311 : ! == K.K., 23.10. 1979 ==
1312 : ! ==------------------------------------------------------------==
1313 : END DO
1314 : ! Remove unused operations
1315 2880 : i = 0
1316 2880 : ni = 13
1317 2880 : IF (ihg < 6) ni = 25
1318 2880 : li = 0
1319 86736 : DO n = 1, nc
1320 83856 : l = ib(n)
1321 83856 : IF (iis(l) == 0) CYCLE
1322 59764 : i = i + 1
1323 59764 : ib(i) = ib(n)
1324 59764 : IF (ib(i) == ni) li = i
1325 239628 : DO k = 1, nat
1326 260840 : f0(i, k) = f0(n, k)
1327 : END DO
1328 : END DO
1329 : ! ==--------------------------------------------------------------==
1330 2880 : nc = i
1331 2880 : vs = 0._dp
1332 62644 : DO n = 1, nc
1333 62644 : vs = vs + ABS(v(1, n)) + ABS(v(2, n)) + ABS(v(3, n))
1334 : END DO
1335 : ! THE ORIGINAL VALUE DELTA=0.0001 WAS MODIFIED
1336 : ! BY K.K. , SEPTEMBER 1979 TO 0.0005
1337 : ! AND RETURNED TO 0.0001 BY RJN OCT 1987
1338 2880 : IF (vs > delta) THEN
1339 616 : isy = 0
1340 : ELSE
1341 2264 : isy = 1
1342 : END IF
1343 : ! ==--------------------------------------------------------------==
1344 : ! Determination of the point group
1345 : ! (Thierry Deutsch - 1998 [Maybe not complete!!])
1346 2880 : IF (ihg < 6) THEN
1347 2670 : IF (nc == 0) THEN
1348 0 : IF (iout > 0) THEN
1349 0 : WRITE (iout, '(" ATFTM1! IHG=",A," NC=",I2)') icst(ihg), nC
1350 : END IF
1351 0 : CPABORT('ATFTM1: NUMBER OF ROTATION NULL')
1352 : ! Triclinic system
1353 2670 : ELSE IF (nc == 1) THEN
1354 : ! IB=1
1355 368 : indpg = 1 ! 1 (c1)
1356 2302 : ELSE IF (nc == 2 .AND. ib(2) == 25) THEN
1357 : ! IB=125
1358 32 : indpg = 2 ! <1>(ci)
1359 2270 : ELSE IF (nc == 2 .AND. ( &
1360 : ib(2) == 4 .OR. & ! 2[001]
1361 : ib(2) == 2 .OR. & ! 2[100]
1362 : ib(2) == 3)) THEN ! 2[010]
1363 : ! Monoclinic system
1364 : ! IB=14 (z-axis) OR
1365 : ! IB=12 (x-axis) OR
1366 : ! IB=13 (y-axis)
1367 16 : indpg = 3 ! 2 (c2)
1368 2254 : ELSE IF (nc == 2 .AND. ( &
1369 : ib(2) == 28 .OR. &
1370 : ib(2) == 26 .OR. &
1371 : ib(2) == 27)) THEN
1372 : ! IB=128 (z-axis) OR
1373 : ! IB=126 (x-axis) OR
1374 : ! IB=127 (y-axis)
1375 120 : indpg = 4 ! m (c1h)
1376 2134 : ELSE IF (nc == 4 .AND. ( &
1377 : ib(4) == 28 .OR. & ! 2[001]
1378 : ib(4) == 27 .OR. & ! 2[010]
1379 : ib(4) == 26 .OR. & ! 2[100]
1380 : ib(4) == 37 .OR. & ! -2[-110]
1381 : ib(4) == 40)) THEN ! 2[110]
1382 : ! IB=1 425 28 (z-axis) OR
1383 : ! IB=1 225 26 (x-axis) OR
1384 : ! IB=1 325 27 (y-axis) OR
1385 : ! IB=113 2537 (-xy-axis)OR
1386 : ! IB=116 2540 (xy-axis)
1387 330 : indpg = 5 ! 2/m(c2h)
1388 1804 : ELSE IF (nc == 4 .AND. ( &
1389 : ib(4) == 15 .OR. &
1390 : ib(4) == 20 .OR. &
1391 : ib(4) == 24)) THEN
1392 : ! Tetragonal system
1393 : ! IB=14 1415 (z-axis) OR
1394 : ! IB=12 1920 (x-axis) OR
1395 : ! IB=13 2224 (y-axis)
1396 4 : indpg = 11 ! 4 (c4)
1397 1800 : ELSE IF (nc == 4 .AND. ( &
1398 : ib(4) == 39 .OR. &
1399 : ib(4) == 44 .OR. &
1400 : ib(4) == 48)) THEN
1401 : ! IB=14 3839 (z-axis) OR
1402 : ! IB=12 4344 (x-axis) OR
1403 : ! IB=13 4648 (y-axis)
1404 0 : indpg = 12 ! <4>(s4)
1405 1800 : ELSE IF (nc == 8 .AND. ( &
1406 : (ib(3) == 14 .AND. ib(8) == 39) .OR. &
1407 : (ib(3) == 19 .AND. ib(8) == 44) .OR. &
1408 : (ib(3) == 22 .AND. ib(8) == 48))) THEN
1409 : ! IB=14 1415 2825 3839 (z-axis) OR
1410 : ! IB=12 1920 2625 4344 (x-axis) OR
1411 : ! IB=13 2224 2725 4648 (y-axis)
1412 0 : indpg = 13 ! 422(d4)
1413 1800 : ELSE IF (nc == 8 .AND. ib(4) == 4 .AND. ( &
1414 : ib(8) == 16 .OR. &
1415 : ib(8) == 20 .OR. &
1416 : ib(8) == 24)) THEN
1417 : ! IB=12 3 413 1415 16 (z-axis) OR
1418 : ! IB=12 3 417 1920 18 (x-axis) OR
1419 : ! IB=12 3 421 2224 23 (y-axis)
1420 24 : indpg = 14 ! 4/m(c4h)
1421 1776 : ELSE IF (nc == 8 .AND. ( &
1422 : ib(8) == 40 .OR. &
1423 : ib(8) == 42 .OR. &
1424 : ib(8) == 47)) THEN
1425 : ! IB=14 1415 2627 3740 (z-axis) OR
1426 : ! IB=12 1920 2827 4142 (x-axis) OR
1427 : ! IB=13 2224 2628 4547 (y-axis)
1428 4 : indpg = 15 ! 4mm(c4v)
1429 1772 : ELSE IF (nc == 8 .AND. ( &
1430 : (ib(3) == 13 .AND. ib(8) == 39) .OR. &
1431 : (ib(3) == 17 .AND. ib(8) == 44) .OR. &
1432 : (ib(3) == 21 .AND. ib(8) == 48))) THEN
1433 : ! IB=14 1316 2627 3839 (z-axis) OR
1434 : ! IB=12 1718 2827 4344 (x-axis) OR
1435 : ! IB=13 2123 2628 4648 (y-axis)
1436 0 : indpg = 16 ! <4>2m(d2d)
1437 1772 : ELSE IF (nc == 16 .AND. ( &
1438 : ib(16) == 40 .OR. &
1439 : ib(16) == 44 .OR. &
1440 : ib(16) == 48)) THEN
1441 : ! IB=12 3 413 1415 1625 2627 2837 3839 40 (z-axis) OR
1442 : ! IB=12 3 417 1920 1825 2627 2841 4344 42 (x-axis) OR
1443 : ! IB=12 3 421 2224 2325 2627 2845 4648 47 (y-axis)
1444 484 : indpg = 17 ! 4/mmm(d4h)
1445 1288 : ELSE IF (nc == 4 .AND. (ib(4) == 4)) THEN
1446 : ! Orthorhombic system
1447 : ! IB=12 3 4
1448 0 : indpg = 25 ! 222(d2)
1449 1288 : ELSE IF (nc == 4 .AND. ( &
1450 : ib(4) == 27 .OR. &
1451 : ib(4) == 28)) THEN
1452 : ! IB=13 2627 (z-axis) OR
1453 : ! IB=12 2728 (x-axis) OR
1454 : ! IB=14 2628 (y-axis) OR
1455 0 : indpg = 26 ! mm2(c2v)
1456 1288 : ELSE IF (nc == 8) THEN
1457 : ! IB=12 3 425 2627 28
1458 42 : indpg = 27 ! mmm(d2h)
1459 1246 : ELSE IF (nc == 12 .AND. ( &
1460 : ib(12) == 12 .OR. &
1461 : ib(12) == 47 .OR. &
1462 : ib(12) == 45)) THEN
1463 : ! Cubic system
1464 : ! IB=12 3 4 5 6 7 8 910 1112 OR
1465 : ! IB=15 1113 1823 2530 3537 4247 OR
1466 : ! IB=18 1016 1821 2532 3440 4245
1467 0 : indpg = 28 ! 23 (t)
1468 1246 : ELSE IF (nc == 24 .AND. ib(24) == 36) THEN
1469 : ! IB= 1 2 3 4 5 6 7 8 910 1112
1470 : ! 2526 2728 2930 3132 3334 3536
1471 0 : indpg = 29 ! m3 (th)
1472 1246 : ELSE IF (nc == 24 .AND. ib(24) == 24) THEN
1473 : ! IB=12 3 45 6 78 9 1011 12
1474 : ! 1314 1516 1718 1920 2122 2324
1475 0 : indpg = 30 ! 432 (o)
1476 1246 : ELSE IF (nc == 24 .AND. ib(24) == 48) THEN
1477 : ! IB=12 3 45 6 78 9 1011 12
1478 : ! 3738 3940 4142 4345 4647 48
1479 0 : indpg = 31 ! <4>3m(td)
1480 1246 : ELSE IF (nc == 48) THEN
1481 : ! IB=1..48
1482 954 : indpg = 32 ! m3m(oh)
1483 : ELSE
1484 : ! WRITE(6,'(" ATFTM1! IHG=",A," NC=",I2)') ICST(IHG),NC
1485 : ! WRITE(6,'(" ATFTM1!",19I3)') (IB(I),I=1,NC)
1486 : ! WRITE(6,'(" ATFTM1! THIS CASE IS UNKNOWN IN THE DATABASE")')
1487 : ! Probably a sub-group of 32
1488 292 : indpg = -32
1489 : END IF
1490 : ELSE IF (ihg >= 6) THEN
1491 210 : IF (nc == 0) THEN
1492 0 : IF (iout > 0) THEN
1493 0 : WRITE (iout, '(" ATFTM1! IHG=",A," NC=",I2)') icst(ihg), nC
1494 : END IF
1495 0 : CPABORT('ATFTM1: NUMBER OF ROTATION NULL')
1496 : ! Triclinic system
1497 210 : ELSE IF (nc == 1) THEN
1498 : ! IB=1
1499 0 : indpg = 1 ! 1 (c1)
1500 210 : ELSE IF (nc == 2 .AND. ib(2) == 13) THEN
1501 : ! IB=113
1502 0 : indpg = 2 ! <1>(ci)
1503 210 : ELSE IF (nc == 2 .AND. ( &
1504 : ib(2) == 4)) THEN ! 2[001]
1505 : ! Monoclinic system
1506 : ! IB=1 4
1507 0 : indpg = 3 ! 2 (c2)
1508 210 : ELSE IF (nc == 2 .AND. ( &
1509 : ib(2) == 16)) THEN
1510 : ! IB=116
1511 0 : indpg = 4 ! m (c1h)
1512 210 : ELSE IF (nc == 4 .AND. ( &
1513 : ib(4) == 24 .OR. &
1514 : ib(4) == 20)) THEN
1515 : ! IB=112 1324 OR
1516 : ! IB=1 813 20
1517 0 : indpg = 5 ! 2/m(c2h)
1518 210 : ELSE IF (nc == 3 .AND. ib(3) == 5) THEN
1519 : ! Trigonal system
1520 : ! IB=13 5
1521 8 : indpg = 6 ! 3 (c3)
1522 202 : ELSE IF (nc == 6 .AND. ib(6) == 17) THEN
1523 : ! IB=113 1517 35
1524 0 : indpg = 7 ! <3>(c3i)
1525 202 : ELSE IF (nc == 6 .AND. ib(6) == 11) THEN
1526 : ! IB=17 9 1135
1527 0 : indpg = 8 ! 32 (d3)
1528 202 : ELSE IF (nc == 6 .AND. ib(6) == 23) THEN
1529 : ! IB=13 5 1921 23
1530 0 : indpg = 9 ! 3m (c3v)
1531 202 : ELSE IF (nc == 12 .AND. ib(12) == 23) THEN
1532 : ! IB=13 5 79 1113 1517 1921 23
1533 194 : indpg = 10 ! <3>m(d3d)
1534 8 : ELSE IF (nc == 6 .AND. ib(6) == 6) THEN
1535 : ! Hexagonal system
1536 : ! IB=12 3 45 6
1537 0 : indpg = 18 ! 6 (c6)
1538 8 : ELSE IF (nc == 6 .AND. ib(6) == 18) THEN
1539 : ! IB=13 5 1416 18
1540 0 : indpg = 19 ! <6>(c3h)
1541 8 : ELSE IF (nc == 12 .AND. ib(12) == 18) THEN
1542 : ! IB=12 3 45 6 1314 1516 1718
1543 0 : indpg = 20 ! 6/m(c6h)
1544 8 : ELSE IF (nc == 12 .AND. ib(12) == 12) THEN
1545 : ! IB=12 3 45 6 78 9 1011 12
1546 0 : indpg = 21 ! 622(d6)
1547 8 : ELSE IF (nc == 12 .AND. ib(2) == 2 .AND. ib(12) == 24) THEN
1548 : ! IB=12 3 45 6 1920 2122 2324
1549 0 : indpg = 22 ! 6mm(c6v)
1550 8 : ELSE IF (nc == 12 .AND. ib(2) == 3 .AND. ib(12) == 24) THEN
1551 : ! IB=13 5 79 1114 1618 2022 24
1552 0 : indpg = 23 ! <6>m2(d3h)
1553 8 : ELSE IF (nc == 24) THEN
1554 : ! IB=1..24
1555 8 : indpg = 24 ! 6/mmm(d6h)
1556 : ELSE
1557 : ! Probably a sub-group of 24
1558 : ! WRITE(6,'(" ATFTM1! IHG=",A," NC=",I2)') ICST(IHG),NC
1559 : ! WRITE(6,'(" ATFTM1!",48I3)') (IB(I),I=1,NC)
1560 : ! WRITE(6,'(" ATFTM1! THIS CASE IS UNKNOWN IN THE DATABASE")')
1561 0 : indpg = -24
1562 : END IF
1563 : END IF
1564 : ! ==--------------------------------------------------------------==
1565 : ! == Determination if the space group is symmorphic or not ==
1566 : ! ==--------------------------------------------------------------==
1567 2880 : IF (isy /= 1) THEN
1568 : ! Transform V in cartesian coordinates
1569 15852 : DO n = 1, nc
1570 15236 : vc(1, n) = a(1, 1)*v(1, n) + a(1, 2)*v(2, n) + a(1, 3)*v(3, n)
1571 15236 : vc(2, n) = a(2, 1)*v(1, n) + a(2, 2)*v(2, n) + a(2, 3)*v(3, n)
1572 15852 : vc(3, n) = a(3, 1)*v(1, n) + a(3, 2)*v(2, n) + a(3, 3)*v(3, n)
1573 : END DO
1574 616 : CALL symmorphic(nc, ib, r, vc, ai, info, origin, delta)
1575 616 : IF (info == 1) THEN
1576 40 : CALL rlv3(ai, origin, xb, il, delta)
1577 : ! !!!RLV3 determines -XB in crystal coordinates
1578 : ! !!We want between 0.0 and 1.0
1579 160 : DO i = 1, 3
1580 160 : IF (-xb(i) >= 0._dp) THEN
1581 72 : origin(i) = -xb(i)
1582 : ELSE
1583 48 : origin(i) = 1._dp - xb(i)
1584 : END IF
1585 : END DO
1586 160 : DO i = 1, 3
1587 160 : xb(i) = a(i, 1)*origin(1) + a(i, 2)*origin(2) + a(i, 3)*origin(3)
1588 : END DO
1589 40 : isy = -1
1590 576 : ELSE IF (info == 0) THEN
1591 576 : isy = 0
1592 : ELSE
1593 0 : isy = -2
1594 : END IF
1595 : ELSE
1596 9056 : DO i = 1, 3
1597 9056 : origin(i) = 0._dp
1598 : END DO
1599 : END IF
1600 : ! ==--------------------------------------------------------------==
1601 : ! == Output ==
1602 : ! ==--------------------------------------------------------------==
1603 2880 : IF (iout > 0) THEN
1604 : IF (iout > 0) THEN
1605 670 : WRITE (iout, *)
1606 : END IF
1607 670 : CALL xstring(icst(ihg), i, j)
1608 670 : IF ((ihg == 7 .AND. nc == 24) .OR. &
1609 : (ihg == 5 .AND. nc == 48)) THEN
1610 185 : IF (iout > 0) THEN
1611 : WRITE (iout, '(A,A,A)') &
1612 185 : ' KPSYM| THE POINT GROUP OF THE CRYSTAL IS THE FULL ', &
1613 185 : icst(ihg) (i:j), &
1614 370 : ' GROUP'
1615 : END IF
1616 : ELSE
1617 485 : IF (iout > 0) THEN
1618 : WRITE (iout, '(A,A,A,I2,A)') &
1619 485 : ' KPSYM| THE CRYSTAL SYSTEM IS ', &
1620 485 : icst(ihg) (i:j), &
1621 970 : ' WITH ', nc, ' OPERATIONS:'
1622 : END IF
1623 485 : IF (ihc == 0) THEN
1624 18 : IF (iout > 0) THEN
1625 225 : WRITE (iout, '( 5(5(A13),/))') (rname_hexai(ib(i)), i=1, nc)
1626 : END IF
1627 : ELSE
1628 467 : IF (iout > 0) THEN
1629 4302 : WRITE (iout, '(10(5(A13),/))') (rname_cubic(ib(i)), i=1, nc)
1630 : END IF
1631 : END IF
1632 : END IF
1633 : ! ==------------------------------------------------------------==
1634 670 : IF (isy == 1) THEN
1635 566 : IF (iout > 0) THEN
1636 : WRITE (iout, '(A)') &
1637 566 : ' KPSYM| THE SPACE GROUP OF THE CRYSTAL IS SYMMORPHIC'
1638 : END IF
1639 104 : ELSE IF (isy == -1) THEN
1640 8 : IF (iout > 0) THEN
1641 : WRITE (iout, '(A)') &
1642 8 : ' KPSYM| THE SPACE GROUP OF THE CRYSTAL IS SYMMORPHIC'
1643 : END IF
1644 8 : IF (iout > 0) THEN
1645 : WRITE (iout, '(A,A,/,T3,3F10.6,3X,3F10.6)') &
1646 8 : ' KPSYM| THE STANDARD ORIGIN OF COORDINATES IS: ', &
1647 16 : '[CARTESIAN] [CRYSTAL]', xb, origin
1648 : END IF
1649 96 : ELSE IF (isy == 0) THEN
1650 96 : IF (iout > 0) THEN
1651 : WRITE (iout, '(A,/,3X,A,F15.6,A)') &
1652 96 : ' KPSYM| THE SPACE GROUP IS NON-SYMMORPHIC,', &
1653 192 : ' (SUM OF TRANSLATION VECTORS=', vs, ')'
1654 : END IF
1655 0 : ELSE IF (isy == -2) THEN
1656 0 : IF (iout > 0) THEN
1657 : WRITE (iout, '(A,A)') &
1658 0 : ' KPSYM| CANNOT DETERMINE IF THE SPACE GROUP IS', &
1659 0 : ' SYMMORPHIC OR NOT'
1660 : END IF
1661 0 : IF (iout > 0) THEN
1662 : WRITE (iout, '(A,/,A,/,3X,A,F15.6,A)') &
1663 0 : ' KPSYM| THE SPACE GROUP IS NON-SYMMORPHIC,', &
1664 0 : ' KPSYM| OR ELSE A NON STANDARD ORIGIN OF COORDINATES WAS USED.', &
1665 0 : ' KPSYM| (SUM OF TRANSLATION VECTORS=', vs, ')'
1666 : END IF
1667 : END IF
1668 670 : IF (indpg > 0) THEN
1669 623 : CALL xstring(pgrp(indpg), i, j)
1670 623 : CALL xstring(pgrd(indpg), k, l)
1671 623 : IF (iout > 0) THEN
1672 : WRITE (iout, '(A,A,"(",A,")",T56,"[INDEX=",I2,"]")') &
1673 623 : ' KPSYM| THE POINT GROUP OF THE CRYSTAL IS ', pgrp(indpg) (i:j), &
1674 1246 : pgrd(indpg) (k:l), indpg
1675 : END IF
1676 : ELSE
1677 47 : CALL xstring(pgrp(-indpg), i, j)
1678 47 : CALL xstring(pgrd(-indpg), k, l)
1679 47 : IF (iout > 0) THEN
1680 : WRITE (iout, '(A,I2,A,A,"(",A,")",T56,"[INDEX=",I2,"]")') &
1681 47 : ' KPSYM| POINT GROUP: GROUP ORDER=', nc, &
1682 47 : ' SUBGROUP OF ', pgrp(-indpg) (i:j), &
1683 94 : pgrd(-indpg) (k:l), -indpg
1684 : END IF
1685 : END IF
1686 670 : IF (ntvec == 1) THEN
1687 610 : IF (iout > 0) THEN
1688 : WRITE (iout, '(A,T60,I6)') &
1689 610 : ' KPSYM| NUMBER OF PRIMITIVE CELL:', ntvec
1690 : END IF
1691 : ELSE
1692 60 : IF (iout > 0) THEN
1693 : WRITE (iout, '(A,T60,I6)') &
1694 60 : ' KPSYM| NUMBER OF PRIMITIVE CELLS:', ntvec
1695 : END IF
1696 : END IF
1697 : END IF
1698 :
1699 2880 : END SUBROUTINE atftm1
1700 :
1701 : ! **************************************************************************************************
1702 : !> \brief ...
1703 : !> \param n ...
1704 : !> \param nat ...
1705 : !> \param ty ...
1706 : !> \param rx ...
1707 : !> \param x ...
1708 : !> \param vr ...
1709 : !> \param f0 ...
1710 : !> \param ai ...
1711 : !> \param isc ...
1712 : !> \param nodupli ...
1713 : !> \param oksym ...
1714 : !> \param delta ...
1715 : ! **************************************************************************************************
1716 223340 : SUBROUTINE checkrlv3(n, nat, ty, rx, x, vr, f0, ai, isc, &
1717 : nodupli, oksym, delta)
1718 : ! ==--------------------------------------------------------------==
1719 : ! == WRITTEN IN MAY 14TH, 1998 (T.D.) ==
1720 : ! == CHECK IF RX+VR GIVES THE SAME LATTICE AS X ==
1721 : ! == BUILD THE ATOM TRANSFORMATION TABLE ==
1722 : ! ==--------------------------------------------------------------==
1723 : ! == INPUT: ==
1724 : ! == N ROTATION NUMBER (INDEX USED IN F0 BETWEEN 1 AND 48) ==
1725 : ! == NAT NUMBER OF ATOMS ==
1726 : ! == TY(1:NAT) TYPE OF ATOMS ==
1727 : ! == RX(1:3,1:NAT) ATOMIC COORDINATES FROM Nth ROTATION (CART.) ==
1728 : ! == X(1:3,1:NAT) ATOMIC COORDINATES (CARTESIAN) ==
1729 : ! == VR(1:3) TRANSLATION VECTOR (CRYSTAL COOR.) ==
1730 : ! == AI(1:3,1:3) LATTICE RECIPROCAL VECTORS ==
1731 : ! == NODUPLI .TRUE., THE CELL IS A PRIMITIVE ONE ==
1732 : ! == WE CAN SPEED UP ==
1733 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
1734 : ! == OUTPUT: ==
1735 : ! == F0(1:49,1:NAT) ATOM TRANSFORMATION TABLE ==
1736 : ! == F0 IS THE FUNCTION DEFINED IN MARADUDIN AND VOSK0 ==
1737 : ! == BY EQ.(2.35). ==
1738 : ! == IT DEFINES THE ATOM TRANSFORMATION TABLE ==
1739 : ! == OKSYM TRUE IF RX+VR = X ==
1740 : ! == ISC(1:NAT) SCRATCH ARRAY ==
1741 : ! == USED TO SPEED UP THE ROUTINE ==
1742 : ! == EACH ATOM IS ONLY ONCE AN IMAGE ==
1743 : ! == IF NO DUPLICATION OF THE CELL ==
1744 : ! ==--------------------------------------------------------------==
1745 : INTEGER :: n, nat, ty(nat)
1746 : REAL(dp) :: rx(3, nat), x(3, nat), vr(3)
1747 : INTEGER :: f0(49, nat)
1748 : REAL(dp) :: ai(3, 3)
1749 : INTEGER :: isc(nat)
1750 : LOGICAL :: nodupli, oksym
1751 : REAL(dp) :: delta
1752 :
1753 : INTEGER :: ia, ib, il
1754 : REAL(dp) :: tol, vt(3), xb(3)
1755 :
1756 : ! Fractional residuals at the tolerance boundary accumulate Cartesian
1757 : ! rotation and lattice-conversion roundoff. Account for that roundoff so
1758 : ! equivalent operations are not accepted or rejected by a few ulps.
1759 : tol = delta + 32.0_dp*EPSILON(1.0_dp)*MAX(1.0_dp, &
1760 15827624 : MAXVAL(SUM(ABS(ai), DIM=2))*MAX(MAXVAL(ABS(rx)), MAXVAL(ABS(x))))
1761 :
1762 1810948 : DO ia = 1, nat
1763 1810948 : isc(ia) = 0
1764 : END DO
1765 : ! Now we check if ROT(N)+VR gives a correct symmetry.
1766 669980 : atom: DO ia = 1, nat
1767 3241556 : DO ib = 1, nat
1768 3241556 : IF (ty(ia) == ty(ib) .AND. isc(ib) == 0) THEN
1769 2482640 : xb(1) = rx(1, ia) - x(1, ib)
1770 2482640 : xb(2) = rx(2, ia) - x(2, ib)
1771 2482640 : xb(3) = rx(3, ia) - x(3, ib)
1772 2482640 : CALL rlv3(ai, xb, vt, il, delta)
1773 : ! VT STANDS FOR V-TEST
1774 4103832 : oksym = ALL(ABS((vr - vt) - ANINT(vr - vt)) <= tol)
1775 2482640 : IF (oksym) THEN
1776 446640 : IF (nodupli) isc(ib) = 1
1777 446640 : f0(n, ia) = ib
1778 : ! IR+VR is the good one: another symmetry operation
1779 : ! Next atom
1780 : CYCLE atom
1781 : END IF
1782 : END IF
1783 : END DO
1784 : ! VR is not the correct translation vector
1785 223340 : RETURN
1786 : END DO atom
1787 : END SUBROUTINE checkrlv3
1788 : ! ==================================================================
1789 : ! **************************************************************************************************
1790 : !> \brief ...
1791 : !> \param nc ...
1792 : !> \param ib ...
1793 : !> \param r ...
1794 : !> \param v ...
1795 : !> \param ai ...
1796 : !> \param info ...
1797 : !> \param origin ...
1798 : !> \param delta ...
1799 : ! **************************************************************************************************
1800 616 : SUBROUTINE symmorphic(nc, ib, r, v, ai, info, origin, delta)
1801 : ! ==--------------------------------------------------------------==
1802 : ! == Check if the group is symmorphic with a non-standard origin ==
1803 : ! == WARNING: If there are equivalent atoms, this routine could ==
1804 : ! == not determine if the space group is symmorphic ==
1805 : ! == So you have to check if the solution V=0 works (see ATFTM1) ==
1806 : ! ==--------------------------------------------------------------==
1807 : ! == INPUT: ==
1808 : ! == NC Number of operations ==
1809 : ! == IB(NC) Index of operation in R ==
1810 : ! == R(3,3,48) Rotations ==
1811 : ! == V(3,NC) Fractional translations related to R(3,3,IB(NC)) ==
1812 : ! == R AND V ARE IN CARTESIAN COORDINATES ==
1813 : ! == AI(I,J) ARE THE RECIPROCAL LATTICE VECTORS, ==
1814 : ! == B(I) = AI(I,J),J=1,2,3 ==
1815 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
1816 : ! == ==
1817 : ! == OUTPUT: ==
1818 : ! == ORIGIN(1:3) Give standard origin (cartesian coordinates) ==
1819 : ! == Give the standard origin with smallest coordinates==
1820 : ! == if NTVEC /= 1 ==
1821 : ! == INFO = 1 The group is symmorphic ==
1822 : ! == INFO = 0 The group is not symmorphic ==
1823 : ! == INFO =-1 The routine cannot determine ==
1824 : ! ==--------------------------------------------------------------==
1825 : INTEGER :: nc, ib(nc)
1826 : REAL(dp) :: r(3, 3, 48), v(3, nc), ai(3, 3)
1827 : INTEGER :: info
1828 : REAL(dp) :: origin(3), delta
1829 :
1830 : INTEGER :: i, i1, ierror, igood(3), il, imissing2, &
1831 : imissing3, iok(3), ionly, ir, j, j1
1832 : REAL(dp) :: diag, dif, r2(2, 2), r3(3, 3), vr(3), &
1833 : xb(3)
1834 :
1835 : ! Variables
1836 : ! ==--------------------------------------------------------------==
1837 : ! Find a point A / V_R = (1-R).OA
1838 :
1839 2464 : DO i = 1, 3
1840 2464 : iok(i) = 0
1841 : END DO
1842 2464 : DO i = 1, 3
1843 2464 : origin(i) = 0._dp
1844 : END DO
1845 5268 : DO ir = 1, nc
1846 5212 : dif = v(1, ir)*v(1, ir) + v(2, ir)*v(2, ir) + v(3, ir)*v(3, ir)
1847 5268 : IF (dif > delta*delta) THEN
1848 5696 : DO i = 1, 3
1849 5696 : igood(i) = 1
1850 : END DO
1851 : ! V is non-zero. Construct matrix 1-R
1852 5696 : DO i = 1, 3
1853 17088 : DO j = 1, 3
1854 17088 : r3(i, j) = -r(i, j, ib(ir))
1855 : END DO
1856 5696 : r3(i, i) = 1 + r3(i, i)
1857 : END DO
1858 1424 : CALL invmat(r3, ierror)
1859 1424 : IF (ierror == 0) THEN
1860 : ! The matrix 3x3 has an inverse.
1861 912 : DO i = 1, 3
1862 : vr(i) = r3(i, 1)*v(1, ir) &
1863 : + r3(i, 2)*v(2, ir) &
1864 912 : + r3(i, 3)*v(3, ir)
1865 : END DO
1866 : ELSE
1867 : ! IERROR gives the column which causes some trouble
1868 : ! Construct matrix 1-R with 2x2
1869 1196 : igood(ierror) = 0
1870 1196 : imissing3 = ierror
1871 1196 : i1 = 0
1872 4784 : DO i = 1, 3
1873 4784 : IF (i /= ierror) THEN
1874 2392 : i1 = i1 + 1
1875 2392 : j1 = 0
1876 9568 : DO j = 1, 3
1877 9568 : IF (j /= ierror) THEN
1878 4784 : j1 = j1 + 1
1879 4784 : r2(i1, j1) = -r(i, j, ib(ir))
1880 : END IF
1881 : END DO
1882 2392 : r2(i1, i1) = 1 + r2(i1, i1)
1883 : END IF
1884 : END DO
1885 1196 : CALL invmat(r2, ierror)
1886 1196 : IF (ierror == 0) THEN
1887 : ! The matrix 2X2 has an inverse.
1888 : ! Solve Vxy = (1-R).OAxy + OAz R3z (z is IMISSING3)
1889 : i1 = 0
1890 4080 : DO i = 1, 3
1891 4080 : IF (igood(i) == 1) THEN
1892 2040 : i1 = i1 + 1
1893 2040 : vr(i) = 0._dp
1894 2040 : j1 = 0
1895 8160 : DO j = 1, 3
1896 8160 : IF (igood(j) == 1) THEN
1897 4080 : j1 = j1 + 1
1898 : vr(i) = vr(i) + r2(i1, j1)*(v(j, ir) + &
1899 4080 : origin(imissing3)*r(j, imissing3, ib(ir)))
1900 : END IF
1901 : END DO
1902 : ELSE
1903 1020 : vr(i) = origin(i)
1904 : END IF
1905 : END DO
1906 : ELSE
1907 : ! Construct matrix 1-R with 1x1
1908 : i1 = 0
1909 704 : DO i = 1, 3
1910 704 : IF (i /= imissing3) THEN
1911 352 : i1 = i1 + 1
1912 352 : IF (i1 == ierror) THEN
1913 176 : igood(i) = 0
1914 176 : imissing2 = i
1915 : ELSE
1916 : ionly = i
1917 : END IF
1918 : END IF
1919 : END DO
1920 176 : diag = (1 - r(ionly, ionly, ib(ir)))
1921 176 : IF (ABS(diag) > delta) THEN
1922 : vr(ionly) = 1._dp/diag*(v(ionly, ir) + &
1923 : origin(imissing3)*r(ionly, imissing3, ib(ir)) + &
1924 176 : origin(imissing2)*r(ionly, imissing2, ib(ir)))
1925 : ELSE
1926 0 : vr(ionly) = origin(ionly)
1927 0 : igood(ionly) = 0
1928 : END IF
1929 176 : vr(imissing3) = origin(imissing3)
1930 176 : vr(imissing2) = origin(imissing2)
1931 : END IF
1932 : END IF
1933 : ! ==----------------------------------------------------------==
1934 : ! Compare VR with ORIGIN
1935 1424 : dif = 0._dp
1936 : ! If NTVEC /=1 there are NTVEC possible standard origins
1937 5696 : DO i = 1, 3
1938 5696 : IF (iok(i) == 1) THEN
1939 1832 : dif = dif + ABS(origin(i) - vr(i))
1940 : END IF
1941 : END DO
1942 1424 : IF (dif > delta) THEN
1943 : ! Non-symmorphic
1944 560 : info = 0
1945 560 : RETURN
1946 : ELSE
1947 3456 : DO i = 1, 3
1948 3456 : IF (iok(i) /= 1 .AND. igood(i) == 1) THEN
1949 1464 : iok(i) = 1
1950 1464 : origin(i) = vr(i)
1951 : END IF
1952 : END DO
1953 : END IF
1954 : END IF
1955 : END DO
1956 : ! ==--------------------------------------------------------------==
1957 56 : IF (iok(1) == 0 .AND. iok(2) == 0 .AND. iok(3) == 0) THEN
1958 : ! Cannot not determine
1959 0 : info = -1
1960 0 : RETURN
1961 : END IF
1962 : ! The group is symmorphic
1963 56 : info = 1
1964 : ! Check
1965 168 : DO ir = 1, nc
1966 512 : DO i = 1, 3
1967 : vr(i) = r(i, 1, ib(ir))*origin(1) &
1968 : + r(i, 2, ib(ir))*origin(2) &
1969 384 : + r(i, 3, ib(ir))*origin(3)
1970 512 : vr(i) = (origin(i) - vr(i)) - v(i, ir)
1971 : END DO
1972 128 : CALL rlv3(ai, vr, xb, il, delta)
1973 128 : dif = ABS(xb(1)) + ABS(xb(2)) + ABS(xb(3))
1974 168 : IF (dif > delta) THEN
1975 : ! Non-symmorphic
1976 16 : info = 0
1977 16 : RETURN
1978 : END IF
1979 : END DO
1980 : ! ==--------------------------------------------------------------==
1981 : RETURN
1982 : END SUBROUTINE symmorphic
1983 : ! ==================================================================
1984 : ! **************************************************************************************************
1985 : !> \brief ...
1986 : !> \param ihc ...
1987 : !> \param r ...
1988 : ! **************************************************************************************************
1989 6294 : SUBROUTINE rot1(ihc, r)
1990 : ! ==--------------------------------------------------------------==
1991 : ! == WRITTEN ON FEBRUARY 17TH, 1976 ==
1992 : ! == GENERATION OF THE X,Y,Z-TRANSFORMATION MATRICES 3X3 ==
1993 : ! == FOR HEXAGONAL AND CUBIC GROUPS ==
1994 : ! == SUBROUTINES NEEDED -- NONE ==
1995 : ! ==--------------------------------------------------------------==
1996 : ! == THIS IS IDENTICAL WITH THE SUBROUTINE ROT OF WORLTON-WARREN ==
1997 : ! == (IN THE AC-COMPLEX), ONLY THE WAY OF TRANSFERRING THE DATA ==
1998 : ! == WAS CHANGED ==
1999 : ! ==--------------------------------------------------------------==
2000 : ! == INPUT DATA: ==
2001 : ! == IHC SWITCH DETERMINING IF WE DESIRE ==
2002 : ! == THE HEXAGONAL GROUP(IHC=0) OR THE CUBIC GROUP (IHC=1) ==
2003 : ! == OUTPUT DATA: ==
2004 : ! == R...THE 3X3 MATRICES OF THE DESIRED COORDINATE REPRESENTATION==
2005 : ! == THEIR NUMBERING CORRESPONDS TO THE SYMMETRY ELEMENTS AS ==
2006 : ! == LISTE IN WORLTON-WARREN ==
2007 : ! == (COMPUT. PHYS. COMM. 3(1972) 88--117) ==
2008 : ! == FOR IHC=0 THE FIRST 24 MATRICES OF THE ARRAY R REPRESENT ==
2009 : ! == THE FULL HEXAGONAL GROUP D(6H) ==
2010 : ! == FOR IHC=1 THE FIRST 48 MATRICES OF THE ARRAY R REPRESENT ==
2011 : ! == THE FULL CUBIC GROUP O(H) ==
2012 : ! ==--------------------------------------------------------------==
2013 : INTEGER :: ihc
2014 : REAL(dp) :: r(3, 3, 48)
2015 :
2016 : INTEGER :: i, j, k, n, nv
2017 : REAL(dp) :: c, s
2018 :
2019 25176 : DO j = 1, 3
2020 81822 : DO i = 1, 3
2021 2794536 : DO n = 1, 48
2022 2775654 : r(i, j, n) = 0._dp
2023 : END DO
2024 : END DO
2025 : END DO
2026 6294 : IF (ihc == 0) THEN
2027 : ! ==------------------------------------------------------------==
2028 : ! DEFINE THE GENERATORS FOR THE ROTATION MATRICES--HEXAGONAL GROUP
2029 : ! ==------------------------------------------------------------==
2030 3252 : c = 0.5_dp
2031 3252 : s = 0.5_dp*SQRT(3.0_dp)
2032 3252 : r(1, 1, 2) = c
2033 3252 : r(1, 2, 2) = -s
2034 3252 : r(2, 1, 2) = s
2035 3252 : r(2, 2, 2) = c
2036 3252 : r(1, 1, 7) = -c
2037 3252 : r(1, 2, 7) = -s
2038 3252 : r(2, 1, 7) = -s
2039 3252 : r(2, 2, 7) = c
2040 22764 : DO n = 1, 6
2041 19512 : r(3, 3, n) = 1._dp
2042 19512 : r(3, 3, n + 18) = 1._dp
2043 19512 : r(3, 3, n + 6) = -1._dp
2044 22764 : r(3, 3, n + 12) = -1._dp
2045 : END DO
2046 : ! ==------------------------------------------------------------==
2047 : ! == GENERATE THE REST OF THE ROTATION MATRICES ==
2048 : ! ==------------------------------------------------------------==
2049 9756 : DO i = 1, 2
2050 6504 : r(i, i, 1) = 1._dp
2051 22764 : DO j = 1, 2
2052 13008 : r(i, j, 6) = r(j, i, 2)
2053 45528 : DO k = 1, 2
2054 26016 : r(i, j, 3) = r(i, j, 3) + r(i, k, 2)*r(k, j, 2)
2055 26016 : r(i, j, 8) = r(i, j, 8) + r(i, k, 2)*r(k, j, 7)
2056 39024 : r(i, j, 12) = r(i, j, 12) + r(i, k, 7)*r(k, j, 2)
2057 : END DO
2058 : END DO
2059 : END DO
2060 9756 : DO i = 1, 2
2061 22764 : DO j = 1, 2
2062 13008 : r(i, j, 5) = r(j, i, 3)
2063 45528 : DO k = 1, 2
2064 26016 : r(i, j, 4) = r(i, j, 4) + r(i, k, 2)*r(k, j, 3)
2065 26016 : r(i, j, 9) = r(i, j, 9) + r(i, k, 2)*r(k, j, 8)
2066 26016 : r(i, j, 10) = r(i, j, 10) + r(i, k, 12)*r(k, j, 3)
2067 39024 : r(i, j, 11) = r(i, j, 11) + r(i, k, 12)*r(k, j, 2)
2068 : END DO
2069 : END DO
2070 : END DO
2071 42276 : DO n = 1, 12
2072 39024 : nv = n + 12
2073 120324 : DO i = 1, 2
2074 273168 : DO j = 1, 2
2075 234144 : r(i, j, nv) = -r(i, j, n)
2076 : END DO
2077 : END DO
2078 : END DO
2079 : ELSE
2080 : ! ==------------------------------------------------------------==
2081 : ! == DEFINE THE GENERATORS FOR THE ROTATION MATRICES-CUBIC GROUP==
2082 : ! ==------------------------------------------------------------==
2083 3042 : r(1, 3, 9) = 1._dp
2084 3042 : r(2, 1, 9) = 1._dp
2085 3042 : r(3, 2, 9) = 1._dp
2086 3042 : r(1, 1, 19) = 1._dp
2087 3042 : r(2, 3, 19) = -1._dp
2088 3042 : r(3, 2, 19) = 1._dp
2089 12168 : DO i = 1, 3
2090 9126 : r(i, i, 1) = 1._dp
2091 39546 : DO j = 1, 3
2092 27378 : r(i, j, 20) = r(j, i, 19)
2093 27378 : r(i, j, 5) = r(j, i, 9)
2094 118638 : DO k = 1, 3
2095 82134 : r(i, j, 2) = r(i, j, 2) + r(i, k, 19)*r(k, j, 19)
2096 82134 : r(i, j, 16) = r(i, j, 16) + r(i, k, 9)*r(k, j, 19)
2097 109512 : r(i, j, 23) = r(i, j, 23) + r(i, k, 19)*r(k, j, 9)
2098 : END DO
2099 : END DO
2100 : END DO
2101 12168 : DO i = 1, 3
2102 39546 : DO j = 1, 3
2103 118638 : DO k = 1, 3
2104 82134 : r(i, j, 6) = r(i, j, 6) + r(i, k, 2)*r(k, j, 5)
2105 82134 : r(i, j, 7) = r(i, j, 7) + r(i, k, 16)*r(k, j, 23)
2106 82134 : r(i, j, 8) = r(i, j, 8) + r(i, k, 5)*r(k, j, 2)
2107 82134 : r(i, j, 10) = r(i, j, 10) + r(i, k, 2)*r(k, j, 9)
2108 82134 : r(i, j, 11) = r(i, j, 11) + r(i, k, 9)*r(k, j, 2)
2109 82134 : r(i, j, 12) = r(i, j, 12) + r(i, k, 23)*r(k, j, 16)
2110 82134 : r(i, j, 14) = r(i, j, 14) + r(i, k, 16)*r(k, j, 2)
2111 82134 : r(i, j, 15) = r(i, j, 15) + r(i, k, 2)*r(k, j, 16)
2112 82134 : r(i, j, 22) = r(i, j, 22) + r(i, k, 23)*r(k, j, 2)
2113 109512 : r(i, j, 24) = r(i, j, 24) + r(i, k, 2)*r(k, j, 23)
2114 : END DO
2115 : END DO
2116 : END DO
2117 12168 : DO i = 1, 3
2118 39546 : DO j = 1, 3
2119 118638 : DO k = 1, 3
2120 82134 : r(i, j, 3) = r(i, j, 3) + r(i, k, 5)*r(k, j, 12)
2121 82134 : r(i, j, 4) = r(i, j, 4) + r(i, k, 5)*r(k, j, 10)
2122 82134 : r(i, j, 13) = r(i, j, 13) + r(i, k, 23)*r(k, j, 11)
2123 82134 : r(i, j, 17) = r(i, j, 17) + r(i, k, 16)*r(k, j, 12)
2124 82134 : r(i, j, 18) = r(i, j, 18) + r(i, k, 16)*r(k, j, 10)
2125 109512 : r(i, j, 21) = r(i, j, 21) + r(i, k, 12)*r(k, j, 15)
2126 : END DO
2127 : END DO
2128 : END DO
2129 76050 : DO n = 1, 24
2130 73008 : nv = n + 24
2131 73008 : r(1, 1, nv) = -r(1, 1, n)
2132 73008 : r(1, 2, nv) = -r(1, 2, n)
2133 73008 : r(1, 3, nv) = -r(1, 3, n)
2134 73008 : r(2, 1, nv) = -r(2, 1, n)
2135 73008 : r(2, 2, nv) = -r(2, 2, n)
2136 73008 : r(2, 3, nv) = -r(2, 3, n)
2137 73008 : r(3, 1, nv) = -r(3, 1, n)
2138 73008 : r(3, 2, nv) = -r(3, 2, n)
2139 76050 : r(3, 3, nv) = -r(3, 3, n)
2140 : END DO
2141 : END IF
2142 : ! ==--------------------------------------------------------------==
2143 6294 : RETURN
2144 : END SUBROUTINE rot1
2145 : ! ==================================================================
2146 : ! **************************************************************************************************
2147 : !> \brief ...
2148 : !> \param iout ...
2149 : !> \param iq1 ...
2150 : !> \param iq2 ...
2151 : !> \param iq3 ...
2152 : !> \param wvk0 ...
2153 : !> \param nkpoint ...
2154 : !> \param a1 ...
2155 : !> \param a2 ...
2156 : !> \param a3 ...
2157 : !> \param b1 ...
2158 : !> \param b2 ...
2159 : !> \param b3 ...
2160 : !> \param inv ...
2161 : !> \param nc ...
2162 : !> \param ib ...
2163 : !> \param r ...
2164 : !> \param ntot ...
2165 : !> \param wvkl ...
2166 : !> \param lwght ...
2167 : !> \param lrot ...
2168 : !> \param ncbrav ...
2169 : !> \param ibrav ...
2170 : !> \param istriz ...
2171 : !> \param nhash ...
2172 : !> \param includ ...
2173 : !> \param list ...
2174 : !> \param rlist ...
2175 : !> \param delta ...
2176 : ! **************************************************************************************************
2177 960 : SUBROUTINE sppt2(iout, iq1, iq2, iq3, wvk0, nkpoint, &
2178 : a1, a2, a3, b1, b2, b3, &
2179 960 : inv, nc, ib, r, ntot, wvkl, lwght, lrot, &
2180 : ncbrav, ibrav, istriz, &
2181 960 : nhash, includ, list, rlist, delta)
2182 : ! ==--------------------------------------------------------------==
2183 : ! == WRITTEN ON SEPTEMBER 12-20TH, 1979 BY K.K. ==
2184 : ! == MODIFIED 26-MAY-82 BY OLE HOLM NIELSEN ==
2185 : ! == GENERATION OF SPECIAL POINTS FOR AN ARBITRARY LATTICE, ==
2186 : ! == FOLLOWING THE METHOD MONKHORST,PACK, ==
2187 : ! == PHYS. REV. B13 (1976) 5188 ==
2188 : ! == MODIFIED BY MACDONALD, PHYS. REV. B18 (1978) 5897 ==
2189 : ! == THE SUBROUTINE IS WRITTEN ASSUMING THAT THE POINTS ARE ==
2190 : ! == GENERATED IN THE RECIPROCAL SPACE. ==
2191 : ! == IF, HOWEVER, THE B1,B2,B3 ARE REPLACED BY A1,A2,A3, THEN ==
2192 : ! == SPECIAL POINTS IN THE DIRECT SPACE CAN BE PRODUCED, AS WELL. ==
2193 : ! == (NO MULTIPLICATION BY 2PI IS THEN NECESSARY.) ==
2194 : ! == IN THE CASE OF NONSYMMORPHIC GROUPS, THE APPLICATION IN THE ==
2195 : ! == DIRECT SPACE WOULD PROBABLY REQUIRE A CERTAIN CAUTION. ==
2196 : ! == SUBROUTINES NEEDED: BZDEFI,BZRDUC,INBZ,MESH ==
2197 : ! == IN THE CASES WHERE THE POINT GROUP OF THE CRYSTAL DOES NOT ==
2198 : ! == CONTAIN INVERSION. THE LATTER MAY BE ADDED IF WE WISH ==
2199 : ! == (SEE COMMENT TO THE SWITCH INV). ==
2200 : ! == REDUCTION TO THE 1ST BRILLOUIN ZONE IS DONE ==
2201 : ! == BY ADDING G-VECTORS TO FIND THE SHORTEST WAVE-VECTOR. ==
2202 : ! == THE ROTATIONS OF THE BRAVAIS LATTICE ARE APPLIED TO THE ==
2203 : ! == MONKHORST/PACK MESH IN ORDER TO FIND ALL K-POINTS ==
2204 : ! == THAT ARE RELATED BY SYMMETRY. (OLE HOLM NIELSEN) ==
2205 : ! ==--------------------------------------------------------------==
2206 : ! == INPUT DATA: ==
2207 : ! == IOUT: LOGICAL UNIT FOR OUTPUT ==
2208 : ! == IF (IOUT<=0) NO MESSAGE ==
2209 : ! == IQ1,IQ2,IQ3 .. PARAMETER Q OF MONKHORST AND PACK, ==
2210 : ! == GENERALIZED AND DIFFERENT FOR THE 3 DIRECTIONS B1, ==
2211 : ! == B2 AND B3 ==
2212 : ! == WVK0 ... THE 'ARBITRARY' SHIFT OF THE WHOLE MESH, DENOTED K0 ==
2213 : ! == IN MACDONALD. WVK0 = 0 CORRESPONDS TO THE ORIGINAL ==
2214 : ! == SCHEME OF MONKHORST AND PACK. ==
2215 : ! == UNITS: 2PI/(UNITS OF LENGTH USED IN A1, A2, A3), ==
2216 : ! == I.E. THE SAME UNITS AS THE GENERATED SPECIAL POINTS==
2217 : ! == NKPOINT .. VARIABLE DIMENSION OF THE (OUTPUT) ARRAYS WVKL, ==
2218 : ! == LWGHT,LROT, I.E. SPACE RESERVED FOR THE SPECIAL ==
2219 : ! == POINTS AND ACCESSORIES. ==
2220 : ! == NKPOINT HAS TO BE >= NTOT (TOTAL NUMBER OF SPECIAL==
2221 : ! == POINTS. THIS IS CHECKED BY THE SUBROUTINE. ==
2222 : ! == ISTRIZ . INDICATES WHETHER ADDITIONAL MESH POINTS SHOULD BE ==
2223 : ! == GENERATED BY APPLYING GROUP OPERATIONS TO THE MESH. ==
2224 : ! == ISTRIZ=+1 MEANS SYMMETRIZE ==
2225 : ! == ISTRIZ=-1 MEANS DO NOT SYMMETRIZE ==
2226 : ! == THE FOLLOWING INPUT DATA MAY BE OBTAINED FROM THE SBRT. ==
2227 : ! == B1,B2,B3 .. RECIPROCAL LATTICE VECTORS, NOT MULTIPLIED BY ==
2228 : ! == GROUP1: ANY 2PI (IN UNITS RECIPROCAL TO THOSE ==
2229 : ! == OF A1,A2,A3) ==
2230 : ! == INV .... CODE INDICATING WHETHER WE WISH TO ADD THE INVERSION==
2231 : ! == TO THE POINT GROUP OF THE CRYSTAL OR NOT (IN THE ==
2232 : ! == CASE THAT THE POINT GROUP DOES NOT CONTAIN ANY). ==
2233 : ! == INV=0 MEANS: DO NOT ADD INVERSION ==
2234 : ! == INV/=0 MEANS: ADD THE INVERSION ==
2235 : ! == INV/=0 SHOULD BE THE STANDARD CHOICE WHEN SPPT2 ==
2236 : ! == IS USED IN RECIPROCAL SPACE - IN ORDER TO MAKE ==
2237 : ! == USE OF THE HERMITICITY OF HAMILTONIAN. ==
2238 : ! == WHEN USED IN DIRECT SPACE, THE RIGHT CHOICE OF INV ==
2239 : ! == WILL DEPEND ON THE NATURE OF THE PHYSICAL PROBLEM. ==
2240 : ! == IN THE CASES WHERE THE INVERSION IS ADDED BY THE ==
2241 : ! == SWITCH INV, THE LIST IB WILL NOT BE MODIFIED BUT IN ==
2242 : ! == THE OUTPUT LIST LROT SOME OF THE OPERATIONS WILL ==
2243 : ! == APPEAR WITH NEGATIVE SIGN; THIS MEANS THAT THEY HAVE==
2244 : ! == TO BE APPLIED MULTIPLIED BY INVERSION. ==
2245 : ! == NC ..... TOTAL NUMBER OF ELEMENTS IN THE POINT GROUP OF THE ==
2246 : ! == CRYSTAL ==
2247 : ! == IB ..... LIST OF THE ROTATIONS CONSTITUTING THE POINT GROUP ==
2248 : ! == OF THE CRYSTAL. THE NUMBERING IS THAT DEFINED IN ==
2249 : ! == WORLTON AND WARREN, I.E. THE ONE MATERIALIZED IN THE==
2250 : ! == ARRAY R (SEE BELOW) ==
2251 : ! == ONLY THE FIRST NC ELEMENTS OF THE ARRAY IB ARE ==
2252 : ! == MEANINGFUL ==
2253 : ! == R ...... LIST OF THE 3 X 3 ROTATION MATRICES ==
2254 : ! == (XYZ REPRESENTATION OF THE O(H) OR D(6)H GROUPS) ==
2255 : ! == ALL 48 OR 24 MATRICES ARE LISTED. ==
2256 : ! == NCBRAV . TOTAL NUMBER OF ELEMENTS IN RBRAV ==
2257 : ! == IBRAV .. LIST OF NCBRAV OPERATIONS OF THE BRAVAIS LATTICE ==
2258 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
2259 : ! ==--------------------------------------------------------------==
2260 : ! == OUTPUT DATA: ==
2261 : ! == NTOT ... TOTAL NUMBER OF SPECIAL POINTS ==
2262 : ! == IF NTOT APPEARS NEGATIVE, THIS IS AN ERROR SIGNAL ==
2263 : ! == WHICH MEANS THAT THE DIMENSION NKPOINT WAS CHOSEN ==
2264 : ! == TOO SMALL SO THAT THE ARRAYS WVKL ETC. CANNOT ==
2265 : ! == ACCOMODATE ALL THE GENERATED SPECIAL POINTS. ==
2266 : ! == IN THIS CASE THE ARRAYS WILL BE FILLED UP TO NKPOINT==
2267 : ! == AND FURTHER GENERATION OF NEW POINTS WILL BE ==
2268 : ! == INTERRUPTED. ==
2269 : ! == WVKL ... LIST OF SPECIAL POINTS. ==
2270 : ! == CARTESIAN COORDINATES AND NOT MULTIPLIED BY 2*PI. ==
2271 : ! == ONLY THE FIRST NTOT VECTORS ARE MEANINGFUL ==
2272 : ! == ALTHOUGH NO 2 POINTS FROM THE LIST ARE EQUIVALENT ==
2273 : ! == BY SYMMETRY, THIS SUBROUTINE STILL HAS A KIND OF ==
2274 : ! == 'BEAUTY DEFECT': THE POINTS FINALLY ==
2275 : ! == SELECTED ARE NOT NECESSARILY SITUATED IN A ==
2276 : ! == 'COMPACT' IRREDUCIBLE BRILL.ZONE; THEY MIGHT LIE IN ==
2277 : ! == DIFFERENT IRREDUCIBLE PARTS OF THE B.Z. - BUT THEY ==
2278 : ! == DO REPRESENT AN IRREDUCIBLE SET FOR INTEGRATION ==
2279 : ! == OVER THE ENTIRE B.Z. ==
2280 : ! == LWGHT ... THE LIST OF WEIGHTS OF THE CORRESPONDING POINTS. ==
2281 : ! == THESE WEIGHTS ARE NOT NORMALIZED (JUST INTEGERS) ==
2282 : ! == LROT ... FOR EACH SPECIAL POINT THE 'UNFOLDING ROTATIONS' ==
2283 : ! == ARE LISTED. IF E.G. THE WEIGHT OF THE I-TH SPECIAL ==
2284 : ! == POINT IS LWGHT(I), THEN THE ROTATIONS WITH NUMBERS ==
2285 : ! == LROT(J,I), J=1,2,...,LWGHT(I) WILL 'SPREAD' THIS ==
2286 : ! == SINGLE POINT FROM THE IRREDUCIBLE PART OF B.Z. INTO ==
2287 : ! == SEVERAL POINTS IN AN ELEMENTARY UNIT CELL ==
2288 : ! == (PARALLELOPIPED) OF THE RECIPROCAL SPACE. ==
2289 : ! == SOME OPERATION NUMBERS IN THE LIST LROT MAY APPEAR ==
2290 : ! == NEGATIVE, THIS MEANS THAT THE CORRESPONDING ROTATION==
2291 : ! == HAS TO BE APPLIED WITH INVERSION (THE LATTER HAVING ==
2292 : ! == BEEN ARTIFICIALLY ADDED AS SYMMETRY OPERATION IN ==
2293 : ! == CASE INV/=0).NO OTHER EFFORT WAS TAKEN,TO RENUMBER==
2294 : ! == THE ROTATIONS WITH MINUS SIGN OR TO EXTEND THE ==
2295 : ! == LIST OF THE POINT-GROUP OPERATIONS IN THE LIST NB. ==
2296 : ! == INCLUD ... INTEGER ARRAY USED BY SPPT2 INCLUD(NKPOINT) ==
2297 : ! == THE FIRST BIT (0) IS USED BY THE ROUTINE. ==
2298 : ! == THE OTHER BITS GIVE THE K-POINT INDEX IN ==
2299 : ! == THE SPECIAL K-POINT TABLE. ==
2300 : ! ==--------------------------------------------------------------==
2301 : ! == NHASH USED BY MESH ROUTINE ==
2302 : ! == LIST INTEGER ARRAY USED BY MESH LIST(NHASH+NKPOINT) ==
2303 : ! == RLIST real(8) :: ARRAY USED BY MESH RLIST(3,NKPOINT) ==
2304 : ! ==--------------------------------------------------------------==
2305 : ! == Use bit manipulations functions ==
2306 : ! == IBSET(I,POS) sets the bit POS to 1 in I integer ==
2307 : ! == IBCLR(I,POS) clears the bit POS to 1 in I integer ==
2308 : ! == BTEST(I,POS) .TRUE. if bit POS is 1 in I integer ==
2309 : ! ==--------------------------------------------------------------==
2310 : INTEGER :: iout, iq1, iq2, iq3
2311 : REAL(dp) :: wvk0(3)
2312 : INTEGER :: nkpoint
2313 : REAL(dp) :: a1(3), a2(3), a3(3), b1(3), b2(3), b3(3)
2314 : INTEGER :: inv, nc, ib(48)
2315 : REAL(dp) :: r(3, 3, 48)
2316 : INTEGER :: ntot
2317 : REAL(dp) :: wvkl(3, nkpoint)
2318 : INTEGER :: lwght(nkpoint), lrot(48, nkpoint), &
2319 : ncbrav, ibrav(48), istriz, nhash, &
2320 : includ(nkpoint), list(nkpoint + nhash)
2321 : REAL(dp) :: rlist(3, nkpoint), delta
2322 :
2323 : INTEGER, PARAMETER :: no = 0, nrsdir = 100
2324 :
2325 : INTEGER :: i, i1, i2, i3, ibsign, igarb0, igarbage, &
2326 : igarbg, ii, imesh, iop, iplace, &
2327 : iremov, iwvk, j, jplace, k, n, nplane
2328 : REAL(dp) :: diff, proja(3), projb(3), &
2329 : rsdir(4, nrsdir), ur1, ur2, ur3, &
2330 : wva(3), wvk(3)
2331 :
2332 : ! ==--------------------------------------------------------------==
2333 :
2334 960 : ntot = 0
2335 135820 : DO i = 1, nkpoint
2336 134860 : lrot(1, i) = 1
2337 6474240 : DO j = 2, 48
2338 6473280 : lrot(j, i) = 0
2339 : END DO
2340 : END DO
2341 135820 : DO i = 1, nkpoint
2342 135820 : includ(i) = no
2343 : END DO
2344 3840 : DO i = 1, 3
2345 3840 : wva(i) = 0._dp
2346 : END DO
2347 : ! ==--------------------------------------------------------------==
2348 : ! == DEFINE THE 1ST BRILLOUIN ZONE ==
2349 : ! ==--------------------------------------------------------------==
2350 960 : CALL bzdefine(iout, b1, b2, b3, rsdir, nplane, delta)
2351 : ! ==--------------------------------------------------------------==
2352 : ! == Generation of the mesh (they are not multiplied by 2*pi) by ==
2353 : ! == the Monkhorst/Pack algorithm, supplemented by all rotations ==
2354 : ! ==--------------------------------------------------------------==
2355 : ! Initialize the list of vectors
2356 960 : iplace = -2
2357 : CALL mesh(iout, wva, iplace, igarb0, igarbg, nkpoint, nhash, &
2358 960 : list, rlist, delta)
2359 960 : imesh = 0
2360 3250 : DO i1 = 1, iq1
2361 8180 : DO i2 = 1, iq2
2362 20706 : DO i3 = 1, iq3
2363 13486 : ur1 = REAL(1 + iq1 - 2*i1, kind=dp)/REAL(2*iq1, kind=dp)
2364 13486 : ur2 = REAL(1 + iq2 - 2*i2, kind=dp)/REAL(2*iq2, kind=dp)
2365 13486 : ur3 = REAL(1 + iq3 - 2*i3, kind=dp)/REAL(2*iq3, kind=dp)
2366 53944 : DO i = 1, 3
2367 53944 : wvk(i) = ur1*b1(i) + ur2*b2(i) + ur3*b3(i) + wvk0(i)
2368 : END DO
2369 : ! Reduce WVK to the 1st Brillouin zone
2370 : CALL bzrduc(wvk, a1, a2, a3, b1, b2, b3, rsdir, &
2371 13486 : nrsdir, nplane, delta)
2372 18416 : IF (istriz == 1) THEN
2373 : ! Symmetrization of the k-points mesh.
2374 : ! Apply all the Bravais lattice operations to WVK
2375 521454 : DO iop = 1, ncbrav
2376 2031872 : DO i = 1, 3
2377 1523904 : wva(i) = 0._dp
2378 6603584 : DO j = 1, 3
2379 6095616 : wva(i) = wva(i) + r(i, j, ibrav(iop))*wvk(j)
2380 : END DO
2381 : END DO
2382 : ! Check that WVA is inside the 1 Bz.
2383 507968 : IF (.NOT. inside_bz(wva, rsdir, nplane, delta)) THEN
2384 0 : IF (iout > 0) THEN
2385 0 : WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
2386 : END IF
2387 0 : IF (iout > 0) THEN
2388 : WRITE (iout, '(A,3F10.4,/,A,3F10.4,A,/,A,I3,A)') &
2389 0 : ' THE VECTOR ', wva, &
2390 0 : ' GENERATED FROM ', wvk, ' IN THE BASIC MESH', &
2391 0 : ' BY ROTATION NO. ', ibrav(iop), ' IS OUTSIDE THE 1BZ'
2392 : END IF
2393 0 : CPABORT('SPPT2: VECTOR OUTSIDE THE 1BZ')
2394 : END IF
2395 : ! Place WVA in list
2396 507968 : iplace = 0
2397 : CALL mesh(iout, wva, iplace, igarb0, igarbg, &
2398 507968 : nkpoint, nhash, list, rlist, delta)
2399 : ! If WVA was new (and therefore inserted),
2400 : ! IPLACE is the number.
2401 507968 : IF (iplace > 0) imesh = iplace
2402 521454 : IF (iplace > nkpoint) THEN
2403 0 : IF (iout > 0) THEN
2404 0 : WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
2405 : END IF
2406 0 : IF (iout > 0) THEN
2407 0 : WRITE (iout, *) 'MESH SIZE EXCEEDS NKPOINT=', nkpoint
2408 : END IF
2409 0 : CPABORT('SPPT2: MESH SIZE EXCEEDED')
2410 : END IF
2411 : END DO
2412 : ELSE
2413 : ! Place WVK in list
2414 0 : iplace = 0
2415 : CALL mesh(iout, wvk, iplace, igarb0, igarbg, &
2416 0 : nkpoint, nhash, list, rlist, delta)
2417 0 : imesh = iplace
2418 0 : IF (iplace > nkpoint) THEN
2419 0 : IF (iout > 0) THEN
2420 0 : WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
2421 : END IF
2422 0 : IF (iout > 0) THEN
2423 0 : WRITE (iout, *) 'MESH SIZE EXCEEDS NKPOINT=', nkpoint
2424 : END IF
2425 0 : CPABORT('SPPT2: MESH SIZE EXCEEDED')
2426 : END IF
2427 : END IF
2428 : END DO
2429 : END DO
2430 : END DO
2431 : !deb
2432 : !deb get full mesh
2433 : !deb
2434 960 : IF (iout > 0) THEN
2435 : ! IMESH: Number of k points in the mesh.
2436 : WRITE (iout, &
2437 335 : '(" KPSYM| THE WAVEVECTOR MESH CONTAINS ",I5," POINTS")') imesh
2438 335 : WRITE (iout, '(" KPSYM| THE POINTS ARE:")')
2439 4903 : DO ii = 1, imesh
2440 4568 : i = ii
2441 : CALL mesh(iout, wva, i, igarb0, igarbg, nkpoint, nhash, &
2442 4568 : list, rlist, delta)
2443 4903 : IF (MOD(i, 2) == 1) THEN
2444 2306 : WRITE (iout, '(1X,I5,3F10.4)', advance="no") i, wva
2445 : ELSE
2446 2262 : WRITE (iout, '(1X,I5,3F10.4)') i, wva
2447 : END IF
2448 : END DO
2449 335 : WRITE (iout, *)
2450 : END IF
2451 : ! ==--------------------------------------------------------------==
2452 960 : IF (istriz == 1) THEN
2453 : ! Now figure out if any special point difference (K - K'') is an
2454 : ! integral multiple of a reciprocal-space vector
2455 960 : iremov = 0
2456 14458 : DO i = 1, (imesh - 1)
2457 13498 : iplace = i
2458 : CALL mesh(iout, wva, iplace, igarb0, igarbg, &
2459 13498 : nkpoint, nhash, list, rlist, delta)
2460 : ! Project WVA onto B1,2,3:
2461 13498 : proja(1) = 0._dp
2462 13498 : proja(2) = 0._dp
2463 13498 : proja(3) = 0._dp
2464 53992 : DO k = 1, 3
2465 40494 : proja(1) = proja(1) + wva(k)*a1(k)
2466 40494 : proja(2) = proja(2) + wva(k)*a2(k)
2467 53992 : proja(3) = proja(3) + wva(k)*a3(k)
2468 : END DO
2469 : ! Now loop over all the rest of the mesh points
2470 359258 : loop_mesh: DO j = (i + 1), imesh
2471 344800 : jplace = j
2472 : CALL mesh(iout, wvk, jplace, igarb0, igarbg, &
2473 344800 : nkpoint, nhash, list, rlist, delta)
2474 : ! Project WVK onto B1,2,3:
2475 344800 : projb(1) = 0._dp
2476 344800 : projb(2) = 0._dp
2477 344800 : projb(3) = 0._dp
2478 1379200 : DO k = 1, 3
2479 1034400 : projb(1) = projb(1) + wvk(k)*a1(k)
2480 1034400 : projb(2) = projb(2) + wvk(k)*a2(k)
2481 1379200 : projb(3) = projb(3) + wvk(k)*a3(k)
2482 : END DO
2483 : ! Check (PROJA - PROJB): Is it integral ?
2484 438190 : DO k = 1, 3
2485 438142 : diff = proja(k) - projb(k)
2486 438190 : IF (ABS(REAL(NINT(diff), kind=dp) - diff) > delta) CYCLE loop_mesh
2487 : END DO
2488 : ! DIFF is integral: remove WVK from mesh:
2489 : CALL remove(wvk, jplace, igarb0, igarbg, &
2490 48 : nkpoint, nhash, list, rlist, delta)
2491 : ! If WVK actually removed, increment IREMOV
2492 13546 : IF (jplace > 0) iremov = iremov + 1
2493 : END DO loop_mesh
2494 : END DO
2495 960 : IF (iremov > 0 .AND. iout > 0) THEN
2496 : WRITE (iout, '(A,A,/,A,1X,I6,A,/)') &
2497 1 : ' KPSYM| SOME OF THESE MESH POINTS ARE RELATED BY LATTICE ', &
2498 1 : 'TRANSLATION VECTORS', &
2499 2 : ' KPSYM|', iremov, ' OF THE MESH POINTS REMOVED.'
2500 : END IF
2501 : END IF
2502 : ! ==--------------------------------------------------------------==
2503 : ! == IN THE MESH OF WAVEVECTORS, NOW SEARCH FOR EQUIVALENT POINTS:==
2504 : ! == THE INVERSION (TIME REVERSAL !) MAY BE USED. ==
2505 : ! ==--------------------------------------------------------------==
2506 15418 : DO iwvk = 1, imesh
2507 : ! IF(INCLUD(IWVK) == YES) CYCLE
2508 14458 : IF (BTEST(includ(iwvk), 0)) CYCLE
2509 : ! IWVK has not been encountered previously: new special point,
2510 : ! (only if WVK is not a garbage vector, however.)
2511 : ! INCLUD(IWVK) = YES
2512 2558 : includ(iwvk) = IBSET(includ(iwvk), 0)
2513 2558 : iplace = iwvk
2514 : CALL mesh(iout, wvk, iplace, igarb0, igarbg, &
2515 2558 : nkpoint, nhash, list, rlist, delta)
2516 : ! Find out whether Wvk is in the garbage list
2517 : CALL garbag(wvk, igarbage, igarb0, &
2518 2558 : nkpoint, nhash, list, rlist, delta)
2519 2558 : IF (igarbage > 0) CYCLE
2520 2510 : ntot = ntot + 1
2521 : ! Give the index in the special k points table.
2522 2510 : includ(iwvk) = includ(iwvk) + ntot*2
2523 10040 : DO i = 1, 3
2524 10040 : wvkl(i, ntot) = wvk(i)
2525 : END DO
2526 2510 : lwght(ntot) = 1
2527 : ! ==-----------------------------------------------------------==
2528 : ! Find all the equivalent points (symmetry given by atoms)
2529 46838 : equivalent_points: DO n = 1, nc
2530 : ! Rotate:
2531 173472 : DO i = 1, 3
2532 130104 : wva(i) = 0._dp
2533 563784 : DO j = 1, 3
2534 520416 : wva(i) = wva(i) + r(i, j, ib(n))*wvk(j)
2535 : END DO
2536 : END DO
2537 : ibsign = +1
2538 14458 : DO
2539 : ! Find WVA in the list
2540 47108 : iplace = -1
2541 : CALL mesh(iout, wva, iplace, igarb0, igarbg, &
2542 47108 : nkpoint, nhash, list, rlist, delta)
2543 47108 : IF (iplace == 0) THEN
2544 48 : IF (istriz /= -1) THEN
2545 : ! Find out whether WVA is in the garbage list
2546 : CALL garbag(wva, igarbage, igarb0, &
2547 48 : nkpoint, nhash, list, rlist, delta)
2548 48 : IF (igarbage == 0) THEN
2549 : ! I think this case is impossible (NC <= NCBRAV)
2550 : ! Error message
2551 0 : IF (iout > 0) THEN
2552 0 : WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
2553 : END IF
2554 0 : IF (iout > 0) THEN
2555 : WRITE (iout, '(A,3F10.4,/,A,3F10.4,A,/,A,I3,A)') &
2556 0 : ' THE VECTOR ', wva, &
2557 0 : ' GENERATED FROM ', wvk, ' IN THE BASIC MESH', &
2558 0 : ' BY ROTATION NO. ', ib(n), ' IS NOT IN THE LIST'
2559 : END IF
2560 0 : CPABORT('SPPT2: VECTOR NOT IN THE LIST')
2561 : END IF
2562 : END IF
2563 : END IF
2564 47108 : IF (iplace /= 0 .OR. istriz /= -1) THEN
2565 : ! Find out whether WVA is in the garbage list
2566 : CALL garbag(wva, igarbage, igarb0, &
2567 47108 : nkpoint, nhash, list, rlist, delta)
2568 47108 : IF (igarbage > 0) CYCLE equivalent_points
2569 : ! Was WVA encountered before ?
2570 47060 : IF (.NOT. BTEST(includ(iplace), 0)) THEN
2571 : ! Increment weight.
2572 11900 : lwght(ntot) = lwght(ntot) + 1
2573 11900 : lrot(lwght(ntot), ntot) = ib(n)*ibsign
2574 : ! INCLUD(IPLACE) = YES
2575 11900 : includ(iplace) = IBSET(includ(iplace), 0)
2576 : ! This k-point is an image of a special k-point.
2577 : ! Put the index of the special k-point.
2578 11900 : includ(iplace) = includ(iplace) + ntot*2
2579 : END IF
2580 : END IF
2581 47060 : IF (ibsign == -1 .OR. inv == 0) CYCLE equivalent_points
2582 : ! The case where we also apply the inversion to WVA
2583 : ! Repeat the search, but for -WVA
2584 3740 : ibsign = -1
2585 58328 : DO i = 1, 3
2586 14960 : wva(i) = -wva(i)
2587 : END DO
2588 : END DO
2589 : END DO equivalent_points
2590 : END DO
2591 : ! ==--------------------------------------------------------------==
2592 : ! == TOTAL NUMBER OF SPECIAL POINTS: NTOT ==
2593 : ! == BEFORE USING THE LIST WVKL AS WAVE VECTORS, THEY HAVE TO BE ==
2594 : ! == MULTIPLIED BY 2*PI ==
2595 : ! == THE LIST OF WEIGHTS LWGHT IS NOT NORMALIZED ==
2596 : ! ==--------------------------------------------------------------==
2597 960 : IF (ntot > nkpoint) THEN
2598 0 : IF (iout > 0) THEN
2599 0 : WRITE (iout, *) 'IN SPPT2 NUMBER OF SPECIAL POINTS = ', ntot
2600 : END IF
2601 0 : IF (iout > 0) THEN
2602 0 : WRITE (iout, *) 'BUT NKPOINT = ', nkpoint
2603 : END IF
2604 0 : ntot = -1
2605 : END IF
2606 960 : IF (iout > 0) THEN
2607 : ! Write the index table relating k points in the mesh
2608 : ! with special k points
2609 : IF (iout > 0) THEN
2610 : WRITE (iout, '(/,A,4X,A)') &
2611 335 : ' KPSYM|', 'CROSS TABLE RELATING MESH POINTS WITH SPECIAL POINTS:'
2612 : END IF
2613 335 : IF (iout > 0) THEN
2614 335 : WRITE (iout, '(5(4X,"IK -> SK"))')
2615 : END IF
2616 4903 : DO i = 1, imesh
2617 4568 : iplace = includ(i)/2
2618 4568 : IF (iout > 0) THEN
2619 4568 : WRITE (iout, '(1X,I5,1X,I5)', advance="no") i, iplace
2620 : END IF
2621 4903 : IF ((MOD(i, 5) == 0) .AND. iout > 0) THEN
2622 721 : WRITE (iout, *)
2623 : END IF
2624 : END DO
2625 335 : IF ((MOD(j - 1, 5) /= 0) .AND. iout > 0) THEN
2626 335 : WRITE (iout, *)
2627 : END IF
2628 : END IF
2629 960 : END SUBROUTINE sppt2
2630 : ! **************************************************************************************************
2631 : !> \brief ...
2632 : !> \param iout ...
2633 : !> \param wvk ...
2634 : !> \param iplace ...
2635 : !> \param igarb0 ...
2636 : !> \param igarbg ...
2637 : !> \param nmesh ...
2638 : !> \param nhash ...
2639 : !> \param list ...
2640 : !> \param rlist ...
2641 : !> \param delta ...
2642 : ! **************************************************************************************************
2643 921460 : SUBROUTINE mesh(iout, wvk, iplace, igarb0, igarbg, &
2644 921460 : nmesh, nhash, list, rlist, delta)
2645 : ! ==--------------------------------------------------------------==
2646 : ! == MESH MAINTAINS A LIST OF VECTORS FOR PLACEMENT AND/OR LOOKUP ==
2647 : ! == ==
2648 : ! == ADDITIONAL ENTRY POINTS: REMOVE .... REMOVE VECTOR FROM LIST ==
2649 : ! == GARBAG .... WAS VECTOR REMOVED ? ==
2650 : ! == ==
2651 : ! == WVK ....... VECTOR ==
2652 : ! == IPLACE .... ON INPUT: -2 MEANS: INITIALIZE THE LIST ==
2653 : ! == (AND RETURN) ==
2654 : ! == -1 MEANS: FIND WVK IN THE LIST ==
2655 : ! == 0 MEANS: ADD WVK TO THE LIST ==
2656 : ! == >0 MEANS: RETURN WVK NO. IPLACE ==
2657 : ! == ON OUTPUT: THE POSITION ASSIGNED TO WVK ==
2658 : ! == (=0 IF WVK IS NOT IN THE LIST) ==
2659 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
2660 : ! ==--------------------------------------------------------------==
2661 : INTEGER :: iout
2662 : REAL(dp) :: wvk(3)
2663 : INTEGER :: iplace, igarb0, igarbg, nmesh, nhash, &
2664 : list(nhash + nmesh)
2665 : REAL(dp) :: rlist(3, nmesh), delta
2666 :
2667 : INTEGER, PARAMETER :: nil = 0
2668 :
2669 : INTEGER :: i, ihash, ipoint
2670 : INTEGER, SAVE :: istore
2671 : REAL(dp) :: delta1, rhash
2672 :
2673 : ! ==--------------------------------------------------------------==
2674 : ! == Initialization ==
2675 : ! ==--------------------------------------------------------------==
2676 :
2677 921460 : delta1 = 10._dp*delta
2678 921460 : IF (iplace <= -2) THEN
2679 1100460 : DO i = 1, nhash + nmesh
2680 1100460 : list(i) = nil
2681 : END DO
2682 960 : istore = 1
2683 : ! IGARB0 points to a linked list of removed WVKS (the garbage).
2684 960 : igarb0 = 0
2685 960 : igarbg = 0
2686 960 : RETURN
2687 : ! ==--------------------------------------------------------------==
2688 920500 : ELSE IF ((iplace > -2) .AND. (iplace <= 0)) THEN
2689 : ! The particular HASH function used in this case:
2690 : rhash = 0.7890_dp*wvk(1) &
2691 : + 0.6810_dp*wvk(2) &
2692 555076 : + 0.5811_dp*wvk(3) + delta
2693 555076 : ihash = INT(ABS(rhash)*REAL(nhash, kind=dp))
2694 555076 : ihash = MOD(ihash, nhash) + nmesh + 1
2695 : ! Search for WVK in linked list
2696 555076 : ipoint = list(ihash)
2697 803582 : DO i = 1, 100
2698 : ! List exhausted
2699 803582 : IF (ipoint == nil) EXIT
2700 : ! Compare WVK with this element
2701 2439106 : IF (ALL(ABS(wvk(:) - rlist(:, ipoint)) <= delta1)) THEN
2702 : ! WVK located
2703 540570 : IF (iplace == 0) RETURN
2704 : ! IPLACE=-1
2705 47060 : iplace = ipoint
2706 47060 : RETURN
2707 : END IF
2708 : ! Next element of list
2709 248506 : ihash = ipoint
2710 263012 : ipoint = list(ihash)
2711 : END DO
2712 14506 : IF (ipoint /= nil) THEN
2713 : ! List too long
2714 0 : IF (iout > 0) THEN
2715 : WRITE (iout, '(2A,/,A)') &
2716 0 : ' SUBROUTINE MESH *** FATAL ERROR *** LINKED LIST', &
2717 0 : ' TOO LONG ***', ' CHOOSE A BETTER HASH-FUNCTION'
2718 : END IF
2719 0 : CPABORT('MESH: WARNING')
2720 : END IF
2721 : ! WVK was not found
2722 14506 : IF (iplace == -1) THEN
2723 : ! IPLACE=-1 : search for WVK unsuccessful
2724 48 : iplace = 0
2725 48 : RETURN
2726 : ELSE
2727 : ! IPLACE=0: add WVK to the list
2728 14458 : list(ihash) = istore
2729 14458 : IF (istore > nmesh) THEN
2730 0 : IF (iout > 0) THEN
2731 0 : WRITE (iout, '(A)') 'SUBROUTINE MESH *** FATAL ERROR ***'
2732 : END IF
2733 0 : IF (iout > 0) THEN
2734 : WRITE (iout, '(A,I10,A,/,A,3F10.5)') &
2735 0 : ' ISTORE=', istore, ' EXCEEDS DIMENSIONS', &
2736 0 : ' WVK = ', wvk
2737 : END IF
2738 0 : CPABORT('MESH: WARNING')
2739 : END IF
2740 14458 : list(istore) = nil
2741 57832 : DO i = 1, 3
2742 57832 : rlist(i, istore) = wvk(i)
2743 : END DO
2744 14458 : istore = istore + 1
2745 14458 : iplace = istore - 1
2746 14458 : RETURN
2747 : END IF
2748 : ! WVK was found
2749 : ELSE
2750 : ! ==--------------------------------------------------------------==
2751 : ! == Return a wavevector (IPLACE > 0) ==
2752 : ! ==--------------------------------------------------------------==
2753 365424 : ipoint = iplace
2754 365424 : IF (ipoint >= istore) THEN
2755 0 : IF (iout > 0) THEN
2756 : WRITE (iout, '(A,/,A,I5,A,/)') &
2757 0 : ' SUBROUTINE MESH *** WARNING ***', &
2758 0 : ' IPLACE = ', iplace, &
2759 0 : ' IS BEYOND THE LISTS - WVK SET TO 1.0E38'
2760 : END IF
2761 0 : DO i = 1, 3
2762 0 : wvk(i) = 1.0e38_dp
2763 : END DO
2764 : END IF
2765 1461696 : DO i = 1, 3
2766 1461696 : wvk(i) = rlist(i, ipoint)
2767 : END DO
2768 : END IF
2769 : END SUBROUTINE mesh
2770 : ! **************************************************************************************************
2771 : !> \brief ...
2772 : !> \param wvk ...
2773 : !> \param iplace ...
2774 : !> \param igarb0 ...
2775 : !> \param igarbg ...
2776 : !> \param nmesh ...
2777 : !> \param nhash ...
2778 : !> \param list ...
2779 : !> \param rlist ...
2780 : !> \param delta ...
2781 : ! **************************************************************************************************
2782 48 : SUBROUTINE remove(wvk, iplace, igarb0, igarbg, &
2783 48 : nmesh, nhash, list, rlist, delta)
2784 : ! ==--------------------------------------------------------------==
2785 : ! == ENTRY POINT FOR REMOVING A WAVEVECTOR ==
2786 : ! == ==
2787 : ! == INPUT: ==
2788 : ! == WVK(3) ==
2789 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
2790 : ! == OUTPUT:
2791 : ! == IPLACE .....1 IF WVK WAS REMOVED ==
2792 : ! == 0 IF WVK WAS NOT REMOVED ==
2793 : ! == (WVK NOT IN THE LINKED LISTS) ==
2794 : ! ==--------------------------------------------------------------==
2795 : REAL(dp) :: wvk(3)
2796 : INTEGER :: iplace, igarb0, igarbg, nmesh, nhash, &
2797 : list(nhash + nmesh)
2798 : REAL(dp) :: rlist(3, nmesh), delta
2799 :
2800 : INTEGER, PARAMETER :: nil = 0
2801 :
2802 : INTEGER :: i, ihash, ipoint
2803 : REAL(dp) :: delta1, rhash
2804 :
2805 : ! ==--------------------------------------------------------------==
2806 : ! Variables
2807 : ! ==--------------------------------------------------------------==
2808 :
2809 48 : delta1 = 10._dp*delta
2810 : ! The particular hash function used in this case:
2811 : rhash = 0.7890_dp*wvk(1) &
2812 : + 0.6810_dp*wvk(2) &
2813 48 : + 0.5811_dp*wvk(3) + delta
2814 48 : ihash = INT(ABS(rhash)*REAL(nhash, kind=dp))
2815 48 : ihash = MOD(ihash, nhash) + nmesh + 1
2816 : ! Search for WVK in linked list
2817 48 : ipoint = list(ihash)
2818 96 : DO i = 1, 100
2819 : ! List exhausted
2820 96 : IF (ipoint == nil) THEN
2821 : ! WVK was not found in the mesh:
2822 0 : iplace = 0
2823 0 : RETURN
2824 : END IF
2825 : ! Compare WVK with this element
2826 240 : IF (.NOT. ANY(ABS(wvk(:) - rlist(:, ipoint)) > delta1)) THEN
2827 : ! WVK located, now remove it from the list:
2828 48 : list(ihash) = list(ipoint)
2829 : ! LIST(IHASH) now points to the next element in the list,
2830 : ! and the present WVK has become garbage.
2831 : ! Add WVK to the list of garbage:
2832 48 : IF (igarb0 == 0) THEN
2833 : ! Start up the garbage list:
2834 4 : igarb0 = ipoint
2835 : ELSE
2836 44 : list(igarbg) = ipoint
2837 : END IF
2838 48 : igarbg = ipoint
2839 48 : list(igarbg) = nil
2840 48 : iplace = 1
2841 48 : RETURN
2842 : END IF
2843 : ! Next element of list
2844 48 : ihash = ipoint
2845 48 : ipoint = list(ihash)
2846 : END DO
2847 : ! List too long
2848 0 : CPABORT('MESH: LIST TOO LONG')
2849 : END SUBROUTINE remove
2850 : ! **************************************************************************************************
2851 : !> \brief ...
2852 : !> \param wvk ...
2853 : !> \param iplace ...
2854 : !> \param igarb0 ...
2855 : !> \param nmesh ...
2856 : !> \param nhash ...
2857 : !> \param list ...
2858 : !> \param rlist ...
2859 : !> \param delta ...
2860 : ! **************************************************************************************************
2861 49714 : SUBROUTINE garbag(wvk, iplace, igarb0, &
2862 49714 : nmesh, nhash, list, rlist, delta)
2863 : ! ==--------------------------------------------------------------==
2864 : ! == ENTRY POINT FOR CHECKING IF A WAVEVECTOR ==
2865 : ! == IS IN THE GARBAGE LIST ==
2866 : ! == INPUT: ==
2867 : ! == WVK(3) ==
2868 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
2869 : ! == ==
2870 : ! == OUTPUT: ==
2871 : ! == IPLACE ..... I > 0 IS THE PLACE IN THE GARBAGE LIST ==
2872 : ! == 0 IF WVK NOT AMONG THE GARBAGE ==
2873 : ! ==--------------------------------------------------------------==
2874 : REAL(dp) :: wvk(3)
2875 : INTEGER :: iplace, igarb0, nmesh, nhash, &
2876 : list(nhash + nmesh)
2877 : REAL(dp) :: rlist(3, nmesh), delta
2878 :
2879 : INTEGER, PARAMETER :: nil = 0
2880 :
2881 : INTEGER :: i, ihash, ipoint
2882 : REAL(dp) :: delta1
2883 :
2884 : ! ==--------------------------------------------------------------==
2885 : ! Variables
2886 : ! ==--------------------------------------------------------------==
2887 :
2888 49714 : delta1 = 10._dp*delta
2889 : ! Search for WVK in linked list
2890 : ! Point to the garbage list
2891 49714 : ipoint = igarb0
2892 55978 : DO i = 1, nmesh
2893 : ! LIST EXHAUSTED
2894 55978 : IF (ipoint == nil) THEN
2895 : ! WVK was not found in the mesh:
2896 49570 : iplace = 0
2897 49570 : RETURN
2898 : END IF
2899 : ! Compare WVK with this element
2900 7432 : IF (.NOT. ANY(ABS(wvk(:) - rlist(:, ipoint)) > delta1)) THEN
2901 : ! WVK was located in the garbage list
2902 144 : iplace = i
2903 144 : RETURN
2904 : END IF
2905 : ! Next element of list
2906 6264 : ihash = ipoint
2907 6264 : ipoint = list(ihash)
2908 : END DO
2909 : ! List too long
2910 0 : CPABORT('GARBAG: LIST TOO LONG')
2911 : END SUBROUTINE garbag
2912 :
2913 : ! **************************************************************************************************
2914 : !> \brief ...
2915 : !> \param wvk ...
2916 : !> \param a1 ...
2917 : !> \param a2 ...
2918 : !> \param a3 ...
2919 : !> \param b1 ...
2920 : !> \param b2 ...
2921 : !> \param b3 ...
2922 : !> \param rsdir ...
2923 : !> \param nrsdir ...
2924 : !> \param nplane ...
2925 : !> \param delta ...
2926 : ! **************************************************************************************************
2927 13486 : SUBROUTINE bzrduc(wvk, a1, a2, a3, b1, b2, b3, rsdir, nrsdir, nplane, delta)
2928 : ! ==--------------------------------------------------------------==
2929 : ! == REDUCE WVK TO LIE ENTIRELY WITHIN THE 1ST BRILLOUIN ZONE ==
2930 : ! == BY ADDING B-VECTORS ==
2931 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
2932 : ! ==--------------------------------------------------------------==
2933 : REAL(dp) :: wvk(3), a1(3), a2(3), a3(3), b1(3), &
2934 : b2(3), b3(3)
2935 : INTEGER :: nrsdir
2936 : REAL(dp) :: rsdir(4, nrsdir)
2937 : INTEGER :: nplane
2938 : REAL(dp) :: delta
2939 :
2940 : INTEGER, PARAMETER :: nzones = 4, nnn = 2*nzones + 1, &
2941 : nn = nzones + 1
2942 :
2943 : INTEGER :: i, i1, i2, i3, n1, n2, n3, nn1, nn2, nn3
2944 : LOGICAL :: inside
2945 : REAL(dp) :: wb(3), wva(3)
2946 :
2947 : ! ==--------------------------------------------------------------==
2948 : ! Variables
2949 : ! Look around +/- "NZONES" to locate vector
2950 : ! NZONES may need to be increased for very anisotropic zones
2951 : ! ==--------------------------------------------------------------==
2952 :
2953 13486 : IF (.NOT. inside_bz(wvk, rsdir, nplane, delta)) THEN
2954 40 : inside = .FALSE.
2955 : ! Express WVK in the basis of B1,2,3.
2956 : ! This permits an estimate of how far WVK is from the 1Bz.
2957 40 : wb(1) = wvk(1)*a1(1) + wvk(2)*a1(2) + wvk(3)*a1(3)
2958 40 : wb(2) = wvk(1)*a2(1) + wvk(2)*a2(2) + wvk(3)*a2(3)
2959 40 : wb(3) = wvk(1)*a3(1) + wvk(2)*a3(2) + wvk(3)*a3(3)
2960 40 : nn1 = NINT(wb(1))
2961 40 : nn2 = NINT(wb(2))
2962 40 : nn3 = NINT(wb(3))
2963 : ! Look around the estimated vector for the one truly inside the 1Bz
2964 192 : n1_loop: DO n1 = 1, nnn
2965 192 : i1 = nn - n1 - nn1
2966 1768 : DO n2 = 1, nnn
2967 1576 : i2 = nn - n2 - nn2
2968 15712 : DO n3 = 1, nnn
2969 14024 : i3 = nn - n3 - nn3
2970 56096 : DO i = 1, 3
2971 : wva(i) = wvk(i) + REAL(i1, kind=dp)*b1(i) + REAL(i2, kind=dp)*b2(i) + &
2972 56096 : REAL(i3, kind=dp)*b3(i)
2973 : END DO
2974 14024 : inside = inside_bz(wva, rsdir, nplane, delta)
2975 15560 : IF (inside) EXIT n1_loop
2976 : END DO
2977 : END DO
2978 : END DO n1_loop
2979 40 : CPASSERT(inside)
2980 40 : wvk(1:3) = wva(1:3)
2981 : END IF
2982 :
2983 13486 : END SUBROUTINE bzrduc
2984 :
2985 : ! **************************************************************************************************
2986 : !> \brief Is wvk in the 1st Brillouin zone ?
2987 : !> Check whether wvk lies inside all the planes that define the 1bz.
2988 : !> \param wvk ...
2989 : !> \param rsdir ...
2990 : !> \param nplane ...
2991 : !> \param delta ...
2992 : !> \return ...
2993 : ! **************************************************************************************************
2994 535478 : FUNCTION inside_bz(wvk, rsdir, nplane, delta) RESULT(inbz)
2995 : REAL(KIND=dp), DIMENSION(3) :: wvk
2996 : REAL(KIND=dp), DIMENSION(:, :) :: rsdir
2997 : INTEGER :: nplane
2998 : REAL(KIND=dp) :: delta
2999 : LOGICAL :: inbz
3000 :
3001 : INTEGER :: n
3002 : REAL(KIND=dp) :: projct
3003 :
3004 535478 : inbz = .TRUE.
3005 2125652 : DO n = 1, nplane
3006 1604198 : projct = (rsdir(1, n)*wvk(1) + rsdir(2, n)*wvk(2) + rsdir(3, n)*wvk(3))/rsdir(4, n)
3007 2125652 : IF (ABS(projct) > 0.5_dp + delta) THEN
3008 : inbz = .FALSE.
3009 : EXIT
3010 : END IF
3011 : END DO
3012 :
3013 535478 : END FUNCTION inside_bz
3014 :
3015 : ! **************************************************************************************************
3016 : !> \brief Find the vectors whose halves define the 1st Brillouin zone
3017 : !> Output:
3018 : !> nplane -- How many elements of rsdir contain normal vectors defining the planes
3019 : !> Method:
3020 : !> Starting with the parallelopiped spanned by b1,2,3 around the origin,
3021 : !> vectors inside a sufficiently large sphere are tested to see whether
3022 : !> the planes at 1/2*b will further confine the 1bz.
3023 : !> The resulting vectors are not cleaned to avoid redundant planes
3024 : !> \param iout ...
3025 : !> \param b1 ...
3026 : !> \param b2 ...
3027 : !> \param b3 ...
3028 : !> \param rsdir ...
3029 : !> \param nplane ...
3030 : !> \param delta ...
3031 : ! **************************************************************************************************
3032 960 : SUBROUTINE bzdefine(iout, b1, b2, b3, rsdir, nplane, delta)
3033 : INTEGER :: iout
3034 : REAL(KIND=dp), DIMENSION(3) :: b1, b2, b3
3035 : REAL(KIND=dp), DIMENSION(:, :) :: rsdir
3036 : INTEGER :: nplane
3037 : REAL(KIND=dp) :: delta
3038 :
3039 : INTEGER :: i, i1, i2, i3, n, n1, n2, n3, nb1, nb2, &
3040 : nb3, nnb1, nnb2, nnb3, nrsdir
3041 : REAL(KIND=dp) :: b1len, b2len, b3len, bmax, projct
3042 : REAL(KIND=dp), DIMENSION(3) :: bvec
3043 :
3044 960 : nrsdir = SIZE(rsdir, 2)
3045 :
3046 960 : b1len = b1(1)**2 + b1(2)**2 + b1(3)**2
3047 960 : b2len = b2(1)**2 + b2(2)**2 + b2(3)**2
3048 960 : b3len = b3(1)**2 + b3(2)**2 + b3(3)**2
3049 : ! Lattice containing entirely the Brillouin zone
3050 960 : bmax = b1len + b2len + b3len
3051 960 : nb1 = INT(SQRT(bmax/b1len) + delta) + 1
3052 960 : nb2 = INT(SQRT(bmax/b2len) + delta) + 1
3053 960 : nb3 = INT(SQRT(bmax/b3len) + delta) + 1
3054 480960 : rsdir(:, :) = 0._dp
3055 : ! 1Bz is certainly confined inside the 1/2(B1,B2,B3) parallelopiped
3056 3840 : rsdir(1:3, 1) = b1(1:3)
3057 3840 : rsdir(1:3, 2) = b2(1:3)
3058 3840 : rsdir(1:3, 3) = b3(1:3)
3059 960 : rsdir(4, 1) = b1len
3060 960 : rsdir(4, 2) = b2len
3061 960 : rsdir(4, 3) = b3len
3062 : ! Starting confinement: 3 planes
3063 960 : nplane = 3
3064 960 : nnb1 = 2*nb1 + 1
3065 960 : nnb2 = 2*nb2 + 1
3066 960 : nnb3 = 2*nb3 + 1
3067 :
3068 5760 : DO n1 = 1, nnb1
3069 4800 : i1 = nb1 + 1 - n1
3070 37080 : DO n2 = 1, nnb2
3071 31320 : i2 = nb2 + 1 - n2
3072 222780 : inner_loop: DO n3 = 1, nnb3
3073 186660 : i3 = nb3 + 1 - n3
3074 186660 : IF (i1 == 0 .AND. i2 == 0 .AND. i3 == 0) CYCLE inner_loop
3075 742800 : DO i = 1, 3
3076 : bvec(i) = REAL(i1, kind=dp)*b1(i) + REAL(i2, kind=dp)*b2(i) + &
3077 742800 : REAL(i3, kind=dp)*b3(i)
3078 : END DO
3079 : ! Does the plane of 1/2*BVEC narrow down the 1Bz ?
3080 232156 : DO n = 1, nplane
3081 : projct = 0.5_dp*(rsdir(1, n)*bvec(1) + rsdir(2, n)*bvec(2) &
3082 231780 : + rsdir(3, n)*bvec(3))/rsdir(4, n)
3083 : ! 1/2*BVEC is outside the Bz - skip this direction
3084 : ! The 1.e-6_dp takes care of single points touching the Bz,
3085 : ! and of the -(plane)
3086 232156 : IF (ABS(projct) > 0.5_dp - delta) CYCLE inner_loop
3087 : END DO
3088 : ! 1/2*BVEC further confines the 1Bz - include into RSDIR
3089 376 : nplane = nplane + 1
3090 376 : CPASSERT(nplane <= nrsdir)
3091 1504 : DO i = 1, 3
3092 1504 : rsdir(i, nplane) = bvec(i)
3093 : END DO
3094 : ! Length squared
3095 217980 : rsdir(4, nplane) = bvec(1)**2 + bvec(2)**2 + bvec(3)**2
3096 : END DO inner_loop
3097 : END DO
3098 : END DO
3099 :
3100 960 : IF (iout > 0) THEN
3101 : WRITE (iout, '(A,I3,A,/,A,/,100(" KPSYM|",1X,3F10.4,/))') &
3102 335 : ' KPSYM| The 1st Brillouin zone is confined by (at most)', &
3103 335 : nplane, ' planes', &
3104 335 : ' KPSYM| as defined by the +/- halves of the vectors:', &
3105 670 : ((rsdir(i, n), i=1, 3), n=1, nplane)
3106 : END IF
3107 :
3108 960 : END SUBROUTINE bzdefine
3109 :
3110 : END MODULE kpsym
|
