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