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_LIGHT_SPEED    = 0
20PTR_COMMON_ORIG    = 1
21PTR_SHIELDING_ORIG = 4
22PTR_RINV_ORIG      = 4
23PTR_RINV_ZETA      = 7
24PTR_ENV_START      = 20
25
26CHARGE_OF  = 0
27PTR_COORD  = 1
28NUC_MOD_OF = 2
29PTR_ZETA   = 3
30RAD_GRIDS  = 4
31ANG_GRIDS  = 5
32ATM_SLOTS  = 6
33
34ATOM_OF   = 0
35ANG_OF    = 1
36NPRIM_OF  = 2
37NCTR_OF   = 3
38KAPPA_OF  = 4
39PTR_EXP   = 5
40PTR_COEFF = 6
41BAS_SLOTS = 8
42
43natm = 4
44nbas = 0
45atm = numpy.zeros((natm,ATM_SLOTS), dtype=numpy.int32)
46bas = numpy.zeros((1000,BAS_SLOTS), dtype=numpy.int32)
47env = numpy.zeros(10000)
48off = PTR_ENV_START
49for i in range(natm):
50    atm[i, CHARGE_OF] = (i+1)*2
51    atm[i, PTR_COORD] = off
52    env[off+0] = .2 * (i+1)
53    env[off+1] = .3 + (i+1) * .5
54    env[off+2] = .1 - (i+1) * .5
55    off += 3
56off0 = off
57
58# basis with kappa > 0
59nh = 0
60
61bas[nh,ATOM_OF ]  = 0
62bas[nh,ANG_OF  ]  = 1
63bas[nh,KAPPA_OF]  = 1
64bas[nh,NPRIM_OF]  = 1
65bas[nh,NCTR_OF ]  = 1
66bas[nh,PTR_EXP]   = off
67env[off+0] = 1
68bas[nh,PTR_COEFF] = off + 1
69env[off+1] = 1
70off += 2
71nh += 1
72
73bas[nh,ATOM_OF ]  = 1
74bas[nh,ANG_OF  ]  = 2
75bas[nh,KAPPA_OF]  = 2
76bas[nh,NPRIM_OF]  = 2
77bas[nh,NCTR_OF ]  = 2
78bas[nh,PTR_EXP]   = off
79env[off+0] = 5
80env[off+1] = 3
81bas[nh,PTR_COEFF] = off + 2
82env[off+2] = 1
83env[off+3] = 2
84env[off+4] = 4
85env[off+5] = 1
86off += 6
87nh += 1
88
89bas[nh,ATOM_OF ]  = 2
90bas[nh,ANG_OF  ]  = 3
91bas[nh,KAPPA_OF]  = 3
92bas[nh,NPRIM_OF]  = 1
93bas[nh,NCTR_OF ]  = 1
94bas[nh,PTR_EXP ]  = off
95env[off+0] = 1
96bas[nh,PTR_COEFF] = off + 1
97env[off+1] = 1
98off += 2
99nh += 1
100
101bas[nh,ATOM_OF ]  = 3
102bas[nh,ANG_OF  ]  = 4
103bas[nh,KAPPA_OF]  = 4
104bas[nh,NPRIM_OF]  = 1
105bas[nh,NCTR_OF ]  = 1
106bas[nh,PTR_EXP ]  = off
107env[off+0] = .5
108bas[nh,PTR_COEFF] = off + 1
109env[off+1] = 1.
110off = off + 2
111nh += 1
112
113nbas = nh
114
115# basis with kappa < 0
116n = off - off0
117for i in range(n):
118    env[off+i] = env[off0+i]
119
120for i in range(nh):
121        bas[i+nh,ATOM_OF ] = bas[i,ATOM_OF ]
122        bas[i+nh,ANG_OF  ] = bas[i,ANG_OF  ] - 1
123        bas[i+nh,KAPPA_OF] =-bas[i,KAPPA_OF]
124        bas[i+nh,NPRIM_OF] = bas[i,NPRIM_OF]
125        bas[i+nh,NCTR_OF ] = bas[i,NCTR_OF ]
126        bas[i+nh,PTR_EXP ] = bas[i,PTR_EXP ]  + n
127        bas[i+nh,PTR_COEFF]= bas[i,PTR_COEFF] + n
128        env[bas[i+nh,PTR_COEFF]] /= 2 * env[bas[i,PTR_EXP]]
129
130env[bas[5,PTR_COEFF]+0] = env[bas[1,PTR_COEFF]+0] / (2 * env[bas[1,PTR_EXP]+0])
131env[bas[5,PTR_COEFF]+1] = env[bas[1,PTR_COEFF]+1] / (2 * env[bas[1,PTR_EXP]+1])
132env[bas[5,PTR_COEFF]+2] = env[bas[1,PTR_COEFF]+2] / (2 * env[bas[1,PTR_EXP]+0])
133env[bas[5,PTR_COEFF]+3] = env[bas[1,PTR_COEFF]+3] / (2 * env[bas[1,PTR_EXP]+1])
134
135natm = ctypes.c_int(natm)
136nbas = ctypes.c_int(nbas)
137c_atm = atm.ctypes.data_as(ctypes.c_void_p)
138c_bas = bas.ctypes.data_as(ctypes.c_void_p)
139c_env = env.ctypes.data_as(ctypes.c_void_p)
140
141opt = ctypes.POINTER(ctypes.c_void_p)()
142_cint.CINTlen_spinor.restype = ctypes.c_int
143
144
145from pyscf import gto
146mol = gto.M()
147mol._atm = atm[:natm.value]
148mol._bas = bas[:nbas.value]
149mol._env = env
150coords = mol.atom_coords()
151ao = mol.eval_gto('GTOval_sph', coords)
152
153
154def test_int2c1e_sph():
155    fnpp1 = _cint.cint1e_ipiprinv_sph
156    fnp1p = _cint.cint1e_iprinvip_sph
157    nullptr = ctypes.POINTER(ctypes.c_void_p)()
158    def by_pp(shls, shape):
159        buf = numpy.empty(shape+(9,), order='F')
160        fnpp1(buf.ctypes.data_as(ctypes.c_void_p), (ctypes.c_int*4)(*shls),
161              c_atm, natm, c_bas, nbas, c_env, nullptr)
162        ref = buf[:,:,0] + buf[:,:,4] + buf[:,:,8]
163        fnp1p(buf.ctypes.data_as(ctypes.c_void_p), (ctypes.c_int*4)(*shls),
164              c_atm, natm, c_bas, nbas, c_env, nullptr)
165        ref+=(buf[:,:,0] + buf[:,:,4] + buf[:,:,8])*2
166        shls = (shls[1], shls[0])
167        shape = (shape[1], shape[0]) + (9,)
168        buf = numpy.empty(shape, order='F')
169        fnpp1(buf.ctypes.data_as(ctypes.c_void_p), (ctypes.c_int*4)(*shls),
170              c_atm, natm, c_bas, nbas, c_env, nullptr)
171        ref+= (buf[:,:,0] + buf[:,:,4] + buf[:,:,8]).transpose(1,0)
172        return ref * (-.25/numpy.pi)
173
174    #intor = _cint.cint4c1e_sph
175    ao_loc = mol.ao_loc_nr()
176    for nucid in range(mol.natm):
177        mol.set_rinv_orig(coords[nucid])
178        for j in range(nbas.value):
179            j0 = ao_loc[j]
180            j1 = ao_loc[j+1]
181            for i in range(j+1):
182                di = (bas[i,ANG_OF] * 2 + 1) * bas[i,NCTR_OF]
183                dj = (bas[j,ANG_OF] * 2 + 1) * bas[j,NCTR_OF]
184                shls = (i, j)
185                i0 = ao_loc[i]
186                i1 = ao_loc[i+1]
187                buf = numpy.einsum('i,j->ij', ao[nucid,i0:i1], ao[nucid,j0:j1])
188                ref = by_pp(shls, (di,dj))
189                dd = abs(ref - buf).sum()
190                if dd > 1e-8:
191                    print("* FAIL: cint2c1e", "  shell:", i, j, "err:", dd)
192                    return
193    print('cint1e_ipiprinv_sph cint1e_iprinvip_sph pass')
194
195def test_int4c1e_sph():
196    fnpp1 = _cint.cint2e_ipip1_sph
197    fnp1p = _cint.cint2e_ipvip1_sph
198    nullptr = ctypes.POINTER(ctypes.c_void_p)()
199    def by_pp(shls, shape):
200        buf = numpy.empty(shape+(9,), order='F')
201        fnpp1(buf.ctypes.data_as(ctypes.c_void_p), (ctypes.c_int*4)(*shls),
202              c_atm, natm, c_bas, nbas, c_env, nullptr)
203        ref = buf[:,:,:,:,0] + buf[:,:,:,:,4] + buf[:,:,:,:,8]
204        fnp1p(buf.ctypes.data_as(ctypes.c_void_p), (ctypes.c_int*4)(*shls),
205              c_atm, natm, c_bas, nbas, c_env, nullptr)
206        ref+=(buf[:,:,:,:,0] + buf[:,:,:,:,4] + buf[:,:,:,:,8])*2
207        shls = (shls[1], shls[0]) + shls[2:]
208        shape = (shape[1], shape[0]) + shape[2:] + (9,)
209        buf = numpy.empty(shape, order='F')
210        fnpp1(buf.ctypes.data_as(ctypes.c_void_p), (ctypes.c_int*4)(*shls),
211              c_atm, natm, c_bas, nbas, c_env, nullptr)
212        ref+= (buf[:,:,:,:,0] + buf[:,:,:,:,4] + buf[:,:,:,:,8]).transpose(1,0,2,3)
213        return ref * (-.25/numpy.pi)
214
215    intor = _cint.cint4c1e_sph
216    for l in range(nbas.value):
217        for k in range(l+1):
218            for j in range(nbas.value):
219                for i in range(j+1):
220                    di = (bas[i,ANG_OF] * 2 + 1) * bas[i,NCTR_OF]
221                    dj = (bas[j,ANG_OF] * 2 + 1) * bas[j,NCTR_OF]
222                    dk = (bas[k,ANG_OF] * 2 + 1) * bas[k,NCTR_OF]
223                    dl = (bas[l,ANG_OF] * 2 + 1) * bas[l,NCTR_OF]
224                    shls = (i, j, k, l)
225                    buf = numpy.empty((di,dj,dk,dl), order='F')
226                    intor(buf.ctypes.data_as(ctypes.c_void_p), (ctypes.c_int*4)(*shls),
227                          c_atm, natm, c_bas, nbas, c_env, nullptr)
228                    ref = by_pp(shls, (di,dj,dk,dl))
229                    dd = abs(ref - buf).max()
230                    if dd > 1e-6:
231                        print("* FAIL: cint4c1e", "  shell:", i, j, k, l, "err:", dd)
232                        return
233    print('cint4c1e_sph pass')
234
235test_int2c1e_sph()
236#test_int4c1e_sph()
237