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
22! 1 = A0
23! 2 = B1
24! 3 = B2
25! 4 = B3
26! 5 = B4
27! 6 = B5
28! 7 = R1
29! 8 = RCUT
30! 9 = TAIL1
31! 10 = TAIL2
32! 11 = TAIL3
33! 12 = TAIL4
34! 13 = TAIL5
35! 14 = TAIL6
36
37FUNCTION DUNIVSCALE(I, J, L1, L2, MP, R, WHICHINT)
38
39  USE CONSTANTS_MOD
40  USE SETUPARRAY
41  USE UNIVARRAY
42  USE MYPRECISION
43
44  IMPLICIT NONE
45
46  INTEGER :: I, J, K, L1, L2, IP1, IP2, MP, IC, MYINTEGRAL
47  INTEGER :: BREAKLOOP
48  INTEGER :: KLO, KHI
49  REAL(LATTEPREC) :: DUNIVSCALE
50  REAL(LATTEPREC) :: SA, SB, DX
51  REAL(LATTEPREC) :: A(14), R, RMINUSR1,DPOLYNOM, RMOD
52  CHARACTER(LEN=1) :: WHICHINT
53  CHARACTER(LEN=3) :: IGLTYPE
54  IF (EXISTERROR) RETURN
55
56  ! can't test directly on L values because basis strings always list
57  ! lower L values first
58
59  IF (L1 .GT. L2) THEN
60     IP1 = L2
61     IP2 = L1
62  ELSE
63     IP1 = L1
64     IP2 = L2
65  ENDIF
66
67  ! build basis string from L and M values - pure hackery
68
69  SELECT CASE(IP1)
70  CASE(0)
71     IGLTYPE = "s"
72  CASE(1)
73     IGLTYPE = "p"
74  CASE(2)
75     IGLTYPE = "d"
76  CASE(3)
77     IGLTYPE = "f"
78  END SELECT
79
80  SELECT CASE(IP2)
81  CASE(0)
82     IGLTYPE = TRIM(IGLTYPE)//"s"
83  CASE(1)
84     IGLTYPE = TRIM(IGLTYPE)//"p"
85  CASE(2)
86     IGLTYPE = TRIM(IGLTYPE)//"d"
87  CASE(3)
88     IGLTYPE = TRIM(IGLTYPE)//"f"
89  END SELECT
90
91  SELECT CASE(MP)
92  CASE(0)
93     IGLTYPE = TRIM(IGLTYPE)//"s"
94  CASE(1)
95     IGLTYPE = TRIM(IGLTYPE)//"p"
96  CASE(2)
97     IGLTYPE = TRIM(IGLTYPE)//"d"
98  CASE(3)
99     IGLTYPE = TRIM(IGLTYPE)//"f"
100  END SELECT
101
102  ! It makes a difference if our atoms are of the species or not...
103
104  ! Easier case first ATELE(I) = ATELE(J)
105
106  MYINTEGRAL = 0
107
108  IF (ATELE(I) .EQ. ATELE(J)) THEN
109
110     BREAKLOOP = 0
111     IC = 0
112     DO WHILE (BREAKLOOP .EQ. 0)
113
114        IC = IC + 1
115        IF (ATELE(I) .EQ. ELE1(IC) .AND. ATELE(J) .EQ. ELE2(IC) .AND. &
116             IGLTYPE .EQ. BTYPE(IC)) THEN
117
118           ! Now we've ID'ed our bond integral
119
120           MYINTEGRAL = IC
121
122           BREAKLOOP = 1
123
124        ENDIF
125     ENDDO
126
127  ELSE
128
129     ! Elements are different - care must be taken with p-s, s-p etc.
130
131     IF (L1 .EQ. L2) THEN ! This is a special case sss, pps, ppp etc.
132
133        BREAKLOOP = 0
134        IC = 0
135        DO WHILE (BREAKLOOP .EQ. 0)
136           IC = IC + 1
137           IF (((ATELE(I) .EQ. ELE1(IC) .AND. ATELE(J) .EQ. ELE2(IC)) .OR. &
138                (ATELE(I) .EQ. ELE2(IC) .AND. ATELE(J) .EQ. ELE1(IC))) .AND. &
139                IGLTYPE .EQ. BTYPE(IC)) THEN
140
141              ! Now we've ID'ed our bond integral
142
143              MYINTEGRAL = IC
144
145              BREAKLOOP = 1
146
147           ENDIF
148        ENDDO
149
150     ELSE  ! L1 .NE. L2
151
152        IF (L1 .LT. L2) THEN
153
154           DO IC = 1, NOINT
155
156              IF ((ATELE(I) .EQ. ELE1(IC) .AND. ATELE(J) .EQ. ELE2(IC)) .AND. &
157                   IGLTYPE .EQ. BTYPE(IC)) THEN
158
159                 ! Now we've ID'ed our bond integral
160
161                 MYINTEGRAL = IC
162
163              ENDIF
164           ENDDO
165
166        ELSE
167
168           DO IC = 1, NOINT
169
170              IF ((ATELE(I) .EQ. ELE2(IC) .AND. ATELE(J) .EQ. ELE1(IC)) .AND. &
171                   IGLTYPE .EQ. BTYPE(IC)) THEN
172
173                 ! Now we've ID'ed our bond integral
174
175                 MYINTEGRAL = IC
176
177              ENDIF
178           ENDDO
179
180        ENDIF
181     ENDIF
182
183  ENDIF
184
185  IF (SCLTYPE .EQ. "EXP") THEN
186
187     SELECT CASE(WHICHINT)
188     CASE("H") ! We're doing the H matrix build
189        A = BOND(:,MYINTEGRAL)
190     CASE("S") ! We're doing the S matrix build
191        A = OVERL(:,MYINTEGRAL)
192     END SELECT
193
194     IF (R .LE. A(7)) THEN
195
196        RMOD = R - A(6)
197
198        DPOLYNOM = A(2) + RMOD*(TWO*A(3) + RMOD*(THREE*A(4) + FOUR*A(5)*RMOD))
199
200        DUNIVSCALE = DPOLYNOM * EXP( RMOD*(A(2) + &
201             RMOD*(A(3) + RMOD*(A(4) + A(5)*RMOD))) )
202
203     ELSEIF (R .GT. A(7) .AND. R .LT. A(8)) THEN
204
205        RMINUSR1 = R - A(7)
206
207        DUNIVSCALE = A(10) + RMINUSR1*(TWO*A(11) + &
208             RMINUSR1*(THREE*A(12) + RMINUSR1*(FOUR*A(13) + &
209             RMINUSR1*FIVE*A(14))))
210
211     ELSE
212
213        DUNIVSCALE = ZERO
214
215     END IF
216
217     DUNIVSCALE = -A(1)*DUNIVSCALE
218
219  ELSEIF (SCLTYPE .EQ. "TABLE") THEN
220
221     KLO = 1
222     KHI = LENTABINT(MYINTEGRAL)
223
224     DO WHILE (KHI - KLO .GT. 1)
225
226        K = (KHI + KLO)/2
227
228        IF (TABR(K,MYINTEGRAL) .GT. R) THEN
229           KHI = K
230        ELSE
231           KLO = K
232        ENDIF
233
234     ENDDO
235
236     DX = TABR(KHI, MYINTEGRAL) - TABR(KLO,MYINTEGRAL)
237
238     SA = (TABR(KHI, MYINTEGRAL) - R)/DX
239     SB = (R - TABR(KLO, MYINTEGRAL))/DX
240
241     ! Negative gradient
242
243     IF (WHICHINT .EQ. "H") THEN
244
245        DUNIVSCALE = -((TABH(KHI,MYINTEGRAL) - TABH(KLO,MYINTEGRAL))/DX + &
246             ((ONE - THREE*SA*SA)*HSPL(KLO,MYINTEGRAL) + &
247             (THREE*SB*SB - ONE)*HSPL(KHI,MYINTEGRAL))*(DX/SIX))
248
249        IF (R .GT. HCUT(MYINTEGRAL)) DUNIVSCALE = ZERO
250
251     ELSEIF (WHICHINT .EQ. "S") THEN
252
253        DUNIVSCALE = -((TABS(KHI,MYINTEGRAL) - TABS(KLO,MYINTEGRAL))/DX + &
254             ((ONE - THREE*SA*SA)*SSPL(KLO,MYINTEGRAL) + &
255             (THREE*SB*SB - ONE)*SSPL(KHI,MYINTEGRAL))*(DX/SIX))
256
257        IF (R .GT. SCUT(MYINTEGRAL)) DUNIVSCALE = ZERO
258
259     ENDIF
260
261  ENDIF
262
263  ! permutation symmetry
264
265  IF (L1 .GT. L2 .AND. MOD(L1 + L2, 2) .NE. 0) DUNIVSCALE = -DUNIVSCALE
266
267  !  PRINT*, UNIVSCALE
268
269  RETURN
270
271END FUNCTION DUNIVSCALE
272
273