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 GRADH
23
24  USE CONSTANTS_MOD
25  USE SETUPARRAY
26  USE NEBLISTARRAY
27  USE UNIVARRAY
28  USE SPINARRAY
29  USE VIRIALARRAY
30  USE MYPRECISION
31
32  IMPLICIT NONE
33
34  INTEGER :: I, J, K, L, M, N, KK, INDI, INDJ
35  INTEGER :: LBRA, MBRA, LKET, MKET
36  INTEGER :: PREVJ, NEWJ
37  INTEGER :: PBCI, PBCJ, PBCK
38  INTEGER :: BASISI(5), BASISJ(5), LBRAINC, LKETINC
39  REAL(LATTEPREC) :: ALPHA, BETA, PHI, RHO, COSBETA
40  REAL(LATTEPREC) :: RIJ(3), DC(3)
41  REAL(LATTEPREC) :: MAGR, MAGR2, MAGRP, MAGRP2, FTMP(3)
42  REAL(LATTEPREC) :: MAXRCUT, MAXRCUT2
43  REAL(LATTEPREC) :: MYDFDA, MYDFDB, MYDFDR, RCUTTB
44  REAL(LATTEPREC), EXTERNAL :: DFDA, DFDB, DFDR
45  LOGICAL PATH
46  IF (EXISTERROR) RETURN
47
48  F = ZERO
49  VIRBOND = ZERO
50
51!$OMP PARALLEL DO DEFAULT (NONE) &
52!$OMP SHARED(NATS, BASIS, ELEMPOINTER, TOTNEBTB, NEBTB) &
53!$OMP SHARED(CR, BOX, BO, RHOUP, RHODOWN, SPINON, NOINT, ATELE, ELE1, ELE2) &
54!$OMP SHARED(HCUT, SCUT, MATINDLIST, BASISTYPE) &
55!$OMP SHARED(F) &
56!$OMP PRIVATE(I, J, K, NEWJ, BASISI, BASISJ, INDI, INDJ, PBCI, PBCJ, PBCK) &
57!$OMP PRIVATE(RIJ, MAGR2, MAGR, MAGRP2, MAGRP, PATH, PHI, ALPHA, BETA, COSBETA, FTMP) &
58!$OMP PRIVATE(DC, LBRAINC, LBRA, MBRA, L, LKETINC, LKET, MKET, RHO) &
59!$OMP PRIVATE(MYDFDA, MYDFDB, MYDFDR, RCUTTB) &
60!$OMP REDUCTION(+:VIRBOND)
61
62  DO I = 1, NATS
63
64     ! Build list of orbitals on atom I
65     SELECT CASE(BASIS(ELEMPOINTER(I)))
66
67     CASE("s")
68        BASISI(1) = 0
69        BASISI(2) = -1
70     CASE("p")
71        BASISI(1) = 1
72        BASISI(2) = -1
73     CASE("d")
74        BASISI(1) = 2
75        BASISI(2) = -1
76     CASE("f")
77        BASISI(1) = 3
78        BASISI(2) = -1
79     CASE("sp")
80        BASISI(1) = 0
81        BASISI(2) = 1
82        BASISI(3) = -1
83     CASE("sd")
84        BASISI(1) = 0
85        BASISI(2) = 2
86        BASISI(3) = -1
87     CASE("sf")
88        BASISI(1) = 0
89        BASISI(2) = 3
90        BASISI(3) = -1
91     CASE("pd")
92        BASISI(1) = 1
93        BASISI(2) = 2
94        BASISI(3) = -1
95     CASE("pf")
96        BASISI(1) = 1
97        BASISI(2) = 3
98        BASISI(3) = -1
99     CASE("df")
100        BASISI(1) = 2
101        BASISI(2) = 3
102        BASISI(3) = -1
103     CASE("spd")
104        BASISI(1) = 0
105        BASISI(2) = 1
106        BASISI(3) = 2
107        BASISI(4) = -1
108     CASE("spf")
109        BASISI(1) = 0
110        BASISI(2) = 1
111        BASISI(3) = 3
112        BASISI(4) = -1
113     CASE("sdf")
114        BASISI(1) = 0
115        BASISI(2) = 2
116        BASISI(3) = 3
117        BASISI(4) = -1
118     CASE("pdf")
119        BASISI(1) = 1
120        BASISI(2) = 2
121        BASISI(3) = 3
122        BASISI(4) = -1
123     CASE("spdf")
124        BASISI(1) = 0
125        BASISI(2) = 1
126        BASISI(3) = 2
127        BASISI(4) = 3
128        BASISI(5) = -1
129     END SELECT
130
131     INDI = MATINDLIST(I)
132
133     DO NEWJ = 1, TOTNEBTB(I)
134
135        J = NEBTB(1, NEWJ, I)
136        PBCI = NEBTB(2, NEWJ, I)
137        PBCJ = NEBTB(3, NEWJ, I)
138        PBCK = NEBTB(4, NEWJ, I)
139
140        RIJ(1) = CR(1,J) + REAL(PBCI)*BOX(1,1) + REAL(PBCJ)*BOX(2,1) + &
141             REAL(PBCK)*BOX(3,1) - CR(1,I)
142
143        RIJ(2) = CR(2,J) + REAL(PBCI)*BOX(1,2) + REAL(PBCJ)*BOX(2,2) + &
144             REAL(PBCK)*BOX(3,2) - CR(2,I)
145
146        RIJ(3) = CR(3,J) + REAL(PBCI)*BOX(1,3) + REAL(PBCJ)*BOX(2,3) + &
147             REAL(PBCK)*BOX(3,3) - CR(3,I)
148
149        MAGR2 = RIJ(1)*RIJ(1) + RIJ(2)*RIJ(2) + RIJ(3)*RIJ(3)
150
151        ! Building the forces is expensive - use the cut-off
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
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           !        PRINT*, ALPHA, BETA
275
276           DC = RIJ/MAGR
277
278           ! build forces using PRB 72 165107 eq. (12) - the sign of the
279           ! dfda contribution seems to be wrong, but gives the right
280           ! answer(?)
281
282           FTMP = ZERO
283           K = INDI
284
285           LBRAINC = 1
286           DO WHILE (BASISI(LBRAINC) .NE. -1)
287
288              LBRA = BASISI(LBRAINC)
289              LBRAINC = LBRAINC + 1
290
291              DO MBRA = -LBRA, LBRA
292
293                 K = K + 1
294                 L = INDJ
295
296                 LKETINC = 1
297                 DO WHILE (BASISJ(LKETINC) .NE. -1)
298
299                    LKET = BASISJ(LKETINC)
300                    LKETINC = LKETINC + 1
301
302                    DO MKET = -LKET, LKET
303
304                       L = L + 1
305
306                       SELECT CASE(SPINON)
307                       CASE(0)
308                          RHO = BO(L, K)
309                       CASE(1)
310                          RHO = RHOUP(L, K) + RHODOWN(L, K)
311                       END SELECT
312
313                       IF (.NOT. PATH) THEN
314
315                          ! Unroll loops and pre-compute
316
317                          MYDFDA = DFDA(I, J, LBRA, LKET, MBRA, &
318                               MKET, MAGR, ALPHA, COSBETA, "H")
319
320                          MYDFDB = DFDB(I, J, LBRA, LKET, MBRA, &
321                               MKET, MAGR, ALPHA, COSBETA, "H")
322
323                          MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, &
324                               MKET, MAGR, ALPHA, COSBETA, "H")
325
326                          !
327                          ! d/d_alpha
328                          !
329
330                          FTMP(1) = FTMP(1) + RHO * &
331                               (-RIJ(2) / MAGRP2 * MYDFDA)
332
333                          FTMP(2) = FTMP(2) + RHO * &
334                               (RIJ(1)/ MAGRP2 * MYDFDA)
335
336                          !
337                          ! d/d_beta
338                          !
339
340                          FTMP(1) = FTMP(1) + RHO * &
341                               (((((RIJ(3) * RIJ(1)) / &
342                               MAGR2)) / MAGRP) * MYDFDB)
343
344                          FTMP(2) = FTMP(2) + RHO * &
345                               (((((RIJ(3) * RIJ(2)) / &
346                               MAGR2)) / MAGRP) * MYDFDB)
347
348                          FTMP(3) = FTMP(3) - RHO * &
349                               (((ONE - ((RIJ(3) * RIJ(3)) / &
350                               MAGR2)) / MAGRP) * MYDFDB)
351
352                          !
353                          ! d/dR
354                          !
355
356                          FTMP(1) = FTMP(1) - RHO * DC(1) * &
357                               MYDFDR
358
359                          FTMP(2) = FTMP(2) - RHO * DC(2) * &
360                               MYDFDR
361
362                          FTMP(3) = FTMP(3) - RHO * DC(3) * &
363                               MYDFDR
364
365
366                       ELSE
367
368                          ! pathological configuration in which beta=0
369                          ! or pi => alpha undefined
370
371                          ! fixed: MJC 12/17/13
372
373                          MYDFDB = DFDB(I, J, LBRA, LKET, &
374                               MBRA, MKET, MAGR, ZERO, COSBETA, "H") / MAGR
375
376                          FTMP(1) = FTMP(1) - RHO * (COSBETA * MYDFDB)
377
378                          MYDFDB = DFDB(I, J, LBRA, LKET, &
379                               MBRA, MKET, MAGR, PI/TWO, COSBETA, "H") / MAGR
380
381                          FTMP(2) = FTMP(2) - RHO * (COSBETA * MYDFDB)
382
383                          MYDFDR = DFDR(I, J, LBRA, LKET, MBRA, &
384                               MKET, MAGR, ZERO, COSBETA, "H")
385
386                          FTMP(3) = FTMP(3) - RHO * COSBETA * MYDFDR
387
388                       ENDIF
389
390                    ENDDO
391                 ENDDO
392              ENDDO
393           ENDDO
394
395           F(1,I) = F(1,I) + FTMP(1)
396           F(2,I) = F(2,I) + FTMP(2)
397           F(3,I) = F(3,I) + FTMP(3)
398
399           VIRBOND(1) = VIRBOND(1) + RIJ(1) * FTMP(1)
400           VIRBOND(2) = VIRBOND(2) + RIJ(2) * FTMP(2)
401           VIRBOND(3) = VIRBOND(3) + RIJ(3) * FTMP(3)
402           VIRBOND(4) = VIRBOND(4) + RIJ(1) * FTMP(2)
403           VIRBOND(5) = VIRBOND(5) + RIJ(2) * FTMP(3)
404           VIRBOND(6) = VIRBOND(6) + RIJ(3) * FTMP(1)
405
406        ENDIF
407
408     ENDDO
409
410     !     INDI = INDI + NORBI
411
412  ENDDO
413
414!$OMP END PARALLEL DO
415
416  RETURN
417
418END SUBROUTINE GRADH
419