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