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