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