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 FLCNNONO
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 MYPRECISION
33
34  IMPLICIT NONE
35
36  INTEGER :: I, J, K, L, M, N, KK, INDI, INDJ
37  INTEGER :: LBRA, MBRA, LKET, MKET
38  INTEGER :: PREVJ, NEWJ
39  INTEGER :: PBCI, PBCJ, PBCK
40  INTEGER :: BASISI(5), BASISJ(5), LBRAINC, LKETINC
41  REAL(LATTEPREC) :: ALPHA, BETA, PHI, RHO, COSBETA
42  REAL(LATTEPREC) :: RIJ(3), DC(3)
43  REAL(LATTEPREC) :: MAGR, MAGR2, MAGRP, MAGRP2, FTMP(3)
44  REAL(LATTEPREC) :: MAXRCUT, MAXRCUT2
45  REAL(LATTEPREC) :: MYDFDA, MYDFDB, MYDFDR, RCUTTB
46  REAL(LATTEPREC), EXTERNAL :: DFDA, DFDB, DFDR
47  LOGICAL PATH
48
49
50  FSLCN = ZERO
51  VIRSLCN = ZERO
52
53!$OMP PARALLEL DO DEFAULT (NONE) &
54!$OMP SHARED(NATS, BASIS, ELEMPOINTER, TOTNEBTB, NEBTB) &
55!$OMP SHARED(CR, BOX, BO, RHOUP, RHODOWN, SPINON, NOINT, ATELE, ELE1, ELE2) &
56!$OMP SHARED(HCUT, SCUT, MATINDLIST, BASISTYPE) &
57!$OMP SHARED(LCNSHIFT) &
58!$OMP PRIVATE(I, J, K, NEWJ, BASISI, BASISJ, INDI, INDJ, PBCI, PBCJ, PBCK) &
59!$OMP PRIVATE(RIJ, MAGR2, MAGR, MAGRP2, MAGRP, PATH, PHI, ALPHA, BETA, COSBETA, FTMP) &
60!$OMP PRIVATE(DC, LBRAINC, LBRA, MBRA, L, LKETINC, LKET, MKET, RHO) &
61!$OMP PRIVATE(MYDFDA, MYDFDB, MYDFDR, RCUTTB) &
62!$OMP REDUCTION(+:FSLCN, VIRSLCN)
63
64  DO I = 1, NATS
65
66     ! Build list of orbitals on atom I
67     SELECT CASE(BASIS(ELEMPOINTER(I)))
68
69     CASE("s")
70        BASISI(1) = 0
71        BASISI(2) = -1
72     CASE("p")
73        BASISI(1) = 1
74        BASISI(2) = -1
75     CASE("d")
76        BASISI(1) = 2
77        BASISI(2) = -1
78     CASE("f")
79        BASISI(1) = 3
80        BASISI(2) = -1
81     CASE("sp")
82        BASISI(1) = 0
83        BASISI(2) = 1
84        BASISI(3) = -1
85     CASE("sd")
86        BASISI(1) = 0
87        BASISI(2) = 2
88        BASISI(3) = -1
89     CASE("sf")
90        BASISI(1) = 0
91        BASISI(2) = 3
92        BASISI(3) = -1
93     CASE("pd")
94        BASISI(1) = 1
95        BASISI(2) = 2
96        BASISI(3) = -1
97     CASE("pf")
98        BASISI(1) = 1
99        BASISI(2) = 3
100        BASISI(3) = -1
101     CASE("df")
102        BASISI(1) = 2
103        BASISI(2) = 3
104        BASISI(3) = -1
105     CASE("spd")
106        BASISI(1) = 0
107        BASISI(2) = 1
108        BASISI(3) = 2
109        BASISI(4) = -1
110     CASE("spf")
111        BASISI(1) = 0
112        BASISI(2) = 1
113        BASISI(3) = 3
114        BASISI(4) = -1
115     CASE("sdf")
116        BASISI(1) = 0
117        BASISI(2) = 2
118        BASISI(3) = 3
119        BASISI(4) = -1
120     CASE("pdf")
121        BASISI(1) = 1
122        BASISI(2) = 2
123        BASISI(3) = 3
124        BASISI(4) = -1
125     CASE("spdf")
126        BASISI(1) = 0
127        BASISI(2) = 1
128        BASISI(3) = 2
129        BASISI(4) = 3
130        BASISI(5) = -1
131     END SELECT
132
133     INDI = MATINDLIST(I)
134
135     DO NEWJ = 1, TOTNEBTB(I)
136
137        J = NEBTB(1, NEWJ, I)
138        PBCI = NEBTB(2, NEWJ, I)
139        PBCJ = NEBTB(3, NEWJ, I)
140        PBCK = NEBTB(4, NEWJ, I)
141
142        RIJ(1) = CR(1,J) + REAL(PBCI)*BOX(1,1) + REAL(PBCJ)*BOX(2,1) + &
143             REAL(PBCK)*BOX(3,1) - CR(1,I)
144
145        RIJ(2) = CR(2,J) + REAL(PBCI)*BOX(1,2) + REAL(PBCJ)*BOX(2,2) + &
146             REAL(PBCK)*BOX(3,2) - CR(2,I)
147
148        RIJ(3) = CR(3,J) + REAL(PBCI)*BOX(1,3) + REAL(PBCJ)*BOX(2,3) + &
149             REAL(PBCK)*BOX(3,3) - CR(3,I)
150
151        MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3)
152
153        RCUTTB = ZERO
154        DO K = 1, NOINT
155
156           IF ( (ATELE(I) .EQ. ELE1(K) .AND. &
157                ATELE(J) .EQ. ELE2(K)) .OR. &
158                (ATELE(J) .EQ. ELE1(K) .AND. &
159                ATELE(I) .EQ. ELE2(K) )) THEN
160
161              IF (HCUT(K) .GT. RCUTTB ) RCUTTB = HCUT(K)
162
163              IF (BASISTYPE .EQ. "NONORTHO") THEN
164                 IF (SCUT(K) .GT. RCUTTB ) RCUTTB = SCUT(K)
165              ENDIF
166
167           ENDIF
168
169        ENDDO
170
171        IF (MAGR2 .LT. RCUTTB*RCUTTB) THEN
172
173           MAGR = SQRT(MAGR2)
174
175           ! Build list of orbitals on atom J
176
177           SELECT CASE(BASIS(ELEMPOINTER(J)))
178           CASE("s")
179              BASISJ(1) = 0
180              BASISJ(2) = -1
181           CASE("p")
182              BASISJ(1) = 1
183              BASISJ(2) = -1
184           CASE("d")
185              BASISJ(1) = 2
186              BASISJ(2) = -1
187           CASE("f")
188              BASISJ(1) = 3
189              BASISJ(2) = -1
190           CASE("sp")
191              BASISJ(1) = 0
192              BASISJ(2) = 1
193              BASISJ(3) = -1
194           CASE("sd")
195              BASISJ(1) = 0
196              BASISJ(2) = 2
197              BASISJ(3) = -1
198           CASE("sf")
199              BASISJ(1) = 0
200              BASISJ(2) = 3
201              BASISJ(3) = -1
202           CASE("pd")
203              BASISJ(1) = 1
204              BASISJ(2) = 2
205              BASISJ(3) = -1
206           CASE("pf")
207              BASISJ(1) = 1
208              BASISJ(2) = 3
209              BASISJ(3) = -1
210           CASE("df")
211              BASISJ(1) = 2
212              BASISJ(2) = 3
213              BASISJ(3) = -1
214           CASE("spd")
215              BASISJ(1) = 0
216              BASISJ(2) = 1
217              BASISJ(3) = 2
218              BASISJ(4) = -1
219           CASE("spf")
220              BASISJ(1) = 0
221              BASISJ(2) = 1
222              BASISJ(3) = 3
223              BASISJ(4) = -1
224           CASE("sdf")
225              BASISJ(1) = 0
226              BASISJ(2) = 2
227              BASISJ(3) = 3
228              BASISJ(4) = -1
229           CASE("pdf")
230              BASISJ(1) = 1
231              BASISJ(2) = 2
232              BASISJ(3) = 3
233              BASISJ(4) = -1
234           CASE("spdf")
235              BASISJ(1) = 0
236              BASISJ(2) = 1
237              BASISJ(3) = 2
238              BASISJ(4) = 3
239              BASISJ(5) = -1
240           END SELECT
241
242           INDJ = MATINDLIST(J)
243
244           MAGRP2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2)
245           MAGRP = SQRT(MAGRP2)
246
247           ! transform to system in which z-axis is aligned with RIJ
248
249           PATH = .FALSE.
250           IF (ABS(RIJ(1)) .GT. 1E-12) THEN
251              IF (RIJ(1) .GT. ZERO .AND. RIJ(2) .GE. ZERO) THEN
252                 PHI = ZERO
253              ELSEIF (RIJ(1) .GT. ZERO .AND. RIJ(2) .LT. ZERO) THEN
254                 PHI = TWO * PI
255              ELSE
256                 PHI = PI
257              ENDIF
258              ALPHA = ATAN(RIJ(2) / RIJ(1)) + PHI
259           ELSEIF (ABS(RIJ(2)) .GT. 1E-12) THEN
260              IF (RIJ(2) .GT. 1E-12) THEN
261                 ALPHA = PI / TWO
262              ELSE
263                 ALPHA = THREE * PI / TWO
264              ENDIF
265           ELSE
266              ! pathological case: pole in alpha at beta=0
267              PATH = .TRUE.
268           ENDIF
269
270           COSBETA = RIJ(3)/MAGR
271           BETA = ACOS(RIJ(3) / MAGR)
272
273           DC = RIJ/MAGR
274
275           ! build forces using PRB 72 165107 eq. (12) - the sign of the
276           ! dfda contribution seems to be wrong, but gives the right
277           ! answer(?)
278
279           FTMP = ZERO
280           K = INDI
281
282           LBRAINC = 1
283           DO WHILE (BASISI(LBRAINC) .NE. -1)
284
285              LBRA = BASISI(LBRAINC)
286              LBRAINC = LBRAINC + 1
287
288              DO MBRA = -LBRA, LBRA
289
290                 K = K + 1
291                 L = INDJ
292
293                 LKETINC = 1
294                 DO WHILE (BASISJ(LKETINC) .NE. -1)
295
296                    LKET = BASISJ(LKETINC)
297                    LKETINC = LKETINC + 1
298
299                    DO MKET = -LKET, LKET
300
301                       L = L + 1
302
303                       SELECT CASE(SPINON)
304                       CASE(0)
305                          RHO = BO(L, K)
306                       CASE(1)
307                          RHO = RHOUP(L, K) + RHODOWN(L, K)
308                       END SELECT
309
310                       IF (.NOT. PATH) THEN
311
312                          ! Unroll loops and pre-compute
313
314                          MYDFDA = DFDA(I, J, LBRA, LKET, MBRA, &
315                               MKET, MAGR, ALPHA, COSBETA, "S")
316
317                          MYDFDB = DFDB(I, J, LBRA, LKET, MBRA, &
318                               MKET, MAGR, ALPHA, COSBETA, "S")
319
320                          MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, &
321                               MKET, MAGR, ALPHA, COSBETA, "S")
322
323                          !
324                          ! d/d_alpha
325                          !
326
327                          FTMP(1) = FTMP(1) + RHO * &
328                               (-RIJ(2) / MAGRP2 * MYDFDA)
329
330                          FTMP(2) = FTMP(2) + RHO * &
331                               (RIJ(1)/ MAGRP2 * MYDFDA)
332
333                          !
334                          ! d/d_beta
335                          !
336
337                          FTMP(1) = FTMP(1) + RHO * &
338                               (((((RIJ(3) * RIJ(1)) / &
339                               MAGR2)) / MAGRP) * MYDFDB)
340
341                          FTMP(2) = FTMP(2) + RHO * &
342                               (((((RIJ(3) * RIJ(2)) / &
343                               MAGR2)) / MAGRP) * MYDFDB)
344
345                          FTMP(3) = FTMP(3) - RHO * &
346                               (((ONE - ((RIJ(3) * RIJ(3)) / &
347                               MAGR2)) / MAGRP) * MYDFDB)
348
349                          !
350                          ! d/dR
351                          !
352
353                          FTMP(1) = FTMP(1) - RHO * DC(1) * &
354                               MYDFDR
355
356                          FTMP(2) = FTMP(2) - RHO * DC(2) * &
357                               MYDFDR
358
359                          FTMP(3) = FTMP(3) - RHO * DC(3) * &
360                               MYDFDR
361
362
363                       ELSE
364
365                          ! pathological configuration in which beta=0
366                          ! or pi => alpha undefined
367
368                          ! fixed: MJC 12/17/13
369
370                          MYDFDB = DFDB(I, J, LBRA, LKET, &
371                               MBRA, MKET, MAGR, ZERO, COSBETA, "S") / MAGR
372
373                          FTMP(1) = FTMP(1) - RHO * (COSBETA * MYDFDB)
374
375                          MYDFDB = DFDB(I, J, LBRA, LKET, &
376                               MBRA, MKET, MAGR, PI/TWO, COSBETA, "S") / MAGR
377
378                          FTMP(2) = FTMP(2) - RHO * (COSBETA * MYDFDB)
379
380                          MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, &
381                               MKET, MAGR, ZERO, COSBETA, "S")
382
383                          FTMP(3) = FTMP(3) - RHO * COSBETA * MYDFDR
384
385                       ENDIF
386
387                    ENDDO
388                 ENDDO
389              ENDDO
390           ENDDO
391
392           FTMP = FTMP * (LCNSHIFT(I) + LCNSHIFT(J))
393!           FTMP = FTMP * ( HUBBARDU(ELEMPOINTER(J))*DELTAQ(J) + COULOMBV(J) &
394!                +HUBBARDU(ELEMPOINTER(I))*DELTAQ(I) + COULOMBV(I))
395
396
397           FSLCN(1,I) = FSLCN(1,I) + FTMP(1)
398           FSLCN(2,I) = FSLCN(2,I) + FTMP(2)
399           FSLCN(3,I) = FSLCN(3,I) + FTMP(3)
400
401           ! with the factor of 2...
402
403           VIRSLCN(1) = VIRSLCN(1) + RIJ(1)*FTMP(1)
404           VIRSLCN(2) = VIRSLCN(2) + RIJ(2)*FTMP(2)
405           VIRSLCN(3) = VIRSLCN(3) + RIJ(3)*FTMP(3)
406           VIRSLCN(4) = VIRSLCN(4) + RIJ(1)*FTMP(2)
407           VIRSLCN(5) = VIRSLCN(5) + RIJ(2)*FTMP(3)
408           VIRSLCN(6) = VIRSLCN(6) + RIJ(3)*FTMP(1)
409
410
411        ENDIF
412     ENDDO
413
414  ENDDO
415
416!$OMP END PARALLEL DO
417
418  VIRSLCN = VIRSLCN/TWO
419
420!  DO I = 1, NATS
421!     WRITE(6,10) I, FSLCN(1,I), FSLCN(2,I), FSLCN(3,I)
422!  ENDDO
423
424!10 FORMAT(I4, 3F12.6)
425
426  RETURN
427
428END SUBROUTINE FLCNNONO
429