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