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