1from math import *
2import numpy
3import ctypes
4
5# general contracted DZ basis [3s1p/2s1p] for H2
6#     exponents    contract-coeff
7# S   6.0          0.7               0.4
8#     2.0          0.6               0.3
9#     0.8          0.5               0.2
10# P   0.9          1.
11
12def gto_norm(n, a):
13    # normalization factor of function r^n e^{-a r^2}
14    s = 2**(2*n+3) * factorial(n+1) * (2*a)**(n+1.5) \
15            / (factorial(2*n+2) * sqrt(pi))
16    return sqrt(s)
17
18CHARGE_OF  = 1
19PTR_COORD  = 2
20NUC_MOD_OF = 3
21PTR_ZETA   = 4
22ATM_SLOTS  = 6
23
24ATOM_OF    = 1
25ANG_OF     = 2
26NPRIM_OF   = 3
27NCTR_OF    = 4
28KAPPA_OF   = 5
29PTR_EXP    = 6
30PTR_COEFF  = 7
31BAS_SLOTS  = 8
32
33ptr_env = 20
34atm = []
35bas = []
36env = [0] * ptr_env
37
38i = 0
39#           CHARGE_OF,PTR_COORD
40atm.append([1,        ptr_env,  0, 0, 0, 0])
41#           x  y  z (Bohr)
42env.extend([0, 0, -0.8])
43ptr_env += 3
44atm.append([1,        ptr_env,  0, 0, 0, 0])
45env.extend([0, 0,  0.8])
46ptr_env += 3
47
48# basis for atom #0
49# 3s -> 2s
50env.extend([6., 2., .8])
51env.append(.7 * gto_norm(0, 6.))
52env.append(.6 * gto_norm(0, 2.))
53env.append(.5 * gto_norm(0, .8))
54env.append(.4 * gto_norm(0, 6.))
55env.append(.3 * gto_norm(0, 2.))
56env.append(.2 * gto_norm(0, .8))
57#           ATOM_OF, ANG_OF, NPRIM_OF, NCTR_OF, KAPPA_OF, PTR_EXP, PTR_COEFF
58bas.append([0,       0,      3,        2,       0,        ptr_env, ptr_env+3, 0])
59ptr_env += 9
60env.extend([.9])
61env.append(1. * gto_norm(1, 0.9))
62bas.append([0,       1,      1,        1,       0,        ptr_env, ptr_env+1, 0])
63ptr_env += 2
64
65# basis functions for atom #1, they are the same to the parameters for atom #0
66bas.extend(bas[-2:])
67
68# note the integer type
69atm = numpy.array(atm, dtype=numpy.intc)
70bas = numpy.array(bas, dtype=numpy.intc)
71env = numpy.array(env, dtype=numpy.double)
72
73_cint = ctypes.cdll.LoadLibrary('/path/to/libcint.so')
74
75_cint.CINTcgto_spheric.restype = ctypes.c_int
76_cint.CINTcgto_spheric.argtypes = [ctypes.c_int, numpy.ctypeslib.ndpointer(dtype=numpy.intc, ndim=2)]
77di = _cint.CINTcgto_spheric(0, bas)
78dj = _cint.CINTcgto_spheric(1, bas)
79
80_cint.cint1e_ipnuc_sph.argtypes = [
81    numpy.ctypeslib.ndpointer(dtype=numpy.double, ndim=3),
82    (ctypes.c_int * 2),
83    numpy.ctypeslib.ndpointer(dtype=numpy.intc, ndim=2),
84    ctypes.c_int,
85    numpy.ctypeslib.ndpointer(dtype=numpy.intc, ndim=2),
86    ctypes.c_int,
87    numpy.ctypeslib.ndpointer(dtype=numpy.double, ndim=1)
88]
89buf = numpy.empty((di, dj, 3), order='F')
90_cint.cint1e_ipnuc_sph(buf, (ctypes.c_int * 2)(0, 1), atm, atm.shape[0], bas, bas.shape[1], env)
91print(buf)
92