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 728 : SUBROUTINE k290s(iout, nat, nkpoint, nsp, iq1, iq2, iq3, istriz, &
79 728 : a1, a2, a3, alat, strain, xkapa, rx, tvec, &
80 728 : ty, isc, f0, ntvec, wvk0, wvkl, lwght, lrot, &
81 728 : 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 728 : f00 = 0
181 728 : x0 = 0._dp
182 728 : v = 0._dp
183 728 : v0 = 0._dp
184 : ! ==--------------------------------------------------------------==
185 : ! READ IN LATTICE STRUCTURE
186 : ! ==--------------------------------------------------------------==
187 2912 : DO i = 1, 3
188 2184 : a01(i) = a1(i)/alat
189 2184 : a02(i) = a2(i)/alat
190 2912 : a03(i) = a3(i)/alat
191 : END DO
192 728 : IF (iout > 0) THEN
193 292 : WRITE (iout, '(" KPSYM| NUMBER OF ATOMS (STRUCT):",I6)') nat
194 : END IF
195 728 : IF (iout > 0) THEN
196 292 : WRITE (iout, '(" KPSYM|",10X,"K TYPE",14X,"X(K)")')
197 : END IF
198 728 : itype = 0
199 4940 : DO i = 1, nat
200 : ! Assign an atomic type (for internal purposes)
201 4212 : located_type = .FALSE.
202 4212 : IF (i /= 1) THEN
203 4600 : DO j = 1, (i - 1)
204 4600 : 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 3484 : IF (.NOT. located_type) THEN
213 940 : itype = itype + 1
214 940 : 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 4940 : IF (iout > 0) THEN
228 : WRITE (iout, '(" KPSYM|",6X,I5,I6,3F10.5)') &
229 1626 : i, ty(i), (xkapa(j, i), j=1, 3)
230 : END IF
231 : END DO
232 : ! ==--------------------------------------------------------------==
233 : ! IS THE STRAIN SIGNIFICANT ?
234 : ! ==--------------------------------------------------------------==
235 728 : dtotstr = delta*delta
236 728 : totstr = 0._dp
237 728 : istrin = 0
238 5096 : DO i = 1, 6
239 5096 : totstr = totstr + ABS(strain(i))
240 : END DO
241 728 : 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 728 : A2(3)*A3(2)*A1(1) - A3(3)*A1(2)*A2(1)
247 728 : volum = ABS(volum)
248 728 : b1(1) = (a2(2)*a3(3) - a2(3)*a3(2))/volum
249 728 : b1(2) = (a2(3)*a3(1) - a2(1)*a3(3))/volum
250 728 : b1(3) = (a2(1)*a3(2) - a2(2)*a3(1))/volum
251 728 : b2(1) = (a3(2)*a1(3) - a3(3)*a1(2))/volum
252 728 : b2(2) = (a3(3)*a1(1) - a3(1)*a1(3))/volum
253 728 : b2(3) = (a3(1)*a1(2) - a3(2)*a1(1))/volum
254 728 : b3(1) = (a1(2)*a2(3) - a1(3)*a2(2))/volum
255 728 : b3(2) = (a1(3)*a2(1) - a1(1)*a2(3))/volum
256 728 : b3(3) = (a1(1)*a2(2) - a1(2)*a2(1))/volum
257 : ! ==--------------------------------------------------------------==
258 2912 : DO i = 1, 3
259 2184 : b01(i) = b1(i)*alat
260 2184 : b02(i) = b2(i)*alat
261 2912 : 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 728 : v, f0, r, tvec, origin, rx, isc, delta)
269 : ! ==--------------------------------------------------------------==
270 26952 : DO n = nc + 1, 48
271 26952 : ib(n) = 0
272 : END DO
273 : ! ==--------------------------------------------------------------==
274 728 : invadd = 0
275 728 : IF (li == 0) THEN
276 436 : IF (iout > 0) THEN
277 : WRITE (iout, '(A,/,A,/,A)') &
278 187 : ' KPSYM| ALTHOUGH THE POINT GROUP OF THE CRYSTAL DOES NOT', &
279 187 : ' KPSYM| CONTAIN INVERSION, THE SPECIAL POINT GENERATION ALGORITHM', &
280 374 : ' KPSYM| WILL CONSIDER IT AS A SYMMETRY OPERATION'
281 : END IF
282 436 : invadd = 1
283 : END IF
284 : ! ==--------------------------------------------------------------==
285 : ! == CRYSTALLOGRAPHIC DATA ==
286 : ! ==--------------------------------------------------------------==
287 728 : IF (iout > 0) THEN
288 292 : WRITE (iout, '(/," KPSYM| CRYSTALLOGRAPHIC DATA:")')
289 292 : WRITE (iout, '(4X,"A1",3F10.5,10X,"B1",3F10.5)') a1, b1
290 292 : WRITE (iout, '(4X,"A2",3F10.5,10X,"B2",3F10.5)') a2, b2
291 292 : WRITE (iout, '(4X,"A3",3F10.5,10X,"B3",3F10.5)') a3, b3
292 : END IF
293 : ! ==--------------------------------------------------------------==
294 : ! == GROUP-THEORETICAL INFORMATION ==
295 : ! ==--------------------------------------------------------------==
296 728 : IF (iout > 0) THEN
297 292 : WRITE (iout, '(/," KPSYM| GROUP-THEORETICAL INFORMATION:")')
298 : END IF
299 : ! IHG .... Point group of the primitive lattice, holohedral
300 728 : IF (iout > 0) THEN
301 : WRITE (iout, &
302 : '(" KPSYM| POINT GROUP OF THE PRIMITIVE LATTICE: ",A," SYSTEM")') &
303 292 : 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 728 : IF (isy == 0) THEN
308 244 : IF (iout > 0) THEN
309 91 : WRITE (iout, '(" KPSYM|",4X,"NONSYMMORPHIC GROUP")')
310 : END IF
311 484 : ELSE IF (isy == 1) THEN
312 464 : IF (iout > 0) THEN
313 193 : 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 728 : IF (li == 0) THEN
326 436 : IF (iout > 0) THEN
327 187 : WRITE (iout, '(" KPSYM|",4X,"NO INVERSION SYMMETRY")')
328 : END IF
329 292 : ELSE IF (li > 0) THEN
330 292 : IF (iout > 0) THEN
331 105 : WRITE (iout, '(" KPSYM|",4X,"INVERSION SYMMETRY")')
332 : END IF
333 : END IF
334 : ! NC ..... Total number of elements in the point group
335 728 : IF (iout > 0) THEN
336 : WRITE (iout, &
337 292 : '(" KPSYM|",4X,"TOTAL NUMBER OF ELEMENTS IN THE POINT GROUP:",I3)') nc
338 : END IF
339 728 : IF (iout > 0) THEN
340 : WRITE (iout, '(" KPSYM|",4X,"TO SUM UP: (",I1,5I3,")")') &
341 292 : ihg, ihc, isy, li, nc, indpg
342 : END IF
343 : ! IB ..... List of the rotations constituting the point group
344 728 : IF (iout > 0) THEN
345 292 : WRITE (iout, '(/," KPSYM|",4X,"LIST OF THE ROTATIONS:")')
346 : END IF
347 728 : IF (iout > 0) THEN
348 292 : WRITE (iout, '(7X,12I4)') (ib(i), i=1, nc)
349 : END IF
350 : ! V ...... Nonprimitive translations (for nonsymmorphic groups)
351 728 : IF (isy <= 0) THEN
352 264 : IF (iout > 0) THEN
353 99 : WRITE (iout, '(/," KPSYM|",4X,"NONPRIMITIVE TRANSLATIONS:")')
354 : END IF
355 264 : IF (iout > 0) THEN
356 : WRITE (iout, '(A,A)') &
357 99 : ' ROT V IN THE BASIS A1, A2, A3 ', &
358 198 : 'V IN CARTESIAN COORDINATES'
359 : END IF
360 : ! Cartesian components of nonprimitive translation.
361 6814 : DO i = 1, nc
362 26200 : DO j = 1, 3
363 26200 : vv0(j) = v(1, i)*a1(j) + v(2, i)*a2(j) + v(3, i)*a3(j)
364 : END DO
365 6814 : IF (iout > 0) THEN
366 : WRITE (iout, '(1X,I3,3F10.5,3X,3F10.5)') &
367 1948 : 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 728 : IF (iout > 0) THEN
374 : WRITE (iout, &
375 292 : '(/," KPSYM|",4X,"ATOM TRANSFORMATION TABLE (MARADUDIN,VOSKO):")')
376 : END IF
377 728 : IF (iout > 0) THEN
378 292 : WRITE (iout, '(5(4X,"R AT->AT"))')
379 : END IF
380 728 : IF (iout > 0) THEN
381 292 : WRITE (iout, '(I5," [Identity]")') 1
382 : END IF
383 8720 : DO k = 2, nc
384 61676 : DO j = 1, nat
385 53684 : IF (iout > 0) THEN
386 15540 : WRITE (iout, '(I5,2I4)', advance="no") ib(k), j, f0(k, j)
387 : END IF
388 61676 : IF ((MOD(j, 5) == 0) .AND. iout > 0) THEN
389 1823 : WRITE (iout, *)
390 : END IF
391 : END DO
392 8720 : IF ((MOD(j - 1, 5) /= 0) .AND. iout > 0) THEN
393 2301 : WRITE (iout, *)
394 : END IF
395 : END DO
396 : ! R ...... List of the 3 x 3 rotation matrices
397 728 : IF (iout > 0) THEN
398 292 : WRITE (iout, '(/," KPSYM|",4X,"LIST OF THE 3 X 3 ROTATION MATRICES:")')
399 : END IF
400 728 : IF (ihc == 0) THEN
401 172 : DO k = 1, nc
402 172 : IF (iout > 0) THEN
403 : WRITE (iout, &
404 : '(4X,I3," (",I2,": ",A11,")",2(3F14.6,/,25X),3F14.6)') &
405 507 : 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 9276 : DO k = 1, nc
410 9276 : IF (iout > 0) THEN
411 : WRITE (iout, &
412 : '(4X,I3," (",I2,": ",A10,") ",2(3F14.6,/,25X),3F14.6)') &
413 33202 : 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 728 : 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 728 : IF (iout > 0) THEN
429 : WRITE (iout, '(/,1X,19("*"),A,25("*"))') &
430 292 : ' GENERATION OF SPECIAL POINTS '
431 : END IF
432 : ! Parameter Q of Monkhorst and Pack, generalized for 3 axes B1,2,3
433 728 : IF (iout > 0) THEN
434 : WRITE (iout, '(A,/,1X,3I5)') &
435 292 : ' KPSYM| MONKHORST-PACK PARAMETERS (GENERALIZED) IQ1,IQ2,IQ3:', &
436 584 : iq1, iq2, iq3
437 : END IF
438 : ! WVK0 is the shift of the whole mesh (see Macdonald)
439 728 : IF (iout > 0) THEN
440 : WRITE (iout, '(A,/,1X,3F10.5)') &
441 292 : ' KPSYM| CONSTANT VECTOR SHIFT (MACDONALD) OF THIS MESH:', wvk0
442 : END IF
443 728 : IF (ABS(iq1) + ABS(iq2) + ABS(iq3) == 0) RETURN
444 728 : 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 728 : IF (iout > 0) THEN
454 292 : WRITE (iout, '(" KPSYM| SYMMETRIZATION SWITCH: ",I3)', advance="no") istriz
455 : END IF
456 728 : IF (istriz == 1) THEN
457 728 : IF (iout > 0) THEN
458 292 : 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 101748 : DO i = 1, nkpoint
467 101748 : 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 728 : 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 728 : nhash, includ, list, rlist, delta)
516 : ! ==--------------------------------------------------------------==
517 : ! == Check on error signals ==
518 : ! ==--------------------------------------------------------------==
519 728 : IF (iout > 0) THEN
520 292 : WRITE (iout, '(/," KPSYM|",1X,I5," SPECIAL POINTS GENERATED")') ntot
521 : END IF
522 728 : IF (ntot == 0) THEN
523 : RETURN
524 728 : 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 728 : iswght = 0
536 2578 : DO i = 1, ntot
537 2578 : iswght = iswght + lwght(i)
538 : END DO
539 728 : IF (iout > 0) THEN
540 : WRITE (iout, '(8X,A,T33,A,4X,A)') &
541 292 : 'WAVEVECTOR K', 'WEIGHT', 'UNFOLDING ROTATIONS'
542 : END IF
543 : ! Set near-zeroes equal to zero:
544 2578 : DO l = 1, ntot
545 7400 : DO i = 1, 3
546 7400 : IF (ABS(wvkl(i, l)) < delta) wvkl(i, l) = 0._dp
547 : END DO
548 1850 : 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 1850 : lmax = lwght(l)
563 1850 : IF (iout > 0) THEN
564 : WRITE (iout, fmt='(1X,I5,3F8.4,I8,T42,12I3)') &
565 715 : l, (wvkl(i, l), i=1, 3), lwght(l), (lrot(i, l), i=1, MIN(lmax, 12))
566 : END IF
567 2758 : DO j = 13, lmax, 12
568 2030 : 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 728 : IF (iout > 0) THEN
575 292 : 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 2184 : SUBROUTINE group1s(iout, a1, a2, a3, nat, ty, x, b1, b2, b3, &
610 : ihg, ihc, isy, li, nc, indpg, ib, ntvec, &
611 2184 : 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 8736 : DO i = 1, 3
731 6552 : a(i, 1) = a1(i)
732 6552 : a(i, 2) = a2(i)
733 8736 : 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 2184 : CALL calbrec(a, ai)
746 8736 : DO i = 1, 3
747 6552 : b1(i) = ai(1, i)
748 6552 : b2(i) = ai(2, i)
749 8736 : 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 2184 : 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 2184 : CALL pgl1(ap, api, ihc, nc, ib, ihg, r, delta)
759 2184 : 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 340 : ncprim = nc
763 : ! The hexagonal system is found if the z axis is the sixfold axis
764 340 : CALL pgl1(a, ai, ihc, nc, ib, ihg, r, delta)
765 340 : 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 2184 : nc, indpg, ntvec, a, ai, li, isy, isc, delta)
774 :
775 2184 : IF (iout > 0) THEN
776 584 : IF (li > 0) THEN
777 : IF (iout > 0) THEN
778 : WRITE (iout, '(1X,A)') &
779 397 : 'KPSYM| THE POINT GROUP OF THE CRYSTAL CONTAINS THE INVERSION'
780 : END IF
781 : END IF
782 584 : IF (iout > 0) THEN
783 584 : WRITE (iout, *)
784 : END IF
785 : END IF
786 :
787 2184 : END SUBROUTINE group1s
788 : ! **************************************************************************************************
789 : !> \brief ...
790 : !> \param a ...
791 : !> \param ai ...
792 : ! **************************************************************************************************
793 2848 : 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 2848 : A(2, 1)*A(1, 2)*A(3, 3) - A(3, 1)*A(1, 3)*A(2, 2)
811 2848 : det = 1._dp/det
812 11392 : DO i = 1, 3
813 8544 : il = 1
814 8544 : iu = 3
815 8544 : IF (i == 1) il = 2
816 5696 : IF (i == 3) iu = 2
817 37024 : DO j = 1, 3
818 25632 : jl = 1
819 25632 : ju = 3
820 25632 : IF (j == 1) jl = 2
821 17088 : IF (j == 3) ju = 2
822 : ai(j, i) = (-1._dp)**(i + j)*det* &
823 34176 : (A(IL, JL)*A(IU, JU) - A(IL, JU)*A(IU, JL))
824 : END DO
825 : END DO
826 : ! ==--------------------------------------------------------------==
827 2848 : 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 2184 : 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 2184 : ntvec = 1
890 2184 : tvec(1, 1) = 0._dp
891 2184 : tvec(2, 1) = 0._dp
892 2184 : tvec(3, 1) = 0._dp
893 11336 : DO i = 1, nat
894 11336 : f0(49, i) = i
895 : END DO
896 9152 : DO k2 = 2, nat
897 6968 : IF (ty(1) /= ty(k2)) CYCLE
898 24864 : DO i = 1, 3
899 24864 : xb(i) = x(i, k2) - x(i, 1)
900 : END DO
901 : ! A fractional translation vector VR is defined.
902 6216 : CALL rlv3(ai, xb, vr, il, delta)
903 6216 : CALL checkrlv3(1, nat, ty, x, x, vr, f0, ai, isc, .TRUE., oksym, delta)
904 8400 : IF (oksym) THEN
905 : ! A fractional translational vector is found
906 988 : ntvec = ntvec + 1
907 : ! F0(49,1:NAT) gives number of equivalent atoms
908 : ! and has atom indexes of inequivalent atoms (for translation)
909 8796 : DO i = 1, nat
910 8796 : IF (f0(49, i) > f0(1, i)) f0(49, i) = f0(1, i)
911 : END DO
912 3952 : DO i = 1, 3
913 3952 : tvec(i, ntvec) = vr(i)
914 : END DO
915 : END IF
916 : END DO
917 : ! ==-------------------------------------------------------------==
918 8736 : DO i = 1, 3
919 6552 : ap(1, i) = a(1, i)
920 6552 : ap(2, i) = a(2, i)
921 6552 : ap(3, i) = a(3, i)
922 6552 : api(1, i) = ai(1, i)
923 6552 : api(2, i) = ai(2, i)
924 8736 : api(3, i) = ai(3, i)
925 : END DO
926 2184 : 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 1328 : DO iv = 2, ntvec
933 : ! TVEC in cartesian coordinates
934 3952 : DO i = 1, 3
935 : xb(i) = tvec(1, iv)*a(i, 1) &
936 : + TVEC(2, IV)*A(I, 2) &
937 3952 : + TVEC(3, IV)*A(I, 3)
938 : END DO
939 : ! We calculare TVEC in AP basis
940 988 : CALL rlv3(api, xb, vr, il, delta)
941 2624 : DO i = 1, 3
942 2284 : IF (ABS(vr(i)) > delta) THEN
943 664 : il = NINT(1._dp/ABS(vr(i)))
944 664 : IF (il > 1) THEN
945 : ! We replace AP(1:3,I) by TVEC(1:3,IV)
946 2656 : DO j = 1, 3
947 2656 : ap(j, i) = xb(j)
948 : END DO
949 : ! Calculate new API
950 664 : CALL calbrec(ap, api)
951 664 : EXIT
952 : END IF
953 : END IF
954 : END DO
955 : END DO
956 : END IF
957 : ! ==--------------------------------------------------------------==
958 2184 : 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 2524 : 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 5000 : DO ihc = 0, 1
1017 : ! IHC is 0 for hexagonal groups and 1 for cubic groups.
1018 5000 : IF (ihc == 0) THEN
1019 : nr = 24
1020 : ELSE
1021 2476 : nr = 48
1022 : END IF
1023 5000 : nc = 0
1024 : ! Constructs rotation operations.
1025 5000 : CALL rot1(ihc, r)
1026 184424 : loop_rotation: DO n = 1, nr
1027 179424 : ib(n) = 0
1028 : ! Rotate the A1,2,3 vectors by rotation No. N
1029 467216 : DO k = 1, 3
1030 1505536 : DO i = 1, 3
1031 1129152 : xa(i) = 0._dp
1032 4892992 : DO j = 1, 3
1033 4516608 : xa(i) = xa(i) + r(i, j, n)*a(j, k)
1034 : END DO
1035 : END DO
1036 376384 : CALL rlv3(ai, xa, vr, lx, delta)
1037 376384 : tr = 0._dp
1038 1505536 : DO i = 1, 3
1039 1505536 : tr = tr + ABS(vr(i))
1040 : END DO
1041 : ! If VR.ne.0, then XA cannot be a multiple of a lattice vector
1042 467216 : IF (tr > delta) CYCLE loop_rotation
1043 : END DO
1044 90832 : nc = nc + 1
1045 184424 : ib(nc) = n
1046 : END DO loop_rotation
1047 : ! ==------------------------------------------------------------==
1048 : ! IHG stands for holohedral group number.
1049 5000 : IF (ihc == 0) THEN
1050 : ! Hexagonal group:
1051 2524 : IF (nc == 12) ihg = 6
1052 2524 : IF (nc > 12) ihg = 7
1053 2524 : IF (nc >= 12) RETURN
1054 : ! Too few operations, try cubic group: (IHC=1,NR=48)
1055 : ELSE
1056 : ! Cubic group:
1057 2476 : IF (nc < 4) ihg = 1
1058 2476 : IF (nc == 4) ihg = 2
1059 2476 : IF (nc > 4) ihg = 3
1060 2476 : IF (nc == 16) ihg = 4
1061 2476 : IF (nc > 16) ihg = 5
1062 2476 : 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 2511992 : 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 2511992 : il = 0
1106 10047968 : DO i = 1, 3
1107 10047968 : vr(i) = 0._dp
1108 : END DO
1109 2511992 : ts = ABS(xb(1)) + ABS(xb(2)) + ABS(xb(3))
1110 2511992 : IF (ts <= delta) RETURN
1111 9274800 : DO i = 1, 3
1112 6956100 : vr(i) = vr(i) + ai(i, 1)*xb(1) + ai(i, 2)*xb(2) + ai(i, 3)*xb(3)
1113 6956100 : 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 9274800 : 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 2184 : SUBROUTINE atftm1(iout, r, v, x, f0, origin, ib, ty, nat, ihg, ihc, &
1148 2184 : 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 2184 : nodupli = ntvec == 1
1255 2184 : nca = 0
1256 107016 : DO n = 1, 48
1257 107016 : iis(n) = 0
1258 : END DO
1259 : ! Calculate translational vector for each operation
1260 : ! and atom transformation table.
1261 60816 : DO n = 1, nc
1262 58632 : l = ib(n)
1263 58632 : iis(l) = 1
1264 313440 : DO k = 1, nat
1265 1077864 : DO i = 1, 3
1266 1019232 : 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 234528 : DO k = 1, 3
1270 234528 : 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 58632 : CALL checkrlv3(n, nat, ty, rx, x, vr, f0, ai, isc, nodupli, oksym, delta)
1275 58632 : IF (.NOT. oksym) THEN
1276 : ! Now we try other possible VR
1277 : ! F0(49,1:NAT) has only inequivalent atom indexes for translation
1278 172860 : DO k2 = 1, nat
1279 151212 : IF (f0(49, k2) < k2) CYCLE
1280 132780 : IF (ty(1) /= ty(k2)) CYCLE
1281 465264 : DO i = 1, 3
1282 465264 : xb(i) = rx(i, 1) - x(i, k2)
1283 : END DO
1284 : ! A translation vector VR is defined.
1285 116316 : 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 116316 : CALL checkrlv3(n, nat, ty, rx, x, vr, f0, ai, isc, nodupli, oksym, delta)
1295 137964 : IF (oksym) EXIT
1296 : END DO
1297 28348 : IF (.NOT. oksym) THEN
1298 21648 : iis(l) = 0
1299 21648 : CYCLE
1300 : END IF
1301 : END IF
1302 36984 : nca = nca + 1
1303 150120 : DO i = 1, 3
1304 147936 : 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 2184 : i = 0
1316 2184 : ni = 13
1317 2184 : IF (ihg < 6) ni = 25
1318 2184 : li = 0
1319 60816 : DO n = 1, nc
1320 58632 : l = ib(n)
1321 58632 : IF (iis(l) == 0) CYCLE
1322 36984 : i = i + 1
1323 36984 : ib(i) = ib(n)
1324 36984 : IF (ib(i) == ni) li = i
1325 174504 : DO k = 1, nat
1326 193968 : f0(i, k) = f0(n, k)
1327 : END DO
1328 : END DO
1329 : ! ==--------------------------------------------------------------==
1330 2184 : nc = i
1331 2184 : vs = 0._dp
1332 39168 : DO n = 1, nc
1333 39168 : 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 2184 : IF (vs > delta) THEN
1339 528 : isy = 0
1340 : ELSE
1341 1656 : isy = 1
1342 : END IF
1343 : ! ==--------------------------------------------------------------==
1344 : ! Determination of the point group
1345 : ! (Thierry Deutsch - 1998 [Maybe not complete!!])
1346 2184 : IF (ihg < 6) THEN
1347 2136 : 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 2136 : ELSE IF (nc == 1) THEN
1354 : ! IB=1
1355 356 : indpg = 1 ! 1 (c1)
1356 1780 : ELSE IF (nc == 2 .AND. ib(2) == 25) THEN
1357 : ! IB=125
1358 32 : indpg = 2 ! <1>(ci)
1359 1748 : 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 1732 : 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 1612 : 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 318 : indpg = 5 ! 2/m(c2h)
1388 1294 : 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 1290 : 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 1290 : 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 1290 : 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 0 : indpg = 14 ! 4/m(c4h)
1421 1290 : 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 1286 : 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 1286 : 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 450 : indpg = 17 ! 4/mmm(d4h)
1445 836 : ELSE IF (nc == 4 .AND. (ib(4) == 4)) THEN
1446 : ! Orthorhombic system
1447 : ! IB=12 3 4
1448 0 : indpg = 25 ! 222(d2)
1449 836 : 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 836 : ELSE IF (nc == 8) THEN
1457 : ! IB=12 3 425 2627 28
1458 38 : indpg = 27 ! mmm(d2h)
1459 798 : 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 798 : 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 798 : 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 798 : 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 798 : ELSE IF (nc == 48) THEN
1481 : ! IB=1..48
1482 542 : 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 256 : indpg = -32
1489 : END IF
1490 : ELSE IF (ihg >= 6) THEN
1491 48 : 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 48 : ELSE IF (nc == 1) THEN
1498 : ! IB=1
1499 0 : indpg = 1 ! 1 (c1)
1500 48 : ELSE IF (nc == 2 .AND. ib(2) == 13) THEN
1501 : ! IB=113
1502 0 : indpg = 2 ! <1>(ci)
1503 48 : 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 48 : ELSE IF (nc == 2 .AND. ( &
1509 : ib(2) == 16)) THEN
1510 : ! IB=116
1511 0 : indpg = 4 ! m (c1h)
1512 48 : 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 48 : ELSE IF (nc == 3 .AND. ib(3) == 5) THEN
1519 : ! Trigonal system
1520 : ! IB=13 5
1521 8 : indpg = 6 ! 3 (c3)
1522 40 : ELSE IF (nc == 6 .AND. ib(6) == 17) THEN
1523 : ! IB=113 1517 35
1524 0 : indpg = 7 ! <3>(c3i)
1525 40 : ELSE IF (nc == 6 .AND. ib(6) == 11) THEN
1526 : ! IB=17 9 1135
1527 0 : indpg = 8 ! 32 (d3)
1528 40 : ELSE IF (nc == 6 .AND. ib(6) == 23) THEN
1529 : ! IB=13 5 1921 23
1530 0 : indpg = 9 ! 3m (c3v)
1531 40 : ELSE IF (nc == 12 .AND. ib(12) == 23) THEN
1532 : ! IB=13 5 79 1113 1517 1921 23
1533 32 : 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 2184 : IF (isy /= 1) THEN
1568 : ! Transform V in cartesian coordinates
1569 13628 : DO n = 1, nc
1570 13100 : vc(1, n) = a(1, 1)*v(1, n) + a(1, 2)*v(2, n) + a(1, 3)*v(3, n)
1571 13100 : vc(2, n) = a(2, 1)*v(1, n) + a(2, 2)*v(2, n) + a(2, 3)*v(3, n)
1572 13628 : 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 528 : CALL symmorphic(nc, ib, r, vc, ai, info, origin, delta)
1575 528 : 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 488 : ELSE IF (info == 0) THEN
1591 488 : isy = 0
1592 : ELSE
1593 0 : isy = -2
1594 : END IF
1595 : ELSE
1596 6624 : DO i = 1, 3
1597 6624 : origin(i) = 0._dp
1598 : END DO
1599 : END IF
1600 : ! ==--------------------------------------------------------------==
1601 : ! == Output ==
1602 : ! ==--------------------------------------------------------------==
1603 2184 : IF (iout > 0) THEN
1604 : IF (iout > 0) THEN
1605 584 : WRITE (iout, *)
1606 : END IF
1607 584 : CALL xstring(icst(ihg), i, j)
1608 584 : IF ((ihg == 7 .AND. nc == 24) .OR. &
1609 : (ihg == 5 .AND. nc == 48)) THEN
1610 115 : IF (iout > 0) THEN
1611 : WRITE (iout, '(A,A,A)') &
1612 115 : ' KPSYM| THE POINT GROUP OF THE CRYSTAL IS THE FULL ', &
1613 115 : icst(ihg) (i:j), &
1614 230 : ' GROUP'
1615 : END IF
1616 : ELSE
1617 469 : IF (iout > 0) THEN
1618 : WRITE (iout, '(A,A,A,I2,A)') &
1619 469 : ' KPSYM| THE CRYSTAL SYSTEM IS ', &
1620 469 : icst(ihg) (i:j), &
1621 938 : ' WITH ', nc, ' OPERATIONS:'
1622 : END IF
1623 469 : IF (ihc == 0) THEN
1624 6 : IF (iout > 0) THEN
1625 69 : WRITE (iout, '( 5(5(A13),/))') (rname_hexai(ib(i)), i=1, nc)
1626 : END IF
1627 : ELSE
1628 463 : IF (iout > 0) THEN
1629 4265 : WRITE (iout, '(10(5(A13),/))') (rname_cubic(ib(i)), i=1, nc)
1630 : END IF
1631 : END IF
1632 : END IF
1633 : ! ==------------------------------------------------------------==
1634 584 : IF (isy == 1) THEN
1635 485 : IF (iout > 0) THEN
1636 : WRITE (iout, '(A)') &
1637 485 : ' KPSYM| THE SPACE GROUP OF THE CRYSTAL IS SYMMORPHIC'
1638 : END IF
1639 99 : 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 91 : ELSE IF (isy == 0) THEN
1650 91 : IF (iout > 0) THEN
1651 : WRITE (iout, '(A,/,3X,A,F15.6,A)') &
1652 91 : ' KPSYM| THE SPACE GROUP IS NON-SYMMORPHIC,', &
1653 182 : ' (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 584 : IF (indpg > 0) THEN
1669 537 : CALL xstring(pgrp(indpg), i, j)
1670 537 : CALL xstring(pgrd(indpg), k, l)
1671 537 : IF (iout > 0) THEN
1672 : WRITE (iout, '(A,A,"(",A,")",T56,"[INDEX=",I2,"]")') &
1673 537 : ' KPSYM| THE POINT GROUP OF THE CRYSTAL IS ', pgrp(indpg) (i:j), &
1674 1074 : 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 584 : IF (ntvec == 1) THEN
1687 528 : IF (iout > 0) THEN
1688 : WRITE (iout, '(A,T60,I6)') &
1689 528 : ' KPSYM| NUMBER OF PRIMITIVE CELL:', ntvec
1690 : END IF
1691 : ELSE
1692 56 : IF (iout > 0) THEN
1693 : WRITE (iout, '(A,T60,I6)') &
1694 56 : ' KPSYM| NUMBER OF PRIMITIVE CELLS:', ntvec
1695 : END IF
1696 : END IF
1697 : END IF
1698 :
1699 2184 : 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 181164 : 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) :: vt(3), xb(3)
1755 :
1756 1375692 : DO ia = 1, nat
1757 1375692 : isc(ia) = 0
1758 : END DO
1759 : ! Now we check if ROT(N)+VR gives a correct symmetry.
1760 563432 : atom: DO ia = 1, nat
1761 2523132 : DO ib = 1, nat
1762 2523132 : IF (ty(ia) == ty(ib) .AND. isc(ib) == 0) THEN
1763 2011908 : xb(1) = rx(1, ia) - x(1, ib)
1764 2011908 : xb(2) = rx(2, ia) - x(2, ib)
1765 2011908 : xb(3) = rx(3, ia) - x(3, ib)
1766 2011908 : CALL rlv3(ai, xb, vt, il, delta)
1767 : ! VT STANDS FOR V-TEST
1768 : oksym = (ABS((vr(1) - vt(1)) - NINT(vr(1) - vt(1))) < delta) .AND. &
1769 : (ABS((vr(2) - vt(2)) - NINT(vr(2) - vt(2))) < delta) .AND. &
1770 2011908 : (ABS((vr(3) - vt(3)) - NINT(vr(3) - vt(3))) < delta)
1771 2011908 : IF (oksym) THEN
1772 382268 : IF (nodupli) isc(ib) = 1
1773 382268 : f0(n, ia) = ib
1774 : ! IR+VR is the good one: another symmetry operation
1775 : ! Next atom
1776 : CYCLE atom
1777 : END IF
1778 : END IF
1779 : END DO
1780 : ! VR is not the correct translation vector
1781 181164 : RETURN
1782 : END DO atom
1783 : END SUBROUTINE checkrlv3
1784 : ! ==================================================================
1785 : ! **************************************************************************************************
1786 : !> \brief ...
1787 : !> \param nc ...
1788 : !> \param ib ...
1789 : !> \param r ...
1790 : !> \param v ...
1791 : !> \param ai ...
1792 : !> \param info ...
1793 : !> \param origin ...
1794 : !> \param delta ...
1795 : ! **************************************************************************************************
1796 528 : SUBROUTINE symmorphic(nc, ib, r, v, ai, info, origin, delta)
1797 : ! ==--------------------------------------------------------------==
1798 : ! == Check if the group is symmorphic with a non-standard origin ==
1799 : ! == WARNING: If there are equivalent atoms, this routine could ==
1800 : ! == not determine if the space group is symmorphic ==
1801 : ! == So you have to check if the solution V=0 works (see ATFTM1) ==
1802 : ! ==--------------------------------------------------------------==
1803 : ! == INPUT: ==
1804 : ! == NC Number of operations ==
1805 : ! == IB(NC) Index of operation in R ==
1806 : ! == R(3,3,48) Rotations ==
1807 : ! == V(3,NC) Fractional translations related to R(3,3,IB(NC)) ==
1808 : ! == R AND V ARE IN CARTESIAN COORDINATES ==
1809 : ! == AI(I,J) ARE THE RECIPROCAL LATTICE VECTORS, ==
1810 : ! == B(I) = AI(I,J),J=1,2,3 ==
1811 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
1812 : ! == ==
1813 : ! == OUTPUT: ==
1814 : ! == ORIGIN(1:3) Give standard origin (cartesian coordinates) ==
1815 : ! == Give the standard origin with smallest coordinates==
1816 : ! == if NTVEC /= 1 ==
1817 : ! == INFO = 1 The group is symmorphic ==
1818 : ! == INFO = 0 The group is not symmorphic ==
1819 : ! == INFO =-1 The routine cannot determine ==
1820 : ! ==--------------------------------------------------------------==
1821 : INTEGER :: nc, ib(nc)
1822 : REAL(dp) :: r(3, 3, 48), v(3, nc), ai(3, 3)
1823 : INTEGER :: info
1824 : REAL(dp) :: origin(3), delta
1825 :
1826 : INTEGER :: i, i1, ierror, igood(3), il, imissing2, &
1827 : imissing3, iok(3), ionly, ir, j, j1
1828 : REAL(dp) :: diag, dif, r2(2, 2), r3(3, 3), vr(3), &
1829 : xb(3)
1830 :
1831 : ! Variables
1832 : ! ==--------------------------------------------------------------==
1833 : ! Find a point A / V_R = (1-R).OA
1834 :
1835 2112 : DO i = 1, 3
1836 2112 : iok(i) = 0
1837 : END DO
1838 2112 : DO i = 1, 3
1839 2112 : origin(i) = 0._dp
1840 : END DO
1841 4584 : DO ir = 1, nc
1842 4524 : dif = v(1, ir)*v(1, ir) + v(2, ir)*v(2, ir) + v(3, ir)*v(3, ir)
1843 4584 : IF (dif > delta*delta) THEN
1844 4816 : DO i = 1, 3
1845 4816 : igood(i) = 1
1846 : END DO
1847 : ! V is non-zero. Construct matrix 1-R
1848 4816 : DO i = 1, 3
1849 14448 : DO j = 1, 3
1850 14448 : r3(i, j) = -r(i, j, ib(ir))
1851 : END DO
1852 4816 : r3(i, i) = 1 + r3(i, i)
1853 : END DO
1854 1204 : CALL invmat(r3, ierror)
1855 1204 : IF (ierror == 0) THEN
1856 : ! The matrix 3x3 has an inverse.
1857 864 : DO i = 1, 3
1858 : vr(i) = r3(i, 1)*v(1, ir) &
1859 : + r3(i, 2)*v(2, ir) &
1860 864 : + r3(i, 3)*v(3, ir)
1861 : END DO
1862 : ELSE
1863 : ! IERROR gives the column which causes some trouble
1864 : ! Construct matrix 1-R with 2x2
1865 988 : igood(ierror) = 0
1866 988 : imissing3 = ierror
1867 988 : i1 = 0
1868 3952 : DO i = 1, 3
1869 3952 : IF (i /= ierror) THEN
1870 1976 : i1 = i1 + 1
1871 1976 : j1 = 0
1872 7904 : DO j = 1, 3
1873 7904 : IF (j /= ierror) THEN
1874 3952 : j1 = j1 + 1
1875 3952 : r2(i1, j1) = -r(i, j, ib(ir))
1876 : END IF
1877 : END DO
1878 1976 : r2(i1, i1) = 1 + r2(i1, i1)
1879 : END IF
1880 : END DO
1881 988 : CALL invmat(r2, ierror)
1882 988 : IF (ierror == 0) THEN
1883 : ! The matrix 2X2 has an inverse.
1884 : ! Solve Vxy = (1-R).OAxy + OAz R3z (z is IMISSING3)
1885 : i1 = 0
1886 3248 : DO i = 1, 3
1887 3248 : IF (igood(i) == 1) THEN
1888 1624 : i1 = i1 + 1
1889 1624 : vr(i) = 0._dp
1890 1624 : j1 = 0
1891 6496 : DO j = 1, 3
1892 6496 : IF (igood(j) == 1) THEN
1893 3248 : j1 = j1 + 1
1894 : vr(i) = vr(i) + r2(i1, j1)*(v(j, ir) + &
1895 3248 : origin(imissing3)*r(j, imissing3, ib(ir)))
1896 : END IF
1897 : END DO
1898 : ELSE
1899 812 : vr(i) = origin(i)
1900 : END IF
1901 : END DO
1902 : ELSE
1903 : ! Construct matrix 1-R with 1x1
1904 : i1 = 0
1905 704 : DO i = 1, 3
1906 704 : IF (i /= imissing3) THEN
1907 352 : i1 = i1 + 1
1908 352 : IF (i1 == ierror) THEN
1909 176 : igood(i) = 0
1910 176 : imissing2 = i
1911 : ELSE
1912 : ionly = i
1913 : END IF
1914 : END IF
1915 : END DO
1916 176 : diag = (1 - r(ionly, ionly, ib(ir)))
1917 176 : IF (ABS(diag) > delta) THEN
1918 : vr(ionly) = 1._dp/diag*(v(ionly, ir) + &
1919 : origin(imissing3)*r(ionly, imissing3, ib(ir)) + &
1920 176 : origin(imissing2)*r(ionly, imissing2, ib(ir)))
1921 : ELSE
1922 0 : vr(ionly) = origin(ionly)
1923 0 : igood(ionly) = 0
1924 : END IF
1925 176 : vr(imissing3) = origin(imissing3)
1926 176 : vr(imissing2) = origin(imissing2)
1927 : END IF
1928 : END IF
1929 : ! ==----------------------------------------------------------==
1930 : ! Compare VR with ORIGIN
1931 1204 : dif = 0._dp
1932 : ! If NTVEC /=1 there are NTVEC possible standard origins
1933 4816 : DO i = 1, 3
1934 4816 : IF (iok(i) == 1) THEN
1935 1528 : dif = dif + ABS(origin(i) - vr(i))
1936 : END IF
1937 : END DO
1938 1204 : IF (dif > delta) THEN
1939 : ! Non-symmorphic
1940 468 : info = 0
1941 468 : RETURN
1942 : ELSE
1943 2944 : DO i = 1, 3
1944 2944 : IF (iok(i) /= 1 .AND. igood(i) == 1) THEN
1945 1264 : iok(i) = 1
1946 1264 : origin(i) = vr(i)
1947 : END IF
1948 : END DO
1949 : END IF
1950 : END IF
1951 : END DO
1952 : ! ==--------------------------------------------------------------==
1953 60 : IF (iok(1) == 0 .AND. iok(2) == 0 .AND. iok(3) == 0) THEN
1954 : ! Cannot not determine
1955 0 : info = -1
1956 0 : RETURN
1957 : END IF
1958 : ! The group is symmorphic
1959 60 : info = 1
1960 : ! Check
1961 180 : DO ir = 1, nc
1962 560 : DO i = 1, 3
1963 : vr(i) = r(i, 1, ib(ir))*origin(1) &
1964 : + r(i, 2, ib(ir))*origin(2) &
1965 420 : + r(i, 3, ib(ir))*origin(3)
1966 560 : vr(i) = (origin(i) - vr(i)) - v(i, ir)
1967 : END DO
1968 140 : CALL rlv3(ai, vr, xb, il, delta)
1969 140 : dif = ABS(xb(1)) + ABS(xb(2)) + ABS(xb(3))
1970 180 : IF (dif > delta) THEN
1971 : ! Non-symmorphic
1972 20 : info = 0
1973 20 : RETURN
1974 : END IF
1975 : END DO
1976 : ! ==--------------------------------------------------------------==
1977 : RETURN
1978 : END SUBROUTINE symmorphic
1979 : ! ==================================================================
1980 : ! **************************************************************************************************
1981 : !> \brief ...
1982 : !> \param ihc ...
1983 : !> \param r ...
1984 : ! **************************************************************************************************
1985 5000 : SUBROUTINE rot1(ihc, r)
1986 : ! ==--------------------------------------------------------------==
1987 : ! == WRITTEN ON FEBRUARY 17TH, 1976 ==
1988 : ! == GENERATION OF THE X,Y,Z-TRANSFORMATION MATRICES 3X3 ==
1989 : ! == FOR HEXAGONAL AND CUBIC GROUPS ==
1990 : ! == SUBROUTINES NEEDED -- NONE ==
1991 : ! ==--------------------------------------------------------------==
1992 : ! == THIS IS IDENTICAL WITH THE SUBROUTINE ROT OF WORLTON-WARREN ==
1993 : ! == (IN THE AC-COMPLEX), ONLY THE WAY OF TRANSFERRING THE DATA ==
1994 : ! == WAS CHANGED ==
1995 : ! ==--------------------------------------------------------------==
1996 : ! == INPUT DATA: ==
1997 : ! == IHC SWITCH DETERMINING IF WE DESIRE ==
1998 : ! == THE HEXAGONAL GROUP(IHC=0) OR THE CUBIC GROUP (IHC=1) ==
1999 : ! == OUTPUT DATA: ==
2000 : ! == R...THE 3X3 MATRICES OF THE DESIRED COORDINATE REPRESENTATION==
2001 : ! == THEIR NUMBERING CORRESPONDS TO THE SYMMETRY ELEMENTS AS ==
2002 : ! == LISTE IN WORLTON-WARREN ==
2003 : ! == (COMPUT. PHYS. COMM. 3(1972) 88--117) ==
2004 : ! == FOR IHC=0 THE FIRST 24 MATRICES OF THE ARRAY R REPRESENT ==
2005 : ! == THE FULL HEXAGONAL GROUP D(6H) ==
2006 : ! == FOR IHC=1 THE FIRST 48 MATRICES OF THE ARRAY R REPRESENT ==
2007 : ! == THE FULL CUBIC GROUP O(H) ==
2008 : ! ==--------------------------------------------------------------==
2009 : INTEGER :: ihc
2010 : REAL(dp) :: r(3, 3, 48)
2011 :
2012 : INTEGER :: i, j, k, n, nv
2013 : REAL(dp) :: c, s
2014 :
2015 20000 : DO j = 1, 3
2016 65000 : DO i = 1, 3
2017 2220000 : DO n = 1, 48
2018 2205000 : r(i, j, n) = 0._dp
2019 : END DO
2020 : END DO
2021 : END DO
2022 5000 : IF (ihc == 0) THEN
2023 : ! ==------------------------------------------------------------==
2024 : ! DEFINE THE GENERATORS FOR THE ROTATION MATRICES--HEXAGONAL GROUP
2025 : ! ==------------------------------------------------------------==
2026 2524 : c = 0.5_dp
2027 2524 : s = 0.5_dp*SQRT(3.0_dp)
2028 2524 : r(1, 1, 2) = c
2029 2524 : r(1, 2, 2) = -s
2030 2524 : r(2, 1, 2) = s
2031 2524 : r(2, 2, 2) = c
2032 2524 : r(1, 1, 7) = -c
2033 2524 : r(1, 2, 7) = -s
2034 2524 : r(2, 1, 7) = -s
2035 2524 : r(2, 2, 7) = c
2036 17668 : DO n = 1, 6
2037 15144 : r(3, 3, n) = 1._dp
2038 15144 : r(3, 3, n + 18) = 1._dp
2039 15144 : r(3, 3, n + 6) = -1._dp
2040 17668 : r(3, 3, n + 12) = -1._dp
2041 : END DO
2042 : ! ==------------------------------------------------------------==
2043 : ! == GENERATE THE REST OF THE ROTATION MATRICES ==
2044 : ! ==------------------------------------------------------------==
2045 7572 : DO i = 1, 2
2046 5048 : r(i, i, 1) = 1._dp
2047 17668 : DO j = 1, 2
2048 10096 : r(i, j, 6) = r(j, i, 2)
2049 35336 : DO k = 1, 2
2050 20192 : r(i, j, 3) = r(i, j, 3) + r(i, k, 2)*r(k, j, 2)
2051 20192 : r(i, j, 8) = r(i, j, 8) + r(i, k, 2)*r(k, j, 7)
2052 30288 : r(i, j, 12) = r(i, j, 12) + r(i, k, 7)*r(k, j, 2)
2053 : END DO
2054 : END DO
2055 : END DO
2056 7572 : DO i = 1, 2
2057 17668 : DO j = 1, 2
2058 10096 : r(i, j, 5) = r(j, i, 3)
2059 35336 : DO k = 1, 2
2060 20192 : r(i, j, 4) = r(i, j, 4) + r(i, k, 2)*r(k, j, 3)
2061 20192 : r(i, j, 9) = r(i, j, 9) + r(i, k, 2)*r(k, j, 8)
2062 20192 : r(i, j, 10) = r(i, j, 10) + r(i, k, 12)*r(k, j, 3)
2063 30288 : r(i, j, 11) = r(i, j, 11) + r(i, k, 12)*r(k, j, 2)
2064 : END DO
2065 : END DO
2066 : END DO
2067 32812 : DO n = 1, 12
2068 30288 : nv = n + 12
2069 93388 : DO i = 1, 2
2070 212016 : DO j = 1, 2
2071 181728 : r(i, j, nv) = -r(i, j, n)
2072 : END DO
2073 : END DO
2074 : END DO
2075 : ELSE
2076 : ! ==------------------------------------------------------------==
2077 : ! == DEFINE THE GENERATORS FOR THE ROTATION MATRICES-CUBIC GROUP==
2078 : ! ==------------------------------------------------------------==
2079 2476 : r(1, 3, 9) = 1._dp
2080 2476 : r(2, 1, 9) = 1._dp
2081 2476 : r(3, 2, 9) = 1._dp
2082 2476 : r(1, 1, 19) = 1._dp
2083 2476 : r(2, 3, 19) = -1._dp
2084 2476 : r(3, 2, 19) = 1._dp
2085 9904 : DO i = 1, 3
2086 7428 : r(i, i, 1) = 1._dp
2087 32188 : DO j = 1, 3
2088 22284 : r(i, j, 20) = r(j, i, 19)
2089 22284 : r(i, j, 5) = r(j, i, 9)
2090 96564 : DO k = 1, 3
2091 66852 : r(i, j, 2) = r(i, j, 2) + r(i, k, 19)*r(k, j, 19)
2092 66852 : r(i, j, 16) = r(i, j, 16) + r(i, k, 9)*r(k, j, 19)
2093 89136 : r(i, j, 23) = r(i, j, 23) + r(i, k, 19)*r(k, j, 9)
2094 : END DO
2095 : END DO
2096 : END DO
2097 9904 : DO i = 1, 3
2098 32188 : DO j = 1, 3
2099 96564 : DO k = 1, 3
2100 66852 : r(i, j, 6) = r(i, j, 6) + r(i, k, 2)*r(k, j, 5)
2101 66852 : r(i, j, 7) = r(i, j, 7) + r(i, k, 16)*r(k, j, 23)
2102 66852 : r(i, j, 8) = r(i, j, 8) + r(i, k, 5)*r(k, j, 2)
2103 66852 : r(i, j, 10) = r(i, j, 10) + r(i, k, 2)*r(k, j, 9)
2104 66852 : r(i, j, 11) = r(i, j, 11) + r(i, k, 9)*r(k, j, 2)
2105 66852 : r(i, j, 12) = r(i, j, 12) + r(i, k, 23)*r(k, j, 16)
2106 66852 : r(i, j, 14) = r(i, j, 14) + r(i, k, 16)*r(k, j, 2)
2107 66852 : r(i, j, 15) = r(i, j, 15) + r(i, k, 2)*r(k, j, 16)
2108 66852 : r(i, j, 22) = r(i, j, 22) + r(i, k, 23)*r(k, j, 2)
2109 89136 : r(i, j, 24) = r(i, j, 24) + r(i, k, 2)*r(k, j, 23)
2110 : END DO
2111 : END DO
2112 : END DO
2113 9904 : DO i = 1, 3
2114 32188 : DO j = 1, 3
2115 96564 : DO k = 1, 3
2116 66852 : r(i, j, 3) = r(i, j, 3) + r(i, k, 5)*r(k, j, 12)
2117 66852 : r(i, j, 4) = r(i, j, 4) + r(i, k, 5)*r(k, j, 10)
2118 66852 : r(i, j, 13) = r(i, j, 13) + r(i, k, 23)*r(k, j, 11)
2119 66852 : r(i, j, 17) = r(i, j, 17) + r(i, k, 16)*r(k, j, 12)
2120 66852 : r(i, j, 18) = r(i, j, 18) + r(i, k, 16)*r(k, j, 10)
2121 89136 : r(i, j, 21) = r(i, j, 21) + r(i, k, 12)*r(k, j, 15)
2122 : END DO
2123 : END DO
2124 : END DO
2125 61900 : DO n = 1, 24
2126 59424 : nv = n + 24
2127 59424 : r(1, 1, nv) = -r(1, 1, n)
2128 59424 : r(1, 2, nv) = -r(1, 2, n)
2129 59424 : r(1, 3, nv) = -r(1, 3, n)
2130 59424 : r(2, 1, nv) = -r(2, 1, n)
2131 59424 : r(2, 2, nv) = -r(2, 2, n)
2132 59424 : r(2, 3, nv) = -r(2, 3, n)
2133 59424 : r(3, 1, nv) = -r(3, 1, n)
2134 59424 : r(3, 2, nv) = -r(3, 2, n)
2135 61900 : r(3, 3, nv) = -r(3, 3, n)
2136 : END DO
2137 : END IF
2138 : ! ==--------------------------------------------------------------==
2139 5000 : RETURN
2140 : END SUBROUTINE rot1
2141 : ! ==================================================================
2142 : ! **************************************************************************************************
2143 : !> \brief ...
2144 : !> \param iout ...
2145 : !> \param iq1 ...
2146 : !> \param iq2 ...
2147 : !> \param iq3 ...
2148 : !> \param wvk0 ...
2149 : !> \param nkpoint ...
2150 : !> \param a1 ...
2151 : !> \param a2 ...
2152 : !> \param a3 ...
2153 : !> \param b1 ...
2154 : !> \param b2 ...
2155 : !> \param b3 ...
2156 : !> \param inv ...
2157 : !> \param nc ...
2158 : !> \param ib ...
2159 : !> \param r ...
2160 : !> \param ntot ...
2161 : !> \param wvkl ...
2162 : !> \param lwght ...
2163 : !> \param lrot ...
2164 : !> \param ncbrav ...
2165 : !> \param ibrav ...
2166 : !> \param istriz ...
2167 : !> \param nhash ...
2168 : !> \param includ ...
2169 : !> \param list ...
2170 : !> \param rlist ...
2171 : !> \param delta ...
2172 : ! **************************************************************************************************
2173 728 : SUBROUTINE sppt2(iout, iq1, iq2, iq3, wvk0, nkpoint, &
2174 : a1, a2, a3, b1, b2, b3, &
2175 728 : inv, nc, ib, r, ntot, wvkl, lwght, lrot, &
2176 : ncbrav, ibrav, istriz, &
2177 728 : nhash, includ, list, rlist, delta)
2178 : ! ==--------------------------------------------------------------==
2179 : ! == WRITTEN ON SEPTEMBER 12-20TH, 1979 BY K.K. ==
2180 : ! == MODIFIED 26-MAY-82 BY OLE HOLM NIELSEN ==
2181 : ! == GENERATION OF SPECIAL POINTS FOR AN ARBITRARY LATTICE, ==
2182 : ! == FOLLOWING THE METHOD MONKHORST,PACK, ==
2183 : ! == PHYS. REV. B13 (1976) 5188 ==
2184 : ! == MODIFIED BY MACDONALD, PHYS. REV. B18 (1978) 5897 ==
2185 : ! == THE SUBROUTINE IS WRITTEN ASSUMING THAT THE POINTS ARE ==
2186 : ! == GENERATED IN THE RECIPROCAL SPACE. ==
2187 : ! == IF, HOWEVER, THE B1,B2,B3 ARE REPLACED BY A1,A2,A3, THEN ==
2188 : ! == SPECIAL POINTS IN THE DIRECT SPACE CAN BE PRODUCED, AS WELL. ==
2189 : ! == (NO MULTIPLICATION BY 2PI IS THEN NECESSARY.) ==
2190 : ! == IN THE CASE OF NONSYMMORPHIC GROUPS, THE APPLICATION IN THE ==
2191 : ! == DIRECT SPACE WOULD PROBABLY REQUIRE A CERTAIN CAUTION. ==
2192 : ! == SUBROUTINES NEEDED: BZDEFI,BZRDUC,INBZ,MESH ==
2193 : ! == IN THE CASES WHERE THE POINT GROUP OF THE CRYSTAL DOES NOT ==
2194 : ! == CONTAIN INVERSION. THE LATTER MAY BE ADDED IF WE WISH ==
2195 : ! == (SEE COMMENT TO THE SWITCH INV). ==
2196 : ! == REDUCTION TO THE 1ST BRILLOUIN ZONE IS DONE ==
2197 : ! == BY ADDING G-VECTORS TO FIND THE SHORTEST WAVE-VECTOR. ==
2198 : ! == THE ROTATIONS OF THE BRAVAIS LATTICE ARE APPLIED TO THE ==
2199 : ! == MONKHORST/PACK MESH IN ORDER TO FIND ALL K-POINTS ==
2200 : ! == THAT ARE RELATED BY SYMMETRY. (OLE HOLM NIELSEN) ==
2201 : ! ==--------------------------------------------------------------==
2202 : ! == INPUT DATA: ==
2203 : ! == IOUT: LOGICAL UNIT FOR OUTPUT ==
2204 : ! == IF (IOUT<=0) NO MESSAGE ==
2205 : ! == IQ1,IQ2,IQ3 .. PARAMETER Q OF MONKHORST AND PACK, ==
2206 : ! == GENERALIZED AND DIFFERENT FOR THE 3 DIRECTIONS B1, ==
2207 : ! == B2 AND B3 ==
2208 : ! == WVK0 ... THE 'ARBITRARY' SHIFT OF THE WHOLE MESH, DENOTED K0 ==
2209 : ! == IN MACDONALD. WVK0 = 0 CORRESPONDS TO THE ORIGINAL ==
2210 : ! == SCHEME OF MONKHORST AND PACK. ==
2211 : ! == UNITS: 2PI/(UNITS OF LENGTH USED IN A1, A2, A3), ==
2212 : ! == I.E. THE SAME UNITS AS THE GENERATED SPECIAL POINTS==
2213 : ! == NKPOINT .. VARIABLE DIMENSION OF THE (OUTPUT) ARRAYS WVKL, ==
2214 : ! == LWGHT,LROT, I.E. SPACE RESERVED FOR THE SPECIAL ==
2215 : ! == POINTS AND ACCESSORIES. ==
2216 : ! == NKPOINT HAS TO BE >= NTOT (TOTAL NUMBER OF SPECIAL==
2217 : ! == POINTS. THIS IS CHECKED BY THE SUBROUTINE. ==
2218 : ! == ISTRIZ . INDICATES WHETHER ADDITIONAL MESH POINTS SHOULD BE ==
2219 : ! == GENERATED BY APPLYING GROUP OPERATIONS TO THE MESH. ==
2220 : ! == ISTRIZ=+1 MEANS SYMMETRIZE ==
2221 : ! == ISTRIZ=-1 MEANS DO NOT SYMMETRIZE ==
2222 : ! == THE FOLLOWING INPUT DATA MAY BE OBTAINED FROM THE SBRT. ==
2223 : ! == B1,B2,B3 .. RECIPROCAL LATTICE VECTORS, NOT MULTIPLIED BY ==
2224 : ! == GROUP1: ANY 2PI (IN UNITS RECIPROCAL TO THOSE ==
2225 : ! == OF A1,A2,A3) ==
2226 : ! == INV .... CODE INDICATING WHETHER WE WISH TO ADD THE INVERSION==
2227 : ! == TO THE POINT GROUP OF THE CRYSTAL OR NOT (IN THE ==
2228 : ! == CASE THAT THE POINT GROUP DOES NOT CONTAIN ANY). ==
2229 : ! == INV=0 MEANS: DO NOT ADD INVERSION ==
2230 : ! == INV/=0 MEANS: ADD THE INVERSION ==
2231 : ! == INV/=0 SHOULD BE THE STANDARD CHOICE WHEN SPPT2 ==
2232 : ! == IS USED IN RECIPROCAL SPACE - IN ORDER TO MAKE ==
2233 : ! == USE OF THE HERMITICITY OF HAMILTONIAN. ==
2234 : ! == WHEN USED IN DIRECT SPACE, THE RIGHT CHOICE OF INV ==
2235 : ! == WILL DEPEND ON THE NATURE OF THE PHYSICAL PROBLEM. ==
2236 : ! == IN THE CASES WHERE THE INVERSION IS ADDED BY THE ==
2237 : ! == SWITCH INV, THE LIST IB WILL NOT BE MODIFIED BUT IN ==
2238 : ! == THE OUTPUT LIST LROT SOME OF THE OPERATIONS WILL ==
2239 : ! == APPEAR WITH NEGATIVE SIGN; THIS MEANS THAT THEY HAVE==
2240 : ! == TO BE APPLIED MULTIPLIED BY INVERSION. ==
2241 : ! == NC ..... TOTAL NUMBER OF ELEMENTS IN THE POINT GROUP OF THE ==
2242 : ! == CRYSTAL ==
2243 : ! == IB ..... LIST OF THE ROTATIONS CONSTITUTING THE POINT GROUP ==
2244 : ! == OF THE CRYSTAL. THE NUMBERING IS THAT DEFINED IN ==
2245 : ! == WORLTON AND WARREN, I.E. THE ONE MATERIALIZED IN THE==
2246 : ! == ARRAY R (SEE BELOW) ==
2247 : ! == ONLY THE FIRST NC ELEMENTS OF THE ARRAY IB ARE ==
2248 : ! == MEANINGFUL ==
2249 : ! == R ...... LIST OF THE 3 X 3 ROTATION MATRICES ==
2250 : ! == (XYZ REPRESENTATION OF THE O(H) OR D(6)H GROUPS) ==
2251 : ! == ALL 48 OR 24 MATRICES ARE LISTED. ==
2252 : ! == NCBRAV . TOTAL NUMBER OF ELEMENTS IN RBRAV ==
2253 : ! == IBRAV .. LIST OF NCBRAV OPERATIONS OF THE BRAVAIS LATTICE ==
2254 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
2255 : ! ==--------------------------------------------------------------==
2256 : ! == OUTPUT DATA: ==
2257 : ! == NTOT ... TOTAL NUMBER OF SPECIAL POINTS ==
2258 : ! == IF NTOT APPEARS NEGATIVE, THIS IS AN ERROR SIGNAL ==
2259 : ! == WHICH MEANS THAT THE DIMENSION NKPOINT WAS CHOSEN ==
2260 : ! == TOO SMALL SO THAT THE ARRAYS WVKL ETC. CANNOT ==
2261 : ! == ACCOMODATE ALL THE GENERATED SPECIAL POINTS. ==
2262 : ! == IN THIS CASE THE ARRAYS WILL BE FILLED UP TO NKPOINT==
2263 : ! == AND FURTHER GENERATION OF NEW POINTS WILL BE ==
2264 : ! == INTERRUPTED. ==
2265 : ! == WVKL ... LIST OF SPECIAL POINTS. ==
2266 : ! == CARTESIAN COORDINATES AND NOT MULTIPLIED BY 2*PI. ==
2267 : ! == ONLY THE FIRST NTOT VECTORS ARE MEANINGFUL ==
2268 : ! == ALTHOUGH NO 2 POINTS FROM THE LIST ARE EQUIVALENT ==
2269 : ! == BY SYMMETRY, THIS SUBROUTINE STILL HAS A KIND OF ==
2270 : ! == 'BEAUTY DEFECT': THE POINTS FINALLY ==
2271 : ! == SELECTED ARE NOT NECESSARILY SITUATED IN A ==
2272 : ! == 'COMPACT' IRREDUCIBLE BRILL.ZONE; THEY MIGHT LIE IN ==
2273 : ! == DIFFERENT IRREDUCIBLE PARTS OF THE B.Z. - BUT THEY ==
2274 : ! == DO REPRESENT AN IRREDUCIBLE SET FOR INTEGRATION ==
2275 : ! == OVER THE ENTIRE B.Z. ==
2276 : ! == LWGHT ... THE LIST OF WEIGHTS OF THE CORRESPONDING POINTS. ==
2277 : ! == THESE WEIGHTS ARE NOT NORMALIZED (JUST INTEGERS) ==
2278 : ! == LROT ... FOR EACH SPECIAL POINT THE 'UNFOLDING ROTATIONS' ==
2279 : ! == ARE LISTED. IF E.G. THE WEIGHT OF THE I-TH SPECIAL ==
2280 : ! == POINT IS LWGHT(I), THEN THE ROTATIONS WITH NUMBERS ==
2281 : ! == LROT(J,I), J=1,2,...,LWGHT(I) WILL 'SPREAD' THIS ==
2282 : ! == SINGLE POINT FROM THE IRREDUCIBLE PART OF B.Z. INTO ==
2283 : ! == SEVERAL POINTS IN AN ELEMENTARY UNIT CELL ==
2284 : ! == (PARALLELOPIPED) OF THE RECIPROCAL SPACE. ==
2285 : ! == SOME OPERATION NUMBERS IN THE LIST LROT MAY APPEAR ==
2286 : ! == NEGATIVE, THIS MEANS THAT THE CORRESPONDING ROTATION==
2287 : ! == HAS TO BE APPLIED WITH INVERSION (THE LATTER HAVING ==
2288 : ! == BEEN ARTIFICIALLY ADDED AS SYMMETRY OPERATION IN ==
2289 : ! == CASE INV/=0).NO OTHER EFFORT WAS TAKEN,TO RENUMBER==
2290 : ! == THE ROTATIONS WITH MINUS SIGN OR TO EXTEND THE ==
2291 : ! == LIST OF THE POINT-GROUP OPERATIONS IN THE LIST NB. ==
2292 : ! == INCLUD ... INTEGER ARRAY USED BY SPPT2 INCLUD(NKPOINT) ==
2293 : ! == THE FIRST BIT (0) IS USED BY THE ROUTINE. ==
2294 : ! == THE OTHER BITS GIVE THE K-POINT INDEX IN ==
2295 : ! == THE SPECIAL K-POINT TABLE. ==
2296 : ! ==--------------------------------------------------------------==
2297 : ! == NHASH USED BY MESH ROUTINE ==
2298 : ! == LIST INTEGER ARRAY USED BY MESH LIST(NHASH+NKPOINT) ==
2299 : ! == RLIST real(8) :: ARRAY USED BY MESH RLIST(3,NKPOINT) ==
2300 : ! ==--------------------------------------------------------------==
2301 : ! == Use bit manipulations functions ==
2302 : ! == IBSET(I,POS) sets the bit POS to 1 in I integer ==
2303 : ! == IBCLR(I,POS) clears the bit POS to 1 in I integer ==
2304 : ! == BTEST(I,POS) .TRUE. if bit POS is 1 in I integer ==
2305 : ! ==--------------------------------------------------------------==
2306 : INTEGER :: iout, iq1, iq2, iq3
2307 : REAL(dp) :: wvk0(3)
2308 : INTEGER :: nkpoint
2309 : REAL(dp) :: a1(3), a2(3), a3(3), b1(3), b2(3), b3(3)
2310 : INTEGER :: inv, nc, ib(48)
2311 : REAL(dp) :: r(3, 3, 48)
2312 : INTEGER :: ntot
2313 : REAL(dp) :: wvkl(3, nkpoint)
2314 : INTEGER :: lwght(nkpoint), lrot(48, nkpoint), &
2315 : ncbrav, ibrav(48), istriz, nhash, &
2316 : includ(nkpoint), list(nkpoint + nhash)
2317 : REAL(dp) :: rlist(3, nkpoint), delta
2318 :
2319 : INTEGER, PARAMETER :: no = 0, nrsdir = 100
2320 :
2321 : INTEGER :: i, i1, i2, i3, ibsign, igarb0, igarbage, &
2322 : igarbg, ii, imesh, iop, iplace, &
2323 : iremov, iwvk, j, jplace, k, n, nplane
2324 : REAL(dp) :: diff, proja(3), projb(3), &
2325 : rsdir(4, nrsdir), ur1, ur2, ur3, &
2326 : wva(3), wvk(3)
2327 :
2328 : ! ==--------------------------------------------------------------==
2329 :
2330 728 : ntot = 0
2331 101748 : DO i = 1, nkpoint
2332 101020 : lrot(1, i) = 1
2333 4849688 : DO j = 2, 48
2334 4848960 : lrot(j, i) = 0
2335 : END DO
2336 : END DO
2337 101748 : DO i = 1, nkpoint
2338 101748 : includ(i) = no
2339 : END DO
2340 2912 : DO i = 1, 3
2341 2912 : wva(i) = 0._dp
2342 : END DO
2343 : ! ==--------------------------------------------------------------==
2344 : ! == DEFINE THE 1ST BRILLOUIN ZONE ==
2345 : ! ==--------------------------------------------------------------==
2346 728 : CALL bzdefine(iout, b1, b2, b3, rsdir, nplane, delta)
2347 : ! ==--------------------------------------------------------------==
2348 : ! == Generation of the mesh (they are not multiplied by 2*pi) by ==
2349 : ! == the Monkhorst/Pack algorithm, supplemented by all rotations ==
2350 : ! ==--------------------------------------------------------------==
2351 : ! Initialize the list of vectors
2352 728 : iplace = -2
2353 : CALL mesh(iout, wva, iplace, igarb0, igarbg, nkpoint, nhash, &
2354 728 : list, rlist, delta)
2355 728 : imesh = 0
2356 2394 : DO i1 = 1, iq1
2357 5992 : DO i2 = 1, iq2
2358 15366 : DO i3 = 1, iq3
2359 10102 : ur1 = REAL(1 + iq1 - 2*i1, kind=dp)/REAL(2*iq1, kind=dp)
2360 10102 : ur2 = REAL(1 + iq2 - 2*i2, kind=dp)/REAL(2*iq2, kind=dp)
2361 10102 : ur3 = REAL(1 + iq3 - 2*i3, kind=dp)/REAL(2*iq3, kind=dp)
2362 40408 : DO i = 1, 3
2363 40408 : wvk(i) = ur1*b1(i) + ur2*b2(i) + ur3*b3(i) + wvk0(i)
2364 : END DO
2365 : ! Reduce WVK to the 1st Brillouin zone
2366 : CALL bzrduc(wvk, a1, a2, a3, b1, b2, b3, rsdir, &
2367 10102 : nrsdir, nplane, delta)
2368 13700 : IF (istriz == 1) THEN
2369 : ! Symmetrization of the k-points mesh.
2370 : ! Apply all the Bravais lattice operations to WVK
2371 373734 : DO iop = 1, ncbrav
2372 1454528 : DO i = 1, 3
2373 1090896 : wva(i) = 0._dp
2374 4727216 : DO j = 1, 3
2375 4363584 : wva(i) = wva(i) + r(i, j, ibrav(iop))*wvk(j)
2376 : END DO
2377 : END DO
2378 : ! Check that WVA is inside the 1 Bz.
2379 363632 : IF (.NOT. inside_bz(wva, rsdir, nplane, delta)) THEN
2380 0 : IF (iout > 0) THEN
2381 0 : WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
2382 : END IF
2383 0 : IF (iout > 0) THEN
2384 : WRITE (iout, '(A,3F10.4,/,A,3F10.4,A,/,A,I3,A)') &
2385 0 : ' THE VECTOR ', wva, &
2386 0 : ' GENERATED FROM ', wvk, ' IN THE BASIC MESH', &
2387 0 : ' BY ROTATION NO. ', ibrav(iop), ' IS OUTSIDE THE 1BZ'
2388 : END IF
2389 0 : CPABORT('SPPT2: VECTOR OUTSIDE THE 1BZ')
2390 : END IF
2391 : ! Place WVA in list
2392 363632 : iplace = 0
2393 : CALL mesh(iout, wva, iplace, igarb0, igarbg, &
2394 363632 : nkpoint, nhash, list, rlist, delta)
2395 : ! If WVA was new (and therefore inserted),
2396 : ! IPLACE is the number.
2397 363632 : IF (iplace > 0) imesh = iplace
2398 373734 : IF (iplace > nkpoint) THEN
2399 0 : IF (iout > 0) THEN
2400 0 : WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
2401 : END IF
2402 0 : IF (iout > 0) THEN
2403 0 : WRITE (iout, *) 'MESH SIZE EXCEEDS NKPOINT=', nkpoint
2404 : END IF
2405 0 : CPABORT('SPPT2: MESH SIZE EXCEEDED')
2406 : END IF
2407 : END DO
2408 : ELSE
2409 : ! Place WVK in list
2410 0 : iplace = 0
2411 : CALL mesh(iout, wvk, iplace, igarb0, igarbg, &
2412 0 : nkpoint, nhash, list, rlist, delta)
2413 0 : imesh = iplace
2414 0 : IF (iplace > nkpoint) THEN
2415 0 : IF (iout > 0) THEN
2416 0 : WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
2417 : END IF
2418 0 : IF (iout > 0) THEN
2419 0 : WRITE (iout, *) 'MESH SIZE EXCEEDS NKPOINT=', nkpoint
2420 : END IF
2421 0 : CPABORT('SPPT2: MESH SIZE EXCEEDED')
2422 : END IF
2423 : END IF
2424 : END DO
2425 : END DO
2426 : END DO
2427 : !deb
2428 : !deb get full mesh
2429 : !deb
2430 728 : IF (iout > 0) THEN
2431 : ! IMESH: Number of k points in the mesh.
2432 : WRITE (iout, &
2433 292 : '(" KPSYM| THE WAVEVECTOR MESH CONTAINS ",I5," POINTS")') imesh
2434 292 : WRITE (iout, '(" KPSYM| THE POINTS ARE:")')
2435 4016 : DO ii = 1, imesh
2436 3724 : i = ii
2437 : CALL mesh(iout, wva, i, igarb0, igarbg, nkpoint, nhash, &
2438 3724 : list, rlist, delta)
2439 4016 : IF (MOD(i, 2) == 1) THEN
2440 1868 : WRITE (iout, '(1X,I5,3F10.4)', advance="no") i, wva
2441 : ELSE
2442 1856 : WRITE (iout, '(1X,I5,3F10.4)') i, wva
2443 : END IF
2444 : END DO
2445 292 : WRITE (iout, *)
2446 : END IF
2447 : ! ==--------------------------------------------------------------==
2448 728 : IF (istriz == 1) THEN
2449 : ! Now figure out if any special point difference (K - K'') is an
2450 : ! integral multiple of a reciprocal-space vector
2451 728 : iremov = 0
2452 10738 : DO i = 1, (imesh - 1)
2453 10010 : iplace = i
2454 : CALL mesh(iout, wva, iplace, igarb0, igarbg, &
2455 10010 : nkpoint, nhash, list, rlist, delta)
2456 : ! Project WVA onto B1,2,3:
2457 10010 : proja(1) = 0._dp
2458 10010 : proja(2) = 0._dp
2459 10010 : proja(3) = 0._dp
2460 40040 : DO k = 1, 3
2461 30030 : proja(1) = proja(1) + wva(k)*a1(k)
2462 30030 : proja(2) = proja(2) + wva(k)*a2(k)
2463 40040 : proja(3) = proja(3) + wva(k)*a3(k)
2464 : END DO
2465 : ! Now loop over all the rest of the mesh points
2466 315742 : loop_mesh: DO j = (i + 1), imesh
2467 305004 : jplace = j
2468 : CALL mesh(iout, wvk, jplace, igarb0, igarbg, &
2469 305004 : nkpoint, nhash, list, rlist, delta)
2470 : ! Project WVK onto B1,2,3:
2471 305004 : projb(1) = 0._dp
2472 305004 : projb(2) = 0._dp
2473 305004 : projb(3) = 0._dp
2474 1220016 : DO k = 1, 3
2475 915012 : projb(1) = projb(1) + wvk(k)*a1(k)
2476 915012 : projb(2) = projb(2) + wvk(k)*a2(k)
2477 1220016 : projb(3) = projb(3) + wvk(k)*a3(k)
2478 : END DO
2479 : ! Check (PROJA - PROJB): Is it integral ?
2480 382970 : DO k = 1, 3
2481 382922 : diff = proja(k) - projb(k)
2482 382970 : IF (ABS(REAL(NINT(diff), kind=dp) - diff) > delta) CYCLE loop_mesh
2483 : END DO
2484 : ! DIFF is integral: remove WVK from mesh:
2485 : CALL remove(wvk, jplace, igarb0, igarbg, &
2486 48 : nkpoint, nhash, list, rlist, delta)
2487 : ! If WVK actually removed, increment IREMOV
2488 10058 : IF (jplace > 0) iremov = iremov + 1
2489 : END DO loop_mesh
2490 : END DO
2491 728 : IF (iremov > 0 .AND. iout > 0) THEN
2492 : WRITE (iout, '(A,A,/,A,1X,I6,A,/)') &
2493 1 : ' KPSYM| SOME OF THESE MESH POINTS ARE RELATED BY LATTICE ', &
2494 1 : 'TRANSLATION VECTORS', &
2495 2 : ' KPSYM|', iremov, ' OF THE MESH POINTS REMOVED.'
2496 : END IF
2497 : END IF
2498 : ! ==--------------------------------------------------------------==
2499 : ! == IN THE MESH OF WAVEVECTORS, NOW SEARCH FOR EQUIVALENT POINTS:==
2500 : ! == THE INVERSION (TIME REVERSAL !) MAY BE USED. ==
2501 : ! ==--------------------------------------------------------------==
2502 11466 : DO iwvk = 1, imesh
2503 : ! IF(INCLUD(IWVK) == YES) CYCLE
2504 10738 : IF (BTEST(includ(iwvk), 0)) CYCLE
2505 : ! IWVK has not been encountered previously: new special point,
2506 : ! (only if WVK is not a garbage vector, however.)
2507 : ! INCLUD(IWVK) = YES
2508 1898 : includ(iwvk) = IBSET(includ(iwvk), 0)
2509 1898 : iplace = iwvk
2510 : CALL mesh(iout, wvk, iplace, igarb0, igarbg, &
2511 1898 : nkpoint, nhash, list, rlist, delta)
2512 : ! Find out whether Wvk is in the garbage list
2513 : CALL garbag(wvk, igarbage, igarb0, &
2514 1898 : nkpoint, nhash, list, rlist, delta)
2515 1898 : IF (igarbage > 0) CYCLE
2516 1850 : ntot = ntot + 1
2517 : ! Give the index in the special k points table.
2518 1850 : includ(iwvk) = includ(iwvk) + ntot*2
2519 7400 : DO i = 1, 3
2520 7400 : wvkl(i, ntot) = wvk(i)
2521 : END DO
2522 1850 : lwght(ntot) = 1
2523 : ! ==-----------------------------------------------------------==
2524 : ! Find all the equivalent points (symmetry given by atoms)
2525 24224 : equivalent_points: DO n = 1, nc
2526 : ! Rotate:
2527 86584 : DO i = 1, 3
2528 64938 : wva(i) = 0._dp
2529 281398 : DO j = 1, 3
2530 259752 : wva(i) = wva(i) + r(i, j, ib(n))*wvk(j)
2531 : END DO
2532 : END DO
2533 : ibsign = +1
2534 10738 : DO
2535 : ! Find WVA in the list
2536 24604 : iplace = -1
2537 : CALL mesh(iout, wva, iplace, igarb0, igarbg, &
2538 24604 : nkpoint, nhash, list, rlist, delta)
2539 24604 : IF (iplace == 0) THEN
2540 48 : IF (istriz /= -1) THEN
2541 : ! Find out whether WVA is in the garbage list
2542 : CALL garbag(wva, igarbage, igarb0, &
2543 48 : nkpoint, nhash, list, rlist, delta)
2544 48 : IF (igarbage == 0) THEN
2545 : ! I think this case is impossible (NC <= NCBRAV)
2546 : ! Error message
2547 0 : IF (iout > 0) THEN
2548 0 : WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
2549 : END IF
2550 0 : IF (iout > 0) THEN
2551 : WRITE (iout, '(A,3F10.4,/,A,3F10.4,A,/,A,I3,A)') &
2552 0 : ' THE VECTOR ', wva, &
2553 0 : ' GENERATED FROM ', wvk, ' IN THE BASIC MESH', &
2554 0 : ' BY ROTATION NO. ', ib(n), ' IS NOT IN THE LIST'
2555 : END IF
2556 0 : CPABORT('SPPT2: VECTOR NOT IN THE LIST')
2557 : END IF
2558 : END IF
2559 : END IF
2560 24604 : IF (iplace /= 0 .OR. istriz /= -1) THEN
2561 : ! Find out whether WVA is in the garbage list
2562 : CALL garbag(wva, igarbage, igarb0, &
2563 24604 : nkpoint, nhash, list, rlist, delta)
2564 24604 : IF (igarbage > 0) CYCLE equivalent_points
2565 : ! Was WVA encountered before ?
2566 24556 : IF (.NOT. BTEST(includ(iplace), 0)) THEN
2567 : ! Increment weight.
2568 8840 : lwght(ntot) = lwght(ntot) + 1
2569 8840 : lrot(lwght(ntot), ntot) = ib(n)*ibsign
2570 : ! INCLUD(IPLACE) = YES
2571 8840 : includ(iplace) = IBSET(includ(iplace), 0)
2572 : ! This k-point is an image of a special k-point.
2573 : ! Put the index of the special k-point.
2574 8840 : includ(iplace) = includ(iplace) + ntot*2
2575 : END IF
2576 : END IF
2577 24556 : IF (ibsign == -1 .OR. inv == 0) CYCLE equivalent_points
2578 : ! The case where we also apply the inversion to WVA
2579 : ! Repeat the search, but for -WVA
2580 2958 : ibsign = -1
2581 33478 : DO i = 1, 3
2582 11832 : wva(i) = -wva(i)
2583 : END DO
2584 : END DO
2585 : END DO equivalent_points
2586 : END DO
2587 : ! ==--------------------------------------------------------------==
2588 : ! == TOTAL NUMBER OF SPECIAL POINTS: NTOT ==
2589 : ! == BEFORE USING THE LIST WVKL AS WAVE VECTORS, THEY HAVE TO BE ==
2590 : ! == MULTIPLIED BY 2*PI ==
2591 : ! == THE LIST OF WEIGHTS LWGHT IS NOT NORMALIZED ==
2592 : ! ==--------------------------------------------------------------==
2593 728 : IF (ntot > nkpoint) THEN
2594 0 : IF (iout > 0) THEN
2595 0 : WRITE (iout, *) 'IN SPPT2 NUMBER OF SPECIAL POINTS = ', ntot
2596 : END IF
2597 0 : IF (iout > 0) THEN
2598 0 : WRITE (iout, *) 'BUT NKPOINT = ', nkpoint
2599 : END IF
2600 0 : ntot = -1
2601 : END IF
2602 728 : IF (iout > 0) THEN
2603 : ! Write the index table relating k points in the mesh
2604 : ! with special k points
2605 : IF (iout > 0) THEN
2606 : WRITE (iout, '(/,A,4X,A)') &
2607 292 : ' KPSYM|', 'CROSS TABLE RELATING MESH POINTS WITH SPECIAL POINTS:'
2608 : END IF
2609 292 : IF (iout > 0) THEN
2610 292 : WRITE (iout, '(5(4X,"IK -> SK"))')
2611 : END IF
2612 4016 : DO i = 1, imesh
2613 3724 : iplace = includ(i)/2
2614 3724 : IF (iout > 0) THEN
2615 3724 : WRITE (iout, '(1X,I5,1X,I5)', advance="no") i, iplace
2616 : END IF
2617 4016 : IF ((MOD(i, 5) == 0) .AND. iout > 0) THEN
2618 571 : WRITE (iout, *)
2619 : END IF
2620 : END DO
2621 292 : IF ((MOD(j - 1, 5) /= 0) .AND. iout > 0) THEN
2622 292 : WRITE (iout, *)
2623 : END IF
2624 : END IF
2625 728 : END SUBROUTINE sppt2
2626 : ! **************************************************************************************************
2627 : !> \brief ...
2628 : !> \param iout ...
2629 : !> \param wvk ...
2630 : !> \param iplace ...
2631 : !> \param igarb0 ...
2632 : !> \param igarbg ...
2633 : !> \param nmesh ...
2634 : !> \param nhash ...
2635 : !> \param list ...
2636 : !> \param rlist ...
2637 : !> \param delta ...
2638 : ! **************************************************************************************************
2639 709600 : SUBROUTINE mesh(iout, wvk, iplace, igarb0, igarbg, &
2640 709600 : nmesh, nhash, list, rlist, delta)
2641 : ! ==--------------------------------------------------------------==
2642 : ! == MESH MAINTAINS A LIST OF VECTORS FOR PLACEMENT AND/OR LOOKUP ==
2643 : ! == ==
2644 : ! == ADDITIONAL ENTRY POINTS: REMOVE .... REMOVE VECTOR FROM LIST ==
2645 : ! == GARBAG .... WAS VECTOR REMOVED ? ==
2646 : ! == ==
2647 : ! == WVK ....... VECTOR ==
2648 : ! == IPLACE .... ON INPUT: -2 MEANS: INITIALIZE THE LIST ==
2649 : ! == (AND RETURN) ==
2650 : ! == -1 MEANS: FIND WVK IN THE LIST ==
2651 : ! == 0 MEANS: ADD WVK TO THE LIST ==
2652 : ! == >0 MEANS: RETURN WVK NO. IPLACE ==
2653 : ! == ON OUTPUT: THE POSITION ASSIGNED TO WVK ==
2654 : ! == (=0 IF WVK IS NOT IN THE LIST) ==
2655 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
2656 : ! ==--------------------------------------------------------------==
2657 : INTEGER :: iout
2658 : REAL(dp) :: wvk(3)
2659 : INTEGER :: iplace, igarb0, igarbg, nmesh, nhash, &
2660 : list(nhash + nmesh)
2661 : REAL(dp) :: rlist(3, nmesh), delta
2662 :
2663 : INTEGER, PARAMETER :: nil = 0
2664 :
2665 : INTEGER :: i, ihash, ipoint
2666 : INTEGER, SAVE :: istore
2667 : REAL(dp) :: delta1, rhash
2668 :
2669 : ! ==--------------------------------------------------------------==
2670 : ! == Initialization ==
2671 : ! ==--------------------------------------------------------------==
2672 :
2673 709600 : delta1 = 10._dp*delta
2674 709600 : IF (iplace <= -2) THEN
2675 834388 : DO i = 1, nhash + nmesh
2676 834388 : list(i) = nil
2677 : END DO
2678 728 : istore = 1
2679 : ! IGARB0 points to a linked list of removed WVKS (the garbage).
2680 728 : igarb0 = 0
2681 728 : igarbg = 0
2682 728 : RETURN
2683 : ! ==--------------------------------------------------------------==
2684 708872 : ELSE IF ((iplace > -2) .AND. (iplace <= 0)) THEN
2685 : ! The particular HASH function used in this case:
2686 : rhash = 0.7890_dp*wvk(1) &
2687 : + 0.6810_dp*wvk(2) &
2688 388236 : + 0.5811_dp*wvk(3) + delta
2689 388236 : ihash = INT(ABS(rhash)*REAL(nhash, kind=dp))
2690 388236 : ihash = MOD(ihash, nhash) + nmesh + 1
2691 : ! Search for WVK in linked list
2692 388236 : ipoint = list(ihash)
2693 580326 : DO i = 1, 100
2694 : ! List exhausted
2695 580326 : IF (ipoint == nil) EXIT
2696 : ! Compare WVK with this element
2697 1707542 : IF (ALL(ABS(wvk(:) - rlist(:, ipoint)) <= delta1)) THEN
2698 : ! WVK located
2699 377450 : IF (iplace == 0) RETURN
2700 : ! IPLACE=-1
2701 24556 : iplace = ipoint
2702 24556 : RETURN
2703 : END IF
2704 : ! Next element of list
2705 192090 : ihash = ipoint
2706 202876 : ipoint = list(ihash)
2707 : END DO
2708 10786 : IF (ipoint /= nil) THEN
2709 : ! List too long
2710 0 : IF (iout > 0) THEN
2711 : WRITE (iout, '(2A,/,A)') &
2712 0 : ' SUBROUTINE MESH *** FATAL ERROR *** LINKED LIST', &
2713 0 : ' TOO LONG ***', ' CHOOSE A BETTER HASH-FUNCTION'
2714 : END IF
2715 0 : CPABORT('MESH: WARNING')
2716 : END IF
2717 : ! WVK was not found
2718 10786 : IF (iplace == -1) THEN
2719 : ! IPLACE=-1 : search for WVK unsuccessful
2720 48 : iplace = 0
2721 48 : RETURN
2722 : ELSE
2723 : ! IPLACE=0: add WVK to the list
2724 10738 : list(ihash) = istore
2725 10738 : IF (istore > nmesh) THEN
2726 0 : IF (iout > 0) THEN
2727 0 : WRITE (iout, '(A)') 'SUBROUTINE MESH *** FATAL ERROR ***'
2728 : END IF
2729 0 : IF (iout > 0) THEN
2730 : WRITE (iout, '(A,I10,A,/,A,3F10.5)') &
2731 0 : ' ISTORE=', istore, ' EXCEEDS DIMENSIONS', &
2732 0 : ' WVK = ', wvk
2733 : END IF
2734 0 : CPABORT('MESH: WARNING')
2735 : END IF
2736 10738 : list(istore) = nil
2737 42952 : DO i = 1, 3
2738 42952 : rlist(i, istore) = wvk(i)
2739 : END DO
2740 10738 : istore = istore + 1
2741 10738 : iplace = istore - 1
2742 10738 : RETURN
2743 : END IF
2744 : ! WVK was found
2745 : ELSE
2746 : ! ==--------------------------------------------------------------==
2747 : ! == Return a wavevector (IPLACE > 0) ==
2748 : ! ==--------------------------------------------------------------==
2749 320636 : ipoint = iplace
2750 320636 : IF (ipoint >= istore) THEN
2751 0 : IF (iout > 0) THEN
2752 : WRITE (iout, '(A,/,A,I5,A,/)') &
2753 0 : ' SUBROUTINE MESH *** WARNING ***', &
2754 0 : ' IPLACE = ', iplace, &
2755 0 : ' IS BEYOND THE LISTS - WVK SET TO 1.0E38'
2756 : END IF
2757 0 : DO i = 1, 3
2758 0 : wvk(i) = 1.0e38_dp
2759 : END DO
2760 : END IF
2761 1282544 : DO i = 1, 3
2762 1282544 : wvk(i) = rlist(i, ipoint)
2763 : END DO
2764 : END IF
2765 : END SUBROUTINE mesh
2766 : ! **************************************************************************************************
2767 : !> \brief ...
2768 : !> \param wvk ...
2769 : !> \param iplace ...
2770 : !> \param igarb0 ...
2771 : !> \param igarbg ...
2772 : !> \param nmesh ...
2773 : !> \param nhash ...
2774 : !> \param list ...
2775 : !> \param rlist ...
2776 : !> \param delta ...
2777 : ! **************************************************************************************************
2778 48 : SUBROUTINE remove(wvk, iplace, igarb0, igarbg, &
2779 48 : nmesh, nhash, list, rlist, delta)
2780 : ! ==--------------------------------------------------------------==
2781 : ! == ENTRY POINT FOR REMOVING A WAVEVECTOR ==
2782 : ! == ==
2783 : ! == INPUT: ==
2784 : ! == WVK(3) ==
2785 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
2786 : ! == OUTPUT:
2787 : ! == IPLACE .....1 IF WVK WAS REMOVED ==
2788 : ! == 0 IF WVK WAS NOT REMOVED ==
2789 : ! == (WVK NOT IN THE LINKED LISTS) ==
2790 : ! ==--------------------------------------------------------------==
2791 : REAL(dp) :: wvk(3)
2792 : INTEGER :: iplace, igarb0, igarbg, nmesh, nhash, &
2793 : list(nhash + nmesh)
2794 : REAL(dp) :: rlist(3, nmesh), delta
2795 :
2796 : INTEGER, PARAMETER :: nil = 0
2797 :
2798 : INTEGER :: i, ihash, ipoint
2799 : REAL(dp) :: delta1, rhash
2800 :
2801 : ! ==--------------------------------------------------------------==
2802 : ! Variables
2803 : ! ==--------------------------------------------------------------==
2804 :
2805 48 : delta1 = 10._dp*delta
2806 : ! The particular hash function used in this case:
2807 : rhash = 0.7890_dp*wvk(1) &
2808 : + 0.6810_dp*wvk(2) &
2809 48 : + 0.5811_dp*wvk(3) + delta
2810 48 : ihash = INT(ABS(rhash)*REAL(nhash, kind=dp))
2811 48 : ihash = MOD(ihash, nhash) + nmesh + 1
2812 : ! Search for WVK in linked list
2813 48 : ipoint = list(ihash)
2814 96 : DO i = 1, 100
2815 : ! List exhausted
2816 96 : IF (ipoint == nil) THEN
2817 : ! WVK was not found in the mesh:
2818 0 : iplace = 0
2819 0 : RETURN
2820 : END IF
2821 : ! Compare WVK with this element
2822 240 : IF (.NOT. ANY(ABS(wvk(:) - rlist(:, ipoint)) > delta1)) THEN
2823 : ! WVK located, now remove it from the list:
2824 48 : list(ihash) = list(ipoint)
2825 : ! LIST(IHASH) now points to the next element in the list,
2826 : ! and the present WVK has become garbage.
2827 : ! Add WVK to the list of garbage:
2828 48 : IF (igarb0 == 0) THEN
2829 : ! Start up the garbage list:
2830 4 : igarb0 = ipoint
2831 : ELSE
2832 44 : list(igarbg) = ipoint
2833 : END IF
2834 48 : igarbg = ipoint
2835 48 : list(igarbg) = nil
2836 48 : iplace = 1
2837 48 : RETURN
2838 : END IF
2839 : ! Next element of list
2840 48 : ihash = ipoint
2841 48 : ipoint = list(ihash)
2842 : END DO
2843 : ! List too long
2844 0 : CPABORT('MESH: LIST TOO LONG')
2845 : END SUBROUTINE remove
2846 : ! **************************************************************************************************
2847 : !> \brief ...
2848 : !> \param wvk ...
2849 : !> \param iplace ...
2850 : !> \param igarb0 ...
2851 : !> \param nmesh ...
2852 : !> \param nhash ...
2853 : !> \param list ...
2854 : !> \param rlist ...
2855 : !> \param delta ...
2856 : ! **************************************************************************************************
2857 26550 : SUBROUTINE garbag(wvk, iplace, igarb0, &
2858 26550 : nmesh, nhash, list, rlist, delta)
2859 : ! ==--------------------------------------------------------------==
2860 : ! == ENTRY POINT FOR CHECKING IF A WAVEVECTOR ==
2861 : ! == IS IN THE GARBAGE LIST ==
2862 : ! == INPUT: ==
2863 : ! == WVK(3) ==
2864 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
2865 : ! == ==
2866 : ! == OUTPUT: ==
2867 : ! == IPLACE ..... I > 0 IS THE PLACE IN THE GARBAGE LIST ==
2868 : ! == 0 IF WVK NOT AMONG THE GARBAGE ==
2869 : ! ==--------------------------------------------------------------==
2870 : REAL(dp) :: wvk(3)
2871 : INTEGER :: iplace, igarb0, nmesh, nhash, &
2872 : list(nhash + nmesh)
2873 : REAL(dp) :: rlist(3, nmesh), delta
2874 :
2875 : INTEGER, PARAMETER :: nil = 0
2876 :
2877 : INTEGER :: i, ihash, ipoint
2878 : REAL(dp) :: delta1
2879 :
2880 : ! ==--------------------------------------------------------------==
2881 : ! Variables
2882 : ! ==--------------------------------------------------------------==
2883 :
2884 26550 : delta1 = 10._dp*delta
2885 : ! Search for WVK in linked list
2886 : ! Point to the garbage list
2887 26550 : ipoint = igarb0
2888 32814 : DO i = 1, nmesh
2889 : ! LIST EXHAUSTED
2890 32814 : IF (ipoint == nil) THEN
2891 : ! WVK was not found in the mesh:
2892 26406 : iplace = 0
2893 26406 : RETURN
2894 : END IF
2895 : ! Compare WVK with this element
2896 7432 : IF (.NOT. ANY(ABS(wvk(:) - rlist(:, ipoint)) > delta1)) THEN
2897 : ! WVK was located in the garbage list
2898 144 : iplace = i
2899 144 : RETURN
2900 : END IF
2901 : ! Next element of list
2902 6264 : ihash = ipoint
2903 6264 : ipoint = list(ihash)
2904 : END DO
2905 : ! List too long
2906 0 : CPABORT('GARBAG: LIST TOO LONG')
2907 : END SUBROUTINE garbag
2908 :
2909 : ! **************************************************************************************************
2910 : !> \brief ...
2911 : !> \param wvk ...
2912 : !> \param a1 ...
2913 : !> \param a2 ...
2914 : !> \param a3 ...
2915 : !> \param b1 ...
2916 : !> \param b2 ...
2917 : !> \param b3 ...
2918 : !> \param rsdir ...
2919 : !> \param nrsdir ...
2920 : !> \param nplane ...
2921 : !> \param delta ...
2922 : ! **************************************************************************************************
2923 10102 : SUBROUTINE bzrduc(wvk, a1, a2, a3, b1, b2, b3, rsdir, nrsdir, nplane, delta)
2924 : ! ==--------------------------------------------------------------==
2925 : ! == REDUCE WVK TO LIE ENTIRELY WITHIN THE 1ST BRILLOUIN ZONE ==
2926 : ! == BY ADDING B-VECTORS ==
2927 : ! == DELTA REQUIRED ACCURACY (1.e-6_dp IS A GOOD VALUE) ==
2928 : ! ==--------------------------------------------------------------==
2929 : REAL(dp) :: wvk(3), a1(3), a2(3), a3(3), b1(3), &
2930 : b2(3), b3(3)
2931 : INTEGER :: nrsdir
2932 : REAL(dp) :: rsdir(4, nrsdir)
2933 : INTEGER :: nplane
2934 : REAL(dp) :: delta
2935 :
2936 : INTEGER, PARAMETER :: nzones = 4, nnn = 2*nzones + 1, &
2937 : nn = nzones + 1
2938 :
2939 : INTEGER :: i, i1, i2, i3, n1, n2, n3, nn1, nn2, nn3
2940 : LOGICAL :: inside
2941 : REAL(dp) :: wb(3), wva(3)
2942 :
2943 : ! ==--------------------------------------------------------------==
2944 : ! Variables
2945 : ! Look around +/- "NZONES" to locate vector
2946 : ! NZONES may need to be increased for very anisotropic zones
2947 : ! ==--------------------------------------------------------------==
2948 :
2949 10102 : IF (.NOT. inside_bz(wvk, rsdir, nplane, delta)) THEN
2950 40 : inside = .FALSE.
2951 : ! Express WVK in the basis of B1,2,3.
2952 : ! This permits an estimate of how far WVK is from the 1Bz.
2953 40 : wb(1) = wvk(1)*a1(1) + wvk(2)*a1(2) + wvk(3)*a1(3)
2954 40 : wb(2) = wvk(1)*a2(1) + wvk(2)*a2(2) + wvk(3)*a2(3)
2955 40 : wb(3) = wvk(1)*a3(1) + wvk(2)*a3(2) + wvk(3)*a3(3)
2956 40 : nn1 = NINT(wb(1))
2957 40 : nn2 = NINT(wb(2))
2958 40 : nn3 = NINT(wb(3))
2959 : ! Look around the estimated vector for the one truly inside the 1Bz
2960 192 : n1_loop: DO n1 = 1, nnn
2961 192 : i1 = nn - n1 - nn1
2962 1768 : DO n2 = 1, nnn
2963 1576 : i2 = nn - n2 - nn2
2964 15712 : DO n3 = 1, nnn
2965 14024 : i3 = nn - n3 - nn3
2966 56096 : DO i = 1, 3
2967 : wva(i) = wvk(i) + REAL(i1, kind=dp)*b1(i) + REAL(i2, kind=dp)*b2(i) + &
2968 56096 : REAL(i3, kind=dp)*b3(i)
2969 : END DO
2970 14024 : inside = inside_bz(wva, rsdir, nplane, delta)
2971 15560 : IF (inside) EXIT n1_loop
2972 : END DO
2973 : END DO
2974 : END DO n1_loop
2975 40 : CPASSERT(inside)
2976 40 : wvk(1:3) = wva(1:3)
2977 : END IF
2978 :
2979 10102 : END SUBROUTINE bzrduc
2980 :
2981 : ! **************************************************************************************************
2982 : !> \brief Is wvk in the 1st Brillouin zone ?
2983 : !> Check whether wvk lies inside all the planes that define the 1bz.
2984 : !> \param wvk ...
2985 : !> \param rsdir ...
2986 : !> \param nplane ...
2987 : !> \param delta ...
2988 : !> \return ...
2989 : ! **************************************************************************************************
2990 387758 : FUNCTION inside_bz(wvk, rsdir, nplane, delta) RESULT(inbz)
2991 : REAL(KIND=dp), DIMENSION(3) :: wvk
2992 : REAL(KIND=dp), DIMENSION(:, :) :: rsdir
2993 : INTEGER :: nplane
2994 : REAL(KIND=dp) :: delta
2995 : LOGICAL :: inbz
2996 :
2997 : INTEGER :: n
2998 : REAL(KIND=dp) :: projct
2999 :
3000 387758 : inbz = .TRUE.
3001 1525032 : DO n = 1, nplane
3002 1151298 : projct = (rsdir(1, n)*wvk(1) + rsdir(2, n)*wvk(2) + rsdir(3, n)*wvk(3))/rsdir(4, n)
3003 1525032 : IF (ABS(projct) > 0.5_dp + delta) THEN
3004 : inbz = .FALSE.
3005 : EXIT
3006 : END IF
3007 : END DO
3008 :
3009 387758 : END FUNCTION inside_bz
3010 :
3011 : ! **************************************************************************************************
3012 : !> \brief Find the vectors whose halves define the 1st Brillouin zone
3013 : !> Output:
3014 : !> nplane -- How many elements of rsdir contain normal vectors defining the planes
3015 : !> Method:
3016 : !> Starting with the parallelopiped spanned by b1,2,3 around the origin,
3017 : !> vectors inside a sufficiently large sphere are tested to see whether
3018 : !> the planes at 1/2*b will further confine the 1bz.
3019 : !> The resulting vectors are not cleaned to avoid redundant planes
3020 : !> \param iout ...
3021 : !> \param b1 ...
3022 : !> \param b2 ...
3023 : !> \param b3 ...
3024 : !> \param rsdir ...
3025 : !> \param nplane ...
3026 : !> \param delta ...
3027 : ! **************************************************************************************************
3028 728 : SUBROUTINE bzdefine(iout, b1, b2, b3, rsdir, nplane, delta)
3029 : INTEGER :: iout
3030 : REAL(KIND=dp), DIMENSION(3) :: b1, b2, b3
3031 : REAL(KIND=dp), DIMENSION(:, :) :: rsdir
3032 : INTEGER :: nplane
3033 : REAL(KIND=dp) :: delta
3034 :
3035 : INTEGER :: i, i1, i2, i3, n, n1, n2, n3, nb1, nb2, &
3036 : nb3, nnb1, nnb2, nnb3, nrsdir
3037 : REAL(KIND=dp) :: b1len, b2len, b3len, bmax, projct
3038 : REAL(KIND=dp), DIMENSION(3) :: bvec
3039 :
3040 728 : nrsdir = SIZE(rsdir, 2)
3041 :
3042 728 : b1len = b1(1)**2 + b1(2)**2 + b1(3)**2
3043 728 : b2len = b2(1)**2 + b2(2)**2 + b2(3)**2
3044 728 : b3len = b3(1)**2 + b3(2)**2 + b3(3)**2
3045 : ! Lattice containing entirely the Brillouin zone
3046 728 : bmax = b1len + b2len + b3len
3047 728 : nb1 = INT(SQRT(bmax/b1len) + delta) + 1
3048 728 : nb2 = INT(SQRT(bmax/b2len) + delta) + 1
3049 728 : nb3 = INT(SQRT(bmax/b3len) + delta) + 1
3050 364728 : rsdir(:, :) = 0._dp
3051 : ! 1Bz is certainly confined inside the 1/2(B1,B2,B3) parallelopiped
3052 2912 : rsdir(1:3, 1) = b1(1:3)
3053 2912 : rsdir(1:3, 2) = b2(1:3)
3054 2912 : rsdir(1:3, 3) = b3(1:3)
3055 728 : rsdir(4, 1) = b1len
3056 728 : rsdir(4, 2) = b2len
3057 728 : rsdir(4, 3) = b3len
3058 : ! Starting confinement: 3 planes
3059 728 : nplane = 3
3060 728 : nnb1 = 2*nb1 + 1
3061 728 : nnb2 = 2*nb2 + 1
3062 728 : nnb3 = 2*nb3 + 1
3063 :
3064 4368 : DO n1 = 1, nnb1
3065 3640 : i1 = nb1 + 1 - n1
3066 29888 : DO n2 = 1, nnb2
3067 25520 : i2 = nb2 + 1 - n2
3068 186820 : inner_loop: DO n3 = 1, nnb3
3069 157660 : i3 = nb3 + 1 - n3
3070 157660 : IF (i1 == 0 .AND. i2 == 0 .AND. i3 == 0) CYCLE inner_loop
3071 627728 : DO i = 1, 3
3072 : bvec(i) = REAL(i1, kind=dp)*b1(i) + REAL(i2, kind=dp)*b2(i) + &
3073 627728 : REAL(i3, kind=dp)*b3(i)
3074 : END DO
3075 : ! Does the plane of 1/2*BVEC narrow down the 1Bz ?
3076 193624 : DO n = 1, nplane
3077 : projct = 0.5_dp*(rsdir(1, n)*bvec(1) + rsdir(2, n)*bvec(2) &
3078 193530 : + rsdir(3, n)*bvec(3))/rsdir(4, n)
3079 : ! 1/2*BVEC is outside the Bz - skip this direction
3080 : ! The 1.e-6_dp takes care of single points touching the Bz,
3081 : ! and of the -(plane)
3082 193624 : IF (ABS(projct) > 0.5_dp - delta) CYCLE inner_loop
3083 : END DO
3084 : ! 1/2*BVEC further confines the 1Bz - include into RSDIR
3085 94 : nplane = nplane + 1
3086 94 : CPASSERT(nplane <= nrsdir)
3087 376 : DO i = 1, 3
3088 376 : rsdir(i, nplane) = bvec(i)
3089 : END DO
3090 : ! Length squared
3091 183180 : rsdir(4, nplane) = bvec(1)**2 + bvec(2)**2 + bvec(3)**2
3092 : END DO inner_loop
3093 : END DO
3094 : END DO
3095 :
3096 728 : IF (iout > 0) THEN
3097 : WRITE (iout, '(A,I3,A,/,A,/,100(" KPSYM|",1X,3F10.4,/))') &
3098 292 : ' KPSYM| The 1st Brillouin zone is confined by (at most)', &
3099 292 : nplane, ' planes', &
3100 292 : ' KPSYM| as defined by the +/- halves of the vectors:', &
3101 584 : ((rsdir(i, n), i=1, 3), n=1, nplane)
3102 : END IF
3103 :
3104 728 : END SUBROUTINE bzdefine
3105 :
3106 : END MODULE kpsym
|
