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