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