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 BLDNEWHS
23
24  USE CONSTANTS_MOD
25  USE SETUPARRAY
26  USE NEBLISTARRAY
27  USE XBOARRAY
28  USE NONOARRAY
29  USE UNIVARRAY
30  USE MYPRECISION
31
32#ifdef PROGRESSON
33  USE GENXPROGRESS
34#endif
35
36  IMPLICIT NONE
37
38  INTEGER :: I, J, NEWJ, K, L, II, JJ, KK, MM, MP, NN, SUBI
39  INTEGER :: IBRA, IKET, LBRA, LKET, MBRA, MKET
40  INTEGER :: INDEX, INDI, INDJ
41  INTEGER :: SWITCH, PREVJ
42  INTEGER :: PBCI, PBCJ, PBCK
43  INTEGER :: BASISI(5), BASISJ(5), LBRAINC, LKETINC
44  REAL(LATTEPREC) :: ALPHA, BETA, COSBETA, PHI, TMP, PERM
45  REAL(LATTEPREC) :: RIJ(3), MAGR2, MAGR, MAGRP, RCUTTB
46  REAL(LATTEPREC) :: MAXRCUT, MAXRCUT2
47  REAL(LATTEPREC) :: ANGFACTOR, AMMBRA, WIGLBRAMBRA
48  !  REAL(LATTEPREC), ALLOCATABLE :: SAVEHDIAG(:)
49  REAL(LATTEPREC), EXTERNAL :: UNIVSCALE, WIGNERD, SLMMP, TLMMP, AM, BM
50  IF (EXISTERROR) RETURN
51
52
53  H = ZERO
54
55  IF (BASISTYPE .EQ. "NONORTHO") THEN
56
57     SMAT = ZERO
58     DO I = 1, HDIM
59        SMAT(I,I) = ONE
60     ENDDO
61
62  ENDIF
63
64  INDEX = 0
65
66  ! Build diagonal elements
67  DO I = 1, NATS
68
69     K = ELEMPOINTER(I)
70
71     SELECT CASE(BASIS(K))
72
73     CASE("s")
74
75        INDEX = INDEX + 1
76        H(INDEX, INDEX) = HES(K)
77
78     CASE("p")
79
80        DO SUBI = 1, 3
81           INDEX = INDEX + 1
82           H(INDEX,INDEX) = HEP(K)
83        ENDDO
84
85     CASE("d")
86
87        DO SUBI = 1, 5
88           INDEX = INDEX + 1
89           H(INDEX,INDEX) = HED(K)
90        ENDDO
91
92     CASE("f")
93
94        DO SUBI = 1, 7
95           INDEX = INDEX + 1
96           H(INDEX,INDEX) = HEF(K)
97        ENDDO
98
99     CASE("sp")
100
101        DO SUBI = 1, 4
102
103           INDEX = INDEX + 1
104           IF (SUBI .EQ. 1) THEN
105              H(INDEX,INDEX) = HES(K)
106           ELSE
107              H(INDEX,INDEX) = HEP(K)
108           ENDIF
109
110        ENDDO
111
112     CASE("sd")
113
114        DO SUBI = 1, 6
115
116           INDEX = INDEX + 1
117           IF (SUBI .EQ. 1) THEN
118              H(INDEX,INDEX) = HES(K)
119           ELSE
120              H(INDEX,INDEX) = HED(K)
121           ENDIF
122
123        ENDDO
124
125     CASE("sf")
126
127        DO SUBI = 1, 8
128
129           INDEX = INDEX + 1
130           IF (SUBI .EQ. 1) THEN
131              H(INDEX,INDEX) = HES(K)
132           ELSE
133              H(INDEX,INDEX) = HEF(K)
134           ENDIF
135
136        ENDDO
137
138     CASE("pd")
139
140        DO SUBI = 1, 8
141
142           INDEX = INDEX + 1
143           IF (SUBI .LE. 3) THEN
144              H(INDEX,INDEX) = HEP(K)
145           ELSE
146              H(INDEX,INDEX) = HED(K)
147           ENDIF
148
149        ENDDO
150
151     CASE("pf")
152
153        DO SUBI = 1, 10
154
155           INDEX = INDEX + 1
156           IF (SUBI .LE. 3) THEN
157              H(INDEX,INDEX) = HEP(K)
158           ELSE
159              H(INDEX,INDEX) = HEF(K)
160           ENDIF
161
162        ENDDO
163
164     CASE("df")
165
166        DO SUBI = 1, 12
167
168           INDEX = INDEX + 1
169           IF (SUBI .LE. 5) THEN
170              H(INDEX,INDEX) = HED(K)
171           ELSE
172              H(INDEX,INDEX) = HEF(K)
173           ENDIF
174
175        ENDDO
176
177     CASE("spd")
178
179        DO SUBI = 1, 9
180
181           INDEX = INDEX + 1
182           IF (SUBI .EQ. 1) THEN
183              H(INDEX,INDEX) = HES(K)
184           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN
185              H(INDEX,INDEX) = HEP(K)
186           ELSE
187              H(INDEX,INDEX) = HED(K)
188           ENDIF
189
190        ENDDO
191
192     CASE("spf")
193
194        DO SUBI = 1, 11
195
196           INDEX = INDEX + 1
197           IF (SUBI .EQ. 1) THEN
198              H(INDEX,INDEX) = HES(K)
199           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN
200              H(INDEX,INDEX) = HEP(K)
201           ELSE
202              H(INDEX,INDEX) = HEF(K)
203           ENDIF
204
205        ENDDO
206
207     CASE("sdf")
208
209        DO SUBI = 1, 13
210
211           INDEX = INDEX + 1
212           IF (SUBI .EQ. 1) THEN
213              H(INDEX,INDEX) = HES(K)
214           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 6) THEN
215              H(INDEX,INDEX) = HED(K)
216           ELSE
217              H(INDEX,INDEX) = HEF(K)
218           ENDIF
219
220        ENDDO
221
222
223     CASE("pdf")
224
225        DO SUBI = 1, 15
226
227           INDEX = INDEX + 1
228           IF (SUBI .LE. 3) THEN
229              H(INDEX,INDEX) = HEP(K)
230           ELSEIF (SUBI .GT. 3 .AND. SUBI .LE. 8) THEN
231              H(INDEX,INDEX) = HED(K)
232           ELSE
233              H(INDEX,INDEX) = HEF(K)
234           ENDIF
235
236        ENDDO
237
238     CASE("spdf")
239
240        DO SUBI = 1, 16
241
242           INDEX = INDEX + 1
243           IF (SUBI .EQ. 1) THEN
244              H(INDEX, INDEX) = HES(K)
245           ELSEIF (SUBI .GT. 1 .AND. SUBI .LE. 4) THEN
246              H(INDEX, INDEX) = HEP(K)
247           ELSEIF (SUBI .GT. 4 .AND. SUBI .LE. 9) THEN
248              H(INDEX, INDEX) = HED(K)
249           ELSE
250              H(INDEX, INDEX) = HEF(K)
251           ENDIF
252
253        ENDDO
254
255     END SELECT
256
257  ENDDO
258
259!$OMP PARALLEL DO DEFAULT (NONE) &
260!$OMP SHARED(NATS, BASIS, ELEMPOINTER, TOTNEBTB, NEBTB) &
261!$OMP SHARED(CR, BOX, H, SMAT, NOINT, ATELE, ELE1, ELE2) &
262!$OMP SHARED(HCUT, SCUT, MATINDLIST, BASISTYPE) &
263!$OMP PRIVATE(I, J, K, NEWJ, BASISI, BASISJ, INDI, INDJ, PBCI, PBCJ, PBCK) &
264!$OMP PRIVATE(RIJ, MAGR2, MAGR, MAGRP, PHI, ALPHA, BETA, COSBETA) &
265!$OMP PRIVATE(LBRAINC, LBRA, MBRA, L, LKETINC, LKET, MKET) &
266!$OMP PRIVATE(RCUTTB, IBRA, IKET, AMMBRA, WIGLBRAMBRA, ANGFACTOR, MP)
267
268  ! open loop over atoms I in system
269  DO I = 1, NATS
270
271     ! Build the lists of orbitals on each atom
272
273     SELECT CASE(BASIS(ELEMPOINTER(I)))
274
275     CASE("s")
276        BASISI(1) = 0
277        BASISI(2) = -1
278     CASE("p")
279        BASISI(1) = 1
280        BASISI(2) = -1
281     CASE("d")
282        BASISI(1) = 2
283        BASISI(2) = -1
284     CASE("f")
285        BASISI(1) = 3
286        BASISI(2) = -1
287     CASE("sp")
288        BASISI(1) = 0
289        BASISI(2) = 1
290        BASISI(3) = -1
291     CASE("sd")
292        BASISI(1) = 0
293        BASISI(2) = 2
294        BASISI(3) = -1
295     CASE("sf")
296        BASISI(1) = 0
297        BASISI(2) = 3
298        BASISI(3) = -1
299     CASE("pd")
300        BASISI(1) = 1
301        BASISI(2) = 2
302        BASISI(3) = -1
303     CASE("pf")
304        BASISI(1) = 1
305        BASISI(2) = 3
306        BASISI(3) = -1
307     CASE("df")
308        BASISI(1) = 2
309        BASISI(2) = 3
310        BASISI(3) = -1
311     CASE("spd")
312        BASISI(1) = 0
313        BASISI(2) = 1
314        BASISI(3) = 2
315        BASISI(4) = -1
316     CASE("spf")
317        BASISI(1) = 0
318        BASISI(2) = 1
319        BASISI(3) = 3
320        BASISI(4) = -1
321     CASE("sdf")
322        BASISI(1) = 0
323        BASISI(2) = 2
324        BASISI(3) = 3
325        BASISI(4) = -1
326     CASE("pdf")
327        BASISI(1) = 1
328        BASISI(2) = 2
329        BASISI(3) = 3
330        BASISI(4) = -1
331     CASE("spdf")
332        BASISI(1) = 0
333        BASISI(2) = 1
334        BASISI(3) = 2
335        BASISI(4) = 3
336        BASISI(5) = -1
337     END SELECT
338
339     INDI = MATINDLIST(I)
340
341     ! open loop over neighbors J of atom I
342     DO NEWJ = 1, TOTNEBTB(I)
343
344        J = NEBTB(1, NEWJ, I)
345
346        PBCI = NEBTB(2, NEWJ, I)
347        PBCJ = NEBTB(3, NEWJ, I)
348        PBCK = NEBTB(4, NEWJ, I)
349
350        RIJ(1) = CR(1,J) + REAL(PBCI)*BOX(1,1) + REAL(PBCJ)*BOX(2,1) + &
351             REAL(PBCK)*BOX(3,1) - CR(1,I)
352
353        RIJ(2) = CR(2,J) + REAL(PBCI)*BOX(1,2) + REAL(PBCJ)*BOX(2,2) + &
354             REAL(PBCK)*BOX(3,2) - CR(2,I)
355
356        RIJ(3) = CR(3,J) + REAL(PBCI)*BOX(1,3) + REAL(PBCJ)*BOX(2,3) + &
357             REAL(PBCK)*BOX(3,3) - CR(3,I)
358
359        MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3)
360
361        RCUTTB = ZERO
362        DO K = 1, NOINT
363
364           IF ( (ATELE(I) .EQ. ELE1(K) .AND. &
365                ATELE(J) .EQ. ELE2(K)) .OR. &
366                (ATELE(J) .EQ. ELE1(K) .AND. &
367                ATELE(I) .EQ. ELE2(K) )) THEN
368
369              IF (HCUT(K) .GT. RCUTTB ) RCUTTB = HCUT(K)
370
371              IF (BASISTYPE .EQ. "NONORTHO") THEN
372                 IF (SCUT(K) .GT. RCUTTB ) RCUTTB = SCUT(K)
373              ENDIF
374
375           ENDIF
376
377        ENDDO
378
379        IF (MAGR2 .LT. RCUTTB*RCUTTB) THEN
380
381           MAGR = SQRT(MAGR2)
382
383           SELECT CASE(BASIS(ELEMPOINTER(J)))
384           CASE("s")
385              BASISJ(1) = 0
386              BASISJ(2) = -1
387           CASE("p")
388              BASISJ(1) = 1
389              BASISJ(2) = -1
390           CASE("d")
391              BASISJ(1) = 2
392              BASISJ(2) = -1
393           CASE("f")
394              BASISJ(1) = 3
395              BASISJ(2) = -1
396           CASE("sp")
397              BASISJ(1) = 0
398              BASISJ(2) = 1
399              BASISJ(3) = -1
400           CASE("sd")
401              BASISJ(1) = 0
402              BASISJ(2) = 2
403              BASISJ(3) = -1
404           CASE("sf")
405              BASISJ(1) = 0
406              BASISJ(2) = 3
407              BASISJ(3) = -1
408           CASE("pd")
409              BASISJ(1) = 1
410              BASISJ(2) = 2
411              BASISJ(3) = -1
412           CASE("pf")
413              BASISJ(1) = 1
414              BASISJ(2) = 3
415              BASISJ(3) = -1
416           CASE("df")
417              BASISJ(1) = 2
418              BASISJ(2) = 3
419              BASISJ(3) = -1
420           CASE("spd")
421              BASISJ(1) = 0
422              BASISJ(2) = 1
423              BASISJ(3) = 2
424              BASISJ(4) = -1
425           CASE("spf")
426              BASISJ(1) = 0
427              BASISJ(2) = 1
428              BASISJ(3) = 3
429              BASISJ(4) = -1
430           CASE("sdf")
431              BASISJ(1) = 0
432              BASISJ(2) = 2
433              BASISJ(3) = 3
434              BASISJ(4) = -1
435           CASE("pdf")
436              BASISJ(1) = 1
437              BASISJ(2) = 2
438              BASISJ(3) = 3
439              BASISJ(4) = -1
440           CASE("spdf")
441              BASISJ(1) = 0
442              BASISJ(2) = 1
443              BASISJ(3) = 2
444              BASISJ(4) = 3
445              BASISJ(5) = -1
446           END SELECT
447
448           INDJ = MATINDLIST(J)
449
450           MAGRP = SQRT(RIJ(1) * RIJ(1) + RIJ(2) * RIJ(2))
451
452           ! transform to system in which z-axis is aligned with RIJ,
453           IF (ABS(RIJ(1)) .GT. 1.0E-12) THEN
454
455              IF (RIJ(1) .GT. ZERO .AND. RIJ(2) .GE. ZERO) THEN
456                 PHI = ZERO
457              ELSEIF (RIJ(1) .GT. ZERO .AND. RIJ(2) .LT. ZERO) THEN
458                 PHI = TWO * PI
459              ELSE
460                 PHI = PI
461              ENDIF
462              ALPHA = ATAN(RIJ(2) / RIJ(1)) + PHI
463
464           ELSEIF (ABS(RIJ(2)) .GT. 1.0E-12) THEN
465
466              IF (RIJ(2) .GT. 1.0E-12) THEN
467                 ALPHA = PI / TWO
468              ELSE
469                 ALPHA = THREE * PI / TWO
470              ENDIF
471
472           ELSE
473              ! pathological case: beta=0 and alpha undefined, but
474              ! this doesn't matter for matrix elements
475
476              ALPHA = ZERO
477
478           ENDIF
479
480           COSBETA = RIJ(3)/MAGR
481           BETA = ACOS(RIJ(3) / MAGR)
482
483           ! Build matrix elements using eqns (1)-(9) in PRB 72 165107
484
485           ! The loops over LBRA and LKET need to take into account
486           ! the orbitals assigned to each atom, e.g., sd rather than
487           ! spd...
488
489           IBRA = INDI + 1
490
491           LBRAINC = 1
492           DO WHILE (BASISI(LBRAINC) .NE. -1)
493
494              LBRA = BASISI(LBRAINC)
495              LBRAINC = LBRAINC + 1
496
497              DO MBRA = -LBRA, LBRA
498
499                 ! We can calculate these two outside the
500                 ! MKET loop...
501
502                 AMMBRA = AM(MBRA, ALPHA)
503                 WIGLBRAMBRA = WIGNERD(LBRA, ABS(MBRA), 0, COSBETA)
504
505                 IKET = INDJ + 1
506
507                 LKETINC = 1
508                 DO WHILE (BASISJ(LKETINC) .NE. -1)
509
510                    LKET = BASISJ(LKETINC)
511                    LKETINC = LKETINC + 1
512
513                    DO MKET = -LKET, LKET
514
515                       ! This is the sigma bonds (mp = 0)
516
517                       ! Hamiltonian build
518
519                       ! Pre-compute the angular part so we can use it
520                          ! again later if we're building the S matrix too
521
522                       ANGFACTOR = TWO * AMMBRA * &
523                            AM(MKET, ALPHA) * &
524                            WIGLBRAMBRA  * &
525                            WIGNERD(LKET, ABS(MKET), 0, COSBETA)
526
527                       H(IBRA, IKET) = H(IBRA, IKET) + ANGFACTOR * &
528                            UNIVSCALE(I, J, LBRA, LKET, 0, MAGR, "H")
529
530                       ! Overlap matrix build
531
532                       IF (BASISTYPE .EQ. "NONORTHO") THEN
533
534                          SMAT(IBRA, IKET) = SMAT(IBRA, IKET) + ANGFACTOR * &
535                               UNIVSCALE(I, J, LBRA, LKET, 0, MAGR, "S")
536
537                       ENDIF
538
539                       ! everything else
540
541                       DO MP = 1, MIN(LBRA, LKET)
542
543                          ANGFACTOR = SLMMP(LBRA, MBRA, MP, ALPHA, COSBETA) * &
544                               SLMMP(LKET, MKET, MP, ALPHA, COSBETA) + &
545                               TLMMP(LBRA, MBRA, MP, ALPHA, COSBETA) * &
546                               TLMMP(LKET, MKET, MP, ALPHA, COSBETA)
547
548                          H(IBRA, IKET) = H(IBRA, IKET) + ANGFACTOR * &
549                               UNIVSCALE(I, J, LBRA, LKET, MP, MAGR, "H")
550
551                          IF (BASISTYPE .EQ. "NONORTHO") THEN
552
553                             SMAT(IBRA, IKET) = SMAT(IBRA, IKET) + &
554                                  ANGFACTOR * &
555                                  UNIVSCALE(I, J, LBRA, LKET, MP, MAGR, "S")
556
557                          ENDIF
558
559                       ENDDO
560
561
562                       IKET = IKET + 1
563
564                    ENDDO
565
566                 ENDDO
567
568                 IBRA = IBRA + 1
569
570              ENDDO
571           ENDDO
572        ENDIF
573     ENDDO
574
575  ENDDO
576
577!$OMP END PARALLEL DO
578
579  DO I = 1, HDIM
580     HDIAG(I) = H(I,I)
581  ENDDO
582
583
584  IF (BASISTYPE .EQ. "NONORTHO") THEN
585
586     H0 = H
587
588#ifdef PROGRESSON
589     IF(LATTEINEXISTS)THEN
590        CALL GENXBML
591     ELSE
592        CALL GENX
593     ENDIF
594#else
595     CALL GENX
596#endif
597
598
599     IF (DEBUGON .EQ. 1) THEN
600
601        OPEN(UNIT=30, STATUS="UNKNOWN", FILE="myS.dat")
602        OPEN(UNIT=31, STATUS="UNKNOWN", FILE="myH0.dat")
603
604        PRINT*, "Caution - the Slater-Koster H and overlap matrices are being written to file"
605
606        DO I = 1, HDIM
607           WRITE(30,10) (SMAT(I,J), J = 1, HDIM)
608           WRITE(31,10) (H0(I,J), J = 1, HDIM)
609        ENDDO
610
611        CLOSE(30)
612        CLOSE(31)
613
61410      FORMAT(100F12.6)
615
616     ENDIF
617
618  ENDIF
619
620  RETURN
621
622END SUBROUTINE BLDNEWHS
623