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