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