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  REAL(LATTEPREC) :: B1(3), B2(3), B3(3), A1A2XA3, K0(3)
51  COMPLEX(LATTEPREC) :: BLOCH, KHTMP, KSTMP
52  REAL(LATTEPREC), EXTERNAL :: UNIVSCALE, WIGNERD, SLMMP, TLMMP, AM, BM
53  IF (EXISTERROR) RETURN
54
55  HK = CMPLX(ZERO)
56
57  ! Computing the reciprocal lattice vectors
58
59  B1(1) = BOX(2,2)*BOX(3,3) - BOX(3,2)*BOX(2,3)
60  B1(2) = BOX(3,1)*BOX(2,3) - BOX(2,1)*BOX(3,3)
61  B1(3) = BOX(2,1)*BOX(3,2) - BOX(3,1)*BOX(2,2)
62
63  A1A2XA3 = BOX(1,1)*B1(1) + BOX(1,2)*B1(2) + BOX(1,3)*B1(3)
64
65  ! B1 = (A2 X A3)/(A1.(A2 X A3))
66
67  B1 = B1/A1A2XA3
68
69  ! B2 = (A3 x A1)/(A1(A2 X A3))
70
71  B2(1) = (BOX(3,2)*BOX(1,3) - BOX(1,2)*BOX(3,3))/A1A2XA3
72  B2(2) = (BOX(1,1)*BOX(3,3) - BOX(3,1)*BOX(1,3))/A1A2XA3
73  B2(3) = (BOX(3,1)*BOX(1,2) - BOX(1,1)*BOX(3,2))/A1A2XA3
74
75  ! B3 = (A1 x A2)/(A1(A2 X A3))
76
77  B3(1) = (BOX(1,2)*BOX(2,3) - BOX(2,2)*BOX(1,3))/A1A2XA3
78  B3(2) = (BOX(2,1)*BOX(1,3) - BOX(1,1)*BOX(2,3))/A1A2XA3
79  B3(3) = (BOX(1,1)*BOX(2,2) - BOX(2,1)*BOX(1,2))/A1A2XA3
80
81  INDEX = 0
82
83  ! Build diagonal elements
84  DO I = 1, NATS
85
86     SELECT CASE(BASIS(ELEMPOINTER(I)))
87
88     CASE("s")
89
90        INDEX = INDEX + 1
91        HK(INDEX, INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
92
93     CASE("p")
94
95        DO SUBI = 1, 3
96           INDEX = INDEX + 1
97           HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
98        ENDDO
99
100     CASE("d")
101
102        DO SUBI = 1, 5
103           INDEX = INDEX + 1
104           HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
105        ENDDO
106
107     CASE("f")
108
109        DO SUBI = 1, 7
110           INDEX = INDEX + 1
111           HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
112        ENDDO
113
114     CASE("sp")
115
116        DO SUBI = 1, 4
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(HEP(ELEMPOINTER(I)))
123           ENDIF
124
125        ENDDO
126
127     CASE("sd")
128
129        DO SUBI = 1, 6
130
131           INDEX = INDEX + 1
132           IF (SUBI .EQ. 1) THEN
133              HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
134           ELSE
135              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
136           ENDIF
137
138        ENDDO
139
140     CASE("sf")
141
142        DO SUBI = 1, 8
143
144           INDEX = INDEX + 1
145           IF (SUBI .EQ. 1) THEN
146              HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
147           ELSE
148              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
149           ENDIF
150
151        ENDDO
152
153     CASE("pd")
154
155        DO SUBI = 1, 8
156
157           INDEX = INDEX + 1
158           IF (SUBI .LE. 3) THEN
159              HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
160           ELSE
161              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
162           ENDIF
163
164        ENDDO
165
166     CASE("pf")
167
168        DO SUBI = 1, 10
169
170           INDEX = INDEX + 1
171           IF (SUBI .LE. 3) THEN
172              HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
173           ELSE
174              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
175           ENDIF
176
177        ENDDO
178
179     CASE("df")
180
181        DO SUBI = 1, 12
182
183           INDEX = INDEX + 1
184           IF (SUBI .LE. 5) THEN
185              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
186           ELSE
187              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
188           ENDIF
189
190        ENDDO
191
192     CASE("spd")
193
194        DO SUBI = 1, 9
195
196           INDEX = INDEX + 1
197           IF (SUBI .EQ. 1) THEN
198              HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
199           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN
200              HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
201           ELSE
202              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
203           ENDIF
204
205        ENDDO
206
207     CASE("spf")
208
209        DO SUBI = 1, 11
210
211           INDEX = INDEX + 1
212           IF (SUBI .EQ. 1) THEN
213              HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
214           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN
215              HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
216           ELSE
217              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
218           ENDIF
219
220        ENDDO
221
222     CASE("sdf")
223
224        DO SUBI = 1, 13
225
226           INDEX = INDEX + 1
227           IF (SUBI .EQ. 1) THEN
228              HK(INDEX,INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
229           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 6) THEN
230              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
231           ELSE
232              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
233           ENDIF
234
235        ENDDO
236
237
238     CASE("pdf")
239
240        DO SUBI = 1, 15
241
242           INDEX = INDEX + 1
243           IF (SUBI .LE. 3) THEN
244              HK(INDEX,INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
245           ELSEIF (SUBI .GT. 3 .AND. SUBI .LE. 8) THEN
246              HK(INDEX,INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
247           ELSE
248              HK(INDEX,INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
249           ENDIF
250
251        ENDDO
252
253     CASE("spdf")
254
255        DO SUBI = 1, 16
256
257           INDEX = INDEX + 1
258           IF (SUBI .EQ. 1) THEN
259              HK(INDEX, INDEX, 1) = CMPLX(HES(ELEMPOINTER(I)))
260           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN
261              HK(INDEX, INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
262           ELSEIF (SUBI .GT.  1 .AND. SUBI .LE. 4) THEN
263              HK(INDEX, INDEX, 1) = CMPLX(HEP(ELEMPOINTER(I)))
264           ELSEIF (SUBI .GT. 4 .AND. SUBI .LE. 9) THEN
265              HK(INDEX, INDEX, 1) = CMPLX(HED(ELEMPOINTER(I)))
266           ELSE
267              HK(INDEX, INDEX, 1) = CMPLX(HEF(ELEMPOINTER(I)))
268           ENDIF
269
270        ENDDO
271
272     END SELECT
273
274  ENDDO
275
276  DO I = 2, NKTOT
277     DO J  = 1, HDIM
278        HK(J,J,I) = HK(J,J,1)
279     ENDDO
280  ENDDO
281
282  ! We assign the diagonal elements in ADDQDEP
283
284
285  IF (BASISTYPE .EQ. "NONORTHO") THEN
286
287     SK = CMPLX(ZERO)
288
289     DO I = 1, NKTOT
290        DO J = 1, HDIM
291           SK(J,J,I) = CMPLX(ONE,ZERO)
292        ENDDO
293     ENDDO
294
295  ENDIF
296
297  K0 = PI*(ONE - REAL(NKX))/(REAL(NKX))*B1 + &
298       PI*(ONE - REAL(NKY))/(REAL(NKY))*B2 + &
299       PI*(ONE - REAL(NKZ))/(REAL(NKZ))*B3 - PI*KSHIFT
300
301!$OMP PARALLEL DO DEFAULT (NONE) &
302!$OMP SHARED(NATS, BASIS, ELEMPOINTER, TOTNEBTB, NEBTB) &
303!$OMP SHARED(CR, BOX, B1, B2, B3, NOINT, ATELE, ELE1, ELE2, HK, SK) &
304!$OMP SHARED(HCUT, SCUT, MATINDLIST, BASISTYPE, K0, NKX, NKY, NKZ) &
305!$OMP PRIVATE(I, J, K, NEWJ, BASISI, BASISJ, INDI, INDJ, PBCI, PBCJ, PBCK) &
306!$OMP PRIVATE(RIJ, MAGR2, MAGR, MAGRP, PHI, ALPHA, BETA, COSBETA) &
307!$OMP PRIVATE(LBRAINC, LBRA, MBRA, L, LKETINC, LKET, MKET) &
308!$OMP PRIVATE(BLOCH, KDOTL, KPOINT, KCOUNT, KHTMP, KSTMP)&
309!$OMP PRIVATE(RCUTTB, IBRA, IKET, AMMBRA, WIGLBRAMBRA, ANGFACTOR, MP)
310!!$OMP REDUCTION(+:HK)
311
312  DO I = 1, NATS
313
314     ! Build the lists of orbitals on each atom
315
316     SELECT CASE(BASIS(ELEMPOINTER(I)))
317
318     CASE("s")
319        BASISI(1) = 0
320        BASISI(2) = -1
321     CASE("p")
322        BASISI(1) = 1
323        BASISI(2) = -1
324     CASE("d")
325        BASISI(1) = 2
326        BASISI(2) = -1
327     CASE("f")
328        BASISI(1) = 3
329        BASISI(2) = -1
330     CASE("sp")
331        BASISI(1) = 0
332        BASISI(2) = 1
333        BASISI(3) = -1
334     CASE("sd")
335        BASISI(1) = 0
336        BASISI(2) = 2
337        BASISI(3) = -1
338     CASE("sf")
339        BASISI(1) = 0
340        BASISI(2) = 3
341        BASISI(3) = -1
342     CASE("pd")
343        BASISI(1) = 1
344        BASISI(2) = 2
345        BASISI(3) = -1
346     CASE("pf")
347        BASISI(1) = 1
348        BASISI(2) = 3
349        BASISI(3) = -1
350     CASE("df")
351        BASISI(1) = 2
352        BASISI(2) = 3
353        BASISI(3) = -1
354     CASE("spd")
355        BASISI(1) = 0
356        BASISI(2) = 1
357        BASISI(3) = 2
358        BASISI(4) = -1
359     CASE("spf")
360        BASISI(1) = 0
361        BASISI(2) = 1
362        BASISI(3) = 3
363        BASISI(4) = -1
364     CASE("sdf")
365        BASISI(1) = 0
366        BASISI(2) = 2
367        BASISI(3) = 3
368        BASISI(4) = -1
369     CASE("pdf")
370        BASISI(1) = 1
371        BASISI(2) = 2
372        BASISI(3) = 3
373        BASISI(4) = -1
374     CASE("spdf")
375        BASISI(1) = 0
376        BASISI(2) = 1
377        BASISI(3) = 2
378        BASISI(4) = 3
379        BASISI(5) = -1
380     END SELECT
381
382     INDI = MATINDLIST(I)
383
384     ! open loop over neighbors J of atom I
385     DO NEWJ = 1, TOTNEBTB(I)
386
387        J = NEBTB(1, NEWJ, I)
388
389        PBCI = NEBTB(2, NEWJ, I)
390        PBCJ = NEBTB(3, NEWJ, I)
391        PBCK = NEBTB(4, NEWJ, I)
392
393        RIJ(1) = CR(1,J) + REAL(PBCI)*BOX(1,1) + REAL(PBCJ)*BOX(2,1) + &
394             REAL(PBCK)*BOX(3,1) - CR(1,I)
395
396        RIJ(2) = CR(2,J) + REAL(PBCI)*BOX(1,2) + REAL(PBCJ)*BOX(2,2) + &
397             REAL(PBCK)*BOX(3,2) - CR(2,I)
398
399        RIJ(3) = CR(3,J) + REAL(PBCI)*BOX(1,3) + REAL(PBCJ)*BOX(2,3) + &
400             REAL(PBCK)*BOX(3,3) - CR(3,I)
401
402        MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3)
403
404        RCUTTB = ZERO
405
406        DO K = 1, NOINT
407
408           IF ( (ATELE(I) .EQ. ELE1(K) .AND. &
409                ATELE(J) .EQ. ELE2(K)) .OR. &
410                (ATELE(J) .EQ. ELE1(K) .AND. &
411                ATELE(I) .EQ. ELE2(K) )) THEN
412
413              IF (HCUT(K) .GT. RCUTTB ) RCUTTB = HCUT(K)
414
415              IF (BASISTYPE .EQ. "NONORTHO") THEN
416                 IF (SCUT(K) .GT. RCUTTB ) RCUTTB = SCUT(K)
417              ENDIF
418
419           ENDIF
420
421        ENDDO
422
423        IF (MAGR2 .LT. RCUTTB*RCUTTB) THEN
424
425           MAGR = SQRT(MAGR2)
426
427           SELECT CASE(BASIS(ELEMPOINTER(J)))
428           CASE("s")
429              BASISJ(1) = 0
430              BASISJ(2) = -1
431           CASE("p")
432              BASISJ(1) = 1
433              BASISJ(2) = -1
434           CASE("d")
435              BASISJ(1) = 2
436              BASISJ(2) = -1
437           CASE("f")
438              BASISJ(1) = 3
439              BASISJ(2) = -1
440           CASE("sp")
441              BASISJ(1) = 0
442              BASISJ(2) = 1
443              BASISJ(3) = -1
444           CASE("sd")
445              BASISJ(1) = 0
446              BASISJ(2) = 2
447              BASISJ(3) = -1
448           CASE("sf")
449              BASISJ(1) = 0
450              BASISJ(2) = 3
451              BASISJ(3) = -1
452           CASE("pd")
453              BASISJ(1) = 1
454              BASISJ(2) = 2
455              BASISJ(3) = -1
456           CASE("pf")
457              BASISJ(1) = 1
458              BASISJ(2) = 3
459              BASISJ(3) = -1
460           CASE("df")
461              BASISJ(1) = 2
462              BASISJ(2) = 3
463              BASISJ(3) = -1
464           CASE("spd")
465              BASISJ(1) = 0
466              BASISJ(2) = 1
467              BASISJ(3) = 2
468              BASISJ(4) = -1
469           CASE("spf")
470              BASISJ(1) = 0
471              BASISJ(2) = 1
472              BASISJ(3) = 3
473              BASISJ(4) = -1
474           CASE("sdf")
475              BASISJ(1) = 0
476              BASISJ(2) = 2
477              BASISJ(3) = 3
478              BASISJ(4) = -1
479           CASE("pdf")
480              BASISJ(1) = 1
481              BASISJ(2) = 2
482              BASISJ(3) = 3
483              BASISJ(4) = -1
484           CASE("spdf")
485              BASISJ(1) = 0
486              BASISJ(2) = 1
487              BASISJ(3) = 2
488              BASISJ(4) = 3
489              BASISJ(5) = -1
490           END SELECT
491
492           INDJ = MATINDLIST(J)
493
494           MAGRP = SQRT(RIJ(1) * RIJ(1) + RIJ(2) * RIJ(2))
495
496           ! transform to system in which z-axis is aligned with RIJ,
497           IF (ABS(RIJ(1)) .GT. 1.0E-12) THEN
498
499              IF (RIJ(1) .GT. ZERO .AND. RIJ(2) .GE. ZERO) THEN
500                 PHI = ZERO
501              ELSEIF (RIJ(1) .GT. ZERO .AND. RIJ(2) .LT. ZERO) THEN
502                 PHI = TWO * PI
503              ELSE
504                 PHI = PI
505              ENDIF
506              ALPHA = ATAN(RIJ(2) / RIJ(1)) + PHI
507
508           ELSEIF (ABS(RIJ(2)) .GT. 1.0E-12) THEN
509
510              IF (RIJ(2) .GT. 1.0E-12) THEN
511                 ALPHA = PI / TWO
512              ELSE
513                 ALPHA = THREE * PI / TWO
514              ENDIF
515
516           ELSE
517
518              ! pathological case: beta=0 and alpha undefined, but
519              ! this doesn't matter for matrix elements
520
521              ALPHA = ZERO
522
523           ENDIF
524
525           COSBETA = RIJ(3)/MAGR
526           BETA = ACOS( RIJ(3) / MAGR )
527
528           ! Build matrix elements using eqns (1)-(9) in PRB 72 165107
529
530           ! The loops over LBRA and LKET need to take into account
531           ! the orbitals assigned to each atom, e.g., sd rather than
532           ! spd...
533
534           IBRA = INDI + 1
535
536           LBRAINC = 1
537           DO WHILE (BASISI(LBRAINC) .NE. -1)
538
539              LBRA = BASISI(LBRAINC)
540              LBRAINC = LBRAINC + 1
541
542              DO MBRA = -LBRA, LBRA
543
544                 ! We can calculate these two outside the
545                 ! MKET loop...
546
547                 AMMBRA = AM(MBRA, ALPHA)
548                 WIGLBRAMBRA = WIGNERD(LBRA, ABS(MBRA), 0, COSBETA)
549
550                 IKET = INDJ + 1
551
552                 LKETINC = 1
553                 DO WHILE (BASISJ(LKETINC) .NE. -1)
554
555                    LKET = BASISJ(LKETINC)
556                    LKETINC = LKETINC + 1
557
558                    DO MKET = -LKET, LKET
559
560                       ! This is the sigma bonds (mp = 0)
561
562                       ! Hamiltonian build
563
564                       ! Pre-compute the angular part so we can use it
565                       ! again later if we're building the S matrix too
566
567                       ANGFACTOR = TWO * AMMBRA * &
568                            AM(MKET, ALPHA) * &
569                            WIGLBRAMBRA * &
570                            WIGNERD(LKET, ABS(MKET), 0, COSBETA)
571
572                       KHTMP = CMPLX(ANGFACTOR * &
573                            UNIVSCALE(I, J, LBRA, LKET, &
574                            0, MAGR, "H"))
575
576                       IF (BASISTYPE .EQ. "NONORTHO") THEN
577                          KSTMP = CMPLX(ANGFACTOR * &
578                               UNIVSCALE(I, J, LBRA, LKET, &
579                               0, MAGR, "S"))
580                       ENDIF
581
582                       KCOUNT = 0
583
584                       !Loop over k-points
585
586                       DO KX = 1, NKX
587
588                          DO KY = 1, NKY
589
590                             DO KZ = 1, NKZ
591
592                                KPOINT = TWO*PI*(REAL(KX-1)*B1/REAL(NKX) + &
593                                     REAL(KY-1)*B2/REAL(NKY) + &
594                                     REAL(KZ-1)*B3/REAL(NKZ)) + K0
595
596                                KCOUNT = KCOUNT+1
597
598                                KDOTL = KPOINT(1)*RIJ(1) + KPOINT(2)*RIJ(2) + &
599                                     KPOINT(3)*RIJ(3)
600
601                                BLOCH = EXP(CMPLX(ZERO,KDOTL))
602
603                                HK(IBRA, IKET, KCOUNT) = &
604                                     HK(IBRA, IKET, KCOUNT) + &
605                                     BLOCH*KHTMP
606
607                                IF (BASISTYPE .EQ. "NONORTHO") THEN
608                                   SK(IBRA, IKET, KCOUNT) = &
609                                        SK(IBRA, IKET, KCOUNT) + &
610                                        BLOCH*KSTMP
611                                ENDIF
612
613                             ENDDO
614                          ENDDO
615                       ENDDO
616
617
618                       ! everything else
619
620                       DO MP = 1, MIN(LBRA, LKET)
621
622                          ANGFACTOR = &
623                               SLMMP(LBRA, MBRA, MP, ALPHA, COSBETA)* &
624                               SLMMP(LKET, MKET, MP, ALPHA, COSBETA)+ &
625                               TLMMP(LBRA, MBRA, MP, ALPHA, COSBETA)* &
626                               TLMMP(LKET, MKET, MP, ALPHA, COSBETA)
627
628                          KHTMP = CMPLX(ANGFACTOR * &
629                               UNIVSCALE(I, J, LBRA, LKET, &
630                               MP, MAGR, "H"))
631
632                          IF (BASISTYPE .EQ. "NONORTHO") THEN
633                             KSTMP = CMPLX(ANGFACTOR * &
634                                  UNIVSCALE(I, J, LBRA, LKET, &
635                                  MP, MAGR, "S"))
636                          ENDIF
637
638                          KCOUNT = 0
639
640                          DO KX = 1, NKX
641                             DO KY = 1, NKY
642                                DO KZ = 1, NKZ
643
644                                   KPOINT = TWO*PI*(REAL(KX-1)*B1/REAL(NKX) + &
645                                        REAL(KY-1)*B2/REAL(NKY) + &
646                                        REAL(KZ-1)*B3/REAL(NKZ)) + K0
647
648                                   KCOUNT = KCOUNT+1
649
650                                   KDOTL = KPOINT(1)*RIJ(1) + &
651                                        KPOINT(2)*RIJ(2) + KPOINT(3)*RIJ(3)
652
653                                   BLOCH = EXP(CMPLX(ZERO,KDOTL))
654
655                                   HK(IBRA, IKET, KCOUNT) = &
656                                        HK(IBRA, IKET, KCOUNT) + &
657                                        BLOCH*KHTMP
658
659                                   IF (BASISTYPE .EQ. "NONORTHO") THEN
660                                      SK(IBRA, IKET, KCOUNT) = &
661                                           SK(IBRA, IKET, KCOUNT) + &
662                                           BLOCH*KSTMP
663                                   ENDIF
664
665                                ENDDO
666                             ENDDO
667                          ENDDO
668
669                       ENDDO
670
671                       IKET = IKET + 1
672
673                    ENDDO
674
675                 ENDDO
676
677                 IBRA = IBRA + 1
678
679              ENDDO
680           ENDDO
681        ENDIF
682     ENDDO
683  ENDDO
684
685!$OMP END PARALLEL DO
686
687  ! Save the diagonal elements: it will help a lot when we add in the partial charges
688
689  IF (BASISTYPE .EQ. "ORTHO") THEN
690
691     DO I = 1, NKTOT
692        DO J = 1,HDIM
693           HKDIAG(J,I) = HK(J,J,I)
694        ENDDO
695     ENDDO
696
697  ELSEIF (BASISTYPE .EQ. "NONORTHO") THEN
698
699     ! keep a copy of the charge-independent Hamiltonian
700
701     HK0 = HK
702
703     CALL KGENX ! Get the factors for the Lowdin transform
704
705  ENDIF
706
707
708  IF (DEBUGON .EQ. 1 ) THEN
709
710     OPEN(UNIT=31, STATUS="UNKNOWN", FILE="myS0_k_R.dat")
711
712     OPEN(UNIT=32, STATUS="UNKNOWN", FILE="myS0_k_I.dat")
713
714     DO K = 1, NKTOT
715        WRITE(31,*) K
716        WRITE(32,*) K
717        DO I = 1, HDIM
718           WRITE(31,10) (REAL(SK(I,J,K)), J = 1, HDIM)
719           WRITE(32,10) (AIMAG(SK(I,J,K)), J = 1, HDIM)
720        ENDDO
721     ENDDO
722
723     CLOSE(31)
724     CLOSE(32)
725
72610   FORMAT(648F12.6)
727
728  ENDIF
729
730  RETURN
731
732END SUBROUTINE KBLDNEWH
733