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 KFCOULNONO
23
24  USE CONSTANTS_MOD
25  USE SETUPARRAY
26  USE UNIVARRAY
27  USE NONOARRAY
28  USE COULOMBARRAY
29  USE NEBLISTARRAY
30  USE SPINARRAY
31  USE VIRIALARRAY
32  USE KSPACEARRAY
33  USE MYPRECISION
34
35  IMPLICIT NONE
36
37  INTEGER :: I, J, K, L, M, N, KK, INDI, INDJ
38  INTEGER :: LBRA, MBRA, LKET, MKET
39  INTEGER :: PREVJ, NEWJ
40  INTEGER :: PBCI, PBCJ, PBCK
41  INTEGER :: BASISI(5), BASISJ(5), LBRAINC, LKETINC
42  INTEGER :: KX, KY, KZ, KCOUNT
43  REAL(LATTEPREC) :: ALPHA, BETA, PHI, COSBETA
44  REAL(LATTEPREC) :: RIJ(3), DC(3)
45  REAL(LATTEPREC) :: MAGR, MAGR2, MAGRP, MAGRP2
46  REAL(LATTEPREC) :: MAXRCUT, MAXRCUT2
47  REAL(LATTEPREC) :: MYDFDA, MYDFDB, MYDFDR, RCUTTB
48  REAL(LATTEPREC) :: KPOINT(3), KX0, KY0, KZ0, KDOTL
49  REAL(LATTEPREC) :: B1(3), B2(3), B3(3), MAG1, MAG2, MAG3, A1A2XA3, K0(3)
50  REAL(LATTEPREC), EXTERNAL :: DFDA, DFDB, DFDR
51  COMPLEX(LATTEPREC) :: FTMP(3), RHO, CONJGBLOCH
52  LOGICAL PATH
53  IF (EXISTERROR) RETURN
54
55  ! These were allocated elsewhere. We'll use them to accumulate the complex forces
56
57  IF (SPINON .EQ. 1) THEN
58     CALL ERRORS("kfcoulnono","Non-ortho k-space and spin polarization not yet implemented")
59  ENDIF
60
61  KF = CMPLX(ZERO)
62  VIRBONDK = CMPLX(ZERO)
63
64  ! Computing the reciprocal lattice vectors
65
66  B1(1) = BOX(2,2)*BOX(3,3) - BOX(3,2)*BOX(2,3)
67  B1(2) = BOX(3,1)*BOX(2,3) - BOX(2,1)*BOX(3,3)
68  B1(3) = BOX(2,1)*BOX(3,2) - BOX(3,1)*BOX(2,2)
69
70  A1A2XA3 = BOX(1,1)*B1(1) + BOX(1,2)*B1(2) + BOX(1,3)*B1(3)
71
72  ! B1 = 2*PI*(A2 X A3)/(A1.(A2 X A3))
73
74  B1 = B1/A1A2XA3
75
76  ! B2 = 2*PI*(A3 x A1)/(A1(A2 X A3))
77
78  B2(1) = (BOX(3,2)*BOX(1,3) - BOX(1,2)*BOX(3,3))/A1A2XA3
79  B2(2) = (BOX(1,1)*BOX(3,3) - BOX(3,1)*BOX(1,3))/A1A2XA3
80  B2(3) = (BOX(3,1)*BOX(1,2) - BOX(1,1)*BOX(3,2))/A1A2XA3
81
82  ! B3 = 2*PI*(A1 x A2)/(A1(A2 X A3))
83
84  B3(1) = (BOX(1,2)*BOX(2,3) - BOX(2,2)*BOX(1,3))/A1A2XA3
85  B3(2) = (BOX(2,1)*BOX(1,3) - BOX(1,1)*BOX(2,3))/A1A2XA3
86  B3(3) = (BOX(1,1)*BOX(2,2) - BOX(2,1)*BOX(1,2))/A1A2XA3
87
88  K0 = PI*(ONE - REAL(NKX))/(REAL(NKX))*B1 + &
89       PI*(ONE - REAL(NKY))/(REAL(NKY))*B2 + &
90       PI*(ONE - REAL(NKZ))/(REAL(NKZ))*B3 - PI*KSHIFT
91
92
93!$OMP PARALLEL DO DEFAULT (NONE) &
94!$OMP SHARED(NATS, BASIS, ELEMPOINTER, TOTNEBTB, NEBTB) &
95!$OMP SHARED(CR, BOX, KBO, RHOUP, RHODOWN, SPINON, NOINT, ATELE, ELE1, ELE2) &
96!$OMP SHARED(HCUT, SCUT, MATINDLIST, BASISTYPE) &
97!$OMP SHARED(HUBBARDU, DELTAQ, COULOMBV) &
98!$OMP SHARED(K0, B1, B2, B3, NKX, NKY, NKZ, KF) &
99!$OMP PRIVATE(I, J, K, NEWJ, BASISI, BASISJ, INDI, INDJ, PBCI, PBCJ, PBCK) &
100!$OMP PRIVATE(RIJ, MAGR2, MAGR, MAGRP2, MAGRP, PATH, PHI, ALPHA, BETA, COSBETA, FTMP) &
101!$OMP PRIVATE(DC, LBRAINC, LBRA, MBRA, L, LKETINC, LKET, MKET, RHO) &
102!$OMP PRIVATE(MYDFDA, MYDFDB, MYDFDR, RCUTTB, CONJGBLOCH, KDOTL) &
103!$OMP PRIVATE(KPOINT, KCOUNT) &
104!$OMP REDUCTION(+: VIRBONDK)
105
106
107  DO I = 1, NATS
108
109     ! Build list of orbitals on atom I
110     SELECT CASE(BASIS(ELEMPOINTER(I)))
111
112     CASE("s")
113        BASISI(1) = 0
114        BASISI(2) = -1
115     CASE("p")
116        BASISI(1) = 1
117        BASISI(2) = -1
118     CASE("d")
119        BASISI(1) = 2
120        BASISI(2) = -1
121     CASE("f")
122        BASISI(1) = 3
123        BASISI(2) = -1
124     CASE("sp")
125        BASISI(1) = 0
126        BASISI(2) = 1
127        BASISI(3) = -1
128     CASE("sd")
129        BASISI(1) = 0
130        BASISI(2) = 2
131        BASISI(3) = -1
132     CASE("sf")
133        BASISI(1) = 0
134        BASISI(2) = 3
135        BASISI(3) = -1
136     CASE("pd")
137        BASISI(1) = 1
138        BASISI(2) = 2
139        BASISI(3) = -1
140     CASE("pf")
141        BASISI(1) = 1
142        BASISI(2) = 3
143        BASISI(3) = -1
144     CASE("df")
145        BASISI(1) = 2
146        BASISI(2) = 3
147        BASISI(3) = -1
148     CASE("spd")
149        BASISI(1) = 0
150        BASISI(2) = 1
151        BASISI(3) = 2
152        BASISI(4) = -1
153     CASE("spf")
154        BASISI(1) = 0
155        BASISI(2) = 1
156        BASISI(3) = 3
157        BASISI(4) = -1
158     CASE("sdf")
159        BASISI(1) = 0
160        BASISI(2) = 2
161        BASISI(3) = 3
162        BASISI(4) = -1
163     CASE("pdf")
164        BASISI(1) = 1
165        BASISI(2) = 2
166        BASISI(3) = 3
167        BASISI(4) = -1
168     CASE("spdf")
169        BASISI(1) = 0
170        BASISI(2) = 1
171        BASISI(3) = 2
172        BASISI(4) = 3
173        BASISI(5) = -1
174     END SELECT
175
176     INDI = MATINDLIST(I)
177
178     DO NEWJ = 1, TOTNEBTB(I)
179
180        J = NEBTB(1, NEWJ, I)
181        PBCI = NEBTB(2, NEWJ, I)
182        PBCJ = NEBTB(3, NEWJ, I)
183        PBCK = NEBTB(4, NEWJ, I)
184
185        RIJ(1) = CR(1,J) + REAL(PBCI)*BOX(1,1) + REAL(PBCJ)*BOX(2,1) + &
186             REAL(PBCK)*BOX(3,1) - CR(1,I)
187
188        RIJ(2) = CR(2,J) + REAL(PBCI)*BOX(1,2) + REAL(PBCJ)*BOX(2,2) + &
189             REAL(PBCK)*BOX(3,2) - CR(2,I)
190
191        RIJ(3) = CR(3,J) + REAL(PBCI)*BOX(1,3) + REAL(PBCJ)*BOX(2,3) + &
192             REAL(PBCK)*BOX(3,3) - CR(3,I)
193
194        MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3)
195
196        RCUTTB = ZERO
197        DO K = 1, NOINT
198
199           IF ( (ATELE(I) .EQ. ELE1(K) .AND. &
200                ATELE(J) .EQ. ELE2(K)) .OR. &
201                (ATELE(J) .EQ. ELE1(K) .AND. &
202                ATELE(I) .EQ. ELE2(K) )) THEN
203
204              IF (SCUT(K) .GT. RCUTTB ) RCUTTB = HCUT(K)
205
206              IF (BASISTYPE .EQ. "NONORTHO") THEN
207                 IF (SCUT(K) .GT. RCUTTB ) RCUTTB = SCUT(K)
208              ENDIF
209
210           ENDIF
211
212        ENDDO
213
214        IF (MAGR2 .LT. RCUTTB*RCUTTB) THEN
215
216           MAGR = SQRT(MAGR2)
217
218           ! Build list of orbitals on atom J
219
220           SELECT CASE(BASIS(ELEMPOINTER(J)))
221           CASE("s")
222              BASISJ(1) = 0
223              BASISJ(2) = -1
224           CASE("p")
225              BASISJ(1) = 1
226              BASISJ(2) = -1
227           CASE("d")
228              BASISJ(1) = 2
229              BASISJ(2) = -1
230           CASE("f")
231              BASISJ(1) = 3
232              BASISJ(2) = -1
233           CASE("sp")
234              BASISJ(1) = 0
235              BASISJ(2) = 1
236              BASISJ(3) = -1
237           CASE("sd")
238              BASISJ(1) = 0
239              BASISJ(2) = 2
240              BASISJ(3) = -1
241           CASE("sf")
242              BASISJ(1) = 0
243              BASISJ(2) = 3
244              BASISJ(3) = -1
245           CASE("pd")
246              BASISJ(1) = 1
247              BASISJ(2) = 2
248              BASISJ(3) = -1
249           CASE("pf")
250              BASISJ(1) = 1
251              BASISJ(2) = 3
252              BASISJ(3) = -1
253           CASE("df")
254              BASISJ(1) = 2
255              BASISJ(2) = 3
256              BASISJ(3) = -1
257           CASE("spd")
258              BASISJ(1) = 0
259              BASISJ(2) = 1
260              BASISJ(3) = 2
261              BASISJ(4) = -1
262           CASE("spf")
263              BASISJ(1) = 0
264              BASISJ(2) = 1
265              BASISJ(3) = 3
266              BASISJ(4) = -1
267           CASE("sdf")
268              BASISJ(1) = 0
269              BASISJ(2) = 2
270              BASISJ(3) = 3
271              BASISJ(4) = -1
272           CASE("pdf")
273              BASISJ(1) = 1
274              BASISJ(2) = 2
275              BASISJ(3) = 3
276              BASISJ(4) = -1
277           CASE("spdf")
278              BASISJ(1) = 0
279              BASISJ(2) = 1
280              BASISJ(3) = 2
281              BASISJ(4) = 3
282              BASISJ(5) = -1
283           END SELECT
284
285           INDJ = MATINDLIST(J)
286
287           MAGRP2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2)
288           MAGRP = SQRT(MAGRP2)
289
290           ! transform to system in which z-axis is aligned with RIJ
291
292           PATH = .FALSE.
293           IF (ABS(RIJ(1)) .GT. 1E-12) THEN
294              IF (RIJ(1) .GT. ZERO .AND. RIJ(2) .GE. ZERO) THEN
295                 PHI = ZERO
296              ELSEIF (RIJ(1) .GT. ZERO .AND. RIJ(2) .LT. ZERO) THEN
297                 PHI = TWO * PI
298              ELSE
299                 PHI = PI
300              ENDIF
301              ALPHA = ATAN(RIJ(2) / RIJ(1)) + PHI
302           ELSEIF (ABS(RIJ(2)) .GT. 1E-12) THEN
303              IF (RIJ(2) .GT. 1E-12) THEN
304                 ALPHA = PI / TWO
305              ELSE
306                 ALPHA = THREE * PI / TWO
307              ENDIF
308           ELSE
309              ! pathological case: pole in alpha at beta=0
310              PATH = .TRUE.
311           ENDIF
312
313           COSBETA = RIJ(3)/MAGR
314           BETA = ACOS(RIJ(3) / MAGR)
315
316           DC = RIJ/MAGR
317
318           ! build forces using PRB 72 165107 eq. (12) - the sign of the
319           ! dfda contribution seems to be wrong, but gives the right
320           ! answer(?)
321
322           FTMP = ZERO
323           K = INDI
324
325           LBRAINC = 1
326           DO WHILE (BASISI(LBRAINC) .NE. -1)
327
328              LBRA = BASISI(LBRAINC)
329              LBRAINC = LBRAINC + 1
330
331              DO MBRA = -LBRA, LBRA
332
333                 K = K + 1
334                 L = INDJ
335
336                 LKETINC = 1
337                 DO WHILE (BASISJ(LKETINC) .NE. -1)
338
339                    LKET = BASISJ(LKETINC)
340                    LKETINC = LKETINC + 1
341
342                    DO MKET = -LKET, LKET
343
344                       L = L + 1
345
346                       !                       SELECT CASE(SPINON)
347                       !                       CASE(0)
348                       !                          RHO = BO(L, K)
349                       !                       CASE(1)
350                       !                          RHO = RHOUP(L, K) + RHODOWN(L, K)
351                       !                       END SELECT
352
353                       IF (.NOT. PATH) THEN
354
355                          ! Unroll loops and pre-compute
356
357                          MYDFDA = DFDA(I, J, LBRA, LKET, MBRA, &
358                               MKET, MAGR, ALPHA, COSBETA, "S")
359
360                          MYDFDB = DFDB(I, J, LBRA, LKET, MBRA, &
361                               MKET, MAGR, ALPHA, COSBETA, "S")
362
363                          MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, &
364                               MKET, MAGR, ALPHA, COSBETA, "S")
365
366                          KCOUNT = 0
367
368                          DO KX = 1, NKX
369
370                             DO KY = 1, NKY
371
372                                DO KZ = 1, NKZ
373
374                                   KPOINT = TWO*PI*(REAL(KX-1)*B1/REAL(NKX) + &
375                                        REAL(KY-1)*B2/REAL(NKY) + &
376                                        REAL(KZ-1)*B3/REAL(NKZ)) + K0
377
378                                   KCOUNT = KCOUNT+1
379
380                                   KDOTL = KPOINT(1)*RIJ(1) + KPOINT(2)*RIJ(2) + &
381                                        KPOINT(3)*RIJ(3)
382
383                                   CONJGBLOCH = EXP(CMPLX(ZERO,-KDOTL))
384
385                                   RHO = KBO(K,L,KCOUNT)*CONJGBLOCH
386
387                                   !
388                                   ! d/d_alpha
389                                   !
390
391                                   FTMP(1) = FTMP(1) + RHO * &
392                                        (-RIJ(2) / MAGRP2 * MYDFDA)
393
394                                   FTMP(2) = FTMP(2) + RHO * &
395                                        (RIJ(1)/ MAGRP2 * MYDFDA)
396
397                                   !
398                                   ! d/d_beta
399                                   !
400
401                                   FTMP(1) = FTMP(1) + RHO * &
402                                        (((((RIJ(3) * RIJ(1)) / &
403                                        MAGR2)) / MAGRP) * MYDFDB)
404
405                                   FTMP(2) = FTMP(2) + RHO * &
406                                        (((((RIJ(3) * RIJ(2)) / &
407                                        MAGR2)) / MAGRP) * MYDFDB)
408
409                                   FTMP(3) = FTMP(3) - RHO * &
410                                        (((ONE - ((RIJ(3) * RIJ(3)) / &
411                                        MAGR2)) / MAGRP) * MYDFDB)
412
413                                   !
414                                   ! d/dR
415                                   !
416
417                                   FTMP(1) = FTMP(1) - RHO * DC(1) * &
418                                        MYDFDR
419
420                                   FTMP(2) = FTMP(2) - RHO * DC(2) * &
421                                        MYDFDR
422
423                                   FTMP(3) = FTMP(3) - RHO * DC(3) * &
424                                        MYDFDR
425
426                                ENDDO
427                             ENDDO
428                          ENDDO
429
430                       ELSE
431
432                          ! pathological configuration in which beta=0
433                          ! or pi => alpha undefined
434
435                          ! fixed: MJC 12/17/13
436
437                          MYDFDB = DFDB(I, J, LBRA, LKET, &
438                               MBRA, MKET, MAGR, ZERO, COSBETA, "S") / MAGR
439
440                          MYDFDB = DFDB(I, J, LBRA, LKET, &
441                               MBRA, MKET, MAGR, PI/TWO, COSBETA, "S") / MAGR
442
443                          MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, &
444                               MKET, MAGR, ZERO, COSBETA, "S")
445
446                          KCOUNT = 0
447
448                          DO KX = 1, NKX
449                             DO KY = 1, NKY
450                                DO KZ = 1, NKZ
451
452                                   KPOINT = TWO*PI*(REAL(KX-1)*B1/REAL(NKX) + &
453                                        REAL(KY-1)*B2/REAL(NKY) + &
454                                        REAL(KZ-1)*B3/REAL(NKZ)) + K0
455
456                                   KCOUNT = KCOUNT+1
457
458                                   KDOTL = KPOINT(1)*RIJ(1) + KPOINT(2)*RIJ(2) + &
459                                        KPOINT(3)*RIJ(3)
460
461                                   CONJGBLOCH = EXP(CMPLX(ZERO,-KDOTL))
462
463                                   RHO = KBO(K,L,KCOUNT)*CONJGBLOCH
464
465                                   FTMP(1) = FTMP(1) - RHO * COSBETA * MYDFDB
466                                   FTMP(2) = FTMP(2) - RHO * COSBETA * MYDFDB
467                                   FTMP(3) = FTMP(3) - RHO * COSBETA * MYDFDR
468
469                                ENDDO
470                             ENDDO
471                          ENDDO
472
473                       ENDIF
474
475                    ENDDO
476                 ENDDO
477              ENDDO
478           ENDDO
479
480           FTMP = FTMP * ( HUBBARDU(ELEMPOINTER(J))*DELTAQ(J) + COULOMBV(J) &
481                +HUBBARDU(ELEMPOINTER(I))*DELTAQ(I) + COULOMBV(I))
482
483
484           KF(1,I) = KF(1,I) + FTMP(1)
485           KF(2,I) = KF(2,I) + FTMP(2)
486           KF(3,I) = KF(3,I) + FTMP(3)
487
488           ! with the factor of 2...
489
490           VIRBONDK(1) = VIRBONDK(1) + RIJ(1)*FTMP(1)
491           VIRBONDK(2) = VIRBONDK(2) + RIJ(2)*FTMP(2)
492           VIRBONDK(3) = VIRBONDK(3) + RIJ(3)*FTMP(3)
493           VIRBONDK(4) = VIRBONDK(4) + RIJ(1)*FTMP(2)
494           VIRBONDK(5) = VIRBONDK(5) + RIJ(2)*FTMP(3)
495           VIRBONDK(6) = VIRBONDK(6) + RIJ(3)*FTMP(1)
496
497
498        ENDIF
499     ENDDO
500
501  ENDDO
502
503  !$OMP END PARALLEL DO
504
505  VIRBONDK = VIRBONDK/TWO
506
507  !  PRINT*, KF(1,1), KF(2,1), KF(3,1)
508
509  FSCOUL = REAL(KF)/REAL(NKTOT)
510  VIRSCOUL = REAL(VIRBONDK)/REAL(NKTOT)
511
512  !  DO I = 1, NATS
513  !     WRITE(6,10) I, FSCOUL(1,I), FSCOUL(2,I), FSCOUL(3,I)
514  !  ENDDO
515
516  !10 FORMAT(I4, 3F12.6)
517
518  RETURN
519
520END SUBROUTINE KFCOULNONO
521