1*
2* routines from basisP.F that use blas calls
3*
4* $Id$
5*
6*.....................................................................
7      logical function bas_add_ucnt(basis, tag, l_value, ngen, nprim,
8     $     rex, expnt, coeffs, ldc, stdtag, oshell_is_rel)
9cc AJL/Begin: This is a wrapper call to the updated version of this function
10cc SPIN-POLARISED ECPS
11      implicit none
12c#include "errquit.fh"
13c#include "mafdecls.fh"
14c#include "basdeclsP.fh"
15c#include "nwc_const.fh"
16c#include "basP.fh"
17c#include "bas_exndcf_dec.fh"
18c#include "ecpso_decP.fh"
19c#include "stdio.fh"
20      integer basis ! [input] basis handle
21      character*(*) tag ! [input] tag on which to add contraction
22      character*(*) stdtag ! [input] standard basis set name associated
23*                                    with the tag
24      integer l_value ! [input] type of contraction (s, p, ..., sp)
25      integer ngen ! [input] no. of contractions
26      integer nprim ! [input] no. of primitives
27      integer ldc ! [input] leading dimension of coeffs
28      double precision expnt(nprim) ! [input] exponents
29      double precision coeffs(ldc, 1:*) ! [input] coefficients
30      double precision rex(nprim) ! [input] gaussian R exponents
31*. . . . . . . . . . . . . . . . .          (offset by 2)
32      logical oshell_is_rel ! [input] flag for relativistic shell
33
34      logical  bas_add_ucnt0
35      external bas_add_ucnt0
36
37      bas_add_ucnt =  bas_add_ucnt0(basis, tag, l_value,ngen, nprim,
38     &                rex, expnt, coeffs, ldc, stdtag, oshell_is_rel, 0)
39
40      end
41c
42      logical function bas_add_ucnt0(basis, tag, l_value, ngen, nprim,
43     $     rex, expnt, coeffs, ldc, stdtag, oshell_is_rel, channel)
44cc AJL/End
45      implicit none
46#include "errquit.fh"
47#include "mafdecls.fh"
48#include "basdeclsP.fh"
49#include "nwc_const.fh"
50#include "basP.fh"
51#include "bas_exndcf_dec.fh"
52#include "ecpso_decP.fh"
53#include "stdio.fh"
54      integer basis ! [input] basis handle
55      character*(*) tag ! [input] tag on which to add contraction
56      character*(*) stdtag ! [input] standard basis set name associated
57*                                    with the tag
58      integer l_value ! [input] type of contraction (s, p, ..., sp)
59      integer ngen ! [input] no. of contractions
60      integer nprim ! [input] no. of primitives
61      integer ldc ! [input] leading dimension of coeffs
62      double precision expnt(nprim) ! [input] exponents
63      double precision coeffs(ldc, 1:*) ! [input] coefficients
64      double precision rex(nprim) ! [input] gaussian R exponents
65*. . . . . . . . . . . . . . . . .          (offset by 2)
66      logical oshell_is_rel ! [input] flag for relativistic shell
67c
68      integer size_add ! amount to be added to exndcf array
69      integer ind ! Index into basis function structures
70      integer free ! Free space pointer
71c
72      integer i, itag, jtag, iu_cont, ntags ! Locals
73      integer s_old, s_new, k_old, k_new, h_old, h_new
74c
75      logical oIs_ecp ! flag saying if this is ecp
76      logical oIs_so ! flag saying if this is so potential
77      logical bas_add_utag
78      external bas_add_utag
79c
80      logical bas_check_handle
81      external bas_check_handle
82cc AJL/Begin/SPIN-POLARISED ECPS
83      integer channel           ! [input] ecp channel [both/alpha/beta]
84cc AJL/End
85c
86#include "bas_exndcf_sfn.fh"
87#include "ecpso_sfnP.fh"
88c
89      oIs_ecp = Is_ECP_in(basis)
90      oIs_so  = Is_SO_in(basis)
91c
92*Ul with rexp=0 :      if (oIs_ecp.and.l_value.eq.(-1)) then
93*Ul with rexp=0 :        do i = 1,nprim
94*Ul with rexp=0 :          if (int(rex(i)).eq.0) then
95*Ul with rexp=0 :            write(luout,*)'This version of nwchem does not support ',
96*Ul with rexp=0 :     &          'local potentials with an r-exponent of 0'
97*Ul with rexp=0 :            call util_flush(luout)
98*Ul with rexp=0 :            write(luout,*)'This should be fixed in a future release',
99*Ul with rexp=0 :     &          ' of NWChem'
100*Ul with rexp=0 :            call util_flush(luout)
101*Ul with rexp=0 :            call errquit(' {ecp}bas_input fatal error ',911)
102*Ul with rexp=0 :          endif
103*Ul with rexp=0 :        enddo
104*Ul with rexp=0 :      endif
105c
106c     adds a new general contraction on the specified tag.  If the
107c     tag is not present it will also add that by calling bas_add_utag
108c
109CC AJL/Begin/SPIN POLARISED ECPS
110c      bas_add_ucnt = bas_check_handle(basis,'bas_add_ucnt')
111c      if (.not. bas_add_ucnt) return
112      bas_add_ucnt0 = bas_check_handle(basis,'bas_add_ucnt0')
113      if (.not. bas_add_ucnt0) return
114      ind = basis + BASIS_HANDLE_OFFSET
115c
116c     Make sure that the tag is in the list
117c
118c      bas_add_ucnt = bas_add_utag(basis, tag, stdtag, itag)
119c      if (.not. bas_add_ucnt) return
120      bas_add_ucnt0 = bas_add_utag(basis, tag, stdtag, itag)
121      if (.not. bas_add_ucnt0) return
122CC AJL/End
123c
124c     Update header information about all unique contractions on all
125c     tags.  Free points to next free word in the exndcf
126c
127*old: free = infbs_head(HEAD_NPRIM,ind)+infbs_head(HEAD_NCOEF,ind)+1
128      free = infbs_head(HEAD_EXCFPTR,ind) + 1
129      s_old = exndcf(SZ_exndcf,ind)
130      size_add = nprim*ngen + nprim
131      if (oIs_ecp.or.oIs_so) size_add = size_add + nprim
132      if ((free+size_add-1) .gt. s_old) then
133        h_old = exndcf(H_exndcf,ind)
134        k_old = exndcf(K_exndcf,ind)
135        s_new = free+size_add-1
136        if (.not.ma_alloc_get(
137     &      mt_dbl,s_new,' input for basis heap ',
138     &      h_new, k_new)) then
139          write(LuOut,*)'bas_add_ucnt0: too many prims/coeffs'
140          write(LuOut,*)' allocated size for input is :',
141     &        exndcf(SZ_exndcf,ind)
142          write(LuOut,*)' size requested here         :',
143     &        (free+size_add-1)
144          bas_add_ucnt0 = .false.
145          return
146        endif
147        call dfill(s_new,0.0d00,dbl_mb(k_new),1)
148        exndcf(H_exndcf,ind) = h_new
149        exndcf(K_exndcf,ind) = k_new
150        exndcf(SZ_exndcf,ind) = s_new
151        call dcopy(s_old,dbl_mb(k_old),1,dbl_mb(k_new),1)
152        if (.not.ma_free_heap(h_old)) call errquit
153     &      ('bas_add_ucnt: error freeing old exponents',911,
154     &       BASIS_ERR)
155      endif
156*      if (infbs_head(HEAD_NCONT,ind)+1 .gt. nucont_bsmx) then
157      if (infbs_head(HEAD_NCONT,ind) .gt. nucont_bsmx) then
158         write(LuOut,*) 'bas_add_ucnt0: too many contractions '
159         bas_add_ucnt0 = .false.
160         return
161      endif
162c
163      infbs_head(HEAD_NCONT,ind) = infbs_head(HEAD_NCONT,ind) + 1
164      infbs_head(HEAD_NPRIM,ind) = infbs_head(HEAD_NPRIM,ind) + nprim
165      infbs_head(HEAD_NCOEF,ind) = infbs_head(HEAD_NCOEF,ind) +
166     $     ngen*nprim
167      infbs_head(HEAD_EXCFPTR,ind) =  infbs_head(HEAD_EXCFPTR,ind) +
168     &      size_add
169c
170      ntags = infbs_head(HEAD_NTAGS,ind)
171      if (itag .ne. ntags) then
172         do jtag = ntags, itag+1, -1
173c
174c     Shuffle data+pointers for following tags up one contraction
175c
176            do iu_cont = infbs_tags(TAG_LCONT,jtag,ind),
177     $           infbs_tags(TAG_FCONT,jtag,ind), -1
178               do i = 1, ndbs_ucont
179                  infbs_cont(i,iu_cont+1,ind) =
180     $                 infbs_cont(i,iu_cont,ind)
181               enddo
182            enddo
183c
184c     Increment first and last contractions on following tags
185c
186            infbs_tags(TAG_FCONT,jtag,ind) =
187     $           infbs_tags(TAG_FCONT,jtag,ind) + 1
188            infbs_tags(TAG_LCONT,jtag,ind) =
189     $           infbs_tags(TAG_LCONT,jtag,ind) + 1
190         enddo
191      endif
192c
193c     Increment basis info on this tag
194c
195      infbs_tags(Tag_High_Ang,itag,ind) =
196     &      max(infbs_tags(Tag_High_Ang,itag,ind),abs(l_value))
197      infbs_tags(TAG_NCONT,itag,ind) = infbs_tags(TAG_NCONT,itag,ind)
198     $     + 1
199      infbs_tags(TAG_NPRIM,itag,ind) = infbs_tags(TAG_NPRIM,itag,ind)
200     $     + nprim
201      infbs_tags(TAG_NCOEF,itag,ind) = infbs_tags(TAG_NCOEF,itag,ind)
202     $     + nprim*ngen
203      if (infbs_tags(TAG_FCONT,itag,ind).eq.0) then
204         if (itag .ne. ntags) call errquit
205     $        ('bas_add_ucnt0: tag error', itag, BASIS_ERR)
206         infbs_tags(TAG_FCONT,itag,ind) = infbs_head(HEAD_NCONT,ind)
207         infbs_tags(TAG_LCONT,itag,ind) = infbs_head(HEAD_NCONT,ind)
208      else
209         infbs_tags(TAG_LCONT,itag,ind) =
210     &      infbs_tags(TAG_LCONT,itag,ind) + 1
211      endif
212c
213*. . . . . . . . . . . . . . . . . . . . . ! Index of new contraction
214      iu_cont = infbs_tags(TAG_LCONT,itag,ind)
215c
216      infbs_cont(CONT_TYPE, iu_cont,ind) = l_value
217      infbs_cont(CONT_NPRIM,iu_cont,ind) = nprim
218      infbs_cont(CONT_NGEN, iu_cont,ind) = ngen
219      infbs_cont(CONT_IEXP, iu_cont,ind) = free
220      infbs_cont(CONT_ICFP, iu_cont,ind) = free + nprim
221cc AJL/Begin/SPIN-POLARISED ECPs
222      infbs_cont(CONT_CHANNEL, iu_cont, ind) = channel
223cc AJL/End
224      if (oIs_ecp.or.oIs_so) then
225        infbs_cont(Cont_Irexp, iu_cont, ind) =
226     &      free + nprim + nprim*ngen
227      else
228*. . . . . . . . . . . . . . . . . . ! point to exponents for safety?
229        infbs_cont(Cont_Irexp, iu_cont,ind) = free
230      endif
231      if (oshell_is_rel) then
232        infbs_cont(CONT_RELLS, iu_cont,ind) = 1
233      else
234        infbs_cont(CONT_RELLS, iu_cont,ind) = 0
235      end if
236c
237c     Copy real data over
238c
239      call dcopy(nprim, expnt, 1, dbl_mb(mb_exndcf(free,ind)), 1)
240      free = free + nprim
241      do i = 1, ngen
242         call dcopy
243     &      (nprim, coeffs(1,i), 1, dbl_mb(mb_exndcf(free,ind)), 1)
244         free = free + nprim
245      enddo
246      if (oIs_ecp.or.oIs_so)
247     &    call dcopy(nprim, rex, 1, dbl_mb(mb_exndcf(free,ind)), 1)
248c
249c     Done
250c
251      end
252*.....................................................................
253      logical function bas_num_uce(basisin,nucent)
254      implicit none
255#include "mafdecls.fh"
256#include "nwc_const.fh"
257#include "basP.fh"
258#include "geobasmapP.fh"
259#include "basdeclsP.fh"
260#include "bas_exndcf_dec.fh"
261#include "stdio.fh"
262c::function
263      logical bas_check_handle
264      external bas_check_handle
265c::passed
266      integer basisin, nucent
267c::local
268      integer basis
269c
270#include "bas_exndcf_sfn.fh"
271c
272      bas_num_uce = bas_check_handle(basisin,'bas_getu_coeff')
273      if (.not.bas_num_uce) return
274
275      basis = basisin + BASIS_HANDLE_OFFSET
276c
277      nucent = infbs_head(HEAD_NTAGS,basis)
278c
279      bas_num_uce = .true.
280c
281      return
282      end
283*.....................................................................
284      logical function bas_uce2cnr(basisin,ucenter,ifirst,ilast)
285      implicit none
286#include "mafdecls.fh"
287#include "nwc_const.fh"
288#include "basP.fh"
289#include "geobasmapP.fh"
290#include "basdeclsP.fh"
291#include "bas_exndcf_dec.fh"
292#include "stdio.fh"
293c::function
294      logical bas_check_handle
295      external bas_check_handle
296c::passed
297      integer basisin, ucenter, ifirst, ilast
298c::local
299      integer basis
300c
301#include "bas_exndcf_sfn.fh"
302c
303      bas_uce2cnr = bas_check_handle(basisin,'bas_getu_coeff')
304      if (.not.bas_uce2cnr) return
305
306      basis = basisin + BASIS_HANDLE_OFFSET
307c
308      ifirst = infbs_tags(TAG_FCONT,ucenter,basis)
309      ilast = infbs_tags(TAG_LCONT,ucenter,basis)
310c
311      bas_uce2cnr = .true.
312c
313      return
314      end
315*.....................................................................
316      logical function bas_uce_tag(basisin,ucent,tagout)
317      implicit none
318#include "nwc_const.fh"
319#include "basP.fh"
320#include "stdio.fh"
321#include "mafdecls.fh"
322#include "geobasmapP.fh"
323#include "bas_ibs_dec.fh"
324c::-function
325      logical bas_check_handle
326      external bas_check_handle
327c::-passed
328      integer basisin
329      integer ucent
330      character*(*) tagout
331c::-local
332      integer basis
333      integer len_tagout
334#include "bas_ibs_sfn.fh"
335c
336      bas_uce_tag = bas_check_handle(basisin,'bas_cont_tag')
337      if (.not.bas_uce_tag) return
338
339      basis = basisin + Basis_Handle_Offset
340
341      len_tagout = len(tagout)
342      tagout = bs_tags(ucent,basis)(1:len_tagout)
343      bas_uce_tag = .true.
344      end
345*.....................................................................
346      logical function bas_getu_coeff(basisin,icont,coeff)
347      implicit none
348#include "mafdecls.fh"
349#include "nwc_const.fh"
350#include "basP.fh"
351#include "geobasmapP.fh"
352#include "basdeclsP.fh"
353#include "bas_exndcf_dec.fh"
354#include "stdio.fh"
355c::function
356      logical bas_check_handle
357      external bas_check_handle
358c:blas
359c     dcopy
360c::passed
361      integer basisin, icont
362      double precision coeff(*)
363c::local
364      integer basis, myucont, icontmax
365      integer mycoeffptr, myprim, mygen
366c
367#include "bas_exndcf_sfn.fh"
368c
369      bas_getu_coeff = bas_check_handle(basisin,'bas_getu_coeff')
370      if (.not.bas_getu_coeff) return
371
372      basis = basisin + BASIS_HANDLE_OFFSET
373c
374      icontmax = infbs_head(HEAD_NCONT,basis)
375      myucont  = icont
376c
377      bas_getu_coeff = icont.gt.0.and.icont.le.icontmax
378      if (.not.(bas_getu_coeff)) then
379        write(LuOut,*)' bas_getu_coeff: ERROR '
380        write(LuOut,*)' contraction range for basis is 1:',
381     &         icontmax
382        write(LuOut,*)' information requested for contraction:',icont
383        return
384      endif
385c
386      mycoeffptr = infbs_cont(CONT_ICFP,myucont,basis)
387      myprim  = infbs_cont(CONT_NPRIM,myucont,basis)
388      mygen   = infbs_cont(CONT_NGEN,myucont,basis)
389      call dcopy ((myprim*mygen),
390     &    dbl_mb(mb_exndcf(mycoeffptr,basis)),1,coeff,1)
391c
392      bas_getu_coeff = .true.
393c
394      return
395      end
396*.....................................................................
397      logical function bas_getu_exponent(basisin,icont,exp)
398      implicit none
399#include "mafdecls.fh"
400#include "nwc_const.fh"
401#include "basP.fh"
402#include "geobasmapP.fh"
403#include "basdeclsP.fh"
404#include "bas_exndcf_dec.fh"
405#include "stdio.fh"
406c::function
407      logical bas_check_handle
408      external bas_check_handle
409c:blas
410c     dcopy
411c::passed
412      integer basisin, icont
413      double precision exp(*)
414c::local
415      integer basis, myucont, icontmax
416      integer myprim,myexptr
417c
418#include "bas_exndcf_sfn.fh"
419c
420      bas_getu_exponent =
421     &    bas_check_handle(basisin,'bas_getu_exponent')
422
423      if (.not.bas_getu_exponent) return
424
425      basis = basisin + BASIS_HANDLE_OFFSET
426
427      icontmax = infbs_head(HEAD_NCONT,basis)
428      myucont = icont
429
430      bas_getu_exponent = icont.gt.0.and.icont.le.icontmax
431      if (.not.(bas_getu_exponent)) then
432        write(LuOut,*)' bas_getu_exponent: ERROR '
433        write(LuOut,*)' contraction range for basis is 1:',
434     &         icontmax
435        write(LuOut,*)' information requested for contraction:',icont
436        return
437      endif
438c
439      myexptr = infbs_cont(CONT_IEXP,myucont,basis)
440      myprim  = infbs_cont(CONT_NPRIM,myucont,basis)
441      call dcopy(myprim,dbl_mb(mb_exndcf(myexptr,basis)),1,exp,1)
442c
443      bas_getu_exponent = .true.
444c
445      return
446      end
447*.....................................................................
448      logical function bas_setu_coeff(basisin,icont,coeff,ncoeff)
449      implicit none
450#include "mafdecls.fh"
451#include "nwc_const.fh"
452#include "basP.fh"
453#include "geobasmapP.fh"
454#include "basdeclsP.fh"
455#include "bas_exndcf_dec.fh"
456#include "stdio.fh"
457c::function
458      logical bas_check_handle
459      external bas_check_handle
460c:blas
461c     dcopy
462c::passed
463      integer basisin, icont, ncoeff
464      double precision coeff(ncoeff)
465c::local
466      integer basis, myucont, icontmax
467      integer mycoeffptr, myprim, mygen
468c
469#include "bas_exndcf_sfn.fh"
470c
471      bas_setu_coeff = bas_check_handle(basisin,'bas_setu_coeff')
472      if (.not.bas_setu_coeff) return
473
474      basis = basisin + BASIS_HANDLE_OFFSET
475c
476      icontmax = infbs_head(HEAD_NCONT,basis)
477      myucont  = icont
478c
479      bas_setu_coeff = icont.gt.0.and.icont.le.icontmax
480      if (.not.(bas_setu_coeff)) then
481        write(LuOut,*)' bas_setu_coeff: ERROR '
482        write(LuOut,*)' contraction range for basis is 1:',
483     &         icontmax
484        write(LuOut,*)' information requested for contraction:',icont
485        return
486      endif
487c
488      mycoeffptr = infbs_cont(CONT_ICFP,myucont,basis)
489      myprim  = infbs_cont(CONT_NPRIM,myucont,basis)
490      mygen   = infbs_cont(CONT_NGEN,myucont,basis)
491c
492      bas_setu_coeff = ncoeff .eq. (myprim*mygen)
493      if(.not.bas_setu_coeff) then
494        write(LuOut,*)' bas_setu_coeff: ERROR '
495        write(LuOut,*)' input and stored number of coefficients ',
496     &         '(nprim*ngen) differ '
497        write(LuOut,*)' input  nprim*ngen: ',ncoeff
498        write(LuOut,*)' stored nprim*ngen: ',(myprim*mygen)
499        return
500      endif
501      call dcopy(ncoeff,coeff,1,
502     &    dbl_mb(mb_exndcf(mycoeffptr,basis)),1)
503c
504      bas_setu_coeff = .true.
505c
506      return
507      end
508*.....................................................................
509      logical function bas_setu_exponent(basisin,icont,exp,nexp)
510      implicit none
511#include "nwc_const.fh"
512#include "basP.fh"
513#include "geobasmapP.fh"
514#include "basdeclsP.fh"
515#include "mafdecls.fh"
516#include "bas_exndcf_dec.fh"
517#include "stdio.fh"
518c::function
519      logical bas_check_handle
520      external bas_check_handle
521c:blas
522c     dcopy
523c::passed
524      integer basisin, icont, nexp
525      double precision exp(nexp)
526c::local
527      integer basis, myucont, icontmax
528      integer myprim,myexptr
529c
530#include "bas_exndcf_sfn.fh"
531c
532      bas_setu_exponent =
533     &    bas_check_handle(basisin,'bas_setu_exponent')
534
535      if (.not.bas_setu_exponent) return
536
537      basis = basisin + BASIS_HANDLE_OFFSET
538
539      icontmax = infbs_head(HEAD_NCONT,basis)
540      myucont = icont
541
542      bas_setu_exponent = icont.gt.0.and.icont.le.icontmax
543      if (.not.(bas_setu_exponent)) then
544        write(LuOut,*)' bas_setu_exponent: ERROR '
545        write(LuOut,*)' contraction range for basis is 1:',
546     &         icontmax
547        write(LuOut,*)' information requested for contraction:',icont
548        return
549      endif
550c
551      myexptr = infbs_cont(CONT_IEXP,myucont,basis)
552      myprim  = infbs_cont(CONT_NPRIM,myucont,basis)
553      bas_setu_exponent = myprim.eq.nexp
554      if (.not.bas_setu_exponent) then
555        write(LuOut,*)' bas_setu_exponent: ERROR '
556        write(LuOut,*)' input and stored number of exponents ',
557     &         '(nprim) differ '
558        write(LuOut,*)' input  nprim: ',nexp
559        write(LuOut,*)' stored nprim: ',myprim
560        return
561      endif
562c
563      call dcopy(nexp,exp,1,dbl_mb(mb_exndcf(myexptr,basis)),1)
564c
565      bas_setu_exponent = .true.
566*.....................................................................
567      end
568