1#!/usr/bin/env python 2# $Id$ 3# -*- coding: utf-8 4 5''' 6test libcint 7''' 8 9__author__ = "Qiming Sun <osirpt.sun@gmail.com>" 10 11import sys 12import os 13import ctypes 14import numpy 15 16_cint = numpy.ctypeslib.load_library('libcint', '/opengrok/src/dports/science/libcint/.build') 17 18 19PTR_EXPCUTOFF = 0 20PTR_COMMON_ORIG = 1 21PTR_SHIELDING_ORIG = 4 22PTR_RINV_ORIG = 4 23PTR_RINV_ZETA = 7 24PTR_RANGE_OMEGA = 8 25PTR_ENV_START = 20 26NGRIDS = 11 27PTR_GRIDS = 12 28 29CHARGE_OF = 0 30PTR_COORD = 1 31NUC_MOD_OF = 2 32PTR_ZETA = 3 33RAD_GRIDS = 4 34ANG_GRIDS = 5 35ATM_SLOTS = 6 36 37ATOM_OF = 0 38ANG_OF = 1 39NPRIM_OF = 2 40NCTR_OF = 3 41KAPPA_OF = 4 42PTR_EXP = 5 43PTR_COEFF = 6 44BAS_SLOTS = 8 45 46natm = 4 47nbas = 0 48atm = numpy.zeros((natm,ATM_SLOTS), dtype=numpy.int32) 49bas = numpy.zeros((1000,BAS_SLOTS), dtype=numpy.int32) 50env = numpy.zeros(10000) 51off = PTR_ENV_START 52for i in range(natm): 53 atm[i, CHARGE_OF] = (i+1)*2 54 atm[i, PTR_COORD] = off 55 env[off+0] = .2 * (i+1) 56 env[off+1] = .3 + (i+1) * .5 57 env[off+2] = .1 - (i+1) * .5 58 off += 3 59off0 = off 60 61# basis with kappa > 0 62nh = 0 63 64bas[nh,ATOM_OF ] = 0 65bas[nh,ANG_OF ] = 1 66bas[nh,KAPPA_OF] = 1 67bas[nh,NPRIM_OF] = 1 68bas[nh,NCTR_OF ] = 1 69bas[nh,PTR_EXP] = off 70env[off+0] = 1 71bas[nh,PTR_COEFF] = off + 1 72env[off+1] = 1 73off += 2 74nh += 1 75 76bas[nh,ATOM_OF ] = 1 77bas[nh,ANG_OF ] = 2 78bas[nh,KAPPA_OF] = 2 79bas[nh,NPRIM_OF] = 2 80bas[nh,NCTR_OF ] = 2 81bas[nh,PTR_EXP] = off 82env[off+0] = 5 83env[off+1] = 3 84bas[nh,PTR_COEFF] = off + 2 85env[off+2] = 1 86env[off+3] = 2 87env[off+4] = 4 88env[off+5] = 1 89off += 6 90nh += 1 91 92bas[nh,ATOM_OF ] = 2 93bas[nh,ANG_OF ] = 3 94bas[nh,KAPPA_OF] = 3 95bas[nh,NPRIM_OF] = 1 96bas[nh,NCTR_OF ] = 1 97bas[nh,PTR_EXP ] = off 98env[off+0] = 1 99bas[nh,PTR_COEFF] = off + 1 100env[off+1] = 1 101off += 2 102nh += 1 103 104bas[nh,ATOM_OF ] = 3 105bas[nh,ANG_OF ] = 4 106bas[nh,KAPPA_OF] = 4 107bas[nh,NPRIM_OF] = 1 108bas[nh,NCTR_OF ] = 1 109bas[nh,PTR_EXP ] = off 110env[off+0] = .5 111bas[nh,PTR_COEFF] = off + 1 112env[off+1] = 1. 113off = off + 2 114nh += 1 115 116nbas = nh 117 118# basis with kappa < 0 119n = off - off0 120for i in range(n): 121 env[off+i] = env[off0+i] 122 123for i in range(nh): 124 bas[i+nh,ATOM_OF ] = bas[i,ATOM_OF ] 125 bas[i+nh,ANG_OF ] = bas[i,ANG_OF ] - 1 126 bas[i+nh,KAPPA_OF] =-bas[i,KAPPA_OF] 127 bas[i+nh,NPRIM_OF] = bas[i,NPRIM_OF] 128 bas[i+nh,NCTR_OF ] = bas[i,NCTR_OF ] 129 bas[i+nh,PTR_EXP ] = bas[i,PTR_EXP ] + n 130 bas[i+nh,PTR_COEFF]= bas[i,PTR_COEFF] + n 131 env[bas[i+nh,PTR_COEFF]] /= 2 * env[bas[i,PTR_EXP]] 132 133env[bas[5,PTR_COEFF]+0] = env[bas[1,PTR_COEFF]+0] / (2 * env[bas[1,PTR_EXP]+0]) 134env[bas[5,PTR_COEFF]+1] = env[bas[1,PTR_COEFF]+1] / (2 * env[bas[1,PTR_EXP]+1]) 135env[bas[5,PTR_COEFF]+2] = env[bas[1,PTR_COEFF]+2] / (2 * env[bas[1,PTR_EXP]+0]) 136env[bas[5,PTR_COEFF]+3] = env[bas[1,PTR_COEFF]+3] / (2 * env[bas[1,PTR_EXP]+1]) 137 138natm = ctypes.c_int(natm) 139nbas = ctypes.c_int(nbas) 140c_atm = atm.ctypes.data_as(ctypes.c_void_p) 141c_bas = bas.ctypes.data_as(ctypes.c_void_p) 142c_env = env.ctypes.data_as(ctypes.c_void_p) 143 144opt = ctypes.POINTER(ctypes.c_void_p)() 145_cint.CINTlen_spinor.restype = ctypes.c_int 146 147 148def close(v1, vref, count, place): 149 return round(abs(v1-vref)/count**.5, place) == 0 150 151def test_int1e_sph(name, vref, dim, place): 152 intor = getattr(_cint, name) 153 intor.restype = ctypes.c_void_p 154 op = (ctypes.c_double * (10000 * dim))() 155 v1 = 0 156 cnt = 0 157 for j in range(nbas.value*2): 158 for i in range(j+1): 159 di = (bas[i,ANG_OF] * 2 + 1) * bas[i,NCTR_OF] 160 dj = (bas[j,ANG_OF] * 2 + 1) * bas[j,NCTR_OF] 161 shls = (ctypes.c_int * 2)(i, j) 162 intor(op, shls, c_atm, natm, c_bas, nbas, c_env) 163 v1 += abs(numpy.array(op[:di*dj*dim])).sum() 164 cnt += di*dj*dim 165 if close(v1, vref, cnt, place): 166 print("pass: ", name) 167 else: 168 print("* FAIL: ", name, ". err:", '%.16g' % abs(v1-vref), "/", vref) 169 170def cdouble_to_cmplx(arr): 171 return numpy.array(arr)[0::2] + numpy.array(arr)[1::2] * 1j 172 173def test_int1e_spinor(name, vref, dim, place): 174 intor = getattr(_cint, name) 175 intor.restype = ctypes.c_void_p 176 op = (ctypes.c_double * (20000 * dim))() 177 v1 = 0 178 cnt = 0 179 for j in range(nbas.value*2): 180 for i in range(j+1): 181 di = _cint.CINTlen_spinor(i, c_bas, nbas) * bas[i,NCTR_OF] 182 dj = _cint.CINTlen_spinor(j, c_bas, nbas) * bas[j,NCTR_OF] 183 shls = (ctypes.c_int * 2)(i, j) 184 intor(op, shls, c_atm, natm, c_bas, nbas, c_env) 185 v1 += abs(cdouble_to_cmplx(op[:di*dj*dim*2])).sum() 186 cnt += di*dj*dim*2 187 if close(v1, vref, cnt, place): 188 print("pass: ", name) 189 else: 190 print("* FAIL: ", name, ". err:", '%.16g' % abs(v1-vref), "/", vref) 191 192def max_loc(arr): 193 loc = [] 194 maxi = arr.argmax() 195 n = maxi 196 for i in arr.shape: 197 loc.append(n % i) 198 n /= i 199 loc.reverse() 200 return maxi, loc 201 202def test_comp1e_spinor(name1, name_ref, shift, dim, place): 203 intor = getattr(_cint, name1) 204 intor.restype = ctypes.c_void_p 205 intor_ref = getattr(_cint, name_ref) 206 intor_ref.restype = ctypes.c_void_p 207 op = (ctypes.c_double * (20000 * dim))() 208 op_ref = (ctypes.c_double * (20000 * dim))() 209 210 pfac = 1 211 if shift[0] > 0: 212 pfac *= -1j 213 if shift[1] > 0: 214 pfac *= 1j 215 216 for j in range(nbas.value*2 - shift[1]): 217 for i in range(min(nbas.value*2-shift[0],j+1)): 218 di = _cint.CINTlen_spinor(i, c_bas, nbas) * bas[i,NCTR_OF] 219 dj = _cint.CINTlen_spinor(j, c_bas, nbas) * bas[j,NCTR_OF] 220 shls = (ctypes.c_int * 2)(i+shift[0], j+shift[1]) 221 intor(op, shls, c_atm, natm, c_bas, nbas, c_env) 222 shls_ref = (ctypes.c_int * 2)(i, j) 223 intor_ref(op_ref, shls_ref, c_atm, natm, c_bas, nbas, c_env) 224 dd = abs(pfac * cdouble_to_cmplx(op[:di*dj*dim*2]).reshape(di,dj,dim) 225 - cdouble_to_cmplx(op_ref[:di*dj*dim*2]).reshape(di,dj,dim)) 226 if numpy.round(dd, place).sum(): 227 maxi = dd.argmax() 228 print("* FAIL: ", name1, "/", name_ref, ". shell:", i, j, \ 229 "err:", dd.flatten()[maxi], \ 230 "/", op_ref[maxi*2]+op_ref[maxi*2+1]*1j) 231 return 232 print("pass: ", name1, "/", name_ref) 233 234#################### 235def test_int2e_sph(name, vref, dim, place): 236 intor = getattr(_cint, name) 237 intor.restype = ctypes.c_void_p 238 op = (ctypes.c_double * (1000000 * dim))() 239 v1 = 0 240 cnt = 0 241 for l in range(nbas.value*2): 242 for k in range(l+1): 243 for j in range(nbas.value*2): 244 for i in range(j+1): 245 di = (bas[i,ANG_OF] * 2 + 1) * bas[i,NCTR_OF] 246 dj = (bas[j,ANG_OF] * 2 + 1) * bas[j,NCTR_OF] 247 dk = (bas[k,ANG_OF] * 2 + 1) * bas[k,NCTR_OF] 248 dl = (bas[l,ANG_OF] * 2 + 1) * bas[l,NCTR_OF] 249 shls = (ctypes.c_int * 4)(i, j, k, l) 250 intor(op, shls, c_atm, natm, c_bas, nbas, c_env, opt) 251 v1 += abs(numpy.array(op[:di*dj*dk*dl*dim])).sum() 252 cnt += di*dj*dk*dl*dim 253 if close(v1, vref, cnt, place): 254 print("pass: ", name) 255 else: 256 print("* FAIL: ", name, ". err:", '%.16g' % abs(v1-vref), "/", vref) 257 258def test_erf(name, omega, place): 259 intor = getattr(_cint, name) 260 intor.restype = ctypes.c_void_p 261 env_lr = env.copy() 262 env_lr[PTR_RANGE_OMEGA] = omega 263 env_sr = env.copy() 264 env_sr[PTR_RANGE_OMEGA] = -omega 265 c_env_lr = env_lr.ctypes.data_as(ctypes.c_void_p) 266 c_env_sr = env_sr.ctypes.data_as(ctypes.c_void_p) 267 268 op_lr = numpy.empty(10000) 269 op_sr = numpy.empty(10000) 270 op = numpy.empty(10000) 271 v1 = 0 272 cnt = 0 273 shls = (ctypes.c_int * 4)(3,3,3,3) 274 intor(op.ctypes.data_as(ctypes.c_void_p), shls, c_atm, natm, c_bas, nbas, c_env, opt) 275 intor(op_lr.ctypes.data_as(ctypes.c_void_p), shls, c_atm, natm, c_bas, nbas, c_env_lr, opt) 276 intor(op_sr.ctypes.data_as(ctypes.c_void_p), shls, c_atm, natm, c_bas, nbas, c_env_sr, opt) 277 278 for l in range(nbas.value): 279 for k in range(l+1): 280 for j in range(l+1): 281 for i in range(j+1): 282 di = (bas[i,ANG_OF] * 2 + 1) * bas[i,NCTR_OF] 283 dj = (bas[j,ANG_OF] * 2 + 1) * bas[j,NCTR_OF] 284 dk = (bas[k,ANG_OF] * 2 + 1) * bas[k,NCTR_OF] 285 dl = (bas[l,ANG_OF] * 2 + 1) * bas[l,NCTR_OF] 286 shls = (ctypes.c_int * 4)(i, j, k, l) 287 intor(op.ctypes.data_as(ctypes.c_void_p), shls, c_atm, natm, c_bas, nbas, c_env, opt) 288 intor(op_lr.ctypes.data_as(ctypes.c_void_p), shls, c_atm, natm, c_bas, nbas, c_env_lr, opt) 289 intor(op_sr.ctypes.data_as(ctypes.c_void_p), shls, c_atm, natm, c_bas, nbas, c_env_sr, opt) 290 cnt = di * dj * dk * dl 291 dd = abs(op_lr[:cnt] + op_sr[:cnt] - op[:cnt]) 292 293 with numpy.errstate(invalid='ignore'): 294 dd[dd > 1e-6] / op[:cnt][dd > 1e-6] 295 if numpy.round(dd.max(), place) > 0: 296 maxi = dd.argmax() 297 print("* FAIL: erf+erfc", name, "omega=", omega, 298 " shell:", i, j, k, l, 299 "erf:", op_lr[maxi], "erfc:", op_sr[maxi], 300 "regular:", op[maxi], "err:", dd[maxi]) 301 return 302 print("pass: erf+erfc", name, "omega=", omega) 303 304def test_int2e_spinor(name, vref, dim, place): 305 intor = getattr(_cint, name) 306 intor.restype = ctypes.c_void_p 307 op = (ctypes.c_double * (2000000 * dim))() 308 v1 = 0 309 cnt = 0 310 for l in range(nbas.value*2): 311 for k in range(l+1): 312 for j in range(nbas.value*2): 313 for i in range(j+1): 314 di = _cint.CINTlen_spinor(i, c_bas, nbas) * bas[i,NCTR_OF] 315 dj = _cint.CINTlen_spinor(j, c_bas, nbas) * bas[j,NCTR_OF] 316 dk = _cint.CINTlen_spinor(k, c_bas, nbas) * bas[k,NCTR_OF] 317 dl = _cint.CINTlen_spinor(l, c_bas, nbas) * bas[l,NCTR_OF] 318 shls = (ctypes.c_int * 4)(i, j, k, l) 319 intor(op, shls, c_atm, natm, c_bas, nbas, c_env, opt) 320 v1 += abs(cdouble_to_cmplx(op[:di*dj*dk*dl*dim*2])).sum() 321 cnt += di*dj*dk*dl*dim*2 322 if close(v1, vref, cnt, place): 323 print("pass: ", name) 324 else: 325 print("* FAIL: ", name, ". err:", '%.16g' % abs(v1-vref), "/", vref) 326 327def test_int1e_grids_sph(name, vref, dim, place): 328 intor = getattr(_cint, name) 329 intor.restype = ctypes.c_void_p 330 ngrids = 141 331 numpy.random.seed(12) 332 grids = numpy.random.random((ngrids, 3)) - 5.2 333 env_g = numpy.append(env, grids.ravel()) 334 env_g[NGRIDS] = ngrids 335 env_g[PTR_GRIDS] = env.size 336 op = (ctypes.c_double * (1000000 * dim))() 337 v1 = 0 338 cnt = 0 339 for j in range(nbas.value*2): 340 for i in range(j+1): 341 di = (bas[i,ANG_OF] * 2 + 1) * bas[i,NCTR_OF] 342 dj = (bas[j,ANG_OF] * 2 + 1) * bas[j,NCTR_OF] 343 shls = (ctypes.c_int * 4)(i, j, 0, ngrids) 344 intor(op, shls, c_atm, natm, c_bas, nbas, 345 env_g.ctypes.data_as(ctypes.c_void_p), opt) 346 v1 += abs(numpy.array(op[:ngrids*di*dj*dim])).sum() 347 cnt += ngrids*di*dj*dim 348 if close(v1, vref, cnt, place): 349 print("pass: ", name) 350 else: 351 print("* FAIL: ", name, ". err:", '%.16g' % abs(v1-vref), "/", vref) 352 353def test_int1e_grids_spinor(name, vref, dim, place): 354 intor = getattr(_cint, name) 355 intor.restype = ctypes.c_void_p 356 ngrids = 141 357 numpy.random.seed(12) 358 grids = numpy.random.random((ngrids, 3)) - 5.2 359 env_g = numpy.append(env, grids.ravel()) 360 env_g[NGRIDS] = ngrids 361 env_g[PTR_GRIDS] = env.size 362 op = (ctypes.c_double * (1000000 * dim))() 363 v1 = 0 364 cnt = 0 365 for j in range(nbas.value*2): 366 for i in range(j+1): 367 di = _cint.CINTlen_spinor(i, c_bas, nbas) * bas[i,NCTR_OF] 368 dj = _cint.CINTlen_spinor(j, c_bas, nbas) * bas[j,NCTR_OF] 369 shls = (ctypes.c_int * 4)(i, j, 0, ngrids) 370 intor(op, shls, c_atm, natm, c_bas, nbas, c_env, opt) 371 v1 += abs(cdouble_to_cmplx(op[:ngrids*di*dj*dim*2])).sum() 372 cnt += ngrids*di*dj*dim*2 373 if close(v1, vref, cnt, place): 374 print("pass: ", name) 375 else: 376 print("* FAIL: ", name, ". err:", '%.16g' % abs(v1-vref), "/", vref) 377 378def test_comp2e_spinor(name1, name_ref, shift, dim, place): 379 intor = getattr(_cint, name1) 380 intor.restype = ctypes.c_void_p 381 intor_ref = getattr(_cint, name_ref) 382 intor_ref.restype = ctypes.c_void_p 383 op = (ctypes.c_double * (2000000 * dim))() 384 op_ref = (ctypes.c_double * (2000000 * dim))() 385 386 pfac = 1 387 if shift[0] > 0: 388 pfac *= -1j 389 if shift[1] > 0: 390 pfac *= 1j 391 if shift[2] > 0: 392 pfac *= -1j 393 if shift[3] > 0: 394 pfac *= 1j 395 396 for l in range(nbas.value*2 - shift[3]): 397 for k in range(min(nbas.value*2-shift[0],l+1)): 398 for j in range(nbas.value*2 - shift[1]): 399 for i in range(min(nbas.value*2-shift[0],j+1)): 400 di = _cint.CINTlen_spinor(i, c_bas, nbas) * bas[i,NCTR_OF] 401 dj = _cint.CINTlen_spinor(j, c_bas, nbas) * bas[j,NCTR_OF] 402 dk = _cint.CINTlen_spinor(k, c_bas, nbas) * bas[k,NCTR_OF] 403 dl = _cint.CINTlen_spinor(l, c_bas, nbas) * bas[l,NCTR_OF] 404 shls = (ctypes.c_int * 4)(i+shift[0], j+shift[1], k+shift[2], l+shift[3]) 405 intor(op, shls, c_atm, natm, c_bas, nbas, c_env, opt) 406 shls_ref = (ctypes.c_int * 4)(i, j, k, l) 407 intor_ref(op_ref, shls_ref, c_atm, natm, c_bas, nbas, c_env, opt) 408 dd = abs(pfac * cdouble_to_cmplx(op[:di*dj*dk*dl*dim*2]).reshape(di,dj,dk,dl,dim) 409 - cdouble_to_cmplx(op_ref[:di*dj*dk*dl*dim*2]).reshape(di,dj,dk,dl,dim)) 410 if numpy.round(dd, place).sum(): 411 maxi = dd.argmax() 412 print("* FAIL: ", name1, "/", name_ref, ". shell:", i, j, k, l, \ 413 "err:", dd.flatten()[maxi], \ 414 "/", op_ref[maxi*2]+op_ref[maxi*2+1]*1j) 415 return 416 print("pass: ", name1, "/", name_ref) 417 418 419if __name__ == "__main__": 420 if "--high-prec" in sys.argv: 421 def close(v1, vref, count, place): 422 return round(abs(v1-vref), place) == 0 423 424 for f in (('cint1e_ovlp_sph' , 320.9470780962389, 1, 11), 425 ('cint1e_nuc_sph' , 3664.898206036863, 1, 10), 426 ('cint1e_kin_sph' , 887.2525599069498, 1, 11), 427 ('cint1e_ia01p_sph' , 210.475021425001 , 3, 12), 428 ('cint1e_cg_irxp_sph', 3464.41761486531, 3, 10), 429 ('cint1e_giao_irjxp_sph', 2529.89787038728, 3, 10), 430 ('cint1e_igkin_sph' , 107.6417224161130, 3, 11), 431 ('cint1e_igovlp_sph', 37.94968860771099, 3, 12), 432 ('cint1e_ignuc_sph' , 478.5594827282386, 3, 11), 433 ('cint1e_ipovlp_sph', 429.5284222008585, 3, 11), 434 ('cint1e_ipkin_sph' , 1307.395170673386, 3, 10), 435 ('cint1e_ipnuc_sph' , 8358.422626593954, 3, 10), 436 ('cint1e_iprinv_sph', 385.1108471512923, 3, 11), 437 ('cint1e_prinvxp_sph',210.475021425001, 3, 11), 438 ('cint1e_z_sph' , 651.8811101988866, 1, 11), 439 ('cint1e_zz_sph' , 1881.075059037941, 1, 10), 440 ('cint1e_r_sph' , 1803.13043674652 , 3, 10), 441 ('cint1e_rr_sph' , 13379.47937680471, 9, 9), 442 ('cint1e_r2_sph' , 5237.899221349136, 1, 10), 443 ): 444 test_int1e_sph(*f) 445 446 for f in (('cint1e_ovlp' , 284.4528456839759, 1, 11), 447 ('cint1e_nuc' , 3025.766689838620, 1, 10), 448 ('cint1e_gnuc' , 296.6160944673867, 3, 11), 449 ('cint1e_srsr' , 1430.424389624617, 1, 10), 450 ('cint1e_sr' , 240.4064385362524, 1, 11), 451 ('cint1e_srsp' , 1022.805155947573, 1, 10), 452 ('cint1e_spsp' , 1554.251610462129, 1, 10), 453 ('cint1e_sp' , 265.2448605537422, 1, 11), 454 ('cint1e_spspsp' , 1551.856568558924, 1, 10), 455 ('cint1e_spnuc' , 3905.024862120781, 1, 10), 456 ('cint1e_spnucsp' , 20689.33926165072, 1, 9 ), 457 ('cint1e_srnucsr' , 13408.06084488522, 1, 9 ), 458 ('cint1e_cg_sa10sa01', 319.6545034966355, 9, 11), 459 ('cint1e_cg_sa10sp' , 1705.563585675829, 3, 10), 460 ('cint1e_cg_sa10nucsp',16506.04502697362, 3, 9 ), 461 ('cint1e_giao_sa10sa01' , 358.7833729392868, 9, 11), 462 ('cint1e_giao_sa10sp' , 1070.550400465705, 3, 10), 463 ('cint1e_giao_sa10nucsp', 12819.05472701636, 3, 9 ), 464 ('cint1e_govlp' , 23.2674074483772, 3, 12), 465 ('cint1e_sa01sp' , 218.244203172625, 3, 11), 466 ('cint1e_spgsp' , 96.9346217249868, 3, 12), 467 ('cint1e_spgnucsp' , 1659.37670007911, 3, 10), 468 ('cint1e_spgsa01' , 37.8884662927634, 9, 12), 469 ('cint1e_ipovlp' , 153.860148521121, 3, 12), 470 ('cint1e_ipkin' , 497.249399637873, 3, 11), 471 ('cint1e_ipnuc' , 4506.61348255897, 3, 10), 472 ('cint1e_iprinv' , 240.036283917245, 3, 11), 473 ('cint1e_ipspnucsp', 35059.4071347107, 3, 9), 474 ('cint1e_ipsprinvsp',1166.20850563398, 3, 10), 475 ): 476 test_int1e_spinor(*f) 477 478 for f in (# rys_roots for i,j,k,l=3,3,3,3 has round-off error ~ 1e-5 479 ('cint2e_sph' , 56243.88080655417, 1, 8 ), 480 ('cint2e_ip1_sph', 115489.8647398112, 3, 8 ), 481 ): 482 test_int2e_sph(*f) 483 484 for f in (('cint1e_grids_sph', 5262.154357565998, 1, 9), 485 ('cint1e_grids_ip_sph', 3121.659392012421, 1, 9), 486 ): 487 test_int1e_grids_sph(*f) 488 489 for f in (('cint1e_grids_spvsp', 85792.74665645276, 1, 9), 490 ): 491 test_int1e_grids_spinor(*f) 492 493 test_erf('cint2e_sph', 0.2, 9) 494 test_erf('cint2e_sph', 0.5, 9) 495 test_erf('cint2e_sph', 0.8, 9) 496 497 if "--quick" not in sys.argv: 498 for f in (('cint2e_ip1_sph', 115489.8647398112, 3, 8 ), 499 ('cint2e_p1vxp1_sph', 89014.88690617065, 3, 9), 500 ): 501 test_int2e_sph(*f) 502 for f in (('cint2e' , 37737.11178943688, 1, 8), 503 ('cint2e_spsp1' , 221528.4801800184, 1, 8), 504 ('cint2e_spsp1spsp2' , 1391716.710120983, 1, 7), 505 ('cint2e_srsr1' , 178572.7403565846, 1, 8), 506 ('cint2e_srsr1srsr2' , 860883.5184736011, 1, 8), 507 ('cint2e_cg_sa10sp1' , 241519.2142324595, 3, 8), 508 ('cint2e_cg_sa10sp1spsp2' , 1419443.36668730, 3, 7), 509 ('cint2e_giao_sa10sp1' , 153861.923093945, 3, 8), 510 ('cint2e_giao_sa10sp1spsp2', 918284.847389833, 3, 8), 511 ('cint2e_g1' , 3755.251962011825, 3, 10), 512 ('cint2e_spgsp1' , 16626.99245038913, 3, 9 ), 513 ('cint2e_g1spsp2' , 22186.56802347518, 3, 9 ), 514 ('cint2e_spgsp1spsp2' , 107110.2316610047, 3, 8 ), 515 ('cint2e_ip1' , 34912.85673865274, 3, 9 ), 516 ('cint2e_ipspsp1' , 221092.5698965603, 3, 8 ), 517 ('cint2e_ip1spsp2' , 212447.11625827 , 3, 8 ), 518 ('cint2e_ipspsp1spsp2', 1443972.847900875, 3, 7 ), 519 ): 520 test_int2e_spinor(*f) 521 522 test_comp2e_spinor('cint2e_spsp1', 'cint2e', (4,4,0,0), 1, 11) 523 test_comp2e_spinor('cint2e_spsp1spsp2', 'cint2e', (4,4,4,4), 1, 11) 524 test_comp2e_spinor('cint2e_spsp1spsp2', 'cint2e_spsp1', (0,0,4,4), 1, 11) 525 test_comp2e_spinor('cint2e_spgsp1', 'cint2e_g1', (4,4,0,0), 3, 11) 526 test_comp2e_spinor('cint2e_g1spsp2', 'cint2e_g1', (0,0,4,4), 3, 11) 527 test_comp2e_spinor('cint2e_spgsp1spsp2', 'cint2e_g1', (4,4,4,4), 3, 11) 528 test_comp2e_spinor('cint2e_ipspsp1', 'cint2e_ip1', (4,4,0,0), 3, 11) 529 test_comp2e_spinor('cint2e_ip1spsp2', 'cint2e_ip1', (0,0,4,4), 3, 11) 530 test_comp2e_spinor('cint2e_ipspsp1spsp2', 'cint2e_ip1', (4,4,4,4), 3, 11) 531 532 fz = getattr(_cint, 'cint1e_z_sph') 533 fzz = getattr(_cint, 'cint1e_zz_sph') 534 fr = getattr(_cint, 'cint1e_r_sph') 535 fr2 = getattr(_cint, 'cint1e_r2_sph') 536 frr = getattr(_cint, 'cint1e_rr_sph') 537 v1 = 0 538 for j in range(nbas.value*2): 539 for i in range(j+1): 540 di = (bas[i,ANG_OF] * 2 + 1) * bas[i,NCTR_OF] 541 dj = (bas[j,ANG_OF] * 2 + 1) * bas[j,NCTR_OF] 542 opz = numpy.empty((di,dj) , order='F') 543 opzz = numpy.empty((di,dj) , order='F') 544 opr = numpy.empty((di,dj,3), order='F') 545 opr2 = numpy.empty((di,dj) , order='F') 546 oprr = numpy.empty((di,dj,9), order='F') 547 shls = (ctypes.c_int * 2)(i, j) 548 fz ( opz.ctypes.data_as(ctypes.c_void_p), shls, c_atm, natm, c_bas, nbas, c_env) 549 fzz(opzz.ctypes.data_as(ctypes.c_void_p), shls, c_atm, natm, c_bas, nbas, c_env) 550 fr ( opr.ctypes.data_as(ctypes.c_void_p), shls, c_atm, natm, c_bas, nbas, c_env) 551 fr2(opr2.ctypes.data_as(ctypes.c_void_p), shls, c_atm, natm, c_bas, nbas, c_env) 552 frr(oprr.ctypes.data_as(ctypes.c_void_p), shls, c_atm, natm, c_bas, nbas, c_env) 553 v1 = abs(opz-opr[:,:,2]).sum() 554 v1 += abs(opzz-oprr[:,:,8]).sum() 555 v1 += abs(opr2-oprr[:,:,0]-oprr[:,:,4]-oprr[:,:,8]).sum() 556 if round(v1/(di*dj), 13): 557 print("* FAIL: ", i, j, v1) 558