1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2! Copyright 2010.  Los Alamos National Security, LLC. This material was    !
3! produced under U.S. Government contract DE-AC52-06NA25396 for Los Alamos !
4! National Laboratory (LANL), which is operated by Los Alamos National     !
5! Security, LLC for the U.S. Department of Energy. The U.S. Government has !
6! rights to use, reproduce, and distribute this software.  NEITHER THE     !
7! GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY,     !
8! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS         !
9! SOFTWARE.  If software is modified to produce derivative works, such     !
10! modified software should be clearly marked, so as not to confuse it      !
11! with the version available from LANL.                                    !
12!                                                                          !
13! Additionally, this program is free software; you can redistribute it     !
14! and/or modify it under the terms of the GNU General Public License as    !
15! published by the Free Software Foundation; version 2.0 of the License.   !
16! Accordingly, this program is distributed in the hope that it will be     !
17! useful, but WITHOUT ANY WARRANTY; without even the implied warranty of   !
18! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General !
19! Public License for more details.                                         !
20!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
21
22SUBROUTINE KBLDNEWH
23
24  USE CONSTANTS_MOD
25  USE SETUPARRAY
26  USE NEBLISTARRAY
27  USE XBOARRAY
28  USE NONOARRAY
29  USE UNIVARRAY
30  USE KSPACEARRAY
31  USE MYPRECISION
32
33  IMPLICIT NONE
34
35  INTEGER :: I, J, NEWJ, K, L, II, JJ, KK, MM, MP, NN, SUBI
36  INTEGER :: IBRA, IKET, LBRA, LKET, MBRA, MKET
37  INTEGER :: INDEX, INDI, INDJ
38  INTEGER :: SWITCH, PREVJ
39  INTEGER :: PBCI, PBCJ, PBCK
40  INTEGER :: BASISI(5), BASISJ(5), LBRAINC, LKETINC
41  INTEGER :: NORBI, KCOUNT
42  INTEGER :: KX, KY, KZ
43  REAL(LATTEPREC) :: KX0, KY0, KZ0
44  REAL(LATTEPREC) :: ALPHA, BETA, COSBETA, PHI, TMP, PERM
45  REAL(LATTEPREC) :: RIJ(3), MAGR2, MAGR, MAGRP, RCUTTB
46  REAL(LATTEPREC) :: MAXARRAY(20), MAXRCUT, MAXRCUT2
47  REAL(LATTEPREC) :: ANGFACTOR
48  REAL(LATTEPREC) :: KPOINT(3), KDOTL
49  REAL(LATTEPREC) :: AMMBRA, WIGLBRAMBRA
50  COMPLEX(LATTEPREC) :: BLOCH, KHTMP
51  REAL(LATTEPREC), EXTERNAL :: UNIVSCALE, WIGNERD, SLMMP, TLMMP, AM, BM
52
53  HK = CMPLX(ZERO)
54
55  INDEX = 0
56
57  ! Build diagonal elements
58  DO I = 1, NATS
59
60     SELECT CASE(BASIS(ELEMPOINTER(I)))
61
62     CASE("s")
63
64        INDEX = INDEX + 1
65        HK(INDEX, INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
66
67     CASE("p")
68
69        DO SUBI = 1, 3
70           INDEX = INDEX + 1
71           HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
72        ENDDO
73
74     CASE("d")
75
76        DO SUBI = 1, 5
77           INDEX = INDEX + 1
78           HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
79        ENDDO
80
81     CASE("f")
82
83        DO SUBI = 1, 7
84           INDEX = INDEX + 1
85           HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
86        ENDDO
87
88     CASE("sp")
89
90        DO SUBI = 1, 4
91
92           INDEX = INDEX + 1
93           IF (SUBI .EQ. 1) THEN
94              HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
95           ELSE
96              HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
97           ENDIF
98
99        ENDDO
100
101     CASE("sd")
102
103        DO SUBI = 1, 6
104
105           INDEX = INDEX + 1
106           IF (SUBI .EQ. 1) THEN
107              HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
108           ELSE
109              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
110           ENDIF
111
112        ENDDO
113
114     CASE("sf")
115
116        DO SUBI = 1, 8
117
118           INDEX = INDEX + 1
119           IF (SUBI .EQ. 1) THEN
120              HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
121           ELSE
122              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
123           ENDIF
124
125        ENDDO
126
127     CASE("pd")
128
129        DO SUBI = 1, 8
130
131           INDEX = INDEX + 1
132           IF (SUBI .LE. 3) THEN
133              HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
134           ELSE
135              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
136           ENDIF
137
138        ENDDO
139
140     CASE("pf")
141
142        DO SUBI = 1, 10
143
144           INDEX = INDEX + 1
145           IF (SUBI .LE. 3) THEN
146              HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
147           ELSE
148              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
149           ENDIF
150
151        ENDDO
152
153     CASE("df")
154
155        DO SUBI = 1, 12
156
157           INDEX = INDEX + 1
158           IF (SUBI .LE. 5) THEN
159              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
160           ELSE
161              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
162           ENDIF
163
164        ENDDO
165
166     CASE("spd")
167
168        DO SUBI = 1, 9
169
170           INDEX = INDEX + 1
171           IF (SUBI .EQ. 1) THEN
172              HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
173           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN
174              HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
175           ELSE
176              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
177           ENDIF
178
179        ENDDO
180
181     CASE("spf")
182
183        DO SUBI = 1, 11
184
185           INDEX = INDEX + 1
186           IF (SUBI .EQ. 1) THEN
187              HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
188           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN
189              HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
190           ELSE
191              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
192           ENDIF
193
194        ENDDO
195
196     CASE("sdf")
197
198        DO SUBI = 1, 13
199
200           INDEX = INDEX + 1
201           IF (SUBI .EQ. 1) THEN
202              HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
203           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 6) THEN
204              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
205           ELSE
206              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
207           ENDIF
208
209        ENDDO
210
211
212     CASE("pdf")
213
214        DO SUBI = 1, 15
215
216           INDEX = INDEX + 1
217           IF (SUBI .LE. 3) THEN
218              HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
219           ELSEIF (SUBI .GT. 3 .AND. SUBI .LE. 8) THEN
220              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
221           ELSE
222              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
223           ENDIF
224
225        ENDDO
226
227     CASE("spdf")
228
229        DO SUBI = 1, 16
230
231           INDEX = INDEX + 1
232           IF (SUBI .EQ. 1) THEN
233              HK(INDEX, INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
234           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN
235              HK(INDEX, INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
236           ELSEIF (SUBI .GT.  1 .AND. SUBI .LE. 4) THEN
237              HK(INDEX, INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
238           ELSEIF (SUBI .GT. 4 .AND. SUBI .LE. 9) THEN
239              HK(INDEX, INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
240           ELSE
241              HK(INDEX, INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
242           ENDIF
243
244        ENDDO
245
246     END SELECT
247
248  ENDDO
249
250  DO I = 2, NKTOT
251     DO J  = 1, HDIM
252        HK(J,J,I) = HK(J,J,1)
253     ENDDO
254  ENDDO
255
256  ! We assign the diagonal elements in ADDQDEP
257
258  !  IF (BASISTYPE .EQ. "NONORTHO") THEN
259
260  !     SMAT = ZERO
261
262  !     DO I = 1, HDIM
263  !        SMAT(I,I) = ONE
264  !     ENDDO
265
266  !  ENDIF
267
268  ! Loop over Kpoints
269
270
271  KX0 = PI*(ONE - REAL(NKX))/(REAL(NKX)*BOXDIMS(1)) - PI*KSHIFT(1)
272  KY0 = PI*(ONE - REAL(NKY))/(REAL(NKY)*BOXDIMS(2)) - PI*KSHIFT(2)
273  KZ0 = PI*(ONE - REAL(NKZ))/(REAL(NKZ)*BOXDIMS(3)) - PI*KSHIFT(3)
274
275  !  KCOUNT = 0
276
277  !  DO KX = 1, NKX
278
279  !     KPOINT(1) = KX0 + TWO*PI*REAL(KX-1)/(BOXDIMS(1)*REAL(NKX))
280
281  !     DO KY = 1, NKY
282
283  !        KPOINT(2) = KY0 + TWO*PI*REAL(KY-1)/(BOXDIMS(2)*REAL(NKY))
284
285  !        DO KZ = 1, NKZ
286
287  !           KPOINT(3) = KZ0 + TWO*PI*REAL(KZ-1)/(BOXDIMS(3)*REAL(NKZ))
288
289  !           KCOUNT = KCOUNT+1
290
291  !$OMP PARALLEL DO DEFAULT (NONE) &
292  !$OMP SHARED(KPOINT, KCOUNT, NATS, BASIS, ELEMPOINTER, TOTNEBTB, NEBTB) &
293  !$OMP SHARED(CR, BOXDIMS, HK, NOINT, ATELE, ELE1, ELE2) &
294  !$OMP SHARED(BOND, OVERL, MATINDLIST, BASISTYPE) &
295  !$OMP PRIVATE(I, J, K, NEWJ, BASISI, BASISJ, INDI, INDJ, PBCI, PBCJ, PBCK) &
296  !$OMP PRIVATE(RIJ, MAGR2, MAGR, MAGRP, PHI, ALPHA, BETA, COSBETA) &
297  !$OMP PRIVATE(LBRAINC, LBRA, MBRA, L, LKETINC, LKET, MKET) &
298  !$OMP PRIVATE(BLOCH, KDOTL) &
299  !$OMP PRIVATE(RCUTTB, IBRA, IKET, AMMBRA, WIGLBRAMBRA, ANGFACTOR, MP)
300
301  DO I = 1, NATS
302
303     ! Build the lists of orbitals on each atom
304
305     SELECT CASE(BASIS(ELEMPOINTER(I)))
306
307     CASE("s")
308        BASISI(1) = 0
309        BASISI(2) = -1
310     CASE("p")
311        BASISI(1) = 1
312        BASISI(2) = -1
313     CASE("d")
314        BASISI(1) = 2
315        BASISI(2) = -1
316     CASE("f")
317        BASISI(1) = 3
318        BASISI(2) = -1
319     CASE("sp")
320        BASISI(1) = 0
321        BASISI(2) = 1
322        BASISI(3) = -1
323     CASE("sd")
324        BASISI(1) = 0
325        BASISI(2) = 2
326        BASISI(3) = -1
327     CASE("sf")
328        BASISI(1) = 0
329        BASISI(2) = 3
330        BASISI(3) = -1
331     CASE("pd")
332        BASISI(1) = 1
333        BASISI(2) = 2
334        BASISI(3) = -1
335     CASE("pf")
336        BASISI(1) = 1
337        BASISI(2) = 3
338        BASISI(3) = -1
339     CASE("df")
340        BASISI(1) = 2
341        BASISI(2) = 3
342        BASISI(3) = -1
343     CASE("spd")
344        BASISI(1) = 0
345        BASISI(2) = 1
346        BASISI(3) = 2
347        BASISI(4) = -1
348     CASE("spf")
349        BASISI(1) = 0
350        BASISI(2) = 1
351        BASISI(3) = 3
352        BASISI(4) = -1
353     CASE("sdf")
354        BASISI(1) = 0
355        BASISI(2) = 2
356        BASISI(3) = 3
357        BASISI(4) = -1
358     CASE("pdf")
359        BASISI(1) = 1
360        BASISI(2) = 2
361        BASISI(3) = 3
362        BASISI(4) = -1
363     CASE("spdf")
364        BASISI(1) = 0
365        BASISI(2) = 1
366        BASISI(3) = 2
367        BASISI(4) = 3
368        BASISI(5) = -1
369     END SELECT
370
371     INDI = MATINDLIST(I)
372
373     ! open loop over neighbors J of atom I
374     DO NEWJ = 1, TOTNEBTB(I)
375
376        J = NEBTB(1, NEWJ, I)
377
378        !                                  IF ( J .GE. I ) THEN
379
380        PBCI = NEBTB(2, NEWJ, I)
381        PBCJ = NEBTB(3, NEWJ, I)
382        PBCK = NEBTB(4, NEWJ, I)
383
384        RIJ(1) = CR(1,J) + REAL(PBCI) * BOXDIMS(1) &
385             - CR(1,I)
386        RIJ(2) = CR(2,J) + REAL(PBCJ) * BOXDIMS(2) &
387             - CR(2,I)
388        RIJ(3) = CR(3,J) + REAL(PBCK) * BOXDIMS(3) &
389             - CR(3,I)
390
391        MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3)
392
393        RCUTTB = ZERO
394
395        DO K = 1, NOINT
396
397           IF ( (ATELE(I) .EQ. ELE1(K) .AND. &
398                ATELE(J) .EQ. ELE2(K)) .OR. &
399                (ATELE(J) .EQ. ELE1(K) .AND. &
400                ATELE(I) .EQ. ELE2(K) )) THEN
401
402              IF (BOND(8,K) .GT. RCUTTB ) RCUTTB = BOND(8,K)
403
404              IF (BASISTYPE .EQ. "NONORTHO") THEN
405                 IF (OVERL(8,K) .GT. RCUTTB ) RCUTTB = OVERL(8,K)
406              ENDIF
407
408           ENDIF
409
410        ENDDO
411
412        IF (MAGR2 .LT. RCUTTB*RCUTTB) THEN
413
414           !                    KDOTL = KPOINT(1)*RIJ(1) + KPOINT(2)*RIJ(2) + &
415           !                         KPOINT(3)*RIJ(3)
416
417           !                    BLOCH = EXP(CMPLX(ZERO,KDOTL))
418
419           MAGR = SQRT(MAGR2)
420
421           SELECT CASE(BASIS(ELEMPOINTER(J)))
422           CASE("s")
423              BASISJ(1) = 0
424              BASISJ(2) = -1
425           CASE("p")
426              BASISJ(1) = 1
427              BASISJ(2) = -1
428           CASE("d")
429              BASISJ(1) = 2
430              BASISJ(2) = -1
431           CASE("f")
432              BASISJ(1) = 3
433              BASISJ(2) = -1
434           CASE("sp")
435              BASISJ(1) = 0
436              BASISJ(2) = 1
437              BASISJ(3) = -1
438           CASE("sd")
439              BASISJ(1) = 0
440              BASISJ(2) = 2
441              BASISJ(3) = -1
442           CASE("sf")
443              BASISJ(1) = 0
444              BASISJ(2) = 3
445              BASISJ(3) = -1
446           CASE("pd")
447              BASISJ(1) = 1
448              BASISJ(2) = 2
449              BASISJ(3) = -1
450           CASE("pf")
451              BASISJ(1) = 1
452              BASISJ(2) = 3
453              BASISJ(3) = -1
454           CASE("df")
455              BASISJ(1) = 2
456              BASISJ(2) = 3
457              BASISJ(3) = -1
458           CASE("spd")
459              BASISJ(1) = 0
460              BASISJ(2) = 1
461              BASISJ(3) = 2
462              BASISJ(4) = -1
463           CASE("spf")
464              BASISJ(1) = 0
465              BASISJ(2) = 1
466              BASISJ(3) = 3
467              BASISJ(4) = -1
468           CASE("sdf")
469              BASISJ(1) = 0
470              BASISJ(2) = 2
471              BASISJ(3) = 3
472              BASISJ(4) = -1
473           CASE("pdf")
474              BASISJ(1) = 1
475              BASISJ(2) = 2
476              BASISJ(3) = 3
477              BASISJ(4) = -1
478           CASE("spdf")
479              BASISJ(1) = 0
480              BASISJ(2) = 1
481              BASISJ(3) = 2
482              BASISJ(4) = 3
483              BASISJ(5) = -1
484           END SELECT
485
486           INDJ = MATINDLIST(J)
487
488           MAGRP = SQRT(RIJ(1) * RIJ(1) + RIJ(2) * RIJ(2))
489
490           ! transform to system in which z-axis is aligned with RIJ,
491           IF (ABS(RIJ(1)) .GT. 1E-12) THEN
492
493              IF (RIJ(1) .GT. ZERO .AND. RIJ(2) .GE. ZERO) THEN
494                 PHI = ZERO
495              ELSEIF (RIJ(1) .GT. ZERO .AND. RIJ(2) .LT. ZERO) THEN
496                 PHI = TWO * PI
497              ELSE
498                 PHI = PI
499              ENDIF
500              ALPHA = ATAN(RIJ(2) / RIJ(1)) + PHI
501
502           ELSEIF (ABS(RIJ(2)) .GT. 1E-12) THEN
503
504              IF (RIJ(2) .GT. 1E-12) THEN
505                 ALPHA = PI / TWO
506              ELSE
507                 ALPHA = THREE * PI / TWO
508              ENDIF
509
510           ELSE
511              ! pathological case: beta=0 and alpha undefined, but
512              ! this doesn't matter for matrix elements
513
514              ALPHA = ZERO
515
516           ENDIF
517
518           COSBETA = RIJ(3)/MAGR
519           BETA = ACOS(RIJ(3) / MAGR)
520
521           ! Build matrix elements using eqns (1)-(9) in PRB 72 165107
522
523           ! The loops over LBRA and LKET need to take into account
524           ! the orbitals assigned to each atom, e.g., sd rather than
525           ! spd...
526
527           IBRA = INDI + 1
528
529           LBRAINC = 1
530           DO WHILE (BASISI(LBRAINC) .NE. -1)
531
532              LBRA = BASISI(LBRAINC)
533              LBRAINC = LBRAINC + 1
534
535              DO MBRA = -LBRA, LBRA
536
537                 ! We can calculate these two outside the
538                 ! MKET loop...
539
540                 AMMBRA = AM(MBRA, ALPHA)
541                 WIGLBRAMBRA = WIGNERD(LBRA, ABS(MBRA), 0, COSBETA)
542
543                 IKET = INDJ + 1
544
545                 LKETINC = 1
546                 DO WHILE (BASISJ(LKETINC) .NE. -1)
547
548                    LKET = BASISJ(LKETINC)
549                    LKETINC = LKETINC + 1
550
551                    DO MKET = -LKET, LKET
552
553                       ! This is the sigma bonds (mp = 0)
554
555                       ! Hamiltonian build
556
557                       ! Pre-compute the angular part so we can use it
558                       ! again later if we're building the S matrix too
559
560                       ANGFACTOR = TWO * AMMBRA * &
561                            AM(MKET, ALPHA) * &
562                            WIGLBRAMBRA * &
563                            WIGNERD(LKET, ABS(MKET), 0, COSBETA)
564
565                       KHTMP = CMPLX(ANGFACTOR * &
566                            UNIVSCALE(I, J, LBRA, LKET, &
567                            0, MAGR, "H"))
568
569                       KCOUNT = 0
570
571                       DO KX = 1, NKX
572                          KPOINT(1) = KX0 + TWO*PI*REAL(KX-1)/(BOXDIMS(1)*REAL(NKX))
573                          DO KY = 1, NKY
574                             KPOINT(2) = KY0 + TWO*PI*REAL(KY-1)/(BOXDIMS(2)*REAL(NKY))
575                             DO KZ = 1, NKZ
576                                KPOINT(3) = KZ0 + TWO*PI*REAL(KZ-1)/(BOXDIMS(3)*REAL(NKZ))
577
578                                KCOUNT = KCOUNT+1
579
580                                KDOTL = KPOINT(1)*RIJ(1) + KPOINT(2)*RIJ(2) + &
581                                     KPOINT(3)*RIJ(3)
582
583                                BLOCH = EXP(CMPLX(ZERO,KDOTL))
584
585                                HK(IBRA, IKET, KCOUNT) = &
586                                     HK(IBRA, IKET, KCOUNT) + &
587                                     BLOCH*KHTMP
588
589                             ENDDO
590                          ENDDO
591                       ENDDO
592
593
594                       ! Overlap matrix build
595
596                       !                          IF (BASISTYPE .EQ. "NONORTHO") THEN
597
598                       !                             SMAT(IBRA, IKET) = SMAT(IBRA, IKET) + &
599                       !                                  BLOCH*ANGFACTOR * &
600                       !                                  UNIVSCALE(I, J, LBRA, LKET, 0, MAGR, "S")
601
602                       !                          ENDIF
603
604                       ! everything else
605
606                       DO MP = 1, MIN(LBRA, LKET)
607
608                          ANGFACTOR = &
609                               SLMMP(LBRA, MBRA, MP, ALPHA, COSBETA)* &
610                               SLMMP(LKET, MKET, MP, ALPHA, COSBETA)+ &
611                               TLMMP(LBRA, MBRA, MP, ALPHA, COSBETA)* &
612                               TLMMP(LKET, MKET, MP, ALPHA, COSBETA)
613
614                          KHTMP = CMPLX(ANGFACTOR * &
615                               UNIVSCALE(I, J, LBRA, LKET, &
616                               MP, MAGR, "H"))
617
618                          KCOUNT = 0
619                          DO KX = 1, NKX
620                             KPOINT(1) = KX0 + TWO*PI*REAL(KX-1)/(BOXDIMS(1)*REAL(NKX))
621                             DO KY = 1, NKY
622                                KPOINT(2) = KY0 + TWO*PI*REAL(KY-1)/(BOXDIMS(2)*REAL(NKY))
623                                DO KZ = 1, NKZ
624                                   KPOINT(3) = KZ0 + TWO*PI*REAL(KZ-1)/(BOXDIMS(3)*REAL(NKZ))
625
626                                   KCOUNT = KCOUNT+1
627
628                                   KDOTL = KPOINT(1)*RIJ(1) + KPOINT(2)*RIJ(2) + &
629                                        KPOINT(3)*RIJ(3)
630
631                                   BLOCH = EXP(CMPLX(ZERO,KDOTL))
632
633                                   HK(IBRA, IKET, KCOUNT) = &
634                                        HK(IBRA, IKET, KCOUNT) + &
635                                        BLOCH*KHTMP
636
637                                ENDDO
638                             ENDDO
639                          ENDDO
640
641
642
643                          !   HK(IBRA, IKET, KCOUNT) = &
644                          !        HK(IBRA, IKET, KCOUNT) + &
645                          !        BLOCH*CMPLX(ANGFACTOR * &
646                          !        UNIVSCALE(I, J, LBRA, LKET, &
647                          !        MP, MAGR, "H"))
648
649                          !                             IF (BASISTYPE .EQ. "NONORTHO") THEN
650
651                          !                                SMAT(IBRA, IKET) = SMAT(IBRA, IKET) + &
652                          !                                     BLOCH*ANGFACTOR * &
653                          !                                     UNIVSCALE(I, J, LBRA, LKET, MP, MAGR, "S")
654
655                          !                             ENDIF
656
657                       ENDDO
658
659                       !                         HK(IKET, IBRA, KCOUNT) = CONJG(HK(IBRA, IKET, KCOUNT))
660                       ! IF (BASISTYPE .EQ. "NONORTHO") &
661                       !      SMAT(IKET, IBRA) = SMAT(IBRA, IKET)
662
663                       IKET = IKET + 1
664
665                    ENDDO
666
667                 ENDDO
668
669                 IBRA = IBRA + 1
670
671              ENDDO
672           ENDDO
673        ENDIF
674        !              ENDIF
675     ENDDO
676     !              INDI = INDI + NORBI
677
678  ENDDO
679  !$OMP END PARALLEL DO
680  !        ENDDO
681  !     ENDDO
682  !  ENDDO
683
684  ! Save the diagonal elements: it will help a lot when we add in the partial charges
685
686  DO I = 1, NKTOT
687     DO J = 1,HDIM
688        HKDIAG(J,I) = HK(J,J,I)
689     ENDDO
690  ENDDO
691
692
693  !  IF (BASISTYPE .EQ. "NONORTHO") THEN
694
695  !     H0 = H
696  !     CALL GENX
697
698  !     IF (DEBUGON .EQ. 1) THEN
699
700  !        OPEN(UNIT=30, STATUS="UNKNOWN", FILE="myS.dat")
701  !        OPEN(UNIT=31, STATUS="UNKNOWN", FILE="myH0.dat")
702
703  !        PRINT*, "Caution - the Slater-Koster H and overlap matrices are being written to file"
704
705  !        DO I = 1, HDIM
706  !           WRITE(30,10) (SMAT(I,J), J = 1, HDIM)
707  !           WRITE(31,10) (H0(I,J), J = 1, HDIM)
708  !        ENDDO
709
710  !        CLOSE(30)
711  !        CLOSE(31)
712
713  !10      FORMAT(100F12.6)
714
715  !     ENDIF
716
717  !  ENDIF
718
719  RETURN
720
721END SUBROUTINE KBLDNEWH
722