1      program testbasis
2c $Id$
3      implicit none
4#include "errquit.fh"
5c
6#include "nwc_const.fh"
7#include "bas.fh"
8#include "rtdb.fh"
9#include "mafdecls.fh"
10#include "geom.fh"
11#include "inp.fh"
12#include "stdio.fh"
13c
14c::friends functions
15      logical bas_ucontinfo
16      logical bas_getu_exponent
17      logical bas_getu_coeff
18      logical bas_setu_exponent
19      logical bas_setu_coeff
20      external bas_ucontinfo
21      external bas_getu_exponent
22      external bas_getu_coeff
23      external bas_setu_exponent
24      external bas_setu_coeff
25c
26      integer rtdb, geom, basis, ngen,nprim, iang
27*      integer ncenters
28      integer sphcart, i,j
29      integer nbf,ibf,iat,icont, type
30      character*256 namename, nametrans
31*      character*16 drivtags(20)
32*      double precision coords(3,20), charge(20)
33      integer lexcf
34      parameter (lexcf=400)
35      double precision exp(lexcf), coeff(lexcf)
36      logical status
37      double precision expnt_new(3), coeff_new(4,3)
38      logical bas_add_ucnt, bas_geobas_build
39      external bas_add_ucnt, bas_geobas_build
40      data expnt_new/1., 2., 3./
41      data coeff_new/-1., -2., -3., 0.0, -4., -5., -6., 0.0, -7., -8.,
42     $     -9., 0.0/
43c
44      if (.not. ma_init(MT_DBL, 400 000, 400 000))
45     &      call errquit('testbasis: ma_init failed',911, MA_ERR)
46
47      status = rtdb_open('testbasis.rtdb','empty', rtdb)
48      if (.not.status) call errquit
49     &      ('testbasis rtdb open failed',911, RTDB_ERR)
50c
51c      write(luout,*)' rtdb handle ', rtdb
52*old      status = bas_321g_load(rtdb)
53c
54      open(unit=luin,file='testbasis.nw',
55     &      form='formatted',
56     &      access='sequential',
57     &      status='old',
58     &      err=99565)
59      call input_parse(rtdb)
60c
61      if(.not.geom_create(geom,'geometry')) then
62        write(LuOut,*)' error getting geometry handle '
63        stop ' error '
64      endif
65c
66      if(.not.geom_rtdb_load(rtdb,geom,'geometry'))
67     &      call errquit('error loading geometry',911, GEOM_ERR)
68      status=geom_print(geom)
69c
70      basis = 0
71      if(.not.bas_create(basis,'xy basis')) then
72        write(LuOut,*)' error getting basis handle '
73        stop ' error '
74      endif
75      basis = 0
76      if(.not.bas_create(basis,'wz basis')) then
77        write(LuOut,*)' error getting second basis handle '
78        stop ' error '
79      endif
80      basis = 0
81      if(.not.bas_create(basis,'ao basis')) then
82        write(LuOut,*)' error getting third basis handle '
83        stop ' error '
84      endif
85c
86c
87      status = bas_rtdb_load(rtdb,geom,basis,'ao basis')
88      status = bas_print(basis)
89      status = gbs_map_print(basis)
90c
91      status = bas_continfo(basis,1,iang,nprim,ngen,sphcart)
92      write(LuOut,*)' user:query: status  cont 1 ',status
93      write(LuOut,*)' user:query: type    cont 1 ',iang
94      write(LuOut,*)' user:query: nprim   cont 1 ',nprim
95      write(LuOut,*)' user:query: ngen    cont 1 ',ngen
96      write(LuOut,*)' user:query: sphcart cont 1 ',sphcart
97      if (nprim*ngen.gt.lexcf) call errquit
98     &      (' lexcf too small ',(nprim*ngen), BASIS_ERR)
99      call dfill(lexcf,0.0d00,exp  ,1)
100      call dfill(lexcf,0.0d00,coeff,1)
101      status = bas_get_exponent(basis,1,exp)
102      status = bas_get_coeff(basis,1,coeff)
103      write(LuOut,*)' user exponenents and coefficients '
104      do 00100 i=1,nprim
105        write(LuOut,*)exp(i),(coeff(i+(j-1)*nprim),j=1,ngen)
10600100 continue
107c
108      status = bas_ucontinfo(basis,1,type,nprim,ngen,sphcart)
109      write(LuOut,*)' unique:query: status  cont 1 ',status
110      write(LuOut,*)' unique:query: type    cont 1 ',type
111      write(LuOut,*)' unique:query: nprim   cont 1 ',nprim
112      if (nprim*ngen.gt.lexcf) call errquit
113     &      (' lexcf too small ',(nprim*ngen), BASIS_ERR)
114      write(LuOut,*)' unique:query: ngen    cont 1 ',ngen
115      write(LuOut,*)' unique:query: sphcart cont 1 ',sphcart
116      call dfill(lexcf,0.0d00,exp  ,1)
117      call dfill(lexcf,0.0d00,coeff,1)
118      status = bas_getu_exponent(basis,1,exp)
119      status = bas_getu_coeff(basis,1,coeff)
120      write(LuOut,*)' unique exponenents and coefficients '
121      do 00200 i=1,nprim
122        write(LuOut,*)exp(i),(coeff(i+(j-1)*nprim),j=1,ngen)
12300200 continue
124c
125      exp(1) = 565.6589
126      coeff(1) = 6.021023
127      status = bas_setu_exponent(basis,1,exp,(nprim+1))
128      status = bas_setu_coeff(basis,1,coeff,(nprim*ngen)+1)
129      status = bas_setu_exponent(basis,1,exp,nprim)
130      status = bas_setu_coeff(basis,1,coeff,nprim*ngen)
131c
132c
133      status = bas_ucontinfo(basis,1,type,nprim,ngen,sphcart)
134      write(LuOut,*)' modified '
135      write(LuOut,*)' unique:query: status  cont 1 ',status
136      write(LuOut,*)' unique:query: type    cont 1 ',type
137      write(LuOut,*)' unique:query: nprim   cont 1 ',nprim
138      write(LuOut,*)' unique:query: ngen    cont 1 ',ngen
139      write(LuOut,*)' unique:query: sphcart cont 1 ',sphcart
140      if (nprim*ngen.gt.lexcf) call errquit
141     &      (' lexcf too small ',(nprim*ngen), BASIS_ERR)
142      status = bas_getu_exponent(basis,1,exp)
143      status = bas_getu_coeff(basis,1,coeff)
144      write(LuOut,*)' exponenents and coefficients '
145      do 00300 i=1,nprim
146        write(LuOut,*)exp(i),(coeff(i+(j-1)*nprim),j=1,ngen)
14700300 continue
148c
149c
150c     Try adding new contractions on an existing center
151c
152      write(LuOut,*) ' adding 3*3 d function on H'
153      if (.not. bas_add_ucnt(basis, 'H', 2, 3, 3, expnt_new,
154     &    expnt_new, coeff_new, 4, 'none', .false.))
155     &    write(LuOut,*) ' basis_add_ucnt failed'
156      if (.not. bas_print(basis)) write(LuOut,*) ' print ?'
157c
158c     Try adding new contractions on a new center
159c
160      write(LuOut,*) ' adding 2*3 g function on Cl'
161      if (.not. bas_add_ucnt(basis, 'Cl', 4, 2, 3, expnt_new,
162     $     expnt_new,coeff_new, 4, 'none', .false.))
163     &    write(LuOut,*) ' basis_add_ucnt failed'
164      if (.not. bas_print(basis)) write(LuOut,*) ' print ?'
165c
166c
167      write(LuOut,'(///,1x,a)')' bas_print_all '
168      status = bas_print_all()
169      call bas_err_info(' testbasis bas_err_info info ')
170c
171      status = bas_high_angular(basis,iang)
172      write(LuOut,*)' high angular momentum ', iang
173      status = bas_version()
174      namename = ' '
175      nametrans = ' '
176      status = bas_name(basis,namename, nametrans)
177      write(LuOut,*)' handle         :  ',basis, status
178      i = inp_strlen(namename)
179      write(LuOut,'(a,a,a)')' name           : <',
180     &       namename(1:i),'>'
181      i = inp_strlen(nametrans)
182      write(LuOut,*)' translated name: <',
183     &       nametrans(1:i),'>'
184c
185      write(LuOut,*)' basis function to ce/cn check '
186      status = bas_numbf(basis,nbf)
187      do 00400 ibf=1,nbf
188        status = bas_bf2ce(basis,ibf,iat)
189        status = bas_bf2cn(basis,ibf,icont)
190        write(LuOut,*)' basis function ',ibf,
191     &            ' is on center ',iat,
192     &         ' and in contraction',icont
19300400 continue
194c
195      status = bas_geobas_build(basis)
196      write(LuOut,*)' basis handle ',basis
197      call bas_geomap_check(basis)
198c
199      status = bas_nbf_cn_max(basis,nbf)
200      write(LuOut,*)' max block of nbf on a shell for basis is:',nbf
201c
202      status = bas_nbf_ce_max(basis,nbf)
203      write(LuOut,*)' max block of nbf on an atom for basis is:',nbf
204c
205      status = bas_ncoef_cn_max(basis,nbf)
206      write(LuOut,*)' max num coefs in any contraction ',nbf
207c
208      status = bas_nprim_cn_max(basis,nbf)
209      write(LuOut,*)' max num prims in any contraction ',nbf
210c
211c
212      stop ' testbasis done '
21399565 continue
214      stop ' error opening testbasis.nw'
215      end
216      subroutine bas_geomap_check(basisin)
217      implicit none
218#include "nwc_const.fh"
219#include "bas.fh"
220#include "basP.fh"
221#include "geobasmapP.fh"
222#include "basdeclsP.fh"
223#include "mafdecls.fh"
224#include "bas_ibs_dec.fh"
225#include "geom.fh"
226#include "stdio.fh"
227c
228      integer basisin
229c
230      logical status
231      integer basis, geom
232      integer ncont, icont, ucont, ucont2
233      integer cent_api, cent_sf, cent_mb
234      integer inat, nat, atom_sf, atom_mb
235      integer bflo_api, bflo_sf, bflo_mb
236      integer bfhi_api, bfhi_sf, bfhi_mb
237      integer cnlo_api, cnlo_sf, cnlo_mb
238      integer cnhi_api, cnhi_sf, cnhi_mb
239c
240#include "bas_ibs_sfn.fh"
241c
242      status = bas_geom(basisin,geom)
243      basis = basisin + BASIS_HANDLE_OFFSET
244      write(LuOut,*)' '
245      write(LuOut,*)' '
246      write(LuOut,*)' '
247      write(LuOut,*)' basis geometry mapping checks '
248      write(LuOut,*)' basis handle is    :',basisin
249      write(LuOut,*)' geometry handle is :',geom
250      write(LuOut,*)' '
251      write(LuOut,*)' current memory heap pointers'
252      write(LuOut,*)' ibs_cn2ucn ma handle :',ibs_cn2ucn(H_ibs,basis)
253      write(LuOut,*)' ibs_cn2ucn ma pointer:',ibs_cn2ucn(K_ibs,basis)
254      write(LuOut,*)' ibs_cn2ucn size      :',ibs_cn2ucn(SZ_ibs,basis)
255      write(LuOut,*)' '
256      write(LuOut,*)' ibs_cn2ce ma handle  :',ibs_cn2ce(H_ibs,basis)
257      write(LuOut,*)' ibs_cn2ce ma pointer :',ibs_cn2ce(K_ibs,basis)
258      write(LuOut,*)' ibs_cn2ce size       :',ibs_cn2ce(SZ_ibs,basis)
259      write(LuOut,*)' '
260      write(LuOut,*)' ibs_ce2uce ma handle :',ibs_ce2uce(H_ibs,basis)
261      write(LuOut,*)' ibs_ce2uce ma pointer:',ibs_ce2uce(K_ibs,basis)
262      write(LuOut,*)' ibs_ce2uce size      :',ibs_ce2uce(SZ_ibs,basis)
263      write(LuOut,*)' '
264      write(LuOut,*)' ibs_cn2bfr ma handle :',ibs_cn2bfr(H_ibs,basis)
265      write(LuOut,*)' ibs_cn2bfr ma pointer:',ibs_cn2bfr(K_ibs,basis)
266      write(LuOut,*)' ibs_cn2bfr size      :',ibs_cn2bfr(SZ_ibs,basis)
267      write(LuOut,*)' '
268      write(LuOut,*)' ibs_ce2cnr ma handle :',ibs_ce2cnr(H_ibs,basis)
269      write(LuOut,*)' ibs_ce2cnr ma pointer:',ibs_ce2cnr(K_ibs,basis)
270      write(LuOut,*)' ibs_ce2cnr size      :',ibs_ce2cnr(SZ_ibs,basis)
271      write(LuOut,*)' '
272      write(LuOut,*)' '
273      write(LuOut,*)' '
274c
275      status = bas_numcont(basisin,ncont)
276      status = geom_ncent(geom,nat)
277      write(LuOut,*)' total number of contractions is ',ncont
278      write(LuOut,*)' total number of centers      is ',nat
279      write(LuOut,*)' '
280      write(LuOut,*)' '
281      write(LuOut,*)' '
282      write(LuOut,*)' cont 2 uniq cont check '
283      do icont = 0,ncont
284        ucont = 23234
285        ucont2 = ucont
286        ucont = sf_ibs_cn2ucn(icont,basis)
287        ucont2 = int_mb(mb_ibs_cn2ucn(icont,basis))
288        write(LuOut,'(a,i5,a,i5,i5)')
289     &        ' contraction ',icont,
290     &        ' maps to unique contraction ',ucont,ucont2
291      enddo
292
293      write(LuOut,*)' '
294      write(LuOut,*)' '
295      write(LuOut,*)' cont 2 center check '
296      do icont = 0,ncont
297        status = bas_cn2ce(basisin,icont,cent_api)
298        cent_sf = sf_ibs_cn2ce(icont,basis)
299        cent_mb = int_mb(mb_ibs_cn2ce(icont,basis))
300        write(LuOut,'(a,i5,a,3i5,a)')
301     &        ' contraction ',icont,' is on center ',
302     &        cent_api,cent_sf,cent_mb,
303     &        ' ... (sould be 3 identical numbers)'
304      enddo
305c
306      write(LuOut,*)' '
307      write(LuOut,*)' '
308      write(LuOut,*)' center 2 unique center check '
309      do inat = 1,nat
310        atom_sf = sf_ibs_ce2uce(inat,basis)
311        atom_mb = int_mb(mb_ibs_ce2uce(inat,basis))
312        write(LuOut,'(a,i5,a,2i5,a)')
313     &        ' center ',inat,
314     &        ' maps to unique center/tag ',atom_sf,atom_mb,
315     &        ' ... (sould be 2 identical numbers)'
316      enddo
317c
318      write(LuOut,*)' '
319      write(LuOut,*)' '
320      write(LuOut,*)' contraction 2 basis function range '
321      do icont = 1,ncont
322        status = bas_cn2bfr(basisin,icont,bflo_api,bfhi_api)
323        bflo_sf = sf_ibs_cn2bfr(1,icont,basis)
324        bfhi_sf = sf_ibs_cn2bfr(2,icont,basis)
325        bflo_mb = int_mb(mb_ibs_cn2bfr(1,icont,basis))
326        bfhi_mb = int_mb(mb_ibs_cn2bfr(2,icont,basis))
327        write(LuOut,'(a,i5,a,6i5,a)')
328     &        ' contraction ',icont,
329     &        ' maps to range of basis functions ',
330     &        bflo_api,bfhi_api,bflo_sf,bfhi_sf,bflo_mb,bfhi_mb,
331     &        ' ... (sould be 3 sets identical numbers)'
332      enddo
333c
334      write(LuOut,*)' '
335      write(LuOut,*)' '
336      write(LuOut,*)' center 2 contraction range '
337      do inat = 1,nat
338        status = bas_ce2cnr(basisin,inat,cnlo_api,cnhi_api)
339        cnlo_sf = sf_ibs_ce2cnr(1,inat,basis)
340        cnhi_sf = sf_ibs_ce2cnr(2,inat,basis)
341        cnlo_mb = int_mb(mb_ibs_ce2cnr(1,inat,basis))
342        cnhi_mb = int_mb(mb_ibs_ce2cnr(2,inat,basis))
343        write(LuOut,'(a,i5,a,6i5,a)')
344     &        ' center ',inat,
345     &        ' maps to range of contractons',
346     &        cnlo_api,cnhi_api,cnlo_sf,cnhi_sf,cnlo_mb,cnhi_mb,
347     &        ' ... (sould be 3 sets identical numbers)'
348      enddo
349c
350      end
351